EmcRecSplitWeighted Class Reference

#include <EmcRecSplitWeighted.h>

Inheritance diagram for EmcRecSplitWeighted:

EmcRecSplitAbs List of all members.

Public Member Functions

 EmcRecSplitWeighted ()
 ~EmcRecSplitWeighted ()
virtual void Split (RecEmcCluster &aCluster, const RecEmcIDVector &aMaxVec, RecEmcShowerMap &aShowerMap)

Public Attributes

EmcRecShowerEnergy * fShowerE
EmcRecShowerPosAbs * fShowerPos
EmcRecShowerShape * fShowerShape

Detailed Description

Definition at line 17 of file EmcRecSplitWeighted.h.

Constructor & Destructor Documentation

EmcRecSplitWeighted::EmcRecSplitWeighted (  ) 

Definition at line 23 of file EmcRecSplitWeighted.cxx.

References fShowerE, fShowerPos, fShowerShape, EmcRecParameter::GetInstance(), and EmcRecParameter::PositionMode1().

00024 {
00025   fShowerE=new EmcRecShowerEnergy;
00026   fShowerShape=new EmcRecShowerShape;
00027   
00028   EmcRecParameter& Para=EmcRecParameter::GetInstance();
00029 //  if(Para.PositionMode1()=="lin") {
00030 //    fShowerPos=new EmcRecShowerPosLin;
00031 //  } else {
00032 //    fShowerPos=new EmcRecShowerPosLog;
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 }

EmcRecSplitWeighted::~EmcRecSplitWeighted (  ) 

Definition at line 50 of file EmcRecSplitWeighted.cxx.

References fShowerE, fShowerPos, and fShowerShape.

00051 {
00052    delete fShowerE;
00053    delete fShowerPos;
00054    delete fShowerShape;
00055 }

Member Function Documentation

void EmcRecSplitWeighted::Split ( RecEmcCluster &  aCluster,
const RecEmcIDVector &  aMaxVec,
RecEmcShowerMap &  aShowerMap 
) [virtual]

Implements EmcRecSplitAbs.

Definition at line 57 of file EmcRecSplitWeighted.cxx.

References RecEmcCluster::Begin(), EmcRecShowerShape::CalculateMoment(), RecEmcShower::Clear(), RecEmcShower::ClusterId(), cut, check_raw_filter::dist, EmcRecParameter::ElectronicsNoiseLevel(), RecEmcCluster::End(), EmcRecShowerEnergy::Energy(), EmcRecParameter::EThresholdCluster(), exp(), RecEmcCluster::Find(), RecEmcFraction::Fraction(), fShowerE, fShowerPos, fShowerShape, RecEmcHit::getCellId(), IEmcRecGeoSvc::GetCFrontCenter(), RecEmcCluster::getClusterId(), RecEmcHit::getEnergy(), RecEmcCluster::getEnergy(), RecEmcFraction::getFraction(), EmcRecParameter::GetInstance(), EmcID::getPHI_BARREL_MAX(), RecEmcCluster::getSecondMoment(), RecEmcShower::getShowerId(), RecEmcShower::Insert(), RecEmcCluster::InsertShowerId(), EmcRecParameter::LateralProfile(), EmcRecParameter::MoliereRadius(), phi2, EmcID::phi_module(), RecEmcShower::PhiGap(), DstEmcShower::position(), EmcRecShowerPosAbs::Position(), EvtCyclic3::second(), DstEmcShower::setStatus(), DstEmcShower::setTime(), RecEmcShower::ShowerId(), EmcRecParameter::SmCut(), EmcID::theta_module(), RecEmcShower::ThetaGap(), EmcRecParameter::TimeMax(), EmcRecParameter::TimeMin(), and weight.

00060 {
00061    //A exponential function will be used to describe the shape of a shower.
00062    //Energy fall off with distance from centre of shower is exp(-a*dist/r)
00063    //dist: distance from center
00064    //r: moliere radius
00065    //a: a parameter, needed to be optimised
00066   EmcRecParameter& Para=EmcRecParameter::GetInstance();
00067 
00068   //Calculate the cut of second moment
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);  //constant
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     // the cluster is abandoned
00084   }
00085    
00086   // only one seed or second moment belows cut
00087   else if(aMaxVec.size()==1||aMaxVec.size()>20||aCluster.getSecondMoment()<cut) {
00088     aShower.Clear();
00089     aShower.setStatus(1);
00090     // ID
00091     aShower.ShowerId(aMaxVec[0]);
00092     aShower.ClusterId(aCluster.getClusterId());
00093     aCluster.InsertShowerId(aMaxVec[0]);
00094     // each fraction
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         //if((aHit.time()>=Para.TimeMin()) &&
00102            //(aHit.time()<=Para.TimeMax())){
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       // push it into map
00116       aShowerMap[aMaxVec[0]]=aShower;
00117     }
00118   }
00119 
00120 
00121   // two or more seeds and second moment beyonds cut
00122   else if(aMaxVec.size()>1&&aMaxVec.size()<=20&&aCluster.getSecondMoment()>cut) {
00123     //cout<<"In EmcRecSplitWeighted: ";
00124     //cout<<aMaxVec.size();
00125     //cout<<" seeds in a cluster"<<endl;      
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;  //iteration time
00136     double centroidShift;   //center shift between 2 iterations
00137     map<RecEmcID,HepPoint3D,less<RecEmcID> > showerCentroid;  //shower center
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       //cout<<"Error: Can't get EmcRecGeoSvc"<<endl;
00145     }
00146 
00147     //Preparation
00148     for(ciMax=aMaxVec.begin();
00149         ciMax!=aMaxVec.end();
00150         ++ciMax) {
00151       //initialize: seed's front center
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         // figure out the weight for each cell for each seed
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           //if((aFrac.time()>=Para.TimeMin() &&
00201               //aFrac.time()<=Para.TimeMax()) ||
00202               //isSeed) {
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       //cout<<"iterations: "<<iterations<<"\tcentroidShift: "<<centroidShift<<endl;
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     // shower attributes
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             // min theta, phi gap
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             //cout<<"+++++"<<id<<" "<<thetagap<<" "<<phigap<<endl;
00275             i_ShowerMap->second.NearestSeed(nearest);
00276             i_ShowerMap->second.ThetaGap(thetagap);
00277             i_ShowerMap->second.PhiGap(phigap);
00278          
00279             // save it
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 }

Member Data Documentation

EmcRecShowerEnergy* EmcRecSplitWeighted::fShowerE

Definition at line 29 of file EmcRecSplitWeighted.h.

Referenced by EmcRecSplitWeighted(), Split(), and ~EmcRecSplitWeighted().

EmcRecShowerPosAbs* EmcRecSplitWeighted::fShowerPos

Definition at line 30 of file EmcRecSplitWeighted.h.

Referenced by EmcRecSplitWeighted(), Split(), and ~EmcRecSplitWeighted().

EmcRecShowerShape* EmcRecSplitWeighted::fShowerShape

Definition at line 31 of file EmcRecSplitWeighted.h.

Referenced by EmcRecSplitWeighted(), Split(), and ~EmcRecSplitWeighted().

Generated on Tue Nov 29 23:18:46 2016 for BOSS_7.0.2 by  doxygen 1.4.7