Fork me on GitHub

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

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

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