Fork me on GitHub

source: svn/trunk/modules/Calorimeter.cc@ 1245

Last change on this file since 1245 was 1242, checked in by Pavel Demin, 11 years ago

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