EmcRecShowerPosLoglin Class Reference

#include <EmcRecShowerPosLoglin.h>

Inheritance diagram for EmcRecShowerPosLoglin:

EmcRecShowerPosAbs List of all members.

Public Member Functions

 EmcRecShowerPosLoglin ()
virtual ~EmcRecShowerPosLoglin ()
virtual void Position (RecEmcShower &aShower)

Detailed Description

Definition at line 17 of file EmcRecShowerPosLoglin.h.


Constructor & Destructor Documentation

EmcRecShowerPosLoglin::EmcRecShowerPosLoglin (  )  [inline]

Definition at line 21 of file EmcRecShowerPosLoglin.h.

00021 {}

virtual EmcRecShowerPosLoglin::~EmcRecShowerPosLoglin (  )  [inline, virtual]

Definition at line 22 of file EmcRecShowerPosLoglin.h.

00022 {}


Member Function Documentation

void EmcRecShowerPosLoglin::Position ( RecEmcShower aShower  )  [virtual]

Implements EmcRecShowerPosAbs.

Definition at line 19 of file EmcRecShowerPosLoglin.cxx.

References a0, EmcID::barrel_ec(), EmcRecParameter::BarrLoglinPhiPara(), EmcRecParameter::BarrLoglinThetaPara(), RecEmcShower::Begin(), cos(), EmcID::crystal_id(), DstEmcShower::e3x3(), DstEmcShower::e5x5(), RecEmcShower::End(), DstEmcShower::energy(), exp(), IEmcRecGeoSvc::GetBarrelR(), IEmcRecGeoSvc::GetCFrontCenter(), IEmcRecGeoSvc::GetCrystalPoint(), RecEmcShower::getEAll(), RecEmcShower::getFractionMap3x3(), RecEmcShower::getFractionMap5x5(), EmcRecParameter::GetInstance(), RecEmcShower::getShowerId(), genRecEmupikp::i, ganga-rec::j, num, 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(), and w.

