[u/mrichter/AliRoot.git] / PWG2 / FLOW / macros / runAliAnalysisTaskFlow.C

// from CreateESDChain.C - instead of  #include "CreateESDChain.C"
TChain* CreateESDChain(const char* aDataDir = "ESDfiles.txt", Int_t aRuns = 20, Int_t offset = 0) ;
void LookupWrite(TChain* chain, const char* target) ;

//SETTING THE ANALYSIS

//Flow analysis method can be:
// SP    = Scalar Product
// LYZ1  = Lee Yang Zeroes first run
// LYZ2  = Lee Yang Zeroes second run
// LYZEP = Lee Yang Zeroes Event Plane
// CUM   = Cumulants
const TString method = "SP";

//Type of analysis can be:
// ESD, AOD, MC, ESDMC0, ESDMC1
const TString type = "ESDMC1";


//SETTING THE CUTS

//for integrated flow
const Double_t ptmin1 = 0.0;
const Double_t ptmax1 = 1000.0;
const Double_t ymin1  = -2.;
const Double_t ymax1  = 2.;
const Int_t mintrackrefsTPC1 = 2;
const Int_t mintrackrefsITS1 = 3;
const Int_t charge1 = 1;
const Int_t PDG1 = 211;
const Int_t minclustersTPC1 = 50;
const Int_t maxnsigmatovertex1 = 3;

//for differential flow
const Double_t ptmin2 = 0.0;
const Double_t ptmax2 = 1000.0;
const Double_t ymin2  = -2.;
const Double_t ymax2  = 2.;
const Int_t mintrackrefsTPC2 = 2;
const Int_t mintrackrefsITS2 = 3;
const Int_t charge2 = 1;
const Int_t PDG2 = 321;
const Int_t minclustersTPC2 = 50;
const Int_t maxnsigmatovertex2 = 3;


void runAliAnalysisTaskFlow(Int_t nRuns = 10, const Char_t* dataDir="/data/alice1/kolk/TherminatorFIX", Int_t offset = 0) 

