00001
00002
00003
00004
00005
00006 #include "EmcRec/EmcRecShowerPosLin.h"
00007
00008 #include <iostream>
00009 #include <fstream>
00010
00011 #include "EmcRecGeoSvc/EmcRecGeoSvc.h"
00012 #include "GaudiKernel/Bootstrap.h"
00013 #include "GaudiKernel/ISvcLocator.h"
00014
00015 #include "EmcRec/EmcRecParameter.h"
00016
00017 void EmcRecShowerPosLin::Position(RecEmcShower& aShower)
00018 {
00019 RecEmcFractionMap::const_iterator cit;
00020 HepPoint3D pos(0,0,0);
00021 HepPoint3D possum(0,0,0);
00022 double w,wsum;
00023 double etot;
00024
00025
00026 IEmcRecGeoSvc* iGeoSvc;
00027 ISvcLocator* svcLocator = Gaudi::svcLocator();
00028 StatusCode sc = svcLocator->service("EmcRecGeoSvc",iGeoSvc);
00029 if(sc!=StatusCode::SUCCESS) {
00030
00031 }
00032
00033 EmcRecParameter& Para=EmcRecParameter::GetInstance();
00034
00035 double e5x5 = aShower.e5x5();
00036
00037 etot=0;
00038
00039 for(cit=aShower.Begin();cit!=aShower.End();cit++){
00040 etot+=(cit->second.getEnergy()*cit->second.getFraction());
00041 }
00042 wsum=0;
00043 for(cit=aShower.Begin();cit!=aShower.End();cit++){
00044 w=(cit->second.getEnergy()*cit->second.getFraction());
00045 wsum+=w;
00046 pos=iGeoSvc->GetCFrontCenter(cit->second.getCellId());
00047 possum+=pos*w;
00048 }
00049
00050 if(wsum<=0) {
00051 cout<<"[[[[[[[[[[[[[[["<<wsum;
00052 }
00053
00054 pos=possum/wsum;
00055
00056
00057
00058
00059 if(Para.PosCorr()==0){ aShower.setPosition(pos);}
00060
00061
00062
00063
00064 RecEmcID id = aShower.getShowerId();
00065 const unsigned int module = EmcID::barrel_ec(id);;
00066 const unsigned int thetaModule = EmcID::theta_module(id);
00067 const unsigned int phiModule = EmcID::phi_module(id);
00068 HepPoint3D posCorr(pos);
00069
00070 if(module==1) {
00071 unsigned int theModule;
00072 if(thetaModule>21) {
00073 theModule = 43 - thetaModule;
00074 id = EmcID::crystal_id(module,theModule,phiModule);
00075 posCorr.setZ(-posCorr.z());
00076 } else {
00077 theModule = thetaModule;
00078 }
00079
00080
00081 EmcRecParameter& Para=EmcRecParameter::GetInstance();
00082 double IRShift, parTheta[5],parPhi[5];
00083
00084
00085 if(theModule==21) {
00086 IRShift = 0;
00087 } else if(theModule==20) {
00088 IRShift = 2.5;
00089 } else {
00090 IRShift = 5.;
00091 }
00092
00093
00094 if(Para.MethodMode()==0){
00095
00096
00097 for(int i=0;i<5;i++){
00098 if(theModule==21) {
00099 double par[5] = { 1.698, -1.553, 0.9384, 0.1193, -0.01916 };
00100 parTheta[i] = par[i];
00101 }else if(theModule==20){
00102 double par[5] = { 1.593, -1.582, 0.8881, 0.3298, -0.08856 };
00103 parTheta[i] = par[i];
00104 }else{
00105 double par[5] = { 1.605, -1.702, 0.8433, 0.3072, -0.1809 };
00106 parTheta[i] = par[i];
00107 }
00108 }
00109
00110 } else if(Para.MethodMode()==1){
00111
00112 for(int i=0;i<5;i++){
00113
00114 if(e5x5>1.0) parTheta[i] = Para.BarrLinThetaPara(theModule,i);
00115 else if(e5x5<=1.0 && e5x5>0.5) parTheta[i] = Para.BarrLinThetaPara(theModule+22,i);
00116 else if(e5x5<=0.5) parTheta[i] = Para.BarrLinThetaPara(theModule+44,i);
00117
00118 if(e5x5>1.0) parPhi[i] = Para.BarrLinPhiPara(0,i);
00119 else if(e5x5<=1.0 && e5x5>0.5) parPhi[i] = Para.BarrLinPhiPara(1,i);
00120 else if(e5x5<=0.5) parPhi[i] = Para.BarrLinPhiPara(2,i);
00121
00122 }
00123 }
00124
00125
00126
00127 HepPoint3D IR(0,0,IRShift);
00128 HepPoint3D center01, center23;
00129 center01 = (iGeoSvc->GetCrystalPoint(id,0)+iGeoSvc->GetCrystalPoint(id,1))/2;
00130 center23 = (iGeoSvc->GetCrystalPoint(id,2)+iGeoSvc->GetCrystalPoint(id,3))/2;
00131
00132
00133 double theta01,theta23,length,d;
00134 theta01 = (center01-IR).theta();
00135 theta23 = (center23-IR).theta();
00136 length = (center01-IR).mag();
00137 d = length*tan(theta01-theta23);
00138
00139
00140
00141 double xIn,xOut,deltaTheta;
00142 xIn = length*tan(theta01-(posCorr-IR).theta())-d/2;
00143 xOut = xIn-(parTheta[0]*atan(xIn*parTheta[1]-parTheta[4])+parTheta[2]*xIn+parTheta[3]);
00144 deltaTheta = atan((xOut+d/2)/length)-atan((xIn+d/2)/length);
00145
00146
00147 double parPhi1[5] = { 0.7185, -2.539, 0.6759, 0.3472, -0.3108 };
00148 double parPhi2[5] = { 0.5781, -2.917, 0.5516, -0.5745, 0.3777 };
00149
00150
00151
00152
00153
00154
00155
00156 double yIn,yOut,deltaPhi;
00157
00158
00159 if(phiModule==0&&posCorr.phi()<0){
00160
00161 yIn = posCorr.phi()*180./CLHEP::pi+360.-119*3-1.843;
00162
00163 }else if (phiModule==119&&posCorr.phi()>0){
00164
00165 yIn = posCorr.phi()*180./CLHEP::pi-1.843;
00166
00167 }else {
00168
00169 yIn = posCorr.phi() < 0 ? posCorr.phi()*180./CLHEP::pi+360.-phiModule*3.-1.843
00170 : posCorr.phi()*180./CLHEP::pi-phiModule*3.-1.843;
00171
00172 }
00173
00174 yOut = parPhi[0]*atan(yIn*parPhi[1]-parPhi[4])+parPhi[2]*yIn+parPhi[3];
00175 deltaPhi = yOut*CLHEP::pi/180.;
00176
00177
00178 HepPoint3D rotateVector(0,1,0);
00179 rotateVector.rotateZ(posCorr.phi());
00180 posCorr.setZ(posCorr.z()-IRShift);
00181 posCorr.rotate(-deltaTheta,rotateVector);
00182 posCorr.setZ(posCorr.z()+IRShift);
00183 if(Para.MethodMode()==0){
00184 posCorr.rotateZ(-0.002994);
00185 }else if(Para.MethodMode()==1){
00186 posCorr.rotateZ(-deltaPhi);
00187 }
00188
00189 if(thetaModule>21) {
00190 posCorr.setZ(-posCorr.z());
00191 }
00192 } else {
00193
00194 EmcRecParameter& Para=EmcRecParameter::GetInstance();
00195
00196 if(Para.MethodMode()==1){
00197
00198 double IRShift = 0;
00199 HepPoint3D IR(0,0,IRShift);
00200 double theta12=-9999,theta03=-9999,theta23=-9999,theta14=-9999,delta=-9999;
00201 double phi=-9999,phi01=-9999,phi23=-9999,phi12=-9999,phi34=-9999,deltaphi=-9999;
00202 double xIn,xOut,deltaTheta,yIn,yOut,deltaPhi;
00203 double posphi=posCorr.phi();
00204
00205 double parTheta[5],parPhi[5];
00206
00207 HepPoint3D point0 ,point1 ,point2 ,point3 ,point4 ;
00208 point0 = iGeoSvc->GetCrystalPoint(id,0);
00209 point1 = iGeoSvc->GetCrystalPoint(id,1);
00210 point2 = iGeoSvc->GetCrystalPoint(id,2);
00211 point3 = iGeoSvc->GetCrystalPoint(id,3);
00212 point4 = iGeoSvc->GetCrystalPoint(id,4);
00213
00214 HepPoint3D center,center01, center23, center12,center03,center34,center14;
00215 center01 = (point0+point1)/2;
00216 center23 = (point2+point3)/2;
00217 center12 = (point1+point2)/2;
00218 center34 = (point3+point4)/2;
00219 center03 = (point0+point3)/2;
00220 center14 = (point1+point4)/2;
00221
00222 phi01 = center01.phi();
00223 phi23 = center23.phi();
00224 phi12 = center12.phi();
00225 phi34 = center34.phi();
00226
00227
00228 if( (thetaModule==1&&((phiModule+3)%4 == 0||(phiModule+2)%4 == 0))
00229 ||(thetaModule==3&&((phiModule+4)%5 == 0||(phiModule+3)%5 == 0||(phiModule+2)%5 == 0)))
00230 {
00231 center12 = center23;
00232 center03 = center14;
00233
00234 center01 = center12;
00235 center23 = center34;
00236
00237 phi01 = phi12;
00238 phi23 = phi34;
00239 }
00240
00241 if(phiModule==0){
00242 posphi = posphi;
00243 }else if( ((thetaModule==0||thetaModule==1)&&phiModule==63)
00244 ||((thetaModule==2||thetaModule==3)&&phiModule==79)
00245 ||((thetaModule==4||thetaModule==5)&&phiModule==95)
00246 ){
00247
00248 posphi = posphi+2*CLHEP::pi;
00249
00250 }else {
00251 posphi = posphi < 0 ? posphi+2*CLHEP::pi : posphi;
00252 }
00253
00254 phi01 = phi01 <0 ? phi01+2*CLHEP::pi : phi01;
00255 phi23 = phi23 <0 ? phi23+2*CLHEP::pi : phi23;
00256
00257 if(module==0){
00258
00259 IRShift = 10;
00260 center = center23;
00261 phi = phi23;
00262
00263 for(int i=0;i<5;i++) {
00264 parTheta[i] = Para.EastLinThetaPara(thetaModule,i);
00265 parPhi[i] = Para.EastLinPhiPara(0,i);
00266 }
00267
00268 }else if (module==2){
00269
00270 IRShift = -10;
00271 center = center01;
00272 phi = phi01;
00273
00274 for(int i=0;i<5;i++){
00275 parTheta[i] = Para.WestLinThetaPara(thetaModule,i);
00276 parPhi[i] = Para.WestLinPhiPara(0,i);
00277 }
00278
00279 }
00280
00281 theta12 = (center12 - IR).theta();
00282 theta03 = (center03 - IR).theta();
00283 delta = theta03 - theta12;
00284 xIn = (((posCorr-IR).theta() - theta12) - delta*0.5)/delta;
00285 xOut = parTheta[0]*atan(xIn*parTheta[1]-parTheta[4]) + parTheta[2]*xIn + parTheta[3] ;
00286 deltaTheta = xOut*delta;
00287
00288 deltaphi = fabs(phi23 - phi01);
00289 yIn = ((posphi - phi) - deltaphi*0.5)/deltaphi;
00290 yOut = parPhi[0]*atan(yIn*parPhi[1]-parPhi[4]) + parPhi[2]*yIn + parPhi[3] ;
00291 deltaPhi = yOut*deltaphi;
00292
00293 HepPoint3D rotateVector(0,1,0);
00294 rotateVector.rotateZ(center.phi());
00295 posCorr.setZ(posCorr.z()-IRShift);
00296 posCorr.rotate(-deltaTheta,rotateVector);
00297 posCorr.setZ(posCorr.z()+IRShift);
00298 posCorr.rotateZ(-deltaPhi);
00299
00300 }
00301 }
00302
00303
00304 if(Para.PosCorr()==1){ aShower.setPosition(posCorr);}
00305
00306
00307 if(aShower.energy()<1e-5) return;
00308
00309 double r,theta,phi;
00310 double dtheta,dphi,dx,dy,dz;
00311
00312 r = posCorr.mag();
00313 theta = posCorr.theta();
00314 phi = posCorr.phi();
00315
00316
00317
00318 if(aShower.energy()>0) {
00319 dtheta = Para.SigTheta(0)/sqrt(aShower.energy())+Para.SigTheta(1);
00320 dphi = Para.SigPhi(0)/sqrt(aShower.energy())+Para.SigPhi(1);
00321 }
00322 else {
00323 dtheta = 0;
00324 dphi = 0;
00325 }
00326
00327
00328 dx = fabs(iGeoSvc->GetBarrelR()*sin(phi)*dphi);
00329 dy = fabs(iGeoSvc->GetBarrelR()*cos(phi)*dphi);
00330 if(theta>0)
00331 dz = fabs(iGeoSvc->GetBarrelR()*dtheta/(sin(theta)*sin(theta)));
00332 else
00333 dz = 0;
00334
00335
00336 aShower.setDtheta(dtheta);
00337 aShower.setDphi(dphi);
00338
00339
00340
00341 HepSymMatrix matrix(3);
00342 matrix[0][0]=0;
00343 matrix[1][1]=dtheta*dtheta;
00344 matrix[2][2]=dphi*dphi;
00345 matrix[0][1]=0;
00346 matrix[0][2]=0;
00347 matrix[1][2]=0;
00348
00349
00350 HepMatrix matrix1(3,3),matrix2(3,3);
00351 matrix1[0][0]=sin(theta)*cos(phi);
00352 matrix1[0][1]=r*cos(theta)*cos(phi);
00353 matrix1[0][2]=-r*sin(theta)*sin(phi);
00354 matrix1[1][0]=sin(theta)*sin(phi);
00355 matrix1[1][1]=r*cos(theta)*sin(phi);
00356 matrix1[1][2]=r*sin(theta)*cos(phi);
00357 matrix1[2][0]=cos(theta);
00358 matrix1[2][1]=-r*sin(theta);
00359 matrix1[2][2]=0;
00360
00361
00362 for(int i=0;i<3;i++) {
00363 for(int j=0;j<3;j++) {
00364 matrix2[i][j]=matrix1[j][i];
00365 }
00366 }
00367
00368
00369 HepMatrix matrix4(3,3);
00370 matrix4=matrix1*matrix*matrix2;
00371
00372
00373 HepSymMatrix matrix5(3);
00374 for(int i=0;i<3;i++) {
00375 for(int j=0;j<=i;j++) {
00376 matrix5[i][j]=matrix4[i][j];
00377 }
00378 }
00379
00380
00381 aShower.setErrorMatrix(matrix5);
00382 }