Fork me on GitHub

source: svn/trunk/src/CaloUtil.cc@ 898

Last change on this file since 898 was 551, checked in by Xavier Rouby, 15 years ago

EtRatio for muon and electron; apparently bug free
File size: 9.1 KB
Line 
1/***********************************************************************
2** **
3** /----------------------------------------------\ **
4** | Delphes, a framework for the fast simulation | **
5** | of a generic collider experiment | **
6** \------------- arXiv:0903.2225v1 ------------/ **
7** **
8** **
9** This package uses: **
10** ------------------ **
11** ROOT: Nucl. Inst. & Meth. in Phys. Res. A389 (1997) 81-86 **
12** FastJet algorithm: Phys. Lett. B641 (2006) [hep-ph/0512210] **
13** Hector: JINST 2:P09005 (2007) [physics.acc-ph:0707.1198v2] **
14** FROG: [hep-ex/0901.2718v1] **
15** HepMC: Comput. Phys. Commun.134 (2001) 41 **
16** **
17** ------------------------------------------------------------------ **
18** **
19** Main authors: **
20** ------------- **
21** **
22** Severine Ovyn Xavier Rouby **
23** severine.ovyn@uclouvain.be xavier.rouby@cern **
24** **
25** Center for Particle Physics and Phenomenology (CP3) **
26** Universite catholique de Louvain (UCL) **
27** Louvain-la-Neuve, Belgium **
28** **
29** Copyright (C) 2008-2009, **
30** All rights reserved. **
31** **
32***********************************************************************/
33
34/// \file D_CaloUtil.cc
35/// \brief Implementation of D_CaloElement, D_CaloResolution and D_CaloList classes
36
37#include "CaloUtil.h"
38#include <cmath>
39#include <algorithm>
40#include <iostream>
41#include "TRandom.h"
42
43D_CaloElement& D_CaloElement::operator=(const D_CaloElement& calo) {
44 if(this == &calo) return *this;
45 _name = calo._name;
46 _etamin = calo._etamin;
47 _etamax = calo._etamax;
48 _resol_em = calo._resol_em;
49 _resol_had= calo._resol_had;
50 return *this;
51}
52
53D_CaloResolution& D_CaloResolution::operator=(const D_CaloResolution& resolution) {
54 if(this == &resolution) return *this;
55 _c = resolution._c;
56 _n = resolution._n;
57 _s = resolution._s;
58 return *this;
59}
60
61void D_CaloResolution::print() const {
62 std::cout << "C: " << _c << " N:" << _n << " S:" << _s << std::endl;
63}
64
65const D_CaloElement& D_CaloList::getElement(const float eta) const {
66 for (unsigned int i=0; i<_calorimeters.size(); i++) {
67 if(_calorimeters[i].getEtamin() < eta && eta < _calorimeters[i].getEtamax()) return _calorimeters[i];
68 }
69 // if not found
70 return dummyCalo;
71}
72
73void D_CaloList::addElement(const D_CaloElement & calo) {
74 _sorted = false;
75 _calorimeters.push_back(calo);
76}
77
78void D_CaloList::sortElements() {
79 // here: it sorts the elements with respect to their eta
80 sort(_calorimeters.begin(),_calorimeters.end(),ordering());
81
82 for (unsigned int i=0; i<_calorimeters.size()-1; i++) {
83 if(_calorimeters[i].getEtamax() > _calorimeters[i+1].getEtamin() ) {
84 _etamin = UNDEFINED;
85 _etamax = UNDEFINED;
86 cerr << "** Error: the calorimeter list should never contain overlapping **"<< endl;
87 cerr << "** calorimeters. **"<< endl;
88 return;
89 }
90 }
91 _etamin = _calorimeters[0].getEtamin();
92 _etamax = _calorimeters[_calorimeters.size()-1].getEtamax();
93 _sorted = true;
94}
95
96D_CaloList& D_CaloList::operator=(const D_CaloList& list) {
97 if(this == &list) return *this;
98 _calorimeters = list._calorimeters;
99 _etamin = list._etamin;
100 _etamax = list._etamax;
101 _sorted = list._sorted;
102 return *this;
103}
104
105void D_CaloList::print() const {
106 if(_sorted) cout << "The list has been sorted\n";
107 else cout << "The list has not been sorted\n";
108
109 for (unsigned int i=0; i<_calorimeters.size(); i++) {
110 std::cout << _calorimeters[i].getName() << " " << _calorimeters[i].getEtamin() <<
111 " " << _calorimeters[i].getEtamax() << std::endl;
112 _calorimeters[i].getElectromagneticResolution().print();
113 _calorimeters[i].getHadronicResolution().print();
114 }
115
116}
117
118
119
120
121
122
123
124
125D_CaloTower::D_CaloTower(const float iEta, const float iPhi) :
126 _eta(iEta), _phi(iPhi), _Eem(UNDEFINED), _Ehad(UNDEFINED), _E(UNDEFINED), _ET(UNDEFINED) {}
127
128D_CaloTower::D_CaloTower(const D_CaloTower& tower) :
129 _eta(tower._eta), _phi(tower._phi), _Eem(tower._Eem),
130 _Ehad(tower._Ehad), _E(tower._E), _ET(tower._ET) {}
131
132D_CaloTower& D_CaloTower::operator=(const D_CaloTower& tower) {
133 if(this==&tower) return *this;
134 _eta = tower._eta; _phi = tower._phi;
135 _Eem = tower._Eem; _Ehad = tower._Ehad;
136 _E = tower._E; _ET = tower._ET;
137 return *this;
138}
139
140
141const float D_CaloResolution::Smear(const float energy) const {
142 float sE = gRandom->Gaus(energy,
143 sqrt( pow(_n,2) +
144 pow(_c*energy,2) +
145 pow(_s*sqrt(energy),2) ));
146 return (sE>0)? sE : 0;
147}
148
149const bool D_CaloTower::operator==(const D_CaloTower& el) const { if (this == &el) return true; return this->_eta==el._eta && this->_phi == el._phi;}
150
151void D_CaloTowerList::addTower(const D_CaloTower& tower) {
152 _calotowers.push_back(tower);
153 _sorted = false;
154 _merged = false;
155}
156
157D_CaloTowerList& D_CaloTowerList::operator=(const D_CaloTowerList& list) {
158 if(this==&list) return *this;
159 _calotowers = list._calotowers;
160 _sorted = list._sorted;
161 _merged = list._merged;
162 return *this;
163}
164
165void D_CaloTowerList::sortElements() {
166 _sorted = true;
167 if(_calotowers.size()==0) return;
168 sort(_calotowers.begin(),_calotowers.end(),ordering());
169}
170
171void D_CaloTowerList::mergeDuplicates() {
172 if(!_sorted) sortElements();
173 _merged = true;
174 if (_calotowers.size() ==0) { return;}
175
176 vector<D_CaloTower> temp;
177 temp.push_back(_calotowers[0]);
178 for (unsigned int i=1; i<_calotowers.size(); i++) {
179
180 if(temp[temp.size()-1] == _calotowers[i]) {
181 //i.e. if towers have the same eta/phi
182 float em = temp[temp.size()-1].getEem() + _calotowers[i].getEem();
183 float had= temp[temp.size()-1].getEhad()+ _calotowers[i].getEhad();
184 temp[temp.size()-1].Set_Eem_Ehad_E_ET(em,had);
185 }
186 else temp.push_back(_calotowers[i]);
187
188 }
189 _calotowers = temp;
190}
191
192void D_CaloTowerList::print() const {
193 for (unsigned int i=0; i<_calotowers.size(); i++) {
194 std::cout << _calotowers[i].getEta() << " " <<
195 _calotowers[i].getPhi() << " " <<
196 _calotowers[i].getEem() << " " <<
197 _calotowers[i].getEhad() << std::endl;
198 }
199}
200
201void D_CaloTowerList::smearTowers(const D_CaloList& list_of_calorimeters) {
202 for (unsigned int i=0; i<_calotowers.size(); i++) {
203 // first, finds in which calo the tower is, in order to have its resolution em /had
204 D_CaloElement currentCalo = list_of_calorimeters.getElement(_calotowers[i].getEta() );
205 // then, smears separately the em and had parts
206 float sEem = currentCalo.getElectromagneticResolution().Smear( _calotowers[i].getEem() );
207 float sEhad = currentCalo.getHadronicResolution().Smear( _calotowers[i].getEhad() );
208 // finally, sets the tower smeared energy
209 _calotowers[i].Set_Eem_Ehad_E_ET(sEem,sEhad);
210 }
211}
212
213
214const bool D_CaloTower::operator>(const D_CaloTower& tocomp) const {
215 if(_eta == tocomp._eta) {
216 if(_phi > tocomp._phi) return true;
217 else return false;
218 }
219 else if(_eta > tocomp._eta) return true;
220 else return false;
221};
222
223const bool D_CaloTower::operator<(const D_CaloTower& tocomp) const {
224 if(_eta == tocomp._eta) {
225 if(_phi < tocomp._phi) return true;
226 else return false;
227 }
228 else if(_eta < tocomp._eta) return true;
229 else return false;
230};
231
232const D_CaloTower& D_CaloTowerList::getElement(const float eta, const float phi) {
233 if(!_sorted) sortElements();
234 for (unsigned int i=0; i< _calotowers.size(); i++) {
235 if(_calotowers[i].getEta() == eta &&
236 _calotowers[i].getPhi() == phi )
237 return _calotowers[i];
238 }
239 return dummyTower; // if not found
240}
241
242const D_CaloTower& D_CaloTowerList::getElement(const float eta, const float phi) const {
243 if(!_sorted) cout << "Warning: unsorted list!";
244 for (unsigned int i=0; i< _calotowers.size(); i++) {
245 if(_calotowers[i].getEta() == eta &&
246 _calotowers[i].getPhi() == phi )
247 return _calotowers[i];
248 }
249 return dummyTower; // if not found
250}
251
Note: See TracBrowser for help on using the repository browser.