Fork me on GitHub

source: git/modules/EnergyLoss.cc@ e70228d

Timing
Last change on this file since e70228d was b982244, checked in by Michele Selvaggi <michele.selvaggi@…>, 5 years ago

adding energy loss module

  • Property mode set to 100644
File size: 7.0 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/** \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. The active material
23 * is assumed to be uniformly distributed in the detector volume. The actual active volume
24 * get normalized by multiplying the path length by the parameter ActiveFraction.
25 *
26 *
27 * \author M. Selvaggi - CERN
28 *
29 */
30
31#include "modules/EnergyLoss.h"
32
33#include "classes/DelphesClasses.h"
34#include "classes/DelphesFactory.h"
35#include "classes/DelphesFormula.h"
36
37#include "ExRootAnalysis/ExRootClassifier.h"
38#include "ExRootAnalysis/ExRootFilter.h"
39#include "ExRootAnalysis/ExRootResult.h"
40
41#include "TDatabasePDG.h"
42#include "TFormula.h"
43#include "TLorentzVector.h"
44#include "TMath.h"
45#include "TObjArray.h"
46#include "TRandom3.h"
47#include "TString.h"
48
49#include <algorithm>
50#include <iostream>
51#include <sstream>
52#include <stdexcept>
53
54using namespace std;
55
56//------------------------------------------------------------------------------
57
58EnergyLoss::EnergyLoss()
59{
60}
61
62//------------------------------------------------------------------------------
63
64EnergyLoss::~EnergyLoss()
65{
66}
67
68//------------------------------------------------------------------------------
69
70void EnergyLoss::Init()
71{
72
73 fActiveFraction = GetDouble("ActiveFraction", 0.013); // fraction of active material that measures the deposited charge
74 fChargeCollectionEfficiency = GetDouble("ChargeCollectionEfficiency", 0.75); // this number shifts Landau to the left
75
76 // fixme: this number should probably be charge/energy dependent
77 fResolution = GetDouble("Resolution", 0.15); // 0 - perfect Landau energy loss (0.15 gives good agreement with CMS pixel detector)
78
79 // active material properties (cf. http://pdg.lbl.gov/2014/AtomicNuclearProperties/properties8.dat)
80 fZ = GetDouble("Z", 14.);
81 fA = GetDouble("A", 28.0855); // in g/mol
82 frho = GetDouble("rho", 2.329); // in g/cm3
83 fa = GetDouble("a", 0.1492);
84 fm = GetDouble("m", 3.2546);
85 fx0 = GetDouble("x0", 0.2015);
86 fx1 = GetDouble("x1", 2.8716);
87 fI = GetDouble("I", 173.0); // mean excitation potential in (eV)
88 fc0 = GetDouble("c0", 4.4355);
89
90 // import arrays with output from other modules
91
92 ExRootConfParam param = GetParam("InputArray");
93 Long_t i, size;
94 const TObjArray *array;
95 TIterator *iterator;
96
97 size = param.GetSize();
98 for(i = 0; i < size; ++i)
99 {
100 array = ImportArray(param[i].GetString());
101 iterator = array->MakeIterator();
102
103 fInputList.push_back(iterator);
104 }
105
106}
107
108//------------------------------------------------------------------------------
109
110void EnergyLoss::Finish()
111{
112 vector<TIterator *>::iterator itInputList;
113 TIterator *iterator;
114
115 for(itInputList = fInputList.begin(); itInputList != fInputList.end(); ++itInputList)
116 {
117 iterator = *itInputList;
118 if(iterator) delete iterator;
119 }
120
121}
122
123//------------------------------------------------------------------------------
124
125void EnergyLoss::Process()
126{
127 Candidate *candidate, *particle, *particleTest;
128 vector<TIterator *>::iterator itInputList;
129 TIterator *iterator;
130 TObjArray *array;
131
132 Double_t beta, gamma, charge, x;
133 Double_t kappa, chi, me, I, Wmax, delta, avdE, dP, dx, dE, dEdx;
134
135 //cout<<"---------------- new event -------------------"<<endl;
136
137
138 // loop over all input arrays
139 for(itInputList = fInputList.begin(); itInputList != fInputList.end(); ++itInputList)
140 {
141 iterator = *itInputList;
142
143 // loop over all candidates
144 iterator->Reset();
145 while((candidate = static_cast<Candidate *>(iterator->Next())))
146 {
147 //cout<<" ---------------- new candidate -------------------"<<endl;
148 const TLorentzVector &candidateMomentum = candidate->Momentum;
149
150 beta = candidateMomentum.Beta();
151 gamma = candidateMomentum.Gamma();
152
153 charge = TMath::Abs(candidate->Charge);
154
155 // length of the track normalized by the fraction of active material and the charge collection efficiency in the tracker (in cm)
156 dx = candidate->L * fActiveFraction * 0.1;
157 x = dx * fChargeCollectionEfficiency;
158
159 kappa = 2*0.1535*TMath::Abs(charge)*TMath::Abs(charge)*fZ*frho*x/(fA*beta*beta); //energy loss in MeV
160
161 chi = 0.5*kappa;
162 me = 0.510998; // electron mass in MeV, need
163 I = fI*1e-6; // convert I in MeV
164
165 // fixme: max energy transfer wrong for electrons
166 Wmax = 2*me*beta*beta*gamma*gamma; // this is not valid for electrons
167
168 delta = Deltaf(fc0, fa, fm, fx0, fx1, beta, gamma);
169
170 // Bethe-Bloch energy loss in MeV (not used here)
171 avdE = kappa*( TMath::Log(Wmax/I) - beta*beta - delta/2);
172
173 // most probable energy (MPV) loss for Landau
174 dP = chi*( TMath::Log(Wmax/I) + TMath::Log(chi/I) + 0.2 - beta*beta - delta);
175
176 //cout<<"x: "<<x<<", Beta: "<< beta<<", Gamma: "<<gamma <<", Charge: "<<charge<<endl;
177
178 //cout<<x<<","<<kappa<<endl;
179 //cout<<" Wmax: "<<Wmax<<", Chi: "<<chi<<", delta: "<<delta<<", DeDx: "<<avdE<<", DeltaP: "<<dP<<endl;
180
181 // compute total energy loss in MeV predicted by a Landau
182 dE = gRandom->Landau(dP,chi); // this is the total energy loss in MeV predicted by a Landau
183
184 // apply additionnal gaussian smearing to simulate finite resolution in charge measurement
185 dE = gRandom->Gaus(dE,fResolution*dP);
186
187 dEdx = dx > 0 ? dE/dx : -1. ;
188
189 // store computed dEdx in MeV/cm
190 candidate->DeDx = dEdx;
191
192 // add dedx also in Muons in electrons classes in treeWriter
193 // fix electrons here
194 // think whether any relevance for hits
195
196
197 //cout<<" eloss: "<<dE<<", dx: "<<dx<<", dEdx: "<<dEdx<<endl;
198 }
199 }
200
201}
202
203
204//------------------------------------------------------------------------------
205
206// formula Taken from Leo (2.30) pg. 26
207Double_t EnergyLoss::Deltaf(Double_t c0, Double_t a, Double_t m, Double_t x0, Double_t x1, Double_t beta, Double_t gamma)
208{
209 Double_t x= TMath::Log10(beta*gamma);
210 Double_t delta = 0.;
211
212 //cout<<x<<","<<x0<<","<<x1<<","<<endl;
213
214 if (x < x0)
215 delta = 0.;
216 if (x >= x0 && x< x1)
217 delta = 4.6052*x - c0 + a*TMath::Power(x1 - x,m);
218 if (x> x1)
219 delta = 4.6052*x - c0;
220
221 return delta;
222}
Note: See TracBrowser for help on using the repository browser.