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

#include "TStopwatch.h"
#include "TObjArray"
#include "Riostream.h"
#include "TFile.h"

//--------------------------------------------------------------------------------------
// RUN SETTINGS
//flow analysis method can be: (set to kTRUE or kFALSE)
Bool_t SP    = kTRUE;
Bool_t LYZ1  = kTRUE;
Bool_t LYZ2  = kFALSE;  
Bool_t LYZEP = kFALSE; 
Bool_t GFC   = kTRUE;
Bool_t QC    = kTRUE;
Bool_t FQD   = kTRUE;
Bool_t MCEP  = kFALSE; //does not work yet 24/12/08
//--------------------------------------------------------------------------------------

//--------------------------------------------------------------------------------------
// CUT SETTINGS
//integrated selection
Double_t ptMaxInt  = 10.;
Double_t ptMinInt  = 0.;
Double_t etaMaxInt = 1.;
Double_t etaMinInt = -1.;
Double_t phiMaxInt = 7.5;
Double_t phiMinInt = 0.;
Int_t PIDInt       = 211;

//differential selection
Double_t ptMaxDiff  = 10.;
Double_t ptMinDiff  = 0.;
Double_t etaMaxDiff = 1.;
Double_t etaMinDiff = -1.;
Double_t phiMaxDiff = 7.5;
Double_t phiMinDiff = 0.;
Int_t PIDDiff       = 211;
//--------------------------------------------------------------------------------------

//--------------------------------------------------------------------------------------
// WEIGHTS SETTINGS 
//to use or not to use the weights - that is a question!
Bool_t usePhiWeights = kFALSE; //Phi (correction for non-uniform azimuthal acceptance)
Bool_t usePtWeights  = kFALSE; //v'(pt) (differential flow in pt)
Bool_t useEtaWeights = kFALSE; //v'(eta) (differential flow in eta)
//--------------------------------------------------------------------------------------

Int_t offset = 0;
                                          
