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