Fork me on GitHub

source: git/modules/OldCalorimeter.cc@ 0852ba95

Last change on this file since 0852ba95 was 341014c, checked in by Pavel Demin <pavel-demin@…>, 6 years ago

apply .clang-format to all .h, .cc and .cpp files
  • Property mode set to 100644
File size: 17.2 KB
Line 
1
2/** \class OldCalorimeter
3 *
4 * Fills calorimeter towers, performs calorimeter resolution smearing,
5 * preselects towers hit by photons and creates energy flow objects.
6 *
7 * $Date$
8 * $Revision$
9 *
10 *
11 * \author P. Demin - UCL, Louvain-la-Neuve
12 *
13 */
14
15#include "modules/OldCalorimeter.h"
16
17#include "classes/DelphesClasses.h"
18#include "classes/DelphesFactory.h"
19#include "classes/DelphesFormula.h"
20
21#include "ExRootAnalysis/ExRootClassifier.h"
22#include "ExRootAnalysis/ExRootFilter.h"
23#include "ExRootAnalysis/ExRootResult.h"
24
25#include "TDatabasePDG.h"
26#include "TFormula.h"
27#include "TLorentzVector.h"
28#include "TMath.h"
29#include "TObjArray.h"
30#include "TRandom3.h"
31#include "TString.h"
32
33#include <algorithm>
34#include <iostream>
35#include <sstream>
36#include <stdexcept>
37
38using namespace std;
39
40//------------------------------------------------------------------------------
41
42OldCalorimeter::OldCalorimeter() :
43 fECalResolutionFormula(0), fHCalResolutionFormula(0),
44 fItParticleInputArray(0), fItTrackInputArray(0),
45 fTowerECalArray(0), fItTowerECalArray(0),
46 fTowerHCalArray(0), fItTowerHCalArray(0),
47 fTowerTrackArray(0), fItTowerTrackArray(0),
48 fTowerECalTrackArray(0), fItTowerECalTrackArray(0),
49 fTowerHCalTrackArray(0), fItTowerHCalTrackArray(0)
50{
51 fECalResolutionFormula = new DelphesFormula;
52 fHCalResolutionFormula = new DelphesFormula;
53
54 fTowerECalArray = new TObjArray;
55 fItTowerECalArray = fTowerECalArray->MakeIterator();
56 fTowerHCalArray = new TObjArray;
57 fItTowerHCalArray = fTowerHCalArray->MakeIterator();
58
59 fTowerTrackArray = new TObjArray;
60 fItTowerTrackArray = fTowerTrackArray->MakeIterator();
61 fTowerECalTrackArray = new TObjArray;
62 fItTowerECalTrackArray = fTowerECalTrackArray->MakeIterator();
63 fTowerHCalTrackArray = new TObjArray;
64 fItTowerHCalTrackArray = fTowerHCalTrackArray->MakeIterator();
65}
66
67//------------------------------------------------------------------------------
68
69OldCalorimeter::~OldCalorimeter()
70{
71 if(fECalResolutionFormula) delete fECalResolutionFormula;
72 if(fHCalResolutionFormula) delete fHCalResolutionFormula;
73
74 if(fTowerECalArray) delete fTowerECalArray;
75 if(fItTowerECalArray) delete fItTowerECalArray;
76 if(fTowerHCalArray) delete fTowerHCalArray;
77 if(fItTowerHCalArray) delete fItTowerHCalArray;
78
79 if(fTowerTrackArray) delete fTowerTrackArray;
80 if(fItTowerTrackArray) delete fItTowerTrackArray;
81 if(fTowerECalTrackArray) delete fTowerECalTrackArray;
82 if(fItTowerECalTrackArray) delete fItTowerECalTrackArray;
83 if(fTowerHCalTrackArray) delete fTowerHCalTrackArray;
84 if(fItTowerHCalTrackArray) delete fItTowerHCalTrackArray;
85}
86
87//------------------------------------------------------------------------------
88
89void OldCalorimeter::Init()
90{
91 ExRootConfParam param, paramEtaBins, paramPhiBins, paramFractions;
92 Long_t i, j, k, size, sizeEtaBins, sizePhiBins, sizeFractions;
93 Double_t ecalFraction, hcalFraction;
94 TBinMap::iterator itEtaBin;
95 set<Double_t>::iterator itPhiBin;
96 vector<Double_t> *phiBins;
97
98 // read eta and phi bins
99 param = GetParam("EtaPhiBins");
100 size = param.GetSize();
101 fBinMap.clear();
102 fEtaBins.clear();
103 fPhiBins.clear();
104 for(i = 0; i < size / 2; ++i)
105 {
106 paramEtaBins = param[i * 2];
107 sizeEtaBins = paramEtaBins.GetSize();
108 paramPhiBins = param[i * 2 + 1];
109 sizePhiBins = paramPhiBins.GetSize();
110
111 for(j = 0; j < sizeEtaBins; ++j)
112 {
113 for(k = 0; k < sizePhiBins; ++k)
114 {
115 fBinMap[paramEtaBins[j].GetDouble()].insert(paramPhiBins[k].GetDouble());
116 }
117 }
118 }
119
120 // for better performance we transform map of sets to parallel vectors:
121 // vector< double > and vector< vector< double >* >
122 for(itEtaBin = fBinMap.begin(); itEtaBin != fBinMap.end(); ++itEtaBin)
123 {
124 fEtaBins.push_back(itEtaBin->first);
125 phiBins = new vector<double>(itEtaBin->second.size());
126 fPhiBins.push_back(phiBins);
127 phiBins->clear();
128 for(itPhiBin = itEtaBin->second.begin(); itPhiBin != itEtaBin->second.end(); ++itPhiBin)
129 {
130 phiBins->push_back(*itPhiBin);
131 }
132 }
133
134 // read energy fractions for different particles
135 param = GetParam("EnergyFraction");
136 size = param.GetSize();
137
138 // set default energy fractions values
139 fFractionMap.clear();
140 fFractionMap[0] = make_pair(0.0, 1.0);
141
142 for(i = 0; i < size / 2; ++i)
143 {
144 paramFractions = param[i * 2 + 1];
145 sizeFractions = paramFractions.GetSize();
146
147 ecalFraction = paramFractions[0].GetDouble();
148 hcalFraction = paramFractions[1].GetDouble();
149
150 fFractionMap[param[i * 2].GetInt()] = make_pair(ecalFraction, hcalFraction);
151 }
152 /*
153 TFractionMap::iterator itFractionMap;
154 for(itFractionMap = fFractionMap.begin(); itFractionMap != fFractionMap.end(); ++itFractionMap)
155 {
156 cout << itFractionMap->first << " " << itFractionMap->second.first << " " << itFractionMap->second.second << endl;
157 }
158*/
159 // read resolution formulas
160 fECalResolutionFormula->Compile(GetString("ECalResolutionFormula", "0"));
161 fHCalResolutionFormula->Compile(GetString("HCalResolutionFormula", "0"));
162
163 // import array with output from other modules
164 fParticleInputArray = ImportArray(GetString("ParticleInputArray", "ParticlePropagator/particles"));
165 fItParticleInputArray = fParticleInputArray->MakeIterator();
166
167 fTrackInputArray = ImportArray(GetString("TrackInputArray", "ParticlePropagator/tracks"));
168 fItTrackInputArray = fTrackInputArray->MakeIterator();
169
170 // create output arrays
171 fTowerOutputArray = ExportArray(GetString("TowerOutputArray", "towers"));
172 fPhotonOutputArray = ExportArray(GetString("PhotonOutputArray", "photons"));
173
174 fEFlowTrackOutputArray = ExportArray(GetString("EFlowTrackOutputArray", "eflowTracks"));
175 fEFlowTowerOutputArray = ExportArray(GetString("EFlowTowerOutputArray", "eflowTowers"));
176}
177
178//------------------------------------------------------------------------------
179
180void OldCalorimeter::Finish()
181{
182 vector<vector<Double_t> *>::iterator itPhiBin;
183 if(fItParticleInputArray) delete fItParticleInputArray;
184 if(fItTrackInputArray) delete fItTrackInputArray;
185 for(itPhiBin = fPhiBins.begin(); itPhiBin != fPhiBins.end(); ++itPhiBin)
186 {
187 delete *itPhiBin;
188 }
189}
190
191//------------------------------------------------------------------------------
192
193void OldCalorimeter::Process()
194{
195 Candidate *particle, *track;
196 TLorentzVector position, momentum;
197 Short_t etaBin, phiBin, flags;
198 Int_t number;
199 Long64_t towerHit, towerEtaPhi, hitEtaPhi;
200 Double_t ecalFraction, hcalFraction;
201 Double_t ecalEnergy, hcalEnergy;
202 Int_t pdgCode;
203
204 TFractionMap::iterator itFractionMap;
205
206 vector<Double_t>::iterator itEtaBin;
207 vector<Double_t>::iterator itPhiBin;
208 vector<Double_t> *phiBins;
209
210 vector<Long64_t>::iterator itTowerHits;
211
212 DelphesFactory *factory = GetFactory();
213 fTowerHits.clear();
214 fECalFractions.clear();
215 fHCalFractions.clear();
216
217 // loop over all particles
218 fItParticleInputArray->Reset();
219 number = -1;
220 while((particle = static_cast<Candidate *>(fItParticleInputArray->Next())))
221 {
222 const TLorentzVector &particlePosition = particle->Position;
223 ++number;
224
225 pdgCode = TMath::Abs(particle->PID);
226
227 itFractionMap = fFractionMap.find(pdgCode);
228 if(itFractionMap == fFractionMap.end())
229 {
230 itFractionMap = fFractionMap.find(0);
231 }
232
233 ecalFraction = itFractionMap->second.first;
234 hcalFraction = itFractionMap->second.second;
235
236 fECalFractions.push_back(ecalFraction);
237 fHCalFractions.push_back(hcalFraction);
238
239 if(ecalFraction < 1.0E-9 && hcalFraction < 1.0E-9) continue;
240
241 // find eta bin [1, fEtaBins.size - 1]
242 itEtaBin = lower_bound(fEtaBins.begin(), fEtaBins.end(), particlePosition.Eta());
243 if(itEtaBin == fEtaBins.begin() || itEtaBin == fEtaBins.end()) continue;
244 etaBin = distance(fEtaBins.begin(), itEtaBin);
245
246 // phi bins for given eta bin
247 phiBins = fPhiBins[etaBin];
248
249 // find phi bin [1, phiBins.size - 1]
250 itPhiBin = lower_bound(phiBins->begin(), phiBins->end(), particlePosition.Phi());
251 if(itPhiBin == phiBins->begin() || itPhiBin == phiBins->end()) continue;
252 phiBin = distance(phiBins->begin(), itPhiBin);
253
254 flags = 0;
255 flags |= (ecalFraction >= 1.0E-9) << 1;
256 flags |= (hcalFraction >= 1.0E-9) << 2;
257 flags |= (pdgCode == 11 || pdgCode == 22) << 3;
258
259 // make tower hit {16-bits for eta bin number, 16-bits for phi bin number, 8-bits for flags, 24-bits for particle number}
260 towerHit = (Long64_t(etaBin) << 48) | (Long64_t(phiBin) << 32) | (Long64_t(flags) << 24) | Long64_t(number);
261
262 fTowerHits.push_back(towerHit);
263 }
264
265 // loop over all tracks
266 fItTrackInputArray->Reset();
267 number = -1;
268 while((track = static_cast<Candidate *>(fItTrackInputArray->Next())))
269 {
270 const TLorentzVector &trackPosition = track->Position;
271 ++number;
272
273 pdgCode = TMath::Abs(track->PID);
274
275 itFractionMap = fFractionMap.find(pdgCode);
276 if(itFractionMap == fFractionMap.end())
277 {
278 itFractionMap = fFractionMap.find(0);
279 }
280
281 ecalFraction = itFractionMap->second.first;
282 hcalFraction = itFractionMap->second.second;
283
284 // find eta bin [1, fEtaBins.size - 1]
285 itEtaBin = lower_bound(fEtaBins.begin(), fEtaBins.end(), trackPosition.Eta());
286 if(itEtaBin == fEtaBins.begin() || itEtaBin == fEtaBins.end()) continue;
287 etaBin = distance(fEtaBins.begin(), itEtaBin);
288
289 // phi bins for given eta bin
290 phiBins = fPhiBins[etaBin];
291
292 // find phi bin [1, phiBins.size - 1]
293 itPhiBin = lower_bound(phiBins->begin(), phiBins->end(), trackPosition.Phi());
294 if(itPhiBin == phiBins->begin() || itPhiBin == phiBins->end()) continue;
295 phiBin = distance(phiBins->begin(), itPhiBin);
296
297 flags = 1;
298 flags |= (ecalFraction >= 1.0E-9) << 1;
299 flags |= (hcalFraction >= 1.0E-9) << 2;
300
301 // make tower hit {16-bits for eta bin number, 16-bits for phi bin number, 8-bits for flags, 24-bits for track number}
302 towerHit = (Long64_t(etaBin) << 48) | (Long64_t(phiBin) << 32) | (Long64_t(flags) << 24) | Long64_t(number);
303
304 fTowerHits.push_back(towerHit);
305 }
306
307 // all hits are sorted first by eta bin number, then by phi bin number,
308 // then by flags and then by particle or track number
309 sort(fTowerHits.begin(), fTowerHits.end());
310
311 // loop over all hits
312 towerEtaPhi = 0;
313 fTower = 0;
314 for(itTowerHits = fTowerHits.begin(); itTowerHits != fTowerHits.end(); ++itTowerHits)
315 {
316 towerHit = (*itTowerHits);
317 flags = (towerHit >> 24) & 0x00000000000000FFLL;
318 number = (towerHit)&0x0000000000FFFFFFLL;
319 hitEtaPhi = towerHit >> 32;
320
321 if(towerEtaPhi != hitEtaPhi)
322 {
323 // switch to next tower
324 towerEtaPhi = hitEtaPhi;
325
326 // finalize previous tower
327 FinalizeTower();
328
329 // create new tower
330 fTower = factory->NewCandidate();
331
332 phiBin = (towerHit >> 32) & 0x000000000000FFFFLL;
333 etaBin = (towerHit >> 48) & 0x000000000000FFFFLL;
334
335 // phi bins for given eta bin
336 phiBins = fPhiBins[etaBin];
337
338 // calculate eta and phi of the tower's center
339 fTowerEta = 0.5 * (fEtaBins[etaBin - 1] + fEtaBins[etaBin]);
340 fTowerPhi = 0.5 * ((*phiBins)[phiBin - 1] + (*phiBins)[phiBin]);
341
342 fTowerEdges[0] = fEtaBins[etaBin - 1];
343 fTowerEdges[1] = fEtaBins[etaBin];
344 fTowerEdges[2] = (*phiBins)[phiBin - 1];
345 fTowerEdges[3] = (*phiBins)[phiBin];
346
347 fTowerECalEnergy = 0.0;
348 fTowerHCalEnergy = 0.0;
349
350 fTowerPhotonHits = 0;
351
352 fTowerAllHits = 0;
353 fTowerECalHits = 0;
354 fTowerHCalHits = 0;
355
356 fTowerTrackAllHits = 0;
357 fTowerECalTrackHits = 0;
358 fTowerHCalTrackHits = 0;
359
360 fTowerECalArray->Clear();
361 fTowerHCalArray->Clear();
362
363 fTowerTrackArray->Clear();
364 fTowerECalTrackArray->Clear();
365 fTowerHCalTrackArray->Clear();
366 }
367
368 // check for track hits
369 if(flags & 1)
370 {
371 track = static_cast<Candidate *>(fTrackInputArray->At(number));
372
373 ++fTowerTrackAllHits;
374 fTowerTrackArray->Add(track);
375
376 // check for track ECAL hits
377 if(flags & 2)
378 {
379 ++fTowerECalTrackHits;
380 fTowerECalTrackArray->Add(track);
381 }
382
383 // check for track HCAL hits
384 if(flags & 4)
385 {
386 ++fTowerHCalTrackHits;
387 fTowerHCalTrackArray->Add(track);
388 }
389 continue;
390 }
391
392 ++fTowerAllHits;
393
394 // check for ECAL hits
395 if(flags & 2)
396 {
397 ++fTowerECalHits;
398 fTowerECalArray->Add(particle);
399 }
400
401 // check for HCAL hits
402 if(flags & 4)
403 {
404 ++fTowerHCalHits;
405 fTowerHCalArray->Add(particle);
406 }
407
408 // check for photon and electron hits in current tower
409 if(flags & 8) ++fTowerPhotonHits;
410
411 particle = static_cast<Candidate *>(fParticleInputArray->At(number));
412 momentum = particle->Momentum;
413
414 // fill current tower
415 ecalEnergy = momentum.E() * fECalFractions[number];
416 hcalEnergy = momentum.E() * fHCalFractions[number];
417
418 fTowerECalEnergy += ecalEnergy;
419 fTowerHCalEnergy += hcalEnergy;
420
421 fTower->AddCandidate(particle);
422 }
423
424 // finalize last tower
425 FinalizeTower();
426}
427
428//------------------------------------------------------------------------------
429
430void OldCalorimeter::FinalizeTower()
431{
432 Candidate *particle, *track, *tower;
433 Double_t energy, pt, eta, phi;
434 Double_t ecalEnergy, hcalEnergy;
435 TIterator *itTowerTrackArray;
436
437 if(!fTower) return;
438
439 // ecalEnergy = gRandom->Gaus(fTowerECalEnergy, fECalResolutionFormula->Eval(0.0, fTowerEta, 0.0, fTowerECalEnergy));
440 // if(ecalEnergy < 0.0) ecalEnergy = 0.0;
441
442 ecalEnergy = LogNormal(fTowerECalEnergy, fECalResolutionFormula->Eval(0.0, fTowerEta, 0.0, fTowerECalEnergy));
443
444 // hcalEnergy = gRandom->Gaus(fTowerHCalEnergy, fHCalResolutionFormula->Eval(0.0, fTowerEta, 0.0, fTowerHCalEnergy));
445 // if(hcalEnergy < 0.0) hcalEnergy = 0.0;
446
447 hcalEnergy = LogNormal(fTowerHCalEnergy, fHCalResolutionFormula->Eval(0.0, fTowerEta, 0.0, fTowerHCalEnergy));
448
449 energy = ecalEnergy + hcalEnergy;
450
451 // eta = fTowerEta;
452 // phi = fTowerPhi;
453
454 eta = gRandom->Uniform(fTowerEdges[0], fTowerEdges[1]);
455 phi = gRandom->Uniform(fTowerEdges[2], fTowerEdges[3]);
456
457 pt = energy / TMath::CosH(eta);
458
459 fTower->Position.SetPtEtaPhiE(1.0, eta, phi, 0.0);
460 fTower->Momentum.SetPtEtaPhiE(pt, eta, phi, energy);
461 fTower->Eem = ecalEnergy;
462 fTower->Ehad = hcalEnergy;
463
464 fTower->Edges[0] = fTowerEdges[0];
465 fTower->Edges[1] = fTowerEdges[1];
466 fTower->Edges[2] = fTowerEdges[2];
467 fTower->Edges[3] = fTowerEdges[3];
468
469 // fill calorimeter towers and photon candidates
470 if(energy > 0.0)
471 {
472 if(fTowerPhotonHits > 0 && fTowerTrackAllHits == 0)
473 {
474 fPhotonOutputArray->Add(fTower);
475 }
476
477 fTowerOutputArray->Add(fTower);
478 }
479
480 // fill energy flow candidates
481 if(fTowerTrackAllHits == fTowerAllHits)
482 {
483 fItTowerTrackArray->Reset();
484 while((track = static_cast<Candidate *>(fItTowerTrackArray->Next())))
485 {
486 fEFlowTrackOutputArray->Add(track);
487 }
488 }
489 else if(fTowerTrackAllHits > 0 && fTowerECalHits + fTowerHCalHits == fTowerAllHits)
490 {
491 if(fTowerECalHits == fTowerECalTrackHits && fTowerHCalHits == fTowerHCalTrackHits)
492 {
493 itTowerTrackArray = fItTowerTrackArray;
494 }
495 else if(fTowerECalHits == fTowerECalTrackHits)
496 {
497 itTowerTrackArray = fItTowerECalTrackArray;
498
499 if(hcalEnergy > 0.0)
500 {
501 DelphesFactory *factory = GetFactory();
502
503 // create new tower
504 tower = factory->NewCandidate();
505
506 fItTowerHCalArray->Reset();
507 while((particle = static_cast<Candidate *>(fItTowerHCalArray->Next())))
508 {
509 tower->AddCandidate(particle);
510 }
511
512 pt = hcalEnergy / TMath::CosH(eta);
513
514 tower->Position.SetPtEtaPhiE(1.0, eta, phi, 0.0);
515 tower->Momentum.SetPtEtaPhiE(pt, eta, phi, hcalEnergy);
516 tower->Eem = 0.0;
517 tower->Ehad = hcalEnergy;
518
519 tower->Edges[0] = fTowerEdges[0];
520 tower->Edges[1] = fTowerEdges[1];
521 tower->Edges[2] = fTowerEdges[2];
522 tower->Edges[3] = fTowerEdges[3];
523
524 fEFlowTowerOutputArray->Add(tower);
525 }
526 }
527 else if(fTowerHCalHits == fTowerHCalTrackHits)
528 {
529 itTowerTrackArray = fItTowerHCalTrackArray;
530
531 if(ecalEnergy > 0.0)
532 {
533 DelphesFactory *factory = GetFactory();
534
535 // create new tower
536 tower = factory->NewCandidate();
537
538 fItTowerECalArray->Reset();
539 while((particle = static_cast<Candidate *>(fItTowerECalArray->Next())))
540 {
541 tower->AddCandidate(particle);
542 }
543
544 pt = ecalEnergy / TMath::CosH(eta);
545
546 tower->Position.SetPtEtaPhiE(1.0, eta, phi, 0.0);
547 tower->Momentum.SetPtEtaPhiE(pt, eta, phi, ecalEnergy);
548 tower->Eem = ecalEnergy;
549 tower->Ehad = 0.0;
550
551 tower->Edges[0] = fTowerEdges[0];
552 tower->Edges[1] = fTowerEdges[1];
553 tower->Edges[2] = fTowerEdges[2];
554 tower->Edges[3] = fTowerEdges[3];
555
556 fEFlowTowerOutputArray->Add(tower);
557 }
558 }
559 else
560 {
561 itTowerTrackArray = 0;
562 fEFlowTowerOutputArray->Add(fTower);
563 }
564
565 if(itTowerTrackArray)
566 {
567 itTowerTrackArray->Reset();
568 while((track = static_cast<Candidate *>(itTowerTrackArray->Next())))
569 {
570 fEFlowTrackOutputArray->Add(track);
571 }
572 }
573 }
574 else if(energy > 0.0)
575 {
576 fEFlowTowerOutputArray->Add(fTower);
577 }
578}
579
580//------------------------------------------------------------------------------
581
582Double_t OldCalorimeter::LogNormal(Double_t mean, Double_t sigma)
583{
584 Double_t a, b;
585
586 if(mean > 0.0)
587 {
588 b = TMath::Sqrt(TMath::Log((1.0 + (sigma * sigma) / (mean * mean))));
589 a = TMath::Log(mean) - 0.5 * b * b;
590
591 return TMath::Exp(a + b * gRandom->Gaus(0, 1));
592 }
593 else
594 {
595 return 0.0;
596 }
597}
Note: See TracBrowser for help on using the repository browser.