int runFlowAnalysis(Int_t aRuns = 100, const char* 
		    // 			  dir="/data/alice1/kolk/KineOnly3/")
		    dir="/Users/snelling/alice_data/KineOnly3/")
{
  TStopwatch timer;
  timer.Start();
  
  if (LYZ1 && LYZ2) {cout<<"WARNING: you cannot run LYZ1 and LYZ2 at the same time! LYZ2 needs the output from LYZ1."<<endl; exit(); }
  if (LYZ2 && LYZEP) {cout<<"WARNING: you cannot run LYZ2 and LYZEP at the same time! LYZEP needs the output from LYZ2."<<endl; exit(); }
  if (LYZ1 && LYZEP) {cout<<"WARNING: you cannot run LYZ1 and LYZEP at the same time! LYZEP needs the output from LYZ2."<<endl; exit(); }


  cout<<endl;
  cout<<" ---- BEGIN ANALYSIS ---- "<<endl;
  cout<<endl;
  

  gSystem->AddIncludePath("-I$ROOTSYS/include");

  // In AliRoot
  gSystem->AddIncludePath("-I$ALICE_ROOT/include");

  gSystem->Load("libANALYSIS.so");
  gSystem->Load("libPWG2flowCommon.so");
  cerr<<"libPWG2flowCommon.so loaded ..."<<endl;
  gSystem->Load("libPWG2flowTasks.so");
  cerr<<"libPWG2flowTasks.so loaded ..."<<endl;
  cout<<endl; 
  
  // Flow event in AliRoot
  AliFlowEventSimpleMaker* fEventMaker = new AliFlowEventSimpleMaker();
   
  // In root

  /*    
  // Constants  
  gROOT->LoadMacro("AliFlowCommon/AliFlowCommonConstants.cxx+");
  gROOT->LoadMacro("AliFlowCommon/AliFlowLYZConstants.cxx+");
  gROOT->LoadMacro("AliFlowCommon/AliFlowCumuConstants.cxx+");

  // Flow event
  gROOT->LoadMacro("AliFlowCommon/AliFlowVector.cxx+"); 
  gROOT->LoadMacro("AliFlowCommon/AliFlowTrackSimple.cxx+");    
  gROOT->LoadMacro("AliFlowCommon/AliFlowEventSimple.cxx+");

  // Cuts
  gROOT->LoadMacro("AliFlowCommon/AliFlowTrackSimpleCuts.cxx+");    

  // Output histosgrams
  gROOT->LoadMacro("AliFlowCommon/AliFlowCommonHist.cxx+");
  gROOT->LoadMacro("AliFlowCommon/AliFlowCommonHistResults.cxx+");
  gROOT->LoadMacro("AliFlowCommon/AliFlowLYZHist1.cxx+");
  gROOT->LoadMacro("AliFlowCommon/AliFlowLYZHist2.cxx+");

  // Functions needed for various methods
  gROOT->LoadMacro("AliFlowCommon/AliCumulantsFunctions.cxx+");
  gROOT->LoadMacro("AliFlowCommon/AliQCumulantsFunctions.cxx+");
  gROOT->LoadMacro("AliFlowCommon/AliFittingFunctionsForQDistribution.cxx+");
  gROOT->LoadMacro("AliFlowCommon/AliFlowLYZEventPlane.cxx+");

  // Flow Analysis code for various methods
  gROOT->LoadMacro("AliFlowCommon/AliFlowAnalysisWithMCEventPlane.cxx+"); 
  gROOT->LoadMacro("AliFlowCommon/AliFlowAnalysisWithScalarProduct.cxx+");
  gROOT->LoadMacro("AliFlowCommon/AliFlowAnalysisWithLYZEventPlane.cxx+");
  gROOT->LoadMacro("AliFlowCommon/AliFlowAnalysisWithLeeYangZeros.cxx+");
  gROOT->LoadMacro("AliFlowCommon/AliFlowAnalysisWithCumulants.cxx+");
  gROOT->LoadMacro("AliFlowCommon/AliFlowAnalysisWithQCumulants.cxx+"); 
  gROOT->LoadMacro("AliFlowCommon/AliFittingQDistribution.cxx+");

  // Class to fill the FlowEvent without aliroot dependence
  // can be found in the directory FlowEventMakers
  gROOT->LoadMacro("FlowEventMakers/FlowEventSimpleMaker.cxx+");   

  cout << "finished loading macros!" << endl;  

  //flow event
  FlowEventSimpleMaker* fEventMaker = new FlowEventSimpleMaker(); 
  //------------------------------------------------------------------------
  */

  //load needed libraries
  gSystem->Load("libTree.so");

  //------------------------------------------------------------------------
  //cuts:
  AliFlowTrackSimpleCuts* cutsInt = new AliFlowTrackSimpleCuts();
  cutsInt->SetPtMax(ptMaxInt);
  cutsInt->SetPtMin(ptMinInt);
  cutsInt->SetEtaMax(etaMaxInt);
  cutsInt->SetEtaMin(etaMinInt);
  cutsInt->SetPhiMax(phiMaxInt);
  cutsInt->SetPhiMin(phiMinInt);
  cutsInt->SetPID(PIDInt);

  AliFlowTrackSimpleCuts* cutsDiff = new AliFlowTrackSimpleCuts();
  cutsDiff->SetPtMax(ptMaxDiff);
  cutsDiff->SetPtMin(ptMinDiff);
  cutsDiff->SetEtaMax(etaMaxDiff);
  cutsDiff->SetEtaMin(etaMinDiff);
  cutsDiff->SetPhiMax(phiMaxDiff);
  cutsDiff->SetPhiMin(phiMinDiff);
  cutsDiff->SetPID(PIDDiff);

  //if the weights are used: 
  TFile *fileWithWeights = NULL;
  TList *listWithWeights = NULL;
  
  if(usePhiWeights||usePtWeights||useEtaWeights)
  {
   fileWithWeights = TFile::Open("weights.root","READ");
   if(fileWithWeights)
   {
    listWithWeights = (TList*)fileWithWeights->Get("weights");
   }else{cout<<" WARNING: the file <weights.root> with weights from the previous run was not accessed."<<endl;}    
  }

  //flow methods:  
  AliFlowAnalysisWithQCumulants    *qc    = NULL;
  AliFlowAnalysisWithCumulants     *gfc   = NULL;
  AliFittingQDistribution          *fqd   = NULL;
  AliFlowAnalysisWithLeeYangZeros  *lyz1  = NULL;
  AliFlowAnalysisWithLeeYangZeros  *lyz2  = NULL;
  AliFlowAnalysisWithLYZEventPlane *lyzep = NULL;
  AliFlowAnalysisWithScalarProduct *sp    = NULL;
  AliFlowAnalysisWithMCEventPlane  *mcep  = NULL;   

  //MCEP = monte carlo event plane
  if (MCEP) {
    AliFlowAnalysisWithMCEventPlane *mcep = new AliFlowAnalysisWithMCEventPlane();
    mcep->Init();
  }

  //QC = Q-cumulants  
  if(QC) { 
    AliFlowAnalysisWithQCumulants* qc = new AliFlowAnalysisWithQCumulants();
    qc->Init();
    if(listWithWeights) qc->SetWeightsList(listWithWeights);
    if(usePhiWeights) qc->SetUsePhiWeights(usePhiWeights);
    if(usePtWeights) qc->SetUsePtWeights(usePtWeights);
    if(useEtaWeights) qc->SetUseEtaWeights(useEtaWeights);
  }
  
  //GFC = Generating Function Cumulants 
  if(GFC) {
    AliFlowAnalysisWithCumulants* gfc = new AliFlowAnalysisWithCumulants();
    gfc->CreateOutputObjects();
  }
  
  //FQD = Fitting q-distribution 
  if(FQD) {
    AliFittingQDistribution* fqd = new AliFittingQDistribution();
    fqd->CreateOutputObjects();
  }

  //SP = Scalar Product 
  if(SP) {
    AliFlowAnalysisWithScalarProduct* sp = new AliFlowAnalysisWithScalarProduct();
    sp->Init();
  }

  //LYZ1 = Lee-Yang Zeroes first run
  if(LYZ1) {
    AliFlowAnalysisWithLeeYangZeros* lyz1 = new AliFlowAnalysisWithLeeYangZeros();
    lyz1->SetFirstRun(kTRUE);
    lyz1->SetUseSum(kTRUE);
    lyz1->Init();
  }

  //LYZ2 = Lee-Yang Zeroes second run
  if(LYZ2) {
    AliFlowAnalysisWithLeeYangZeros* lyz2 = new AliFlowAnalysisWithLeeYangZeros();
    // read the input file from the first run 
    TString inputFileNameLYZ2 = "outputLYZ1analysis.root" ;
    TFile* inputFileLYZ2 = new TFile(inputFileNameLYZ2.Data(),"READ");
    if(!inputFileLYZ2 || inputFileLYZ2->IsZombie()) { 
      cerr << " ERROR: NO First Run file... " << endl ; 
    }
    else { 
      TList* inputListLYZ2 = (TList*)inputFileLYZ2->Get("cobjLYZ1");  
      if (!inputListLYZ2) {cout<<"list is NULL pointer!"<<endl;}
      else {
	cout<<"LYZ2 input file/list read..."<<endl;
	lyz2->SetFirstRunList(inputListLYZ2);
	lyz2->SetFirstRun(kFALSE);
	lyz2->SetUseSum(kTRUE);
	lyz2->Init();
      }
    }
  }
  
 //LYZEP = Lee-Yang Zeroes event plane
  if(LYZEP) {
    AliFlowLYZEventPlane* ep = new AliFlowLYZEventPlane() ;
    AliFlowAnalysisWithLYZEventPlane* lyzep = new AliFlowAnalysisWithLYZEventPlane();
    // read the input file from the second lyz run 
    TString inputFileNameLYZEP = "outputLYZ2analysis.root" ;
    TFile* inputFileLYZEP = new TFile(inputFileNameLYZEP.Data(),"READ");
    if(!inputFileLYZEP || inputFileLYZEP->IsZombie()) { 
      cerr << " ERROR: NO Second Run file... " << endl ; }
    else { 
      TList* inputListLYZEP = (TList*)inputFileLYZEP->Get("cobjLYZ2");  
      if (!inputListLYZEP) {cout<<"list is NULL pointer!"<<endl;}
      else {
	cout<<"LYZEP input file/list read..."<<endl;
	ep   ->SetSecondRunList(inputListLYZEP);
	lyzep->SetSecondRunList(inputListLYZEP);
	ep   ->Init();
	lyzep->Init();
      }
    }
  }
  //------------------------------------------------------------------------
  
  
  //standard code
  Int_t fCount = 0;
  TString execDir(gSystem->pwd());
  TString targetDir(dir);
  TSystemDirectory* baseDir = new TSystemDirectory(".", dir);  
  TList* dirList 	   = baseDir->GetListOfFiles();
  if (!dirList) {
    cout << endl << "No input files in: " << targetDir.Data() << endl;
    break;
  }
  Int_t nDirs		   = dirList->GetEntries();
  cout<<endl;
  cout<<"Int_t nDirs = "<<nDirs<<endl;
  gSystem->cd(execDir);
  
  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) 
	{
	  cout << endl; 
	  cout << "Directory (" << iDir << "): " << presentDir->GetName() << 
	    " - Skipping ... " << endl;
	  continue ;   
	}
      
      if(offset > 0) { --offset ; continue ; }
      if((aRuns > 0) && (fCount >= aRuns)) { break ; }
      
      TString presentDirName(dir); // aDataDir
      presentDirName += presentDir->GetName();
      presentDirName += "/";
      //cerr<<" presentDirName = "<<presentDirName<<endl;
      
      TString fileName = presentDirName; 
      fileName += "galice.root";
      Long_t *id, *size, *flags, *modtime;
      if(gSystem->GetPathInfo(fileName.Data(),id,size,flags,modtime)) 
	{ 
	  cout << " File : " << fileName << " does NOT exist ! - Skipping ... " 
	       << endl; 
	  continue; 
	}
      //      cout << endl ; cout << "Directory (" << iDir << "): " << presentDirName << " ... " << endl;
      
      //loop (simulations in the present dir) 
      TSystemDirectory* evtsDir = new TSystemDirectory(".", 
						       presentDirName.Data());
      TList* fileList 	    = evtsDir->GetListOfFiles();
      Int_t nFiles		    = fileList->GetEntries();
      //cout<<" Int_t nFiles = "<<nFiles<<endl;
      gSystem->cd(execDir);      
      for(Int_t iFiles=0; iFiles<nFiles; ++iFiles)
	{
	  TSystemFile* presentFile = (TSystemFile*) fileList->At(iFiles);
	  TString presentFileName(presentDirName);
	  presentFileName += presentFile->GetName();
	  
	  if(!(presentFileName.Contains("Kinematics") && 
	       presentFileName.Contains("root"))) { continue ; }
	  
	  //cout << " found: " << presentFileName.Data() << endl; 
	  
	  TFile* kineFile = new TFile(presentFileName.Data(), "READ"); 
	  // kineFile->ls();
	  Int_t nEvts = kineFile->GetNkeys() ; 
	  //cout << "  . found: " << nEvts << " KineTree(s) in " << presentFileName.Data() << endl;
	  TList* kineEventsList = (TList*)kineFile->GetListOfKeys(); 
	  TTree* kTree;
	  TIter next(kineEventsList); 
	  TKey* key;
	  
	  //loop over the events
	  while( key=(TKey *)next() ) 
	    {
	      TDirectory* tDir = (TDirectory*)key->ReadObj();
	      if(!tDir) break;
	      
	      TString evtDir(tDir->GetName()); 
	      //cout << "  . . found: " << tDir->GetName() << endl;
	      
	      kTree = (TTree *)tDir->Get("TreeK");
	      if(!kTree) break;
	      
	      Int_t nPart = kTree->GetEntries();
	      //cout << "  . . . kTree " << fCount << " has " << nPart << " particles " << endl; 
	      
	      //-----------------------------------------------------------
	      //fill and save the flow event	      
	      AliFlowEventSimple *fEvent = fEventMaker->FillTracks(kTree, cutsInt, cutsDiff); 
	                    
	      //pass the flow event to flow methods for analysis 
	      Double_t fEP = 0.; // temporary number need true value
	      if(MCEP) mcep->Make(fEvent,fEP);  //fix fEP
	      if(QC) qc->Make(fEvent);
	      if(GFC) gfc->Make(fEvent);
	      if(FQD) fqd->Make(fEvent);
	      if(LYZ1) lyz1->Make(fEvent);
	      if(LYZ2) lyz2->Make(fEvent);
	      if(LYZEP) lyzep->Make(fEvent,ep);
	      if(SP) sp->Make(fEvent);
	      //-----------------------------------------------------------
	      
	      fCount++;
	      cout << "# " << fCount << " events processed" << endl;
	      delete kTree;
	      delete fEvent;
	    }
	  delete kineFile ;
	}
      delete evtsDir ;
    }

  //--------------------------------------------------------------
  //calculating and storing the final results of flow analysis
  if(MCEP) {mcep->Finish(); mcep->WriteHistograms("outputMCEPanalysis.root");}
  if(QC) {qc->Finish(); qc->WriteHistograms("outputQCanalysis.root");}
  if(GFC) {gfc->Finish(); gfc->WriteHistograms("outputGFCanalysis.root");}
  if(FQD) {fqd->Finish(); fqd->WriteHistograms("outputFQDanalysis.root");}
  if(LYZ1) {lyz1->Finish(); lyz1->WriteHistograms("outputLYZ1analysis.root");}
  if(LYZ2) {lyz2->Finish(); lyz2->WriteHistograms("outputLYZ2analysis.root");}
  if(LYZEP) {lyzep->Finish(); lyzep->WriteHistograms("outputLYZEPanalysis.root");}
  if(SP) {sp->Finish(); sp->WriteHistograms("outputSPanalysis.root");}
  //--------------------------------------------------------------
  
  cout << endl;
  cout << " Finished ... " << endl;
  cout << endl;
  
  timer.Stop();
  cout << endl;
  timer.Print();
}
Commit	Line	Data
334e3256	1	#include "TStopwatch.h"
	2	#include "TObjArray"
	3	#include "Riostream.h"
