[b982244] | 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 |
|
---|
[6777565] | 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 | */
|
---|
[b982244] | 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 |
|
---|
[6777565] | 54 | EnergyLoss::EnergyLoss()
|
---|
[b982244] | 55 | {
|
---|
| 56 | }
|
---|
| 57 |
|
---|
| 58 | //------------------------------------------------------------------------------
|
---|
| 59 |
|
---|
| 60 | EnergyLoss::~EnergyLoss()
|
---|
| 61 | {
|
---|
| 62 | }
|
---|
| 63 |
|
---|
| 64 | //------------------------------------------------------------------------------
|
---|
| 65 |
|
---|
| 66 | void EnergyLoss::Init()
|
---|
| 67 | {
|
---|
| 68 |
|
---|
[6777565] | 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
|
---|
[b982244] | 73 |
|
---|
| 74 | // active material properties (cf. http://pdg.lbl.gov/2014/AtomicNuclearProperties/properties8.dat)
|
---|
[6777565] | 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)
|
---|
[fec809d] | 83 | fC0 = GetDouble("c0", 4.4355);
|
---|
[b982244] | 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 | {
|
---|
[fec809d] | 122 | Candidate *candidate;
|
---|
[b982244] | 123 | vector<TIterator *>::iterator itInputList;
|
---|
| 124 | TIterator *iterator;
|
---|
| 125 |
|
---|
[6777565] | 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;
|
---|
[b982244] | 132 |
|
---|
[fec809d] | 133 | //cout<<"---------------- new event -------------------"<<endl;
|
---|
[b982244] | 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 | {
|
---|
[fec809d] | 145 | //cout<<" ---------------- new candidate -------------------"<<endl;
|
---|
[b982244] | 146 | const TLorentzVector &candidateMomentum = candidate->Momentum;
|
---|
| 147 |
|
---|
[6777565] | 148 | beta = candidateMomentum.Beta();
|
---|
| 149 | gamma = candidateMomentum.Gamma();
|
---|
[b982244] | 150 | charge = TMath::Abs(candidate->Charge);
|
---|
| 151 |
|
---|
[fec809d] | 152 |
|
---|
[b982244] | 153 | // length of the track normalized by the fraction of active material and the charge collection efficiency in the tracker (in cm)
|
---|
[6777565] | 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 |
|
---|
[fec809d] | 164 |
|
---|
| 165 | //beta = 0.999945;
|
---|
| 166 | //gamma = 95.6446;
|
---|
| 167 | //charge = 1.;
|
---|
| 168 | //nhits = 100;
|
---|
| 169 |
|
---|
[6777565] | 170 | //cout<<L<<","<<fActiveFraction<<","<<dx<<","<<nhits<<endl;
|
---|
[b982244] | 171 |
|
---|
[6777565] | 172 | kappa = 2*0.1535*TMath::Abs(charge)*TMath::Abs(charge)*fZ*fRho*dx/(fA*beta*beta); //energy loss in MeV
|
---|
[b982244] | 173 |
|
---|
| 174 | chi = 0.5*kappa;
|
---|
[6777565] | 175 | me = 0.510998; // electron mass in MeV, need
|
---|
[b982244] | 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 |
|
---|
[6777565] | 181 | delta = Deltaf(fC0, fAa, fM, fX0, fX1, beta, gamma);
|
---|
[b982244] | 182 |
|
---|
| 183 | // Bethe-Bloch energy loss in MeV (not used here)
|
---|
| 184 | avdE = kappa*( TMath::Log(Wmax/I) - beta*beta - delta/2);
|
---|
| 185 |
|
---|
[6777565] | 186 | // most probable energy (MPV) loss for Landau in a single layer
|
---|
[b982244] | 187 | dP = chi*( TMath::Log(Wmax/I) + TMath::Log(chi/I) + 0.2 - beta*beta - delta);
|
---|
| 188 |
|
---|
[fec809d] | 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;
|
---|
[b982244] | 191 | //cout<<x<<","<<kappa<<endl;
|
---|
[fec809d] | 192 | //cout<<" Wmax: "<<Wmax<<", Chi: "<<chi<<", delta: "<<delta<<", DeDx: "<<avdE<<", DeltaP: "<<dP<<endl;
|
---|
[6777565] | 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;
|
---|
[b982244] | 202 |
|
---|
[6777565] | 203 | // apply additionnal gaussian smearing
|
---|
| 204 | dE = gRandom->Gaus(dE,res);
|
---|
| 205 | elosses.push_back(dE);
|
---|
| 206 | }
|
---|
[b982244] | 207 |
|
---|
[6777565] | 208 | sort (elosses.begin(), elosses.end());
|
---|
| 209 | eloss_truncmean = TruncatedMean(elosses, fTruncatedMeanFraction);
|
---|
| 210 |
|
---|
| 211 |
|
---|
| 212 | dEdx = dx > 0 ? eloss_truncmean/dx : -1. ;
|
---|
[b982244] | 213 |
|
---|
| 214 | // store computed dEdx in MeV/cm
|
---|
[6777565] | 215 | candidate->DeDx = dEdx;
|
---|
| 216 |
|
---|
[b982244] | 217 | // add dedx also in Muons in electrons classes in treeWriter
|
---|
| 218 | // fix electrons here
|
---|
| 219 | // think whether any relevance for hits
|
---|
| 220 |
|
---|
[6777565] | 221 |
|
---|
[b982244] | 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.;
|
---|
[6777565] | 235 |
|
---|
[b982244] | 236 | //cout<<x<<","<<x0<<","<<x1<<","<<endl;
|
---|
[6777565] | 237 |
|
---|
[b982244] | 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;
|
---|
[6777565] | 244 |
|
---|
[b982244] | 245 | return delta;
|
---|
| 246 | }
|
---|
[6777565] | 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 | }
|
---|