00020 {
00021   RecEmcFractionMap::const_iterator cit;
00022   HepPoint3D pos(0,0,0);
00023   HepPoint3D possum(0,0,0);
00024   double w,w1,w2,wsum=0;
00025   double num,num0;
00026 
00027   // cout<<"EmcRecShowerPosLoglin::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 
00039   //  double offset=Para.LogPosOffset();
00040   double a0=4;
00041   double offset;
00042   //cout<<offset<<endl;
00043   double e5x5 = -99999;
00044   if(Para.PositionMode2()=="all") {
00045     double etot=aShower.getEAll();
00046     for(cit=aShower.Begin();cit!=aShower.End();cit++){
00047       w=offset+log((cit->second.getEnergy()/etot)*cit->second.getFraction());
00048       if(w>0){
00049         wsum+=w;
00050         pos=iGeoSvc->GetCFrontCenter(cit->second.getCellId());
00051         possum+=pos*w;
00052       }
00053     }
00054   } else if(Para.PositionMode2()=="3x3") {
00055     double e3x3=aShower.e3x3();
00056     RecEmcFractionMap fracMap3x3=aShower.getFractionMap3x3();
00057     for(cit=fracMap3x3.begin();cit!=fracMap3x3.end();cit++) {
00058       w=offset+log((cit->second.getEnergy()/e3x3)*cit->second.getFraction());
00059       if(w>0){
00060         wsum+=w;
00061         pos=iGeoSvc->GetCFrontCenter(cit->second.getCellId());
00062         possum+=pos*w;
00063       }
00064     }
00065   } else {
00066     e5x5=aShower.e5x5();
00067     offset=a0-1.594*exp(-2.543*e5x5);
00068     RecEmcFractionMap fracMap5x5=aShower.getFractionMap5x5();
00069 
00070     num0=0;
00071     num=0;
00072     for(cit=fracMap5x5.begin();cit!=fracMap5x5.end();cit++) {
00073       w1=offset+log((cit->second.getEnergy()/e5x5)*cit->second.getFraction());      
00074       w2=(cit->second.getEnergy()/e5x5)*cit->second.getFraction();
00075       
00076       num0++;
00077       
00078       if(w1>0){
00079         num++;
00080       }
00081       
00082     }
00083 
00084     // cout<<"num0=="<<num0<<endl;
00085     //    cout<<"num==="<<num<<endl;
00086     for(cit=fracMap5x5.begin();cit!=fracMap5x5.end();cit++) {
00087       w1=offset+log((cit->second.getEnergy()/e5x5)*cit->second.getFraction());      
00088       w2=(cit->second.getEnergy()/e5x5)*cit->second.getFraction();
00089       
00090       num0++;
00091      
00092       if(w1>0){
00093         if(num>1){
00094           //      cout<<"w1 be used"<<endl;
00095           wsum+=w1;
00096           pos=iGeoSvc->GetCFrontCenter(cit->second.getCellId());
00097           possum+=pos*w1;
00098         } else{
00099           //      cout<<"w2 be used"<<endl;
00100           wsum+=w2;
00101           pos=iGeoSvc->GetCFrontCenter(cit->second.getCellId());
00102           possum+=pos*w2;
00103         }
00104       }           
00105     }
00106 
00107     
00108   }
00109   if(wsum<=0) {
00110     //    cout<<"[[[[[[[[[[[[[[["<<wsum<<endl;
00111   }
00112   pos=possum/wsum;
00113 
00114   //PosCorr=0 without position correction
00115   //PosCorr=1 with position correction
00116   if(Para.PosCorr()==0){        
00117 
00118     aShower.setPosition(pos);
00119   }
00120   
00121   //----------------------------------------------------------------------
00122   //position correction
00123   RecEmcID id = aShower.getShowerId();
00124   const unsigned int module = EmcID::barrel_ec(id);;
00125   const unsigned int thetaModule = EmcID::theta_module(id);
00126   const unsigned int phiModule = EmcID::phi_module(id);
00127   HepPoint3D posCorr(pos);
00128   
00129   if(module==1) {   //barrel
00130     unsigned int theModule;
00131     if(thetaModule>21) {
00132       theModule = 43 - thetaModule;
00133       id = EmcID::crystal_id(module,theModule,phiModule);
00134       posCorr.setZ(-posCorr.z());
00135     } else {
00136       theModule = thetaModule;
00137     }
00138     
00139     double IRShift=0, parTheta[5],parPhi[5];
00140     if(theModule==21) {
00141       IRShift = 0;
00142     } else if(theModule==20) {
00143       IRShift = 2.5;
00144     } else {
00145       IRShift = 5.;
00146     }  
00147 
00148  
00149       EmcRecParameter& Para=EmcRecParameter::GetInstance();
00150       for(int i=0;i<5;i++){
00151         
00152         parTheta[i] = Para.BarrLoglinThetaPara(theModule,i);
00153          
00154         parPhi[i] = Para.BarrLoglinPhiPara(0,i);
00155 
00156 //      cout<<"Para.BarrLoglinThetaPara"<<i<<"===="<<Para.BarrLoglinThetaPara(theModule,i)<<endl;
00157 //      cout<<"Para.BarrLoglinPhiPara"<<i<<"===="<<Para.BarrLoglinPhiPara(0,i)<<endl;
00158 
00159 //      cout<<"theModule===="<<theModule<<endl;
00160 //      cout<<"partheta"<<i<<"===="<<parTheta[i]<<endl;
00161 //      cout<<"parphi"<<i<<"===="<<parPhi[i]<<endl;
00162 
00163         
00164       }
00165     
00166     HepPoint3D IR(0,0,IRShift);
00167     HepPoint3D center01, center23;  //center of point0,1 and point2,3 in crystal
00168     center01 = (iGeoSvc->GetCrystalPoint(id,0)+iGeoSvc->GetCrystalPoint(id,1))/2;
00169     center23 = (iGeoSvc->GetCrystalPoint(id,2)+iGeoSvc->GetCrystalPoint(id,3))/2;
00170 
00171     double theta01,theta23,length,d;
00172     theta01 = (center01-IR).theta();
00173     theta23 = (center23-IR).theta();
00174     length = (center01-IR).mag();
00175     d = length*tan(theta01-theta23);  //length in x direction
00176     
00177     double xIn,xOut,deltaTheta;
00178     xIn = length*tan(theta01-(posCorr-IR).theta())-d/2;
00179     xOut = xIn-( parTheta[0]*atan(xIn*parTheta[1]-parTheta[4])+ parTheta[2]*xIn + parTheta[3] );
00180     deltaTheta = atan((xOut+d/2)/length)-atan((xIn+d/2)/length);
00181     //cout<<"xIn: "<<xIn<<" xOut: "<<xOut<<" deltaTheta: "<<deltaTheta<<endl;
00182     //cout<<"dx: "<<parTheta[0]*atan(xIn*parTheta[1]-parTheta[4])
00183     //      + parTheta[2]*xIn + parTheta[3] <<endl;
00184     
00185     //obtained by photon, not used, just for comparison
00186     //    double parPhi[5] = { 51.1, -0.1499, 7.62, 0.1301, 0.005581 };
00187     double yIn,yOut,deltaPhi;
00188     //    yIn = posCorr.phi() < 0 ? posCorr.phi()*180./CLHEP::pi+360.-phiModule*3.-1.1537
00189     //      : posCorr.phi()*180./CLHEP::pi-phiModule*3.-1.1537;
00190     //    yOut = parPhi[0]*atan(yIn*parPhi[1]-parPhi[4])+parPhi[2]*yIn+parPhi[3];
00191     
00192     yIn = posCorr.phi() < 0 ? posCorr.phi()*180./CLHEP::pi+360.-phiModule*3.-1.843
00193       : posCorr.phi()*180./CLHEP::pi-phiModule*3.-1.843;
00194 
00195     if(phiModule==0&&posCorr.phi()<0){
00196       yIn = posCorr.phi()*180./CLHEP::pi+360.-119*3-1.843;
00197     }else if (phiModule==119&&posCorr.phi()>0){
00198       yIn = posCorr.phi()*180./CLHEP::pi-1.843;
00199     }else {
00200       yIn = posCorr.phi() < 0 ? posCorr.phi()*180./CLHEP::pi+360.-phiModule*3.-1.843 : posCorr.phi()*180./CLHEP::pi-phiModule*3.-1.843;
00201     }
00202     yOut = parPhi[0]*atan(yIn*parPhi[1]-parPhi[4])+parPhi[2]*yIn+parPhi[3];
00203     deltaPhi = yOut*CLHEP::pi/180.;
00204    
00205     HepPoint3D rotateVector(0,1,0);   //y axis
00206     rotateVector.rotateZ(center01.phi());
00207     posCorr.setZ(posCorr.z()-IRShift);
00208     posCorr.rotate(-deltaTheta,rotateVector);
00209     posCorr.setZ(posCorr.z()+IRShift);
00210     //phi parameter is gotten by the average of e+ e- peak
00211    
00212       posCorr.rotateZ(-deltaPhi);
00213     
00214     if(thetaModule>21) {
00215       posCorr.setZ(-posCorr.z());
00216     }
00217   }  
00218   if(Para.PosCorr()==1){  
00219         
00220     aShower.setPosition(posCorr);
00221   }
00223   if(aShower.energy()<1e-5) return;
00224 
00225   double r,theta,phi;
00226   double dtheta,dphi,dx,dy,dz;
00227 
00228   r = posCorr.mag();
00229   theta = posCorr.theta();
00230   phi = posCorr.phi();
00231   //cout<<"energy: "<<aShower.energy()<<" theta: "<<theta<<" phi: "<<phi<<endl;
00232   
00233   if(aShower.energy()>0) {
00234     dtheta = Para.SigTheta(0)/sqrt(aShower.energy())+Para.SigTheta(1);
00235     dphi = Para.SigPhi(0)/sqrt(aShower.energy())+Para.SigPhi(1);
00236   }
00237   else {
00238     dtheta = 0;
00239     dphi = 0;
00240   }
00241   //cout<<"dtheta: "<<dtheta<<" dphi: "<<dphi<<endl;
00242   
00243   dx = fabs(iGeoSvc->GetBarrelR()*sin(phi)*dphi);
00244   dy = fabs(iGeoSvc->GetBarrelR()*cos(phi)*dphi);
00245   if(theta>0)
00246     dz = fabs(iGeoSvc->GetBarrelR()*dtheta/(sin(theta)*sin(theta)));
00247   else
00248     dz = 0;
00249   //cout<<" dx: "<<dx<<" dy: "<<dy<<" dz: "<<dz<<endl;
00250   
00251   aShower.setDtheta(dtheta);
00252   aShower.setDphi(dphi);
00253 
00254   //----------------------------------------------------------------------
00255   // Error matrix
00256   HepSymMatrix matrix(3); //error matrix of r, theta, phi
00257   matrix[0][0]=0;
00258   matrix[1][1]=dtheta*dtheta;
00259   matrix[2][2]=dphi*dphi;
00260   matrix[0][1]=0;
00261   matrix[0][2]=0;
00262   matrix[1][2]=0;
00263   //cout<<"matrix:\n"<<matrix<<endl;
00264 
00265   HepMatrix matrix1(3,3),matrix2(3,3);
00266   matrix1[0][0]=sin(theta)*cos(phi);
00267   matrix1[0][1]=r*cos(theta)*cos(phi);
00268   matrix1[0][2]=-r*sin(theta)*sin(phi);
00269   matrix1[1][0]=sin(theta)*sin(phi);
00270   matrix1[1][1]=r*cos(theta)*sin(phi);
00271   matrix1[1][2]=r*sin(theta)*cos(phi);
00272   matrix1[2][0]=cos(theta);
00273   matrix1[2][1]=-r*sin(theta);
00274   matrix1[2][2]=0;
00275   //cout<<"matrix1:\n"<<matrix1<<endl;
00276 
00277   for(int i=0;i<3;i++) {
00278     for(int j=0;j<3;j++) {
00279       matrix2[i][j]=matrix1[j][i];
00280     }
00281   }
00282   //cout<<"matrix2:\n"<<matrix2<<endl;
00283   
00284   HepMatrix matrix4(3,3);
00285   matrix4=matrix1*matrix*matrix2;
00286   //cout<<"matrix4:\n"<<matrix4<<endl;
00287   
00288   HepSymMatrix matrix5(3); //error matrix of x, y, z
00289   for(int i=0;i<3;i++) {
00290     for(int j=0;j<=i;j++) {
00291       matrix5[i][j]=matrix4[i][j];
00292     }
00293   }
00294   //cout<<"matrix5:\n"<<matrix5<<endl;
00295 
00296   aShower.setErrorMatrix(matrix5);
00297 
00298   if(matrix5[0][0]>0) dx=sqrt(matrix5[0][0]);
00299   if(matrix5[1][1]>0) dy=sqrt(matrix5[1][1]);
00300   if(matrix5[2][2]>0) dz=sqrt(matrix5[2][2]);
00301 }


Generated on Tue Nov 29 23:18:45 2016 for BOSS_7.0.2 by  doxygen 1.4.7