Fork me on GitHub

source: git/modules/SimpleCalorimeter.cc@ 28027d5

ImprovedOutputFile Timing dual_readout llp
Last change on this file since 28027d5 was fc6300d, checked in by Michele Selvaggi <selvaggi@…>, 9 years ago

added isEcal flag to simple calorimeter and fill Eem and Ehad variables accordingly (answer to #639)

  • Property mode set to 100644
File size: 13.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 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 fTowerTrackArray(0), fItTowerTrackArray(0)
60{
61 fResolutionFormula = new DelphesFormula;
62
63 fTowerTrackArray = new TObjArray;
64 fItTowerTrackArray = fTowerTrackArray->MakeIterator();
65}
66
67//------------------------------------------------------------------------------
68
69SimpleCalorimeter::~SimpleCalorimeter()
70{
71 if(fResolutionFormula) delete fResolutionFormula;
72
73 if(fTowerTrackArray) delete fTowerTrackArray;
74 if(fItTowerTrackArray) delete fItTowerTrackArray;
75}
76
77//------------------------------------------------------------------------------
78
79void SimpleCalorimeter::Init()
80{
81 ExRootConfParam param, paramEtaBins, paramPhiBins, paramFractions;
82 Long_t i, j, k, size, sizeEtaBins, sizePhiBins;
83 Double_t fraction;
84 TBinMap::iterator itEtaBin;
85 set< Double_t >::iterator itPhiBin;
86 vector< Double_t > *phiBins;
87
88 // read eta and phi bins
89 param = GetParam("EtaPhiBins");
90 size = param.GetSize();
91 fBinMap.clear();
92 fEtaBins.clear();
93 fPhiBins.clear();
94 for(i = 0; i < size/2; ++i)
95 {
96 paramEtaBins = param[i*2];
97 sizeEtaBins = paramEtaBins.GetSize();
98 paramPhiBins = param[i*2 + 1];
99 sizePhiBins = paramPhiBins.GetSize();
100
101 for(j = 0; j < sizeEtaBins; ++j)
102 {
103 for(k = 0; k < sizePhiBins; ++k)
104 {
105 fBinMap[paramEtaBins[j].GetDouble()].insert(paramPhiBins[k].GetDouble());
106 }
107 }
108 }
109
110 // for better performance we transform map of sets to parallel vectors:
111 // vector< double > and vector< vector< double >* >
112 for(itEtaBin = fBinMap.begin(); itEtaBin != fBinMap.end(); ++itEtaBin)
113 {
114 fEtaBins.push_back(itEtaBin->first);
115 phiBins = new vector< double >(itEtaBin->second.size());
116 fPhiBins.push_back(phiBins);
117 phiBins->clear();
118 for(itPhiBin = itEtaBin->second.begin(); itPhiBin != itEtaBin->second.end(); ++itPhiBin)
119 {
120 phiBins->push_back(*itPhiBin);
121 }
122 }
123
124 // read energy fractions for different particles
125 param = GetParam("EnergyFraction");
126 size = param.GetSize();
127
128 // set default energy fractions values
129 fFractionMap.clear();
130 fFractionMap[0] = 1.0;
131
132 for(i = 0; i < size/2; ++i)
133 {
134 paramFractions = param[i*2 + 1];
135 fraction = paramFractions[0].GetDouble();
136 fFractionMap[param[i*2].GetInt()] = fraction;
137 }
138
139/*
140 TFractionMap::iterator itFractionMap;
141 for(itFractionMap = fFractionMap.begin(); itFractionMap != fFractionMap.end(); ++itFractionMap)
142 {
143 cout << itFractionMap->first << " " << itFractionMap->second.first << " " << itFractionMap->second.second << endl;
144 }
145*/
146
147 // read min E value for towers to be saved
148 fEnergyMin = GetDouble("EnergyMin", 0.0);
149
150 fEnergySignificanceMin = GetDouble("EnergySignificanceMin", 0.0);
151
152 // flag that says if current calo is Ecal of Hcal (will then fill correct values of Eem and Ehad)
153 fIsEcal = GetBool("IsEcal", false);
154
155 // switch on or off the dithering of the center of calorimeter towers
156 fSmearTowerCenter = GetBool("SmearTowerCenter", true);
157
158 // read resolution formulas
159 fResolutionFormula->Compile(GetString("ResolutionFormula", "0"));
160
161 // import array with output from other modules
162 fParticleInputArray = ImportArray(GetString("ParticleInputArray", "ParticlePropagator/particles"));
163 fItParticleInputArray = fParticleInputArray->MakeIterator();
164
165 fTrackInputArray = ImportArray(GetString("TrackInputArray", "ParticlePropagator/tracks"));
166 fItTrackInputArray = fTrackInputArray->MakeIterator();
167
168 // create output arrays
169 fTowerOutputArray = ExportArray(GetString("TowerOutputArray", "towers"));
170 fEFlowTowerOutputArray = ExportArray(GetString("EFlowTowerOutputArray", "eflowTowers"));
171}
172
173//------------------------------------------------------------------------------
174
175void SimpleCalorimeter::Finish()
176{
177 vector< vector< Double_t >* >::iterator itPhiBin;
178 if(fItParticleInputArray) delete fItParticleInputArray;
179 if(fItTrackInputArray) delete fItTrackInputArray;
180 for(itPhiBin = fPhiBins.begin(); itPhiBin != fPhiBins.end(); ++itPhiBin)
181 {
182 delete *itPhiBin;
183 }
184}
185
186//------------------------------------------------------------------------------
187
188void SimpleCalorimeter::Process()
189{
190 Candidate *particle, *track;
191 TLorentzVector position, momentum;
192 Short_t etaBin, phiBin, flags;
193 Int_t number;
194 Long64_t towerHit, towerEtaPhi, hitEtaPhi;
195 Double_t fraction;
196 Double_t energy;
197 Int_t pdgCode;
198
199 TFractionMap::iterator itFractionMap;
200
201 vector< Double_t >::iterator itEtaBin;
202 vector< Double_t >::iterator itPhiBin;
203 vector< Double_t > *phiBins;
204
205 vector< Long64_t >::iterator itTowerHits;
206
207 DelphesFactory *factory = GetFactory();
208 fTowerHits.clear();
209 fTowerFractions.clear();
210 fTrackFractions.clear();
211
212 // loop over all particles
213 fItParticleInputArray->Reset();
214 number = -1;
215 while((particle = static_cast<Candidate*>(fItParticleInputArray->Next())))
216 {
217 const TLorentzVector &particlePosition = particle->Position;
218 ++number;
219
220 pdgCode = TMath::Abs(particle->PID);
221
222 itFractionMap = fFractionMap.find(pdgCode);
223 if(itFractionMap == fFractionMap.end())
224 {
225 itFractionMap = fFractionMap.find(0);
226 }
227
228 fraction = itFractionMap->second;
229 fTowerFractions.push_back(fraction);
230
231 if(fraction < 1.0E-9) continue;
232
233 // find eta bin [1, fEtaBins.size - 1]
234 itEtaBin = lower_bound(fEtaBins.begin(), fEtaBins.end(), particlePosition.Eta());
235 if(itEtaBin == fEtaBins.begin() || itEtaBin == fEtaBins.end()) continue;
236 etaBin = distance(fEtaBins.begin(), itEtaBin);
237
238 // phi bins for given eta bin
239 phiBins = fPhiBins[etaBin];
240
241 // find phi bin [1, phiBins.size - 1]
242 itPhiBin = lower_bound(phiBins->begin(), phiBins->end(), particlePosition.Phi());
243 if(itPhiBin == phiBins->begin() || itPhiBin == phiBins->end()) continue;
244 phiBin = distance(phiBins->begin(), itPhiBin);
245
246 flags = 0;
247 flags |= (pdgCode == 11 || pdgCode == 22) << 1;
248
249 // make tower hit {16-bits for eta bin number, 16-bits for phi bin number, 8-bits for flags, 24-bits for particle number}
250 towerHit = (Long64_t(etaBin) << 48) | (Long64_t(phiBin) << 32) | (Long64_t(flags) << 24) | Long64_t(number);
251
252 fTowerHits.push_back(towerHit);
253 }
254
255 // loop over all tracks
256 fItTrackInputArray->Reset();
257 number = -1;
258 while((track = static_cast<Candidate*>(fItTrackInputArray->Next())))
259 {
260 const TLorentzVector &trackPosition = track->Position;
261 ++number;
262
263 pdgCode = TMath::Abs(track->PID);
264
265 itFractionMap = fFractionMap.find(pdgCode);
266 if(itFractionMap == fFractionMap.end())
267 {
268 itFractionMap = fFractionMap.find(0);
269 }
270
271 fraction = itFractionMap->second;
272
273 fTrackFractions.push_back(fraction);
274
275 // find eta bin [1, fEtaBins.size - 1]
276 itEtaBin = lower_bound(fEtaBins.begin(), fEtaBins.end(), trackPosition.Eta());
277 if(itEtaBin == fEtaBins.begin() || itEtaBin == fEtaBins.end()) continue;
278 etaBin = distance(fEtaBins.begin(), itEtaBin);
279
280 // phi bins for given eta bin
281 phiBins = fPhiBins[etaBin];
282
283 // find phi bin [1, phiBins.size - 1]
284 itPhiBin = lower_bound(phiBins->begin(), phiBins->end(), trackPosition.Phi());
285 if(itPhiBin == phiBins->begin() || itPhiBin == phiBins->end()) continue;
286 phiBin = distance(phiBins->begin(), itPhiBin);
287
288 flags = 1;
289
290 // make tower hit {16-bits for eta bin number, 16-bits for phi bin number, 8-bits for flags, 24-bits for track number}
291 towerHit = (Long64_t(etaBin) << 48) | (Long64_t(phiBin) << 32) | (Long64_t(flags) << 24) | Long64_t(number);
292
293 fTowerHits.push_back(towerHit);
294 }
295
296 // all hits are sorted first by eta bin number, then by phi bin number,
297 // then by flags and then by particle or track number
298 sort(fTowerHits.begin(), fTowerHits.end());
299
300 // loop over all hits
301 towerEtaPhi = 0;
302 fTower = 0;
303 for(itTowerHits = fTowerHits.begin(); itTowerHits != fTowerHits.end(); ++itTowerHits)
304 {
305 towerHit = (*itTowerHits);
306 flags = (towerHit >> 24) & 0x00000000000000FFLL;
307 number = (towerHit) & 0x0000000000FFFFFFLL;
308 hitEtaPhi = towerHit >> 32;
309
310 if(towerEtaPhi != hitEtaPhi)
311 {
312 // switch to next tower
313 towerEtaPhi = hitEtaPhi;
314
315 // finalize previous tower
316 FinalizeTower();
317
318 // create new tower
319 fTower = factory->NewCandidate();
320
321 phiBin = (towerHit >> 32) & 0x000000000000FFFFLL;
322 etaBin = (towerHit >> 48) & 0x000000000000FFFFLL;
323
324 // phi bins for given eta bin
325 phiBins = fPhiBins[etaBin];
326
327 // calculate eta and phi of the tower's center
328 fTowerEta = 0.5*(fEtaBins[etaBin - 1] + fEtaBins[etaBin]);
329 fTowerPhi = 0.5*((*phiBins)[phiBin - 1] + (*phiBins)[phiBin]);
330
331 fTowerEdges[0] = fEtaBins[etaBin - 1];
332 fTowerEdges[1] = fEtaBins[etaBin];
333 fTowerEdges[2] = (*phiBins)[phiBin - 1];
334 fTowerEdges[3] = (*phiBins)[phiBin];
335
336 fTowerEnergy = 0.0;
337 fTrackEnergy = 0.0;
338
339 fTowerTime = 0.0;
340 fTrackTime = 0.0;
341
342 fTowerTimeWeight = 0.0;
343
344 fTowerTrackHits = 0;
345 fTowerPhotonHits = 0;
346
347 fTowerTrackArray->Clear();
348 }
349
350 // check for track hits
351 if(flags & 1)
352 {
353 ++fTowerTrackHits;
354
355 track = static_cast<Candidate*>(fTrackInputArray->At(number));
356 momentum = track->Momentum;
357 position = track->Position;
358
359 energy = momentum.E() * fTrackFractions[number];
360
361 fTrackEnergy += energy;
362
363 fTrackTime += TMath::Sqrt(energy)*position.T();
364 fTrackTimeWeight += TMath::Sqrt(energy);
365
366 fTowerTrackArray->Add(track);
367
368 continue;
369 }
370
371 // check for photon and electron hits in current tower
372 if(flags & 2) ++fTowerPhotonHits;
373
374 particle = static_cast<Candidate*>(fParticleInputArray->At(number));
375 momentum = particle->Momentum;
376 position = particle->Position;
377
378 // fill current tower
379 energy = momentum.E() * fTowerFractions[number];
380
381 fTowerEnergy += energy;
382
383 fTowerTime += TMath::Sqrt(energy)*position.T();
384 fTowerTimeWeight += TMath::Sqrt(energy);
385
386 fTower->AddCandidate(particle);
387 }
388
389 // finalize last tower
390 FinalizeTower();
391}
392
393//------------------------------------------------------------------------------
394
395void SimpleCalorimeter::FinalizeTower()
396{
397 Candidate *tower;
398 Double_t energy, pt, eta, phi;
399 Double_t sigma;
400 Double_t time;
401
402 if(!fTower) return;
403
404 sigma = fResolutionFormula->Eval(0.0, fTowerEta, 0.0, fTowerEnergy);
405
406 energy = LogNormal(fTowerEnergy, sigma);
407
408 time = (fTowerTimeWeight < 1.0E-09 ) ? 0.0 : fTowerTime/fTowerTimeWeight;
409
410 sigma = fResolutionFormula->Eval(0.0, fTowerEta, 0.0, energy);
411
412 if(energy < fEnergyMin || energy < fEnergySignificanceMin*sigma) energy = 0.0;
413
414 if(fSmearTowerCenter)
415 {
416 eta = gRandom->Uniform(fTowerEdges[0], fTowerEdges[1]);
417 phi = gRandom->Uniform(fTowerEdges[2], fTowerEdges[3]);
418 }
419 else
420 {
421 eta = fTowerEta;
422 phi = fTowerPhi;
423 }
424
425 pt = energy / TMath::CosH(eta);
426
427 fTower->Position.SetPtEtaPhiE(1.0, eta, phi, time);
428 fTower->Momentum.SetPtEtaPhiE(pt, eta, phi, energy);
429
430 fTower->Eem = (!fIsEcal) ? 0 : energy;
431 fTower->Ehad = (fIsEcal) ? 0 : energy;
432
433 fTower->Edges[0] = fTowerEdges[0];
434 fTower->Edges[1] = fTowerEdges[1];
435 fTower->Edges[2] = fTowerEdges[2];
436 fTower->Edges[3] = fTowerEdges[3];
437
438 // fill SimpleCalorimeter towers
439 if(energy > 0.0) fTowerOutputArray->Add(fTower);
440
441 // fill energy flow candidates
442 energy -= fTrackEnergy;
443
444 sigma = fResolutionFormula->Eval(0.0, fTowerEta, 0.0, energy);
445
446 if(energy < fEnergyMin || energy < fEnergySignificanceMin*sigma) energy = 0.0;
447
448 // save energy excess as an energy flow tower
449 if(energy > 0.0)
450 {
451 // create new photon tower
452 tower = static_cast<Candidate*>(fTower->Clone());
453 pt = energy / TMath::CosH(eta);
454
455 tower->Eem = (!fIsEcal) ? 0 : energy;
456 tower->Ehad = (fIsEcal) ? 0 : energy;
457
458 tower->Momentum.SetPtEtaPhiE(pt, eta, phi, energy);
459 fEFlowTowerOutputArray->Add(tower);
460 }
461}
462
463//------------------------------------------------------------------------------
464
465Double_t SimpleCalorimeter::LogNormal(Double_t mean, Double_t sigma)
466{
467 Double_t a, b;
468
469 if(mean > 0.0)
470 {
471 b = TMath::Sqrt(TMath::Log((1.0 + (sigma*sigma)/(mean*mean))));
472 a = TMath::Log(mean) - 0.5*b*b;
473
474 return TMath::Exp(a + b*gRandom->Gaus(0.0, 1.0));
475 }
476 else
477 {
478 return 0.0;
479 }
480}
Note: See TracBrowser for help on using the repository browser.