Fork me on GitHub

source: svn/trunk/modules/SimpleCalorimeter.cc

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

added simplecalo

File size: 12.2 KB
Line 
1
2/** \class SimpleCalorimeter
3 *
4 * Fills SimpleCalorimeter towers, performs SimpleCalorimeter resolution smearing,
5 * and creates energy flow objects (tracks, photons, and neutral hadrons).
6 *
7 * $Date: 2014-04-16 15:29:31 +0200 (Wed, 16 Apr 2014) $
8 * $Revision: 1364 $
9 *
10 *
11 * \author P. Demin - UCL, Louvain-la-Neuve
12 *
13 */
14
15#include "modules/SimpleCalorimeter.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
42SimpleCalorimeter::SimpleCalorimeter() :
43 fResolutionFormula(0),
44 fItParticleInputArray(0), fItTrackInputArray(0),
45 fTowerTrackArray(0), fItTowerTrackArray(0)
46{
47 fResolutionFormula = new DelphesFormula;
48
49 fTowerTrackArray = new TObjArray;
50 fItTowerTrackArray = fTowerTrackArray->MakeIterator();
51}
52
53//------------------------------------------------------------------------------
54
55SimpleCalorimeter::~SimpleCalorimeter()
56{
57 if(fResolutionFormula) delete fResolutionFormula;
58
59 if(fTowerTrackArray) delete fTowerTrackArray;
60 if(fItTowerTrackArray) delete fItTowerTrackArray;
61}
62
63//------------------------------------------------------------------------------
64
65void SimpleCalorimeter::Init()
66{
67 ExRootConfParam param, paramEtaBins, paramPhiBins, paramFractions;
68 Long_t i, j, k, size, sizeEtaBins, sizePhiBins;
69 Double_t fraction;
70 TBinMap::iterator itEtaBin;
71 set< Double_t >::iterator itPhiBin;
72 vector< Double_t > *phiBins;
73
74 // read eta and phi bins
75 param = GetParam("EtaPhiBins");
76 size = param.GetSize();
77 fBinMap.clear();
78 fEtaBins.clear();
79 fPhiBins.clear();
80 for(i = 0; i < size/2; ++i)
81 {
82 paramEtaBins = param[i*2];
83 sizeEtaBins = paramEtaBins.GetSize();
84 paramPhiBins = param[i*2 + 1];
85 sizePhiBins = paramPhiBins.GetSize();
86
87 for(j = 0; j < sizeEtaBins; ++j)
88 {
89 for(k = 0; k < sizePhiBins; ++k)
90 {
91 fBinMap[paramEtaBins[j].GetDouble()].insert(paramPhiBins[k].GetDouble());
92 }
93 }
94 }
95
96 // for better performance we transform map of sets to parallel vectors:
97 // vector< double > and vector< vector< double >* >
98 for(itEtaBin = fBinMap.begin(); itEtaBin != fBinMap.end(); ++itEtaBin)
99 {
100 fEtaBins.push_back(itEtaBin->first);
101 phiBins = new vector< double >(itEtaBin->second.size());
102 fPhiBins.push_back(phiBins);
103 phiBins->clear();
104 for(itPhiBin = itEtaBin->second.begin(); itPhiBin != itEtaBin->second.end(); ++itPhiBin)
105 {
106 phiBins->push_back(*itPhiBin);
107 }
108 }
109
110 // read energy fractions for different particles
111 param = GetParam("EnergyFraction");
112 size = param.GetSize();
113
114 // set default energy fractions values
115 fFractionMap.clear();
116 fFractionMap[0] = 1.0;
117
118 for(i = 0; i < size/2; ++i)
119 {
120 paramFractions = param[i*2 + 1];
121 fraction = paramFractions[0].GetDouble();
122 fFractionMap[param[i*2].GetInt()] = fraction;
123 }
124/*
125 TFractionMap::iterator itFractionMap;
126 for(itFractionMap = fFractionMap.begin(); itFractionMap != fFractionMap.end(); ++itFractionMap)
127 {
128 cout << itFractionMap->first << " " << itFractionMap->second.first << " " << itFractionMap->second.second << endl;
129 }
130*/
131 // read resolution formulas
132 fResolutionFormula->Compile(GetString("ResolutionFormula", "0"));
133
134 // import array with output from other modules
135 fParticleInputArray = ImportArray(GetString("ParticleInputArray", "ParticlePropagator/particles"));
136 fItParticleInputArray = fParticleInputArray->MakeIterator();
137
138 fTrackInputArray = ImportArray(GetString("TrackInputArray", "ParticlePropagator/tracks"));
139 fItTrackInputArray = fTrackInputArray->MakeIterator();
140
141 // create output arrays
142 fTowerOutputArray = ExportArray(GetString("TowerOutputArray", "towers"));
143 fEFlowTowerOutputArray = ExportArray(GetString("EFlowTowerOutputArray", "eflowTowers"));
144
145}
146
147//------------------------------------------------------------------------------
148
149void SimpleCalorimeter::Finish()
150{
151 vector< vector< Double_t >* >::iterator itPhiBin;
152 if(fItParticleInputArray) delete fItParticleInputArray;
153 if(fItTrackInputArray) delete fItTrackInputArray;
154 for(itPhiBin = fPhiBins.begin(); itPhiBin != fPhiBins.end(); ++itPhiBin)
155 {
156 delete *itPhiBin;
157 }
158}
159
160//------------------------------------------------------------------------------
161
162void SimpleCalorimeter::Process()
163{
164 Candidate *particle, *track;
165 TLorentzVector position, momentum;
166 Short_t etaBin, phiBin, flags;
167 Int_t number;
168 Long64_t towerHit, towerEtaPhi, hitEtaPhi;
169 Double_t fraction;
170 Double_t energy;
171 Int_t pdgCode;
172
173 TFractionMap::iterator itFractionMap;
174
175 vector< Double_t >::iterator itEtaBin;
176 vector< Double_t >::iterator itPhiBin;
177 vector< Double_t > *phiBins;
178
179 vector< Long64_t >::iterator itTowerHits;
180
181 DelphesFactory *factory = GetFactory();
182 fTowerHits.clear();
183 fTowerFractions.clear();
184 fTrackFractions.clear();
185
186 // loop over all particles
187 fItParticleInputArray->Reset();
188 number = -1;
189 while((particle = static_cast<Candidate*>(fItParticleInputArray->Next())))
190 {
191 const TLorentzVector &particlePosition = particle->Position;
192 ++number;
193
194 pdgCode = TMath::Abs(particle->PID);
195
196 itFractionMap = fFractionMap.find(pdgCode);
197 if(itFractionMap == fFractionMap.end())
198 {
199 itFractionMap = fFractionMap.find(0);
200 }
201
202 fraction = itFractionMap->second;
203 fTowerFractions.push_back(fraction);
204
205 if(fraction < 1.0E-9) continue;
206
207 // find eta bin [1, fEtaBins.size - 1]
208 itEtaBin = lower_bound(fEtaBins.begin(), fEtaBins.end(), particlePosition.Eta());
209 if(itEtaBin == fEtaBins.begin() || itEtaBin == fEtaBins.end()) continue;
210 etaBin = distance(fEtaBins.begin(), itEtaBin);
211
212 // phi bins for given eta bin
213 phiBins = fPhiBins[etaBin];
214
215 // find phi bin [1, phiBins.size - 1]
216 itPhiBin = lower_bound(phiBins->begin(), phiBins->end(), particlePosition.Phi());
217 if(itPhiBin == phiBins->begin() || itPhiBin == phiBins->end()) continue;
218 phiBin = distance(phiBins->begin(), itPhiBin);
219
220 flags = 0;
221 flags |= (pdgCode == 11 || pdgCode == 22) << 1;
222
223 // make tower hit {16-bits for eta bin number, 16-bits for phi bin number, 8-bits for flags, 24-bits for particle number}
224 towerHit = (Long64_t(etaBin) << 48) | (Long64_t(phiBin) << 32) | (Long64_t(flags) << 24) | Long64_t(number);
225
226 fTowerHits.push_back(towerHit);
227 }
228
229 // loop over all tracks
230 fItTrackInputArray->Reset();
231 number = -1;
232 while((track = static_cast<Candidate*>(fItTrackInputArray->Next())))
233 {
234 const TLorentzVector &trackPosition = track->Position;
235 ++number;
236
237 pdgCode = TMath::Abs(track->PID);
238
239 itFractionMap = fFractionMap.find(pdgCode);
240 if(itFractionMap == fFractionMap.end())
241 {
242 itFractionMap = fFractionMap.find(0);
243 }
244
245 fraction = itFractionMap->second;
246
247 fTrackFractions.push_back(fraction);
248
249 // find eta bin [1, fEtaBins.size - 1]
250 itEtaBin = lower_bound(fEtaBins.begin(), fEtaBins.end(), trackPosition.Eta());
251 if(itEtaBin == fEtaBins.begin() || itEtaBin == fEtaBins.end()) continue;
252 etaBin = distance(fEtaBins.begin(), itEtaBin);
253
254 // phi bins for given eta bin
255 phiBins = fPhiBins[etaBin];
256
257 // find phi bin [1, phiBins.size - 1]
258 itPhiBin = lower_bound(phiBins->begin(), phiBins->end(), trackPosition.Phi());
259 if(itPhiBin == phiBins->begin() || itPhiBin == phiBins->end()) continue;
260 phiBin = distance(phiBins->begin(), itPhiBin);
261
262 flags = 1;
263
264 // make tower hit {16-bits for eta bin number, 16-bits for phi bin number, 8-bits for flags, 24-bits for track number}
265 towerHit = (Long64_t(etaBin) << 48) | (Long64_t(phiBin) << 32) | (Long64_t(flags) << 24) | Long64_t(number);
266
267 fTowerHits.push_back(towerHit);
268 }
269
270 // all hits are sorted first by eta bin number, then by phi bin number,
271 // then by flags and then by particle or track number
272 sort(fTowerHits.begin(), fTowerHits.end());
273
274 // loop over all hits
275 towerEtaPhi = 0;
276 fTower = 0;
277 for(itTowerHits = fTowerHits.begin(); itTowerHits != fTowerHits.end(); ++itTowerHits)
278 {
279 towerHit = (*itTowerHits);
280 flags = (towerHit >> 24) & 0x00000000000000FFLL;
281 number = (towerHit) & 0x0000000000FFFFFFLL;
282 hitEtaPhi = towerHit >> 32;
283
284 if(towerEtaPhi != hitEtaPhi)
285 {
286 // switch to next tower
287 towerEtaPhi = hitEtaPhi;
288
289 // finalize previous tower
290 FinalizeTower();
291
292 // create new tower
293 fTower = factory->NewCandidate();
294
295 phiBin = (towerHit >> 32) & 0x000000000000FFFFLL;
296 etaBin = (towerHit >> 48) & 0x000000000000FFFFLL;
297
298 // phi bins for given eta bin
299 phiBins = fPhiBins[etaBin];
300
301 // calculate eta and phi of the tower's center
302 fTowerEta = 0.5*(fEtaBins[etaBin - 1] + fEtaBins[etaBin]);
303 fTowerPhi = 0.5*((*phiBins)[phiBin - 1] + (*phiBins)[phiBin]);
304
305 fTowerEdges[0] = fEtaBins[etaBin - 1];
306 fTowerEdges[1] = fEtaBins[etaBin];
307 fTowerEdges[2] = (*phiBins)[phiBin - 1];
308 fTowerEdges[3] = (*phiBins)[phiBin];
309
310 fTowerEnergy = 0.0;
311 fTrackEnergy = 0.0;
312
313 fTowerTime = 0.0;
314 fTrackTime = 0.0;
315
316 fTowerWeightTime = 0.0;
317
318 fTowerTrackHits = 0;
319 fTowerPhotonHits = 0;
320
321 fTowerTrackArray->Clear();
322 }
323
324 // check for track hits
325 if(flags & 1)
326 {
327 ++fTowerTrackHits;
328
329 track = static_cast<Candidate*>(fTrackInputArray->At(number));
330 momentum = track->Momentum;
331 position = track->Position;
332
333 energy = momentum.E() * fTrackFractions[number];
334
335 fTrackEnergy += energy;
336
337 fTrackTime += TMath::Sqrt(energy)*position.T();
338 fTrackWeightTime += TMath::Sqrt(energy);
339
340 fTowerTrackArray->Add(track);
341
342 continue;
343 }
344
345 // check for photon and electron hits in current tower
346 if(flags & 2) ++fTowerPhotonHits;
347
348 particle = static_cast<Candidate*>(fParticleInputArray->At(number));
349 momentum = particle->Momentum;
350 position = particle->Position;
351
352 // fill current tower
353 energy = momentum.E() * fTowerFractions[number];
354
355 fTowerEnergy += energy;
356
357 fTowerTime += TMath::Sqrt(energy)*position.T();
358 fTowerWeightTime += TMath::Sqrt(energy);
359
360 fTower->AddCandidate(particle);
361 }
362
363 // finalize last tower
364 FinalizeTower();
365}
366
367//------------------------------------------------------------------------------
368
369void SimpleCalorimeter::FinalizeTower()
370{
371 Candidate *tower;
372 Double_t energy, pt, eta, phi;
373 Double_t sigma;
374 Double_t time;
375
376 if(!fTower) return;
377
378 sigma = fResolutionFormula->Eval(0.0, fTowerEta, 0.0, fTowerEnergy);
379
380// energy = gRandom->Gaus(fTowerEnergy, sigma);
381// if(energy < 0.0) energy = 0.0;
382
383 energy = LogNormal(fTowerEnergy, sigma);
384 time = (fTowerWeightTime < 1.0E-09 ) ? 0 : fTowerTime/fTowerWeightTime;
385
386 eta = gRandom->Uniform(fTowerEdges[0], fTowerEdges[1]);
387 phi = gRandom->Uniform(fTowerEdges[2], fTowerEdges[3]);
388
389 pt = energy / TMath::CosH(eta);
390
391 // fTower->Position.SetXYZT(-time, 0.0, 0.0, time);
392 fTower->Position.SetPtEtaPhiE(1.0, eta, phi, time);
393 fTower->Momentum.SetPtEtaPhiE(pt, eta, phi, energy);
394
395 fTower->Edges[0] = fTowerEdges[0];
396 fTower->Edges[1] = fTowerEdges[1];
397 fTower->Edges[2] = fTowerEdges[2];
398 fTower->Edges[3] = fTowerEdges[3];
399
400
401 // fill SimpleCalorimeter towers
402 if(energy > 0.0) fTowerOutputArray->Add(fTower);
403
404
405 // fill energy flow candidates
406 energy -= fTrackEnergy;
407 if(energy < 0.0) energy = 0.0;
408
409 // save energy excess as an energy flow tower
410 if(energy > 0.0)
411 {
412 // create new photon tower
413 tower = static_cast<Candidate*>(fTower->Clone());
414 pt = energy / TMath::CosH(eta);
415
416 tower->Momentum.SetPtEtaPhiE(pt, eta, phi, energy);
417 fEFlowTowerOutputArray->Add(tower);
418 }
419
420}
421
422//------------------------------------------------------------------------------
423
424Double_t SimpleCalorimeter::LogNormal(Double_t mean, Double_t sigma)
425{
426 Double_t a, b;
427
428 if(mean > 0.0)
429 {
430 b = TMath::Sqrt(TMath::Log((1.0 + (sigma*sigma)/(mean*mean))));
431 a = TMath::Log(mean) - 0.5*b*b;
432
433 return TMath::Exp(a + b*gRandom->Gaus(0, 1));
434 }
435 else
436 {
437 return 0.0;
438 }
439}
Note: See TracBrowser for help on using the repository browser.