Fork me on GitHub

Changeset 6777565 in git


Ignore:
Timestamp:
Jan 26, 2020, 6:25:30 PM (5 years ago)
Author:
michele <michele.selvaggi@…>
Branches:
Timing
Children:
fec809d
Parents:
77249aa
Message:

truncated mean implementation of dedx

Files:
3 edited

Legend:

Unmodified
Added
Removed

  • cards/FCC/FCChh_PileUpVtx.tcl

    r77249aa r6777565  
    256256  add InputArray MuonMomentumSmearing/muons
    257257
    258   set ActiveFraction 0.013
    259   set ChargeCollectionEfficiency 0.75
    260 
    261   # fixme: this number should probably be charge/energy dependent, or absolute number in MeV/cm?
    262   set Resolution 0.15
    263 
    264   # active material properties (cf. http://pdg.lbl.gov/2014/AtomicNuclearProperties/properties8.dat)
     258
     259  # absolute resolution per measurement (normalized in MeV/cm)
     260  set Resolution 0.4
     261
     262  # fraction of measurements to ignore when computing truncated mean
     263  set TruncatedMeanFraction 0.5
     264
     265  # detector properties (active fraction = nhits*thickness/L)
     266  set Thickness 200E-6
     267  set ActiveFraction 0.00133
     268
     269  # Silicon properties, for other materials:
     270  # cf. http://pdg.lbl.gov/2014/AtomicNuclearProperties/properties8.dat
     271
    265272  set Z 14.
    266273  set A 28.0855

  • modules/EnergyLoss.cc

    r77249aa r6777565  
    1717 */
    1818
    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 
     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  */
    3127#include "modules/EnergyLoss.h"
    3228
    … …  
    5652//------------------------------------------------------------------------------
    5753
    58 EnergyLoss::EnergyLoss() 
     54EnergyLoss::EnergyLoss()
    5955{
    6056}
    … …  
    7167{
    7268
    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)
     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
    7873
    7974  // 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);
     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  fX0   =  GetDouble("c0", 4.4355);
    8984
    9085  // import arrays with output from other modules
    … …  
    125120void EnergyLoss::Process()
    126121{
    127   Candidate *candidate, *particle, *particleTest;
     122  Candidate *candidate, *particle;
    128123  vector<TIterator *>::iterator itInputList;
    129124  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;
     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;
    136134
    137135
    … …  
    145143    while((candidate = static_cast<Candidate *>(iterator->Next())))
    146144    {
    147       //cout<<"    ---------------- new candidate -------------------"<<endl;
     145      cout<<"    ---------------- new candidate -------------------"<<endl;
    148146      const TLorentzVector &candidateMomentum = candidate->Momentum;
    149147
    150       beta         = candidateMomentum.Beta();
    151       gamma        = candidateMomentum.Gamma();
    152 
     148      beta      = candidateMomentum.Beta();
     149      gamma     = candidateMomentum.Gamma();
    153150      charge    = TMath::Abs(candidate->Charge);
    154151
    155152      // 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
     153      //dx = candidate->L * fActiveFraction * 0.1;
     154      // amount of material in one sensor (converted in cm)
     155      dx = fThickness * 100.;
     156      // path length in cm
     157      L  = candidate->L * 0.1;
     158
     159
     160     // compute number of hits as path length over active length
     161      nhits = Int_t(L*fActiveFraction/dx);
     162
     163      //cout<<L<<","<<fActiveFraction<<","<<dx<<","<<nhits<<endl;
     164
     165      kappa = 2*0.1535*TMath::Abs(charge)*TMath::Abs(charge)*fZ*fRho*dx/(fA*beta*beta); //energy loss in MeV
    160166
    161167      chi = 0.5*kappa;
    162       me = 0.510998; // electron mass in MeV, need 
     168      me = 0.510998; // electron mass in MeV, need
    163169      I = fI*1e-6; // convert I in MeV
    164170
    … …  
    166172      Wmax = 2*me*beta*beta*gamma*gamma; // this is not valid for electrons
    167173
    168       delta = Deltaf(fc0, fa, fm, fx0, fx1, beta, gamma);
     174      delta = Deltaf(fC0, fAa, fM, fX0, fX1, beta, gamma);
    169175
    170176      // Bethe-Bloch energy loss in MeV (not used here)
    171177      avdE  = kappa*( TMath::Log(Wmax/I) - beta*beta - delta/2);
    172178
    173       // most probable energy (MPV) loss for Landau
     179      // most probable energy (MPV) loss for Landau in a single layer
    174180      dP = chi*( TMath::Log(Wmax/I) + TMath::Log(chi/I) + 0.2 - beta*beta - delta);
    175181
    176       //cout<<"x: "<<x<<", Beta: "<< beta<<", Gamma: "<<gamma <<", Charge: "<<charge<<endl;
    177 
     182      if (candidateMomentum.Pt() > 5) {
     183        cout<<"Nhits: "<<nhits<<", dx: "<<dx<<", Charge: "<<charge<<", Beta: "<< beta<<", Gamma: "<<gamma<<", PT: "<<candidateMomentum.Pt()<<endl;
    178184      //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. ;
     185        cout<<"    Wmax: "<<Wmax<<", Chi: "<<chi<<", delta: "<<delta<<", DeDx: "<<avdE<<", DeltaP: "<<dP<<endl;
     186      }
     187
     188      // simulate Nhits energy loss measurements
     189      elosses.clear();
     190      for (Int_t j=0; j<nhits; j++){
     191        // compute total energy loss in MeV predicted by a Landau
     192        dE = gRandom->Landau(dP,chi); // this is the total energy loss in MeV predicted by a Landau
     193
     194        // convert resolution given in Mev/cm into absolute for this sensor
     195        res = fResolution*dx;
     196
     197        // apply additionnal gaussian smearing
     198        dE = gRandom->Gaus(dE,res);
     199        elosses.push_back(dE);
     200      }
     201
     202      sort (elosses.begin(), elosses.end());
     203      eloss_truncmean = TruncatedMean(elosses, fTruncatedMeanFraction);
     204
     205
     206      dEdx = dx > 0 ? eloss_truncmean/dx : -1. ;
    188207
    189208      // store computed dEdx in MeV/cm
    190       candidate->DeDx = dEdx; 
    191  
     209      candidate->DeDx = dEdx;
     210
    192211      // add dedx also in Muons in electrons classes in treeWriter
    193212      // fix electrons here
    194213      // think whether any relevance for hits
    195214
    196  
     215
    197216      //cout<<"    eloss: "<<dE<<", dx: "<<dx<<", dEdx: "<<dEdx<<endl;
    198217    }
    … …  
    201220}
    202221
    203 
    204222//------------------------------------------------------------------------------
    205223
    … …  
    209227   Double_t x= TMath::Log10(beta*gamma);
    210228   Double_t delta = 0.;
    211    
     229
    212230   //cout<<x<<","<<x0<<","<<x1<<","<<endl;
    213    
     231
    214232   if (x < x0)
    215233       delta = 0.;
    … …  
    218236   if  (x> x1)
    219237       delta = 4.6052*x - c0;
    220        
     238
    221239   return delta;
    222240}
     241
     242//------------------------------------------------------------------------------
     243Double_t EnergyLoss::TruncatedMean(std::vector<Double_t> elosses, Double_t truncFrac)
     244{
     245     Int_t new_size = Int_t( elosses.size() * (1 - truncFrac));
     246
     247     // remove outliers and re-compute mean
     248     elosses.resize(new_size);
     249     return accumulate( elosses.begin(), elosses.end(), 0.0)/elosses.size();
     250}

  • modules/EnergyLoss.h

    r77249aa r6777565  
    2222/** \class EnergyLoss
    2323 *
    24  *  Subtract pile-up contribution from tracks.
     24 *  This module computes the charged energy loss according to the active material properties.
     25 *  The energy loss is simulated with a Landau convoluted by a Gaussian.
    2526 *
    26  *  \author P. Demin - UCL, Louvain-la-Neuve
     27 *  \author M. Selvaggi - CERN
    2728 *
    2829 */
    … …  
    4748
    4849private:
    49   Double_t fActiveFraction;
    50   Double_t fChargeCollectionEfficiency;
    51   Double_t fResolution; 
     50
     51  Double_t   fActiveFraction;
     52  Double_t   fThickness;
     53
     54  Double_t   fResolution;
     55  Double_t   fTruncatedMeanFraction;
    5256
    5357  // material parameters
    54   Double_t   fZ;   
    55   Double_t   fA;   
    56   Double_t   frho;
     58  Double_t   fZ;
     59  Double_t   fA;
     60  Double_t   fRho;
    5761
    58   Double_t   fa;   
    59   Double_t   fm;   
    60   Double_t   fx0; 
    61   Double_t   fx1; 
    62   Double_t   fI; 
    63   Double_t   fc0; 
     62  Double_t   fAa;
     63  Double_t   fM;
     64  Double_t   fX0;
     65  Double_t   fX1;
     66  Double_t   fI;
     67  Double_t   fC0;
    6468
    6569  // this function computes corrections due to polarisation of the material
    6670  Double_t Deltaf(Double_t c0, Double_t a, Double_t m, Double_t x0, Double_t x1, Double_t beta, Double_t gamma);
     71
     72  // this function computes truncated of dEdx measurements
     73  Double_t TruncatedMean(std::vector<Double_t> elosses, Double_t truncFrac);
    6774
    6875  std::vector<TIterator *> fInputList; //!
Note: See TracChangeset for help on using the changeset viewer.