00001
00002
00003
00004
00005
00006 #include "EmcRec/EmcRecShowerPosLinShMax.h"
00007
00008 #include <iostream>
00009 #include <fstream>
00010 #include "EmcRecGeoSvc/EmcRecGeoSvc.h"
00011 #include "GaudiKernel/Bootstrap.h"
00012 #include "GaudiKernel/ISvcLocator.h"
00013 #include "EmcRec/EmcRecParameter.h"
00014 #include "CLHEP/Vector/LorentzVector.h"
00015
00016 using namespace std;
00017
00018 void EmcRecShowerPosLinShMax::Position(RecEmcShower& aShower)
00019 {
00020 RecEmcFractionMap::const_iterator cit;
00021 HepPoint3D pos(0,0,0);
00022 HepPoint3D posMax(0,0,0);
00023 HepPoint3D possum(0,0,0);
00024 double w,wsum=0;
00025 double etot;
00026
00027
00028
00029 IEmcRecGeoSvc* iGeoSvc;
00030 ISvcLocator* svcLocator = Gaudi::svcLocator();
00031 StatusCode sc = svcLocator->service("EmcRecGeoSvc",iGeoSvc);
00032 if(sc!=StatusCode::SUCCESS) {
00033
00034 }
00035
00036 EmcRecParameter& Para=EmcRecParameter::GetInstance();
00037
00038 for(cit=aShower.Begin();cit!=aShower.End();cit++){
00039 etot+=(cit->second.getEnergy()*cit->second.getFraction());
00040
00041 }
00042 wsum=0;
00043
00044 double etota=aShower.getEAll();
00045 double ltot=(log(etota/0.0145)+0.5)*1.86;
00046 double e55=aShower.e5x5();
00047 double len55=(log(etota/0.0145)+0.5)*1.86;
00048 Identifier seedID = aShower.getShowerId();
00049 HepPoint3D seedPoint(iGeoSvc->GetCFrontCenter(seedID));
00050
00051
00052 for(cit=aShower.Begin();cit!=aShower.End();cit++){
00053 w=(cit->second.getEnergy()*cit->second.getFraction());
00054
00055 HepPoint3D hitPoint(iGeoSvc->GetCFrontCenter(cit->second.getCellId()));
00056
00057
00058 double theDistance=0;
00059 if(cit->second.getCellId()==seedID) {
00060
00061 theDistance=(seedPoint.mag()+ltot);
00062 } else {
00063 theDistance=(seedPoint.mag()+ltot)/cos(seedPoint.angle(hitPoint));
00064
00065 }
00066
00067
00068 wsum+=w;
00069 pos=iGeoSvc->GetCFrontCenter(cit->second.getCellId());
00070
00071 posMax = (theDistance/pos.mag())*pos;
00072
00073 possum+=posMax*w;
00074 }
00075
00076 if(wsum<=0) {
00077 cout<<"[[[[[[[[[[[[[[["<<wsum;
00078 }
00079 pos=possum/wsum;
00080
00081
00082 if(Para.PosCorr()==0){ aShower.setPosition(pos);}
00083
00084
00085
00086 RecEmcID id = aShower.getShowerId();
00087 const unsigned int module = EmcID::barrel_ec(id);;
00088 const unsigned int thetaModule = EmcID::theta_module(id);
00089 const unsigned int phiModule = EmcID::phi_module(id);
00090 HepPoint3D posCorr(pos);
00091
00092 if(module==1) {
00093 unsigned int theModule;
00094 if(thetaModule>21) {
00095 theModule = 43 - thetaModule;
00096 id = EmcID::crystal_id(module,theModule,phiModule);
00097 posCorr.setZ(-posCorr.z());
00098 } else {
00099 theModule = thetaModule;
00100 }
00101
00102 double IRShift, parTheta[5],parPhi[5];
00103 if(theModule==21) {
00104 IRShift = 0;
00105 } else if(theModule==20) {
00106 IRShift = 2.5;
00107 } else {
00108 IRShift = 5.;
00109 }
00110
00111 EmcRecParameter& Para=EmcRecParameter::GetInstance();
00112
00113 for(int i=0;i<5;i++){
00114
00115 if(e55>1.0) parTheta[i] = Para.BarrShLinThetaPara(theModule,i);
00116 else if(e55<=1.0 && e55>0.5) parTheta[i] = Para.BarrShLinThetaPara(theModule+22,i);
00117 else if(e55<=0.5) parTheta[i] = Para.BarrShLinThetaPara(theModule+44,i);
00118
00119
00120 if(e55>1.0) parPhi[i] = Para.BarrShLinPhiPara(0,i);
00121 else if(e55<=1.0 && e55>0.5) parPhi[i] = Para.BarrShLinPhiPara(1,i);
00122 else if(e55<=0.5) parPhi[i] = Para.BarrShLinPhiPara(2,i);
00123
00124 }
00125
00126 HepPoint3D IR(0,0,IRShift);
00127 HepPoint3D center01, center23;
00128 center01 = (iGeoSvc->GetCrystalPoint(id,0)+iGeoSvc->GetCrystalPoint(id,1))/2;
00129 center23 = (iGeoSvc->GetCrystalPoint(id,2)+iGeoSvc->GetCrystalPoint(id,3))/2;
00130
00131 double theta01,theta23,length,d;
00132 theta01 = (center01-IR).theta();
00133 theta23 = (center23-IR).theta();
00134 length = (center01-IR).mag();
00135 d = length*tan(theta01-theta23);
00136
00137 double xIn,xOut,deltaTheta;
00138 xIn = length*tan(theta01-(posCorr-IR).theta())-d/2;
00139 xOut = xIn-( parTheta[0]*atan(xIn*parTheta[1]-parTheta[4])
00140 + parTheta[2]*xIn + parTheta[3] );
00141 deltaTheta = atan((xOut+d/2)/length)-atan((xIn+d/2)/length);
00142
00143
00144
00146 double yIn,yOut,deltaPhi;
00147
00148 if(phiModule==0&&posCorr.phi()<0){
00149
00150 yIn = posCorr.phi()*180./CLHEP::pi+360.-119*3-1.843;
00151
00152 }else if (phiModule==119&&posCorr.phi()>0){
00153
00154 yIn = posCorr.phi()*180./CLHEP::pi-1.843;
00155
00156 }else {
00157
00158 yIn = posCorr.phi() < 0 ? posCorr.phi()*180./CLHEP::pi+360.-phiModule*3.-1.843
00159 : posCorr.phi()*180./CLHEP::pi-phiModule*3.-1.843;
00160
00161 }
00162 yOut = parPhi[0]*atan(yIn*parPhi[1]-parPhi[4])+parPhi[2]*yIn+parPhi[3];
00163
00164
00165
00166 deltaPhi = yOut*CLHEP::pi/180.;
00167
00168 HepPoint3D rotateVector(0,1,0);
00169 rotateVector.rotateZ(center01.phi());
00170 posCorr.setZ(posCorr.z()-IRShift);
00171 posCorr.rotate(-deltaTheta,rotateVector);
00172 posCorr.setZ(posCorr.z()+IRShift);
00173
00174
00175 posCorr.rotateZ(-deltaPhi);
00176
00177 if(thetaModule>21) {
00178 posCorr.setZ(-posCorr.z());
00179 }
00180 }
00181 posCorr=((posCorr.mag()-len55)/posCorr.mag())*posCorr;
00182
00183 if(Para.PosCorr()==1){ aShower.setPosition(posCorr);}
00184
00186 if(aShower.energy()<1e-5) return;
00187
00188 double r,theta,phi;
00189 double dtheta,dphi,dx,dy,dz;
00190
00191 r = posCorr.mag();
00192 theta = posCorr.theta();
00193 phi = posCorr.phi();
00194
00195
00196 if(aShower.energy()>0) {
00197 dtheta = Para.SigTheta(0)/sqrt(aShower.energy())+Para.SigTheta(1);
00198 dphi = Para.SigPhi(0)/sqrt(aShower.energy())+Para.SigPhi(1);
00199 }
00200 else {
00201 dtheta = 0;
00202 dphi = 0;
00203 }
00204
00205
00206 dx = fabs(iGeoSvc->GetBarrelR()*sin(phi)*dphi);
00207 dy = fabs(iGeoSvc->GetBarrelR()*cos(phi)*dphi);
00208 if(theta>0)
00209 dz = fabs(iGeoSvc->GetBarrelR()*dtheta/(sin(theta)*sin(theta)));
00210 else
00211 dz = 0;
00212
00213
00214 aShower.setDtheta(dtheta);
00215 aShower.setDphi(dphi);
00216
00217
00218
00219 HepSymMatrix matrix(3);
00220 matrix[0][0]=0;
00221 matrix[1][1]=dtheta*dtheta;
00222 matrix[2][2]=dphi*dphi;
00223 matrix[0][1]=0;
00224 matrix[0][2]=0;
00225 matrix[1][2]=0;
00226
00227
00228 HepMatrix matrix1(3,3),matrix2(3,3);
00229 matrix1[0][0]=sin(theta)*cos(phi);
00230 matrix1[0][1]=r*cos(theta)*cos(phi);
00231 matrix1[0][2]=-r*sin(theta)*sin(phi);
00232 matrix1[1][0]=sin(theta)*sin(phi);
00233 matrix1[1][1]=r*cos(theta)*sin(phi);
00234 matrix1[1][2]=r*sin(theta)*cos(phi);
00235 matrix1[2][0]=cos(theta);
00236 matrix1[2][1]=-r*sin(theta);
00237 matrix1[2][2]=0;
00238
00239
00240 for(int i=0;i<3;i++) {
00241 for(int j=0;j<3;j++) {
00242 matrix2[i][j]=matrix1[j][i];
00243 }
00244 }
00245
00246
00247 HepMatrix matrix4(3,3);
00248 matrix4=matrix1*matrix*matrix2;
00249
00250
00251 HepSymMatrix matrix5(3);
00252 for(int i=0;i<3;i++) {
00253 for(int j=0;j<=i;j++) {
00254 matrix5[i][j]=matrix4[i][j];
00255 }
00256 }
00257
00258
00259 aShower.setErrorMatrix(matrix5);
00260
00261 if(matrix5[0][0]>0) dx=sqrt(matrix5[0][0]);
00262 if(matrix5[1][1]>0) dy=sqrt(matrix5[1][1]);
00263 if(matrix5[2][2]>0) dz=sqrt(matrix5[2][2]);
00264 }
00265