Fork me on GitHub

source: svn/trunk/modules/Calorimeter.cc

Last change on this file was 1369, checked in by Michele Selvaggi, 11 years ago

added simplecalo
  • Property svn:keywords set to Id Revision Date
File size: 15.3 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: 2014-04-16 15:17:35 +0000 (Wed, 16 Apr 2014) $
8 * $Revision: 1369 $
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 // read resolution formulas
138 fECalResolutionFormula->Compile(GetString("ECalResolutionFormula", "0"));
139 fHCalResolutionFormula->Compile(GetString("HCalResolutionFormula", "0"));
140
141 // import array with output from other modules
142 fParticleInputArray = ImportArray(GetString("ParticleInputArray", "ParticlePropagator/particles"));
143 fItParticleInputArray = fParticleInputArray->MakeIterator();
144
145 fTrackInputArray = ImportArray(GetString("TrackInputArray", "ParticlePropagator/tracks"));
146 fItTrackInputArray = fTrackInputArray->MakeIterator();
147
148 // create output arrays
149 fTowerOutputArray = ExportArray(GetString("TowerOutputArray", "towers"));
150 fPhotonOutputArray = ExportArray(GetString("PhotonOutputArray", "photons"));
151
152 fEFlowTrackOutputArray = ExportArray(GetString("EFlowTrackOutputArray", "eflowTracks"));
153 fEFlowPhotonOutputArray = ExportArray(GetString("EFlowPhotonOutputArray", "eflowPhotons"));
154 fEFlowNeutralHadronOutputArray = ExportArray(GetString("EFlowNeutralHadronOutputArray", "eflowNeutralHadrons"));
155
156
157}
158
159//------------------------------------------------------------------------------
160
161void Calorimeter::Finish()
162{
163 vector< vector< Double_t >* >::iterator itPhiBin;
164 if(fItParticleInputArray) delete fItParticleInputArray;
165 if(fItTrackInputArray) delete fItTrackInputArray;
166 for(itPhiBin = fPhiBins.begin(); itPhiBin != fPhiBins.end(); ++itPhiBin)
167 {
168 delete *itPhiBin;
169 }
170}
171
172//------------------------------------------------------------------------------
173
174void Calorimeter::Process()
175{
176 Candidate *particle, *track;
177 TLorentzVector position, momentum;
178 Short_t etaBin, phiBin, flags;
179 Int_t number;
180 Long64_t towerHit, towerEtaPhi, hitEtaPhi;
181 Double_t ecalFraction, hcalFraction;
182 Double_t ecalEnergy, hcalEnergy;
183 Int_t pdgCode;
184
185 TFractionMap::iterator itFractionMap;
186
187 vector< Double_t >::iterator itEtaBin;
188 vector< Double_t >::iterator itPhiBin;
189 vector< Double_t > *phiBins;
190
191 vector< Long64_t >::iterator itTowerHits;
192
193 DelphesFactory *factory = GetFactory();
194 fTowerHits.clear();
195 fTowerECalFractions.clear();
196 fTowerHCalFractions.clear();
197 fTrackECalFractions.clear();
198 fTrackHCalFractions.clear();
199
200 // loop over all particles
201 fItParticleInputArray->Reset();
202 number = -1;
203 while((particle = static_cast<Candidate*>(fItParticleInputArray->Next())))
204 {
205 const TLorentzVector &particlePosition = particle->Position;
206 ++number;
207
208 pdgCode = TMath::Abs(particle->PID);
209
210 itFractionMap = fFractionMap.find(pdgCode);
211 if(itFractionMap == fFractionMap.end())
212 {
213 itFractionMap = fFractionMap.find(0);
214 }
215
216 ecalFraction = itFractionMap->second.first;
217 hcalFraction = itFractionMap->second.second;
218
219 fTowerECalFractions.push_back(ecalFraction);
220 fTowerHCalFractions.push_back(hcalFraction);
221
222 if(ecalFraction < 1.0E-9 && hcalFraction < 1.0E-9) continue;
223
224 // find eta bin [1, fEtaBins.size - 1]
225 itEtaBin = lower_bound(fEtaBins.begin(), fEtaBins.end(), particlePosition.Eta());
226 if(itEtaBin == fEtaBins.begin() || itEtaBin == fEtaBins.end()) continue;
227 etaBin = distance(fEtaBins.begin(), itEtaBin);
228
229 // phi bins for given eta bin
230 phiBins = fPhiBins[etaBin];
231
232 // find phi bin [1, phiBins.size - 1]
233 itPhiBin = lower_bound(phiBins->begin(), phiBins->end(), particlePosition.Phi());
234 if(itPhiBin == phiBins->begin() || itPhiBin == phiBins->end()) continue;
235 phiBin = distance(phiBins->begin(), itPhiBin);
236
237 flags = 0;
238 flags |= (pdgCode == 11 || pdgCode == 22) << 1;
239
240 // make tower hit {16-bits for eta bin number, 16-bits for phi bin number, 8-bits for flags, 24-bits for particle number}
241 towerHit = (Long64_t(etaBin) << 48) | (Long64_t(phiBin) << 32) | (Long64_t(flags) << 24) | Long64_t(number);
242
243 fTowerHits.push_back(towerHit);
244 }
245
246 // loop over all tracks
247 fItTrackInputArray->Reset();
248 number = -1;
249 while((track = static_cast<Candidate*>(fItTrackInputArray->Next())))
250 {
251 const TLorentzVector &trackPosition = track->Position;
252 ++number;
253
254 pdgCode = TMath::Abs(track->PID);
255
256 itFractionMap = fFractionMap.find(pdgCode);
257 if(itFractionMap == fFractionMap.end())
258 {
259 itFractionMap = fFractionMap.find(0);
260 }
261
262 ecalFraction = itFractionMap->second.first;
263 hcalFraction = itFractionMap->second.second;
264
265 fTrackECalFractions.push_back(ecalFraction);
266 fTrackHCalFractions.push_back(hcalFraction);
267
268 // find eta bin [1, fEtaBins.size - 1]
269 itEtaBin = lower_bound(fEtaBins.begin(), fEtaBins.end(), trackPosition.Eta());
270 if(itEtaBin == fEtaBins.begin() || itEtaBin == fEtaBins.end()) continue;
271 etaBin = distance(fEtaBins.begin(), itEtaBin);
272
273 // phi bins for given eta bin
274 phiBins = fPhiBins[etaBin];
275
276 // find phi bin [1, phiBins.size - 1]
277 itPhiBin = lower_bound(phiBins->begin(), phiBins->end(), trackPosition.Phi());
278 if(itPhiBin == phiBins->begin() || itPhiBin == phiBins->end()) continue;
279 phiBin = distance(phiBins->begin(), itPhiBin);
280
281 flags = 1;
282
283 // make tower hit {16-bits for eta bin number, 16-bits for phi bin number, 8-bits for flags, 24-bits for track number}
284 towerHit = (Long64_t(etaBin) << 48) | (Long64_t(phiBin) << 32) | (Long64_t(flags) << 24) | Long64_t(number);
285
286 fTowerHits.push_back(towerHit);
287 }
288
289 // all hits are sorted first by eta bin number, then by phi bin number,
290 // then by flags and then by particle or track number
291 sort(fTowerHits.begin(), fTowerHits.end());
292
293 // loop over all hits
294 towerEtaPhi = 0;
295 fTower = 0;
296 for(itTowerHits = fTowerHits.begin(); itTowerHits != fTowerHits.end(); ++itTowerHits)
297 {
298 towerHit = (*itTowerHits);
299 flags = (towerHit >> 24) & 0x00000000000000FFLL;
300 number = (towerHit) & 0x0000000000FFFFFFLL;
301 hitEtaPhi = towerHit >> 32;
302
303 if(towerEtaPhi != hitEtaPhi)
304 {
305 // switch to next tower
306 towerEtaPhi = hitEtaPhi;
307
308 // finalize previous tower
309 FinalizeTower();
310
311 // create new tower
312 fTower = factory->NewCandidate();
313
314 phiBin = (towerHit >> 32) & 0x000000000000FFFFLL;
315 etaBin = (towerHit >> 48) & 0x000000000000FFFFLL;
316
317 // phi bins for given eta bin
318 phiBins = fPhiBins[etaBin];
319
320 // calculate eta and phi of the tower's center
321 fTowerEta = 0.5*(fEtaBins[etaBin - 1] + fEtaBins[etaBin]);
322 fTowerPhi = 0.5*((*phiBins)[phiBin - 1] + (*phiBins)[phiBin]);
323
324 fTowerEdges[0] = fEtaBins[etaBin - 1];
325 fTowerEdges[1] = fEtaBins[etaBin];
326 fTowerEdges[2] = (*phiBins)[phiBin - 1];
327 fTowerEdges[3] = (*phiBins)[phiBin];
328
329 fTowerECalEnergy = 0.0;
330 fTowerHCalEnergy = 0.0;
331
332 fTrackECalEnergy = 0.0;
333 fTrackHCalEnergy = 0.0;
334
335 fTowerECalTime = 0.0;
336 fTowerHCalTime = 0.0;
337
338 fTrackECalTime = 0.0;
339 fTrackHCalTime = 0.0;
340
341 fTowerECalWeightTime = 0.0;
342 fTowerHCalWeightTime = 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
360 ecalEnergy = momentum.E() * fTrackECalFractions[number];
361 hcalEnergy = momentum.E() * fTrackHCalFractions[number];
362
363 fTrackECalEnergy += ecalEnergy;
364 fTrackHCalEnergy += hcalEnergy;
365
366 fTrackECalTime += TMath::Sqrt(ecalEnergy)*position.T();
367 fTrackHCalTime += TMath::Sqrt(hcalEnergy)*position.T();
368
369 fTrackECalWeightTime += TMath::Sqrt(ecalEnergy);
370 fTrackHCalWeightTime += TMath::Sqrt(hcalEnergy);
371
372 fTowerTrackArray->Add(track);
373
374 continue;
375 }
376
377 // check for photon and electron hits in current tower
378 if(flags & 2) ++fTowerPhotonHits;
379
380 particle = static_cast<Candidate*>(fParticleInputArray->At(number));
381 momentum = particle->Momentum;
382 position = particle->Position;
383
384 // fill current tower
385 ecalEnergy = momentum.E() * fTowerECalFractions[number];
386 hcalEnergy = momentum.E() * fTowerHCalFractions[number];
387
388 fTowerECalEnergy += ecalEnergy;
389 fTowerHCalEnergy += hcalEnergy;
390
391 fTowerECalTime += TMath::Sqrt(ecalEnergy)*position.T();
392 fTowerHCalTime += TMath::Sqrt(hcalEnergy)*position.T();
393
394 fTowerECalWeightTime += TMath::Sqrt(ecalEnergy);
395 fTowerHCalWeightTime += TMath::Sqrt(hcalEnergy);
396
397
398 fTower->AddCandidate(particle);
399 }
400
401 // finalize last tower
402 FinalizeTower();
403}
404
405//------------------------------------------------------------------------------
406
407void Calorimeter::FinalizeTower()
408{
409 Candidate *track, *tower;
410 Double_t energy, pt, eta, phi;
411 Double_t ecalEnergy, hcalEnergy;
412 Double_t ecalSigma, hcalSigma;
413 Double_t ecalTime, hcalTime, time;
414
415 if(!fTower) return;
416
417 ecalSigma = fECalResolutionFormula->Eval(0.0, fTowerEta, 0.0, fTowerECalEnergy);
418
419// ecalEnergy = gRandom->Gaus(fTowerECalEnergy, ecalSigma);
420// if(ecalEnergy < 0.0) ecalEnergy = 0.0;
421
422 ecalEnergy = LogNormal(fTowerECalEnergy, ecalSigma);
423 ecalTime = (fTowerECalWeightTime < 1.0E-09 ) ? 0 : fTowerECalTime/fTowerECalWeightTime;
424
425 hcalSigma = fHCalResolutionFormula->Eval(0.0, fTowerEta, 0.0, fTowerHCalEnergy);
426
427// hcalEnergy = gRandom->Gaus(fTowerHCalEnergy, hcalSigma);
428// if(hcalEnergy < 0.0) hcalEnergy = 0.0;
429
430 hcalEnergy = LogNormal(fTowerHCalEnergy, hcalSigma);
431 hcalTime = (fTowerHCalWeightTime < 1.0E-09 ) ? 0 : fTowerHCalTime/fTowerHCalWeightTime;
432
433 energy = ecalEnergy + hcalEnergy;
434 time = (TMath::Sqrt(ecalEnergy)*ecalTime + TMath::Sqrt(hcalEnergy)*hcalTime)/(TMath::Sqrt(ecalEnergy) + TMath::Sqrt(hcalEnergy));
435
436// eta = fTowerEta;
437// phi = fTowerPhi;
438
439 eta = gRandom->Uniform(fTowerEdges[0], fTowerEdges[1]);
440 phi = gRandom->Uniform(fTowerEdges[2], fTowerEdges[3]);
441
442 pt = energy / TMath::CosH(eta);
443
444 // fTower->Position.SetXYZT(-time, 0.0, 0.0, time);
445 fTower->Position.SetPtEtaPhiE(1.0, eta, phi, time);
446 fTower->Momentum.SetPtEtaPhiE(pt, eta, phi, energy);
447 fTower->Eem = ecalEnergy;
448 fTower->Ehad = hcalEnergy;
449
450 fTower->Edges[0] = fTowerEdges[0];
451 fTower->Edges[1] = fTowerEdges[1];
452 fTower->Edges[2] = fTowerEdges[2];
453 fTower->Edges[3] = fTowerEdges[3];
454
455
456 // fill calorimeter towers
457 if(energy > 0.0)
458 {
459 if(fTowerPhotonHits > 0 && fTowerTrackHits == 0)
460 {
461 fPhotonOutputArray->Add(fTower);
462 }
463
464 fTowerOutputArray->Add(fTower);
465 }
466
467 // fill energy flow candidates
468
469 // save all the tracks as energy flow tracks
470 fItTowerTrackArray->Reset();
471 while((track = static_cast<Candidate*>(fItTowerTrackArray->Next())))
472 {
473 fEFlowTrackOutputArray->Add(track);
474 }
475
476 ecalEnergy -= fTrackECalEnergy;
477 if(ecalEnergy < 0.0) ecalEnergy = 0.0;
478
479 hcalEnergy -= fTrackHCalEnergy;
480 if(hcalEnergy < 0.0) hcalEnergy = 0.0;
481
482 energy = ecalEnergy + hcalEnergy;
483
484
485 // save ECAL and/or HCAL energy excess as an energy flow tower
486 if(ecalEnergy > 0.0)
487 {
488 // create new photon tower
489 tower = static_cast<Candidate*>(fTower->Clone());
490
491 pt = ecalEnergy / TMath::CosH(eta);
492
493 tower->Momentum.SetPtEtaPhiE(pt, eta, phi, ecalEnergy);
494 tower->Eem = ecalEnergy;
495 tower->Ehad = 0;
496
497 fEFlowPhotonOutputArray->Add(tower);
498 }
499
500 if(hcalEnergy > 0.0)
501 {
502 // create new neutral hadron tower
503 tower = static_cast<Candidate*>(fTower->Clone());
504
505 pt = hcalEnergy / TMath::CosH(eta);
506
507 tower->Momentum.SetPtEtaPhiE(pt, eta, phi, hcalEnergy);
508 tower->Eem = 0;
509 tower->Ehad = hcalEnergy;
510
511 fEFlowNeutralHadronOutputArray->Add(tower);
512 }
513
514
515
516
517}
518
519//------------------------------------------------------------------------------
520
521Double_t Calorimeter::LogNormal(Double_t mean, Double_t sigma)
522{
523 Double_t a, b;
524
525 if(mean > 0.0)
526 {
527 b = TMath::Sqrt(TMath::Log((1.0 + (sigma*sigma)/(mean*mean))));
528 a = TMath::Log(mean) - 0.5*b*b;
529
530 return TMath::Exp(a + b*gRandom->Gaus(0, 1));
531 }
532 else
533 {
534 return 0.0;
535 }
536}
Note: See TracBrowser for help on using the repository browser.