00001
00002
00003
00004
00005
00006 #include "EmcRec/EmcRecShowerPosLoglin.h"
00007
00008 #include <iostream>
00009 #include <fstream>
00010
00011
00012 #include "EmcRecGeoSvc/EmcRecGeoSvc.h"
00013 #include "GaudiKernel/Bootstrap.h"
00014 #include "GaudiKernel/ISvcLocator.h"
00015
00016 #include "EmcRec/EmcRecParameter.h"
00017
00018
00019 void EmcRecShowerPosLoglin::Position(RecEmcShower& aShower)
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
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
00039
00040 double a0=4;
00041 double offset;
00042
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
00085
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
00095 wsum+=w1;
00096 pos=iGeoSvc->GetCFrontCenter(cit->second.getCellId());
00097 possum+=pos*w1;
00098 } else{
00099
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
00111 }
00112 pos=possum/wsum;
00113
00114
00115
00116 if(Para.PosCorr()==0){
00117
00118 aShower.setPosition(pos);
00119 }
00120
00121
00122
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) {
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
00157
00158
00159
00160
00161
00162
00163
00164 }
00165
00166 HepPoint3D IR(0,0,IRShift);
00167 HepPoint3D center01, center23;
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);
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
00182
00183
00184
00185
00186
00187 double yIn,yOut,deltaPhi;
00188
00189
00190
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);
00206 rotateVector.rotateZ(center01.phi());
00207 posCorr.setZ(posCorr.z()-IRShift);
00208 posCorr.rotate(-deltaTheta,rotateVector);
00209 posCorr.setZ(posCorr.z()+IRShift);
00210
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
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
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
00250
00251 aShower.setDtheta(dtheta);
00252 aShower.setDphi(dphi);
00253
00254
00255
00256 HepSymMatrix matrix(3);
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
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
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
00283
00284 HepMatrix matrix4(3,3);
00285 matrix4=matrix1*matrix*matrix2;
00286
00287
00288 HepSymMatrix matrix5(3);
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
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 }
00302