00001
00002
00003
00004
00005
00006 #include "EmcRec/EmcRecShowerPosLogShMax.h"
00007
00008 #include <iostream>
00009 #include <fstream>
00010 #include "EmcRecGeoSvc/EmcRecGeoSvc.h"
00011 #include "GaudiKernel/Bootstrap.h"
00012 #include "GaudiKernel/ISvcLocator.h"
00013
00014 #include "EmcRec/EmcRecParameter.h"
00015 using namespace std;
00016
00017 void EmcRecShowerPosLogShMax::Position(RecEmcShower& aShower)
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
00026
00027 IEmcRecGeoSvc* iGeoSvc;
00028 ISvcLocator* svcLocator = Gaudi::svcLocator();
00029 StatusCode sc = svcLocator->service("EmcRecGeoSvc",iGeoSvc);
00030 if(sc!=StatusCode::SUCCESS) {
00031
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
00045
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
00052
00053 double theDistance;
00054 if(cit->second.getCellId()==seedID) {
00055
00056 theDistance=(seedPoint.mag()+ltot);
00057
00058 } else {
00059
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
00075
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
00084
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
00104
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
00114
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
00128 }
00129 }
00130 }
00131
00132 if(wsum<=0) {
00133 cout<<"[[[[[[[[[[[[[[["<<wsum;
00134 }
00135 pos=possum/wsum;
00136
00137
00138
00139
00140 if(Para.PosCorr()==0){ aShower.setPosition(pos);}
00141
00142
00143
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) {
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;
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);
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
00199
00200
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
00221 deltaPhi = yOut*CLHEP::pi/180.;
00222
00223 HepPoint3D rotateVector(0,1,0);
00224 rotateVector.rotateZ(center01.phi());
00225 posCorr.setZ(posCorr.z()-IRShift);
00226 posCorr.rotate(-deltaTheta,rotateVector);
00227 posCorr.setZ(posCorr.z()+IRShift);
00228
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
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
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
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
00267
00268 aShower.setDtheta(dtheta);
00269 aShower.setDphi(dphi);
00270
00271
00272
00273 HepSymMatrix matrix(3);
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
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
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
00300
00301 HepMatrix matrix4(3,3);
00302 matrix4=matrix1*matrix*matrix2;
00303
00304
00305 HepSymMatrix matrix5(3);
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
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 }
00319