26c4cbb9	4	#include "TFile.h"
334e3256	5
03a02aca	6	//--------------------------------------------------------------------------------------
03a02aca	7	// RUN SETTINGS
334e3256	8	//flow analysis method can be: (set to kTRUE or kFALSE)
29b61d43	9	Bool_t SP = kTRUE;
29b61d43	10	Bool_t LYZ1 = kTRUE;
334e3256	11	Bool_t LYZ2 = kFALSE;
	12	Bool_t LYZEP = kFALSE;
	13	Bool_t GFC = kTRUE;
	14	Bool_t QC = kTRUE;
29b61d43	15	Bool_t FQD = kTRUE;
334e3256	16	Bool_t MCEP = kFALSE; //does not work yet 24/12/08
03a02aca	17	//--------------------------------------------------------------------------------------
334e3256	18
03a02aca	19	//--------------------------------------------------------------------------------------
03a02aca	20	// CUT SETTINGS
334e3256	21	//integrated selection
	22	Double_t ptMaxInt = 10.;
	23	Double_t ptMinInt = 0.;
	24	Double_t etaMaxInt = 1.;
	25	Double_t etaMinInt = -1.;
	26	Double_t phiMaxInt = 7.5;
	27	Double_t phiMinInt = 0.;
	28	Int_t PIDInt = 211;
	29
	30	//differential selection
	31	Double_t ptMaxDiff = 10.;
	32	Double_t ptMinDiff = 0.;
	33	Double_t etaMaxDiff = 1.;
	34	Double_t etaMinDiff = -1.;
	35	Double_t phiMaxDiff = 7.5;
	36	Double_t phiMinDiff = 0.;
	37	Int_t PIDDiff = 211;
