#include <EmcRecShowerPosLoglin.h>
Inheritance diagram for EmcRecShowerPosLoglin:
Public Member Functions | |
EmcRecShowerPosLoglin () | |
virtual | ~EmcRecShowerPosLoglin () |
virtual void | Position (RecEmcShower &aShower) |
Definition at line 17 of file EmcRecShowerPosLoglin.h.
EmcRecShowerPosLoglin::EmcRecShowerPosLoglin | ( | ) | [inline] |
virtual EmcRecShowerPosLoglin::~EmcRecShowerPosLoglin | ( | ) | [inline, virtual] |
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 }