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