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 | /** \class EnergyLoss
|
---|
20 | *
|
---|
21 | * This module computes the charged energy loss according to the active material properties.
|
---|
22 | * The energy loss is simulated with a Landau convoluted by a Gaussian.
|
---|
23 | *
|
---|
24 | * \author M. Selvaggi - CERN
|
---|
25 | *
|
---|
26 | */
|
---|
27 | #include "modules/EnergyLoss.h"
|
---|
28 |
|
---|
29 | #include "classes/DelphesClasses.h"
|
---|
30 | #include "classes/DelphesFactory.h"
|
---|
31 | #include "classes/DelphesFormula.h"
|
---|
32 |
|
---|
33 | #include "ExRootAnalysis/ExRootClassifier.h"
|
---|
34 | #include "ExRootAnalysis/ExRootFilter.h"
|
---|
35 | #include "ExRootAnalysis/ExRootResult.h"
|
---|
36 |
|
---|
37 | #include "TDatabasePDG.h"
|
---|
38 | #include "TFormula.h"
|
---|
39 | #include "TLorentzVector.h"
|
---|
40 | #include "TMath.h"
|
---|
41 | #include "TObjArray.h"
|
---|
42 | #include "TRandom3.h"
|
---|
43 | #include "TString.h"
|
---|
44 |
|
---|
45 | #include <algorithm>
|
---|
46 | #include <iostream>
|
---|
47 | #include <sstream>
|
---|
48 | #include <stdexcept>
|
---|
49 | #include <numeric>
|
---|
50 |
|
---|
51 | using namespace std;
|
---|
52 |
|
---|
53 | //------------------------------------------------------------------------------
|
---|
54 |
|
---|
55 | EnergyLoss::EnergyLoss()
|
---|
56 | {
|
---|
57 | }
|
---|
58 |
|
---|
59 | //------------------------------------------------------------------------------
|
---|
60 |
|
---|
61 | EnergyLoss::~EnergyLoss()
|
---|
62 | {
|
---|
63 | }
|
---|
64 |
|
---|
65 | //------------------------------------------------------------------------------
|
---|
66 |
|
---|
67 | void EnergyLoss::Init()
|
---|
68 | {
|
---|
69 |
|
---|
70 | fActiveFraction = GetDouble("ActiveFraction", 0.002); // active fraction of the detector
|
---|
71 | fThickness = GetDouble("Thickness", 200E-6); // active detector thickness
|
---|
72 | fResolution = GetDouble("Resolution", 0.4); // 0 - perfect Landau energy loss (0.15 gives good agreement with CMS pixel detector)
|
---|
73 | fTruncatedMeanFraction = GetDouble("TruncatedMeanFraction", 0.5); // fraction of measurements to ignore when computing mean
|
---|
74 |
|
---|
75 | // active material properties (cf. http://pdg.lbl.gov/2014/AtomicNuclearProperties/properties8.dat)
|
---|
76 | fZ = GetDouble("Z", 14.);
|
---|
77 | fA = GetDouble("A", 28.0855); // in g/mol
|
---|
78 | fRho = GetDouble("rho", 2.329); // in g/cm3
|
---|
79 | fAa = GetDouble("a", 0.1492);
|
---|
80 | fM = GetDouble("m", 3.2546);
|
---|
81 | fX0 = GetDouble("x0", 0.2015);
|
---|
82 | fX1 = GetDouble("x1", 2.8716);
|
---|
83 | fI = GetDouble("I", 173.0); // mean excitation potential in (eV)
|
---|
84 | fC0 = GetDouble("c0", 4.4355);
|
---|
85 |
|
---|
86 | // import arrays with output from other modules
|
---|
87 |
|
---|
88 | ExRootConfParam param = GetParam("InputArray");
|
---|
89 | Long_t i, size;
|
---|
90 | const TObjArray *array;
|
---|
91 | TIterator *iterator;
|
---|
92 |
|
---|
93 | size = param.GetSize();
|
---|
94 | for(i = 0; i < size; ++i)
|
---|
95 | {
|
---|
96 | array = ImportArray(param[i].GetString());
|
---|
97 | iterator = array->MakeIterator();
|
---|
98 |
|
---|
99 | fInputList.push_back(iterator);
|
---|
100 | }
|
---|
101 |
|
---|
102 | }
|
---|
103 |
|
---|
104 | //------------------------------------------------------------------------------
|
---|
105 |
|
---|
106 | void EnergyLoss::Finish()
|
---|
107 | {
|
---|
108 | vector<TIterator *>::iterator itInputList;
|
---|
109 | TIterator *iterator;
|
---|
110 |
|
---|
111 | for(itInputList = fInputList.begin(); itInputList != fInputList.end(); ++itInputList)
|
---|
112 | {
|
---|
113 | iterator = *itInputList;
|
---|
114 | if(iterator) delete iterator;
|
---|
115 | }
|
---|
116 |
|
---|
117 | }
|
---|
118 |
|
---|
119 | //------------------------------------------------------------------------------
|
---|
120 |
|
---|
121 | void EnergyLoss::Process()
|
---|
122 | {
|
---|
123 | Candidate *candidate;
|
---|
124 | vector<TIterator *>::iterator itInputList;
|
---|
125 | TIterator *iterator;
|
---|
126 |
|
---|
127 | Double_t beta, gamma, charge;
|
---|
128 | Double_t kappa, chi, me, I, Wmax, delta, avdE, dP, dx, L, dE, dEdx, res;
|
---|
129 | Double_t eloss_truncmean;
|
---|
130 |
|
---|
131 | Int_t nhits;
|
---|
132 | vector<Double_t> elosses;
|
---|
133 |
|
---|
134 | //cout<<"---------------- new event -------------------"<<endl;
|
---|
135 |
|
---|
136 |
|
---|
137 | // loop over all input arrays
|
---|
138 | for(itInputList = fInputList.begin(); itInputList != fInputList.end(); ++itInputList)
|
---|
139 | {
|
---|
140 | iterator = *itInputList;
|
---|
141 |
|
---|
142 | // loop over all candidates
|
---|
143 | iterator->Reset();
|
---|
144 | while((candidate = static_cast<Candidate *>(iterator->Next())))
|
---|
145 | {
|
---|
146 | //cout<<" ---------------- new candidate -------------------"<<endl;
|
---|
147 | const TLorentzVector &candidateMomentum = candidate->Momentum;
|
---|
148 |
|
---|
149 | beta = candidateMomentum.Beta();
|
---|
150 | gamma = candidateMomentum.Gamma();
|
---|
151 | charge = TMath::Abs(candidate->Charge);
|
---|
152 |
|
---|
153 |
|
---|
154 | // length of the track normalized by the fraction of active material and the charge collection efficiency in the tracker (in cm)
|
---|
155 | //dx = candidate->L * fActiveFraction * 0.1;
|
---|
156 | // amount of material in one sensor (converted in cm)
|
---|
157 | dx = fThickness * 100.;
|
---|
158 | // path length in cm
|
---|
159 | L = candidate->L * 0.1;
|
---|
160 |
|
---|
161 |
|
---|
162 | // compute number of hits as path length over active length
|
---|
163 | nhits = Int_t(L*fActiveFraction/dx);
|
---|
164 |
|
---|
165 |
|
---|
166 | //beta = 0.999945;
|
---|
167 | //gamma = 95.6446;
|
---|
168 | //charge = 1.;
|
---|
169 | //nhits = 100;
|
---|
170 |
|
---|
171 | //cout<<L<<","<<fActiveFraction<<","<<dx<<","<<nhits<<endl;
|
---|
172 |
|
---|
173 | kappa = 2*0.1535*TMath::Abs(charge)*TMath::Abs(charge)*fZ*fRho*dx/(fA*beta*beta); //energy loss in MeV
|
---|
174 |
|
---|
175 | chi = 0.5*kappa;
|
---|
176 | me = 0.510998; // electron mass in MeV, need
|
---|
177 | I = fI*1e-6; // convert I in MeV
|
---|
178 |
|
---|
179 | // fixme: max energy transfer wrong for electrons
|
---|
180 | Wmax = 2*me*beta*beta*gamma*gamma; // this is not valid for electrons
|
---|
181 |
|
---|
182 | delta = Deltaf(fC0, fAa, fM, fX0, fX1, beta, gamma);
|
---|
183 |
|
---|
184 | // Bethe-Bloch energy loss in MeV (not used here)
|
---|
185 | avdE = kappa*( TMath::Log(Wmax/I) - beta*beta - delta/2);
|
---|
186 |
|
---|
187 | // most probable energy (MPV) loss for Landau in a single layer
|
---|
188 | dP = chi*( TMath::Log(Wmax/I) + TMath::Log(chi/I) + 0.2 - beta*beta - delta);
|
---|
189 |
|
---|
190 | //cout<<"L: "<<L<<", PT: "<<candidateMomentum.Pt()<<", Eta: "<<candidateMomentum.Eta()<<", Phi: "<< candidateMomentum.Phi()<<endl;
|
---|
191 | //cout<<"Nhits: "<<nhits<<", dx: "<<dx<<", Charge: "<<charge<<", Beta: "<< beta<<", Gamma: "<<gamma<<", PT: "<<candidateMomentum.Pt()<<endl;
|
---|
192 | //cout<<x<<","<<kappa<<endl;
|
---|
193 | //cout<<" Wmax: "<<Wmax<<", Chi: "<<chi<<", delta: "<<delta<<", DeDx: "<<avdE<<", DeltaP: "<<dP<<endl;
|
---|
194 |
|
---|
195 | // simulate Nhits energy loss measurements
|
---|
196 | elosses.clear();
|
---|
197 | for (Int_t j=0; j<nhits; j++){
|
---|
198 | // compute total energy loss in MeV predicted by a Landau
|
---|
199 | dE = gRandom->Landau(dP,chi); // this is the total energy loss in MeV predicted by a Landau
|
---|
200 |
|
---|
201 | // convert resolution given in Mev/cm into absolute for this sensor
|
---|
202 | res = fResolution*dx;
|
---|
203 |
|
---|
204 | // apply additionnal gaussian smearing
|
---|
205 | dE = gRandom->Gaus(dE,res);
|
---|
206 | elosses.push_back(dE);
|
---|
207 | }
|
---|
208 |
|
---|
209 | sort (elosses.begin(), elosses.end());
|
---|
210 | eloss_truncmean = TruncatedMean(elosses, fTruncatedMeanFraction);
|
---|
211 |
|
---|
212 |
|
---|
213 | dEdx = dx > 0 ? eloss_truncmean/dx : -1. ;
|
---|
214 |
|
---|
215 | // store computed dEdx in MeV/cm
|
---|
216 | candidate->DeDx = dEdx;
|
---|
217 |
|
---|
218 | // add dedx also in Muons in electrons classes in treeWriter
|
---|
219 | // fix electrons here
|
---|
220 | // think whether any relevance for hits
|
---|
221 |
|
---|
222 |
|
---|
223 | //cout<<" eloss: "<<dE<<", dx: "<<dx<<", dEdx: "<<dEdx<<endl;
|
---|
224 | }
|
---|
225 | }
|
---|
226 |
|
---|
227 | }
|
---|
228 |
|
---|
229 | //------------------------------------------------------------------------------
|
---|
230 |
|
---|
231 | // formula Taken from Leo (2.30) pg. 26
|
---|
232 | Double_t EnergyLoss::Deltaf(Double_t c0, Double_t a, Double_t m, Double_t x0, Double_t x1, Double_t beta, Double_t gamma)
|
---|
233 | {
|
---|
234 | Double_t x= TMath::Log10(beta*gamma);
|
---|
235 | Double_t delta = 0.;
|
---|
236 |
|
---|
237 | //cout<<x<<","<<x0<<","<<x1<<","<<endl;
|
---|
238 |
|
---|
239 | if (x < x0)
|
---|
240 | delta = 0.;
|
---|
241 | if (x >= x0 && x< x1)
|
---|
242 | delta = 4.6052*x - c0 + a*TMath::Power(x1 - x,m);
|
---|
243 | if (x> x1)
|
---|
244 | delta = 4.6052*x - c0;
|
---|
245 |
|
---|
246 | return delta;
|
---|
247 | }
|
---|
248 |
|
---|
249 | //------------------------------------------------------------------------------
|
---|
250 | Double_t EnergyLoss::TruncatedMean(std::vector<Double_t> elosses, Double_t truncFrac)
|
---|
251 | {
|
---|
252 | Int_t new_size = Int_t( elosses.size() * (1 - truncFrac));
|
---|
253 |
|
---|
254 | // remove outliers and re-compute mean
|
---|
255 | elosses.resize(new_size);
|
---|
256 | return accumulate( elosses.begin(), elosses.end(), 0.0)/elosses.size();
|
---|
257 | }
|
---|