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 nEvts=440, Int_t mode=mLocal)
 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->SetHarmonic(2); // default is v2
 214    mcep->Init();
 215  }
 216
 217  // QC = Q-cumulants
 218  if(QC) {
 219    AliFlowAnalysisWithQCumulants* qc = new AliFlowAnalysisWithQCumulants();
 220    if(listWithWeights) qc->SetWeightsList(listWithWeights);
 221    if(usePhiWeights) qc->SetUsePhiWeights(usePhiWeights);
 222    if(usePtWeights) qc->SetUsePtWeights(usePtWeights);
 223    if(useEtaWeights) qc->SetUseEtaWeights(useEtaWeights);
 224    // qc->SetHarmonic(2); // default is v2
 225    // qc->SetApplyCorrectionForNUA(kTRUE); // default
 226    // qc->SetCalculate2DFlow(kFALSE); // default
 227    qc->Init();
 228  }
 229
 230  // GFC = Generating Function Cumulants
 231  if(GFC) {
 232    AliFlowAnalysisWithCumulants* gfc = new AliFlowAnalysisWithCumulants();
 233    if(listWithWeights) gfc->SetWeightsList(listWithWeights);
 234    if(usePhiWeights) gfc->SetUsePhiWeights(usePhiWeights);
 235    if(usePtWeights) gfc->SetUsePtWeights(usePtWeights);
 236    if(useEtaWeights) gfc->SetUseEtaWeights(useEtaWeights);
 237    gfc->Init();
 238  }
 239
 240  // FQD = Fitting q-distribution
 241  if(FQD) {
 242    AliFlowAnalysisWithFittingQDistribution* fqd = new AliFlowAnalysisWithFittingQDistribution();
 243    if(listWithWeights) fqd->SetWeightsList(listWithWeights);
 244    if(usePhiWeights) fqd->SetUsePhiWeights(usePhiWeights);
 245    fqd->Init();
 246  }
 247
 248  // SP = Scalar Product
 249  if(SP) {
 250    AliFlowAnalysisWithScalarProduct* sp = new AliFlowAnalysisWithScalarProduct();
 251    sp->Init();
 252  }
 253
 254  // LYZ1 = Lee-Yang Zeroes first run
 255  if(LYZ1SUM) {
 256    AliFlowAnalysisWithLeeYangZeros* lyz1sum = new AliFlowAnalysisWithLeeYangZeros();
 257    lyz1sum->SetFirstRun(kTRUE);
 258    lyz1sum->SetUseSum(kTRUE);
 259    lyz1sum->Init();
 260  }
 261  if(LYZ1PROD) {
 262    AliFlowAnalysisWithLeeYangZeros* lyz1prod = new AliFlowAnalysisWithLeeYangZeros();
 263    lyz1prod->SetFirstRun(kTRUE);
 264    lyz1prod->SetUseSum(kFALSE);
 265    lyz1prod->Init();
 266  }
 267  // LYZ2 = Lee-Yang Zeroes second run
 268  if(LYZ2SUM) {
 269    AliFlowAnalysisWithLeeYangZeros* lyz2sum = new AliFlowAnalysisWithLeeYangZeros();
 270    // read the input file from the first run
 271    TString inputFileNameLYZ2SUM = "outputLYZ1SUManalysis.root" ;
 272    TFile* inputFileLYZ2SUM = new TFile(inputFileNameLYZ2SUM.Data(),"READ");
 273    if(!inputFileLYZ2SUM || inputFileLYZ2SUM->IsZombie()) {
 274      cerr << " ERROR: To run LYZ2SUM you need the output file from LYZ1SUM. This file is not there! Please run LYZ1SUM first." << endl ;
 275      break;
 276    }
 277    else {
 278      TList* inputListLYZ2SUM = (TList*)inputFileLYZ2SUM->Get("cobjLYZ1SUM");
 279      if (!inputListLYZ2SUM) {cout<<"Input list LYZ2SUM is NULL pointer!"<<endl; break;}
 280      else {
 281        cout<<"LYZ2SUM input file/list read..."<<endl;
 282        lyz2sum->SetFirstRunList(inputListLYZ2SUM);
 283        lyz2sum->SetFirstRun(kFALSE);
 284        lyz2sum->SetUseSum(kTRUE);
 285        lyz2sum->Init();
 286      }
 287    }
 288  }
 289  if(LYZ2PROD) {
 290    AliFlowAnalysisWithLeeYangZeros* lyz2prod = new AliFlowAnalysisWithLeeYangZeros();
 291    // read the input file from the first run
 292    TString inputFileNameLYZ2PROD = "outputLYZ1PRODanalysis.root" ;
 293    TFile* inputFileLYZ2PROD = new TFile(inputFileNameLYZ2PROD.Data(),"READ");
 294    if(!inputFileLYZ2PROD || inputFileLYZ2PROD->IsZombie()) {
 295      cerr << " ERROR: To run LYZ2PROD you need the output file from LYZ1PROD. This file is not there! Please run LYZ1PROD first." << endl ;
 296      break;
 297    }
 298    else {
 299      TList* inputListLYZ2PROD = (TList*)inputFileLYZ2PROD->Get("cobjLYZ1PROD");
 300      if (!inputListLYZ2PROD) {cout<<"Input list LYZ2PROD is NULL pointer!"<<endl; break;}
 301      else {
 302        cout<<"LYZ2PROD input file/list read..."<<endl;
 303        lyz2prod->SetFirstRunList(inputListLYZ2PROD);
 304        lyz2prod->SetFirstRun(kFALSE);
 305        lyz2prod->SetUseSum(kFALSE);
 306        lyz2prod->Init();
 307      }
 308    }
 309  }
 310
 311  // LYZEP = Lee-Yang Zeroes event plane
 312  if(LYZEP) {
 313    AliFlowLYZEventPlane* ep = new AliFlowLYZEventPlane() ;
 314    AliFlowAnalysisWithLYZEventPlane* lyzep = new AliFlowAnalysisWithLYZEventPlane();
 315    // read the input file from the second lyz run
 316    TString inputFileNameLYZEP = "outputLYZ2SUManalysis.root" ;
 317    TFile* inputFileLYZEP = new TFile(inputFileNameLYZEP.Data(),"READ");
 318    if(!inputFileLYZEP || inputFileLYZEP->IsZombie()) {
 319      cerr << " ERROR: To run LYZEP you need the output file from LYZ2SUM. This file is not there! Please run LYZ2SUM first." << endl ;
 320      break;
 321    }
 322    else {
 323      TList* inputListLYZEP = (TList*)inputFileLYZEP->Get("cobjLYZ2SUM");
 324      if (!inputListLYZEP) {cout<<"Input list LYZEP is NULL pointer!"<<endl; break;}
 325      else {
 326        cout<<"LYZEP input file/list read..."<<endl;
 327        ep   ->SetSecondRunList(inputListLYZEP);
 328        lyzep->SetSecondRunList(inputListLYZEP);
 329        ep   ->Init();
 330        lyzep->Init();
 331      }
 332    }
 333  }
 334  //---------------------------------------------------------------------------------------
 335
 336  //---------------------------------------------------------------------------------------
 337  // Initialize all the flow methods for additional analysis with different settings/aims
 338  // Label each setting/aim with different label !!!!
 339
 340  // GFC:
 341  TString gfcDefaultName = "outputGFCanalysis";
 342  // 1.) GFC analysis for elliptic flow with r0 = 1.5:
 343  AliFlowAnalysisWithCumulants *gfc_1;
 344  TString gfcAnalysisLabels_1 = "_r0_1.5"; // all histograms and output file name will have this label
 345  TString gfcOutputFileName_1;
 346  gfcOutputFileName_1 = gfcDefaultName.Data();
 347  gfcOutputFileName_1 += gfcAnalysisLabels_1.Data();
 348  gfcOutputFileName_1 += ".root";
 349  if(GFC_Additional_Analysis)
 350  {
 351   gfc_1 = new AliFlowAnalysisWithCumulants();
 352   //gfc_1->SetAnalysisLabel(gfcAnalysisLabels_1.Data());
 353   if(listWithWeights) gfc_1->SetWeightsList(listWithWeights);
 354   if(usePhiWeights) gfc_1->SetUsePhiWeights(usePhiWeights);
 355   if(usePtWeights) gfc_1->SetUsePtWeights(usePtWeights);
 356   if(useEtaWeights) gfc_1->SetUseEtaWeights(useEtaWeights);
 357   gfc_1->Init();
 358  }
 359
 360  // QC:
 361  TString qcDefaultName = "outputQCanalysis";
 362  // 1.) QC analysis for directed flow:
 363  AliFlowAnalysisWithQCumulants *qc_1;
 364  TString qcAnalysisLabels_1 = "_v1"; // all histograms and output file name will have this label
 365  TString qcOutputFileName_1;
 366  qcOutputFileName_1 = qcDefaultName.Data();
 367  qcOutputFileName_1 += qcAnalysisLabels_1.Data();
 368  qcOutputFileName_1 += ".root";
 369  if(QC_Additional_Analysis)
 370  {
 371   qc_1 = new AliFlowAnalysisWithQCumulants();
 372   //qc_1->SetAnalysisLabel(qcAnalysisLabels_1->Data());
 373   if(listWithWeights) qc_1->SetWeightsList(listWithWeights);
 374   if(usePhiWeights) qc_1->SetUsePhiWeights(usePhiWeights);
 375   qc_1->Init();
 376  }
 377
 378  // FQD:
 379  TString fqdDefaultName = "outputFQDanalysis";
 380  // 1.) FQD fitting with fixed sigma:
 381  AliFlowAnalysisWithFittingQDistribution *fqd_1;
 382  TString fqdAnalysisLabels_1 = "_fixedSigma"; // all histograms and output file name will have this label
 383  TString fqdOutputFileName_1;
 384  fqdOutputFileName_1 = fqdDefaultName.Data();
 385  fqdOutputFileName_1 += fqdAnalysisLabels_1.Data();
 386  fqdOutputFileName_1 += ".root";
 387  if(FQD_Additional_Analysis)
 388  {
 389   fqd_1 = new AliFlowAnalysisWithFittingQDistribution();
 390   //fqd_1->SetAnalysisLabel(fqdAnalysisLabels_1->Data());
 391   if(listWithWeights) fqd_1->SetWeightsList(listWithWeights);
 392   if(usePhiWeights) fqd_1->SetUsePhiWeights(usePhiWeights);
 393   fqd_1->Init();
 394  }
 395  //---------------------------------------------------------------------------------------
 396
 397  // set the global event parameters:
 398  eventMakerOnTheFly->SetNoOfLoops(iLoops);
 399
 400  if(bMultDistrOfRPsIsGauss)
 401  {
 402   eventMakerOnTheFly->SetMultDistrOfRPsIsGauss(bMultDistrOfRPsIsGauss);
 403   eventMakerOnTheFly->SetMultiplicityOfRP(iMultiplicityOfRP);
 404   eventMakerOnTheFly->SetMultiplicitySpreadOfRP(dMultiplicitySpreadOfRP);
 405  } else
 406    {
 407     eventMakerOnTheFly->SetMultDistrOfRPsIsGauss(bMultDistrOfRPsIsGauss);
 408     eventMakerOnTheFly->SetMinMultOfRP(iMinMultOfRP);
 409     eventMakerOnTheFly->SetMaxMultOfRP(iMaxMultOfRP);
 410    }
 411
 412  eventMakerOnTheFly->SetTemperatureOfRP(dTemperatureOfRP);
 413
 414  eventMakerOnTheFly->SetV1RP(dV1RP);
 415  eventMakerOnTheFly->SetV1SpreadRP(dV1SpreadRP);
 416  eventMakerOnTheFly->SetV4RP(dV4RP);
 417  eventMakerOnTheFly->SetV4SpreadRP(dV4SpreadRP);
 418
 419  // constant harmonic V2:
 420  if(bConstantHarmonics)
 421  {
 422   eventMakerOnTheFly->SetUseConstantHarmonics(bConstantHarmonics);
 423   if(bV2DistrOfRPsIsGauss)
 424   {
 425    eventMakerOnTheFly->SetV2DistrOfRPsIsGauss(bV2DistrOfRPsIsGauss);
 426    eventMakerOnTheFly->SetV2RP(dV2RP);
 427    eventMakerOnTheFly->SetV2SpreadRP(dV2SpreadRP);
 428   } else
 429     {
 430      eventMakerOnTheFly->SetV2DistrOfRPsIsGauss(bV2DistrOfRPsIsGauss);
 431      eventMakerOnTheFly->SetMinV2RP(dMinV2RP);
 432      eventMakerOnTheFly->SetMaxV2RP(dMaxV2RP);
 433     }
 434  }
 435
 436  // (pt,eta) dependent harmonic V2:
 437  if(!bConstantHarmonics)
 438  {
 439   eventMakerOnTheFly->SetUseConstantHarmonics(bConstantHarmonics);
 440   eventMakerOnTheFly->SetV2RPMax(dV2RPMax);
 441   eventMakerOnTheFly->SetPtCutOff(dPtCutOff);
 442  }
 443
 444  // non-uniform acceptance:
 445  if(!uniformAcceptance)
 446  {
 447   eventMakerOnTheFly->SetFirstSectorPhiMin(phimin1);
 448   eventMakerOnTheFly->SetFirstSectorPhiMax(phimax1);
 449   eventMakerOnTheFly->SetFirstSectorProbability(p1);
 450   eventMakerOnTheFly->SetSecondSectorPhiMin(phimin2);
 451   eventMakerOnTheFly->SetSecondSectorPhiMax(phimax2);
 452   eventMakerOnTheFly->SetSecondSectorProbability(p2);
 453  }
 454
 455  //---------------------------------------------------------------------------------------
 456  // create and analyze events 'on the fly':
 457
 458  for(Int_t i=0;i<nEvts;i++) {
 459    // creating the event with above settings:
 460    AliFlowEventSimple *event = eventMakerOnTheFly->CreateEventOnTheFly();
 461
 462    // analyzing the created event 'on the fly':
 463    // do flow analysis for various methods:
 464    if(MCEP)    mcep->Make(event);
 465    if(QC)      qc->Make(event);
 466    if(GFC)     gfc->Make(event);
 467    if(FQD)     fqd->Make(event);
 468    if(LYZ1SUM) lyz1sum->Make(event);
 469    if(LYZ1PROD)lyz1prod->Make(event);
 470    if(LYZ2SUM) lyz2sum->Make(event);
 471    if(LYZ2PROD)lyz2prod->Make(event);
 472    if(LYZEP)   lyzep->Make(event,ep);
 473    if(SP)      sp->Make(event);
 474
 475    if(GFC_Additional_Analysis)
 476    {
 477     // r0 = 1.5:
 478     gfc_1->Make(event);
 479    }
 480    if(QC_Additional_Analysis)
 481    {
 482     // v1:
 483     qc_1->Make(event);
 484    }
 485    if(FQD_Additional_Analysis)
 486    {
 487     // fixed sigma:
 488     fqd_1->Make(event);
 489    }
 490
 491    delete event;
 492  } // end of for(Int_t i=0;i<nEvts;i++)
 493  //---------------------------------------------------------------------------------------
 494
 495
 496
 497  //---------------------------------------------------------------------------------------
 498  // calculating and storing the final results of default flow analysis:
 499  if(MCEP)    {mcep->Finish();    mcep->WriteHistograms("outputMCEPanalysis.root");}
 500  if(SP)      {sp->Finish();      sp->WriteHistograms("outputSPanalysis.root");}
 501  if(QC)      {qc->Finish();      qc->WriteHistograms("outputQCanalysis.root");}
 502  if(GFC)     {gfc->Finish();     gfc->WriteHistograms("outputGFCanalysis.root");}
 503  if(FQD)     {fqd->Finish();     fqd->WriteHistograms("outputFQDanalysis.root");}
 504  if(LYZ1SUM) {lyz1sum->Finish(); lyz1sum->WriteHistograms("outputLYZ1SUManalysis.root");}
 505  if(LYZ1PROD){lyz1prod->Finish();lyz1prod->WriteHistograms("outputLYZ1PRODanalysis.root");}
 506  if(LYZ2SUM) {lyz2sum->Finish(); lyz2sum->WriteHistograms("outputLYZ2SUManalysis.root");}
 507  if(LYZ2PROD){lyz2prod->Finish();lyz2prod->WriteHistograms("outputLYZ2PRODanalysis.root");}
 508  if(LYZEP)   {lyzep->Finish();   lyzep->WriteHistograms("outputLYZEPanalysis.root");}
 509  //---------------------------------------------------------------------------------------
 510
 511  //---------------------------------------------------------------------------------------
 512  // calculating and storing the final results of flow analysis with different settings/aims:
 513  if(GFC_Additional_Analysis)
 514  {
 515   // r0 = 1.5:
 516   gfc_1->Finish();
 517   gfc_1->WriteHistograms(gfcOutputFileName_1.Data());
 518  }
 519  if(QC_Additional_Analysis)
 520  {
 521   // v1:
 522   qc_1->Finish();
 523   qc_1->WriteHistograms(qcOutputFileName_1.Data());
 524  }
 525  if(FQD_Additional_Analysis)
 526  {
 527   // fixed sigma:
 528   fqd_1->Finish();
 529   fqd_1->WriteHistograms(fqdOutputFileName_1.Data());
 530  }
 531  //---------------------------------------------------------------------------------------
 532
 533  cout<<endl;
 534  cout<<endl;
 535  cout<<" ---- LANDED SUCCESSFULLY ---- "<<endl;
 536  cout<<endl;
 537
 538  timer.Stop();
 539  cout << endl;
 540  timer.Print();
 541 }
 542
 543 void SetupPar(char* pararchivename)
 544 {
 545   //Load par files, create analysis libraries
 546   //For testing, if par file already decompressed and modified
 547   //classes then do not decompress.
 548
 549   TString cdir(Form("%s", gSystem->WorkingDirectory() )) ;
 550   TString parpar(Form("%s.par", pararchivename)) ;
 551   if ( gSystem->AccessPathName(parpar.Data()) ) {
 552     gSystem->ChangeDirectory(gSystem->Getenv("ALICE_ROOT")) ;
 553     TString processline(Form(".! make %s", parpar.Data())) ;
 554     gROOT->ProcessLine(processline.Data()) ;
 555     gSystem->ChangeDirectory(cdir) ;
 556     processline = Form(".! mv /tmp/%s .", parpar.Data()) ;
 557     gROOT->ProcessLine(processline.Data()) ;
 558   }
 559   if ( gSystem->AccessPathName(pararchivename) ) {
 560     TString processline = Form(".! tar xvzf %s",parpar.Data()) ;
 561     gROOT->ProcessLine(processline.Data());
 562   }
 563
 564   TString ocwd = gSystem->WorkingDirectory();
 565   gSystem->ChangeDirectory(pararchivename);
 566
 567   // check for BUILD.sh and execute
 568   if (!gSystem->AccessPathName("PROOF-INF/BUILD.sh")) {
 569     printf("*******************************\n");
 570     printf("*** Building PAR archive    ***\n");
 571     cout<<pararchivename<<endl;
 572     printf("*******************************\n");
 573
 574     if (gSystem->Exec("PROOF-INF/BUILD.sh")) {
 575       Error("runProcess","Cannot Build the PAR Archive! - Abort!");
 576       return -1;
 577     }
 578   }
 579   // check for SETUP.C and execute
 580   if (!gSystem->AccessPathName("PROOF-INF/SETUP.C")) {
 581     printf("*******************************\n");
 582     printf("*** Setup PAR archive       ***\n");
 583     cout<<pararchivename<<endl;
 584     printf("*******************************\n");
 585     gROOT->Macro("PROOF-INF/SETUP.C");
 586   }
 587
 588   gSystem->ChangeDirectory(ocwd.Data());
 589   printf("Current dir: %s\n", ocwd.Data());
 590 }
 591
 592 void LoadLibraries(const anaModes mode) {
 593
 594   //--------------------------------------
 595   // Load the needed libraries most of them already loaded by aliroot
 596   //--------------------------------------
 597   gSystem->Load("libTree");
 598   gSystem->Load("libGeom");
 599   gSystem->Load("libVMC");
 600   gSystem->Load("libXMLIO");
 601   gSystem->Load("libPhysics");
 602
 603   //----------------------------------------------------------
 604   // >>>>>>>>>>> Local mode <<<<<<<<<<<<<<
 605   //----------------------------------------------------------
 606   if (mode==mLocal) {
 607     //--------------------------------------------------------
 608     // If you want to use already compiled libraries
 609     // in the aliroot distribution
 610     //--------------------------------------------------------
 611     gSystem->Load("libSTEERBase");
 612     gSystem->Load("libESD");
 613     gSystem->Load("libAOD");
 614     gSystem->Load("libANALYSIS");
 615     gSystem->Load("libANALYSISalice");
 616     gSystem->Load("libCORRFW");
 617     cerr<<"libCORRFW loaded..."<<endl;
 618     gSystem->Load("libPWG2flowCommon");
 619     cerr<<"libPWG2flowCommon loaded..."<<endl;
 620     gSystem->Load("libPWG2flowTasks");
 621     cerr<<"libPWG2flowTasks loaded..."<<endl;
 622   }
 623
 624   else if (mode == mLocalPAR) {
 625     //--------------------------------------------------------
 626     //If you want to use root and par files from aliroot
 627     //--------------------------------------------------------
 628      //If you want to use root and par files from aliroot
 629     //--------------------------------------------------------
 630     SetupPar("STEERBase");
 631     SetupPar("ESD");
 632     SetupPar("AOD");
 633     SetupPar("ANALYSIS");
 634     SetupPar("ANALYSISalice");
 635     SetupPar("PWG2AOD");
 636     SetupPar("CORRFW");
 637     SetupPar("PWG2flowCommon");
 638     cerr<<"PWG2flowCommon.par loaded..."<<endl;
 639     SetupPar("PWG2flowTasks");
 640     cerr<<"PWG2flowTasks.par loaded..."<<endl;
 641   }
 642
 643   //---------------------------------------------------------
 644   // <<<<<<<<<< Source mode >>>>>>>>>>>>
 645   //---------------------------------------------------------
 646   else if (mode==mLocalSource) {
 647
 648     // In root inline compile
 649
 650
 651     // Constants
 652     gROOT->LoadMacro("AliFlowCommon/AliFlowCommonConstants.cxx+");
 653     gROOT->LoadMacro("AliFlowCommon/AliFlowLYZConstants.cxx+");
 654     gROOT->LoadMacro("AliFlowCommon/AliFlowCumuConstants.cxx+");
 655
 656     // Flow event
 657     gROOT->LoadMacro("AliFlowCommon/AliFlowVector.cxx+");
 658     gROOT->LoadMacro("AliFlowCommon/AliFlowTrackSimple.cxx+");
 659     gROOT->LoadMacro("AliFlowCommon/AliFlowEventSimple.cxx+");
 660
 661     // Cuts
 662     gROOT->LoadMacro("AliFlowCommon/AliFlowTrackSimpleCuts.cxx+");
 663
 664     // Output histosgrams
 665     gROOT->LoadMacro("AliFlowCommon/AliFlowCommonHist.cxx+");
 666     gROOT->LoadMacro("AliFlowCommon/AliFlowCommonHistResults.cxx+");
 667     gROOT->LoadMacro("AliFlowCommon/AliFlowLYZHist1.cxx+");
 668     gROOT->LoadMacro("AliFlowCommon/AliFlowLYZHist2.cxx+");
 669
 670     // Functions needed for various methods
 671     gROOT->LoadMacro("AliFlowCommon/AliCumulantsFunctions.cxx+");
 672     gROOT->LoadMacro("AliFlowCommon/AliFlowLYZEventPlane.cxx+");
 673
 674     // Flow Analysis code for various methods
 675     gROOT->LoadMacro("AliFlowCommon/AliFlowAnalysisWithMCEventPlane.cxx+");
 676     gROOT->LoadMacro("AliFlowCommon/AliFlowAnalysisWithScalarProduct.cxx+");
 677     gROOT->LoadMacro("AliFlowCommon/AliFlowAnalysisWithLYZEventPlane.cxx+");
 678     gROOT->LoadMacro("AliFlowCommon/AliFlowAnalysisWithLeeYangZeros.cxx+");
 679     gROOT->LoadMacro("AliFlowCommon/AliFlowAnalysisWithCumulants.cxx+");
 680     gROOT->LoadMacro("AliFlowCommon/AliFlowAnalysisWithQCumulants.cxx+");
 681     gROOT->LoadMacro("AliFlowCommon/AliFlowAnalysisWithFittingQDistribution.cxx+");
 682
 683     // Class to fill the FlowEvent on the fly (generate Monte Carlo events)
 684     gROOT->LoadMacro("AliFlowCommon/AliFlowEventSimpleMakerOnTheFly.cxx+");
 685
 686     cout << "finished loading macros!" << endl;
 687
 688   }
 689
 690 }
 691
 692