New class for ITS dE/dx samples
[u/mrichter/AliRoot.git] / 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 }