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
  42   for(Int_t i=0; i<fgkMaxSamples; i++) fdEdxSamples[i]=0.;
  43 }
  44
  45 //______________________________________________________________________
  46 AliITSdEdxSamples::AliITSdEdxSamples(Int_t nSamples,Double_t* samples, Double_t mom, Int_t specie) :
  47   TObject(),
  48   fNSamples(nSamples),
  49   fP(mom),
  50   fParticleSpecie(specie)
  51 {
  52   // Standard constructor
  53
  54   SetSamples(nSamples,samples);
  55 }
  56 //______________________________________________________________________
  57 void AliITSdEdxSamples::SetSamples(Int_t nSamples, Double_t* samples){
  58   // Set the samples
  59
  60   if(nSamples>fgkMaxSamples){
  61     AliWarning(Form("Too many dE/dx samples,only first %d will be used",fgkMaxSamples));
  62     fNSamples=fgkMaxSamples;
  63   }else{
  64     fNSamples=nSamples;
  65   }
  66   for(Int_t i=0; i<fNSamples; i++) fdEdxSamples[i]=samples[i];
  67   for(Int_t i=fNSamples; i<fgkMaxSamples; i++) fdEdxSamples[i]=0.;
  68   return;
  69 }
  70 //______________________________________________________________________
  71 void AliITSdEdxSamples::SetSamplesAndMomenta(Int_t nSamples, Double_t* samples, Double_t* mom){
  72   // Set the samples
  73   SetSamples(nSamples,samples);
  74   for(Int_t i=0; i<fNSamples; i++) fPAtSample[i]=mom[i];
  75   for(Int_t i=fNSamples; i<fgkMaxSamples; i++) fPAtSample[i]=0.;
  76   return;
  77 }
  78 //______________________________________________________________________
  79 Double_t AliITSdEdxSamples::GetTruncatedMean(Double_t frac, Double_t mindedx) const {
  80   // compute truncated mean
  81
  82   Int_t nc=0;
  83   Double_t dedx[fgkMaxSamples];
  84   for (Int_t il=0; il<fNSamples; il++) { // count good (>0) dE/dx values
  85     if(fdEdxSamples[il]>mindedx){
  86       dedx[nc]= fdEdxSamples[il];
  87       nc++;
  88     }
  89   }
  90   if(nc<1) return 0.;
  91
  92   Int_t swap; // sort in ascending order
  93   do {
  94     swap=0;
  95     for (Int_t i=0; i<nc-1; i++) {
  96       if (dedx[i]<=dedx[i+1]) continue;
  97       Double_t tmp=dedx[i];
  98       dedx[i]=dedx[i+1];
  99       dedx[i+1]=tmp;
 100       swap++;
 101     }
 102   } while (swap);
 103
 104   Double_t sumamp=0,sumweight=0;
 105   Double_t weight[fgkMaxSamples];
 106   for(Int_t iw=0; iw<fgkMaxSamples; iw++) weight[iw]=0.;
 107   Int_t lastUsed=(Int_t)(frac*nc+0.00001);
 108   if(lastUsed==0) lastUsed=1;
 109   if(lastUsed>fgkMaxSamples) lastUsed=fgkMaxSamples;
 110   for(Int_t iw=0; iw<lastUsed; iw++) weight[iw]=1.;
 111   if((frac*nc-lastUsed)>0.4 && lastUsed<fgkMaxSamples) weight[lastUsed]=0.5;
 112   for (Int_t i=0; i<nc; i++) {
 113     AliDebug(5,Form("dE/dx %f   weight %f",dedx[i],weight[i]));
 114     sumamp+= dedx[i]*weight[i];
 115     sumweight+=weight[i];
 116   }
 117   if(sumweight>0.) return sumamp/sumweight;
 118   else return 0.;
 119
 120 }
 121 //______________________________________________________________________
 122 Double_t AliITSdEdxSamples::GetWeightedMean(Double_t mindedx) const {
 123   // compute generalized mean with k=-2 (used by CMS)
 124
 125   Int_t nc=0;
 126   Double_t dedx[fgkMaxSamples];
 127   for (Int_t il=0; il<fNSamples; il++) { // count good (>0) dE/dx values
 128     if(fdEdxSamples[il]>mindedx){
 129       dedx[nc]= fdEdxSamples[il];
 130       nc++;
 131     }
 132   }
 133   if(nc<1) return 0.;
 134
 135   Double_t weiSum = 0.;
 136   for (Int_t i=0; i<nc; i++) {
 137     weiSum+=TMath::Power(dedx[i],-2);
 138   }
 139   Double_t wMean=0.;
 140   if(weiSum>0.) wMean= TMath::Power(weiSum/nc,-0.5);
 141   return wMean;
 142
 143 }
 144 //______________________________________________________________________
 145 void  AliITSdEdxSamples::GetConditionalProbabilities(AliITSPidParams* pars, Double_t condprob[AliPID::kSPECIES], Double_t mindedx) const {
 146   // compute conditional probablilities
 147   const Int_t nPart = 3;
 148   Double_t itsProb[nPart] = {1,1,1}; // p, K, pi
 149
 150   for(Int_t iS=0; iS<fNSamples; iS++){
 151     Int_t iLayer=iS+3; // to match with present ITS
 152     if(iLayer>6) iLayer=6; // all extra points are treated as SSD
 153     if(fdEdxSamples[iS]<mindedx) continue;
 154     Float_t dedx = fdEdxSamples[iS];
 155     Float_t layProb = pars->GetLandauGausNorm(dedx,AliPID::kProton,fP,iLayer);
 156     itsProb[0] *= layProb;
 157
 158     layProb = pars->GetLandauGausNorm(dedx,AliPID::kKaon,fP,iLayer);
 159     if (fP < 0.16) layProb=0.00001;
 160     itsProb[1] *= layProb;
 161
 162     layProb = pars->GetLandauGausNorm(dedx,AliPID::kPion,fP,iLayer);
 163     itsProb[2] *= layProb;
 164   }
 165
 166   // Normalise probabilities
 167   Double_t sumProb = 0;
 168   for (Int_t iPart = 0; iPart < nPart; iPart++) {
 169     sumProb += itsProb[iPart];
 170   }
 171   sumProb += 2*itsProb[2]; // muon and electron cannot be distinguished from pions
 172
 173   for (Int_t iPart = 0; iPart < nPart; iPart++) {
 174     itsProb[iPart]/=sumProb;
 175   }
 176
 177   condprob[AliPID::kElectron] = itsProb[2];
 178   condprob[AliPID::kMuon] = itsProb[2];
 179   condprob[AliPID::kPion] = itsProb[2];
 180   condprob[AliPID::kKaon] = itsProb[1];
 181   condprob[AliPID::kProton] = itsProb[0];
 182   return;
 183
 184 }