AddTaskFemto for train update
[u/mrichter/AliRoot.git] / ITS / AliITSdEdxSamples.cxx
CommitLineData
9189f38b 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
32ClassImp(AliITSdEdxSamples)
33
34//______________________________________________________________________
35AliITSdEdxSamples::AliITSdEdxSamples():TObject(),
36 fNSamples(0),
0d66557a 37 fClusterMap(0),
9189f38b 38 fP(0.),
0d66557a 39 fParticleSpecie(0),
40 fLayersForPid(0xFFFF)
9189f38b 41{
42 // Default constructor
92329a66 43 for(Int_t i=0; i<kMaxSamples; i++){
44 fdESamples[i]=0.;
45 fdxSamples[i]=0.;
cb0a52e0 46 fPAtSample[i]=0.;
47 }
9189f38b 48}
49
50//______________________________________________________________________
92329a66 51AliITSdEdxSamples::AliITSdEdxSamples(Int_t nSamples, Double_t* esamples, Double_t* xsamples, Double_t mom, Int_t specie) :
9189f38b 52 TObject(),
53 fNSamples(nSamples),
0d66557a 54 fClusterMap(0),
9189f38b 55 fP(mom),
0d66557a 56 fParticleSpecie(specie),
57 fLayersForPid(0xFFFF)
9189f38b 58{
59 // Standard constructor
92329a66 60 SetdESamples(nSamples,esamples);
61 SetdxSamples(nSamples,xsamples);
0d66557a 62 SetClusterMapFromdE();
92329a66 63}
64
65//______________________________________________________________________
66AliITSdEdxSamples::AliITSdEdxSamples(const AliITSdEdxSamples& source) :
67 TObject(),
68 fNSamples(source.fNSamples),
0d66557a 69 fClusterMap(source.fClusterMap),
92329a66 70 fP(source.fP),
0d66557a 71 fParticleSpecie(source.fParticleSpecie),
72 fLayersForPid(source.fLayersForPid)
92329a66 73{
74 // Copy constructor
8d91dc44 75 for(Int_t i=0; i<kMaxSamples; i++){
92329a66 76 fdESamples[i]=source.GetdESample(i);
77 fdxSamples[i]=source.GetdxSample(i);
78 fPAtSample[i]=source.GetMomentumAtSample(i);
79 }
9189f38b 80}
81//______________________________________________________________________
92329a66 82void AliITSdEdxSamples::SetdESamples(Int_t nSamples, Double_t* samples){
9189f38b 83 // Set the samples
84
92329a66 85 if(nSamples>kMaxSamples){
86 AliWarning(Form("Too many dE samples,only first %d will be used",kMaxSamples));
87 fNSamples=kMaxSamples;
9189f38b 88 }else{
89 fNSamples=nSamples;
90 }
92329a66 91 for(Int_t i=0; i<fNSamples; i++) fdESamples[i]=samples[i];
92 for(Int_t i=fNSamples; i<kMaxSamples; i++) fdESamples[i]=0.;
9189f38b 93 return;
94}
95//______________________________________________________________________
92329a66 96void AliITSdEdxSamples::SetdxSamples(Int_t nSamples, Double_t* samples){
9189f38b 97 // Set the samples
92329a66 98
99 if(nSamples>kMaxSamples){
13242232 100 AliWarning(Form("Too many dx samples,only first %d will be used",kMaxSamples));
92329a66 101 fNSamples=kMaxSamples;
102 }else{
103 fNSamples=nSamples;
104 }
105 for(Int_t i=0; i<fNSamples; i++) fdxSamples[i]=samples[i];
106 for(Int_t i=fNSamples; i<kMaxSamples; i++) fdxSamples[i]=0.;
107 return;
108}
109
110//______________________________________________________________________
111void AliITSdEdxSamples::SetSamplesAndMomenta(Int_t nSamples, Double_t* esamples, Double_t* xsamples, Double_t* mom){
112 // Set the samples
113 SetdESamples(nSamples,esamples);
114 SetdxSamples(nSamples,xsamples);
9189f38b 115 for(Int_t i=0; i<fNSamples; i++) fPAtSample[i]=mom[i];
92329a66 116 for(Int_t i=fNSamples; i<kMaxSamples; i++) fPAtSample[i]=0.;
9189f38b 117 return;
118}
119//______________________________________________________________________
0d66557a 120void AliITSdEdxSamples::SetLayerSample(Int_t iLayer, Bool_t haspoint, Double_t dE, Double_t dx, Double_t p){
121 // set info from single layer
122 if(haspoint){
123 SetPointOnLayer(iLayer);
124 fdESamples[iLayer]=dE;
125 fdxSamples[iLayer]=dx;
126 fPAtSample[iLayer]=p;
127 }else{
128 if(HasPointOnLayer(iLayer)) fClusterMap-=(1<<iLayer);
129 fdESamples[iLayer]=0.;
130 fdxSamples[iLayer]=0.;
131 fPAtSample[iLayer]=0.;
132
133 }
134}
135//______________________________________________________________________
9189f38b 136Double_t AliITSdEdxSamples::GetTruncatedMean(Double_t frac, Double_t mindedx) const {
137 // compute truncated mean
138
139 Int_t nc=0;
92329a66 140 Double_t dedx[kMaxSamples];
9189f38b 141 for (Int_t il=0; il<fNSamples; il++) { // count good (>0) dE/dx values
92329a66 142 Double_t dedxsamp=GetdEdxSample(il);
0d66557a 143 if(HasPointOnLayer(il) && UseLayerForPid(il) && dedxsamp>mindedx){
92329a66 144 dedx[nc]= dedxsamp;
9189f38b 145 nc++;
92329a66 146 }
9189f38b 147 }
148 if(nc<1) return 0.;
0d66557a 149
9189f38b 150 Int_t swap; // sort in ascending order
151 do {
152 swap=0;
153 for (Int_t i=0; i<nc-1; i++) {
154 if (dedx[i]<=dedx[i+1]) continue;
155 Double_t tmp=dedx[i];
156 dedx[i]=dedx[i+1];
157 dedx[i+1]=tmp;
158 swap++;
159 }
160 } while (swap);
161
162 Double_t sumamp=0,sumweight=0;
92329a66 163 Double_t weight[kMaxSamples];
164 for(Int_t iw=0; iw<kMaxSamples; iw++) weight[iw]=0.;
9189f38b 165 Int_t lastUsed=(Int_t)(frac*nc+0.00001);
166 if(lastUsed==0) lastUsed=1;
92329a66 167 if(lastUsed>kMaxSamples) lastUsed=kMaxSamples;
9189f38b 168 for(Int_t iw=0; iw<lastUsed; iw++) weight[iw]=1.;
92329a66 169 if((frac*nc-lastUsed)>0.4 && lastUsed<kMaxSamples) weight[lastUsed]=0.5;
9189f38b 170 for (Int_t i=0; i<nc; i++) {
92329a66 171 // AliDebug(5,Form("dE/dx %f weight %f",dedx[i],weight[i]));
9189f38b 172 sumamp+= dedx[i]*weight[i];
173 sumweight+=weight[i];
174 }
175 if(sumweight>0.) return sumamp/sumweight;
176 else return 0.;
177
178}
179//______________________________________________________________________
180Double_t AliITSdEdxSamples::GetWeightedMean(Double_t mindedx) const {
181 // compute generalized mean with k=-2 (used by CMS)
9189f38b 182 Int_t nc=0;
92329a66 183 Double_t dedx[kMaxSamples];
9189f38b 184 for (Int_t il=0; il<fNSamples; il++) { // count good (>0) dE/dx values
92329a66 185 Double_t dedxsamp=GetdEdxSample(il);
0d66557a 186 if(HasPointOnLayer(il) && UseLayerForPid(il) && dedxsamp>mindedx){
92329a66 187 dedx[nc]= dedxsamp;
188 nc++;
9189f38b 189 }
190 }
191 if(nc<1) return 0.;
192
193 Double_t weiSum = 0.;
194 for (Int_t i=0; i<nc; i++) {
195 weiSum+=TMath::Power(dedx[i],-2);
196 }
197 Double_t wMean=0.;
198 if(weiSum>0.) wMean= TMath::Power(weiSum/nc,-0.5);
199 return wMean;
200
201}
202//______________________________________________________________________
203void AliITSdEdxSamples::GetConditionalProbabilities(AliITSPidParams* pars, Double_t condprob[AliPID::kSPECIES], Double_t mindedx) const {
204 // compute conditional probablilities
205 const Int_t nPart = 3;
206 Double_t itsProb[nPart] = {1,1,1}; // p, K, pi
207
208 for(Int_t iS=0; iS<fNSamples; iS++){
0d66557a 209 if(!HasPointOnLayer(iS)) continue;
210 if(!UseLayerForPid(iS)) continue;
9189f38b 211 Int_t iLayer=iS+3; // to match with present ITS
212 if(iLayer>6) iLayer=6; // all extra points are treated as SSD
92329a66 213 Float_t dedx = GetdEdxSample(iS);
214 if(dedx<mindedx) continue;
9189f38b 215 Float_t layProb = pars->GetLandauGausNorm(dedx,AliPID::kProton,fP,iLayer);
216 itsProb[0] *= layProb;
217
218 layProb = pars->GetLandauGausNorm(dedx,AliPID::kKaon,fP,iLayer);
219 if (fP < 0.16) layProb=0.00001;
220 itsProb[1] *= layProb;
221
222 layProb = pars->GetLandauGausNorm(dedx,AliPID::kPion,fP,iLayer);
223 itsProb[2] *= layProb;
224 }
225
226 // Normalise probabilities
227 Double_t sumProb = 0;
228 for (Int_t iPart = 0; iPart < nPart; iPart++) {
229 sumProb += itsProb[iPart];
230 }
231 sumProb += 2*itsProb[2]; // muon and electron cannot be distinguished from pions
232
233 for (Int_t iPart = 0; iPart < nPart; iPart++) {
234 itsProb[iPart]/=sumProb;
235 }
236
237 condprob[AliPID::kElectron] = itsProb[2];
238 condprob[AliPID::kMuon] = itsProb[2];
239 condprob[AliPID::kPion] = itsProb[2];
240 condprob[AliPID::kKaon] = itsProb[1];
241 condprob[AliPID::kProton] = itsProb[0];
242 return;
243
244}
0d66557a 245//______________________________________________________________________
246void AliITSdEdxSamples::PrintAll() const{
247 // print all the infos
248 printf("Particle %d momentum %f GeV/c, number of points %d\n",
249 GetParticleSpecieMC(),
250 fP,
251 GetNumberOfEffectiveSamples());
252 for(Int_t iLay=0; iLay<fNSamples; iLay++){
253 printf(" Layer %d Point %d dE %f keV dx %f cm mom %f GeV/c\n",iLay,
254 HasPointOnLayer(iLay),
255 GetdESample(iLay),
256 GetdxSample(iLay),
257 GetMomentumAtSample(iLay));
258 }
259
260 printf("Layers used for PID:\n");
261 printf("Layer ");
262 for(Int_t iLay=0; iLay<fNSamples; iLay++){
263 printf("%d ",iLay);
264 }
265 printf("\n");
266
267 printf("Use ");
268 for(Int_t iLay=0; iLay<fNSamples; iLay++){
269 printf("%d ",UseLayerForPid(iLay));
270 }
271 printf("\n");
272 printf("Truncated mean = %f\n",GetTruncatedMean());
273}
274//______________________________________________________________________
275void AliITSdEdxSamples::PrintClusterMap() const{
276 // print the cluster map
277
278 printf("Layer ");
279 for(Int_t iLay=0; iLay<fNSamples; iLay++){
280 printf("%d ",iLay);
281 }
282 printf("\n");
283
284 printf("Point ");
285 for(Int_t iLay=0; iLay<fNSamples; iLay++){
286 printf("%d ",HasPointOnLayer(iLay));
287 }
288 printf("\n");
289}