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