[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
//--------------------------------------------------------------------------------------

// Weights 
// Use weights for Q vector
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)

//--------------------------------------------------------------------------------------
// 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;
//--------------------------------------------------------------------------------------

enum anaModes {mLocal,mLocalSource,mLocalPAR,};
//mLocal: Analyze data on your computer using aliroot
//mLocalPAR: Analyze data on your computer using root + PAR files
//mLocalSource: Analyze data on your computer using root + source files


Int_t offset = 0;
                                          
int runFlowAnalysis(Int_t mode=mLocal, Int_t aRuns = 1000, 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;
  
  LoadLibraries(mode);

  if (mode == mLocal || mode == mLocalPAR) {
    // AliFlow event in aliroot or with pars
    AliFlowEventSimpleMaker* fEventMaker = new AliFlowEventSimpleMaker();
  }
  else if (mode == mLocalSource) {
    // flow event in source mode
    FlowEventSimpleMaker* fEventMaker = new FlowEventSimpleMaker(); 
  }
  else{
    cout << "No supported running mode selected!" << endl;
    break;
  }


  //------------------------------------------------------------------------
  //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 found."<<endl;
	break;
      }    
  }

  //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->Init();
    if(listWithWeights) gfc->SetWeightsList(listWithWeights);
    if(usePhiWeights) gfc->SetUsePhiWeights(usePhiWeights);
    if(usePtWeights) gfc->SetUsePtWeights(usePtWeights);
    if(useEtaWeights) gfc->SetUseEtaWeights(useEtaWeights);
  }
  
  //FQD = Fitting q-distribution 
  if(FQD) {
    AliFittingQDistribution* fqd = new AliFittingQDistribution();
    fqd->Init();
    if(listWithWeights) fqd->SetWeightsList(listWithWeights);
    if(usePhiWeights) fqd->SetUsePhiWeights(usePhiWeights);
  }

  //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 ;
      break; 
    }
    else { 
      TList* inputListLYZ2 = (TList*)inputFileLYZ2->Get("cobjLYZ1");  
      if (!inputListLYZ2) {cout<<"Input list is NULL pointer!"<<endl; break;}
      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 ; 
      break;
    }
    else { 
      TList* inputListLYZEP = (TList*)inputFileLYZEP->Get("cobjLYZ2");  
      if (!inputListLYZEP) {cout<<"Input list is NULL pointer!"<<endl; break;}
      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); 
	                    
	      if(MCEP) mcep->Make(fEvent);
	      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(SP) {sp->Finish(); sp->WriteHistograms("outputSPanalysis.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");}

  //--------------------------------------------------------------
  
  cout << endl;
  cout << " Fini ... " << endl;
  cout << endl;
  
  timer.Stop();
  cout << endl;
  timer.Print();
}

void SetupPar(char* pararchivename)
{
  //Load par files, create analysis libraries
  //For testing, if par file already decompressed and modified
  //classes then do not decompress.
 
  TString cdir(Form("%s", gSystem->WorkingDirectory() )) ; 
  TString parpar(Form("%s.par", pararchivename)) ; 
  if ( gSystem->AccessPathName(parpar.Data()) ) {
    gSystem->ChangeDirectory(gSystem->Getenv("ALICE_ROOT")) ;
    TString processline(Form(".! make %s", parpar.Data())) ; 
    gROOT->ProcessLine(processline.Data()) ;
    gSystem->ChangeDirectory(cdir) ; 
    processline = Form(".! mv /tmp/%s .", parpar.Data()) ;
    gROOT->ProcessLine(processline.Data()) ;
  } 
  if ( gSystem->AccessPathName(pararchivename) ) {  
    TString processline = Form(".! tar xvzf %s",parpar.Data()) ;
    gROOT->ProcessLine(processline.Data());
  }
  
  TString ocwd = gSystem->WorkingDirectory();
  gSystem->ChangeDirectory(pararchivename);
  
  // check for BUILD.sh and execute
  if (!gSystem->AccessPathName("PROOF-INF/BUILD.sh")) {
    printf("*******************************\n");
    printf("*** Building PAR archive    ***\n");
    cout<<pararchivename<<endl;
    printf("*******************************\n");
    
    if (gSystem->Exec("PROOF-INF/BUILD.sh")) {
      Error("runProcess","Cannot Build the PAR Archive! - Abort!");
      return -1;
    }
  }
  // check for SETUP.C and execute
  if (!gSystem->AccessPathName("PROOF-INF/SETUP.C")) {
    printf("*******************************\n");
    printf("*** Setup PAR archive       ***\n");
    cout<<pararchivename<<endl;
    printf("*******************************\n");
    gROOT->Macro("PROOF-INF/SETUP.C");
  }
  
  gSystem->ChangeDirectory(ocwd.Data());
  printf("Current dir: %s\n", ocwd.Data());
}

