00001
00002 #include "EmcRec/EmcRecShowerShape.h"
00003 #include "EmcRecEventModel/RecEmcEventModel.h"
00004
00005 #include "EmcRec/EmcRecParameter.h"
00006
00007 #include <vector>
00008 #include <complex>
00009
00010 void EmcRecShowerShape::CalculateMoment(RecEmcShower& aShower) const
00011 {
00012
00013 SecondMoment(aShower);
00014 LatMoment(aShower);
00015 A20Moment(aShower);
00016 A42Moment(aShower);
00017 }
00018
00019 void EmcRecShowerShape::SecondMoment(RecEmcShower& aShower) const
00020 {
00021 double etot=0;
00022 double sum=0;
00023 HepPoint3D center(aShower.position());
00024 RecEmcFractionMap fracMap=aShower.getFractionMap5x5();
00025
00026 RecEmcFractionMap::const_iterator it;
00027 for(it=fracMap.begin();
00028 it!=fracMap.end();
00029 it++){
00030 HepPoint3D pos(it->second.getFrontCenter());
00031
00033 EmcRecParameter& Para=EmcRecParameter::GetInstance();
00034 if(Para.DataMode()==1) {
00035
00036 RecEmcID id = it->second.getCellId();
00037 const unsigned int module = EmcID::barrel_ec(id);
00038 const unsigned int thetaModule = EmcID::theta_module(id);
00039 const unsigned int phiModule = EmcID::phi_module(id);
00040 double lengthemc;
00041 lengthemc = abs(pos.z()/cos(pos.theta()));
00042 double posDataCorPar;
00043 if(module==1) posDataCorPar=Para.BarrPosDataCor(thetaModule,phiModule);
00044 if(module==0) posDataCorPar=Para.EastPosDataCor(thetaModule,phiModule);
00045 if(module==2) posDataCorPar=Para.WestPosDataCor(thetaModule,phiModule);
00046
00047 pos.setTheta(pos.theta()-posDataCorPar/lengthemc);
00048
00049 }
00050
00052
00053 etot+=it->second.getEnergy()*it->second.getFraction();
00054 sum+=it->second.getEnergy()*it->second.getFraction()*pos.distance2(center);
00055 }
00056
00057 if(etot>0) {
00058 sum/=etot;
00059 }
00060
00061 aShower.setSecondMoment(sum);
00062 }
00063
00064 void EmcRecShowerShape::LatMoment(RecEmcShower& aShower) const
00065 {
00066 RecEmcFractionMap fracMap=aShower.getFractionMap5x5();
00067 if(fracMap.size()<2) {
00068 aShower.setLatMoment(0);
00069 return;
00070 }
00071
00072 vector<RecEmcFraction> aFractionVec;
00073 RecEmcFractionMap::const_iterator it;
00074 for(it=fracMap.begin();
00075 it!=fracMap.end();
00076 it++){
00077 aFractionVec.push_back(it->second);
00078 }
00079
00080
00081 partial_sort(aFractionVec.begin(),
00082 aFractionVec.begin()+2,
00083 aFractionVec.end(),
00084 greater<RecEmcFraction>());
00085
00086
00087 vector<RecEmcFraction>::iterator iVec;
00088 double numerator=0;
00089 double denominator=0;
00090 int n=0;
00091 for(iVec=aFractionVec.begin();
00092 iVec!=aFractionVec.end();
00093 iVec++) {
00094 n++;
00095 HepPoint3D pos((*iVec).getFrontCenter());
00096
00098 EmcRecParameter& Para=EmcRecParameter::GetInstance();
00099 if(Para.DataMode()==1) {
00100
00101 RecEmcID id = (*iVec).getCellId();
00102 const unsigned int module = EmcID::barrel_ec(id);
00103 const unsigned int thetaModule = EmcID::theta_module(id);
00104 const unsigned int phiModule = EmcID::phi_module(id);
00105 double lengthemc;
00106 lengthemc = abs(pos.z()/cos(pos.theta()));
00107 double posDataCorPar;
00108 if(module==1) posDataCorPar=Para.BarrPosDataCor(thetaModule,phiModule);
00109 if(module==0) posDataCorPar=Para.EastPosDataCor(thetaModule,phiModule);
00110 if(module==2) posDataCorPar=Para.WestPosDataCor(thetaModule,phiModule);
00111
00112 pos.setTheta(pos.theta()-posDataCorPar/lengthemc);
00113
00114 }
00115
00117
00118 double r=pos.mag()*sin(aShower.position().angle(pos));
00119
00120 double energy=(*iVec).getEnergy()*(*iVec).getFraction();
00121 if(n<3) {
00122 denominator+=5.2*5.2*energy;
00123 } else {
00124 numerator+=r*r*energy;
00125 denominator+=r*r*energy;
00126 }
00127 }
00128 double lat=-99;
00129 if(denominator>0) lat=numerator/denominator;
00130
00131 aShower.setLatMoment(lat);
00132
00133 }
00134
00135 void EmcRecShowerShape::A20Moment(RecEmcShower& aShower) const
00136 {
00137 double a20=0;
00138 const double R0=15.6;
00139 Hep3Vector r0(aShower.position());
00140 RecEmcFractionMap fracMap=aShower.getFractionMap5x5();
00141 RecEmcFractionMap::const_iterator it;
00142 for(it=fracMap.begin();
00143 it!=fracMap.end();
00144 it++){
00145 double energy=it->second.getEnergy()*it->second.getFraction();
00146 HepPoint3D pos(it->second.getFrontCenter());
00147
00148 Hep3Vector r=pos-r0;
00149 r=r-r.dot(r0)*r0/(r0.mag()*r0.mag());
00150
00151 a20+=(energy/aShower.e5x5())*(2*pow(r.mag()/R0,2.)-1);
00152 }
00153
00154 aShower.setA20Moment(fabs(a20));
00155 }
00156
00157 void EmcRecShowerShape::A42Moment(RecEmcShower& aShower) const
00158 {
00159 complex<double> a42(0.,0.);
00160 const double R0=15.6;
00161 Hep3Vector r0(aShower.position());
00162 RecEmcFractionMap fracMap=aShower.getFractionMap5x5();
00163 RecEmcFractionMap::const_iterator it;
00164 for(it=fracMap.begin();
00165 it!=fracMap.end();
00166 it++){
00167 double energy=it->second.getEnergy()*it->second.getFraction();
00168 HepPoint3D pos(it->second.getFrontCenter());
00169
00170 Hep3Vector r=pos-r0;
00171 r=r-r.dot(r0)*r0/(r0.mag()*r0.mag());
00172
00173 complex<double> a(0.,2.*r.phi());
00174
00175 a42+=(energy/aShower.e5x5())*(4.*pow(r.mag()/R0,4.)-3.*pow(r.mag()/R0,2.))*exp(a);
00176 }
00177
00178 aShower.setA42Moment(abs(a42));
00179 }