Fork me on GitHub

source: git/modules/Calorimeter.cc@ 122e1e5

ImprovedOutputFile Timing dual_readout llp
Last change on this file since 122e1e5 was 839deb7, checked in by Pavel Demin <pavel.demin@…>, 9 years ago

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