Fork me on GitHub

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

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

new energy flow that compiles
  • Property svn:keywords set to Id Revision Date
File size: 16.3 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-22 20:00:23 +0000 (Mon, 22 Jul 2013) $
8 * $Revision: 1235 $
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 fTower->AddCandidate(particle);
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
422 // finalize last tower
423 FinalizeTower();
424}
425
426//------------------------------------------------------------------------------
427
428void Calorimeter::FinalizeTower()
429{
430 Candidate *particle, *track, *tower;
431 Double_t energy, pt, eta, phi;
432 Double_t ecalEnergy, hcalEnergy;
433 TIterator *iterator;
434
435 if(!fTower) return;
436
437// ecalEnergy = gRandom->Gaus(fTowerECalEnergy, fECalResolutionFormula->Eval(0.0, fTowerEta, 0.0, fTowerECalEnergy));
438// if(ecalEnergy < 0.0) ecalEnergy = 0.0;
439
440 ecalEnergy = LogNormal(fTowerECalEnergy, fECalResolutionFormula->Eval(0.0, fTowerEta, 0.0, fTowerECalEnergy));
441
442// hcalEnergy = gRandom->Gaus(fTowerHCalEnergy, fHCalResolutionFormula->Eval(0.0, fTowerEta, 0.0, fTowerHCalEnergy));
443// if(hcalEnergy < 0.0) hcalEnergy = 0.0;
444
445 hcalEnergy = LogNormal(fTowerHCalEnergy, fHCalResolutionFormula->Eval(0.0, fTowerEta, 0.0, fTowerHCalEnergy));
446
447 energy = ecalEnergy + hcalEnergy;
448
449// eta = fTowerEta;
450// phi = fTowerPhi;
451
452 eta = gRandom->Uniform(fTowerEdges[0], fTowerEdges[1]);
453 phi = gRandom->Uniform(fTowerEdges[2], fTowerEdges[3]);
454
455 pt = energy / TMath::CosH(eta);
456
457 fTower->Position.SetPtEtaPhiE(1.0, eta, phi, 0.0);
458 fTower->Momentum.SetPtEtaPhiE(pt, eta, phi, energy);
459 fTower->Eem = ecalEnergy;
460 fTower->Ehad = hcalEnergy;
461
462 fTower->Edges[0] = fTowerEdges[0];
463 fTower->Edges[1] = fTowerEdges[1];
464 fTower->Edges[2] = fTowerEdges[2];
465 fTower->Edges[3] = fTowerEdges[3];
466
467 // fill calorimeter towers and photon candidates
468 if(energy > 0.0)
469 {
470 if(fTowerPhotonHits > 0 && fTowerTrackAllHits == 0)
471 {
472 fPhotonOutputArray->Add(fTower);
473 }
474
475 fTowerOutputArray->Add(fTower);
476 }
477
478 // fill energy flow candidates
479 if(fTowerTrackAllHits == fTowerAllHits)
480 {
481 fItTowerTrackArray->Reset();
482 while((track = static_cast<Candidate*>(fItTowerTrackArray->Next())))
483 {
484 fEFlowTrackOutputArray->Add(track);
485 }
486 }
487 else if(fTowerTrackAllHits > 0 &&
488 fTowerECalHits + fTowerHCalHits == fTowerAllHits &&
489 (fTowerECalHits == fTowerECalTrackHits || fTowerHCalHits == fTowerHCalTrackHits))
490 {
491 if(fTowerECalHits == fTowerECalTrackHits)
492 {
493 fItTowerECalTrackArray->Reset();
494 while((track = static_cast<Candidate*>(fItTowerECalTrackArray->Next())))
495 {
496 fEFlowTrackOutputArray->Add(track);
497 }
498 energy = hcalEnergy;
499 iterator = fItTowerHCalArray;
500 }
501
502 if(fTowerHCalHits == fTowerHCalTrackHits)
503 {
504 fItTowerHCalTrackArray->Reset();
505 while((track = static_cast<Candidate*>(fItTowerHCalTrackArray->Next())))
506 {
507 fEFlowTrackOutputArray->Add(track);
508 }
509 energy = ecalEnergy;
510 iterator = fItTowerECalArray;
511 }
512
513 if(energy > 0.0 && iterator)
514 {
515 DelphesFactory *factory = GetFactory();
516
517 // create new tower
518 tower = factory->NewCandidate();
519
520 iterator->Reset();
521 while((particle = static_cast<Candidate*>(iterator->Next())))
522 {
523 tower->AddCandidate(particle);
524 }
525
526 pt = energy / TMath::CosH(eta);
527
528 tower->Position.SetPtEtaPhiE(1.0, eta, phi, 0.0);
529 tower->Momentum.SetPtEtaPhiE(pt, eta, phi, energy);
530 tower->Eem = energy;
531 tower->Ehad = 0.0;
532
533 tower->Edges[0] = fTowerEdges[0];
534 tower->Edges[1] = fTowerEdges[1];
535 tower->Edges[2] = fTowerEdges[2];
536 tower->Edges[3] = fTowerEdges[3];
537
538 fEFlowTowerOutputArray->Add(tower);
539 }
540 }
541 else if(energy > 0.0)
542 {
543 fEFlowTowerOutputArray->Add(fTower);
544 }
545}
546
547//------------------------------------------------------------------------------
548
549Double_t Calorimeter::LogNormal(Double_t mean, Double_t sigma)
550{
551 Double_t a, b;
552
553 if(mean > 0.0)
554 {
555 b = TMath::Sqrt(TMath::Log((1.0 + (sigma*sigma)/(mean*mean))));
556 a = TMath::Log(mean) - 0.5*b*b;
557
558 return TMath::Exp(a + b*gRandom->Gaus(0, 1));
559 }
560 else
561 {
562 return 0.0;
563 }
564}
Note: See TracBrowser for help on using the repository browser.