Fork me on GitHub

source: git/modules/SimpleCalorimeter.cc@ 734b267

ImprovedOutputFile Timing dual_readout llp
Last change on this file since 734b267 was e0f8f99, checked in by Michele Selvaggi <michele.selvaggi@…>, 8 years ago

put back energy weighted time mean

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