PWG2/FLOW/macros/runFlowAnalysisOnTheFly.C

   1 #include "TStopwatch.h"
   2 #include "TObjArray"
   3 #include "Riostream.h"
   4 #include "TFile.h"
   5
   6 //--------------------------------------------------------------------------------------
   7 // RUN SETTINGS
   8 // flow analysis method can be: (set to kTRUE or kFALSE)
   9 Bool_t SP       = kTRUE;
  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 GFC      = kTRUE;
  16 Bool_t QC       = kTRUE;
  17 Bool_t FQD      = kTRUE;
  18 Bool_t MCEP     = kTRUE;
  19 //--------------------------------------------------------------------------------------
  20
  21 // Weights
  22 // use weights for Q-vector:
  23 Bool_t usePhiWeights = kFALSE; // phi weights (correction for non-uniform azimuthal acceptance)
  24 Bool_t usePtWeights  = kFALSE; // pt weights
  25 Bool_t useEtaWeights = kFALSE; // eta weights
  26
  27 // Run same flow analysis method but with different settings/aims
  28 // You will have to label each setting/aim with your own label (see examples bellow):
  29 Bool_t GFC_Additional_Analysis = kFALSE;
  30 Bool_t QC_Additional_Analysis  = kFALSE;
  31 Bool_t FQD_Additional_Analysis = kFALSE;
  32
  33 // Parameters for the simulation of events 'on the fly':
  34 Bool_t bSameSeed = kFALSE; // use always the same seed for random generators.
  35                            // usage of same seed (kTRUE) is relevant in two cases:
  36                            // 1.) If you want to use LYZ method to calcualte differential flow;
  37                            // 2.) If you want to use phi weights for GFC, QC and FQD
  38
  39 Bool_t bConstantHarmonics = kTRUE; // harmonics V1, V2, V4... are constant (kTRUE) or functions of pt and eta (kFALSE)
  40
  41 Int_t iLoops = 1; // number of times to use each track (to simulate nonflow)
  42
  43 Bool_t bMultDistrOfRPsIsGauss = kTRUE; // 1.) if kTRUE  = multiplicitiy of RPs is sampled e-b-e from Gaussian distribution with
  44                                         //                 mean = iMultiplicityOfRP and spread = dMultiplicitySpreadOfRP
  45                                         // 2.) if kFALSE = multiplicitiy of RPs is sampled e-b-e uniformly from
  46                                         //                 interval [iMinMultOfRP,iMaxMultOfRP]
  47                                         // 3.) for a fixed multiplicity use Gaussian with zero spread or use uniform with iMinMult=iMaxMult
  48
  49 Bool_t bV2DistrOfRPsIsGauss = kTRUE; // 1.) if kTRUE  = elliptic flow of RPs is sampled e-b-e from Gaussian distribution with
  50                                       //                 mean = dV2RP and spread = dV2SpreadRP
  51                                       // 2.) if kFALSE = elliptic flow of RPs is sampled e-b-e uniformly from
  52                                       //                 interval [dMinV2RP,dMaxV2RP]
  53                                       // 3.) for a fixed elliptic flow use Gaussian with zero spread or use uniform with dMinV2RP=dMaxV2RP
  54
  55 Bool_t uniformAcceptance = kTRUE; // 1.) if kTRUE = detectors has uniform azimuthal acceptance
  56                                   // 2.) if kFALSE = you will simulate detector with non-uniform acceptance in one or two sectors.
  57                                   //                 For each of two sectors you specify phi_min, phi_max and probability p. Then all particles
  58                                   //                 going in direction phi_min < phi < phi_max will be taken with probability p. If p = 0, that
  59                                   //                 sector is blocked. Set bellow phimin1, phimax1, p1 for the first sector and phimin2, phimax2, p2
  60                                   //                 for the second sector. If you set phimin2 = phimax2 = p2 = 0, only first non-uniform sector is
  61                                   //                 simulated.
  62
  63 Int_t iMultiplicityOfRP = 500;        // mean multiplicity of RPs (if sampled from Gaussian)
  64 Double_t dMultiplicitySpreadOfRP = 0; // multiplicity spread of RPs (if sampled from Gaussian)
  65 Int_t iMinMultOfRP = 400;             // minimal multiplicity of RPs (if sampled uniformly)
  66 Int_t iMaxMultOfRP = 600;             // maximal multiplicity of RPs (if sampled uniformly)
  67
  68 Double_t dTemperatureOfRP = 0.44; // 'temperature' of RPs in GeV/c (increase this parameter to get more high pt RPs)
  69
  70 //......................................................................................
  71 // if you use (pt,eta) dependent harmonics (bConstantHarmonics = kFALSE):
  72 Double_t dPtCutOff = 2.0; // V2(pt) is linear up to pt = 2 GeV and for pt > 2 GeV it is constant: V2(pt) = dVRPMax
  73 Double_t dV2RPMax = 0.20; // maximum value of V2(pt) for pt >= 2GeV
  74 //......................................................................................
  75
  76 //......................................................................................
  77 // if you use constant harmonics (bConstantHarmonics = kTRUE) (i.e. no pt dependence):
  78 Double_t dV2RP = 0.05;       // elliptic flow of RPs (if sampled from Gaussian)
  79 Double_t dV2SpreadRP = 0.0;  // elliptic flow spread of RPs (if sampled from Gaussian)
  80 Double_t dMinV2RP = 0.04;    // minimal elliptic flow of RPs (if sampled uniformly)
  81 Double_t dMaxV2RP = 0.06;    // maximal elliptic flow of RPs (if sampled uniformly)
  82
  83 Double_t dV1RP = 0.0; // directed flow of RPs
  84 Double_t dV1SpreadRP = 0.0; // directed flow spread of RPs
  85
  86 Double_t dV4RP = 0.0; // harmonic V4 of RPs (to be improved: name needed)
  87 Double_t dV4SpreadRP = 0.0; // harmonic V4's spread of RPs (to be improved: name needed)
  88 //......................................................................................
  89
  90 //......................................................................................
  91 // settings for non-uniform acceptance:
  92 // Remark: set the angles in degrees from interval [0,360] and probability from interval [0,1]
  93
  94 // 1st non-uniform sector:
  95 Double_t phimin1 = 60;  // first non-uniform sector starts at this azimuth
  96 Double_t phimax1 = 120; // first non-uniform sector ends at this azimuth
  97 Double_t p1 = 0.33;        // e.g. if p1 = 0 all particles emitted in phimin1 < phi < phimax1 are blocked
  98                         // e.g. if p1 = 0.5 half of the particles emitted in phimin1 < phi < phimax1 are blocked
  99
 100 // 2nd non-uniform sector (Remark: if you do NOT want to simulate this sector, set phimin2 = phimax2 = p2 = 0):
 101 Double_t phimin2 = 0.0; // second non-uniform sector starts at this azimuth (make sure phimin2 > phimax1 !!!!)
 102 Double_t phimax2 = 0.0; // second non-uniform sector ends at this azimuth
 103 Double_t p2 = 0.0;
 104 //......................................................................................
 105
 106 enum anaModes {mLocal,mLocalSource,mLocalPAR};
 107 // mLocal: Analyze data on your computer using aliroot
 108 // mLocalPAR: Analyze data on your computer using root + PAR files
 109 // mLocalSource: Analyze data on your computer using root + source files
 110
 111 int runFlowAnalysisOnTheFly(Int_t mode=mLocal, Int_t nEvts=440)
 112 {
 113  TStopwatch timer;
 114  timer.Start();
 115
 116  if (LYZ1SUM && LYZ2SUM)  {cout<<"WARNING: you cannot run LYZ1 and LYZ2 at the same time! LYZ2 needs the output from LYZ1.  "<<endl; exit(); }
 117  if (LYZ1PROD && LYZ2PROD)  {cout<<"WARNING: you cannot run LYZ1 and LYZ2 at the same time! LYZ2 needs the output from LYZ1.  "<<endl; exit(); }
 118  if (LYZ2SUM && LYZEP) {cout<<"WARNING: you cannot run LYZ2 and LYZEP at the same time! LYZEP needs the output from LYZ2."<<endl; exit(); }
 119  if (LYZ1SUM && LYZEP) {cout<<"WARNING: you cannot run LYZ1 and LYZEP at the same time! LYZEP needs the output from LYZ2."<<endl; exit(); }
 120
 121  if(!uniformAcceptance && phimin1 > phimax1)
 122  {
 123   cout<<"WARNING: you must have phimin1 < phimax1 !!!!"<<endl;
 124   break;
 125  }
 126
 127  if (!uniformAcceptance && !((phimin2 == 0.) && (phimax2 == 0.) && (p2 == 0.)) && (phimin2 < phimax1 || phimin2 > phimax2))
 128  {
 129   cout<<"WARNING: you must have phimin2 > phimax1 and phimin2 < phimax2 !!!!"<<endl;
 130   break;
 131  }
 132
 133  if((phimin1 < 0 || phimin1 > 360) || (phimax1 < 0 || phimax1 > 360) ||
 134     (phimin2 < 0 || phimin2 > 360) || (phimax2 < 0 || phimax2 > 360) )
 135  {
 136   cout<<"WARNING: you must take azimuthal angles from interval [0,360] !!!!"<<endl;
 137   break;
 138  }
 139
 140  if((p1 < 0 || p1 > 1) || (p2 < 0 || p2 > 1))
 141  {
 142   cout<<"WARNING: you must take p1 and p2 from interval [0,1] !!!!"<<endl;
 143   break;
 144  }
 145
 146  cout<<endl;
 147  cout<<endl;
 148  cout<<"      ---- ARE YOU READY TO FLY ? ----      "<<endl;
 149  cout<<endl;
 150
 151  cout<<endl;
 152  cout<<" ---- BEGIN FLOW ANALYSIS 'ON THE FLY' ---- "<<endl;
 153  cout<<endl;
 154  cout<<endl;
 155
 156  LoadLibraries(mode);
 157
 158  // Initialize the seed for random generator
 159  UInt_t sseed = 0;
 160
 161  if(bSameSeed)
 162  {
 163   sseed = 44; // the default constant value for seed for random generators
 164  }
 165
 166  if(!bSameSeed)
 167  {
 168   TTimeStamp dt;
 169   sseed = dt.GetNanoSec()/1000;
 170  }
 171
 172  cout<<endl;
 173  cout<<"Seed for the random generators is "<<sseed<<endl;
 174  cout<<endl;
 175
 176  //---------------------------------------------------------------------------------------
 177  // If the weights are used:
 178  TFile *fileWithWeights = NULL;
 179  TList *listWithWeights = NULL;
 180
 181  if(usePhiWeights||usePtWeights||useEtaWeights) {
 182    fileWithWeights = TFile::Open("weights.root","READ");
 183    if(fileWithWeights) {
 184      listWithWeights = (TList*)fileWithWeights->Get("weights");
 185    }
 186    else
 187      {cout << " WARNING: the file <weights.root> with weights from the previous run was not found."<<endl;
 188       break;
 189      }
 190  }
 191
 192  //---------------------------------------------------------------------------------------
 193  // Initialize the flowevent maker
 194  AliFlowEventSimpleMakerOnTheFly* eventMakerOnTheFly = new AliFlowEventSimpleMakerOnTheFly(sseed);
 195  eventMakerOnTheFly->Init();
 196
 197  //---------------------------------------------------------------------------------------
 198  // Initialize all the flow methods for default analysis:
 199  AliFlowAnalysisWithQCumulants *qc = NULL;
 200  AliFlowAnalysisWithCumulants *gfc = NULL;
 201  AliFlowAnalysisWithFittingQDistribution *fqd = NULL;
 202  AliFlowAnalysisWithLeeYangZeros *lyz1sum  = NULL;
 203  AliFlowAnalysisWithLeeYangZeros *lyz1prod = NULL;
 204  AliFlowAnalysisWithLeeYangZeros *lyz2sum  = NULL;
 205  AliFlowAnalysisWithLeeYangZeros *lyz2prod = NULL;
 206  AliFlowAnalysisWithLYZEventPlane *lyzep = NULL;
 207  AliFlowAnalysisWithScalarProduct *sp = NULL;
 208  AliFlowAnalysisWithMCEventPlane *mcep = NULL;
 209
 210  // MCEP = monte carlo event plane
 211  if (MCEP) {
 212    AliFlowAnalysisWithMCEventPlane *mcep = new AliFlowAnalysisWithMCEventPlane();
 213    mcep->Init();
 214  }
 215
 216  // QC = Q-cumulants
 217  if(QC) {
 218    AliFlowAnalysisWithQCumulants* qc = new AliFlowAnalysisWithQCumulants();
 219    if(listWithWeights) qc->SetWeightsList(listWithWeights);
 220    if(usePhiWeights) qc->SetUsePhiWeights(usePhiWeights);
 221    if(usePtWeights) qc->SetUsePtWeights(usePtWeights);
 222    if(useEtaWeights) qc->SetUseEtaWeights(useEtaWeights);
 223    qc->SetCalculate2DFlow(kTRUE);
 224    qc->SetEvaluateNestedLoopsForIntFlow(kFALSE);
 225    qc->SetEvaluateNestedLoopsForDiffFlow(kFALSE);
 226    qc->Init();
 227  }
 228
 229  // GFC = Generating Function Cumulants
 230  if(GFC) {
 231    AliFlowAnalysisWithCumulants* gfc = new AliFlowAnalysisWithCumulants();
 232    if(listWithWeights) gfc->SetWeightsList(listWithWeights);
 233    if(usePhiWeights) gfc->SetUsePhiWeights(usePhiWeights);
 234    if(usePtWeights) gfc->SetUsePtWeights(usePtWeights);
 235    if(useEtaWeights) gfc->SetUseEtaWeights(useEtaWeights);
 236    gfc->Init();
 237  }
 238
 239  // FQD = Fitting q-distribution
 240  if(FQD) {
 241    AliFlowAnalysisWithFittingQDistribution* fqd = new AliFlowAnalysisWithFittingQDistribution();
 242    if(listWithWeights) fqd->SetWeightsList(listWithWeights);
 243    if(usePhiWeights) fqd->SetUsePhiWeights(usePhiWeights);
 244    fqd->Init();
 245  }
 246
 247  // SP = Scalar Product
 248  if(SP) {
 249    AliFlowAnalysisWithScalarProduct* sp = new AliFlowAnalysisWithScalarProduct();
 250    sp->Init();
 251  }
 252
 253  // LYZ1 = Lee-Yang Zeroes first run
 254  if(LYZ1SUM) {
 255    AliFlowAnalysisWithLeeYangZeros* lyz1sum = new AliFlowAnalysisWithLeeYangZeros();
 256    lyz1sum->SetFirstRun(kTRUE);
 257    lyz1sum->SetUseSum(kTRUE);
 258    lyz1sum->Init();
 259  }
 260  if(LYZ1PROD) {
 261    AliFlowAnalysisWithLeeYangZeros* lyz1prod = new AliFlowAnalysisWithLeeYangZeros();
 262    lyz1prod->SetFirstRun(kTRUE);
 263    lyz1prod->SetUseSum(kFALSE);
 264    lyz1prod->Init();
 265  }
 266  // LYZ2 = Lee-Yang Zeroes second run
 267  if(LYZ2SUM) {
 268    AliFlowAnalysisWithLeeYangZeros* lyz2sum = new AliFlowAnalysisWithLeeYangZeros();
 269    // read the input file from the first run
 270    TString inputFileNameLYZ2SUM = "outputLYZ1SUManalysis.root" ;
 271    TFile* inputFileLYZ2SUM = new TFile(inputFileNameLYZ2SUM.Data(),"READ");
 272    if(!inputFileLYZ2SUM || inputFileLYZ2SUM->IsZombie()) {
 273      cerr << " ERROR: To run LYZ2SUM you need the output file from LYZ1SUM. This file is not there! Please run LYZ1SUM first." << endl ;
 274      break;
 275    }
 276    else {
 277      TList* inputListLYZ2SUM = (TList*)inputFileLYZ2SUM->Get("cobjLYZ1SUM");
 278      if (!inputListLYZ2SUM) {cout<<"Input list LYZ2SUM is NULL pointer!"<<endl; break;}
 279      else {
 280        cout<<"LYZ2SUM input file/list read..."<<endl;
 281        lyz2sum->SetFirstRunList(inputListLYZ2SUM);
 282        lyz2sum->SetFirstRun(kFALSE);
 283        lyz2sum->SetUseSum(kTRUE);
 284        lyz2sum->Init();
 285      }
 286    }
 287  }
 288  if(LYZ2PROD) {
 289    AliFlowAnalysisWithLeeYangZeros* lyz2prod = new AliFlowAnalysisWithLeeYangZeros();
 290    // read the input file from the first run
 291    TString inputFileNameLYZ2PROD = "outputLYZ1PRODanalysis.root" ;
 292    TFile* inputFileLYZ2PROD = new TFile(inputFileNameLYZ2PROD.Data(),"READ");
 293    if(!inputFileLYZ2PROD || inputFileLYZ2PROD->IsZombie()) {
 294      cerr << " ERROR: To run LYZ2PROD you need the output file from LYZ1PROD. This file is not there! Please run LYZ1PROD first." << endl ;
 295      break;
 296    }
 297    else {
 298      TList* inputListLYZ2PROD = (TList*)inputFileLYZ2PROD->Get("cobjLYZ1PROD");
 299      if (!inputListLYZ2PROD) {cout<<"Input list LYZ2PROD is NULL pointer!"<<endl; break;}
 300      else {
 301        cout<<"LYZ2PROD input file/list read..."<<endl;
 302        lyz2prod->SetFirstRunList(inputListLYZ2PROD);
 303        lyz2prod->SetFirstRun(kFALSE);
 304        lyz2prod->SetUseSum(kFALSE);
 305        lyz2prod->Init();
 306      }
 307    }
 308  }
 309
 310  // LYZEP = Lee-Yang Zeroes event plane
 311  if(LYZEP) {
 312    AliFlowLYZEventPlane* ep = new AliFlowLYZEventPlane() ;
 313    AliFlowAnalysisWithLYZEventPlane* lyzep = new AliFlowAnalysisWithLYZEventPlane();
 314    // read the input file from the second lyz run
 315    TString inputFileNameLYZEP = "outputLYZ2SUManalysis.root" ;
 316    TFile* inputFileLYZEP = new TFile(inputFileNameLYZEP.Data(),"READ");
 317    if(!inputFileLYZEP || inputFileLYZEP->IsZombie()) {
 318      cerr << " ERROR: To run LYZEP you need the output file from LYZ2SUM. This file is not there! Please run LYZ2SUM first." << endl ;
 319      break;
 320    }
 321    else {
 322      TList* inputListLYZEP = (TList*)inputFileLYZEP->Get("cobjLYZ2SUM");
 323      if (!inputListLYZEP) {cout<<"Input list LYZEP is NULL pointer!"<<endl; break;}
 324      else {
 325        cout<<"LYZEP input file/list read..."<<endl;
 326        ep   ->SetSecondRunList(inputListLYZEP);
 327        lyzep->SetSecondRunList(inputListLYZEP);
 328        ep   ->Init();
 329        lyzep->Init();
 330      }
 331    }
 332  }
 333  //---------------------------------------------------------------------------------------
 334
 335  //---------------------------------------------------------------------------------------
 336  // Initialize all the flow methods for additional analysis with different settings/aims
 337  // Label each setting/aim with different label !!!!
 338
 339  // GFC:
 340  TString gfcDefaultName = "outputGFCanalysis";
 341  // 1.) GFC analysis for elliptic flow with r0 = 1.5:
 342  AliFlowAnalysisWithCumulants *gfc_1;
 343  TString gfcAnalysisLabels_1 = "_r0_1.5"; // all histograms and output file name will have this label
 344  TString gfcOutputFileName_1;
 345  gfcOutputFileName_1 = gfcDefaultName.Data();
 346  gfcOutputFileName_1 += gfcAnalysisLabels_1.Data();
 347  gfcOutputFileName_1 += ".root";
 348  if(GFC_Additional_Analysis)
 349  {
 350   gfc_1 = new AliFlowAnalysisWithCumulants();
 351   //gfc_1->SetAnalysisLabel(gfcAnalysisLabels_1.Data());
 352   if(listWithWeights) gfc_1->SetWeightsList(listWithWeights);
 353   if(usePhiWeights) gfc_1->SetUsePhiWeights(usePhiWeights);
 354   if(usePtWeights) gfc_1->SetUsePtWeights(usePtWeights);
 355   if(useEtaWeights) gfc_1->SetUseEtaWeights(useEtaWeights);
 356   gfc_1->Init();
 357  }
 358
 359  // QC:
 360  TString qcDefaultName = "outputQCanalysis";
 361  // 1.) QC analysis for directed flow:
 362  AliFlowAnalysisWithQCumulants *qc_1;
 363  TString qcAnalysisLabels_1 = "_v1"; // all histograms and output file name will have this label
 364  TString qcOutputFileName_1;
 365  qcOutputFileName_1 = qcDefaultName.Data();
 366  qcOutputFileName_1 += qcAnalysisLabels_1.Data();
 367  qcOutputFileName_1 += ".root";
 368  if(QC_Additional_Analysis)
 369  {
 370   qc_1 = new AliFlowAnalysisWithQCumulants();
 371   //qc_1->SetAnalysisLabel(qcAnalysisLabels_1->Data());
 372   if(listWithWeights) qc_1->SetWeightsList(listWithWeights);
 373   if(usePhiWeights) qc_1->SetUsePhiWeights(usePhiWeights);
 374   qc_1->Init();
 375  }
 376
 377  // FQD:
 378  TString fqdDefaultName = "outputFQDanalysis";
 379  // 1.) FQD fitting with fixed sigma:
 380  AliFlowAnalysisWithFittingQDistribution *fqd_1;
 381  TString fqdAnalysisLabels_1 = "_fixedSigma"; // all histograms and output file name will have this label
 382  TString fqdOutputFileName_1;
 383  fqdOutputFileName_1 = fqdDefaultName.Data();
 384  fqdOutputFileName_1 += fqdAnalysisLabels_1.Data();
 385  fqdOutputFileName_1 += ".root";
 386  if(FQD_Additional_Analysis)
 387  {
 388   fqd_1 = new AliFlowAnalysisWithFittingQDistribution();
 389   //fqd_1->SetAnalysisLabel(fqdAnalysisLabels_1->Data());
 390   if(listWithWeights) fqd_1->SetWeightsList(listWithWeights);
 391   if(usePhiWeights) fqd_1->SetUsePhiWeights(usePhiWeights);
 392   fqd_1->Init();
 393  }
 394  //---------------------------------------------------------------------------------------
 395
 396  // set the global event parameters:
 397  eventMakerOnTheFly->SetNoOfLoops(iLoops);
 398
 399  if(bMultDistrOfRPsIsGauss)
 400  {
 401   eventMakerOnTheFly->SetMultDistrOfRPsIsGauss(bMultDistrOfRPsIsGauss);
 402   eventMakerOnTheFly->SetMultiplicityOfRP(iMultiplicityOfRP);
 403   eventMakerOnTheFly->SetMultiplicitySpreadOfRP(dMultiplicitySpreadOfRP);
 404  } else
 405    {
 406     eventMakerOnTheFly->SetMultDistrOfRPsIsGauss(bMultDistrOfRPsIsGauss);
 407     eventMakerOnTheFly->SetMinMultOfRP(iMinMultOfRP);
 408     eventMakerOnTheFly->SetMaxMultOfRP(iMaxMultOfRP);
 409    }
 410
 411  eventMakerOnTheFly->SetTemperatureOfRP(dTemperatureOfRP);
 412
 413  eventMakerOnTheFly->SetV1RP(dV1RP);
 414  eventMakerOnTheFly->SetV1SpreadRP(dV1SpreadRP);
 415  eventMakerOnTheFly->SetV4RP(dV4RP);
 416  eventMakerOnTheFly->SetV4SpreadRP(dV4SpreadRP);
 417
 418  // constant harmonic V2:
 419  if(bConstantHarmonics)
 420  {
 421   eventMakerOnTheFly->SetUseConstantHarmonics(bConstantHarmonics);
 422   if(bV2DistrOfRPsIsGauss)
 423   {
 424    eventMakerOnTheFly->SetV2DistrOfRPsIsGauss(bV2DistrOfRPsIsGauss);
 425    eventMakerOnTheFly->SetV2RP(dV2RP);
 426    eventMakerOnTheFly->SetV2SpreadRP(dV2SpreadRP);
 427   } else
 428     {
 429      eventMakerOnTheFly->SetV2DistrOfRPsIsGauss(bV2DistrOfRPsIsGauss);
 430      eventMakerOnTheFly->SetMinV2RP(dMinV2RP);
 431      eventMakerOnTheFly->SetMaxV2RP(dMaxV2RP);
 432     }
 433  }
 434
 435  // (pt,eta) dependent harmonic V2:
 436  if(!bConstantHarmonics)
 437  {
 438   eventMakerOnTheFly->SetUseConstantHarmonics(bConstantHarmonics);
 439   eventMakerOnTheFly->SetV2RPMax(dV2RPMax);
 440   eventMakerOnTheFly->SetPtCutOff(dPtCutOff);
 441  }
 442
 443  // non-uniform acceptance:
 444  if(!uniformAcceptance)
 445  {
 446   eventMakerOnTheFly->SetFirstSectorPhiMin(phimin1);
 447   eventMakerOnTheFly->SetFirstSectorPhiMax(phimax1);
 448   eventMakerOnTheFly->SetFirstSectorProbability(p1);
 449   eventMakerOnTheFly->SetSecondSectorPhiMin(phimin2);
 450   eventMakerOnTheFly->SetSecondSectorPhiMax(phimax2);
 451   eventMakerOnTheFly->SetSecondSectorProbability(p2);
 452  }
 453
 454  //---------------------------------------------------------------------------------------
 455  // create and analyze events 'on the fly':
 456
 457  for(Int_t i=0;i<nEvts;i++) {
 458    // creating the event with above settings:
 459    AliFlowEventSimple *event = eventMakerOnTheFly->CreateEventOnTheFly();
 460
 461    // analyzing the created event 'on the fly':
 462    // do flow analysis for various methods:
 463    if(MCEP)    mcep->Make(event);
 464    if(QC)      qc->Make(event);
 465    if(GFC)     gfc->Make(event);
 466    if(FQD)     fqd->Make(event);
 467    if(LYZ1SUM) lyz1sum->Make(event);
 468    if(LYZ1PROD)lyz1prod->Make(event);
 469    if(LYZ2SUM) lyz2sum->Make(event);
 470    if(LYZ2PROD)lyz2prod->Make(event);
 471    if(LYZEP)   lyzep->Make(event,ep);
 472    if(SP)      sp->Make(event);
 473
 474    if(GFC_Additional_Analysis)
 475    {
 476     // r0 = 1.5:
 477     gfc_1->Make(event);
 478    }
 479    if(QC_Additional_Analysis)
 480    {
 481     // v1:
 482     qc_1->Make(event);
 483    }
 484    if(FQD_Additional_Analysis)
 485    {
 486     // fixed sigma:
 487     fqd_1->Make(event);
 488    }
 489
 490    delete event;
 491  } // end of for(Int_t i=0;i<nEvts;i++)
 492  //---------------------------------------------------------------------------------------
 493
 494
 495
 496  //---------------------------------------------------------------------------------------
 497  // calculating and storing the final results of default flow analysis:
 498  if(MCEP)    {mcep->Finish();    mcep->WriteHistograms("outputMCEPanalysis.root");}
 499  if(SP)      {sp->Finish();      sp->WriteHistograms("outputSPanalysis.root");}
 500  if(QC)      {qc->Finish();      qc->WriteHistograms("outputQCanalysis.root");}
 501  if(GFC)     {gfc->Finish();     gfc->WriteHistograms("outputGFCanalysis.root");}
 502  if(FQD)     {fqd->Finish();     fqd->WriteHistograms("outputFQDanalysis.root");}
 503  if(LYZ1SUM) {lyz1sum->Finish(); lyz1sum->WriteHistograms("outputLYZ1SUManalysis.root");}
 504  if(LYZ1PROD){lyz1prod->Finish();lyz1prod->WriteHistograms("outputLYZ1PRODanalysis.root");}
 505  if(LYZ2SUM) {lyz2sum->Finish(); lyz2sum->WriteHistograms("outputLYZ2SUManalysis.root");}
 506  if(LYZ2PROD){lyz2prod->Finish();lyz2prod->WriteHistograms("outputLYZ2PRODanalysis.root");}
 507  if(LYZEP)   {lyzep->Finish();   lyzep->WriteHistograms("outputLYZEPanalysis.root");}
 508  //---------------------------------------------------------------------------------------
 509
 510  //---------------------------------------------------------------------------------------
 511  // calculating and storing the final results of flow analysis with different settings/aims:
 512  if(GFC_Additional_Analysis)
 513  {
 514   // r0 = 1.5:
 515   gfc_1->Finish();
 516   gfc_1->WriteHistograms(gfcOutputFileName_1.Data());
 517  }
 518  if(QC_Additional_Analysis)
 519  {
 520   // v1:
 521   qc_1->Finish();
 522   qc_1->WriteHistograms(qcOutputFileName_1.Data());
 523  }
 524  if(FQD_Additional_Analysis)
 525  {
 526   // fixed sigma:
 527   fqd_1->Finish();
 528   fqd_1->WriteHistograms(fqdOutputFileName_1.Data());
 529  }
 530  //---------------------------------------------------------------------------------------
 531
 532  cout<<endl;
 533  cout<<endl;
 534  cout<<" ---- LANDED SUCCESSFULLY ---- "<<endl;
 535  cout<<endl;
 536
 537  timer.Stop();
 538  cout << endl;
 539  timer.Print();
 540 }
 541
 542 void SetupPar(char* pararchivename)
 543 {
 544   //Load par files, create analysis libraries
 545   //For testing, if par file already decompressed and modified
 546   //classes then do not decompress.
 547
 548   TString cdir(Form("%s", gSystem->WorkingDirectory() )) ;
 549   TString parpar(Form("%s.par", pararchivename)) ;
 550   if ( gSystem->AccessPathName(parpar.Data()) ) {
 551     gSystem->ChangeDirectory(gSystem->Getenv("ALICE_ROOT")) ;
 552     TString processline(Form(".! make %s", parpar.Data())) ;
 553     gROOT->ProcessLine(processline.Data()) ;
 554     gSystem->ChangeDirectory(cdir) ;
 555     processline = Form(".! mv /tmp/%s .", parpar.Data()) ;
 556     gROOT->ProcessLine(processline.Data()) ;
 557   }
 558   if ( gSystem->AccessPathName(pararchivename) ) {
 559     TString processline = Form(".! tar xvzf %s",parpar.Data()) ;
 560     gROOT->ProcessLine(processline.Data());
 561   }
 562
 563   TString ocwd = gSystem->WorkingDirectory();
 564   gSystem->ChangeDirectory(pararchivename);
 565
 566   // check for BUILD.sh and execute
 567   if (!gSystem->AccessPathName("PROOF-INF/BUILD.sh")) {
 568     printf("*******************************\n");
 569     printf("*** Building PAR archive    ***\n");
 570     cout<<pararchivename<<endl;
 571     printf("*******************************\n");
 572
 573     if (gSystem->Exec("PROOF-INF/BUILD.sh")) {
 574       Error("runProcess","Cannot Build the PAR Archive! - Abort!");
 575       return -1;
 576     }
 577   }
 578   // check for SETUP.C and execute
 579   if (!gSystem->AccessPathName("PROOF-INF/SETUP.C")) {
 580     printf("*******************************\n");
 581     printf("*** Setup PAR archive       ***\n");
 582     cout<<pararchivename<<endl;
 583     printf("*******************************\n");
 584     gROOT->Macro("PROOF-INF/SETUP.C");
 585   }
 586
 587   gSystem->ChangeDirectory(ocwd.Data());
 588   printf("Current dir: %s\n", ocwd.Data());
 589 }
 590
 591 void LoadLibraries(const anaModes mode) {
 592
 593   //--------------------------------------
 594   // Load the needed libraries most of them already loaded by aliroot
 595   //--------------------------------------
 596   gSystem->Load("libTree.so");
 597   gSystem->Load("libGeom.so");
 598   gSystem->Load("libVMC.so");
 599   gSystem->Load("libXMLIO.so");
 600   gSystem->Load("libPhysics.so");
 601
 602   //----------------------------------------------------------
 603   // >>>>>>>>>>> Local mode <<<<<<<<<<<<<<
 604   //----------------------------------------------------------
 605   if (mode==mLocal) {
 606     //--------------------------------------------------------
 607     // If you want to use already compiled libraries
 608     // in the aliroot distribution
 609     //--------------------------------------------------------
 610     gSystem->Load("libSTEERBase");
 611     gSystem->Load("libESD");
 612     gSystem->Load("libAOD");
 613     gSystem->Load("libANALYSIS");
 614     gSystem->Load("libANALYSISalice");
 615     gSystem->Load("libCORRFW.so");
 616     cerr<<"libCORRFW.so loaded..."<<endl;
 617     gSystem->Load("libPWG2flowCommon.so");
 618     cerr<<"libPWG2flowCommon.so loaded..."<<endl;
 619     gSystem->Load("libPWG2flowTasks.so");
 620     cerr<<"libPWG2flowTasks.so loaded..."<<endl;
 621   }
 622
 623   else if (mode == mLocalPAR) {
 624     //--------------------------------------------------------
 625     //If you want to use root and par files from aliroot
 626     //--------------------------------------------------------
 627      //If you want to use root and par files from aliroot
 628     //--------------------------------------------------------
 629     SetupPar("STEERBase");
 630     SetupPar("ESD");
 631     SetupPar("AOD");
 632     SetupPar("ANALYSIS");
 633     SetupPar("ANALYSISalice");
 634     SetupPar("PWG2AOD");
 635     SetupPar("CORRFW");
 636     SetupPar("PWG2flowCommon");
 637     cerr<<"PWG2flowCommon.par loaded..."<<endl;
 638     SetupPar("PWG2flowTasks");
 639     cerr<<"PWG2flowTasks.par loaded..."<<endl;
 640   }
 641
 642   //---------------------------------------------------------
 643   // <<<<<<<<<< Source mode >>>>>>>>>>>>
 644   //---------------------------------------------------------
 645   else if (mode==mLocalSource) {
 646
 647     // In root inline compile
 648
 649
 650     // Constants
 651     gROOT->LoadMacro("AliFlowCommon/AliFlowCommonConstants.cxx+");
 652     gROOT->LoadMacro("AliFlowCommon/AliFlowLYZConstants.cxx+");
 653     gROOT->LoadMacro("AliFlowCommon/AliFlowCumuConstants.cxx+");
 654
 655     // Flow event
 656     gROOT->LoadMacro("AliFlowCommon/AliFlowVector.cxx+");
 657     gROOT->LoadMacro("AliFlowCommon/AliFlowTrackSimple.cxx+");
 658     gROOT->LoadMacro("AliFlowCommon/AliFlowEventSimple.cxx+");
 659
 660     // Cuts
 661     gROOT->LoadMacro("AliFlowCommon/AliFlowTrackSimpleCuts.cxx+");
 662
 663     // Output histosgrams
 664     gROOT->LoadMacro("AliFlowCommon/AliFlowCommonHist.cxx+");
 665     gROOT->LoadMacro("AliFlowCommon/AliFlowCommonHistResults.cxx+");
 666     gROOT->LoadMacro("AliFlowCommon/AliFlowLYZHist1.cxx+");
 667     gROOT->LoadMacro("AliFlowCommon/AliFlowLYZHist2.cxx+");
 668
 669     // Functions needed for various methods
 670     gROOT->LoadMacro("AliFlowCommon/AliCumulantsFunctions.cxx+");
 671     gROOT->LoadMacro("AliFlowCommon/AliFlowLYZEventPlane.cxx+");
 672
 673     // Flow Analysis code for various methods
 674     gROOT->LoadMacro("AliFlowCommon/AliFlowAnalysisWithMCEventPlane.cxx+");
 675     gROOT->LoadMacro("AliFlowCommon/AliFlowAnalysisWithScalarProduct.cxx+");
 676     gROOT->LoadMacro("AliFlowCommon/AliFlowAnalysisWithLYZEventPlane.cxx+");
 677     gROOT->LoadMacro("AliFlowCommon/AliFlowAnalysisWithLeeYangZeros.cxx+");
 678     gROOT->LoadMacro("AliFlowCommon/AliFlowAnalysisWithCumulants.cxx+");
 679     gROOT->LoadMacro("AliFlowCommon/AliFlowAnalysisWithQCumulants.cxx+");
 680     gROOT->LoadMacro("AliFlowCommon/AliFlowAnalysisWithFittingQDistribution.cxx+");
 681
 682     // Class to fill the FlowEvent on the fly (generate Monte Carlo events)
 683     gROOT->LoadMacro("AliFlowCommon/AliFlowEventSimpleMakerOnTheFly.cxx+");
 684
 685     cout << "finished loading macros!" << endl;
 686
 687   }
 688
 689 }
 690
 691