added pid example macros and QA task for the PMD
authorsnelling <snelling@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 30 Mar 2011 13:27:11 +0000 (13:27 +0000)
committersnelling <snelling@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 30 Mar 2011 13:27:11 +0000 (13:27 +0000)
PWG2/CMakelibPWG2flowTasks.pkg
PWG2/FLOW/AliFlowTasks/AliAnalysisTaskQAPmdflow.cxx [new file with mode: 0644]
PWG2/FLOW/AliFlowTasks/AliAnalysisTaskQAPmdflow.h [new file with mode: 0644]
PWG2/FLOW/macros/AddTaskFlowCentralityPID.C [new file with mode: 0644]
PWG2/FLOW/macros/runFlowTaskCentralityPIDTrain.C [new file with mode: 0644]
PWG2/PWG2flowTasksLinkDef.h

index 171152d..c612ef9 100644 (file)
@@ -25,7 +25,7 @@
 # SHLIBS - Shared Libraries and objects for linking (Executables only)           #
 #--------------------------------------------------------------------------------#
 
-set ( SRCS  FLOW/AliFlowTasks/AliFlowEventSimpleMaker.cxx FLOW/AliFlowTasks/AliFlowEvent.cxx FLOW/AliFlowTasks/AliFlowEventCuts.cxx FLOW/AliFlowTasks/AliFlowTrack.cxx FLOW/AliFlowTasks/AliFlowCandidateTrack.cxx FLOW/AliFlowTasks/AliFlowTrackCuts.cxx FLOW/AliFlowTasks/AliAnalysisTaskScalarProduct.cxx FLOW/AliFlowTasks/AliAnalysisTaskMCEventPlane.cxx FLOW/AliFlowTasks/AliAnalysisTaskLYZEventPlane.cxx FLOW/AliFlowTasks/AliAnalysisTaskCumulants.cxx FLOW/AliFlowTasks/AliAnalysisTaskQCumulants.cxx FLOW/AliFlowTasks/AliAnalysisTaskLeeYangZeros.cxx FLOW/AliFlowTasks/AliAnalysisTaskFittingQDistribution.cxx FLOW/AliFlowTasks/AliAnalysisTaskFlowEvent.cxx FLOW/AliFlowTasks/AliAnalysisTaskFlowK0Candidates.cxx FLOW/AliFlowTasks/AliAnalysisTaskFlowEventforRP.cxx FLOW/AliFlowTasks/AliAnalysisTaskMixedHarmonics.cxx FLOW/AliFlowTasks/AliAnalysisTaskNestedLoops.cxx FLOW/AliFlowTasks/AliAnalysisTaskQAflow.cxx FLOW/AliFlowTasks/AliAnalysisTaskPIDflowQA.cxx)
+set ( SRCS  FLOW/AliFlowTasks/AliFlowEventSimpleMaker.cxx FLOW/AliFlowTasks/AliFlowEvent.cxx FLOW/AliFlowTasks/AliFlowEventCuts.cxx FLOW/AliFlowTasks/AliFlowTrack.cxx FLOW/AliFlowTasks/AliFlowCandidateTrack.cxx FLOW/AliFlowTasks/AliFlowTrackCuts.cxx FLOW/AliFlowTasks/AliAnalysisTaskScalarProduct.cxx FLOW/AliFlowTasks/AliAnalysisTaskMCEventPlane.cxx FLOW/AliFlowTasks/AliAnalysisTaskLYZEventPlane.cxx FLOW/AliFlowTasks/AliAnalysisTaskCumulants.cxx FLOW/AliFlowTasks/AliAnalysisTaskQCumulants.cxx FLOW/AliFlowTasks/AliAnalysisTaskLeeYangZeros.cxx FLOW/AliFlowTasks/AliAnalysisTaskFittingQDistribution.cxx FLOW/AliFlowTasks/AliAnalysisTaskFlowEvent.cxx FLOW/AliFlowTasks/AliAnalysisTaskFlowK0Candidates.cxx FLOW/AliFlowTasks/AliAnalysisTaskFlowEventforRP.cxx FLOW/AliFlowTasks/AliAnalysisTaskMixedHarmonics.cxx FLOW/AliFlowTasks/AliAnalysisTaskNestedLoops.cxx FLOW/AliFlowTasks/AliAnalysisTaskQAflow.cxx FLOW/AliFlowTasks/AliAnalysisTaskPIDflowQA.cxx FLOW/AliFlowTasks/AliAnalysisTaskQAPmdflow.cxx)
 
 string ( REPLACE ".cxx" ".h" HDRS "${SRCS}" )
 
