Fork me on GitHub

source: git/modules/Calorimeter.cc@ 1d1f6a4

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

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