Fork me on GitHub

source: git/modules/Calorimeter.cc@ 22dc7fd

ImprovedOutputFile Timing dual_readout llp
Last change on this file since 22dc7fd was 22dc7fd, checked in by mselvaggi <mselvaggi@…>, 11 years ago

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