#include <EmcRecShowerPosLog.h>
Inheritance diagram for EmcRecShowerPosLog:
Public Member Functions | |
| EmcRecShowerPosLog () | |
| EmcRecShowerPosLog () | |
| virtual void | Position (RecEmcShower &aShower) |
| virtual void | Position (RecEmcShower &aShower) |
| virtual | ~EmcRecShowerPosLog () |
| virtual | ~EmcRecShowerPosLog () |
|
|
00018 {}
|
|
|
00019 {}
|
|
|
00018 {}
|
|
|
00019 {}
|
|
|
Implements EmcRecShowerPosAbs. |
|
|
Implements EmcRecShowerPosAbs. 00017 {
00018 RecEmcFractionMap::const_iterator cit;
00019 HepPoint3D pos(0,0,0);
00020 HepPoint3D possum(0,0,0);
00021 double w,wsum=0;
00022
00023 //cout<<"EmcRecShowerPosLog::Position Begin"<<endl;
00024
00025 IEmcRecGeoSvc* iGeoSvc;
00026 ISvcLocator* svcLocator = Gaudi::svcLocator();
00027 StatusCode sc = svcLocator->service("EmcRecGeoSvc",iGeoSvc);
00028 if(sc!=StatusCode::SUCCESS) {
00029 cout<<"Error: Can't get EmcRecGeoSvc"<<endl;
00030 }
00031
00032 EmcRecParameter& Para=EmcRecParameter::GetInstance();
00033 double offset=Para.LogPosOffset();
00034 //cout<<offset<<endl;
00035
00036
00037 if(Para.PositionMode2()=="all") {
00038 double etot=aShower.getEAll();
00039 for(cit=aShower.Begin();cit!=aShower.End();cit++){
00040 w=offset+log((cit->second.getEnergy()/etot)*cit->second.getFraction());
00041 if(w>0){
00042 wsum+=w;
00043 pos=iGeoSvc->GetCFrontCenter(cit->second.getCellId());
00044 possum+=pos*w;
00045 }
00046 }
00047 } else if(Para.PositionMode2()=="3x3") {
00048 double e3x3=aShower.e3x3();
00049 RecEmcFractionMap fracMap3x3=aShower.getFractionMap3x3();
00050 for(cit=fracMap3x3.begin();
00051 cit!=fracMap3x3.end();
00052 cit++) {
00053 w=offset+log((cit->second.getEnergy()/e3x3)*cit->second.getFraction());
00054 if(w>0){
00055 wsum+=w;
00056 pos=iGeoSvc->GetCFrontCenter(cit->second.getCellId());
00057 possum+=pos*w;
00058 }
00059 }
00060 } else {
00061 double e5x5=aShower.e5x5();
00062 RecEmcFractionMap fracMap5x5=aShower.getFractionMap5x5();
00063 for(cit=fracMap5x5.begin();
00064 cit!=fracMap5x5.end();
00065 cit++) {
00066 w=offset+log((cit->second.getEnergy()/e5x5)*cit->second.getFraction());
00067 if(w>0){
00068 wsum+=w;
00069 pos=iGeoSvc->GetCFrontCenter(cit->second.getCellId());
00070 possum+=pos*w;
00071 }
00072 }
00073 }
00074 if(wsum<=0) {
00075 cout<<"[[[[[[[[[[[[[[["<<wsum<<endl;
00076 }
00077 pos=possum/wsum;
00078 //aShower.setPosition(pos);
00079
00080 //----------------------------------------------------------------------
00081 //position correction
00082 RecEmcID id = aShower.getShowerId();
00083 const unsigned int module = EmcID::barrel_ec(id);;
00084 const unsigned int thetaModule = EmcID::theta_module(id);
00085 const unsigned int phiModule = EmcID::phi_module(id);
00086 HepPoint3D posCorr(pos);
00087
00088 if(module==1) { //barrel
00089 unsigned int theModule;
00090 if(thetaModule>21) {
00091 theModule = 43 - thetaModule;
00092 id = EmcID::crystal_id(module,theModule,phiModule);
00093 posCorr.setZ(-posCorr.z());
00094 } else {
00095 theModule = thetaModule;
00096 }
00097
00098 double IRShift, parTheta[5];
00099 if(theModule==21) {
00100 IRShift = 0;
00101 double par[5] = { 51.97, -0.1209, 6.175, -0.1431, -0.005497 };
00102 for(int i=0;i<5;i++) {
00103 parTheta[i] = par[i];
00104 }
00105 } else if(theModule==20) {
00106 IRShift = 2.5;
00107 double par[5] = { 56.31, -0.1135, 6.312, -1.216, -0.02978 };
00108 for(int i=0;i<5;i++) {
00109 parTheta[i] = par[i];
00110 }
00111 } else {
00112 IRShift = 5.;
00113 double par[5] = { 10.65, -0.2257, 2.216, -0.1562, -0.05552 };
00114 for(int i=0;i<5;i++) {
00115 parTheta[i] = par[i];
00116 }
00117 }
00118
00119 HepPoint3D IR(0,0,IRShift);
00120 HepPoint3D center01, center23; //center of point0,1 and point2,3 in crystal
00121 center01 = (iGeoSvc->GetCrystalPoint(id,0)+iGeoSvc->GetCrystalPoint(id,1))/2;
00122 center23 = (iGeoSvc->GetCrystalPoint(id,2)+iGeoSvc->GetCrystalPoint(id,3))/2;
00123
00124 double theta01,theta23,length,d;
00125 theta01 = (center01-IR).theta();
00126 theta23 = (center23-IR).theta();
00127 length = (center01-IR).mag();
00128 d = length*tan(theta01-theta23); //length in x direction
00129
00130 double xIn,xOut,deltaTheta;
00131 xIn = length*tan(theta01-(posCorr-IR).theta())-d/2;
00132 xOut = xIn-( parTheta[0]*atan(xIn*parTheta[1]-parTheta[4])
00133 + parTheta[2]*xIn + parTheta[3] );
00134 deltaTheta = atan((xOut+d/2)/length)-atan((xIn+d/2)/length);
00135 //cout<<"xIn: "<<xIn<<" xOut: "<<xOut<<" deltaTheta: "<<deltaTheta<<endl;
00136 //cout<<"dx: "<<parTheta[0]*atan(xIn*parTheta[1]-parTheta[4])
00137 // + parTheta[2]*xIn + parTheta[3] <<endl;
00138
00139 //obtained by photon, not used, just for comparison
00140 double parPhi[5] = { 51.1, -0.1499, 7.62, 0.1301, 0.005581 };
00141 double yIn,yOut,deltaPhi;
00142 yIn = posCorr.phi() < 0 ? posCorr.phi()*180./CLHEP::pi+360.-phiModule*3.-1.1537
00143 : posCorr.phi()*180./CLHEP::pi-phiModule*3.-1.1537;
00144 yOut = parPhi[0]*atan(yIn*parPhi[1]-parPhi[4])+parPhi[2]*yIn+parPhi[3];
00145 deltaPhi = yOut*CLHEP::pi/180.;
00146 //
00147
00148 HepPoint3D rotateVector(0,1,0); //y axis
00149 rotateVector.rotateZ(center01.phi());
00150 posCorr.setZ(posCorr.z()-IRShift);
00151 posCorr.rotate(-deltaTheta,rotateVector);
00152 posCorr.setZ(posCorr.z()+IRShift);
00153 //phi parameter is gotten by the average of e+ e- peak
00154 posCorr.rotateZ(-0.002994);
00155 //posCorr.rotateZ(-deltaPhi);
00156
00157 if(thetaModule>21) {
00158 posCorr.setZ(-posCorr.z());
00159 }
00160 }
00161
00162 aShower.setPosition(posCorr);
00164 if(aShower.energy()<1e-5) return;
00165
00166 double r,theta,phi;
00167 double dtheta,dphi,dx,dy,dz;
00168
00169 r = posCorr.mag();
00170 theta = posCorr.theta();
00171 phi = posCorr.phi();
00172 //cout<<"energy: "<<aShower.energy()<<" theta: "<<theta<<" phi: "<<phi<<endl;
00173
00174 if(aShower.energy()>0) {
00175 dtheta = Para.SigTheta(0)/sqrt(aShower.energy())+Para.SigTheta(1);
00176 dphi = Para.SigPhi(0)/sqrt(aShower.energy())+Para.SigPhi(1);
00177 }
00178 else {
00179 dtheta = 0;
00180 dphi = 0;
00181 }
00182 //cout<<"dtheta: "<<dtheta<<" dphi: "<<dphi<<endl;
00183
00184 dx = fabs(iGeoSvc->GetBarrelR()*sin(phi)*dphi);
00185 dy = fabs(iGeoSvc->GetBarrelR()*cos(phi)*dphi);
00186 if(theta>0)
00187 dz = fabs(iGeoSvc->GetBarrelR()*dtheta/(sin(theta)*sin(theta)));
00188 else
00189 dz = 0;
00190 //cout<<" dx: "<<dx<<" dy: "<<dy<<" dz: "<<dz<<endl;
00191
00192 aShower.setDtheta(dtheta);
00193 aShower.setDphi(dphi);
00194
00195 //----------------------------------------------------------------------
00196 // Error matrix
00197 HepSymMatrix matrix(3); //error matrix of r, theta, phi
00198 matrix[0][0]=0;
00199 matrix[1][1]=dtheta*dtheta;
00200 matrix[2][2]=dphi*dphi;
00201 matrix[0][1]=0;
00202 matrix[0][2]=0;
00203 matrix[1][2]=0;
00204 //cout<<"matrix:\n"<<matrix<<endl;
00205
00206 HepMatrix matrix1(3,3),matrix2(3,3);
00207 matrix1[0][0]=sin(theta)*cos(phi);
00208 matrix1[0][1]=r*cos(theta)*cos(phi);
00209 matrix1[0][2]=-r*sin(theta)*sin(phi);
00210 matrix1[1][0]=sin(theta)*sin(phi);
00211 matrix1[1][1]=r*cos(theta)*sin(phi);
00212 matrix1[1][2]=r*sin(theta)*cos(phi);
00213 matrix1[2][0]=cos(theta);
00214 matrix1[2][1]=-r*sin(theta);
00215 matrix1[2][2]=0;
00216 //cout<<"matrix1:\n"<<matrix1<<endl;
00217
00218 for(int i=0;i<3;i++) {
00219 for(int j=0;j<3;j++) {
00220 matrix2[i][j]=matrix1[j][i];
00221 }
00222 }
00223 //cout<<"matrix2:\n"<<matrix2<<endl;
00224
00225 HepMatrix matrix4(3,3);
00226 matrix4=matrix1*matrix*matrix2;
00227 //cout<<"matrix4:\n"<<matrix4<<endl;
00228
00229 HepSymMatrix matrix5(3); //error matrix of x, y, z
00230 for(int i=0;i<3;i++) {
00231 for(int j=0;j<=i;j++) {
00232 matrix5[i][j]=matrix4[i][j];
00233 }
00234 }
00235 //cout<<"matrix5:\n"<<matrix5<<endl;
00236
00237 aShower.setErrorMatrix(matrix5);
00238
00239 if(matrix5[0][0]>0) dx=sqrt(matrix5[0][0]);
00240 if(matrix5[1][1]>0) dy=sqrt(matrix5[1][1]);
00241 if(matrix5[2][2]>0) dz=sqrt(matrix5[2][2]);
00242 }
|
1.3.9.1