{
  TStopwatch timer;
  timer.Start();

  // include path (to find the .h files when compiling)
  gSystem->AddIncludePath("-I$ALICE_ROOT/include") ;
  gSystem->AddIncludePath("-I$ROOTSYS/include") ;

  // load needed libraries
  gSystem->Load("libTree.so");
  gSystem->Load("libESD.so");
  cerr<<"libESD loaded..."<<endl;
  gSystem->Load("libANALYSIS.so");
  cerr<<"libANALYSIS.so loaded..."<<endl;
  gSystem->Load("libANALYSISRL.so");
  cerr<<"libANALYSISRL.so loaded..."<<endl;
  gSystem->Load("libCORRFW.so");
  cerr<<"libCORRFW.so loaded..."<<endl;
  gSystem->Load("libPWG2flow.so");
  cerr<<"libPWG2flow.so loaded..."<<endl;

  // create the TChain. CreateESDChain() is defined in CreateESDChain.C
  TChain* chain = CreateESDChain(dataDir, nRuns, offset);
  cout<<"chain ("<<chain<<")"<<endl;
 
  //____________________________________________//
  //Create cuts using correction framework

  //############# cuts on MC
  AliCFTrackKineCuts* mcKineCuts1 = new AliCFTrackKineCuts("mcKineCuts1","MC-level kinematic cuts");
  mcKineCuts1->SetPtRange(ptmin1,ptmax1);
  mcKineCuts1->SetRapidityRange(ymin1,ymax1);
  mcKineCuts1->SetChargeMC(charge1);
  
  AliCFTrackKineCuts* mcKineCuts2 = new AliCFTrackKineCuts("mcKineCuts2","MC-level kinematic cuts");
  mcKineCuts2->SetPtRange(ptmin2,ptmax2);
  mcKineCuts2->SetRapidityRange(ymin2,ymax2);
  mcKineCuts2->SetChargeMC(charge2);

  AliCFParticleGenCuts* mcGenCuts1 = new AliCFParticleGenCuts("mcGenCuts1","MC particle generation cuts");
  mcGenCuts1->SetRequireIsPrimary();
  mcGenCuts1->SetRequirePdgCode(PDG1);

  AliCFParticleGenCuts* mcGenCuts2 = new AliCFParticleGenCuts("mcGenCuts2","MC particle generation cuts");
  mcGenCuts2->SetRequireIsPrimary();
  mcGenCuts2->SetRequirePdgCode(PDG2);

  //############# Acceptance Cuts  ????????
  AliCFAcceptanceCuts *mcAccCuts = new AliCFAcceptanceCuts("mcAccCuts","MC acceptance cuts");
  mcAccCuts->SetMinNHitITS(mintrackrefsITS1);
  mcAccCuts->SetMinNHitTPC(mintrackrefsTPC1);
  
  //############# Rec-Level kinematic cuts
  AliCFTrackKineCuts *recKineCuts1 = new AliCFTrackKineCuts("recKineCuts1","rec-level kine cuts");
  recKineCuts1->SetPtRange(ptmin1,ptmax1);
  recKineCuts1->SetRapidityRange(ymin1,ymax1);
  recKineCuts1->SetChargeRec(charge1);

  AliCFTrackKineCuts *recKineCuts2 = new AliCFTrackKineCuts("recKineCuts2","rec-level kine cuts");
  recKineCuts2->SetPtRange(ptmin2,ptmax2);
  recKineCuts2->SetRapidityRange(ymin2,ymax2);
  recKineCuts2->SetChargeRec(charge2);

  AliCFTrackQualityCuts *recQualityCuts = new AliCFTrackQualityCuts("recQualityCuts","rec-level quality cuts");
  recQualityCuts->SetMinNClusterTPC(minclustersTPC1);
  recQualityCuts->SetRequireITSRefit(kTRUE);

  AliCFTrackIsPrimaryCuts *recIsPrimaryCuts = new AliCFTrackIsPrimaryCuts("recIsPrimaryCuts","rec-level isPrimary cuts");
  recIsPrimaryCuts->SetMaxNSigmaToVertex(maxnsigmatovertex1);
  
  AliCFTrackCutPid* cutPID1 = new AliCFTrackCutPid("cutPID1","ESD_PID") ;
  AliCFTrackCutPid* cutPID2 = new AliCFTrackCutPid("cutPID2","ESD_PID") ;
  int n_species = AliPID::kSPECIES ;
  Double_t* prior = new Double_t[n_species];

  prior[0] = 0.0244519 ;
  prior[1] = 0.0143988 ;
  prior[2] = 0.805747  ;
  prior[3] = 0.0928785 ;
  prior[4] = 0.0625243 ;
  
  cutPID1->SetPriors(prior);
  cutPID1->SetProbabilityCut(0.0);
  cutPID1->SetDetectors("TPC TOF");
  switch(TMath::Abs(PDG1)) {
  case 11   : cutPID1->SetParticleType(AliPID::kElectron, kTRUE); break;
  case 13   : cutPID1->SetParticleType(AliPID::kMuon    , kTRUE); break;
  case 211  : cutPID1->SetParticleType(AliPID::kPion    , kTRUE); break;
  case 321  : cutPID1->SetParticleType(AliPID::kKaon    , kTRUE); break;
  case 2212 : cutPID1->SetParticleType(AliPID::kProton  , kTRUE); break;
  default   : printf("UNDEFINED PID\n"); break;
  }

  cutPID2->SetPriors(prior);
  cutPID2->SetProbabilityCut(0.0);
  cutPID2->SetDetectors("TPC TOF");
  switch(TMath::Abs(PDG2)) {
  case 11   : cutPID2->SetParticleType(AliPID::kElectron, kTRUE); break;
  case 13   : cutPID2->SetParticleType(AliPID::kMuon    , kTRUE); break;
  case 211  : cutPID2->SetParticleType(AliPID::kPion    , kTRUE); break;
  case 321  : cutPID2->SetParticleType(AliPID::kKaon    , kTRUE); break;
  case 2212 : cutPID2->SetParticleType(AliPID::kProton  , kTRUE); break;
  default   : printf("UNDEFINED PID\n"); break;
  }


  printf("CREATE MC KINE CUTS\n");
  TObjArray* mcList1 = new TObjArray(0);
  mcList1->AddLast(mcKineCuts1);
  mcList1->AddLast(mcGenCuts1);
  
  TObjArray* mcList2 = new TObjArray(0);
  mcList2->AddLast(mcKineCuts2);
  mcList2->AddLast(mcGenCuts2);

  printf("CREATE ACCEPTANCE CUTS\n");
  TObjArray* accList = new TObjArray(0) ;
  accList->AddLast(mcAccCuts);

  printf("CREATE RECONSTRUCTION CUTS\n");
  TObjArray* recList1 = new TObjArray(0) ;
  recList1->AddLast(recKineCuts1);
  recList1->AddLast(recQualityCuts);
  recList1->AddLast(recIsPrimaryCuts);

  TObjArray* recList2 = new TObjArray(0) ;
  recList2->AddLast(recKineCuts2);
  recList2->AddLast(recQualityCuts);
  recList2->AddLast(recIsPrimaryCuts);

  printf("CREATE PID CUTS\n");
  TObjArray* fPIDCutList1 = new TObjArray(0) ;
  fPIDCutList1->AddLast(cutPID1);

  TObjArray* fPIDCutList2 = new TObjArray(0) ;
  fPIDCutList2->AddLast(cutPID2);

  printf("CREATE INTERFACE AND CUTS\n");
  AliCFManager* cfmgr1 = new AliCFManager();
  cfmgr1->SetParticleCutsList(AliCFManager::kPartGenCuts,mcList1);
  //cfmgr1->SetParticleCutsList(AliCFManager::kPartAccCuts,accList);
  cfmgr1->SetParticleCutsList(AliCFManager::kPartRecCuts,recList1);
  cfmgr1->SetParticleCutsList(AliCFManager::kPartSelCuts,fPIDCutList1);

  AliCFManager* cfmgr2 = new AliCFManager();
  cfmgr2->SetParticleCutsList(AliCFManager::kPartGenCuts,mcList2);
  //cfmgr2->SetParticleCutsList(AliCFManager::kPartAccCuts,accList);
  cfmgr2->SetParticleCutsList(AliCFManager::kPartRecCuts,recList2);
  cfmgr2->SetParticleCutsList(AliCFManager::kPartSelCuts,fPIDCutList2);

 
  if (method == "LYZ2"){  
    // read the input files from the first run 
    TString firstRunFileName = "outputLYZ1analysis" ;
    firstRunFileName += type;
    firstRunFileName += "_firstrun.root";
    cout<<"The input file is "<<firstRunFileName.Data()<<endl;
    TFile* fFirstRunFile = new TFile(firstRunFileName.Data(),"READ");
    if(!fFirstRunFile || fFirstRunFile->IsZombie()) { cerr << " ERROR: NO First Run file... " << endl ; }
    cout<<"input files read..."<<endl;
  }


  //____________________________________________//
  // Make the analysis manager
  AliAnalysisManager *mgr = new AliAnalysisManager("FlowAnalysisManager");
  
  if (type == "ESD"){
    AliVEventHandler* esdH = new AliESDInputHandler;
    mgr->SetInputEventHandler(esdH); }
   
  if (type == "AOD"){
    AliVEventHandler* aodH = new AliAODInputHandler;
    mgr->SetInputEventHandler(aodH); }
  
  if (type == "MC" || type == "ESDMC0" || type == "ESDMC1"){
    AliVEventHandler* esdH = new AliESDInputHandler;
    mgr->SetInputEventHandler(esdH);

    AliMCEventHandler *mc = new AliMCEventHandler();
    mgr->SetMCtruthEventHandler(mc); }

  //____________________________________________//
  // 1st MC EP task
  if (method == "SP"){
    AliAnalysisTaskScalarProduct *task1 = new AliAnalysisTaskScalarProduct("TaskScalarProduct");
    task1->SetAnalysisType(type);
    task1->SetCFManager1(cfmgr1);
    task1->SetCFManager2(cfmgr2);
    mgr->AddTask(task1);
  }
  else if (method == "LYZ1"){
    AliAnalysisTaskLeeYangZeros *task1 = new AliAnalysisTaskLeeYangZeros("TaskLeeYangZeros",kTRUE);
    task1->SetAnalysisType(type);
    task1->SetFirstRunLYZ(kTRUE);
    task1->SetUseSumLYZ(kTRUE);
    task1->SetCFManager1(cfmgr1);
    task1->SetCFManager2(cfmgr2);
    mgr->AddTask(task1);
  }
  else if (method == "LYZ2"){
    AliAnalysisTaskLeeYangZeros *task1 = new AliAnalysisTaskLeeYangZeros("TaskLeeYangZeros",kFALSE);
    task1->SetAnalysisType(type);
    task1->SetFirstRunLYZ(kFALSE);
    task1->SetUseSumLYZ(kTRUE);
    task1->SetCFManager1(cfmgr1);
    task1->SetCFManager2(cfmgr2);
    mgr->AddTask(task1);
  }
  else if (method == "LYZEP"){
    AliAnalysisTaskLYZEventPlane *task1 = new AliAnalysisTaskLYZEventPlane("TaskLYZEventPlane");
    task1->SetAnalysisType(type);
    //task1->SetCFManager1(cfmgr1);
    //task1->SetCFManager2(cfmgr2);
    mgr->AddTask(task1);
  }
  else if (method == "CUM"){
    AliAnalysisTaskCumulants *task1 = new AliAnalysisTaskCumulants("TaskCumulants");
    task1->SetAnalysisType(type);
    //task1->SetCFManager1(cfmgr1);
    //task1->SetCFManager2(cfmgr2);
    mgr->AddTask(task1);
  }


  // Create containers for input/output
  AliAnalysisDataContainer *cinput1 = 
    mgr->CreateContainer("cchain1",TChain::Class(),AliAnalysisManager::kInputContainer);

  if (method == "LYZ2"){ AliAnalysisDataContainer *cinput2 = 
		    mgr->CreateContainer("cobj2",TList::Class(),AliAnalysisManager::kInputContainer); } 


  TString outputName = "output";
  outputName+= method;
  outputName+= "analysis";
  outputName+= type;
  if (method == "LYZ1") outputName+= "_firstrun";
  if (method == "LYZ2") outputName+= "_secondrun";
  outputName+= ".root";
  AliAnalysisDataContainer *coutput1 = 
    mgr->CreateContainer("cobj1", TList::Class(),AliAnalysisManager::kOutputContainer,outputName);

  //____________________________________________//
  mgr->ConnectInput(task1,0,cinput1);
  if (method == "LYZ2") { mgr->ConnectInput(task1,1,cinput2); }
  mgr->ConnectOutput(task1,0,coutput1);

  if (method == "LYZ2"){ cinput2->SetData(fFirstRunFile);}

  if (!mgr->InitAnalysis()) return;
  mgr->PrintStatus();
  mgr->StartAnalysis("local",chain);

  timer.Stop();
  timer.Print();
}


