00001
00002
00003
00004
00005
00006
00007 #include <iostream>
00008
00009 #include "EmcRec/EmcRecSplitWeighted.h"
00010 #include "EmcRec/EmcRecShowerEnergy.h"
00011 #include "EmcRec/EmcRecShowerPosLin.h"
00012 #include "EmcRec/EmcRecShowerPosLog.h"
00013 #include "EmcRec/EmcRecShowerPosLoglin.h"
00014 #include "EmcRec/EmcRecShowerPosLinShMax.h"
00015 #include "EmcRec/EmcRecShowerPosLogShMax.h"
00016 #include "EmcRec/EmcRecShowerShape.h"
00017 #include "EmcRec/EmcRecParameter.h"
00018 #include "EmcRecGeoSvc/EmcRecGeoSvc.h"
00019 #include "GaudiKernel/Bootstrap.h"
00020 #include "GaudiKernel/ISvcLocator.h"
00021
00022
00023 EmcRecSplitWeighted::EmcRecSplitWeighted()
00024 {
00025 fShowerE=new EmcRecShowerEnergy;
00026 fShowerShape=new EmcRecShowerShape;
00027
00028 EmcRecParameter& Para=EmcRecParameter::GetInstance();
00029
00030
00031
00032
00033
00034
00035 if(Para.PositionMode1()=="lin") {
00036 fShowerPos=new EmcRecShowerPosLin;
00037 } else if(Para.PositionMode1()=="log"){
00038 fShowerPos=new EmcRecShowerPosLog;
00039 } else if(Para.PositionMode1()=="loglin"){
00040 fShowerPos=new EmcRecShowerPosLoglin;
00041 } else if(Para.PositionMode1()=="linsh"){
00042 fShowerPos=new EmcRecShowerPosLinShMax;
00043 } else if(Para.PositionMode1()=="logsh"){
00044 fShowerPos=new EmcRecShowerPosLogShMax;
00045 }
00046
00047 }
00048
00049
00050 EmcRecSplitWeighted::~EmcRecSplitWeighted()
00051 {
00052 delete fShowerE;
00053 delete fShowerPos;
00054 delete fShowerShape;
00055 }
00056
00057 void EmcRecSplitWeighted::Split(RecEmcCluster& aCluster,
00058 const RecEmcIDVector& aMaxVec,
00059 RecEmcShowerMap& aShowerMap)
00060 {
00061
00062
00063
00064
00065
00066 EmcRecParameter& Para=EmcRecParameter::GetInstance();
00067
00068
00069 double eCluster=aCluster.getEnergy();
00070 double cut=0;
00071 if(eCluster>Para.SmCut(3)) {
00072 cut=Para.SmCut(0)*exp(-(eCluster)/Para.SmCut(1))+Para.SmCut(2);
00073 } else {
00074 cut=Para.SmCut(0)*exp(-(Para.SmCut(3))/Para.SmCut(1))+Para.SmCut(2);
00075 }
00076 cut/=10;
00077
00078 RecEmcShower aShower;
00079 RecEmcHit aHit;
00080 RecEmcFraction aFraction;
00081 RecEmcHitMap::const_iterator ciHitMap;
00082 if(aMaxVec.size()==0) {
00083
00084 }
00085
00086
00087 else if(aMaxVec.size()==1||aMaxVec.size()>20||aCluster.getSecondMoment()<cut) {
00088 aShower.Clear();
00089 aShower.setStatus(1);
00090
00091 aShower.ShowerId(aMaxVec[0]);
00092 aShower.ClusterId(aCluster.getClusterId());
00093 aCluster.InsertShowerId(aMaxVec[0]);
00094
00095 double time=aCluster.Find(aMaxVec[0])->second.getTime();
00096 if(time>=Para.TimeMin()&&time<=Para.TimeMax()) {
00097 for(ciHitMap=aCluster.Begin();
00098 ciHitMap!=aCluster.End();
00099 ++ciHitMap) {
00100 aHit=ciHitMap->second;
00101
00102
00103 aFraction=aHit;
00104 aFraction.Fraction(1);
00105 aShower.Insert(aFraction);
00106
00107 }
00108
00109 aShower.setTime(time);
00110 fShowerE->Energy(aShower);
00111 fShowerPos->Position(aShower);
00112 fShowerShape->CalculateMoment(aShower);
00113 aShower.ThetaGap(9);
00114 aShower.PhiGap(9);
00115
00116 aShowerMap[aMaxVec[0]]=aShower;
00117 }
00118 }
00119
00120
00121
00122 else if(aMaxVec.size()>1&&aMaxVec.size()<=20&&aCluster.getSecondMoment()>cut) {
00123
00124
00125
00126
00127 RecEmcCluster tmpCluster=aCluster;
00128 RecEmcHitMap::const_iterator ci_HitMap;
00129 ci_RecEmcIDVector ciMax, ciMax1;;
00130 RecEmcShower aShower;
00131 RecEmcShowerMap tmpShowerMap;
00132 RecEmcShowerMap::iterator i_ShowerMap,i2_ShowerMap;
00133
00134 RecEmcFraction aFrac;
00135 unsigned int iterations=0;
00136 double centroidShift;
00137 map<RecEmcID,HepPoint3D,less<RecEmcID> > showerCentroid;
00138 map<RecEmcID,HepPoint3D,less<RecEmcID> >::const_iterator ci_showerCentroid;
00139
00140 IEmcRecGeoSvc* iGeoSvc;
00141 ISvcLocator* svcLocator = Gaudi::svcLocator();
00142 StatusCode sc = svcLocator->service("EmcRecGeoSvc",iGeoSvc);
00143 if(sc!=StatusCode::SUCCESS) {
00144
00145 }
00146
00147
00148 for(ciMax=aMaxVec.begin();
00149 ciMax!=aMaxVec.end();
00150 ++ciMax) {
00151
00152 showerCentroid[*ciMax]=iGeoSvc->GetCFrontCenter(*ciMax);
00153 }
00154 do {
00155 centroidShift=0.;
00156 tmpShowerMap.clear();
00157 for(ciMax=aMaxVec.begin();
00158 ciMax!=aMaxVec.end();
00159 ++ciMax) {
00160 double weightTotal=0.,weight=0.;
00161
00162 aShower.Clear();
00163 aShower.ShowerId(*ciMax);
00164
00165
00166 for(ci_HitMap=tmpCluster.Begin();
00167 ci_HitMap!=tmpCluster.End();
00168 ++ci_HitMap) {
00169
00170 aFrac=ci_HitMap->second;
00171 double aDistance=0;
00172 RecEmcEnergy aESeed=0;
00173
00174 bool isSeed=false;
00175 for(ciMax1=aMaxVec.begin();
00176 ciMax1!=aMaxVec.end();
00177 ++ciMax1) {
00178 HepPoint3D seedPoint(showerCentroid[*ciMax1]);
00179 HepPoint3D hitPoint(iGeoSvc->GetCFrontCenter(aFrac.getCellId()));
00180
00181 RecEmcEnergy theESeed=tmpCluster.Find(*ciMax1)->second.getEnergy();
00182 double theDistance;
00183 if(*ciMax1==aFrac.getCellId()) {
00184 isSeed=true;
00185 theDistance=0.;
00186 } else {
00187 theDistance=(hitPoint-seedPoint).mag();
00188 }
00189
00190 if(*ciMax1==*ciMax) {
00191 aDistance=theDistance;
00192 aESeed=theESeed;
00193 }
00194
00195 weightTotal+=theESeed*exp(-Para.LateralProfile()*theDistance/Para.MoliereRadius());
00196 }
00197
00198 weight=aESeed*exp(-Para.LateralProfile()*aDistance/Para.MoliereRadius())/weightTotal;
00199 aFrac.Fraction(weight);
00200
00201
00202
00203 if(aFrac.getEnergy()*aFrac.getFraction()>Para.ElectronicsNoiseLevel()) {
00204 aShower.Insert(aFrac);
00205 }
00206
00207 weightTotal=0;
00208 }
00209
00210 fShowerE->Energy(aShower);
00211 fShowerPos->Position(aShower);
00212 HepPoint3D newCentroid(aShower.position());
00213 HepPoint3D oldCentroid(showerCentroid[*ciMax]);
00214 centroidShift+=(newCentroid-oldCentroid).mag();
00215 tmpShowerMap[aShower.getShowerId()]=aShower;
00216 showerCentroid[*ciMax]=newCentroid;
00217 }
00218
00219 centroidShift/=(double)aMaxVec.size();
00220 for(ci_showerCentroid=showerCentroid.begin();
00221 ci_showerCentroid!=showerCentroid.end();
00222 ci_showerCentroid++){
00223 showerCentroid[ci_showerCentroid->first]
00224 =tmpShowerMap[ci_showerCentroid->first].position();
00225 }
00226 iterations++;
00227
00228 }
00229 while((iterations<8)&&(centroidShift>0.01));
00230
00231 unsigned int tht,phi;
00232 unsigned int tht2,phi2;
00233 unsigned int dtht,dphi;
00234 unsigned int thetagap=0,phigap=0;
00235 double distmin,dist;
00236 RecEmcID id,id2,nearest;
00237
00238 for(i_ShowerMap=tmpShowerMap.begin();
00239 i_ShowerMap!=tmpShowerMap.end();
00240 ++i_ShowerMap) {
00241 aCluster.InsertShowerId(i_ShowerMap->second.getShowerId());
00242 i_ShowerMap->second.ClusterId(aCluster.getClusterId());
00243 i_ShowerMap->second.setStatus(2);
00244 fShowerE->Energy(i_ShowerMap->second);
00245 fShowerPos->Position(i_ShowerMap->second);
00246 fShowerShape->CalculateMoment(i_ShowerMap->second);
00247
00248
00249 id=(i_ShowerMap->second).getShowerId();
00250 tht=EmcID::theta_module(id);
00251 phi=EmcID::phi_module(id);
00252 distmin=999999;
00253 for(i2_ShowerMap=tmpShowerMap.begin();
00254 i2_ShowerMap!=tmpShowerMap.end();
00255 ++i2_ShowerMap) {
00256 id2=(i2_ShowerMap->second).getShowerId();
00257 tht2=EmcID::theta_module(id2);
00258 phi2=EmcID::phi_module(id2);
00259 if(id!=id2) {
00260 dtht=tht>tht2 ? tht-tht2 : tht2-tht;
00261 dphi=phi>phi2 ? phi-phi2 : phi2-phi;
00262 if(dphi>(EmcID::getPHI_BARREL_MAX()+1)/2) dphi=EmcID::getPHI_BARREL_MAX()+1-dphi;
00263 dist=sqrt(double(dtht*dtht+dphi*dphi));
00264 if(dist<distmin) {
00265 distmin=dist;
00266 nearest=id2;
00267 if(dtht>6) dtht=6;
00268 if(dphi>6) dphi=6;
00269 thetagap=dtht;
00270 phigap=dphi;
00271 }
00272 }
00273 }
00274
00275 i_ShowerMap->second.NearestSeed(nearest);
00276 i_ShowerMap->second.ThetaGap(thetagap);
00277 i_ShowerMap->second.PhiGap(phigap);
00278
00279
00280 double time=
00281 i_ShowerMap->second.Find(i_ShowerMap->second.getShowerId())->second.getTime();
00282 if(time>=Para.TimeMin() && time<=Para.TimeMax() &&
00283 (i_ShowerMap->second.getEAll()>Para.EThresholdCluster())) {
00284 i_ShowerMap->second.setTime(time);
00285 aShowerMap[i_ShowerMap->first]=i_ShowerMap->second;
00286 }
00287 }
00288 tmpShowerMap.clear();
00289 }
00290
00291 }