Fork me on GitHub

source: git/modules/Calorimeter.cc@ 8624f58

ImprovedOutputFile Timing dual_readout llp
Last change on this file since 8624f58 was 8624f58, checked in by Pavel Demin <pavel.demin@…>, 10 years ago

replace Ecal with ECal and Hcal with HCal
  • 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 * $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;
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
142 ecalFraction = paramFractions[0].GetDouble();
143 hcalFraction = paramFractions[1].GetDouble();
144
145 fFractionMap[param[i*2].GetInt()] = make_pair(ecalFraction, hcalFraction);
146 }
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("ECalMinEnergy", 0.0);
158 fHCalEnergyMin = GetDouble("HCalMinEnergy", 0.0);
159
160 fECalSigmaMin = GetDouble("ECalMinSignificance", 0.0);
161 fHCalSigmaMin = GetDouble("HCalMinSignificance", 0.0);
162
163 // switch on or off the dithering of the center of calorimeter towers
164 fDitherTowerCenter = GetBool("DitherTowerCenter", true);
165
166 // read resolution formulas
167 fECalResolutionFormula->Compile(GetString("ECalResolutionFormula", "0"));
168 fHCalResolutionFormula->Compile(GetString("HCalResolutionFormula", "0"));
169
170 // import array with output from other modules
171 fParticleInputArray = ImportArray(GetString("ParticleInputArray", "ParticlePropagator/particles"));
172 fItParticleInputArray = fParticleInputArray->MakeIterator();
173
174 fTrackInputArray = ImportArray(GetString("TrackInputArray", "ParticlePropagator/tracks"));
175 fItTrackInputArray = fTrackInputArray->MakeIterator();
176
177 // create output arrays
178 fTowerOutputArray = ExportArray(GetString("TowerOutputArray", "towers"));
179 fPhotonOutputArray = ExportArray(GetString("PhotonOutputArray", "photons"));
180
181 fEFlowTrackOutputArray = ExportArray(GetString("EFlowTrackOutputArray", "eflowTracks"));
182 fEFlowPhotonOutputArray = ExportArray(GetString("EFlowPhotonOutputArray", "eflowPhotons"));
183 fEFlowNeutralHadronOutputArray = ExportArray(GetString("EFlowNeutralHadronOutputArray", "eflowNeutralHadrons"));
184
185}
186
187//------------------------------------------------------------------------------
188
189void Calorimeter::Finish()
190{
191 vector< vector< Double_t >* >::iterator itPhiBin;
192 if(fItParticleInputArray) delete fItParticleInputArray;
193 if(fItTrackInputArray) delete fItTrackInputArray;
194 for(itPhiBin = fPhiBins.begin(); itPhiBin != fPhiBins.end(); ++itPhiBin)
195 {
196 delete *itPhiBin;
197 }
198}
199
200//------------------------------------------------------------------------------
201
202void Calorimeter::Process()
203{
204 Candidate *particle, *track;
205 TLorentzVector position, momentum;
206 Short_t etaBin, phiBin, flags;
207 Int_t number;
208 Long64_t towerHit, towerEtaPhi, hitEtaPhi;
209 Double_t ecalFraction, hcalFraction;
210 Double_t ecalEnergy, hcalEnergy;
211 Int_t pdgCode;
212
213 TFractionMap::iterator itFractionMap;
214
215 vector< Double_t >::iterator itEtaBin;
216 vector< Double_t >::iterator itPhiBin;
217 vector< Double_t > *phiBins;
218
219 vector< Long64_t >::iterator itTowerHits;
220
221 DelphesFactory *factory = GetFactory();
222 fTowerHits.clear();
223 fTowerECalFractions.clear();
224 fTowerHCalFractions.clear();
225 fTrackECalFractions.clear();
226 fTrackHCalFractions.clear();
227
228 // loop over all particles
229 fItParticleInputArray->Reset();
230 number = -1;
231 while((particle = static_cast<Candidate*>(fItParticleInputArray->Next())))
232 {
233 const TLorentzVector &particlePosition = particle->Position;
234 ++number;
235
236 pdgCode = TMath::Abs(particle->PID);
237
238 itFractionMap = fFractionMap.find(pdgCode);
239 if(itFractionMap == fFractionMap.end())
240 {
241 itFractionMap = fFractionMap.find(0);
242 }
243
244 ecalFraction = itFractionMap->second.first;
245 hcalFraction = itFractionMap->second.second;
246
247 fTowerECalFractions.push_back(ecalFraction);
248 fTowerHCalFractions.push_back(hcalFraction);
249
250 if(ecalFraction < 1.0E-9 && hcalFraction < 1.0E-9) continue;
251
252 // find eta bin [1, fEtaBins.size - 1]
253 itEtaBin = lower_bound(fEtaBins.begin(), fEtaBins.end(), particlePosition.Eta());
254 if(itEtaBin == fEtaBins.begin() || itEtaBin == fEtaBins.end()) continue;
255 etaBin = distance(fEtaBins.begin(), itEtaBin);
256
257 // phi bins for given eta bin
258 phiBins = fPhiBins[etaBin];
259
260 // find phi bin [1, phiBins.size - 1]
261 itPhiBin = lower_bound(phiBins->begin(), phiBins->end(), particlePosition.Phi());
262 if(itPhiBin == phiBins->begin() || itPhiBin == phiBins->end()) continue;
263 phiBin = distance(phiBins->begin(), itPhiBin);
264
265 flags = 0;
266 flags |= (pdgCode == 11 || pdgCode == 22) << 1;
267
268 // make tower hit {16-bits for eta bin number, 16-bits for phi bin number, 8-bits for flags, 24-bits for particle number}
269 towerHit = (Long64_t(etaBin) << 48) | (Long64_t(phiBin) << 32) | (Long64_t(flags) << 24) | Long64_t(number);
270
271 fTowerHits.push_back(towerHit);
272 }
273
274 // loop over all tracks
275 fItTrackInputArray->Reset();
276 number = -1;
277 while((track = static_cast<Candidate*>(fItTrackInputArray->Next())))
278 {
279 const TLorentzVector &trackPosition = track->Position;
280 ++number;
281
282 pdgCode = TMath::Abs(track->PID);
283
284 itFractionMap = fFractionMap.find(pdgCode);
285 if(itFractionMap == fFractionMap.end())
286 {
287 itFractionMap = fFractionMap.find(0);
288 }
289
290 ecalFraction = itFractionMap->second.first;
291 hcalFraction = itFractionMap->second.second;
292
293 fTrackECalFractions.push_back(ecalFraction);
294 fTrackHCalFractions.push_back(hcalFraction);
295
296 // find eta bin [1, fEtaBins.size - 1]
297 itEtaBin = lower_bound(fEtaBins.begin(), fEtaBins.end(), trackPosition.Eta());
298 if(itEtaBin == fEtaBins.begin() || itEtaBin == fEtaBins.end()) continue;
299 etaBin = distance(fEtaBins.begin(), itEtaBin);
300
301 // phi bins for given eta bin
302 phiBins = fPhiBins[etaBin];
303
304 // find phi bin [1, phiBins.size - 1]
305 itPhiBin = lower_bound(phiBins->begin(), phiBins->end(), trackPosition.Phi());
306 if(itPhiBin == phiBins->begin() || itPhiBin == phiBins->end()) continue;
307 phiBin = distance(phiBins->begin(), itPhiBin);
308
309 flags = 1;
310
311 // make tower hit {16-bits for eta bin number, 16-bits for phi bin number, 8-bits for flags, 24-bits for track number}
312 towerHit = (Long64_t(etaBin) << 48) | (Long64_t(phiBin) << 32) | (Long64_t(flags) << 24) | Long64_t(number);
313
314 fTowerHits.push_back(towerHit);
315 }
316
317 // all hits are sorted first by eta bin number, then by phi bin number,
318 // then by flags and then by particle or track number
319 sort(fTowerHits.begin(), fTowerHits.end());
320
321 // loop over all hits
322 towerEtaPhi = 0;
323 fTower = 0;
324 for(itTowerHits = fTowerHits.begin(); itTowerHits != fTowerHits.end(); ++itTowerHits)
325 {
326 towerHit = (*itTowerHits);
327 flags = (towerHit >> 24) & 0x00000000000000FFLL;
328 number = (towerHit) & 0x0000000000FFFFFFLL;
329 hitEtaPhi = towerHit >> 32;
330
331 if(towerEtaPhi != hitEtaPhi)
332 {
333 // switch to next tower
334 towerEtaPhi = hitEtaPhi;
335
336 // finalize previous tower
337 FinalizeTower();
338
339 // create new tower
340 fTower = factory->NewCandidate();
341
342 phiBin = (towerHit >> 32) & 0x000000000000FFFFLL;
343 etaBin = (towerHit >> 48) & 0x000000000000FFFFLL;
344
345 // phi bins for given eta bin
346 phiBins = fPhiBins[etaBin];
347
348 // calculate eta and phi of the tower's center
349 fTowerEta = 0.5*(fEtaBins[etaBin - 1] + fEtaBins[etaBin]);
350 fTowerPhi = 0.5*((*phiBins)[phiBin - 1] + (*phiBins)[phiBin]);
351
352 fTowerEdges[0] = fEtaBins[etaBin - 1];
353 fTowerEdges[1] = fEtaBins[etaBin];
354 fTowerEdges[2] = (*phiBins)[phiBin - 1];
355 fTowerEdges[3] = (*phiBins)[phiBin];
356
357 fTowerECalEnergy = 0.0;
358 fTowerHCalEnergy = 0.0;
359
360 fTrackECalEnergy = 0.0;
361 fTrackHCalEnergy = 0.0;
362
363 fTowerECalTime = 0.0;
364 fTowerHCalTime = 0.0;
365
366 fTrackECalTime = 0.0;
367 fTrackHCalTime = 0.0;
368
369 fTowerECalWeightTime = 0.0;
370 fTowerHCalWeightTime = 0.0;
371
372 fTowerTrackHits = 0;
373 fTowerPhotonHits = 0;
374
375 fTowerTrackArray->Clear();
376 }
377
378 // check for track hits
379 if(flags & 1)
380 {
381 ++fTowerTrackHits;
382
383 track = static_cast<Candidate*>(fTrackInputArray->At(number));
384 momentum = track->Momentum;
385 position = track->Position;
386
387
388 ecalEnergy = momentum.E() * fTrackECalFractions[number];
389 hcalEnergy = momentum.E() * fTrackHCalFractions[number];
390
391 fTrackECalEnergy += ecalEnergy;
392 fTrackHCalEnergy += hcalEnergy;
393
394 fTrackECalTime += TMath::Sqrt(ecalEnergy)*position.T();
395 fTrackHCalTime += TMath::Sqrt(hcalEnergy)*position.T();
396
397 fTrackECalWeightTime += TMath::Sqrt(ecalEnergy);
398 fTrackHCalWeightTime += TMath::Sqrt(hcalEnergy);
399
400 fTowerTrackArray->Add(track);
401
402 continue;
403 }
404
405 // check for photon and electron hits in current tower
406 if(flags & 2) ++fTowerPhotonHits;
407
408 particle = static_cast<Candidate*>(fParticleInputArray->At(number));
409 momentum = particle->Momentum;
410 position = particle->Position;
411
412 // fill current tower
413 ecalEnergy = momentum.E() * fTowerECalFractions[number];
414 hcalEnergy = momentum.E() * fTowerHCalFractions[number];
415
416 fTowerECalEnergy += ecalEnergy;
417 fTowerHCalEnergy += hcalEnergy;
418
419 fTowerECalTime += TMath::Sqrt(ecalEnergy)*position.T();
420 fTowerHCalTime += TMath::Sqrt(hcalEnergy)*position.T();
421
422 fTowerECalWeightTime += TMath::Sqrt(ecalEnergy);
423 fTowerHCalWeightTime += TMath::Sqrt(hcalEnergy);
424
425
426 fTower->AddCandidate(particle);
427 }
428
429 // finalize last tower
430 FinalizeTower();
431}
432
433//------------------------------------------------------------------------------
434
435void Calorimeter::FinalizeTower()
436{
437 Candidate *track, *tower;
438 Double_t energy, pt, eta, phi;
439 Double_t ecalEnergy, hcalEnergy;
440 Double_t ecalSigma, hcalSigma;
441 Double_t ecalTime, hcalTime, time;
442
443 if(!fTower) return;
444
445 ecalSigma = fECalResolutionFormula->Eval(0.0, fTowerEta, 0.0, fTowerECalEnergy);
446
447 ecalEnergy = LogNormal(fTowerECalEnergy, ecalSigma);
448 ecalTime = (fTowerECalWeightTime < 1.0E-09 ) ? 0 : fTowerECalTime/fTowerECalWeightTime;
449
450 hcalSigma = fHCalResolutionFormula->Eval(0.0, fTowerEta, 0.0, fTowerHCalEnergy);
451
452 hcalEnergy = LogNormal(fTowerHCalEnergy, hcalSigma);
453 hcalTime = (fTowerHCalWeightTime < 1.0E-09 ) ? 0 : fTowerHCalTime/fTowerHCalWeightTime;
454
455
456 ecalSigma = fECalResolutionFormula->Eval(0.0, fTowerEta, 0.0, ecalEnergy);
457 hcalSigma = fHCalResolutionFormula->Eval(0.0, fTowerEta, 0.0, hcalEnergy);
458
459 ecalEnergy = (ecalEnergy < fECalEnergyMin || ecalEnergy < fECalSigmaMin*ecalSigma) ? 0 : ecalEnergy;
460 hcalEnergy = (hcalEnergy < fHCalEnergyMin || hcalEnergy < fHCalSigmaMin*hcalSigma) ? 0 : hcalEnergy;
461
462 energy = ecalEnergy + hcalEnergy;
463 time = (TMath::Sqrt(ecalEnergy)*ecalTime + TMath::Sqrt(hcalEnergy)*hcalTime)/(TMath::Sqrt(ecalEnergy) + TMath::Sqrt(hcalEnergy));
464
465 if(fDitherTowerCenter)
466 {
467 eta = gRandom->Uniform(fTowerEdges[0], fTowerEdges[1]);
468 phi = gRandom->Uniform(fTowerEdges[2], fTowerEdges[3]);
469 }
470 else
471 {
472 eta = fTowerEta;
473 phi = fTowerPhi;
474 }
475
476 pt = energy / TMath::CosH(eta);
477
478 fTower->Position.SetPtEtaPhiE(1.0, eta, phi, time);
479 fTower->Momentum.SetPtEtaPhiE(pt, eta, phi, energy);
480 fTower->Eem = ecalEnergy;
481 fTower->Ehad = hcalEnergy;
482
483 fTower->Edges[0] = fTowerEdges[0];
484 fTower->Edges[1] = fTowerEdges[1];
485 fTower->Edges[2] = fTowerEdges[2];
486 fTower->Edges[3] = fTowerEdges[3];
487
488 if( energy > 0.0 )
489 {
490 if(fTowerPhotonHits > 0 && fTowerTrackHits == 0)
491 {
492 fPhotonOutputArray->Add(fTower);
493 }
494
495 fTowerOutputArray->Add(fTower);
496 }
497
498 // fill energy flow candidates
499
500 // save all the tracks as energy flow tracks
501 fItTowerTrackArray->Reset();
502 while((track = static_cast<Candidate*>(fItTowerTrackArray->Next())))
503 {
504 fEFlowTrackOutputArray->Add(track);
505 }
506
507 ecalEnergy -= fTrackECalEnergy;
508 if(ecalEnergy < fECalEnergyMin || ecalEnergy < fECalSigmaMin*fECalResolutionFormula->Eval(0.0, fTowerEta, 0.0, ecalEnergy)) ecalEnergy = 0.0;
509
510 hcalEnergy -= fTrackHCalEnergy;
511 if(hcalEnergy < fHCalEnergyMin || hcalEnergy < fHCalSigmaMin*fHCalResolutionFormula->Eval(0.0, fTowerEta, 0.0, hcalEnergy)) hcalEnergy = 0.0;
512
513 energy = ecalEnergy + hcalEnergy;
514
515 if(ecalEnergy > 0.0)
516 {
517 // create new photon tower
518 tower = static_cast<Candidate*>(fTower->Clone());
519
520 pt = ecalEnergy / TMath::CosH(eta);
521
522 tower->Momentum.SetPtEtaPhiE(pt, eta, phi, ecalEnergy);
523 tower->Eem = ecalEnergy;
524 tower->Ehad = 0;
525
526 fEFlowPhotonOutputArray->Add(tower);
527 }
528 if(hcalEnergy > 0.0)
529 {
530 // create new neutral hadron tower
531 tower = static_cast<Candidate*>(fTower->Clone());
532
533 pt = hcalEnergy / TMath::CosH(eta);
534
535 tower->Momentum.SetPtEtaPhiE(pt, eta, phi, hcalEnergy);
536 tower->Eem = 0;
537 tower->Ehad = hcalEnergy;
538
539 fEFlowNeutralHadronOutputArray->Add(tower);
540 }
541}
542
543//------------------------------------------------------------------------------
544
545Double_t Calorimeter::LogNormal(Double_t mean, Double_t sigma)
546{
547 Double_t a, b;
548
549 if(mean > 0.0)
550 {
551 b = TMath::Sqrt(TMath::Log((1.0 + (sigma*sigma)/(mean*mean))));
552 a = TMath::Log(mean) - 0.5*b*b;
553
554 return TMath::Exp(a + b*gRandom->Gaus(0, 1));
555 }
556 else
557 {
558 return 0.0;
559 }
560}
Note: See TracBrowser for help on using the repository browser.