BOSS_7.0.2: EmcRecShowerPosLinShMax Class Reference

void EmcRecShowerPosLinShMax::Position ( RecEmcShower & aShower ) [virtual]

Definition at line 18 of file EmcRecShowerPosLinShMax.cxx.

References EmcID::barrel_ec(), EmcRecParameter::BarrShLinPhiPara(), EmcRecParameter::BarrShLinThetaPara(), RecEmcShower::Begin(), cos(), EmcID::crystal_id(), DstEmcShower::e5x5(), RecEmcShower::End(), DstEmcShower::energy(), IEmcRecGeoSvc::GetBarrelR(), IEmcRecGeoSvc::GetCFrontCenter(), IEmcRecGeoSvc::GetCrystalPoint(), RecEmcShower::getEAll(), EmcRecParameter::GetInstance(), RecEmcShower::getShowerId(), genRecEmupikp::i, ganga-rec::j, EmcID::phi_module(), pi, boss::pos, EmcRecParameter::PosCorr(), DstEmcShower::setDphi(), DstEmcShower::setDtheta(), DstEmcShower::setErrorMatrix(), DstEmcShower::setPosition(), EmcRecParameter::SigPhi(), EmcRecParameter::SigTheta(), sin(), tan(), EmcID::theta_module(), and w.

