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
Error occurred while calculating annotation data.
Log Message:
*** empty log message ***

File Contents

# Content
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 #include "TROOT.h"
28 #include "TFile.h"
29 #include "TTree.h"
30 #include "TBranch.h"
31 #include "TLorentzVector.h"
32 #include "TClonesArray.h"
33
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 }