[u/mrichter/AliRoot.git] / PWGCF / FLOW / macros / runStarFlowAnalysis.C

//example script on what to do with the star events
//run e.g. like this:
//                    root runStarFlowAnalysis.C
// flow analysis method can be: (set to kTRUE or kFALSE)
Bool_t MCEP     = kTRUE;
Bool_t SP       = kTRUE;
Bool_t GFC      = kTRUE;
Bool_t QC       = kTRUE;
Bool_t FQD      = kTRUE;
Bool_t LYZ1SUM  = kTRUE;
Bool_t LYZ1PROD = kTRUE;
Bool_t LYZ2SUM  = kFALSE;
Bool_t LYZ2PROD = kFALSE;
Bool_t LYZEP    = kFALSE;
Bool_t MH       = kFALSE; // mixed harmonics 
Bool_t NL       = kFALSE; // nested loops
Bool_t MCEP_AH  = kFALSE; // MCEP in another harmonic 
//--------------------------------------------------------------------------------------

// Weights 
// use weights for Q-vector:
Bool_t usePhiWeights = kFALSE; // phi weights (correction for non-uniform azimuthal acceptance)
Bool_t usePtWeights  = kFALSE; // pt weights 
Bool_t useEtaWeights = kFALSE; // eta weights

// Define the range for eta subevents
Double_t minA = -0.9;
Double_t maxA = -0.01;
Double_t minB = 0.01;
Double_t maxB = 0.9;


