BOSS_7.0.2: EmcRecShowerPosLogShMax Class Reference

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

Definition at line 17 of file EmcRecShowerPosLogShMax.cxx.

References EmcID::barrel_ec(), EmcRecParameter::BarrShLogPhiPara(), EmcRecParameter::BarrShLogThetaPara(), RecEmcShower::Begin(), cos(), EmcID::crystal_id(), DstEmcShower::e3x3(), DstEmcShower::e5x5(), 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(), 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.

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


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

EmcRecShowerPosLogShMax Class Reference

Public Member Functions

Detailed Description

Constructor & Destructor Documentation

Member Function Documentation