// Helper macros for creating chains
// from: CreateESDChain.C,v 1.10 jgrosseo Exp

TChain* CreateESDChain(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>/AliESDs.root
  // <aDataDir>/<dir1>/AliESDs.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("esdTree");
  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 + "/AliESDs.root/esdTree");
	  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 esdfile;
      while(in.good()) {
	in >> esdfile;
	if (!esdfile.Contains("root")) continue; // protection
	
	if (offset > 0)
	  {
	    --offset;
	    continue;
	  }
	
	if (count++ == aRuns)
	  break;
	
        // add esd file
	chain->Add(esdfile);
      }
      
      in.close();
    }
  
  return chain;
}

void LookupWrite(TChain* chain, const char* target)
{
  // looks up the chain and writes the remaining files to the text file target
  
  chain->Lookup();
  
  TObjArray* list = chain->GetListOfFiles();
  TIterator* iter = list->MakeIterator();
  TObject* obj = 0;
  
  ofstream outfile;
  outfile.open(target);
  
  while ((obj = iter->Next()))
    outfile << obj->GetTitle() << "#AliESDs.root" << endl;
  
  outfile.close();
  
  delete iter;
}
Commit	Line	Data
448e8856	1	// from CreateESDChain.C - instead of #include "CreateESDChain.C"
	2	TChain* CreateESDChain(const char* aDataDir = "ESDfiles.txt", Int_t aRuns = 20, Int_t offset = 0) ;
	3	void LookupWrite(TChain* chain, const char* target) ;
	4
