Fork me on GitHub

Context Navigation

source: git/modules/Calorimeter.cc@ 3db5282

Visit:

ImprovedOutputFile Timing dual_readout llp

Last change on this file since 3db5282 was 3db5282, checked in by Michele Selvaggi <michele.selvaggi@…>, 10 years ago
included timing information in Calorimeter
Property mode set to `100644`
File size: 18.2 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 Calorimeter
21	*
22	* Fills calorimeter towers, performs calorimeter 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/Calorimeter.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
52	using namespace std;
53
54	//------------------------------------------------------------------------------
55
56	Calorimeter::Calorimeter() :
57	fECalResolutionFormula(0), fHCalResolutionFormula(0),
58	fItParticleInputArray(0), fItTrackInputArray(0),
59	fTowerTrackArray(0), fItTowerTrackArray(0)
60	{
61	fECalResolutionFormula = new DelphesFormula;
62	fHCalResolutionFormula = new DelphesFormula;
63
64	fTowerTrackArray = new TObjArray;
65	fItTowerTrackArray = fTowerTrackArray->MakeIterator();
66	}
67
68	//------------------------------------------------------------------------------
69
70	Calorimeter::~Calorimeter()
71	{
72	if(fECalResolutionFormula) delete fECalResolutionFormula;
73	if(fHCalResolutionFormula) delete fHCalResolutionFormula;
74
75	if(fTowerTrackArray) delete fTowerTrackArray;
76	if(fItTowerTrackArray) delete fItTowerTrackArray;
77	}
78
79	//------------------------------------------------------------------------------
80
81	void Calorimeter::Init()
82	{
83	ExRootConfParam param, paramEtaBins, paramPhiBins, paramFractions;
84	Long_t i, j, k, size, sizeEtaBins, sizePhiBins;
85	Double_t ecalFraction, hcalFraction;
86	TBinMap::iterator itEtaBin;
87	set< Double_t >::iterator itPhiBin;
88	vector< Double_t > *phiBins;
89
90	// read eta and phi bins
91	param = GetParam("EtaPhiBins");
92	size = param.GetSize();
93	fBinMap.clear();
94	fEtaBins.clear();
95	fPhiBins.clear();
96	for(i = 0; i < size/2; ++i)
97	{
98	paramEtaBins = param[i*2];
99	sizeEtaBins = paramEtaBins.GetSize();
100	paramPhiBins = param[i*2 + 1];
101	sizePhiBins = paramPhiBins.GetSize();
102
103	for(j = 0; j < sizeEtaBins; ++j)
104	{
105	for(k = 0; k < sizePhiBins; ++k)
106	{
107	fBinMap[paramEtaBins[j].GetDouble()].insert(paramPhiBins[k].GetDouble());
108	}
109	}
110	}
111
112	// for better performance we transform map of sets to parallel vectors:
113	// vector< double > and vector< vector< double >* >
114	for(itEtaBin = fBinMap.begin(); itEtaBin != fBinMap.end(); ++itEtaBin)
115	{
116	fEtaBins.push_back(itEtaBin->first);
117	phiBins = new vector< double >(itEtaBin->second.size());
118	fPhiBins.push_back(phiBins);
119	phiBins->clear();
120	for(itPhiBin = itEtaBin->second.begin(); itPhiBin != itEtaBin->second.end(); ++itPhiBin)
121	{
122	phiBins->push_back(*itPhiBin);
123	}
124	}
125
126	// read energy fractions for different particles
127	param = GetParam("EnergyFraction");
128	size = param.GetSize();
129
130	// set default energy fractions values
131	fFractionMap.clear();
132	fFractionMap[0] = make_pair(0.0, 1.0);
133
134	for(i = 0; i < size/2; ++i)
135	{
136	paramFractions = param[i*2 + 1];
137
138	ecalFraction = paramFractions[0].GetDouble();
139	hcalFraction = paramFractions[1].GetDouble();
140
141	fFractionMap[param[i*2].GetInt()] = make_pair(ecalFraction, hcalFraction);
142	}
143
144	/*
145	TFractionMap::iterator itFractionMap;
146	for(itFractionMap = fFractionMap.begin(); itFractionMap != fFractionMap.end(); ++itFractionMap)
147	{
148	cout << itFractionMap->first << " " << itFractionMap->second.first << " " << itFractionMap->second.second << endl;
149	}
150	*/
151
152	// read min E value for timing measurement in ECAL
153	fTimingEMin = GetDouble("TimingEMin",4.);
154	// For timing
155	// So far this flag needs to be false
156	// Curved extrapolation not supported
157	fElectronsFromTrack = false;
158
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 calorimeter 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
192	void Calorimeter::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
205	void Calorimeter::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	Int_t pdgCode;
215
216	TFractionMap::iterator itFractionMap;
217
218	vector< Double_t >::iterator itEtaBin;
219	vector< Double_t >::iterator itPhiBin;
220	vector< Double_t > *phiBins;
221
222	vector< Long64_t >::iterator itTowerHits;
223
224	DelphesFactory *factory = GetFactory();
225	fTowerHits.clear();
226	fTowerECalFractions.clear();
227	fTowerHCalFractions.clear();
228	fTrackECalFractions.clear();
229	fTrackHCalFractions.clear();
230
231	// loop over all particles
232	fItParticleInputArray->Reset();
233	number = -1;
234	while((particle = static_cast<Candidate*>(fItParticleInputArray->Next())))
235	{
236	const TLorentzVector &particlePosition = particle->Position;
237	++number;
238
239	pdgCode = TMath::Abs(particle->PID);
240
241	itFractionMap = fFractionMap.find(pdgCode);
242	if(itFractionMap == fFractionMap.end())
243	{
244	itFractionMap = fFractionMap.find(0);
245	}
246
247	ecalFraction = itFractionMap->second.first;
248	hcalFraction = itFractionMap->second.second;
249
250	fTowerECalFractions.push_back(ecalFraction);
251	fTowerHCalFractions.push_back(hcalFraction);
252
253	if(ecalFraction < 1.0E-9 && hcalFraction < 1.0E-9) continue;
254
255	// find eta bin [1, fEtaBins.size - 1]
256	itEtaBin = lower_bound(fEtaBins.begin(), fEtaBins.end(), particlePosition.Eta());
257	if(itEtaBin == fEtaBins.begin() \|\| itEtaBin == fEtaBins.end()) continue;
258	etaBin = distance(fEtaBins.begin(), itEtaBin);
259
260	// phi bins for given eta bin
261	phiBins = fPhiBins[etaBin];
262
263	// find phi bin [1, phiBins.size - 1]
264	itPhiBin = lower_bound(phiBins->begin(), phiBins->end(), particlePosition.Phi());
265	if(itPhiBin == phiBins->begin() \|\| itPhiBin == phiBins->end()) continue;
266	phiBin = distance(phiBins->begin(), itPhiBin);
267
268	flags = 0;
269	flags \|= (pdgCode == 11 \|\| pdgCode == 22) << 1;
270
271	// make tower hit {16-bits for eta bin number, 16-bits for phi bin number, 8-bits for flags, 24-bits for particle number}
272	towerHit = (Long64_t(etaBin) << 48) \| (Long64_t(phiBin) << 32) \| (Long64_t(flags) << 24) \| Long64_t(number);
273
274	fTowerHits.push_back(towerHit);
275	}
276
277	// loop over all tracks
278	fItTrackInputArray->Reset();
279	number = -1;
280	while((track = static_cast<Candidate*>(fItTrackInputArray->Next())))
281	{
282	const TLorentzVector &trackPosition = track->Position;
283	++number;
284
285	pdgCode = TMath::Abs(track->PID);
286
287	itFractionMap = fFractionMap.find(pdgCode);
288	if(itFractionMap == fFractionMap.end())
289	{
290	itFractionMap = fFractionMap.find(0);
291	}
292
293	ecalFraction = itFractionMap->second.first;
294	hcalFraction = itFractionMap->second.second;
295
296	fTrackECalFractions.push_back(ecalFraction);
297	fTrackHCalFractions.push_back(hcalFraction);
298
299	// find eta bin [1, fEtaBins.size - 1]
300	itEtaBin = lower_bound(fEtaBins.begin(), fEtaBins.end(), trackPosition.Eta());
301	if(itEtaBin == fEtaBins.begin() \|\| itEtaBin == fEtaBins.end()) continue;
302	etaBin = distance(fEtaBins.begin(), itEtaBin);
303
304	// phi bins for given eta bin
305	phiBins = fPhiBins[etaBin];
306
307	// find phi bin [1, phiBins.size - 1]
308	itPhiBin = lower_bound(phiBins->begin(), phiBins->end(), trackPosition.Phi());
309	if(itPhiBin == phiBins->begin() \|\| itPhiBin == phiBins->end()) continue;
310	phiBin = distance(phiBins->begin(), itPhiBin);
311
312	flags = 1;
313
314	// make tower hit {16-bits for eta bin number, 16-bits for phi bin number, 8-bits for flags, 24-bits for track number}
315	towerHit = (Long64_t(etaBin) << 48) \| (Long64_t(phiBin) << 32) \| (Long64_t(flags) << 24) \| Long64_t(number);
316
317	fTowerHits.push_back(towerHit);
318	}
319
320	// all hits are sorted first by eta bin number, then by phi bin number,
321	// then by flags and then by particle or track number
322	sort(fTowerHits.begin(), fTowerHits.end());
323
324	// loop over all hits
325	towerEtaPhi = 0;
326	fTower = 0;
327	for(itTowerHits = fTowerHits.begin(); itTowerHits != fTowerHits.end(); ++itTowerHits)
328	{
329	towerHit = (*itTowerHits);
330	flags = (towerHit >> 24) & 0x00000000000000FFLL;
331	number = (towerHit) & 0x0000000000FFFFFFLL;
332	hitEtaPhi = towerHit >> 32;
333
334	if(towerEtaPhi != hitEtaPhi)
335	{
336	// switch to next tower
337	towerEtaPhi = hitEtaPhi;
338
339	// finalize previous tower
340	FinalizeTower();
341
342	// create new tower
343	fTower = factory->NewCandidate();
344
345	phiBin = (towerHit >> 32) & 0x000000000000FFFFLL;
346	etaBin = (towerHit >> 48) & 0x000000000000FFFFLL;
347
348	// phi bins for given eta bin
349	phiBins = fPhiBins[etaBin];
350
351	// calculate eta and phi of the tower's center
352	fTowerEta = 0.5*(fEtaBins[etaBin - 1] + fEtaBins[etaBin]);
353	fTowerPhi = 0.5((phiBins)[phiBin - 1] + (*phiBins)[phiBin]);
354
355	fTowerEdges[0] = fEtaBins[etaBin - 1];
356	fTowerEdges[1] = fEtaBins[etaBin];
357	fTowerEdges[2] = (*phiBins)[phiBin - 1];
358	fTowerEdges[3] = (*phiBins)[phiBin];
359
360	fTowerECalEnergy = 0.0;
361	fTowerHCalEnergy = 0.0;
362
363	fTrackECalEnergy = 0.0;
364	fTrackHCalEnergy = 0.0;
365
366	fTowerTrackHits = 0;
367	fTowerPhotonHits = 0;
368
369	fTowerTrackArray->Clear();
370	}
371
372	// check for track hits
373	if(flags & 1)
374	{
375	++fTowerTrackHits;
376
377	track = static_cast<Candidate*>(fTrackInputArray->At(number));
378	momentum = track->Momentum;
379	position = track->Position;
380
381	ecalEnergy = momentum.E() * fTrackECalFractions[number];
382	hcalEnergy = momentum.E() * fTrackHCalFractions[number];
383
384	fTrackECalEnergy += ecalEnergy;
385	fTrackHCalEnergy += hcalEnergy;
386
387	bool dbg_scz = false;
388	if (dbg_scz) {
389	cout << " Calorimeter input track has x y z t " << track->Position.X() << " " << track->Position.Y() << " " << track->Position.Z() << " " << track->Position.T()
390	<< endl;
391	Candidate prt = static_cast<Candidate>(track->GetCandidates()->Last());
392	const TLorentzVector &ini = prt->Position;
393
394	cout << " and parent has x y z t " << ini.X() << " " << ini.Y() << " " << ini.Z() << " " << ini.T();
395
396	}
397
398	if (ecalEnergy > fTimingEMin && fTower) {
399	if (fElectronsFromTrack) {
400	// cout << " SCZ Debug pushing back track hit E=" << ecalEnergy << " T=" << track->Position.T() << " isPU=" << track->IsPU << " isRecoPU=" << track->IsRecoPU
401	// << " PID=" << track->PID << endl;
402	fTower->Ecal_E_t.push_back(std::make_pair<float,float>(ecalEnergy,track->Position.T()));
403	} else {
404	// cout << " Skipping track hit E=" << ecalEnergy << " T=" << track->Position.T() << " isPU=" << track->IsPU << " isRecoPU=" << track->IsRecoPU
405	// << " PID=" << track->PID << endl;
406	}
407	}
408
409
410	fTowerTrackArray->Add(track);
411
412	continue;
413	}
414
415	// check for photon and electron hits in current tower
416	if(flags & 2) ++fTowerPhotonHits;
417
418	particle = static_cast<Candidate*>(fParticleInputArray->At(number));
419	momentum = particle->Momentum;
420	position = particle->Position;
421
422	// fill current tower
423	ecalEnergy = momentum.E() * fTowerECalFractions[number];
424	hcalEnergy = momentum.E() * fTowerHCalFractions[number];
425
426	fTowerECalEnergy += ecalEnergy;
427	fTowerHCalEnergy += hcalEnergy;
428
429	if (ecalEnergy > fTimingEMin && fTower) {
430	if (abs(particle->PID) != 11 \|\| !fElectronsFromTrack) {
431	// cout << " SCZ Debug About to push back particle hit E=" << ecalEnergy << " T=" << particle->Position.T() << " isPU=" << particle->IsPU
432	// << " PID=" << particle->PID << endl;
433	fTower->Ecal_E_t.push_back(std::make_pair<Float_t,Float_t>(ecalEnergy,particle->Position.T()));
434	} else {
435
436	// N.B. Only charged particle set to leave ecal energy is the electrons
437	// cout << " SCZ Debug To avoid double-counting, skipping particle hit E=" << ecalEnergy << " T=" << particle->Position.T() << " isPU=" << particle->IsPU
438	// << " PID=" << particle->PID << endl;
439
440	}
441	}
442
443	fTower->AddCandidate(particle);
444	}
445
446	// finalize last tower
447	FinalizeTower();
448	}
449
450	//------------------------------------------------------------------------------
451
452	void Calorimeter::FinalizeTower()
453	{
454	Candidate track, tower;
455	Double_t energy, pt, eta, phi;
456	Double_t ecalEnergy, hcalEnergy;
457	Double_t ecalSigma, hcalSigma;
458
459	if(!fTower) return;
460
461	ecalSigma = fECalResolutionFormula->Eval(0.0, fTowerEta, 0.0, fTowerECalEnergy);
462	hcalSigma = fHCalResolutionFormula->Eval(0.0, fTowerEta, 0.0, fTowerHCalEnergy);
463
464	ecalEnergy = LogNormal(fTowerECalEnergy, ecalSigma);
465	hcalEnergy = LogNormal(fTowerHCalEnergy, hcalSigma);
466
467	ecalSigma = fECalResolutionFormula->Eval(0.0, fTowerEta, 0.0, ecalEnergy);
468	hcalSigma = fHCalResolutionFormula->Eval(0.0, fTowerEta, 0.0, hcalEnergy);
469
470	if(ecalEnergy < fECalEnergyMin \|\| ecalEnergy < fECalEnergySignificanceMin*ecalSigma) ecalEnergy = 0.0;
471	if(hcalEnergy < fHCalEnergyMin \|\| hcalEnergy < fHCalEnergySignificanceMin*hcalSigma) hcalEnergy = 0.0;
472
473	energy = ecalEnergy + hcalEnergy;
474
475	if(fSmearTowerCenter)
476	{
477	eta = gRandom->Uniform(fTowerEdges[0], fTowerEdges[1]);
478	phi = gRandom->Uniform(fTowerEdges[2], fTowerEdges[3]);
479	}
480	else
481	{
482	eta = fTowerEta;
483	phi = fTowerPhi;
484	}
485
486	pt = energy / TMath::CosH(eta);
487
488	// Time calculation for tower
489	fTower->Ntimes = 0;
490	Float_t tow_sumT = 0;
491	Float_t tow_sumW = 0;
492
493	for (Int_t i = 0 ; i < fTower->Ecal_E_t.size() ; i++)
494	{
495	Float_t w = TMath::Sqrt(fTower->Ecal_E_t[i].first);
496	tow_sumT += w*fTower->Ecal_E_t[i].second;
497	tow_sumW += w;
498	fTower->Ntimes++;
499	}
500
501	if (tow_sumW > 0.) {
502	fTower->Position.SetPtEtaPhiE(1.0, eta, phi,tow_sumT/tow_sumW);
503	} else {
504	fTower->Position.SetPtEtaPhiE(1.0,eta,phi,999999.);
505	}
506
507
508	fTower->Momentum.SetPtEtaPhiE(pt, eta, phi, energy);
509	fTower->Eem = ecalEnergy;
510	fTower->Ehad = hcalEnergy;
511
512	fTower->Edges[0] = fTowerEdges[0];
513	fTower->Edges[1] = fTowerEdges[1];
514	fTower->Edges[2] = fTowerEdges[2];
515	fTower->Edges[3] = fTowerEdges[3];
516
517	if(energy > 0.0)
518	{
519	if(fTowerPhotonHits > 0 && fTowerTrackHits == 0)
520	{
521	fPhotonOutputArray->Add(fTower);
522	}
523
524	fTowerOutputArray->Add(fTower);
525	}
526
527	// fill energy flow candidates
528
529	// save all the tracks as energy flow tracks
530	fItTowerTrackArray->Reset();
531	while((track = static_cast<Candidate*>(fItTowerTrackArray->Next())))
532	{
533	fEFlowTrackOutputArray->Add(track);
534	}
535
536	ecalEnergy -= fTrackECalEnergy;
537	hcalEnergy -= fTrackHCalEnergy;
538
539	ecalSigma = fECalResolutionFormula->Eval(0.0, fTowerEta, 0.0, ecalEnergy);
540	hcalSigma = fHCalResolutionFormula->Eval(0.0, fTowerEta, 0.0, hcalEnergy);
541
542	if(ecalEnergy < fECalEnergyMin \|\| ecalEnergy < fECalEnergySignificanceMin*ecalSigma) ecalEnergy = 0.0;
543	if(hcalEnergy < fHCalEnergyMin \|\| hcalEnergy < fHCalEnergySignificanceMin*hcalSigma) hcalEnergy = 0.0;
544
545	energy = ecalEnergy + hcalEnergy;
546
547	if(ecalEnergy > 0.0)
548	{
549	// create new photon tower
550	tower = static_cast<Candidate*>(fTower->Clone());
551
552	pt = ecalEnergy / TMath::CosH(eta);
553
554	tower->Momentum.SetPtEtaPhiE(pt, eta, phi, ecalEnergy);
555	tower->Eem = ecalEnergy;
556	tower->Ehad = 0.0;
557
558	fEFlowPhotonOutputArray->Add(tower);
559	}
560	if(hcalEnergy > 0.0)
561	{
562	// create new neutral hadron tower
563	tower = static_cast<Candidate*>(fTower->Clone());
564
565	pt = hcalEnergy / TMath::CosH(eta);
566
567	tower->Momentum.SetPtEtaPhiE(pt, eta, phi, hcalEnergy);
568	tower->Eem = 0.0;
569	tower->Ehad = hcalEnergy;
570
571	fEFlowNeutralHadronOutputArray->Add(tower);
572	}
573	}
574
575	//------------------------------------------------------------------------------
576
577	Double_t Calorimeter::LogNormal(Double_t mean, Double_t sigma)
578	{
579	Double_t a, b;
580
581	if(mean > 0.0)
582	{
583	b = TMath::Sqrt(TMath::Log((1.0 + (sigmasigma)/(meanmean))));
584	a = TMath::Log(mean) - 0.5bb;
585
586	return TMath::Exp(a + b*gRandom->Gaus(0.0, 1.0));
587	}
588	else
589	{
590	return 0.0;
591	}
592	}

Note: See TracBrowser for help on using the repository browser.

Download in other formats: