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 |
|
---|
38 | using namespace std;
|
---|
39 |
|
---|
40 | //------------------------------------------------------------------------------
|
---|
41 |
|
---|
42 | SimpleCalorimeter::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 |
|
---|
55 | SimpleCalorimeter::~SimpleCalorimeter()
|
---|
56 | {
|
---|
57 | if(fResolutionFormula) delete fResolutionFormula;
|
---|
58 |
|
---|
59 | if(fTowerTrackArray) delete fTowerTrackArray;
|
---|
60 | if(fItTowerTrackArray) delete fItTowerTrackArray;
|
---|
61 | }
|
---|
62 |
|
---|
63 | //------------------------------------------------------------------------------
|
---|
64 |
|
---|
65 | void 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 |
|
---|
149 | void 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 |
|
---|
162 | void 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 |
|
---|
369 | void 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 |
|
---|
424 | Double_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 | }
|
---|