Adding the base class for the bar{p}-p analysis and adopting the new scheme
authorpchrist <pchrist@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 2 Apr 2009 09:38:50 +0000 (09:38 +0000)
committerpchrist <pchrist@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 2 Apr 2009 09:38:50 +0000 (09:38 +0000)
PWG2/AliAnalysisTaskProtons.cxx
PWG2/AliAnalysisTaskProtons.h
PWG2/SPECTRA/AliProtonAnalysis.cxx
PWG2/SPECTRA/AliProtonAnalysis.h
PWG2/SPECTRA/AliProtonAnalysisBase.cxx [new file with mode: 0644]
PWG2/SPECTRA/AliProtonAnalysisBase.h [new file with mode: 0644]
PWG2/configProtonAnalysis.C [new file with mode: 0644]
PWG2/runProtonAnalysis.C

index 703f9af..53a163c 100644 (file)
@@ -1,7 +1,9 @@
+#include "TStyle.h"
 #include "TChain.h"
 #include "TTree.h"
 #include "TString.h"
 #include "TList.h"
+#include "TFile.h"
 #include "TH2F.h"
 #include "TH1I.h"
 #include "TF1.h"
 #include "AliMCEventHandler.h"
 #include "AliMCEvent.h"
 #include "AliStack.h"
-#include "AliCFContainer.h"
 #include "AliVertexerTracks.h"
 #include "AliESDVertex.h"
 
 #include "AliProtonAnalysis.h"
+#include "AliProtonAnalysisBase.h"
 #include "AliAnalysisTaskProtons.h"
 
-// Analysis task creating a the 2d y-p_t spectrum of p and antip
+// Analysis task to run the \bar{p}/p analysis
 // Author: Panos Cristakoglou
 
 ClassImp(AliAnalysisTaskProtons)
-
+  
 //________________________________________________________________________ 
 AliAnalysisTaskProtons::AliAnalysisTaskProtons()
