Fork me on GitHub

source: git/modules/DualReadoutCalorimeter.cc@ a1b19ea

ImprovedOutputFile Timing dual_readout
Last change on this file since a1b19ea was a1b19ea, checked in by Michele Selvaggi <michele.selvaggi@…>, 6 years ago

first attempt at energy flow with Dual Readout

  • Property mode set to 100644
File size: 26.9 KB
Line 
1/*
2 * Delphes: a framework for fast simulation of a generic collider experiment
3 * Copyright (C) 2012-2014 Universite catholique de Louvain (UCL), Belgium
4 *
5 * This program is free software: you can redistribute it and/or modify
6 * it under the terms of the GNU General Public License as published by
7 * the Free Software Foundation, either version 3 of the License, or
8 * (at your option) any later version.
9 *
10 * This program is distributed in the hope that it will be useful,
11 * but WITHOUT ANY WARRANTY; without even the implied warranty of
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 * GNU General Public License for more details.
14 *
15 * You should have received a copy of the GNU General Public License
16 * along with this program. If not, see <http://www.gnu.org/licenses/>.
17 */
18
19
20/** \class DualReadoutCalorimeter
21 *
22
23 // ============ TODO =========================================================
24 // This implementation of dual calorimetry relies on several approximations:
25 // - If hadronic energy is found in the tower the energy resolution then the full tower enrgy is smeared according to hadronic resolution (pessimistic for (e,n) or (pi+,gamma))
26 // - While e+ vs pi+ (or gamma vs n) separation is in principle possible for single particles (using C/S, PMT timing, lateral shower profile) it is not obvious it can be done overlapping particles.
27 // Now we assume that regarless of the number of particle hits per tower we can always distinguish e+ vs pi+, which is probably not true in the case (e+,n) vs (pi+,gamma) without longitudinal segmentation.
28
29 *
30 * \author M. Selvaggi - CERN
31 *
32 */
33
34#include "modules/DualReadoutCalorimeter.h"
35
36#include "classes/DelphesClasses.h"
37#include "classes/DelphesFactory.h"
38#include "classes/DelphesFormula.h"
39
40#include "ExRootAnalysis/ExRootResult.h"
41#include "ExRootAnalysis/ExRootFilter.h"
42#include "ExRootAnalysis/ExRootClassifier.h"
43
44#include "TMath.h"
45#include "TString.h"
46#include "TFormula.h"
47#include "TRandom3.h"
48#include "TObjArray.h"
49#include "TDatabasePDG.h"
50#include "TLorentzVector.h"
51
52#include <algorithm>
53#include <stdexcept>
54#include <iostream>
55#include <sstream>
56
57using namespace std;
58
59//------------------------------------------------------------------------------
60
61DualReadoutCalorimeter::DualReadoutCalorimeter() :
62 fECalResolutionFormula(0), fHCalResolutionFormula(0),
63 fItParticleInputArray(0), fItTrackInputArray(0)
64{
65
66 fECalResolutionFormula = new DelphesFormula;
67 fHCalResolutionFormula = new DelphesFormula;
68
69 fECalTowerTrackArray = new TObjArray;
70 fItECalTowerTrackArray = fECalTowerTrackArray->MakeIterator();
71
72 fHCalTowerTrackArray = new TObjArray;
73 fItHCalTowerTrackArray = fHCalTowerTrackArray->MakeIterator();
74
75 fTowerTrackArray = new TObjArray;
76 fItTowerTrackArray = fTowerTrackArray->MakeIterator();
77
78}
79
80//------------------------------------------------------------------------------
81
82DualReadoutCalorimeter::~DualReadoutCalorimeter()
83{
84
85 if(fECalResolutionFormula) delete fECalResolutionFormula;
86 if(fHCalResolutionFormula) delete fHCalResolutionFormula;
87
88 if(fECalTowerTrackArray) delete fECalTowerTrackArray;
89 if(fItECalTowerTrackArray) delete fItECalTowerTrackArray;
90
91 if(fHCalTowerTrackArray) delete fHCalTowerTrackArray;
92 if(fItHCalTowerTrackArray) delete fItHCalTowerTrackArray;
93
94 if(fTowerTrackArray) delete fTowerTrackArray;
95 if(fItTowerTrackArray) delete fItTowerTrackArray;
96
97}
98
99//------------------------------------------------------------------------------
100
101void DualReadoutCalorimeter::Init()
102{
103 ExRootConfParam param, paramEtaBins, paramPhiBins, paramFractions;
104 Long_t i, j, k, size, sizeEtaBins, sizePhiBins;
105 Double_t ecalFraction, hcalFraction;
106 TBinMap::iterator itEtaBin;
107 set< Double_t >::iterator itPhiBin;
108 vector< Double_t > *phiBins;
109
110 // read eta and phi bins
111 param = GetParam("EtaPhiBins");
112 size = param.GetSize();
113 fBinMap.clear();
114 fEtaBins.clear();
115 fPhiBins.clear();
116 for(i = 0; i < size/2; ++i)
117 {
118 paramEtaBins = param[i*2];
119 sizeEtaBins = paramEtaBins.GetSize();
120 paramPhiBins = param[i*2 + 1];
121 sizePhiBins = paramPhiBins.GetSize();
122
123 for(j = 0; j < sizeEtaBins; ++j)
124 {
125 for(k = 0; k < sizePhiBins; ++k)
126 {
127 fBinMap[paramEtaBins[j].GetDouble()].insert(paramPhiBins[k].GetDouble());
128 }
129 }
130 }
131
132 // for better performance we transform map of sets to parallel vectors:
133 // vector< double > and vector< vector< double >* >
134 for(itEtaBin = fBinMap.begin(); itEtaBin != fBinMap.end(); ++itEtaBin)
135 {
136 fEtaBins.push_back(itEtaBin->first);
137 phiBins = new vector< double >(itEtaBin->second.size());
138 fPhiBins.push_back(phiBins);
139 phiBins->clear();
140 for(itPhiBin = itEtaBin->second.begin(); itPhiBin != itEtaBin->second.end(); ++itPhiBin)
141 {
142 phiBins->push_back(*itPhiBin);
143 }
144 }
145
146 // read energy fractions for different particles
147 param = GetParam("EnergyFraction");
148 size = param.GetSize();
149
150 // set default energy fractions values
151 fFractionMap.clear();
152 fFractionMap[0] = make_pair(0.0, 1.0);
153
154 for(i = 0; i < size/2; ++i)
155 {
156 paramFractions = param[i*2 + 1];
157
158 ecalFraction = paramFractions[0].GetDouble();
159 hcalFraction = paramFractions[1].GetDouble();
160
161 fFractionMap[param[i*2].GetInt()] = make_pair(ecalFraction, hcalFraction);
162 }
163
164 // read min E value for timing measurement in ECAL
165 fTimingEnergyMin = GetDouble("TimingEnergyMin",4.);
166 // For timing
167 // So far this flag needs to be false
168 // Curved extrapolation not supported
169 fElectronsFromTrack = false;
170
171 // read min E value for towers to be saved
172 fECalEnergyMin = GetDouble("ECalEnergyMin", 0.0);
173 fHCalEnergyMin = GetDouble("HCalEnergyMin", 0.0);
174
175 fECalEnergySignificanceMin = GetDouble("ECalEnergySignificanceMin", 0.0);
176 fHCalEnergySignificanceMin = GetDouble("HCalEnergySignificanceMin", 0.0);
177
178 // switch on or off the dithering of the center of DualReadoutCalorimeter towers
179 fSmearTowerCenter = GetBool("SmearTowerCenter", true);
180
181 // read resolution formulas
182 fECalResolutionFormula->Compile(GetString("ECalResolutionFormula", "0"));
183 fHCalResolutionFormula->Compile(GetString("HCalResolutionFormula", "0"));
184
185 // import array with output from other modules
186 fParticleInputArray = ImportArray(GetString("ParticleInputArray", "ParticlePropagator/particles"));
187 fItParticleInputArray = fParticleInputArray->MakeIterator();
188
189 fTrackInputArray = ImportArray(GetString("TrackInputArray", "ParticlePropagator/tracks"));
190 fItTrackInputArray = fTrackInputArray->MakeIterator();
191
192 // create output arrays
193 fTowerOutputArray = ExportArray(GetString("TowerOutputArray", "towers"));
194 fPhotonOutputArray = ExportArray(GetString("PhotonOutputArray", "photons"));
195
196 fEFlowTrackOutputArray = ExportArray(GetString("EFlowTrackOutputArray", "eflowTracks"));
197 fEFlowPhotonOutputArray = ExportArray(GetString("EFlowPhotonOutputArray", "eflowPhotons"));
198 fEFlowNeutralHadronOutputArray = ExportArray(GetString("EFlowNeutralHadronOutputArray", "eflowNeutralHadrons"));
199}
200
201//------------------------------------------------------------------------------
202
203void DualReadoutCalorimeter::Finish()
204{
205 vector< vector< Double_t >* >::iterator itPhiBin;
206 if(fItParticleInputArray) delete fItParticleInputArray;
207 if(fItTrackInputArray) delete fItTrackInputArray;
208 for(itPhiBin = fPhiBins.begin(); itPhiBin != fPhiBins.end(); ++itPhiBin)
209 {
210 delete *itPhiBin;
211 }
212}
213
214//------------------------------------------------------------------------------
215
216void DualReadoutCalorimeter::Process()
217{
218 Candidate *particle, *track;
219 TLorentzVector position, momentum;
220 Short_t etaBin, phiBin, flags;
221 Int_t number;
222 Long64_t towerHit, towerEtaPhi, hitEtaPhi;
223 Double_t ecalFraction, hcalFraction;
224 Double_t ecalEnergy, hcalEnergy;
225 Double_t ecalSigma, hcalSigma, sigma;
226 Double_t energyGuess, energy;
227 Int_t pdgCode;
228
229 TFractionMap::iterator itFractionMap;
230
231 vector< Double_t >::iterator itEtaBin;
232 vector< Double_t >::iterator itPhiBin;
233 vector< Double_t > *phiBins;
234
235 vector< Long64_t >::iterator itTowerHits;
236
237 DelphesFactory *factory = GetFactory();
238 fTowerHits.clear();
239 fECalTowerFractions.clear();
240 fHCalTowerFractions.clear();
241 fECalTrackFractions.clear();
242 fHCalTrackFractions.clear();
243
244 // loop over all particles
245 fItParticleInputArray->Reset();
246 number = -1;
247 while((particle = static_cast<Candidate*>(fItParticleInputArray->Next())))
248 {
249 const TLorentzVector &particlePosition = particle->Position;
250 ++number;
251
252 pdgCode = TMath::Abs(particle->PID);
253
254 itFractionMap = fFractionMap.find(pdgCode);
255 if(itFractionMap == fFractionMap.end())
256 {
257 itFractionMap = fFractionMap.find(0);
258 }
259
260 ecalFraction = itFractionMap->second.first;
261 hcalFraction = itFractionMap->second.second;
262
263 fECalTowerFractions.push_back(ecalFraction);
264 fHCalTowerFractions.push_back(hcalFraction);
265
266 if(ecalFraction < 1.0E-9 && hcalFraction < 1.0E-9) continue;
267
268 // find eta bin [1, fEtaBins.size - 1]
269 itEtaBin = lower_bound(fEtaBins.begin(), fEtaBins.end(), particlePosition.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(), particlePosition.Phi());
278 if(itPhiBin == phiBins->begin() || itPhiBin == phiBins->end()) continue;
279 phiBin = distance(phiBins->begin(), itPhiBin);
280
281 flags = 0;
282 flags |= (pdgCode == 11 || pdgCode == 22) << 1;
283
284 // make tower hit {16-bits for eta bin number, 16-bits for phi bin number, 8-bits for flags, 24-bits for particle number}
285 towerHit = (Long64_t(etaBin) << 48) | (Long64_t(phiBin) << 32) | (Long64_t(flags) << 24) | Long64_t(number);
286
287 fTowerHits.push_back(towerHit);
288 }
289
290 // loop over all tracks
291 fItTrackInputArray->Reset();
292 number = -1;
293 while((track = static_cast<Candidate*>(fItTrackInputArray->Next())))
294 {
295 const TLorentzVector &trackPosition = track->Position;
296 ++number;
297
298 pdgCode = TMath::Abs(track->PID);
299
300 itFractionMap = fFractionMap.find(pdgCode);
301 if(itFractionMap == fFractionMap.end())
302 {
303 itFractionMap = fFractionMap.find(0);
304 }
305
306 ecalFraction = itFractionMap->second.first;
307 hcalFraction = itFractionMap->second.second;
308
309 fECalTrackFractions.push_back(ecalFraction);
310 fHCalTrackFractions.push_back(hcalFraction);
311
312 // find eta bin [1, fEtaBins.size - 1]
313 itEtaBin = lower_bound(fEtaBins.begin(), fEtaBins.end(), trackPosition.Eta());
314 if(itEtaBin == fEtaBins.begin() || itEtaBin == fEtaBins.end()) continue;
315 etaBin = distance(fEtaBins.begin(), itEtaBin);
316
317 // phi bins for given eta bin
318 phiBins = fPhiBins[etaBin];
319
320 // find phi bin [1, phiBins.size - 1]
321 itPhiBin = lower_bound(phiBins->begin(), phiBins->end(), trackPosition.Phi());
322 if(itPhiBin == phiBins->begin() || itPhiBin == phiBins->end()) continue;
323 phiBin = distance(phiBins->begin(), itPhiBin);
324
325 flags = 1;
326
327 // make tower hit {16-bits for eta bin number, 16-bits for phi bin number, 8-bits for flags, 24-bits for track number}
328 towerHit = (Long64_t(etaBin) << 48) | (Long64_t(phiBin) << 32) | (Long64_t(flags) << 24) | Long64_t(number);
329
330 fTowerHits.push_back(towerHit);
331 }
332
333 // all hits are sorted first by eta bin number, then by phi bin number,
334 // then by flags and then by particle or track number
335 sort(fTowerHits.begin(), fTowerHits.end());
336
337 // loop over all hits
338 towerEtaPhi = 0;
339 fTower = 0;
340 for(itTowerHits = fTowerHits.begin(); itTowerHits != fTowerHits.end(); ++itTowerHits)
341 {
342 towerHit = (*itTowerHits);
343 flags = (towerHit >> 24) & 0x00000000000000FFLL;
344 number = (towerHit) & 0x0000000000FFFFFFLL;
345 hitEtaPhi = towerHit >> 32;
346
347 if(towerEtaPhi != hitEtaPhi)
348 {
349 // switch to next tower
350 towerEtaPhi = hitEtaPhi;
351
352 // finalize previous tower
353 FinalizeTower();
354
355 // create new tower
356 fTower = factory->NewCandidate();
357
358 phiBin = (towerHit >> 32) & 0x000000000000FFFFLL;
359 etaBin = (towerHit >> 48) & 0x000000000000FFFFLL;
360
361 // phi bins for given eta bin
362 phiBins = fPhiBins[etaBin];
363
364 // calculate eta and phi of the tower's center
365 fTowerEta = 0.5*(fEtaBins[etaBin - 1] + fEtaBins[etaBin]);
366 fTowerPhi = 0.5*((*phiBins)[phiBin - 1] + (*phiBins)[phiBin]);
367
368 fTowerEdges[0] = fEtaBins[etaBin - 1];
369 fTowerEdges[1] = fEtaBins[etaBin];
370 fTowerEdges[2] = (*phiBins)[phiBin - 1];
371 fTowerEdges[3] = (*phiBins)[phiBin];
372
373 fECalTowerEnergy = 0.0;
374 fHCalTowerEnergy = 0.0;
375
376 fECalTrackEnergy = 0.0;
377 fHCalTrackEnergy = 0.0;
378 fTrackEnergy = 0.0;
379
380 fECalTrackSigma = 0.0;
381 fHCalTrackSigma = 0.0;
382 fTrackSigma = 0.0;
383
384 fTowerTrackHits = 0;
385 fTowerPhotonHits = 0;
386
387 fECalTowerTrackArray->Clear();
388 fHCalTowerTrackArray->Clear();
389 fTowerTrackArray->Clear();
390
391 }
392
393 // check for track hits
394 if(flags & 1)
395 {
396 ++fTowerTrackHits;
397
398 track = static_cast<Candidate*>(fTrackInputArray->At(number));
399 momentum = track->Momentum;
400 position = track->Position;
401
402 ecalEnergy = momentum.E() * fECalTrackFractions[number];
403 hcalEnergy = momentum.E() * fHCalTrackFractions[number];
404 energy = ecalEnergy + hcalEnergy;
405
406 if(ecalEnergy > fTimingEnergyMin && fTower)
407 {
408 if(fElectronsFromTrack)
409 {
410 fTower->ECalEnergyTimePairs.push_back(make_pair<Float_t, Float_t>(ecalEnergy, track->Position.T()));
411 }
412 }
413
414
415 /*
416 if(fECalTrackFractions[number] > 1.0E-9 && fHCalTrackFractions[number] < 1.0E-9)
417 {
418 fECalTrackEnergy += ecalEnergy;
419 ecalSigma = fECalResolutionFormula->Eval(0.0, fTowerEta, 0.0, momentum.E());
420 if(ecalSigma/momentum.E() < track->TrackResolution) energyGuess = ecalEnergy;
421 else energyGuess = momentum.E();
422
423 fECalTrackSigma += (track->TrackResolution)*energyGuess*(track->TrackResolution)*energyGuess;
424 fECalTowerTrackArray->Add(track);
425 }
426
427 else if(fECalTrackFractions[number] < 1.0E-9 && fHCalTrackFractions[number] > 1.0E-9)
428 {
429 fHCalTrackEnergy += hcalEnergy;
430 hcalSigma = fHCalResolutionFormula->Eval(0.0, fTowerEta, 0.0, momentum.E());
431 if(hcalSigma/momentum.E() < track->TrackResolution) energyGuess = hcalEnergy;
432 else energyGuess = momentum.E();
433
434 fHCalTrackSigma += (track->TrackResolution)*energyGuess*(track->TrackResolution)*energyGuess;
435 fHCalTowerTrackArray->Add(track);
436 }
437
438 // muons
439 else if(fECalTrackFractions[number] < 1.0E-9 && fHCalTrackFractions[number] < 1.0E-9)
440 {
441 fEFlowTrackOutputArray->Add(track);
442 }
443 */
444
445 // in Dual Readout we do not care if tracks are ECAL of HCAL
446 if(fECalTrackFractions[number] > 1.0E-9 || fHCalTrackFractions[number] > 1.0E-9)
447 {
448 fTrackEnergy += energy;
449 // this sigma will be used to determine whether neutral excess is significant. We choose the resolution according to bthe higest deposited fraction (in practice had for charged hadrons and em for electrons)
450 sigma = 0.0;
451 if(fHCalTrackFractions[number] > 0)
452 sigma = fHCalResolutionFormula->Eval(0.0, fTowerEta, 0.0, momentum.E());
453 else
454 sigma = fECalResolutionFormula->Eval(0.0, fTowerEta, 0.0, momentum.E());
455
456 if(sigma/momentum.E() < track->TrackResolution)
457 energyGuess = ecalEnergy + hcalEnergy;
458 else
459 energyGuess = momentum.E();
460
461 fTrackSigma += (track->TrackResolution)*energyGuess*(track->TrackResolution)*energyGuess;
462 fTowerTrackArray->Add(track);
463
464 }
465 else
466 {
467 fEFlowTrackOutputArray->Add(track);
468 }
469
470 continue;
471 }
472
473 // check for photon and electron hits in current tower
474 if(flags & 2) ++fTowerPhotonHits;
475
476 particle = static_cast<Candidate*>(fParticleInputArray->At(number));
477 momentum = particle->Momentum;
478 position = particle->Position;
479
480 // fill current tower
481 ecalEnergy = momentum.E() * fECalTowerFractions[number];
482 hcalEnergy = momentum.E() * fHCalTowerFractions[number];
483
484 fECalTowerEnergy += ecalEnergy;
485 fHCalTowerEnergy += hcalEnergy;
486
487 if(ecalEnergy > fTimingEnergyMin && fTower)
488 {
489 if (abs(particle->PID) != 11 || !fElectronsFromTrack)
490 {
491 fTower->ECalEnergyTimePairs.push_back(make_pair<Float_t, Float_t>(ecalEnergy, particle->Position.T()));
492 }
493 }
494
495 fTower->AddCandidate(particle);
496 }
497
498 // finalize last tower
499 FinalizeTower();
500}
501
502//------------------------------------------------------------------------------
503
504void DualReadoutCalorimeter::FinalizeTower()
505{
506
507 Candidate *track, *tower, *mother;
508 Double_t energy, pt, eta, phi;
509 Double_t ecalEnergy, hcalEnergy;
510 Double_t ecalNeutralEnergy, hcalNeutralEnergy, neutralEnergy;
511
512 Double_t ecalSigma, hcalSigma, sigma;
513 Double_t ecalNeutralSigma, hcalNeutralSigma, neutralSigma;
514
515 Double_t weightTrack, weightCalo, bestEnergyEstimate, rescaleFactor;
516
517 TLorentzVector momentum;
518 TFractionMap::iterator itFractionMap;
519
520 Float_t weight, sumWeightedTime, sumWeight;
521
522 if(!fTower) return;
523
524
525 //if (fHCalTowerEnergy < 30 && fECalTowerEnergy < 30) return;
526 //cout<<"----------- New tower ---------"<<endl;
527
528
529 // here we change behaviour w.r.t to standard calorimeter. Start with worse case scenario. If fHCalTowerEnergy > 0, assume total energy smeared by HCAL resolution.
530 // For example, if overlapping charged pions and photons take hadronic resolution as combined measurement
531
532 // if no hadronic fraction at all, then use ECAL resolution
533
534 //cout<<"fECalTowerEnergy: "<<fECalTowerEnergy<<", fHCalTowerEnergy: "<<fHCalTowerEnergy<<", Eta: "<<fTowerEta<<endl;
535
536 // if no hadronic energy, use ECAL resolution
537 if (fHCalTowerEnergy <= fHCalEnergyMin)
538 {
539 energy = fECalTowerEnergy;
540 sigma = fECalResolutionFormula->Eval(0.0, fTowerEta, 0.0, energy);
541 //cout<<"em energy"<<energy<<", sigma: "<<sigma<<endl;
542 }
543
544 // if hadronic fraction > 0, use HCAL resolution
545 else
546 {
547 energy = fECalTowerEnergy + fHCalTowerEnergy;
548 sigma = fHCalResolutionFormula->Eval(0.0, fTowerEta, 0.0, energy);
549 //cout<<"had energy: "<<energy<<", sigma: "<<sigma<<endl;
550 }
551
552 energy = LogNormal(energy, sigma);
553 //cout<<energy<<","<<ecalEnergy<<","<<hcalEnergy<<endl;
554
555 if(energy < fEnergyMin || energy < fEnergySignificanceMin*sigma) energy = 0.0;
556
557 // for now keep this the same
558 ecalSigma = fECalResolutionFormula->Eval(0.0, fTowerEta, 0.0, fECalTowerEnergy);
559 hcalSigma = fHCalResolutionFormula->Eval(0.0, fTowerEta, 0.0, fHCalTowerEnergy);
560
561 ecalEnergy = LogNormal(fECalTowerEnergy, ecalSigma);
562 hcalEnergy = LogNormal(fHCalTowerEnergy, hcalSigma);
563
564 ecalSigma = fECalResolutionFormula->Eval(0.0, fTowerEta, 0.0, ecalEnergy);
565 hcalSigma = fHCalResolutionFormula->Eval(0.0, fTowerEta, 0.0, hcalEnergy);
566
567 if(ecalEnergy < fECalEnergyMin || ecalEnergy < fECalEnergySignificanceMin*ecalSigma) ecalEnergy = 0.0;
568 if(hcalEnergy < fHCalEnergyMin || hcalEnergy < fHCalEnergySignificanceMin*hcalSigma) hcalEnergy = 0.0;
569
570 //cout<<"Measured energy: "<<energy<<endl;
571
572 if(fSmearTowerCenter)
573 {
574 eta = gRandom->Uniform(fTowerEdges[0], fTowerEdges[1]);
575 phi = gRandom->Uniform(fTowerEdges[2], fTowerEdges[3]);
576 }
577 else
578 {
579 eta = fTowerEta;
580 phi = fTowerPhi;
581 }
582
583 pt = energy / TMath::CosH(eta);
584
585 // Time calculation for tower
586 fTower->NTimeHits = 0;
587 sumWeightedTime = 0.0;
588 sumWeight = 0.0;
589
590 for(size_t i = 0; i < fTower->ECalEnergyTimePairs.size(); ++i)
591 {
592 weight = TMath::Sqrt(fTower->ECalEnergyTimePairs[i].first);
593 sumWeightedTime += weight * fTower->ECalEnergyTimePairs[i].second;
594 sumWeight += weight;
595 fTower->NTimeHits++;
596 }
597
598 if(sumWeight > 0.0)
599 {
600 fTower->Position.SetPtEtaPhiE(1.0, eta, phi, sumWeightedTime/sumWeight);
601 }
602 else
603 {
604 fTower->Position.SetPtEtaPhiE(1.0, eta, phi, 999999.9);
605 }
606
607 fTower->Momentum.SetPtEtaPhiE(pt, eta, phi, energy);
608 fTower->Eem = ecalEnergy;
609 fTower->Ehad = hcalEnergy;
610
611 fTower->Edges[0] = fTowerEdges[0];
612 fTower->Edges[1] = fTowerEdges[1];
613 fTower->Edges[2] = fTowerEdges[2];
614 fTower->Edges[3] = fTowerEdges[3];
615
616 if(energy > 0.0)
617 {
618 if(fTowerPhotonHits > 0 && fTowerTrackHits == 0)
619 {
620 fPhotonOutputArray->Add(fTower);
621 }
622 fTowerOutputArray->Add(fTower);
623 }
624
625 // fill energy flow candidates
626
627 fTrackSigma = TMath::Sqrt(fTrackSigma);
628 neutralEnergy = max( (energy - fTrackEnergy) , 0.0);
629 neutralSigma = neutralEnergy / TMath::Sqrt(fTrackSigma*fTrackSigma + sigma*sigma);
630
631 //cout<<"trackEnergy: "<<fTrackEnergy<<", trackSigma: "<<fTrackSigma<<", Ntracks: "<<fTowerTrackArray->GetEntries()<<endl;
632
633 //cout<<"neutralEnergy: "<<neutralEnergy<<", neutralSigma: "<<neutralSigma<<endl;
634
635 // For now, if neutral excess is significant, simply create neutral EflowPhoton tower and clone each track into eflowtrack !!! -> Creating only photons !! EFlowNeutralHadron collection will be empy!!! TO BE FIXED
636 if(neutralEnergy > fEnergyMin && neutralSigma > fEnergySignificanceMin)
637 {
638
639 //cout<<"significant neutral excess found:"<<endl;
640 // create new photon tower
641 tower = static_cast<Candidate*>(fTower->Clone());
642 pt = neutralEnergy / TMath::CosH(eta);
643 //cout<<"Creating tower with Pt, Eta, Phi, Energy: "<<pt<<","<<eta<<","<<phi<<","<<neutralEnergy<<endl;
644 tower->Momentum.SetPtEtaPhiE(pt, eta, phi, neutralEnergy);
645 tower->Eem = neutralEnergy;
646 tower->Ehad = 0.0;
647 tower->PID = 22;
648
649 fEFlowPhotonOutputArray->Add(tower);
650
651
652 //clone tracks
653 fItTowerTrackArray->Reset();
654 while((track = static_cast<Candidate*>(fItTowerTrackArray->Next())))
655 {
656 //cout<<"looping over tracks"<<endl;
657 mother = track;
658 track = static_cast<Candidate*>(track->Clone());
659 track->AddCandidate(mother);
660 fEFlowTrackOutputArray->Add(track);
661 }
662 }
663
664
665 // if neutral excess is not significant, rescale eflow tracks, such that the total charged equals the best measurement given by the DualReadoutCalorimeter and tracking
666 else if(fTrackEnergy > 0.0)
667 {
668 //cout<<"no significant neutral excess found:"<<endl;
669 weightTrack = (fTrackSigma > 0.0) ? 1 / (fTrackSigma*fTrackSigma) : 0.0;
670 weightCalo = (sigma > 0.0) ? 1 / (sigma*sigma) : 0.0;
671
672 bestEnergyEstimate = (weightTrack*fTrackEnergy + weightCalo*energy) / (weightTrack + weightCalo);
673 rescaleFactor = bestEnergyEstimate/fTrackEnergy;
674
675 //rescale tracks
676 fItTowerTrackArray->Reset();
677 while((track = static_cast<Candidate*>(fItTowerTrackArray->Next())))
678 {
679 mother = track;
680 track = static_cast<Candidate*>(track->Clone());
681 track->AddCandidate(mother);
682
683 track->Momentum *= rescaleFactor;
684
685 fEFlowTrackOutputArray->Add(track);
686 }
687 }
688
689
690 /*
691 // fill energy flow candidates
692 fECalTrackSigma = TMath::Sqrt(fECalTrackSigma);
693 fHCalTrackSigma = TMath::Sqrt(fHCalTrackSigma);
694
695 //compute neutral excesses
696 ecalNeutralEnergy = max( (ecalEnergy - fECalTrackEnergy) , 0.0);
697 hcalNeutralEnergy = max( (hcalEnergy - fHCalTrackEnergy) , 0.0);
698
699 ecalNeutralSigma = ecalNeutralEnergy / TMath::Sqrt(fECalTrackSigma*fECalTrackSigma + ecalSigma*ecalSigma);
700 hcalNeutralSigma = hcalNeutralEnergy / TMath::Sqrt(fHCalTrackSigma*fHCalTrackSigma + hcalSigma*hcalSigma);
701
702 // if ecal neutral excess is significant, simply create neutral EflowPhoton tower and clone each track into eflowtrack
703 if(ecalNeutralEnergy > fECalEnergyMin && ecalNeutralSigma > fECalEnergySignificanceMin)
704 {
705 // create new photon tower
706 tower = static_cast<Candidate*>(fTower->Clone());
707 pt = ecalNeutralEnergy / TMath::CosH(eta);
708
709 tower->Momentum.SetPtEtaPhiE(pt, eta, phi, ecalNeutralEnergy);
710 tower->Eem = ecalNeutralEnergy;
711 tower->Ehad = 0.0;
712 tower->PID = 22;
713
714 fEFlowPhotonOutputArray->Add(tower);
715
716 //clone tracks
717 fItECalTowerTrackArray->Reset();
718 while((track = static_cast<Candidate*>(fItECalTowerTrackArray->Next())))
719 {
720 mother = track;
721 track = static_cast<Candidate*>(track->Clone());
722 track->AddCandidate(mother);
723
724 fEFlowTrackOutputArray->Add(track);
725 }
726
727 }
728
729 // if neutral excess is not significant, rescale eflow tracks, such that the total charged equals the best measurement given by the DualReadoutCalorimeter and tracking
730 else if(fECalTrackEnergy > 0.0)
731 {
732 weightTrack = (fECalTrackSigma > 0.0) ? 1 / (fECalTrackSigma*fECalTrackSigma) : 0.0;
733 weightCalo = (ecalSigma > 0.0) ? 1 / (ecalSigma*ecalSigma) : 0.0;
734
735 bestEnergyEstimate = (weightTrack*fECalTrackEnergy + weightCalo*ecalEnergy) / (weightTrack + weightCalo);
736 rescaleFactor = bestEnergyEstimate/fECalTrackEnergy;
737
738 //rescale tracks
739 fItECalTowerTrackArray->Reset();
740 while((track = static_cast<Candidate*>(fItECalTowerTrackArray->Next())))
741 {
742 mother = track;
743 track = static_cast<Candidate*>(track->Clone());
744 track->AddCandidate(mother);
745
746 track->Momentum *= rescaleFactor;
747
748 fEFlowTrackOutputArray->Add(track);
749 }
750 }
751
752
753 // if hcal neutral excess is significant, simply create neutral EflowNeutralHadron tower and clone each track into eflowtrack
754 if(hcalNeutralEnergy > fHCalEnergyMin && hcalNeutralSigma > fHCalEnergySignificanceMin)
755 {
756 // create new photon tower
757 tower = static_cast<Candidate*>(fTower->Clone());
758 pt = hcalNeutralEnergy / TMath::CosH(eta);
759
760 tower->Momentum.SetPtEtaPhiE(pt, eta, phi, hcalNeutralEnergy);
761 tower->Ehad = hcalNeutralEnergy;
762 tower->Eem = 0.0;
763
764 fEFlowNeutralHadronOutputArray->Add(tower);
765
766 //clone tracks
767 fItHCalTowerTrackArray->Reset();
768 while((track = static_cast<Candidate*>(fItHCalTowerTrackArray->Next())))
769 {
770 mother = track;
771 track = static_cast<Candidate*>(track->Clone());
772 track->AddCandidate(mother);
773
774 fEFlowTrackOutputArray->Add(track);
775 }
776
777 }
778
779 // if neutral excess is not significant, rescale eflow tracks, such that the total charged equals the best measurement given by the DualReadoutCalorimeter and tracking
780 else if(fHCalTrackEnergy > 0.0)
781 {
782 weightTrack = (fHCalTrackSigma > 0.0) ? 1 / (fHCalTrackSigma*fHCalTrackSigma) : 0.0;
783 weightCalo = (hcalSigma > 0.0) ? 1 / (hcalSigma*hcalSigma) : 0.0;
784
785 bestEnergyEstimate = (weightTrack*fHCalTrackEnergy + weightCalo*hcalEnergy) / (weightTrack + weightCalo);
786 rescaleFactor = bestEnergyEstimate / fHCalTrackEnergy;
787
788 //rescale tracks
789 fItHCalTowerTrackArray->Reset();
790 while((track = static_cast<Candidate*>(fItHCalTowerTrackArray->Next())))
791 {
792 mother = track;
793 track = static_cast<Candidate*>(track->Clone());
794 track->AddCandidate(mother);
795
796 track->Momentum *= rescaleFactor;
797
798 fEFlowTrackOutputArray->Add(track);
799 }
800 }
801
802 */
803}
804
805//------------------------------------------------------------------------------
806
807Double_t DualReadoutCalorimeter::LogNormal(Double_t mean, Double_t sigma)
808{
809 Double_t a, b;
810
811 if(mean > 0.0)
812 {
813 b = TMath::Sqrt(TMath::Log((1.0 + (sigma*sigma)/(mean*mean))));
814 a = TMath::Log(mean) - 0.5*b*b;
815
816 return TMath::Exp(a + b*gRandom->Gaus(0.0, 1.0));
817 }
818 else
819 {
820 return 0.0;
821 }
822}
Note: See TracBrowser for help on using the repository browser.