ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/yangyong/ECALG4SIM/CaloSD.cc
Revision: 1.2
Committed: Thu Apr 4 08:48:25 2013 UTC (12 years, 1 month ago) by yangyong
Content type: text/plain
Branch: MAIN
CVS Tags: HEAD
Changes since 1.1: +6 -0 lines
Log Message: 
*** empty log message ***

File Contents

# User Rev Content
1 yangyong 1.1 ///////////////////////////////////////////////////////////////////////////////
2     // File: CaloSD.cc
3     // Description: Sensitive Detector class for calorimeters
4     ///////////////////////////////////////////////////////////////////////////////
5    
6     #include "SimG4CMS/Calo/interface/CaloSD.h"
7     #include "SimDataFormats/SimHitMaker/interface/CaloSlaveSD.h"
8     #include "SimG4Core/Notification/interface/TrackInformation.h"
9     #include "SimG4Core/Application/interface/EventAction.h"
10    
11     #include "G4EventManager.hh"
12     #include "G4SDManager.hh"
13     #include "G4Step.hh"
14     #include "G4Track.hh"
15     #include "G4VProcess.hh"
16     #include "G4GFlashSpot.hh"
17     #include "G4ParticleTable.hh"
18    
19    
20     #include "DataFormats/DetId/interface/DetId.h"
21     #include "DataFormats/EcalDetId/interface/EBDetId.h"
22     #include "DataFormats/EcalDetId/interface/EEDetId.h"
23     #include "Geometry/CaloGeometry/interface/CaloGenericDetId.h"
24     #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
25     #include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
26    
27 yangyong 1.2 #include "TROOT.h"
28     #include "TFile.h"
29     #include "TTree.h"
30     #include "TBranch.h"
31     #include "TLorentzVector.h"
32     #include "TClonesArray.h"
33 yangyong 1.1
34    
35     TFile *g4f;
36     TTree *g4t;
37    
38     bool saveTree;
39    
40     G4String csdname;
41    
42     using namespace std;
43    
44     int ncounter1;
45     int ncounter2;
46     std::map<int,std::map<int, std::vector<float> > > estepg4EB;
47     std::map<int,std::map<int, std::vector<float> > > tstepg4EB;
48     std::map<int,std::map<int, std::vector<int> > > idstepg4EB;
49     std::map<int,std::map<int, std::vector<int> > > pidstepg4EB;
50     std::map<int,std::map<int, std::vector<int> > > parentidstepg4EB;
51    
52     std::map<int,std::map<int, std::vector<int> > > enterstepg4EB;
53     std::map<int,std::map<int, std::vector<int> > > leavestepg4EB;
54    
55     std::map<int,std::map<int, std::vector<float> > > postxstepg4EB;
56     std::map<int,std::map<int, std::vector<float> > > postystepg4EB;
57     std::map<int,std::map<int, std::vector<float> > > postzstepg4EB;
58    
59     std::map<int,std::map<int, std::vector<float> > > pretstepg4EB;
60     std::map<int,std::map<int, std::vector<float> > > prexstepg4EB;
61     std::map<int,std::map<int, std::vector<float> > > preystepg4EB;
62     std::map<int,std::map<int, std::vector<float> > > prezstepg4EB;
63    
64     std::map<int,std::map<int, std::vector<float> > > preestepg4EB;
65    
66     //std::map<int,std::map<int, std::vector<float> > > prepxstepg4EB;
67     //std::map<int,std::map<int, std::vector<float> > > prepystepg4EB;
68     //std::map<int,std::map<int, std::vector<float> > > prepzstepg4EB;
69     std::map<int,std::map<int, std::vector<int> > > prechazstepg4EB;
70     std::map<int,std::map<int, std::vector<float> > > premzstepg4EB;
71    
72    
73     std::map<int,std::map<int, std::vector<float> > > estepg4EEm;
74     std::map<int,std::map<int, std::vector<float> > > tstepg4EEm;
75     std::map<int,std::map<int, std::vector<float> > > estepg4EEp;
76     std::map<int,std::map<int, std::vector<float> > > tstepg4EEp;
77     std::map<int,std::map<int, std::vector<int> > > idstepg4EEm;
78     std::map<int,std::map<int, std::vector<int> > > pidstepg4EEm;
79     std::map<int,std::map<int, std::vector<int> > > parentidstepg4EEm;
80     std::map<int,std::map<int, std::vector<int> > > enterstepg4EEm;
81     std::map<int,std::map<int, std::vector<int> > > leavestepg4EEm;
82    
83    
84     std::map<int,std::map<int, std::vector<float> > > postxstepg4EEm;
85     std::map<int,std::map<int, std::vector<float> > > postystepg4EEm;
86     std::map<int,std::map<int, std::vector<float> > > postzstepg4EEm;
87    
88     std::map<int,std::map<int, std::vector<float> > > pretstepg4EEm;
89    
90     std::map<int,std::map<int, std::vector<float> > > prexstepg4EEm;
91     std::map<int,std::map<int, std::vector<float> > > preystepg4EEm;
92     std::map<int,std::map<int, std::vector<float> > > prezstepg4EEm;
93    
94    
95     std::map<int,std::map<int, std::vector<float> > > preestepg4EEm;
96    
97     std::map<int,std::map<int, std::vector<int> > > prechazstepg4EEm;
98     std::map<int,std::map<int, std::vector<float> > > premzstepg4EEm;
99    
100    
101    
102     std::map<int,std::map<int, std::vector<int> > > idstepg4EEp;
103     std::map<int,std::map<int, std::vector<int> > > pidstepg4EEp;
104     std::map<int,std::map<int, std::vector<int> > > parentidstepg4EEp;
105    
106     std::map<int,std::map<int, std::vector<float> > > postxstepg4EEp;
107     std::map<int,std::map<int, std::vector<float> > > postystepg4EEp;
108     std::map<int,std::map<int, std::vector<float> > > postzstepg4EEp;
109     std::map<int,std::map<int, std::vector<int> > > enterstepg4EEp;
110     std::map<int,std::map<int, std::vector<int> > > leavestepg4EEp;
111     std::map<int,std::map<int, std::vector<float> > > pretstepg4EEp;
112     std::map<int,std::map<int, std::vector<float> > > prexstepg4EEp;
113     std::map<int,std::map<int, std::vector<float> > > preystepg4EEp;
114     std::map<int,std::map<int, std::vector<float> > > prezstepg4EEp;
115    
116     std::map<int,std::map<int, std::vector<float> > > preestepg4EEp;
117    
118     std::map<int,std::map<int, std::vector<int> > > prechazstepg4EEp;
119     std::map<int,std::map<int, std::vector<float> > > premzstepg4EEp;
120    
121    
122    
123    
124     int ng4EB;
125     int ietag4EB[61200];
126     int iphig4EB[61200];
127     float esumg4EB[61200];
128     float tming4EB[61200];
129     float xtming4EB[61200];
130     float ytming4EB[61200];
131     float ztming4EB[61200];
132    
133     std::vector<std::vector<float> > *eg4EB;
134     std::vector<std::vector<float> > *tg4EB;
135     std::vector<std::vector<int> > *idg4EB;
136     std::vector<std::vector<int> > *pidg4EB;
137     std::vector<std::vector<int> > *parentidg4EB;
138     std::vector<std::vector<int> > *enterg4EB;
139     std::vector<std::vector<int> > *leaveg4EB;
140     std::vector<std::vector<float> > *preeg4EB;
141    
142    
143    
144     std::vector<std::vector<float> > *postxg4EB;
145     std::vector<std::vector<float> > *postyg4EB;
146     std::vector<std::vector<float> > *postzg4EB;
147    
148     std::vector<std::vector<float> > *prexg4EB;
149     std::vector<std::vector<float> > *preyg4EB;
150     std::vector<std::vector<float> > *prezg4EB;
151     std::vector<std::vector<float> > *pretg4EB;
152     std::vector<std::vector<int> > *prechag4EB;
153     std::vector<std::vector<float> > *premg4EB;
154    
155    
156     int ng4EE;
157     int ixg4EE[14648];
158     int iyg4EE[14648];
159     int izg4EE[14648];
160     float esumg4EE[14648];
161     float tming4EE[14648];
162     float xtming4EE[14648];
163     float ytming4EE[14648];
164     float ztming4EE[14648];
165    
166    
167     std::vector<std::vector<float> > *eg4EE;
168     std::vector<std::vector<float> > *tg4EE;
169     std::vector<std::vector<int> > *idg4EE;
170     std::vector<std::vector<int> > *pidg4EE;
171     std::vector<std::vector<int> > *parentidg4EE;
172    
173     std::vector<std::vector<int> > *enterg4EE;
174     std::vector<std::vector<int> > *leaveg4EE;
175    
176     std::vector<std::vector<float> > *preeg4EE;
177    
178    
179     std::vector<std::vector<float> > *postxg4EE;
180     std::vector<std::vector<float> > *postyg4EE;
181     std::vector<std::vector<float> > *postzg4EE;
182    
183     std::vector<std::vector<float> > *prexg4EE;
184     std::vector<std::vector<float> > *preyg4EE;
185     std::vector<std::vector<float> > *prezg4EE;
186     std::vector<std::vector<float> > *pretg4EE;
187     std::vector<std::vector<int> > *prechag4EE;
188     std::vector<std::vector<float> > *premg4EE;
189    
190    
191     int debugCaloSD = 0;
192    
193     int treeConfig = 1;
194     int saveOnlyEnterCur = 1;
195     //int saveOnlyEnterCur = 0;
196    
197     //#define DebugLog
198    
199     CaloSD::CaloSD(G4String name, const DDCompactView & cpv,
200     SensitiveDetectorCatalog & clg,
201     edm::ParameterSet const & p, const SimTrackManager* manager,
202     int tSlice, bool ignoreTkID) :
203     SensitiveCaloDetector(name, cpv, clg, p),
204     G4VGFlashSensitiveDetector(), theTrack(0), preStepPoint(0), eminHit(0),
205     eminHitD(0), m_trackManager(manager), currentHit(0), runInit(false),
206     timeSlice(tSlice), ignoreTrackID(ignoreTkID), hcID(-1), theHC(0),
207     meanResponse(0) {
208     //Add Hcal Sentitive Detector Names
209    
210     collectionName.insert(name);
211    
212     csdname = name;
213    
214     if(debugCaloSD) cout<<"CaloSD::CaloSD " << name <<" tslice " << tSlice <<endl;
215    
216     if( "EcalHitsEB" == name){
217     if(debugCaloSD) cout<<"creat newtree" <<endl;
218    
219     ncounter1 = 0;
220     ncounter2 = 0;
221    
222    
223     g4f = new TFile("g4simhitsEcal.root","RECREATE");
224     g4t = new TTree("G4SIM","g4sim");
225    
226     eg4EB = new std::vector<std::vector<float> >; eg4EB->clear();
227     tg4EB = new std::vector<std::vector<float> >; tg4EB->clear();
228     idg4EB = new std::vector<std::vector<int> >; idg4EB->clear();
229     pidg4EB = new std::vector<std::vector<int> >; pidg4EB->clear();
230     parentidg4EB = new std::vector<std::vector<int> >; parentidg4EB->clear();
231     enterg4EB = new std::vector<std::vector<int> >; enterg4EB->clear();
232     leaveg4EB = new std::vector<std::vector<int> >; leaveg4EB->clear();
233    
234     preeg4EB = new std::vector<std::vector<float> >; preeg4EB->clear();
235     prechag4EB = new std::vector<std::vector<int> >; prechag4EB->clear();
236     premg4EB = new std::vector<std::vector<float> >; premg4EB->clear();
237    
238    
239    
240     postxg4EB = new std::vector<std::vector<float> >; postxg4EB->clear();
241     postyg4EB = new std::vector<std::vector<float> >; postyg4EB->clear();
242     postzg4EB = new std::vector<std::vector<float> >; postzg4EB->clear();
243     prexg4EB = new std::vector<std::vector<float> >; prexg4EB->clear();
244     preyg4EB = new std::vector<std::vector<float> >; preyg4EB->clear();
245     prezg4EB = new std::vector<std::vector<float> >; prezg4EB->clear();
246     pretg4EB = new std::vector<std::vector<float> >; pretg4EB->clear();
247     prexg4EE = new std::vector<std::vector<float> >; prexg4EE->clear();
248     preyg4EE = new std::vector<std::vector<float> >; preyg4EE->clear();
249     prezg4EE = new std::vector<std::vector<float> >; prezg4EE->clear();
250     pretg4EE = new std::vector<std::vector<float> >; pretg4EE->clear();
251    
252     eg4EE = new std::vector<std::vector<float> >; eg4EE->clear();
253     tg4EE = new std::vector<std::vector<float> >; tg4EE->clear();
254     idg4EE = new std::vector<std::vector<int> >; idg4EE->clear();
255     pidg4EE = new std::vector<std::vector<int> >; pidg4EE->clear();
256     parentidg4EE = new std::vector<std::vector<int> >; parentidg4EE->clear();
257     enterg4EE = new std::vector<std::vector<int> >; enterg4EE->clear();
258     leaveg4EE = new std::vector<std::vector<int> >; leaveg4EE->clear();
259    
260    
261     postxg4EE = new std::vector<std::vector<float> >; postxg4EE->clear();
262     postyg4EE = new std::vector<std::vector<float> >; postyg4EE->clear();
263     postzg4EE = new std::vector<std::vector<float> >; postzg4EE->clear();
264    
265     preeg4EE = new std::vector<std::vector<float> >; preeg4EE->clear();
266    
267    
268    
269     g4t->Branch("ng4EB", &ng4EB, "ng4EB/I");
270     g4t->Branch("ietag4EB", ietag4EB, "ietag4EB[ng4EB]/I");
271     g4t->Branch("iphig4EB", iphig4EB, "iphig4EB[ng4EB]/I");
272     g4t->Branch("esumg4EB", esumg4EB, "esumg4EB[ng4EB]/F");
273     g4t->Branch("tming4EB", tming4EB, "tming4EB[ng4EB]/F");
274     g4t->Branch("xtming4EB", xtming4EB, "xtming4EB[ng4EB]/F");
275     g4t->Branch("ytming4EB", ytming4EB, "ytming4EB[ng4EB]/F");
276     g4t->Branch("ztming4EB", ztming4EB, "ztming4EB[ng4EB]/F");
277    
278    
279     g4t->Branch("ng4EE", &ng4EE, "ng4EE/I");
280     g4t->Branch("ixg4EE", ixg4EE, "ixg4EE[ng4EE]/I");
281     g4t->Branch("iyg4EE", iyg4EE, "iyg4EE[ng4EE]/I");
282     g4t->Branch("izg4EE", izg4EE, "izg4EE[ng4EE]/I");
283     g4t->Branch("esumg4EE", esumg4EE, "esumg4EE[ng4EE]/F");
284     g4t->Branch("tming4EE", tming4EE, "tming4EE[ng4EE]/F");
285     g4t->Branch("xtming4EE", xtming4EE, "xtming4EE[ng4EE]/F");
286     g4t->Branch("ytming4EE", ytming4EE, "ytming4EE[ng4EE]/F");
287     g4t->Branch("ztming4EE", ztming4EE, "ztming4EE[ng4EE]/F");
288    
289    
290     if(treeConfig==1){
291     g4t->Branch("eg4EB", "std::vector<std::vector<float> >", &eg4EB);
292     g4t->Branch("tg4EB", "std::vector<std::vector<float> >", &tg4EB);
293     g4t->Branch("idg4EB", "std::vector<std::vector<int> >", &idg4EB);
294     g4t->Branch("pidg4EB", "std::vector<std::vector<int> >", &pidg4EB);
295     g4t->Branch("parentidg4EB", "std::vector<std::vector<int> >", &parentidg4EB);
296     g4t->Branch("enterg4EB", "std::vector<std::vector<int> >", &enterg4EB);
297     g4t->Branch("leaveg4EB", "std::vector<std::vector<int> >", &leaveg4EB);
298    
299     g4t->Branch("postxg4EB", "std::vector<std::vector<float> >", &postxg4EB);
300     g4t->Branch("postyg4EB", "std::vector<std::vector<float> >", &postyg4EB);
301     g4t->Branch("postzg4EB", "std::vector<std::vector<float> >", &postzg4EB);
302    
303     g4t->Branch("prexg4EB", "std::vector<std::vector<float> >", &prexg4EB);
304     g4t->Branch("preyg4EB", "std::vector<std::vector<float> >", &preyg4EB);
305     g4t->Branch("prezg4EB", "std::vector<std::vector<float> >", &prezg4EB);
306     g4t->Branch("pretg4EB", "std::vector<std::vector<float> >", &pretg4EB);
307    
308     g4t->Branch("preeg4EB", "std::vector<std::vector<float> >", &preeg4EB);
309    
310    
311    
312     g4t->Branch("eg4EE", "std::vector<std::vector<float> >", &eg4EE);
313     g4t->Branch("tg4EE", "std::vector<std::vector<float> >", &tg4EE);
314     g4t->Branch("idg4EE", "std::vector<std::vector<int> >", &idg4EE);
315     g4t->Branch("pidg4EE", "std::vector<std::vector<int> >", &pidg4EE);
316     g4t->Branch("parentidg4EE", "std::vector<std::vector<int> >", &parentidg4EE);
317    
318     g4t->Branch("enterg4EE", "std::vector<std::vector<int> >", &enterg4EE);
319     g4t->Branch("leaveg4EE", "std::vector<std::vector<int> >", &leaveg4EE);
320    
321    
322     g4t->Branch("postxg4EE", "std::vector<std::vector<float> >", &postxg4EE);
323     g4t->Branch("postyg4EE", "std::vector<std::vector<float> >", &postyg4EE);
324     g4t->Branch("postzg4EE", "std::vector<std::vector<float> >", &postzg4EE);
325     g4t->Branch("prexg4EE", "std::vector<std::vector<float> >", &prexg4EE);
326     g4t->Branch("preyg4EE", "std::vector<std::vector<float> >", &preyg4EE);
327     g4t->Branch("prezg4EE", "std::vector<std::vector<float> >", &prezg4EE);
328     g4t->Branch("pretg4EE", "std::vector<std::vector<float> >", &pretg4EE);
329    
330     g4t->Branch("preeg4EE", "std::vector<std::vector<float> >", &preeg4EE);
331    
332    
333    
334     }
335    
336     if(debugCaloSD) cout<<"tree branch set " <<endl;
337    
338     saveTree = true;
339     }
340    
341    
342    
343    
344    
345    
346     //Parameters
347     edm::ParameterSet m_CaloSD = p.getParameter<edm::ParameterSet>("CaloSD");
348     energyCut = m_CaloSD.getParameter<double>("EminTrack")*GeV;
349     tmaxHit = m_CaloSD.getParameter<double>("TmaxHit")*ns;
350     std::vector<double> eminHits = m_CaloSD.getParameter<std::vector<double> >("EminHits");
351     std::vector<double> tmaxHits = m_CaloSD.getParameter<std::vector<double> >("TmaxHits");
352     std::vector<std::string> hcn = m_CaloSD.getParameter<std::vector<std::string> >("HCNames");
353     std::vector<int> useResMap = m_CaloSD.getParameter<std::vector<int> >("UseResponseTables");
354     std::vector<double> eminHitX = m_CaloSD.getParameter<std::vector<double> >("EminHitsDepth");
355     suppressHeavy= m_CaloSD.getParameter<bool>("SuppressHeavy");
356     kmaxIon = m_CaloSD.getParameter<double>("IonThreshold")*MeV;
357     kmaxProton = m_CaloSD.getParameter<double>("ProtonThreshold")*MeV;
358     kmaxNeutron = m_CaloSD.getParameter<double>("NeutronThreshold")*MeV;
359     checkHits = m_CaloSD.getUntrackedParameter<int>("CheckHits", 25);
360     useMap = m_CaloSD.getUntrackedParameter<bool>("UseMap", true);
361     int verbn = m_CaloSD.getUntrackedParameter<int>("Verbosity", 0);
362     corrTOFBeam = m_CaloSD.getParameter<bool>("CorrectTOFBeam");
363     double beamZ = m_CaloSD.getParameter<double>("BeamPosition")*cm;
364     correctT = beamZ/c_light/nanosecond;
365    
366     SetVerboseLevel(verbn);
367     for (unsigned int k=0; k<hcn.size(); ++k) {
368     if (name == (G4String)(hcn[k])) {
369     if (k < eminHits.size()) eminHit = eminHits[k]*MeV;
370     if (k < eminHitX.size()) eminHitD= eminHitX[k]*MeV;
371     if (k < tmaxHits.size()) tmaxHit = tmaxHits[k]*ns;
372     if (k < useResMap.size() && useResMap[k] > 0) meanResponse = new CaloMeanResponse(p);
373     break;
374     }
375     }
376     #ifdef DebugLog
377     LogDebug("CaloSim") << "***************************************************"
378     << "\n"
379     << "* *"
380     << "\n"
381     << "* Constructing a CaloSD with name " << GetName()
382     << "\n"
383     << "* *"
384     << "\n"
385     << "***************************************************";
386     #endif
387     slave = new CaloSlaveSD(name);
388     currentID = CaloHitID(timeSlice, ignoreTrackID);
389     previousID = CaloHitID(timeSlice, ignoreTrackID);
390    
391     primAncestor = 0;
392     cleanIndex = 0;
393     totalHits = 0;
394     forceSave = false;
395    
396     //
397     // Now attach the right detectors (LogicalVolumes) to me
398     //
399     std::vector<std::string> lvNames = clg.logicalNames(name);
400     this->Register();
401     for (std::vector<std::string>::iterator it=lvNames.begin(); it !=lvNames.end(); ++it) {
402     this->AssignSD(*it);
403     #ifdef DebugLog
404     LogDebug("CaloSim") << "CaloSD : Assigns SD to LV " << (*it);
405     #endif
406     }
407    
408     edm::LogInfo("CaloSim") << "CaloSD: Minimum energy of track for saving it "
409     << energyCut/GeV << " GeV" << "\n"
410     << " Use of HitID Map " << useMap << "\n"
411     << " Check last " << checkHits
412     << " before saving the hit\n"
413     << " Correct TOF globally by " << correctT
414     << " ns (Flag =" << corrTOFBeam << ")\n"
415     << " Save hits recorded before " << tmaxHit
416     << " ns and if energy is above " << eminHit/MeV
417     << " MeV (for depth 0) or " << eminHitD/MeV
418     << " MeV (for nonzero depths); Time Slice Unit "
419     << timeSlice << " Ignore TrackID Flag " << ignoreTrackID;
420     }
421    
422     CaloSD::~CaloSD() {
423    
424     if(debugCaloSD) cout<<"save g4t "<<endl;
425    
426     if(saveTree){
427     g4f->cd();
428     g4t->Write();
429     g4f->Write() ;
430     g4f->Close() ;
431     }
432     saveTree = false;
433    
434    
435     if (slave) delete slave;
436     if (theHC) delete theHC;
437     if (meanResponse) delete meanResponse;
438     }
439    
440     bool CaloSD::ProcessHits(G4Step * aStep, G4TouchableHistory * ) {
441    
442     NaNTrap( aStep ) ;
443    
444     if (aStep == NULL) {
445     return true;
446     } else {
447     if (getStepInfo(aStep)) {
448     if (hitExists() == false && edepositEM+edepositHAD>0.)
449     currentHit = createNewHit();
450     }
451     }
452     return true;
453     }
454    
455     bool CaloSD::ProcessHits(G4GFlashSpot* aSpot, G4TouchableHistory*) {
456    
457     if (aSpot != NULL) {
458     theTrack = const_cast<G4Track *>(aSpot->GetOriginatorTrack()->GetPrimaryTrack());
459     G4int particleCode = theTrack->GetDefinition()->GetPDGEncoding();
460    
461     if (particleCode == emPDG ||
462     particleCode == epPDG ||
463     particleCode == gammaPDG ) {
464     edepositEM = aSpot->GetEnergySpot()->GetEnergy();
465     edepositHAD = 0.;
466     } else {
467     edepositEM = 0.;
468     edepositHAD = 0.;
469     }
470    
471     if (edepositEM>0.) {
472     G4Step * fFakeStep = new G4Step();
473     preStepPoint = fFakeStep->GetPreStepPoint();
474     G4StepPoint * fFakePostStepPoint = fFakeStep->GetPostStepPoint();
475     preStepPoint->SetPosition(aSpot->GetPosition());
476     fFakePostStepPoint->SetPosition(aSpot->GetPosition());
477    
478     G4TouchableHandle fTouchableHandle = aSpot->GetTouchableHandle();
479     preStepPoint->SetTouchableHandle(fTouchableHandle);
480     fFakeStep->SetTotalEnergyDeposit(aSpot->GetEnergySpot()->GetEnergy());
481    
482     double time = 0;
483     unsigned int unitID = setDetUnitId(fFakeStep);
484     int primaryID = getTrackID(theTrack);
485     uint16_t depth = getDepth(fFakeStep);
486    
487     if (unitID > 0) {
488     currentID.setID(unitID, time, primaryID, depth);
489     #ifdef DebugLog
490     LogDebug("CaloSim") << "CaloSD:: GetSpotInfo for"
491     << " Unit 0x" << std::hex << currentID.unitID()
492     << std::dec << " Edeposit = " << edepositEM << " "
493     << edepositHAD;
494     #endif
495     // Update if in the same detector, time-slice and for same track
496     if (currentID == previousID) {
497     updateHit(currentHit);
498     } else {
499     posGlobal = aSpot->GetEnergySpot()->GetPosition();
500     // Reset entry point for new primary
501     if (currentID.trackID() != previousID.trackID()) {
502     entrancePoint = aSpot->GetPosition();
503     entranceLocal = aSpot->GetTouchableHandle()->GetHistory()->
504     GetTopTransform().TransformPoint(entrancePoint);
505     incidentEnergy = theTrack->GetKineticEnergy();
506     #ifdef DebugLog
507     LogDebug("CaloSim") << "CaloSD: Incident energy "
508     << incidentEnergy/GeV << " GeV and"
509     << " entrance point " << entrancePoint
510     << " (Global) " << entranceLocal << " (Local)";
511     #endif
512     }
513    
514     if (checkHit() == false) currentHit = createNewHit();
515     }
516     }
517     delete fFakeStep;
518     }
519     return true;
520     }
521     return false;
522     }
523    
524     double CaloSD::getEnergyDeposit(G4Step* aStep) {
525     return aStep->GetTotalEnergyDeposit();
526     }
527    
528     void CaloSD::Initialize(G4HCofThisEvent * HCE) {
529     totalHits = 0;
530    
531     #ifdef DebugLog
532     edm::LogInfo("CaloSim") << "CaloSD : Initialize called for " << GetName();
533     #endif
534    
535     //This initialization is performed at the beginning of an event
536     //------------------------------------------------------------
537     theHC = new CaloG4HitCollection(GetName(), collectionName[0]);
538    
539     if (hcID<0) hcID = G4SDManager::GetSDMpointer()->GetCollectionID(collectionName[0]);
540     HCE->AddHitsCollection(hcID, theHC);
541     }
542    
543     void CaloSD::EndOfEvent(G4HCofThisEvent* ) {
544     // clean the hits for the last tracks
545    
546    
547     if(debugCaloSD) cout<<"CaloSD::EndOfEvent " << csdname<<" "<<ncounter2 <<endl;
548    
549     if(ncounter2%6==0){
550     if(debugCaloSD) cout<<"filltree " <<estepg4EEm[85][58].size() <<endl;
551    
552     ng4EB = 0;
553     for(int j=-85; j<= 85; j++){
554     if(j==0) continue;
555     for(int k = 1; k<= 360; k++){
556     int nhit = estepg4EB[j][k].size();
557     if(nhit<1) continue;
558     vector<float> temp;
559     vector<float> temp1;
560     vector<int> temp2;
561     vector<int> temp2a;
562     vector<int> temp2b;
563    
564     vector<float> temp3;
565     vector<float> temp4;
566     vector<float> temp5;
567     vector<int> temp6;
568     vector<int> temp7;
569    
570     vector<float> temp8;
571     vector<float> temp9;
572     vector<float> temp10;
573    
574     vector<float> temp11;
575    
576     vector<float> temp12;
577     vector<float> temp13;
578     vector<float> temp14;
579     vector<float> temp15;
580    
581    
582    
583     double esum = 0;
584     double tmin = 1E9;
585     int indtmin = 0;
586     for(int n=0; n<nhit; n++){
587    
588     esum += estepg4EB[j][k][n];
589    
590     if(tmin>tstepg4EB[j][k][n]){
591     tmin = tstepg4EB[j][k][n];
592     indtmin = n;
593     }
594    
595     if( saveOnlyEnterCur==1){
596     if(enterstepg4EB[j][k][n]==0) continue;
597     }
598    
599     temp.push_back(estepg4EB[j][k][n]);
600     temp1.push_back(tstepg4EB[j][k][n]);
601     temp2.push_back(idstepg4EB[j][k][n]);
602     temp2a.push_back(pidstepg4EB[j][k][n]);
603     temp2b.push_back(parentidstepg4EB[j][k][n]);
604     temp3.push_back(postxstepg4EB[j][k][n]);
605     temp4.push_back(postystepg4EB[j][k][n]);
606     temp5.push_back(postzstepg4EB[j][k][n]);
607    
608     temp6.push_back(enterstepg4EB[j][k][n]);
609     temp7.push_back(leavestepg4EB[j][k][n]);
610    
611     temp8.push_back(prexstepg4EB[j][k][n]);
612     temp9.push_back(preystepg4EB[j][k][n]);
613     temp10.push_back(prezstepg4EB[j][k][n]);
614     temp11.push_back(pretstepg4EB[j][k][n]);
615    
616     temp12.push_back(preestepg4EB[j][k][n]);
617    
618    
619    
620     }
621    
622     ietag4EB[ng4EB] = j;
623     iphig4EB[ng4EB] = k;
624     esumg4EB[ng4EB] = esum/1000;
625     tming4EB[ng4EB] = tmin;
626     xtming4EB[ng4EB] = postxstepg4EB[j][k][indtmin];
627     ytming4EB[ng4EB] = postystepg4EB[j][k][indtmin];
628     ztming4EB[ng4EB] = postzstepg4EB[j][k][indtmin];
629    
630    
631     eg4EB->push_back(temp);
632     tg4EB->push_back(temp1);
633     idg4EB->push_back(temp2);
634     pidg4EB->push_back(temp2a);
635     parentidg4EB->push_back(temp2b);
636     postxg4EB->push_back(temp3);
637     postyg4EB->push_back(temp4);
638     postzg4EB->push_back(temp5);
639    
640     enterg4EB->push_back(temp6);
641     leaveg4EB->push_back(temp7);
642    
643     prexg4EB->push_back(temp8);
644     preyg4EB->push_back(temp9);
645     prezg4EB->push_back(temp10);
646     pretg4EB->push_back(temp11);
647    
648     preeg4EB->push_back(temp12);
649    
650    
651    
652     ng4EB++;
653     }
654     }
655    
656     ng4EE = 0;
657     for(int j=1;j<=100; j++){
658     for(int k=1; k<=100; k++){
659     int nhit = estepg4EEm[j][k].size();
660     if(nhit<1) continue;
661     vector<float> temp;
662     vector<float> temp1;
663     vector<int> temp2;
664     vector<int> temp2a;
665     vector<int> temp2b;
666    
667     vector<float> temp3;
668     vector<float> temp4;
669     vector<float> temp5;
670    
671     vector<int> temp6;
672     vector<int> temp7;
673    
674     vector<float> temp8;
675     vector<float> temp9;
676     vector<float> temp10;
677    
678     vector<float> temp11;
679    
680     vector<float> temp12;
681     vector<float> temp13;
682     vector<float> temp14;
683     vector<float> temp15;
684    
685    
686     double esum = 0;
687     double tmin = 1E9;
688     int indtmin = 0;
689     for(int n=0; n<nhit; n++){
690    
691     esum += estepg4EEm[j][k][n];
692     if( tmin > tstepg4EEm[j][k][n]){
693     tmin = tstepg4EEm[j][k][n];
694     indtmin = n;
695     }
696    
697    
698     if( saveOnlyEnterCur==1){
699     if(enterstepg4EEm[j][k][n]==0) continue;
700     }
701    
702     temp.push_back(estepg4EEm[j][k][n]);
703     temp1.push_back(tstepg4EEm[j][k][n]);
704     temp2.push_back(idstepg4EEm[j][k][n]);
705     temp2a.push_back(pidstepg4EEm[j][k][n]);
706     temp2b.push_back(parentidstepg4EEm[j][k][n]);
707     temp3.push_back(postxstepg4EEm[j][k][n]);
708     temp4.push_back(postystepg4EEm[j][k][n]);
709     temp5.push_back(postzstepg4EEm[j][k][n]);
710    
711     temp6.push_back(enterstepg4EEm[j][k][n]);
712     temp7.push_back(leavestepg4EEm[j][k][n]);
713    
714     temp8.push_back(prexstepg4EEm[j][k][n]);
715     temp9.push_back(preystepg4EEm[j][k][n]);
716     temp10.push_back(prezstepg4EEm[j][k][n]);
717     temp11.push_back(pretstepg4EEm[j][k][n]);
718    
719     temp12.push_back(preestepg4EEm[j][k][n]);
720    
721    
722    
723     }
724    
725     ixg4EE[ng4EE] = j;
726     iyg4EE[ng4EE] = k;
727     izg4EE[ng4EE] = -1;
728     esumg4EE[ng4EE] = esum/1000;
729     tming4EE[ng4EE] = tmin;
730     xtming4EE[ng4EE] = postxstepg4EEm[j][k][indtmin];
731     ytming4EE[ng4EE] = postystepg4EEm[j][k][indtmin];
732     ztming4EE[ng4EE] = postzstepg4EEm[j][k][indtmin];
733    
734    
735     eg4EE->push_back(temp);
736     tg4EE->push_back(temp1);
737     idg4EE->push_back(temp2);
738     pidg4EE->push_back(temp2a);
739     parentidg4EE->push_back(temp2b);
740     postxg4EE->push_back(temp3);
741     postyg4EE->push_back(temp4);
742     postzg4EE->push_back(temp5);
743    
744     enterg4EE->push_back(temp6);
745     leaveg4EE->push_back(temp7);
746     prexg4EE->push_back(temp8);
747     preyg4EE->push_back(temp9);
748     prezg4EE->push_back(temp10);
749     pretg4EE->push_back(temp11);
750    
751     preeg4EE->push_back(temp12);
752    
753    
754    
755     ng4EE++;
756     }
757     }
758     for(int j=1;j<=100; j++){
759     for(int k=1; k<=100; k++){
760     int nhit = estepg4EEp[j][k].size();
761     if(nhit<1) continue;
762     vector<float> temp;
763     vector<float> temp1;
764     vector<int> temp2;
765     vector<int> temp2a;
766     vector<int> temp2b;
767    
768     vector<float> temp3;
769     vector<float> temp4;
770     vector<float> temp5;
771    
772     vector<int> temp6;
773     vector<int> temp7;
774    
775     vector<float> temp8;
776     vector<float> temp9;
777     vector<float> temp10;
778    
779     vector<float> temp11;
780    
781     vector<float> temp12;
782     vector<float> temp13;
783     vector<float> temp14;
784     vector<float> temp15;
785    
786    
787     double esum = 0;
788     double tmin = 1E9;
789     int indtmin = 0;
790     for(int n=0; n<nhit; n++){
791    
792     esum += estepg4EEp[j][k][n];
793     if(tmin > tstepg4EEp[j][k][n]){
794     tmin = tstepg4EEp[j][k][n];
795     indtmin = n;
796     }
797    
798     if( saveOnlyEnterCur==1){
799     if(enterstepg4EEp[j][k][n]==0) continue;
800     }
801    
802     temp.push_back(estepg4EEp[j][k][n]);
803     temp1.push_back(tstepg4EEp[j][k][n]);
804     temp2.push_back(idstepg4EEp[j][k][n]);
805     temp2a.push_back(pidstepg4EEp[j][k][n]);
806     temp2b.push_back(parentidstepg4EEp[j][k][n]);
807     temp3.push_back(postxstepg4EEp[j][k][n]);
808     temp4.push_back(postystepg4EEp[j][k][n]);
809     temp5.push_back(postzstepg4EEp[j][k][n]);
810     temp6.push_back(enterstepg4EEp[j][k][n]);
811     temp7.push_back(leavestepg4EEp[j][k][n]);
812     temp8.push_back(prexstepg4EEp[j][k][n]);
813     temp9.push_back(preystepg4EEp[j][k][n]);
814     temp10.push_back(prezstepg4EEp[j][k][n]);
815     temp11.push_back(pretstepg4EEp[j][k][n]);
816    
817    
818     temp12.push_back(preestepg4EEp[j][k][n]);
819    
820    
821    
822     }
823     ixg4EE[ng4EE] = j;
824     iyg4EE[ng4EE] = k;
825     izg4EE[ng4EE] = 1;
826     esumg4EE[ng4EE] = esum/1000;
827     tming4EE[ng4EE] = tmin;
828     xtming4EE[ng4EE] = postxstepg4EEp[j][k][indtmin];
829     ytming4EE[ng4EE] = postystepg4EEp[j][k][indtmin];
830     ztming4EE[ng4EE] = postzstepg4EEp[j][k][indtmin];
831    
832     eg4EE->push_back(temp);
833     tg4EE->push_back(temp1);
834     idg4EE->push_back(temp2);
835     pidg4EE->push_back(temp2a);
836     parentidg4EE->push_back(temp2b);
837     postxg4EE->push_back(temp3);
838     postyg4EE->push_back(temp4);
839     postzg4EE->push_back(temp5);
840    
841     enterg4EE->push_back(temp6);
842     leaveg4EE->push_back(temp7);
843     prexg4EE->push_back(temp8);
844     preyg4EE->push_back(temp9);
845     prezg4EE->push_back(temp10);
846     pretg4EE->push_back(temp11);
847    
848     preeg4EE->push_back(temp12);
849    
850    
851    
852    
853     ng4EE++;
854     }
855     }
856    
857     if(debugCaloSD) cout<<"filling tree " << ng4EB <<" "<< eg4EB->size() <<" "<< ng4EE<<" "<<eg4EE->size()<<endl;
858     g4t->Fill();
859    
860    
861     }
862     ncounter2 ++;
863     if(ncounter2 ==6){
864     ncounter2 = 0;
865     }
866    
867    
868     cleanHitCollection();
869    
870     #ifdef DebugLog
871     edm::LogInfo("CaloSim") << "CaloSD: EndofEvent entered with " << theHC->entries()
872     << " entries";
873     #endif
874     // TimeMe("CaloSD:sortAndMergeHits",false);
875     }
876    
877     void CaloSD::clear() {}
878    
879     void CaloSD::DrawAll() {}
880    
881     void CaloSD::PrintAll() {
882     #ifdef DebugLog
883     edm::LogInfo("CaloSim") << "CaloSD: Collection " << theHC->GetName();
884     #endif
885     theHC->PrintAllHits();
886     }
887    
888     void CaloSD::fillHits(edm::PCaloHitContainer& c, std::string n) {
889     if (slave->name() == n) c=slave->hits();
890     slave->Clean();
891     }
892    
893     bool CaloSD::getStepInfo(G4Step* aStep) {
894    
895     preStepPoint = aStep->GetPreStepPoint();
896     theTrack = aStep->GetTrack();
897    
898     G4int particleCode = theTrack->GetDefinition()->GetPDGEncoding();
899     if (particleCode == emPDG ||
900     particleCode == epPDG ||
901     particleCode == gammaPDG ) {
902     edepositEM = getEnergyDeposit(aStep);
903     edepositHAD = 0.;
904     } else {
905     edepositEM = 0.;
906     edepositHAD = getEnergyDeposit(aStep);
907     }
908    
909     double time = (aStep->GetPostStepPoint()->GetGlobalTime())/nanosecond;
910     unsigned int unitID= setDetUnitId(aStep);
911     uint16_t depth = getDepth(aStep);
912     int primaryID = getTrackID(theTrack);
913     int parentID = theTrack->GetParentID();
914    
915    
916     bool flag = (unitID > 0);
917    
918     double timepre = (aStep->GetPreStepPoint()->GetGlobalTime())/nanosecond;
919     bool debug = true;
920    
921     if(flag && debug){
922    
923     G4ThreeVector pos1 = preStepPoint->GetPosition();
924     G4ThreeVector pos2 = aStep->GetPostStepPoint()->GetPosition();
925    
926     const DetId detId(unitID);
927     if(!detId.null() && detId.det() == DetId::Ecal){
928    
929     G4TouchableHandle touch1 = preStepPoint->GetTouchableHandle();
930     G4VPhysicalVolume* volume = touch1->GetVolume();
931     G4String name = volume->GetName();
932     int enterCurrent = 0;
933     G4String name2 = preStepPoint->GetPhysicalVolume()->GetLogicalVolume()->GetName();
934     G4Material* material = preStepPoint ->GetMaterial();
935     G4String matname = material->GetName();
936    
937     if (preStepPoint->GetStepStatus() == fGeomBoundary) {
938     enterCurrent = 1;
939     }
940     int leaveCurrent = 0;
941     if( aStep->GetPostStepPoint()->GetStepStatus() == fGeomBoundary){
942     leaveCurrent = 1;
943     }
944    
945    
946     G4VPhysicalVolume* mother = touch1->GetVolume(depth=1);
947     G4String mothername = mother->GetName();
948    
949     if( EcalBarrel == detId.subdetId()){
950     EBDetId eb(detId);
951     int x = eb.ieta();
952     int y = eb.iphi();
953     estepg4EB[x][y].push_back(edepositEM);
954     tstepg4EB[x][y].push_back(time);
955     idstepg4EB[x][y].push_back(primaryID);
956     pidstepg4EB[x][y].push_back(particleCode);
957     parentidstepg4EB[x][y].push_back(parentID);
958     enterstepg4EB[x][y].push_back(enterCurrent);
959     leavestepg4EB[x][y].push_back(leaveCurrent);
960    
961     postxstepg4EB[x][y].push_back(pos2.x());
962     postystepg4EB[x][y].push_back(pos2.y());
963     postzstepg4EB[x][y].push_back(pos2.z());
964     prexstepg4EB[x][y].push_back(pos1.x());
965     preystepg4EB[x][y].push_back(pos1.y());
966     prezstepg4EB[x][y].push_back(pos1.z());
967     pretstepg4EB[x][y].push_back(timepre);
968    
969     preestepg4EB[x][y].push_back(aStep->GetPreStepPoint()->GetTotalEnergy());
970    
971    
972     if(debugCaloSD && enterCurrent ){
973     if(x==8 && y==9){
974     cout<<"eb1xtal "<<enterCurrent<<" "<<aStep->GetPreStepPoint()->GetTotalEnergy()<<" "<< name <<" "<<name2<<" "
975     <<matname<<" "<<mothername<<" ch/m "<< preStepPoint->GetCharge()<<" "<<preStepPoint->GetMass()<<" "<<particleCode <<" "<<parentID<<endl;
976     }
977     if(x==9 && y==9){
978     cout<<"eb2xtal "<<enterCurrent<<" "<<aStep->GetPreStepPoint()->GetTotalEnergy()<<" "<< name <<" "<<name2<<" "<<" "<<matname<<" "<<mothername<<endl;
979     }
980    
981     }
982    
983    
984     }else if(EcalEndcap == detId.subdetId() ) {
985     EEDetId ee(detId);
986     int x = ee.ix();
987     int y = ee.iy();
988    
989     if( ee.zside() == -1 ){
990     estepg4EEm[x][y].push_back(edepositEM);
991     tstepg4EEm[x][y].push_back(time);
992     idstepg4EEm[x][y].push_back(primaryID);
993     pidstepg4EEm[x][y].push_back(particleCode);
994     enterstepg4EEm[x][y].push_back(enterCurrent);
995     leavestepg4EEm[x][y].push_back(leaveCurrent);
996     parentidstepg4EEm[x][y].push_back(parentID);
997    
998     postxstepg4EEm[x][y].push_back(pos2.x());
999     postystepg4EEm[x][y].push_back(pos2.y());
1000     postzstepg4EEm[x][y].push_back(pos2.z());
1001     prexstepg4EEm[x][y].push_back(pos1.x());
1002     preystepg4EEm[x][y].push_back(pos1.y());
1003     prezstepg4EEm[x][y].push_back(pos1.z());
1004     pretstepg4EEm[x][y].push_back(timepre);
1005    
1006     preestepg4EEm[x][y].push_back(aStep->GetPreStepPoint()->GetTotalEnergy());
1007    
1008    
1009    
1010     if(debugCaloSD && ee.ix() == 85 && ee.iy() == 58 && ee.zside() == -1 ) {
1011     std::cout<<"ee_step_xtal1: "<< edepositEM<<" "<<time<<" "<<primaryID<<" "<<pos1.x()<<" "<<pos1.y()<<" "<<pos1.z()<<" "<<pos2.x()<<" "<<pos2.y()<<" "<<pos2.z()<<" "<<name<<endl;
1012     if(enterCurrent==1){ //theTrack->pos is the same as pos2 and theTrack->GetGlobaltime is the same as time
1013     /// theTrack->GetTotalEnergy() == aStep->GetPostStepPoint()->GetTotalEnergy();
1014     /// aStep->GetPreStepPoint()->GetTotalEnergy() is larger, the delta is ~ the deposted energy of this step
1015     cout<<"xtal1_entercur "<< pos1.x()<<" "<<pos1.y()<<" "<<pos1.z()<<" "<< pos2.x()<<" "<<pos2.y()<<" "<<pos2.z()<<" "<<name<<" en "<<theTrack->GetTotalEnergy()<<" "<<theTrack->GetKineticEnergy()<<" trkt "<<theTrack->GetGlobalTime()<<" trkpos "<< theTrack->GetPosition().x()<<" "<<theTrack->GetPosition().y()<<" "<<theTrack->GetPosition().z()<<" etot "<< aStep->GetPostStepPoint()->GetTotalEnergy()<<" "<<aStep->GetPreStepPoint()->GetTotalEnergy()<<" "<<endl;
1016     }
1017     if(leaveCurrent==1){
1018     cout<<"xtal1_leavecur "<< pos1.x()<<" "<<pos1.y()<<" "<<pos1.z()<<" "<< pos2.x()<<" "<<pos2.y()<<" "<<pos2.z()<<" "<<name<<endl;
1019     }
1020    
1021     }
1022    
1023    
1024     if(debugCaloSD && ee.ix() == 85 && ee.iy() == 59 && ee.zside() == -1 ) std::cout<<"ee_step_xtal2: "<< edepositEM<<" "<<time<<" "<<primaryID<<" "<<pos1.x()<<" "<<pos1.y()<<" "<<pos1.z()<<" "<<pos2.x()<<" "<<pos2.y()<<" "<<pos2.z()<<" "<<name<<endl;
1025    
1026     }else{
1027     estepg4EEp[x][y].push_back(edepositEM);
1028     tstepg4EEp[x][y].push_back(time);
1029     idstepg4EEp[x][y].push_back(primaryID);
1030     pidstepg4EEp[x][y].push_back(particleCode);
1031     enterstepg4EEp[x][y].push_back(enterCurrent);
1032     leavestepg4EEp[x][y].push_back(leaveCurrent);
1033     parentidstepg4EEp[x][y].push_back(parentID);
1034    
1035     postxstepg4EEp[x][y].push_back(pos2.x());
1036     postystepg4EEp[x][y].push_back(pos2.y());
1037     postzstepg4EEp[x][y].push_back(pos2.z());
1038     prexstepg4EEp[x][y].push_back(pos1.x());
1039     preystepg4EEp[x][y].push_back(pos1.y());
1040     prezstepg4EEp[x][y].push_back(pos1.z());
1041     pretstepg4EEp[x][y].push_back(timepre);
1042     preestepg4EEp[x][y].push_back(aStep->GetPreStepPoint()->GetTotalEnergy());
1043    
1044    
1045     }
1046    
1047    
1048     }
1049     }
1050     }
1051    
1052    
1053    
1054     if (flag) {
1055     currentID.setID(unitID, time, primaryID, depth);
1056     #ifdef DebugLog
1057     G4TouchableHistory* touch =(G4TouchableHistory*)(theTrack->GetTouchable());
1058     LogDebug("CaloSim") << "CaloSD:: GetStepInfo for"
1059     << " PV " << touch->GetVolume(0)->GetName()
1060     << " PVid = " << touch->GetReplicaNumber(0)
1061     << " MVid = " << touch->GetReplicaNumber(1)
1062     << " Unit " << currentID.unitID()
1063     << " Edeposit = " << edepositEM << " " << edepositHAD;
1064     } else {
1065     G4TouchableHistory* touch =(G4TouchableHistory*)(theTrack->GetTouchable());
1066     LogDebug("CaloSim") << "CaloSD:: GetStepInfo for"
1067     << " PV " << touch->GetVolume(0)->GetName()
1068     << " PVid = " << touch->GetReplicaNumber(0)
1069     << " MVid = " << touch->GetReplicaNumber(1)
1070     << " Unit " << std::hex << unitID << std::dec
1071     << " Edeposit = " << edepositEM << " " << edepositHAD;
1072     #endif
1073     }
1074     return flag;
1075     }
1076    
1077     G4ThreeVector CaloSD::setToLocal(G4ThreeVector global, const G4VTouchable* touch) {
1078    
1079     G4ThreeVector localPoint = touch->GetHistory()->GetTopTransform().TransformPoint(global);
1080    
1081     return localPoint;
1082     }
1083    
1084     G4bool CaloSD::hitExists() {
1085     #ifdef DebugLog
1086     if (currentID.trackID()<1)
1087     edm::LogWarning("CaloSim") << "***** CaloSD error: primaryID = "
1088     << currentID.trackID()
1089     << " maybe detector name changed";
1090     #endif
1091     // Update if in the same detector, time-slice and for same track
1092     if (currentID == previousID) {
1093     updateHit(currentHit);
1094     return true;
1095     }
1096    
1097     // Reset entry point for new primary
1098     posGlobal = preStepPoint->GetPosition();
1099     if (currentID.trackID() != previousID.trackID())
1100     resetForNewPrimary(preStepPoint->GetPosition(), preStepPoint->GetKineticEnergy());
1101    
1102     return checkHit();
1103     }
1104    
1105     G4bool CaloSD::checkHit() {
1106     //look in the HitContainer whether a hit with the same ID already exists:
1107     bool found = false;
1108     if (useMap) {
1109     std::map<CaloHitID,CaloG4Hit*>::const_iterator it = hitMap.find(currentID);
1110     if (it != hitMap.end()) {
1111     currentHit = it->second;
1112     found = true;
1113     }
1114     } else {
1115     if (checkHits <= 0) return false;
1116     int minhit= (theHC->entries()>checkHits ? theHC->entries()-checkHits : 0);
1117     int maxhit= theHC->entries()-1;
1118    
1119     for (int j=maxhit; j>minhit&&!found; --j) {
1120     if ((*theHC)[j]->getID() == currentID) {
1121     currentHit = (*theHC)[j];
1122     found = true;
1123     }
1124     }
1125     }
1126    
1127     if (found) {
1128     updateHit(currentHit);
1129     return true;
1130     } else {
1131     return false;
1132     }
1133     }
1134    
1135     int CaloSD::getNumberOfHits() { return theHC->entries(); }
1136    
1137     CaloG4Hit* CaloSD::createNewHit() {
1138     #ifdef DebugLog
1139     LogDebug("CaloSim") << "CaloSD::CreateNewHit for"
1140     << " Unit " << currentID.unitID()
1141     << " " << currentID.depth()
1142     << " Edeposit = " << edepositEM << " " << edepositHAD;
1143     LogDebug("CaloSim") << " primary " << currentID.trackID()
1144     << " time slice " << currentID.timeSliceID()
1145     << " For Track " << theTrack->GetTrackID()
1146     << " which is a " <<theTrack->GetDefinition()->GetParticleName()
1147     << " of energy " << theTrack->GetKineticEnergy()/GeV
1148     << " " << theTrack->GetMomentum().mag()/GeV
1149     << " daughter of part. " << theTrack->GetParentID()
1150     << " and created by " ;
1151    
1152     if (theTrack->GetCreatorProcess()!=NULL)
1153     LogDebug("CaloSim") << theTrack->GetCreatorProcess()->GetProcessName() ;
1154     else
1155     LogDebug("CaloSim") << "NO process";
1156     #endif
1157    
1158     CaloG4Hit* aHit;
1159     if (reusehit.size() > 0) {
1160     aHit = reusehit[0];
1161     aHit->setEM(0.);
1162     aHit->setHadr(0.);
1163     reusehit.erase(reusehit.begin());
1164     } else {
1165     aHit = new CaloG4Hit;
1166     }
1167    
1168     aHit->setID(currentID);
1169     aHit->setEntry(entrancePoint.x(),entrancePoint.y(),entrancePoint.z());
1170     aHit->setEntryLocal(entranceLocal.x(),entranceLocal.y(),entranceLocal.z());
1171     aHit->setPosition(posGlobal.x(),posGlobal.y(),posGlobal.z());
1172     aHit->setIncidentEnergy(incidentEnergy);
1173     updateHit(aHit);
1174    
1175     storeHit(aHit);
1176     double etrack = 0;
1177     if (currentID.trackID() == primIDSaved) { // The track is saved; nothing to be done
1178     } else if (currentID.trackID() == theTrack->GetTrackID()) {
1179     etrack= theTrack->GetKineticEnergy();
1180     //edm::LogInfo("CaloSim") << "CaloSD: set save the track " << currentID.trackID()
1181     // << " etrack " << etrack << " eCut " << energyCut << " flag " << forceSave;
1182     if (etrack >= energyCut || forceSave) {
1183     TrackInformation* trkInfo = (TrackInformation *)(theTrack->GetUserInformation());
1184     trkInfo->storeTrack(true);
1185     trkInfo->putInHistory();
1186     // trkInfo->setAncestor();
1187     #ifdef DebugLog
1188     LogDebug("CaloSim") << "CaloSD: set save the track " << currentID.trackID()
1189     << " with Hit";
1190     #endif
1191     }
1192     } else {
1193     TrackWithHistory * trkh = tkMap[currentID.trackID()];
1194     #ifdef DebugLog
1195     LogDebug("CaloSim") << "CaloSD : TrackwithHistory pointer for "
1196     << currentID.trackID() << " is " << trkh;
1197     #endif
1198     if (trkh != NULL) {
1199     etrack = sqrt(trkh->momentum().Mag2());
1200     if (etrack >= energyCut) {
1201     trkh->save();
1202     #ifdef DebugLog
1203     LogDebug("CaloSim") << "CaloSD: set save the track "
1204     << currentID.trackID() << " with Hit";
1205     #endif
1206     }
1207     }
1208     }
1209     primIDSaved = currentID.trackID();
1210     if (useMap) totalHits++;
1211     return aHit;
1212     }
1213    
1214     void CaloSD::updateHit(CaloG4Hit* aHit) {
1215     if (edepositEM+edepositHAD != 0) {
1216     aHit->addEnergyDeposit(edepositEM,edepositHAD);
1217     #ifdef DebugLog
1218     LogDebug("CaloSim") << "CaloSD: Add energy deposit in " << currentID
1219     << " em " << edepositEM/MeV << " hadronic "
1220     << edepositHAD/MeV << " MeV";
1221     #endif
1222     }
1223    
1224     // buffer for next steps:
1225     previousID = currentID;
1226     }
1227    
1228     void CaloSD::resetForNewPrimary(G4ThreeVector point, double energy) {
1229     entrancePoint = point;
1230     entranceLocal = setToLocal(entrancePoint, preStepPoint->GetTouchable());
1231     incidentEnergy = energy;
1232     #ifdef DebugLog
1233     LogDebug("CaloSim") << "CaloSD: Incident energy " << incidentEnergy/GeV
1234     << " GeV and" << " entrance point " << entrancePoint
1235     << " (Global) " << entranceLocal << " (Local)";
1236     #endif
1237     }
1238    
1239     double CaloSD::getAttenuation(G4Step* aStep, double birk1, double birk2, double birk3) {
1240     double weight = 1.;
1241     double charge = aStep->GetPreStepPoint()->GetCharge();
1242    
1243     if (charge != 0. && aStep->GetStepLength() > 0) {
1244     G4Material* mat = aStep->GetPreStepPoint()->GetMaterial();
1245     double density = mat->GetDensity();
1246     double dedx = aStep->GetTotalEnergyDeposit()/aStep->GetStepLength();
1247     double rkb = birk1/density;
1248     double c = birk2*rkb*rkb;
1249     if (std::abs(charge) >= 2.) rkb /= birk3; // based on alpha particle data
1250     weight = 1./(1.+rkb*dedx+c*dedx*dedx);
1251     #ifdef DebugLog
1252     LogDebug("CaloSim") << "CaloSD::getAttenuation in " << mat->GetName()
1253     << " Charge " << charge << " dE/dx " << dedx
1254     << " Birk Const " << rkb << ", " << c << " Weight = "
1255     << weight << " dE " << aStep->GetTotalEnergyDeposit();
1256     #endif
1257     }
1258     return weight;
1259     }
1260    
1261     void CaloSD::update(const BeginOfRun *) {
1262    
1263    
1264     if(debugCaloSD) cout<<"CaloSD::update BeginOfRun" << csdname << endl;
1265    
1266    
1267     G4ParticleTable * theParticleTable = G4ParticleTable::GetParticleTable();
1268     G4String particleName;
1269     emPDG = theParticleTable->FindParticle(particleName="e-")->GetPDGEncoding();
1270     epPDG = theParticleTable->FindParticle(particleName="e+")->GetPDGEncoding();
1271     gammaPDG = theParticleTable->FindParticle(particleName="gamma")->GetPDGEncoding();
1272     #ifdef DebugLog
1273     LogDebug("CaloSim") << "CaloSD: Particle code for e- = " << emPDG
1274     << " for e+ = " << epPDG << " for gamma = " << gammaPDG;
1275     #endif
1276     initRun();
1277     runInit = true;
1278     }
1279    
1280     void CaloSD::update(const BeginOfEvent *) {
1281    
1282     if(debugCaloSD) cout<<"CaloSD::update BeginOfEvent " << csdname<< " "<<ncounter1 <<endl;
1283    
1284     if(ncounter1%6==0){
1285     if(debugCaloSD)cout<<"here_to_cleartree "<< ncounter1 <<endl;
1286    
1287     for(int j=-85; j<= 85; j++){
1288     if(j==0) continue;
1289     for(int k = 1; k<= 360; k++){
1290     estepg4EB[j][k].clear();
1291     tstepg4EB[j][k].clear();
1292     idstepg4EB[j][k].clear();
1293     pidstepg4EB[j][k].clear();
1294     parentidstepg4EB[j][k].clear();
1295     enterstepg4EB[j][k].clear();
1296     leavestepg4EB[j][k].clear();
1297    
1298    
1299     postxstepg4EB[j][k].clear();
1300     postystepg4EB[j][k].clear();
1301     postzstepg4EB[j][k].clear();
1302    
1303     prexstepg4EB[j][k].clear();
1304     preystepg4EB[j][k].clear();
1305     prezstepg4EB[j][k].clear();
1306     pretstepg4EB[j][k].clear();
1307    
1308     preestepg4EB[j][k].clear();
1309    
1310    
1311     }
1312     }
1313    
1314     for(int j=1; j<=100; j++){
1315     for(int k=1; k<=100; k++){
1316     estepg4EEm[j][k].clear();
1317     tstepg4EEm[j][k].clear();
1318     estepg4EEp[j][k].clear();
1319     tstepg4EEp[j][k].clear();
1320    
1321     idstepg4EEm[j][k].clear();
1322     pidstepg4EEm[j][k].clear();
1323     parentidstepg4EEm[j][k].clear();
1324     enterstepg4EEm[j][k].clear();
1325     leavestepg4EEm[j][k].clear();
1326    
1327     prexstepg4EEm[j][k].clear();
1328     preystepg4EEm[j][k].clear();
1329     prezstepg4EEm[j][k].clear();
1330     pretstepg4EEm[j][k].clear();
1331    
1332     preestepg4EEm[j][k].clear();
1333    
1334    
1335    
1336     idstepg4EEp[j][k].clear();
1337     pidstepg4EEp[j][k].clear();
1338     parentidstepg4EEp[j][k].clear();
1339     enterstepg4EEp[j][k].clear();
1340     leavestepg4EEp[j][k].clear();
1341    
1342     postxstepg4EEm[j][k].clear();
1343     postystepg4EEm[j][k].clear();
1344     postzstepg4EEm[j][k].clear();
1345    
1346     postxstepg4EEp[j][k].clear();
1347     postystepg4EEp[j][k].clear();
1348     postzstepg4EEp[j][k].clear();
1349    
1350     prexstepg4EEp[j][k].clear();
1351     preystepg4EEp[j][k].clear();
1352     prezstepg4EEp[j][k].clear();
1353     pretstepg4EEp[j][k].clear();
1354    
1355     preestepg4EEp[j][k].clear();
1356    
1357    
1358    
1359     }
1360     }
1361    
1362     eg4EB->clear();
1363     tg4EB->clear();
1364     idg4EB->clear();
1365     pidg4EB->clear();
1366     parentidg4EB->clear();
1367     enterg4EB->clear();
1368     leaveg4EB->clear();
1369    
1370     postxg4EB->clear();
1371     postyg4EB->clear();
1372     postzg4EB->clear();
1373     prexg4EB->clear();
1374     preyg4EB->clear();
1375     prezg4EB->clear();
1376     pretg4EB->clear();
1377    
1378     preeg4EB->clear();
1379    
1380    
1381     eg4EE->clear();
1382     tg4EE->clear();
1383     idg4EE->clear();
1384     pidg4EE->clear();
1385     parentidg4EE->clear();
1386     postxg4EE->clear();
1387     postyg4EE->clear();
1388     postzg4EE->clear();
1389     enterg4EE->clear();
1390     leaveg4EE->clear();
1391     prexg4EE->clear();
1392     preyg4EE->clear();
1393     prezg4EE->clear();
1394     pretg4EE->clear();
1395    
1396     preeg4EE->clear();
1397    
1398    
1399    
1400     }
1401    
1402     ncounter1 ++;
1403     if(ncounter1==6){
1404     ncounter1 = 0;
1405     }
1406    
1407    
1408     #ifdef DebugLog
1409     LogDebug("CaloSim") << "CaloSD: Dispatched BeginOfEvent for " << GetName()
1410     << " !" ;
1411     #endif
1412     clearHits();
1413     }
1414    
1415     void CaloSD::update(const EndOfTrack * trk) {
1416     int id = (*trk)()->GetTrackID();
1417     TrackInformation *trkI =(TrackInformation *)((*trk)()->GetUserInformation());
1418     int lastTrackID = -1;
1419     if (trkI) lastTrackID = trkI->getIDonCaloSurface();
1420     if (id == lastTrackID) {
1421     const TrackContainer * trksForThisEvent = m_trackManager->trackContainer();
1422     if (trksForThisEvent != NULL) {
1423     int it = (int)(trksForThisEvent->size()) - 1;
1424     if (it >= 0) {
1425     TrackWithHistory * trkH = (*trksForThisEvent)[it];
1426     if (trkH->trackID() == (unsigned int)(id)) tkMap[id] = trkH;
1427     #ifdef DebugLog
1428     LogDebug("CaloSim") << "CaloSD: get track " << it << " from "
1429     << "Container of size " << trksForThisEvent->size()
1430     << " with ID " << trkH->trackID();
1431     } else {
1432     LogDebug("CaloSim") << "CaloSD: get track " << it << " from "
1433     << "Container of size " << trksForThisEvent->size()
1434     << " with no ID";
1435     #endif
1436     }
1437     }
1438     }
1439     }
1440    
1441     void CaloSD::update(const ::EndOfEvent * ) {
1442     int count = 0, wrong = 0;
1443     bool ok;
1444    
1445     slave->ReserveMemory(theHC->entries());
1446    
1447     for (int i=0; i<theHC->entries(); ++i) {
1448     ok = saveHit((*theHC)[i]);
1449     ++count;
1450     if (!ok) ++wrong;
1451     }
1452    
1453     edm::LogInfo("CaloSim") << "CaloSD: " << GetName() << " store " << count
1454     << " hits recorded with " << wrong
1455     << " track IDs not given properly and "
1456     << totalHits-count << " hits not passing cuts";
1457     summarize();
1458    
1459     tkMap.erase (tkMap.begin(), tkMap.end());
1460     }
1461    
1462     void CaloSD::clearHits() {
1463     if (useMap) hitMap.erase (hitMap.begin(), hitMap.end());
1464     for (unsigned int i = 0; i<reusehit.size(); ++i) delete reusehit[i];
1465     std::vector<CaloG4Hit*>().swap(reusehit);
1466     cleanIndex = 0;
1467     previousID.reset();
1468     primIDSaved = -99;
1469     #ifdef DebugLog
1470     LogDebug("CaloSim") << "CaloSD: Clears hit vector for " << GetName() << " " << slave;
1471     #endif
1472     slave->Initialize();
1473     #ifdef DebugLog
1474     LogDebug("CaloSim") << "CaloSD: Initialises slave SD for " << GetName();
1475     #endif
1476     }
1477    
1478     void CaloSD::initRun() {}
1479    
1480     int CaloSD::getTrackID(G4Track* aTrack) {
1481     int primaryID = 0;
1482     forceSave = false;
1483     TrackInformation* trkInfo=(TrackInformation *)(aTrack->GetUserInformation());
1484     if (trkInfo) {
1485     primaryID = trkInfo->getIDonCaloSurface();
1486     #ifdef DebugLog
1487     LogDebug("CaloSim") << "CaloSD: hit update from track Id on Calo Surface "
1488     << trkInfo->getIDonCaloSurface();
1489     #endif
1490     } else {
1491     primaryID = aTrack->GetTrackID();
1492     #ifdef DebugLog
1493     edm::LogWarning("CaloSim") << "CaloSD: Problem with primaryID **** set by "
1494     << "force to TkID **** " << primaryID << " in "
1495     << preStepPoint->GetTouchable()->GetVolume(0)->GetName();
1496     #endif
1497     }
1498     return primaryID;
1499     }
1500    
1501     uint16_t CaloSD::getDepth(G4Step*) { return 0; }
1502    
1503     bool CaloSD::filterHit(CaloG4Hit* hit, double time) {
1504     double emin(eminHit);
1505     if (hit->getDepth() > 0) emin = eminHitD;
1506     #ifdef DebugLog
1507     LogDebug("CaloSim") << "Depth " << hit->getDepth() << " Emin = " << emin << " ("
1508     << eminHit << ", " << eminHitD << ")";
1509     #endif
1510     return ((time <= tmaxHit) && (hit->getEnergyDeposit() > emin));
1511     }
1512    
1513     double CaloSD::getResponseWt(G4Track* aTrack) {
1514     if (meanResponse) {
1515     TrackInformation * trkInfo = (TrackInformation *)(aTrack->GetUserInformation());
1516     return meanResponse->getWeight(trkInfo->genParticlePID(), trkInfo->genParticleP());
1517     } else {
1518     return 1;
1519     }
1520     }
1521    
1522     void CaloSD::storeHit(CaloG4Hit* hit) {
1523     if (previousID.trackID()<0) return;
1524     if (hit == 0) {
1525     edm::LogWarning("CaloSim") << "CaloSD: hit to be stored is NULL !!";
1526     return;
1527     }
1528    
1529     theHC->insert(hit);
1530     if (useMap) hitMap.insert(std::pair<CaloHitID,CaloG4Hit*>(previousID,hit));
1531     }
1532    
1533     bool CaloSD::saveHit(CaloG4Hit* aHit) {
1534     int tkID;
1535     bool ok = true;
1536     if (m_trackManager) {
1537     tkID = m_trackManager->giveMotherNeeded(aHit->getTrackID());
1538     if (tkID == 0) {
1539     if (m_trackManager->trackExists(aHit->getTrackID())) tkID = (aHit->getTrackID());
1540     else {
1541     ok = false;
1542     }
1543     }
1544     } else {
1545     tkID = aHit->getTrackID();
1546     ok = false;
1547     }
1548     // edm::LogInfo("CaloSim") << "CalosD: Track ID " << aHit->getTrackID() << " changed to " << tkID << " by SimTrackManager" << " Status " << ok;
1549     #ifdef DebugLog
1550     LogDebug("CaloSim") << "CalosD: Track ID " << aHit->getTrackID()
1551     << " changed to " << tkID << " by SimTrackManager"
1552     << " Status " << ok;
1553     #endif
1554     double time = aHit->getTimeSlice();
1555     if (corrTOFBeam) time += correctT;
1556     slave->processHits(aHit->getUnitID(), aHit->getEM()/GeV,
1557     aHit->getHadr()/GeV, time, tkID, aHit->getDepth());
1558     #ifdef DebugLog
1559     LogDebug("CaloSim") << "CaloSD: Store Hit at " << std::hex
1560     << aHit->getUnitID() << std::dec << " "
1561     << aHit->getDepth() << " due to " << tkID
1562     << " in time " << time << " of energy "
1563     << aHit->getEM()/GeV << " GeV (EM) and "
1564     << aHit->getHadr()/GeV << " GeV (Hadr)";
1565     #endif
1566     return ok;
1567     }
1568    
1569     void CaloSD::summarize() {}
1570    
1571     void CaloSD::update(const BeginOfTrack * trk) {
1572     int primary = -1;
1573     TrackInformation * trkInfo = (TrackInformation *)((*trk)()->GetUserInformation());
1574     if ( trkInfo->isPrimary() ) primary = (*trk)()->GetTrackID();
1575    
1576     #ifdef DebugLog
1577     LogDebug("CaloSim") << "New track: isPrimary " << trkInfo->isPrimary()
1578     << " primary ID = " << primary
1579     << " primary ancestor ID " << primAncestor;
1580     #endif
1581    
1582     // update the information if a different primary track ID
1583    
1584     if (primary > 0 && primary != primAncestor) {
1585     primAncestor = primary;
1586    
1587     // clean the hits information
1588    
1589     if (theHC->entries()>0) cleanHitCollection();
1590    
1591     }
1592     }
1593    
1594     void CaloSD::cleanHitCollection() {
1595     std::vector<CaloG4Hit*>* theCollection = theHC->GetVector();
1596    
1597     #ifdef DebugLog
1598     LogDebug("CaloSim") << "CaloSD: collection before merging, size = " << theHC->entries();
1599     #endif
1600    
1601     selIndex.reserve(theHC->entries()-cleanIndex);
1602     if ( reusehit.size() == 0 ) reusehit.reserve(theHC->entries()-cleanIndex);
1603    
1604     // if no map used, merge before hits to have the save situation as a map
1605     if ( !useMap ) {
1606     hitvec.swap(*theCollection);
1607     sort((hitvec.begin()+cleanIndex), hitvec.end(), CaloG4HitLess());
1608     #ifdef DebugLog
1609     LogDebug("CaloSim") << "CaloSD::cleanHitCollection: sort hits in buffer "
1610     << "starting from element = " << cleanIndex;
1611     for (unsigned int i = 0; i<hitvec.size(); ++i)
1612     LogDebug("CaloSim")<<i<<" "<<*hitvec[i];
1613     #endif
1614     unsigned int i, j;
1615     CaloG4HitEqual equal;
1616     for (i=cleanIndex; i<hitvec.size(); ++i) {
1617     selIndex.push_back(i-cleanIndex);
1618     int jump = 0;
1619     for (j = i+1; j <hitvec.size() && equal(hitvec[i], hitvec[j]); ++j) {
1620     ++jump;
1621     // merge j to i
1622     (*hitvec[i]).addEnergyDeposit(*hitvec[j]);
1623     (*hitvec[j]).setEM(0.);
1624     (*hitvec[j]).setHadr(0.);
1625     reusehit.push_back(hitvec[j]);
1626     }
1627     i+=jump;
1628     }
1629     #ifdef DebugLog
1630     LogDebug("CaloSim") << "CaloSD: cleanHitCollection merge the hits in buffer ";
1631     for (unsigned int i = 0; i<hitvec.size(); ++i)
1632     LogDebug("CaloSim")<<i<<" "<<*hitvec[i];
1633     #endif
1634     for ( unsigned int i = cleanIndex; i < cleanIndex+selIndex.size(); ++i ) {
1635     hitvec[i] = hitvec[selIndex[i-cleanIndex]+cleanIndex];
1636     }
1637     hitvec.resize(cleanIndex+selIndex.size());
1638     #ifdef DebugLog
1639     LogDebug("CaloSim") << "CaloSD::cleanHitCollection: remove the merged hits in buffer,"
1640     << " new size = " << hitvec.size();
1641     for (unsigned int i = 0; i<hitvec.size(); ++i)
1642     LogDebug("CaloSim")<<i<<" "<<*hitvec[i];
1643     #endif
1644     hitvec.swap(*theCollection);
1645     std::vector<CaloG4Hit*>().swap(hitvec);
1646     selIndex.clear();
1647     totalHits = theHC->entries();
1648     }
1649    
1650     #ifdef DebugLog
1651     LogDebug("CaloSim") << "CaloSD: collection after merging, size = " << theHC->entries();
1652     #endif
1653    
1654     int addhit = 0;
1655    
1656     #ifdef DebugLog
1657     LogDebug("CaloSim") << "CaloSD: Size of reusehit after merge = " << reusehit.size();
1658     LogDebug("CaloSim") << "CaloSD: Starting hit selection from index = " << cleanIndex;
1659     #endif
1660    
1661     selIndex.reserve(theCollection->size()-cleanIndex);
1662     for (unsigned int i = cleanIndex; i<theCollection->size(); ++i) {
1663     CaloG4Hit* aHit((*theCollection)[i]);
1664    
1665     // selection
1666    
1667     double time = aHit->getTimeSlice();
1668     if (corrTOFBeam) time += correctT;
1669     if (!filterHit(aHit,time)) {
1670     #ifdef DebugLog
1671     LogDebug("CaloSim") << "CaloSD: dropped CaloG4Hit " << " " << *aHit;
1672     #endif
1673    
1674     // create the list of hits to be reused
1675    
1676     reusehit.push_back((*theCollection)[i]);
1677     ++addhit;
1678     } else {
1679     selIndex.push_back(i-cleanIndex);
1680     }
1681     }
1682    
1683     #ifdef DebugLog
1684     LogDebug("CaloSim") << "CaloSD: Size of reusehit after selection = " << reusehit.size()
1685     << " Number of added hit = " << addhit;
1686     #endif
1687     if (useMap) {
1688     if ( addhit>0 ) {
1689     int offset = reusehit.size()-addhit;
1690     for (int ii = addhit-1; ii>=0; --ii) {
1691     CaloHitID theID = reusehit[offset+ii]->getID();
1692     hitMap.erase(theID);
1693     }
1694     }
1695     }
1696     for (unsigned int j = 0; j<selIndex.size(); ++j) {
1697     (*theCollection)[cleanIndex+j] = (*theCollection)[cleanIndex+selIndex[j]];
1698     }
1699    
1700     theCollection->resize(cleanIndex+selIndex.size());
1701     std::vector<unsigned int>().swap(selIndex);
1702    
1703     #ifdef DebugLog
1704     LogDebug("CaloSim") << "CaloSD: hit collection after selection, size = "
1705     << theHC->entries();
1706     theHC->PrintAllHits();
1707     #endif
1708    
1709     cleanIndex = theHC->entries();
1710     }