-  : AliAnalysisTask(), fESD(0), fAOD(0), fMC(0), fAnalysisType("ESD"),
-    fList(0), fProtonAnalysis(0),
-    fElectronFunction(0), fMuonFunction(0),
-    fPionFunction(0), fKaonFunction(0), fProtonFunction(0),
-    fFunctionUsed(kFALSE) {
+  : AliAnalysisTask(), fESD(0), fAOD(0), fMC(0),
+    fList(0), fProtonAnalysis(0) {
   //Dummy constructor
-                                                                                                   
+  
 }
 
 //________________________________________________________________________
 AliAnalysisTaskProtons::AliAnalysisTaskProtons(const char *name) 
-: AliAnalysisTask(name, ""), fESD(0), fAOD(0), fMC(0), fAnalysisType("ESD"),
-  fList(0), fProtonAnalysis(0), 
-  fElectronFunction(0), fMuonFunction(0),
-  fPionFunction(0), fKaonFunction(0), fProtonFunction(0),
-  fFunctionUsed(kFALSE),
-  fTriggerMode(kMB2), fProtonAnalysisMode(kTPC),
-  fVxMax(0), fVyMax(0), fVzMax(0) { 
+  : AliAnalysisTask(name, ""), fESD(0), fAOD(0), fMC(0),
+    fList(0), fProtonAnalysis(0) {
   // Constructor
-
+  
   // Define input and output slots here
   // Input slot #0 works with a TChain
   DefineInput(0, TChain::Class());
@@ -63,12 +57,13 @@ AliAnalysisTaskProtons::AliAnalysisTaskProtons(const char *name)
 void AliAnalysisTaskProtons::ConnectInputData(Option_t *) {
   // Connect ESD or AOD here
   // Called once
+  TString gAnalysisLevel = dynamic_cast<AliProtonAnalysisBase*>(fProtonAnalysis->GetProtonAnalysisBaseObject())->GetAnalysisLevel(); 
 
   TTree* tree = dynamic_cast<TTree*> (GetInputData(0));
   if (!tree) {
     Printf("ERROR: Could not read chain from input slot 0");
   } else {
-    if(fAnalysisType == "ESD") {
+    if(gAnalysisLevel == "ESD") {
       AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
      
       if (!esdH) {
@@ -76,7 +71,7 @@ void AliAnalysisTaskProtons::ConnectInputData(Option_t *) {
       } else
        fESD = esdH->GetEvent();
     }
-    else if(fAnalysisType == "AOD") {
+    else if(gAnalysisLevel == "AOD") {
      AliAODInputHandler *aodH = dynamic_cast<AliAODInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
      
      if (!aodH) {
@@ -84,7 +79,7 @@ void AliAnalysisTaskProtons::ConnectInputData(Option_t *) {
      } else
        fAOD = aodH->GetEvent();
     }
-    else if(fAnalysisType == "MC") {
+    else if(gAnalysisLevel == "MC") {
      AliMCEventHandler* mcH = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
      if (!mcH) {
        Printf("ERROR: Could not retrieve MC event handler");
@@ -99,91 +94,8 @@ void AliAnalysisTaskProtons::ConnectInputData(Option_t *) {
 
 //________________________________________________________________________
 void AliAnalysisTaskProtons::CreateOutputObjects() {
-  // Create histograms
+  // Create output objects
   // Called once
-  //Prior probabilities
-  Double_t partFrac[5] = {0.01, 0.01, 0.85, 0.10, 0.05};
-  
-  //proton analysis object
-  fProtonAnalysis = new AliProtonAnalysis();
-
-  if(fAnalysisType == "ESD") {
-    //Use of TPConly tracks
-    if(fProtonAnalysisMode == kTPC) {
-      fProtonAnalysis->InitAnalysisHistograms(10, -0.5, 0.5, 16, 0.5, 0.9);
-      fProtonAnalysis->UseTPCOnly();
-      //fProtonAnalysis->SetTPCpid();
-      fProtonAnalysis->SetMinTPCClusters(100);
-      fProtonAnalysis->SetMaxChi2PerTPCCluster(2.2);
-      fProtonAnalysis->SetMaxCov11(0.5);
-      fProtonAnalysis->SetMaxCov22(0.5);
-      fProtonAnalysis->SetMaxCov33(0.5);
-      fProtonAnalysis->SetMaxCov44(0.5);
-      fProtonAnalysis->SetMaxCov55(0.5);
-      fProtonAnalysis->SetMaxSigmaToVertexTPC(2.0);
-      //fProtonAnalysis->SetMaxDCAXYTPC(1.5);
-      //fProtonAnalysis->SetMaxDCAZTPC(1.5);
-    }
-    //Use of HybridTPC tracks
-    else if(fProtonAnalysisMode == kHybrid) {
-      fProtonAnalysis->InitAnalysisHistograms(10, -0.5, 0.5, 16, 0.5, 0.9);
-      fProtonAnalysis->UseHybridTPC();
-      fProtonAnalysis->SetTPCpid();
-      fProtonAnalysis->SetMinTPCClusters(110);
-      fProtonAnalysis->SetMaxChi2PerTPCCluster(2.2);
-      fProtonAnalysis->SetMaxCov11(0.5);
-      fProtonAnalysis->SetMaxCov22(0.5);
-      fProtonAnalysis->SetMaxCov33(0.5);
-      fProtonAnalysis->SetMaxCov44(0.5);
-      fProtonAnalysis->SetMaxCov55(0.5);
-      fProtonAnalysis->SetMaxSigmaToVertex(2.0);
-      /*fProtonAnalysis->SetMaxDCAXY(1.5);
-       fProtonAnalysis->SetMaxDCAZ(1.5);*/
-      fProtonAnalysis->SetPointOnITSLayer6();
-      fProtonAnalysis->SetPointOnITSLayer5();
-      //fProtonAnalysis->SetPointOnITSLayer4();
-      //fProtonAnalysis->SetPointOnITSLayer3();
-      fProtonAnalysis->SetPointOnITSLayer2();
-      fProtonAnalysis->SetPointOnITSLayer1();
-      fProtonAnalysis->SetMinITSClusters(5);
-    }
-    //Combined tracking
-    else if(fProtonAnalysisMode == kGlobal) {
-      fProtonAnalysis->InitAnalysisHistograms(20, -1.0, 1.0, 48, 0.3, 1.5);
-      fProtonAnalysis->SetMinTPCClusters(110);
-      fProtonAnalysis->SetMaxChi2PerTPCCluster(2.2);
-      fProtonAnalysis->SetMaxCov11(0.5);
-      fProtonAnalysis->SetMaxCov22(0.5);
-      fProtonAnalysis->SetMaxCov33(0.5);
-      fProtonAnalysis->SetMaxCov44(0.5);
-      fProtonAnalysis->SetMaxCov55(0.5);
-      fProtonAnalysis->SetMaxSigmaToVertex(2.0);
-      //fProtonAnalysis->SetMaxDCAXY(2.0);
-      //fProtonAnalysis->SetMaxDCAZ(2.0);
-      fProtonAnalysis->SetTPCRefit();
-      fProtonAnalysis->SetPointOnITSLayer1();
-      fProtonAnalysis->SetPointOnITSLayer2();
-      //fProtonAnalysis->SetPointOnITSLayer3();
-      //fProtonAnalysis->SetPointOnITSLayer4();
-      fProtonAnalysis->SetPointOnITSLayer5();
-      fProtonAnalysis->SetPointOnITSLayer6();
-      fProtonAnalysis->SetMinITSClusters(5);
-      fProtonAnalysis->SetITSRefit();
-      fProtonAnalysis->SetESDpid();
-    }
-
-    if(fFunctionUsed)
-      fProtonAnalysis->SetPriorProbabilityFunctions(fElectronFunction,
-                                             fMuonFunction,
-                                             fPionFunction,
-                                             fKaonFunction,
-                                             fProtonFunction);
-    else
-      fProtonAnalysis->SetPriorProbabilities(partFrac);
-  }//ESD analysis
-  else if(fAnalysisType == "MC") 
-    fProtonAnalysis->InitAnalysisHistograms(10, -0.5, 0.5, 16, 0.5, 0.9);
-  
   fList = new TList();
   fList->Add(fProtonAnalysis->GetProtonYPtHistogram());
   fList->Add(fProtonAnalysis->GetAntiProtonYPtHistogram());
@@ -196,16 +108,16 @@ void AliAnalysisTaskProtons::CreateOutputObjects() {
 void AliAnalysisTaskProtons::Exec(Option_t *) {
   // Main loop
   // Called for each event
+  TString gAnalysisLevel = dynamic_cast<AliProtonAnalysisBase*>(fProtonAnalysis->GetProtonAnalysisBaseObject())->GetAnalysisLevel(); 
   
-  if(fAnalysisType == "ESD") {
+  if(gAnalysisLevel == "ESD") {
     if (!fESD) {
       Printf("ERROR: fESD not available");
       return;
     }
 
-    if(IsEventTriggered(fESD,fTriggerMode)) {
-      const AliESDVertex *vertex = GetVertex(fESD,fProtonAnalysisMode,
-                                            fVxMax,fVyMax,fVzMax);
+    if(dynamic_cast<AliProtonAnalysisBase*>(fProtonAnalysis->GetProtonAnalysisBaseObject())->IsEventTriggered(fESD,dynamic_cast<AliProtonAnalysisBase*>(fProtonAnalysis->GetProtonAnalysisBaseObject())->GetTriggerMode())) {
+      const AliESDVertex *vertex = dynamic_cast<AliProtonAnalysisBase*>(fProtonAnalysis->GetProtonAnalysisBaseObject())->GetVertex(fESD,dynamic_cast<AliProtonAnalysisBase*>(fProtonAnalysis->GetProtonAnalysisBaseObject())->GetAnalysisMode(),dynamic_cast<AliProtonAnalysisBase*>(fProtonAnalysis->GetProtonAnalysisBaseObject())->GetVxMax(),dynamic_cast<AliProtonAnalysisBase*>(fProtonAnalysis->GetProtonAnalysisBaseObject())->GetVyMax(),dynamic_cast<AliProtonAnalysisBase*>(fProtonAnalysis->GetProtonAnalysisBaseObject())->GetVzMax());
       if(vertex) {
        Printf("Proton ESD analysis task: There are %d tracks in this event", fESD->GetNumberOfTracks());
        fProtonAnalysis->Analyze(fESD,vertex);
@@ -213,7 +125,7 @@ void AliAnalysisTaskProtons::Exec(Option_t *) {
     }//triggered event
   }//ESD analysis              
   
-  else if(fAnalysisType == "AOD") {
+  else if(gAnalysisLevel == "AOD") {
     if (!fAOD) {
       Printf("ERROR: fAOD not available");
       return;
@@ -221,9 +133,9 @@ void AliAnalysisTaskProtons::Exec(Option_t *) {
     
     Printf("Proton AOD analysis task: There are %d tracks in this event", fAOD->GetNumberOfTracks());
     fProtonAnalysis->Analyze(fAOD);
-  }//AOD analysis                                                                                                                                                                
+  }//AOD analysis
   
-  else if(fAnalysisType == "MC") {
+  else if(gAnalysisLevel == "MC") {
     if (!fMC) {
       Printf("ERROR: Could not retrieve MC event");
       return;
@@ -246,7 +158,8 @@ void AliAnalysisTaskProtons::Exec(Option_t *) {
 void AliAnalysisTaskProtons::Terminate(Option_t *) {
   // Draw result to the screen
   // Called once at the end of the query
-  
+  gStyle->SetPalette(1,0);
+
   fList = dynamic_cast<TList*> (GetOutputData(0));
   if (!fList) {
     Printf("ERROR: fList not available");
@@ -263,85 +176,10 @@ void AliAnalysisTaskProtons::Terminate(Option_t *) {
   if (fHistYPtProtons) fHistYPtProtons->DrawCopy("colz");
   c1->cd(2)->SetLeftMargin(0.15); c1->cd(2)->SetBottomMargin(0.15);  
   if (fHistYPtAntiProtons) fHistYPtAntiProtons->DrawCopy("colz");
-}
-
-//________________________________________________________________________
-Bool_t AliAnalysisTaskProtons::IsEventTriggered(const AliESDEvent *esd, 
-                                               TriggerMode trigger) {
-  // check if the event was triggered
-  ULong64_t triggerMask = esd->GetTriggerMask();
-  
-  // definitions from p-p.cfg
-  ULong64_t spdFO = (1 << 14);
-  ULong64_t v0left = (1 << 11);
-  ULong64_t v0right = (1 << 12);
 
-  switch (trigger) {
-  case kMB1: {
-    if (triggerMask & spdFO || ((triggerMask & v0left) || (triggerMask & v0right)))
-      return kTRUE;
-    break;
-  }
-  case kMB2: {
-    if (triggerMask & spdFO && ((triggerMask & v0left) || (triggerMask & v0right)))
-      return kTRUE;
-    break;
-  }
-  case kSPDFASTOR: {
-    if (triggerMask & spdFO)
-      return kTRUE;
-    break;
-  }
-  }//switch
-
-  return kFALSE;
+  TCanvas *c2 = dynamic_cast<AliProtonAnalysisBase*>(fProtonAnalysis->GetProtonAnalysisBaseObject())->GetListOfCuts();
+  TFile *flocal = TFile::Open("ListOfCuts.root","recreate");
+  c2->Write();
+  flocal->Close();
 }
 
-//________________________________________________________________________
-const AliESDVertex* AliAnalysisTaskProtons::GetVertex(AliESDEvent* esd, 
-                                                     AnalysisMode mode,
-                                                     Double_t gVxMax,
-                                                     Double_t gVyMax,
-                                                     Double_t gVzMax) {
-  // Get the vertex from the ESD and returns it if the vertex is valid
-  // Second argument decides which vertex is used (this selects
-  // also the quality criteria that are applied)
-  const AliESDVertex* vertex = 0;
-  if(mode == kHybrid)
-    vertex = esd->GetPrimaryVertexSPD();
-  else if(mode == kTPC){
-    Double_t kBz = esd->GetMagneticField();
-    AliVertexerTracks vertexer(kBz);
-    vertexer.SetTPCMode();
-    AliESDVertex *vTPC = vertexer.FindPrimaryVertex(esd);
-    esd->SetPrimaryVertexTPC(vTPC);
-    for (Int_t i=0; i<esd->GetNumberOfTracks(); i++) {
-      AliESDtrack *t = esd->GetTrack(i);
-      t->RelateToVertexTPC(vTPC, kBz, kVeryBig);
-    }
-    delete vTPC;
-    vertex = esd->GetPrimaryVertexTPC();
-  }
-  else if(mode == kGlobal)
-    vertex = esd->GetPrimaryVertex();
-  else
-    Printf("GetVertex: ERROR: Invalid second argument %d", mode);
-  
-  if(!vertex) return 0;
-  
-  // check Ncontributors
-  if(vertex->GetNContributors() <= 0) return 0;
-  
-  // check resolution
-  Double_t zRes = vertex->GetZRes();
-  if(zRes == 0) return 0;
-  
-  //check position
-  if(TMath::Abs(vertex->GetXv()) > gVxMax) return 0;
-  if(TMath::Abs(vertex->GetYv()) > gVyMax) return 0;
-  if(TMath::Abs(vertex->GetZv()) > gVzMax) return 0;
-  
-  return vertex;
-}
-
-
index 5fe4c8c..473e371 100644 (file)
@@ -1,23 +1,19 @@
 #ifndef AliAnalysisTaskProtons_cxx
 #define AliAnalysisTaskProtons_cxx
 
-// Analysis task creating a the 2d y-p_t spectrum of p and antip
+// Analysis task to run the \bar{p}/p analysis
 // Author: Panos Cristakoglou
-class TString;
+//class TString;
 class TList;
 class AliESDEvent;
 class AliAODEvent;
 class AliMCEvent;
 class AliProtonAnalysis;
-class TF1;
 
 #include "AliAnalysisTask.h"
 
 class AliAnalysisTaskProtons : public AliAnalysisTask {
  public:
-  enum TriggerMode { kMB1 = 0, kMB2, kSPDFASTOR }; 
-  enum AnalysisMode { kInvalid = -1, kTPC = 0, kHybrid, kGlobal };
-  
   AliAnalysisTaskProtons();
   AliAnalysisTaskProtons(const char *name);
   virtual ~AliAnalysisTaskProtons() {}
@@ -27,56 +23,18 @@ class AliAnalysisTaskProtons : public AliAnalysisTask {
   virtual void   Exec(Option_t *option);
   virtual void   Terminate(Option_t *);
 
-  void SetType(const char* type) {fAnalysisType = type;}
-  void SetPriorProbabilityFunctions(TF1 *felectrons,
-                                    TF1 *fmuons,
-                                    TF1 *fpions,
-                                    TF1 *fkaons,
-                                    TF1 *fprotons) {
-    fFunctionUsed = kTRUE;
-    fElectronFunction = felectrons;
-    fMuonFunction = fmuons;
-    fPionFunction = fpions;
-    fKaonFunction = fkaons;
-    fProtonFunction = fprotons;
-  }
-
-  void SetTriggerMode(TriggerMode triggermode) {fTriggerMode = triggermode;}
-  void SetAnalysisMode(AnalysisMode analysismode) {fProtonAnalysisMode = analysismode;}
-  void SetAcceptedVertexDiamond(Double_t gVx, Double_t gVy, Double_t gVz) {
-    fVxMax = gVx; fVyMax = gVy; fVzMax = gVz;}
-
-  static Bool_t IsEventTriggered(const AliESDEvent *esd, 
-                                TriggerMode trigger = kMB2);
-  static const  AliESDVertex *GetVertex(AliESDEvent *esd, 
-                                       AnalysisMode mode = kTPC,
-                                       Double_t gVx = 100.,
-                                       Double_t gVy = 100.,
-                                       Double_t gVz = 100.);
+  void SetAnalysisObject(AliProtonAnalysis *analysis) {
+    fProtonAnalysis = analysis;}
   
  private:
   AliESDEvent *fESD;    //ESD object 
   AliAODEvent *fAOD;    //AOD object
   AliMCEvent  *fMC;     //MC object 
   
-  TString fAnalysisType;//"ESD", "AOD" or "MC"
-
   TList  *fList; //TList output object 
   
   AliProtonAnalysis *fProtonAnalysis; //analysis object 
   
-  TF1 *fElectronFunction; //TF1 for e
-  TF1 *fMuonFunction; //TF1 for mu
-  TF1 *fPionFunction; //TF1 for pi
-  TF1 *fKaonFunction; //TF1 for K
-  TF1 *fProtonFunction; //TF1 for p
-  
-  Bool_t fFunctionUsed; //kTRUE if Functions are used
-  
-  TriggerMode   fTriggerMode; //Trigger mode
-  AnalysisMode  fProtonAnalysisMode; //Analysis mode
-  Double_t      fVxMax, fVyMax, fVzMax; //vertex diamond constrain
-
   AliAnalysisTaskProtons(const AliAnalysisTaskProtons&); // not implemented
   AliAnalysisTaskProtons& operator=(const AliAnalysisTaskProtons&); // not implemented
   
index 523a12d..f897823 100644 (file)
@@ -28,6 +28,7 @@
 #include <TParticle.h>
 
 #include "AliProtonAnalysis.h"
+#include "AliProtonAnalysisBase.h"
 
 #include <AliExternalTrackParam.h>
 #include <AliAODEvent.h>
@@ -44,33 +45,9 @@ ClassImp(AliProtonAnalysis)
 
 //____________________________________________________________________//
 AliProtonAnalysis::AliProtonAnalysis() : 
-  TObject(), fAnalysisEtaMode(kFALSE),
+  TObject(), fProtonAnalysisBase(0),
   fNBinsY(0), fMinY(0), fMaxY(0),
   fNBinsPt(0), fMinPt(0), fMaxPt(0),
-  fMinTPCClusters(0), fMinITSClusters(0),
-  fMaxChi2PerTPCCluster(0), fMaxChi2PerITSCluster(0),
-  fMaxCov11(0), fMaxCov22(0), fMaxCov33(0), fMaxCov44(0), fMaxCov55(0),
-  fMaxSigmaToVertex(0), fMaxSigmaToVertexTPC(0),
-  fMaxDCAXY(0), fMaxDCAXYTPC(0),
-  fMaxDCAZ(0), fMaxDCAZTPC(0),
-  fMaxConstrainChi2(0),
-  fMinTPCClustersFlag(kFALSE), fMinITSClustersFlag(kFALSE),
-  fMaxChi2PerTPCClusterFlag(kFALSE), fMaxChi2PerITSClusterFlag(kFALSE),
-  fMaxCov11Flag(kFALSE), fMaxCov22Flag(kFALSE), 
-  fMaxCov33Flag(kFALSE), fMaxCov44Flag(kFALSE), fMaxCov55Flag(kFALSE),
-  fMaxSigmaToVertexFlag(kFALSE), fMaxSigmaToVertexTPCFlag(kFALSE),
-  fMaxDCAXYFlag(kFALSE), fMaxDCAXYTPCFlag(kFALSE),
-  fMaxDCAZFlag(kFALSE), fMaxDCAZTPCFlag(kFALSE),
-  fMaxConstrainChi2Flag(kFALSE),
-  fITSRefitFlag(kFALSE), fTPCRefitFlag(kFALSE),
-  fESDpidFlag(kFALSE), fTPCpidFlag(kFALSE),
-  fPointOnITSLayer1Flag(0), fPointOnITSLayer2Flag(0),
-  fPointOnITSLayer3Flag(0), fPointOnITSLayer4Flag(0),
-  fPointOnITSLayer5Flag(0), fPointOnITSLayer6Flag(0),
-  fFunctionProbabilityFlag(kFALSE), 
-  fElectronFunction(0), fMuonFunction(0),
-  fPionFunction(0), fKaonFunction(0), fProtonFunction(0),
-  fUseTPCOnly(kFALSE), fUseHybridTPC(kFALSE), 
   fProtonContainer(0), fAntiProtonContainer(0),
   fHistEvents(0), fHistYPtProtons(0), fHistYPtAntiProtons(0),
   fEffGridListProtons(0), fCorrectionListProtons2D(0), 
@@ -79,38 +56,16 @@ AliProtonAnalysis::AliProtonAnalysis() :
   fEfficiencyListAntiProtons1D(0), fCorrectionListAntiProtons1D(0),
   fCorrectProtons(0), fCorrectAntiProtons(0) {
   //Default constructor
-  for(Int_t i = 0; i < 5; i++) fPartFrac[i] = 0.0;
 }
 
 //____________________________________________________________________//
-AliProtonAnalysis::AliProtonAnalysis(Int_t nbinsY, Float_t fLowY, Float_t fHighY,Int_t nbinsPt, Float_t fLowPt, Float_t fHighPt) : 
-  TObject(), fAnalysisEtaMode(kFALSE),
+AliProtonAnalysis::AliProtonAnalysis(Int_t nbinsY, 
+                                    Float_t fLowY, Float_t fHighY,
+                                    Int_t nbinsPt, 
+                                    Float_t fLowPt, Float_t fHighPt) : 
+  TObject(), fProtonAnalysisBase(0),
   fNBinsY(nbinsY), fMinY(fLowY), fMaxY(fHighY),
   fNBinsPt(nbinsPt), fMinPt(fLowPt), fMaxPt(fHighPt),
-  fMinTPCClusters(0), fMinITSClusters(0),
-  fMaxChi2PerTPCCluster(0), fMaxChi2PerITSCluster(0),
-  fMaxCov11(0), fMaxCov22(0), fMaxCov33(0), fMaxCov44(0), fMaxCov55(0),
-  fMaxSigmaToVertex(0), fMaxSigmaToVertexTPC(0),
-  fMaxDCAXY(0), fMaxDCAXYTPC(0),
-  fMaxDCAZ(0), fMaxDCAZTPC(0),
-  fMaxConstrainChi2(0),
-  fMinTPCClustersFlag(kFALSE), fMinITSClustersFlag(kFALSE),
-  fMaxChi2PerTPCClusterFlag(kFALSE), fMaxChi2PerITSClusterFlag(kFALSE),
-  fMaxCov11Flag(kFALSE), fMaxCov22Flag(kFALSE), 
-  fMaxCov33Flag(kFALSE), fMaxCov44Flag(kFALSE), fMaxCov55Flag(kFALSE),
-  fMaxSigmaToVertexFlag(kFALSE), fMaxSigmaToVertexTPCFlag(kFALSE),
-  fMaxDCAXYFlag(kFALSE), fMaxDCAXYTPCFlag(kFALSE),
-  fMaxDCAZFlag(kFALSE), fMaxDCAZTPCFlag(kFALSE),
-  fMaxConstrainChi2Flag(kFALSE),
-  fITSRefitFlag(kFALSE), fTPCRefitFlag(kFALSE),
-  fESDpidFlag(kFALSE), fTPCpidFlag(kFALSE),
-  fPointOnITSLayer1Flag(0), fPointOnITSLayer2Flag(0),
-  fPointOnITSLayer3Flag(0), fPointOnITSLayer4Flag(0),
-  fPointOnITSLayer5Flag(0), fPointOnITSLayer6Flag(0),
-  fFunctionProbabilityFlag(kFALSE), 
-  fElectronFunction(0), fMuonFunction(0),
-  fPionFunction(0), fKaonFunction(0), fProtonFunction(0),
-  fUseTPCOnly(kFALSE), fUseHybridTPC(kFALSE),   
   fProtonContainer(0), fAntiProtonContainer(0),
   fHistEvents(0), fHistYPtProtons(0), fHistYPtAntiProtons(0),
   fEffGridListProtons(0), fCorrectionListProtons2D(0), 
@@ -121,18 +76,26 @@ AliProtonAnalysis::AliProtonAnalysis(Int_t nbinsY, Float_t fLowY, Float_t fHighY
   //Default constructor
   fHistEvents = new TH1I("fHistEvents","Analyzed events",1,0,1);
 
-  fHistYPtProtons = new TH2D("fHistYPtProtons","y-Pt Protons",
-                            fNBinsY,fMinY,fMaxY,fNBinsPt,fMinPt,fMaxPt);
+  fHistYPtProtons = new TH2D("fHistYPtProtons","Protons",
+                            fNBinsY,fMinY,fMaxY,
+                            fNBinsPt,fMinPt,fMaxPt);
   fHistYPtProtons->SetStats(kTRUE);
-  fHistYPtProtons->GetYaxis()->SetTitle("P_{T} [GeV]");
-  fHistYPtProtons->GetXaxis()->SetTitle("y");
+  fHistYPtProtons->GetYaxis()->SetTitle("P_{T} [GeV/c]");
+  if(fProtonAnalysisBase->GetEtaMode())
+    fHistYPtProtons->GetXaxis()->SetTitle("#eta");
+  else
+    fHistYPtProtons->GetXaxis()->SetTitle("y");
   fHistYPtProtons->GetXaxis()->SetTitleColor(1);
 
-  fHistYPtAntiProtons = new TH2D("fHistYPtAntiProtons","y-Pt Antiprotons",
-                                fNBinsY,fMinY,fMaxY,fNBinsPt,fMinPt,fMaxPt);
+  fHistYPtAntiProtons = new TH2D("fHistYPtAntiProtons","Antiprotons",
+                                fNBinsY,fMinY,fMaxY,
+                                fNBinsPt,fMinPt,fMaxPt);
   fHistYPtAntiProtons->SetStats(kTRUE);
-  fHistYPtAntiProtons->GetYaxis()->SetTitle("P_{T} [GeV]");
-  fHistYPtAntiProtons->GetXaxis()->SetTitle("y");
+  fHistYPtAntiProtons->GetYaxis()->SetTitle("P_{T} [GeV/c]");
+  if(fProtonAnalysisBase->GetEtaMode())
+    fHistYPtAntiProtons->GetXaxis()->SetTitle("#eta");
+  else
+    fHistYPtAntiProtons->GetXaxis()->SetTitle("y");
   fHistYPtAntiProtons->GetXaxis()->SetTitleColor(1);
 
   //setting up the containers
@@ -150,18 +113,20 @@ AliProtonAnalysis::AliProtonAnalysis(Int_t nbinsY, Float_t fLowY, Float_t fHighY
   fProtonContainer = new AliCFContainer("containerProtons",
                                        "container for protons",
                                        1,2,iBin);
-  fProtonContainer->SetBinLimits(0,binLimY); //rapidity
+  fProtonContainer->SetBinLimits(0,binLimY); //rapidity or eta
   fProtonContainer->SetBinLimits(1,binLimPt); //pT
   fAntiProtonContainer = new AliCFContainer("containerAntiProtons",
                                            "container for antiprotons",
                                            1,2,iBin);
-  fAntiProtonContainer->SetBinLimits(0,binLimY); //rapidity
+  fAntiProtonContainer->SetBinLimits(0,binLimY); //rapidity or eta
   fAntiProtonContainer->SetBinLimits(1,binLimPt); //pT
 } 
 
 //____________________________________________________________________//
 AliProtonAnalysis::~AliProtonAnalysis() {
   //Default destructor
+  if(fProtonAnalysisBase) delete fProtonAnalysisBase;
+
   if(fHistEvents) delete fHistEvents;
   if(fHistYPtProtons) delete fHistYPtProtons;
   if(fHistYPtAntiProtons) delete fHistYPtAntiProtons;
@@ -193,20 +158,28 @@ void AliProtonAnalysis::InitAnalysisHistograms(Int_t nbinsY,
   fMinPt = fLowPt;
   fMaxPt = fHighPt;
 
-  fHistEvents = new TH1I("fHistEvents","Anallyzed events",1,0,1);
+  fHistEvents = new TH1I("fHistEvents","Analyzed events",1,0,1);
 
-  fHistYPtProtons = new TH2D("fHistYPtProtons","y-Pt Protons",
-                            fNBinsY,fMinY,fMaxY,fNBinsPt,fMinPt,fMaxPt);
+  fHistYPtProtons = new TH2D("fHistYPtProtons","Protons",
+                            fNBinsY,fMinY,fMaxY,
+                            fNBinsPt,fMinPt,fMaxPt);
   fHistYPtProtons->SetStats(kTRUE);
-  fHistYPtProtons->GetYaxis()->SetTitle("P_{T} [GeV]");
-  fHistYPtProtons->GetXaxis()->SetTitle("y");
+  fHistYPtProtons->GetYaxis()->SetTitle("P_{T} [GeV/c]");
+  if(fProtonAnalysisBase->GetEtaMode())
+    fHistYPtProtons->GetXaxis()->SetTitle("#eta");
+  else
+    fHistYPtProtons->GetXaxis()->SetTitle("y");
   fHistYPtProtons->GetXaxis()->SetTitleColor(1);
 
-  fHistYPtAntiProtons = new TH2D("fHistYPtAntiProtons","y-Pt Antiprotons",
-                                fNBinsY,fMinY,fMaxY,fNBinsPt,fMinPt,fMaxPt);
+  fHistYPtAntiProtons = new TH2D("fHistYPtAntiProtons","Antiprotons",
+                                fNBinsY,fMinY,fMaxY,
+                                fNBinsPt,fMinPt,fMaxPt);
   fHistYPtAntiProtons->SetStats(kTRUE);
-  fHistYPtAntiProtons->GetYaxis()->SetTitle("P_{T} [GeV]");
-  fHistYPtAntiProtons->GetXaxis()->SetTitle("y");
+  fHistYPtAntiProtons->GetYaxis()->SetTitle("P_{T} [GeV/c]");
+  if(fProtonAnalysisBase->GetEtaMode())
+    fHistYPtAntiProtons->GetXaxis()->SetTitle("#eta");
+  else
+    fHistYPtAntiProtons->GetXaxis()->SetTitle("y");
   fHistYPtAntiProtons->GetXaxis()->SetTitleColor(1);
 
   //setting up the containers
@@ -244,7 +217,7 @@ Bool_t AliProtonAnalysis::ReadFromFile(const char* filename) {
     status = kFALSE;
   }
 
-  TList *list = (TList *)file->Get("outputList1");
+  TList *list = (TList *)file->Get("outputList");
   if(list) {
     cout<<"Retrieving objects from the list "<<list->GetName()<<"..."<<endl; 
     fHistYPtProtons = (TH2D *)list->At(0);
@@ -561,36 +534,23 @@ TH1D *AliProtonAnalysis::GetPtAsymmetryHistogram() {
 }
 
 //____________________________________________________________________//
-Double_t AliProtonAnalysis::GetParticleFraction(Int_t i, Double_t p) {
-  //Return the a priori probs
-  Double_t partFrac=0;
-  if(fFunctionProbabilityFlag) {
-    if(i == 0) partFrac = fElectronFunction->Eval(p);
-    if(i == 1) partFrac = fMuonFunction->Eval(p);
-    if(i == 2) partFrac = fPionFunction->Eval(p);
-    if(i == 3) partFrac = fKaonFunction->Eval(p);
-    if(i == 4) partFrac = fProtonFunction->Eval(p);
-  }
-  else partFrac = fPartFrac[i];
-
-  return partFrac;
-}
-
-//____________________________________________________________________//
 void AliProtonAnalysis::Analyze(AliESDEvent* esd,
                                const AliESDVertex *vertex) {
   //Main analysis part - ESD
+  Int_t nTracks = 0;
+  Int_t nIdentifiedProtons = 0, nIdentifiedAntiProtons = 0;
+  Int_t nSurvivedProtons = 0, nSurvivedAntiProtons = 0;
+
   fHistEvents->Fill(0); //number of analyzed events
   Double_t containerInput[2] ;
   Double_t gPt = 0.0, gP = 0.0;
-  Int_t nGoodTracks = esd->GetNumberOfTracks();
-  for(Int_t iTracks = 0; iTracks < nGoodTracks; iTracks++) {
+  nTracks = esd->GetNumberOfTracks();
+  for(Int_t iTracks = 0; iTracks < nTracks; iTracks++) {
     AliESDtrack* track = esd->GetTrack(iTracks);
-    Double_t probability[5];
     AliESDtrack trackTPC;
 
     //in case it's a TPC only track relate it to the proper vertex
-    if((fUseTPCOnly)&&(!fUseHybridTPC)) {
+    /*if(fProtonAnalysisBase->GetAnalysisMode()==AliProtonAnalysisBase::kTPC) {
       Float_t p[2],cov[3];
       track->GetImpactParametersTPC(p,cov);
       if (p[0]==0 && p[1]==0)  
@@ -599,104 +559,128 @@ void AliProtonAnalysis::Analyze(AliESDEvent* esd,
        continue;
       }
       track = &trackTPC ;
-    }
+      }*/
 
-    if(IsAccepted(esd,vertex,track)) { 
-      if(fUseTPCOnly) {
-       AliExternalTrackParam *tpcTrack = (AliExternalTrackParam *)track->GetTPCInnerParam();
-       if(!tpcTrack) continue;
-       gPt = tpcTrack->Pt();
-       gP = tpcTrack->P();
-       
-       //pid
-       track->GetTPCpid(probability);
-       Double_t rcc = 0.0;
-       for(Int_t i = 0; i < AliPID::kSPECIES; i++) 
-         rcc += probability[i]*GetParticleFraction(i,gP);
-       if(rcc == 0.0) continue;
-       Double_t w[5];
-       for(Int_t i = 0; i < AliPID::kSPECIES; i++) 
-         w[i] = probability[i]*GetParticleFraction(i,gP)/rcc;
-       Long64_t fParticleType = TMath::LocMax(AliPID::kSPECIES,w);
-       if(fParticleType == 4) {
+    if((fProtonAnalysisBase->GetAnalysisMode()==AliProtonAnalysisBase::kTPC)||(fProtonAnalysisBase->GetAnalysisMode()==AliProtonAnalysisBase::kHybrid)) {
+      AliExternalTrackParam *tpcTrack = (AliExternalTrackParam *)track->GetTPCInnerParam();
+      if(!tpcTrack) continue;
+      gPt = tpcTrack->Pt();
+      gP = tpcTrack->P();
+      
+      if(fProtonAnalysisBase->IsProton(track)) {
          if(tpcTrack->Charge() > 0) {
-           fHistYPtProtons->Fill(Rapidity(tpcTrack->Px(),
-                                          tpcTrack->Py(),
-                                          tpcTrack->Pz()),
-                                 gPt);
-           //fill the container
-           containerInput[0] = Rapidity(tpcTrack->Px(),
-                                        tpcTrack->Py(),
-                                        tpcTrack->Pz());
+           nIdentifiedProtons += 1;
+           if(!fProtonAnalysisBase->IsAccepted(esd,vertex,track)) continue;//track cuts
+           if(!fProtonAnalysisBase->IsInPhaseSpace(track)) continue; //track outside the analyzed y-Pt
+           nSurvivedProtons += 1;
+           if(fProtonAnalysisBase->GetEtaMode()) {
+             fHistYPtProtons->Fill(tpcTrack->Eta(),
+                                   gPt);
+             containerInput[0] = tpcTrack->Eta();
+           }
+           else {
+             fHistYPtProtons->Fill(fProtonAnalysisBase->Rapidity(tpcTrack->Px(),
+                                                                 tpcTrack->Py(),
+                                                                 tpcTrack->Pz()),
+                                   gPt);
+             //fill the container
+             containerInput[0] = fProtonAnalysisBase->Rapidity(tpcTrack->Px(),
+                                                               tpcTrack->Py(),
+                                                               tpcTrack->Pz());
+           }
            containerInput[1] = gPt;
            fProtonContainer->Fill(containerInput,0);   
          }//protons
          else if(tpcTrack->Charge() < 0) {
-           fHistYPtAntiProtons->Fill(Rapidity(tpcTrack->Px(),
-                                              tpcTrack->Py(),
-                                              tpcTrack->Pz()),
-                                     gPt);
-           //fill the container
-           containerInput[0] = Rapidity(tpcTrack->Px(),
-                                        tpcTrack->Py(),
-                                        tpcTrack->Pz());
+           nIdentifiedAntiProtons += 1;
+           if(!fProtonAnalysisBase->IsAccepted(esd,vertex,track)) continue;//track cuts
+           if(!fProtonAnalysisBase->IsInPhaseSpace(track)) continue; //track outside the analyzed y-Pt
+           nSurvivedAntiProtons += 1;
+           if(fProtonAnalysisBase->GetEtaMode()) {
+             fHistYPtAntiProtons->Fill(tpcTrack->Eta(),
+                                       gPt);
+             containerInput[0] = tpcTrack->Eta();
+           }
+           else {
+             fHistYPtAntiProtons->Fill(fProtonAnalysisBase->Rapidity(tpcTrack->Px(),
+                                                                     tpcTrack->Py(),
+                                                                     tpcTrack->Pz()),
+                                       gPt);
+             //fill the container
+             containerInput[0] = fProtonAnalysisBase->Rapidity(tpcTrack->Px(),
+                                                               tpcTrack->Py(),
+                                                               tpcTrack->Pz());
+           }
            containerInput[1] = gPt;
            fAntiProtonContainer->Fill(containerInput,0);
          }//antiprotons   
-       }//proton check
-      }//TPC only tracks
-      else if(!fUseTPCOnly) {
-       gPt = track->Pt();
-       gP = track->P();
-       
-       //pid
-       track->GetESDpid(probability);
-       Double_t rcc = 0.0;
-       for(Int_t i = 0; i < AliPID::kSPECIES; i++) 
-         rcc += probability[i]*GetParticleFraction(i,gP);
-       if(rcc == 0.0) continue;
-       Double_t w[5];
-       for(Int_t i = 0; i < AliPID::kSPECIES; i++) 
-         w[i] = probability[i]*GetParticleFraction(i,gP)/rcc;
-       Long64_t fParticleType = TMath::LocMax(AliPID::kSPECIES,w);
-       if(fParticleType == 4) {
-         //cout<<"(Anti)protons found..."<<endl;
-         if(track->Charge() > 0) {
-           fHistYPtProtons->Fill(Rapidity(track->Px(),
-                                          track->Py(),
-                                          track->Pz()),
+      }//proton check
+    }//TPC only tracks
+    else if(fProtonAnalysisBase->GetAnalysisMode() == AliProtonAnalysisBase::kGlobal) {
+      gPt = track->Pt();
+      gP = track->P();
+      
+      if(fProtonAnalysisBase->IsProton(track)) {
+       if(track->Charge() > 0) {
+         nIdentifiedProtons += 1;
+         if(!fProtonAnalysisBase->IsAccepted(esd,vertex,track)) continue;//track cuts
+         if(!fProtonAnalysisBase->IsInPhaseSpace(track)) continue; //track outside the analyzed y-Pt
+         nSurvivedProtons += 1;
+         if(fProtonAnalysisBase->GetEtaMode()) {
+           fHistYPtProtons->Fill(track->Eta(),
+                                 gPt);
+           containerInput[0] = track->Eta();
+         }
+         else {
+           fHistYPtProtons->Fill(fProtonAnalysisBase->Rapidity(track->Px(),
+                                                               track->Py(),
+                                                               track->Pz()),
                                  gPt);
            //fill the container
-           containerInput[0] = Rapidity(track->Px(),
-                                        track->Py(),
-                                        track->Pz());
-           containerInput[1] = gPt;
-           fProtonContainer->Fill(containerInput,0);   
-         }//protons
-         else if(track->Charge() < 0) {
-           fHistYPtAntiProtons->Fill(Rapidity(track->Px(),
-                                              track->Py(),
-                                              track->Pz()),
+           containerInput[0] = fProtonAnalysisBase->Rapidity(track->Px(),
+                                                             track->Py(),
+                                                             track->Pz());
+         }
+         containerInput[1] = gPt;
+         fProtonContainer->Fill(containerInput,0);   
+       }//protons
+       else if(track->Charge() < 0) {
+         nIdentifiedAntiProtons += 1;
+         if(!fProtonAnalysisBase->IsAccepted(esd,vertex,track)) continue;//track cuts
+         if(!fProtonAnalysisBase->IsInPhaseSpace(track)) continue; //track outside the analyzed y-Pt
+         nSurvivedAntiProtons += 1;
+         if(fProtonAnalysisBase->GetEtaMode()) {
+           fHistYPtAntiProtons->Fill(track->Eta(),
+                                     gPt);
+           containerInput[0] = track->Eta();
+         }
+         else {
+           fHistYPtAntiProtons->Fill(fProtonAnalysisBase->Rapidity(track->Px(),
+                                                                   track->Py(),
+                                                                   track->Pz()),
                                      gPt);
            //fill the container
-           containerInput[0] = Rapidity(track->Px(),
-                                        track->Py(),
-                                        track->Pz());
-           containerInput[1] = gPt;
-           fAntiProtonContainer->Fill(containerInput,0);   
-         }//antiprotons
-       }//proton check 
-      }//combined tracking
-    }//cuts
+           containerInput[0] = fProtonAnalysisBase->Rapidity(track->Px(),
+                                                             track->Py(),
+                                                             track->Pz());
+         }
+         containerInput[1] = gPt;
+         fAntiProtonContainer->Fill(containerInput,0);   
+       }//antiprotons
+      }//proton check 
+    }//combined tracking
   }//track loop 
+  
+  if(fProtonAnalysisBase->GetDebugMode())
+    Printf("Initial number of tracks: %d | Identified (anti)protons: %d - %d | Survived (anti)protons: %d - %d",nTracks,nIdentifiedProtons,nIdentifiedAntiProtons,nSurvivedProtons,nSurvivedAntiProtons);
 }
 
 //____________________________________________________________________//
 void AliProtonAnalysis::Analyze(AliAODEvent* const fAOD) {
   //Main analysis part - AOD
   fHistEvents->Fill(0); //number of analyzed events
-  Int_t nGoodTracks = fAOD->GetNumberOfTracks();
-  for(Int_t iTracks = 0; iTracks < nGoodTracks; iTracks++) {
+  Int_t nTracks = fAOD->GetNumberOfTracks();
+  for(Int_t iTracks = 0; iTracks < nTracks; iTracks++) {
     AliAODTrack* track = fAOD->GetTrack(iTracks);
     Double_t gPt = track->Pt();
     Double_t gP = track->P();
@@ -705,10 +689,10 @@ void AliProtonAnalysis::Analyze(AliAODEvent* const fAOD) {
     Double_t probability[10];
     track->GetPID(probability);
     Double_t rcc = 0.0;
-    for(Int_t i = 0; i < AliPID::kSPECIESN; i++) rcc += probability[i]*GetParticleFraction(i,gP);
+    for(Int_t i = 0; i < AliPID::kSPECIESN; i++) rcc += probability[i]*fProtonAnalysisBase->GetParticleFraction(i,gP);
     if(rcc == 0.0) continue;
     Double_t w[10];
-    for(Int_t i = 0; i < AliPID::kSPECIESN; i++) w[i] = probability[i]*GetParticleFraction(i,gP)/rcc;
+    for(Int_t i = 0; i < AliPID::kSPECIESN; i++) w[i] = probability[i]*fProtonAnalysisBase->GetParticleFraction(i,gP)/rcc;
     Long64_t fParticleType = TMath::LocMax(AliPID::kSPECIESN,w);
     if(fParticleType == 4) {
       if(track->Charge() > 0) 
@@ -739,231 +723,21 @@ void AliProtonAnalysis::Analyze(AliStack* const stack,
 
     if(TMath::Abs(particle->Eta()) > 1.0) continue;
     if((particle->Pt() > fMaxPt)||(particle->Pt() < fMinPt)) continue;
-    if((Rapidity(particle->Px(),particle->Py(),particle->Pz()) > fMaxY)||(Rapidity(particle->Px(),particle->Py(),particle->Pz()) < fMinY)) continue;
+    if((fProtonAnalysisBase->Rapidity(particle->Px(),particle->Py(),particle->Pz()) > fMaxY)||(fProtonAnalysisBase->Rapidity(particle->Px(),particle->Py(),particle->Pz()) < fMinY)) continue;
 
     Int_t pdgcode = particle->GetPdgCode();
-    if(pdgcode == 2212) fHistYPtProtons->Fill(Rapidity(particle->Px(),
-                                                      particle->Py(),
-                                                      particle->Pz()),
+    if(pdgcode == 2212) fHistYPtProtons->Fill(fProtonAnalysisBase->Rapidity(particle->Px(),
+                                                                           particle->Py(),
+                                                                           particle->Pz()),
                                              particle->Pt());
-    if(pdgcode == -2212) fHistYPtAntiProtons->Fill(Rapidity(particle->Px(),
-                                                           particle->Py(),
-                                                           particle->Pz()),
+    if(pdgcode == -2212) fHistYPtAntiProtons->Fill(fProtonAnalysisBase->Rapidity(particle->Px(),
+                                                                                particle->Py(),
+                                                                                particle->Pz()),
                                                   particle->Pt());
   }//particle loop
 }
 
 //____________________________________________________________________//
-Bool_t AliProtonAnalysis::IsInPhaseSpace(AliESDtrack* const track) {
-  // Checks if the track is outside the analyzed y-Pt phase space
-  Double_t gPt = 0.0, gPx = 0.0, gPy = 0.0, gPz = 0.0;
-  Double_t eta = 0.0;
-
-  if(fUseTPCOnly) {
-    AliExternalTrackParam *tpcTrack = (AliExternalTrackParam *)track->GetTPCInnerParam();
-    if(!tpcTrack) {
-      gPt = 0.0; gPx = 0.0; gPy = 0.0; gPz = 0.0; eta = -10.0;
-    }
-    else {
-      gPt = tpcTrack->Pt();
-      gPx = tpcTrack->Px();
-      gPy = tpcTrack->Py();
-      gPz = tpcTrack->Pz();
-      eta = tpcTrack->Eta();
-    }
-  }
-  else {
-    gPt = track->Pt();
-    gPx = track->Px();
-    gPy = track->Py();
-    gPz = track->Pz();
-    eta = track->Eta();
-  }
-  
-  if((gPt < fMinPt) || (gPt > fMaxPt)) return kFALSE;
-  if(fAnalysisEtaMode) {
-    if((eta < fMinY) || (eta > fMaxY)) 
-      return kFALSE;
-  }
-  else {
-    if((Rapidity(gPx,gPy,gPz) < fMinY) || (Rapidity(gPx,gPy,gPz) > fMaxY)) 
-      return kFALSE;
-  }
-
-  return kTRUE;
-}
-
-//____________________________________________________________________//
-Bool_t AliProtonAnalysis::IsAccepted(AliESDEvent *esd,
-                                    const AliESDVertex *vertex, 
-                                    AliESDtrack* track) {
-  // Checks if the track is excluded from the cuts
-  Double_t gPt = 0.0, gPx = 0.0, gPy = 0.0, gPz = 0.0;
-  Double_t dca[2] = {0.0,0.0}, cov[3] = {0.0,0.0,0.0};  //The impact parameters and their covariance.
-  
-  if((fUseTPCOnly)&&(!fUseHybridTPC)) {
-    AliExternalTrackParam *tpcTrack = (AliExternalTrackParam *)track->GetTPCInnerParam();
-    if(!tpcTrack) {
-      gPt = 0.0; gPx = 0.0; gPy = 0.0; gPz = 0.0;
-      dca[0] = -100.; dca[1] = -100.;
-      cov[0] = -100.; cov[1] = -100.; cov[2] = -100.;
-    }
-    else {
-      gPt = tpcTrack->Pt();
-      gPx = tpcTrack->Px();
-      gPy = tpcTrack->Py();
-      gPz = tpcTrack->Pz();
-      tpcTrack->PropagateToDCA(vertex,
-                              esd->GetMagneticField(),
-                              100.,dca,cov);
-    }
-  }
-  else if(fUseHybridTPC) {
-     AliExternalTrackParam *tpcTrack = (AliExternalTrackParam *)track->GetTPCInnerParam();
-    if(!tpcTrack) {
-      gPt = 0.0; gPx = 0.0; gPy = 0.0; gPz = 0.0;
-      dca[0] = -100.; dca[1] = -100.;
-      cov[0] = -100.; cov[1] = -100.; cov[2] = -100.;
-    }
-    else {
-      gPt = tpcTrack->Pt();
-      gPx = tpcTrack->Px();
-      gPy = tpcTrack->Py();
-      gPz = tpcTrack->Pz();
-      tpcTrack->PropagateToDCA(vertex,
-                              esd->GetMagneticField(),
-                              100.,dca,cov);
-    }
-  }
-  else{
-    gPt = track->Pt();
-    gPx = track->Px();
-    gPy = track->Py();
-    gPz = track->Pz();
-    track->PropagateToDCA(vertex,
-                         esd->GetMagneticField(),
-                         100.,dca,cov);
-  }
-     
-  Int_t  fIdxInt[200];
-  Int_t nClustersITS = track->GetITSclusters(fIdxInt);
-  Int_t nClustersTPC = track->GetTPCclusters(fIdxInt);
-
-  Float_t chi2PerClusterITS = -1;
-  if (nClustersITS!=0)
-    chi2PerClusterITS = track->GetITSchi2()/Float_t(nClustersITS);
-  Float_t chi2PerClusterTPC = -1;
-  if (nClustersTPC!=0)
-    chi2PerClusterTPC = track->GetTPCchi2()/Float_t(nClustersTPC);
-
-  Double_t extCov[15];
-  track->GetExternalCovariance(extCov);
-
-  if(fPointOnITSLayer1Flag)
-    if(!track->HasPointOnITSLayer(0)) return kFALSE;
-  if(fPointOnITSLayer2Flag)
-    if(!track->HasPointOnITSLayer(1)) return kFALSE;
-  if(fPointOnITSLayer3Flag)
-    if(!track->HasPointOnITSLayer(2)) return kFALSE;
-  if(fPointOnITSLayer4Flag)
-    if(!track->HasPointOnITSLayer(3)) return kFALSE;
-  if(fPointOnITSLayer5Flag)
-    if(!track->HasPointOnITSLayer(4)) return kFALSE;
-  if(fPointOnITSLayer6Flag)
-    if(!track->HasPointOnITSLayer(5)) return kFALSE;
-  if(fMinITSClustersFlag)
-    if(nClustersITS < fMinITSClusters) return kFALSE;
-  if(fMaxChi2PerITSClusterFlag)
-    if(chi2PerClusterITS > fMaxChi2PerITSCluster) return kFALSE; 
-  if(fMinTPCClustersFlag)
-    if(nClustersTPC < fMinTPCClusters) return kFALSE;
-  if(fMaxChi2PerTPCClusterFlag)
-    if(chi2PerClusterTPC > fMaxChi2PerTPCCluster) return kFALSE; 
-  if(fMaxCov11Flag)
-    if(extCov[0] > fMaxCov11) return kFALSE;
-  if(fMaxCov22Flag)
-    if(extCov[2] > fMaxCov22) return kFALSE;
-  if(fMaxCov33Flag)
-    if(extCov[5] > fMaxCov33) return kFALSE;
-  if(fMaxCov44Flag)
-    if(extCov[9] > fMaxCov44) return kFALSE;
-  if(fMaxCov55Flag)
-    if(extCov[14] > fMaxCov55) return kFALSE;
-  if(fMaxSigmaToVertexFlag)
-    if(GetSigmaToVertex(track) > fMaxSigmaToVertex) return kFALSE;
-  if(fMaxSigmaToVertexTPCFlag)
-    if(GetSigmaToVertex(track) > fMaxSigmaToVertexTPC) return kFALSE;
-  if(fMaxDCAXYFlag) 
-    if(TMath::Abs(dca[0]) > fMaxDCAXY) return kFALSE;
-  if(fMaxDCAXYTPCFlag) 
-    if(TMath::Abs(dca[0]) > fMaxDCAXYTPC) return kFALSE;
-    if(fMaxDCAZFlag) 
-    if(TMath::Abs(dca[1]) > fMaxDCAZ) return kFALSE;
-  if(fMaxDCAZTPCFlag) 
-    if(TMath::Abs(dca[1]) > fMaxDCAZTPC) return kFALSE;
-  if(fMaxConstrainChi2Flag) {
-    if(track->GetConstrainedChi2() > 0) 
-      if(TMath::Log(track->GetConstrainedChi2()) > fMaxConstrainChi2) return kFALSE;
-  }
-  if(fITSRefitFlag)
-    if ((track->GetStatus() & AliESDtrack::kITSrefit) == 0) return kFALSE;
-  if(fTPCRefitFlag)
-    if ((track->GetStatus() & AliESDtrack::kTPCrefit) == 0) return kFALSE;
-  if(fESDpidFlag)
-    if ((track->GetStatus() & AliESDtrack::kESDpid) == 0) return kFALSE;
-  if(fTPCpidFlag)
-    if ((track->GetStatus() & AliESDtrack::kTPCpid) == 0) return kFALSE;
-
-  return kTRUE;
-}
-
-//____________________________________________________________________//
-Float_t AliProtonAnalysis::GetSigmaToVertex(AliESDtrack* esdTrack) const {
-  // Calculates the number of sigma to the vertex.
-  
-  Float_t b[2];
-  Float_t bRes[2];
-  Float_t bCov[3];
-  if((fUseTPCOnly)&&(!fUseHybridTPC))
-    esdTrack->GetImpactParametersTPC(b,bCov);
-  else
-    esdTrack->GetImpactParameters(b,bCov);
-  
-  if (bCov[0]<=0 || bCov[2]<=0) {
-    //AliDebug(1, "Estimated b resolution lower or equal zero!");
-    bCov[0]=0; bCov[2]=0;
-  }
-  bRes[0] = TMath::Sqrt(bCov[0]);
-  bRes[1] = TMath::Sqrt(bCov[2]);
-  
-  if (bRes[0] == 0 || bRes[1] ==0) return -1;
-  
-  Float_t d = TMath::Sqrt(TMath::Power(b[0]/bRes[0],2) + TMath::Power(b[1]/bRes[1],2));
-  
-  if (TMath::Exp(-d * d / 2) < 1e-10) return 1000;
-  
-  d = TMath::ErfInverse(1 - TMath::Exp(-d * d / 2)) * TMath::Sqrt(2);
-  
-  return d;
-}
-
-//____________________________________________________________________//
-Double_t AliProtonAnalysis::Rapidity(Double_t gPx, Double_t gPy, Double_t gPz) const {
-  //returns the rapidity of the proton - to be removed
-  Double_t fMass = 9.38270000000000048e-01;
-  
-  Double_t gP = TMath::Sqrt(TMath::Power(gPx,2) + 
-                           TMath::Power(gPy,2) + 
-                          TMath::Power(gPz,2));
-  Double_t energy = TMath::Sqrt(gP*gP + fMass*fMass);
-  Double_t y = -999;
-  if(energy != gPz) 
-    y = 0.5*TMath::Log((energy + gPz)/(energy - gPz));
-
-  return y;
-}
-
-//____________________________________________________________________//
 Bool_t AliProtonAnalysis::PrintMean(TH1 *hist, Double_t edge) {
   //calculates the mean value of the ratio/asymmetry within \pm edge
   Double_t sum = 0.0;
@@ -1037,11 +811,13 @@ void AliProtonAnalysis::Correct(Int_t step) {
   fCorrectProtons = new AliCFDataGrid("correctProtons",
                                      "corrected data",
                                      *fProtonContainer);
+  fCorrectProtons->SetMeasured(0);
   fCorrectProtons->ApplyEffCorrection(*(AliCFEffGrid *)fEffGridListProtons->At(step));
 
   fCorrectAntiProtons = new AliCFDataGrid("correctAntiProtons",
                                          "corrected data",
                                          *fAntiProtonContainer);
+  fCorrectAntiProtons->SetMeasured(0);
   fCorrectAntiProtons->ApplyEffCorrection(*(AliCFEffGrid *)fEffGridListAntiProtons->At(step));
 }
 
index d2f3360..3a485ef 100644 (file)
@@ -31,6 +31,7 @@ class AliESDtrack;
 class AliExternalTrackParam;
 class AliStack;
 class AliESDVertex;
+class AliProtonAnalysisBase;
 
 class AliProtonAnalysis : public TObject {
  public:
@@ -39,11 +40,11 @@ class AliProtonAnalysis : public TObject {
                    Int_t nbinsPt, Float_t fLowPt, Float_t fHighPt);
   virtual ~AliProtonAnalysis();
 
-  void SetEtaMode() {fAnalysisEtaMode = kTRUE;}
+  void SetBaseAnalysis(AliProtonAnalysisBase *baseAnalysis) {
+    fProtonAnalysisBase = baseAnalysis;}
+  AliProtonAnalysisBase *GetProtonAnalysisBaseObject() {
+    return fProtonAnalysisBase;}
 
-  void UseTPCOnly() {fUseTPCOnly = kTRUE;}
-  void UseHybridTPC() {fUseHybridTPC = kTRUE;}
-  
   void InitAnalysisHistograms(Int_t nbinsY, Float_t fLowY, Float_t fHighY,
                              Int_t nbinsPt, Float_t fLowPt, Float_t fHighPt);
   Bool_t ReadFromFile(const char* filename);
@@ -81,89 +82,6 @@ class AliProtonAnalysis : public TObject {
   Bool_t  PrintMean(TH1 *hist, Double_t edge);
   Bool_t  PrintYields(TH1 *hist, Double_t edge); 
 
-  //Cut functions
-  void    SetPointOnITSLayer1() {fPointOnITSLayer1Flag = kTRUE;}
-  void    SetPointOnITSLayer2() {fPointOnITSLayer2Flag = kTRUE;}
-  void    SetPointOnITSLayer3() {fPointOnITSLayer3Flag = kTRUE;}
-  void    SetPointOnITSLayer4() {fPointOnITSLayer4Flag = kTRUE;}
-  void    SetPointOnITSLayer5() {fPointOnITSLayer5Flag = kTRUE;}
-  void    SetPointOnITSLayer6() {fPointOnITSLayer6Flag = kTRUE;}
-  void    SetMinITSClusters(Int_t minITSClusters) {
-    fMinITSClusters = minITSClusters;
-    fMinITSClustersFlag = kTRUE;
-  }
-  void    SetMaxChi2PerITSCluster(Double_t maxChi2PerITSCluster) {
-    fMaxChi2PerITSCluster = maxChi2PerITSCluster;
-    fMaxChi2PerITSClusterFlag = kTRUE;
-  }
-  void    SetMinTPCClusters(Int_t minTPCClusters) {
-    fMinTPCClusters = minTPCClusters;
-    fMinTPCClustersFlag = kTRUE;
-  }
-  void    SetMaxChi2PerTPCCluster(Double_t maxChi2PerTPCCluster) {
-    fMaxChi2PerTPCCluster = maxChi2PerTPCCluster;
-    fMaxChi2PerTPCClusterFlag = kTRUE;
-  }
-  void    SetMaxCov11(Double_t maxCov11) {
-    fMaxCov11 = maxCov11; fMaxCov11Flag = kTRUE;}
-  void    SetMaxCov22(Double_t maxCov22) {
-    fMaxCov22 = maxCov22; fMaxCov22Flag = kTRUE;}
-  void    SetMaxCov33(Double_t maxCov33) {
-    fMaxCov33 = maxCov33; fMaxCov33Flag = kTRUE;}
-  void    SetMaxCov44(Double_t maxCov44) {
-    fMaxCov44 = maxCov44; fMaxCov44Flag = kTRUE;}
-  void    SetMaxCov55(Double_t maxCov55) {
-    fMaxCov55 = maxCov55; fMaxCov55Flag = kTRUE;}
-  void    SetMaxSigmaToVertex(Double_t maxSigmaToVertex) {
-    fMaxSigmaToVertex = maxSigmaToVertex;
-    fMaxSigmaToVertexFlag = kTRUE;
-  }
-  void    SetMaxSigmaToVertexTPC(Double_t maxSigmaToVertex) {
-    fMaxSigmaToVertexTPC = maxSigmaToVertex;
-    fMaxSigmaToVertexTPCFlag = kTRUE;
-  }
-  void    SetMaxDCAXY(Double_t maxDCAXY) {
-    fMaxDCAXY = maxDCAXY;
-    fMaxDCAXYFlag = kTRUE;
-  }
-  void    SetMaxDCAXYTPC(Double_t maxDCAXY) {
-    fMaxDCAXYTPC = maxDCAXY;
-    fMaxDCAXYTPCFlag = kTRUE;
-  }
-  void    SetMaxDCAZ(Double_t maxDCAZ) {
-    fMaxDCAZ = maxDCAZ;
-    fMaxDCAZFlag = kTRUE;
-  }
-  void    SetMaxDCAZTPC(Double_t maxDCAZ) {
-    fMaxDCAZTPC = maxDCAZ;
-    fMaxDCAZTPCFlag = kTRUE;
-  }
-  void    SetMaxConstrainChi2(Double_t maxConstrainChi2) {
-    fMaxConstrainChi2 = maxConstrainChi2;
-    fMaxConstrainChi2Flag = kTRUE;
-  }
-  void    SetITSRefit() {fITSRefitFlag = kTRUE;}
-  void    SetTPCRefit() {fTPCRefitFlag = kTRUE;}
-  void    SetESDpid() {fESDpidFlag = kTRUE;}
-  void    SetTPCpid() {fTPCpidFlag = kTRUE;}
-
-  //Prior probabilities
-  void SetPriorProbabilities(Double_t * const partFrac) {
-    for(Int_t i = 0; i < AliPID::kSPECIESN; i++) fPartFrac[i] = partFrac[i];} 
-  void SetPriorProbabilityFunctions(TF1 *const felectron, 
-                                   TF1 *const fmuon, 
-                                   TF1 *const fpion, 
-                                   TF1 *const fkaon, 
-                                   TF1 *const fproton) {
-    fFunctionProbabilityFlag = kTRUE;
-    fElectronFunction = felectron; 
-    fMuonFunction = fmuon; 
-    fPionFunction = fpion;
-    fKaonFunction = fkaon;
-    fProtonFunction = fproton;
-  } 
-  Double_t GetParticleFraction(Int_t i, Double_t p);
-
   //interface to the correction framework
   void Correct(Int_t step);
   Bool_t ReadCorrectionContainer(const char* filename);
@@ -183,58 +101,14 @@ class AliProtonAnalysis : public TObject {
  private:
   AliProtonAnalysis(const AliProtonAnalysis&); // Not implemented
   AliProtonAnalysis& operator=(const AliProtonAnalysis&); // Not implemented
-
-  Bool_t   IsAccepted(AliESDEvent *esd,
-                     const AliESDVertex *vertex, 
-                     AliESDtrack *track);
-  Bool_t   IsInPhaseSpace(AliESDtrack *track);
-
-  Float_t  GetSigmaToVertex(AliESDtrack* esdTrack) const; 
-  Double_t Rapidity(Double_t Px, Double_t Py, Double_t Pz) const;
   
-  Bool_t fAnalysisEtaMode; //run the analysis in eta or y
+  AliProtonAnalysisBase *fProtonAnalysisBase;//base analysis object
 
-  Int_t fNBinsY; //number of bins in y
-  Float_t fMinY, fMaxY; //min & max value of y
+  Int_t fNBinsY; //number of bins in y or eta
+  Double_t fMinY, fMaxY; //min & max value of y or eta
   Int_t fNBinsPt;  //number of bins in pT
-  Float_t fMinPt, fMaxPt; //min & max value of pT
-  
-  //cuts
-  Int_t fMinTPCClusters, fMinITSClusters; //min TPC & ITS clusters
-  Double_t fMaxChi2PerTPCCluster, fMaxChi2PerITSCluster; //max chi2 per TPC & ITS cluster
-  Double_t fMaxCov11, fMaxCov22, fMaxCov33, fMaxCov44, fMaxCov55; //max values of cov. matrix
-  Double_t fMaxSigmaToVertex; //max sigma to vertex cut
-  Double_t fMaxSigmaToVertexTPC; //max sigma to vertex cut
-  Double_t fMaxDCAXY, fMaxDCAXYTPC; //max DCA xy
-  Double_t fMaxDCAZ, fMaxDCAZTPC; //max DCA z
-  Double_t fMaxConstrainChi2; //max constrain chi2 - vertex
-  Bool_t fMinTPCClustersFlag, fMinITSClustersFlag; //shows if this cut is used or not
-  Bool_t fMaxChi2PerTPCClusterFlag, fMaxChi2PerITSClusterFlag; //shows if this cut is used or not
-  Bool_t fMaxCov11Flag, fMaxCov22Flag, fMaxCov33Flag, fMaxCov44Flag, fMaxCov55Flag; //shows if this cut is used or not
-  Bool_t fMaxSigmaToVertexFlag; //shows if this cut is used or not
-  Bool_t fMaxSigmaToVertexTPCFlag; //shows if this cut is used or not
-  Bool_t fMaxDCAXYFlag, fMaxDCAXYTPCFlag; //shows if this cut is used or not
-  Bool_t fMaxDCAZFlag, fMaxDCAZTPCFlag; //shows if this cut is used or not
-  Bool_t fMaxConstrainChi2Flag; //shows if this cut is used or not
-  Bool_t fITSRefitFlag, fTPCRefitFlag; //shows if this cut is used or not
-  Bool_t fESDpidFlag, fTPCpidFlag; //shows if this cut is used or not
-  Bool_t fPointOnITSLayer1Flag, fPointOnITSLayer2Flag; //shows if this cut is used or not
-  Bool_t fPointOnITSLayer3Flag, fPointOnITSLayer4Flag; //shows if this cut is used or not
-  Bool_t fPointOnITSLayer5Flag, fPointOnITSLayer6Flag; //shows if this cut is used or not
+  Double_t fMinPt, fMaxPt; //min & max value of pT
   
-  //pid
-  Bool_t fFunctionProbabilityFlag; //flag: kTRUE if functions used
-  Double_t fPartFrac[10]; //prior probabilities
-  TF1  *fElectronFunction; //momentum dependence of the prior probs
-  TF1  *fMuonFunction; //momentum dependence of the prior probs
-  TF1  *fPionFunction; //momentum dependence of the prior probs
-  TF1  *fKaonFunction; //momentum dependence of the prior probs
-  TF1  *fProtonFunction; //momentum dependence of the prior probs
-
-  //Detectors
-  Bool_t fUseTPCOnly; //kTRUE if TPC only information is used
-  Bool_t fUseHybridTPC; //kTRUE if TPC info is used for momentum - PID and ITS for vertex & points
-
   //Analysis containers
   AliCFContainer *fProtonContainer; //container for protons
   AliCFContainer *fAntiProtonContainer; //container for antiprotons
diff --git a/PWG2/SPECTRA/AliProtonAnalysisBase.cxx b/PWG2/SPECTRA/AliProtonAnalysisBase.cxx
new file mode 100644 (file)
index 0000000..23f6ae0
--- /dev/null
@@ -0,0 +1,803 @@
+/**************************************************************************
+ * Author: Panos Christakoglou.                                           *
+ * Contributors are mentioned in the code where appropriate.              *
+ *                                                                        *
+ * Permission to use, copy, modify and distribute this software and its   *
+ * documentation strictly for non-commercial purposes is hereby granted   *
+ * without fee, provided that the above copyright notice appears in all   *
+ * copies and that both the copyright notice and this permission notice   *
+ * appear in the supporting documentation. The authors make no claims     *
+ * about the suitability of this software for any purpose. It is          *
+ * provided "as is" without express or implied warranty.                  *
+ **************************************************************************/
+
+/* $Id: AliProtonAnalysisBase.cxx 31056 2009-02-16 14:31:41Z pchrist $ */
+
+//-----------------------------------------------------------------
+//                 AliProtonAnalysisBase class
+//   This is the class to deal with the proton analysis
+//   Origin: Panos Christakoglou | Panos.Christakoglou@cern.ch
+//-----------------------------------------------------------------
+#include <Riostream.h>
+#include <TFile.h>
+#include <TSystem.h>
+#include <TF1.h>
+#include <TH2D.h>
+#include <TH1D.h>
+#include <TH1I.h>
+#include <TParticle.h>
+#include <TCanvas.h>
+#include <TLatex.h>
+
+#include "AliProtonAnalysisBase.h"
+
+#include <AliExternalTrackParam.h>
+#include <AliESDEvent.h>
+#include <AliLog.h>
+#include <AliPID.h>
+#include <AliESDVertex.h>
+#include <AliVertexerTracks.h>
+
+ClassImp(AliProtonAnalysisBase)
+
+//____________________________________________________________________//
+AliProtonAnalysisBase::AliProtonAnalysisBase() : 
+  TObject(),  fProtonAnalysisLevel("ESD"),
+  fTriggerMode(kMB2), fProtonAnalysisMode(kTPC), fProtonPIDMode(kBayesian),
+  fAnalysisEtaMode(kFALSE),
+  fVxMax(100.), fVyMax(100.), fVzMax(100.),
+  fNBinsX(0), fMinX(0), fMaxX(0),
+  fNBinsY(0), fMinY(0), fMaxY(0),
+  fMinTPCClusters(0), fMinITSClusters(0),
+  fMaxChi2PerTPCCluster(0), fMaxChi2PerITSCluster(0),
+  fMaxCov11(0), fMaxCov22(0), fMaxCov33(0), fMaxCov44(0), fMaxCov55(0),
+  fMaxSigmaToVertex(0), fMaxSigmaToVertexTPC(0),
+  fMaxDCAXY(0), fMaxDCAXYTPC(0),
+  fMaxDCAZ(0), fMaxDCAZTPC(0),
+  fMaxDCA3D(0), fMaxDCA3DTPC(0),
+  fMaxConstrainChi2(0),
+  fMinTPCClustersFlag(kFALSE), fMinITSClustersFlag(kFALSE),
+  fMaxChi2PerTPCClusterFlag(kFALSE), fMaxChi2PerITSClusterFlag(kFALSE),
+  fMaxCov11Flag(kFALSE), fMaxCov22Flag(kFALSE), 
+  fMaxCov33Flag(kFALSE), fMaxCov44Flag(kFALSE), fMaxCov55Flag(kFALSE),
+  fMaxSigmaToVertexFlag(kFALSE), fMaxSigmaToVertexTPCFlag(kFALSE),
+  fMaxDCAXYFlag(kFALSE), fMaxDCAXYTPCFlag(kFALSE),
+  fMaxDCAZFlag(kFALSE), fMaxDCAZTPCFlag(kFALSE),
+  fMaxDCA3DFlag(kFALSE), fMaxDCA3DTPCFlag(kFALSE),
+  fMaxConstrainChi2Flag(kFALSE),
+  fITSRefitFlag(kFALSE), fTPCRefitFlag(kFALSE),
+  fESDpidFlag(kFALSE), fTPCpidFlag(kFALSE),
+  fPointOnITSLayer1Flag(0), fPointOnITSLayer2Flag(0),
+  fPointOnITSLayer3Flag(0), fPointOnITSLayer4Flag(0),
+  fPointOnITSLayer5Flag(0), fPointOnITSLayer6Flag(0),
+  fFunctionProbabilityFlag(kFALSE), 
+  fElectronFunction(0), fMuonFunction(0),
+  fPionFunction(0), fKaonFunction(0), fProtonFunction(0),
+  fDebugMode(kFALSE) {
+  //Default constructor
+  for(Int_t i = 0; i < 5; i++) fPartFrac[i] = 0.0;
+}
+
+//____________________________________________________________________//
+AliProtonAnalysisBase::~AliProtonAnalysisBase() {
+  //Default destructor
+  if(fElectronFunction) delete fElectronFunction;
+  if(fMuonFunction) delete fMuonFunction;
+  if(fPionFunction) delete fPionFunction;
+  if(fKaonFunction) delete fKaonFunction;
+  if(fProtonFunction) delete fProtonFunction;
+}
+
+//____________________________________________________________________//
+Double_t AliProtonAnalysisBase::GetParticleFraction(Int_t i, Double_t p) {
+  //Return the a priori probs
+  Double_t partFrac=0;
+  if(fFunctionProbabilityFlag) {
+    if(i == 0) partFrac = fElectronFunction->Eval(p);
+    if(i == 1) partFrac = fMuonFunction->Eval(p);
+    if(i == 2) partFrac = fPionFunction->Eval(p);
+    if(i == 3) partFrac = fKaonFunction->Eval(p);
+    if(i == 4) partFrac = fProtonFunction->Eval(p);
+  }
+  else partFrac = fPartFrac[i];
+
+  return partFrac;
+}
+
+//____________________________________________________________________//
+Bool_t AliProtonAnalysisBase::IsInPhaseSpace(AliESDtrack* const track) {
+  // Checks if the track is outside the analyzed y-Pt phase space
+  Double_t gPt = 0.0, gPx = 0.0, gPy = 0.0, gPz = 0.0;
+  Double_t eta = 0.0;
+
+  if(kTPC) {
+    AliExternalTrackParam *tpcTrack = (AliExternalTrackParam *)track->GetTPCInnerParam();
+    if(!tpcTrack) {
+      gPt = 0.0; gPx = 0.0; gPy = 0.0; gPz = 0.0; eta = -10.0;
+    }
+    else {
+      gPt = tpcTrack->Pt();
+      gPx = tpcTrack->Px();
+      gPy = tpcTrack->Py();
+      gPz = tpcTrack->Pz();
+      eta = tpcTrack->Eta();
+    }
+  }
+  else {
+    gPt = track->Pt();
+    gPx = track->Px();
+    gPy = track->Py();
+    gPz = track->Pz();
+    eta = track->Eta();
+  }
+  
+  if((gPt < fMinY) || (gPt > fMaxY)) {
+      if(fDebugMode)
+       Printf("IsInPhaseSpace: Track rejected because it has a Pt value of %lf (accepted interval: %lf - %lf)",gPt,fMinY,fMaxY);
+      return kFALSE;
+  }
+  if(fAnalysisEtaMode) {
+    if((eta < fMinX) || (eta > fMaxX)) {
+      if(fDebugMode)
+       Printf("IsInPhaseSpace: Track rejected because it has an eta value of %lf (accepted interval: %lf - %lf)",eta,fMinX,fMaxX);
+      return kFALSE;
+    }
+  }
+  else {
+    if((Rapidity(gPx,gPy,gPz) < fMinX) || (Rapidity(gPx,gPy,gPz) > fMaxX)) {
+      if(fDebugMode)
+       Printf("IsInPhaseSpace: Track rejected because it has a y value of %lf (accepted interval: %lf - %lf)",Rapidity(gPx,gPy,gPz),fMinX,fMaxX);
+      return kFALSE;
+    }
+  }
+
+  return kTRUE;
+}
+
+//____________________________________________________________________//
+Bool_t AliProtonAnalysisBase::IsAccepted(AliESDEvent *esd,
+                                        const AliESDVertex *vertex, 
+                                        AliESDtrack* track) {
+  // Checks if the track is excluded from the cuts
+  Double_t gPt = 0.0, gPx = 0.0, gPy = 0.0, gPz = 0.0;
+  Double_t dca[2] = {0.0,0.0}, cov[3] = {0.0,0.0,0.0};  //The impact parameters and their covariance.
+  Double_t dca3D = 0.0;
+  
+  if((kTPC)&&(!kHybrid)) {
+    AliExternalTrackParam *tpcTrack = (AliExternalTrackParam *)track->GetTPCInnerParam();
+    if(!tpcTrack) {
+      gPt = 0.0; gPx = 0.0; gPy = 0.0; gPz = 0.0;
+      dca[0] = -100.; dca[1] = -100.; dca3D = -100.;
+      cov[0] = -100.; cov[1] = -100.; cov[2] = -100.;
+    }
+    else {
+      gPt = tpcTrack->Pt();
+      gPx = tpcTrack->Px();
+      gPy = tpcTrack->Py();
+      gPz = tpcTrack->Pz();
+      tpcTrack->PropagateToDCA(vertex,
+                              esd->GetMagneticField(),
+                              100.,dca,cov);
+    }
+  }
+  else if(kHybrid) {
+     AliExternalTrackParam *tpcTrack = (AliExternalTrackParam *)track->GetTPCInnerParam();
+    if(!tpcTrack) {
+      gPt = 0.0; gPx = 0.0; gPy = 0.0; gPz = 0.0;
+      dca[0] = -100.; dca[1] = -100.;
+      cov[0] = -100.; cov[1] = -100.; cov[2] = -100.;
+    }
+    else {
+      gPt = tpcTrack->Pt();
+      gPx = tpcTrack->Px();
+      gPy = tpcTrack->Py();
+      gPz = tpcTrack->Pz();
+      tpcTrack->PropagateToDCA(vertex,
+                              esd->GetMagneticField(),
+                              100.,dca,cov);
+    }
+  }
+  else{
+    gPt = track->Pt();
+    gPx = track->Px();
+    gPy = track->Py();
+    gPz = track->Pz();
+    track->PropagateToDCA(vertex,
+                         esd->GetMagneticField(),
+                         100.,dca,cov);
+  }
+  dca3D = TMath::Sqrt(TMath::Power(dca[0],2) +
+                     TMath::Power(dca[1],2));
+     
+  Int_t  fIdxInt[200];
+  Int_t nClustersITS = track->GetITSclusters(fIdxInt);
+  Int_t nClustersTPC = track->GetTPCclusters(fIdxInt);
+
+  Float_t chi2PerClusterITS = -1;
+  if (nClustersITS!=0)
+    chi2PerClusterITS = track->GetITSchi2()/Float_t(nClustersITS);
+  Float_t chi2PerClusterTPC = -1;
+  if (nClustersTPC!=0)
+    chi2PerClusterTPC = track->GetTPCchi2()/Float_t(nClustersTPC);
+
+  Double_t extCov[15];
+  track->GetExternalCovariance(extCov);
+
+  if(fPointOnITSLayer1Flag) {
+    if(!track->HasPointOnITSLayer(0)) {
+      if(fDebugMode)
+       Printf("IsAccepted: Track rejected because it doesn't have a point on the 1st ITS layer");
+      return kFALSE;
+    }
+  }
+  if(fPointOnITSLayer2Flag) {
+    if(!track->HasPointOnITSLayer(1)) {
+      if(fDebugMode)
+       Printf("IsAccepted: Track rejected because it doesn't have a point on the 2nd ITS layer");
+      return kFALSE;
+    }
+  }
+  if(fPointOnITSLayer3Flag) {
+    if(!track->HasPointOnITSLayer(2)) {
+      if(fDebugMode)
+       Printf("IsAccepted: Track rejected because it doesn't have a point on the 3rd ITS layer");
+      return kFALSE;
+    }
+  }
+  if(fPointOnITSLayer4Flag) {
+    if(!track->HasPointOnITSLayer(3)) {
+      if(fDebugMode)
+       Printf("IsAccepted: Track rejected because it doesn't have a point on the 4th ITS layer");
+      return kFALSE;
+    }
+  }
+  if(fPointOnITSLayer5Flag) {
+    if(!track->HasPointOnITSLayer(4)) {
+      if(fDebugMode)
+       Printf("IsAccepted: Track rejected because it doesn't have a point on the 5th ITS layer");
+      return kFALSE;
+    }
+  }
+  if(fPointOnITSLayer6Flag) {
+    if(!track->HasPointOnITSLayer(5)) {
+      if(fDebugMode)
+       Printf("IsAccepted: Track rejected because it doesn't have a point on the 6th ITS layer");
+      return kFALSE;
+    }
+  }
+  if(fMinITSClustersFlag) {
+    if(nClustersITS < fMinITSClusters) {
+      if(fDebugMode)
+       Printf("IsAccepted: Track rejected because it has %d ITS points (min. requested: %d)",nClustersITS,fMinITSClusters);
+      return kFALSE;
+    }
+  }
+  if(fMaxChi2PerITSClusterFlag) {
+    if(chi2PerClusterITS > fMaxChi2PerITSCluster) {
+      if(fDebugMode)
+       Printf("IsAccepted: Track rejected because it has a chi2 per ITS cluster %lf (max. requested: %lf)",chi2PerClusterITS,fMaxChi2PerITSCluster);
+      return kFALSE; 
+    }
+  }
+  if(fMinTPCClustersFlag) {
+    if(nClustersTPC < fMinTPCClusters) {
+      if(fDebugMode)
+       Printf("IsAccepted: Track rejected because it has %d TPC clusters (min. requested: %d)",nClustersTPC,fMinTPCClusters);
+      return kFALSE;
+    }
+  }
+  if(fMaxChi2PerTPCClusterFlag) {
+    if(chi2PerClusterTPC > fMaxChi2PerTPCCluster) {
+      if(fDebugMode)
+       Printf("IsAccepted: Track rejected because it has a chi2 per TPC cluster %lf (max. requested: %lf)",chi2PerClusterTPC,fMaxChi2PerTPCCluster);
+      return kFALSE; 
+    }
+  }
+  if(fMaxCov11Flag) {
+    if(extCov[0] > fMaxCov11) {
+      if(fDebugMode)
+       Printf("IsAccepted: Track rejected because it has a cov11 value of %lf (max. requested: %lf)",extCov[0],fMaxCov11);
+      return kFALSE;
+    }
+  }
+  if(fMaxCov22Flag) {
+    if(extCov[2] > fMaxCov22) {
+      if(fDebugMode)
+       Printf("IsAccepted: Track rejected because it has a cov22 value of %lf (max. requested: %lf)",extCov[2],fMaxCov22);
+      return kFALSE;
+    }
+  }
+  if(fMaxCov33Flag) {
+    if(extCov[5] > fMaxCov33) {
+      if(fDebugMode)
+       Printf("IsAccepted: Track rejected because it has a cov33 value of %lf (max. requested: %lf)",extCov[5],fMaxCov33);
+      return kFALSE;
+    }
+  }
+  if(fMaxCov44Flag) {
+    if(extCov[9] > fMaxCov44) {
+      if(fDebugMode)
+       Printf("IsAccepted: Track rejected because it has a cov44 value of %lf (max. requested: %lf)",extCov[9],fMaxCov44);
+      return kFALSE;
+    }
+  }
+  if(fMaxCov55Flag) {
+    if(extCov[14] > fMaxCov55) {
+      if(fDebugMode)
+       Printf("IsAccepted: Track rejected because it has a cov55 value of %lf (max. requested: %lf)",extCov[14],fMaxCov55);
+      return kFALSE;
+    }
+  }
+  if(fMaxSigmaToVertexFlag) {
+    if(GetSigmaToVertex(track) > fMaxSigmaToVertex) {
+      if(fDebugMode)
+       Printf("IsAccepted: Track rejected because it has a %lf sigmas to vertex (max. requested: %lf)",GetSigmaToVertex(track),fMaxSigmaToVertex);
+      return kFALSE;
+    }
+  }
+  if(fMaxSigmaToVertexTPCFlag) {
+    if(GetSigmaToVertex(track) > fMaxSigmaToVertexTPC) {
+      if(fDebugMode)
+       Printf("IsAccepted: Track rejected because it has a %lf sigmas to vertex TPC (max. requested: %lf)",GetSigmaToVertex(track),fMaxSigmaToVertexTPC);
+      return kFALSE;
+    }
+  }
+  if(fMaxDCAXYFlag) { 
+    if(TMath::Abs(dca[0]) > fMaxDCAXY) {
+      if(fDebugMode)
+       Printf("IsAccepted: Track rejected because it has a value of dca(xy) of %lf (max. requested: %lf)",TMath::Abs(dca[0]),fMaxDCAXY);
+      return kFALSE;
+    }
+  }
+  if(fMaxDCAXYTPCFlag) { 
+    if(TMath::Abs(dca[0]) > fMaxDCAXYTPC) {
+      if(fDebugMode)
+       Printf("IsAccepted: Track rejected because it has a value of dca(xy) (TPC) of %lf (max. requested: %lf)",TMath::Abs(dca[0]),fMaxDCAXYTPC);
+      return kFALSE;
+    }
+  }
+  if(fMaxDCAZFlag) { 
+    if(TMath::Abs(dca[1]) > fMaxDCAZ) {
+      if(fDebugMode)
+       Printf("IsAccepted: Track rejected because it has a value of dca(z) of %lf (max. requested: %lf)",TMath::Abs(dca[1]),fMaxDCAZ);
+      return kFALSE;
+    }
+  }
+  if(fMaxDCAZTPCFlag) { 
+    if(TMath::Abs(dca[1]) > fMaxDCAZTPC) {
+      if(fDebugMode)
+       Printf("IsAccepted: Track rejected because it has a value of dca(z) (TPC) of %lf (max. requested: %lf)",TMath::Abs(dca[1]),fMaxDCAZTPC);
+      return kFALSE;
+    }
+  }
+  if(fMaxDCA3DFlag) { 
+    if(TMath::Abs(dca3D) > fMaxDCA3D) {
+      if(fDebugMode)
+       Printf("IsAccepted: Track rejected because it has a value of dca(3D) of %lf (max. requested: %lf)",TMath::Abs(dca3D),fMaxDCA3D);
+      return kFALSE;
+    }
+  }
+  if(fMaxDCA3DTPCFlag) { 
+    if(TMath::Abs(dca3D) > fMaxDCA3DTPC)  {
+      if(fDebugMode)
+       Printf("IsAccepted: Track rejected because it has a value of dca(3D) (TPC) of %lf (max. requested: %lf)",TMath::Abs(dca3D),fMaxDCA3DTPC);
+      return kFALSE;
+    }
+  }
+  if(fMaxConstrainChi2Flag) {
+    if(track->GetConstrainedChi2() > 0) 
+      if(TMath::Log(track->GetConstrainedChi2()) > fMaxConstrainChi2)  {
+      if(fDebugMode)
+       Printf("IsAccepted: Track rejected because it has a value of the constrained chi2 to the vertex of %lf (max. requested: %lf)",TMath::Log(track->GetConstrainedChi2()),fMaxConstrainChi2);
+      return kFALSE;
+      }
+  }
+  if(fITSRefitFlag) {
+    if ((track->GetStatus() & AliESDtrack::kITSrefit) == 0) {
+      if(fDebugMode)
+       Printf("IsAccepted: Track rejected because it has no ITS refit flag");
+      return kFALSE;
+    }
+  }
+  if(fTPCRefitFlag) {
+    if ((track->GetStatus() & AliESDtrack::kTPCrefit) == 0) {
+      if(fDebugMode)
+       Printf("IsAccepted: Track rejected because it has no TPC refit flag");
+      return kFALSE;
+    }
+  }
+  if(fESDpidFlag) {
+    if ((track->GetStatus() & AliESDtrack::kESDpid) == 0) {
+      if(fDebugMode)
+       Printf("IsAccepted: Track rejected because it has no ESD pid flag");
+      return kFALSE;
+    }
+  }
+  if(fTPCpidFlag) {
+    if ((track->GetStatus() & AliESDtrack::kTPCpid) == 0) {
+      if(fDebugMode)
+       Printf("IsAccepted: Track rejected because it has no TPC pid flag");
+      return kFALSE;
+    }
+  }
+
+  return kTRUE;
+}
+
+//____________________________________________________________________//
+Float_t AliProtonAnalysisBase::GetSigmaToVertex(AliESDtrack* esdTrack) const {
+  // Calculates the number of sigma to the vertex.
+  
+  Float_t b[2];
+  Float_t bRes[2];
+  Float_t bCov[3];
+  if((kTPC)&&(!kHybrid))
+    esdTrack->GetImpactParametersTPC(b,bCov);
+  else
+    esdTrack->GetImpactParameters(b,bCov);
+  
+  if (bCov[0]<=0 || bCov[2]<=0) {
+    //AliDebug(1, "Estimated b resolution lower or equal zero!");
+    bCov[0]=0; bCov[2]=0;
+  }
+  bRes[0] = TMath::Sqrt(bCov[0]);
+  bRes[1] = TMath::Sqrt(bCov[2]);
+  
+  if (bRes[0] == 0 || bRes[1] ==0) return -1;
+  
+  Float_t d = TMath::Sqrt(TMath::Power(b[0]/bRes[0],2) + TMath::Power(b[1]/bRes[1],2));
+  
+  if (TMath::Exp(-d * d / 2) < 1e-10) return 1000;
+  
+  d = TMath::ErfInverse(1 - TMath::Exp(-d * d / 2)) * TMath::Sqrt(2);
+  
+  return d;
+}
+
+//____________________________________________________________________//
+Double_t AliProtonAnalysisBase::Rapidity(Double_t gPx, 
+                                        Double_t gPy, 
+                                        Double_t gPz) const {
+  //returns the rapidity of the proton - to be removed
+  Double_t fMass = 9.38270000000000048e-01;
+  
+  Double_t gP = TMath::Sqrt(TMath::Power(gPx,2) + 
+                           TMath::Power(gPy,2) + 
+                          TMath::Power(gPz,2));
+  Double_t energy = TMath::Sqrt(gP*gP + fMass*fMass);
+  Double_t y = -999;
+  if(energy != gPz) 
+    y = 0.5*TMath::Log((energy + gPz)/(energy - gPz));
+
+  return y;
+}
+
+//________________________________________________________________________
+const AliESDVertex* AliProtonAnalysisBase::GetVertex(AliESDEvent* esd, 
+                                                    AnalysisMode mode,
+                                                    Double_t gVxMax,
+                                                    Double_t gVyMax,
+                                                    Double_t gVzMax) {
+  // Get the vertex from the ESD and returns it if the vertex is valid
+  // Second argument decides which vertex is used (this selects
+  // also the quality criteria that are applied)
+  const AliESDVertex* vertex = 0;
+  if(mode == kHybrid)
+    vertex = esd->GetPrimaryVertexSPD();
+  else if(mode == kTPC){
+    Double_t kBz = esd->GetMagneticField();
+    AliVertexerTracks vertexer(kBz);
+    vertexer.SetTPCMode();
+    AliESDVertex *vTPC = vertexer.FindPrimaryVertex(esd);
+    esd->SetPrimaryVertexTPC(vTPC);
+    for (Int_t i=0; i<esd->GetNumberOfTracks(); i++) {
+      AliESDtrack *t = esd->GetTrack(i);
+      t->RelateToVertexTPC(vTPC, kBz, kVeryBig);
+    }
+    delete vTPC;
+    vertex = esd->GetPrimaryVertexTPC();
+  }
+  else if(mode == kGlobal)
+    vertex = esd->GetPrimaryVertex();
+  else
+    Printf("GetVertex: ERROR: Invalid second argument %d", mode);
+  
+  if(!vertex) {
+    if(fDebugMode)
+      Printf("GetVertex: Event rejected because there is no valid vertex object");
+    return 0;
+  }
+  
+  // check Ncontributors
+  if(vertex->GetNContributors() <= 0) {
+    if(fDebugMode)
+      Printf("GetVertex: Event rejected because the number of contributors for the vertex determination is <= 0");
+    return 0;
+  }
+  
+  // check resolution
+  Double_t zRes = vertex->GetZRes();
+  if(zRes == 0) {
+    if(fDebugMode)
+      Printf("GetVertex: Event rejected because the value of the vertex resolution in z is 0");
+    return 0;
+  }
+  
+  //check position
+  if(TMath::Abs(vertex->GetXv()) > gVxMax) {
+    if(fDebugMode)
+      Printf("GetVertex: Event rejected because it has a Vx value of %lf cm (accepted interval: -%lf - %lf)",TMath::Abs(vertex->GetXv()),gVxMax,gVxMax);
+    return 0;
+  }
+  if(TMath::Abs(vertex->GetYv()) > gVyMax)  {
+    if(fDebugMode)
+      Printf("GetVertex: Event rejected because it has a Vy value of %lf cm (accepted interval: -%lf - %lf)",TMath::Abs(vertex->GetYv()),gVyMax,gVyMax);
+    return 0;
+  }
+  if(TMath::Abs(vertex->GetZv()) > gVzMax)  {
+    if(fDebugMode)
+      Printf("GetVertex: Event rejected because it has a Vz value of %lf cm (accepted interval: -%lf - %lf)",TMath::Abs(vertex->GetZv()),gVzMax,gVzMax);
+    return 0;
+  }
+  
+  return vertex;
+}
+
+//________________________________________________________________________
+Bool_t AliProtonAnalysisBase::IsEventTriggered(const AliESDEvent *esd, 
+                                              TriggerMode trigger) {
+  // check if the event was triggered
+  ULong64_t triggerMask = esd->GetTriggerMask();
+  
+  // definitions from p-p.cfg
+  ULong64_t spdFO = (1 << 14);
+  ULong64_t v0left = (1 << 11);
+  ULong64_t v0right = (1 << 12);
+
+  switch (trigger) {
+  case kMB1: {
+    if (triggerMask & spdFO || ((triggerMask & v0left) || (triggerMask & v0right)))
+      return kTRUE;
+    break;
+  }
+  case kMB2: {
+    if (triggerMask & spdFO && ((triggerMask & v0left) || (triggerMask & v0right)))
+      return kTRUE;
+    break;
+  }
+  case kSPDFASTOR: {
+    if (triggerMask & spdFO)
+      return kTRUE;
+    break;
+  }
+  }//switch
+
+  return kFALSE;
+}
+
+//________________________________________________________________________
+TCanvas *AliProtonAnalysisBase::GetListOfCuts() {
+  // return the list of cuts and their cut values for reference
+  TLatex l;
+  l.SetTextAlign(12);
+  l.SetTextSize(0.04);
+
+  TString listOfCuts;
+
+  TCanvas *c = new TCanvas("cListOfCuts","List of cuts",0,0,900,600);
+  c->SetFillColor(10); c->SetHighLightColor(41);
+  c->Divide(3,2);
+
+  c->cd(1)->SetFillColor(10);
+  l.DrawLatex(0.3,0.9,"Analysis details: List of cuts\n\n");
+
+  listOfCuts = "Analysis level: "; listOfCuts += fProtonAnalysisLevel;
+  l.DrawLatex(0.1,0.82,listOfCuts.Data());
+  listOfCuts = "Analysis mode: ";
+  if(fProtonAnalysisMode == kTPC) listOfCuts += "TPC standalone"; 
+  if(fProtonAnalysisMode == kHybrid) listOfCuts += "Hybrid TPC"; 
+  if(fProtonAnalysisMode == kGlobal) listOfCuts += "Global tracking"; 
+  l.DrawLatex(0.1,0.74,listOfCuts.Data());
+  listOfCuts = "Trigger mode: "; 
+  if(fTriggerMode == kMB1) listOfCuts += "Minimum bias 1"; 
+  if(fTriggerMode == kMB2) listOfCuts += "Minimum bias 2"; 
+  if(fTriggerMode == kSPDFASTOR) listOfCuts += "FastOR (SPD)"; 
+  l.DrawLatex(0.1,0.66,listOfCuts.Data());
+  listOfCuts = "PID mode: "; 
+  if(fProtonPIDMode == kBayesian) listOfCuts += "Bayesian PID";
+  if(fProtonPIDMode == kRatio) listOfCuts += "Z = ln((dE/dx)_{exp.}/(dE/dx)_{theor.})"; 
+  if(fProtonPIDMode == kSigma) listOfCuts += "N__{#sigma} area"; 
+  l.DrawLatex(0.1,0.58,listOfCuts.Data());
+  listOfCuts = "Accepted vertex diamond: "; 
+  l.DrawLatex(0.1,0.5,listOfCuts.Data());
+  listOfCuts = "|V_{x}| < "; listOfCuts += fVxMax; listOfCuts += " [cm]";
+  l.DrawLatex(0.6,0.5,listOfCuts.Data());
+  listOfCuts = "|V_{y}| < "; listOfCuts += fVyMax; listOfCuts += " [cm]";
+  l.DrawLatex(0.6,0.4,listOfCuts.Data());
+  listOfCuts = "|V_{z}| < "; listOfCuts += fVzMax; listOfCuts += " [cm]";
+  l.DrawLatex(0.6,0.3,listOfCuts.Data());
+  listOfCuts = "Phase space: "; 
+  l.DrawLatex(0.1,0.2,listOfCuts.Data());
+  if(fAnalysisEtaMode) listOfCuts = "|#eta| < ";  
+  else listOfCuts = "|y| < ";
+  listOfCuts += TMath::Abs(fMinX); 
+  l.DrawLatex(0.35,0.2,listOfCuts.Data());
+  listOfCuts = "N_{bins} = ";
+  listOfCuts += fNBinsX; listOfCuts += " (binning: ";
+  listOfCuts += (fMaxX-fMinX)/fNBinsX; listOfCuts += ")";
+  l.DrawLatex(0.35,0.15,listOfCuts.Data());
+  listOfCuts = ""; listOfCuts += fMinY; listOfCuts += " < P_{T} < ";
+  listOfCuts += fMaxY; listOfCuts += "GeV/c"; 
+  l.DrawLatex(0.35,0.1,listOfCuts.Data());
+  listOfCuts = "N_{bins} = ";
+  listOfCuts += fNBinsY; listOfCuts += " (binning: ";
+  listOfCuts += (fMaxY-fMinY)/fNBinsY; listOfCuts += ")";
+  l.DrawLatex(0.35,0.05,listOfCuts.Data());
+
+  c->cd(2)->SetFillColor(2);
+  l.DrawLatex(0.3,0.9,"ITS related cuts");
+  listOfCuts = "Request a cluster on SPD1: "; 
+  listOfCuts += fPointOnITSLayer1Flag;
+  l.DrawLatex(0.1,0.8,listOfCuts.Data());
+  listOfCuts = "Request a cluster on SPD2: "; 
+  listOfCuts += fPointOnITSLayer2Flag;
+  l.DrawLatex(0.1,0.7,listOfCuts.Data());
+  listOfCuts = "Request a cluster on SDD1: "; 
+  listOfCuts += fPointOnITSLayer3Flag;
+  l.DrawLatex(0.1,0.6,listOfCuts.Data());
+  listOfCuts = "Request a cluster on SDD2: "; 
+  listOfCuts += fPointOnITSLayer4Flag;
+  l.DrawLatex(0.1,0.5,listOfCuts.Data());
+  listOfCuts = "Request a cluster on SSD1: "; 
+  listOfCuts += fPointOnITSLayer5Flag;
+  l.DrawLatex(0.1,0.4,listOfCuts.Data());
+  listOfCuts = "Request a cluster on SSD2: "; 
+  listOfCuts += fPointOnITSLayer6Flag; 
+  l.DrawLatex(0.1,0.3,listOfCuts.Data());  
+  listOfCuts = "Minimum number of ITS clusters: ";
+  if(fMinITSClustersFlag) listOfCuts += fMinITSClusters;
+  else listOfCuts += "Not used";
+  l.DrawLatex(0.1,0.2,listOfCuts.Data());
+  listOfCuts = "Maximum #chi^{2} per ITS cluster: ";
+  if(fMaxChi2PerITSClusterFlag) listOfCuts += fMaxChi2PerITSCluster; 
+  else listOfCuts += "Not used";
+  l.DrawLatex(0.1,0.1,listOfCuts.Data());
+
+  c->cd(3)->SetFillColor(3);
+  l.DrawLatex(0.3,0.9,"TPC related cuts");
+  listOfCuts = "Minimum number of TPC clusters: ";
+  if(fMinTPCClustersFlag) listOfCuts += fMinTPCClusters;
+  else listOfCuts += "Not used";
+  l.DrawLatex(0.1,0.7,listOfCuts.Data());
+  listOfCuts = "Maximum #chi^{2} per TPC cluster: ";
+  if(fMaxChi2PerTPCClusterFlag) listOfCuts += fMaxChi2PerTPCCluster;
+  else listOfCuts += "Not used";
+  l.DrawLatex(0.1,0.5,listOfCuts.Data());
+
+  c->cd(4)->SetFillColor(4);
+  l.DrawLatex(0.3,0.9,"Tracking related cuts");
+  listOfCuts = "Maximum cov11: ";
+  if(fMaxCov11Flag) listOfCuts += fMaxCov11;
+  else listOfCuts += "Not used";
+  l.DrawLatex(0.1,0.75,listOfCuts.Data());
+  listOfCuts = "Maximum cov22: ";
+  if(fMaxCov22Flag) listOfCuts += fMaxCov22;
+  else listOfCuts += "Not used";
+  l.DrawLatex(0.1,0.6,listOfCuts.Data());
+  listOfCuts = "Maximum cov33: ";
+  if(fMaxCov33Flag) listOfCuts += fMaxCov33;
+  else listOfCuts += "Not used";
+  l.DrawLatex(0.1,0.45,listOfCuts.Data());
+  listOfCuts = "Maximum cov44: ";
+  if(fMaxCov44Flag) listOfCuts += fMaxCov44;
+  else listOfCuts += "Not used";
+  l.DrawLatex(0.1,0.3,listOfCuts.Data());
+  listOfCuts = "Maximum cov55: ";
+  if(fMaxCov55Flag) listOfCuts += fMaxCov55;
+  else listOfCuts += "Not used";
+  l.DrawLatex(0.1,0.15,listOfCuts.Data());
+
+  c->cd(5)->SetFillColor(5);
+  l.DrawLatex(0.3,0.9,"DCA related cuts");
+  listOfCuts = "Maximum sigma to vertex: ";
+  if(fMaxSigmaToVertexFlag) listOfCuts += fMaxSigmaToVertex;
+  else listOfCuts += "Not used";
+  l.DrawLatex(0.1,0.8,listOfCuts.Data());
+  listOfCuts = "Maximum sigma to vertex (TPC): ";
+  if(fMaxSigmaToVertexTPCFlag) listOfCuts += fMaxSigmaToVertexTPC;
+  else listOfCuts += "Not used";
+  l.DrawLatex(0.1,0.72,listOfCuts.Data());
+  listOfCuts = "Maximum DCA in xy: ";
+  if(fMaxDCAXYFlag) listOfCuts += fMaxDCAXY;
+  else listOfCuts += "Not used";
+  l.DrawLatex(0.1,0.64,listOfCuts.Data());
+  listOfCuts = "Maximum DCA in z: ";
+  if(fMaxDCAZFlag) listOfCuts += fMaxDCAZ;
+  else listOfCuts += "Not used";
+  l.DrawLatex(0.1,0.56,listOfCuts.Data());
+  listOfCuts = "Maximum DCA in 3D: ";
+  if(fMaxDCA3DFlag) listOfCuts += fMaxDCA3D;
+  else listOfCuts += "Not used";
+  l.DrawLatex(0.1,0.48,listOfCuts.Data());
+  listOfCuts = "Maximum DCA in xy (TPC): ";
+  if(fMaxDCAXYFlag) listOfCuts += fMaxDCAXYTPC;
+  else listOfCuts += "Not used";
+  l.DrawLatex(0.1,0.4,listOfCuts.Data());
+  listOfCuts = "Maximum DCA in z (TPC): ";
+  if(fMaxDCAZFlag) listOfCuts += fMaxDCAZTPC;
+  else listOfCuts += "Not used";
+  l.DrawLatex(0.1,0.32,listOfCuts.Data());
+  listOfCuts = "Maximum DCA in 3D (TPC): ";
+  if(fMaxDCA3DFlag) listOfCuts += fMaxDCA3DTPC;
+  else listOfCuts += "Not used";
+  l.DrawLatex(0.1,0.24,listOfCuts.Data());
+  listOfCuts = "Maximum constrained #chi^{2} (vertex): ";
+  if(fMaxConstrainChi2Flag) listOfCuts += fMaxConstrainChi2;
+  else listOfCuts += "Not used";
+  l.DrawLatex(0.1,0.16,listOfCuts.Data());
+
+  c->cd(6)->SetFillColor(6);
+  l.DrawLatex(0.3,0.9,"Tracking bits related cuts");
+  listOfCuts = "Request the ITS refit bit: "; listOfCuts += fITSRefitFlag;
+  l.DrawLatex(0.1,0.7,listOfCuts.Data());
+  listOfCuts = "Request the TPC refit bit: "; listOfCuts += fTPCRefitFlag; 
+  l.DrawLatex(0.1,0.5,listOfCuts.Data());
+  listOfCuts = "Request the TPC pid bit: "; listOfCuts += fTPCpidFlag;
+  l.DrawLatex(0.1,0.3,listOfCuts.Data());
+  listOfCuts = "Request the ESD pid bit: "; listOfCuts += fESDpidFlag;
+  l.DrawLatex(0.1,0.1,listOfCuts.Data());
+
+  return c;
+}
+
+//________________________________________________________________________
+Bool_t AliProtonAnalysisBase::IsProton(AliESDtrack *track) {
+  //Function that checks if a track is a proton
+  Double_t probability[5];
+  Double_t gPt = 0.0, gP = 0.0;
+  Long64_t fParticleType = 0;
+  //Bayesian approach for the PID
+  if(fProtonPIDMode == kBayesian) {
+    if((fProtonAnalysisMode == kTPC)||(fProtonAnalysisMode == kHybrid)) {
+      AliExternalTrackParam *tpcTrack = (AliExternalTrackParam *)track->GetTPCInnerParam();
+      if(tpcTrack) {
+       gPt = tpcTrack->Pt();
+       gP = tpcTrack->P();
+       track->GetTPCpid(probability);
+      }
+    }//TPC standalone or Hybrid TPC
+    else if(fProtonAnalysisMode == kGlobal) {
+      gPt = track->Pt();
+      gP = track->P();
+      track->GetESDpid(probability);
+    }//Global tracking    
+    
+    Double_t rcc = 0.0;
+    for(Int_t i = 0; i < AliPID::kSPECIES; i++) 
+      rcc += probability[i]*GetParticleFraction(i,gP);
+    if(rcc != 0.0) {
+      Double_t w[5];
+      for(Int_t i = 0; i < AliPID::kSPECIES; i++) 
+       w[i] = probability[i]*GetParticleFraction(i,gP)/rcc;
+      fParticleType = TMath::LocMax(AliPID::kSPECIES,w);
+    }
+    if(fParticleType == 4)
+      return kTRUE;
+  }
+  //Ratio of the measured over the theoretical dE/dx a la STAR
+  else if(fProtonPIDMode == kRatio) {
+    Printf("The kRatio mode is not implemented yet!!!");
+    return kFALSE;
+  }
+  //Definition of an N-sigma area around the dE/dx vs P band
+  else if(fProtonPIDMode == kSigma) {
+    Printf("The kSigma mode is not implemented yet!!!");
+    return kFALSE;
+  }
+
+  return kFALSE;
+}
+
+
+
diff --git a/PWG2/SPECTRA/AliProtonAnalysisBase.h b/PWG2/SPECTRA/AliProtonAnalysisBase.h
new file mode 100644 (file)
index 0000000..91ddfc1
--- /dev/null
@@ -0,0 +1,235 @@
+#ifndef ALIPROTONANALYSISBASE_H
+#define ALIPROTONANALYSISBASE_H
+
+/*  See cxx source for full Copyright notice */
+
+
+/* $Id: AliProtonAnalysisBase.h 31056 2009-02-16 14:31:41Z pchrist $ */
+
+//-------------------------------------------------------------------------
+//                       Class AliProtonAnalysisBase
+//   This is the base class for the baryon (proton) analysis
+//
+//    Origin: Panos Christakoglou | Panos.Christakoglou@cern.ch
+//-------------------------------------------------------------------------
+
+#include "TObject.h"
+#include "TString.h"
+class TF1;
+class TCanvas;
+
+#include "AliPID.h"
+class AliESDEvent;
+class AliESDtrack;
+class AliESDVertex;
+
+class AliProtonAnalysisBase : public TObject {
+ public:
+  enum TriggerMode { kMB1 = 0, kMB2, kSPDFASTOR };
+  enum AnalysisMode { kInvalid = -1, kTPC = 0, kHybrid, kGlobal };
+  enum PIDMode { kBayesian = 0, kRatio, kSigma };
+
+  AliProtonAnalysisBase();
+  virtual ~AliProtonAnalysisBase();
+
+  void SetAnalysisLevel(const char* type) {fProtonAnalysisLevel = type;}
+  void SetAnalysisMode(AnalysisMode analysismode) {fProtonAnalysisMode = analysismode;}
+  void SetEtaMode() {fAnalysisEtaMode = kTRUE;}
+  void SetTriggerMode(TriggerMode triggermode) {fTriggerMode = triggermode;}
+  void SetPIDMode(PIDMode pidmode) {fProtonPIDMode = pidmode;}
+
+  const char *GetAnalysisLevel() {return fProtonAnalysisLevel.Data();}
+  AnalysisMode GetAnalysisMode() {return fProtonAnalysisMode;}
+  Bool_t GetEtaMode() {return fAnalysisEtaMode;}
+  TriggerMode GetTriggerMode() {return fTriggerMode;}
+  PIDMode GetPIDMode() {return fProtonPIDMode;}
+
+  const  AliESDVertex *GetVertex(AliESDEvent *esd,
+                                AnalysisMode mode,
+                                Double_t gVx = 100.,
+                                Double_t gVy = 100.,
+                                Double_t gVz = 100.);
+  void SetAcceptedVertexDiamond(Double_t gVx, Double_t gVy, Double_t gVz) {
+    fVxMax = gVx; fVyMax = gVy; fVzMax = gVz;}
+  Double_t GetVxMax() {return fVxMax;}
+  Double_t GetVyMax() {return fVyMax;}
+  Double_t GetVzMax() {return fVzMax;}
+
+  void SetPhaseSpace(Int_t nBinsX, Double_t gXmin, Double_t gXmax,
+                    Int_t nBinsY, Double_t gYmin, Double_t gYmax) {
+    fNBinsX = nBinsX; fMinX = gXmin; fMaxX = gXmax;
+    fNBinsY = nBinsY; fMinY = gYmin; fMaxY = gYmax;
+  }
+  Int_t GetNBinsX() {return fNBinsX;}
+  Int_t GetNBinsY() {return fNBinsY;}
+  Double_t GetMinX() {return fMinX;}
+  Double_t GetMinY() {return fMinY;}
+  Double_t GetMaxX() {return fMaxX;}
+  Double_t GetMaxY() {return fMaxY;}
+
+  static Bool_t IsEventTriggered(const AliESDEvent *esd,
+                                 TriggerMode trigger = kMB2);
+  Bool_t IsAccepted(AliESDEvent *esd,
+                   const AliESDVertex *vertex, 
+                   AliESDtrack *track);
+  Bool_t IsInPhaseSpace(AliESDtrack *track);
+
+  Float_t GetSigmaToVertex(AliESDtrack* esdTrack) const; 
+  Double_t Rapidity(Double_t Px, Double_t Py, Double_t Pz) const;
+  
+  //Cut functions
+  void    SetPointOnITSLayer1() {fPointOnITSLayer1Flag = kTRUE;}
+  void    SetPointOnITSLayer2() {fPointOnITSLayer2Flag = kTRUE;}
+  void    SetPointOnITSLayer3() {fPointOnITSLayer3Flag = kTRUE;}
+  void    SetPointOnITSLayer4() {fPointOnITSLayer4Flag = kTRUE;}
+  void    SetPointOnITSLayer5() {fPointOnITSLayer5Flag = kTRUE;}
+  void    SetPointOnITSLayer6() {fPointOnITSLayer6Flag = kTRUE;}
+  void    SetMinITSClusters(Int_t minITSClusters) {
+    fMinITSClusters = minITSClusters;
+    fMinITSClustersFlag = kTRUE;
+  }
+  void    SetMaxChi2PerITSCluster(Double_t maxChi2PerITSCluster) {
+    fMaxChi2PerITSCluster = maxChi2PerITSCluster;
+    fMaxChi2PerITSClusterFlag = kTRUE;
+  }
+  void    SetMinTPCClusters(Int_t minTPCClusters) {
+    fMinTPCClusters = minTPCClusters;
+    fMinTPCClustersFlag = kTRUE;
+  }
+  void    SetMaxChi2PerTPCCluster(Double_t maxChi2PerTPCCluster) {
+    fMaxChi2PerTPCCluster = maxChi2PerTPCCluster;
+    fMaxChi2PerTPCClusterFlag = kTRUE;
+  }
+  void    SetMaxCov11(Double_t maxCov11) {
+    fMaxCov11 = maxCov11; fMaxCov11Flag = kTRUE;}
+  void    SetMaxCov22(Double_t maxCov22) {
+    fMaxCov22 = maxCov22; fMaxCov22Flag = kTRUE;}
+  void    SetMaxCov33(Double_t maxCov33) {
+    fMaxCov33 = maxCov33; fMaxCov33Flag = kTRUE;}
+  void    SetMaxCov44(Double_t maxCov44) {
+    fMaxCov44 = maxCov44; fMaxCov44Flag = kTRUE;}
+  void    SetMaxCov55(Double_t maxCov55) {
+    fMaxCov55 = maxCov55; fMaxCov55Flag = kTRUE;}
+  void    SetMaxSigmaToVertex(Double_t maxSigmaToVertex) {
+    fMaxSigmaToVertex = maxSigmaToVertex;
+    fMaxSigmaToVertexFlag = kTRUE;
+  }
+  void    SetMaxSigmaToVertexTPC(Double_t maxSigmaToVertex) {
+    fMaxSigmaToVertexTPC = maxSigmaToVertex;
+    fMaxSigmaToVertexTPCFlag = kTRUE;
+  }
+  void    SetMaxDCAXY(Double_t maxDCAXY) {
+    fMaxDCAXY = maxDCAXY;
+    fMaxDCAXYFlag = kTRUE;
+  }
+  void    SetMaxDCAXYTPC(Double_t maxDCAXY) {
+    fMaxDCAXYTPC = maxDCAXY;
+    fMaxDCAXYTPCFlag = kTRUE;
+  }
+  void    SetMaxDCAZ(Double_t maxDCAZ) {
+    fMaxDCAZ = maxDCAZ;
+    fMaxDCAZFlag = kTRUE;
+  }
+  void    SetMaxDCAZTPC(Double_t maxDCAZ) {
+    fMaxDCAZTPC = maxDCAZ;
+    fMaxDCAZTPCFlag = kTRUE;
+  }
+  void    SetMaxDCA3D(Double_t maxDCA3D) {
+    fMaxDCA3D = maxDCA3D;
+    fMaxDCA3DFlag = kTRUE;
+  }
+  void    SetMaxDCA3DTPC(Double_t maxDCA3D) {
+    fMaxDCA3DTPC = maxDCA3D;
+    fMaxDCA3DTPCFlag = kTRUE;
+  }
+  void    SetMaxConstrainChi2(Double_t maxConstrainChi2) {
+    fMaxConstrainChi2 = maxConstrainChi2;
+    fMaxConstrainChi2Flag = kTRUE;
+  }
+  void    SetITSRefit() {fITSRefitFlag = kTRUE;}
+  void    SetTPCRefit() {fTPCRefitFlag = kTRUE;}
+  void    SetESDpid() {fESDpidFlag = kTRUE;}
+  void    SetTPCpid() {fTPCpidFlag = kTRUE;}
+
+  TCanvas *GetListOfCuts();
+
+  //PID related functions
+  Bool_t IsProton(AliESDtrack *track);
+  void SetPriorProbabilities(Double_t * const partFrac) {
+    for(Int_t i = 0; i < AliPID::kSPECIESN; i++) fPartFrac[i] = partFrac[i];} 
+  void SetPriorProbabilityFunctions(TF1 *const felectron, 
+                                   TF1 *const fmuon, 
+                                   TF1 *const fpion, 
+                                   TF1 *const fkaon, 
+                                   TF1 *const fproton) {
+    fFunctionProbabilityFlag = kTRUE;
+    fElectronFunction = felectron; 
+    fMuonFunction = fmuon; 
+    fPionFunction = fpion;
+    fKaonFunction = fkaon;
+    fProtonFunction = fproton;
+  }
+  Bool_t IsPriorProbabilityFunctionUsed() {return fFunctionProbabilityFlag;}
+  Double_t GetParticleFraction(Int_t i, Double_t p);
+
+  void SetDebugMode() {fDebugMode = kTRUE;}
+  Bool_t GetDebugMode() {return fDebugMode;}
+
+ private:
+  AliProtonAnalysisBase(const AliProtonAnalysisBase&); // Not implemented
+  AliProtonAnalysisBase& operator=(const AliProtonAnalysisBase&); // Not implemented
+
+  TString fProtonAnalysisLevel;//"ESD", "AOD" or "MC"
+  TriggerMode fTriggerMode; //Trigger mode
+  AnalysisMode fProtonAnalysisMode; //Analysis mode: TPC-Hybrid-Global
+  PIDMode fProtonPIDMode; //PID mode: Bayesian-dE/dx ratio-Nsigma areas
+  Bool_t fAnalysisEtaMode; //run the analysis in eta or y
+
+  Double_t fVxMax, fVyMax, fVzMax; //vertex diamond constrain 
+
+  Int_t fNBinsX; //number of bins in y or eta
+  Double_t fMinX, fMaxX; //min & max value of y or eta
+  Int_t fNBinsY;  //number of bins in pT
+  Double_t fMinY, fMaxY; //min & max value of pT
+  
+  //cuts
+  Int_t fMinTPCClusters, fMinITSClusters; //min TPC & ITS clusters
+  Double_t fMaxChi2PerTPCCluster, fMaxChi2PerITSCluster; //max chi2 per TPC & ITS cluster
+  Double_t fMaxCov11, fMaxCov22, fMaxCov33, fMaxCov44, fMaxCov55; //max values of cov. matrix
+  Double_t fMaxSigmaToVertex; //max sigma to vertex cut
+  Double_t fMaxSigmaToVertexTPC; //max sigma to vertex cut
+  Double_t fMaxDCAXY, fMaxDCAXYTPC; //max DCA xy
+  Double_t fMaxDCAZ, fMaxDCAZTPC; //max DCA z
+  Double_t fMaxDCA3D, fMaxDCA3DTPC; //max DCA 3D
+  Double_t fMaxConstrainChi2; //max constrain chi2 - vertex
+  Bool_t fMinTPCClustersFlag, fMinITSClustersFlag; //shows if this cut is used or not
+  Bool_t fMaxChi2PerTPCClusterFlag, fMaxChi2PerITSClusterFlag; //shows if this cut is used or not
+  Bool_t fMaxCov11Flag, fMaxCov22Flag, fMaxCov33Flag, fMaxCov44Flag, fMaxCov55Flag; //shows if this cut is used or not
+  Bool_t fMaxSigmaToVertexFlag; //shows if this cut is used or not
+  Bool_t fMaxSigmaToVertexTPCFlag; //shows if this cut is used or not
+  Bool_t fMaxDCAXYFlag, fMaxDCAXYTPCFlag; //shows if this cut is used or not
+  Bool_t fMaxDCAZFlag, fMaxDCAZTPCFlag; //shows if this cut is used or not
+  Bool_t fMaxDCA3DFlag, fMaxDCA3DTPCFlag; //shows if this cut is used or not
+  Bool_t fMaxConstrainChi2Flag; //shows if this cut is used or not
+  Bool_t fITSRefitFlag, fTPCRefitFlag; //shows if this cut is used or not
+  Bool_t fESDpidFlag, fTPCpidFlag; //shows if this cut is used or not
+  Bool_t fPointOnITSLayer1Flag, fPointOnITSLayer2Flag; //shows if this cut is used or not
+  Bool_t fPointOnITSLayer3Flag, fPointOnITSLayer4Flag; //shows if this cut is used or not
+  Bool_t fPointOnITSLayer5Flag, fPointOnITSLayer6Flag; //shows if this cut is used or not
+  
+  //pid
+  Bool_t fFunctionProbabilityFlag; //flag: kTRUE if functions used
+  Double_t fPartFrac[10]; //prior probabilities
+  TF1  *fElectronFunction; //momentum dependence of the prior probs
+  TF1  *fMuonFunction; //momentum dependence of the prior probs
+  TF1  *fPionFunction; //momentum dependence of the prior probs
+  TF1  *fKaonFunction; //momentum dependence of the prior probs
+  TF1  *fProtonFunction; //momentum dependence of the prior probs
+
+  //Debug
+  Bool_t fDebugMode; //Enable the debug mode
+
+  ClassDef(AliProtonAnalysisBase,0);
+};
+
+#endif
diff --git a/PWG2/configProtonAnalysis.C b/PWG2/configProtonAnalysis.C
new file mode 100644 (file)
index 0000000..17c9aeb
--- /dev/null
@@ -0,0 +1,142 @@
+//__________________________________________________//
+AliProtonAnalysis *GetProtonAnalysisObject(const char* analysisLevel = "ESD", 
+                                          const char* esdAnalysisType = "Hybrid", 
+                                          const char* pidMode = "Bayesian") {
+  //Function to setup the AliProtonAnalysis object and return it
+  AliProtonAnalysisBase *baseAnalysis = GetProtonAnalysisBaseObject(analysisLevel,esdAnalysisType,pidMode);
+
+  AliProtonAnalysis *analysis = new AliProtonAnalysis();
+  analysis->SetBaseAnalysis(baseAnalysis);
+  //if(analysisBase->GetEtaMode()) analysis->SetEtaMode();
+  analysis->InitAnalysisHistograms(baseAnalysis->GetNBinsX(),
+                                  baseAnalysis->GetMinX(),
+                                  baseAnalysis->GetMaxX(),
+                                  baseAnalysis->GetNBinsY(),
+                                  baseAnalysis->GetMinY(),
+                                  baseAnalysis->GetMaxY());
+    
+  return analysis;
+}
+
+//__________________________________________________//
+AliProtonQAAnalysis *GetProtonQAAnalysisObject(const char* analysisLevel = "ESD",
+                                              const char* esdAnalysisType = "Hybrid",
+                                              const char* pidMode = "Bayesian") {
+  //Function to setup the AliProtonQAAnalysis object and return it
+  AliProtonQAAnalysis *analysis = new AliProtonQAAnalysis();
+
+  return analysis;
+}
+
+//__________________________________________________//
+AliProtonAnalysisBase *GetProtonAnalysisBaseObject(const char* analysisLevel = "ESD",
+                                                  const char* esdAnalysisType = "Hybrid",
+                                                  const char* pidMode = "Bayesian") {
+  //Function to setup the AliProtonAnalysisBase object and return it
+  AliProtonAnalysisBase *baseAnalysis = new AliProtonAnalysisBase();
+  //baseAnalysis->SetDebugMode();
+  baseAnalysis->SetAnalysisLevel(analysisLevel);
+  if(analysisLevel == "ESD") {  
+    baseAnalysis->SetTriggerMode(AliProtonAnalysisBase::kMB2);
+    switch(esdAnalysisType) {
+    case "TPC":
+      baseAnalysis->SetAnalysisMode(AliProtonAnalysisBase::kTPC);
+      baseAnalysis->SetPhaseSpace(10, -0.5, 0.5, 16, 0.5, 0.9);
+      baseAnalysis->SetTPCpid();
+      baseAnalysis->SetMinTPCClusters(100);
+      baseAnalysis->SetMaxChi2PerTPCCluster(2.2);
+      baseAnalysis->SetMaxCov11(0.5);
+      baseAnalysis->SetMaxCov22(0.5);
+      baseAnalysis->SetMaxCov33(0.5);
+      baseAnalysis->SetMaxCov44(0.5);
+      baseAnalysis->SetMaxCov55(0.5);
+      baseAnalysis->SetMaxSigmaToVertexTPC(2.0);
+      //baseAnalysis->SetMaxDCAXYTPC(1.5);
+      //baseAnalysis->SetMaxDCAZTPC(1.5);
+      break;
+    case "Hybrid":
+      baseAnalysis->SetAnalysisMode(AliProtonAnalysisBase::kHybrid);
+      baseAnalysis->SetPhaseSpace(10, -0.5, 0.5, 16, 0.5, 0.9);
+      baseAnalysis->SetTPCpid();
+      baseAnalysis->SetMinTPCClusters(110);
+      baseAnalysis->SetMaxChi2PerTPCCluster(2.2);
+      baseAnalysis->SetMaxCov11(0.5);
+      baseAnalysis->SetMaxCov22(0.5);
+      baseAnalysis->SetMaxCov33(0.5);
+      baseAnalysis->SetMaxCov44(0.5);
+      baseAnalysis->SetMaxCov55(0.5);
+      baseAnalysis->SetMaxSigmaToVertex(2.0);
+      /*baseAnalysis->SetMaxDCAXY(1.5);
+       baseAnalysis->SetMaxDCAZ(1.5);*/
+      baseAnalysis->SetPointOnITSLayer6();
+      baseAnalysis->SetPointOnITSLayer5();
+      //baseAnalysis->SetPointOnITSLayer4();
+      //baseAnalysis->SetPointOnITSLayer3();
+      baseAnalysis->SetPointOnITSLayer2();
+      baseAnalysis->SetPointOnITSLayer1();
+      baseAnalysis->SetMinITSClusters(5);
+      break;
+    case "Global":
+      baseAnalysis->SetAnalysisMode(AliProtonAnalysisBase::kGlobal);
+      baseAnalysis->SetPhaseSpace(20, -1.0, 1.0, 48, 0.3, 1.5);
+      baseAnalysis->SetMinTPCClusters(110);
+      baseAnalysis->SetMaxChi2PerTPCCluster(2.2);
+      baseAnalysis->SetMaxCov11(0.5);
+      baseAnalysis->SetMaxCov22(0.5);
+      baseAnalysis->SetMaxCov33(0.5);
+      baseAnalysis->SetMaxCov44(0.5);
+      baseAnalysis->SetMaxCov55(0.5);
+      baseAnalysis->SetMaxSigmaToVertex(2.0);
+      //baseAnalysis->SetMaxDCAXY(2.0);
+      //baseAnalysis->SetMaxDCAZ(2.0);
+      baseAnalysis->SetTPCRefit();
+      baseAnalysis->SetPointOnITSLayer1();
+      baseAnalysis->SetPointOnITSLayer2();
+      //baseAnalysis->SetPointOnITSLayer3();
+      //baseAnalysis->SetPointOnITSLayer4();
+      baseAnalysis->SetPointOnITSLayer5();
+      baseAnalysis->SetPointOnITSLayer6();
+      baseAnalysis->SetMinITSClusters(5);
+      baseAnalysis->SetITSRefit();
+      baseAnalysis->SetESDpid();
+      break;
+    default:
+      break;
+    }
+    baseAnalysis->SetAcceptedVertexDiamond(5.,5.,15.);
+    baseAnalysis->SetEtaMode();
+    switch(pidMode) {
+    case "Bayesian":
+      baseAnalysis->SetPIDMode(AliProtonAnalysisBase::kBayesian);
+      //Momentum dependent priors
+      /*TFile *f = TFile::Open("PriorProb/PriorProbabilities.root ");
+       TF1 *fitElectrons = (TF1 *)f->Get("fitElectrons");
+       TF1 *fitMuons = (TF1 *)f->Get("fitMuons");
+       TF1 *fitPions = (TF1 *)f->Get("fitPions");
+       TF1 *fitKaons = (TF1 *)f->Get("fitKaons");
+       TF1 *fitProtons = (TF1 *)f->Get("fitProtons");
+       baseAnalysis->SetPriorProbabilityFunctions(fitElectrons,
+       fitMuons,
+       fitPions,
+       fitKaons,
+       fitProtons);*/
+      //Fixed prior probabilities
+      Double_t partFrac[5] = {0.01, 0.01, 0.85, 0.10, 0.05};
+      if(!baseAnalysis->IsPriorProbabilityFunctionUsed())
+       baseAnalysis->SetPriorProbabilities(partFrac);
+      break;
+    case "Ratio":
+      baseAnalysis->SetPIDMode(AliProtonAnalysisBase::kRatio);
+      break;
+    case "Sigma":
+      baseAnalysis->SetPIDMode(AliProtonAnalysisBase::kSigma);
+      break;
+    default:
+      break;
+    }//PID mode
+  }//ESD
+  if(analysisLevel == "MC") 
+    baseAnalysis->SetPhaseSpace(56, -1.0, 1.0, 16, 0.1, 1.5);
+
+  return baseAnalysis;
+}
index 45bd87b..3e7f0e1 100644 (file)
@@ -1,27 +1,34 @@
-void runProtonAnalysis(const char* esdAnalysisType = "TPC",) {
+void runProtonAnalysis(const char* esdAnalysisType = "Hybrid",
+                      const char* pidMode = "Bayesian") {
   //Macro to run the proton analysis tested for local, proof & GRID.
-  //Local: Takes three arguments, the analysis mode, the type of the ESD 
-  //       analysis and the path where the tag and ESD or AOD files reside.
-  //Interactive: Takes three arguments, the analysis mode, the type of the ESD 
-  //             analysis and the name of the collection of tag files.
-  //Batch: Takes three arguments, the analysis mode, the type of the ESD 
-  //       analysis and the name of the collection file with the event list 
-  //       for each file.
-  //Proof: Takes four arguments, the analysis mode, the number of events,
-  //       the dataset name and the analysis type in case of ESD.
-  
+  //Local: Takes four arguments, the analysis mode, the type of the ESD 
+  //       analysis, the PID mode and the path where the tag and ESD or 
+  //       AOD files reside.
+  //Interactive: Takes four arguments, the analysis mode, the type of the ESD 
+  //             analysis, the PID mode and the name of the collection of tag 
+  //             files.
+  //Batch: Takes four arguments, the analysis mode, the type of the ESD 
+  //       analysis, the PID mode and the name of the collection file with 
+  //       the event list for each file.
+  //Proof: Takes five arguments, the analysis level, the analysis mode in 
+  //       case of ESD, the PID mode, the number of events and the dataset 
+  //       name and .  
   //Analysis mode can be: "MC", "ESD", "AOD"
   //ESD analysis type can be one of the three: "TPC", "Hybrid", "Global"
+  //PID mode can be one of the three: "Bayesian" (standard Bayesian approach) 
+  //   "Ratio" (ratio of measured over expected/theoretical dE/dx a la STAR) 
+  //   "Sigma" (N-sigma area around the fitted dE/dx vs P band)
   TStopwatch timer;
   timer.Start();
   
-  //runLocal("ESD",esdAnalysisType,"/home/pchrist/ALICE/Alien/Tutorial/November2007/Tags");
-  //runInteractive("ESD",esdAnalysisType,"tag.xml");
-  //runBatch("ESD",esdAnalysisType,"wn.xml");  
-  runProof("MC",
-          200000,
-          "/COMMON/COMMON/LHC08c11_10TeV_0.5T",
-          esdAnalysisType);
+  runLocal("ESD",
+          esdAnalysisType,
+          pidMode,
+          "/home/pchrist/ALICE/Baryons/QA/Local");
+  //runInteractive("ESD",esdAnalysisType,pidMode,"tag.xml");
+  //runBatch("ESD",esdAnalysisType,pidMode,"wn.xml");  
+  /*runProof("ESD",esdAnalysisType,pidMode
+    200000,"/COMMON/COMMON/LHC08c11_10TeV_0.5T",);*/
   
   timer.Stop();
   timer.Print();
@@ -30,10 +37,8 @@ void runProtonAnalysis(const char* esdAnalysisType = "TPC",) {
 //_________________________________________________//
 void runLocal(const char* mode = "ESD",
              const char* analysisType = 0x0,
+             const char* pidMode = 0x0,
              const char* path = "/home/pchrist/ALICE/Alien/Tutorial/November2007/Tags") {
-  TStopwatch timer;
-  timer.Start();
-
   TString smode = mode;
   TString outputFilename = "Protons."; outputFilename += mode;
   if(analysisType) {
@@ -55,8 +60,8 @@ void runLocal(const char* mode = "ESD",
   gSystem->Load("libANALYSIS.so");
   setupPar("ANALYSISalice");
   gSystem->Load("libANALYSISalice.so");
-  setupPar->UploadPackage("CORRFW.par");
-  gSystem->EnablePackage("CORRFW");
+  setupPar("CORRFW");
+  gSystem->Load("libCORRFW.so");
   setupPar("PWG2spectra");
   gSystem->Load("libPWG2spectra.so");
   //____________________________________________________//  
@@ -75,7 +80,10 @@ void runLocal(const char* mode = "ESD",
   chain->SetBranchStatus("*Calo*",0);
 
   //____________________________________________//
-  gROOT->LoadMacro("AliAnalysisTaskProtons.cxx++");
+  gROOT->LoadMacro("configProtonAnalysis.C");
+  AliProtonAnalysis *analysis = GetProtonAnalysisObject(mode,
+                                                       analysisType,
+                                                       pidMode);
   //____________________________________________//
   // Make the analysis manager
   AliAnalysisManager *mgr = new AliAnalysisManager("TestManager");
@@ -85,45 +93,20 @@ void runLocal(const char* mode = "ESD",
     AliMCEventHandler *mc = new AliMCEventHandler();
     mgr->SetMCtruthEventHandler(mc);
   }
-
+  
   //____________________________________________//
-  // 1st Proton task
+  //Create the proton task
   AliAnalysisTaskProtons *taskProtons = new AliAnalysisTaskProtons("TaskProtons");
-  taskProtons->SetType(mode);
-  taskProtons->SetTriggerMode(AliAnalysisTaskProtons::kMB2);
-  switch(analysisType) {
-  case "TPC":
-    taskProtons->SetAnalysisMode(AliAnalysisTaskProtons::kTPC);
-    break;
-  case "Hybrid":
-    taskProtons->SetAnalysisMode(AliAnalysisTaskProtons::kHybrid);
-    break;
-  case "Global":
-    taskProtons->SetAnalysisMode(AliAnalysisTaskProtons::kGlobal);
-    break;
-  default:
-    break;
-  }
-  taskProtons->SetAcceptedVertexDiamond(5.,5.,15.);
-  //Momentum dependent priors
-  /*TFile *f = TFile::Open("PriorProb/PriorProbabilities.root ");
-  TF1 *fitElectrons = (TF1 *)f->Get("fitElectrons");
-  TF1 *fitMuons = (TF1 *)f->Get("fitMuons");
-  TF1 *fitPions = (TF1 *)f->Get("fitPions");
-  TF1 *fitKaons = (TF1 *)f->Get("fitKaons");
-  TF1 *fitProtons = (TF1 *)f->Get("fitProtons");
-  taskProtons->SetPriorProbabilityFunctions(fitElectrons,
-                                           fitMuons,
-                                           fitPions,
-                                           fitKaons,
-                                           fitProtons);*/
+  taskProtons->SetAnalysisObject(analysis);
   mgr->AddTask(taskProtons);
 
-  // Create containers for input/output                                                                              
-  AliAnalysisDataContainer *cinput1 = mgr->GetCommonInputContainer();
-  AliAnalysisDataContainer *coutput1 = mgr->CreateContainer("outputList1",
+  // Create containers for input/output
+  AliAnalysisDataContainer *cinput1 = mgr->CreateContainer("dataChain",
+                                                          TChain::Class(),
+                                                          AliAnalysisManager::kInputContainer);
+  AliAnalysisDataContainer *coutput1 = mgr->CreateContainer("outputList",
                                                             TList::Class(),
-                                                           AliAnalysisManager::kOutputCont
+                                                           AliAnalysisManager::kOutputContainer,
                                                             outputFilename.Data());
 
   //____________________________________________//
@@ -132,17 +115,13 @@ void runLocal(const char* mode = "ESD",
   if (!mgr->InitAnalysis()) return;
   mgr->PrintStatus();
   mgr->StartAnalysis("local",chain);
-
-  timer.Stop();
-  timer.Print();
 }
 
 //_________________________________________________//
 void runInteractive(const char* mode = "ESD",
                    const char* analysisType = 0x0,
+                   const char* pidMode = 0x0,
                    const char* collectionName = "tag.xml") {
-  TStopwatch timer;
-  timer.Start();
   gSystem->Load("libProofPlayer.so");
 
   TString smode = mode;
@@ -169,8 +148,8 @@ void runInteractive(const char* mode = "ESD",
   gSystem->Load("libANALYSIS.so");
   setupPar("ANALYSISalice");
   gSystem->Load("libANALYSISalice.so");
-  setupPar->UploadPackage("CORRFW.par");
-  gSystem->EnablePackage("CORRFW");
+  setupPar("CORRFW");
+  gSystem->Load("libCORRFW.so");
   setupPar("PWG2spectra");
   gSystem->Load("libPWG2spectra.so");
   //____________________________________________________//  
@@ -192,7 +171,10 @@ void runInteractive(const char* mode = "ESD",
   chain->SetBranchStatus("*Calo*",0);
   
   //____________________________________________//
-  gROOT->LoadMacro("AliAnalysisTaskProtons.cxx++");
+  gROOT->LoadMacro("configProtonAnalysis.C");
+  AliProtonAnalysis *analysis = GetProtonAnalysisObject(mode,
+                                                       analysisType,
+                                                       pidMode);
   //____________________________________________//
   // Make the analysis manager
   AliAnalysisManager *mgr = new AliAnalysisManager("TestManager");
@@ -204,43 +186,18 @@ void runInteractive(const char* mode = "ESD",
   }
 
   //____________________________________________//
-  // 1st Proton task
+  //Create the proton task
   AliAnalysisTaskProtons *taskProtons = new AliAnalysisTaskProtons("TaskProtons");
-  taskProtons->SetType(mode);
-  taskProtons->SetTriggerMode(AliAnalysisTaskProtons::kMB2);
-  switch(analysisType) {
-  case "TPC":
-    taskProtons->SetAnalysisMode(AliAnalysisTaskProtons::kTPC);
-    break;
-  case "Hybrid":
-    taskProtons->SetAnalysisMode(AliAnalysisTaskProtons::kHybrid);
-    break;
-  case "Global":
-    taskProtons->SetAnalysisMode(AliAnalysisTaskProtons::kGlobal);
-    break;
-  default:
-    break;
-  }
-  taskProtons->SetAcceptedVertexDiamond(5.,5.,15.);
-  //Momentum dependent priors
-  /*TFile *f = TFile::Open("PriorProb/PriorProbabilities.root ");
-  TF1 *fitElectrons = (TF1 *)f->Get("fitElectrons");
-  TF1 *fitMuons = (TF1 *)f->Get("fitMuons");
-  TF1 *fitPions = (TF1 *)f->Get("fitPions");
-  TF1 *fitKaons = (TF1 *)f->Get("fitKaons");
-  TF1 *fitProtons = (TF1 *)f->Get("fitProtons");
-  taskProtons->SetPriorProbabilityFunctions(fitElectrons,
-                                           fitMuons,
-                                           fitPions,
-                                           fitKaons,
-                                           fitProtons);*/
+  taskProtons->SetAnalysisObject(analysis);
   mgr->AddTask(taskProtons);
 
-  // Create containers for input/output                                                                               
-  AliAnalysisDataContainer *cinput1 = mgr->GetCommonInputContainer();
-  AliAnalysisDataContainer *coutput1 = mgr->CreateContainer("outputList1",
+  // Create containers for input/output
+  AliAnalysisDataContainer *cinput1 = mgr->CreateContainer("dataChain",
+                                                          TChain::Class(),
+                                                          AliAnalysisManager::kInputContainer);
+  AliAnalysisDataContainer *coutput1 = mgr->CreateContainer("outputList",
                                                             TList::Class(),
-                                                           AliAnalysisManager::kOutputCont
+                                                           AliAnalysisManager::kOutputContainer,
                                                             outputFilename.Data());
   
   //____________________________________________//
@@ -249,18 +206,13 @@ void runInteractive(const char* mode = "ESD",
   if (!mgr->InitAnalysis()) return;
   mgr->PrintStatus();
   mgr->StartAnalysis("local",chain);
-
-  timer.Stop();
-  timer.Print();
 }
 
 //_________________________________________________//
 void runBatch(const char* mode = "ESD",
              const char* analysisType = 0x0,
+             const char* pidMode = 0x0,
              const char *collectionfile = "wn.xml") {
-  TStopwatch timer;
-  timer.Start();
-
   TString smode = mode;
   TString outputFilename = "Protons."; outputFilename += mode;
   if(analysisType) {
@@ -286,8 +238,8 @@ void runBatch(const char* mode = "ESD",
   gSystem->Load("libANALYSIS.so");
   setupPar("ANALYSISalice");
   gSystem->Load("libANALYSISalice.so");
-  setupPar->UploadPackage("CORRFW.par");
-  gSystem->EnablePackage("CORRFW");
+  setupPar("CORRFW");
+  gSystem->Load("libCORRFW.so");
   setupPar("PWG2spectra");
   gSystem->Load("libPWG2spectra.so");
   //____________________________________________________//  
@@ -300,7 +252,10 @@ void runBatch(const char* mode = "ESD",
   chain->SetBranchStatus("*Calo*",0);
 
   //____________________________________________//
-  gROOT->LoadMacro("AliAnalysisTaskProtons.cxx++");
+  gROOT->LoadMacro("configProtonAnalysis.C");
+  AliProtonAnalysis *analysis = GetProtonAnalysisObject(mode,
+                                                       analysisType,
+                                                       pidMode);
   //____________________________________________//
   // Make the analysis manager
   AliAnalysisManager *mgr = new AliAnalysisManager("TestManager");
@@ -312,42 +267,18 @@ void runBatch(const char* mode = "ESD",
   }
   
   //____________________________________________//
-  // 1st Proton task
+  //Create the proton task
   AliAnalysisTaskProtons *taskProtons = new AliAnalysisTaskProtons("TaskProtons");
-  taskProtons->SetType(mode);
-  taskProtons->SetTriggerMode(AliAnalysisTaskProtons::kMB2);
-  switch(analysisType) {
-  case "TPC":
-    taskProtons->SetAnalysisMode(AliAnalysisTaskProtons::kTPC);
-    break;
-  case "Hybrid":
-    taskProtons->SetAnalysisMode(AliAnalysisTaskProtons::kHybrid);
-    break;
-  case "Global":
-    taskProtons->SetAnalysisMode(AliAnalysisTaskProtons::kGlobal);
-    break;
-  default:
-    break;
-  }
-  taskProtons->SetAcceptedVertexDiamond(5.,5.,15.);
-  //Momentum dependent priors
-  /*TFile *f = TFile::Open("PriorProb/PriorProbabilities.root ");
-  TF1 *fitElectrons = (TF1 *)f->Get("fitElectrons");
-  TF1 *fitMuons = (TF1 *)f->Get("fitMuons");
-  TF1 *fitPions = (TF1 *)f->Get("fitPions");
-  TF1 *fitKaons = (TF1 *)f->Get("fitKaons");
-  TF1 *fitProtons = (TF1 *)f->Get("fitProtons");
-  taskProtons->SetPriorProbabilityFunctions(fitElectrons,
-                                           fitMuons,
-                                           fitPions,
-                                           fitKaons,
-                                           fitProtons);*/
+  taskProtons->SetAnalysisObject(analysis);
   mgr->AddTask(taskProtons);
 
-  // Create containers for input/output                                                                               
-  AliAnalysisDataContainer *cinput1 = mgr->GetCommonInputContainer();
-  AliAnalysisDataContainer *coutput1 = mgr->CreateContainer("outputList1",
-                                                            TList::Class(),AliAnalysisManager::kOutputCont
+  // Create containers for input/output
+  AliAnalysisDataContainer *cinput1 = mgr->CreateContainer("dataChain",
+                                                          TChain::Class(),
+                                                          AliAnalysisManager::kInputContainer);
+  AliAnalysisDataContainer *coutput1 = mgr->CreateContainer("outputList",
+                                                            TList::Class(),
+                                                           AliAnalysisManager::kOutputContainer,
                                                             outputFilename.Data());
 
   //____________________________________________//
@@ -356,19 +287,14 @@ void runBatch(const char* mode = "ESD",
   if (!mgr->InitAnalysis()) return;
   mgr->PrintStatus();
   mgr->StartAnalysis("grid",chain);
-
-  timer.Stop();
-  timer.Print();
 }
 
 //_________________________________________________//
 void runProof(const char* mode = "ESD",
+             const char* analysisType = 0x0,
+             const char* pidMode = 0x0,
              Int_t stats = 0, 
-             const char* dataset = 0x0,
-             const char* analysisType = 0x0) {
-  TStopwatch timer;
-  timer.Start();
-  
+             const char* dataset = 0x0) {  
   TString smode = mode;
   TString outputFilename = "Protons."; outputFilename += mode;
   if(analysisType) {
@@ -397,7 +323,10 @@ void runProof(const char* mode = "ESD",
   gProof->EnablePackage("PWG2spectra");
   
   //____________________________________________//
-  gProof->Load("AliAnalysisTaskProtons.cxx++");
+  gROOT->LoadMacro("configProtonAnalysis.C");
+  AliProtonAnalysis *analysis = GetProtonAnalysisObject(mode,
+                                                       analysisType,
+                                                       pidMode);
   //____________________________________________//
 
   //____________________________________________//
@@ -410,44 +339,19 @@ void runProof(const char* mode = "ESD",
     mgr->SetMCtruthEventHandler(mc);
   }
   //____________________________________________//
-  // 1st Proton task
+  //Create the proton task
   AliAnalysisTaskProtons *taskProtons = new AliAnalysisTaskProtons("TaskProtons");
-  taskProtons->SetType(mode);
-  taskProtons->SetTriggerMode(AliAnalysisTaskProtons::kMB2);
-  switch(analysisType) {
-  case "TPC":
-    taskProtons->SetAnalysisMode(AliAnalysisTaskProtons::kTPC);
-    break;
-  case "Hybrid":
-    taskProtons->SetAnalysisMode(AliAnalysisTaskProtons::kHybrid);
-    break;
-  case "Global":
-    taskProtons->SetAnalysisMode(AliAnalysisTaskProtons::kGlobal);
-    break;
-  default:
-    break;
-  }
-  taskProtons->SetAcceptedVertexDiamond(5.,5.,15.);
-  //Momentum dependent priors
-  /*TFile *f = TFile::Open("PriorProb/PriorProbabilities.root ");
-  TF1 *fitElectrons = (TF1 *)f->Get("fitElectrons");
-  TF1 *fitMuons = (TF1 *)f->Get("fitMuons");
-  TF1 *fitPions = (TF1 *)f->Get("fitPions");
-  TF1 *fitKaons = (TF1 *)f->Get("fitKaons");
-  TF1 *fitProtons = (TF1 *)f->Get("fitProtons");
-  taskProtons->SetPriorProbabilityFunctions(fitElectrons,
-                                           fitMuons,
-                                           fitPions,
-                                           fitKaons,
-                                           fitProtons);*/
+  taskProtons->SetAnalysisObject(analysis);
   mgr->AddTask(taskProtons);
 
   // Create containers for input/output
-  AliAnalysisDataContainer *cinput1 = mgr->GetCommonInputContainer();
-  AliAnalysisDataContainer *coutput1 = mgr->CreateContainer("outputList1", 
-                                                           TList::Class(),
+  AliAnalysisDataContainer *cinput1 = mgr->CreateContainer("dataChain",
+                                                          TChain::Class(),
+                                                          AliAnalysisManager::kInputContainer);
+  AliAnalysisDataContainer *coutput1 = mgr->CreateContainer("outputList",
+                                                            TList::Class(),
                                                            AliAnalysisManager::kOutputContainer,
-                                                           outputFilename.Data());
+                                                            outputFilename.Data());
 
   //____________________________________________//
   mgr->ConnectInput(taskProtons,0,cinput1);
@@ -467,9 +371,6 @@ void runProof(const char* mode = "ESD",
 
     mgr->StartAnalysis("proof",chain);
   }
-
-  timer.Stop();
-  timer.Print();
 }
 
 //_________________________________________________//