00019 {
00020   RecEmcFractionMap::const_iterator cit;
00021   HepPoint3D pos(0,0,0);
00022   HepPoint3D posMax(0,0,0);
00023   HepPoint3D possum(0,0,0);
00024   double w,wsum=0;
00025   double etot;
00026  
00027 //cout<<"EmcRecShowerPosLinShMax::Position Begin"<<endl;
00028  
00029   IEmcRecGeoSvc* iGeoSvc;
00030   ISvcLocator* svcLocator = Gaudi::svcLocator();
00031   StatusCode sc = svcLocator->service("EmcRecGeoSvc",iGeoSvc);
00032   if(sc!=StatusCode::SUCCESS) {
00033     //cout<<"Error: Can't get EmcRecGeoSvc"<<endl;
00034   } 
00035 
00036   EmcRecParameter& Para=EmcRecParameter::GetInstance();   
00037 
00038   for(cit=aShower.Begin();cit!=aShower.End();cit++){
00039     etot+=(cit->second.getEnergy()*cit->second.getFraction());
00040     //cout<<etot<<endl;
00041   }
00042   wsum=0;
00043 
00044     double etota=aShower.getEAll();
00045     double ltot=(log(etota/0.0145)+0.5)*1.86;
00046     double e55=aShower.e5x5();
00047     double len55=(log(etota/0.0145)+0.5)*1.86;
00048     Identifier   seedID = aShower.getShowerId();
00049     HepPoint3D seedPoint(iGeoSvc->GetCFrontCenter(seedID));
00050     //Rs is the distance from orgin to showerCentroid, unit cm
00051 //    HepDouble Rs=seedPoint.mag();   
00052   for(cit=aShower.Begin();cit!=aShower.End();cit++){
00053     w=(cit->second.getEnergy()*cit->second.getFraction());
00054 
00055           HepPoint3D hitPoint(iGeoSvc->GetCFrontCenter(cit->second.getCellId()));
00056 //          HepDouble Ang=seedPoint.angle(hitPoint);
00057           // theDistance is the distance from orgin to showerMax
00058           double theDistance=0;
00059             if(cit->second.getCellId()==seedID) {
00060 //              theDistance=(Rs+ltot);
00061                 theDistance=(seedPoint.mag()+ltot);
00062             } else {
00063               theDistance=(seedPoint.mag()+ltot)/cos(seedPoint.angle(hitPoint));
00064 //              theDistance=(Rs+ltot)/cos(Ang);
00065             }  
00066 
00067   //cout<<w<<endl;
00068     wsum+=w;
00069     pos=iGeoSvc->GetCFrontCenter(cit->second.getCellId());
00070 
00071     posMax = (theDistance/pos.mag())*pos;
00072 
00073     possum+=posMax*w;
00074   }
00075 
00076   if(wsum<=0) {
00077     cout<<"[[[[[[[[[[[[[[["<<wsum;
00078   }
00079   pos=possum/wsum;
00080 //  pos=((pos.mag()-len55)/pos.mag())*pos;
00081 //  aShower.setPosition(pos);
00082   if(Para.PosCorr()==0){        aShower.setPosition(pos);}
00083   
00084   //----------------------------------------------------------------------
00085   //position correction
00086   RecEmcID id = aShower.getShowerId();
00087   const unsigned int module = EmcID::barrel_ec(id);;
00088   const unsigned int thetaModule = EmcID::theta_module(id);
00089   const unsigned int phiModule = EmcID::phi_module(id);
00090   HepPoint3D posCorr(pos);
00091 
00092   if(module==1) {   //barrel
00093     unsigned int theModule;
00094     if(thetaModule>21) {
00095       theModule = 43 - thetaModule;
00096       id = EmcID::crystal_id(module,theModule,phiModule);
00097       posCorr.setZ(-posCorr.z());
00098     } else {
00099       theModule = thetaModule;
00100     }
00101 
00102     double IRShift, parTheta[5],parPhi[5];
00103     if(theModule==21) {
00104       IRShift = 0;
00105     } else if(theModule==20) {
00106       IRShift = 2.5;
00107     } else {
00108       IRShift = 5.;
00109     }
00110 
00111     EmcRecParameter& Para=EmcRecParameter::GetInstance();
00112 
00113     for(int i=0;i<5;i++){
00114 
00115       if(e55>1.0)                     parTheta[i] = Para.BarrShLinThetaPara(theModule,i);
00116       else if(e55<=1.0 && e55>0.5)    parTheta[i] = Para.BarrShLinThetaPara(theModule+22,i);
00117       else if(e55<=0.5)               parTheta[i] = Para.BarrShLinThetaPara(theModule+44,i);
00118 
00119 
00120       if(e55>1.0)                      parPhi[i] = Para.BarrShLinPhiPara(0,i);
00121       else if(e55<=1.0 && e55>0.5)     parPhi[i] = Para.BarrShLinPhiPara(1,i);
00122       else if(e55<=0.5)                parPhi[i] = Para.BarrShLinPhiPara(2,i);
00123       
00124     }
00125 
00126     HepPoint3D IR(0,0,IRShift);
00127     HepPoint3D center01, center23;  //center of point0,1 and point2,3 in crystal
00128     center01 = (iGeoSvc->GetCrystalPoint(id,0)+iGeoSvc->GetCrystalPoint(id,1))/2;
00129     center23 = (iGeoSvc->GetCrystalPoint(id,2)+iGeoSvc->GetCrystalPoint(id,3))/2;
00130 
00131     double theta01,theta23,length,d;
00132     theta01 = (center01-IR).theta();
00133     theta23 = (center23-IR).theta();
00134     length = (center01-IR).mag();
00135     d = length*tan(theta01-theta23);  //length in x direction
00136 
00137     double xIn,xOut,deltaTheta;
00138     xIn = length*tan(theta01-(posCorr-IR).theta())-d/2;
00139     xOut = xIn-( parTheta[0]*atan(xIn*parTheta[1]-parTheta[4])
00140         + parTheta[2]*xIn + parTheta[3] );
00141     deltaTheta = atan((xOut+d/2)/length)-atan((xIn+d/2)/length);
00142     //cout<<"xIn: "<<xIn<<" xOut: "<<xOut<<" deltaTheta: "<<deltaTheta<<endl;
00143     //cout<<"dx: "<<parTheta[0]*atan(xIn*parTheta[1]-parTheta[4])
00144     //      + parTheta[2]*xIn + parTheta[3] <<endl;
00146     double yIn,yOut,deltaPhi;
00147 
00148     if(phiModule==0&&posCorr.phi()<0){
00149     
00150       yIn = posCorr.phi()*180./CLHEP::pi+360.-119*3-1.843;
00151   
00152     }else if (phiModule==119&&posCorr.phi()>0){
00153     
00154       yIn = posCorr.phi()*180./CLHEP::pi-1.843;
00155 
00156     }else {
00157 
00158       yIn = posCorr.phi() < 0 ? posCorr.phi()*180./CLHEP::pi+360.-phiModule*3.-1.843
00159       : posCorr.phi()*180./CLHEP::pi-phiModule*3.-1.843;
00160 
00161     }
00162     yOut = parPhi[0]*atan(yIn*parPhi[1]-parPhi[4])+parPhi[2]*yIn+parPhi[3];
00163 
00164 //    yOut = parPhi[0]+parPhi[1]*yIn+parPhi[2]*yIn*yIn+parPhi[3]*yIn*yIn*yIn+parPhi[4]*yIn*yIn*yIn*yIn+parPhi[5]*yIn*yIn*yIn*yIn*yIn;
00165 
00166     deltaPhi = yOut*CLHEP::pi/180.;
00167 
00168     HepPoint3D rotateVector(0,1,0);   //y axis
00169     rotateVector.rotateZ(center01.phi());
00170     posCorr.setZ(posCorr.z()-IRShift);
00171     posCorr.rotate(-deltaTheta,rotateVector);
00172     posCorr.setZ(posCorr.z()+IRShift);
00173     //phi parameter is gotten by the average of e+ e- peak
00174 //    posCorr.rotateZ(-0.002994);
00175     posCorr.rotateZ(-deltaPhi);
00176 
00177     if(thetaModule>21) {
00178       posCorr.setZ(-posCorr.z());
00179     }
00180   }
00181   posCorr=((posCorr.mag()-len55)/posCorr.mag())*posCorr;
00182   //  aShower.setPosition(posCorr);
00183   if(Para.PosCorr()==1){  aShower.setPosition(posCorr);}
00184   
00186   if(aShower.energy()<1e-5) return;
00187 
00188   double r,theta,phi;
00189   double dtheta,dphi,dx,dy,dz;
00190 
00191   r = posCorr.mag();
00192   theta = posCorr.theta();
00193   phi = posCorr.phi();
00194   //cout<<"energy: "<<aShower.energy()<<" theta: "<<theta<<" phi: "<<phi<<endl;
00195 
00196   if(aShower.energy()>0) {
00197     dtheta = Para.SigTheta(0)/sqrt(aShower.energy())+Para.SigTheta(1);
00198     dphi = Para.SigPhi(0)/sqrt(aShower.energy())+Para.SigPhi(1);
00199   }
00200   else {
00201     dtheta = 0;
00202     dphi = 0;
00203   }
00204   //cout<<"dtheta: "<<dtheta<<" dphi: "<<dphi<<endl;
00205 
00206   dx = fabs(iGeoSvc->GetBarrelR()*sin(phi)*dphi);
00207   dy = fabs(iGeoSvc->GetBarrelR()*cos(phi)*dphi);
00208   if(theta>0)
00209     dz = fabs(iGeoSvc->GetBarrelR()*dtheta/(sin(theta)*sin(theta)));
00210   else
00211     dz = 0;
00212   //cout<<" dx: "<<dx<<" dy: "<<dy<<" dz: "<<dz<<endl;
00213 
00214   aShower.setDtheta(dtheta);
00215   aShower.setDphi(dphi);
00216 
00217   //----------------------------------------------------------------------
00218   // Error matrix
00219   HepSymMatrix matrix(3); //error matrix of r, theta, phi
00220   matrix[0][0]=0;
00221   matrix[1][1]=dtheta*dtheta;
00222   matrix[2][2]=dphi*dphi;
00223   matrix[0][1]=0;
00224   matrix[0][2]=0;
00225   matrix[1][2]=0;
00226   //cout<<"matrix:\n"<<matrix<<endl;
00227 
00228   HepMatrix matrix1(3,3),matrix2(3,3);
00229   matrix1[0][0]=sin(theta)*cos(phi);
00230   matrix1[0][1]=r*cos(theta)*cos(phi);
00231   matrix1[0][2]=-r*sin(theta)*sin(phi);
00232   matrix1[1][0]=sin(theta)*sin(phi);
00233   matrix1[1][1]=r*cos(theta)*sin(phi);
00234   matrix1[1][2]=r*sin(theta)*cos(phi);
00235   matrix1[2][0]=cos(theta);
00236   matrix1[2][1]=-r*sin(theta);
00237   matrix1[2][2]=0;
00238   //cout<<"matrix1:\n"<<matrix1<<endl;
00239 
00240   for(int i=0;i<3;i++) {
00241     for(int j=0;j<3;j++) {
00242       matrix2[i][j]=matrix1[j][i];
00243     }
00244   }
00245   //cout<<"matrix2:\n"<<matrix2<<endl;
00246 
00247   HepMatrix matrix4(3,3);
00248   matrix4=matrix1*matrix*matrix2;
00249   //cout<<"matrix4:\n"<<matrix4<<endl;
00250 
00251   HepSymMatrix matrix5(3); //error matrix of x, y, z
00252   for(int i=0;i<3;i++) {
00253     for(int j=0;j<=i;j++) {
00254       matrix5[i][j]=matrix4[i][j];
00255     }
00256   }
00257   //cout<<"matrix5:\n"<<matrix5<<endl;
00258 
00259   aShower.setErrorMatrix(matrix5);
00260 
00261   if(matrix5[0][0]>0) dx=sqrt(matrix5[0][0]);
00262   if(matrix5[1][1]>0) dy=sqrt(matrix5[1][1]);
00263   if(matrix5[2][2]>0) dz=sqrt(matrix5[2][2]);
00264 }


Public Member Functions
	EmcRecShowerPosLinShMax ()
virtual	~EmcRecShowerPosLinShMax ()
virtual void	Position (RecEmcShower &aShower)

EmcRecShowerPosLinShMax Class Reference

Public Member Functions

Detailed Description

Constructor & Destructor Documentation

Member Function Documentation