Updated class to store 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   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<fNSamples; 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 }