void  runStarFlowAnalysis(const char* inputDataFiles="/Users/snelling/alice_data/jthomas/testData/")
{
  gSystem->Load("libTree");
  gSystem->Load("libVMC");
  gSystem->Load("libPhysics");
  gSystem->Load("libPWGflowBase");

  Int_t maxNumberOfEvents = 1000000;

  TStopwatch timer;
  timer.Start();
  
  if (LYZ1SUM && LYZ2SUM) {cout<<"WARNING: you cannot run LYZ1 and LYZ2 at the same time! LYZ2 needs the output from LYZ1."<<endl; exit(1); }
  if (LYZ1PROD && LYZ2PROD) {cout<<"WARNING: you cannot run LYZ1 and LYZ2 at the same time! LYZ2 needs the output from LYZ1."<<endl; exit(1); }
  if (LYZ2SUM && LYZEP) {cout<<"WARNING: you cannot run LYZ2 and LYZEP at the same time! LYZEP needs the output from LYZ2."<<endl; exit(1); }
  if (LYZ1SUM && LYZEP) {cout<<"WARNING: you cannot run LYZ1 and LYZEP at the same time! LYZEP needs the output from LYZ2."<<endl; exit(1); }


  //define reference particles
  AliStarTrackCuts* rpCuts = AliStarTrackCuts::StandardCuts();
  rpCuts->SetPtMin(0.05);
  rpCuts->SetPtMax(10.);

  //define particles of interest
  AliStarTrackCuts* poiCuts = AliStarTrackCuts::StandardCuts();
  poiCuts->SetPtMin(0.05);
  poiCuts->SetPtMax(10.);

  //define event cuts
  AliStarEventCuts* starEventCuts = AliStarEventCuts::StandardCuts();
  starEventCuts-> SetCentralityIDMax(3);
  starEventCuts-> SetCentralityIDMin(3);

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

  //---------------------------------------------------------------------------------------
  // Initialize all the flow methods for default analysis:  
  AliFlowAnalysisWithQCumulants *qc = NULL;
  AliFlowAnalysisWithCumulants *gfc = NULL;
  AliFlowAnalysisWithFittingQDistribution *fqd = NULL;
  AliFlowAnalysisWithLeeYangZeros *lyz1sum  = NULL;
  AliFlowAnalysisWithLeeYangZeros *lyz1prod = NULL;
  AliFlowAnalysisWithLeeYangZeros *lyz2sum  = NULL;
  AliFlowAnalysisWithLeeYangZeros *lyz2prod = NULL;
  AliFlowAnalysisWithLYZEventPlane *lyzep = NULL;
  AliFlowAnalysisWithScalarProduct *sp = NULL;
  AliFlowAnalysisWithMixedHarmonics *mh = NULL;
  AliFlowAnalysisWithNestedLoops *nl = NULL;
  AliFlowAnalysisWithMCEventPlane *mcep = NULL;   
  AliFlowAnalysisWithMCEventPlane *mcep_ah = NULL;   

// MCEP = monte carlo event plane
 if (MCEP) {
   AliFlowAnalysisWithMCEventPlane *mcep = new AliFlowAnalysisWithMCEventPlane();
   //mcep->SetHarmonic(2); // default is v2
   //mcep->SetFlowOfResonances(kTRUE);
   mcep->Init();
 }

 // MCEP = monte carlo event plane in another harmonic: 
 if(MCEP_AH)
 {
  AliFlowAnalysisWithMCEventPlane *mcep_ah = new AliFlowAnalysisWithMCEventPlane();
  mcep_ah->SetHarmonic(1);
  mcep_ah->Init();
 }
 
 // Mixed harmonics:
 if(MH) 
 {
  AliFlowAnalysisWithMixedHarmonics *mh = new AliFlowAnalysisWithMixedHarmonics();
  mh->SetHarmonic(1); // integer n in expression cos[n(2phi1-phi2-phi3)] = v2n*vn^2
  mh->SetMinMultiplicity(100); 
  mh->SetNoOfMultipicityBins(5);  
  mh->SetMultipicityBinWidth(200);   
  mh->Init(); 
 }
 
 // NL = nested loops:
 if(NL) {
   AliFlowAnalysisWithNestedLoops *nl = new AliFlowAnalysisWithNestedLoops();
   nl->Init();
 }

 // QC = Q-cumulants  
 if(QC) { 
   AliFlowAnalysisWithQCumulants* qc = new AliFlowAnalysisWithQCumulants();
   if(listWithWeights) qc->SetWeightsList(listWithWeights);
   if(usePhiWeights) qc->SetUsePhiWeights(usePhiWeights);
   if(usePtWeights) qc->SetUsePtWeights(usePtWeights);
   if(useEtaWeights) qc->SetUseEtaWeights(useEtaWeights);
   // qc->SetHarmonic(2); // default is v2
   // qc->SetApplyCorrectionForNUA(kTRUE); // default
   // qc->SetCalculate2DFlow(kFALSE); // default
   // qc->SetMultiplicityWeight("combinations"); // default
   // qc->SetMultiplicityWeight("unit");
   // qc->SetMultiplicityWeight("multiplicity");  
   qc->SetnBinsMult(10000);
   qc->SetMinMult(0);
   qc->SetMaxMult(10000);      
   qc->Init();  
 }
  
 // GFC = Generating Function Cumulants 
 if(GFC) {
   AliFlowAnalysisWithCumulants* gfc = new AliFlowAnalysisWithCumulants();
   if(listWithWeights) gfc->SetWeightsList(listWithWeights);
   if(usePhiWeights) gfc->SetUsePhiWeights(usePhiWeights);
   if(usePtWeights) gfc->SetUsePtWeights(usePtWeights);
   if(useEtaWeights) gfc->SetUseEtaWeights(useEtaWeights);
   // calculation vs multiplicity:
   gfc->SetCalculateVsMultiplicity(kFALSE);   
   gfc->SetnBinsMult(10000);
   gfc->SetMinMult(0);
   gfc->SetMaxMult(10000);   
   // tuning of interpolating parameters:
   gfc->SetTuneParameters(kFALSE);
   Double_t r0[10] = {1.8,1.9,2.0,2.1,2.2,2.3,2.4,2.5,2.6,2.7}; // up to 10 values allowed
   for(Int_t r=0;r<10;r++) {gfc->SetTuningR0(r0[r],r);}
   gfc->Init();
 }
 
 // FQD = Fitting q-distribution 
 if(FQD) {
   AliFlowAnalysisWithFittingQDistribution* fqd = new AliFlowAnalysisWithFittingQDistribution();
   if(listWithWeights) fqd->SetWeightsList(listWithWeights);
   if(usePhiWeights) fqd->SetUsePhiWeights(usePhiWeights);  
   fqd->Init();
 }
 
 // SP = Scalar Product 
 if(SP) {
   AliFlowAnalysisWithScalarProduct* sp = new AliFlowAnalysisWithScalarProduct();
   if(listWithWeights) sp->SetWeightsList(listWithWeights);
   sp->SetUsePhiWeights(usePhiWeights);  
   sp->Init();
 }

 // LYZ1 = Lee-Yang Zeroes first run
 if(LYZ1SUM) {
   AliFlowAnalysisWithLeeYangZeros* lyz1sum = new AliFlowAnalysisWithLeeYangZeros();
   lyz1sum->SetFirstRun(kTRUE);
   lyz1sum->SetUseSum(kTRUE);
   lyz1sum->Init();
 }
 if(LYZ1PROD) {
   AliFlowAnalysisWithLeeYangZeros* lyz1prod = new AliFlowAnalysisWithLeeYangZeros();
   lyz1prod->SetFirstRun(kTRUE);
   lyz1prod->SetUseSum(kFALSE);
   lyz1prod->Init();
 }
 // LYZ2 = Lee-Yang Zeroes second run
 if(LYZ2SUM) {
   AliFlowAnalysisWithLeeYangZeros* lyz2sum = new AliFlowAnalysisWithLeeYangZeros();
   // read the input file from the first run 
   TString inputFileNameLYZ2SUM = "outputLYZ1SUManalysis.root" ;
   TFile* inputFileLYZ2SUM = new TFile(inputFileNameLYZ2SUM.Data(),"READ");
   if(!inputFileLYZ2SUM || inputFileLYZ2SUM->IsZombie()) { 
     cerr << " ERROR: To run LYZ2SUM you need the output file from LYZ1SUM. This file is not there! Please run LYZ1SUM first." << endl ;
     break; 
   }
   else { 
     TList* inputListLYZ2SUM = (TList*)inputFileLYZ2SUM->Get("cobjLYZ1SUM");  
     if (!inputListLYZ2SUM) {cout<<"Input list LYZ2SUM is NULL pointer!"<<endl; break;}
     else {
       cout<<"LYZ2SUM input file/list read..."<<endl;
       lyz2sum->SetFirstRunList(inputListLYZ2SUM);
       lyz2sum->SetFirstRun(kFALSE);
       lyz2sum->SetUseSum(kTRUE);
       lyz2sum->Init();
     }
   }
 }
 if(LYZ2PROD) {
   AliFlowAnalysisWithLeeYangZeros* lyz2prod = new AliFlowAnalysisWithLeeYangZeros();
   // read the input file from the first run 
   TString inputFileNameLYZ2PROD = "outputLYZ1PRODanalysis.root" ;
   TFile* inputFileLYZ2PROD = new TFile(inputFileNameLYZ2PROD.Data(),"READ");
   if(!inputFileLYZ2PROD || inputFileLYZ2PROD->IsZombie()) { 
     cerr << " ERROR: To run LYZ2PROD you need the output file from LYZ1PROD. This file is not there! Please run LYZ1PROD first." << endl ;
     break; 
   }
   else { 
     TList* inputListLYZ2PROD = (TList*)inputFileLYZ2PROD->Get("cobjLYZ1PROD");  
     if (!inputListLYZ2PROD) {cout<<"Input list LYZ2PROD is NULL pointer!"<<endl; break;}
     else {
       cout<<"LYZ2PROD input file/list read..."<<endl;
       lyz2prod->SetFirstRunList(inputListLYZ2PROD);
       lyz2prod->SetFirstRun(kFALSE);
       lyz2prod->SetUseSum(kFALSE);
       lyz2prod->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 = "outputLYZ2SUManalysis.root" ;
   TFile* inputFileLYZEP = new TFile(inputFileNameLYZEP.Data(),"READ");
   if(!inputFileLYZEP || inputFileLYZEP->IsZombie()) { 
     cerr << " ERROR: To run LYZEP you need the output file from LYZ2SUM. This file is not there! Please run LYZ2SUM first." << endl ; 
     break;
   }
   else { 
     TList* inputListLYZEP = (TList*)inputFileLYZEP->Get("cobjLYZ2SUM");  
     if (!inputListLYZEP) {cout<<"Input list LYZEP is NULL pointer!"<<endl; break;}
     else {
       cout<<"LYZEP input file/list read..."<<endl;
       ep   ->SetSecondRunList(inputListLYZEP);
       lyzep->SetSecondRunList(inputListLYZEP);
       ep   ->Init();
       lyzep->Init();
     }
   }
 }
 //---------------------------------------------------------------------------------------


  Int_t i=0;
  AliStarEventReader starReader(inputDataFiles) ;
  while ( starReader.GetNextEvent() )                                // Get next event
  {
    AliStarEvent* starEvent = starReader.GetEvent();
    if ( !starEventCuts->PassesCuts(starEvent) ) continue;              // Test if the event is good

    AliFlowEventSimple* flowEvent = new AliFlowEventStar(starEvent,rpCuts,poiCuts);  // make a flow event from a star event (aka "the magic")
    flowEvent->TagSubeventsInEta(minA, maxA, minB, maxB );

    /////analysis here////////////////

    // do flow analysis for various methods
    if(MCEP)    mcep->Make(flowEvent);
    if(QC)      qc->Make(flowEvent);
    if(GFC)     gfc->Make(flowEvent);
    if(FQD)     fqd->Make(flowEvent);
    if(LYZ1SUM) lyz1sum->Make(flowEvent);
    if(LYZ1PROD)lyz1prod->Make(flowEvent);
    if(LYZ2SUM) lyz2sum->Make(flowEvent);
    if(LYZ2PROD)lyz2prod->Make(flowEvent);
    if(LYZEP)   lyzep->Make(flowEvent,ep);
    if(SP)      sp->Make(flowEvent);	      
    if(MH)      mh->Make(flowEvent);
    if(NL)      nl->Make(flowEvent);
    if(MCEP_AH) mcep_ah->Make(flowEvent);

    //////////////////////////////////

    flowEvent->Print();

    delete flowEvent;

    i++;
    if (i>maxNumberOfEvents) break;
  }

  //---------------------------------------------------------------------------------------  
  // create a new file which will hold the final results of all methods:
  TString outputFileName = "AnalysisResults.root";  
  TFile *outputFile = new TFile(outputFileName.Data(),"RECREATE");
  // create a new file for each method wich will hold list with final results:
  const Int_t nMethods = 13;
  TString method[nMethods] = {"MCEP","SP","GFC","QC","FQD","LYZ1SUM","LYZ1PROD","LYZ2SUM","LYZ2PROD","LYZEP","MH","NL","MCEP_AH"};
  TDirectoryFile *dirFileFinal[nMethods] = {NULL};
  TString fileName[nMethods]; 
  for(Int_t i=0;i<nMethods;i++)
    {
      // form a file name for each method:
      fileName[i]+="output";
      fileName[i]+=method[i].Data();
      fileName[i]+="analysis";
      dirFileFinal[i] = new TDirectoryFile(fileName[i].Data(),fileName[i].Data());
    } 
  
  // calculating and storing the final results of default flow analysis:
  if(MCEP)    {mcep->Finish();    mcep->WriteHistograms(dirFileFinal[0]);}
  if(SP)      {sp->Finish();      sp->WriteHistograms(dirFileFinal[1]);}
  if(GFC)     {gfc->Finish();     gfc->WriteHistograms(dirFileFinal[2]);}
  if(QC)      {qc->Finish();      qc->WriteHistograms(dirFileFinal[3]);}
  if(FQD)     {fqd->Finish();     fqd->WriteHistograms(dirFileFinal[4]);}
  if(LYZ1SUM) {lyz1sum->Finish(); lyz1sum->WriteHistograms(dirFileFinal[5]);}
  if(LYZ1PROD){lyz1prod->Finish();lyz1prod->WriteHistograms(dirFileFinal[6]);}
  if(LYZ2SUM) {lyz2sum->Finish(); lyz2sum->WriteHistograms(dirFileFinal[7]);}
  if(LYZ2PROD){lyz2prod->Finish();lyz2prod->WriteHistograms(dirFileFinal[8]);}
  if(LYZEP)   {lyzep->Finish();   lyzep->WriteHistograms(dirFileFinal[9]);}
  if(MH)      {mh->Finish();      mh->WriteHistograms(dirFileFinal[10]);}
  if(NL)      {nl->Finish();      nl->WriteHistograms(dirFileFinal[11]);}
  if(MCEP_AH) {mcep_ah->Finish(); mcep_ah->WriteHistograms(dirFileFinal[12]);}
  //---------------------------------------------------------------------------------------  
  
  outputFile->Close();
  delete outputFile;
  
  cout<<endl;
  cout<<endl;
  cout<<" ---- LANDED SUCCESSFULLY ---- "<<endl;
  cout<<endl; 
  
  timer.Stop();
  cout << endl;
  timer.Print();
  
  delete rpCuts;
  delete poiCuts;
  delete starEventCuts;
}
Commit	Line	Data
6167f59c	1	//example script on what to do with the star events
	2	//run e.g. like this:
	3	// root runStarFlowAnalysis.C
b80d67e1	4	// flow analysis method can be: (set to kTRUE or kFALSE)
b80d67e1	5	Bool_t MCEP = kTRUE;
6167f59c	6	Bool_t SP = kTRUE;
6167f59c	7	Bool_t GFC = kTRUE;
	8	Bool_t QC = kTRUE;
	9	Bool_t FQD = kTRUE;
b80d67e1	10	Bool_t LYZ1SUM = kTRUE;
	11	Bool_t LYZ1PROD = kTRUE;
	12	Bool_t LYZ2SUM = kFALSE;
	13	Bool_t LYZ2PROD = kFALSE;
	14	Bool_t LYZEP = kFALSE;
	15	Bool_t MH = kFALSE; // mixed harmonics
	16	Bool_t NL = kFALSE; // nested loops
	17	Bool_t MCEP_AH = kFALSE; // MCEP in another harmonic
6167f59c	18	//--------------------------------------------------------------------------------------
	19
	20	// Weights
b80d67e1	21	// use weights for Q-vector:
	22	Bool_t usePhiWeights = kFALSE; // phi weights (correction for non-uniform azimuthal acceptance)
	23	Bool_t usePtWeights = kFALSE; // pt weights
	24	Bool_t useEtaWeights = kFALSE; // eta weights
	25
	26	// Define the range for eta subevents
	27	Double_t minA = -0.9;
	28	Double_t maxA = -0.01;
	29	Double_t minB = 0.01;
	30	Double_t maxB = 0.9;
6167f59c	31
6167f59c	32
afc49a4b	33	void runStarFlowAnalysis(const char* inputDataFiles="/Users/snelling/alice_data/jthomas/testData/")
6167f59c	34	{
4070f709	35	gSystem->Load("libTree");
	36	gSystem->Load("libVMC");
	37	gSystem->Load("libPhysics");
2e311896	38	gSystem->Load("libPWGflowBase");
6167f59c	39
	40	Int_t maxNumberOfEvents = 1000000;
	41
	42	TStopwatch timer;
	43	timer.Start();
	44
	45	if (LYZ1SUM && LYZ2SUM) {cout<<"WARNING: you cannot run LYZ1 and LYZ2 at the same time! LYZ2 needs the output from LYZ1."<<endl; exit(1); }
	46	if (LYZ1PROD && LYZ2PROD) {cout<<"WARNING: you cannot run LYZ1 and LYZ2 at the same time! LYZ2 needs the output from LYZ1."<<endl; exit(1); }
	47	if (LYZ2SUM && LYZEP) {cout<<"WARNING: you cannot run LYZ2 and LYZEP at the same time! LYZEP needs the output from LYZ2."<<endl; exit(1); }
	48	if (LYZ1SUM && LYZEP) {cout<<"WARNING: you cannot run LYZ1 and LYZEP at the same time! LYZEP needs the output from LYZ2."<<endl; exit(1); }
	49
	50
	51
	52	//define reference particles
	53	AliStarTrackCuts* rpCuts = AliStarTrackCuts::StandardCuts();
cdfe7ebe	54	rpCuts->SetPtMin(0.05);
cdfe7ebe	55	rpCuts->SetPtMax(10.);
6167f59c	56
	57	//define particles of interest
	58	AliStarTrackCuts* poiCuts = AliStarTrackCuts::StandardCuts();
	59	poiCuts->SetPtMin(0.05);
cdfe7ebe	60	poiCuts->SetPtMax(10.);
6167f59c	61
	62	//define event cuts
	63	AliStarEventCuts* starEventCuts = AliStarEventCuts::StandardCuts();
cdfe7ebe	64	starEventCuts-> SetCentralityIDMax(3);
cdfe7ebe	65	starEventCuts-> SetCentralityIDMin(3);
6167f59c	66
	67	//if the weights are used:
	68	TFile *fileWithWeights = NULL;
	69	TList *listWithWeights = NULL;
	70
	71	if(usePhiWeights\|\|usePtWeights\|\|useEtaWeights) {
	72	fileWithWeights = TFile::Open("weights.root","READ");
	73	if(fileWithWeights) {
	74	listWithWeights = (TList*)fileWithWeights->Get("weights");
	75	}
	76	else
	77	{cout << " WARNING: the file <weights.root> with weights from the previous run was not found."<<endl;
	78	break;
	79	}
	80	}
	81
b80d67e1	82	//---------------------------------------------------------------------------------------
b80d67e1	83	// Initialize all the flow methods for default analysis:
6167f59c	84	AliFlowAnalysisWithQCumulants *qc = NULL;
	85	AliFlowAnalysisWithCumulants *gfc = NULL;
	86	AliFlowAnalysisWithFittingQDistribution *fqd = NULL;
b80d67e1	87	AliFlowAnalysisWithLeeYangZeros *lyz1sum = NULL;
6167f59c	88	AliFlowAnalysisWithLeeYangZeros *lyz1prod = NULL;
b80d67e1	89	AliFlowAnalysisWithLeeYangZeros *lyz2sum = NULL;
6167f59c	90	AliFlowAnalysisWithLeeYangZeros *lyz2prod = NULL;
	91	AliFlowAnalysisWithLYZEventPlane *lyzep = NULL;
	92	AliFlowAnalysisWithScalarProduct *sp = NULL;
6167f59c	93	AliFlowAnalysisWithMixedHarmonics *mh = NULL;
6167f59c	94	AliFlowAnalysisWithNestedLoops *nl = NULL;
b80d67e1	95	AliFlowAnalysisWithMCEventPlane *mcep = NULL;
b80d67e1	96	AliFlowAnalysisWithMCEventPlane *mcep_ah = NULL;
6167f59c	97
b80d67e1	98	// MCEP = monte carlo event plane
	99	if (MCEP) {
	100	AliFlowAnalysisWithMCEventPlane *mcep = new AliFlowAnalysisWithMCEventPlane();
	101	//mcep->SetHarmonic(2); // default is v2
	102	//mcep->SetFlowOfResonances(kTRUE);
	103	mcep->Init();
	104	}
6167f59c	105
b80d67e1	106	// MCEP = monte carlo event plane in another harmonic:
	107	if(MCEP_AH)
	108	{
	109	AliFlowAnalysisWithMCEventPlane *mcep_ah = new AliFlowAnalysisWithMCEventPlane();
	110	mcep_ah->SetHarmonic(1);
	111	mcep_ah->Init();
	112	}
	113
	114	// Mixed harmonics:
	115	if(MH)
	116	{
	117	AliFlowAnalysisWithMixedHarmonics *mh = new AliFlowAnalysisWithMixedHarmonics();
	118	mh->SetHarmonic(1); // integer n in expression cos[n(2phi1-phi2-phi3)] = v2n*vn^2
	119	mh->SetMinMultiplicity(100);
	120	mh->SetNoOfMultipicityBins(5);
	121	mh->SetMultipicityBinWidth(200);
	122	mh->Init();
	123	}
	124
	125	// NL = nested loops:
	126	if(NL) {
	127	AliFlowAnalysisWithNestedLoops *nl = new AliFlowAnalysisWithNestedLoops();
	128	nl->Init();
	129	}
	130
	131	// QC = Q-cumulants
	132	if(QC) {
	133	AliFlowAnalysisWithQCumulants* qc = new AliFlowAnalysisWithQCumulants();
	134	if(listWithWeights) qc->SetWeightsList(listWithWeights);
	135	if(usePhiWeights) qc->SetUsePhiWeights(usePhiWeights);
	136	if(usePtWeights) qc->SetUsePtWeights(usePtWeights);
	137	if(useEtaWeights) qc->SetUseEtaWeights(useEtaWeights);
	138	// qc->SetHarmonic(2); // default is v2
	139	// qc->SetApplyCorrectionForNUA(kTRUE); // default
	140	// qc->SetCalculate2DFlow(kFALSE); // default
	141	// qc->SetMultiplicityWeight("combinations"); // default
	142	// qc->SetMultiplicityWeight("unit");
	143	// qc->SetMultiplicityWeight("multiplicity");
	144	qc->SetnBinsMult(10000);
	145	qc->SetMinMult(0);
	146	qc->SetMaxMult(10000);
	147	qc->Init();
	148	}
6167f59c	149
b80d67e1	150	// GFC = Generating Function Cumulants
	151	if(GFC) {
	152	AliFlowAnalysisWithCumulants* gfc = new AliFlowAnalysisWithCumulants();
	153	if(listWithWeights) gfc->SetWeightsList(listWithWeights);
	154	if(usePhiWeights) gfc->SetUsePhiWeights(usePhiWeights);
	155	if(usePtWeights) gfc->SetUsePtWeights(usePtWeights);
	156	if(useEtaWeights) gfc->SetUseEtaWeights(useEtaWeights);
	157	// calculation vs multiplicity:
	158	gfc->SetCalculateVsMultiplicity(kFALSE);
	159	gfc->SetnBinsMult(10000);
	160	gfc->SetMinMult(0);
	161	gfc->SetMaxMult(10000);
	162	// tuning of interpolating parameters:
	163	gfc->SetTuneParameters(kFALSE);
	164	Double_t r0[10] = {1.8,1.9,2.0,2.1,2.2,2.3,2.4,2.5,2.6,2.7}; // up to 10 values allowed
	165	for(Int_t r=0;r<10;r++) {gfc->SetTuningR0(r0[r],r);}
	166	gfc->Init();
	167	}
	168
	169	// FQD = Fitting q-distribution
	170	if(FQD) {
	171	AliFlowAnalysisWithFittingQDistribution* fqd = new AliFlowAnalysisWithFittingQDistribution();
	172	if(listWithWeights) fqd->SetWeightsList(listWithWeights);
	173	if(usePhiWeights) fqd->SetUsePhiWeights(usePhiWeights);
	174	fqd->Init();
	175	}
	176
	177	// SP = Scalar Product
	178	if(SP) {
	179	AliFlowAnalysisWithScalarProduct* sp = new AliFlowAnalysisWithScalarProduct();
	180	if(listWithWeights) sp->SetWeightsList(listWithWeights);
	181	sp->SetUsePhiWeights(usePhiWeights);
	182	sp->Init();
	183	}
6167f59c	184
b80d67e1	185	// LYZ1 = Lee-Yang Zeroes first run
	186	if(LYZ1SUM) {
	187	AliFlowAnalysisWithLeeYangZeros* lyz1sum = new AliFlowAnalysisWithLeeYangZeros();
	188	lyz1sum->SetFirstRun(kTRUE);
	189	lyz1sum->SetUseSum(kTRUE);
	190	lyz1sum->Init();
	191	}
	192	if(LYZ1PROD) {
	193	AliFlowAnalysisWithLeeYangZeros* lyz1prod = new AliFlowAnalysisWithLeeYangZeros();
	194	lyz1prod->SetFirstRun(kTRUE);
	195	lyz1prod->SetUseSum(kFALSE);
	196	lyz1prod->Init();
	197	}
	198	// LYZ2 = Lee-Yang Zeroes second run
	199	if(LYZ2SUM) {
	200	AliFlowAnalysisWithLeeYangZeros* lyz2sum = new AliFlowAnalysisWithLeeYangZeros();
	201	// read the input file from the first run
	202	TString inputFileNameLYZ2SUM = "outputLYZ1SUManalysis.root" ;
	203	TFile* inputFileLYZ2SUM = new TFile(inputFileNameLYZ2SUM.Data(),"READ");
	204	if(!inputFileLYZ2SUM \|\| inputFileLYZ2SUM->IsZombie()) {
	205	cerr << " ERROR: To run LYZ2SUM you need the output file from LYZ1SUM. This file is not there! Please run LYZ1SUM first." << endl ;
	206	break;
	207	}
	208	else {
	209	TList* inputListLYZ2SUM = (TList*)inputFileLYZ2SUM->Get("cobjLYZ1SUM");
	210	if (!inputListLYZ2SUM) {cout<<"Input list LYZ2SUM is NULL pointer!"<<endl; break;}
	211	else {
	212	cout<<"LYZ2SUM input file/list read..."<<endl;
	213	lyz2sum->SetFirstRunList(inputListLYZ2SUM);
	214	lyz2sum->SetFirstRun(kFALSE);
	215	lyz2sum->SetUseSum(kTRUE);
	216	lyz2sum->Init();
	217	}
	218	}
	219	}
	220	if(LYZ2PROD) {
	221	AliFlowAnalysisWithLeeYangZeros* lyz2prod = new AliFlowAnalysisWithLeeYangZeros();
	222	// read the input file from the first run
	223	TString inputFileNameLYZ2PROD = "outputLYZ1PRODanalysis.root" ;
	224	TFile* inputFileLYZ2PROD = new TFile(inputFileNameLYZ2PROD.Data(),"READ");
	225	if(!inputFileLYZ2PROD \|\| inputFileLYZ2PROD->IsZombie()) {
	226	cerr << " ERROR: To run LYZ2PROD you need the output file from LYZ1PROD. This file is not there! Please run LYZ1PROD first." << endl ;
	227	break;
	228	}
	229	else {
	230	TList* inputListLYZ2PROD = (TList*)inputFileLYZ2PROD->Get("cobjLYZ1PROD");
	231	if (!inputListLYZ2PROD) {cout<<"Input list LYZ2PROD is NULL pointer!"<<endl; break;}
	232	else {
	233	cout<<"LYZ2PROD input file/list read..."<<endl;
	234	lyz2prod->SetFirstRunList(inputListLYZ2PROD);
	235	lyz2prod->SetFirstRun(kFALSE);
	236	lyz2prod->SetUseSum(kFALSE);
	237	lyz2prod->Init();
	238	}
	239	}
	240	}
6167f59c	241
b80d67e1	242	// LYZEP = Lee-Yang Zeroes event plane
	243	if(LYZEP) {
	244	AliFlowLYZEventPlane* ep = new AliFlowLYZEventPlane() ;
	245	AliFlowAnalysisWithLYZEventPlane* lyzep = new AliFlowAnalysisWithLYZEventPlane();
	246	// read the input file from the second lyz run
	247	TString inputFileNameLYZEP = "outputLYZ2SUManalysis.root" ;
	248	TFile* inputFileLYZEP = new TFile(inputFileNameLYZEP.Data(),"READ");
	249	if(!inputFileLYZEP \|\| inputFileLYZEP->IsZombie()) {
	250	cerr << " ERROR: To run LYZEP you need the output file from LYZ2SUM. This file is not there! Please run LYZ2SUM first." << endl ;
	251	break;
	252	}
	253	else {
	254	TList* inputListLYZEP = (TList*)inputFileLYZEP->Get("cobjLYZ2SUM");
	255	if (!inputListLYZEP) {cout<<"Input list LYZEP is NULL pointer!"<<endl; break;}
	256	else {
	257	cout<<"LYZEP input file/list read..."<<endl;
	258	ep ->SetSecondRunList(inputListLYZEP);
	259	lyzep->SetSecondRunList(inputListLYZEP);
	260	ep ->Init();
	261	lyzep->Init();
	262	}
	263	}
	264	}
	265	//---------------------------------------------------------------------------------------
6167f59c	266
	267
	268	Int_t i=0;
afc49a4b	269	AliStarEventReader starReader(inputDataFiles) ;
6167f59c	270	while ( starReader.GetNextEvent() ) // Get next event
	271	{
	272	AliStarEvent* starEvent = starReader.GetEvent();
	273	if ( !starEventCuts->PassesCuts(starEvent) ) continue; // Test if the event is good
	274
	275	AliFlowEventSimple* flowEvent = new AliFlowEventStar(starEvent,rpCuts,poiCuts); // make a flow event from a star event (aka "the magic")
b80d67e1	276	flowEvent->TagSubeventsInEta(minA, maxA, minB, maxB );
6167f59c	277
	278	/////analysis here////////////////
	279
	280	// do flow analysis for various methods
	281	if(MCEP) mcep->Make(flowEvent);
	282	if(QC) qc->Make(flowEvent);
	283	if(GFC) gfc->Make(flowEvent);
	284	if(FQD) fqd->Make(flowEvent);
	285	if(LYZ1SUM) lyz1sum->Make(flowEvent);
	286	if(LYZ1PROD)lyz1prod->Make(flowEvent);
	287	if(LYZ2SUM) lyz2sum->Make(flowEvent);
	288	if(LYZ2PROD)lyz2prod->Make(flowEvent);
	289	if(LYZEP) lyzep->Make(flowEvent,ep);
	290	if(SP) sp->Make(flowEvent);
	291	if(MH) mh->Make(flowEvent);
	292	if(NL) nl->Make(flowEvent);
b80d67e1	293	if(MCEP_AH) mcep_ah->Make(flowEvent);
6167f59c	294
	295	//////////////////////////////////
	296
b80d67e1	297	flowEvent->Print();
6167f59c	298
	299	delete flowEvent;
	300
	301	i++;
	302	if (i>maxNumberOfEvents) break;
	303	}
	304
b80d67e1	305	//---------------------------------------------------------------------------------------
6167f59c	306	// create a new file which will hold the final results of all methods:
	307	TString outputFileName = "AnalysisResults.root";
	308	TFile *outputFile = new TFile(outputFileName.Data(),"RECREATE");
	309	// create a new file for each method wich will hold list with final results:
b80d67e1	310	const Int_t nMethods = 13;
b80d67e1	311	TString method[nMethods] = {"MCEP","SP","GFC","QC","FQD","LYZ1SUM","LYZ1PROD","LYZ2SUM","LYZ2PROD","LYZEP","MH","NL","MCEP_AH"};
6167f59c	312	TDirectoryFile *dirFileFinal[nMethods] = {NULL};
b80d67e1	313	TString fileName[nMethods];
6167f59c	314	for(Int_t i=0;i<nMethods;i++)
	315	{
	316	// form a file name for each method:
b80d67e1	317	fileName[i]+="output";
	318	fileName[i]+=method[i].Data();
	319	fileName[i]+="analysis";
	320	dirFileFinal[i] = new TDirectoryFile(fileName[i].Data(),fileName[i].Data());
6167f59c	321	}
	322
	323	// calculating and storing the final results of default flow analysis:
	324	if(MCEP) {mcep->Finish(); mcep->WriteHistograms(dirFileFinal[0]);}
	325	if(SP) {sp->Finish(); sp->WriteHistograms(dirFileFinal[1]);}
	326	if(GFC) {gfc->Finish(); gfc->WriteHistograms(dirFileFinal[2]);}
	327	if(QC) {qc->Finish(); qc->WriteHistograms(dirFileFinal[3]);}
	328	if(FQD) {fqd->Finish(); fqd->WriteHistograms(dirFileFinal[4]);}
	329	if(LYZ1SUM) {lyz1sum->Finish(); lyz1sum->WriteHistograms(dirFileFinal[5]);}
	330	if(LYZ1PROD){lyz1prod->Finish();lyz1prod->WriteHistograms(dirFileFinal[6]);}
	331	if(LYZ2SUM) {lyz2sum->Finish(); lyz2sum->WriteHistograms(dirFileFinal[7]);}
	332	if(LYZ2PROD){lyz2prod->Finish();lyz2prod->WriteHistograms(dirFileFinal[8]);}
	333	if(LYZEP) {lyzep->Finish(); lyzep->WriteHistograms(dirFileFinal[9]);}
	334	if(MH) {mh->Finish(); mh->WriteHistograms(dirFileFinal[10]);}
	335	if(NL) {nl->Finish(); nl->WriteHistograms(dirFileFinal[11]);}
b80d67e1	336	if(MCEP_AH) {mcep_ah->Finish(); mcep_ah->WriteHistograms(dirFileFinal[12]);}
6167f59c	337	//---------------------------------------------------------------------------------------
	338
	339	outputFile->Close();
	340	delete outputFile;
	341
b80d67e1	342	cout<<endl;
	343	cout<<endl;
	344	cout<<" ---- LANDED SUCCESSFULLY ---- "<<endl;
	345	cout<<endl;
6167f59c	346
	347	timer.Stop();
	348	cout << endl;
	349	timer.Print();
	350
	351	delete rpCuts;
	352	delete poiCuts;
	353	delete starEventCuts;
	354	}