Fork me on GitHub

source: git/modules/Calorimeter.cc@ 4b9a2dc

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

added EnergyMin and SigmaMin parameters

  • Property mode set to 100644
File size: 16.2 KB
Line 
1
2/** \class Calorimeter
3 *
4 * Fills calorimeter towers, performs calorimeter resolution smearing,
5 * and creates energy flow objects (tracks, photons, and neutral hadrons).
6 *
7 * $Date$
8 * $Revision$
9 *
10 *
11 * \author P. Demin - UCL, Louvain-la-Neuve
12 *
13 */
14
15#include "modules/Calorimeter.h"
16
17#include "classes/DelphesClasses.h"
18#include "classes/DelphesFactory.h"
19#include "classes/DelphesFormula.h"
20
21#include "ExRootAnalysis/ExRootResult.h"
22#include "ExRootAnalysis/ExRootFilter.h"
23#include "ExRootAnalysis/ExRootClassifier.h"
24
25#include "TMath.h"
26#include "TString.h"
27#include "TFormula.h"
28#include "TRandom3.h"
29#include "TObjArray.h"
30#include "TDatabasePDG.h"
31#include "TLorentzVector.h"
32
33#include <algorithm>
34#include <stdexcept>
35#include <iostream>
36#include <sstream>
37
38using namespace std;
39
40//------------------------------------------------------------------------------
41
42Calorimeter::Calorimeter() :
43 fECalResolutionFormula(0), fHCalResolutionFormula(0),
44 fItParticleInputArray(0), fItTrackInputArray(0),
45 fTowerTrackArray(0), fItTowerTrackArray(0)
46{
47 fECalResolutionFormula = new DelphesFormula;
48 fHCalResolutionFormula = new DelphesFormula;
49
50 fTowerTrackArray = new TObjArray;
51 fItTowerTrackArray = fTowerTrackArray->MakeIterator();
52}
53
54//------------------------------------------------------------------------------
55
56Calorimeter::~Calorimeter()
57{
58 if(fECalResolutionFormula) delete fECalResolutionFormula;
59 if(fHCalResolutionFormula) delete fHCalResolutionFormula;
60
61 if(fTowerTrackArray) delete fTowerTrackArray;
62 if(fItTowerTrackArray) delete fItTowerTrackArray;
63}
64
65//------------------------------------------------------------------------------
66
67void Calorimeter::Init()
68{
69 ExRootConfParam param, paramEtaBins, paramPhiBins, paramFractions;
70 Long_t i, j, k, size, sizeEtaBins, sizePhiBins, sizeFractions;
71 Double_t ecalFraction, hcalFraction;
72 TBinMap::iterator itEtaBin;
73 set< Double_t >::iterator itPhiBin;
74 vector< Double_t > *phiBins;
75
76 // read eta and phi bins
77 param = GetParam("EtaPhiBins");
78 size = param.GetSize();
79 fBinMap.clear();
80 fEtaBins.clear();
81 fPhiBins.clear();
82 for(i = 0; i < size/2; ++i)
83 {
84 paramEtaBins = param[i*2];
85 sizeEtaBins = paramEtaBins.GetSize();
86 paramPhiBins = param[i*2 + 1];
87 sizePhiBins = paramPhiBins.GetSize();
88
89 for(j = 0; j < sizeEtaBins; ++j)
90 {
91 for(k = 0; k < sizePhiBins; ++k)
92 {
93 fBinMap[paramEtaBins[j].GetDouble()].insert(paramPhiBins[k].GetDouble());
94 }
95 }
96 }
97
98 // for better performance we transform map of sets to parallel vectors:
99 // vector< double > and vector< vector< double >* >
100 for(itEtaBin = fBinMap.begin(); itEtaBin != fBinMap.end(); ++itEtaBin)
101 {
102 fEtaBins.push_back(itEtaBin->first);
103 phiBins = new vector< double >(itEtaBin->second.size());
104 fPhiBins.push_back(phiBins);
105 phiBins->clear();
106 for(itPhiBin = itEtaBin->second.begin(); itPhiBin != itEtaBin->second.end(); ++itPhiBin)
107 {
108 phiBins->push_back(*itPhiBin);
109 }
110 }
111
112 // read energy fractions for different particles
113 param = GetParam("EnergyFraction");
114 size = param.GetSize();
115
116 // set default energy fractions values
117 fFractionMap.clear();
118 fFractionMap[0] = make_pair(0.0, 1.0);
119
120 for(i = 0; i < size/2; ++i)
121 {
122 paramFractions = param[i*2 + 1];
123 sizeFractions = paramFractions.GetSize();
124
125 ecalFraction = paramFractions[0].GetDouble();
126 hcalFraction = paramFractions[1].GetDouble();
127
128 fFractionMap[param[i*2].GetInt()] = make_pair(ecalFraction, hcalFraction);
129 }
130/*
131 TFractionMap::iterator itFractionMap;
132 for(itFractionMap = fFractionMap.begin(); itFractionMap != fFractionMap.end(); ++itFractionMap)
133 {
134 cout << itFractionMap->first << " " << itFractionMap->second.first << " " << itFractionMap->second.second << endl;
135 }
136*/
137
138 // read min E value for towers to be saved
139 fEcalEnergyMin = GetDouble("EcalTowerMinEnergy", 0.0);
140 fHcalEnergyMin = GetDouble("HcalTowerMinEnergy", 0.0);
141
142 fEcalSigmaMin = GetDouble("EcalTowerMinSignificance", 0.0);
143 fHcalSigmaMin = GetDouble("HcalTowerMinSignificance", 0.0);
144
145
146 // read resolution formulas
147 fECalResolutionFormula->Compile(GetString("ECalResolutionFormula", "0"));
148 fHCalResolutionFormula->Compile(GetString("HCalResolutionFormula", "0"));
149
150 // import array with output from other modules
151 fParticleInputArray = ImportArray(GetString("ParticleInputArray", "ParticlePropagator/particles"));
152 fItParticleInputArray = fParticleInputArray->MakeIterator();
153
154 fTrackInputArray = ImportArray(GetString("TrackInputArray", "ParticlePropagator/tracks"));
155 fItTrackInputArray = fTrackInputArray->MakeIterator();
156
157 // create output arrays
158 fTowerOutputArray = ExportArray(GetString("TowerOutputArray", "towers"));
159 fPhotonOutputArray = ExportArray(GetString("PhotonOutputArray", "photons"));
160
161 fEFlowTrackOutputArray = ExportArray(GetString("EFlowTrackOutputArray", "eflowTracks"));
162 fEFlowPhotonOutputArray = ExportArray(GetString("EFlowPhotonOutputArray", "eflowPhotons"));
163 fEFlowNeutralHadronOutputArray = ExportArray(GetString("EFlowNeutralHadronOutputArray", "eflowNeutralHadrons"));
164
165
166}
167
168//------------------------------------------------------------------------------
169
170void Calorimeter::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 Calorimeter::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 ecalFraction, hcalFraction;
191 Double_t ecalEnergy, hcalEnergy;
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 fTowerECalFractions.clear();
205 fTowerHCalFractions.clear();
206 fTrackECalFractions.clear();
207 fTrackHCalFractions.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 ecalFraction = itFractionMap->second.first;
226 hcalFraction = itFractionMap->second.second;
227
228 fTowerECalFractions.push_back(ecalFraction);
229 fTowerHCalFractions.push_back(hcalFraction);
230
231 if(ecalFraction < 1.0E-9 && hcalFraction < 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 ecalFraction = itFractionMap->second.first;
272 hcalFraction = itFractionMap->second.second;
273
274 fTrackECalFractions.push_back(ecalFraction);
275 fTrackHCalFractions.push_back(hcalFraction);
276
277 // find eta bin [1, fEtaBins.size - 1]
278 itEtaBin = lower_bound(fEtaBins.begin(), fEtaBins.end(), trackPosition.Eta());
279 if(itEtaBin == fEtaBins.begin() || itEtaBin == fEtaBins.end()) continue;
280 etaBin = distance(fEtaBins.begin(), itEtaBin);
281
282 // phi bins for given eta bin
283 phiBins = fPhiBins[etaBin];
284
285 // find phi bin [1, phiBins.size - 1]
286 itPhiBin = lower_bound(phiBins->begin(), phiBins->end(), trackPosition.Phi());
287 if(itPhiBin == phiBins->begin() || itPhiBin == phiBins->end()) continue;
288 phiBin = distance(phiBins->begin(), itPhiBin);
289
290 flags = 1;
291
292 // make tower hit {16-bits for eta bin number, 16-bits for phi bin number, 8-bits for flags, 24-bits for track number}
293 towerHit = (Long64_t(etaBin) << 48) | (Long64_t(phiBin) << 32) | (Long64_t(flags) << 24) | Long64_t(number);
294
295 fTowerHits.push_back(towerHit);
296 }
297
298 // all hits are sorted first by eta bin number, then by phi bin number,
299 // then by flags and then by particle or track number
300 sort(fTowerHits.begin(), fTowerHits.end());
301
302 // loop over all hits
303 towerEtaPhi = 0;
304 fTower = 0;
305 for(itTowerHits = fTowerHits.begin(); itTowerHits != fTowerHits.end(); ++itTowerHits)
306 {
307 towerHit = (*itTowerHits);
308 flags = (towerHit >> 24) & 0x00000000000000FFLL;
309 number = (towerHit) & 0x0000000000FFFFFFLL;
310 hitEtaPhi = towerHit >> 32;
311
312 if(towerEtaPhi != hitEtaPhi)
313 {
314 // switch to next tower
315 towerEtaPhi = hitEtaPhi;
316
317 // finalize previous tower
318 FinalizeTower();
319
320 // create new tower
321 fTower = factory->NewCandidate();
322
323 phiBin = (towerHit >> 32) & 0x000000000000FFFFLL;
324 etaBin = (towerHit >> 48) & 0x000000000000FFFFLL;
325
326 // phi bins for given eta bin
327 phiBins = fPhiBins[etaBin];
328
329 // calculate eta and phi of the tower's center
330 fTowerEta = 0.5*(fEtaBins[etaBin - 1] + fEtaBins[etaBin]);
331 fTowerPhi = 0.5*((*phiBins)[phiBin - 1] + (*phiBins)[phiBin]);
332
333 fTowerEdges[0] = fEtaBins[etaBin - 1];
334 fTowerEdges[1] = fEtaBins[etaBin];
335 fTowerEdges[2] = (*phiBins)[phiBin - 1];
336 fTowerEdges[3] = (*phiBins)[phiBin];
337
338 fTowerECalEnergy = 0.0;
339 fTowerHCalEnergy = 0.0;
340
341 fTrackECalEnergy = 0.0;
342 fTrackHCalEnergy = 0.0;
343
344 fTowerECalTime = 0.0;
345 fTowerHCalTime = 0.0;
346
347 fTrackECalTime = 0.0;
348 fTrackHCalTime = 0.0;
349
350 fTowerECalWeightTime = 0.0;
351 fTowerHCalWeightTime = 0.0;
352
353 fTowerTrackHits = 0;
354 fTowerPhotonHits = 0;
355
356 fTowerTrackArray->Clear();
357 }
358
359 // check for track hits
360 if(flags & 1)
361 {
362 ++fTowerTrackHits;
363
364 track = static_cast<Candidate*>(fTrackInputArray->At(number));
365 momentum = track->Momentum;
366 position = track->Position;
367
368
369 ecalEnergy = momentum.E() * fTrackECalFractions[number];
370 hcalEnergy = momentum.E() * fTrackHCalFractions[number];
371
372 fTrackECalEnergy += ecalEnergy;
373 fTrackHCalEnergy += hcalEnergy;
374
375 fTrackECalTime += TMath::Sqrt(ecalEnergy)*position.T();
376 fTrackHCalTime += TMath::Sqrt(hcalEnergy)*position.T();
377
378 fTrackECalWeightTime += TMath::Sqrt(ecalEnergy);
379 fTrackHCalWeightTime += TMath::Sqrt(hcalEnergy);
380
381 fTowerTrackArray->Add(track);
382
383 continue;
384 }
385
386 // check for photon and electron hits in current tower
387 if(flags & 2) ++fTowerPhotonHits;
388
389 particle = static_cast<Candidate*>(fParticleInputArray->At(number));
390 momentum = particle->Momentum;
391 position = particle->Position;
392
393 // fill current tower
394 ecalEnergy = momentum.E() * fTowerECalFractions[number];
395 hcalEnergy = momentum.E() * fTowerHCalFractions[number];
396
397 fTowerECalEnergy += ecalEnergy;
398 fTowerHCalEnergy += hcalEnergy;
399
400 fTowerECalTime += TMath::Sqrt(ecalEnergy)*position.T();
401 fTowerHCalTime += TMath::Sqrt(hcalEnergy)*position.T();
402
403 fTowerECalWeightTime += TMath::Sqrt(ecalEnergy);
404 fTowerHCalWeightTime += TMath::Sqrt(hcalEnergy);
405
406
407 fTower->AddCandidate(particle);
408 }
409
410 // finalize last tower
411 FinalizeTower();
412}
413
414//------------------------------------------------------------------------------
415
416void Calorimeter::FinalizeTower()
417{
418 Candidate *track, *tower;
419 Double_t energy, pt, eta, phi;
420 Double_t ecalEnergy, hcalEnergy;
421 Double_t ecalSigma, hcalSigma;
422 Double_t ecalTime, hcalTime, time;
423
424 if(!fTower) return;
425// cout<<"----------------------"<<endl;
426// cout<<"Finalize Tower"<<endl;
427// cout<<""<<endl;
428
429
430 ecalSigma = fECalResolutionFormula->Eval(0.0, fTowerEta, 0.0, fTowerECalEnergy);
431
432// ecalEnergy = gRandom->Gaus(fTowerECalEnergy, ecalSigma);
433// if(ecalEnergy < 0.0) ecalEnergy = 0.0;
434
435 ecalEnergy = LogNormal(fTowerECalEnergy, ecalSigma);
436 ecalTime = (fTowerECalWeightTime < 1.0E-09 ) ? 0 : fTowerECalTime/fTowerECalWeightTime;
437
438 hcalSigma = fHCalResolutionFormula->Eval(0.0, fTowerEta, 0.0, fTowerHCalEnergy);
439
440// hcalEnergy = gRandom->Gaus(fTowerHCalEnergy, hcalSigma);
441// if(hcalEnergy < 0.0) hcalEnergy = 0.0;
442
443 hcalEnergy = LogNormal(fTowerHCalEnergy, hcalSigma);
444 hcalTime = (fTowerHCalWeightTime < 1.0E-09 ) ? 0 : fTowerHCalTime/fTowerHCalWeightTime;
445
446
447 ecalSigma = fECalResolutionFormula->Eval(0.0, fTowerEta, 0.0, ecalEnergy);
448 hcalSigma = fHCalResolutionFormula->Eval(0.0, fTowerEta, 0.0, hcalEnergy);
449
450 ecalEnergy = (ecalEnergy < fEcalEnergyMin || ecalEnergy < fEcalSigmaMin*ecalSigma) ? 0 : ecalEnergy;
451 hcalEnergy = (hcalEnergy < fHcalEnergyMin || hcalEnergy < fHcalSigmaMin*hcalSigma) ? 0 : hcalEnergy;
452
453 energy = ecalEnergy + hcalEnergy;
454 time = (TMath::Sqrt(ecalEnergy)*ecalTime + TMath::Sqrt(hcalEnergy)*hcalTime)/(TMath::Sqrt(ecalEnergy) + TMath::Sqrt(hcalEnergy));
455
456// eta = fTowerEta;
457// phi = fTowerPhi;
458
459 eta = gRandom->Uniform(fTowerEdges[0], fTowerEdges[1]);
460 phi = gRandom->Uniform(fTowerEdges[2], fTowerEdges[3]);
461
462 pt = energy / TMath::CosH(eta);
463
464 // fTower->Position.SetXYZT(-time, 0.0, 0.0, time);
465 fTower->Position.SetPtEtaPhiE(1.0, eta, phi, time);
466 fTower->Momentum.SetPtEtaPhiE(pt, eta, phi, energy);
467 fTower->Eem = ecalEnergy;
468 fTower->Ehad = hcalEnergy;
469
470 fTower->Edges[0] = fTowerEdges[0];
471 fTower->Edges[1] = fTowerEdges[1];
472 fTower->Edges[2] = fTowerEdges[2];
473 fTower->Edges[3] = fTowerEdges[3];
474
475 if( energy > 0.0 )
476 {
477 if(fTowerPhotonHits > 0 && fTowerTrackHits == 0)
478 {
479 fPhotonOutputArray->Add(fTower);
480 }
481
482 fTowerOutputArray->Add(fTower);
483 }
484
485 // fill energy flow candidates
486
487 // save all the tracks as energy flow tracks
488 fItTowerTrackArray->Reset();
489 while((track = static_cast<Candidate*>(fItTowerTrackArray->Next())))
490 {
491 fEFlowTrackOutputArray->Add(track);
492 }
493
494 ecalEnergy -= fTrackECalEnergy;
495 if(ecalEnergy < fEcalEnergyMin || ecalEnergy < fEcalSigmaMin*fECalResolutionFormula->Eval(0.0, fTowerEta, 0.0, ecalEnergy)) ecalEnergy = 0.0;
496
497 hcalEnergy -= fTrackHCalEnergy;
498 if(hcalEnergy < fHcalEnergyMin || hcalEnergy < fHcalSigmaMin*fHCalResolutionFormula->Eval(0.0, fTowerEta, 0.0, hcalEnergy)) hcalEnergy = 0.0;
499
500 energy = ecalEnergy + hcalEnergy;
501
502 if(ecalEnergy > 0.0)
503 {
504 // create new photon tower
505 tower = static_cast<Candidate*>(fTower->Clone());
506
507 pt = ecalEnergy / TMath::CosH(eta);
508
509 tower->Momentum.SetPtEtaPhiE(pt, eta, phi, ecalEnergy);
510 tower->Eem = ecalEnergy;
511 tower->Ehad = 0;
512
513 fEFlowPhotonOutputArray->Add(tower);
514 }
515 if(hcalEnergy > 0.0)
516 {
517 // create new neutral hadron tower
518 tower = static_cast<Candidate*>(fTower->Clone());
519
520 pt = hcalEnergy / TMath::CosH(eta);
521
522 tower->Momentum.SetPtEtaPhiE(pt, eta, phi, hcalEnergy);
523 tower->Eem = 0;
524 tower->Ehad = hcalEnergy;
525
526 fEFlowNeutralHadronOutputArray->Add(tower);
527 }
528
529
530
531
532}
533
534//------------------------------------------------------------------------------
535
536Double_t Calorimeter::LogNormal(Double_t mean, Double_t sigma)
537{
538 Double_t a, b;
539
540 if(mean > 0.0)
541 {
542 b = TMath::Sqrt(TMath::Log((1.0 + (sigma*sigma)/(mean*mean))));
543 a = TMath::Log(mean) - 0.5*b*b;
544
545 return TMath::Exp(a + b*gRandom->Gaus(0, 1));
546 }
547 else
548 {
549 return 0.0;
550 }
551}
Note: See TracBrowser for help on using the repository browser.