#include <EmcRecShowerPosLin.h>
Inheritance diagram for EmcRecShowerPosLin:
Public Member Functions | |
| EmcRecShowerPosLin () | |
| EmcRecShowerPosLin () | |
| virtual void | Position (RecEmcShower &aShower) |
| virtual void | Position (RecEmcShower &aShower) |
| virtual | ~EmcRecShowerPosLin () |
| virtual | ~EmcRecShowerPosLin () |
|
|
00018 {}
|
|
|
00019 {}
|
|
|
00018 {}
|
|
|
00019 {}
|
|
|
Implements EmcRecShowerPosAbs. |
|
|
Implements EmcRecShowerPosAbs. 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 //cout<<"EmcRecShowerPosLin::Position Begin"<<endl;
00025
00026 IEmcRecGeoSvc* iGeoSvc;
00027 ISvcLocator* svcLocator = Gaudi::svcLocator();
00028 StatusCode sc = svcLocator->service("EmcRecGeoSvc",iGeoSvc);
00029 if(sc!=StatusCode::SUCCESS) {
00030 cout<<"Error: Can't get EmcRecGeoSvc"<<endl;
00031 }
00032 etot=0;
00033 //cout<<etot<<endl;
00034 for(cit=aShower.Begin();cit!=aShower.End();cit++){
00035 etot+=(cit->second.getEnergy()*cit->second.getFraction());
00036 //cout<<etot<<endl;
00037 }
00038 wsum=0;
00039 for(cit=aShower.Begin();cit!=aShower.End();cit++){
00040 w=(cit->second.getEnergy()*cit->second.getFraction());
00041 //cout<<w<<endl;
00042 wsum+=w;
00043 pos=iGeoSvc->GetCFrontCenter(cit->second.getCellId());
00044 possum+=pos*w;
00045 }
00046
00047 if(wsum<=0) {
00048 cout<<"[[[[[[[[[[[[[[["<<wsum;
00049 }
00050
00051 pos=possum/wsum;
00052 //aShower.setPosition(pos);
00053
00054 //----------------------------------------------------------------------
00055 //position correction
00056 RecEmcID id = aShower.getShowerId();
00057 const unsigned int module = EmcID::barrel_ec(id);;
00058 const unsigned int thetaModule = EmcID::theta_module(id);
00059 const unsigned int phiModule = EmcID::phi_module(id);
00060 HepPoint3D posCorr(pos);
00061
00062 if(module==1) { //barrel
00063 unsigned int theModule;
00064 if(thetaModule>21) {
00065 theModule = 43 - thetaModule;
00066 id = EmcID::crystal_id(module,theModule,phiModule);
00067 posCorr.setZ(-posCorr.z());
00068 } else {
00069 theModule = thetaModule;
00070 }
00071 //cout<<"EmcID: "<<id<<" theta: "<<theModule<<endl;
00072
00073 double IRShift, parTheta[5];
00074 if(theModule==21) {
00075 IRShift = 0;
00076 double par[5] = { 1.698, -1.553, 0.9384, 0.1193, -0.01916 };
00077 for(int i=0;i<5;i++) {
00078 parTheta[i] = par[i];
00079 }
00080 } else if(theModule==20) {
00081 IRShift = 2.5;
00082 double par[5] = { 1.593, -1.582, 0.8881, 0.3298, -0.08856 };
00083 for(int i=0;i<5;i++) {
00084 parTheta[i] = par[i];
00085 }
00086 } else {
00087 IRShift = 5.;
00088 double par[5] = { 1.605, -1.702, 0.8433, 0.3072, -0.1809 };
00089 for(int i=0;i<5;i++) {
00090 parTheta[i] = par[i];
00091 }
00092 }
00093
00094 HepPoint3D IR(0,0,IRShift);
00095 HepPoint3D center01, center23; //center of point0,1 and point2,3 in crystal
00096 center01 = (iGeoSvc->GetCrystalPoint(id,0)+iGeoSvc->GetCrystalPoint(id,1))/2;
00097 center23 = (iGeoSvc->GetCrystalPoint(id,2)+iGeoSvc->GetCrystalPoint(id,3))/2;
00098 //cout<<"before correction: "<<(posCorr-IR).theta()<<"\t"<<posCorr.phi()<<endl;
00099
00100 double theta01,theta23,length,d;
00101 theta01 = (center01-IR).theta();
00102 theta23 = (center23-IR).theta();
00103 length = (center01-IR).mag();
00104 d = length*tan(theta01-theta23); //length in x direction
00105 //cout<<"theta01: "<<theta01<<"\t"<<" theta23: "<<theta23<<"\t"
00106 //<<" length: "<<length<<" d: "<<d<<endl;
00107
00108 double xIn,xOut,deltaTheta;
00109 xIn = length*tan(theta01-(posCorr-IR).theta())-d/2;
00110 xOut = xIn-(parTheta[0]*atan(xIn*parTheta[1]-parTheta[4])+parTheta[2]*xIn+parTheta[3]);
00111 deltaTheta = atan((xOut+d/2)/length)-atan((xIn+d/2)/length);
00112 //cout<<"xIn: "<<xIn<<" xOut: "<<xOut<<" deltaTheta: "<<deltaTheta<<endl;
00113
00114 double parPhi1[5] = { 0.7185, -2.539, 0.6759, 0.3472, -0.3108 }; //e-
00115 double parPhi2[5] = { 0.5781, -2.917, 0.5516, -0.5745, 0.3777 }; //e+
00116 double parPhi[5];// = { 0.7723, -3.1, 0.6992, -0.1441, -0.01012 }; //by photon, not used
00117 for(int i=0;i<5;i++) {
00118 parPhi[i] = (parPhi1[i]+parPhi2[i])/2; //average of e- and e+
00119 }
00120
00121 double yIn,yOut,deltaPhi;
00122 yIn = posCorr.phi() < 0 ? posCorr.phi()*180./CLHEP::pi+360.-phiModule*3.-1.1537
00123 : posCorr.phi()*180./CLHEP::pi-phiModule*3.-1.1537;
00124 yOut = parPhi[0]*atan(yIn*parPhi[1]-parPhi[4])+parPhi[2]*yIn+parPhi[3];
00125 deltaPhi = yOut*CLHEP::pi/180.;
00126 //cout<<"yIn="<<yIn<<"\tyOut="<<yOut<<"\tdeltaPhi="<<deltaPhi<<endl;
00127
00128 HepPoint3D rotateVector(0,1,0); //y axis
00129 rotateVector.rotateZ(posCorr.phi());
00130 posCorr.setZ(posCorr.z()-IRShift);
00131 posCorr.rotate(-deltaTheta,rotateVector);
00132 posCorr.setZ(posCorr.z()+IRShift);
00133 //posCorr.rotateZ(-deltaPhi);
00134 posCorr.rotateZ(-0.002994);
00135 //cout<<"after correction: "<<(posCorr-IR).theta()<<"\t"<<posCorr.phi()<<endl;
00136
00137 if(thetaModule>21) {
00138 posCorr.setZ(-posCorr.z());
00139 }
00140 }
00141
00142 aShower.setPosition(posCorr);
00143
00144 //----------------------------------------------------------------------
00145 if(aShower.energy()<1e-5) return;
00146
00147 double r,theta,phi;
00148 double dtheta,dphi,dx,dy,dz;
00149
00150 r = posCorr.mag();
00151 theta = posCorr.theta();
00152 phi = posCorr.phi();
00153 //cout<<"theta: "<<theta<<" phi: "<<phi<<endl;
00154
00155 EmcRecParameter& Para=EmcRecParameter::GetInstance();
00156 if(aShower.energy()>0) {
00157 dtheta = Para.SigTheta(0)/sqrt(aShower.energy())+Para.SigTheta(1);
00158 dphi = Para.SigPhi(0)/sqrt(aShower.energy())+Para.SigPhi(1);
00159 }
00160 else {
00161 dtheta = 0;
00162 dphi = 0;
00163 }
00164 //cout<<"dtheta: "<<dtheta<<" dphi: "<<dphi<<endl;
00165
00166 dx = fabs(iGeoSvc->GetBarrelR()*sin(phi)*dphi);
00167 dy = fabs(iGeoSvc->GetBarrelR()*cos(phi)*dphi);
00168 if(theta>0)
00169 dz = fabs(iGeoSvc->GetBarrelR()*dtheta/(sin(theta)*sin(theta)));
00170 else
00171 dz = 0;
00172 //cout<<" dx: "<<dx<<" dy: "<<dy<<" dz: "<<dz<<endl;
00173
00174 aShower.setDtheta(dtheta);
00175 aShower.setDphi(dphi);
00176
00177 //----------------------------------------------------------------------
00178 // Error matrix
00179 HepSymMatrix matrix(3); //error matrix of r, theta, phi
00180 matrix[0][0]=0;
00181 matrix[1][1]=dtheta*dtheta;
00182 matrix[2][2]=dphi*dphi;
00183 matrix[0][1]=0;
00184 matrix[0][2]=0;
00185 matrix[1][2]=0;
00186 //cout<<"matrix: \n"<<matrix;
00187
00188 HepMatrix matrix1(3,3),matrix2(3,3);
00189 matrix1[0][0]=sin(theta)*cos(phi);
00190 matrix1[0][1]=r*cos(theta)*cos(phi);
00191 matrix1[0][2]=-r*sin(theta)*sin(phi);
00192 matrix1[1][0]=sin(theta)*sin(phi);
00193 matrix1[1][1]=r*cos(theta)*sin(phi);
00194 matrix1[1][2]=r*sin(theta)*cos(phi);
00195 matrix1[2][0]=cos(theta);
00196 matrix1[2][1]=-r*sin(theta);
00197 matrix1[2][2]=0;
00198 //cout<<"matrix1: \n"<<matrix1;
00199
00200 for(int i=0;i<3;i++) {
00201 for(int j=0;j<3;j++) {
00202 matrix2[i][j]=matrix1[j][i];
00203 }
00204 }
00205 //cout<<"matrix2: \n"<<matrix2;
00206
00207 HepMatrix matrix4(3,3);
00208 matrix4=matrix1*matrix*matrix2;
00209 //cout<<"matrix4: \n"<<matrix4;
00210
00211 HepSymMatrix matrix5(3); //error matrix of x, y, z
00212 for(int i=0;i<3;i++) {
00213 for(int j=0;j<=i;j++) {
00214 matrix5[i][j]=matrix4[i][j];
00215 }
00216 }
00217 //cout<<"matrix5: \n"<<matrix5;
00218
00219 aShower.setErrorMatrix(matrix5);
00220 }
|
1.3.9.1