1 /**************************************************************************
2 * Copyright(c) 2009-2012, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
19 ///////////////////////////////////////////////////////////////////
21 // Class to store information for PID with ITS //
22 // and truncated mean computation methods //
23 // Origin: F.Prino, Torino, prino@to.infn.it //
25 ///////////////////////////////////////////////////////////////////
27 #include "AliITSPidParams.h"
28 #include "AliITSdEdxSamples.h"
32 ClassImp(AliITSdEdxSamples)
34 //______________________________________________________________________
35 AliITSdEdxSamples::AliITSdEdxSamples():TObject(),
40 // Default constructor
42 for(Int_t i=0; i<fgkMaxSamples; i++) {
48 //______________________________________________________________________
49 AliITSdEdxSamples::AliITSdEdxSamples(Int_t nSamples,Double_t* samples, Double_t mom, Int_t specie) :
53 fParticleSpecie(specie)
55 // Standard constructor
56 for(Int_t i=0; i<fgkMaxSamples; i++)fPAtSample[i]=0.;
57 SetSamples(nSamples,samples);
59 //______________________________________________________________________
60 void AliITSdEdxSamples::SetSamples(Int_t nSamples, Double_t* samples){
63 if(nSamples>fgkMaxSamples){
64 AliWarning(Form("Too many dE/dx samples,only first %d will be used",fgkMaxSamples));
65 fNSamples=fgkMaxSamples;
69 for(Int_t i=0; i<fNSamples; i++) fdEdxSamples[i]=samples[i];
70 for(Int_t i=fNSamples; i<fgkMaxSamples; i++) fdEdxSamples[i]=0.;
73 //______________________________________________________________________
74 void AliITSdEdxSamples::SetSamplesAndMomenta(Int_t nSamples, Double_t* samples, Double_t* mom){
76 SetSamples(nSamples,samples);
77 for(Int_t i=0; i<fNSamples; i++) fPAtSample[i]=mom[i];
78 for(Int_t i=fNSamples; i<fgkMaxSamples; i++) fPAtSample[i]=0.;
81 //______________________________________________________________________
82 Double_t AliITSdEdxSamples::GetTruncatedMean(Double_t frac, Double_t mindedx) const {
83 // compute truncated mean
86 Double_t dedx[fgkMaxSamples];
87 for (Int_t il=0; il<fNSamples; il++) { // count good (>0) dE/dx values
88 if(fdEdxSamples[il]>mindedx){
89 dedx[nc]= fdEdxSamples[il];
95 Int_t swap; // sort in ascending order
98 for (Int_t i=0; i<nc-1; i++) {
99 if (dedx[i]<=dedx[i+1]) continue;
100 Double_t tmp=dedx[i];
107 Double_t sumamp=0,sumweight=0;
108 Double_t weight[fgkMaxSamples];
109 for(Int_t iw=0; iw<fgkMaxSamples; iw++) weight[iw]=0.;
110 Int_t lastUsed=(Int_t)(frac*nc+0.00001);
111 if(lastUsed==0) lastUsed=1;
112 if(lastUsed>fgkMaxSamples) lastUsed=fgkMaxSamples;
113 for(Int_t iw=0; iw<lastUsed; iw++) weight[iw]=1.;
114 if((frac*nc-lastUsed)>0.4 && lastUsed<fgkMaxSamples) weight[lastUsed]=0.5;
115 for (Int_t i=0; i<nc; i++) {
116 AliDebug(5,Form("dE/dx %f weight %f",dedx[i],weight[i]));
117 sumamp+= dedx[i]*weight[i];
118 sumweight+=weight[i];
120 if(sumweight>0.) return sumamp/sumweight;
124 //______________________________________________________________________
125 Double_t AliITSdEdxSamples::GetWeightedMean(Double_t mindedx) const {
126 // compute generalized mean with k=-2 (used by CMS)
129 Double_t dedx[fgkMaxSamples];
130 for (Int_t il=0; il<fNSamples; il++) { // count good (>0) dE/dx values
131 if(fdEdxSamples[il]>mindedx){
132 dedx[nc]= fdEdxSamples[il];
138 Double_t weiSum = 0.;
139 for (Int_t i=0; i<nc; i++) {
140 weiSum+=TMath::Power(dedx[i],-2);
143 if(weiSum>0.) wMean= TMath::Power(weiSum/nc,-0.5);
147 //______________________________________________________________________
148 void AliITSdEdxSamples::GetConditionalProbabilities(AliITSPidParams* pars, Double_t condprob[AliPID::kSPECIES], Double_t mindedx) const {
149 // compute conditional probablilities
150 const Int_t nPart = 3;
151 Double_t itsProb[nPart] = {1,1,1}; // p, K, pi
153 for(Int_t iS=0; iS<fNSamples; iS++){
154 Int_t iLayer=iS+3; // to match with present ITS
155 if(iLayer>6) iLayer=6; // all extra points are treated as SSD
156 if(fdEdxSamples[iS]<mindedx) continue;
157 Float_t dedx = fdEdxSamples[iS];
158 Float_t layProb = pars->GetLandauGausNorm(dedx,AliPID::kProton,fP,iLayer);
159 itsProb[0] *= layProb;
161 layProb = pars->GetLandauGausNorm(dedx,AliPID::kKaon,fP,iLayer);
162 if (fP < 0.16) layProb=0.00001;
163 itsProb[1] *= layProb;
165 layProb = pars->GetLandauGausNorm(dedx,AliPID::kPion,fP,iLayer);
166 itsProb[2] *= layProb;
169 // Normalise probabilities
170 Double_t sumProb = 0;
171 for (Int_t iPart = 0; iPart < nPart; iPart++) {
172 sumProb += itsProb[iPart];
174 sumProb += 2*itsProb[2]; // muon and electron cannot be distinguished from pions
176 for (Int_t iPart = 0; iPart < nPart; iPart++) {
177 itsProb[iPart]/=sumProb;
180 condprob[AliPID::kElectron] = itsProb[2];
181 condprob[AliPID::kMuon] = itsProb[2];
182 condprob[AliPID::kPion] = itsProb[2];
183 condprob[AliPID::kKaon] = itsProb[1];
184 condprob[AliPID::kProton] = itsProb[0];