Fork me on GitHub

source: git/modules/Calorimeter.cc@ 0c871e5

ImprovedOutputFile Timing dual_readout llp
Last change on this file since 0c871e5 was 01f457a, checked in by Michele <michele.selvaggi@…>, 10 years ago

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