03a02aca	38	//--------------------------------------------------------------------------------------
26c4cbb9	39
03a02aca	40	//--------------------------------------------------------------------------------------
03a02aca	41	// WEIGHTS SETTINGS
26c4cbb9	42	//to use or not to use the weights - that is a question!
03a02aca	43	Bool_t usePhiWeights = kFALSE; //Phi (correction for non-uniform azimuthal acceptance)
	44	Bool_t usePtWeights = kFALSE; //v'(pt) (differential flow in pt)
	45	Bool_t useEtaWeights = kFALSE; //v'(eta) (differential flow in eta)
	46	//--------------------------------------------------------------------------------------
334e3256	47
26c4cbb9	48	Int_t offset = 0;
26c4cbb9	49
29b61d43	50	int runFlowAnalysis(Int_t aRuns = 100, const char*
	51	// dir="/data/alice1/kolk/KineOnly3/")
	52	dir="/Users/snelling/alice_data/KineOnly3/")
334e3256	53	{
	54	TStopwatch timer;
	55	timer.Start();
	56
	57	if (LYZ1 && LYZ2) {cout<<"WARNING: you cannot run LYZ1 and LYZ2 at the same time! LYZ2 needs the output from LYZ1."<<endl; exit(); }
334e3256	58	if (LYZ2 && LYZEP) {cout<<"WARNING: you cannot run LYZ2 and LYZEP at the same time! LYZEP needs the output from LYZ2."<<endl; exit(); }
334e3256	59	if (LYZ1 && LYZEP) {cout<<"WARNING: you cannot run LYZ1 and LYZEP at the same time! LYZEP needs the output from LYZ2."<<endl; exit(); }
	60
	61
	62	cout<<endl;
	63	cout<<" ---- BEGIN ANALYSIS ---- "<<endl;
	64	cout<<endl;
	65
2a413443	66
334e3256	67	gSystem->AddIncludePath("-I$ROOTSYS/include");
2a413443	68
	69	// In AliRoot
	70	gSystem->AddIncludePath("-I$ALICE_ROOT/include");
	71
	72	gSystem->Load("libANALYSIS.so");
29b61d43	73	gSystem->Load("libPWG2flowCommon.so");
	74	cerr<<"libPWG2flowCommon.so loaded ..."<<endl;
	75	gSystem->Load("libPWG2flowTasks.so");
	76	cerr<<"libPWG2flowTasks.so loaded ..."<<endl;
2a413443	77	cout<<endl;
26c4cbb9	78
29b61d43	79	// Flow event in AliRoot
26c4cbb9	80	AliFlowEventSimpleMaker* fEventMaker = new AliFlowEventSimpleMaker();
03a02aca	81
2a413443	82	// In root
	83
	84	/*
29b61d43	85	// Constants
	86	gROOT->LoadMacro("AliFlowCommon/AliFlowCommonConstants.cxx+");
	87	gROOT->LoadMacro("AliFlowCommon/AliFlowLYZConstants.cxx+");
	88	gROOT->LoadMacro("AliFlowCommon/AliFlowCumuConstants.cxx+");
	89
	90	// Flow event
	91	gROOT->LoadMacro("AliFlowCommon/AliFlowVector.cxx+");
	92	gROOT->LoadMacro("AliFlowCommon/AliFlowTrackSimple.cxx+");
	93	gROOT->LoadMacro("AliFlowCommon/AliFlowEventSimple.cxx+");
	94
	95	// Cuts
	96	gROOT->LoadMacro("AliFlowCommon/AliFlowTrackSimpleCuts.cxx+");
	97
	98	// Output histosgrams
	99	gROOT->LoadMacro("AliFlowCommon/AliFlowCommonHist.cxx+");
	100	gROOT->LoadMacro("AliFlowCommon/AliFlowCommonHistResults.cxx+");
	101	gROOT->LoadMacro("AliFlowCommon/AliFlowLYZHist1.cxx+");
	102	gROOT->LoadMacro("AliFlowCommon/AliFlowLYZHist2.cxx+");
	103
	104	// Functions needed for various methods
	105	gROOT->LoadMacro("AliFlowCommon/AliCumulantsFunctions.cxx+");
	106	gROOT->LoadMacro("AliFlowCommon/AliQCumulantsFunctions.cxx+");
	107	gROOT->LoadMacro("AliFlowCommon/AliFittingFunctionsForQDistribution.cxx+");
	108	gROOT->LoadMacro("AliFlowCommon/AliFlowLYZEventPlane.cxx+");
334e3256	109
334e3256	110	// Flow Analysis code for various methods
29b61d43	111	gROOT->LoadMacro("AliFlowCommon/AliFlowAnalysisWithMCEventPlane.cxx+");
	112	gROOT->LoadMacro("AliFlowCommon/AliFlowAnalysisWithScalarProduct.cxx+");
	113	gROOT->LoadMacro("AliFlowCommon/AliFlowAnalysisWithLYZEventPlane.cxx+");
	114	gROOT->LoadMacro("AliFlowCommon/AliFlowAnalysisWithLeeYangZeros.cxx+");
	115	gROOT->LoadMacro("AliFlowCommon/AliFlowAnalysisWithCumulants.cxx+");
	116	gROOT->LoadMacro("AliFlowCommon/AliFlowAnalysisWithQCumulants.cxx+");
	117	gROOT->LoadMacro("AliFlowCommon/AliFittingQDistribution.cxx+");
334e3256	118
2a413443	119	// Class to fill the FlowEvent without aliroot dependence
29b61d43	120	// can be found in the directory FlowEventMakers
29b61d43	121	gROOT->LoadMacro("FlowEventMakers/FlowEventSimpleMaker.cxx+");
334e3256	122
2a413443	123	cout << "finished loading macros!" << endl;
	124
	125	//flow event
	126	FlowEventSimpleMaker* fEventMaker = new FlowEventSimpleMaker();
	127	//------------------------------------------------------------------------
	128	*/
	129
334e3256	130	//load needed libraries
	131	gSystem->Load("libTree.so");
	132
334e3256	133	//------------------------------------------------------------------------
e085f1a9	134	//cuts:
334e3256	135	AliFlowTrackSimpleCuts* cutsInt = new AliFlowTrackSimpleCuts();
	136	cutsInt->SetPtMax(ptMaxInt);
	137	cutsInt->SetPtMin(ptMinInt);
	138	cutsInt->SetEtaMax(etaMaxInt);
	139	cutsInt->SetEtaMin(etaMinInt);
	140	cutsInt->SetPhiMax(phiMaxInt);
	141	cutsInt->SetPhiMin(phiMinInt);
	142	cutsInt->SetPID(PIDInt);
	143
	144	AliFlowTrackSimpleCuts* cutsDiff = new AliFlowTrackSimpleCuts();
	145	cutsDiff->SetPtMax(ptMaxDiff);
	146	cutsDiff->SetPtMin(ptMinDiff);
	147	cutsDiff->SetEtaMax(etaMaxDiff);
	148	cutsDiff->SetEtaMin(etaMinDiff);
	149	cutsDiff->SetPhiMax(phiMaxDiff);
	150	cutsDiff->SetPhiMin(phiMinDiff);
	151	cutsDiff->SetPID(PIDDiff);
	152
03a02aca	153	//if the weights are used:
	154	TFile *fileWithWeights = NULL;
	155	TList *listWithWeights = NULL;
334e3256	156
03a02aca	157	if(usePhiWeights\|\|usePtWeights\|\|useEtaWeights)
e085f1a9	158	{
03a02aca	159	fileWithWeights = TFile::Open("weights.root","READ");
	160	if(fileWithWeights)
	161	{
	162	listWithWeights = (TList*)fileWithWeights->Get("weights");
	163	}else{cout<<" WARNING: the file <weights.root> with weights from the previous run was not accessed."<<endl;}
e085f1a9	164	}
	165
	166	//flow methods:
334e3256	167	AliFlowAnalysisWithQCumulants *qc = NULL;
	168	AliFlowAnalysisWithCumulants *gfc = NULL;
	169	AliFittingQDistribution *fqd = NULL;
	170	AliFlowAnalysisWithLeeYangZeros *lyz1 = NULL;
	171	AliFlowAnalysisWithLeeYangZeros *lyz2 = NULL;
	172	AliFlowAnalysisWithLYZEventPlane *lyzep = NULL;
	173	AliFlowAnalysisWithScalarProduct *sp = NULL;
e085f1a9	174	AliFlowAnalysisWithMCEventPlane *mcep = NULL;
deebd72f	175
334e3256	176	//MCEP = monte carlo event plane
deebd72f	177	if (MCEP) {
	178	AliFlowAnalysisWithMCEventPlane *mcep = new AliFlowAnalysisWithMCEventPlane();
	179	mcep->Init();
	180	}
334e3256	181
334e3256	182	//QC = Q-cumulants
deebd72f	183	if(QC) {
deebd72f	184	AliFlowAnalysisWithQCumulants* qc = new AliFlowAnalysisWithQCumulants();
e085f1a9	185	qc->Init();
03a02aca	186	if(listWithWeights) qc->SetWeightsList(listWithWeights);
	187	if(usePhiWeights) qc->SetUsePhiWeights(usePhiWeights);
	188	if(usePtWeights) qc->SetUsePtWeights(usePtWeights);
	189	if(useEtaWeights) qc->SetUseEtaWeights(useEtaWeights);
deebd72f	190	}
334e3256	191
334e3256	192	//GFC = Generating Function Cumulants
deebd72f	193	if(GFC) {
	194	AliFlowAnalysisWithCumulants* gfc = new AliFlowAnalysisWithCumulants();
	195	gfc->CreateOutputObjects();
	196	}
334e3256	197
334e3256	198	//FQD = Fitting q-distribution
deebd72f	199	if(FQD) {
	200	AliFittingQDistribution* fqd = new AliFittingQDistribution();
	201	fqd->CreateOutputObjects();
	202	}
60ec66fb	203
	204	//SP = Scalar Product
	205	if(SP) {
	206	AliFlowAnalysisWithScalarProduct* sp = new AliFlowAnalysisWithScalarProduct();
	207	sp->Init();
	208	}
	209
334e3256	210	//LYZ1 = Lee-Yang Zeroes first run
deebd72f	211	if(LYZ1) {
	212	AliFlowAnalysisWithLeeYangZeros* lyz1 = new AliFlowAnalysisWithLeeYangZeros();
	213	lyz1->SetFirstRun(kTRUE);
	214	lyz1->SetUseSum(kTRUE);
	215	lyz1->Init();
	216	}
334e3256	217
334e3256	218	//LYZ2 = Lee-Yang Zeroes second run
deebd72f	219	if(LYZ2) {
	220	AliFlowAnalysisWithLeeYangZeros* lyz2 = new AliFlowAnalysisWithLeeYangZeros();
	221	// read the input file from the first run
	222	TString inputFileNameLYZ2 = "outputLYZ1analysis.root" ;
	223	TFile* inputFileLYZ2 = new TFile(inputFileNameLYZ2.Data(),"READ");
	224	if(!inputFileLYZ2 \|\| inputFileLYZ2->IsZombie()) {
	225	cerr << " ERROR: NO First Run file... " << endl ;
	226	}
	227	else {
	228	TList* inputListLYZ2 = (TList*)inputFileLYZ2->Get("cobjLYZ1");
	229	if (!inputListLYZ2) {cout<<"list is NULL pointer!"<<endl;}
	230	else {
	231	cout<<"LYZ2 input file/list read..."<<endl;
	232	lyz2->SetFirstRunList(inputListLYZ2);
	233	lyz2->SetFirstRun(kFALSE);
	234	lyz2->SetUseSum(kTRUE);
	235	lyz2->Init();
334e3256	236	}
334e3256	237	}
deebd72f	238	}
deebd72f	239
334e3256	240	//LYZEP = Lee-Yang Zeroes event plane
deebd72f	241	if(LYZEP) {
	242	AliFlowLYZEventPlane* ep = new AliFlowLYZEventPlane() ;
	243	AliFlowAnalysisWithLYZEventPlane* lyzep = new AliFlowAnalysisWithLYZEventPlane();
	244	// read the input file from the second lyz run
	245	TString inputFileNameLYZEP = "outputLYZ2analysis.root" ;
	246	TFile* inputFileLYZEP = new TFile(inputFileNameLYZEP.Data(),"READ");
	247	if(!inputFileLYZEP \|\| inputFileLYZEP->IsZombie()) {
	248	cerr << " ERROR: NO Second Run file... " << endl ; }
	249	else {
	250	TList* inputListLYZEP = (TList*)inputFileLYZEP->Get("cobjLYZ2");
	251	if (!inputListLYZEP) {cout<<"list is NULL pointer!"<<endl;}
	252	else {
	253	cout<<"LYZEP input file/list read..."<<endl;
	254	ep ->SetSecondRunList(inputListLYZEP);
	255	lyzep->SetSecondRunList(inputListLYZEP);
	256	ep ->Init();
	257	lyzep->Init();
334e3256	258	}
334e3256	259	}
deebd72f	260	}
334e3256	261	//------------------------------------------------------------------------
334e3256	262
26c4cbb9	263
334e3256	264	//standard code
	265	Int_t fCount = 0;
	266	TString execDir(gSystem->pwd());
	267	TString targetDir(dir);
	268	TSystemDirectory* baseDir = new TSystemDirectory(".", dir);
	269	TList* dirList = baseDir->GetListOfFiles();
	270	if (!dirList) {
	271	cout << endl << "No input files in: " << targetDir.Data() << endl;
	272	break;
	273	}
	274	Int_t nDirs = dirList->GetEntries();
	275	cout<<endl;
	276	cout<<"Int_t nDirs = "<<nDirs<<endl;
	277	gSystem->cd(execDir);
	278
	279	for(Int_t iDir=0;iDir<nDirs;++iDir)
	280	{
	281	TSystemFile* presentDir = (TSystemFile*)dirList->At(iDir);
	282	if(!presentDir \|\| !presentDir->IsDirectory() \|\|
	283	strcmp(presentDir->GetName(), ".") == 0 \|\|
	284	strcmp(presentDir->GetName(), "..") == 0)
	285	{
	286	cout << endl;
	287	cout << "Directory (" << iDir << "): " << presentDir->GetName() <<
	288	" - Skipping ... " << endl;
	289	continue ;
	290	}
	291
	292	if(offset > 0) { --offset ; continue ; }
	293	if((aRuns > 0) && (fCount >= aRuns)) { break ; }
	294
	295	TString presentDirName(dir); // aDataDir
	296	presentDirName += presentDir->GetName();
	297	presentDirName += "/";
	298	//cerr<<" presentDirName = "<<presentDirName<<endl;
	299
	300	TString fileName = presentDirName;
	301	fileName += "galice.root";
	302	Long_t id, size, flags, modtime;
	303	if(gSystem->GetPathInfo(fileName.Data(),id,size,flags,modtime))
	304	{
	305	cout << " File : " << fileName << " does NOT exist ! - Skipping ... "
	306	<< endl;
	307	continue;
	308	}
deebd72f	309	// cout << endl ; cout << "Directory (" << iDir << "): " << presentDirName << " ... " << endl;
334e3256	310
	311	//loop (simulations in the present dir)
	312	TSystemDirectory* evtsDir = new TSystemDirectory(".",
	313	presentDirName.Data());
	314	TList* fileList = evtsDir->GetListOfFiles();
	315	Int_t nFiles = fileList->GetEntries();
deebd72f	316	//cout<<" Int_t nFiles = "<<nFiles<<endl;
334e3256	317	gSystem->cd(execDir);
	318	for(Int_t iFiles=0; iFiles<nFiles; ++iFiles)
	319	{
	320	TSystemFile* presentFile = (TSystemFile*) fileList->At(iFiles);
	321	TString presentFileName(presentDirName);
	322	presentFileName += presentFile->GetName();
	323
	324	if(!(presentFileName.Contains("Kinematics") &&
	325	presentFileName.Contains("root"))) { continue ; }
	326
deebd72f	327	//cout << " found: " << presentFileName.Data() << endl;
334e3256	328
	329	TFile* kineFile = new TFile(presentFileName.Data(), "READ");
	330	// kineFile->ls();
	331	Int_t nEvts = kineFile->GetNkeys() ;
deebd72f	332	//cout << " . found: " << nEvts << " KineTree(s) in " << presentFileName.Data() << endl;
334e3256	333	TList* kineEventsList = (TList*)kineFile->GetListOfKeys();
	334	TTree* kTree;
	335	TIter next(kineEventsList);
	336	TKey* key;
	337
	338	//loop over the events
	339	while( key=(TKey *)next() )
	340	{
	341	TDirectory* tDir = (TDirectory*)key->ReadObj();
	342	if(!tDir) break;
	343
	344	TString evtDir(tDir->GetName());
deebd72f	345	//cout << " . . found: " << tDir->GetName() << endl;
334e3256	346
	347	kTree = (TTree *)tDir->Get("TreeK");
	348	if(!kTree) break;
	349
	350	Int_t nPart = kTree->GetEntries();
deebd72f	351	//cout << " . . . kTree " << fCount << " has " << nPart << " particles " << endl;
334e3256	352
334e3256	353	//-----------------------------------------------------------
26c4cbb9	354	//fill and save the flow event
	355	AliFlowEventSimple *fEvent = fEventMaker->FillTracks(kTree, cutsInt, cutsDiff);
	356
334e3256	357	//pass the flow event to flow methods for analysis
60ec66fb	358	Double_t fEP = 0.; // temporary number need true value
60ec66fb	359	if(MCEP) mcep->Make(fEvent,fEP); //fix fEP
deebd72f	360	if(QC) qc->Make(fEvent);
	361	if(GFC) gfc->Make(fEvent);
	362	if(FQD) fqd->Make(fEvent);
	363	if(LYZ1) lyz1->Make(fEvent);
	364	if(LYZ2) lyz2->Make(fEvent);
	365	if(LYZEP) lyzep->Make(fEvent,ep);
	366	if(SP) sp->Make(fEvent);
334e3256	367	//-----------------------------------------------------------
	368
	369	fCount++;
deebd72f	370	cout << "# " << fCount << " events processed" << endl;
334e3256	371	delete kTree;
	372	delete fEvent;
	373	}
	374	delete kineFile ;
	375	}
	376	delete evtsDir ;
	377	}
	378
	379	//--------------------------------------------------------------
	380	//calculating and storing the final results of flow analysis
60ec66fb	381	if(MCEP) {mcep->Finish(); mcep->WriteHistograms("outputMCEPanalysis.root");}
deebd72f	382	if(QC) {qc->Finish(); qc->WriteHistograms("outputQCanalysis.root");}
	383	if(GFC) {gfc->Finish(); gfc->WriteHistograms("outputGFCanalysis.root");}
	384	if(FQD) {fqd->Finish(); fqd->WriteHistograms("outputFQDanalysis.root");}
	385	if(LYZ1) {lyz1->Finish(); lyz1->WriteHistograms("outputLYZ1analysis.root");}
	386	if(LYZ2) {lyz2->Finish(); lyz2->WriteHistograms("outputLYZ2analysis.root");}
	387	if(LYZEP) {lyzep->Finish(); lyzep->WriteHistograms("outputLYZEPanalysis.root");}
	388	if(SP) {sp->Finish(); sp->WriteHistograms("outputSPanalysis.root");}
334e3256	389	//--------------------------------------------------------------
	390
	391	cout << endl;
	392	cout << " Finished ... " << endl;
	393	cout << endl;
	394
	395	timer.Stop();
	396	cout << endl;
	397	timer.Print();
	398	}