Fork me on GitHub

source: git/modules/DualReadoutCalorimeter.cc@ ef8a06d

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

first attempt to DR calo
  • Property mode set to 100644
File size: 22.4 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 * Fills DualReadoutCalorimeter towers, performs DualReadoutCalorimeter resolution smearing,
23 * and creates energy flow objects (tracks, photons, and neutral hadrons).
24 *
25 * \author P. Demin - UCL, Louvain-la-Neuve
26 *
27 */
28
29#include "modules/DualReadoutCalorimeter.h"
30
31#include "classes/DelphesClasses.h"
32#include "classes/DelphesFactory.h"
33#include "classes/DelphesFormula.h"
34
35#include "ExRootAnalysis/ExRootResult.h"
36#include "ExRootAnalysis/ExRootFilter.h"
37#include "ExRootAnalysis/ExRootClassifier.h"
38
39#include "TMath.h"
40#include "TString.h"
41#include "TFormula.h"
42#include "TRandom3.h"
43#include "TObjArray.h"
44#include "TDatabasePDG.h"
45#include "TLorentzVector.h"
46
47#include <algorithm>
48#include <stdexcept>
49#include <iostream>
50#include <sstream>
51
52using namespace std;
53
54//------------------------------------------------------------------------------
55
56DualReadoutCalorimeter::DualReadoutCalorimeter() :
57 fECalResolutionFormula(0), fHCalResolutionFormula(0),
58 fItParticleInputArray(0), fItTrackInputArray(0)
59{
60
61 fECalResolutionFormula = new DelphesFormula;
62 fHCalResolutionFormula = new DelphesFormula;
63
64 fECalTowerTrackArray = new TObjArray;
65 fItECalTowerTrackArray = fECalTowerTrackArray->MakeIterator();
66
67 fHCalTowerTrackArray = new TObjArray;
68 fItHCalTowerTrackArray = fHCalTowerTrackArray->MakeIterator();
69
70}
71
72//------------------------------------------------------------------------------
73
74DualReadoutCalorimeter::~DualReadoutCalorimeter()
75{
76
77 if(fECalResolutionFormula) delete fECalResolutionFormula;
78 if(fHCalResolutionFormula) delete fHCalResolutionFormula;
79
80 if(fECalTowerTrackArray) delete fECalTowerTrackArray;
81 if(fItECalTowerTrackArray) delete fItECalTowerTrackArray;
82
83 if(fHCalTowerTrackArray) delete fHCalTowerTrackArray;
84 if(fItHCalTowerTrackArray) delete fItHCalTowerTrackArray;
85
86}
87
88//------------------------------------------------------------------------------
89
90void DualReadoutCalorimeter::Init()
91{
92 ExRootConfParam param, paramEtaBins, paramPhiBins, paramFractions;
93 Long_t i, j, k, size, sizeEtaBins, sizePhiBins;
94 Double_t ecalFraction, hcalFraction;
95 TBinMap::iterator itEtaBin;
96 set< Double_t >::iterator itPhiBin;
97 vector< Double_t > *phiBins;
98
99 // read eta and phi bins
100 param = GetParam("EtaPhiBins");
101 size = param.GetSize();
102 fBinMap.clear();
103 fEtaBins.clear();
104 fPhiBins.clear();
105 for(i = 0; i < size/2; ++i)
106 {
107 paramEtaBins = param[i*2];
108 sizeEtaBins = paramEtaBins.GetSize();
109 paramPhiBins = param[i*2 + 1];
110 sizePhiBins = paramPhiBins.GetSize();
111
112 for(j = 0; j < sizeEtaBins; ++j)
113 {
114 for(k = 0; k < sizePhiBins; ++k)
115 {
116 fBinMap[paramEtaBins[j].GetDouble()].insert(paramPhiBins[k].GetDouble());
117 }
118 }
119 }
120
121 // for better performance we transform map of sets to parallel vectors:
122 // vector< double > and vector< vector< double >* >
123 for(itEtaBin = fBinMap.begin(); itEtaBin != fBinMap.end(); ++itEtaBin)
124 {
125 fEtaBins.push_back(itEtaBin->first);
126 phiBins = new vector< double >(itEtaBin->second.size());
127 fPhiBins.push_back(phiBins);
128 phiBins->clear();
129 for(itPhiBin = itEtaBin->second.begin(); itPhiBin != itEtaBin->second.end(); ++itPhiBin)
130 {
131 phiBins->push_back(*itPhiBin);
132 }
133 }
134
135 // read energy fractions for different particles
136 param = GetParam("EnergyFraction");
137 size = param.GetSize();
138
139 // set default energy fractions values
140 fFractionMap.clear();
141 fFractionMap[0] = make_pair(0.0, 1.0);
142
143 for(i = 0; i < size/2; ++i)
144 {
145 paramFractions = param[i*2 + 1];
146
147 ecalFraction = paramFractions[0].GetDouble();
148 hcalFraction = paramFractions[1].GetDouble();
149
150 fFractionMap[param[i*2].GetInt()] = make_pair(ecalFraction, hcalFraction);
151 }
152
153 // read min E value for timing measurement in ECAL
154 fTimingEnergyMin = GetDouble("TimingEnergyMin",4.);
155 // For timing
156 // So far this flag needs to be false
157 // Curved extrapolation not supported
158 fElectronsFromTrack = false;
159
160 // read min E value for towers to be saved
161 fECalEnergyMin = GetDouble("ECalEnergyMin", 0.0);
162 fHCalEnergyMin = GetDouble("HCalEnergyMin", 0.0);
163
164 fECalEnergySignificanceMin = GetDouble("ECalEnergySignificanceMin", 0.0);
165 fHCalEnergySignificanceMin = GetDouble("HCalEnergySignificanceMin", 0.0);
166
167 // switch on or off the dithering of the center of DualReadoutCalorimeter towers
168 fSmearTowerCenter = GetBool("SmearTowerCenter", true);
169
170 // read resolution formulas
171 fECalResolutionFormula->Compile(GetString("ECalResolutionFormula", "0"));
172 fHCalResolutionFormula->Compile(GetString("HCalResolutionFormula", "0"));
173
174 // import array with output from other modules
175 fParticleInputArray = ImportArray(GetString("ParticleInputArray", "ParticlePropagator/particles"));
176 fItParticleInputArray = fParticleInputArray->MakeIterator();
177
178 fTrackInputArray = ImportArray(GetString("TrackInputArray", "ParticlePropagator/tracks"));
179 fItTrackInputArray = fTrackInputArray->MakeIterator();
180
181 // create output arrays
182 fTowerOutputArray = ExportArray(GetString("TowerOutputArray", "towers"));
183 fPhotonOutputArray = ExportArray(GetString("PhotonOutputArray", "photons"));
184
185 fEFlowTrackOutputArray = ExportArray(GetString("EFlowTrackOutputArray", "eflowTracks"));
186 fEFlowPhotonOutputArray = ExportArray(GetString("EFlowPhotonOutputArray", "eflowPhotons"));
187 fEFlowNeutralHadronOutputArray = ExportArray(GetString("EFlowNeutralHadronOutputArray", "eflowNeutralHadrons"));
188}
189
190//------------------------------------------------------------------------------
191
192void DualReadoutCalorimeter::Finish()
193{
194 vector< vector< Double_t >* >::iterator itPhiBin;
195 if(fItParticleInputArray) delete fItParticleInputArray;
196 if(fItTrackInputArray) delete fItTrackInputArray;
197 for(itPhiBin = fPhiBins.begin(); itPhiBin != fPhiBins.end(); ++itPhiBin)
198 {
199 delete *itPhiBin;
200 }
201}
202
203//------------------------------------------------------------------------------
204
205void DualReadoutCalorimeter::Process()
206{
207 Candidate *particle, *track;
208 TLorentzVector position, momentum;
209 Short_t etaBin, phiBin, flags;
210 Int_t number;
211 Long64_t towerHit, towerEtaPhi, hitEtaPhi;
212 Double_t ecalFraction, hcalFraction;
213 Double_t ecalEnergy, hcalEnergy;
214 Double_t ecalSigma, hcalSigma;
215 Double_t energyGuess;
216 Int_t pdgCode;
217
218 TFractionMap::iterator itFractionMap;
219
220 vector< Double_t >::iterator itEtaBin;
221 vector< Double_t >::iterator itPhiBin;
222 vector< Double_t > *phiBins;
223
224 vector< Long64_t >::iterator itTowerHits;
225
226 DelphesFactory *factory = GetFactory();
227 fTowerHits.clear();
228 fECalTowerFractions.clear();
229 fHCalTowerFractions.clear();
230 fECalTrackFractions.clear();
231 fHCalTrackFractions.clear();
232
233 // loop over all particles
234 fItParticleInputArray->Reset();
235 number = -1;
236 while((particle = static_cast<Candidate*>(fItParticleInputArray->Next())))
237 {
238 const TLorentzVector &particlePosition = particle->Position;
239 ++number;
240
241 pdgCode = TMath::Abs(particle->PID);
242
243 itFractionMap = fFractionMap.find(pdgCode);
244 if(itFractionMap == fFractionMap.end())
245 {
246 itFractionMap = fFractionMap.find(0);
247 }
248
249 ecalFraction = itFractionMap->second.first;
250 hcalFraction = itFractionMap->second.second;
251
252 fECalTowerFractions.push_back(ecalFraction);
253 fHCalTowerFractions.push_back(hcalFraction);
254
255 if(ecalFraction < 1.0E-9 && hcalFraction < 1.0E-9) continue;
256
257 // find eta bin [1, fEtaBins.size - 1]
258 itEtaBin = lower_bound(fEtaBins.begin(), fEtaBins.end(), particlePosition.Eta());
259 if(itEtaBin == fEtaBins.begin() || itEtaBin == fEtaBins.end()) continue;
260 etaBin = distance(fEtaBins.begin(), itEtaBin);
261
262 // phi bins for given eta bin
263 phiBins = fPhiBins[etaBin];
264
265 // find phi bin [1, phiBins.size - 1]
266 itPhiBin = lower_bound(phiBins->begin(), phiBins->end(), particlePosition.Phi());
267 if(itPhiBin == phiBins->begin() || itPhiBin == phiBins->end()) continue;
268 phiBin = distance(phiBins->begin(), itPhiBin);
269
270 flags = 0;
271 flags |= (pdgCode == 11 || pdgCode == 22) << 1;
272
273 // make tower hit {16-bits for eta bin number, 16-bits for phi bin number, 8-bits for flags, 24-bits for particle number}
274 towerHit = (Long64_t(etaBin) << 48) | (Long64_t(phiBin) << 32) | (Long64_t(flags) << 24) | Long64_t(number);
275
276 fTowerHits.push_back(towerHit);
277 }
278
279 // loop over all tracks
280 fItTrackInputArray->Reset();
281 number = -1;
282 while((track = static_cast<Candidate*>(fItTrackInputArray->Next())))
283 {
284 const TLorentzVector &trackPosition = track->Position;
285 ++number;
286
287 pdgCode = TMath::Abs(track->PID);
288
289 itFractionMap = fFractionMap.find(pdgCode);
290 if(itFractionMap == fFractionMap.end())
291 {
292 itFractionMap = fFractionMap.find(0);
293 }
294
295 ecalFraction = itFractionMap->second.first;
296 hcalFraction = itFractionMap->second.second;
297
298 fECalTrackFractions.push_back(ecalFraction);
299 fHCalTrackFractions.push_back(hcalFraction);
300
301 // find eta bin [1, fEtaBins.size - 1]
302 itEtaBin = lower_bound(fEtaBins.begin(), fEtaBins.end(), trackPosition.Eta());
303 if(itEtaBin == fEtaBins.begin() || itEtaBin == fEtaBins.end()) continue;
304 etaBin = distance(fEtaBins.begin(), itEtaBin);
305
306 // phi bins for given eta bin
307 phiBins = fPhiBins[etaBin];
308
309 // find phi bin [1, phiBins.size - 1]
310 itPhiBin = lower_bound(phiBins->begin(), phiBins->end(), trackPosition.Phi());
311 if(itPhiBin == phiBins->begin() || itPhiBin == phiBins->end()) continue;
312 phiBin = distance(phiBins->begin(), itPhiBin);
313
314 flags = 1;
315
316 // make tower hit {16-bits for eta bin number, 16-bits for phi bin number, 8-bits for flags, 24-bits for track number}
317 towerHit = (Long64_t(etaBin) << 48) | (Long64_t(phiBin) << 32) | (Long64_t(flags) << 24) | Long64_t(number);
318
319 fTowerHits.push_back(towerHit);
320 }
321
322 // all hits are sorted first by eta bin number, then by phi bin number,
323 // then by flags and then by particle or track number
324 sort(fTowerHits.begin(), fTowerHits.end());
325
326 // loop over all hits
327 towerEtaPhi = 0;
328 fTower = 0;
329 for(itTowerHits = fTowerHits.begin(); itTowerHits != fTowerHits.end(); ++itTowerHits)
330 {
331 towerHit = (*itTowerHits);
332 flags = (towerHit >> 24) & 0x00000000000000FFLL;
333 number = (towerHit) & 0x0000000000FFFFFFLL;
334 hitEtaPhi = towerHit >> 32;
335
336 if(towerEtaPhi != hitEtaPhi)
337 {
338 // switch to next tower
339 towerEtaPhi = hitEtaPhi;
340
341 // finalize previous tower
342 FinalizeTower();
343
344 // create new tower
345 fTower = factory->NewCandidate();
346
347 phiBin = (towerHit >> 32) & 0x000000000000FFFFLL;
348 etaBin = (towerHit >> 48) & 0x000000000000FFFFLL;
349
350 // phi bins for given eta bin
351 phiBins = fPhiBins[etaBin];
352
353 // calculate eta and phi of the tower's center
354 fTowerEta = 0.5*(fEtaBins[etaBin - 1] + fEtaBins[etaBin]);
355 fTowerPhi = 0.5*((*phiBins)[phiBin - 1] + (*phiBins)[phiBin]);
356
357 fTowerEdges[0] = fEtaBins[etaBin - 1];
358 fTowerEdges[1] = fEtaBins[etaBin];
359 fTowerEdges[2] = (*phiBins)[phiBin - 1];
360 fTowerEdges[3] = (*phiBins)[phiBin];
361
362 fECalTowerEnergy = 0.0;
363 fHCalTowerEnergy = 0.0;
364
365 fECalTrackEnergy = 0.0;
366 fHCalTrackEnergy = 0.0;
367
368 fECalTrackSigma = 0.0;
369 fHCalTrackSigma = 0.0;
370
371 fTowerTrackHits = 0;
372 fTowerPhotonHits = 0;
373
374 fECalTowerTrackArray->Clear();
375 fHCalTowerTrackArray->Clear();
376
377 }
378
379 // check for track hits
380 if(flags & 1)
381 {
382 ++fTowerTrackHits;
383
384 track = static_cast<Candidate*>(fTrackInputArray->At(number));
385 momentum = track->Momentum;
386 position = track->Position;
387
388 ecalEnergy = momentum.E() * fECalTrackFractions[number];
389 hcalEnergy = momentum.E() * fHCalTrackFractions[number];
390
391 if(ecalEnergy > fTimingEnergyMin && fTower)
392 {
393 if(fElectronsFromTrack)
394 {
395 fTower->ECalEnergyTimePairs.push_back(make_pair<Float_t, Float_t>(ecalEnergy, track->Position.T()));
396 }
397 }
398
399 if(fECalTrackFractions[number] > 1.0E-9 && fHCalTrackFractions[number] < 1.0E-9)
400 {
401 fECalTrackEnergy += ecalEnergy;
402 ecalSigma = fECalResolutionFormula->Eval(0.0, fTowerEta, 0.0, momentum.E());
403 if(ecalSigma/momentum.E() < track->TrackResolution) energyGuess = ecalEnergy;
404 else energyGuess = momentum.E();
405
406 fECalTrackSigma += (track->TrackResolution)*energyGuess*(track->TrackResolution)*energyGuess;
407 fECalTowerTrackArray->Add(track);
408 }
409
410 else if(fECalTrackFractions[number] < 1.0E-9 && fHCalTrackFractions[number] > 1.0E-9)
411 {
412 fHCalTrackEnergy += hcalEnergy;
413 hcalSigma = fHCalResolutionFormula->Eval(0.0, fTowerEta, 0.0, momentum.E());
414 if(hcalSigma/momentum.E() < track->TrackResolution) energyGuess = hcalEnergy;
415 else energyGuess = momentum.E();
416
417 fHCalTrackSigma += (track->TrackResolution)*energyGuess*(track->TrackResolution)*energyGuess;
418 fHCalTowerTrackArray->Add(track);
419 }
420
421 else if(fECalTrackFractions[number] < 1.0E-9 && fHCalTrackFractions[number] < 1.0E-9)
422 {
423 fEFlowTrackOutputArray->Add(track);
424 }
425
426 continue;
427 }
428
429 // check for photon and electron hits in current tower
430 if(flags & 2) ++fTowerPhotonHits;
431
432 particle = static_cast<Candidate*>(fParticleInputArray->At(number));
433 momentum = particle->Momentum;
434 position = particle->Position;
435
436 // fill current tower
437 ecalEnergy = momentum.E() * fECalTowerFractions[number];
438 hcalEnergy = momentum.E() * fHCalTowerFractions[number];
439
440 fECalTowerEnergy += ecalEnergy;
441 fHCalTowerEnergy += hcalEnergy;
442
443 if(ecalEnergy > fTimingEnergyMin && fTower)
444 {
445 if (abs(particle->PID) != 11 || !fElectronsFromTrack)
446 {
447 fTower->ECalEnergyTimePairs.push_back(make_pair<Float_t, Float_t>(ecalEnergy, particle->Position.T()));
448 }
449 }
450
451 fTower->AddCandidate(particle);
452 }
453
454 // finalize last tower
455 FinalizeTower();
456}
457
458//------------------------------------------------------------------------------
459
460void DualReadoutCalorimeter::FinalizeTower()
461{
462
463
464
465
466 Candidate *track, *tower, *mother;
467 Double_t energy, pt, eta, phi;
468 Double_t ecalEnergy, hcalEnergy;
469 Double_t ecalNeutralEnergy, hcalNeutralEnergy;
470
471 Double_t ecalSigma, hcalSigma, sigma;
472 Double_t ecalNeutralSigma, hcalNeutralSigma;
473
474 Double_t weightTrack, weightCalo, bestEnergyEstimate, rescaleFactor;
475
476 TLorentzVector momentum;
477 TFractionMap::iterator itFractionMap;
478
479 Float_t weight, sumWeightedTime, sumWeight;
480
481 if(!fTower) return;
482
483
484 //if (fHCalTowerEnergy < 30 && fECalTowerEnergy < 30) return;
485 //cout<<"----------- New tower ---------"<<endl;
486
487
488 // 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.
489 // For example, if overlapping charged pions and photons take hadronic resolution as combined measurement
490
491 // if no hadronic fraction at all, then use ECAL resolution
492
493 //cout<<"fECalTowerEnergy: "<<fECalTowerEnergy<<", fHCalTowerEnergy: "<<fHCalTowerEnergy<<endl;
494
495 if (fHCalTowerEnergy <= fHCalEnergyMin)
496 {
497 energy = fECalTowerEnergy;
498 sigma = fECalResolutionFormula->Eval(0.0, fTowerEta, 0.0, energy);
499 //cout<<"em: "<<energy<<","<<sigma<<endl;
500 }
501
502 // if hadronic fraction > 0, use HCAL resolution
503 else
504 {
505 energy = fECalTowerEnergy + fHCalTowerEnergy;
506 sigma = fHCalResolutionFormula->Eval(0.0, fTowerEta, 0.0, energy);
507 //cout<<"had: "<<energy<<","<<sigma<<endl;
508 }
509
510 energy = LogNormal(energy, sigma);
511 //cout<<energy<<","<<ecalEnergy<<","<<hcalEnergy<<endl;
512
513 if(energy < fEnergyMin || energy < fEnergySignificanceMin*sigma) energy = 0.0;
514
515 // for now keep this the same
516 ecalSigma = fECalResolutionFormula->Eval(0.0, fTowerEta, 0.0, fECalTowerEnergy);
517 hcalSigma = fHCalResolutionFormula->Eval(0.0, fTowerEta, 0.0, fHCalTowerEnergy);
518
519 ecalEnergy = LogNormal(fECalTowerEnergy, ecalSigma);
520 hcalEnergy = LogNormal(fHCalTowerEnergy, hcalSigma);
521
522 ecalSigma = fECalResolutionFormula->Eval(0.0, fTowerEta, 0.0, ecalEnergy);
523 hcalSigma = fHCalResolutionFormula->Eval(0.0, fTowerEta, 0.0, hcalEnergy);
524
525 if(ecalEnergy < fECalEnergyMin || ecalEnergy < fECalEnergySignificanceMin*ecalSigma) ecalEnergy = 0.0;
526 if(hcalEnergy < fHCalEnergyMin || hcalEnergy < fHCalEnergySignificanceMin*hcalSigma) hcalEnergy = 0.0;
527 //cout<<energy<<","<<ecalEnergy<<","<<hcalEnergy<<endl;
528
529
530 if(fSmearTowerCenter)
531 {
532 eta = gRandom->Uniform(fTowerEdges[0], fTowerEdges[1]);
533 phi = gRandom->Uniform(fTowerEdges[2], fTowerEdges[3]);
534 }
535 else
536 {
537 eta = fTowerEta;
538 phi = fTowerPhi;
539 }
540
541 pt = energy / TMath::CosH(eta);
542
543 // Time calculation for tower
544 fTower->NTimeHits = 0;
545 sumWeightedTime = 0.0;
546 sumWeight = 0.0;
547
548 for(size_t i = 0; i < fTower->ECalEnergyTimePairs.size(); ++i)
549 {
550 weight = TMath::Sqrt(fTower->ECalEnergyTimePairs[i].first);
551 sumWeightedTime += weight * fTower->ECalEnergyTimePairs[i].second;
552 sumWeight += weight;
553 fTower->NTimeHits++;
554 }
555
556 if(sumWeight > 0.0)
557 {
558 fTower->Position.SetPtEtaPhiE(1.0, eta, phi, sumWeightedTime/sumWeight);
559 }
560 else
561 {
562 fTower->Position.SetPtEtaPhiE(1.0, eta, phi, 999999.9);
563 }
564
565
566 fTower->Momentum.SetPtEtaPhiE(pt, eta, phi, energy);
567 fTower->Eem = ecalEnergy;
568 fTower->Ehad = hcalEnergy;
569
570 fTower->Edges[0] = fTowerEdges[0];
571 fTower->Edges[1] = fTowerEdges[1];
572 fTower->Edges[2] = fTowerEdges[2];
573 fTower->Edges[3] = fTowerEdges[3];
574
575 if(energy > 0.0)
576 {
577 if(fTowerPhotonHits > 0 && fTowerTrackHits == 0)
578 {
579 fPhotonOutputArray->Add(fTower);
580 }
581
582 fTowerOutputArray->Add(fTower);
583 }
584
585 // fill energy flow candidates
586 fECalTrackSigma = TMath::Sqrt(fECalTrackSigma);
587 fHCalTrackSigma = TMath::Sqrt(fHCalTrackSigma);
588
589 //compute neutral excesses
590 ecalNeutralEnergy = max( (ecalEnergy - fECalTrackEnergy) , 0.0);
591 hcalNeutralEnergy = max( (hcalEnergy - fHCalTrackEnergy) , 0.0);
592
593 ecalNeutralSigma = ecalNeutralEnergy / TMath::Sqrt(fECalTrackSigma*fECalTrackSigma + ecalSigma*ecalSigma);
594 hcalNeutralSigma = hcalNeutralEnergy / TMath::Sqrt(fHCalTrackSigma*fHCalTrackSigma + hcalSigma*hcalSigma);
595
596 // if ecal neutral excess is significant, simply create neutral EflowPhoton tower and clone each track into eflowtrack
597 if(ecalNeutralEnergy > fECalEnergyMin && ecalNeutralSigma > fECalEnergySignificanceMin)
598 {
599 // create new photon tower
600 tower = static_cast<Candidate*>(fTower->Clone());
601 pt = ecalNeutralEnergy / TMath::CosH(eta);
602
603 tower->Momentum.SetPtEtaPhiE(pt, eta, phi, ecalNeutralEnergy);
604 tower->Eem = ecalNeutralEnergy;
605 tower->Ehad = 0.0;
606 tower->PID = 22;
607
608 fEFlowPhotonOutputArray->Add(tower);
609
610 //clone tracks
611 fItECalTowerTrackArray->Reset();
612 while((track = static_cast<Candidate*>(fItECalTowerTrackArray->Next())))
613 {
614 mother = track;
615 track = static_cast<Candidate*>(track->Clone());
616 track->AddCandidate(mother);
617
618 fEFlowTrackOutputArray->Add(track);
619 }
620
621 }
622
623 // if neutral excess is not significant, rescale eflow tracks, such that the total charged equals the best measurement given by the DualReadoutCalorimeter and tracking
624 else if(fECalTrackEnergy > 0.0)
625 {
626 weightTrack = (fECalTrackSigma > 0.0) ? 1 / (fECalTrackSigma*fECalTrackSigma) : 0.0;
627 weightCalo = (ecalSigma > 0.0) ? 1 / (ecalSigma*ecalSigma) : 0.0;
628
629 bestEnergyEstimate = (weightTrack*fECalTrackEnergy + weightCalo*ecalEnergy) / (weightTrack + weightCalo);
630 rescaleFactor = bestEnergyEstimate/fECalTrackEnergy;
631
632 //rescale tracks
633 fItECalTowerTrackArray->Reset();
634 while((track = static_cast<Candidate*>(fItECalTowerTrackArray->Next())))
635 {
636 mother = track;
637 track = static_cast<Candidate*>(track->Clone());
638 track->AddCandidate(mother);
639
640 track->Momentum *= rescaleFactor;
641
642 fEFlowTrackOutputArray->Add(track);
643 }
644 }
645
646
647 // if hcal neutral excess is significant, simply create neutral EflowNeutralHadron tower and clone each track into eflowtrack
648 if(hcalNeutralEnergy > fHCalEnergyMin && hcalNeutralSigma > fHCalEnergySignificanceMin)
649 {
650 // create new photon tower
651 tower = static_cast<Candidate*>(fTower->Clone());
652 pt = hcalNeutralEnergy / TMath::CosH(eta);
653
654 tower->Momentum.SetPtEtaPhiE(pt, eta, phi, hcalNeutralEnergy);
655 tower->Ehad = hcalNeutralEnergy;
656 tower->Eem = 0.0;
657
658 fEFlowNeutralHadronOutputArray->Add(tower);
659
660 //clone tracks
661 fItHCalTowerTrackArray->Reset();
662 while((track = static_cast<Candidate*>(fItHCalTowerTrackArray->Next())))
663 {
664 mother = track;
665 track = static_cast<Candidate*>(track->Clone());
666 track->AddCandidate(mother);
667
668 fEFlowTrackOutputArray->Add(track);
669 }
670
671 }
672
673 // if neutral excess is not significant, rescale eflow tracks, such that the total charged equals the best measurement given by the DualReadoutCalorimeter and tracking
674 else if(fHCalTrackEnergy > 0.0)
675 {
676 weightTrack = (fHCalTrackSigma > 0.0) ? 1 / (fHCalTrackSigma*fHCalTrackSigma) : 0.0;
677 weightCalo = (hcalSigma > 0.0) ? 1 / (hcalSigma*hcalSigma) : 0.0;
678
679 bestEnergyEstimate = (weightTrack*fHCalTrackEnergy + weightCalo*hcalEnergy) / (weightTrack + weightCalo);
680 rescaleFactor = bestEnergyEstimate / fHCalTrackEnergy;
681
682 //rescale tracks
683 fItHCalTowerTrackArray->Reset();
684 while((track = static_cast<Candidate*>(fItHCalTowerTrackArray->Next())))
685 {
686 mother = track;
687 track = static_cast<Candidate*>(track->Clone());
688 track->AddCandidate(mother);
689
690 track->Momentum *= rescaleFactor;
691
692 fEFlowTrackOutputArray->Add(track);
693 }
694 }
695
696
697}
698
699//------------------------------------------------------------------------------
700
701Double_t DualReadoutCalorimeter::LogNormal(Double_t mean, Double_t sigma)
702{
703 Double_t a, b;
704
705 if(mean > 0.0)
706 {
707 b = TMath::Sqrt(TMath::Log((1.0 + (sigma*sigma)/(mean*mean))));
708 a = TMath::Log(mean) - 0.5*b*b;
709
710 return TMath::Exp(a + b*gRandom->Gaus(0.0, 1.0));
711 }
712 else
713 {
714 return 0.0;
715 }
716}
Note: See TracBrowser for help on using the repository browser.