void LoadLibraries(const anaModes mode) {
  
  //--------------------------------------
  // Load the needed libraries most of them already loaded by aliroot
  //--------------------------------------
  gSystem->Load("libTree.so");
  gSystem->Load("libGeom.so");
  gSystem->Load("libVMC.so");
  gSystem->Load("libXMLIO.so");
  gSystem->Load("libPhysics.so");
  
  //----------------------------------------------------------
  // >>>>>>>>>>> Local mode <<<<<<<<<<<<<< 
  //----------------------------------------------------------
  if (mode==mLocal) {
    //--------------------------------------------------------
    // If you want to use already compiled libraries 
    // in the aliroot distribution
    //--------------------------------------------------------
    gSystem->Load("libSTEERBase");
    gSystem->Load("libESD");
    gSystem->Load("libAOD");
    gSystem->Load("libANALYSIS");
    gSystem->Load("libANALYSISalice");
    gSystem->Load("libCORRFW.so");
    cerr<<"libCORRFW.so loaded..."<<endl;
    gSystem->Load("libPWG2flowCommon.so");
    cerr<<"libPWG2flowCommon.so loaded..."<<endl;
    gSystem->Load("libPWG2flowTasks.so");
    cerr<<"libPWG2flowTasks.so loaded..."<<endl;
  }
  
  else if (mode == mLocalPAR) {
    //--------------------------------------------------------
    //If you want to use root and par files from aliroot
    //--------------------------------------------------------  
     //If you want to use root and par files from aliroot
    //--------------------------------------------------------  
    SetupPar("STEERBase");
    SetupPar("ESD");
    SetupPar("AOD");
    SetupPar("ANALYSIS");
    SetupPar("ANALYSISalice");
    SetupPar("PWG2AOD");
    SetupPar("CORRFW");
    SetupPar("PWG2flowCommon");
    cerr<<"PWG2flowCommon.par loaded..."<<endl;
    SetupPar("PWG2flowTasks");
    cerr<<"PWG2flowTasks.par loaded..."<<endl;
  }
  
  //---------------------------------------------------------
  // <<<<<<<<<< Source mode >>>>>>>>>>>>
  //---------------------------------------------------------
  else if (mode==mLocalSource) {
 
    // In root inline compile

   
    // 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;  
    
  }  
  
}
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;
ce78b98c	15	Bool_t FQD = kTRUE;
334e3256	16	Bool_t MCEP = kFALSE; //does not work yet 24/12/08
03a02aca	17	//--------------------------------------------------------------------------------------
334e3256	18
20e52807	19	// Weights
	20	// Use weights for Q vector
	21	Bool_t usePhiWeights = kFALSE; //Phi (correction for non-uniform azimuthal acceptance)
	22	Bool_t usePtWeights = kFALSE; //v'(pt) (differential flow in pt)
	23	Bool_t useEtaWeights = kFALSE; //v'(eta) (differential flow in eta)
	24
03a02aca	25	//--------------------------------------------------------------------------------------
03a02aca	26	// CUT SETTINGS
334e3256	27	//integrated selection
	28	Double_t ptMaxInt = 10.;
	29	Double_t ptMinInt = 0.;
	30	Double_t etaMaxInt = 1.;
	31	Double_t etaMinInt = -1.;
	32	Double_t phiMaxInt = 7.5;
	33	Double_t phiMinInt = 0.;
	34	Int_t PIDInt = 211;
	35
	36	//differential selection
	37	Double_t ptMaxDiff = 10.;
	38	Double_t ptMinDiff = 0.;
	39	Double_t etaMaxDiff = 1.;
	40	Double_t etaMinDiff = -1.;
	41	Double_t phiMaxDiff = 7.5;
	42	Double_t phiMinDiff = 0.;
	43	Int_t PIDDiff = 211;
