1 | |
---|
2 | /** \class Calorimeter |
---|
3 | * |
---|
4 | * Fills calorimeter towers, performs calorimeter resolution smearing, |
---|
5 | * preselects towers hit by photons and creates energy flow objects. |
---|
6 | * |
---|
7 | * $Date: 2013-04-09 18:05:58 +0200 (Tue, 09 Apr 2013) $ |
---|
8 | * $Revision: 1086 $ |
---|
9 | * |
---|
10 | * |
---|
11 | * \author P. Demin - UCL, Louvain-la-Neuve |
---|
12 | * |
---|
13 | */ |
---|
14 | |
---|
15 | #include "modules/Calorimeter.h" |
---|
16 | |
---|
17 | #include "classes/DelphesClasses.h" |
---|
18 | #include "classes/DelphesFactory.h" |
---|
19 | #include "classes/DelphesFormula.h" |
---|
20 | |
---|
21 | #include "ExRootAnalysis/ExRootResult.h" |
---|
22 | #include "ExRootAnalysis/ExRootFilter.h" |
---|
23 | #include "ExRootAnalysis/ExRootClassifier.h" |
---|
24 | |
---|
25 | #include "TMath.h" |
---|
26 | #include "TString.h" |
---|
27 | #include "TFormula.h" |
---|
28 | #include "TRandom3.h" |
---|
29 | #include "TObjArray.h" |
---|
30 | #include "TDatabasePDG.h" |
---|
31 | #include "TLorentzVector.h" |
---|
32 | |
---|
33 | #include <algorithm> |
---|
34 | #include <stdexcept> |
---|
35 | #include <iostream> |
---|
36 | #include <sstream> |
---|
37 | |
---|
38 | using namespace std; |
---|
39 | |
---|
40 | //------------------------------------------------------------------------------ |
---|
41 | |
---|
42 | Calorimeter::Calorimeter() : |
---|
43 | fECalResolutionFormula(0), fHCalResolutionFormula(0), |
---|
44 | fItParticleInputArray(0), fItTrackInputArray(0), |
---|
45 | fTowerTrackArray(0), fItTowerTrackArray(0), |
---|
46 | fTowerPhotonArray(0), fItTowerPhotonArray(0) |
---|
47 | { |
---|
48 | fECalResolutionFormula = new DelphesFormula; |
---|
49 | fHCalResolutionFormula = new DelphesFormula; |
---|
50 | fTowerTrackArray = new TObjArray; |
---|
51 | fItTowerTrackArray = fTowerTrackArray->MakeIterator(); |
---|
52 | fTowerPhotonArray = new TObjArray; |
---|
53 | fItTowerPhotonArray = fTowerPhotonArray->MakeIterator(); |
---|
54 | } |
---|
55 | |
---|
56 | //------------------------------------------------------------------------------ |
---|
57 | |
---|
58 | Calorimeter::~Calorimeter() |
---|
59 | { |
---|
60 | if(fECalResolutionFormula) delete fECalResolutionFormula; |
---|
61 | if(fHCalResolutionFormula) delete fHCalResolutionFormula; |
---|
62 | if(fTowerTrackArray) delete fTowerTrackArray; |
---|
63 | if(fItTowerTrackArray) delete fItTowerTrackArray; |
---|
64 | if(fTowerPhotonArray) delete fTowerPhotonArray; |
---|
65 | if(fItTowerPhotonArray) delete fItTowerPhotonArray; |
---|
66 | } |
---|
67 | |
---|
68 | //------------------------------------------------------------------------------ |
---|
69 | |
---|
70 | void Calorimeter::Init() |
---|
71 | { |
---|
72 | ExRootConfParam param, paramEtaBins, paramPhiBins, paramFractions; |
---|
73 | Long_t i, j, k, size, sizeEtaBins, sizePhiBins, sizeFractions; |
---|
74 | Double_t ecalFraction, hcalFraction; |
---|
75 | TBinMap::iterator itEtaBin; |
---|
76 | set< Double_t >::iterator itPhiBin; |
---|
77 | vector< Double_t > *phiBins; |
---|
78 | |
---|
79 | // read eta and phi bins |
---|
80 | param = GetParam("EtaPhiBins"); |
---|
81 | size = param.GetSize(); |
---|
82 | fBinMap.clear(); |
---|
83 | fEtaBins.clear(); |
---|
84 | fPhiBins.clear(); |
---|
85 | for(i = 0; i < size/2; ++i) |
---|
86 | { |
---|
87 | paramEtaBins = param[i*2]; |
---|
88 | sizeEtaBins = paramEtaBins.GetSize(); |
---|
89 | paramPhiBins = param[i*2 + 1]; |
---|
90 | sizePhiBins = paramPhiBins.GetSize(); |
---|
91 | |
---|
92 | for(j = 0; j < sizeEtaBins; ++j) |
---|
93 | { |
---|
94 | for(k = 0; k < sizePhiBins; ++k) |
---|
95 | { |
---|
96 | fBinMap[paramEtaBins[j].GetDouble()].insert(paramPhiBins[k].GetDouble()); |
---|
97 | } |
---|
98 | } |
---|
99 | } |
---|
100 | |
---|
101 | // for better performance we transform map of sets to parallel vectors: |
---|
102 | // vector< double > and vector< vector< double >* > |
---|
103 | for(itEtaBin = fBinMap.begin(); itEtaBin != fBinMap.end(); ++itEtaBin) |
---|
104 | { |
---|
105 | fEtaBins.push_back(itEtaBin->first); |
---|
106 | phiBins = new vector< double >(itEtaBin->second.size()); |
---|
107 | fPhiBins.push_back(phiBins); |
---|
108 | phiBins->clear(); |
---|
109 | for(itPhiBin = itEtaBin->second.begin(); itPhiBin != itEtaBin->second.end(); ++itPhiBin) |
---|
110 | { |
---|
111 | phiBins->push_back(*itPhiBin); |
---|
112 | } |
---|
113 | } |
---|
114 | |
---|
115 | // read energy fractions for different particles |
---|
116 | param = GetParam("EnergyFraction"); |
---|
117 | size = param.GetSize(); |
---|
118 | |
---|
119 | // set default energy fractions values |
---|
120 | fFractionMap.clear(); |
---|
121 | fFractionMap[0] = make_pair(0.0, 1.0); |
---|
122 | |
---|
123 | for(i = 0; i < size/2; ++i) |
---|
124 | { |
---|
125 | paramFractions = param[i*2 + 1]; |
---|
126 | sizeFractions = paramFractions.GetSize(); |
---|
127 | |
---|
128 | ecalFraction = paramFractions[0].GetDouble(); |
---|
129 | hcalFraction = paramFractions[1].GetDouble(); |
---|
130 | |
---|
131 | fFractionMap[param[i*2].GetInt()] = make_pair(ecalFraction, hcalFraction); |
---|
132 | } |
---|
133 | /* |
---|
134 | TFractionMap::iterator itFractionMap; |
---|
135 | for(itFractionMap = fFractionMap.begin(); itFractionMap != fFractionMap.end(); ++itFractionMap) |
---|
136 | { |
---|
137 | cout << itFractionMap->first << " " << itFractionMap->second.first << " " << itFractionMap->second.second << endl; |
---|
138 | } |
---|
139 | */ |
---|
140 | // read resolution formulas |
---|
141 | fECalResolutionFormula->Compile(GetString("ECalResolutionFormula", "0")); |
---|
142 | fHCalResolutionFormula->Compile(GetString("HCalResolutionFormula", "0")); |
---|
143 | |
---|
144 | // import array with output from other modules |
---|
145 | fParticleInputArray = ImportArray(GetString("ParticleInputArray", "ParticlePropagator/particles")); |
---|
146 | fItParticleInputArray = fParticleInputArray->MakeIterator(); |
---|
147 | |
---|
148 | fTrackInputArray = ImportArray(GetString("TrackInputArray", "ParticlePropagator/tracks")); |
---|
149 | fItTrackInputArray = fTrackInputArray->MakeIterator(); |
---|
150 | |
---|
151 | // create output arrays |
---|
152 | fTowerOutputArray = ExportArray(GetString("TowerOutputArray", "towers")); |
---|
153 | fPhotonOutputArray = ExportArray(GetString("PhotonOutputArray", "photons")); |
---|
154 | |
---|
155 | fEFlowTrackOutputArray = ExportArray(GetString("EFlowTrackOutputArray", "eflowTracks")); |
---|
156 | fEFlowTowerOutputArray = ExportArray(GetString("EFlowTowerOutputArray", "eflowTowers")); |
---|
157 | } |
---|
158 | |
---|
159 | //------------------------------------------------------------------------------ |
---|
160 | |
---|
161 | void Calorimeter::Finish() |
---|
162 | { |
---|
163 | vector< vector< Double_t>* >::iterator itPhiBin; |
---|
164 | if(fItParticleInputArray) delete fItParticleInputArray; |
---|
165 | if(fItTrackInputArray) delete fItTrackInputArray; |
---|
166 | for(itPhiBin = fPhiBins.begin(); itPhiBin != fPhiBins.end(); ++itPhiBin) |
---|
167 | { |
---|
168 | delete *itPhiBin; |
---|
169 | } |
---|
170 | } |
---|
171 | |
---|
172 | //------------------------------------------------------------------------------ |
---|
173 | |
---|
174 | void Calorimeter::Process() |
---|
175 | { |
---|
176 | Candidate *particle, *track; |
---|
177 | TLorentzVector position, momentum; |
---|
178 | Short_t etaBin, phiBin, flags; |
---|
179 | Int_t number; |
---|
180 | Long64_t towerHit, towerEtaPhi, hitEtaPhi; |
---|
181 | Double_t ecalFraction, hcalFraction; |
---|
182 | Double_t ecalEnergy, hcalEnergy; |
---|
183 | Int_t pdgCode; |
---|
184 | |
---|
185 | TFractionMap::iterator itFractionMap; |
---|
186 | |
---|
187 | vector< Double_t >::iterator itEtaBin; |
---|
188 | vector< Double_t >::iterator itPhiBin; |
---|
189 | vector< Double_t > *phiBins; |
---|
190 | |
---|
191 | vector< Long64_t >::iterator itTowerHits; |
---|
192 | |
---|
193 | DelphesFactory *factory = GetFactory(); |
---|
194 | fTowerHits.clear(); |
---|
195 | fECalFractions.clear(); |
---|
196 | fHCalFractions.clear(); |
---|
197 | |
---|
198 | // loop over all particles |
---|
199 | fItParticleInputArray->Reset(); |
---|
200 | number = -1; |
---|
201 | while((particle = static_cast<Candidate*>(fItParticleInputArray->Next()))) |
---|
202 | { |
---|
203 | const TLorentzVector &particlePosition = particle->Position; |
---|
204 | ++number; |
---|
205 | |
---|
206 | pdgCode = TMath::Abs(particle->PID); |
---|
207 | |
---|
208 | itFractionMap = fFractionMap.find(pdgCode); |
---|
209 | if(itFractionMap == fFractionMap.end()) |
---|
210 | { |
---|
211 | itFractionMap = fFractionMap.find(0); |
---|
212 | } |
---|
213 | |
---|
214 | ecalFraction = itFractionMap->second.first; |
---|
215 | hcalFraction = itFractionMap->second.second; |
---|
216 | |
---|
217 | fECalFractions.push_back(ecalFraction); |
---|
218 | fHCalFractions.push_back(hcalFraction); |
---|
219 | |
---|
220 | if(ecalFraction < 1.0E-9 && hcalFraction < 1.0E-9) continue; |
---|
221 | |
---|
222 | // find eta bin [1, fEtaBins.size - 1] |
---|
223 | itEtaBin = lower_bound(fEtaBins.begin(), fEtaBins.end(), particlePosition.Eta()); |
---|
224 | if(itEtaBin == fEtaBins.begin() || itEtaBin == fEtaBins.end()) continue; |
---|
225 | etaBin = distance(fEtaBins.begin(), itEtaBin); |
---|
226 | |
---|
227 | // phi bins for given eta bin |
---|
228 | phiBins = fPhiBins[etaBin]; |
---|
229 | |
---|
230 | // find phi bin [1, phiBins.size - 1] |
---|
231 | itPhiBin = lower_bound(phiBins->begin(), phiBins->end(), particlePosition.Phi()); |
---|
232 | if(itPhiBin == phiBins->begin() || itPhiBin == phiBins->end()) continue; |
---|
233 | phiBin = distance(phiBins->begin(), itPhiBin); |
---|
234 | |
---|
235 | flags = (particle->Charge == 0); |
---|
236 | flags |= (pdgCode == 22) << 1; |
---|
237 | flags |= (pdgCode == 11) << 2; |
---|
238 | |
---|
239 | // make tower hit {16-bits for eta bin number, 16-bits for phi bin number, 8-bits for flags, 24-bits for particle number} |
---|
240 | towerHit = (Long64_t(etaBin) << 48) | (Long64_t(phiBin) << 32) | (Long64_t(flags) << 24) | Long64_t(number); |
---|
241 | |
---|
242 | fTowerHits.push_back(towerHit); |
---|
243 | } |
---|
244 | |
---|
245 | // loop over all tracks |
---|
246 | fItTrackInputArray->Reset(); |
---|
247 | number = -1; |
---|
248 | while((track = static_cast<Candidate*>(fItTrackInputArray->Next()))) |
---|
249 | { |
---|
250 | const TLorentzVector &trackPosition = track->Position; |
---|
251 | ++number; |
---|
252 | |
---|
253 | // find eta bin [1, fEtaBins.size - 1] |
---|
254 | itEtaBin = lower_bound(fEtaBins.begin(), fEtaBins.end(), trackPosition.Eta()); |
---|
255 | if(itEtaBin == fEtaBins.begin() || itEtaBin == fEtaBins.end()) continue; |
---|
256 | etaBin = distance(fEtaBins.begin(), itEtaBin); |
---|
257 | |
---|
258 | // phi bins for given eta bin |
---|
259 | phiBins = fPhiBins[etaBin]; |
---|
260 | |
---|
261 | // find phi bin [1, phiBins.size - 1] |
---|
262 | itPhiBin = lower_bound(phiBins->begin(), phiBins->end(), trackPosition.Phi()); |
---|
263 | if(itPhiBin == phiBins->begin() || itPhiBin == phiBins->end()) continue; |
---|
264 | phiBin = distance(phiBins->begin(), itPhiBin); |
---|
265 | |
---|
266 | // make tower hit {16-bits for eta bin number, 16-bits for phi bin number, 8-bits for flags, 24-bits for track number} |
---|
267 | towerHit = (Long64_t(etaBin) << 48) | (Long64_t(phiBin) << 32) | (Long64_t(1) << 27) | Long64_t(number); |
---|
268 | |
---|
269 | fTowerHits.push_back(towerHit); |
---|
270 | } |
---|
271 | |
---|
272 | // all hits are sorted first by eta bin number, then by phi bin number, |
---|
273 | // then by flags and then by particle or track number |
---|
274 | sort(fTowerHits.begin(), fTowerHits.end()); |
---|
275 | |
---|
276 | // loop over all hits |
---|
277 | towerEtaPhi = 0; |
---|
278 | fTower = 0; |
---|
279 | for(itTowerHits = fTowerHits.begin(); itTowerHits != fTowerHits.end(); ++itTowerHits) |
---|
280 | { |
---|
281 | towerHit = (*itTowerHits); |
---|
282 | flags = (towerHit >> 24) & 0x00000000000000FFLL; |
---|
283 | number = (towerHit) & 0x0000000000FFFFFFLL; |
---|
284 | hitEtaPhi = towerHit >> 32; |
---|
285 | |
---|
286 | if(towerEtaPhi != hitEtaPhi) |
---|
287 | { |
---|
288 | // switch to next tower |
---|
289 | towerEtaPhi = hitEtaPhi; |
---|
290 | |
---|
291 | // finalize previous tower |
---|
292 | FinalizeTower(); |
---|
293 | |
---|
294 | // create new tower |
---|
295 | fTower = factory->NewCandidate(); |
---|
296 | |
---|
297 | phiBin = (towerHit >> 32) & 0x000000000000FFFFLL; |
---|
298 | etaBin = (towerHit >> 48) & 0x000000000000FFFFLL; |
---|
299 | |
---|
300 | // phi bins for given eta bin |
---|
301 | phiBins = fPhiBins[etaBin]; |
---|
302 | |
---|
303 | // calculate eta and phi of the tower's center |
---|
304 | fTowerEta = 0.5*(fEtaBins[etaBin - 1] + fEtaBins[etaBin]); |
---|
305 | fTowerPhi = 0.5*((*phiBins)[phiBin - 1] + (*phiBins)[phiBin]); |
---|
306 | |
---|
307 | fTowerEdges[0] = fEtaBins[etaBin - 1]; |
---|
308 | fTowerEdges[1] = fEtaBins[etaBin]; |
---|
309 | fTowerEdges[2] = (*phiBins)[phiBin - 1]; |
---|
310 | fTowerEdges[3] = (*phiBins)[phiBin]; |
---|
311 | |
---|
312 | fTowerECalEnergy = 0.0; |
---|
313 | fTowerHCalEnergy = 0.0; |
---|
314 | |
---|
315 | fTowerECalNeutralEnergy = 0.0; |
---|
316 | fTowerHCalNeutralEnergy = 0.0; |
---|
317 | |
---|
318 | fTowerNeutralHits = 0; |
---|
319 | fTowerPhotonHits = 0; |
---|
320 | fTowerElectronHits = 0; |
---|
321 | fTowerTrackHits = 0; |
---|
322 | fTowerAllHits = 0; |
---|
323 | |
---|
324 | fTowerTrackArray->Clear(); |
---|
325 | fTowerPhotonArray->Clear(); |
---|
326 | } |
---|
327 | |
---|
328 | // check for track hits |
---|
329 | if(flags & 8) |
---|
330 | { |
---|
331 | ++fTowerTrackHits; |
---|
332 | track = static_cast<Candidate*>(fTrackInputArray->At(number)); |
---|
333 | fTowerTrackArray->Add(track); |
---|
334 | continue; |
---|
335 | } |
---|
336 | |
---|
337 | particle = static_cast<Candidate*>(fParticleInputArray->At(number)); |
---|
338 | momentum = particle->Momentum; |
---|
339 | |
---|
340 | // fill current tower |
---|
341 | ecalEnergy = momentum.E() * fECalFractions[number]; |
---|
342 | hcalEnergy = momentum.E() * fHCalFractions[number]; |
---|
343 | |
---|
344 | fTowerECalEnergy += ecalEnergy; |
---|
345 | fTowerHCalEnergy += hcalEnergy; |
---|
346 | |
---|
347 | ++fTowerAllHits; |
---|
348 | fTower->AddCandidate(particle); |
---|
349 | |
---|
350 | // check for neutral hits in current tower |
---|
351 | if(flags & 1) ++fTowerNeutralHits; |
---|
352 | |
---|
353 | // check for photon hits in current tower |
---|
354 | if(flags & 2) |
---|
355 | { |
---|
356 | ++fTowerPhotonHits; |
---|
357 | fTowerPhotonArray->Add(particle); |
---|
358 | } |
---|
359 | |
---|
360 | // check for electron hits in current tower |
---|
361 | if(flags & 4) ++fTowerElectronHits; |
---|
362 | } |
---|
363 | |
---|
364 | // finalize last tower |
---|
365 | FinalizeTower(); |
---|
366 | } |
---|
367 | |
---|
368 | //------------------------------------------------------------------------------ |
---|
369 | |
---|
370 | void Calorimeter::FinalizeTower() |
---|
371 | { |
---|
372 | Candidate *particle, *track, *tower; |
---|
373 | Double_t energy, pt, eta, phi; |
---|
374 | Double_t ecalEnergy, hcalEnergy; |
---|
375 | |
---|
376 | if(!fTower) return; |
---|
377 | |
---|
378 | ecalEnergy = gRandom->Gaus(fTowerECalEnergy, fECalResolutionFormula->Eval(0.0, fTowerEta, 0.0, fTowerECalEnergy)); |
---|
379 | if(ecalEnergy < 0.0) ecalEnergy = 0.0; |
---|
380 | |
---|
381 | hcalEnergy = gRandom->Gaus(fTowerHCalEnergy, fHCalResolutionFormula->Eval(0.0, fTowerEta, 0.0, fTowerHCalEnergy)); |
---|
382 | if(hcalEnergy < 0.0) hcalEnergy = 0.0; |
---|
383 | |
---|
384 | energy = ecalEnergy + hcalEnergy; |
---|
385 | |
---|
386 | // eta = fTowerEta; |
---|
387 | // phi = fTowerPhi; |
---|
388 | |
---|
389 | eta = gRandom->Uniform(fTowerEdges[0], fTowerEdges[1]); |
---|
390 | phi = gRandom->Uniform(fTowerEdges[2], fTowerEdges[3]); |
---|
391 | |
---|
392 | pt = energy / TMath::CosH(eta); |
---|
393 | |
---|
394 | fTower->Position.SetPtEtaPhiE(1.0, eta, phi, 0.0); |
---|
395 | fTower->Momentum.SetPtEtaPhiE(pt, eta, phi, energy); |
---|
396 | fTower->Eem = ecalEnergy; |
---|
397 | fTower->Ehad = hcalEnergy; |
---|
398 | |
---|
399 | fTower->Edges[0] = fTowerEdges[0]; |
---|
400 | fTower->Edges[1] = fTowerEdges[1]; |
---|
401 | fTower->Edges[2] = fTowerEdges[2]; |
---|
402 | fTower->Edges[3] = fTowerEdges[3]; |
---|
403 | |
---|
404 | // fill calorimeter towers and photon candidates |
---|
405 | if(energy > 0.0) |
---|
406 | { |
---|
407 | if((fTowerPhotonHits > 0 || fTowerElectronHits > 0) && |
---|
408 | fTowerTrackHits == 0) |
---|
409 | { |
---|
410 | fPhotonOutputArray->Add(fTower); |
---|
411 | } |
---|
412 | |
---|
413 | fTowerOutputArray->Add(fTower); |
---|
414 | } |
---|
415 | |
---|
416 | // fill energy flow candidates |
---|
417 | if(fTowerTrackHits == fTowerAllHits) |
---|
418 | { |
---|
419 | fItTowerTrackArray->Reset(); |
---|
420 | while((track = static_cast<Candidate*>(fItTowerTrackArray->Next()))) |
---|
421 | { |
---|
422 | fEFlowTrackOutputArray->Add(track); |
---|
423 | } |
---|
424 | } |
---|
425 | else if(fTowerTrackHits > 0 && |
---|
426 | fTowerElectronHits == 0 && |
---|
427 | fTowerPhotonHits + fTowerTrackHits == fTowerAllHits) |
---|
428 | { |
---|
429 | fItTowerTrackArray->Reset(); |
---|
430 | while((track = static_cast<Candidate*>(fItTowerTrackArray->Next()))) |
---|
431 | { |
---|
432 | fEFlowTrackOutputArray->Add(track); |
---|
433 | } |
---|
434 | |
---|
435 | if(ecalEnergy > 0.0) |
---|
436 | { |
---|
437 | DelphesFactory *factory = GetFactory(); |
---|
438 | |
---|
439 | // create new tower |
---|
440 | tower = factory->NewCandidate(); |
---|
441 | |
---|
442 | fItTowerPhotonArray->Reset(); |
---|
443 | while((particle = static_cast<Candidate*>(fItTowerPhotonArray->Next()))) |
---|
444 | { |
---|
445 | tower->AddCandidate(particle); |
---|
446 | } |
---|
447 | |
---|
448 | pt = ecalEnergy / TMath::CosH(eta); |
---|
449 | |
---|
450 | tower->Position.SetPtEtaPhiE(1.0, eta, phi, 0.0); |
---|
451 | tower->Momentum.SetPtEtaPhiE(pt, eta, phi, ecalEnergy); |
---|
452 | tower->Eem = ecalEnergy; |
---|
453 | tower->Ehad = 0.0; |
---|
454 | |
---|
455 | tower->Edges[0] = fTowerEdges[0]; |
---|
456 | tower->Edges[1] = fTowerEdges[1]; |
---|
457 | tower->Edges[2] = fTowerEdges[2]; |
---|
458 | tower->Edges[3] = fTowerEdges[3]; |
---|
459 | |
---|
460 | fEFlowTowerOutputArray->Add(tower); |
---|
461 | } |
---|
462 | } |
---|
463 | else if(energy > 0.0) |
---|
464 | { |
---|
465 | fEFlowTowerOutputArray->Add(fTower); |
---|
466 | } |
---|
467 | } |
---|
468 | |
---|
469 | //------------------------------------------------------------------------------ |
---|