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 | Panos.Christakoglou@cern.ch
  20 //-----------------------------------------------------------------
  21 #include <Riostream.h>
  22 #include <TFile.h>
  23 #include <TSystem.h>
  24 #include <TF1.h>
  25 #include <TH2D.h>
  26 #include <TH1D.h>
  27 #include <TH1I.h>
  28 #include <TParticle.h>
  29
  30 #include "AliProtonAnalysis.h"
  31
  32 #include <AliExternalTrackParam.h>
  33 #include <AliAODEvent.h>
  34 #include <AliESDEvent.h>
  35 #include <AliLog.h>
  36 #include <AliPID.h>
  37 #include <AliStack.h>
  38 #include <AliCFContainer.h>
  39 #include <AliCFEffGrid.h>
  40 #include <AliCFDataGrid.h>
  41 #include <AliESDVertex.h>
  42
  43 ClassImp(AliProtonAnalysis)
  44
  45 //____________________________________________________________________//
  46 AliProtonAnalysis::AliProtonAnalysis() :
  47   TObject(), fAnalysisEtaMode(kFALSE),
  48   fNBinsY(0), fMinY(0), fMaxY(0),
  49   fNBinsPt(0), fMinPt(0), fMaxPt(0),
  50   fMinTPCClusters(0), fMinITSClusters(0),
  51   fMaxChi2PerTPCCluster(0), fMaxChi2PerITSCluster(0),
  52   fMaxCov11(0), fMaxCov22(0), fMaxCov33(0), fMaxCov44(0), fMaxCov55(0),
  53   fMaxSigmaToVertex(0), fMaxSigmaToVertexTPC(0),
  54   fMaxDCAXY(0), fMaxDCAXYTPC(0),
  55   fMaxDCAZ(0), fMaxDCAZTPC(0),
  56   fMaxConstrainChi2(0),
  57   fMinTPCClustersFlag(kFALSE), fMinITSClustersFlag(kFALSE),
  58   fMaxChi2PerTPCClusterFlag(kFALSE), fMaxChi2PerITSClusterFlag(kFALSE),
  59   fMaxCov11Flag(kFALSE), fMaxCov22Flag(kFALSE),
  60   fMaxCov33Flag(kFALSE), fMaxCov44Flag(kFALSE), fMaxCov55Flag(kFALSE),
  61   fMaxSigmaToVertexFlag(kFALSE), fMaxSigmaToVertexTPCFlag(kFALSE),
  62   fMaxDCAXYFlag(kFALSE), fMaxDCAXYTPCFlag(kFALSE),
  63   fMaxDCAZFlag(kFALSE), fMaxDCAZTPCFlag(kFALSE),
  64   fMaxConstrainChi2Flag(kFALSE),
  65   fITSRefitFlag(kFALSE), fTPCRefitFlag(kFALSE),
  66   fESDpidFlag(kFALSE), fTPCpidFlag(kFALSE),
  67   fPointOnITSLayer1Flag(0), fPointOnITSLayer2Flag(0),
  68   fPointOnITSLayer3Flag(0), fPointOnITSLayer4Flag(0),
  69   fPointOnITSLayer5Flag(0), fPointOnITSLayer6Flag(0),
  70   fFunctionProbabilityFlag(kFALSE),
  71   fElectronFunction(0), fMuonFunction(0),
  72   fPionFunction(0), fKaonFunction(0), fProtonFunction(0),
  73   fUseTPCOnly(kFALSE), fUseHybridTPC(kFALSE),
  74   fProtonContainer(0), fAntiProtonContainer(0),
  75   fHistEvents(0), fHistYPtProtons(0), fHistYPtAntiProtons(0),
  76   fEffGridListProtons(0), fCorrectionListProtons2D(0),
  77   fEfficiencyListProtons1D(0), fCorrectionListProtons1D(0),
  78   fEffGridListAntiProtons(0), fCorrectionListAntiProtons2D(0),
  79   fEfficiencyListAntiProtons1D(0), fCorrectionListAntiProtons1D(0),
  80   fCorrectProtons(0), fCorrectAntiProtons(0) {
  81   //Default constructor
  82   for(Int_t i = 0; i < 5; i++) fPartFrac[i] = 0.0;
  83 }
  84
  85 //____________________________________________________________________//
  86 AliProtonAnalysis::AliProtonAnalysis(Int_t nbinsY, Float_t fLowY, Float_t fHighY,Int_t nbinsPt, Float_t fLowPt, Float_t fHighPt) :
  87   TObject(), fAnalysisEtaMode(kFALSE),
  88   fNBinsY(nbinsY), fMinY(fLowY), fMaxY(fHighY),
  89   fNBinsPt(nbinsPt), fMinPt(fLowPt), fMaxPt(fHighPt),
  90   fMinTPCClusters(0), fMinITSClusters(0),
  91   fMaxChi2PerTPCCluster(0), fMaxChi2PerITSCluster(0),
  92   fMaxCov11(0), fMaxCov22(0), fMaxCov33(0), fMaxCov44(0), fMaxCov55(0),
  93   fMaxSigmaToVertex(0), fMaxSigmaToVertexTPC(0),
  94   fMaxDCAXY(0), fMaxDCAXYTPC(0),
  95   fMaxDCAZ(0), fMaxDCAZTPC(0),
  96   fMaxConstrainChi2(0),
  97   fMinTPCClustersFlag(kFALSE), fMinITSClustersFlag(kFALSE),
  98   fMaxChi2PerTPCClusterFlag(kFALSE), fMaxChi2PerITSClusterFlag(kFALSE),
  99   fMaxCov11Flag(kFALSE), fMaxCov22Flag(kFALSE),
 100   fMaxCov33Flag(kFALSE), fMaxCov44Flag(kFALSE), fMaxCov55Flag(kFALSE),
 101   fMaxSigmaToVertexFlag(kFALSE), fMaxSigmaToVertexTPCFlag(kFALSE),
 102   fMaxDCAXYFlag(kFALSE), fMaxDCAXYTPCFlag(kFALSE),
 103   fMaxDCAZFlag(kFALSE), fMaxDCAZTPCFlag(kFALSE),
 104   fMaxConstrainChi2Flag(kFALSE),
 105   fITSRefitFlag(kFALSE), fTPCRefitFlag(kFALSE),
 106   fESDpidFlag(kFALSE), fTPCpidFlag(kFALSE),
 107   fPointOnITSLayer1Flag(0), fPointOnITSLayer2Flag(0),
 108   fPointOnITSLayer3Flag(0), fPointOnITSLayer4Flag(0),
 109   fPointOnITSLayer5Flag(0), fPointOnITSLayer6Flag(0),
 110   fFunctionProbabilityFlag(kFALSE),
 111   fElectronFunction(0), fMuonFunction(0),
 112   fPionFunction(0), fKaonFunction(0), fProtonFunction(0),
 113   fUseTPCOnly(kFALSE), fUseHybridTPC(kFALSE),
 114   fProtonContainer(0), fAntiProtonContainer(0),
 115   fHistEvents(0), fHistYPtProtons(0), fHistYPtAntiProtons(0),
 116   fEffGridListProtons(0), fCorrectionListProtons2D(0),
 117   fEfficiencyListProtons1D(0), fCorrectionListProtons1D(0),
 118   fEffGridListAntiProtons(0), fCorrectionListAntiProtons2D(0),
 119   fEfficiencyListAntiProtons1D(0), fCorrectionListAntiProtons1D(0),
 120   fCorrectProtons(0), fCorrectAntiProtons(0){
 121   //Default constructor
 122   fHistEvents = new TH1I("fHistEvents","Analyzed events",1,0,1);
 123
 124   fHistYPtProtons = new TH2D("fHistYPtProtons","y-Pt Protons",
 125                              fNBinsY,fMinY,fMaxY,fNBinsPt,fMinPt,fMaxPt);
 126   fHistYPtProtons->SetStats(kTRUE);
 127   fHistYPtProtons->GetYaxis()->SetTitle("P_{T} [GeV]");
 128   fHistYPtProtons->GetXaxis()->SetTitle("y");
 129   fHistYPtProtons->GetXaxis()->SetTitleColor(1);
 130
 131   fHistYPtAntiProtons = new TH2D("fHistYPtAntiProtons","y-Pt Antiprotons",
 132                                  fNBinsY,fMinY,fMaxY,fNBinsPt,fMinPt,fMaxPt);
 133   fHistYPtAntiProtons->SetStats(kTRUE);
 134   fHistYPtAntiProtons->GetYaxis()->SetTitle("P_{T} [GeV]");
 135   fHistYPtAntiProtons->GetXaxis()->SetTitle("y");
 136   fHistYPtAntiProtons->GetXaxis()->SetTitleColor(1);
 137
 138   //setting up the containers
 139   Int_t iBin[2];
 140   iBin[0] = nbinsY;
 141   iBin[1] = nbinsPt;
 142   Double_t *binLimY = new Double_t[nbinsY+1];
 143   Double_t *binLimPt = new Double_t[nbinsPt+1];
 144   //values for bin lower bounds
 145   for(Int_t i = 0; i <= nbinsY; i++)
 146     binLimY[i]=(Double_t)fLowY  + (fHighY - fLowY)  /nbinsY*(Double_t)i;
 147   for(Int_t i = 0; i <= nbinsPt; i++)
 148     binLimPt[i]=(Double_t)fLowPt  + (fHighPt - fLowPt)  /nbinsPt*(Double_t)i;
 149
 150   fProtonContainer = new AliCFContainer("containerProtons",
 151                                         "container for protons",
 152                                         1,2,iBin);
 153   fProtonContainer->SetBinLimits(0,binLimY); //rapidity
 154   fProtonContainer->SetBinLimits(1,binLimPt); //pT
 155   fAntiProtonContainer = new AliCFContainer("containerAntiProtons",
 156                                             "container for antiprotons",
 157                                             1,2,iBin);
 158   fAntiProtonContainer->SetBinLimits(0,binLimY); //rapidity
 159   fAntiProtonContainer->SetBinLimits(1,binLimPt); //pT
 160 }
 161
 162 //____________________________________________________________________//
 163 AliProtonAnalysis::~AliProtonAnalysis() {
 164   //Default destructor
 165   if(fHistEvents) delete fHistEvents;
 166   if(fHistYPtProtons) delete fHistYPtProtons;
 167   if(fHistYPtAntiProtons) delete fHistYPtAntiProtons;
 168   if(fProtonContainer) delete fProtonContainer;
 169   if(fAntiProtonContainer) delete fAntiProtonContainer;
 170
 171   if(fEffGridListProtons) delete fEffGridListProtons;
 172   if(fCorrectionListProtons2D) delete fCorrectionListProtons2D;
 173   if(fEfficiencyListProtons1D) delete fEfficiencyListProtons1D;
 174   if(fCorrectionListProtons1D) delete fCorrectionListProtons1D;
 175   if(fEffGridListAntiProtons) delete fEffGridListAntiProtons;
 176   if(fCorrectionListAntiProtons2D) delete fCorrectionListAntiProtons2D;
 177   if(fEfficiencyListAntiProtons1D) delete fEfficiencyListAntiProtons1D;
 178   if(fCorrectionListAntiProtons1D) delete fCorrectionListAntiProtons1D;
 179   if(fCorrectProtons) delete fCorrectProtons;
 180   if(fCorrectAntiProtons) delete fCorrectAntiProtons;
 181 }
 182
 183 //____________________________________________________________________//
 184 void AliProtonAnalysis::InitAnalysisHistograms(Int_t nbinsY,
 185                                                Float_t fLowY, Float_t fHighY,
 186                                                Int_t nbinsPt,
 187                                                Float_t fLowPt, Float_t fHighPt) {
 188   //Initializes the histograms
 189   fNBinsY = nbinsY;
 190   fMinY = fLowY;
 191   fMaxY = fHighY;
 192   fNBinsPt = nbinsPt;
 193   fMinPt = fLowPt;
 194   fMaxPt = fHighPt;
 195
 196   fHistEvents = new TH1I("fHistEvents","Anallyzed events",1,0,1);
 197
 198   fHistYPtProtons = new TH2D("fHistYPtProtons","y-Pt Protons",
 199                              fNBinsY,fMinY,fMaxY,fNBinsPt,fMinPt,fMaxPt);
 200   fHistYPtProtons->SetStats(kTRUE);
 201   fHistYPtProtons->GetYaxis()->SetTitle("P_{T} [GeV]");
 202   fHistYPtProtons->GetXaxis()->SetTitle("y");
 203   fHistYPtProtons->GetXaxis()->SetTitleColor(1);
 204
 205   fHistYPtAntiProtons = new TH2D("fHistYPtAntiProtons","y-Pt Antiprotons",
 206                                  fNBinsY,fMinY,fMaxY,fNBinsPt,fMinPt,fMaxPt);
 207   fHistYPtAntiProtons->SetStats(kTRUE);
 208   fHistYPtAntiProtons->GetYaxis()->SetTitle("P_{T} [GeV]");
 209   fHistYPtAntiProtons->GetXaxis()->SetTitle("y");
 210   fHistYPtAntiProtons->GetXaxis()->SetTitleColor(1);
 211
 212   //setting up the containers
 213   Int_t iBin[2];
 214   iBin[0] = nbinsY;
 215   iBin[1] = nbinsPt;
 216   Double_t *binLimY = new Double_t[nbinsY+1];
 217   Double_t *binLimPt = new Double_t[nbinsPt+1];
 218   //values for bin lower bounds
 219   for(Int_t i = 0; i <= nbinsY; i++)
 220     binLimY[i]=(Double_t)fLowY  + (fHighY - fLowY)  /nbinsY*(Double_t)i;
 221   for(Int_t i = 0; i <= nbinsPt; i++)
 222     binLimPt[i]=(Double_t)fLowPt  + (fHighPt - fLowPt)  /nbinsPt*(Double_t)i;
 223
 224   fProtonContainer = new AliCFContainer("containerProtons",
 225                                         "container for protons",
 226                                         1,2,iBin);
 227   fProtonContainer->SetBinLimits(0,binLimY); //rapidity
 228   fProtonContainer->SetBinLimits(1,binLimPt); //pT
 229   fAntiProtonContainer = new AliCFContainer("containerAntiProtons",
 230                                             "container for antiprotons",
 231                                             1,2,iBin);
 232   fAntiProtonContainer->SetBinLimits(0,binLimY); //rapidity
 233   fAntiProtonContainer->SetBinLimits(1,binLimPt); //pT
 234 }
 235
 236 //____________________________________________________________________//
 237 Bool_t AliProtonAnalysis::ReadFromFile(const char* filename) {
 238   //Read the containers from the existing file
 239   Bool_t status = kTRUE;
 240
 241   TFile *file = TFile::Open(filename);
 242   if(!file) {
 243     cout<<"Could not find the input file "<<filename<<endl;
 244     status = kFALSE;
 245   }
 246
 247   TList *list = (TList *)file->Get("outputList1");
 248   if(list) {
 249     cout<<"Retrieving objects from the list "<<list->GetName()<<"..."<<endl;
 250     fHistYPtProtons = (TH2D *)list->At(0);
 251     fHistYPtAntiProtons = (TH2D *)list->At(1);
 252     fHistEvents = (TH1I *)list->At(2);
 253     fProtonContainer = (AliCFContainer *)list->At(3);
 254     fAntiProtonContainer = (AliCFContainer *)list->At(4);
 255   }
 256   else if(!list) {
 257     cout<<"Retrieving objects from the file... "<<endl;
 258     fHistYPtProtons = (TH2D *)file->Get("fHistYPtProtons");
 259     fHistYPtAntiProtons = (TH2D *)file->Get("fHistYPtAntiProtons");
 260     fHistEvents = (TH1I *)file->Get("fHistEvents");
 261     fProtonContainer = (AliCFContainer *)file->Get("containerProtons");
 262     fAntiProtonContainer = (AliCFContainer *)file->Get("containerAntiProtons");
 263   }
 264   if((!fHistYPtProtons)||(!fHistYPtAntiProtons)||(!fHistEvents)
 265      ||(!fProtonContainer)||(!fAntiProtonContainer)) {
 266     cout<<"Input containers were not found!!!"<<endl;
 267     status = kFALSE;
 268   }
 269   else {
 270     //fHistYPtProtons = fProtonContainer->ShowProjection(0,1,0);
 271     //fHistYPtAntiProtons = fAntiProtonContainer->ShowProjection(0,1,0);
 272     fHistYPtProtons->Sumw2();
 273     fHistYPtAntiProtons->Sumw2();
 274   }
 275
 276   return status;
 277 }
 278
 279 //____________________________________________________________________//
 280 TH1D *AliProtonAnalysis::GetProtonYHistogram() {
 281   //Get the y histogram for protons
 282   Int_t nAnalyzedEvents = GetNumberOfAnalyzedEvents();
 283
 284   TH1D *fYProtons = (TH1D *)fHistYPtProtons->ProjectionX("fYProtons",0,fHistYPtProtons->GetYaxis()->GetNbins(),"");
 285   //TH1D *fYProtons = fProtonContainer->ShowProjection(0,0); //variable-step
 286
 287   fYProtons->SetStats(kFALSE);
 288   fYProtons->GetYaxis()->SetTitle("(1/N_{events})(dN/dy)");
 289   fYProtons->SetTitle("dN/dy protons");
 290   fYProtons->SetMarkerStyle(kFullCircle);
 291   fYProtons->SetMarkerColor(4);
 292   if(nAnalyzedEvents > 0)
 293   fYProtons->Scale(1./nAnalyzedEvents);
 294
 295   return fYProtons;
 296 }
 297
 298 //____________________________________________________________________//
 299 TH1D *AliProtonAnalysis::GetAntiProtonYHistogram() {
 300   //Get the y histogram for antiprotons
 301   Int_t nAnalyzedEvents = GetNumberOfAnalyzedEvents();
 302
 303   TH1D *fYAntiProtons = (TH1D *)fHistYPtAntiProtons->ProjectionX("fYAntiProtons",0,fHistYPtAntiProtons->GetYaxis()->GetNbins(),"");
 304   //TH1D *fYAntiProtons = fAntiProtonContainer->ShowProjection(0,0);//variable-step
 305
 306   fYAntiProtons->SetStats(kFALSE);
 307   fYAntiProtons->GetYaxis()->SetTitle("(1/N_{events})(dN/dy)");
 308   fYAntiProtons->SetTitle("dN/dy antiprotons");
 309   fYAntiProtons->SetMarkerStyle(kFullCircle);
 310   fYAntiProtons->SetMarkerColor(4);
 311   if(nAnalyzedEvents > 0)
 312     fYAntiProtons->Scale(1./nAnalyzedEvents);
 313
 314   return fYAntiProtons;
 315 }
 316
 317 //____________________________________________________________________//
 318 TH1D *AliProtonAnalysis::GetProtonPtHistogram() {
 319   //Get the Pt histogram for protons
 320   Int_t nAnalyzedEvents = GetNumberOfAnalyzedEvents();
 321
 322   TH1D *fPtProtons = (TH1D *)fHistYPtProtons->ProjectionY("fPtProtons",0,fHistYPtProtons->GetXaxis()->GetNbins(),"");
 323   //TH1D *fPtProtons = fProtonContainer->ShowProjection(1,0); //variable-step
 324
 325   fPtProtons->SetStats(kFALSE);
 326   fPtProtons->GetYaxis()->SetTitle("(1/N_{events})(dN/dP_{T})");
 327   fPtProtons->SetTitle("dN/dPt protons");
 328   fPtProtons->SetMarkerStyle(kFullCircle);
 329   fPtProtons->SetMarkerColor(4);
 330   if(nAnalyzedEvents > 0)
 331     fPtProtons->Scale(1./nAnalyzedEvents);
 332
 333   return fPtProtons;
 334 }
 335
 336 //____________________________________________________________________//
 337 TH1D *AliProtonAnalysis::GetAntiProtonPtHistogram() {
 338   //Get the Pt histogram for antiprotons
 339   Int_t nAnalyzedEvents = GetNumberOfAnalyzedEvents();
 340
 341   TH1D *fPtAntiProtons = (TH1D *)fHistYPtAntiProtons->ProjectionY("fPtAntiProtons",0,fHistYPtProtons->GetXaxis()->GetNbins(),"");
 342   //TH1D *fPtAntiProtons = fAntiProtonContainer->ShowProjection(1,0); //variable-step
 343
 344   fPtAntiProtons->SetStats(kFALSE);
 345   fPtAntiProtons->GetYaxis()->SetTitle("(1/N_{events})(dN/dP_{T})");
 346   fPtAntiProtons->SetTitle("dN/dPt antiprotons");
 347   fPtAntiProtons->SetMarkerStyle(kFullCircle);
 348   fPtAntiProtons->SetMarkerColor(4);
 349   if(nAnalyzedEvents > 0)
 350     fPtAntiProtons->Scale(1./nAnalyzedEvents);
 351
 352   return fPtAntiProtons;
 353 }
 354
 355 //____________________________________________________________________//
 356 TH1D *AliProtonAnalysis::GetProtonCorrectedYHistogram() {
 357   //Get the corrected y histogram for protons
 358   Int_t nAnalyzedEvents = GetNumberOfAnalyzedEvents();
 359
 360   TH1D *fYProtons = fCorrectProtons->Project(0); //0: rapidity
 361
 362   fYProtons->SetStats(kFALSE);
 363   fYProtons->GetYaxis()->SetTitle("(1/N_{events})(dN/dy)");
 364   fYProtons->GetXaxis()->SetTitle("y");
 365   fYProtons->SetTitle("dN/dy protons");
 366   fYProtons->SetMarkerStyle(kFullCircle);
 367   fYProtons->SetMarkerColor(4);
 368   if(nAnalyzedEvents > 0)
 369     fYProtons->Scale(1./nAnalyzedEvents);
 370
 371   return fYProtons;
 372 }
 373
 374 //____________________________________________________________________//
 375 TH1D *AliProtonAnalysis::GetAntiProtonCorrectedYHistogram() {
 376   //Get the corrected y histogram for antiprotons
 377   Int_t nAnalyzedEvents = GetNumberOfAnalyzedEvents();
 378
 379   TH1D *fYAntiProtons = fCorrectAntiProtons->Project(0); //0: rapidity
 380
 381   fYAntiProtons->SetStats(kFALSE);
 382   fYAntiProtons->GetYaxis()->SetTitle("(1/N_{events})(dN/dy)");
 383   fYAntiProtons->GetXaxis()->SetTitle("y");
 384   fYAntiProtons->SetTitle("dN/dy protons");
 385   fYAntiProtons->SetMarkerStyle(kFullCircle);
 386   fYAntiProtons->SetMarkerColor(4);
 387   if(nAnalyzedEvents > 0)
 388     fYAntiProtons->Scale(1./nAnalyzedEvents);
 389
 390   return fYAntiProtons;
 391 }
 392
 393 //____________________________________________________________________//
 394 TH1D *AliProtonAnalysis::GetProtonCorrectedPtHistogram() {
 395   //Get the corrected Pt histogram for protons
 396   Int_t nAnalyzedEvents = GetNumberOfAnalyzedEvents();
 397
 398   TH1D *fPtProtons = fCorrectProtons->Project(0); //0: rapidity
 399
 400   fPtProtons->SetStats(kFALSE);
 401   fPtProtons->GetYaxis()->SetTitle("(1/N_{events})(dN/dP_{T})");
 402   fPtProtons->GetXaxis()->SetTitle("P_{T} [GeV/c]");
 403   fPtProtons->SetTitle("dN/dPt protons");
 404   fPtProtons->SetMarkerStyle(kFullCircle);
 405   fPtProtons->SetMarkerColor(4);
 406   if(nAnalyzedEvents > 0)
 407     fPtProtons->Scale(1./nAnalyzedEvents);
 408
 409   return fPtProtons;
 410 }
 411
 412 //____________________________________________________________________//
 413 TH1D *AliProtonAnalysis::GetAntiProtonCorrectedPtHistogram() {
 414   //Get the corrected Pt histogram for antiprotons
 415   Int_t nAnalyzedEvents = GetNumberOfAnalyzedEvents();
 416
 417   TH1D *fPtAntiProtons = fCorrectAntiProtons->Project(0); //0: rapidity
 418
 419   fPtAntiProtons->SetStats(kFALSE);
 420   fPtAntiProtons->GetYaxis()->SetTitle("(1/N_{events})(dN/dP_{T})");
 421   fPtAntiProtons->GetXaxis()->SetTitle("P_{T} [GeV/c]");
 422   fPtAntiProtons->SetTitle("dN/dPt antiprotons");
 423   fPtAntiProtons->SetMarkerStyle(kFullCircle);
 424   fPtAntiProtons->SetMarkerColor(4);
 425   if(nAnalyzedEvents > 0)
 426     fPtAntiProtons->Scale(1./nAnalyzedEvents);
 427
 428   return fPtAntiProtons;
 429 }
 430
 431 //____________________________________________________________________//
 432 TH1D *AliProtonAnalysis::GetYRatioHistogram() {
 433   //Returns the rapidity dependence of the ratio (uncorrected)
 434   TH1D *fYProtons = GetProtonYHistogram();
 435   TH1D *fYAntiProtons = GetAntiProtonYHistogram();
 436
 437   TH1D *hRatioY = new TH1D("hRatioY","",fYProtons->GetNbinsX(),fYProtons->GetXaxis()->GetXmin(),fYProtons->GetXaxis()->GetXmax());
 438   hRatioY->Divide(fYAntiProtons,fYProtons,1.0,1.0);
 439   hRatioY->SetMarkerStyle(kFullCircle);
 440   hRatioY->SetMarkerColor(4);
 441   hRatioY->GetYaxis()->SetTitle("#bar{p}/p");
 442   hRatioY->GetYaxis()->SetTitleOffset(1.4);
 443   hRatioY->GetXaxis()->SetTitle("y");
 444   hRatioY->GetXaxis()->SetTitleColor(1);
 445   hRatioY->SetStats(kFALSE);
 446
 447   return hRatioY;
 448 }
 449
 450 //____________________________________________________________________//
 451 TH1D *AliProtonAnalysis::GetYRatioCorrectedHistogram(TH2D *gCorrectionProtons,
 452                                                      TH2D *gCorrectionAntiProtons) {
 453   //Returns the rapidity dependence of the ratio (corrected)
 454   fHistYPtProtons->Multiply(gCorrectionProtons);
 455   TH1D *fYProtons = GetProtonYHistogram();
 456   fHistYPtAntiProtons->Multiply(gCorrectionAntiProtons);
 457   TH1D *fYAntiProtons = GetAntiProtonYHistogram();
 458
 459   TH1D *hRatioY = new TH1D("hRatioY","",fYProtons->GetNbinsX(),fYProtons->GetXaxis()->GetXmin(),fYProtons->GetXaxis()->GetXmax());
 460   hRatioY->Divide(fYAntiProtons,fYProtons,1.0,1.0);
 461   hRatioY->SetMarkerStyle(kFullCircle);
 462   hRatioY->SetMarkerColor(4);
 463   hRatioY->GetYaxis()->SetTitle("#bar{p}/p");
 464   hRatioY->GetYaxis()->SetTitleOffset(1.4);
 465   hRatioY->GetXaxis()->SetTitle("y");
 466   hRatioY->GetXaxis()->SetTitleColor(1);
 467   hRatioY->SetStats(kFALSE);
 468
 469   return hRatioY;
 470 }
 471
 472 //____________________________________________________________________//
 473 TH1D *AliProtonAnalysis::GetPtRatioHistogram() {
 474   //Returns the pT dependence of the ratio (uncorrected)
 475   TH1D *fPtProtons = GetProtonPtHistogram();
 476   TH1D *fPtAntiProtons = GetAntiProtonPtHistogram();
 477
 478   TH1D *hRatioPt = new TH1D("hRatioPt","",fPtProtons->GetNbinsX(),fPtProtons->GetXaxis()->GetXmin(),fPtProtons->GetXaxis()->GetXmax());
 479   hRatioPt->Divide(fPtAntiProtons,fPtProtons,1.0,1.0);
 480   hRatioPt->SetMarkerStyle(kFullCircle);
 481   hRatioPt->SetMarkerColor(4);
 482   hRatioPt->GetYaxis()->SetTitle("#bar{p}/p");
 483   hRatioPt->GetYaxis()->SetTitleOffset(1.4);
 484   hRatioPt->GetXaxis()->SetTitle("P_{T} [GeV/c]");
 485   hRatioPt->GetXaxis()->SetTitleColor(1);
 486   hRatioPt->SetStats(kFALSE);
 487
 488   return hRatioPt;
 489 }
 490
 491 //____________________________________________________________________//
 492 TH1D *AliProtonAnalysis::GetPtRatioCorrectedHistogram(TH2D *gCorrectionProtons,
 493                                                       TH2D *gCorrectionAntiProtons) {
 494   //Returns the Pt dependence of the ratio (corrected)
 495   fHistYPtProtons->Multiply(gCorrectionProtons);
 496   TH1D *fPtProtons = GetProtonPtHistogram();
 497   fHistYPtAntiProtons->Multiply(gCorrectionAntiProtons);
 498   TH1D *fPtAntiProtons = GetAntiProtonPtHistogram();
 499
 500   TH1D *hRatioPt = new TH1D("hRatioPt","",fPtProtons->GetNbinsX(),fPtProtons->GetXaxis()->GetXmin(),fPtProtons->GetXaxis()->GetXmax());
 501   hRatioPt->Divide(fPtAntiProtons,fPtProtons,1.0,1.0);
 502   hRatioPt->SetMarkerStyle(kFullCircle);
 503   hRatioPt->SetMarkerColor(4);
 504   hRatioPt->GetYaxis()->SetTitle("#bar{p}/p");
 505   hRatioPt->GetYaxis()->SetTitleOffset(1.4);
 506   hRatioPt->GetXaxis()->SetTitle("y");
 507   hRatioPt->GetXaxis()->SetTitleColor(1);
 508   hRatioPt->SetStats(kFALSE);
 509
 510   return hRatioPt;
 511 }
 512
 513 //____________________________________________________________________//
 514 TH1D *AliProtonAnalysis::GetYAsymmetryHistogram() {
 515   //Returns the rapidity dependence of the asymmetry (uncorrected)
 516   TH1D *fYProtons = GetProtonYHistogram();
 517   TH1D *fYAntiProtons = GetAntiProtonYHistogram();
 518
 519   TH1D *hsum = new TH1D("hsumY","",fYProtons->GetNbinsX(),fYProtons->GetXaxis()->GetXmin(),fYProtons->GetXaxis()->GetXmax());
 520   hsum->Add(fYProtons,fYAntiProtons,1.0,1.0);
 521
 522   TH1D *hdiff = new TH1D("hdiffY","",fYProtons->GetNbinsX(),fYProtons->GetXaxis()->GetXmin(),fYProtons->GetXaxis()->GetXmax());
 523   hdiff->Add(fYProtons,fYAntiProtons,1.0,-1.0);
 524
 525   TH1D *hAsymmetryY = new TH1D("hAsymmetryY","",fYProtons->GetNbinsX(),fYProtons->GetXaxis()->GetXmin(),fYProtons->GetXaxis()->GetXmax());
 526   hAsymmetryY->Divide(hdiff,hsum,2.0,1.);
 527   hAsymmetryY->SetMarkerStyle(kFullCircle);
 528   hAsymmetryY->SetMarkerColor(4);
 529   hAsymmetryY->GetYaxis()->SetTitle("A_{p}");
 530   hAsymmetryY->GetYaxis()->SetTitleOffset(1.4);
 531   hAsymmetryY->GetXaxis()->SetTitle("y");
 532   hAsymmetryY->GetXaxis()->SetTitleColor(1);
 533   hAsymmetryY->SetStats(kFALSE);
 534
 535   return hAsymmetryY;
 536 }
 537
 538 //____________________________________________________________________//
 539 TH1D *AliProtonAnalysis::GetPtAsymmetryHistogram() {
 540   //Returns the pT dependence of the asymmetry (uncorrected)
 541   TH1D *fPtProtons = GetProtonPtHistogram();
 542   TH1D *fPtAntiProtons = GetAntiProtonPtHistogram();
 543
 544   TH1D *hsum = new TH1D("hsumPt","",fPtProtons->GetNbinsX(),fPtProtons->GetXaxis()->GetXmin(),fPtProtons->GetXaxis()->GetXmax());
 545   hsum->Add(fPtProtons,fPtAntiProtons,1.0,1.0);
 546
 547   TH1D *hdiff = new TH1D("hdiffPt","",fPtProtons->GetNbinsX(),fPtProtons->GetXaxis()->GetXmin(),fPtProtons->GetXaxis()->GetXmax());
 548   hdiff->Add(fPtProtons,fPtAntiProtons,1.0,-1.0);
 549
 550   TH1D *hAsymmetryPt = new TH1D("hAsymmetryPt","",fPtProtons->GetNbinsX(),fPtProtons->GetXaxis()->GetXmin(),fPtProtons->GetXaxis()->GetXmax());
 551   hAsymmetryPt->Divide(hdiff,hsum,2.0,1.);
 552   hAsymmetryPt->SetMarkerStyle(kFullCircle);
 553   hAsymmetryPt->SetMarkerColor(4);
 554   hAsymmetryPt->GetYaxis()->SetTitle("A_{p}");
 555   hAsymmetryPt->GetYaxis()->SetTitleOffset(1.4);
 556   hAsymmetryPt->GetXaxis()->SetTitle("P_{T} [GeV/c]");
 557   hAsymmetryPt->GetXaxis()->SetTitleColor(1);
 558   hAsymmetryPt->SetStats(kFALSE);
 559
 560   return hAsymmetryPt;
 561 }
 562
 563 //____________________________________________________________________//
 564 Double_t AliProtonAnalysis::GetParticleFraction(Int_t i, Double_t p) {
 565   //Return the a priori probs
 566   Double_t partFrac=0;
 567   if(fFunctionProbabilityFlag) {
 568     if(i == 0) partFrac = fElectronFunction->Eval(p);
 569     if(i == 1) partFrac = fMuonFunction->Eval(p);
 570     if(i == 2) partFrac = fPionFunction->Eval(p);
 571     if(i == 3) partFrac = fKaonFunction->Eval(p);
 572     if(i == 4) partFrac = fProtonFunction->Eval(p);
 573   }
 574   else partFrac = fPartFrac[i];
 575
 576   return partFrac;
 577 }
 578
 579 //____________________________________________________________________//
 580 void AliProtonAnalysis::Analyze(AliESDEvent* esd,
 581                                 const AliESDVertex *vertex) {
 582   //Main analysis part - ESD
 583   fHistEvents->Fill(0); //number of analyzed events
 584   Double_t containerInput[2] ;
 585   Double_t gPt = 0.0, gP = 0.0;
 586   Int_t nGoodTracks = esd->GetNumberOfTracks();
 587   for(Int_t iTracks = 0; iTracks < nGoodTracks; iTracks++) {
 588     AliESDtrack* track = esd->GetTrack(iTracks);
 589     Double_t probability[5];
 590     AliESDtrack trackTPC;
 591
 592     //in case it's a TPC only track relate it to the proper vertex
 593     if((fUseTPCOnly)&&(!fUseHybridTPC)) {
 594       Float_t p[2],cov[3];
 595       track->GetImpactParametersTPC(p,cov);
 596       if (p[0]==0 && p[1]==0)
 597         track->RelateToVertexTPC(((AliESDEvent*)esd)->GetPrimaryVertexTPC(),esd->GetMagneticField(),kVeryBig);
 598       if (!track->FillTPCOnlyTrack(trackTPC)) {
 599         continue;
 600       }
 601       track = &trackTPC ;
 602     }
 603
 604     if(IsAccepted(esd,vertex,track)) {
 605       if(fUseTPCOnly) {
 606         AliExternalTrackParam *tpcTrack = (AliExternalTrackParam *)track->GetTPCInnerParam();
 607         if(!tpcTrack) continue;
 608         gPt = tpcTrack->Pt();
 609         gP = tpcTrack->P();
 610
 611         //pid
 612         track->GetTPCpid(probability);
 613         Double_t rcc = 0.0;
 614         for(Int_t i = 0; i < AliPID::kSPECIES; i++)
 615           rcc += probability[i]*GetParticleFraction(i,gP);
 616         if(rcc == 0.0) continue;
 617         Double_t w[5];
 618         for(Int_t i = 0; i < AliPID::kSPECIES; i++)
 619           w[i] = probability[i]*GetParticleFraction(i,gP)/rcc;
 620         Long64_t fParticleType = TMath::LocMax(AliPID::kSPECIES,w);
 621         if(fParticleType == 4) {
 622           if(tpcTrack->Charge() > 0) {
 623             fHistYPtProtons->Fill(Rapidity(tpcTrack->Px(),
 624                                            tpcTrack->Py(),
 625                                            tpcTrack->Pz()),
 626                                   gPt);
 627             //fill the container
 628             containerInput[0] = Rapidity(tpcTrack->Px(),
 629                                          tpcTrack->Py(),
 630                                          tpcTrack->Pz());
 631             containerInput[1] = gPt;
 632             fProtonContainer->Fill(containerInput,0);
 633           }//protons
 634           else if(tpcTrack->Charge() < 0) {
 635             fHistYPtAntiProtons->Fill(Rapidity(tpcTrack->Px(),
 636                                                tpcTrack->Py(),
 637                                                tpcTrack->Pz()),
 638                                       gPt);
 639             //fill the container
 640             containerInput[0] = Rapidity(tpcTrack->Px(),
 641                                          tpcTrack->Py(),
 642                                          tpcTrack->Pz());
 643             containerInput[1] = gPt;
 644             fAntiProtonContainer->Fill(containerInput,0);
 645           }//antiprotons
 646         }//proton check
 647       }//TPC only tracks
 648       else if(!fUseTPCOnly) {
 649         gPt = track->Pt();
 650         gP = track->P();
 651
 652         //pid
 653         track->GetESDpid(probability);
 654         Double_t rcc = 0.0;
 655         for(Int_t i = 0; i < AliPID::kSPECIES; i++)
 656           rcc += probability[i]*GetParticleFraction(i,gP);
 657         if(rcc == 0.0) continue;
 658         Double_t w[5];
 659         for(Int_t i = 0; i < AliPID::kSPECIES; i++)
 660           w[i] = probability[i]*GetParticleFraction(i,gP)/rcc;
 661         Long64_t fParticleType = TMath::LocMax(AliPID::kSPECIES,w);
 662         if(fParticleType == 4) {
 663           //cout<<"(Anti)protons found..."<<endl;
 664           if(track->Charge() > 0) {
 665             fHistYPtProtons->Fill(Rapidity(track->Px(),
 666                                            track->Py(),
 667                                            track->Pz()),
 668                                   gPt);
 669             //fill the container
 670             containerInput[0] = Rapidity(track->Px(),
 671                                          track->Py(),
 672                                          track->Pz());
 673             containerInput[1] = gPt;
 674             fProtonContainer->Fill(containerInput,0);
 675           }//protons
 676           else if(track->Charge() < 0) {
 677             fHistYPtAntiProtons->Fill(Rapidity(track->Px(),
 678                                                track->Py(),
 679                                                track->Pz()),
 680                                       gPt);
 681             //fill the container
 682             containerInput[0] = Rapidity(track->Px(),
 683                                          track->Py(),
 684                                          track->Pz());
 685             containerInput[1] = gPt;
 686             fAntiProtonContainer->Fill(containerInput,0);
 687           }//antiprotons
 688         }//proton check
 689       }//combined tracking
 690     }//cuts
 691   }//track loop
 692 }
 693
 694 //____________________________________________________________________//
 695 void AliProtonAnalysis::Analyze(AliAODEvent* const fAOD) {
 696   //Main analysis part - AOD
 697   fHistEvents->Fill(0); //number of analyzed events
 698   Int_t nGoodTracks = fAOD->GetNumberOfTracks();
 699   for(Int_t iTracks = 0; iTracks < nGoodTracks; iTracks++) {
 700     AliAODTrack* track = fAOD->GetTrack(iTracks);
 701     Double_t gPt = track->Pt();
 702     Double_t gP = track->P();
 703
 704     //pid
 705     Double_t probability[10];
 706     track->GetPID(probability);
 707     Double_t rcc = 0.0;
 708     for(Int_t i = 0; i < AliPID::kSPECIESN; i++) rcc += probability[i]*GetParticleFraction(i,gP);
 709     if(rcc == 0.0) continue;
 710     Double_t w[10];
 711     for(Int_t i = 0; i < AliPID::kSPECIESN; i++) w[i] = probability[i]*GetParticleFraction(i,gP)/rcc;
 712     Long64_t fParticleType = TMath::LocMax(AliPID::kSPECIESN,w);
 713     if(fParticleType == 4) {
 714       if(track->Charge() > 0)
 715         fHistYPtProtons->Fill(track->Y(fParticleType),gPt);
 716       else if(track->Charge() < 0)
 717         fHistYPtAntiProtons->Fill(track->Y(fParticleType),gPt);
 718     }//proton check
 719   }//track loop
 720 }
 721
 722 //____________________________________________________________________//
 723 void AliProtonAnalysis::Analyze(AliStack* const stack,
 724                                 Bool_t iInclusive) {
 725   //Main analysis part - MC
 726   fHistEvents->Fill(0); //number of analyzed events
 727
 728   Int_t nParticles = 0;
 729   //inclusive protons -
 730   if(iInclusive) nParticles = stack->GetNtrack();
 731   else nParticles = stack->GetNprimary();
 732
 733   for(Int_t i = 0; i < nParticles; i++) {
 734     TParticle *particle = stack->Particle(i);
 735     if(!particle) continue;
 736
 737     //in case of inclusive protons reject the secondaries from hadronic inter.
 738     if(particle->GetUniqueID() == 13) continue;
 739
 740     if(TMath::Abs(particle->Eta()) > 1.0) continue;
 741     if((particle->Pt() > fMaxPt)||(particle->Pt() < fMinPt)) continue;
 742     if((Rapidity(particle->Px(),particle->Py(),particle->Pz()) > fMaxY)||(Rapidity(particle->Px(),particle->Py(),particle->Pz()) < fMinY)) continue;
 743
 744     Int_t pdgcode = particle->GetPdgCode();
 745     if(pdgcode == 2212) fHistYPtProtons->Fill(Rapidity(particle->Px(),
 746                                                        particle->Py(),
 747                                                        particle->Pz()),
 748                                               particle->Pt());
 749     if(pdgcode == -2212) fHistYPtAntiProtons->Fill(Rapidity(particle->Px(),
 750                                                             particle->Py(),
 751                                                             particle->Pz()),
 752                                                    particle->Pt());
 753   }//particle loop
 754 }
 755
 756 //____________________________________________________________________//
 757 Bool_t AliProtonAnalysis::IsInPhaseSpace(AliESDtrack* const track) {
 758   // Checks if the track is outside the analyzed y-Pt phase space
 759   Double_t gPt = 0.0, gPx = 0.0, gPy = 0.0, gPz = 0.0;
 760   Double_t eta = 0.0;
 761
 762   if(fUseTPCOnly) {
 763     AliExternalTrackParam *tpcTrack = (AliExternalTrackParam *)track->GetTPCInnerParam();
 764     if(!tpcTrack) {
 765       gPt = 0.0; gPx = 0.0; gPy = 0.0; gPz = 0.0; eta = -10.0;
 766     }
 767     else {
 768       gPt = tpcTrack->Pt();
 769       gPx = tpcTrack->Px();
 770       gPy = tpcTrack->Py();
 771       gPz = tpcTrack->Pz();
 772       eta = tpcTrack->Eta();
 773     }
 774   }
 775   else {
 776     gPt = track->Pt();
 777     gPx = track->Px();
 778     gPy = track->Py();
 779     gPz = track->Pz();
 780     eta = track->Eta();
 781   }
 782
 783   if((gPt < fMinPt) || (gPt > fMaxPt)) return kFALSE;
 784   if(fAnalysisEtaMode) {
 785     if((eta < fMinY) || (eta > fMaxY))
 786       return kFALSE;
 787   }
 788   else {
 789     if((Rapidity(gPx,gPy,gPz) < fMinY) || (Rapidity(gPx,gPy,gPz) > fMaxY))
 790       return kFALSE;
 791   }
 792
 793   return kTRUE;
 794 }
 795
 796 //____________________________________________________________________//
 797 Bool_t AliProtonAnalysis::IsAccepted(AliESDEvent *esd,
 798                                      const AliESDVertex *vertex,
 799                                      AliESDtrack* track) {
 800   // Checks if the track is excluded from the cuts
 801   Double_t gPt = 0.0, gPx = 0.0, gPy = 0.0, gPz = 0.0;
 802   Double_t dca[2] = {0.0,0.0}, cov[3] = {0.0,0.0,0.0};  //The impact parameters and their covariance.
 803
 804   if((fUseTPCOnly)&&(!fUseHybridTPC)) {
 805     AliExternalTrackParam *tpcTrack = (AliExternalTrackParam *)track->GetTPCInnerParam();
 806     if(!tpcTrack) {
 807       gPt = 0.0; gPx = 0.0; gPy = 0.0; gPz = 0.0;
 808       dca[0] = -100.; dca[1] = -100.;
 809       cov[0] = -100.; cov[1] = -100.; cov[2] = -100.;
 810     }
 811     else {
 812       gPt = tpcTrack->Pt();
 813       gPx = tpcTrack->Px();
 814       gPy = tpcTrack->Py();
 815       gPz = tpcTrack->Pz();
 816       tpcTrack->PropagateToDCA(vertex,
 817                                esd->GetMagneticField(),
 818                                100.,dca,cov);
 819     }
 820   }
 821   else if(fUseHybridTPC) {
 822      AliExternalTrackParam *tpcTrack = (AliExternalTrackParam *)track->GetTPCInnerParam();
 823     if(!tpcTrack) {
 824       gPt = 0.0; gPx = 0.0; gPy = 0.0; gPz = 0.0;
 825       dca[0] = -100.; dca[1] = -100.;
 826       cov[0] = -100.; cov[1] = -100.; cov[2] = -100.;
 827     }
 828     else {
 829       gPt = tpcTrack->Pt();
 830       gPx = tpcTrack->Px();
 831       gPy = tpcTrack->Py();
 832       gPz = tpcTrack->Pz();
 833       tpcTrack->PropagateToDCA(vertex,
 834                                esd->GetMagneticField(),
 835                                100.,dca,cov);
 836     }
 837   }
 838   else{
 839     gPt = track->Pt();
 840     gPx = track->Px();
 841     gPy = track->Py();
 842     gPz = track->Pz();
 843     track->PropagateToDCA(vertex,
 844                           esd->GetMagneticField(),
 845                           100.,dca,cov);
 846   }
 847
 848   Int_t  fIdxInt[200];
 849   Int_t nClustersITS = track->GetITSclusters(fIdxInt);
 850   Int_t nClustersTPC = track->GetTPCclusters(fIdxInt);
 851
 852   Float_t chi2PerClusterITS = -1;
 853   if (nClustersITS!=0)
 854     chi2PerClusterITS = track->GetITSchi2()/Float_t(nClustersITS);
 855   Float_t chi2PerClusterTPC = -1;
 856   if (nClustersTPC!=0)
 857     chi2PerClusterTPC = track->GetTPCchi2()/Float_t(nClustersTPC);
 858
 859   Double_t extCov[15];
 860   track->GetExternalCovariance(extCov);
 861
 862   if(fPointOnITSLayer1Flag)
 863     if(!track->HasPointOnITSLayer(0)) return kFALSE;
 864   if(fPointOnITSLayer2Flag)
 865     if(!track->HasPointOnITSLayer(1)) return kFALSE;
 866   if(fPointOnITSLayer3Flag)
 867     if(!track->HasPointOnITSLayer(2)) return kFALSE;
 868   if(fPointOnITSLayer4Flag)
 869     if(!track->HasPointOnITSLayer(3)) return kFALSE;
 870   if(fPointOnITSLayer5Flag)
 871     if(!track->HasPointOnITSLayer(4)) return kFALSE;
 872   if(fPointOnITSLayer6Flag)
 873     if(!track->HasPointOnITSLayer(5)) return kFALSE;
 874   if(fMinITSClustersFlag)
 875     if(nClustersITS < fMinITSClusters) return kFALSE;
 876   if(fMaxChi2PerITSClusterFlag)
 877     if(chi2PerClusterITS > fMaxChi2PerITSCluster) return kFALSE;
 878   if(fMinTPCClustersFlag)
 879     if(nClustersTPC < fMinTPCClusters) return kFALSE;
 880   if(fMaxChi2PerTPCClusterFlag)
 881     if(chi2PerClusterTPC > fMaxChi2PerTPCCluster) return kFALSE;
 882   if(fMaxCov11Flag)
 883     if(extCov[0] > fMaxCov11) return kFALSE;
 884   if(fMaxCov22Flag)
 885     if(extCov[2] > fMaxCov22) return kFALSE;
 886   if(fMaxCov33Flag)
 887     if(extCov[5] > fMaxCov33) return kFALSE;
 888   if(fMaxCov44Flag)
 889     if(extCov[9] > fMaxCov44) return kFALSE;
 890   if(fMaxCov55Flag)
 891     if(extCov[14] > fMaxCov55) return kFALSE;
 892   if(fMaxSigmaToVertexFlag)
 893     if(GetSigmaToVertex(track) > fMaxSigmaToVertex) return kFALSE;
 894   if(fMaxSigmaToVertexTPCFlag)
 895     if(GetSigmaToVertex(track) > fMaxSigmaToVertexTPC) return kFALSE;
 896   if(fMaxDCAXYFlag)
 897     if(TMath::Abs(dca[0]) > fMaxDCAXY) return kFALSE;
 898   if(fMaxDCAXYTPCFlag)
 899     if(TMath::Abs(dca[0]) > fMaxDCAXYTPC) return kFALSE;
 900     if(fMaxDCAZFlag)
 901     if(TMath::Abs(dca[1]) > fMaxDCAZ) return kFALSE;
 902   if(fMaxDCAZTPCFlag)
 903     if(TMath::Abs(dca[1]) > fMaxDCAZTPC) return kFALSE;
 904   if(fMaxConstrainChi2Flag) {
 905     if(track->GetConstrainedChi2() > 0)
 906       if(TMath::Log(track->GetConstrainedChi2()) > fMaxConstrainChi2) return kFALSE;
 907   }
 908   if(fITSRefitFlag)
 909     if ((track->GetStatus() & AliESDtrack::kITSrefit) == 0) return kFALSE;
 910   if(fTPCRefitFlag)
 911     if ((track->GetStatus() & AliESDtrack::kTPCrefit) == 0) return kFALSE;
 912   if(fESDpidFlag)
 913     if ((track->GetStatus() & AliESDtrack::kESDpid) == 0) return kFALSE;
 914   if(fTPCpidFlag)
 915     if ((track->GetStatus() & AliESDtrack::kTPCpid) == 0) return kFALSE;
 916
 917   return kTRUE;
 918 }
 919
 920 //____________________________________________________________________//
 921 Float_t AliProtonAnalysis::GetSigmaToVertex(AliESDtrack* esdTrack) const {
 922   // Calculates the number of sigma to the vertex.
 923
 924   Float_t b[2];
 925   Float_t bRes[2];
 926   Float_t bCov[3];
 927   if((fUseTPCOnly)&&(!fUseHybridTPC))
 928     esdTrack->GetImpactParametersTPC(b,bCov);
 929   else
 930     esdTrack->GetImpactParameters(b,bCov);
 931
 932   if (bCov[0]<=0 || bCov[2]<=0) {
 933     //AliDebug(1, "Estimated b resolution lower or equal zero!");
 934     bCov[0]=0; bCov[2]=0;
 935   }
 936   bRes[0] = TMath::Sqrt(bCov[0]);
 937   bRes[1] = TMath::Sqrt(bCov[2]);
 938
 939   if (bRes[0] == 0 || bRes[1] ==0) return -1;
 940
 941   Float_t d = TMath::Sqrt(TMath::Power(b[0]/bRes[0],2) + TMath::Power(b[1]/bRes[1],2));
 942
 943   if (TMath::Exp(-d * d / 2) < 1e-10) return 1000;
 944
 945   d = TMath::ErfInverse(1 - TMath::Exp(-d * d / 2)) * TMath::Sqrt(2);
 946
 947   return d;
 948 }
 949
 950 //____________________________________________________________________//
 951 Double_t AliProtonAnalysis::Rapidity(Double_t gPx, Double_t gPy, Double_t gPz) const {
 952   //returns the rapidity of the proton - to be removed
 953   Double_t fMass = 9.38270000000000048e-01;
 954
 955   Double_t gP = TMath::Sqrt(TMath::Power(gPx,2) +
 956                            TMath::Power(gPy,2) +
 957                            TMath::Power(gPz,2));
 958   Double_t energy = TMath::Sqrt(gP*gP + fMass*fMass);
 959   Double_t y = -999;
 960   if(energy != gPz)
 961     y = 0.5*TMath::Log((energy + gPz)/(energy - gPz));
 962
 963   return y;
 964 }
 965
 966 //____________________________________________________________________//
 967 Bool_t AliProtonAnalysis::PrintMean(TH1 *hist, Double_t edge) {
 968   //calculates the mean value of the ratio/asymmetry within \pm edge
 969   Double_t sum = 0.0;
 970   Int_t nentries = 0;
 971   //calculate the mean
 972   for(Int_t i = 0; i < hist->GetXaxis()->GetNbins(); i++) {
 973     Double_t x = hist->GetBinCenter(i+1);
 974     Double_t y = hist->GetBinContent(i+1);
 975     if(TMath::Abs(x) < edge) {
 976       sum += y;
 977       nentries += 1;
 978     }
 979   }
 980   Double_t mean = 0.0;
 981   if(nentries != 0)
 982     mean = sum/nentries;
 983
 984   //calculate the error
 985   for(Int_t i = 0; i < hist->GetXaxis()->GetNbins(); i++) {
 986     Double_t x = hist->GetBinCenter(i+1);
 987     Double_t y = hist->GetBinContent(i+1);
 988     if(TMath::Abs(x) < edge) {
 989       sum += TMath::Power((mean - y),2);
 990       nentries += 1;
 991     }
 992   }
 993
 994   Double_t error = 0.0;
 995   if(nentries != 0)
 996     error =  TMath::Sqrt(sum)/nentries;
 997
 998   cout<<"========================================="<<endl;
 999   cout<<"Input distribution: "<<hist->GetName()<<endl;
1000   cout<<"Interval used: -"<<edge<<" -> "<<edge<<endl;
1001   cout<<"Mean value :"<<mean<<endl;
1002   cout<<"Error: "<<error<<endl;
1003   cout<<"========================================="<<endl;
1004
1005   return 0;
1006 }
1007
1008 //____________________________________________________________________//
1009 Bool_t AliProtonAnalysis::PrintYields(TH1 *hist, Double_t edge) {
1010   //calculates the (anti)proton yields within the \pm edge
1011   Double_t sum = 0.0, sumerror = 0.0;
1012   Double_t error = 0.0;
1013   for(Int_t i = 0; i < hist->GetXaxis()->GetNbins(); i++) {
1014     Double_t x = hist->GetBinCenter(i+1);
1015     Double_t y = hist->GetBinContent(i+1);
1016     if(TMath::Abs(x) < edge) {
1017       sum += y;
1018       sumerror += TMath::Power(hist->GetBinError(i+1),2);
1019     }
1020   }
1021
1022   error = TMath::Sqrt(sumerror);
1023
1024   cout<<"========================================="<<endl;
1025   cout<<"Input distribution: "<<hist->GetName()<<endl;
1026   cout<<"Interval used: -"<<edge<<" -> "<<edge<<endl;
1027   cout<<"Yields :"<<sum<<endl;
1028   cout<<"Error: "<<error<<endl;
1029   cout<<"========================================="<<endl;
1030
1031   return 0;
1032 }
1033
1034 //____________________________________________________________________//
1035 void AliProtonAnalysis::Correct(Int_t step) {
1036   //Applies the correction maps to the initial containers
1037   fCorrectProtons = new AliCFDataGrid("correctProtons",
1038                                       "corrected data",
1039                                       *fProtonContainer);
1040   fCorrectProtons->ApplyEffCorrection(*(AliCFEffGrid *)fEffGridListProtons->At(step));
1041
1042   fCorrectAntiProtons = new AliCFDataGrid("correctAntiProtons",
1043                                           "corrected data",
1044                                           *fAntiProtonContainer);
1045   fCorrectAntiProtons->ApplyEffCorrection(*(AliCFEffGrid *)fEffGridListAntiProtons->At(step));
1046 }
1047
1048 //____________________________________________________________________//
1049 Bool_t AliProtonAnalysis::ReadCorrectionContainer(const char* filename) {
1050   // Reads the outout of the correction framework task
1051   // Creates the correction maps
1052   // Puts the results in the different TList objects
1053   Bool_t status = kTRUE;
1054
1055   TFile *file = TFile::Open(filename);
1056   if(!file) {
1057     cout<<"Could not find the input CORRFW file "<<filename<<endl;
1058     status = kFALSE;
1059   }
1060
1061   //________________________________________//
1062   //Protons
1063   fEffGridListProtons = new TList();
1064   fCorrectionListProtons2D = new TList();
1065   fEfficiencyListProtons1D = new TList();
1066   fCorrectionListProtons1D = new TList();
1067
1068   AliCFContainer *corrfwContainerProtons = (AliCFContainer*) (file->Get("containerProtons"));
1069   if(!corrfwContainerProtons) {
1070     cout<<"CORRFW container for protons not found!"<<endl;
1071     status = kFALSE;
1072   }
1073
1074   Int_t nSteps = corrfwContainerProtons->GetNStep();
1075   TH2D *gYPt[4];
1076   //currently the GRID is formed by the y-pT parameters
1077   //Add Vz as a next step
1078   Int_t iRap = 0, iPt = 1;
1079   AliCFEffGrid *effProtonsStep0Step1 = new AliCFEffGrid("eff10",
1080                                          "effProtonsStep0Step1",
1081                                          *corrfwContainerProtons);
1082   effProtonsStep0Step1->CalculateEfficiency(1,0); //eff= step1/step0
1083   fEffGridListProtons->Add(effProtonsStep0Step1);
1084   gYPt[0] = effProtonsStep0Step1->Project(iRap,iPt);
1085   fCorrectionListProtons2D->Add(gYPt[0]);
1086
1087   AliCFEffGrid *effProtonsStep0Step2 = new AliCFEffGrid("eff20",
1088                                          "effProtonsStep0Step2",
1089                                          *corrfwContainerProtons);
1090   effProtonsStep0Step2->CalculateEfficiency(2,0); //eff= step2/step0
1091   fEffGridListProtons->Add(effProtonsStep0Step2);
1092   gYPt[1] = effProtonsStep0Step2->Project(iRap,iPt);
1093   fCorrectionListProtons2D->Add(gYPt[1]);
1094
1095   AliCFEffGrid *effProtonsStep0Step3 = new AliCFEffGrid("eff30",
1096                                          "effProtonsStep0Step3",
1097                                          *corrfwContainerProtons);
1098   effProtonsStep0Step3->CalculateEfficiency(3,0); //eff= step1/step0
1099   fEffGridListProtons->Add(effProtonsStep0Step3);
1100   gYPt[2] = effProtonsStep0Step3->Project(iRap,iPt);
1101   fCorrectionListProtons2D->Add(gYPt[2]);
1102
1103   TH1D *gEfficiency[2][3]; //efficiency as a function of pT and of y (raws-[2])
1104   TH1D *gCorrection[2][3]; //efficiency as a function of pT and of y (raws-[2])
1105   TString gTitle = 0;
1106   //Get the projection of the efficiency maps
1107   for(Int_t iParameter = 0; iParameter < 2; iParameter++) {
1108     gEfficiency[iParameter][0] = effProtonsStep0Step1->Project(iParameter);
1109     gTitle = "ProtonsEfficiency_Parameter"; gTitle += iParameter+1;
1110     gTitle += "_Step0_Step1";
1111     gEfficiency[iParameter][0]->SetName(gTitle.Data());
1112     fEfficiencyListProtons1D->Add(gEfficiency[iParameter][0]);
1113     gTitle = "ProtonsCorrection_Parameter"; gTitle += iParameter+1;
1114     gTitle += "_Step0_Step1";
1115     gCorrection[iParameter][0] = new TH1D(gTitle.Data(),
1116                                           gTitle.Data(),
1117                                           gEfficiency[iParameter][0]->GetNbinsX(),
1118                                           gEfficiency[iParameter][0]->GetXaxis()->GetXmin(),
1119                                           gEfficiency[iParameter][0]->GetXaxis()->GetXmax());
1120     //initialisation of the correction
1121     for(Int_t iBin = 1; iBin <= gEfficiency[iParameter][0]->GetNbinsX(); iBin++)
1122       gCorrection[iParameter][0]->SetBinContent(iBin,1.0);
1123
1124     gEfficiency[iParameter][1] = effProtonsStep0Step2->Project(iParameter);
1125     gTitle = "ProtonsEfficiency_Parameter"; gTitle += iParameter+1;
1126     gTitle += "_Step0_Step2";
1127     gEfficiency[iParameter][1]->SetName(gTitle.Data());
1128     fEfficiencyListProtons1D->Add(gEfficiency[iParameter][1]);
1129     gTitle = "ProtonsCorrection_Parameter"; gTitle += iParameter+1;
1130     gTitle += "_Step0_Step2";
1131     gCorrection[iParameter][1] = new TH1D(gTitle.Data(),
1132                                           gTitle.Data(),
1133                                           gEfficiency[iParameter][1]->GetNbinsX(),
1134                                           gEfficiency[iParameter][1]->GetXaxis()->GetXmin(),
1135                                           gEfficiency[iParameter][1]->GetXaxis()->GetXmax());
1136     //initialisation of the correction
1137     for(Int_t iBin = 1; iBin <= gEfficiency[iParameter][1]->GetNbinsX(); iBin++)
1138       gCorrection[iParameter][1]->SetBinContent(iBin,1.0);
1139
1140     gEfficiency[iParameter][2] = effProtonsStep0Step3->Project(iParameter);
1141     gTitle = "ProtonsEfficiency_Parameter"; gTitle += iParameter+1;
1142     gTitle += "_Step0_Step3";
1143     gEfficiency[iParameter][2]->SetName(gTitle.Data());
1144     fEfficiencyListProtons1D->Add(gEfficiency[iParameter][2]);
1145     gTitle = "ProtonsCorrection_Parameter"; gTitle += iParameter+1;
1146     gTitle += "_Step0_Step3";
1147     gCorrection[iParameter][2] = new TH1D(gTitle.Data(),
1148                                           gTitle.Data(),
1149                                           gEfficiency[iParameter][2]->GetNbinsX(),
1150                                           gEfficiency[iParameter][2]->GetXaxis()->GetXmin(),
1151                                           gEfficiency[iParameter][2]->GetXaxis()->GetXmax());
1152     //initialisation of the correction
1153     for(Int_t iBin = 1; iBin <= gEfficiency[iParameter][2]->GetNbinsX(); iBin++)
1154       gCorrection[iParameter][2]->SetBinContent(iBin,1.0);
1155   }//parameter loop
1156   //Calculate the 1D correction parameters as a function of y and pT
1157   for(Int_t iParameter = 0; iParameter < 2; iParameter++) {
1158     for(Int_t iStep = 1; iStep < nSteps; iStep++) {
1159       gCorrection[iParameter][iStep-1]->Divide(gEfficiency[iParameter][iStep-1]);
1160       fCorrectionListProtons1D->Add(gCorrection[iParameter][iStep-1]);
1161     }
1162   }
1163
1164   //________________________________________//
1165   //AntiProtons
1166   fEffGridListAntiProtons = new TList();
1167   fCorrectionListAntiProtons2D = new TList();
1168   fEfficiencyListAntiProtons1D = new TList();
1169   fCorrectionListAntiProtons1D = new TList();
1170
1171   AliCFContainer *corrfwContainerAntiProtons = (AliCFContainer*) (file->Get("containerAntiProtons"));
1172   if(!corrfwContainerAntiProtons) {
1173     cout<<"CORRFW container for antiprotons not found!"<<endl;
1174     status = kFALSE;
1175   }
1176
1177   nSteps = corrfwContainerAntiProtons->GetNStep();
1178   //currently the GRID is formed by the y-pT parameters
1179   //Add Vz as a next step
1180   AliCFEffGrid *effAntiProtonsStep0Step1 = new AliCFEffGrid("eff10",
1181                                          "effAntiProtonsStep0Step1",
1182                                          *corrfwContainerAntiProtons);
1183   effAntiProtonsStep0Step1->CalculateEfficiency(1,0); //eff= step1/step0
1184   fEffGridListAntiProtons->Add(effAntiProtonsStep0Step1);
1185   gYPt[0] = effAntiProtonsStep0Step1->Project(iRap,iPt);
1186   fCorrectionListAntiProtons2D->Add(gYPt[0]);
1187
1188   AliCFEffGrid *effAntiProtonsStep0Step2 = new AliCFEffGrid("eff20",
1189                                          "effAntiProtonsStep0Step2",
1190                                          *corrfwContainerAntiProtons);
1191   effAntiProtonsStep0Step2->CalculateEfficiency(2,0); //eff= step2/step0
1192   fEffGridListAntiProtons->Add(effAntiProtonsStep0Step2);
1193   gYPt[1] = effAntiProtonsStep0Step2->Project(iRap,iPt);
1194   fCorrectionListAntiProtons2D->Add(gYPt[1]);
1195
1196   AliCFEffGrid *effAntiProtonsStep0Step3 = new AliCFEffGrid("eff30",
1197                                          "effAntiProtonsStep0Step3",
1198                                          *corrfwContainerAntiProtons);
1199   effAntiProtonsStep0Step3->CalculateEfficiency(3,0); //eff= step1/step0
1200   fEffGridListAntiProtons->Add(effAntiProtonsStep0Step3);
1201   gYPt[2] = effAntiProtonsStep0Step3->Project(iRap,iPt);
1202   fCorrectionListAntiProtons2D->Add(gYPt[2]);
1203
1204   //Get the projection of the efficiency maps
1205   for(Int_t iParameter = 0; iParameter < 2; iParameter++) {
1206     gEfficiency[iParameter][0] = effAntiProtonsStep0Step1->Project(iParameter);
1207     gTitle = "AntiProtonsEfficiency_Parameter"; gTitle += iParameter+1;
1208     gTitle += "_Step0_Step1";
1209     gEfficiency[iParameter][0]->SetName(gTitle.Data());
1210     fEfficiencyListAntiProtons1D->Add(gEfficiency[iParameter][0]);
1211     gTitle = "AntiProtonsCorrection_Parameter"; gTitle += iParameter+1;
1212     gTitle += "_Step0_Step1";
1213     gCorrection[iParameter][0] = new TH1D(gTitle.Data(),
1214                                           gTitle.Data(),
1215                                           gEfficiency[iParameter][0]->GetNbinsX(),
1216                                           gEfficiency[iParameter][0]->GetXaxis()->GetXmin(),
1217                                           gEfficiency[iParameter][0]->GetXaxis()->GetXmax());
1218     //initialisation of the correction
1219     for(Int_t iBin = 1; iBin <= gEfficiency[iParameter][0]->GetNbinsX(); iBin++)
1220       gCorrection[iParameter][0]->SetBinContent(iBin,1.0);
1221
1222     gEfficiency[iParameter][1] = effAntiProtonsStep0Step2->Project(iParameter);
1223     gTitle = "AntiProtonsEfficiency_Parameter"; gTitle += iParameter+1;
1224     gTitle += "_Step0_Step2";
1225     gEfficiency[iParameter][1]->SetName(gTitle.Data());
1226     fEfficiencyListAntiProtons1D->Add(gEfficiency[iParameter][1]);
1227     gTitle = "AntiProtonsCorrection_Parameter"; gTitle += iParameter+1;
1228     gTitle += "_Step0_Step2";
1229     gCorrection[iParameter][1] = new TH1D(gTitle.Data(),
1230                                           gTitle.Data(),
1231                                           gEfficiency[iParameter][1]->GetNbinsX(),
1232                                           gEfficiency[iParameter][1]->GetXaxis()->GetXmin(),
1233                                           gEfficiency[iParameter][1]->GetXaxis()->GetXmax());
1234     //initialisation of the correction
1235     for(Int_t iBin = 1; iBin <= gEfficiency[iParameter][1]->GetNbinsX(); iBin++)
1236       gCorrection[iParameter][1]->SetBinContent(iBin,1.0);
1237
1238     gEfficiency[iParameter][2] = effAntiProtonsStep0Step3->Project(iParameter);
1239     gTitle = "AntiProtonsEfficiency_Parameter"; gTitle += iParameter+1;
1240     gTitle += "_Step0_Step3";
1241     gEfficiency[iParameter][2]->SetName(gTitle.Data());
1242     fEfficiencyListAntiProtons1D->Add(gEfficiency[iParameter][2]);
1243     gTitle = "AntiProtonsCorrection_Parameter"; gTitle += iParameter+1;
1244     gTitle += "_Step0_Step3";
1245     gCorrection[iParameter][2] = new TH1D(gTitle.Data(),
1246                                           gTitle.Data(),
1247                                           gEfficiency[iParameter][2]->GetNbinsX(),
1248                                           gEfficiency[iParameter][2]->GetXaxis()->GetXmin(),
1249                                           gEfficiency[iParameter][2]->GetXaxis()->GetXmax());
1250     //initialisation of the correction
1251     for(Int_t iBin = 1; iBin <= gEfficiency[iParameter][2]->GetNbinsX(); iBin++)
1252       gCorrection[iParameter][2]->SetBinContent(iBin,1.0);
1253   }//parameter loop
1254   //Calculate the 1D correction parameters as a function of y and pT
1255   for(Int_t iParameter = 0; iParameter < 2; iParameter++) {
1256     for(Int_t iStep = 1; iStep < nSteps; iStep++) {
1257       gCorrection[iParameter][iStep-1]->Divide(gEfficiency[iParameter][iStep-1]);
1258       fCorrectionListAntiProtons1D->Add(gCorrection[iParameter][iStep-1]);
1259     }
1260   }
1261
1262   return status;
1263 }
1264
1265
1266
1267
1268
1269
1270
1271
1272