00001
00002
00003
00004
00005
00006 #include <iostream>
00007
00008 #include "EmcRecEventModel/RecEmcCluster.h"
00009 #include "EmcRecGeoSvc/EmcRecGeoSvc.h"
00010 #include "GaudiKernel/Bootstrap.h"
00011 #include "GaudiKernel/ISvcLocator.h"
00012
00013
00014 RecEmcCluster::RecEmcCluster()
00015 {
00016 Clear();
00017 }
00018
00019 RecEmcCluster::~RecEmcCluster()
00020 {
00021 Clear();
00022 }
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059
00060
00061
00062
00063
00064
00065
00066
00067
00068
00069
00070
00071
00072
00073
00074
00075
00076
00077
00078
00079
00080
00081
00082
00083
00084
00085
00086
00087
00088
00089
00090
00091
00092
00093
00094
00095
00096 void RecEmcCluster::Clear()
00097 {
00098 fClusterId.clear();
00099 fHitMap.clear();
00100 fSeedMap.clear();
00101 fShowerIdVec.clear();
00102 }
00103
00104
00105
00106
00107
00108
00109
00110 void RecEmcCluster::ClusterId(const RecEmcID id)
00111 {
00112 fClusterId=id;
00113 }
00114
00115 RecEmcHitMap::const_iterator RecEmcCluster::Begin() const
00116 {
00117 return fHitMap.begin();
00118 }
00119
00120 RecEmcHitMap::const_iterator RecEmcCluster::End() const
00121 {
00122 return fHitMap.end();
00123 }
00124
00125 RecEmcHitMap::const_iterator RecEmcCluster::Find(const RecEmcID& CellId) const
00126 {
00127
00128 return fHitMap.find(CellId);
00129 }
00130
00131
00132 void RecEmcCluster::Insert(const RecEmcHit& aHit)
00133 {
00134 fHitMap[aHit.getCellId()]=aHit;
00135
00136 return;
00137 }
00138
00139 void RecEmcCluster::Erase(const RecEmcHit& aHit)
00140 {
00141 RecEmcHitMap::const_iterator pHitMap;
00142 pHitMap=fHitMap.find(aHit.getCellId());
00143
00144
00145 if(fHitMap.empty())
00146 {
00147 return;
00148 }
00149
00150
00151 if(pHitMap==End())
00152 {
00153 return;
00154 }
00155
00156
00157 if(pHitMap!=End())
00158 {
00159 fHitMap.erase(pHitMap->first);
00160
00161 if(fHitMap.empty())
00162 {
00163 Clear();
00164 return;
00165 }
00166
00167 else
00168 {
00169 fClusterId=fHitMap.begin()->first;
00170 return;
00171 }
00172 }
00173 }
00174
00175
00176 RecEmcHitMap::const_iterator RecEmcCluster::BeginSeed() const
00177 {
00178 return fSeedMap.begin();
00179 }
00180
00181 RecEmcHitMap::const_iterator RecEmcCluster::EndSeed() const
00182 {
00183 return fSeedMap.end();
00184 }
00185
00186 RecEmcHitMap::const_iterator RecEmcCluster::FindSeed(const RecEmcID& CellId) const
00187 {
00188
00189 return fSeedMap.find(CellId);
00190 }
00191
00192
00193 void RecEmcCluster::InsertSeed(const RecEmcHit& aSeed)
00194 {
00195 fSeedMap[aSeed.getCellId()]=aSeed;
00196 return;
00197 }
00198
00199 int RecEmcCluster::getSeedSize() const
00200 {
00201 return fSeedMap.size();
00202 }
00203
00204 int RecEmcCluster::getShowerSize() const
00205 {
00206 return fShowerIdVec.size();
00207 }
00208
00209
00210 void RecEmcCluster::InsertShowerId(const RecEmcID id)
00211 {
00212 fShowerIdVec.push_back(id);
00213 }
00214
00215
00216 RecEmcEnergy RecEmcCluster::getEnergy() const
00217 {
00218 RecEmcHitMap::const_iterator pHitMap;
00219 double etot=0;
00220
00221 for(pHitMap=fHitMap.begin();
00222 pHitMap!=fHitMap.end();
00223 pHitMap++) {
00224 etot+=pHitMap->second.getEnergy();
00225 }
00226 return etot;
00227 }
00228
00229
00230 HepPoint3D RecEmcCluster::getPosition() const
00231 {
00232 IEmcRecGeoSvc* iGeoSvc;
00233 ISvcLocator* svcLocator = Gaudi::svcLocator();
00234 StatusCode sc = svcLocator->service("EmcRecGeoSvc",iGeoSvc);
00235 if(sc!=StatusCode::SUCCESS) {
00236 cout<<"Error: Can't get EmcRecGeoSvc"<<endl;
00237 }
00238
00239 RecEmcHitMap::const_iterator pHitMap;
00240 HepPoint3D pos(0,0,0);
00241 HepPoint3D possum(0,0,0);
00242 double etot=0;
00243
00244 for(pHitMap=fHitMap.begin();
00245 pHitMap!=fHitMap.end();
00246 pHitMap++) {
00247 etot+=pHitMap->second.getEnergy();
00248 pos=iGeoSvc->GetCFrontCenter(pHitMap->second.getCellId());
00249 possum+=pos*pHitMap->second.getEnergy();
00250 }
00251
00252 if(etot>0) {
00253 possum/=etot;
00254 }
00255 return possum;
00256 }
00257
00258
00259 double RecEmcCluster::getSecondMoment() const
00260 {
00261 IEmcRecGeoSvc* iGeoSvc;
00262 ISvcLocator* svcLocator = Gaudi::svcLocator();
00263 StatusCode sc = svcLocator->service("EmcRecGeoSvc",iGeoSvc);
00264 if(sc!=StatusCode::SUCCESS) {
00265 cout<<"Error: Can't get EmcRecGeoSvc"<<endl;
00266 }
00267
00268 double etot=0;
00269 double sum=0;
00270 HepPoint3D center(getPosition());
00271 RecEmcHitMap::const_iterator pHitMap;
00272
00273 for(pHitMap=fHitMap.begin();
00274 pHitMap!=fHitMap.end();
00275 pHitMap++) {
00276 HepPoint3D pos(pHitMap->second.getFrontCenter());
00277 etot+=pHitMap->second.getEnergy();
00278 sum+=pHitMap->second.getEnergy()*pos.distance2(center);
00279 }
00280
00281 if(etot>0) {
00282 sum/=etot;
00283 }
00284 return sum;
00285 }
00286
00287
00288 void RecEmcCluster::Dump() const
00289 {
00290 RecEmcHitMap::const_iterator pHitMap;
00291
00292 cout<<"EMC Cluster: ";
00293
00294 cout<<"Cluster Id= ";
00295 cout<<fClusterId<<endl;
00296
00297 for(pHitMap=fHitMap.begin();
00298 pHitMap!=fHitMap.end();
00299 pHitMap++)
00300 {
00301 pHitMap->second.Dump();
00302 }
00303 }
00304
00305 ostream& operator<<(ostream & os, const RecEmcCluster& aCluster)
00306 {
00307 RecEmcHitMap::const_iterator pHitMap;
00308
00309 cout<<"EMC Cluster: ";
00310
00311 cout<<"Cluster Id= ";
00312 cout<<aCluster.getClusterId()<<endl;
00313
00314 for(pHitMap=aCluster.Begin();
00315 pHitMap!=aCluster.End();
00316 pHitMap++)
00317 {
00318 os<< (pHitMap->second);
00319 }
00320
00321 if(aCluster.getSeedSize()>0) {
00322 cout<<"Contains "<<aCluster.getSeedSize()<<" Seeds:"<<endl;
00323 RecEmcHitMap::const_iterator pSeedMap;
00324 for(pSeedMap=aCluster.BeginSeed();
00325 pSeedMap!=aCluster.EndSeed();
00326 pSeedMap++)
00327 {
00328 os<< (pSeedMap->second);
00329 }
00330 }
00331
00332 if(aCluster.getShowerSize()>0) {
00333 vector<RecEmcID> aShowerIdVec=aCluster.getShowerIdVec();
00334 vector<RecEmcID>::iterator iShowerId;
00335 os<<"Contains "<<aCluster.getShowerSize()<<" Showers:"<<endl;
00336 for(iShowerId=aShowerIdVec.begin();
00337 iShowerId!=aShowerIdVec.end();
00338 iShowerId++) {
00339 os << *iShowerId <<endl;
00340 }
00341 }
00342
00343 return os;
00344 }