Fork me on GitHub

source: git/modules/SimpleCalorimeter.cc@ 1df571a

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

prepend GPLv3 header

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