Fork me on GitHub

source: git/modules/SimpleCalorimeter.cc@ ebf40fd

Last change on this file since ebf40fd was eee94204, checked in by michele <michele.selvaggi@…>, 3 years ago

fixed bug introduced on calo tower postion

  • Property mode set to 100644
File size: 16.1 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/** \class SimpleCalorimeter
20 *
21 * Fills SimpleCalorimeter towers, performs SimpleCalorimeter resolution smearing,
22 * and creates energy flow objects (tracks, photons, and neutral hadrons).
23 *
24 * \author P. Demin - UCL, Louvain-la-Neuve
25 *
26 */
27
28#include "modules/SimpleCalorimeter.h"
29
30#include "classes/DelphesClasses.h"
31#include "classes/DelphesFactory.h"
32#include "classes/DelphesFormula.h"
33
34#include "ExRootAnalysis/ExRootClassifier.h"
35#include "ExRootAnalysis/ExRootFilter.h"
36#include "ExRootAnalysis/ExRootResult.h"
37
38#include "TDatabasePDG.h"
39#include "TFormula.h"
40#include "TLorentzVector.h"
41#include "TMath.h"
42#include "TObjArray.h"
43#include "TRandom3.h"
44#include "TString.h"
45
46#include <algorithm>
47#include <iostream>
48#include <sstream>
49#include <stdexcept>
50
51using namespace std;
52
53//------------------------------------------------------------------------------
54
55SimpleCalorimeter::SimpleCalorimeter() :
56 fResolutionFormula(0),
57 fItParticleInputArray(0), fItTrackInputArray(0)
58{
59
60 fResolutionFormula = new DelphesFormula;
61 fTowerTrackArray = new TObjArray;
62 fItTowerTrackArray = fTowerTrackArray->MakeIterator();
63}
64
65//------------------------------------------------------------------------------
66
67SimpleCalorimeter::~SimpleCalorimeter()
68{
69
70 if(fResolutionFormula) delete fResolutionFormula;
71 if(fTowerTrackArray) delete fTowerTrackArray;
72 if(fItTowerTrackArray) delete fItTowerTrackArray;
73}
74
75//------------------------------------------------------------------------------
76
77void SimpleCalorimeter::Init()
78{
79 ExRootConfParam param, paramEtaBins, paramPhiBins, paramFractions;
80 Long_t i, j, k, size, sizeEtaBins, sizePhiBins;
81 Double_t fraction;
82 TBinMap::iterator itEtaBin;
83 set<Double_t>::iterator itPhiBin;
84 vector<Double_t> *phiBins;
85
86 // read eta and phi bins
87 param = GetParam("EtaPhiBins");
88 size = param.GetSize();
89 fBinMap.clear();
90 fEtaBins.clear();
91 fPhiBins.clear();
92 for(i = 0; i < size / 2; ++i)
93 {
94 paramEtaBins = param[i * 2];
95 sizeEtaBins = paramEtaBins.GetSize();
96 paramPhiBins = param[i * 2 + 1];
97 sizePhiBins = paramPhiBins.GetSize();
98
99 for(j = 0; j < sizeEtaBins; ++j)
100 {
101 for(k = 0; k < sizePhiBins; ++k)
102 {
103 fBinMap[paramEtaBins[j].GetDouble()].insert(paramPhiBins[k].GetDouble());
104 }
105 }
106 }
107
108 // for better performance we transform map of sets to parallel vectors:
109 // vector< double > and vector< vector< double >* >
110 for(itEtaBin = fBinMap.begin(); itEtaBin != fBinMap.end(); ++itEtaBin)
111 {
112 fEtaBins.push_back(itEtaBin->first);
113 phiBins = new vector<double>(itEtaBin->second.size());
114 fPhiBins.push_back(phiBins);
115 phiBins->clear();
116 for(itPhiBin = itEtaBin->second.begin(); itPhiBin != itEtaBin->second.end(); ++itPhiBin)
117 {
118 phiBins->push_back(*itPhiBin);
119 }
120 }
121
122 // read energy fractions for different particles
123 param = GetParam("EnergyFraction");
124 size = param.GetSize();
125
126 // set default energy fractions values
127 fFractionMap.clear();
128 fFractionMap[0] = 1.0;
129
130 for(i = 0; i < size / 2; ++i)
131 {
132 paramFractions = param[i * 2 + 1];
133 fraction = paramFractions[0].GetDouble();
134 fFractionMap[param[i * 2].GetInt()] = fraction;
135 }
136
137 // read min E value for towers to be saved
138 fEnergyMin = GetDouble("EnergyMin", 0.0);
139
140 fEnergySignificanceMin = GetDouble("EnergySignificanceMin", 0.0);
141
142 // flag that says if current calo is Ecal of Hcal (will then fill correct values of Eem and Ehad)
143 fIsEcal = GetBool("IsEcal", false);
144
145 // switch on or off the dithering of the center of calorimeter towers
146 fSmearTowerCenter = GetBool("SmearTowerCenter", true);
147
148 // read resolution formulas
149 fResolutionFormula->Compile(GetString("ResolutionFormula", "0"));
150
151 // import array with output from other modules
152 fParticleInputArray = ImportArray(GetString("ParticleInputArray", "ParticlePropagator/particles"));
153 fItParticleInputArray = fParticleInputArray->MakeIterator();
154
155 fTrackInputArray = ImportArray(GetString("TrackInputArray", "ParticlePropagator/tracks"));
156 fItTrackInputArray = fTrackInputArray->MakeIterator();
157
158 // create output arrays
159 fTowerOutputArray = ExportArray(GetString("TowerOutputArray", "towers"));
160
161 fEFlowTrackOutputArray = ExportArray(GetString("EFlowTrackOutputArray", "eflowTracks"));
162 fEFlowTowerOutputArray = ExportArray(GetString("EFlowTowerOutputArray", "eflowTowers"));
163}
164
165//------------------------------------------------------------------------------
166
167void SimpleCalorimeter::Finish()
168{
169 vector<vector<Double_t> *>::iterator itPhiBin;
170 if(fItParticleInputArray) delete fItParticleInputArray;
171 if(fItTrackInputArray) delete fItTrackInputArray;
172 for(itPhiBin = fPhiBins.begin(); itPhiBin != fPhiBins.end(); ++itPhiBin)
173 {
174 delete *itPhiBin;
175 }
176}
177
178//------------------------------------------------------------------------------
179
180void SimpleCalorimeter::Process()
181{
182 Candidate *particle, *track;
183 TLorentzVector position, momentum;
184 Short_t etaBin, phiBin, flags;
185 Int_t number;
186 Long64_t towerHit, towerEtaPhi, hitEtaPhi;
187 Double_t fraction;
188 Double_t energy;
189 Double_t sigma;
190 Double_t energyGuess;
191
192 Int_t pdgCode;
193
194 TFractionMap::iterator itFractionMap;
195
196 vector<Double_t>::iterator itEtaBin;
197 vector<Double_t>::iterator itPhiBin;
198 vector<Double_t> *phiBins;
199
200 vector<Long64_t>::iterator itTowerHits;
201
202 DelphesFactory *factory = GetFactory();
203 fTowerHits.clear();
204 fTowerFractions.clear();
205 fTrackFractions.clear();
206
207 // loop over all particles
208 fItParticleInputArray->Reset();
209 number = -1;
210 fTowerRmax=0.;
211 while((particle = static_cast<Candidate *>(fItParticleInputArray->Next())))
212 {
213 const TLorentzVector &particlePosition = particle->Position;
214 ++number;
215
216 // compute maximum radius (needed in FinalizeTower to assess whether barrel or endcap tower)
217 if (particlePosition.Perp() > fTowerRmax)
218 fTowerRmax=particlePosition.Perp();
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
338 fTrackEnergy = 0.0;
339 fTrackSigma = 0.0;
340
341 fTowerTime = 0.0;
342 fTrackTime = 0.0;
343
344 fTowerTimeWeight = 0.0;
345
346 fTowerTrackHits = 0;
347 fTowerPhotonHits = 0;
348
349 fTowerTrackArray->Clear();
350 }
351
352 // check for track hits
353 if(flags & 1)
354 {
355 ++fTowerTrackHits;
356
357 track = static_cast<Candidate *>(fTrackInputArray->At(number));
358 momentum = track->Momentum;
359 position = track->Position;
360
361 energy = momentum.E() * fTrackFractions[number];
362
363 fTrackTime += TMath::Sqrt(energy) * position.T();
364 fTrackTimeWeight += TMath::Sqrt(energy);
365
366 if(fTrackFractions[number] > 1.0E-9)
367 {
368
369 // compute total charged energy
370 fTrackEnergy += energy;
371 sigma = fResolutionFormula->Eval(0.0, fTowerEta, 0.0, momentum.E());
372 if(sigma / momentum.E() < track->TrackResolution)
373 energyGuess = energy;
374 else
375 energyGuess = momentum.E();
376
377 fTrackSigma += ((track->TrackResolution) * energyGuess) * ((track->TrackResolution) * energyGuess);
378 fTowerTrackArray->Add(track);
379 }
380
381 else
382 {
383 fEFlowTrackOutputArray->Add(track);
384 }
385
386 continue;
387 }
388
389 // check for photon and electron hits in current tower
390 if(flags & 2) ++fTowerPhotonHits;
391
392 particle = static_cast<Candidate *>(fParticleInputArray->At(number));
393 momentum = particle->Momentum;
394 position = particle->Position;
395
396 // fill current tower
397 energy = momentum.E() * fTowerFractions[number];
398
399 fTowerEnergy += energy;
400
401 fTowerTime += energy * energy * position.T(); //sigma_t ~ 1/E
402 fTowerTimeWeight += energy * energy;
403
404 fTower->AddCandidate(particle);
405 fTower->Position = position;
406
407 }
408
409 // finalize last tower
410 FinalizeTower();
411}
412
413//------------------------------------------------------------------------------
414
415void SimpleCalorimeter::FinalizeTower()
416{
417 Candidate *tower, *track, *mother;
418 Double_t energy, neutralEnergy, pt, eta, phi, r;
419 Double_t sigma, neutralSigma;
420 Double_t time;
421
422 Double_t weightTrack, weightCalo, bestEnergyEstimate, rescaleFactor;
423
424 TLorentzVector momentum;
425 TFractionMap::iterator itFractionMap;
426
427 if(!fTower) return;
428
429 sigma = fResolutionFormula->Eval(0.0, fTowerEta, 0.0, fTowerEnergy);
430
431 energy = LogNormal(fTowerEnergy, sigma);
432
433 time = (fTowerTimeWeight < 1.0E-09) ? 0.0 : fTowerTime / fTowerTimeWeight;
434
435 sigma = fResolutionFormula->Eval(0.0, fTowerEta, 0.0, energy);
436
437 if(energy < fEnergyMin || energy < fEnergySignificanceMin * sigma) energy = 0.0;
438
439 if(fSmearTowerCenter)
440 {
441 eta = gRandom->Uniform(fTowerEdges[0], fTowerEdges[1]);
442 phi = gRandom->Uniform(fTowerEdges[2], fTowerEdges[3]);
443 }
444 else
445 {
446 eta = fTowerEta;
447 phi = fTowerPhi;
448 }
449
450 pt = energy / TMath::CosH(eta);
451
452 // endcap
453 if (TMath::Abs(fTower->Position.Pt() - fTowerRmax) > 1.e-06 && TMath::Abs(eta) > 0.){
454 r = fTower->Position.Z()/TMath::SinH(eta);
455 }
456 // barrel
457 else {
458 r = fTower->Position.Pt();
459 }
460
461 fTower->Position.SetPtEtaPhiE(r, eta, phi, time);
462 fTower->Momentum.SetPtEtaPhiE(pt, eta, phi, energy);
463 fTower->L = fTower->Position.Vect().Mag();
464
465 fTower->Eem = (!fIsEcal) ? 0 : energy;
466 fTower->Ehad = (fIsEcal) ? 0 : energy;
467 fTower->Etrk = fTrackEnergy;
468
469 fTower->Edges[0] = fTowerEdges[0];
470 fTower->Edges[1] = fTowerEdges[1];
471 fTower->Edges[2] = fTowerEdges[2];
472 fTower->Edges[3] = fTowerEdges[3];
473
474 // fill SimpleCalorimeter towers
475 if(energy > 0.0) fTowerOutputArray->Add(fTower);
476
477 // e-flow candidates
478
479 //compute neutral excess
480
481 fTrackSigma = TMath::Sqrt(fTrackSigma);
482 neutralEnergy = max((energy - fTrackEnergy), 0.0);
483
484 //compute sigma_trk total
485 neutralSigma = neutralEnergy / TMath::Sqrt(fTrackSigma * fTrackSigma + sigma * sigma);
486
487 // if neutral excess is significant, simply create neutral Eflow tower and clone each track into eflowtrack
488 if(neutralEnergy > fEnergyMin && neutralSigma > fEnergySignificanceMin)
489 {
490 // create new photon tower
491 tower = static_cast<Candidate *>(fTower->Clone());
492 pt = neutralEnergy / TMath::CosH(eta);
493
494 tower->Eem = (!fIsEcal) ? 0 : neutralEnergy;
495 tower->Ehad = (fIsEcal) ? 0 : neutralEnergy;
496 tower->PID = (fIsEcal) ? 22 : 0;
497
498 tower->Momentum.SetPtEtaPhiE(pt, eta, phi, neutralEnergy);
499 fEFlowTowerOutputArray->Add(tower);
500
501 fItTowerTrackArray->Reset();
502 while((track = static_cast<Candidate *>(fItTowerTrackArray->Next())))
503 {
504 mother = track;
505 track = static_cast<Candidate *>(track->Clone());
506 track->AddCandidate(mother);
507
508 fEFlowTrackOutputArray->Add(track);
509 }
510 }
511
512 // if neutral excess is not significant, rescale eflow tracks, such that the total charged equals the best measurement given by the calorimeter and tracking
513 else if(fTrackEnergy > 0.0)
514 {
515 weightTrack = (fTrackSigma > 0.0) ? 1 / (fTrackSigma * fTrackSigma) : 0.0;
516 weightCalo = (sigma > 0.0) ? 1 / (sigma * sigma) : 0.0;
517
518 bestEnergyEstimate = (weightTrack * fTrackEnergy + weightCalo * energy) / (weightTrack + weightCalo);
519 rescaleFactor = bestEnergyEstimate / fTrackEnergy;
520
521 fItTowerTrackArray->Reset();
522 while((track = static_cast<Candidate *>(fItTowerTrackArray->Next())))
523 {
524 mother = track;
525 track = static_cast<Candidate *>(track->Clone());
526 track->AddCandidate(mother);
527 track->Momentum.SetPtEtaPhiM(track->Momentum.Pt()*rescaleFactor, track->Momentum.Eta(), track->Momentum.Phi(), track->Momentum.M());
528 fEFlowTrackOutputArray->Add(track);
529 }
530 }
531}
532
533//------------------------------------------------------------------------------
534
535Double_t SimpleCalorimeter::LogNormal(Double_t mean, Double_t sigma)
536{
537 Double_t a, b;
538
539 if(mean > 0.0)
540 {
541 b = TMath::Sqrt(TMath::Log((1.0 + (sigma * sigma) / (mean * mean))));
542 a = TMath::Log(mean) - 0.5 * b * b;
543
544 return TMath::Exp(a + b * gRandom->Gaus(0.0, 1.0));
545 }
546 else
547 {
548 return 0.0;
549 }
550}
Note: See TracBrowser for help on using the repository browser.