03a02aca	44	//--------------------------------------------------------------------------------------
26c4cbb9	45
fdd8c18f	46	enum anaModes {mLocal,mLocalSource,mLocalPAR,};
	47	//mLocal: Analyze data on your computer using aliroot
	48	//mLocalPAR: Analyze data on your computer using root + PAR files
	49	//mLocalSource: Analyze data on your computer using root + source files
	50
	51
334e3256	52
26c4cbb9	53	Int_t offset = 0;
26c4cbb9	54
d7671632	55	int runFlowAnalysis(Int_t mode=mLocal, Int_t aRuns = 1000, const char*
ce78b98c	56	// dir="/data/alice1/kolk/KineOnly3/")
ce78b98c	57	dir="/Users/snelling/alice_data/KineOnly3/")
334e3256	58	{
	59	TStopwatch timer;
	60	timer.Start();
	61
	62	if (LYZ1 && LYZ2) {cout<<"WARNING: you cannot run LYZ1 and LYZ2 at the same time! LYZ2 needs the output from LYZ1."<<endl; exit(); }
334e3256	63	if (LYZ2 && LYZEP) {cout<<"WARNING: you cannot run LYZ2 and LYZEP at the same time! LYZEP needs the output from LYZ2."<<endl; exit(); }
334e3256	64	if (LYZ1 && LYZEP) {cout<<"WARNING: you cannot run LYZ1 and LYZEP at the same time! LYZEP needs the output from LYZ2."<<endl; exit(); }
	65
	66
	67	cout<<endl;
	68	cout<<" ---- BEGIN ANALYSIS ---- "<<endl;
	69	cout<<endl;
	70
fdd8c18f	71	LoadLibraries(mode);
2a413443	72
fdd8c18f	73	if (mode == mLocal \|\| mode == mLocalPAR) {
	74	// AliFlow event in aliroot or with pars
	75	AliFlowEventSimpleMaker* fEventMaker = new AliFlowEventSimpleMaker();
	76	}
	77	else if (mode == mLocalSource) {
	78	// flow event in source mode
	79	FlowEventSimpleMaker* fEventMaker = new FlowEventSimpleMaker();
	80	}
	81	else{
	82	cout << "No supported running mode selected!" << endl;
	83	break;
	84	}
2a413443	85
2a413443	86
334e3256	87	//------------------------------------------------------------------------
e085f1a9	88	//cuts:
334e3256	89	AliFlowTrackSimpleCuts* cutsInt = new AliFlowTrackSimpleCuts();
	90	cutsInt->SetPtMax(ptMaxInt);
	91	cutsInt->SetPtMin(ptMinInt);
	92	cutsInt->SetEtaMax(etaMaxInt);
	93	cutsInt->SetEtaMin(etaMinInt);
	94	cutsInt->SetPhiMax(phiMaxInt);
	95	cutsInt->SetPhiMin(phiMinInt);
	96	cutsInt->SetPID(PIDInt);
	97
	98	AliFlowTrackSimpleCuts* cutsDiff = new AliFlowTrackSimpleCuts();
	99	cutsDiff->SetPtMax(ptMaxDiff);
	100	cutsDiff->SetPtMin(ptMinDiff);
	101	cutsDiff->SetEtaMax(etaMaxDiff);
	102	cutsDiff->SetEtaMin(etaMinDiff);
	103	cutsDiff->SetPhiMax(phiMaxDiff);
	104	cutsDiff->SetPhiMin(phiMinDiff);
	105	cutsDiff->SetPID(PIDDiff);
	106
03a02aca	107	//if the weights are used:
	108	TFile *fileWithWeights = NULL;
	109	TList *listWithWeights = NULL;
334e3256	110
14db9c04	111	if(usePhiWeights\|\|usePtWeights\|\|useEtaWeights) {
	112	fileWithWeights = TFile::Open("weights.root","READ");
	113	if(fileWithWeights) {
	114	listWithWeights = (TList*)fileWithWeights->Get("weights");
	115	}
	116	else
	117	{cout << " WARNING: the file <weights.root> with weights from the previous run was not found."<<endl;
	118	break;
	119	}
e085f1a9	120	}
	121
	122	//flow methods:
334e3256	123	AliFlowAnalysisWithQCumulants *qc = NULL;
	124	AliFlowAnalysisWithCumulants *gfc = NULL;
	125	AliFittingQDistribution *fqd = NULL;
	126	AliFlowAnalysisWithLeeYangZeros *lyz1 = NULL;
	127	AliFlowAnalysisWithLeeYangZeros *lyz2 = NULL;
	128	AliFlowAnalysisWithLYZEventPlane *lyzep = NULL;
	129	AliFlowAnalysisWithScalarProduct *sp = NULL;
e085f1a9	130	AliFlowAnalysisWithMCEventPlane *mcep = NULL;
deebd72f	131
334e3256	132	//MCEP = monte carlo event plane
deebd72f	133	if (MCEP) {
	134	AliFlowAnalysisWithMCEventPlane *mcep = new AliFlowAnalysisWithMCEventPlane();
	135	mcep->Init();
	136	}
334e3256	137
334e3256	138	//QC = Q-cumulants
deebd72f	139	if(QC) {
deebd72f	140	AliFlowAnalysisWithQCumulants* qc = new AliFlowAnalysisWithQCumulants();
e085f1a9	141	qc->Init();
03a02aca	142	if(listWithWeights) qc->SetWeightsList(listWithWeights);
	143	if(usePhiWeights) qc->SetUsePhiWeights(usePhiWeights);
	144	if(usePtWeights) qc->SetUsePtWeights(usePtWeights);
	145	if(useEtaWeights) qc->SetUseEtaWeights(useEtaWeights);
deebd72f	146	}
334e3256	147
334e3256	148	//GFC = Generating Function Cumulants
deebd72f	149	if(GFC) {
deebd72f	150	AliFlowAnalysisWithCumulants* gfc = new AliFlowAnalysisWithCumulants();
e5e75b58	151	gfc->Init();
	152	if(listWithWeights) gfc->SetWeightsList(listWithWeights);
	153	if(usePhiWeights) gfc->SetUsePhiWeights(usePhiWeights);
	154	if(usePtWeights) gfc->SetUsePtWeights(usePtWeights);
	155	if(useEtaWeights) gfc->SetUseEtaWeights(useEtaWeights);
deebd72f	156	}
334e3256	157
334e3256	158	//FQD = Fitting q-distribution
deebd72f	159	if(FQD) {
deebd72f	160	AliFittingQDistribution* fqd = new AliFittingQDistribution();
14db9c04	161	fqd->Init();
	162	if(listWithWeights) fqd->SetWeightsList(listWithWeights);
	163	if(usePhiWeights) fqd->SetUsePhiWeights(usePhiWeights);
deebd72f	164	}
60ec66fb	165
	166	//SP = Scalar Product
	167	if(SP) {
	168	AliFlowAnalysisWithScalarProduct* sp = new AliFlowAnalysisWithScalarProduct();
	169	sp->Init();
	170	}
	171
334e3256	172	//LYZ1 = Lee-Yang Zeroes first run
deebd72f	173	if(LYZ1) {
	174	AliFlowAnalysisWithLeeYangZeros* lyz1 = new AliFlowAnalysisWithLeeYangZeros();
	175	lyz1->SetFirstRun(kTRUE);
	176	lyz1->SetUseSum(kTRUE);
	177	lyz1->Init();
	178	}
334e3256	179
334e3256	180	//LYZ2 = Lee-Yang Zeroes second run
deebd72f	181	if(LYZ2) {
	182	AliFlowAnalysisWithLeeYangZeros* lyz2 = new AliFlowAnalysisWithLeeYangZeros();
	183	// read the input file from the first run
	184	TString inputFileNameLYZ2 = "outputLYZ1analysis.root" ;
	185	TFile* inputFileLYZ2 = new TFile(inputFileNameLYZ2.Data(),"READ");
	186	if(!inputFileLYZ2 \|\| inputFileLYZ2->IsZombie()) {
ce78b98c	187	cerr << " ERROR: NO First Run file... " << endl ;
ce78b98c	188	break;
deebd72f	189	}
	190	else {
	191	TList* inputListLYZ2 = (TList*)inputFileLYZ2->Get("cobjLYZ1");
ace38bad	192	if (!inputListLYZ2) {cout<<"Input list is NULL pointer!"<<endl; break;}
deebd72f	193	else {
	194	cout<<"LYZ2 input file/list read..."<<endl;
	195	lyz2->SetFirstRunList(inputListLYZ2);
	196	lyz2->SetFirstRun(kFALSE);
	197	lyz2->SetUseSum(kTRUE);
	198	lyz2->Init();
334e3256	199	}
334e3256	200	}
deebd72f	201	}
deebd72f	202
334e3256	203	//LYZEP = Lee-Yang Zeroes event plane
deebd72f	204	if(LYZEP) {
	205	AliFlowLYZEventPlane* ep = new AliFlowLYZEventPlane() ;
	206	AliFlowAnalysisWithLYZEventPlane* lyzep = new AliFlowAnalysisWithLYZEventPlane();
	207	// read the input file from the second lyz run
	208	TString inputFileNameLYZEP = "outputLYZ2analysis.root" ;
	209	TFile* inputFileLYZEP = new TFile(inputFileNameLYZEP.Data(),"READ");
	210	if(!inputFileLYZEP \|\| inputFileLYZEP->IsZombie()) {
ce78b98c	211	cerr << " ERROR: NO Second Run file... " << endl ;
	212	break;
	213	}
deebd72f	214	else {
deebd72f	215	TList* inputListLYZEP = (TList*)inputFileLYZEP->Get("cobjLYZ2");
ace38bad	216	if (!inputListLYZEP) {cout<<"Input list is NULL pointer!"<<endl; break;}
deebd72f	217	else {
	218	cout<<"LYZEP input file/list read..."<<endl;
	219	ep ->SetSecondRunList(inputListLYZEP);
	220	lyzep->SetSecondRunList(inputListLYZEP);
	221	ep ->Init();
	222	lyzep->Init();
334e3256	223	}
334e3256	224	}
deebd72f	225	}
334e3256	226	//------------------------------------------------------------------------
334e3256	227
26c4cbb9	228
334e3256	229	//standard code
	230	Int_t fCount = 0;
	231	TString execDir(gSystem->pwd());
	232	TString targetDir(dir);
	233	TSystemDirectory* baseDir = new TSystemDirectory(".", dir);
	234	TList* dirList = baseDir->GetListOfFiles();
	235	if (!dirList) {
	236	cout << endl << "No input files in: " << targetDir.Data() << endl;
	237	break;
	238	}
	239	Int_t nDirs = dirList->GetEntries();
	240	cout<<endl;
	241	cout<<"Int_t nDirs = "<<nDirs<<endl;
	242	gSystem->cd(execDir);
	243
	244	for(Int_t iDir=0;iDir<nDirs;++iDir)
	245	{
	246	TSystemFile* presentDir = (TSystemFile*)dirList->At(iDir);
	247	if(!presentDir \|\| !presentDir->IsDirectory() \|\|
	248	strcmp(presentDir->GetName(), ".") == 0 \|\|
	249	strcmp(presentDir->GetName(), "..") == 0)
	250	{
	251	cout << endl;
	252	cout << "Directory (" << iDir << "): " << presentDir->GetName() <<
	253	" - Skipping ... " << endl;
	254	continue ;
	255	}
	256
	257	if(offset > 0) { --offset ; continue ; }
	258	if((aRuns > 0) && (fCount >= aRuns)) { break ; }
	259
	260	TString presentDirName(dir); // aDataDir
	261	presentDirName += presentDir->GetName();
	262	presentDirName += "/";
	263	//cerr<<" presentDirName = "<<presentDirName<<endl;
	264
	265	TString fileName = presentDirName;
	266	fileName += "galice.root";
	267	Long_t id, size, flags, modtime;
	268	if(gSystem->GetPathInfo(fileName.Data(),id,size,flags,modtime))
	269	{
	270	cout << " File : " << fileName << " does NOT exist ! - Skipping ... "
	271	<< endl;
	272	continue;
	273	}
deebd72f	274	// cout << endl ; cout << "Directory (" << iDir << "): " << presentDirName << " ... " << endl;
334e3256	275
	276	//loop (simulations in the present dir)
	277	TSystemDirectory* evtsDir = new TSystemDirectory(".",
	278	presentDirName.Data());
	279	TList* fileList = evtsDir->GetListOfFiles();
	280	Int_t nFiles = fileList->GetEntries();
deebd72f	281	//cout<<" Int_t nFiles = "<<nFiles<<endl;
334e3256	282	gSystem->cd(execDir);
	283	for(Int_t iFiles=0; iFiles<nFiles; ++iFiles)
	284	{
	285	TSystemFile* presentFile = (TSystemFile*) fileList->At(iFiles);
	286	TString presentFileName(presentDirName);
	287	presentFileName += presentFile->GetName();
	288
	289	if(!(presentFileName.Contains("Kinematics") &&
	290	presentFileName.Contains("root"))) { continue ; }
	291
deebd72f	292	//cout << " found: " << presentFileName.Data() << endl;
334e3256	293
	294	TFile* kineFile = new TFile(presentFileName.Data(), "READ");
	295	// kineFile->ls();
	296	Int_t nEvts = kineFile->GetNkeys() ;
deebd72f	297	//cout << " . found: " << nEvts << " KineTree(s) in " << presentFileName.Data() << endl;
334e3256	298	TList* kineEventsList = (TList*)kineFile->GetListOfKeys();
	299	TTree* kTree;
	300	TIter next(kineEventsList);
	301	TKey* key;
	302
	303	//loop over the events
	304	while( key=(TKey *)next() )
	305	{
	306	TDirectory* tDir = (TDirectory*)key->ReadObj();
	307	if(!tDir) break;
	308
	309	TString evtDir(tDir->GetName());
deebd72f	310	//cout << " . . found: " << tDir->GetName() << endl;
334e3256	311
	312	kTree = (TTree *)tDir->Get("TreeK");
	313	if(!kTree) break;
	314
	315	Int_t nPart = kTree->GetEntries();
deebd72f	316	//cout << " . . . kTree " << fCount << " has " << nPart << " particles " << endl;
334e3256	317
334e3256	318	//-----------------------------------------------------------
26c4cbb9	319	//fill and save the flow event
	320	AliFlowEventSimple *fEvent = fEventMaker->FillTracks(kTree, cutsInt, cutsDiff);
	321
d7671632	322	if(MCEP) mcep->Make(fEvent);
deebd72f	323	if(QC) qc->Make(fEvent);
	324	if(GFC) gfc->Make(fEvent);
	325	if(FQD) fqd->Make(fEvent);
	326	if(LYZ1) lyz1->Make(fEvent);
	327	if(LYZ2) lyz2->Make(fEvent);
	328	if(LYZEP) lyzep->Make(fEvent,ep);
	329	if(SP) sp->Make(fEvent);
334e3256	330	//-----------------------------------------------------------
	331
	332	fCount++;
deebd72f	333	cout << "# " << fCount << " events processed" << endl;
334e3256	334	delete kTree;
	335	delete fEvent;
	336	}
	337	delete kineFile ;
	338	}
	339	delete evtsDir ;
	340	}
	341
	342	//--------------------------------------------------------------
	343	//calculating and storing the final results of flow analysis
60ec66fb	344	if(MCEP) {mcep->Finish(); mcep->WriteHistograms("outputMCEPanalysis.root");}
ace38bad	345	if(SP) {sp->Finish(); sp->WriteHistograms("outputSPanalysis.root");}
deebd72f	346	if(QC) {qc->Finish(); qc->WriteHistograms("outputQCanalysis.root");}
	347	if(GFC) {gfc->Finish(); gfc->WriteHistograms("outputGFCanalysis.root");}
	348	if(FQD) {fqd->Finish(); fqd->WriteHistograms("outputFQDanalysis.root");}
	349	if(LYZ1) {lyz1->Finish(); lyz1->WriteHistograms("outputLYZ1analysis.root");}
	350	if(LYZ2) {lyz2->Finish(); lyz2->WriteHistograms("outputLYZ2analysis.root");}
	351	if(LYZEP) {lyzep->Finish(); lyzep->WriteHistograms("outputLYZEPanalysis.root");}
ace38bad	352
334e3256	353	//--------------------------------------------------------------
	354
	355	cout << endl;
ace38bad	356	cout << " Fini ... " << endl;
334e3256	357	cout << endl;
	358
	359	timer.Stop();
	360	cout << endl;
	361	timer.Print();
	362	}
