BOSS_7.0.2: EmcRecShowerPosLin Class Reference

void EmcRecShowerPosLin::Position ( RecEmcShower & aShower ) [virtual]
Definition at line 17 of file EmcRecShowerPosLin.cxx.
References EmcID::barrel_ec(), EmcRecParameter::BarrLinPhiPara(), EmcRecParameter::BarrLinThetaPara(), RecEmcShower::Begin(), cos(), EmcID::crystal_id(), DstEmcShower::e5x5(), EmcRecParameter::EastLinPhiPara(), EmcRecParameter::EastLinThetaPara(), RecEmcShower::End(), DstEmcShower::energy(), IEmcRecGeoSvc::GetBarrelR(), IEmcRecGeoSvc::GetCFrontCenter(), IEmcRecGeoSvc::GetCrystalPoint(), EmcRecParameter::GetInstance(), RecEmcShower::getShowerId(), genRecEmupikp::i, ganga-rec::j, EmcRecParameter::MethodMode(), 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(), w, EmcRecParameter::WestLinPhiPara(), and EmcRecParameter::WestLinThetaPara().
00018 {
00019   RecEmcFractionMap::const_iterator cit;
00020   HepPoint3D pos(0,0,0);
00021   HepPoint3D possum(0,0,0);
00022   double w,wsum;
00023   double etot;
00024   //cout<<"EmcRecShowerPosLin::Position Begin"<<endl; 
00025   
00026   IEmcRecGeoSvc* iGeoSvc;
00027   ISvcLocator* svcLocator = Gaudi::svcLocator();
00028   StatusCode sc = svcLocator->service("EmcRecGeoSvc",iGeoSvc);
00029   if(sc!=StatusCode::SUCCESS) {
00030 //    cout<<"Error: Can't get EmcRecGeoSvc"<<endl;
00031   }
00032 
00033   EmcRecParameter& Para=EmcRecParameter::GetInstance();
00034 
00035   double e5x5 = aShower.e5x5();
00036   
00037   etot=0;
00038 
00039   for(cit=aShower.Begin();cit!=aShower.End();cit++){      
00040     etot+=(cit->second.getEnergy()*cit->second.getFraction());
00041   }
00042   wsum=0;
00043   for(cit=aShower.Begin();cit!=aShower.End();cit++){
00044     w=(cit->second.getEnergy()*cit->second.getFraction());
00045     wsum+=w;
00046     pos=iGeoSvc->GetCFrontCenter(cit->second.getCellId());
00047     possum+=pos*w;
00048   }
00049 
00050   if(wsum<=0) {
00051     cout<<"[[[[[[[[[[[[[[["<<wsum;
00052   }
00053 
00054   pos=possum/wsum;
00055 //  aShower.setPosition(pos);
00056 
00057   //PosCorr=0 without position correction
00058   //PosCorr=1 with position correction
00059   if(Para.PosCorr()==0){        aShower.setPosition(pos);}
00060 
00061   
00062   //----------------------------------------------------------------------
00063   //position correction
00064   RecEmcID id = aShower.getShowerId();
00065   const unsigned int module = EmcID::barrel_ec(id);;
00066   const unsigned int thetaModule = EmcID::theta_module(id);
00067   const unsigned int phiModule = EmcID::phi_module(id);
00068   HepPoint3D posCorr(pos);
00069 
00070   if(module==1) {   //barrel
00071     unsigned int theModule;
00072     if(thetaModule>21) {
00073       theModule = 43 - thetaModule;
00074       id = EmcID::crystal_id(module,theModule,phiModule);
00075       posCorr.setZ(-posCorr.z());
00076     } else {
00077       theModule = thetaModule;
00078     }
00079     //cout<<"EmcID: "<<id<<" theta: "<<theModule<<endl;
00080 
00081     EmcRecParameter& Para=EmcRecParameter::GetInstance();
00082     double IRShift, parTheta[5],parPhi[5];
00083 
00084 
00085     if(theModule==21) {
00086       IRShift = 0;
00087     } else if(theModule==20) {
00088       IRShift = 2.5;
00089     } else {
00090       IRShift = 5.;
00091     }
00092         
00093     
00094     if(Para.MethodMode()==0){
00095 
00096             
00097       for(int i=0;i<5;i++){
00098         if(theModule==21) {
00099           double par[5] = { 1.698, -1.553, 0.9384, 0.1193, -0.01916 };
00100           parTheta[i] = par[i];
00101         }else if(theModule==20){
00102           double par[5] = { 1.593, -1.582, 0.8881, 0.3298, -0.08856 };
00103           parTheta[i] = par[i];
00104         }else{
00105           double par[5] = { 1.605, -1.702, 0.8433, 0.3072, -0.1809 };
00106           parTheta[i] = par[i];
00107         }
00108       }
00109 
00110     } else if(Para.MethodMode()==1){
00111 
00112       for(int i=0;i<5;i++){
00113       
00114         if(e5x5>1.0)                      parTheta[i] = Para.BarrLinThetaPara(theModule,i);
00115         else if(e5x5<=1.0 && e5x5>0.5)    parTheta[i] = Para.BarrLinThetaPara(theModule+22,i);
00116         else if(e5x5<=0.5)                parTheta[i] = Para.BarrLinThetaPara(theModule+44,i);
00117             
00118         if(e5x5>1.0)                      parPhi[i] = Para.BarrLinPhiPara(0,i);
00119         else if(e5x5<=1.0 && e5x5>0.5)    parPhi[i] = Para.BarrLinPhiPara(1,i);
00120         else if(e5x5<=0.5)                parPhi[i] = Para.BarrLinPhiPara(2,i);
00121             
00122       }
00123     }
00124   
00125 
00126 
00127     HepPoint3D IR(0,0,IRShift);
00128     HepPoint3D center01, center23;  //center of point0,1 and point2,3 in crystal
00129     center01 = (iGeoSvc->GetCrystalPoint(id,0)+iGeoSvc->GetCrystalPoint(id,1))/2;
00130     center23 = (iGeoSvc->GetCrystalPoint(id,2)+iGeoSvc->GetCrystalPoint(id,3))/2;
00131     //cout<<"before correction: "<<(posCorr-IR).theta()<<"\t"<<posCorr.phi()<<endl;
00132 
00133     double theta01,theta23,length,d;
00134     theta01 = (center01-IR).theta();
00135     theta23 = (center23-IR).theta();
00136     length = (center01-IR).mag();
00137     d = length*tan(theta01-theta23);  //length in x direction
00138     //cout<<"theta01: "<<theta01<<"\t"<<" theta23: "<<theta23<<"\t"
00139     //<<" length: "<<length<<" d: "<<d<<endl;
00140 
00141     double xIn,xOut,deltaTheta;
00142     xIn = length*tan(theta01-(posCorr-IR).theta())-d/2;
00143     xOut = xIn-(parTheta[0]*atan(xIn*parTheta[1]-parTheta[4])+parTheta[2]*xIn+parTheta[3]);
00144     deltaTheta = atan((xOut+d/2)/length)-atan((xIn+d/2)/length);
00145     //cout<<"xIn: "<<xIn<<" xOut: "<<xOut<<" deltaTheta: "<<deltaTheta<<endl;
00146 
00147     double parPhi1[5] = { 0.7185, -2.539, 0.6759, 0.3472, -0.3108 };  //e-
00148     double parPhi2[5] = { 0.5781, -2.917, 0.5516, -0.5745, 0.3777 };  //e+
00149 
00150 //    double parPhi[5];// = { 0.7723, -3.1, 0.6992, -0.1441, -0.01012 }; //by photon, not used
00151 
00152 //    for(int i=0;i<5;i++) {
00153 //      parPhi[i] = (parPhi1[i]+parPhi2[i])/2;  //average of e- and e+
00154 //    }
00155 
00156     double yIn,yOut,deltaPhi;
00157 //    yIn = posCorr.phi() < 0 ? posCorr.phi()*180./CLHEP::pi+360.-phiModule*3.-1.1537
00158 //      : posCorr.phi()*180./CLHEP::pi-phiModule*3.-1.1537;
00159     if(phiModule==0&&posCorr.phi()<0){
00160 
00161       yIn = posCorr.phi()*180./CLHEP::pi+360.-119*3-1.843;
00162 
00163     }else if (phiModule==119&&posCorr.phi()>0){
00164 
00165       yIn = posCorr.phi()*180./CLHEP::pi-1.843;
00166 
00167     }else {
00168 
00169       yIn = posCorr.phi() < 0 ? posCorr.phi()*180./CLHEP::pi+360.-phiModule*3.-1.843
00170       : posCorr.phi()*180./CLHEP::pi-phiModule*3.-1.843;
00171 
00172     }
00173 
00174     yOut = parPhi[0]*atan(yIn*parPhi[1]-parPhi[4])+parPhi[2]*yIn+parPhi[3]; 
00175     deltaPhi = yOut*CLHEP::pi/180.;
00176     //cout<<"yIn="<<yIn<<"\tyOut="<<yOut<<"\tdeltaPhi="<<deltaPhi<<endl;
00177 
00178     HepPoint3D rotateVector(0,1,0);   //y axis
00179     rotateVector.rotateZ(posCorr.phi());
00180     posCorr.setZ(posCorr.z()-IRShift);
00181     posCorr.rotate(-deltaTheta,rotateVector);
00182     posCorr.setZ(posCorr.z()+IRShift);
00183     if(Para.MethodMode()==0){
00184         posCorr.rotateZ(-0.002994);
00185     }else if(Para.MethodMode()==1){
00186         posCorr.rotateZ(-deltaPhi);    
00187     }
00188 
00189     if(thetaModule>21) {
00190       posCorr.setZ(-posCorr.z());
00191     }
00192   } else {  //endcap
00193 
00194         EmcRecParameter& Para=EmcRecParameter::GetInstance();
00195 
00196     if(Para.MethodMode()==1){
00197 
00198       double IRShift = 0;
00199       HepPoint3D IR(0,0,IRShift);
00200       double theta12=-9999,theta03=-9999,theta23=-9999,theta14=-9999,delta=-9999;
00201       double phi=-9999,phi01=-9999,phi23=-9999,phi12=-9999,phi34=-9999,deltaphi=-9999; 
00202       double xIn,xOut,deltaTheta,yIn,yOut,deltaPhi;        
00203       double posphi=posCorr.phi();
00204       
00205       double parTheta[5],parPhi[5];     
00206 
00207       HepPoint3D point0 ,point1 ,point2 ,point3 ,point4 ;
00208       point0 = iGeoSvc->GetCrystalPoint(id,0);
00209       point1 = iGeoSvc->GetCrystalPoint(id,1); 
00210       point2 = iGeoSvc->GetCrystalPoint(id,2);
00211       point3 = iGeoSvc->GetCrystalPoint(id,3);
00212       point4 = iGeoSvc->GetCrystalPoint(id,4);
00213 
00214       HepPoint3D center,center01, center23, center12,center03,center34,center14;
00215       center01 = (point0+point1)/2;
00216       center23 = (point2+point3)/2;
00217       center12 = (point1+point2)/2;
00218       center34 = (point3+point4)/2;
00219       center03 = (point0+point3)/2;
00220       center14 = (point1+point4)/2;
00221 
00222       phi01 = center01.phi();
00223       phi23 = center23.phi();
00224       phi12 = center12.phi();
00225       phi34 = center34.phi();
00226 
00227 
00228       if(  (thetaModule==1&&((phiModule+3)%4 == 0||(phiModule+2)%4 == 0))
00229            ||(thetaModule==3&&((phiModule+4)%5 == 0||(phiModule+3)%5 == 0||(phiModule+2)%5 == 0)))
00230         {
00231           center12 = center23;
00232           center03 = center14;
00233           
00234           center01 = center12; 
00235           center23 = center34;
00236 
00237           phi01 = phi12;
00238           phi23 = phi34;
00239         }
00240 
00241       if(phiModule==0){
00242         posphi = posphi;
00243       }else if(  ((thetaModule==0||thetaModule==1)&&phiModule==63)
00244                  ||((thetaModule==2||thetaModule==3)&&phiModule==79)
00245                  ||((thetaModule==4||thetaModule==5)&&phiModule==95)
00246                  ){
00247         
00248         posphi = posphi+2*CLHEP::pi;
00249         
00250       }else {
00251         posphi = posphi < 0 ? posphi+2*CLHEP::pi : posphi;      
00252       }
00253 
00254       phi01 = phi01 <0 ? phi01+2*CLHEP::pi : phi01;
00255       phi23 = phi23 <0 ? phi23+2*CLHEP::pi : phi23;
00256       
00257       if(module==0){  //east endcap     
00258 
00259         IRShift = 10;
00260         center = center23;
00261         phi = phi23;
00262         
00263         for(int i=0;i<5;i++) {
00264           parTheta[i] = Para.EastLinThetaPara(thetaModule,i);
00265           parPhi[i] = Para.EastLinPhiPara(0,i);
00266         }
00267 
00268       }else if (module==2){  //west endcap          
00269 
00270         IRShift = -10;
00271         center = center01;   
00272         phi = phi01;
00273 
00274         for(int i=0;i<5;i++){
00275           parTheta[i] = Para.WestLinThetaPara(thetaModule,i);
00276           parPhi[i] = Para.WestLinPhiPara(0,i);
00277         }
00278 
00279       }
00280 
00281      theta12 = (center12 - IR).theta();
00282      theta03    = (center03 - IR).theta();
00283      delta = theta03 - theta12;
00284      xIn = (((posCorr-IR).theta() - theta12) - delta*0.5)/delta;
00285      xOut = parTheta[0]*atan(xIn*parTheta[1]-parTheta[4]) + parTheta[2]*xIn + parTheta[3] ;
00286      deltaTheta = xOut*delta;
00287       
00288      deltaphi = fabs(phi23 - phi01);
00289      yIn = ((posphi - phi) - deltaphi*0.5)/deltaphi;
00290      yOut = parPhi[0]*atan(yIn*parPhi[1]-parPhi[4]) + parPhi[2]*yIn + parPhi[3] ;
00291      deltaPhi = yOut*deltaphi;
00292      
00293      HepPoint3D rotateVector(0,1,0);   //y axis
00294      rotateVector.rotateZ(center.phi());
00295      posCorr.setZ(posCorr.z()-IRShift);
00296      posCorr.rotate(-deltaTheta,rotateVector);
00297      posCorr.setZ(posCorr.z()+IRShift);
00298      posCorr.rotateZ(-deltaPhi);
00299         
00300     }
00301     }  
00302 
00303   //  aShower.setPosition(posCorr);
00304   if(Para.PosCorr()==1){  aShower.setPosition(posCorr);}
00305 
00306   //----------------------------------------------------------------------
00307   if(aShower.energy()<1e-5) return;
00308 
00309   double r,theta,phi;
00310   double dtheta,dphi,dx,dy,dz;
00311 
00312   r = posCorr.mag();
00313   theta = posCorr.theta();
00314   phi = posCorr.phi();
00315   //cout<<"theta: "<<theta<<" phi: "<<phi<<endl;
00316 
00317 //  EmcRecParameter& Para=EmcRecParameter::GetInstance();   
00318   if(aShower.energy()>0) {
00319     dtheta = Para.SigTheta(0)/sqrt(aShower.energy())+Para.SigTheta(1);
00320     dphi = Para.SigPhi(0)/sqrt(aShower.energy())+Para.SigPhi(1);
00321   }
00322   else {
00323     dtheta = 0;
00324     dphi = 0;
00325   }
00326   //cout<<"dtheta: "<<dtheta<<" dphi: "<<dphi<<endl;
00327 
00328   dx = fabs(iGeoSvc->GetBarrelR()*sin(phi)*dphi);
00329   dy = fabs(iGeoSvc->GetBarrelR()*cos(phi)*dphi);
00330   if(theta>0)
00331     dz = fabs(iGeoSvc->GetBarrelR()*dtheta/(sin(theta)*sin(theta)));
00332   else
00333     dz = 0;
00334   //cout<<" dx: "<<dx<<" dy: "<<dy<<" dz: "<<dz<<endl;
00335 
00336   aShower.setDtheta(dtheta);
00337   aShower.setDphi(dphi);
00338 
00339   //----------------------------------------------------------------------
00340   // Error matrix
00341   HepSymMatrix matrix(3); //error matrix of r, theta, phi
00342   matrix[0][0]=0;
00343   matrix[1][1]=dtheta*dtheta;
00344   matrix[2][2]=dphi*dphi;
00345   matrix[0][1]=0;
00346   matrix[0][2]=0;
00347   matrix[1][2]=0;
00348   //cout<<"matrix: \n"<<matrix;
00349 
00350   HepMatrix matrix1(3,3),matrix2(3,3);
00351   matrix1[0][0]=sin(theta)*cos(phi);
00352   matrix1[0][1]=r*cos(theta)*cos(phi);
00353   matrix1[0][2]=-r*sin(theta)*sin(phi);
00354   matrix1[1][0]=sin(theta)*sin(phi);
00355   matrix1[1][1]=r*cos(theta)*sin(phi);
00356   matrix1[1][2]=r*sin(theta)*cos(phi);
00357   matrix1[2][0]=cos(theta);
00358   matrix1[2][1]=-r*sin(theta);
00359   matrix1[2][2]=0;
00360   //cout<<"matrix1: \n"<<matrix1;
00361 
00362   for(int i=0;i<3;i++) {
00363     for(int j=0;j<3;j++) {
00364       matrix2[i][j]=matrix1[j][i];
00365     }
00366   }
00367   //cout<<"matrix2: \n"<<matrix2;
00368 
00369   HepMatrix matrix4(3,3);
00370   matrix4=matrix1*matrix*matrix2;
00371   //cout<<"matrix4: \n"<<matrix4;
00372 
00373   HepSymMatrix matrix5(3); //error matrix of x, y, z
00374   for(int i=0;i<3;i++) {
00375     for(int j=0;j<=i;j++) {
00376       matrix5[i][j]=matrix4[i][j];
00377     }
00378   }
00379   //cout<<"matrix5: \n"<<matrix5;
00380 
00381   aShower.setErrorMatrix(matrix5);
00382 }

Public Member Functions
	EmcRecShowerPosLin ()
virtual	~EmcRecShowerPosLin ()
virtual void	Position (RecEmcShower &aShower)
EmcRecShowerPosLin Class Reference

Public Member Functions

Detailed Description

Constructor & Destructor Documentation

Member Function Documentation