diff --git a/PWG2/FLOW/AliFlowTasks/AliAnalysisTaskQAPmdflow.cxx b/PWG2/FLOW/AliFlowTasks/AliAnalysisTaskQAPmdflow.cxx
new file mode 100644 (file)
index 0000000..8698238
--- /dev/null
@@ -0,0 +1,170 @@
+#include "TMath.h"
+#include "TH1D.h"
+#include "TH2D.h"
+#include "TSeqCollection.h"
+#include "TObjArray.h"
+#include "TObjArray.h"
+#include "TChain.h"
+#include "TMCProcess.h"
+#include "TLorentzVector.h"
+#include "TDirectory.h"
+#include "TROOT.h"
+#include "TNtuple.h"
+
+#include "AliLog.h"
+#include "AliVParticle.h"
+#include "AliMCParticle.h"
+#include "AliStack.h"
+#include "AliVEvent.h"
+#include "AliESDEvent.h"
+#include "AliMCEvent.h"
+#include "AliESDVZERO.h"
+#include "AliESDZDC.h"
+#include "AliESDtrack.h"
+#include "AliFlowTrackCuts.h"
+#include "AliFlowEventCuts.h"
+#include "AliMultiplicity.h"
+#include "AliESDtrackCuts.h"
+#include "AliVertex.h"
+#include "AliFlowEventSimple.h"
+#include "AliFlowEvent.h"
+#include "AliFlowVector.h"
+#include "AliESDPmdTrack.h"
+#include "AliAnalysisTaskSE.h"
+#include "AliAnalysisManager.h"
+#include "AliInputEventHandler.h"
+
+#include "AliAnalysisTaskQAPmdflow.h"
+
+ClassImp(AliAnalysisTaskQAPmdflow)
+
+//________________________________________________________________________
+AliAnalysisTaskQAPmdflow::AliAnalysisTaskQAPmdflow()
+  : AliAnalysisTaskSE(),
+    fOutput(NULL),
+    fEventCuts(NULL),
+    fRPTrackCuts(NULL),
+    fPOITrackCuts(NULL)
+{
+  // Default constructor
+}
+
+//________________________________________________________________________
+AliAnalysisTaskQAPmdflow::AliAnalysisTaskQAPmdflow(const char* name)
+  : AliAnalysisTaskSE(name),
+    fOutput(NULL),
+    fEventCuts(NULL),
+    fRPTrackCuts(NULL),
+    fPOITrackCuts(NULL)
+{
+  // Constructor
+  DefineOutput(1, TObjArray::Class());
+}
+
+//________________________________________________________________________
+void AliAnalysisTaskQAPmdflow::UserCreateOutputObjects()
+{
+  // Called once at the beginning
+  fOutput=new TObjArray();
+  
+  //define histograms
+  TH1F* histB = new TH1F("PMD ADC cutB","PMD ADC cutB",500,0,10000);
+  fOutput->Add(histB); 
+
+  TH1F* histA = new TH1F("PMD ADC cutA","PMD ADC cutA",500,0,10000);
+  fOutput->Add(histA); 
+
+  TH1F* histCelB = new TH1F("PMD ncell CutB", "PMD ncell CutB",100,0,100);
+  fOutput->Add(histCelB); 
+  TH1F* histCelA = new TH1F("PMD ncell CutA", "PMD ncell CutA",100,0,100);
+  fOutput->Add(histCelA);
+  
+  //post data here as it doesn't change anyway (the pointer to list anyway)
+
+  PostData(1, fOutput);
+}
+
+//________________________________________________________________________
+void AliAnalysisTaskQAPmdflow::UserExec(Option_t *)
+{
+  //get teh input data
+  AliESDEvent* event = dynamic_cast<AliESDEvent*>(InputEvent());
+  if (!event)
+    {
+      AliFatal("no ESD event");
+      return;
+    }
+
+  fRPTrackCuts->SetEvent(event);
+
+  AliFlowTrackCuts::trackParameterType sourceRP = fRPTrackCuts->GetParamType();
+  AliFlowTrackCuts::trackParameterType sourcePOI = fPOITrackCuts->GetParamType();
+  Bool_t PmdTrRp = kFALSE;
+  Bool_t PmdTrPoi = kFALSE;
+  
+  if(sourcePOI == 4) PmdTrPoi = kTRUE;
+  if(sourceRP  == 4) PmdTrRp = kTRUE;
+  if((!PmdTrPoi) && (!PmdTrRp))
+    {
+      printf("Error : PMD track is not used as POI or RP");
+      return;
+    }
+  
+  TH1F* hPmdAdcB = static_cast<TH1F*>(fOutput->At(0));
+  TH1F* hPmdAdcA = static_cast<TH1F*>(fOutput->At(1));
+  TH1F* hPmdNcelB = static_cast<TH1F*>(fOutput->At(2));
+  TH1F* hPmdNcelA = static_cast<TH1F*>(fOutput->At(3));
+  
+  //Bool_t passevent = fEventCuts->IsSelected(event);
+  //Bool_t isSelectedEventSelection = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kMB);
+
+  Int_t ntracks = 0;
+  
+  if(PmdTrRp) ntracks = fRPTrackCuts->GetNumberOfInputObjects();
+  if(PmdTrPoi) ntracks = fPOITrackCuts->GetNumberOfInputObjects();
+    
+  for (Int_t i=0; i < ntracks; i++)
+    {
+      Bool_t pass = kFALSE;
+      TObject* obj = 0x0;
+      if(PmdTrPoi){
+       obj = fPOITrackCuts->GetInputObject(i);
+       if (!obj) continue;
+       pass = fPOITrackCuts->IsSelected(obj,i);
+      }
+      
+      if(PmdTrRp){
+       obj = fRPTrackCuts->GetInputObject(i);
+       if (!obj) continue;
+       pass = fRPTrackCuts->IsSelected(obj,i);
+      }
+      AliESDPmdTrack* trackpmd = dynamic_cast<AliESDPmdTrack*>(obj);
+      if (trackpmd)
+       {
+         Int_t   det   = trackpmd->GetDetector();
+         Float_t adc   = trackpmd->GetClusterADC(); 
+         Float_t ncel  = trackpmd->GetClusterCells(); 
+         if(det == 0){
+         hPmdAdcB->Fill(adc); if(pass) hPmdAdcA->Fill(adc);
+         hPmdNcelB->Fill(ncel); if(pass) hPmdNcelA->Fill(ncel); 
+         }
+       }
+    }
+}
+
+//________________________________________________________________________
+void AliAnalysisTaskQAPmdflow::Terminate(Option_t *)
+{
+  //terminate
+  
+}
+
+//________________________________________________________________________
+AliAnalysisTaskQAPmdflow::~AliAnalysisTaskQAPmdflow()
+{
+  //dtor
+  delete fRPTrackCuts;
+  delete fPOITrackCuts;
+  delete fEventCuts;
+}
+
diff --git a/PWG2/FLOW/AliFlowTasks/AliAnalysisTaskQAPmdflow.h b/PWG2/FLOW/AliFlowTasks/AliAnalysisTaskQAPmdflow.h
new file mode 100644 (file)
index 0000000..2c47bb2
--- /dev/null
@@ -0,0 +1,40 @@
+#ifndef ALIANALYSISTASKQAPMDFLOW_CXX
+#define ALIANALYSISTASKQAPMDFLOW_CXX
+
+class TObjArray;
+class TNtuple;
+class AliESDEvent;
+class AliFlowEventCuts;
+class AliFlowTrackCuts;
+
+#include "AliAnalysisTaskSE.h"
+
+class AliAnalysisTaskQAPmdflow: public AliAnalysisTaskSE
+{
+  public:
+    AliAnalysisTaskQAPmdflow();
+    AliAnalysisTaskQAPmdflow(const char* name);
+    virtual ~AliAnalysisTaskQAPmdflow();
+    
+    virtual void UserCreateOutputObjects();
+    virtual void UserExec(Option_t *option);
+    virtual void Terminate(Option_t *);
+
+
+    void SetTrackCuts(AliFlowTrackCuts* trackcutsrp) {fRPTrackCuts=trackcutsrp;}
+    void SetPOITrackCuts(AliFlowTrackCuts* trackcutspoi) {fPOITrackCuts=trackcutspoi;}
+    void SetEventCuts(AliFlowEventCuts* eventcuts) {fEventCuts=eventcuts;}
+  
+  private:
+    TObjArray* fOutput; //output histograms
+    AliFlowEventCuts* fEventCuts; //AliAnalysisCuts - applied before analysis - for comparing different event classes
+    AliFlowTrackCuts* fRPTrackCuts; //AliFlowTrackCuts go in here
+    AliFlowTrackCuts* fPOITrackCuts; //AliFlowTrackCuts go in here
+
+    AliAnalysisTaskQAPmdflow(const AliAnalysisTaskQAPmdflow&); // not implemented
+    AliAnalysisTaskQAPmdflow& operator=(const AliAnalysisTaskQAPmdflow&); // not implemented
+
+    ClassDef(AliAnalysisTaskQAPmdflow, 1); // example of analysis 
+};
+
+#endif
diff --git a/PWG2/FLOW/macros/AddTaskFlowCentralityPID.C b/PWG2/FLOW/macros/AddTaskFlowCentralityPID.C
new file mode 100644 (file)
index 0000000..a8986eb
--- /dev/null
@@ -0,0 +1,634 @@
+/////////////////////////////////////////////////////////////////////////////////////////////
+//
+// AddTask* macro for flow analysis
+// Creates a Flow Event task and adds it to the analysis manager.
+// Sets the cuts using the correction framework (CORRFW) classes.
+// Also creates Flow Analysis tasks and connects them to the output of the flow event task.
+//
+/////////////////////////////////////////////////////////////////////////////////////////////
+
+// Define the range for eta subevents (for SP method)
+Double_t minA = -0.9;
+Double_t maxA = -0.5;
+Double_t minB = 0.5;
+Double_t maxB = 0.9;
+
+// AFTERBURNER
+Bool_t useAfterBurner=kFALSE;
+Double_t v1=0.0;
+Double_t v2=0.0;
+Double_t v3=0.0;
+Double_t v4=0.0;
+Int_t numberOfTrackClones=0; //non-flow
+
+// Define a range of the detector to exclude
+Bool_t ExcludeRegion = kFALSE;
+Double_t excludeEtaMin = -0.;
+Double_t excludeEtaMax = 0.;
+Double_t excludePhiMin = 0.;
+Double_t excludePhiMax = 0.;
+
+// use physics selection class
+Bool_t  UsePhysicsSelection = kTRUE;
+
+// QA
+Bool_t runQAtask=kFALSE;
+Bool_t FillQAntuple=kFALSE;
+Bool_t DoQAcorrelations=kFALSE;
+
+// RUN SETTINGS
+// Flow analysis method can be:(set to kTRUE or kFALSE)
+Bool_t MCEP     = kFALSE;  // correlation with Monte Carlo reaction plane
+Bool_t SP       = kFALSE;  // scalar product method (similar to eventplane method)
+Bool_t GFC      = kFALSE;  // cumulants based on generating function
+Bool_t QC       = kTRUE;  // cumulants using Q vectors
+Bool_t FQD      = kFALSE;  // fit of the distribution of the Q vector (only integrated v)
+Bool_t LYZ1SUM  = kFALSE;  // Lee Yang Zeroes using sum generating function (integrated v)
+Bool_t LYZ1PROD = kFALSE;  // Lee Yang Zeroes using product generating function (integrated v)
+Bool_t LYZ2SUM  = kFALSE; // Lee Yang Zeroes using sum generating function (second pass differential v)
+Bool_t LYZ2PROD = kFALSE; // Lee Yang Zeroes using product generating function (second pass differential v)
+Bool_t LYZEP    = kFALSE; // Lee Yang Zeroes Event plane using sum generating function (gives eventplane + weight)
+Bool_t MH       = kFALSE;  // azimuthal correlators in mixed harmonics  
+Bool_t NL       = kFALSE;  // nested loops (for instance distribution of phi1-phi2 for all distinct pairs)
+
+Bool_t METHODS[] = {SP,LYZ1SUM,LYZ1PROD,LYZ2SUM,LYZ2PROD,LYZEP,GFC,QC,FQD,MCEP,MH,NL};
+
+// Boolean to use/not use weights for the Q vector
+Bool_t WEIGHTS[] = {kFALSE,kFALSE,kFALSE}; //Phi, v'(pt), v'(eta)
+
+// SETTING THE CUTS
+
+//---------Data selection----------
+//kMC, kGlobal, kESD_TPConly, kESD_SPDtracklet
+AliFlowTrackCuts::trackParameterType rptype = AliFlowTrackCuts::kGlobal;
+AliFlowTrackCuts::trackParameterType poitype = AliFlowTrackCuts::kGlobal;
+
+//---------Parameter mixing--------
+//kPure - no mixing, kTrackWithMCkine, kTrackWithMCPID, kTrackWithMCpt
+AliFlowTrackCuts::trackParameterMix rpmix = AliFlowTrackCuts::kPure;
+AliFlowTrackCuts::trackParameterMix poimix = AliFlowTrackCuts::kPure;
+
+
+const char* rptypestr = AliFlowTrackCuts::GetParamTypeName(rptype);
+const char* poitypestr = AliFlowTrackCuts::GetParamTypeName(poitype);
+
+void AddTaskFlowCentralityPID( Float_t centrMin=0.,
+                               Float_t centrMax=100.,
+                               TString fileNameBase="output",
+                               AliPID::EParticleType particleType=AliPID::kPion,
+                               AliFlowTrackCuts::PIDsource sourcePID = AliFlowTrackCuts::kTOFpid,
+                               Int_t charge=0,
+                               Int_t harmonic=2,
+                               Bool_t doQA=kFALSE )
+{
+  //===========================================================================
+  // EVENTS CUTS:
+  AliFlowEventCuts* cutsEvent = new AliFlowEventCuts("event cuts");
+  cutsEvent->SetCentralityPercentileRange(centrMin,centrMax);
+  cutsEvent->SetCentralityPercentileMethod(AliFlowEventCuts::kV0);
+  cutsEvent->SetRefMultMethod(AliFlowEventCuts::kV0);
+  //cutsEvent->SetCentralityPercentileMethod(AliFlowEventCuts::kSPD1tracklets);
+  cutsEvent->SetNContributorsRange(2);
+  cutsEvent->SetPrimaryVertexZrange(-7.,7.);
+  cutsEvent->SetCutSPDvertexerAnomaly(); //"Francesco's cut"
+  cutsEvent->SetCutZDCtiming();
+  cutsEvent->SetCutTPCmultiplicityOutliers();
+  cutsEvent->SetQA(doQA);
+  
+  // RP TRACK CUTS:
+  AliFlowTrackCuts* cutsRP = new AliFlowTrackCuts("TPConlyRP");
+  cutsRP->SetParamType(rptype);
+  cutsRP->SetParamMix(rpmix);
+  cutsRP->SetPtRange(0.2,5.);
+  cutsRP->SetEtaRange(-0.8,0.8);
+  cutsRP->SetMinNClustersTPC(70);
+  cutsRP->SetMinChi2PerClusterTPC(0.1);
+  cutsRP->SetMaxChi2PerClusterTPC(4.0);
+  cutsRP->SetMaxDCAToVertexXY(3.0);
+  cutsRP->SetMaxDCAToVertexZ(3.0);
+  cutsRP->SetAcceptKinkDaughters(kFALSE);
+  cutsRP->SetMinimalTPCdedx(10.);
+  cutsRP->SetQA(doQA);
+
+  // POI TRACK CUTS:
+  AliFlowTrackCuts* cutsPOI = new AliFlowTrackCuts("GlobalPOI");
+  cutsPOI->SetParamType(poitype);
+  cutsPOI->SetParamMix(poimix);
+  cutsPOI->SetPtRange(0.0,10.);
+  cutsPOI->SetEtaRange(-0.8,0.8);
+  //cutsPOI->SetRequireCharge(kTRUE);
+  //cutsPOI->SetPID(PdgRP);
+  cutsPOI->SetMinNClustersTPC(70);
+  cutsPOI->SetMinChi2PerClusterTPC(0.1);
+  cutsPOI->SetMaxChi2PerClusterTPC(4.0);
+  cutsPOI->SetRequireITSRefit(kTRUE);
+  cutsPOI->SetRequireTPCRefit(kTRUE);
+  cutsPOI->SetMinNClustersITS(2);
+  //cutsPOI->SetMaxChi2PerClusterITS(1.e+09);
+  cutsPOI->SetMaxDCAToVertexXY(0.3);
+  cutsPOI->SetMaxDCAToVertexZ(0.3);
+  //cutsPOI->SetDCAToVertex2D(kTRUE);
+  //cutsPOI->SetMaxNsigmaToVertex(1.e+10);
+  //cutsPOI->SetRequireSigmaToVertex(kFALSE);
+  cutsPOI->SetAcceptKinkDaughters(kFALSE);
+  cutsPOI->SetPID(particleType, sourcePID);
+  if (charge!=0) cutsPOI->SetCharge(charge);
+  cutsPOI->SetAllowTOFmismatch(kFALSE);
+  //iexample: francesco's tunig TPC Bethe Bloch for data:
+  //cutsPOI->GetESDpid().GetTPCResponse().SetBetheBlochParameters(4.36414e-02,1.75977e+01,1.14385e-08,2.27907e+00,3.36699e+00);
+  //cutsPOI->GetESDpid().GetTPCResponse().SetMip(49);
+  cutsPOI->SetMinimalTPCdedx(10.);
+  cutsPOI->SetQA(doQA);
+
+  TString outputSlotName("");
+  outputSlotName+=Form("V%i ",harmonic);
+  outputSlotName+=cutsRP->GetName();
+  outputSlotName+=" ";
+  outputSlotName+=cutsPOI->GetName();
+  outputSlotName+=Form(" %.0f-",centrMin);
+  outputSlotName+=Form("%.0f ",centrMax);
+  outputSlotName+=AliFlowTrackCuts::PIDsourceName(sourcePID);
+  outputSlotName+=" ";
+  outputSlotName+=AliPID::ParticleName(particleType);
+  if (charge<0) outputSlotName+="-";
+  if (charge>0) outputSlotName+="+";
+
+  TString fileName(fileNameBase);
+  fileName.Append(".root");
+
+  Bool_t useWeights  = WEIGHTS[0] || WEIGHTS[1] || WEIGHTS[2];
+  if (useWeights) cout<<"Weights are used"<<endl;
+  else cout<<"Weights are not used"<<endl;
+  
+  // Get the pointer to the existing analysis manager via the static access method.
+  //==============================================================================
+  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
+  if (!mgr) {
+    Error("AddTaskFlowEvent", "No analysis manager to connect to.");
+    return NULL;
+  }
+  
+  // Check the analysis type using the event handlers connected to the analysis
+  // manager. The availability of MC handler can also be checked here.
+  //==============================================================================
+  if (!mgr->GetInputEventHandler()) {
+    ::Error("AddTaskFlowEvent", "This task requires an input event handler");
+    return NULL;
+  }  
+
+  // Open external input files
+  //===========================================================================
+  //weights: 
+  TFile *weightsFile = NULL;
+  TList *weightsList = NULL;
+
+  if(useWeights) {
+    //open the file with the weights:
+    weightsFile = TFile::Open("weights.root","READ");
+    if(weightsFile) {
+      //access the list which holds the histos with weigths:
+      weightsList = (TList*)weightsFile->Get("weights");
+    }
+    else {
+      cout<<" WARNING: the file <weights.root> with weights from the previous run was not available."<<endl;
+      break;
+    } 
+  }
+  
+  //LYZ2
+  if (LYZ2SUM || LYZ2PROD) {
+    //read the outputfile of the first run
+    TString outputFileName = "AnalysisResults1.root";
+    TString pwd(gSystem->pwd());
+    pwd+="/";
+    pwd+=outputFileName.Data();
+    TFile *outputFile = NULL;
+    if(gSystem->AccessPathName(pwd.Data(),kFileExists)) {
+      cout<<"WARNING: You do not have an output file:"<<endl;
+      cout<<"         "<<pwd.Data()<<endl;
+      exit(0);
+    } else { outputFile = TFile::Open(pwd.Data(),"READ");}
+    
+    if (LYZ2SUM){  
+      // read the output directory from LYZ1SUM 
+      TString inputFileNameLYZ2SUM = "outputLYZ1SUManalysis" ;
+      inputFileNameLYZ2SUM += rptypestr;
+      cout<<"The input directory is "<<inputFileNameLYZ2SUM.Data()<<endl;
+      TFile* fInputFileLYZ2SUM = (TFile*)outputFile->FindObjectAny(inputFileNameLYZ2SUM.Data());
+      if(!fInputFileLYZ2SUM || fInputFileLYZ2SUM->IsZombie()) { 
+       cerr << " ERROR: To run LYZ2SUM you need the output file from LYZ1SUM. This file is not there! Please run LYZ1SUM first." << endl ; 
+       break;
+      }
+      else {
+       TList* fInputListLYZ2SUM = (TList*)fInputFileLYZ2SUM->Get("LYZ1SUM");
+       if (!fInputListLYZ2SUM) {cout<<"list is NULL pointer!"<<endl;}
+      }
+      cout<<"LYZ2SUM input file/list read..."<<endl;
+    }
+
+    if (LYZ2PROD){  
+      // read the output directory from LYZ1PROD 
+      TString inputFileNameLYZ2PROD = "outputLYZ1PRODanalysis" ;
+      inputFileNameLYZ2PROD += rptypestr;
+      cout<<"The input directory is "<<inputFileNameLYZ2PROD.Data()<<endl;
+      TFile* fInputFileLYZ2PROD = (TFile*)outputFile->FindObjectAny(inputFileNameLYZ2PROD.Data());
+      if(!fInputFileLYZ2PROD || fInputFileLYZ2PROD->IsZombie()) { 
+       cerr << " ERROR: To run LYZ2PROD you need the output file from LYZ1PROD. This file is not there! Please run LYZ1PROD first." << endl ; 
+       break;
+      }
+      else {
+       TList* fInputListLYZ2PROD = (TList*)fInputFileLYZ2PROD->Get("LYZ1PROD");
+       if (!fInputListLYZ2PROD) {cout<<"list is NULL pointer!"<<endl;}
+      }
+      cout<<"LYZ2PROD input file/list read..."<<endl;
+    }
+  }
+
+  if (LYZEP) {
+    //read the outputfile of the second run
+    TString outputFileName = "AnalysisResults2.root";
+    TString pwd(gSystem->pwd());
+    pwd+="/";
+    pwd+=outputFileName.Data();
+    TFile *outputFile = NULL;
+    if(gSystem->AccessPathName(pwd.Data(),kFileExists)) {
+      cout<<"WARNING: You do not have an output file:"<<endl;
+      cout<<"         "<<pwd.Data()<<endl;
+      exit(0);
+    } else {
+      outputFile = TFile::Open(pwd.Data(),"READ");
+    }
+    
+    // read the output file from LYZ2SUM
+    TString inputFileNameLYZEP = "outputLYZ2SUManalysis" ;
+    inputFileNameLYZEP += rptypestr;
+    cout<<"The input file is "<<inputFileNameLYZEP.Data()<<endl;
+    TFile* fInputFileLYZEP = (TFile*)outputFile->FindObjectAny(inputFileNameLYZEP.Data());
+    if(!fInputFileLYZEP || fInputFileLYZEP->IsZombie()) { 
+      cerr << " ERROR: To run LYZEP you need the output file from LYZ2SUM. This file is not there! Please run LYZ2SUM first." << endl ; 
+      break;
+    }
+    else {
+      TList* fInputListLYZEP = (TList*)fInputFileLYZEP->Get("LYZ2SUM");
+      if (!fInputListLYZEP) {cout<<"list is NULL pointer!"<<endl;}
+    }
+    cout<<"LYZEP input file/list read..."<<endl;
+  }
+  
+  
+  // Create the FMD task and add it to the manager
+  //===========================================================================
+  if (rptypestr == "FMD") {
+    AliFMDAnalysisTaskSE *taskfmd = NULL;
+    if (rptypestr == "FMD") {
+      taskfmd = new AliFMDAnalysisTaskSE("TaskFMD");
+      mgr->AddTask(taskfmd);
+      
+      AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
+      pars->Init();
+      pars->SetProcessPrimary(kTRUE); //for MC only
+      pars->SetProcessHits(kFALSE);
+      
+      //pars->SetRealData(kTRUE); //for real data
+      //pars->SetProcessPrimary(kFALSE); //for real data
+    }
+  }
+  
+  // Create the flow event task, add it to the manager.
+  //===========================================================================
+  AliAnalysisTaskFlowEvent *taskFE = NULL;
+
+  if(useAfterBurner)
+    { 
+      taskFE = new AliAnalysisTaskFlowEvent(Form("TaskFlowEvent %s",outputSlotName.Data()),"",doQA,1);
+      taskFE->SetFlow(v1,v2,v3,v4); 
+      taskFE->SetNonFlowNumberOfTrackClones(numberOfTrackClones);
+      taskFE->SetAfterburnerOn();
+    }
+  else {taskFE = new AliAnalysisTaskFlowEvent(Form("TaskFlowEvent %s",outputSlotName.Data()),"",doQA); }
+  if (ExcludeRegion) {
+    taskFE->DefineDeadZone(excludeEtaMin, excludeEtaMax, excludePhiMin, excludePhiMax); 
+  }
+  taskFE->SetSubeventEtaRange(minA, maxA, minB, maxB);
+  if (UsePhysicsSelection) {
+    taskFE->SelectCollisionCandidates(AliVEvent::kMB);
+    cout<<"Using Physics Selection"<<endl;
+  }
+  mgr->AddTask(taskFE);
+  
+  // Pass cuts for RPs and POIs to the task:
+  taskFE->SetCutsEvent(cutsEvent);
+  taskFE->SetCutsRP(cutsRP);
+  taskFE->SetCutsPOI(cutsPOI);
+
+  // Create the analysis tasks, add them to the manager.
+  //===========================================================================
+  if (SP){
+    AliAnalysisTaskScalarProduct *taskSP = new AliAnalysisTaskScalarProduct(Form("TaskScalarProduct %s",outputSlotName.Data()),WEIGHTS[0]);
+    taskSP->SetRelDiffMsub(1.0);
+    taskSP->SetApplyCorrectionForNUA(kTRUE);
+    mgr->AddTask(taskSP);
+  }
+  if (LYZ1SUM){
+    AliAnalysisTaskLeeYangZeros *taskLYZ1SUM = new AliAnalysisTaskLeeYangZeros(Form("TaskLeeYangZerosSUM %s",outputSlotName.Data()),kTRUE);
+    taskLYZ1SUM->SetFirstRunLYZ(kTRUE);
+    taskLYZ1SUM->SetUseSumLYZ(kTRUE);
+    mgr->AddTask(taskLYZ1SUM);
+  }
+  if (LYZ1PROD){
+    AliAnalysisTaskLeeYangZeros *taskLYZ1PROD = new AliAnalysisTaskLeeYangZeros(Form("TaskLeeYangZerosPROD %s",outputSlotName.Data()),kTRUE);
+    taskLYZ1PROD->SetFirstRunLYZ(kTRUE);
+    taskLYZ1PROD->SetUseSumLYZ(kFALSE);
+    mgr->AddTask(taskLYZ1PROD);
+  }
+  if (LYZ2SUM){
+    AliAnalysisTaskLeeYangZeros *taskLYZ2SUM = new AliAnalysisTaskLeeYangZeros(Form("TaskLeeYangZerosSUM %s",outputSlotName.Data()),kFALSE);
+    taskLYZ2SUM->SetFirstRunLYZ(kFALSE);
+    taskLYZ2SUM->SetUseSumLYZ(kTRUE);
+    mgr->AddTask(taskLYZ2SUM);
+  }
+  if (LYZ2PROD){
+    AliAnalysisTaskLeeYangZeros *taskLYZ2PROD = new AliAnalysisTaskLeeYangZeros(Form("TaskLeeYangZerosPROD %s",outputSlotName.Data()),kFALSE);
+    taskLYZ2PROD->SetFirstRunLYZ(kFALSE);
+    taskLYZ2PROD->SetUseSumLYZ(kFALSE);
+    mgr->AddTask(taskLYZ2PROD);
+  }
+  if (LYZEP){
+    AliAnalysisTaskLYZEventPlane *taskLYZEP = new AliAnalysisTaskLYZEventPlane(Form("TaskLYZEventPlane %s",outputSlotName.Data()));
+    mgr->AddTask(taskLYZEP);
+  }
+  if (GFC){
+    AliAnalysisTaskCumulants *taskGFC = new AliAnalysisTaskCumulants(Form("TaskCumulants %s",outputSlotName.Data()),useWeights);
+    taskGFC->SetUsePhiWeights(WEIGHTS[0]); 
+    taskGFC->SetUsePtWeights(WEIGHTS[1]);
+    taskGFC->SetUseEtaWeights(WEIGHTS[2]); 
+    mgr->AddTask(taskGFC);
+  }
+  if (QC){
+    AliAnalysisTaskQCumulants *taskQC = new AliAnalysisTaskQCumulants(Form("TaskQCumulants %s",outputSlotName.Data()),useWeights);
+    taskQC->SetUsePhiWeights(WEIGHTS[0]); 
+    taskQC->SetUsePtWeights(WEIGHTS[1]);
+    taskQC->SetUseEtaWeights(WEIGHTS[2]); 
+    taskQC->SetCalculateCumulantsVsM(kFALSE);
+    taskQC->SetnBinsMult(10000);
+    taskQC->SetMinMult(0.);
+    taskQC->SetMaxMult(10000.);
+    taskQC->SetHarmonic(harmonic);
+    taskQC->SetApplyCorrectionForNUA(kFALSE);
+    taskQC->SetFillMultipleControlHistograms(kFALSE);     
+    mgr->AddTask(taskQC);
+  }
+  if (FQD){
+    AliAnalysisTaskFittingQDistribution *taskFQD = new AliAnalysisTaskFittingQDistribution(Form("TaskFittingQDistribution %s",outputSlotName.Data()),kFALSE);
+    taskFQD->SetUsePhiWeights(WEIGHTS[0]); 
+    taskFQD->SetqMin(0.);
+    taskFQD->SetqMax(1000.);
+    taskFQD->SetqNbins(10000);
+    mgr->AddTask(taskFQD);
+  }
+  if (MCEP){
+    AliAnalysisTaskMCEventPlane *taskMCEP = new AliAnalysisTaskMCEventPlane(Form("TaskMCEventPlane %s",outputSlotName.Data()));
+    mgr->AddTask(taskMCEP);
+  }
+  if (MH){
+    AliAnalysisTaskMixedHarmonics *taskMH = new AliAnalysisTaskMixedHarmonics(Form("TaskMixedHarmonics %s",outputSlotName.Data()),useWeights);
+    taskMH->SetHarmonic(1); // n in cos[n(phi1+phi2-2phi3)] and cos[n(psi1+psi2-2phi3)]
+    taskMH->SetNoOfMultipicityBins(10000);
+    taskMH->SetMultipicityBinWidth(1.);
+    taskMH->SetMinMultiplicity(1.);
+    taskMH->SetCorrectForDetectorEffects(kTRUE);
+    taskMH->SetEvaluateDifferential3pCorrelator(kFALSE); // evaluate <<cos[n(psi1+psi2-2phi3)]>> (Remark: two nested loops)    
+    taskMH->SetOppositeChargesPOI(kFALSE); // POIs psi1 and psi2 in cos[n(psi1+psi2-2phi3)] will have opposite charges  
+    mgr->AddTask(taskMH);
+  }  
+  if (NL){
+    AliAnalysisTaskNestedLoops *taskNL = new AliAnalysisTaskNestedLoops(Form("TaskNestedLoops %s",outputSlotName.Data()),useWeights);
+    taskNL->SetHarmonic(1); // n in cos[n(phi1+phi2-2phi3)] and cos[n(psi1+psi2-2phi3)]
+    taskNL->SetEvaluateNestedLoopsForRAD(kTRUE); // RAD = Relative Angle Distribution
+    taskNL->SetEvaluateNestedLoopsForMH(kTRUE); // evalaute <<cos[n(phi1+phi2-2phi3)]>> (Remark: three nested loops)   
+    taskNL->SetEvaluateDifferential3pCorrelator(kFALSE); // evaluate <<cos[n(psi1+psi2-2phi3)]>>  (Remark: three nested loops)   
+    taskNL->SetOppositeChargesPOI(kFALSE); // POIs psi1 and psi2 in cos[n(psi1+psi2-2phi3)] will have opposite charges  
+    mgr->AddTask(taskNL);
+  }
+
+  // Create the output container for the data produced by the task
+  // Connect to the input and output containers
+  //===========================================================================
+  AliAnalysisDataContainer *cinput1 = mgr->GetCommonInputContainer();
+  
+  if (rptypestr == "FMD") {
+    AliAnalysisDataContainer *coutputFMD = 
+      mgr->CreateContainer(Form("BackgroundCorrected %s",outputSlotName.Data()), TList::Class(), AliAnalysisManager::kExchangeContainer);
+    //input and output taskFMD     
+    mgr->ConnectInput(taskfmd, 0, cinput1);
+    mgr->ConnectOutput(taskfmd, 1, coutputFMD);
+    //input into taskFE
+    mgr->ConnectInput(taskFE,1,coutputFMD);
+  }
+  
+  AliAnalysisDataContainer *coutputFE = 
+  mgr->CreateContainer(Form("FlowEventSimple %s",outputSlotName.Data()),AliFlowEventSimple::Class(),AliAnalysisManager::kExchangeContainer);
+  mgr->ConnectInput(taskFE,0,cinput1); 
+  mgr->ConnectOutput(taskFE,1,coutputFE);
+  if (taskFE->GetQAOn())
+  {
+    AliAnalysisDataContainer* coutputFEQA = 
+    mgr->CreateContainer(Form("QA %s",outputSlotName.Data()), TList::Class(),AliAnalysisManager::kOutputContainer,Form("%s:QA",fileName));
+    mgr->ConnectOutput(taskFE,2,coutputFEQA);
+  }
+
+  // Create the output containers for the data produced by the analysis tasks
+  // Connect to the input and output containers
+  //===========================================================================
+  if (useWeights) {    
+    AliAnalysisDataContainer *cinputWeights = mgr->CreateContainer(Form("Weights %s",outputSlotName.Data()),
+                                                                  TList::Class(),AliAnalysisManager::kInputContainer); 
+  }
+
+  if(SP) {
+    TString outputSP = fileName;
+    outputSP += ":outputSPanalysis";
+    outputSP+= rptypestr;
+    AliAnalysisDataContainer *coutputSP = mgr->CreateContainer(Form("SP %s",outputSlotName.Data()), 
+                                                              TList::Class(),AliAnalysisManager::kOutputContainer,outputSP); 
+    mgr->ConnectInput(taskSP,0,coutputFE); 
+    mgr->ConnectOutput(taskSP,1,coutputSP); 
+    if (WEIGHTS[0]) {
+      mgr->ConnectInput(taskSP,1,cinputWeights);
+      cinputWeights->SetData(weightsList);
+    }
+  }
+  if(LYZ1SUM) {
+    TString outputLYZ1SUM = fileName;
+    outputLYZ1SUM += ":outputLYZ1SUManalysis";
+    outputLYZ1SUM+= rptypestr;
+    AliAnalysisDataContainer *coutputLYZ1SUM = mgr->CreateContainer(Form("LYZ1SUM %s",outputSlotName.Data()), 
+                                                                   TList::Class(),AliAnalysisManager::kOutputContainer,outputLYZ1SUM); 
+    mgr->ConnectInput(taskLYZ1SUM,0,coutputFE);
+    mgr->ConnectOutput(taskLYZ1SUM,1,coutputLYZ1SUM);
+  }
+  if(LYZ1PROD) {
+    TString outputLYZ1PROD = fileName;
+    outputLYZ1PROD += ":outputLYZ1PRODanalysis";
+    outputLYZ1PROD+= rptypestr;
+    AliAnalysisDataContainer *coutputLYZ1PROD = mgr->CreateContainer(Form("LYZ1PROD %s",outputSlotName.Data()), 
+                                                                    TList::Class(),AliAnalysisManager::kOutputContainer,outputLYZ1PROD); 
+    mgr->ConnectInput(taskLYZ1PROD,0,coutputFE); 
+    mgr->ConnectOutput(taskLYZ1PROD,1,coutputLYZ1PROD);
+  }
+  if(LYZ2SUM) {
+    AliAnalysisDataContainer *cinputLYZ2SUM = mgr->CreateContainer(Form("LYZ2SUMin %s",outputSlotName.Data()),
+                                                                  TList::Class(),AliAnalysisManager::kInputContainer);
+    TString outputLYZ2SUM = fileName;
+    outputLYZ2SUM += ":outputLYZ2SUManalysis";
+    outputLYZ2SUM+= rptypestr;
+    
+    AliAnalysisDataContainer *coutputLYZ2SUM = mgr->CreateContainer(Form("LYZ2SUM %s",outputSlotName.Data()), 
+                                                                   TList::Class(),AliAnalysisManager::kOutputContainer,outputLYZ2SUM); 
+    mgr->ConnectInput(taskLYZ2SUM,0,coutputFE); 
+    mgr->ConnectInput(taskLYZ2SUM,1,cinputLYZ2SUM);
+    mgr->ConnectOutput(taskLYZ2SUM,1,coutputLYZ2SUM); 
+    cinputLYZ2SUM->SetData(fInputListLYZ2SUM);
+  }
+  if(LYZ2PROD) {
+    AliAnalysisDataContainer *cinputLYZ2PROD = mgr->CreateContainer(Form("LYZ2PRODin %s",outputSlotName.Data()),
+                                                                   TList::Class(),AliAnalysisManager::kInputContainer);
+    TString outputLYZ2PROD = fileName;
+    outputLYZ2PROD += ":outputLYZ2PRODanalysis";
+    outputLYZ2PROD+= rptypestr;
+    
+    AliAnalysisDataContainer *coutputLYZ2PROD = mgr->CreateContainer(Form("LYZ2PROD %s",outputSlotName.Data()), 
+                                                                    TList::Class(),AliAnalysisManager::kOutputContainer,outputLYZ2PROD); 
+    mgr->ConnectInput(taskLYZ2PROD,0,coutputFE); 
+    mgr->ConnectInput(taskLYZ2PROD,1,cinputLYZ2PROD);
+    mgr->ConnectOutput(taskLYZ2PROD,1,coutputLYZ2PROD); 
+    cinputLYZ2PROD->SetData(fInputListLYZ2PROD);
+  }
+  if(LYZEP) {
+    AliAnalysisDataContainer *cinputLYZEP = mgr->CreateContainer(Form("LYZEPin %s",outputSlotName.Data()),
+                                                                TList::Class(),AliAnalysisManager::kInputContainer);
+    TString outputLYZEP = fileName;
+    outputLYZEP += ":outputLYZEPanalysis";
+    outputLYZEP+= rptypestr;
+    
+    AliAnalysisDataContainer *coutputLYZEP = mgr->CreateContainer(Form("LYZEP %s",outputSlotName.Data()), 
+                                                                 TList::Class(),AliAnalysisManager::kOutputContainer,outputLYZEP); 
+    mgr->ConnectInput(taskLYZEP,0,coutputFE); 
+    mgr->ConnectInput(taskLYZEP,1,cinputLYZEP);
+    mgr->ConnectOutput(taskLYZEP,1,coutputLYZEP); 
+    cinputLYZEP->SetData(fInputListLYZEP);
+  }
+  if(GFC) {
+    TString outputGFC = fileName;
+    outputGFC += ":outputGFCanalysis";
+    outputGFC+= rptypestr;
+    
+    AliAnalysisDataContainer *coutputGFC = mgr->CreateContainer(Form("GFC %s",outputSlotName.Data()), 
+                                                               TList::Class(),AliAnalysisManager::kOutputContainer,outputGFC); 
+    mgr->ConnectInput(taskGFC,0,coutputFE); 
+    mgr->ConnectOutput(taskGFC,1,coutputGFC);
+    if (useWeights) {
+      mgr->ConnectInput(taskGFC,1,cinputWeights);
+      cinputWeights->SetData(weightsList);
+    } 
+  }
+  if(QC) {
+    TString outputQC = fileName;
+    outputQC += ":outputQCanalysis";
+    outputQC+= rptypestr;
+
+    AliAnalysisDataContainer *coutputQC = mgr->CreateContainer(Form("QC %s",outputSlotName.Data()), 
+                                                              TList::Class(),AliAnalysisManager::kOutputContainer,outputQC); 
+    mgr->ConnectInput(taskQC,0,coutputFE); 
+    mgr->ConnectOutput(taskQC,1,coutputQC);
+    if (useWeights) {
+      mgr->ConnectInput(taskQC,1,cinputWeights);
+      cinputWeights->SetData(weightsList);
+    }
+  }
+  if(FQD) {
+    TString outputFQD = fileName;
+    outputFQD += ":outputFQDanalysis";
+    outputFQD+= rptypestr;
+    
+    AliAnalysisDataContainer *coutputFQD = mgr->CreateContainer(Form("FQD %s",outputSlotName.Data()), 
+                                                               TList::Class(),AliAnalysisManager::kOutputContainer,outputFQD); 
+    mgr->ConnectInput(taskFQD,0,coutputFE); 
+    mgr->ConnectOutput(taskFQD,1,coutputFQD);
+    if(useWeights) {
+      mgr->ConnectInput(taskFQD,1,cinputWeights);
+      cinputWeights->SetData(weightsList);
+    } 
+  }
+  if(MCEP) {
+    TString outputMCEP = fileName;
+    outputMCEP += ":outputMCEPanalysis";
+    outputMCEP+= rptypestr;
+    
+    AliAnalysisDataContainer *coutputMCEP = mgr->CreateContainer(Form("MCEP %s",outputSlotName.Data()), 
+                                                                TList::Class(),AliAnalysisManager::kOutputContainer,outputMCEP); 
+    mgr->ConnectInput(taskMCEP,0,coutputFE);
+    mgr->ConnectOutput(taskMCEP,1,coutputMCEP); 
+  }
+  if(MH) {
+    TString outputMH = fileName;
+    outputMH += ":outputMHanalysis";
+    outputMH += rptypestr;
+        
+    AliAnalysisDataContainer *coutputMH = mgr->CreateContainer(Form("MH %s",outputSlotName.Data()), 
+                                                              TList::Class(),AliAnalysisManager::kOutputContainer,outputMH); 
+    mgr->ConnectInput(taskMH,0,coutputFE); 
+    mgr->ConnectOutput(taskMH,1,coutputMH); 
+    //if (useWeights) {
+    //  mgr->ConnectInput(taskMH,1,cinputWeights);
+    //  cinputWeights->SetData(weightsList);
+    //} 
+  }
+  if(NL) {
+    TString outputNL = fileName;
+    outputNL += ":outputNLanalysis";
+    outputNL += rptypestr;
+
+    AliAnalysisDataContainer *coutputNL = mgr->CreateContainer(Form("NL %s",outputSlotName.Data()), 
+                                                              TList::Class(),AliAnalysisManager::kOutputContainer,outputNL); 
+    mgr->ConnectInput(taskNL,0,coutputFE);
+    mgr->ConnectOutput(taskNL,1,coutputNL);
+    //if (useWeights) {
+    //  mgr->ConnectInput(taskNL,1,cinputWeights);
+    //  cinputWeights->SetData(weightsList);
+    //} 
+  }
+
+  ///////////////////////////////////////////////////////////////////////////////////////////
+  if (runQAtask)
+  {
+    AliAnalysisTaskQAflow* taskQAflow = new AliAnalysisTaskQAflow(Form("TaskQAflow %s",outputSlotName.Data()));
+    taskQAflow->SetEventCuts(cutsEvent);
+    taskQAflow->SetTrackCuts(cutsRP);
+    taskQAflow->SetFillNTuple(FillQAntuple);
+    taskQAflow->SetDoCorrelations(DoQAcorrelations);
+    mgr->AddTask(taskQAflow);
+    
+    Printf("outputSlotName %s",outputSlotName.Data());
+    TString taskQAoutputFileName(fileNameBase);
+    taskQAoutputFileName.Append("_QA.root");
+    AliAnalysisDataContainer* coutputQAtask = mgr->CreateContainer(Form("flowQA %s",outputSlotName.Data()),
+                                              TObjArray::Class(),
+                                              AliAnalysisManager::kOutputContainer,
+                                              taskQAoutputFileName);
+    AliAnalysisDataContainer* coutputQAtaskTree = mgr->CreateContainer(Form("flowQAntuple %s",outputSlotName.Data()),
+                                              TNtuple::Class(),
+                                              AliAnalysisManager::kOutputContainer,
+                                              taskQAoutputFileName);
+    mgr->ConnectInput(taskQAflow,0,mgr->GetCommonInputContainer());
+    mgr->ConnectInput(taskQAflow,1,coutputFE);
+    mgr->ConnectOutput(taskQAflow,1,coutputQAtask);
+    if (FillQAntuple) mgr->ConnectOutput(taskQAflow,2,coutputQAtaskTree);
+  }
+}
+
+
+
+
+
diff --git a/PWG2/FLOW/macros/runFlowTaskCentralityPIDTrain.C b/PWG2/FLOW/macros/runFlowTaskCentralityPIDTrain.C
new file mode 100644 (file)
index 0000000..78ec247
--- /dev/null
@@ -0,0 +1,370 @@
+enum anaModes {mLocal,mPROOF,mGrid};
+//mLocal: Analyze locally files in your computer using aliroot
+//mPROOF: Analyze CAF files with PROOF
+//mGrid: Analyze files on Grid via AliEn plug-in and using precompiled FLOW libraries
+
+// CENTRALITY DEFINITION
+//Int_t binfirst = 4;  //where do we start numbering bins
+//Int_t binlast = 6;  //where do we stop numbering bins
+//const Int_t numberOfCentralityBins = 9;
+Int_t binfirst = 0;  //where do we start numbering bins
+Int_t binlast = 8;  //where do we stop numbering bins
+const Int_t numberOfCentralityBins = 9;
+Float_t centralityArray[numberOfCentralityBins+1] = {0.,5.,10.,20.,30.,40.,50.,60.,70.,80.}; // in centrality percentile
+//Int_t centralityArray[numberOfCentralityBins+1] = {41,80,146,245,384,576,835,1203,1471,10000}; // in terms of TPC only reference multiplicity
+
+TString commonOutputFileName = "outputCentrality"; // e.g.: result for centrality bin 0 will be in the file "outputCentrality0.root", etc
+
+//void runFlowTaskCentralityPIDTrain(Int_t mode=mLocal, Int_t nEvents = 10,
+//Bool_t DATA = kFALSE, const Char_t* dataDir="/Users/snelling/alice_data/Therminator_midcentral", Int_t offset = 0)
+
+void runFlowTaskCentralityPIDTrain( Int_t mode = mGrid,
+                                    Bool_t useFlowParFiles = kTRUE,
+                                    Bool_t DATA = kTRUE,
+                                    //const Char_t* dataDir="/data/alice2/ab/grid/21-16/TEST/data/",
+                                    const Char_t* dataDir="/data/alice2/mikolaj/flowPaper/data/137161/",
+                                    Int_t nEvents = 1e4,
+                                    Int_t offset=0 )
+{
+  // Time:
+  TStopwatch timer;
+  timer.Start();
+  // Cross-check user settings before starting:
+  //  CrossCheckUserSettings(DATA);
+
+  // Load needed libraries:
+  LoadLibraries(mode,useFlowParFiles);
+
+  // Create analysis manager:
+  AliAnalysisManager *mgr = new AliAnalysisManager("FlowAnalysisManager");
+
+  // Chains:
+  if(mode == mLocal)
+  {
+    gROOT->LoadMacro("$ALICE_ROOT/PWG0/CreateESDChain.C");
+    TChain* chain = CreateESDChain(dataDir, nEvents, offset);
+    //TChain* chain = CreateAODChain(dataDir, nEvents, offset);
+  }
+
+  // Connect plug-in to the analysis manager:
+  if(mode == mGrid)
+  {
+    gROOT->LoadMacro("CreateAlienHandler.C");
+    AliAnalysisGrid *alienHandler = CreateAlienHandler(useFlowParFiles);
+    if(!alienHandler) return;
+    mgr->SetGridHandler(alienHandler);
+  }
+
+  // Event handlers:
+  AliVEventHandler* esdH = new AliESDInputHandler;
+  mgr->SetInputEventHandler(esdH);
+  if (!DATA)
+  {
+    AliMCEventHandler *mc = new AliMCEventHandler();
+    mgr->SetMCtruthEventHandler(mc);
+  }
+
+  // Task to check the offline trigger:
+  gROOT->LoadMacro("$ALICE_ROOT/ANALYSIS/macros/AddTaskPhysicsSelection.C");
+  AddTaskPhysicsSelection(!DATA);
+
+  //Add the centrality determination task
+  gROOT->LoadMacro("$ALICE_ROOT/ANALYSIS/macros/AddTaskCentrality.C");
+  AddTaskCentrality();
+
+  //Add the TOF tender
+  gROOT->LoadMacro("$ALICE_ROOT/PWG2/FLOW/macros/AddTaskTenderTOF.C");
+  AddTaskTenderTOF();
+
+  // Setup analysis per centrality bin:
+  gROOT->LoadMacro("AddTaskFlowCentralityPID.C");
+  for (Int_t i=binfirst; i<binlast+1; i++)
+  {
+    Float_t lowCentralityBinEdge = centralityArray[i];
+    Float_t highCentralityBinEdge = centralityArray[i+1];
+    AddTaskFlowCentralityPID( lowCentralityBinEdge,
+                              highCentralityBinEdge,
+                              commonOutputFileName,
+                              AliPID::kPion,
+                              AliFlowTrackCuts::kTOFbetaSimple,
+                              -1,2,kTRUE );
+    AddTaskFlowCentralityPID( lowCentralityBinEdge,
+                              highCentralityBinEdge,
+                              commonOutputFileName,
+                              AliPID::kPion,
+                              AliFlowTrackCuts::kTOFbetaSimple,
+                              1,2,kTRUE );
+    AddTaskFlowCentralityPID( lowCentralityBinEdge,
+                              highCentralityBinEdge,
+                              commonOutputFileName,
+                              AliPID::kKaon,
+                              AliFlowTrackCuts::kTOFbetaSimple,
+                              -1,2,kTRUE );
+    AddTaskFlowCentralityPID( lowCentralityBinEdge,
+                              highCentralityBinEdge,
+                              commonOutputFileName,
+                              AliPID::kKaon,
+                              AliFlowTrackCuts::kTOFbetaSimple,
+                              1,2,kTRUE );
+    AddTaskFlowCentralityPID( lowCentralityBinEdge,
+                              highCentralityBinEdge,
+                              commonOutputFileName,
+                              AliPID::kProton,
+                              AliFlowTrackCuts::kTOFbetaSimple,
+                              -1,2,kTRUE );
+    AddTaskFlowCentralityPID( lowCentralityBinEdge,
+                              highCentralityBinEdge,
+                              commonOutputFileName,
+                              AliPID::kProton,
+                              AliFlowTrackCuts::kTOFbetaSimple,
+                              1,2,kTRUE );
+    AddTaskFlowCentralityPID( lowCentralityBinEdge,
+                              highCentralityBinEdge,
+                              commonOutputFileName,
+                              AliPID::kPion,
+                              AliFlowTrackCuts::kTOFbetaSimple,
+                              -1,3 );
+    AddTaskFlowCentralityPID( lowCentralityBinEdge,
+                              highCentralityBinEdge,
+                              commonOutputFileName,
+                              AliPID::kPion,
+                              AliFlowTrackCuts::kTOFbetaSimple,
+                              1,3 );
+    AddTaskFlowCentralityPID( lowCentralityBinEdge,
+                              highCentralityBinEdge,
+                              commonOutputFileName,
+                              AliPID::kKaon,
+                              AliFlowTrackCuts::kTOFbetaSimple,
+                              -1,3 );
+    AddTaskFlowCentralityPID( lowCentralityBinEdge,
+                              highCentralityBinEdge,
+                              commonOutputFileName,
+                              AliPID::kKaon,
+                              AliFlowTrackCuts::kTOFbetaSimple,
+                              1,3 );
+    AddTaskFlowCentralityPID( lowCentralityBinEdge,
+                              highCentralityBinEdge,
+                              commonOutputFileName,
+                              AliPID::kProton,
+                              AliFlowTrackCuts::kTOFbetaSimple,
+                              -1,3 );
+    AddTaskFlowCentralityPID( lowCentralityBinEdge,
+                              highCentralityBinEdge,
+                              commonOutputFileName,
+                              AliPID::kProton,
+                              AliFlowTrackCuts::kTOFbetaSimple,
+                              1,3 );
+  } // end of for (Int_t i=0; i<numberOfCentralityBins; i++)
+
+  // Enable debug printouts:
+  mgr->SetDebugLevel(2);
+  // Run the analysis:
+  if(!mgr->InitAnalysis()) return;
+  mgr->PrintStatus();
+  if(mode == mLocal)
+  {
+    mgr->StartAnalysis("local",chain);
+  }
+  else if(mode == mPROOF)
+  {
+    mgr->StartAnalysis("proof",dataDir,nEvents,offset);
+  }
+  else if(mode == mGrid)
+  {
+    mgr->StartAnalysis("grid");
+  }
+
+  // Print real and CPU time used for analysis:
+  timer.Stop();
+  timer.Print();
+
+} // end of void runFlowTaskCentralityPIDTrain(...)
+
+//===============================================================================================
+/*
+void CrossCheckUserSettings(Bool_t bData)
+{
+ // Check in this method if the user settings make sense.
+ if(LYZ1SUM && LYZ2SUM) {cout<<" WARNING: You cannot run LYZ1 and LYZ2 at the same time! LYZ2 needs the output from LYZ1 !!!!"<<endl; exit(0); }
+ if(LYZ1PROD && LYZ2PROD) {cout<<" WARNING: You cannot run LYZ1 and LYZ2 at the same time! LYZ2 needs the output from LYZ1 !!!!"<<endl; exit(0); }
+ if(LYZ2SUM && LYZEP) {cout<<" WARNING: You cannot run LYZ2 and LYZEP at the same time! LYZEP needs the output from LYZ2 !!!!"<<endl; exit(0); }
+ if(LYZ1SUM && LYZEP) {cout<<" WARNING: You cannot run LYZ1 and LYZEP at the same time! LYZEP needs the output from LYZ2 !!!!"<<endl; exit(0); }
+} // end of void CrossCheckUserSettings()
+*/
+//===============================================================================================
+
+void LoadLibraries(const anaModes mode, Bool_t useFlowParFiles )
+{
+  //--------------------------------------
+  // Load the needed libraries most of them already loaded by aliroot
+  //--------------------------------------
+
+  gSystem->Load("libCore");
+  gSystem->Load("libTree");
+  gSystem->Load("libGeom");
+  gSystem->Load("libVMC");
+  gSystem->Load("libXMLIO");
+  gSystem->Load("libPhysics");
+  gSystem->Load("libXMLParser");
+  gSystem->Load("libProof");
+  gSystem->Load("libMinuit");
+
+  if (mode==mLocal || mode==mGrid)
+  {
+    gSystem->Load("libSTEERBase");
+    gSystem->Load("libCDB");
+    gSystem->Load("libRAWDatabase");
+    gSystem->Load("libRAWDatarec");
+    gSystem->Load("libESD");
+    gSystem->Load("libAOD");
+    gSystem->Load("libSTEER");
+    gSystem->Load("libANALYSIS");
+    gSystem->Load("libANALYSISalice");
+    gSystem->Load("libTOFbase");
+    gSystem->Load("libTOFrec");
+    gSystem->Load("libTRDbase");
+    gSystem->Load("libVZERObase");
+    gSystem->Load("libVZEROrec");
+    gSystem->Load("libT0base");
+    gSystem->Load("libT0rec");
+    gSystem->Load("libTENDER");
+    gSystem->Load("libTENDERSupplies");
+
+    if (useFlowParFiles)
+    {
+      AliAnalysisAlien::SetupPar("PWG2flowCommon");
+      AliAnalysisAlien::SetupPar("PWG2flowTasks");
+    }
+    else
+    {
+      gSystem->Load("libPWG2flowCommon");
+      gSystem->Load("libPWG2flowTasks");
+    }
+  }
+  else if (mode==mPROOF)
+  {
+    TList* list = new TList();
+    list->Add(new TNamed("ALIROOT_MODE", "ALIROOT"));
+    if (useFlowParFiles)
+      list->Add(new TNamed("ALIROOT_EXTRA_LIBS", "ANALYSIS:ANALYSISalice:TENDER:TENDERSupplies"));
+    else
+      list->Add(new TNamed("ALIROOT_EXTRA_LIBS", "ANALYSIS:ANALYSISalice:TENDER:TENDERSupplies:PWG2flowCommon:PWG2flowTasks"));
+
+    //list->Add(new TNamed("ALIROOT_EXTRA_INCLUDES","PWG2/FLOW/AliFlowCommon:PWG2/FLOW/AliFlowTasks"));
+
+    // Connect to proof
+    printf("*** Connect to PROOF ***\n");
+    gEnv->SetValue("XSec.GSI.DelegProxy","2");
+    //TProof* proof = TProof::Open("alice-caf.cern.ch");
+    TProof* proof = TProof::Open("skaf.saske.sk");
+
+    // list the data available
+    //gProof->ShowDataSets("/*/*");
+    //gProof->ShowDataSets("/alice/sim/"); //for MC Data
+    //gProof->ShowDataSets("/alice/data/"); //for REAL Data
+
+    proof->ClearPackages();
+    proof->EnablePackage("VO_ALICE@AliRoot::v4-21-14-AN",list);
+
+    if (useFlowParFiles)
+    {
+      gProof->UploadPackage("PWG2flowCommon.par");
+      gProof->UploadPackage("PWG2flowTasks.par");
+    }
+
+    // Show enables Packages
+    gProof->ShowEnabledPackages();
+  }
+} // end of void LoadLibraries(const anaModes mode)
+
+//===============================================================================================
+
+TChain* CreateAODChain(const char* aDataDir, Int_t aRuns, Int_t offset)
+{
+  // creates chain of files in a given directory or file containing a list.
+  // In case of directory the structure is expected as:
+  // <aDataDir>/<dir0>/AliAOD.root
+  // <aDataDir>/<dir1>/AliAOD.root
+  // ...
+
+  if (!aDataDir)
+    return 0;
+
+  Long_t id, size, flags, modtime;
+  if (gSystem->GetPathInfo(aDataDir, &id, &size, &flags, &modtime))
+  {
+    printf("%s not found.\n", aDataDir);
+    return 0;
+  }
+
+  TChain* chain = new TChain("aodTree");
+  TChain* chaingAlice = 0;
+
+  if (flags & 2)
+  {
+    TString execDir(gSystem->pwd());
+    TSystemDirectory* baseDir = new TSystemDirectory(".", aDataDir);
+    TList* dirList            = baseDir->GetListOfFiles();
+    Int_t nDirs               = dirList->GetEntries();
+    gSystem->cd(execDir);
+
+    Int_t count = 0;
+
+    for (Int_t iDir=0; iDir<nDirs; ++iDir)
+    {
+      TSystemFile* presentDir = (TSystemFile*) dirList->At(iDir);
+      if (!presentDir || !presentDir->IsDirectory() || strcmp(presentDir->GetName(), ".") == 0 || strcmp(presentDir->GetName(), "..") == 0)
+        continue;
+
+      if (offset > 0)
+      {
+        --offset;
+        continue;
+      }
+
+      if (count++ == aRuns)
+        break;
+
+      TString presentDirName(aDataDir);
+      presentDirName += "/";
+      presentDirName += presentDir->GetName();
+      chain->Add(presentDirName + "/AliAOD.root/aodTree");
+      // cerr<<presentDirName<<endl;
+    }
+
+  }
+  else
+  {
+    // Open the input stream
+    ifstream in;
+    in.open(aDataDir);
+
+    Int_t count = 0;
+
+    // Read the input list of files and add them to the chain
+    TString aodfile;
+    while(in.good())
+    {
+      in >> aodfile;
+      if (!aodfile.Contains("root")) continue; // protection
+
+      if (offset > 0)
+      {
+        --offset;
+        continue;
+      }
+
+      if (count++ == aRuns)
+        break;
+
+      // add aod file
+      chain->Add(aodfile);
+    }
+
+    in.close();
+  }
+
+  return chain;
+
+} // end of TChain* CreateAODChain(const char* aDataDir, Int_t aRuns, Int_t offset)
+
index d22d28d..19d9ea9 100644 (file)
@@ -25,5 +25,6 @@
 #pragma link C++ class AliAnalysisTaskNestedLoops+;
 #pragma link C++ class AliAnalysisTaskQAflow+;
 #pragma link C++ class AliAnalysisTaskPIDflowQA+;
+#pragma link C++ class AliAnalysisTaskQAPmdflow+;
 
 #endif