Fork me on GitHub

source: git/modules/SimpleCalorimeter.cc@ 04290b1

ImprovedOutputFile Timing dual_readout llp
Last change on this file since 04290b1 was 2628d2b, checked in by Michele Selvaggi <michele.selvaggi@…>, 8 years ago

energy weighted time in ECAL

  • Property mode set to 100644
File size: 14.8 KB
Line 
1/*
2 * Delphes: a framework for fast simulation of a generic collider experiment
3 * Copyright (C) 2012-2014 Universite catholique de Louvain (UCL), Belgium
4 *
5 * This program is free software: you can redistribute it and/or modify
6 * it under the terms of the GNU General Public License as published by
7 * the Free Software Foundation, either version 3 of the License, or
8 * (at your option) any later version.
9 *
10 * This program is distributed in the hope that it will be useful,
11 * but WITHOUT ANY WARRANTY; without even the implied warranty of
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 * GNU General Public License for more details.
14 *
15 * You should have received a copy of the GNU General Public License
16 * along with this program. If not, see <http://www.gnu.org/licenses/>.
17 */
18
19
20/** \class SimpleCalorimeter
21 *
22 * Fills SimpleCalorimeter towers, performs SimpleCalorimeter resolution smearing,
23 * and creates energy flow objects (tracks, photons, and neutral hadrons).
24 *
25 * \author P. Demin - UCL, Louvain-la-Neuve
26 *
27 */
28
29#include "modules/SimpleCalorimeter.h"
30
31#include "classes/DelphesClasses.h"
32#include "classes/DelphesFactory.h"
33#include "classes/DelphesFormula.h"
34
35#include "ExRootAnalysis/ExRootResult.h"
36#include "ExRootAnalysis/ExRootFilter.h"
37#include "ExRootAnalysis/ExRootClassifier.h"
38
39#include "TMath.h"
40#include "TString.h"
41#include "TFormula.h"
42#include "TRandom3.h"
43#include "TObjArray.h"
44#include "TDatabasePDG.h"
45#include "TLorentzVector.h"
46
47#include <algorithm>
48#include <stdexcept>
49#include <iostream>
50#include <sstream>
51
52using namespace std;
53
54//------------------------------------------------------------------------------
55
56SimpleCalorimeter::SimpleCalorimeter() :
57 fResolutionFormula(0),
58 fItParticleInputArray(0), fItTrackInputArray(0)
59{
60 Int_t i;
61
62 fResolutionFormula = new DelphesFormula;
63
64 for(i = 0; i < 2; ++i)
65 {
66 fTowerTrackArray[i] = new TObjArray;
67 fItTowerTrackArray[i] = fTowerTrackArray[i]->MakeIterator();
68 }
69}
70
71//------------------------------------------------------------------------------
72
73SimpleCalorimeter::~SimpleCalorimeter()
74{
75 Int_t i;
76
77 if(fResolutionFormula) delete fResolutionFormula;
78
79 for(i = 0; i < 2; ++i)
80 {
81 if(fTowerTrackArray[i]) delete fTowerTrackArray[i];
82 if(fItTowerTrackArray[i]) delete fItTowerTrackArray[i];
83 }
84}
85
86//------------------------------------------------------------------------------
87
88void SimpleCalorimeter::Init()
89{
90 ExRootConfParam param, paramEtaBins, paramPhiBins, paramFractions;
91 Long_t i, j, k, size, sizeEtaBins, sizePhiBins;
92 Double_t fraction;
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] = 1.0;
140
141 for(i = 0; i < size/2; ++i)
142 {
143 paramFractions = param[i*2 + 1];
144 fraction = paramFractions[0].GetDouble();
145 fFractionMap[param[i*2].GetInt()] = fraction;
146 }
147
148 // read min E value for towers to be saved
149 fEnergyMin = GetDouble("EnergyMin", 0.0);
150
151 fEnergySignificanceMin = GetDouble("EnergySignificanceMin", 0.0);
152
153 // flag that says if current calo is Ecal of Hcal (will then fill correct values of Eem and Ehad)
154 fIsEcal = GetBool("IsEcal", false);
155
156 // switch on or off the dithering of the center of calorimeter towers
157 fSmearTowerCenter = GetBool("SmearTowerCenter", true);
158
159 // read resolution formulas
160 fResolutionFormula->Compile(GetString("ResolutionFormula", "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
172 fEFlowTrackOutputArray = ExportArray(GetString("EFlowTrackOutputArray", "eflowTracks"));
173 fEFlowTowerOutputArray = ExportArray(GetString("EFlowTowerOutputArray", "eflowTowers"));
174}
175
176//------------------------------------------------------------------------------
177
178void SimpleCalorimeter::Finish()
179{
180 vector< vector< Double_t >* >::iterator itPhiBin;
181 if(fItParticleInputArray) delete fItParticleInputArray;
182 if(fItTrackInputArray) delete fItTrackInputArray;
183 for(itPhiBin = fPhiBins.begin(); itPhiBin != fPhiBins.end(); ++itPhiBin)
184 {
185 delete *itPhiBin;
186 }
187}
188
189//------------------------------------------------------------------------------
190
191void SimpleCalorimeter::Process()
192{
193 Candidate *particle, *track;
194 TLorentzVector position, momentum;
195 Short_t etaBin, phiBin, flags;
196 Int_t number;
197 Long64_t towerHit, towerEtaPhi, hitEtaPhi;
198 Double_t fraction;
199 Double_t energy;
200 Double_t sigma;
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 fTowerFractions.clear();
214 fTrackFractions.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 fraction = itFractionMap->second;
233 fTowerFractions.push_back(fraction);
234
235 if(fraction < 1.0E-9) continue;
236
237 // find eta bin [1, fEtaBins.size - 1]
238 itEtaBin = lower_bound(fEtaBins.begin(), fEtaBins.end(), particlePosition.Eta());
239 if(itEtaBin == fEtaBins.begin() || itEtaBin == fEtaBins.end()) continue;
240 etaBin = distance(fEtaBins.begin(), itEtaBin);
241
242 // phi bins for given eta bin
243 phiBins = fPhiBins[etaBin];
244
245 // find phi bin [1, phiBins.size - 1]
246 itPhiBin = lower_bound(phiBins->begin(), phiBins->end(), particlePosition.Phi());
247 if(itPhiBin == phiBins->begin() || itPhiBin == phiBins->end()) continue;
248 phiBin = distance(phiBins->begin(), itPhiBin);
249
250 flags = 0;
251 flags |= (pdgCode == 11 || pdgCode == 22) << 1;
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 fraction = itFractionMap->second;
276
277 fTrackFractions.push_back(fraction);
278
279 // find eta bin [1, fEtaBins.size - 1]
280 itEtaBin = lower_bound(fEtaBins.begin(), fEtaBins.end(), trackPosition.Eta());
281 if(itEtaBin == fEtaBins.begin() || itEtaBin == fEtaBins.end()) continue;
282 etaBin = distance(fEtaBins.begin(), itEtaBin);
283
284 // phi bins for given eta bin
285 phiBins = fPhiBins[etaBin];
286
287 // find phi bin [1, phiBins.size - 1]
288 itPhiBin = lower_bound(phiBins->begin(), phiBins->end(), trackPosition.Phi());
289 if(itPhiBin == phiBins->begin() || itPhiBin == phiBins->end()) continue;
290 phiBin = distance(phiBins->begin(), itPhiBin);
291
292 flags = 1;
293
294 // make tower hit {16-bits for eta bin number, 16-bits for phi bin number, 8-bits for flags, 24-bits for track number}
295 towerHit = (Long64_t(etaBin) << 48) | (Long64_t(phiBin) << 32) | (Long64_t(flags) << 24) | Long64_t(number);
296
297 fTowerHits.push_back(towerHit);
298 }
299
300 // all hits are sorted first by eta bin number, then by phi bin number,
301 // then by flags and then by particle or track number
302 sort(fTowerHits.begin(), fTowerHits.end());
303
304 // loop over all hits
305 towerEtaPhi = 0;
306 fTower = 0;
307 for(itTowerHits = fTowerHits.begin(); itTowerHits != fTowerHits.end(); ++itTowerHits)
308 {
309 towerHit = (*itTowerHits);
310 flags = (towerHit >> 24) & 0x00000000000000FFLL;
311 number = (towerHit) & 0x0000000000FFFFFFLL;
312 hitEtaPhi = towerHit >> 32;
313
314 if(towerEtaPhi != hitEtaPhi)
315 {
316 // switch to next tower
317 towerEtaPhi = hitEtaPhi;
318
319 // finalize previous tower
320 FinalizeTower();
321
322 // create new tower
323 fTower = factory->NewCandidate();
324
325 phiBin = (towerHit >> 32) & 0x000000000000FFFFLL;
326 etaBin = (towerHit >> 48) & 0x000000000000FFFFLL;
327
328 // phi bins for given eta bin
329 phiBins = fPhiBins[etaBin];
330
331 // calculate eta and phi of the tower's center
332 fTowerEta = 0.5*(fEtaBins[etaBin - 1] + fEtaBins[etaBin]);
333 fTowerPhi = 0.5*((*phiBins)[phiBin - 1] + (*phiBins)[phiBin]);
334
335 fTowerEdges[0] = fEtaBins[etaBin - 1];
336 fTowerEdges[1] = fEtaBins[etaBin];
337 fTowerEdges[2] = (*phiBins)[phiBin - 1];
338 fTowerEdges[3] = (*phiBins)[phiBin];
339
340 fTowerEnergy = 0.0;
341
342 fTrackEnergy[0] = 0.0;
343 fTrackEnergy[1] = 0.0;
344
345 fTowerTime = 0.0;
346 fTrackTime = 0.0;
347
348 fTowerTimeWeight = 0.0;
349
350 fTowerTrackHits = 0;
351 fTowerPhotonHits = 0;
352
353 fTowerTrackArray[0]->Clear();
354 fTowerTrackArray[1]->Clear();
355 }
356
357 // check for track hits
358 if(flags & 1)
359 {
360 ++fTowerTrackHits;
361
362 track = static_cast<Candidate*>(fTrackInputArray->At(number));
363 momentum = track->Momentum;
364 position = track->Position;
365
366 energy = momentum.E() * fTrackFractions[number];
367
368 fTrackTime += energy*position.T();
369 fTrackTimeWeight += energy;
370
371 if(fTrackFractions[number] > 1.0E-9)
372 {
373 sigma = fResolutionFormula->Eval(0.0, fTowerEta, 0.0, momentum.E());
374 if(sigma/momentum.E() < track->TrackResolution)
375 {
376 fTrackEnergy[0] += energy;
377 fTowerTrackArray[0]->Add(track);
378 }
379 else
380 {
381 fTrackEnergy[1] += energy;
382 fTowerTrackArray[1]->Add(track);
383 }
384 }
385 else
386 {
387 fEFlowTrackOutputArray->Add(track);
388 }
389
390 continue;
391 }
392
393 // check for photon and electron hits in current tower
394 if(flags & 2) ++fTowerPhotonHits;
395
396 particle = static_cast<Candidate*>(fParticleInputArray->At(number));
397 momentum = particle->Momentum;
398 position = particle->Position;
399
400 // fill current tower
401 energy = momentum.E() * fTowerFractions[number];
402
403 fTowerEnergy += energy;
404
405 fTowerTime += energy*position.T();
406 fTowerTimeWeight += energy;
407
408 fTower->AddCandidate(particle);
409 }
410
411 // finalize last tower
412 FinalizeTower();
413}
414
415//------------------------------------------------------------------------------
416
417void SimpleCalorimeter::FinalizeTower()
418{
419 Candidate *tower, *track, *mother;
420 Double_t energy, pt, eta, phi;
421 Double_t sigma;
422 Double_t time;
423
424 TLorentzVector momentum;
425 TFractionMap::iterator itFractionMap;
426
427 if(!fTower) return;
428
429 sigma = fResolutionFormula->Eval(0.0, fTowerEta, 0.0, fTowerEnergy);
430
431 energy = LogNormal(fTowerEnergy, sigma);
432
433 time = (fTowerTimeWeight < 1.0E-09 ) ? 0.0 : fTowerTime/fTowerTimeWeight;
434
435 sigma = fResolutionFormula->Eval(0.0, fTowerEta, 0.0, energy);
436
437 if(energy < fEnergyMin || energy < fEnergySignificanceMin*sigma) energy = 0.0;
438
439 if(fSmearTowerCenter)
440 {
441 eta = gRandom->Uniform(fTowerEdges[0], fTowerEdges[1]);
442 phi = gRandom->Uniform(fTowerEdges[2], fTowerEdges[3]);
443 }
444 else
445 {
446 eta = fTowerEta;
447 phi = fTowerPhi;
448 }
449
450 pt = energy / TMath::CosH(eta);
451
452 fTower->Position.SetPtEtaPhiE(1.0, eta, phi, time);
453 fTower->Momentum.SetPtEtaPhiE(pt, eta, phi, energy);
454
455 fTower->Eem = (!fIsEcal) ? 0 : energy;
456 fTower->Ehad = (fIsEcal) ? 0 : energy;
457
458 fTower->Edges[0] = fTowerEdges[0];
459 fTower->Edges[1] = fTowerEdges[1];
460 fTower->Edges[2] = fTowerEdges[2];
461 fTower->Edges[3] = fTowerEdges[3];
462
463 // fill SimpleCalorimeter towers
464 if(energy > 0.0) fTowerOutputArray->Add(fTower);
465
466 // fill e-flow candidates
467
468 energy -= fTrackEnergy[1];
469
470 fItTowerTrackArray[0]->Reset();
471 while((track = static_cast<Candidate*>(fItTowerTrackArray[0]->Next())))
472 {
473 mother = track;
474 track = static_cast<Candidate*>(track->Clone());
475 track->AddCandidate(mother);
476
477 track->Momentum *= energy/fTrackEnergy[0];
478
479 fEFlowTrackOutputArray->Add(track);
480 }
481
482 fItTowerTrackArray[1]->Reset();
483 while((track = static_cast<Candidate*>(fItTowerTrackArray[1]->Next())))
484 {
485 mother = track;
486 track = static_cast<Candidate*>(track->Clone());
487 track->AddCandidate(mother);
488
489 fEFlowTrackOutputArray->Add(track);
490 }
491
492 if(fTowerTrackArray[0]->GetEntriesFast() > 0) energy = 0.0;
493
494 sigma = fResolutionFormula->Eval(0.0, fTowerEta, 0.0, energy);
495 if(energy < fEnergyMin || energy < fEnergySignificanceMin*sigma) energy = 0.0;
496
497 // save energy excess as an energy flow tower
498 if(energy > 0.0)
499 {
500 // create new photon tower
501 tower = static_cast<Candidate*>(fTower->Clone());
502 pt = energy / TMath::CosH(eta);
503
504 tower->Eem = (!fIsEcal) ? 0 : energy;
505 tower->Ehad = (fIsEcal) ? 0 : energy;
506
507 tower->Momentum.SetPtEtaPhiE(pt, eta, phi, energy);
508
509 tower->PID = (fIsEcal) ? 22 : 0;
510
511 fEFlowTowerOutputArray->Add(tower);
512 }
513}
514
515//------------------------------------------------------------------------------
516
517Double_t SimpleCalorimeter::LogNormal(Double_t mean, Double_t sigma)
518{
519 Double_t a, b;
520
521 if(mean > 0.0)
522 {
523 b = TMath::Sqrt(TMath::Log((1.0 + (sigma*sigma)/(mean*mean))));
524 a = TMath::Log(mean) - 0.5*b*b;
525
526 return TMath::Exp(a + b*gRandom->Gaus(0.0, 1.0));
527 }
528 else
529 {
530 return 0.0;
531 }
532}
Note: See TracBrowser for help on using the repository browser.