BOSS_7.0.2: EmcRecShowerPosLog Class Reference

void EmcRecShowerPosLog::Position ( RecEmcShower & aShower ) [virtual]
Definition at line 18 of file EmcRecShowerPosLog.cxx.
References abs, EmcID::barrel_ec(), EmcRecParameter::BarrLogPhiPara(), EmcRecParameter::BarrLogThetaPara(), EmcRecParameter::BarrPosDataCor(), EmcRecParameter::BarrPosMCCor(), RecEmcShower::Begin(), cos(), EmcID::crystal_id(), EmcRecParameter::DataMode(), DstEmcShower::e3x3(), DstEmcShower::e5x5(), EmcRecParameter::EastLogPhiPara(), EmcRecParameter::EastLogThetaPara(), EmcRecParameter::EastPosDataCor(), RecEmcShower::End(), DstEmcShower::energy(), IEmcRecGeoSvc::GetBarrelR(), IEmcRecGeoSvc::GetCFrontCenter(), IEmcRecGeoSvc::GetCrystalPoint(), RecEmcShower::getEAll(), RecEmcShower::getFractionMap3x3(), RecEmcShower::getFractionMap5x5(), EmcRecParameter::GetInstance(), RecEmcShower::getShowerId(), genRecEmupikp::i, ganga-rec::j, EmcRecParameter::LogPosOffset(), EmcRecParameter::MethodMode(), EmcID::phi_module(), pi, boss::pos, EmcRecParameter::PosCorr(), EmcRecParameter::PositionMode2(), DstEmcShower::setDphi(), DstEmcShower::setDtheta(), DstEmcShower::setErrorMatrix(), DstEmcShower::setPosition(), EmcRecParameter::SigPhi(), EmcRecParameter::SigTheta(), sin(), tan(), EmcID::theta_module(), w, EmcRecParameter::WestLogPhiPara(), EmcRecParameter::WestLogThetaPara(), and EmcRecParameter::WestPosDataCor().
00019 {
00020   RecEmcFractionMap::const_iterator cit;
00021   HepPoint3D pos(0,0,0);
00022   HepPoint3D possum(0,0,0);
00023   double w,wsum=0;
00024 
00025   //  cout<<"EmcRecShowerPosLog::Position Begin"<<endl; 
00026 
00027   IEmcRecGeoSvc* iGeoSvc;
00028   ISvcLocator* svcLocator = Gaudi::svcLocator();
00029   StatusCode sc = svcLocator->service("EmcRecGeoSvc",iGeoSvc);
00030   if(sc!=StatusCode::SUCCESS) {
00031     //    cout<<"Error: Can't get EmcRecGeoSvc"<<endl;
00032   } 
00033 
00034   EmcRecParameter& Para=EmcRecParameter::GetInstance();   
00035 
00036   double offset=Para.LogPosOffset(); 
00037   //cout<<offset<<endl;
00038 
00039   double e5x5 = -99999;
00040   if(Para.PositionMode2()=="all") {
00041     double etot=aShower.getEAll();
00042     for(cit=aShower.Begin();cit!=aShower.End();cit++){
00043       w=offset+log((cit->second.getEnergy()/etot)*cit->second.getFraction());
00044       if(w>0){
00045         wsum+=w;
00046         pos=iGeoSvc->GetCFrontCenter(cit->second.getCellId());
00047         possum+=pos*w;
00048       }
00049     }
00050   } else if(Para.PositionMode2()=="3x3") {
00051     double e3x3=aShower.e3x3();
00052     RecEmcFractionMap fracMap3x3=aShower.getFractionMap3x3();
00053     for(cit=fracMap3x3.begin();cit!=fracMap3x3.end();cit++) {
00054       w=offset+log((cit->second.getEnergy()/e3x3)*cit->second.getFraction());
00055       if(w>0){
00056         wsum+=w;
00057         pos=iGeoSvc->GetCFrontCenter(cit->second.getCellId());
00058         possum+=pos*w;
00059       }
00060     }
00061   } else {
00062     e5x5=aShower.e5x5();
00063     RecEmcFractionMap fracMap5x5=aShower.getFractionMap5x5();
00064     for(cit=fracMap5x5.begin();cit!=fracMap5x5.end();cit++) {
00065       w=offset+log((cit->second.getEnergy()/e5x5)*cit->second.getFraction());
00066       if(w>0){
00067         wsum+=w;
00068         pos=iGeoSvc->GetCFrontCenter(cit->second.getCellId());
00069         possum+=pos*w;
00070       }
00071     }
00072   }
00073   if(wsum<=0) {
00074     //    cout<<"[[[[[[[[[[[[[[["<<wsum<<endl;
00075   }
00076   pos=possum/wsum;
00077 
00078   // aShower.setPositionNoCor(pos);
00079   //----------------------------------------------------------------------
00080   //position correction
00081   RecEmcID id = aShower.getShowerId();
00082   const unsigned int module = EmcID::barrel_ec(id);
00083   const unsigned int thetaModule = EmcID::theta_module(id);
00084   const unsigned int phiModule = EmcID::phi_module(id);
00085   HepPoint3D posCorr(pos);
00086 
00087   if(module==1) {   //barrel
00088 
00089 
00090     double IRShift=0, parTheta[5],parPhi[5]; 
00091  
00092     unsigned int theModule;
00093     if(thetaModule>21) {
00094       theModule = 43 - thetaModule;
00095       id = EmcID::crystal_id(module,theModule,phiModule);
00096       posCorr.setZ(-posCorr.z());
00097     } else {
00098       theModule = thetaModule;
00099     }
00100 
00101  
00102     if(theModule==21) {
00103       IRShift = 0;
00104     } else if(theModule==20) {
00105       IRShift = 2.5;
00106     } else {
00107       IRShift = 5.;
00108     }   
00109     HepPoint3D IR(0,0,IRShift);
00110 
00111     //MethodMode=0 old parameters
00112     //MethodMode=1 new parameters
00113     if(Para.MethodMode()==0){ //using old parameter given by hemiao
00114       
00115       for(int i=0;i<5;i++){
00116         if(theModule==21) {
00117           double par[5] = { 51.97, -0.1209, 6.175, -0.1431, -0.005497 };
00118           parTheta[i] = par[i];
00119         }else if(theModule==20){
00120           double par[5] = { 56.31, -0.1135, 6.312, -1.216, -0.02978 };
00121           parTheta[i] = par[i];
00122         }else{
00123           double par[5] = { 10.65, -0.2257, 2.216, -0.1562, -0.05552 };
00124           parTheta[i] = par[i];
00125         }
00126       }
00127       
00128     }  else  if(Para.MethodMode()==1){ //using new parameter given by liuchunxiu
00129       
00130       for(int i=0;i<5;i++){
00131                 
00132         if(e5x5>1.0)                    parTheta[i] = Para.BarrLogThetaPara(theModule,i);
00133         else if(e5x5<=1.0 && e5x5>0.5)  parTheta[i] = Para.BarrLogThetaPara(theModule+22,i);
00134         else if(e5x5<=0.5)              parTheta[i] = Para.BarrLogThetaPara(theModule+44,i);
00135 
00136         if(e5x5>1.0)                    parPhi[i] = Para.BarrLogPhiPara(theModule,i);
00137         else if(e5x5<=1.0 && e5x5>0.5)  parPhi[i] = Para.BarrLogPhiPara(theModule+22,i);
00138         else if(e5x5<=0.5)              parPhi[i] = Para.BarrLogPhiPara(theModule+44,i);
00139 
00140 
00141       }
00142 
00143     }
00144  
00145     HepPoint3D center01, center23;  //center of point0,1 and point2,3 in crystal
00146     center01 = (iGeoSvc->GetCrystalPoint(id,0)+iGeoSvc->GetCrystalPoint(id,1))/2;
00147     center23 = (iGeoSvc->GetCrystalPoint(id,2)+iGeoSvc->GetCrystalPoint(id,3))/2;
00148 
00149     double theta01,theta23,length,d;
00150     theta01 = (center01-IR).theta();
00151     theta23 = (center23-IR).theta();
00152     length = (center01-IR).mag();
00153     d = length*tan(theta01-theta23);  //length in x direction
00154     
00155     double xIn,xOut,deltaTheta;
00156     xIn = length*tan(theta01-(posCorr-IR).theta())-d/2;
00157     xOut = xIn-( parTheta[0]*atan(xIn*parTheta[1]-parTheta[4])+ parTheta[2]*xIn + parTheta[3] );
00158     deltaTheta = atan((xOut+d/2)/length)-atan((xIn+d/2)/length);
00159     
00160     //obtained by photon, not used, just for comparison
00161     //    double parPhi[5] = { 51.1, -0.1499, 7.62, 0.1301, 0.005581 };
00162     double yIn,deltaY,deltaPhi;
00163 
00164     //1.843=1.5+0.343 degree
00165     if(phiModule==0){
00166       yIn = posCorr.phi()*180./CLHEP::pi -phiModule*3. - 1.843;
00167     }else if (phiModule==119){
00168       yIn = posCorr.phi()*180./CLHEP::pi + 360.-phiModule*3. -1.843;
00169     }else {
00170       yIn = posCorr.phi() < 0 ? posCorr.phi()*180./CLHEP::pi+360.-phiModule*3.-1.843 
00171         : posCorr.phi()*180./CLHEP::pi-phiModule*3.-1.843;
00172     }
00173 
00174     //deltaY=yIn-Ymc => Ymc=yIn-deltaY
00175     deltaY = parPhi[0]*atan(yIn*parPhi[1]-parPhi[4])+parPhi[2]*yIn+parPhi[3];
00176     deltaPhi = deltaY*CLHEP::pi/180.;
00177    
00178     HepPoint3D rotateVector(0,1,0);   //y axis
00179     rotateVector.rotateZ(center01.phi());
00180     posCorr.setZ(posCorr.z()-IRShift);
00181     posCorr.rotate(-deltaTheta,rotateVector);
00182     posCorr.setZ(posCorr.z()+IRShift);
00183     //phi parameter is gotten by the average of e+ e- peak
00184     if(Para.MethodMode()==0){
00185       posCorr.rotateZ(-0.002994);
00186     }else if(Para.MethodMode()==1){
00187       posCorr.rotateZ(-deltaPhi);
00188     }
00189 
00190     if(thetaModule>21) {
00191       posCorr.setZ(-posCorr.z());
00192     }
00193 
00194     double lengthemc;
00195     lengthemc = abs(posCorr.z()/cos(posCorr.theta()));
00196     // for Data
00197     double posDataCorPar;
00198     if(Para.DataMode()==1) {
00199       posDataCorPar=Para.BarrPosDataCor(thetaModule,phiModule); 
00200       posCorr.setTheta(posCorr.theta()-posDataCorPar/lengthemc);
00201     }
00202     
00203     //for MC
00204     double posMCCorPar;
00205     if(Para.DataMode()==0) {  
00206       posMCCorPar=Para.BarrPosMCCor(thetaModule,phiModule);     
00207       posCorr.setTheta(posCorr.theta()-posMCCorPar/lengthemc);
00208     }
00209     
00210   }   else {  //endcap position corretion
00211     if(Para.MethodMode()==1){
00212 
00213       double IRShift=0, parTheta[5],parPhi[5];  
00214       
00215       if(module==0){  //east endcap
00216         
00217         IRShift = 10;
00218         
00219         for(int i=0;i<5;i++) {
00220           
00221           
00222           if(e5x5>1.0)                   parTheta[i] = Para.EastLogThetaPara(thetaModule,i);
00223           else if(e5x5<=1.0 && e5x5>0.5) parTheta[i] = Para.EastLogThetaPara(thetaModule+6,i);
00224           else if (e5x5<=0.5)            parTheta[i] = Para.EastLogThetaPara(thetaModule+6,i);
00225           
00226           if(e5x5>1.0)                    parPhi[i] = Para.EastLogPhiPara(0,i);
00227           else if(e5x5<=1.0 && e5x5>0.5)  parPhi[i] = Para.EastLogPhiPara(1,i);
00228           else if (e5x5<=0.5)             parPhi[i] = Para.EastLogPhiPara(2,i);
00229         }
00230         
00231       }else if (module==2){  //west endcap          
00232         IRShift = -10;
00233         for(int i=0;i<5;i++){
00234           
00235           if(e5x5>1.0)                   parTheta[i] = Para.WestLogThetaPara(thetaModule,i);
00236           else if(e5x5<=1.0 && e5x5>0.5) parTheta[i] = Para.WestLogThetaPara(thetaModule+6,i);
00237           else if (e5x5<=0.5)            parTheta[i] = Para.WestLogThetaPara(thetaModule+6,i);
00238           
00239           if(e5x5>1.0)                    parPhi[i] = Para.WestLogPhiPara(0,i);
00240           else if(e5x5<=1.0 && e5x5>0.5)  parPhi[i] = Para.WestLogPhiPara(1,i);
00241           else if (e5x5<=0.5)             parPhi[i] = Para.WestLogPhiPara(2,i);
00242           
00243         }
00244        
00245       }
00246      
00247       HepPoint3D IR(0,0,IRShift);
00248      
00249       double theta12=-9999,theta03=-9999,theta23=-9999,theta14=-9999,delta=-9999;
00250       double phi01=-9999,phi23=-9999,phi12=-9999,phi34=-9999,deltaphi=-9999; 
00251       double xIn,deltaX,deltaTheta,yIn,deltaY,deltaPhi;  
00252       
00253       double posphi=posCorr.phi();
00254       if(phiModule==0){
00255         posphi = posphi;
00256       }else if(  ((thetaModule==0||thetaModule==1)&&phiModule==63)
00257                  ||((thetaModule==2||thetaModule==3)&&phiModule==79)
00258                  ||((thetaModule==4||thetaModule==5)&&phiModule==95)
00259                  ){     
00260         posphi = posphi+2*CLHEP::pi;    
00261       }else {
00262         posphi = posphi < 0 ? posphi+2*CLHEP::pi : posphi;      
00263       }  
00264       
00265       HepPoint3D center,center01, center23, center12,center03,center34,center14;
00266       
00267       
00268       if(  (thetaModule==1&&((phiModule+3)%4 == 0||(phiModule+2)%4 == 0))
00269            ||(thetaModule==3&&((phiModule+4)%5 == 0||(phiModule+3)%5 == 0||(phiModule+2)%5 == 0))) //pentagon
00270         {
00271           
00272           center23=(iGeoSvc->GetCrystalPoint(id,2)+ iGeoSvc->GetCrystalPoint(id,3))/2.0;
00273           center14=(iGeoSvc->GetCrystalPoint(id,1)+ iGeoSvc->GetCrystalPoint(id,4))/2.0;
00274           center=center23;
00275           theta23 = (center23 - IR).theta();
00276           theta14 = (center14 - IR).theta();
00277           delta = theta14 - theta23;
00278           xIn = (((posCorr-IR).theta() - theta23) - delta*0.5)/delta;
00279           
00280           
00281           center12=(iGeoSvc->GetCrystalPoint(id,1)+ iGeoSvc->GetCrystalPoint(id,2))/2.0;
00282           center34=(iGeoSvc->GetCrystalPoint(id,3)+ iGeoSvc->GetCrystalPoint(id,4))/2.0;
00283           phi12 = center12.phi();
00284           phi34 = center34.phi();
00285           
00286           phi12 = phi12 <0 ? phi12+2*CLHEP::pi : phi12;
00287           phi34 = phi34 <0 ? phi34+2*CLHEP::pi : phi34;
00288           deltaphi = phi34 - phi12;
00289           yIn = ((posphi - phi12) - deltaphi*0.5)/deltaphi;    
00290           
00291         } else {
00292         
00293         center12=(iGeoSvc->GetCrystalPoint(id,1)+ iGeoSvc->GetCrystalPoint(id,2))/2.0;
00294         center03=(iGeoSvc->GetCrystalPoint(id,0)+ iGeoSvc->GetCrystalPoint(id,3))/2.0;
00295         center=center12;
00296         
00297         theta12 = (center12 - IR).theta();
00298         theta03 = (center03 - IR).theta();
00299         delta = theta03 - theta12;
00300         xIn = (((posCorr-IR).theta() - theta12) - delta*0.5)/delta;
00301         
00302         
00303         center01=(iGeoSvc->GetCrystalPoint(id,0)+ iGeoSvc->GetCrystalPoint(id,1))/2.0;
00304         center23=(iGeoSvc->GetCrystalPoint(id,2)+ iGeoSvc->GetCrystalPoint(id,3))/2.0;
00305         phi01 = center01.phi();
00306         phi23 = center23.phi();
00307         phi01 = phi01 <0 ? phi01+2*CLHEP::pi : phi01;
00308         phi23 = phi23 <0 ? phi23+2*CLHEP::pi : phi23;
00309         deltaphi = phi23 - phi01;
00310         yIn = ((posphi - phi01) - deltaphi*0.5)/deltaphi;    
00311         
00312       }
00313       //deltaX=xIn-Xext => Xext=xIn-deltaX
00314       deltaX = parTheta[0]*atan(xIn*parTheta[1]-parTheta[4]) + parTheta[2]*xIn + parTheta[3] ;
00315       deltaTheta = deltaX*delta;   
00316       //deltaY=yIn-Ymc => Ymc=yIn-deltaY
00317       deltaY = parPhi[0]*atan(yIn*parPhi[1]-parPhi[4]) + parPhi[2]*yIn + parPhi[3] ;
00318       deltaPhi = deltaY*deltaphi;
00319       
00320       HepPoint3D rotateVector(0,1,0);   //y axis
00321       rotateVector.rotateZ(center.phi());
00322       posCorr.setZ(posCorr.z()-IRShift);
00323       posCorr.rotate(-deltaTheta,rotateVector);
00324       posCorr.setZ(posCorr.z()+IRShift);
00325       posCorr.rotateZ(-deltaPhi);
00326       
00327       double lengthemc;
00328       lengthemc = abs(posCorr.z()/cos(posCorr.theta()));
00329       double posDataCorPar;
00330       if(Para.DataMode()==1) {
00331         if(module==0) posDataCorPar=Para.EastPosDataCor(thetaModule,phiModule);
00332         if(module==2) posDataCorPar=Para.WestPosDataCor(thetaModule,phiModule); 
00333         posCorr.setTheta(posCorr.theta()-posDataCorPar/lengthemc);
00334       }         
00335     }
00336 
00337   } 
00338 
00339 
00340   //PosCorr=0 without position correction
00341   //PosCorr=1 with position correction
00342   if(Para.PosCorr()==0){ aShower.setPosition(pos);}
00343 
00344   if(Para.PosCorr()==1){
00345     aShower.setPosition(posCorr);
00346   }
00348   if(aShower.energy()<1e-5) return;
00349 
00350   double r,theta,phi;
00351   double dtheta,dphi,dx,dy,dz;
00352 
00353   r = posCorr.mag();
00354   theta = posCorr.theta();
00355   phi = posCorr.phi();
00356   //cout<<"energy: "<<aShower.energy()<<" theta: "<<theta<<" phi: "<<phi<<endl;
00357   
00358   if(aShower.energy()>0) {
00359     dtheta = Para.SigTheta(0)/sqrt(aShower.energy())+Para.SigTheta(1);
00360     dphi = Para.SigPhi(0)/sqrt(aShower.energy())+Para.SigPhi(1);
00361   }
00362   else {
00363     dtheta = 0;
00364     dphi = 0;
00365   }
00366   //cout<<"dtheta: "<<dtheta<<" dphi: "<<dphi<<endl;
00367   
00368   dx = fabs(iGeoSvc->GetBarrelR()*sin(phi)*dphi);
00369   dy = fabs(iGeoSvc->GetBarrelR()*cos(phi)*dphi);
00370   if(theta>0)
00371     dz = fabs(iGeoSvc->GetBarrelR()*dtheta/(sin(theta)*sin(theta)));
00372   else
00373     dz = 0;
00374   //cout<<" dx: "<<dx<<" dy: "<<dy<<" dz: "<<dz<<endl;
00375   
00376   aShower.setDtheta(dtheta);
00377   aShower.setDphi(dphi);
00378 
00379   //----------------------------------------------------------------------
00380   // Error matrix
00381   HepSymMatrix matrix(3); //error matrix of r, theta, phi
00382   matrix[0][0]=0;
00383   matrix[1][1]=dtheta*dtheta;
00384   matrix[2][2]=dphi*dphi;
00385   matrix[0][1]=0;
00386   matrix[0][2]=0;
00387   matrix[1][2]=0;
00388   //cout<<"matrix:\n"<<matrix<<endl;
00389 
00390   HepMatrix matrix1(3,3),matrix2(3,3);
00391   matrix1[0][0]=sin(theta)*cos(phi);
00392   matrix1[0][1]=r*cos(theta)*cos(phi);
00393   matrix1[0][2]=-r*sin(theta)*sin(phi);
00394   matrix1[1][0]=sin(theta)*sin(phi);
00395   matrix1[1][1]=r*cos(theta)*sin(phi);
00396   matrix1[1][2]=r*sin(theta)*cos(phi);
00397   matrix1[2][0]=cos(theta);
00398   matrix1[2][1]=-r*sin(theta);
00399   matrix1[2][2]=0;
00400   //cout<<"matrix1:\n"<<matrix1<<endl;
00401 
00402   for(int i=0;i<3;i++) {
00403     for(int j=0;j<3;j++) {
00404       matrix2[i][j]=matrix1[j][i];
00405     }
00406   }
00407   //cout<<"matrix2:\n"<<matrix2<<endl;
00408   
00409   HepMatrix matrix4(3,3);
00410   matrix4=matrix1*matrix*matrix2;
00411   //cout<<"matrix4:\n"<<matrix4<<endl;
00412   
00413   HepSymMatrix matrix5(3); //error matrix of x, y, z
00414   for(int i=0;i<3;i++) {
00415     for(int j=0;j<=i;j++) {
00416       matrix5[i][j]=matrix4[i][j];
00417     }
00418   }
00419   //cout<<"matrix5:\n"<<matrix5<<endl;
00420 
00421   aShower.setErrorMatrix(matrix5);
00422 
00423   if(matrix5[0][0]>0) dx=sqrt(matrix5[0][0]);
00424   if(matrix5[1][1]>0) dy=sqrt(matrix5[1][1]);
00425   if(matrix5[2][2]>0) dz=sqrt(matrix5[2][2]);
00426 }

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

Public Member Functions

Detailed Description

Constructor & Destructor Documentation

Member Function Documentation
