AOD analysis for protons - analysis train
[u/mrichter/AliRoot.git] / PWG2 / SPECTRA / AliProtonAnalysis.cxx
1 /**************************************************************************
2  * Author: Panos Christakoglou.                                           *
3  * Contributors are mentioned in the code where appropriate.              *
4  *                                                                        *
5  * Permission to use, copy, modify and distribute this software and its   *
6  * documentation strictly for non-commercial purposes is hereby granted   *
7  * without fee, provided that the above copyright notice appears in all   *
8  * copies and that both the copyright notice and this permission notice   *
9  * appear in the supporting documentation. The authors make no claims     *
10  * about the suitability of this software for any purpose. It is          *
11  * provided "as is" without express or implied warranty.                  *
12  **************************************************************************/
13
14 /* $Id$ */
15
16 //-----------------------------------------------------------------
17 //                 AliProtonAnalysis class
18 //   This is the class to deal with the proton analysis
19 //   Origin: Panos Christakoglou, UOA-CERN, Panos.Christakoglou@cern.ch
20 //-----------------------------------------------------------------
21 #include <Riostream.h>
22 #include <TFile.h>
23 #include <TSystem.h>
24 #include <TF1.h>
25 #include <TH2F.h>
26 #include <TH1D.h>
27
28 #include "AliProtonAnalysis.h"
29
30 #include <AliAODEvent.h>
31 #include <AliESDEvent.h>
32 #include <AliLog.h>
33 #include <AliPID.h>
34
35 ClassImp(AliProtonAnalysis)
36
37 //____________________________________________________________________//
38 AliProtonAnalysis::AliProtonAnalysis() : 
39   TObject(), 
40   fNBinsY(0), fMinY(0), fMaxY(0),
41   fNBinsPt(0), fMinPt(0), fMaxPt(0),
42   fMinTPCClusters(0), fMinITSClusters(0),
43   fMaxChi2PerTPCCluster(0), fMaxChi2PerITSCluster(0),
44   fMaxCov11(0), fMaxCov22(0), fMaxCov33(0), fMaxCov44(0), fMaxCov55(0),
45   fMaxSigmaToVertex(0),
46   fMinTPCClustersFlag(kFALSE), fMinITSClustersFlag(kFALSE),
47   fMaxChi2PerTPCClusterFlag(kFALSE), fMaxChi2PerITSClusterFlag(kFALSE),
48   fMaxCov11Flag(kFALSE), fMaxCov22Flag(kFALSE), fMaxCov33Flag(kFALSE), fMaxCov44Flag(kFALSE), fMaxCov55Flag(kFALSE),
49   fMaxSigmaToVertexFlag(kFALSE),
50   fITSRefitFlag(kFALSE), fTPCRefitFlag(kFALSE),
51   fFunctionProbabilityFlag(kFALSE), 
52   fElectronFunction(0), fMuonFunction(0),
53   fPionFunction(0), fKaonFunction(0), fProtonFunction(0),
54   fHistYPtProtons(0), fHistYPtAntiProtons(0) {
55   //Default constructor
56   for(Int_t i = 0; i < 5; i++) fPartFrac[i] = 0.0;
57 }
58
59 //____________________________________________________________________//
60 AliProtonAnalysis::AliProtonAnalysis(Int_t nbinsY, Float_t fLowY, Float_t fHighY,Int_t nbinsPt, Float_t fLowPt, Float_t fHighPt) : 
61   TObject(),
62   fNBinsY(nbinsY), fMinY(fLowY), fMaxY(fHighY),
63   fNBinsPt(nbinsPt), fMinPt(fLowPt), fMaxPt(fHighPt),
64   fMinTPCClusters(0), fMinITSClusters(0),
65   fMaxChi2PerTPCCluster(0), fMaxChi2PerITSCluster(0),
66   fMaxCov11(0), fMaxCov22(0), fMaxCov33(0), fMaxCov44(0), fMaxCov55(0),
67   fMaxSigmaToVertex(0),
68   fMinTPCClustersFlag(kFALSE), fMinITSClustersFlag(kFALSE),
69   fMaxChi2PerTPCClusterFlag(kFALSE), fMaxChi2PerITSClusterFlag(kFALSE),
70   fMaxCov11Flag(kFALSE), fMaxCov22Flag(kFALSE), fMaxCov33Flag(kFALSE), fMaxCov44Flag(kFALSE), fMaxCov55Flag(kFALSE),
71   fMaxSigmaToVertexFlag(kFALSE),
72   fITSRefitFlag(kFALSE), fTPCRefitFlag(kFALSE),
73   fFunctionProbabilityFlag(kFALSE), 
74   fElectronFunction(0), fMuonFunction(0),
75   fPionFunction(0), fKaonFunction(0), fProtonFunction(0),
76   fHistYPtProtons(0), fHistYPtAntiProtons(0) {
77   //Default constructor
78
79   fHistYPtProtons = new TH2F("fHistYPtProtons","y-Pt Protons",fNBinsY,fMinY,fMaxY,fNBinsPt,fMinPt,fMaxPt);
80   fHistYPtProtons->SetStats(kTRUE);
81   fHistYPtProtons->GetYaxis()->SetTitle("P_{T} [GeV]");
82   fHistYPtProtons->GetXaxis()->SetTitle("y");
83   fHistYPtProtons->GetXaxis()->SetTitleColor(1);
84
85   fHistYPtAntiProtons = new TH2F("fHistYPtAntiProtons","y-Pt Antiprotons",fNBinsY,fMinY,fMaxY,fNBinsPt,fMinPt,fMaxPt);
86   fHistYPtAntiProtons->SetStats(kTRUE);
87   fHistYPtAntiProtons->GetYaxis()->SetTitle("P_{T} [GeV]");
88   fHistYPtAntiProtons->GetXaxis()->SetTitle("y");
89   fHistYPtAntiProtons->GetXaxis()->SetTitleColor(1);
90
91
92 //____________________________________________________________________//
93 AliProtonAnalysis::~AliProtonAnalysis() {
94   //Default destructor
95   
96 }
97
98 //____________________________________________________________________//
99 void AliProtonAnalysis::InitHistograms(Int_t nbinsY, Float_t fLowY, Float_t fHighY, Int_t nbinsPt, Float_t fLowPt, Float_t fHighPt) {
100   fNBinsY = nbinsY;
101   fMinY = fLowY;
102   fMaxY = fHighY;
103   fNBinsPt = nbinsPt;
104   fMinPt = fLowPt;
105   fMaxPt = fHighPt;
106
107   fHistYPtProtons = new TH2F("fHistYPtProtons","y-Pt Protons",fNBinsY,fMinY,fMaxY,fNBinsPt,fMinPt,fMaxPt);
108   fHistYPtProtons->SetStats(kTRUE);
109   fHistYPtProtons->GetYaxis()->SetTitle("P_{T} [GeV]");
110   fHistYPtProtons->GetXaxis()->SetTitle("y");
111   fHistYPtProtons->GetXaxis()->SetTitleColor(1);
112
113   fHistYPtAntiProtons = new TH2F("fHistYPtAntiProtons","y-Pt Antiprotons",fNBinsY,fMinY,fMaxY,fNBinsPt,fMinPt,fMaxPt);
114   fHistYPtAntiProtons->SetStats(kTRUE);
115   fHistYPtAntiProtons->GetYaxis()->SetTitle("P_{T} [GeV]");
116   fHistYPtAntiProtons->GetXaxis()->SetTitle("y");
117   fHistYPtAntiProtons->GetXaxis()->SetTitleColor(1);
118 }
119
120 //____________________________________________________________________//
121 void AliProtonAnalysis::ReadFromFile(const char* filename) {
122   TFile *file = TFile::Open(filename);
123   fHistYPtProtons = (TH2F *)file->Get("fHistYPtProtons");
124   fHistYPtAntiProtons = (TH2F *)file->Get("fHistYPtAntiProtons");
125   fHistYPtProtons->Sumw2();
126   fHistYPtAntiProtons->Sumw2();
127 }
128
129 //____________________________________________________________________//
130 TH1D *AliProtonAnalysis::GetProtonYHistogram() {
131   TH1D *fYProtons = (TH1D *)fHistYPtProtons->ProjectionX("fYProtons",0,fHistYPtProtons->GetYaxis()->GetNbins(),"e"); 
132   fYProtons->SetStats(kFALSE);
133   fYProtons->GetYaxis()->SetTitle("dN/dy");
134   fYProtons->SetTitle("dN/dy protons");
135   fYProtons->SetMarkerStyle(kFullCircle);
136   fYProtons->SetMarkerColor(4);
137
138   return fYProtons;
139 }
140
141 //____________________________________________________________________//
142 TH1D *AliProtonAnalysis::GetAntiProtonYHistogram() {
143   TH1D *fYAntiProtons = (TH1D *)fHistYPtAntiProtons->ProjectionX("fYAntiProtons",0,fHistYPtAntiProtons->GetYaxis()->GetNbins(),"e"); 
144   fYAntiProtons->SetStats(kFALSE);
145   fYAntiProtons->GetYaxis()->SetTitle("dN/dy");
146   fYAntiProtons->SetTitle("dN/dy antiprotons");
147   fYAntiProtons->SetMarkerStyle(kFullCircle);
148   fYAntiProtons->SetMarkerColor(4);
149
150   return fYAntiProtons;
151 }
152
153 //____________________________________________________________________//
154 TH1D *AliProtonAnalysis::GetProtonPtHistogram() {
155   TH1D *fPtProtons = (TH1D *)fHistYPtProtons->ProjectionY("fPtProtons",0,fHistYPtProtons->GetXaxis()->GetNbins(),"e"); 
156   fPtProtons->SetStats(kFALSE);
157   fPtProtons->GetYaxis()->SetTitle("dN/dP_{T}");
158   fPtProtons->SetTitle("dN/dPt protons");
159   fPtProtons->SetMarkerStyle(kFullCircle);
160   fPtProtons->SetMarkerColor(4);
161
162   return fPtProtons;
163 }
164
165 //____________________________________________________________________//
166 TH1D *AliProtonAnalysis::GetAntiProtonPtHistogram() {
167   TH1D *fPtAntiProtons = (TH1D *)fHistYPtAntiProtons->ProjectionY("fPtAntiProtons",0,fHistYPtProtons->GetXaxis()->GetNbins(),"e"); 
168   fPtAntiProtons->SetStats(kFALSE);
169   fPtAntiProtons->GetYaxis()->SetTitle("dN/dP_{T}");
170   fPtAntiProtons->SetTitle("dN/dPt antiprotons");
171   fPtAntiProtons->SetMarkerStyle(kFullCircle);
172   fPtAntiProtons->SetMarkerColor(4);
173
174   return fPtAntiProtons;
175 }
176
177 //____________________________________________________________________//
178 TH1D *AliProtonAnalysis::GetYRatioHistogram() {
179   TH1D *fYProtons = GetProtonYHistogram();
180   TH1D *fYAntiProtons = GetAntiProtonYHistogram();
181   
182   TH1D *hRatioY = new TH1D("hRatioY","",fYProtons->GetNbinsX(),fYProtons->GetXaxis()->GetXmin(),fYProtons->GetXaxis()->GetXmax());
183   hRatioY->Divide(fYAntiProtons,fYProtons,1.0,1.0);
184   hRatioY->SetMarkerStyle(kFullCircle);
185   hRatioY->SetMarkerColor(4);
186   hRatioY->GetYaxis()->SetTitle("#bar{p}/p");
187   hRatioY->GetYaxis()->SetTitleOffset(1.4);
188   hRatioY->GetXaxis()->SetTitle("y");
189   hRatioY->GetXaxis()->SetTitleColor(1);
190   hRatioY->SetStats(kFALSE);
191
192   return hRatioY;
193 }
194
195 //____________________________________________________________________//
196 TH1D *AliProtonAnalysis::GetPtRatioHistogram() {
197   TH1D *fPtProtons = GetProtonPtHistogram();
198   TH1D *fPtAntiProtons = GetAntiProtonPtHistogram();
199   
200   TH1D *hRatioPt = new TH1D("hRatioPt","",fPtProtons->GetNbinsX(),fPtProtons->GetXaxis()->GetXmin(),fPtProtons->GetXaxis()->GetXmax());
201   hRatioPt->Divide(fPtAntiProtons,fPtProtons,1.0,1.0);
202   hRatioPt->SetMarkerStyle(kFullCircle);
203   hRatioPt->SetMarkerColor(4);
204   hRatioPt->GetYaxis()->SetTitle("#bar{p}/p");
205   hRatioPt->GetYaxis()->SetTitleOffset(1.4);
206   hRatioPt->GetXaxis()->SetTitle("P_{T} [GeV/c]");
207   hRatioPt->GetXaxis()->SetTitleColor(1);
208   hRatioPt->SetStats(kFALSE);
209
210   return hRatioPt;
211 }
212
213 //____________________________________________________________________//
214 TH1D *AliProtonAnalysis::GetYAsymmetryHistogram() {
215   TH1D *fYProtons = GetProtonYHistogram();
216   TH1D *fYAntiProtons = GetAntiProtonYHistogram();
217   
218   TH1D *hsum = new TH1D("hsumY","",fYProtons->GetNbinsX(),fYProtons->GetXaxis()->GetXmin(),fYProtons->GetXaxis()->GetXmax());
219   hsum->Add(fYProtons,fYAntiProtons,1.0,1.0);
220
221   TH1D *hdiff = new TH1D("hdiffY","",fYProtons->GetNbinsX(),fYProtons->GetXaxis()->GetXmin(),fYProtons->GetXaxis()->GetXmax());
222   hdiff->Add(fYProtons,fYAntiProtons,1.0,-1.0);
223
224   TH1D *hAsymmetryY = new TH1D("hAsymmetryY","",fYProtons->GetNbinsX(),fYProtons->GetXaxis()->GetXmin(),fYProtons->GetXaxis()->GetXmax());
225   hAsymmetryY->Divide(hdiff,hsum,2.0,1.);
226   hAsymmetryY->SetMarkerStyle(kFullCircle);
227   hAsymmetryY->SetMarkerColor(4);
228   hAsymmetryY->GetYaxis()->SetTitle("A_{p}");
229   hAsymmetryY->GetYaxis()->SetTitleOffset(1.4);
230   hAsymmetryY->GetXaxis()->SetTitle("y");
231   hAsymmetryY->GetXaxis()->SetTitleColor(1);
232   hAsymmetryY->SetStats(kFALSE);
233
234   return hAsymmetryY;
235 }
236
237 //____________________________________________________________________//
238 TH1D *AliProtonAnalysis::GetPtAsymmetryHistogram() {
239   TH1D *fPtProtons = GetProtonPtHistogram();
240   TH1D *fPtAntiProtons = GetAntiProtonPtHistogram();
241   
242   TH1D *hsum = new TH1D("hsumPt","",fPtProtons->GetNbinsX(),fPtProtons->GetXaxis()->GetXmin(),fPtProtons->GetXaxis()->GetXmax());
243   hsum->Add(fPtProtons,fPtAntiProtons,1.0,1.0);
244
245   TH1D *hdiff = new TH1D("hdiffPt","",fPtProtons->GetNbinsX(),fPtProtons->GetXaxis()->GetXmin(),fPtProtons->GetXaxis()->GetXmax());
246   hdiff->Add(fPtProtons,fPtAntiProtons,1.0,-1.0);
247
248   TH1D *hAsymmetryPt = new TH1D("hAsymmetryPt","",fPtProtons->GetNbinsX(),fPtProtons->GetXaxis()->GetXmin(),fPtProtons->GetXaxis()->GetXmax());
249   hAsymmetryPt->Divide(hdiff,hsum,2.0,1.);
250   hAsymmetryPt->SetMarkerStyle(kFullCircle);
251   hAsymmetryPt->SetMarkerColor(4);
252   hAsymmetryPt->GetYaxis()->SetTitle("A_{p}");
253   hAsymmetryPt->GetYaxis()->SetTitleOffset(1.4);
254   hAsymmetryPt->GetXaxis()->SetTitle("P_{T} [GeV/c]");
255   hAsymmetryPt->GetXaxis()->SetTitleColor(1);
256   hAsymmetryPt->SetStats(kFALSE);
257
258   return hAsymmetryPt;
259 }
260
261 //____________________________________________________________________//
262 Double_t AliProtonAnalysis::GetParticleFraction(Int_t i, Double_t p) {
263   Double_t partFrac;
264   if(fFunctionProbabilityFlag) {
265     if(i == 0) partFrac = fElectronFunction->Eval(p);
266     if(i == 1) partFrac = fMuonFunction->Eval(p);
267     if(i == 2) partFrac = fPionFunction->Eval(p);
268     if(i == 3) partFrac = fKaonFunction->Eval(p);
269     if(i == 4) partFrac = fProtonFunction->Eval(p);
270   }
271   else partFrac = fPartFrac[i];
272
273   return partFrac;
274 }
275
276 //____________________________________________________________________//
277 void AliProtonAnalysis::Analyze(AliESDEvent* fESD) {
278   //Main analysis part
279   Int_t nGoodTracks = fESD->GetNumberOfTracks();
280   for(Int_t iTracks = 0; iTracks < nGoodTracks; iTracks++) {
281     AliESDtrack* track = fESD->GetTrack(iTracks);
282     if(IsAccepted(track)) {
283       Double_t Pt = track->Pt();
284       Double_t P = track->P();
285         
286       //pid
287       Double_t probability[5];
288           track->GetESDpid(probability);
289       Double_t rcc = 0.0;
290       for(Int_t i = 0; i < AliPID::kSPECIES; i++) rcc += probability[i]*GetParticleFraction(i,P);
291       if(rcc == 0.0) continue;
292       Double_t w[5];
293       for(Int_t i = 0; i < AliPID::kSPECIES; i++) w[i] = probability[i]*GetParticleFraction(i,P)/rcc;
294       Long64_t fParticleType = TMath::LocMax(AliPID::kSPECIES,w);
295       if(fParticleType == 4) {
296         if(track->Charge() > 0) fHistYPtProtons->Fill(Rapidity(track),Pt);
297         else if(track->Charge() < 0) fHistYPtAntiProtons->Fill(Rapidity(track),Pt);
298       }//proton check
299     }//cuts
300   }//track loop 
301 }
302
303 //____________________________________________________________________//
304 void AliProtonAnalysis::Analyze(AliAODEvent* fAOD) {
305   //Main analysis part
306   Int_t nGoodTracks = fAOD->GetNumberOfTracks();
307   for(Int_t iTracks = 0; iTracks < nGoodTracks; iTracks++) {
308     AliAODTrack* track = fAOD->GetTrack(iTracks);
309     Double_t Pt = track->Pt();
310     Double_t P = track->P();
311     
312     //pid
313     Double_t probability[10];
314     track->GetPID(probability);
315     Double_t rcc = 0.0;
316     for(Int_t i = 0; i < AliPID::kSPECIESN; i++) rcc += probability[i]*GetParticleFraction(i,P);
317     if(rcc == 0.0) continue;
318     Double_t w[5];
319     for(Int_t i = 0; i < AliPID::kSPECIESN; i++) w[i] = probability[i]*GetParticleFraction(i,P)/rcc;
320     Long64_t fParticleType = TMath::LocMax(AliPID::kSPECIESN,w);
321     if(fParticleType == 4) {
322       if(track->Charge() > 0) fHistYPtProtons->Fill(track->Y(fParticleType),Pt);
323       else if(track->Charge() < 0) fHistYPtAntiProtons->Fill(track->Y(fParticleType),Pt);
324     }//proton check
325   }//track loop 
326 }
327
328 //____________________________________________________________________//
329 Bool_t AliProtonAnalysis::IsAccepted(AliESDtrack* track) {
330   // Checks if the track is excluded from the cuts
331   Int_t  fIdxInt[200];
332   Int_t nClustersITS = track->GetITSclusters(fIdxInt);
333   Int_t nClustersTPC = track->GetTPCclusters(fIdxInt);
334
335   Float_t chi2PerClusterITS = -1;
336   Float_t chi2PerClusterTPC = -1;
337   if (nClustersTPC!=0)
338     chi2PerClusterTPC = track->GetTPCchi2()/Float_t(nClustersTPC);
339
340   Double_t extCov[15];
341   track->GetExternalCovariance(extCov);
342
343   Double_t Pt = track->Pt();
344
345   if(fMinITSClustersFlag)
346     if(nClustersITS < fMinITSClusters) return kFALSE;
347   if(fMinTPCClustersFlag)
348     if(nClustersTPC < fMinTPCClusters) return kFALSE;
349   if(fMaxChi2PerTPCClusterFlag)
350     if(chi2PerClusterTPC > fMaxChi2PerTPCCluster) return kFALSE; 
351   if(fMaxChi2PerITSClusterFlag)
352     if(chi2PerClusterITS > fMaxChi2PerITSCluster) return kFALSE; 
353   if(fMaxCov11Flag)
354     if(extCov[0] > fMaxCov11) return kFALSE;
355   if(fMaxCov22Flag)
356     if(extCov[2] > fMaxCov22) return kFALSE;
357   if(fMaxCov33Flag)
358     if(extCov[5] > fMaxCov33) return kFALSE;
359   if(fMaxCov44Flag)
360     if(extCov[9] > fMaxCov44) return kFALSE;
361   if(fMaxCov55Flag)
362     if(extCov[14] > fMaxCov55) return kFALSE;
363   if(fMaxSigmaToVertexFlag)
364     if(GetSigmaToVertex(track) > fMaxSigmaToVertex) return kFALSE;
365   if(fITSRefitFlag)
366     if ((track->GetStatus() & AliESDtrack::kITSrefit) == 0) return kFALSE;
367   if(fTPCRefitFlag)
368     if ((track->GetStatus() & AliESDtrack::kTPCrefit) == 0) return kFALSE;
369
370   if((Pt < fMinPt) || (Pt > fMaxPt)) return kFALSE;
371
372   return kTRUE;
373 }
374
375 //____________________________________________________________________//
376 Float_t AliProtonAnalysis::GetSigmaToVertex(AliESDtrack* esdTrack) {
377   // Calculates the number of sigma to the vertex.
378   
379   Float_t b[2];
380   Float_t bRes[2];
381   Float_t bCov[3];
382   esdTrack->GetImpactParameters(b,bCov);
383   if (bCov[0]<=0 || bCov[2]<=0) {
384     //AliDebug(1, "Estimated b resolution lower or equal zero!");
385     bCov[0]=0; bCov[2]=0;
386   }
387   bRes[0] = TMath::Sqrt(bCov[0]);
388   bRes[1] = TMath::Sqrt(bCov[2]);
389   
390   if (bRes[0] == 0 || bRes[1] ==0) return -1;
391   
392   Float_t d = TMath::Sqrt(TMath::Power(b[0]/bRes[0],2) + TMath::Power(b[1]/bRes[1],2));
393   
394   if (TMath::Exp(-d * d / 2) < 1e-10) return 1000;
395   
396   d = TMath::ErfInverse(1 - TMath::Exp(-d * d / 2)) * TMath::Sqrt(2);
397   
398   return d;
399 }
400
401 Double_t AliProtonAnalysis::Rapidity(AliESDtrack *track) {
402   Double_t fMass = 9.38270000000000048e-01;
403   
404   Double_t P = TMath::Sqrt(TMath::Power(track->Px(),2) + 
405                            TMath::Power(track->Py(),2) + 
406                                                    TMath::Power(track->Pz(),2));
407   Double_t energy = TMath::Sqrt(P*P + fMass*fMass);
408   Double_t y = -999;
409   if(energy != track->Pz()) 
410     y = 0.5*TMath::Log((energy + track->Pz())/(energy - track->Pz()));
411
412   return y;
413 }