21e694dd	363
	364	void SetupPar(char* pararchivename)
	365	{
	366	//Load par files, create analysis libraries
	367	//For testing, if par file already decompressed and modified
	368	//classes then do not decompress.
	369
	370	TString cdir(Form("%s", gSystem->WorkingDirectory() )) ;
	371	TString parpar(Form("%s.par", pararchivename)) ;
	372	if ( gSystem->AccessPathName(parpar.Data()) ) {
	373	gSystem->ChangeDirectory(gSystem->Getenv("ALICE_ROOT")) ;
	374	TString processline(Form(".! make %s", parpar.Data())) ;
	375	gROOT->ProcessLine(processline.Data()) ;
	376	gSystem->ChangeDirectory(cdir) ;
	377	processline = Form(".! mv /tmp/%s .", parpar.Data()) ;
	378	gROOT->ProcessLine(processline.Data()) ;
	379	}
	380	if ( gSystem->AccessPathName(pararchivename) ) {
	381	TString processline = Form(".! tar xvzf %s",parpar.Data()) ;
	382	gROOT->ProcessLine(processline.Data());
	383	}
	384
	385	TString ocwd = gSystem->WorkingDirectory();
	386	gSystem->ChangeDirectory(pararchivename);
	387
	388	// check for BUILD.sh and execute
	389	if (!gSystem->AccessPathName("PROOF-INF/BUILD.sh")) {
	390	printf("*******************************\n");
	391	printf("* Building PAR archive *\n");
	392	cout<<pararchivename<<endl;
	393	printf("*******************************\n");
	394
	395	if (gSystem->Exec("PROOF-INF/BUILD.sh")) {
	396	Error("runProcess","Cannot Build the PAR Archive! - Abort!");
	397	return -1;
	398	}
	399	}
	400	// check for SETUP.C and execute
	401	if (!gSystem->AccessPathName("PROOF-INF/SETUP.C")) {
	402	printf("*******************************\n");
	403	printf("* Setup PAR archive *\n");
	404	cout<<pararchivename<<endl;
	405	printf("*******************************\n");
	406	gROOT->Macro("PROOF-INF/SETUP.C");
	407	}
	408
	409	gSystem->ChangeDirectory(ocwd.Data());
	410	printf("Current dir: %s\n", ocwd.Data());
	411	}
	412