88e00a8a	5	//SETTING THE ANALYSIS
448e8856	6
88e00a8a	7	//Flow analysis method can be:
	8	// SP = Scalar Product
	9	// LYZ1 = Lee Yang Zeroes first run
	10	// LYZ2 = Lee Yang Zeroes second run
	11	// LYZEP = Lee Yang Zeroes Event Plane
	12	// CUM = Cumulants
	13	const TString method = "SP";
	14
	15	//Type of analysis can be:
	16	// ESD, AOD, MC, ESDMC0, ESDMC1
	17	const TString type = "ESDMC1";
	18
	19
	20	//SETTING THE CUTS
	21
	22	//for integrated flow
	23	const Double_t ptmin1 = 0.0;
	24	const Double_t ptmax1 = 1000.0;
	25	const Double_t ymin1 = -2.;
	26	const Double_t ymax1 = 2.;
	27	const Int_t mintrackrefsTPC1 = 2;
	28	const Int_t mintrackrefsITS1 = 3;
	29	const Int_t charge1 = 1;
	30	const Int_t PDG1 = 211;
	31	const Int_t minclustersTPC1 = 50;
	32	const Int_t maxnsigmatovertex1 = 3;
	33
	34	//for differential flow
	35	const Double_t ptmin2 = 0.0;
	36	const Double_t ptmax2 = 1000.0;
	37	const Double_t ymin2 = -2.;
	38	const Double_t ymax2 = 2.;
	39	const Int_t mintrackrefsTPC2 = 2;
	40	const Int_t mintrackrefsITS2 = 3;
	41	const Int_t charge2 = 1;
	42	const Int_t PDG2 = 321;
	43	const Int_t minclustersTPC2 = 50;
	44	const Int_t maxnsigmatovertex2 = 3;
	45
	46
	47
	48	void runAliAnalysisTaskFlow(Int_t nRuns = 10, const Char_t* dataDir="/data/alice1/kolk/TherminatorFIX", Int_t offset = 0)
