coverity fix
[u/mrichter/AliRoot.git] / PWGLF / SPECTRA / MultEvShape / AliAnalysisCentralCutESD.cxx
1 /*************************************************************************
2 * Copyright(c) 1998-2008, 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
16 // -----------------------------------------------
17 // particle level cuts for azimuthal isotropic
18 // expansion in highly central collisions analysis
19 // author: Cristian Andrei
20 //         acristian@niham.nipne.ro
21 // ------------------------------------------------
22
23 #include <TF1.h>
24 #include <TFile.h>
25
26 #include "AliAnalysisCentralCutESD.h"
27 #include "AliESDtrack.h"
28
29
30 class TObject;
31 class TParticle;
32
33
34 //____________________________________________________________________
35 ClassImp(AliAnalysisCentralCutESD)
36
37 //____________________________________________________________________
38 AliAnalysisCentralCutESD::AliAnalysisCentralCutESD(const Char_t* name, const Char_t* title) 
39     :AliAnalysisCuts(name,title)
40     ,fReqPID(kFALSE)
41     ,fReqCharge(kFALSE)
42     ,fPartType(kPiPlus)
43     ,fPIDtype("Custom")
44     ,fPriorsFunc(kFALSE)
45     ,fPartPriors()
46     ,fElectronFunction(0)
47     ,fMuonFunction(0)
48     ,fPionFunction(0)
49     ,fKaonFunction(0)
50     ,fProtonFunction(0) 
51 {
52 // Constructor
53 // Initialize the priors
54     fPartPriors[0] = 0.01;
55     fPartPriors[1] = 0.01;
56     fPartPriors[2] = 0.85;
57     fPartPriors[3] = 0.1;
58     fPartPriors[4] = 0.05;
59
60     if(fPriorsFunc){
61         TFile *f = TFile::Open("$ALICE_ROOT/PWG2/data/PriorProbabilities.root ");
62         if(!f){
63             printf("Can't open PWG2 prior probabilities file!\n Exiting ...\n");
64             return;
65         }
66         fElectronFunction = (TF1 *)f->Get("fitElectrons");
67         fMuonFunction = (TF1 *)f->Get("fitMuons");
68         fPionFunction = (TF1 *)f->Get("fitPions");
69         fKaonFunction = (TF1 *)f->Get("fitKaons");
70         fProtonFunction = (TF1 *)f->Get("fitProtons");
71     }
72
73 }
74
75 AliAnalysisCentralCutESD::~AliAnalysisCentralCutESD() {
76 // Destructor
77 // Delete the created priors
78
79
80         if(fElectronFunction) delete fElectronFunction;
81         if(fMuonFunction) delete fMuonFunction;
82         if(fPionFunction) delete fPionFunction;
83         if(fKaonFunction) delete fKaonFunction;
84         if(fProtonFunction) delete fProtonFunction;
85
86 }
87
88
89 Bool_t AliAnalysisCentralCutESD::IsSelected(TObject *obj){
90 // Checks if a particle passes the cuts
91
92     AliESDtrack *track = dynamic_cast<AliESDtrack *>(obj);
93
94     if(!track){
95                 printf("AliAnalysisCentralCutESD:IsSelected ->Can't get track!\n");
96                 return kFALSE;
97     }
98
99     if(fReqCharge){
100                 if(!IsCharged(track)) return kFALSE;
101     } 
102
103     if(fReqPID){
104                 if(!IsA(track, fPartType)) return kFALSE;
105     }
106
107     return kTRUE;
108 }
109
110
111 Double_t AliAnalysisCentralCutESD::GetPriors(Int_t i, Double_t p) {
112 //Return the a priori probs
113
114 Double_t priors=0;
115     if(fPriorsFunc) {
116         if(i == 0) priors = fElectronFunction->Eval(p);
117         if(i == 1) priors = fMuonFunction->Eval(p);
118         if(i == 2) priors = fPionFunction->Eval(p);
119         if(i == 3) priors = fKaonFunction->Eval(p);
120         if(i == 4) priors = fProtonFunction->Eval(p);
121     }
122     else {
123         priors = fPartPriors[i];
124     }
125
126   return priors;
127 }
128
129
130
131 Bool_t AliAnalysisCentralCutESD::IsA(AliESDtrack *track, PDG_t reqPartType){
132 // Determines the type of the particle
133     Int_t charge;
134
135     if(reqPartType < 0){
136       charge = -1;
137     }
138     else{
139       charge = 1;
140     }
141
142     Double_t probability[5] = {0.0,0.0,0.0,0.0,0.0};
143     Double_t w[5] = {0.0,0.0,0.0,0.0,0.0};
144
145     Long64_t partType = 0;
146
147     Double_t p = track->P();
148
149     track->GetESDpid(probability);
150
151     Double_t s = 0.0;
152
153
154     for(Int_t i = 0; i < AliPID::kSPECIES; i++){ 
155                 s += probability[i]*GetPriors(i,p);
156     }
157     
158     if(!s < 0.000001) {
159                 for(Int_t i = 0; i < AliPID::kSPECIES; i++){ 
160                 w[i] = probability[i]*GetPriors(i,p)/s;
161                 }
162     }
163
164
165     if(fPIDtype.Contains("Bayesian")) {
166                 partType = TMath::LocMax(AliPID::kSPECIES,w);
167                 if(partType<0.) return kFALSE;
168     }
169
170     else if(fPIDtype.Contains("Custom")){
171                 for(Int_t i=0;i<AliPID::kSPECIES;i++)   {
172                         if(w[i]>0.9){
173                                 partType = i;
174                         }
175                 }
176     }
177
178     else{
179                 printf("Unknown PID method!\n");
180                 return kFALSE;
181     }
182
183     if(partType<0.) {
184       return kFALSE;
185     }
186     else { // then here partType >=0 so ok for ParticleCode !
187       if((AliPID::ParticleCode(partType)) != reqPartType) return kFALSE;
188     }
189     if(track->Charge() != charge) return kFALSE;
190     
191     return kTRUE;
192
193 }
194
195 Bool_t AliAnalysisCentralCutESD::IsCharged(AliESDtrack* const track) const{
196
197     if(track->Charge() == 0) return kFALSE;
198
199     return kTRUE;
200
201 }