fdd8c18f	413	void LoadLibraries(const anaModes mode) {
	414
	415	//--------------------------------------
	416	// Load the needed libraries most of them already loaded by aliroot
	417	//--------------------------------------
	418	gSystem->Load("libTree.so");
	419	gSystem->Load("libGeom.so");
	420	gSystem->Load("libVMC.so");
	421	gSystem->Load("libXMLIO.so");
	422	gSystem->Load("libPhysics.so");
	423
	424	//----------------------------------------------------------
	425	// >>>>>>>>>>> Local mode <<<<<<<<<<<<<<
	426	//----------------------------------------------------------
	427	if (mode==mLocal) {
	428	//--------------------------------------------------------
	429	// If you want to use already compiled libraries
	430	// in the aliroot distribution
	431	//--------------------------------------------------------
	432	gSystem->Load("libSTEERBase");
	433	gSystem->Load("libESD");
	434	gSystem->Load("libAOD");
	435	gSystem->Load("libANALYSIS");
	436	gSystem->Load("libANALYSISalice");
	437	gSystem->Load("libCORRFW.so");
	438	cerr<<"libCORRFW.so loaded..."<<endl;
	439	gSystem->Load("libPWG2flowCommon.so");
	440	cerr<<"libPWG2flowCommon.so loaded..."<<endl;
	441	gSystem->Load("libPWG2flowTasks.so");
	442	cerr<<"libPWG2flowTasks.so loaded..."<<endl;
	443	}
	444
	445	else if (mode == mLocalPAR) {
	446	//--------------------------------------------------------
	447	//If you want to use root and par files from aliroot
	448	//--------------------------------------------------------
	449	//If you want to use root and par files from aliroot
	450	//--------------------------------------------------------
	451	SetupPar("STEERBase");
	452	SetupPar("ESD");
	453	SetupPar("AOD");
	454	SetupPar("ANALYSIS");
	455	SetupPar("ANALYSISalice");
	456	SetupPar("PWG2AOD");
	457	SetupPar("CORRFW");
	458	SetupPar("PWG2flowCommon");
	459	cerr<<"PWG2flowCommon.par loaded..."<<endl;
	460	SetupPar("PWG2flowTasks");
	461	cerr<<"PWG2flowTasks.par loaded..."<<endl;
	462	}
	463
	464	//---------------------------------------------------------
	465	// <<<<<<<<<< Source mode >>>>>>>>>>>>
	466	//---------------------------------------------------------
	467	else if (mode==mLocalSource) {
	468
	469	// In root inline compile
	470
	471
	472	// Constants
	473	gROOT->LoadMacro("AliFlowCommon/AliFlowCommonConstants.cxx+");
	474	gROOT->LoadMacro("AliFlowCommon/AliFlowLYZConstants.cxx+");
	475	gROOT->LoadMacro("AliFlowCommon/AliFlowCumuConstants.cxx+");
	476
477	// Flow event
478	gROOT->LoadMacro("AliFlowCommon/AliFlowVector.cxx+");
479	gROOT->LoadMacro("AliFlowCommon/AliFlowTrackSimple.cxx+");
480	gROOT->LoadMacro("AliFlowCommon/AliFlowEventSimple.cxx+");
481
482	// Cuts
483	gROOT->LoadMacro("AliFlowCommon/AliFlowTrackSimpleCuts.cxx+");
484
485	// Output histosgrams
486	gROOT->LoadMacro("AliFlowCommon/AliFlowCommonHist.cxx+");
487	gROOT->LoadMacro("AliFlowCommon/AliFlowCommonHistResults.cxx+");
488	gROOT->LoadMacro("AliFlowCommon/AliFlowLYZHist1.cxx+");
489	gROOT->LoadMacro("AliFlowCommon/AliFlowLYZHist2.cxx+");
490
491	// Functions needed for various methods
492	gROOT->LoadMacro("AliFlowCommon/AliCumulantsFunctions.cxx+");
493	gROOT->LoadMacro("AliFlowCommon/AliQCumulantsFunctions.cxx+");
494	gROOT->LoadMacro("AliFlowCommon/AliFittingFunctionsForQDistribution.cxx+");
495	gROOT->LoadMacro("AliFlowCommon/AliFlowLYZEventPlane.cxx+");
496
497	// Flow Analysis code for various methods
498	gROOT->LoadMacro("AliFlowCommon/AliFlowAnalysisWithMCEventPlane.cxx+");
499	gROOT->LoadMacro("AliFlowCommon/AliFlowAnalysisWithScalarProduct.cxx+");
500	gROOT->LoadMacro("AliFlowCommon/AliFlowAnalysisWithLYZEventPlane.cxx+");
501	gROOT->LoadMacro("AliFlowCommon/AliFlowAnalysisWithLeeYangZeros.cxx+");
502	gROOT->LoadMacro("AliFlowCommon/AliFlowAnalysisWithCumulants.cxx+");
503	gROOT->LoadMacro("AliFlowCommon/AliFlowAnalysisWithQCumulants.cxx+");
504	gROOT->LoadMacro("AliFlowCommon/AliFittingQDistribution.cxx+");
505
506	// Class to fill the FlowEvent without aliroot dependence
507	// can be found in the directory FlowEventMakers
508	gROOT->LoadMacro("FlowEventMakers/FlowEventSimpleMaker.cxx+");
509
510	cout << "finished loading macros!" << endl;
511
512	}
513
514	}
515
516