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