00001
00002 #include "EmcRec/EmcRecShowerEnergy.h"
00003 #include "EmcRec/EmcRecNeighbor.h"
00004 #include "EmcRec/EmcRecParameter.h"
00005
00006 void EmcRecShowerEnergy::Energy(RecEmcShower& aShower)
00007 {
00008 RecEmcFractionMap::const_iterator cit;
00009 RecEmcEnergy e1=0;
00010 RecEmcEnergy e9=0;
00011 RecEmcEnergy e25=0;
00012 RecEmcEnergy elepton=0;
00013 RecEmcEnergy eall=0;
00014
00015 EmcRecNeighbor nhb;
00016
00017 RecEmcID CellId=aShower.getShowerId();
00018 int module=EmcID::barrel_ec(CellId);
00019 RecEmcIDVector NearCell=nhb.GetNeighbors(CellId);
00020 RecEmcIDVector NextNearCell=nhb.GetNextNeighbors(CellId);
00021 RecEmcIDVector tmpNearCell;
00022 RecEmcIDVector tmpNextNearCell;
00023 i_RecEmcIDVector pNearCell;
00024 i_RecEmcIDVector pNextNearCell;
00025 vector<RecEmcEnergy> eVec;
00026 vector<RecEmcEnergy>::const_iterator ciVec;
00027
00028 tmpNearCell.push_back(CellId);
00029 tmpNextNearCell.push_back(CellId);
00030
00031 cit=aShower.Find(CellId);
00032
00033 e1=(cit->second.getEnergy())*(cit->second.getFraction());
00034 e9+=(cit->second.getEnergy())*(cit->second.getFraction());
00035 e25+=(cit->second.getEnergy())*(cit->second.getFraction());
00036
00037
00038 for(pNearCell=NearCell.begin();
00039 pNearCell!=NearCell.end();
00040 pNearCell++) {
00041 cit=aShower.Find(*pNearCell);
00042 if(cit!=aShower.End()) {
00043 tmpNearCell.push_back(*pNearCell);
00044 tmpNextNearCell.push_back(*pNearCell);
00045 e9+=cit->second.getEnergy()*cit->second.getFraction();
00046 e25+=cit->second.getEnergy()*cit->second.getFraction();
00047 }
00048 }
00049
00050
00051 for(pNextNearCell=NextNearCell.begin();
00052 pNextNearCell!=NextNearCell.end();
00053 pNextNearCell++) {
00054 cit=aShower.Find(*pNextNearCell);
00055 if(cit!=aShower.End()) {
00056 tmpNextNearCell.push_back(*pNextNearCell);
00057 e25+=cit->second.getEnergy()*cit->second.getFraction();
00058 }
00059 }
00060
00061
00062 for(cit=aShower.Begin();cit!=aShower.End();++cit) {
00063 eall+=(cit->second.getEnergy())*(cit->second.getFraction());
00064 eVec.push_back(cit->second.getEnergy()*cit->second.getFraction());
00065 }
00066
00067
00068 int nHit,n;
00069 EmcRecParameter& Para=EmcRecParameter::GetInstance();
00070 nHit=(int)(Para.HitNb(0)*log(Para.HitNb(1)*e9+Para.HitNb(2)));
00071 n=0;
00072
00073
00074 sort(eVec.begin(), eVec.end(), greater<RecEmcEnergy>());
00075
00076 for(ciVec=eVec.begin();ciVec!=eVec.end();ciVec++) {
00077 if(n<nHit) {
00078 elepton+=*ciVec;
00079 n++;
00080 }
00081 }
00082
00083
00084
00085
00086 int getthetaid = EmcID::theta_module(CellId);
00087 int getmodule = EmcID::barrel_ec(CellId);
00088 if(getthetaid>21)getthetaid=43-getthetaid;
00089 if(getmodule==1)getthetaid=getthetaid+6;
00090 double dthetaid=double(getthetaid);
00091 double eCorr = Para.ECorrMC(e25,dthetaid);
00092
00093
00094 RecEmcEnergy de,de1,de2,de3;
00095 de1 = Para.SigE(0)/eCorr;
00096 de2 = Para.SigE(1)/pow(eCorr,0.25);
00097 de3 = Para.SigE(2);
00098 de = sqrt(de1*de1+de2*de2+de3*de3)*perCent*eCorr;
00099
00100 double err = Para.ErrMC(e25,dthetaid);
00101 if(err>0) de = err*e25;
00102
00103 aShower.setTrackId(-1);
00104 aShower.setCellId(CellId);
00105 aShower.setModule(module);
00106 aShower.setNumHits(aShower.getSize());
00107 aShower.setDE(de);
00108 aShower.CellId3x3(tmpNearCell);
00109 aShower.CellId5x5(tmpNextNearCell);
00110 aShower.setESeed(e1);
00111 aShower.setE3x3(e9);
00112 aShower.setE5x5(e25);
00113 aShower.ELepton(elepton);
00114 aShower.EAll(eall);
00115 aShower.setEnergy(eCorr);
00116
00117
00118
00119
00120
00121
00122
00123 }
00124
00125 RecEmcEnergy EmcRecShowerEnergy::ECorrection(const RecEmcEnergy eIn)
00126 {
00127 if(eIn>3.) return eIn;
00128
00129 EmcRecParameter& Para=EmcRecParameter::GetInstance();
00130
00131 RecEmcEnergy eOut=0;
00132 double par[4];
00133 for(int i=0;i<4;i++) {
00134 par[i]=Para.ECorr(i);
00135 }
00136
00137 eOut = eIn/(par[0]+par[1]*eIn+par[2]*eIn*eIn+par[3]*eIn*eIn*eIn);
00138 return eOut;
00139 }
00140
00141 RecEmcEnergy EmcRecShowerEnergy::ECorrTheta(const RecEmcEnergy eIn, const RecEmcID &id)
00142 {
00143 EmcRecParameter& Para=EmcRecParameter::GetInstance();
00144 RecEmcEnergy eOut=eIn;
00145
00146 unsigned int npart = EmcID::barrel_ec(id);
00147 unsigned int ntheta = EmcID::theta_module(id);
00148
00149 if(npart==1) {
00150 eOut *= 1.843/Para.Peak(ntheta);
00151 } else if(npart==0) {
00152 eOut *= 1.843/Para.Peak(ntheta+44);
00153 } else if(npart==2) {
00154 eOut *= 1.843/Para.Peak(ntheta+50);
00155 }
00156
00157 return eOut;
00158 }