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