PWG2/FLOW/macros/runAliAnalysisTaskFlow.C

   1 /////////////////////////////////////////////////////////////////////////////////
   2 //
   3 // HOW TO USE THIS MACRO:
   4 //
   5 // With this macro several flow analysis can be run.
   6 // SP    = Scalar Product                (for PbPb or pp)
   7 // LYZ1  = Lee Yang Zeroes first run     (for PbPb)
   8 // LYZ2  = Lee Yang Zeroes second run    (for PbPb)
   9 // LYZEP = Lee Yang Zeroes Event Plane   (for PbPb)
  10 // GFC   = Cumulants                     (for PbPb)
  11 // QC    = Q-cumulants                   (for PbPb or pp)
  12 // FQD   = Fitting q-distribution        (for PbPb)
  13 // MCEP  = Flow calculated from the real MC event plane (for PbPb only)
  14 //
  15 // The LYZ analysis should be done in the following order;
  16 // LYZ1 -> LYZ2 -> LYZEP,
  17 // because LYZ2 depends on the outputfile of LYZ1 and LYZEP on the outputfile
  18 // of LYZ2.
  19 //
  20 // The MCEP method is a reference method.
  21 // It can only be run when MC information (kinematics.root & galice.root file)
  22 // is available in which the reaction plane is stored.
  23 //
  24 // One can run on ESD, AOD or MC.
  25 // Additional options are ESDMC0, ESDMC1. In these options the ESD and MC
  26 // information is combined. Tracks are selected in the ESD, the PID information
  27 // is taken from the MC (perfect PID). For ESDMC0 the track kinematics is taken
  28 // from the ESD and for ESDMC1 it is taken from the MC information.
  29 //
  30 // the macro can be used to run local in aliroot or root, on the grid and on caf
  31 ///////////////////////////////////////////////////////////////////////////////////
  32
  33 enum anaModes {mLocal,mLocalPAR,mPROOF,mGRID};
  34 //mLocal: Analyze locally files in your computer using aliroot
  35 //mLocalPAR: Analyze locally files in your computer using root + PAR files
  36 //mPROOF: Analyze CAF files with PROOF
  37
  38 // RUN SETTINGS
  39
  40 // Flow analysis method can be:(set to kTRUE or kFALSE)
  41 Bool_t SP    = kTRUE;
  42 Bool_t LYZ1  = kTRUE;
  43 Bool_t LYZ2  = kFALSE;
  44 Bool_t LYZEP = kFALSE;
  45 Bool_t GFC   = kTRUE;
  46 Bool_t QC    = kTRUE;
  47 Bool_t FQD   = kTRUE;
  48 Bool_t MCEP  = kTRUE;
  49
  50 // Analysis type can be ESD, AOD, MC, ESDMC0, ESDMC1
  51 const TString type = "ESD";
  52
  53 // Boolean to fill/not fill the QA histograms
  54 Bool_t QA = kFALSE;
  55
  56 // Weights
  57 //Use weights for Q vector
  58 Bool_t usePhiWeights = kFALSE; //Phi
  59 Bool_t usePtWeights  = kFALSE; //v'(pt)
  60 Bool_t useEtaWeights = kFALSE; //v'(eta)
  61 Bool_t useWeights = usePhiWeights||usePtWeights||useEtaWeights;
  62
  63
  64 // SETTING THE CUTS
  65
  66 // For integrated flow
  67 const Double_t ptmin1 = 0.0;
  68 const Double_t ptmax1 = 10.0;
  69 const Double_t ymin1  = -1.;
  70 const Double_t ymax1  = 1.;
  71 const Int_t mintrackrefsTPC1 = 2;
  72 const Int_t mintrackrefsITS1 = 3;
  73 const Int_t charge1 = 1;
  74 Bool_t UsePIDIntegratedFlow = kFALSE;
  75 const Int_t PDG1 = 211;
  76 const Int_t minclustersTPC1 = 50;
  77 const Int_t maxnsigmatovertex1 = 3;
  78
  79 // For differential flow
  80 const Double_t ptmin2 = 0.0;
  81 const Double_t ptmax2 = 10.0;
  82 const Double_t ymin2  = -1.;
  83 const Double_t ymax2  = 1.;
  84 const Int_t mintrackrefsTPC2 = 2;
  85 const Int_t mintrackrefsITS2 = 3;
  86 const Int_t charge2 = 1;
  87 Bool_t UsePIDDifferentialFlow = kFALSE;
  88 const Int_t PDG2 = 211;
  89 const Int_t minclustersTPC2 = 50;
  90 const Int_t maxnsigmatovertex2 = 3;
  91
  92 // ESD data on PROOF cluster
  93 //void runAliAnalysisTaskFlow(Int_t mode=mPROOF, const Char_t* data="/PWG2/akisiel/Therminator_c2030", Int_t nRuns=-1, Int_t offset=0)
  94
  95 // AOD data on PROOF cluster
  96 //void runAliAnalysisTaskFlow(Int_t mode=mPROOF, const Char_t* data="/PWG2/nkolk/myDataSet", Int_t nRuns=-1, Int_t offset=0)
  97 //void runAliAnalysisTaskFlow(Int_t mode=mPROOF, const Char_t* data="/PWG2/akisiel/Therminator_midcentral_AOD", Int_t nRuns=44, Int_t offset=0)
  98
  99 // Data at Nikhef
 100 //void runAliAnalysisTaskFlow(Int_t mode=mLocal, Int_t nRuns = 4, const Char_t* dataDir="/data/alice2/kolk/Therminator_midcentral", Int_t offset = 0)
 101 //void runAliAnalysisTaskFlow(Int_t mode=mLocalPAR, Int_t nRuns = 55, const Char_t* dataDir="/data/alice2/kolk/Therminator_midcentral", Int_t offset = 0)
 102
 103 // Data on my mac
 104 void runAliAnalysisTaskFlow(Int_t mode=mLocal, Int_t nRuns = -1, const Char_t* dataDir="/Users/snelling/alice_data/Therminator_midcentral", Int_t offset = 0)
 105 //void runAliAnalysisTaskFlow(Int_t mode=mLocalPAR, Int_t nRuns = 55, const Char_t* dataDir="/Users/snelling/alice_data/Therminator_midcentral", Int_t offset = 0)
 106
 107 {
 108
 109  TStopwatch timer;
 110  timer.Start();
 111
 112  if (LYZ1 && LYZ2) {cout<<"WARNING: you cannot run LYZ1 and LYZ2 at the same time! LYZ2 needs the output from LYZ1."<<endl; exit(); }
 113  if (LYZ2 && LYZEP) {cout<<"WARNING: you cannot run LYZ2 and LYZEP at the same time! LYZEP needs the output from LYZ2."<<endl; exit(); }
 114  if (LYZ1 && LYZEP) {cout<<"WARNING: you cannot run LYZ1 and LYZEP at the same time! LYZEP needs the output from LYZ2."<<endl; exit(); }
 115
 116  LoadLibraries(mode);
 117
 118 if (mode==mLocal || mode == mLocalPAR || mode == mGRID) {
 119   if (type!="AOD") { TChain* chain = CreateESDChain(dataDir, nRuns, offset);}
 120   else { TChain* chain = CreateAODChain(dataDir, nRuns, offset);}
 121 }
 122
 123  //Create cuts using correction framework
 124
 125  //Set TList for the QA histograms
 126  if (QA) {
 127    TList* qaIntFE = new TList(); TList* qaDiffFE = new TList();
 128  }
 129
 130  //############# cuts on MC
 131  AliCFTrackKineCuts* mcKineCuts1 = new AliCFTrackKineCuts("mcKineCuts1","MC-level kinematic cuts");
 132  mcKineCuts1->SetPtRange(ptmin1,ptmax1);
 133  mcKineCuts1->SetRapidityRange(ymin1,ymax1);
 134  mcKineCuts1->SetChargeMC(charge1);
 135  if (QA) {
 136    mcKineCuts1->SetQAOn(qaIntFE);
 137  }
 138
 139  AliCFTrackKineCuts* mcKineCuts2 = new AliCFTrackKineCuts("mcKineCuts2","MC-level kinematic cuts");
 140  mcKineCuts2->SetPtRange(ptmin2,ptmax2);
 141  mcKineCuts2->SetRapidityRange(ymin2,ymax2);
 142  mcKineCuts2->SetChargeMC(charge2);
 143  if (QA) {
 144    mcKineCuts2->SetQAOn(qaDiffFE);
 145  }
 146
 147  AliCFParticleGenCuts* mcGenCuts1 = new AliCFParticleGenCuts("mcGenCuts1","MC particle generation cuts for integrated flow");
 148  mcGenCuts1->SetRequireIsPrimary();
 149  if (UsePIDIntegratedFlow) {mcGenCuts1->SetRequirePdgCode(PDG1);}
 150  if (QA) {
 151    mcGenCuts1->SetQAOn(qaIntFE);
 152  }
 153
 154  AliCFParticleGenCuts* mcGenCuts2 = new AliCFParticleGenCuts("mcGenCuts2","MC particle generation cuts for differential flow");
 155  mcGenCuts2->SetRequireIsPrimary();
 156  if (UsePIDDifferentialFlow) {mcGenCuts2->SetRequirePdgCode(PDG2);}
 157  if (QA) {
 158    mcGenCuts2->SetQAOn(qaDiffFE);
 159  }
 160
 161  //############# Acceptance Cuts
 162  AliCFAcceptanceCuts *mcAccCuts1 = new AliCFAcceptanceCuts("mcAccCuts1","MC acceptance cuts");
 163  mcAccCuts1->SetMinNHitITS(mintrackrefsITS1);
 164  mcAccCuts1->SetMinNHitTPC(mintrackrefsTPC1);
 165  if (QA) {
 166    mcAccCuts1->SetQAOn(qaIntFE);
 167  }
 168
 169  AliCFAcceptanceCuts *mcAccCuts2 = new AliCFAcceptanceCuts("mcAccCuts2","MC acceptance cuts");
 170  mcAccCuts2->SetMinNHitITS(mintrackrefsITS2);
 171  mcAccCuts2->SetMinNHitTPC(mintrackrefsTPC2);
 172  if (QA) {
 173    mcAccCuts2->SetQAOn(qaDiffFE);
 174  }
 175
 176  //############# Rec-Level kinematic cuts
 177  AliCFTrackKineCuts *recKineCuts1 = new AliCFTrackKineCuts("recKineCuts1","rec-level kine cuts");
 178  recKineCuts1->SetPtRange(ptmin1,ptmax1);
 179  recKineCuts1->SetRapidityRange(ymin1,ymax1);
 180  recKineCuts1->SetChargeRec(charge1);
 181  if (QA) {
 182    recKineCuts1->SetQAOn(qaIntFE);
 183  }
 184
 185  AliCFTrackKineCuts *recKineCuts2 = new AliCFTrackKineCuts("recKineCuts2","rec-level kine cuts");
 186  recKineCuts2->SetPtRange(ptmin2,ptmax2);
 187  recKineCuts2->SetRapidityRange(ymin2,ymax2);
 188  recKineCuts2->SetChargeRec(charge2);
 189  if (QA) {
 190    recKineCuts2->SetQAOn(qaDiffFE);
 191  }
 192
 193  AliCFTrackQualityCuts *recQualityCuts1 = new AliCFTrackQualityCuts("recQualityCuts1","rec-level quality cuts");
 194  recQualityCuts1->SetMinNClusterTPC(minclustersTPC1);
 195  recQualityCuts1->SetStatus(AliESDtrack::kITSrefit);
 196  if (QA) {
 197    recQualityCuts1->SetQAOn(qaIntFE);
 198  }
 199
 200  AliCFTrackQualityCuts *recQualityCuts2 = new AliCFTrackQualityCuts("recQualityCuts2","rec-level quality cuts");
 201  recQualityCuts2->SetMinNClusterTPC(minclustersTPC2);
 202  recQualityCuts2->SetStatus(AliESDtrack::kITSrefit);
 203  if (QA) {
 204    recQualityCuts2->SetQAOn(qaDiffFE);
 205  }
 206
 207  AliCFTrackIsPrimaryCuts *recIsPrimaryCuts1 = new AliCFTrackIsPrimaryCuts("recIsPrimaryCuts1","rec-level isPrimary cuts");
 208  recIsPrimaryCuts1->SetMaxNSigmaToVertex(maxnsigmatovertex1);
 209  if (QA) {
 210    recIsPrimaryCuts1->SetQAOn(qaIntFE);
 211  }
 212
 213  AliCFTrackIsPrimaryCuts *recIsPrimaryCuts2 = new AliCFTrackIsPrimaryCuts("recIsPrimaryCuts2","rec-level isPrimary cuts");
 214  recIsPrimaryCuts2->SetMaxNSigmaToVertex(maxnsigmatovertex2);
 215  if (QA) {
 216    recIsPrimaryCuts2->SetQAOn(qaDiffFE);
 217  }
 218
 219  int n_species = AliPID::kSPECIES ;
 220  Double_t* prior = new Double_t[n_species];
 221
 222  prior[0] = 0.0244519 ;
 223  prior[1] = 0.0143988 ;
 224  prior[2] = 0.805747  ;
 225  prior[3] = 0.0928785 ;
 226  prior[4] = 0.0625243 ;
 227
 228  AliCFTrackCutPid* cutPID1 = NULL;
 229  if(UsePIDIntegratedFlow) {
 230    AliCFTrackCutPid* cutPID1 = new AliCFTrackCutPid("cutPID1","ESD_PID for integrated flow") ;
 231    cutPID1->SetPriors(prior);
 232    cutPID1->SetProbabilityCut(0.0);
 233    cutPID1->SetDetectors("TPC TOF");
 234    switch(TMath::Abs(PDG1)) {
 235    case 11   : cutPID1->SetParticleType(AliPID::kElectron, kTRUE); break;
 236    case 13   : cutPID1->SetParticleType(AliPID::kMuon    , kTRUE); break;
 237    case 211  : cutPID1->SetParticleType(AliPID::kPion    , kTRUE); break;
 238    case 321  : cutPID1->SetParticleType(AliPID::kKaon    , kTRUE); break;
 239    case 2212 : cutPID1->SetParticleType(AliPID::kProton  , kTRUE); break;
 240    default   : printf("UNDEFINED PID\n"); break;
 241    }
 242    if (QA) {
 243      cutPID1->SetQAOn(qaIntFE);
 244    }
 245  }
 246
 247  AliCFTrackCutPid* cutPID2 = NULL;
 248  if (UsePIDDifferentialFlow) {
 249    AliCFTrackCutPid* cutPID2 = new AliCFTrackCutPid("cutPID2","ESD_PID for differential flow") ;
 250    cutPID2->SetPriors(prior);
 251    cutPID2->SetProbabilityCut(0.0);
 252    cutPID2->SetDetectors("TPC TOF");
 253    switch(TMath::Abs(PDG2)) {
 254    case 11   : cutPID2->SetParticleType(AliPID::kElectron, kTRUE); break;
 255    case 13   : cutPID2->SetParticleType(AliPID::kMuon    , kTRUE); break;
 256    case 211  : cutPID2->SetParticleType(AliPID::kPion    , kTRUE); break;
 257    case 321  : cutPID2->SetParticleType(AliPID::kKaon    , kTRUE); break;
 258    case 2212 : cutPID2->SetParticleType(AliPID::kProton  , kTRUE); break;
 259    default   : printf("UNDEFINED PID\n"); break;
 260    }
 261    if (QA) {
 262      cutPID2->SetQAOn(qaIntFE);
 263    }
 264  }
 265
 266  printf("CREATE MC KINE CUTS\n");
 267  TObjArray* mcList1 = new TObjArray(0);
 268  mcList1->AddLast(mcKineCuts1);
 269  mcList1->AddLast(mcGenCuts1);
 270
 271  TObjArray* mcList2 = new TObjArray(0);
 272  mcList2->AddLast(mcKineCuts2);
 273  mcList2->AddLast(mcGenCuts2);
 274
 275  printf("CREATE ACCEPTANCE CUTS\n");
 276  TObjArray* accList1 = new TObjArray(0) ;
 277  accList1->AddLast(mcAccCuts1);
 278
 279  TObjArray* accList2 = new TObjArray(0) ;
 280  accList2->AddLast(mcAccCuts2);
 281
 282  printf("CREATE RECONSTRUCTION CUTS\n");
 283  TObjArray* recList1 = new TObjArray(0) ;
 284  recList1->AddLast(recKineCuts1);
 285  recList1->AddLast(recQualityCuts1);
 286  recList1->AddLast(recIsPrimaryCuts1);
 287
 288  TObjArray* recList2 = new TObjArray(0) ;
 289  recList2->AddLast(recKineCuts2);
 290  recList2->AddLast(recQualityCuts2);
 291  recList2->AddLast(recIsPrimaryCuts2);
 292
 293  printf("CREATE PID CUTS\n");
 294  TObjArray* fPIDCutList1 = new TObjArray(0) ;
 295  if(UsePIDIntegratedFlow) {fPIDCutList1->AddLast(cutPID1);}
 296
 297  TObjArray* fPIDCutList2 = new TObjArray(0) ;
 298  if (UsePIDDifferentialFlow)  {fPIDCutList2->AddLast(cutPID2);}
 299
 300  printf("CREATE INTERFACE AND CUTS\n");
 301  AliCFManager* cfmgr1 = new AliCFManager();
 302  cfmgr1->SetNStepParticle(4);
 303  cfmgr1->SetParticleCutsList(AliCFManager::kPartGenCuts,mcList1);
 304  cfmgr1->SetParticleCutsList(AliCFManager::kPartAccCuts,accList1);
 305  cfmgr1->SetParticleCutsList(AliCFManager::kPartRecCuts,recList1);
 306  cfmgr1->SetParticleCutsList(AliCFManager::kPartSelCuts,fPIDCutList1);
 307
 308  AliCFManager* cfmgr2 = new AliCFManager();
 309  cfmgr2->SetNStepParticle(4);
 310  cfmgr2->SetParticleCutsList(AliCFManager::kPartGenCuts,mcList2);
 311  cfmgr2->SetParticleCutsList(AliCFManager::kPartAccCuts,accList2);
 312  cfmgr2->SetParticleCutsList(AliCFManager::kPartRecCuts,recList2);
 313  cfmgr2->SetParticleCutsList(AliCFManager::kPartSelCuts,fPIDCutList2);
 314
 315  //weights:
 316  TFile *weightsFile = NULL;
 317  TList *weightsList = NULL;
 318
 319  if(useWeights) {
 320    //open the file with the weights:
 321    weightsFile = TFile::Open("weights.root","READ");
 322    if(weightsFile) {
 323      //access the list which holds the histos with weigths:
 324      weightsList = (TList*)weightsFile->Get("weights");
 325    }
 326    else {
 327      cout<<" WARNING: the file <weights.root> with weights from the previous run was not available."<<endl;
 328      break;
 329    }
 330  }
 331
 332
 333  if (LYZ2){
 334    // read the input file from the first run
 335    TString inputFileNameLYZ2 = "outputLYZ1analysis" ;
 336    inputFileNameLYZ2 += type;
 337    inputFileNameLYZ2 += ".root";
 338    cout<<"The input file is "<<inputFileNameLYZ2.Data()<<endl;
 339    TFile* fInputFileLYZ2 = new TFile(inputFileNameLYZ2.Data(),"READ");
 340    if(!fInputFileLYZ2 || fInputFileLYZ2->IsZombie()) {
 341      cerr << " ERROR: NO First Run file... " << endl ;
 342      break;
 343    }
 344    else {
 345      TList* fInputListLYZ2 = (TList*)fInputFileLYZ2->Get("cobjLYZ1");
 346      if (!fInputListLYZ2) {cout<<"list is NULL pointer!"<<endl;}
 347    }
 348    cout<<"LYZ2 input file/list read..."<<endl;
 349  }
 350
 351  if (LYZEP) {
 352    // read the input file from the second LYZ run
 353    TString inputFileNameLYZEP = "outputLYZ2analysis" ;
 354    inputFileNameLYZEP += type;
 355    inputFileNameLYZEP += ".root";
 356    cout<<"The input file is "<<inputFileNameLYZEP.Data()<<endl;
 357    TFile* fInputFileLYZEP = new TFile(inputFileNameLYZEP.Data(),"READ");
 358    if(!fInputFileLYZEP || fInputFileLYZEP->IsZombie()) {
 359      cerr << " ERROR: NO First Run file... " << endl ;
 360      break;
 361    }
 362    else {
 363      TList* fInputListLYZEP = (TList*)fInputFileLYZEP->Get("cobjLYZ2");
 364      if (!fInputListLYZEP) {cout<<"list is NULL pointer!"<<endl;}
 365    }
 366    cout<<"LYZEP input file/list read..."<<endl;
 367  }
 368
 369  //____________________________________________//
 370  // Make the analysis manager
 371  AliAnalysisManager *mgr = new AliAnalysisManager("FlowAnalysisManager");
 372
 373  if (type == "ESD"){
 374    AliVEventHandler* esdH = new AliESDInputHandler;
 375    mgr->SetInputEventHandler(esdH);
 376    if (MCEP) {
 377      AliMCEventHandler *mc = new AliMCEventHandler();
 378      mgr->SetMCtruthEventHandler(mc);}  }
 379
 380  if (type == "AOD"){
 381    AliVEventHandler* aodH = new AliAODInputHandler;
 382    mgr->SetInputEventHandler(aodH);
 383    if (MCEP) {
 384      AliMCEventHandler *mc = new AliMCEventHandler();
 385      mgr->SetMCtruthEventHandler(mc);}  }
 386
 387  if (type == "MC" || type == "ESDMC0" || type == "ESDMC1"){
 388    AliVEventHandler* esdH = new AliESDInputHandler;
 389    mgr->SetInputEventHandler(esdH);
 390
 391    AliMCEventHandler *mc = new AliMCEventHandler();
 392    mgr->SetMCtruthEventHandler(mc); }
 393
 394  //____________________________________________//
 395  // tasks
 396  AliAnalysisTaskFlowEvent *taskFE = NULL;
 397  if (QA) {
 398    taskFE = new AliAnalysisTaskFlowEvent("TaskFlowEvent",kTRUE);
 399    taskFE->SetQAList1(qaIntFE);
 400    taskFE->SetQAList2(qaDiffFE);
 401    taskFE->SetAnalysisType(type);
 402    taskFE->SetCFManager1(cfmgr1);
 403    taskFE->SetCFManager2(cfmgr2);
 404    mgr->AddTask(taskFE);
 405  }
 406  else {
 407    taskFE = new AliAnalysisTaskFlowEvent("TaskFlowEvent",kFALSE);
 408    taskFE->SetAnalysisType(type);
 409    taskFE->SetCFManager1(cfmgr1);
 410    taskFE->SetCFManager2(cfmgr2);
 411    mgr->AddTask(taskFE);
 412  }
 413  if (FQD){
 414    AliAnalysisTaskFittingQDistribution *taskFQD = new AliAnalysisTaskFittingQDistribution("TaskFittingQDistribution",useWeights);
 415    taskFQD->SetUsePhiWeights(usePhiWeights);
 416    mgr->AddTask(taskFQD);
 417  }
 418  if (SP){
 419    AliAnalysisTaskScalarProduct *taskSP = new AliAnalysisTaskScalarProduct("TaskScalarProduct",kFALSE);
 420    taskSP->SetAnalysisType(type);
 421    taskSP->SetCFManager1(cfmgr1);
 422    taskSP->SetCFManager2(cfmgr2);
 423    mgr->AddTask(taskSP);
 424  }
 425  if (LYZ1){
 426    AliAnalysisTaskLeeYangZeros *taskLYZ1 = new AliAnalysisTaskLeeYangZeros("TaskLeeYangZeros",kTRUE,kFALSE);
 427    taskLYZ1->SetAnalysisType(type);
 428    taskLYZ1->SetFirstRunLYZ(kTRUE);
 429    taskLYZ1->SetUseSumLYZ(kTRUE);
 430    taskLYZ1->SetCFManager1(cfmgr1);
 431    taskLYZ1->SetCFManager2(cfmgr2);
 432    mgr->AddTask(taskLYZ1);
 433  }
 434  if (LYZ2){
 435    AliAnalysisTaskLeeYangZeros *taskLYZ2 = new AliAnalysisTaskLeeYangZeros("TaskLeeYangZeros",kFALSE,kFALSE);
 436    taskLYZ2->SetAnalysisType(type);
 437    taskLYZ2->SetFirstRunLYZ(kFALSE);
 438    taskLYZ2->SetUseSumLYZ(kTRUE);
 439    taskLYZ2->SetCFManager1(cfmgr1);
 440    taskLYZ2->SetCFManager2(cfmgr2);
 441    mgr->AddTask(taskLYZ2);
 442  }
 443  if (LYZEP){
 444    AliAnalysisTaskLYZEventPlane *taskLYZEP = new AliAnalysisTaskLYZEventPlane("TaskLYZEventPlane",kFALSE);
 445    taskLYZEP->SetAnalysisType(type);
 446    taskLYZEP->SetCFManager1(cfmgr1);
 447    taskLYZEP->SetCFManager2(cfmgr2);
 448    mgr->AddTask(taskLYZEP);
 449  }
 450  if (GFC){
 451    AliAnalysisTaskCumulants *taskGFC = new AliAnalysisTaskCumulants("TaskCumulants",kFALSE,useWeights);
 452    taskGFC->SetAnalysisType(type);
 453    taskGFC->SetUsePhiWeights(usePhiWeights);
 454    taskGFC->SetUsePtWeights(usePtWeights);
 455    taskGFC->SetUseEtaWeights(useEtaWeights);
 456    taskGFC->SetCFManager1(cfmgr1);
 457    taskGFC->SetCFManager2(cfmgr2);
 458    mgr->AddTask(taskGFC);
 459  }
 460  if (QC){
 461    AliAnalysisTaskQCumulants *taskQC = new AliAnalysisTaskQCumulants("TaskQCumulants",kFALSE,useWeights);
 462    taskQC->SetAnalysisType(type);
 463    taskQC->SetUsePhiWeights(usePhiWeights);
 464    taskQC->SetUsePtWeights(usePtWeights);
 465    taskQC->SetUseEtaWeights(useEtaWeights);
 466    taskQC->SetCFManager1(cfmgr1);
 467    taskQC->SetCFManager2(cfmgr2);
 468    mgr->AddTask(taskQC);
 469  }
 470  if (MCEP){
 471    AliAnalysisTaskMCEventPlane *taskMCEP = new AliAnalysisTaskMCEventPlane("TaskMCEventPlane",kFALSE);
 472    taskMCEP->SetAnalysisType(type);
 473    taskMCEP->SetCFManager1(cfmgr1);
 474    taskMCEP->SetCFManager2(cfmgr2);
 475    mgr->AddTask(taskMCEP);
 476  }
 477
 478
 479  // Create containers for input/output
 480  AliAnalysisDataContainer *cinput1 = mgr->GetCommonInputContainer();
 481  // TString outputFE = "outputFlowEvent";
 482  // outputFE+= type;
 483  // outputFE+= ".root";
 484  AliAnalysisDataContainer *coutputFE = mgr->CreateContainer("cobjFlowEventSimple",  AliFlowEventSimple::Class(),AliAnalysisManager::kExchangeContainer);
 485
 486  if (useWeights) {
 487    AliAnalysisDataContainer *cinputWeights = mgr->CreateContainer("cobjWeights",TList::Class(),AliAnalysisManager::kInputContainer);
 488  }
 489
 490  if (LYZ2){
 491    AliAnalysisDataContainer *cinputLYZ2 = mgr->CreateContainer("cobjLYZ2in",TList::Class(),AliAnalysisManager::kInputContainer); }
 492  if (LYZEP){
 493    AliAnalysisDataContainer *cinputLYZEP = mgr->CreateContainer("cobjLYZEPin",TList::Class(),AliAnalysisManager::kInputContainer); }
 494
 495  if(SP) {
 496    TString outputSP = "outputSPanalysis";
 497    outputSP+= type;
 498    outputSP+= ".root";
 499    AliAnalysisDataContainer *coutputSP = mgr->CreateContainer("cobjSP", TList::Class(),AliAnalysisManager::kOutputContainer,outputSP);
 500  }
 501
 502  if(LYZ1) {
 503    TString outputLYZ1 = "outputLYZ1analysis";
 504    outputLYZ1+= type;
 505    outputLYZ1+= ".root";
 506    AliAnalysisDataContainer *coutputLYZ1 = mgr->CreateContainer("cobjLYZ1", TList::Class(),AliAnalysisManager::kOutputContainer,outputLYZ1);
 507  }
 508
 509  if(LYZ2) {
 510    TString outputLYZ2 = "outputLYZ2analysis";
 511    outputLYZ2+= type;
 512    outputLYZ2+= ".root";
 513    AliAnalysisDataContainer *coutputLYZ2 = mgr->CreateContainer("cobjLYZ2", TList::Class(),AliAnalysisManager::kOutputContainer,outputLYZ2);
 514  }
 515
 516  if(LYZEP) {
 517    TString outputLYZEP = "outputLYZEPanalysis";
 518    outputLYZEP+= type;
 519    outputLYZEP+= ".root";
 520    AliAnalysisDataContainer *coutputLYZEP = mgr->CreateContainer("cobjLYZEP", TList::Class(),AliAnalysisManager::kOutputContainer,outputLYZEP);
 521  }
 522
 523  if(GFC) {
 524    TString outputGFC = "outputGFCanalysis";
 525    outputGFC+= type;
 526    outputGFC+= ".root";
 527    AliAnalysisDataContainer *coutputGFC = mgr->CreateContainer("cobjGFC", TList::Class(),AliAnalysisManager::kOutputContainer,outputGFC);
 528  }
 529
 530  if(QC) {
 531    TString outputQC = "outputQCanalysis";
 532    outputQC+= type;
 533    outputQC+= ".root";
 534    AliAnalysisDataContainer *coutputQC = mgr->CreateContainer("cobjQC", TList::Class(),AliAnalysisManager::kOutputContainer,outputQC);
 535  }
 536
 537  if(FQD) {
 538    TString outputFQD = "outputFQDanalysis";
 539    outputFQD+= type;
 540    outputFQD+= ".root";
 541    AliAnalysisDataContainer *coutputFQD = mgr->CreateContainer("cobjFQD", TList::Class(),AliAnalysisManager::kOutputContainer,outputFQD);
 542  }
 543
 544  if(MCEP) {
 545    TString outputMCEP = "outputMCEPanalysis";
 546    outputMCEP+= type;
 547    outputMCEP+= ".root";
 548    AliAnalysisDataContainer *coutputMCEP = mgr->CreateContainer("cobjMCEP", TList::Class(),AliAnalysisManager::kOutputContainer,outputMCEP);
 549  }
 550
 551
 552  if (QA) {
 553    TString qaNameIntFE = "QAforInt_FE_";
 554    qaNameIntFE += type;
 555    qaNameIntFE += ".root";
 556    AliAnalysisDataContainer *coutputQA1FE =
 557      mgr->CreateContainer("QAintFE", TList::Class(),AliAnalysisManager::kOutputContainer,qaNameIntFE);
 558
 559    TString qaNameDiffFE = "QAforDiff_FE_";
 560    qaNameDiffFE += type;
 561    qaNameDiffFE += ".root";
 562    AliAnalysisDataContainer *coutputQA2FE =
 563      mgr->CreateContainer("QAdiffFE", TList::Class(),AliAnalysisManager::kOutputContainer,qaNameDiffFE);
 564  }
 565
 566
 567  //____________________________________________//
 568
 569
 570  // the flow event simple is produced here
 571  mgr->ConnectInput(taskFE,0,cinput1);
 572  mgr->ConnectOutput(taskFE,0,coutputFE);
 573  if (QA) {
 574    mgr->ConnectOutput(taskFE,1,coutputQA1FE);
 575    mgr->ConnectOutput(taskFE,2,coutputQA2FE);
 576  }
 577
 578  if (FQD) {
 579    mgr->ConnectInput(taskFQD,0,coutputFE);
 580    mgr->ConnectOutput(taskFQD,0,coutputFQD);
 581    if(useWeights) {
 582      mgr->ConnectInput(taskFQD,1,cinputWeights);
 583      cinputWeights->SetData(weightsList);
 584    }
 585  }
 586  if (SP) {
 587    mgr->ConnectInput(taskSP,0,cinput1);
 588    mgr->ConnectOutput(taskSP,0,coutputSP);
 589  }
 590  if (LYZ1) {
 591    mgr->ConnectInput(taskLYZ1,0,cinput1);
 592    mgr->ConnectOutput(taskLYZ1,0,coutputLYZ1);
 593  }
 594  if (LYZ2) {
 595    mgr->ConnectInput(taskLYZ2,0,cinput1);
 596    mgr->ConnectInput(taskLYZ2,1,cinputLYZ2);
 597    mgr->ConnectOutput(taskLYZ2,0,coutputLYZ2);
 598    cinputLYZ2->SetData(fInputListLYZ2);
 599  }
 600  if (LYZEP) {
 601    mgr->ConnectInput(taskLYZEP,0,cinput1);
 602    mgr->ConnectInput(taskLYZEP,1,cinputLYZEP);
 603    mgr->ConnectOutput(taskLYZEP,0,coutputLYZEP);
 604    cinputLYZEP->SetData(fInputListLYZEP);
 605  }
 606  if (GFC) {
 607    mgr->ConnectInput(taskGFC,0,cinput1);
 608    mgr->ConnectOutput(taskGFC,0,coutputGFC);
 609    if (useWeights) {
 610      mgr->ConnectInput(taskGFC,1,cinputWeights);
 611      cinputWeights->SetData(weightsList);
 612    }
 613  }
 614  if (QC) {
 615    mgr->ConnectInput(taskQC,0,cinput1);
 616    mgr->ConnectOutput(taskQC,0,coutputQC);
 617    if (useWeights) {
 618      mgr->ConnectInput(taskQC,1,cinputWeights);
 619      cinputWeights->SetData(weightsList);
 620    }
 621  }
 622  if (MCEP) {
 623    mgr->ConnectInput(taskMCEP,0,cinput1);
 624    mgr->ConnectOutput(taskMCEP,0,coutputMCEP);
 625  }
 626
 627
 628
 629  //----------------------------------------------------------
 630  // Run the analysis
 631  //----------------------------------------------------------
 632
 633  if (!mgr->InitAnalysis()) return;
 634  mgr->PrintStatus();
 635
 636  if (mode==mLocal || mode == mLocalPAR) {
 637    mgr->StartAnalysis("local",chain);
 638  }
 639  else if (mode==mPROOF) {
 640    //  mgr->StartAnalysis("proof",chain);
 641    mgr->StartAnalysis("proof",data,nRuns,offset);
 642  }
 643  else if (mode==mGRID) {
 644    mgr->StartAnalysis("local",chain);
 645  }
 646
 647  timer.Stop();
 648  timer.Print();
 649 }
 650
 651
 652 void LoadLibraries(const anaModes mode) {
 653
 654  //--------------------------------------
 655  // Load the needed libraries most of them already loaded by aliroot
 656  //--------------------------------------
 657  gSystem->Load("libTree.so");
 658  gSystem->Load("libGeom.so");
 659  gSystem->Load("libVMC.so");
 660  gSystem->Load("libXMLIO.so");
 661  gSystem->Load("libPhysics.so");
 662
 663  //----------------------------------------------------------
 664  // >>>>>>>>>>> Local mode <<<<<<<<<<<<<<
 665  //----------------------------------------------------------
 666  if (mode==mLocal) {
 667    //--------------------------------------------------------
 668    // If you want to use already compiled libraries
 669    // in the aliroot distribution
 670    //--------------------------------------------------------
 671    gSystem->Load("libSTEERBase");
 672    gSystem->Load("libESD");
 673    gSystem->Load("libAOD");
 674    gSystem->Load("libANALYSIS");
 675    gSystem->Load("libANALYSISalice");
 676    gSystem->Load("libCORRFW.so");
 677    cerr<<"libCORRFW.so loaded..."<<endl;
 678    gSystem->Load("libPWG2flowCommon.so");
 679    cerr<<"libPWG2flowCommon.so loaded..."<<endl;
 680    gSystem->Load("libPWG2flowTasks.so");
 681    cerr<<"libPWG2flowTasks.so loaded..."<<endl;
 682  }
 683
 684  else if (mode == mLocalPAR || mode == mGRID) {
 685    //--------------------------------------------------------
 686    //If you want to use root and par files from aliroot
 687    //--------------------------------------------------------
 688    SetupPar("STEERBase");
 689    SetupPar("ESD");
 690    SetupPar("AOD");
 691    SetupPar("ANALYSIS");
 692    SetupPar("ANALYSISalice");
 693    SetupPar("PWG2AOD");
 694    SetupPar("CORRFW");
 695    SetupPar("PWG2flowCommon");
 696    cerr<<"PWG2flowCommon.par loaded..."<<endl;
 697    SetupPar("PWG2flowTasks");
 698    cerr<<"PWG2flowTasks.par loaded..."<<endl;
 699  }
 700
 701  //---------------------------------------------------------
 702  // <<<<<<<<<< PROOF mode >>>>>>>>>>>>
 703  //---------------------------------------------------------
 704  else if (mode==mPROOF) {
 705    //
 706
 707    //  set to debug root versus if needed
 708    //  TProof::Mgr("alicecaf")->SetROOTVersion("v5-21-01-alice_dbg");
 709    //  TProof::Mgr("alicecaf")->SetROOTVersion("v5-21-01-alice");
 710
 711    // Connect to proof
 712    // Put appropriate username here
 713    // TProof::Reset("proof://snelling@alicecaf.cern.ch");
 714    printf("*** Connect to PROOF ***\n");
 715    //  TProof::Open("abilandz@alicecaf.cern.ch");
 716    TProof::Open("snelling@localhost");
 717
 718    // Enable the STEERBase Package
 719    gProof->ClearPackage("STEERBase.par");
 720    gProof->UploadPackage("STEERBase.par");
 721    gProof->EnablePackage("STEERBase");
 722    // Enable the ESD Package
 723    gProof->ClearPackage("ESD.par");
 724    gProof->UploadPackage("ESD.par");
 725    gProof->EnablePackage("ESD");
 726    // Enable the AOD Package
 727    gProof->ClearPackage("AOD.par");
 728    gProof->UploadPackage("AOD.par");
 729    gProof->EnablePackage("AOD");
 730    // Enable the Analysis Package
 731    gProof->ClearPackage("ANALYSIS.par");
 732    gProof->UploadPackage("ANALYSIS.par");
 733    gProof->EnablePackage("ANALYSIS");
 734    // Enable the Analysis Package alice
 735    gProof->ClearPackage("ANALYSISalice.par");
 736    gProof->UploadPackage("ANALYSISalice.par");
 737    gProof->EnablePackage("ANALYSISalice");
 738    // Load the PWG2 AOD
 739    gProof->ClearPackage("PWG2AOD.par");
 740    gProof->UploadPackage("PWG2AOD.par");
 741    gProof->EnablePackage("PWG2AOD");
 742    // Enable the Correction Framework
 743    gProof->ClearPackage("CORRFW.par");
 744    gProof->UploadPackage("CORRFW.par");
 745    gProof->EnablePackage("CORRFW");
 746    // Enable Flow Analysis
 747    gProof->ClearPackage("PWG2flowCommon");
 748    gProof->UploadPackage("PWG2flowCommon.par");
 749    gProof->EnablePackage("PWG2flowCommon");
 750    gProof->ClearPackage("PWG2flowTasks");
 751    gProof->UploadPackage("PWG2flowTasks.par");
 752    gProof->EnablePackage("PWG2flowTasks");
 753    //
 754    gProof->ShowEnabledPackages();
 755  }
 756
 757 }
 758
 759 void SetupPar(char* pararchivename) {
 760  //Load par files, create analysis libraries
 761  //For testing, if par file already decompressed and modified
 762  //classes then do not decompress.
 763
 764  TString cdir(Form("%s", gSystem->WorkingDirectory() )) ;
 765  TString parpar(Form("%s.par", pararchivename)) ;
 766  if ( gSystem->AccessPathName(parpar.Data()) ) {
 767    gSystem->ChangeDirectory(gSystem->Getenv("ALICE_ROOT")) ;
 768    TString processline(Form(".! make %s", parpar.Data())) ;
 769    gROOT->ProcessLine(processline.Data()) ;
 770    gSystem->ChangeDirectory(cdir) ;
 771    processline = Form(".! mv /tmp/%s .", parpar.Data()) ;
 772    gROOT->ProcessLine(processline.Data()) ;
 773  }
 774  if ( gSystem->AccessPathName(pararchivename) ) {
 775    TString processline = Form(".! tar xvzf %s",parpar.Data()) ;
 776    gROOT->ProcessLine(processline.Data());
 777  }
 778
 779  TString ocwd = gSystem->WorkingDirectory();
 780  gSystem->ChangeDirectory(pararchivename);
 781
 782  // check for BUILD.sh and execute
 783  if (!gSystem->AccessPathName("PROOF-INF/BUILD.sh")) {
 784    printf("*******************************\n");
 785    printf("*** Building PAR archive    ***\n");
 786    cout<<pararchivename<<endl;
 787    printf("*******************************\n");
 788
 789    if (gSystem->Exec("PROOF-INF/BUILD.sh")) {
 790      Error("runProcess","Cannot Build the PAR Archive! - Abort!");
 791      return -1;
 792    }
 793  }
 794  // check for SETUP.C and execute
 795  if (!gSystem->AccessPathName("PROOF-INF/SETUP.C")) {
 796    printf("*******************************\n");
 797    printf("*** Setup PAR archive       ***\n");
 798    cout<<pararchivename<<endl;
 799    printf("*******************************\n");
 800    gROOT->Macro("PROOF-INF/SETUP.C");
 801  }
 802
 803  gSystem->ChangeDirectory(ocwd.Data());
 804  printf("Current dir: %s\n", ocwd.Data());
 805 }
 806
 807
 808 // Helper macros for creating chains
 809 // from: CreateESDChain.C,v 1.10 jgrosseo Exp
 810
 811 TChain* CreateESDChain(const char* aDataDir, Int_t aRuns, Int_t offset)
 812 {
 813  // creates chain of files in a given directory or file containing a list.
 814  // In case of directory the structure is expected as:
 815  // <aDataDir>/<dir0>/AliESDs.root
 816  // <aDataDir>/<dir1>/AliESDs.root
 817  // ...
 818
 819  if (!aDataDir)
 820    return 0;
 821
 822  Long_t id, size, flags, modtime;
 823  if (gSystem->GetPathInfo(aDataDir, &id, &size, &flags, &modtime))
 824    {
 825      printf("%s not found.\n", aDataDir);
 826      return 0;
 827    }
 828
 829  TChain* chain = new TChain("esdTree");
 830  TChain* chaingAlice = 0;
 831
 832  if (flags & 2)
 833    {
 834      TString execDir(gSystem->pwd());
 835      TSystemDirectory* baseDir = new TSystemDirectory(".", aDataDir);
 836      TList* dirList            = baseDir->GetListOfFiles();
 837      Int_t nDirs               = dirList->GetEntries();
 838      gSystem->cd(execDir);
 839
 840      Int_t count = 0;
 841
 842      for (Int_t iDir=0; iDir<nDirs; ++iDir)
 843         {
 844           TSystemFile* presentDir = (TSystemFile*) dirList->At(iDir);
 845           if (!presentDir || !presentDir->IsDirectory() || strcmp(presentDir->GetName(), ".") == 0 || strcmp(presentDir->GetName(), "..") == 0)
 846             continue;
 847
 848           if (offset > 0)
 849             {
 850               --offset;
 851               continue;
 852             }
 853
 854           if (count++ == aRuns)
 855             break;
 856
 857           TString presentDirName(aDataDir);
 858           presentDirName += "/";
 859           presentDirName += presentDir->GetName();
 860           chain->Add(presentDirName + "/AliESDs.root/esdTree");
 861           //  cerr<<presentDirName<<endl;
 862         }
 863
 864    }
 865  else
 866    {
 867      // Open the input stream
 868      ifstream in;
 869      in.open(aDataDir);
 870
 871      Int_t count = 0;
 872
 873      // Read the input list of files and add them to the chain
 874      TString esdfile;
 875      while(in.good()) {
 876         in >> esdfile;
 877         if (!esdfile.Contains("root")) continue; // protection
 878
 879         if (offset > 0)
 880           {
 881             --offset;
 882             continue;
 883           }
 884
 885         if (count++ == aRuns)
 886           break;
 887
 888        // add esd file
 889         chain->Add(esdfile);
 890      }
 891
 892      in.close();
 893    }
 894
 895  return chain;
 896 }
 897
 898 // Helper macros for creating chains
 899 // from: CreateESDChain.C,v 1.10 jgrosseo Exp
 900
 901 TChain* CreateAODChain(const char* aDataDir, Int_t aRuns, Int_t offset)
 902 {
 903  // creates chain of files in a given directory or file containing a list.
 904  // In case of directory the structure is expected as:
 905  // <aDataDir>/<dir0>/AliAOD.root
 906  // <aDataDir>/<dir1>/AliAOD.root
 907  // ...
 908
 909  if (!aDataDir)
 910    return 0;
 911
 912  Long_t id, size, flags, modtime;
 913  if (gSystem->GetPathInfo(aDataDir, &id, &size, &flags, &modtime))
 914    {
 915      printf("%s not found.\n", aDataDir);
 916      return 0;
 917    }
 918
 919  TChain* chain = new TChain("aodTree");
 920  TChain* chaingAlice = 0;
 921
 922  if (flags & 2)
 923    {
 924      TString execDir(gSystem->pwd());
 925      TSystemDirectory* baseDir = new TSystemDirectory(".", aDataDir);
 926      TList* dirList            = baseDir->GetListOfFiles();
 927      Int_t nDirs               = dirList->GetEntries();
 928      gSystem->cd(execDir);
 929
 930      Int_t count = 0;
 931
 932      for (Int_t iDir=0; iDir<nDirs; ++iDir)
 933         {
 934           TSystemFile* presentDir = (TSystemFile*) dirList->At(iDir);
 935           if (!presentDir || !presentDir->IsDirectory() || strcmp(presentDir->GetName(), ".") == 0 || strcmp(presentDir->GetName(), "..") == 0)
 936             continue;
 937
 938           if (offset > 0)
 939             {
 940               --offset;
 941               continue;
 942             }
 943
 944           if (count++ == aRuns)
 945             break;
 946
 947           TString presentDirName(aDataDir);
 948           presentDirName += "/";
 949           presentDirName += presentDir->GetName();
 950           chain->Add(presentDirName + "/AliAOD.root/aodTree");
 951           // cerr<<presentDirName<<endl;
 952         }
 953
 954    }
 955  else
 956    {
 957      // Open the input stream
 958      ifstream in;
 959      in.open(aDataDir);
 960
 961      Int_t count = 0;
 962
 963      // Read the input list of files and add them to the chain
 964      TString aodfile;
 965      while(in.good()) {
 966         in >> aodfile;
 967         if (!aodfile.Contains("root")) continue; // protection
 968
 969         if (offset > 0)
 970           {
 971             --offset;
 972             continue;
 973           }
 974
 975         if (count++ == aRuns)
 976           break;
 977
 978        // add aod file
 979         chain->Add(aodfile);
 980      }
 981
 982      in.close();
 983    }
 984
 985  return chain;
 986 }
 987
 988