448e8856	49
448e8856	50	{
	51	TStopwatch timer;
	52	timer.Start();
	53
88e00a8a	54	// include path (to find the .h files when compiling)
	55	gSystem->AddIncludePath("-I$ALICE_ROOT/include") ;
	56	gSystem->AddIncludePath("-I$ROOTSYS/include") ;
	57
448e8856	58	// load needed libraries
	59	gSystem->Load("libTree.so");
	60	gSystem->Load("libESD.so");
	61	cerr<<"libESD loaded..."<<endl;
	62	gSystem->Load("libANALYSIS.so");
	63	cerr<<"libANALYSIS.so loaded..."<<endl;
88e00a8a	64	gSystem->Load("libANALYSISRL.so");
	65	cerr<<"libANALYSISRL.so loaded..."<<endl;
	66	gSystem->Load("libCORRFW.so");
	67	cerr<<"libCORRFW.so loaded..."<<endl;
448e8856	68	gSystem->Load("libPWG2flow.so");
88e00a8a	69	cerr<<"libPWG2flow.so loaded..."<<endl;
88e00a8a	70
448e8856	71	// create the TChain. CreateESDChain() is defined in CreateESDChain.C
	72	TChain* chain = CreateESDChain(dataDir, nRuns, offset);
	73	cout<<"chain ("<<chain<<")"<<endl;
88e00a8a	74
	75	//____________________________________________//
	76	//Create cuts using correction framework
	77
	78	//############# cuts on MC
	79	AliCFTrackKineCuts* mcKineCuts1 = new AliCFTrackKineCuts("mcKineCuts1","MC-level kinematic cuts");
	80	mcKineCuts1->SetPtRange(ptmin1,ptmax1);
	81	mcKineCuts1->SetRapidityRange(ymin1,ymax1);
	82	mcKineCuts1->SetChargeMC(charge1);
	83
	84	AliCFTrackKineCuts* mcKineCuts2 = new AliCFTrackKineCuts("mcKineCuts2","MC-level kinematic cuts");
	85	mcKineCuts2->SetPtRange(ptmin2,ptmax2);
	86	mcKineCuts2->SetRapidityRange(ymin2,ymax2);
	87	mcKineCuts2->SetChargeMC(charge2);
	88
	89	AliCFParticleGenCuts* mcGenCuts1 = new AliCFParticleGenCuts("mcGenCuts1","MC particle generation cuts");
	90	mcGenCuts1->SetRequireIsPrimary();
	91	mcGenCuts1->SetRequirePdgCode(PDG1);
	92
	93	AliCFParticleGenCuts* mcGenCuts2 = new AliCFParticleGenCuts("mcGenCuts2","MC particle generation cuts");
	94	mcGenCuts2->SetRequireIsPrimary();
	95	mcGenCuts2->SetRequirePdgCode(PDG2);
	96
	97	//############# Acceptance Cuts ????????
	98	AliCFAcceptanceCuts *mcAccCuts = new AliCFAcceptanceCuts("mcAccCuts","MC acceptance cuts");
	99	mcAccCuts->SetMinNHitITS(mintrackrefsITS1);
	100	mcAccCuts->SetMinNHitTPC(mintrackrefsTPC1);
	101
	102	//############# Rec-Level kinematic cuts
	103	AliCFTrackKineCuts *recKineCuts1 = new AliCFTrackKineCuts("recKineCuts1","rec-level kine cuts");
	104	recKineCuts1->SetPtRange(ptmin1,ptmax1);
	105	recKineCuts1->SetRapidityRange(ymin1,ymax1);
	106	recKineCuts1->SetChargeRec(charge1);
	107
	108	AliCFTrackKineCuts *recKineCuts2 = new AliCFTrackKineCuts("recKineCuts2","rec-level kine cuts");
	109	recKineCuts2->SetPtRange(ptmin2,ptmax2);
	110	recKineCuts2->SetRapidityRange(ymin2,ymax2);
	111	recKineCuts2->SetChargeRec(charge2);
	112
	113	AliCFTrackQualityCuts *recQualityCuts = new AliCFTrackQualityCuts("recQualityCuts","rec-level quality cuts");
	114	recQualityCuts->SetMinNClusterTPC(minclustersTPC1);
	115	recQualityCuts->SetRequireITSRefit(kTRUE);
	116
	117	AliCFTrackIsPrimaryCuts *recIsPrimaryCuts = new AliCFTrackIsPrimaryCuts("recIsPrimaryCuts","rec-level isPrimary cuts");
	118	recIsPrimaryCuts->SetMaxNSigmaToVertex(maxnsigmatovertex1);
	119
	120	AliCFTrackCutPid* cutPID1 = new AliCFTrackCutPid("cutPID1","ESD_PID") ;
	121	AliCFTrackCutPid* cutPID2 = new AliCFTrackCutPid("cutPID2","ESD_PID") ;
	122	int n_species = AliPID::kSPECIES ;
	123	Double_t* prior = new Double_t[n_species];
	124
	125	prior[0] = 0.0244519 ;
	126	prior[1] = 0.0143988 ;
	127	prior[2] = 0.805747 ;
	128	prior[3] = 0.0928785 ;
	129	prior[4] = 0.0625243 ;
	130
	131	cutPID1->SetPriors(prior);
	132	cutPID1->SetProbabilityCut(0.0);
	133	cutPID1->SetDetectors("TPC TOF");
	134	switch(TMath::Abs(PDG1)) {
	135	case 11 : cutPID1->SetParticleType(AliPID::kElectron, kTRUE); break;
	136	case 13 : cutPID1->SetParticleType(AliPID::kMuon , kTRUE); break;
	137	case 211 : cutPID1->SetParticleType(AliPID::kPion , kTRUE); break;
138	case 321 : cutPID1->SetParticleType(AliPID::kKaon , kTRUE); break;
139	case 2212 : cutPID1->SetParticleType(AliPID::kProton , kTRUE); break;
140	default : printf("UNDEFINED PID\n"); break;
141	}
142
143	cutPID2->SetPriors(prior);
144	cutPID2->SetProbabilityCut(0.0);
145	cutPID2->SetDetectors("TPC TOF");
146	switch(TMath::Abs(PDG2)) {
147	case 11 : cutPID2->SetParticleType(AliPID::kElectron, kTRUE); break;
148	case 13 : cutPID2->SetParticleType(AliPID::kMuon , kTRUE); break;
149	case 211 : cutPID2->SetParticleType(AliPID::kPion , kTRUE); break;
150	case 321 : cutPID2->SetParticleType(AliPID::kKaon , kTRUE); break;
151	case 2212 : cutPID2->SetParticleType(AliPID::kProton , kTRUE); break;
152	default : printf("UNDEFINED PID\n"); break;
153	}
154
155
156	printf("CREATE MC KINE CUTS\n");
157	TObjArray* mcList1 = new TObjArray(0);
158	mcList1->AddLast(mcKineCuts1);
159	mcList1->AddLast(mcGenCuts1);
160
161	TObjArray* mcList2 = new TObjArray(0);
162	mcList2->AddLast(mcKineCuts2);
163	mcList2->AddLast(mcGenCuts2);
164
165	printf("CREATE ACCEPTANCE CUTS\n");
166	TObjArray* accList = new TObjArray(0) ;
167	accList->AddLast(mcAccCuts);
168
169	printf("CREATE RECONSTRUCTION CUTS\n");
170	TObjArray* recList1 = new TObjArray(0) ;
171	recList1->AddLast(recKineCuts1);
172	recList1->AddLast(recQualityCuts);
173	recList1->AddLast(recIsPrimaryCuts);
174
175	TObjArray* recList2 = new TObjArray(0) ;
176	recList2->AddLast(recKineCuts2);
177	recList2->AddLast(recQualityCuts);
178	recList2->AddLast(recIsPrimaryCuts);
179
180	printf("CREATE PID CUTS\n");
181	TObjArray* fPIDCutList1 = new TObjArray(0) ;
182	fPIDCutList1->AddLast(cutPID1);
183
184	TObjArray* fPIDCutList2 = new TObjArray(0) ;
185	fPIDCutList2->AddLast(cutPID2);
448e8856	186
88e00a8a	187	printf("CREATE INTERFACE AND CUTS\n");
	188	AliCFManager* cfmgr1 = new AliCFManager();
	189	cfmgr1->SetParticleCutsList(AliCFManager::kPartGenCuts,mcList1);
	190	//cfmgr1->SetParticleCutsList(AliCFManager::kPartAccCuts,accList);
	191	cfmgr1->SetParticleCutsList(AliCFManager::kPartRecCuts,recList1);
	192	cfmgr1->SetParticleCutsList(AliCFManager::kPartSelCuts,fPIDCutList1);
	193
	194	AliCFManager* cfmgr2 = new AliCFManager();
	195	cfmgr2->SetParticleCutsList(AliCFManager::kPartGenCuts,mcList2);
	196	//cfmgr2->SetParticleCutsList(AliCFManager::kPartAccCuts,accList);
	197	cfmgr2->SetParticleCutsList(AliCFManager::kPartRecCuts,recList2);
	198	cfmgr2->SetParticleCutsList(AliCFManager::kPartSelCuts,fPIDCutList2);
	199
	200
	201	if (method == "LYZ2"){
	202	// read the input files from the first run
	203	TString firstRunFileName = "outputLYZ1analysis" ;
448e8856	204	firstRunFileName += type;
448e8856	205	firstRunFileName += "_firstrun.root";
88e00a8a	206	cout<<"The input file is "<<firstRunFileName.Data()<<endl;
448e8856	207	TFile* fFirstRunFile = new TFile(firstRunFileName.Data(),"READ");
	208	if(!fFirstRunFile \|\| fFirstRunFile->IsZombie()) { cerr << " ERROR: NO First Run file... " << endl ; }
	209	cout<<"input files read..."<<endl;
	210	}
	211
448e8856	212
	213	//____________________________________________//
	214	// Make the analysis manager
88e00a8a	215	AliAnalysisManager *mgr = new AliAnalysisManager("FlowAnalysisManager");
	216
	217	if (type == "ESD"){
	218	AliVEventHandler* esdH = new AliESDInputHandler;
	219	mgr->SetInputEventHandler(esdH); }
	220
	221	if (type == "AOD"){
	222	AliVEventHandler* aodH = new AliAODInputHandler;
	223	mgr->SetInputEventHandler(aodH); }
	224
	225	if (type == "MC" \|\| type == "ESDMC0" \|\| type == "ESDMC1"){
	226	AliVEventHandler* esdH = new AliESDInputHandler;
	227	mgr->SetInputEventHandler(esdH);
448e8856	228
88e00a8a	229	AliMCEventHandler *mc = new AliMCEventHandler();
	230	mgr->SetMCtruthEventHandler(mc); }
	231
	232	//____________________________________________//
	233	// 1st MC EP task
	234	if (method == "SP"){
	235	AliAnalysisTaskScalarProduct *task1 = new AliAnalysisTaskScalarProduct("TaskScalarProduct");
	236	task1->SetAnalysisType(type);
	237	task1->SetCFManager1(cfmgr1);
	238	task1->SetCFManager2(cfmgr2);
	239	mgr->AddTask(task1);
	240	}
	241	else if (method == "LYZ1"){
	242	AliAnalysisTaskLeeYangZeros *task1 = new AliAnalysisTaskLeeYangZeros("TaskLeeYangZeros",kTRUE);
	243	task1->SetAnalysisType(type);
	244	task1->SetFirstRunLYZ(kTRUE);
	245	task1->SetUseSumLYZ(kTRUE);
	246	task1->SetCFManager1(cfmgr1);
	247	task1->SetCFManager2(cfmgr2);
	248	mgr->AddTask(task1);
	249	}
	250	else if (method == "LYZ2"){
	251	AliAnalysisTaskLeeYangZeros *task1 = new AliAnalysisTaskLeeYangZeros("TaskLeeYangZeros",kFALSE);
	252	task1->SetAnalysisType(type);
	253	task1->SetFirstRunLYZ(kFALSE);
	254	task1->SetUseSumLYZ(kTRUE);
	255	task1->SetCFManager1(cfmgr1);
	256	task1->SetCFManager2(cfmgr2);
	257	mgr->AddTask(task1);
	258	}
	259	else if (method == "LYZEP"){
	260	AliAnalysisTaskLYZEventPlane *task1 = new AliAnalysisTaskLYZEventPlane("TaskLYZEventPlane");
	261	task1->SetAnalysisType(type);
	262	//task1->SetCFManager1(cfmgr1);
	263	//task1->SetCFManager2(cfmgr2);
	264	mgr->AddTask(task1);
	265	}
	266	else if (method == "CUM"){
	267	AliAnalysisTaskCumulants *task1 = new AliAnalysisTaskCumulants("TaskCumulants");
	268	task1->SetAnalysisType(type);
	269	//task1->SetCFManager1(cfmgr1);
	270	//task1->SetCFManager2(cfmgr2);
	271	mgr->AddTask(task1);
	272	}
448e8856	273
448e8856	274
	275	// Create containers for input/output
	276	AliAnalysisDataContainer *cinput1 =
	277	mgr->CreateContainer("cchain1",TChain::Class(),AliAnalysisManager::kInputContainer);
448e8856	278
88e00a8a	279	if (method == "LYZ2"){ AliAnalysisDataContainer *cinput2 =
88e00a8a	280	mgr->CreateContainer("cobj2",TList::Class(),AliAnalysisManager::kInputContainer); }
448e8856	281
88e00a8a	282
	283	TString outputName = "output";
	284	outputName+= method;
	285	outputName+= "analysis";
	286	outputName+= type;
	287	if (method == "LYZ1") outputName+= "_firstrun";
	288	if (method == "LYZ2") outputName+= "_secondrun";
	289	outputName+= ".root";
448e8856	290	AliAnalysisDataContainer *coutput1 =
88e00a8a	291	mgr->CreateContainer("cobj1", TList::Class(),AliAnalysisManager::kOutputContainer,outputName);
448e8856	292
448e8856	293	//____________________________________________//
448e8856	294	mgr->ConnectInput(task1,0,cinput1);
88e00a8a	295	if (method == "LYZ2") { mgr->ConnectInput(task1,1,cinput2); }
448e8856	296	mgr->ConnectOutput(task1,0,coutput1);
448e8856	297
88e00a8a	298	if (method == "LYZ2"){ cinput2->SetData(fFirstRunFile);}
448e8856	299
	300	if (!mgr->InitAnalysis()) return;
	301	mgr->PrintStatus();
	302	mgr->StartAnalysis("local",chain);
	303
	304	timer.Stop();
	305	timer.Print();
	306	}
	307
	308
	309	// Helper macros for creating chains
	310	// from: CreateESDChain.C,v 1.10 jgrosseo Exp
	311
	312	TChain* CreateESDChain(const char* aDataDir, Int_t aRuns, Int_t offset)
	313	{
	314	// creates chain of files in a given directory or file containing a list.
	315	// In case of directory the structure is expected as:
	316	// <aDataDir>/<dir0>/AliESDs.root
	317	// <aDataDir>/<dir1>/AliESDs.root
	318	// ...
	319
	320	if (!aDataDir)
	321	return 0;
	322
	323	Long_t id, size, flags, modtime;
	324	if (gSystem->GetPathInfo(aDataDir, &id, &size, &flags, &modtime))
	325	{
	326	printf("%s not found.\n", aDataDir);
	327	return 0;
	328	}
	329
	330	TChain* chain = new TChain("esdTree");
	331	TChain* chaingAlice = 0;
	332
	333	if (flags & 2)
	334	{
	335	TString execDir(gSystem->pwd());
	336	TSystemDirectory* baseDir = new TSystemDirectory(".", aDataDir);
	337	TList* dirList = baseDir->GetListOfFiles();
	338	Int_t nDirs = dirList->GetEntries();
	339	gSystem->cd(execDir);
	340
	341	Int_t count = 0;
	342
	343	for (Int_t iDir=0; iDir<nDirs; ++iDir)
	344	{
	345	TSystemFile* presentDir = (TSystemFile*) dirList->At(iDir);
	346	if (!presentDir \|\| !presentDir->IsDirectory() \|\| strcmp(presentDir->GetName(), ".") == 0 \|\| strcmp(presentDir->GetName(), "..") == 0)
	347	continue;
	348
	349	if (offset > 0)
	350	{
	351	--offset;
	352	continue;
	353	}
	354
	355	if (count++ == aRuns)
	356	break;
	357
	358	TString presentDirName(aDataDir);
	359	presentDirName += "/";
	360	presentDirName += presentDir->GetName();
	361	chain->Add(presentDirName + "/AliESDs.root/esdTree");
	362	cerr<<presentDirName<<endl;
363	}
364
365	}
366	else
367	{
368	// Open the input stream
369	ifstream in;
370	in.open(aDataDir);
371
372	Int_t count = 0;
373
374	// Read the input list of files and add them to the chain
375	TString esdfile;
376	while(in.good()) {
377	in >> esdfile;
378	if (!esdfile.Contains("root")) continue; // protection
379
380	if (offset > 0)
381	{
382	--offset;
383	continue;
384	}
385
386	if (count++ == aRuns)
387	break;
388
389	// add esd file
390	chain->Add(esdfile);
391	}
392
393	in.close();
394	}
395
396	return chain;
397	}
398
399	void LookupWrite(TChain* chain, const char* target)
400	{
401	// looks up the chain and writes the remaining files to the text file target
402
403	chain->Lookup();
404
405	TObjArray* list = chain->GetListOfFiles();
406	TIterator* iter = list->MakeIterator();
407	TObject* obj = 0;
408
409	ofstream outfile;
410	outfile.open(target);
411
412	while ((obj = iter->Next()))
413	outfile << obj->GetTitle() << "#AliESDs.root" << endl;
414
415	outfile.close();
416
417	delete iter;
418	}