Fork me on GitHub

source: git/modules/Calorimeter.cc@ a817a22

ImprovedOutputFile Timing dual_readout llp
Last change on this file since a817a22 was d7d2da3, checked in by pavel <pavel@…>, 12 years ago

move branches/ModularDelphes to trunk

  • Property mode set to 100644
File size: 13.8 KB
Line 
1
2/** \class Calorimeter
3 *
4 * Fills calorimeter towers, performs calorimeter resolution smearing,
5 * preselects towers hit by photons and creates energy flow objects.
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 fTowerPhotonArray(0), fItTowerPhotonArray(0)
47{
48 fECalResolutionFormula = new DelphesFormula;
49 fHCalResolutionFormula = new DelphesFormula;
50 fTowerTrackArray = new TObjArray;
51 fItTowerTrackArray = fTowerTrackArray->MakeIterator();
52 fTowerPhotonArray = new TObjArray;
53 fItTowerPhotonArray = fTowerPhotonArray->MakeIterator();
54}
55
56//------------------------------------------------------------------------------
57
58Calorimeter::~Calorimeter()
59{
60 if(fECalResolutionFormula) delete fECalResolutionFormula;
61 if(fHCalResolutionFormula) delete fHCalResolutionFormula;
62 if(fTowerTrackArray) delete fTowerTrackArray;
63 if(fItTowerTrackArray) delete fItTowerTrackArray;
64 if(fTowerPhotonArray) delete fTowerPhotonArray;
65 if(fItTowerPhotonArray) delete fItTowerPhotonArray;
66}
67
68//------------------------------------------------------------------------------
69
70void Calorimeter::Init()
71{
72 ExRootConfParam param, paramEtaBins, paramPhiBins, paramFractions;
73 Long_t i, j, k, size, sizeEtaBins, sizePhiBins, sizeFractions;
74 Double_t ecalFraction, hcalFraction;
75 TBinMap::iterator itEtaBin;
76 set< Double_t >::iterator itPhiBin;
77 vector< Double_t > *phiBins;
78
79 // read eta and phi bins
80 param = GetParam("EtaPhiBins");
81 size = param.GetSize();
82 fBinMap.clear();
83 fEtaBins.clear();
84 fPhiBins.clear();
85 for(i = 0; i < size/2; ++i)
86 {
87 paramEtaBins = param[i*2];
88 sizeEtaBins = paramEtaBins.GetSize();
89 paramPhiBins = param[i*2 + 1];
90 sizePhiBins = paramPhiBins.GetSize();
91
92 for(j = 0; j < sizeEtaBins; ++j)
93 {
94 for(k = 0; k < sizePhiBins; ++k)
95 {
96 fBinMap[paramEtaBins[j].GetDouble()].insert(paramPhiBins[k].GetDouble());
97 }
98 }
99 }
100
101 // for better performance we transform map of sets to parallel vectors:
102 // vector< double > and vector< vector< double >* >
103 for(itEtaBin = fBinMap.begin(); itEtaBin != fBinMap.end(); ++itEtaBin)
104 {
105 fEtaBins.push_back(itEtaBin->first);
106 phiBins = new vector< double >(itEtaBin->second.size());
107 fPhiBins.push_back(phiBins);
108 phiBins->clear();
109 for(itPhiBin = itEtaBin->second.begin(); itPhiBin != itEtaBin->second.end(); ++itPhiBin)
110 {
111 phiBins->push_back(*itPhiBin);
112 }
113 }
114
115 // read energy fractions for different particles
116 param = GetParam("EnergyFraction");
117 size = param.GetSize();
118
119 // set default energy fractions values
120 fFractionMap.clear();
121 fFractionMap[0] = make_pair(0.0, 1.0);
122
123 for(i = 0; i < size/2; ++i)
124 {
125 paramFractions = param[i*2 + 1];
126 sizeFractions = paramFractions.GetSize();
127
128 ecalFraction = paramFractions[0].GetDouble();
129 hcalFraction = paramFractions[1].GetDouble();
130
131 fFractionMap[param[i*2].GetInt()] = make_pair(ecalFraction, hcalFraction);
132 }
133/*
134 TFractionMap::iterator itFractionMap;
135 for(itFractionMap = fFractionMap.begin(); itFractionMap != fFractionMap.end(); ++itFractionMap)
136 {
137 cout << itFractionMap->first << " " << itFractionMap->second.first << " " << itFractionMap->second.second << endl;
138 }
139*/
140 // read resolution formulas
141 fECalResolutionFormula->Compile(GetString("ECalResolutionFormula", "0"));
142 fHCalResolutionFormula->Compile(GetString("HCalResolutionFormula", "0"));
143
144 // import array with output from other modules
145 fParticleInputArray = ImportArray(GetString("ParticleInputArray", "ParticlePropagator/particles"));
146 fItParticleInputArray = fParticleInputArray->MakeIterator();
147
148 fTrackInputArray = ImportArray(GetString("TrackInputArray", "ParticlePropagator/tracks"));
149 fItTrackInputArray = fTrackInputArray->MakeIterator();
150
151 // create output arrays
152 fTowerOutputArray = ExportArray(GetString("TowerOutputArray", "towers"));
153 fPhotonOutputArray = ExportArray(GetString("PhotonOutputArray", "photons"));
154
155 fEFlowTrackOutputArray = ExportArray(GetString("EFlowTrackOutputArray", "eflowTracks"));
156 fEFlowTowerOutputArray = ExportArray(GetString("EFlowTowerOutputArray", "eflowTowers"));
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 fECalFractions.clear();
196 fHCalFractions.clear();
197
198 // loop over all particles
199 fItParticleInputArray->Reset();
200 number = -1;
201 while((particle = static_cast<Candidate*>(fItParticleInputArray->Next())))
202 {
203 const TLorentzVector &particlePosition = particle->Position;
204 ++number;
205
206 pdgCode = TMath::Abs(particle->PID);
207
208 itFractionMap = fFractionMap.find(pdgCode);
209 if(itFractionMap == fFractionMap.end())
210 {
211 itFractionMap = fFractionMap.find(0);
212 }
213
214 ecalFraction = itFractionMap->second.first;
215 hcalFraction = itFractionMap->second.second;
216
217 fECalFractions.push_back(ecalFraction);
218 fHCalFractions.push_back(hcalFraction);
219
220 if(ecalFraction < 1.0E-9 && hcalFraction < 1.0E-9) continue;
221
222 // find eta bin [1, fEtaBins.size - 1]
223 itEtaBin = lower_bound(fEtaBins.begin(), fEtaBins.end(), particlePosition.Eta());
224 if(itEtaBin == fEtaBins.begin() || itEtaBin == fEtaBins.end()) continue;
225 etaBin = distance(fEtaBins.begin(), itEtaBin);
226
227 // phi bins for given eta bin
228 phiBins = fPhiBins[etaBin];
229
230 // find phi bin [1, phiBins.size - 1]
231 itPhiBin = lower_bound(phiBins->begin(), phiBins->end(), particlePosition.Phi());
232 if(itPhiBin == phiBins->begin() || itPhiBin == phiBins->end()) continue;
233 phiBin = distance(phiBins->begin(), itPhiBin);
234
235 flags = (particle->Charge == 0);
236 flags |= (pdgCode == 22) << 1;
237 flags |= (pdgCode == 11) << 2;
238
239 // make tower hit {16-bits for eta bin number, 16-bits for phi bin number, 8-bits for flags, 24-bits for particle number}
240 towerHit = (Long64_t(etaBin) << 48) | (Long64_t(phiBin) << 32) | (Long64_t(flags) << 24) | Long64_t(number);
241
242 fTowerHits.push_back(towerHit);
243 }
244
245 // loop over all tracks
246 fItTrackInputArray->Reset();
247 number = -1;
248 while((track = static_cast<Candidate*>(fItTrackInputArray->Next())))
249 {
250 const TLorentzVector &trackPosition = track->Position;
251 ++number;
252
253 // find eta bin [1, fEtaBins.size - 1]
254 itEtaBin = lower_bound(fEtaBins.begin(), fEtaBins.end(), trackPosition.Eta());
255 if(itEtaBin == fEtaBins.begin() || itEtaBin == fEtaBins.end()) continue;
256 etaBin = distance(fEtaBins.begin(), itEtaBin);
257
258 // phi bins for given eta bin
259 phiBins = fPhiBins[etaBin];
260
261 // find phi bin [1, phiBins.size - 1]
262 itPhiBin = lower_bound(phiBins->begin(), phiBins->end(), trackPosition.Phi());
263 if(itPhiBin == phiBins->begin() || itPhiBin == phiBins->end()) continue;
264 phiBin = distance(phiBins->begin(), itPhiBin);
265
266 // make tower hit {16-bits for eta bin number, 16-bits for phi bin number, 8-bits for flags, 24-bits for track number}
267 towerHit = (Long64_t(etaBin) << 48) | (Long64_t(phiBin) << 32) | (Long64_t(1) << 27) | Long64_t(number);
268
269 fTowerHits.push_back(towerHit);
270 }
271
272 // all hits are sorted first by eta bin number, then by phi bin number,
273 // then by flags and then by particle or track number
274 sort(fTowerHits.begin(), fTowerHits.end());
275
276 // loop over all hits
277 towerEtaPhi = 0;
278 fTower = 0;
279 for(itTowerHits = fTowerHits.begin(); itTowerHits != fTowerHits.end(); ++itTowerHits)
280 {
281 towerHit = (*itTowerHits);
282 flags = (towerHit >> 24) & 0x00000000000000FFLL;
283 number = (towerHit) & 0x0000000000FFFFFFLL;
284 hitEtaPhi = towerHit >> 32;
285
286 if(towerEtaPhi != hitEtaPhi)
287 {
288 // switch to next tower
289 towerEtaPhi = hitEtaPhi;
290
291 // finalize previous tower
292 FinalizeTower();
293
294 // create new tower
295 fTower = factory->NewCandidate();
296
297 phiBin = (towerHit >> 32) & 0x000000000000FFFFLL;
298 etaBin = (towerHit >> 48) & 0x000000000000FFFFLL;
299
300 // phi bins for given eta bin
301 phiBins = fPhiBins[etaBin];
302
303 // calculate eta and phi of the tower's center
304 fTowerEta = 0.5*(fEtaBins[etaBin - 1] + fEtaBins[etaBin]);
305 fTowerPhi = 0.5*((*phiBins)[phiBin - 1] + (*phiBins)[phiBin]);
306
307 fTowerEdges[0] = fEtaBins[etaBin - 1];
308 fTowerEdges[1] = fEtaBins[etaBin];
309 fTowerEdges[2] = (*phiBins)[phiBin - 1];
310 fTowerEdges[3] = (*phiBins)[phiBin];
311
312 fTowerECalEnergy = 0.0;
313 fTowerHCalEnergy = 0.0;
314
315 fTowerECalNeutralEnergy = 0.0;
316 fTowerHCalNeutralEnergy = 0.0;
317
318 fTowerNeutralHits = 0;
319 fTowerPhotonHits = 0;
320 fTowerElectronHits = 0;
321 fTowerTrackHits = 0;
322 fTowerAllHits = 0;
323
324 fTowerTrackArray->Clear();
325 fTowerPhotonArray->Clear();
326 }
327
328 // check for track hits
329 if(flags & 8)
330 {
331 ++fTowerTrackHits;
332 track = static_cast<Candidate*>(fTrackInputArray->At(number));
333 fTowerTrackArray->Add(track);
334 continue;
335 }
336
337 particle = static_cast<Candidate*>(fParticleInputArray->At(number));
338 momentum = particle->Momentum;
339
340 // fill current tower
341 ecalEnergy = momentum.E() * fECalFractions[number];
342 hcalEnergy = momentum.E() * fHCalFractions[number];
343
344 fTowerECalEnergy += ecalEnergy;
345 fTowerHCalEnergy += hcalEnergy;
346
347 ++fTowerAllHits;
348 fTower->AddCandidate(particle);
349
350 // check for neutral hits in current tower
351 if(flags & 1) ++fTowerNeutralHits;
352
353 // check for photon hits in current tower
354 if(flags & 2)
355 {
356 ++fTowerPhotonHits;
357 fTowerPhotonArray->Add(particle);
358 }
359
360 // check for electron hits in current tower
361 if(flags & 4) ++fTowerElectronHits;
362 }
363
364 // finalize last tower
365 FinalizeTower();
366}
367
368//------------------------------------------------------------------------------
369
370void Calorimeter::FinalizeTower()
371{
372 Candidate *particle, *track, *tower;
373 Double_t energy, pt, eta, phi;
374 Double_t ecalEnergy, hcalEnergy;
375
376 if(!fTower) return;
377
378 ecalEnergy = gRandom->Gaus(fTowerECalEnergy, fECalResolutionFormula->Eval(0.0, fTowerEta, 0.0, fTowerECalEnergy));
379 if(ecalEnergy < 0.0) ecalEnergy = 0.0;
380
381 hcalEnergy = gRandom->Gaus(fTowerHCalEnergy, fHCalResolutionFormula->Eval(0.0, fTowerEta, 0.0, fTowerHCalEnergy));
382 if(hcalEnergy < 0.0) hcalEnergy = 0.0;
383
384 energy = ecalEnergy + hcalEnergy;
385
386 // eta = fTowerEta;
387 // phi = fTowerPhi;
388
389 eta = gRandom->Uniform(fTowerEdges[0], fTowerEdges[1]);
390 phi = gRandom->Uniform(fTowerEdges[2], fTowerEdges[3]);
391
392 pt = energy / TMath::CosH(eta);
393
394 fTower->Position.SetPtEtaPhiE(1.0, eta, phi, 0.0);
395 fTower->Momentum.SetPtEtaPhiE(pt, eta, phi, energy);
396 fTower->Eem = ecalEnergy;
397 fTower->Ehad = hcalEnergy;
398
399 fTower->Edges[0] = fTowerEdges[0];
400 fTower->Edges[1] = fTowerEdges[1];
401 fTower->Edges[2] = fTowerEdges[2];
402 fTower->Edges[3] = fTowerEdges[3];
403
404 // fill calorimeter towers and photon candidates
405 if(energy > 0.0)
406 {
407 if((fTowerPhotonHits > 0 || fTowerElectronHits > 0) &&
408 fTowerTrackHits == 0)
409 {
410 fPhotonOutputArray->Add(fTower);
411 }
412
413 fTowerOutputArray->Add(fTower);
414 }
415
416 // fill energy flow candidates
417 if(fTowerTrackHits == fTowerAllHits)
418 {
419 fItTowerTrackArray->Reset();
420 while((track = static_cast<Candidate*>(fItTowerTrackArray->Next())))
421 {
422 fEFlowTrackOutputArray->Add(track);
423 }
424 }
425 else if(fTowerTrackHits > 0 &&
426 fTowerElectronHits == 0 &&
427 fTowerPhotonHits + fTowerTrackHits == fTowerAllHits)
428 {
429 fItTowerTrackArray->Reset();
430 while((track = static_cast<Candidate*>(fItTowerTrackArray->Next())))
431 {
432 fEFlowTrackOutputArray->Add(track);
433 }
434
435 if(ecalEnergy > 0.0)
436 {
437 DelphesFactory *factory = GetFactory();
438
439 // create new tower
440 tower = factory->NewCandidate();
441
442 fItTowerPhotonArray->Reset();
443 while((particle = static_cast<Candidate*>(fItTowerPhotonArray->Next())))
444 {
445 tower->AddCandidate(particle);
446 }
447
448 pt = ecalEnergy / TMath::CosH(eta);
449
450 tower->Position.SetPtEtaPhiE(1.0, eta, phi, 0.0);
451 tower->Momentum.SetPtEtaPhiE(pt, eta, phi, ecalEnergy);
452 tower->Eem = ecalEnergy;
453 tower->Ehad = 0.0;
454
455 tower->Edges[0] = fTowerEdges[0];
456 tower->Edges[1] = fTowerEdges[1];
457 tower->Edges[2] = fTowerEdges[2];
458 tower->Edges[3] = fTowerEdges[3];
459
460 fEFlowTowerOutputArray->Add(tower);
461 }
462 }
463 else if(energy > 0.0)
464 {
465 fEFlowTowerOutputArray->Add(fTower);
466 }
467}
468
469//------------------------------------------------------------------------------
Note: See TracBrowser for help on using the repository browser.