ITS/AliITSdEdxSamples.cxx

   1 /**************************************************************************
   2  * Copyright(c) 2009-2012, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15 /* $Id$ */
  16
  17
  18
  19 ///////////////////////////////////////////////////////////////////
  20 //                                                               //
  21 // Class to store information for PID with ITS                   //
  22 // and truncated mean computation methods                        //
  23 // Origin: F.Prino, Torino, prino@to.infn.it                     //
  24 //                                                               //
  25 ///////////////////////////////////////////////////////////////////
  26
  27 #include "AliITSPidParams.h"
  28 #include "AliITSdEdxSamples.h"
  29 #include "AliLog.h"
  30 #include <TMath.h>
  31
  32 ClassImp(AliITSdEdxSamples)
  33
  34 //______________________________________________________________________
  35 AliITSdEdxSamples::AliITSdEdxSamples():TObject(),
  36   fNSamples(0),
  37   fP(0.),
  38   fParticleSpecie(0)
  39 {
  40   // Default constructor
  41   for(Int_t i=0; i<kMaxSamples; i++){
  42     fdESamples[i]=0.;
  43     fdxSamples[i]=0.;
  44     fPAtSample[i]=0.;
  45   }
  46 }
  47
  48 //______________________________________________________________________
  49 AliITSdEdxSamples::AliITSdEdxSamples(Int_t nSamples, Double_t* esamples, Double_t* xsamples, Double_t mom, Int_t specie) :
  50   TObject(),
  51   fNSamples(nSamples),
  52   fP(mom),
  53   fParticleSpecie(specie)
  54 {
  55   // Standard constructor
  56   SetdESamples(nSamples,esamples);
  57   SetdxSamples(nSamples,xsamples);
  58 }
  59
  60 //______________________________________________________________________
  61 AliITSdEdxSamples::AliITSdEdxSamples(const AliITSdEdxSamples& source) :
  62   TObject(),
  63   fNSamples(source.fNSamples),
  64   fP(source.fP),
  65   fParticleSpecie(source.fParticleSpecie)
  66 {
  67   // Copy constructor
  68   for(Int_t i=0; i<kMaxSamples; i++){
  69     fdESamples[i]=source.GetdESample(i);
  70     fdxSamples[i]=source.GetdxSample(i);
  71     fPAtSample[i]=source.GetMomentumAtSample(i);
  72   }
  73 }
  74 //______________________________________________________________________
  75 void AliITSdEdxSamples::SetdESamples(Int_t nSamples, Double_t* samples){
  76   // Set the samples
  77
  78   if(nSamples>kMaxSamples){
  79     AliWarning(Form("Too many dE samples,only first %d will be used",kMaxSamples));
  80     fNSamples=kMaxSamples;
  81   }else{
  82     fNSamples=nSamples;
  83   }
  84   for(Int_t i=0; i<fNSamples; i++) fdESamples[i]=samples[i];
  85   for(Int_t i=fNSamples; i<kMaxSamples; i++) fdESamples[i]=0.;
  86   return;
  87 }
  88 //______________________________________________________________________
  89 void AliITSdEdxSamples::SetdxSamples(Int_t nSamples, Double_t* samples){
  90   // Set the samples
  91
  92   if(nSamples>kMaxSamples){
  93     AliWarning(Form("Too many dx samples,only first %d will be used",kMaxSamples))
  94     fNSamples=kMaxSamples;
  95   }else{
  96     fNSamples=nSamples;
  97   }
  98   for(Int_t i=0; i<fNSamples; i++) fdxSamples[i]=samples[i];
  99   for(Int_t i=fNSamples; i<kMaxSamples; i++) fdxSamples[i]=0.;
 100   return;
 101 }
 102
 103 //______________________________________________________________________
 104 void AliITSdEdxSamples::SetSamplesAndMomenta(Int_t nSamples, Double_t* esamples, Double_t* xsamples, Double_t* mom){
 105   // Set the samples
 106   SetdESamples(nSamples,esamples);
 107   SetdxSamples(nSamples,xsamples);
 108   for(Int_t i=0; i<fNSamples; i++) fPAtSample[i]=mom[i];
 109   for(Int_t i=fNSamples; i<kMaxSamples; i++) fPAtSample[i]=0.;
 110   return;
 111 }
 112 //______________________________________________________________________
 113 Double_t AliITSdEdxSamples::GetTruncatedMean(Double_t frac, Double_t mindedx) const {
 114   // compute truncated mean
 115
 116   Int_t nc=0;
 117   Double_t dedx[kMaxSamples];
 118   for (Int_t il=0; il<fNSamples; il++) { // count good (>0) dE/dx values
 119     Double_t dedxsamp=GetdEdxSample(il);
 120     if(dedxsamp>mindedx){
 121       dedx[nc]= dedxsamp;
 122       nc++;
 123     }
 124   }
 125   if(nc<1) return 0.;
 126
 127   Int_t swap; // sort in ascending order
 128   do {
 129     swap=0;
 130     for (Int_t i=0; i<nc-1; i++) {
 131       if (dedx[i]<=dedx[i+1]) continue;
 132       Double_t tmp=dedx[i];
 133       dedx[i]=dedx[i+1];
 134       dedx[i+1]=tmp;
 135       swap++;
 136     }
 137   } while (swap);
 138
 139   Double_t sumamp=0,sumweight=0;
 140   Double_t weight[kMaxSamples];
 141   for(Int_t iw=0; iw<kMaxSamples; iw++) weight[iw]=0.;
 142   Int_t lastUsed=(Int_t)(frac*nc+0.00001);
 143   if(lastUsed==0) lastUsed=1;
 144   if(lastUsed>kMaxSamples) lastUsed=kMaxSamples;
 145   for(Int_t iw=0; iw<lastUsed; iw++) weight[iw]=1.;
 146   if((frac*nc-lastUsed)>0.4 && lastUsed<kMaxSamples) weight[lastUsed]=0.5;
 147   for (Int_t i=0; i<nc; i++) {
 148     // AliDebug(5,Form("dE/dx %f   weight %f",dedx[i],weight[i]));
 149     sumamp+= dedx[i]*weight[i];
 150     sumweight+=weight[i];
 151   }
 152   if(sumweight>0.) return sumamp/sumweight;
 153   else return 0.;
 154
 155 }
 156 //______________________________________________________________________
 157 Double_t AliITSdEdxSamples::GetWeightedMean(Double_t mindedx) const {
 158   // compute generalized mean with k=-2 (used by CMS)
 159   Int_t nc=0;
 160   Double_t dedx[kMaxSamples];
 161   for (Int_t il=0; il<fNSamples; il++) { // count good (>0) dE/dx values
 162     Double_t dedxsamp=GetdEdxSample(il);
 163     if(dedxsamp>mindedx){
 164       dedx[nc]= dedxsamp;
 165       nc++;
 166     }
 167   }
 168   if(nc<1) return 0.;
 169
 170   Double_t weiSum = 0.;
 171   for (Int_t i=0; i<nc; i++) {
 172     weiSum+=TMath::Power(dedx[i],-2);
 173   }
 174   Double_t wMean=0.;
 175   if(weiSum>0.) wMean= TMath::Power(weiSum/nc,-0.5);
 176   return wMean;
 177
 178 }
 179 //______________________________________________________________________
 180 void  AliITSdEdxSamples::GetConditionalProbabilities(AliITSPidParams* pars, Double_t condprob[AliPID::kSPECIES], Double_t mindedx) const {
 181   // compute conditional probablilities
 182   const Int_t nPart = 3;
 183   Double_t itsProb[nPart] = {1,1,1}; // p, K, pi
 184
 185   for(Int_t iS=0; iS<fNSamples; iS++){
 186     Int_t iLayer=iS+3; // to match with present ITS
 187     if(iLayer>6) iLayer=6; // all extra points are treated as SSD
 188     Float_t dedx = GetdEdxSample(iS);
 189     if(dedx<mindedx) continue;
 190     Float_t layProb = pars->GetLandauGausNorm(dedx,AliPID::kProton,fP,iLayer);
 191     itsProb[0] *= layProb;
 192
 193     layProb = pars->GetLandauGausNorm(dedx,AliPID::kKaon,fP,iLayer);
 194     if (fP < 0.16) layProb=0.00001;
 195     itsProb[1] *= layProb;
 196
 197     layProb = pars->GetLandauGausNorm(dedx,AliPID::kPion,fP,iLayer);
 198     itsProb[2] *= layProb;
 199   }
 200
 201   // Normalise probabilities
 202   Double_t sumProb = 0;
 203   for (Int_t iPart = 0; iPart < nPart; iPart++) {
 204     sumProb += itsProb[iPart];
 205   }
 206   sumProb += 2*itsProb[2]; // muon and electron cannot be distinguished from pions
 207
 208   for (Int_t iPart = 0; iPart < nPart; iPart++) {
 209     itsProb[iPart]/=sumProb;
 210   }
 211
 212   condprob[AliPID::kElectron] = itsProb[2];
 213   condprob[AliPID::kMuon] = itsProb[2];
 214   condprob[AliPID::kPion] = itsProb[2];
 215   condprob[AliPID::kKaon] = itsProb[1];
 216   condprob[AliPID::kProton] = itsProb[0];
 217   return;
 218
 219 }