PWG2/FLOW/macros/runAliAnalysisTaskFlow.C

   1 // from CreateESDChain.C - instead of  #include "CreateESDChain.C"
   2 TChain* CreateESDChain(const char* aDataDir = "ESDfiles.txt", Int_t aRuns = 10, Int_t offset = 0) ;
   3 TChain* CreateAODChain(const char* aDataDir = "ESDfiles.txt", Int_t aRuns = 10, Int_t offset = 0) ;
   4 void LookupWrite(TChain* chain, const char* target) ;
   5
   6 //////////////////////////////////////////////////////////
   7 //
   8 // WARNING: THIS MACRO IS NOT TESTED FOR THE OPTION "CUM"
   9 //
  10 //////////////////////////////////////////////////////////
  11
  12
  13 //SETTING THE ANALYSIS
  14
  15 //Flow analysis method can be:
  16 // SP    = Scalar Product
  17 // LYZ1  = Lee Yang Zeroes first run
  18 // LYZ2  = Lee Yang Zeroes second run
  19 // LYZEP = Lee Yang Zeroes Event Plane
  20 // CUM   = Cumulants
  21 // MCEP  = Flow calculated from the real MC event plane (only for simulated data)
  22 const TString method = "SP";
  23
  24 //Type of analysis can be:
  25 // ESD, AOD, MC, ESDMC0, ESDMC1
  26 const TString type = "AOD";
  27
  28 //Bolean to fill/not fill the QA histograms
  29 Bool_t QA = kTRUE;
  30
  31 //SETTING THE CUTS
  32
  33 //for integrated flow
  34 const Double_t ptmin1 = 0.0;
  35 const Double_t ptmax1 = 1000.0;
  36 const Double_t ymin1  = -2.;
  37 const Double_t ymax1  = 2.;
  38 const Int_t mintrackrefsTPC1 = 2;
  39 const Int_t mintrackrefsITS1 = 3;
  40 const Int_t charge1 = 1;
  41 const Int_t PDG1 = 211;
  42 const Int_t minclustersTPC1 = 50;
  43 const Int_t maxnsigmatovertex1 = 3;
  44
  45 //for differential flow
  46 const Double_t ptmin2 = 0.0;
  47 const Double_t ptmax2 = 1000.0;
  48 const Double_t ymin2  = -2.;
  49 const Double_t ymax2  = 2.;
  50 const Int_t mintrackrefsTPC2 = 2;
  51 const Int_t mintrackrefsITS2 = 3;
  52 const Int_t charge2 = 1;
  53 const Int_t PDG2 = 211;
  54 const Int_t minclustersTPC2 = 50;
  55 const Int_t maxnsigmatovertex2 = 3;
  56
  57
  58
  59 //void runAliAnalysisTaskFlow(Int_t nRuns = 10, const Char_t* dataDir="/data/alice1/kolk/TherminatorFIX", Int_t offset = 0)
  60 void runAliAnalysisTaskFlow(Int_t nRuns = 3, const Char_t* dataDir="/data/alice1/kolk/AOD", Int_t offset = 0)
  61 {
  62   TStopwatch timer;
  63   timer.Start();
  64
  65   // include path (to find the .h files when compiling)
  66   gSystem->AddIncludePath("-I$ALICE_ROOT/include") ;
  67   gSystem->AddIncludePath("-I$ROOTSYS/include") ;
  68
  69   // load needed libraries
  70   gSystem->Load("libTree.so");
  71   gSystem->Load("libESD.so");
  72   cerr<<"libESD loaded..."<<endl;
  73   gSystem->Load("libANALYSIS.so");
  74   cerr<<"libANALYSIS.so loaded..."<<endl;
  75   gSystem->Load("libANALYSISRL.so");
  76   cerr<<"libANALYSISRL.so loaded..."<<endl;
  77   gSystem->Load("libCORRFW.so");
  78   cerr<<"libCORRFW.so loaded..."<<endl;
  79   gSystem->Load("libPWG2flow.so");
  80   cerr<<"libPWG2flow.so loaded..."<<endl;
  81
  82   // create the TChain. CreateESDChain() is defined in CreateESDChain.C
  83   if (type!="AOD") { TChain* chain = CreateESDChain(dataDir, nRuns, offset);
  84   cout<<"chain ("<<chain<<")"<<endl; }
  85   else { TChain* chain = CreateAODChain(dataDir, nRuns, offset);
  86   cout<<"chain ("<<chain<<")"<<endl; }
  87
  88   //____________________________________________//
  89   //Create cuts using correction framework
  90
  91   //Set TList for the QA histograms
  92   if (QA) {
  93     TList* qaInt = new TList();
  94     TList* qaDiff = new TList();
  95   }
  96
  97   //############# cuts on MC
  98   AliCFTrackKineCuts* mcKineCuts1 = new AliCFTrackKineCuts("mcKineCuts1","MC-level kinematic cuts");
  99   mcKineCuts1->SetPtRange(ptmin1,ptmax1);
 100   mcKineCuts1->SetRapidityRange(ymin1,ymax1);
 101   mcKineCuts1->SetChargeMC(charge1);
 102   if (QA) { mcKineCuts1->SetQAOn(qaInt); }
 103
 104   AliCFTrackKineCuts* mcKineCuts2 = new AliCFTrackKineCuts("mcKineCuts2","MC-level kinematic cuts");
 105   mcKineCuts2->SetPtRange(ptmin2,ptmax2);
 106   mcKineCuts2->SetRapidityRange(ymin2,ymax2);
 107   mcKineCuts2->SetChargeMC(charge2);
 108   if (QA) { mcKineCuts2->SetQAOn(qaDiff); }
 109
 110   AliCFParticleGenCuts* mcGenCuts1 = new AliCFParticleGenCuts("mcGenCuts1","MC particle generation cuts");
 111   mcGenCuts1->SetRequireIsPrimary();
 112   mcGenCuts1->SetRequirePdgCode(PDG1);
 113   if (QA) { mcGenCuts1->SetQAOn(qaInt); }
 114
 115   AliCFParticleGenCuts* mcGenCuts2 = new AliCFParticleGenCuts("mcGenCuts2","MC particle generation cuts");
 116   mcGenCuts2->SetRequireIsPrimary();
 117   mcGenCuts2->SetRequirePdgCode(PDG2);
 118   if (QA) { mcGenCuts2->SetQAOn(qaDiff); }
 119
 120   //############# Acceptance Cuts
 121   AliCFAcceptanceCuts *mcAccCuts1 = new AliCFAcceptanceCuts("mcAccCuts1","MC acceptance cuts");
 122   mcAccCuts1->SetMinNHitITS(mintrackrefsITS1);
 123   mcAccCuts1->SetMinNHitTPC(mintrackrefsTPC1);
 124   if (QA) { mcAccCuts1->SetQAOn(qaInt); }
 125
 126   AliCFAcceptanceCuts *mcAccCuts2 = new AliCFAcceptanceCuts("mcAccCuts2","MC acceptance cuts");
 127   mcAccCuts2->SetMinNHitITS(mintrackrefsITS2);
 128   mcAccCuts2->SetMinNHitTPC(mintrackrefsTPC2);
 129   if (QA) { mcAccCuts2->SetQAOn(qaDiff); }
 130
 131   //############# Rec-Level kinematic cuts
 132   AliCFTrackKineCuts *recKineCuts1 = new AliCFTrackKineCuts("recKineCuts1","rec-level kine cuts");
 133   recKineCuts1->SetPtRange(ptmin1,ptmax1);
 134   recKineCuts1->SetRapidityRange(ymin1,ymax1);
 135   recKineCuts1->SetChargeRec(charge1);
 136   if (QA) { recKineCuts1->SetQAOn(qaInt); }
 137
 138   AliCFTrackKineCuts *recKineCuts2 = new AliCFTrackKineCuts("recKineCuts2","rec-level kine cuts");
 139   recKineCuts2->SetPtRange(ptmin2,ptmax2);
 140   recKineCuts2->SetRapidityRange(ymin2,ymax2);
 141   recKineCuts2->SetChargeRec(charge2);
 142   if (QA) { recKineCuts2->SetQAOn(qaDiff); }
 143
 144   AliCFTrackQualityCuts *recQualityCuts1 = new AliCFTrackQualityCuts("recQualityCuts1","rec-level quality cuts");
 145   recQualityCuts1->SetMinNClusterTPC(minclustersTPC1);
 146   recQualityCuts1->SetRequireITSRefit(kTRUE);
 147   if (QA) { recQualityCuts1->SetQAOn(qaInt); }
 148
 149   AliCFTrackQualityCuts *recQualityCuts2 = new AliCFTrackQualityCuts("recQualityCuts2","rec-level quality cuts");
 150   recQualityCuts2->SetMinNClusterTPC(minclustersTPC2);
 151   recQualityCuts2->SetRequireITSRefit(kTRUE);
 152   if (QA) { recQualityCuts2->SetQAOn(qaDiff); }
 153
 154   AliCFTrackIsPrimaryCuts *recIsPrimaryCuts1 = new AliCFTrackIsPrimaryCuts("recIsPrimaryCuts1","rec-level isPrimary cuts");
 155   recIsPrimaryCuts1->SetMaxNSigmaToVertex(maxnsigmatovertex1);
 156   if (QA) { recIsPrimaryCuts1->SetQAOn(qaInt); }
 157
 158   AliCFTrackIsPrimaryCuts *recIsPrimaryCuts2 = new AliCFTrackIsPrimaryCuts("recIsPrimaryCuts2","rec-level isPrimary cuts");
 159   recIsPrimaryCuts2->SetMaxNSigmaToVertex(maxnsigmatovertex2);
 160   if (QA) { recIsPrimaryCuts2->SetQAOn(qaDiff); }
 161
 162   AliCFTrackCutPid* cutPID1 = new AliCFTrackCutPid("cutPID1","ESD_PID") ;
 163   AliCFTrackCutPid* cutPID2 = new AliCFTrackCutPid("cutPID2","ESD_PID") ;
 164   int n_species = AliPID::kSPECIES ;
 165   Double_t* prior = new Double_t[n_species];
 166
 167   prior[0] = 0.0244519 ;
 168   prior[1] = 0.0143988 ;
 169   prior[2] = 0.805747  ;
 170   prior[3] = 0.0928785 ;
 171   prior[4] = 0.0625243 ;
 172
 173   cutPID1->SetPriors(prior);
 174   cutPID1->SetProbabilityCut(0.0);
 175   cutPID1->SetDetectors("TPC TOF");
 176   switch(TMath::Abs(PDG1)) {
 177   case 11   : cutPID1->SetParticleType(AliPID::kElectron, kTRUE); break;
 178   case 13   : cutPID1->SetParticleType(AliPID::kMuon    , kTRUE); break;
 179   case 211  : cutPID1->SetParticleType(AliPID::kPion    , kTRUE); break;
 180   case 321  : cutPID1->SetParticleType(AliPID::kKaon    , kTRUE); break;
 181   case 2212 : cutPID1->SetParticleType(AliPID::kProton  , kTRUE); break;
 182   default   : printf("UNDEFINED PID\n"); break;
 183   }
 184
 185   cutPID2->SetPriors(prior);
 186   cutPID2->SetProbabilityCut(0.0);
 187   cutPID2->SetDetectors("TPC TOF");
 188   switch(TMath::Abs(PDG2)) {
 189   case 11   : cutPID2->SetParticleType(AliPID::kElectron, kTRUE); break;
 190   case 13   : cutPID2->SetParticleType(AliPID::kMuon    , kTRUE); break;
 191   case 211  : cutPID2->SetParticleType(AliPID::kPion    , kTRUE); break;
 192   case 321  : cutPID2->SetParticleType(AliPID::kKaon    , kTRUE); break;
 193   case 2212 : cutPID2->SetParticleType(AliPID::kProton  , kTRUE); break;
 194   default   : printf("UNDEFINED PID\n"); break;
 195   }
 196   if (QA) { cutPID1->SetQAOn(qaInt);
 197   cutPID2->SetQAOn(qaDiff); }
 198
 199   printf("CREATE MC KINE CUTS\n");
 200   TObjArray* mcList1 = new TObjArray(0);
 201   mcList1->AddLast(mcKineCuts1);
 202   mcList1->AddLast(mcGenCuts1);
 203
 204   TObjArray* mcList2 = new TObjArray(0);
 205   mcList2->AddLast(mcKineCuts2);
 206   mcList2->AddLast(mcGenCuts2);
 207
 208   printf("CREATE ACCEPTANCE CUTS\n");
 209   TObjArray* accList1 = new TObjArray(0) ;
 210   accList1->AddLast(mcAccCuts1);
 211
 212   TObjArray* accList2 = new TObjArray(0) ;
 213   accList2->AddLast(mcAccCuts2);
 214
 215   printf("CREATE RECONSTRUCTION CUTS\n");
 216   TObjArray* recList1 = new TObjArray(0) ;
 217   recList1->AddLast(recKineCuts1);
 218   recList1->AddLast(recQualityCuts1);
 219   recList1->AddLast(recIsPrimaryCuts1);
 220
 221   TObjArray* recList2 = new TObjArray(0) ;
 222   recList2->AddLast(recKineCuts2);
 223   recList2->AddLast(recQualityCuts2);
 224   recList2->AddLast(recIsPrimaryCuts2);
 225
 226   printf("CREATE PID CUTS\n");
 227   TObjArray* fPIDCutList1 = new TObjArray(0) ;
 228   fPIDCutList1->AddLast(cutPID1);
 229
 230   TObjArray* fPIDCutList2 = new TObjArray(0) ;
 231   fPIDCutList2->AddLast(cutPID2);
 232
 233   printf("CREATE INTERFACE AND CUTS\n");
 234   AliCFManager* cfmgr1 = new AliCFManager();
 235   cfmgr1->SetParticleCutsList(AliCFManager::kPartGenCuts,mcList1);       //on MC
 236   cfmgr1->SetParticleCutsList(AliCFManager::kPartAccCuts,accList1);
 237   cfmgr1->SetParticleCutsList(AliCFManager::kPartRecCuts,recList1);      //on ESD
 238   cfmgr1->SetParticleCutsList(AliCFManager::kPartSelCuts,fPIDCutList1);  //on ESD
 239
 240   AliCFManager* cfmgr2 = new AliCFManager();
 241   cfmgr2->SetParticleCutsList(AliCFManager::kPartGenCuts,mcList2);
 242   cfmgr2->SetParticleCutsList(AliCFManager::kPartAccCuts,accList2);
 243   cfmgr2->SetParticleCutsList(AliCFManager::kPartRecCuts,recList2);
 244   cfmgr2->SetParticleCutsList(AliCFManager::kPartSelCuts,fPIDCutList2);
 245
 246
 247   if (method == "LYZ2"){
 248     // read the input file from the first run
 249     TString inputFileName = "outputLYZ1analysis" ;
 250     inputFileName += type;
 251     inputFileName += "_firstrun.root";
 252     cout<<"The input file is "<<inputFileName.Data()<<endl;
 253     TFile* fInputFile = new TFile(inputFileName.Data(),"READ");
 254     if(!fInputFile || fInputFile->IsZombie()) {
 255       cerr << " ERROR: NO First Run file... " << endl ; }
 256     else {
 257       TList* fInputList = (TList*)fInputFile->Get("cobj1");
 258       if (!fInputList) {cout<<"list is NULL pointer!"<<endl;}
 259     }
 260     cout<<"input file/list read..."<<endl;
 261   }
 262
 263   if (method == "LYZEP") {
 264     // read the input file from the second LYZ run
 265     TString inputFileName = "outputLYZ2analysis" ;
 266     inputFileName += type;
 267     inputFileName += "_secondrun.root";
 268     cout<<"The input file is "<<inputFileName.Data()<<endl;
 269     TFile* fInputFile = new TFile(inputFileName.Data(),"READ");
 270     if(!fInputFile || fInputFile->IsZombie()) { cerr << " ERROR: NO First Run file... " << endl ; }
 271     else {
 272       TList* fInputList = (TList*)fInputFile->Get("cobj1");
 273       if (!fInputList) {cout<<"list is NULL pointer!"<<endl;}
 274     }
 275     cout<<"input file/list read..."<<endl;
 276   }
 277
 278
 279   //____________________________________________//
 280   // Make the analysis manager
 281   AliAnalysisManager *mgr = new AliAnalysisManager("FlowAnalysisManager");
 282
 283   if (type == "ESD"){
 284     AliVEventHandler* esdH = new AliESDInputHandler;
 285     mgr->SetInputEventHandler(esdH);
 286     if (method == "MCEP") {
 287       AliMCEventHandler *mc = new AliMCEventHandler();
 288       mgr->SetMCtruthEventHandler(mc);}  }
 289
 290   if (type == "AOD"){
 291     AliVEventHandler* aodH = new AliAODInputHandler;
 292     mgr->SetInputEventHandler(aodH);
 293     if (method == "MCEP") {
 294       AliMCEventHandler *mc = new AliMCEventHandler();
 295       mgr->SetMCtruthEventHandler(mc);}  }
 296
 297   if (type == "MC" || type == "ESDMC0" || type == "ESDMC1"){
 298     AliVEventHandler* esdH = new AliESDInputHandler;
 299     mgr->SetInputEventHandler(esdH);
 300
 301     AliMCEventHandler *mc = new AliMCEventHandler();
 302     mgr->SetMCtruthEventHandler(mc); }
 303
 304   //____________________________________________//
 305   // 1st MC EP task
 306   if (method == "SP"){
 307     if (QA) { AliAnalysisTaskScalarProduct *taskSP = new AliAnalysisTaskScalarProduct("TaskScalarProduct",kTRUE); }
 308     else { AliAnalysisTaskScalarProduct *taskSP = new AliAnalysisTaskScalarProduct("TaskScalarProduct",kFALSE); }
 309     taskSP->SetAnalysisType(type);
 310     taskSP->SetCFManager1(cfmgr1);
 311     taskSP->SetCFManager2(cfmgr2);
 312     if (QA) {
 313       taskSP->SetQAList1(qaInt);
 314       taskSP->SetQAList2(qaDiff); }
 315     mgr->AddTask(taskSP);
 316   }
 317   else if (method == "LYZ1"){
 318     if (QA) { AliAnalysisTaskLeeYangZeros *taskLYZ1 = new AliAnalysisTaskLeeYangZeros("TaskLeeYangZeros",kTRUE,kTRUE);}
 319     else { AliAnalysisTaskLeeYangZeros *taskLYZ1 = new AliAnalysisTaskLeeYangZeros("TaskLeeYangZeros",kTRUE,kFALSE);}
 320     taskLYZ1->SetAnalysisType(type);
 321     taskLYZ1->SetFirstRunLYZ(kTRUE);
 322     taskLYZ1->SetUseSumLYZ(kTRUE);
 323     taskLYZ1->SetCFManager1(cfmgr1);
 324     taskLYZ1->SetCFManager2(cfmgr2);
 325     if (QA) {
 326       taskLYZ1->SetQAList1(qaInt);
 327       taskLYZ1->SetQAList2(qaDiff);}
 328     mgr->AddTask(taskLYZ1);
 329   }
 330   else if (method == "LYZ2"){
 331     if (QA) { AliAnalysisTaskLeeYangZeros *taskLYZ2 = new AliAnalysisTaskLeeYangZeros("TaskLeeYangZeros",kFALSE,kTRUE);}
 332     else { AliAnalysisTaskLeeYangZeros *taskLYZ2 = new AliAnalysisTaskLeeYangZeros("TaskLeeYangZeros",kFALSE,kFALSE); }
 333     taskLYZ2->SetAnalysisType(type);
 334     taskLYZ2->SetFirstRunLYZ(kFALSE);
 335     taskLYZ2->SetUseSumLYZ(kTRUE);
 336     taskLYZ2->SetCFManager1(cfmgr1);
 337     taskLYZ2->SetCFManager2(cfmgr2);
 338     if (QA) {
 339       taskLYZ2->SetQAList1(qaInt);
 340       taskLYZ2->SetQAList2(qaDiff); }
 341     mgr->AddTask(taskLYZ2);
 342   }
 343   else if (method == "LYZEP"){
 344     if (QA) { AliAnalysisTaskLYZEventPlane *taskLYZEP = new AliAnalysisTaskLYZEventPlane("TaskLYZEventPlane",kTRUE); }
 345     else { AliAnalysisTaskLYZEventPlane *taskLYZEP = new AliAnalysisTaskLYZEventPlane("TaskLYZEventPlane",kFALSE); }
 346     taskLYZEP->SetAnalysisType(type);
 347     taskLYZEP->SetCFManager1(cfmgr1);
 348     taskLYZEP->SetCFManager2(cfmgr2);
 349     if (QA) {
 350       taskLYZEP->SetQAList1(qaInt);
 351       taskLYZEP->SetQAList2(qaDiff); }
 352     mgr->AddTask(taskLYZEP);
 353   }
 354   else if (method == "CUM"){
 355     if (QA) { AliAnalysisTaskCumulants *taskCUM = new AliAnalysisTaskCumulants("TaskCumulants",kTRUE);}
 356     else { AliAnalysisTaskCumulants *taskCUM = new AliAnalysisTaskCumulants("TaskCumulants",kFALSE);}
 357     taskCUM->SetAnalysisType(type);
 358     taskCUM->SetCFManager1(cfmgr1);
 359     taskCUM->SetCFManager2(cfmgr2);
 360     if (QA) {
 361       taskCUM->SetQAList1(qaInt);
 362       taskCUM->SetQAList2(qaDiff); }
 363     mgr->AddTask(taskCUM);
 364   }
 365   else if (method == "MCEP"){
 366     if (QA) { AliAnalysisTaskMCEventPlane *taskMCEP = new AliAnalysisTaskMCEventPlane("TaskMCEventPlane",kTRUE);}
 367     else { AliAnalysisTaskMCEventPlane *taskMCEP = new AliAnalysisTaskMCEventPlane("TaskMCEventPlane",kFALSE);}
 368     taskMCEP->SetAnalysisType(type);
 369     taskMCEP->SetCFManager1(cfmgr1);
 370     taskMCEP->SetCFManager2(cfmgr2);
 371     if (QA) {
 372       taskMCEP->SetQAList1(qaInt);
 373       taskMCEP->SetQAList2(qaDiff); }
 374     mgr->AddTask(taskMCEP);
 375   }
 376
 377
 378   // Create containers for input/output
 379   AliAnalysisDataContainer *cinput1 =
 380     mgr->CreateContainer("cchain1",TChain::Class(),AliAnalysisManager::kInputContainer);
 381
 382   if (method == "LYZ2" || method == "LYZEP"){
 383     AliAnalysisDataContainer *cinput2 =
 384                     mgr->CreateContainer("cobj2",TList::Class(),AliAnalysisManager::kInputContainer); }
 385
 386
 387   TString outputName = "output";
 388   outputName+= method;
 389   outputName+= "analysis";
 390   outputName+= type;
 391   if (method == "LYZ1") outputName+= "_firstrun";
 392   if (method == "LYZ2") outputName+= "_secondrun";
 393   outputName+= ".root";
 394   AliAnalysisDataContainer *coutput1 =
 395     mgr->CreateContainer("cobj1", TList::Class(),AliAnalysisManager::kOutputContainer,outputName);
 396
 397   if (QA) {
 398     TString qaNameInt = "QAforInt_";
 399     qaNameInt += method;
 400     qaNameInt += "_";
 401     qaNameInt += type;
 402     qaNameInt += ".root";
 403     AliAnalysisDataContainer *coutput2 =
 404       mgr->CreateContainer("QAint", TList::Class(),AliAnalysisManager::kOutputContainer,qaNameInt);
 405
 406     TString qaNameDiff = "QAforDiff_";
 407     qaNameDiff += method;
 408     qaNameDiff += "_";
 409     qaNameDiff += type;
 410     qaNameDiff += ".root";
 411     AliAnalysisDataContainer *coutput3 =
 412       mgr->CreateContainer("QAdiff", TList::Class(),AliAnalysisManager::kOutputContainer,qaNameDiff);
 413   }
 414
 415   //____________________________________________//
 416
 417   if (method == "SP")   {
 418     mgr->ConnectInput(taskSP,0,cinput1);
 419     mgr->ConnectOutput(taskSP,0,coutput1);
 420     if (QA) { mgr->ConnectOutput(taskSP,1,coutput2);
 421     mgr->ConnectOutput(taskSP,2,coutput3); }
 422   }
 423   else if (method == "LYZ1") {
 424     mgr->ConnectInput(taskLYZ1,0,cinput1);
 425     mgr->ConnectOutput(taskLYZ1,0,coutput1);
 426     if (QA) { mgr->ConnectOutput(taskLYZ1,1,coutput2);
 427     mgr->ConnectOutput(taskLYZ1,2,coutput3); }
 428   }
 429   else if (method == "LYZ2") {
 430     mgr->ConnectInput(taskLYZ2,0,cinput1);
 431     mgr->ConnectInput(taskLYZ2,1,cinput2);
 432     mgr->ConnectOutput(taskLYZ2,0,coutput1);
 433     if (QA) { mgr->ConnectOutput(taskLYZ2,1,coutput2);
 434     mgr->ConnectOutput(taskLYZ2,2,coutput3); }
 435     cinput2->SetData(fInputList);
 436   }
 437   else if (method == "LYZEP") {
 438     mgr->ConnectInput(taskLYZEP,0,cinput1);
 439     mgr->ConnectInput(taskLYZEP,1,cinput2);
 440     mgr->ConnectOutput(taskLYZEP,0,coutput1);
 441     if (QA) { mgr->ConnectOutput(taskLYZEP,1,coutput2);
 442     mgr->ConnectOutput(taskLYZEP,2,coutput3); }
 443     cinput2->SetData(fInputList);
 444   }
 445   else if (method == "CUM")   {
 446     mgr->ConnectInput(taskCUM,0,cinput1);
 447     mgr->ConnectOutput(taskCUM,0,coutput1);
 448     if (QA) { mgr->ConnectOutput(taskCUM,1,coutput2);
 449     mgr->ConnectOutput(taskCUM,2,coutput3); }
 450   }
 451   else if (method == "MCEP")  {
 452     mgr->ConnectInput(taskMCEP,0,cinput1);
 453     mgr->ConnectOutput(taskMCEP,0,coutput1);
 454     if (QA) { mgr->ConnectOutput(taskMCEP,1,coutput2);
 455     mgr->ConnectOutput(taskMCEP,2,coutput3); }
 456   }
 457
 458   if (!mgr->InitAnalysis()) return;
 459   mgr->PrintStatus();
 460   mgr->StartAnalysis("local",chain);
 461
 462   timer.Stop();
 463   timer.Print();
 464 }
 465
 466
 467 // Helper macros for creating chains
 468 // from: CreateESDChain.C,v 1.10 jgrosseo Exp
 469
 470 TChain* CreateESDChain(const char* aDataDir, Int_t aRuns, Int_t offset)
 471 {
 472   // creates chain of files in a given directory or file containing a list.
 473   // In case of directory the structure is expected as:
 474   // <aDataDir>/<dir0>/AliESDs.root
 475   // <aDataDir>/<dir1>/AliESDs.root
 476   // ...
 477
 478   if (!aDataDir)
 479     return 0;
 480
 481   Long_t id, size, flags, modtime;
 482   if (gSystem->GetPathInfo(aDataDir, &id, &size, &flags, &modtime))
 483     {
 484       printf("%s not found.\n", aDataDir);
 485       return 0;
 486     }
 487
 488   TChain* chain = new TChain("esdTree");
 489   TChain* chaingAlice = 0;
 490
 491   if (flags & 2)
 492     {
 493       TString execDir(gSystem->pwd());
 494       TSystemDirectory* baseDir = new TSystemDirectory(".", aDataDir);
 495       TList* dirList            = baseDir->GetListOfFiles();
 496       Int_t nDirs               = dirList->GetEntries();
 497       gSystem->cd(execDir);
 498
 499       Int_t count = 0;
 500
 501       for (Int_t iDir=0; iDir<nDirs; ++iDir)
 502         {
 503           TSystemFile* presentDir = (TSystemFile*) dirList->At(iDir);
 504           if (!presentDir || !presentDir->IsDirectory() || strcmp(presentDir->GetName(), ".") == 0 || strcmp(presentDir->GetName(), "..") == 0)
 505             continue;
 506
 507           if (offset > 0)
 508             {
 509               --offset;
 510               continue;
 511             }
 512
 513           if (count++ == aRuns)
 514             break;
 515
 516           TString presentDirName(aDataDir);
 517           presentDirName += "/";
 518           presentDirName += presentDir->GetName();
 519           chain->Add(presentDirName + "/AliESDs.root/esdTree");
 520           cerr<<presentDirName<<endl;
 521         }
 522
 523     }
 524   else
 525     {
 526       // Open the input stream
 527       ifstream in;
 528       in.open(aDataDir);
 529
 530       Int_t count = 0;
 531
 532       // Read the input list of files and add them to the chain
 533       TString esdfile;
 534       while(in.good()) {
 535         in >> esdfile;
 536         if (!esdfile.Contains("root")) continue; // protection
 537
 538         if (offset > 0)
 539           {
 540             --offset;
 541             continue;
 542           }
 543
 544         if (count++ == aRuns)
 545           break;
 546
 547         // add esd file
 548         chain->Add(esdfile);
 549       }
 550
 551       in.close();
 552     }
 553
 554   return chain;
 555 }
 556
 557 // Helper macros for creating chains
 558 // from: CreateESDChain.C,v 1.10 jgrosseo Exp
 559
 560 TChain* CreateAODChain(const char* aDataDir, Int_t aRuns, Int_t offset)
 561 {
 562   // creates chain of files in a given directory or file containing a list.
 563   // In case of directory the structure is expected as:
 564   // <aDataDir>/<dir0>/AliAOD.root
 565   // <aDataDir>/<dir1>/AliAOD.root
 566   // ...
 567
 568   if (!aDataDir)
 569     return 0;
 570
 571   Long_t id, size, flags, modtime;
 572   if (gSystem->GetPathInfo(aDataDir, &id, &size, &flags, &modtime))
 573     {
 574       printf("%s not found.\n", aDataDir);
 575       return 0;
 576     }
 577
 578   TChain* chain = new TChain("aodTree");
 579   TChain* chaingAlice = 0;
 580
 581   if (flags & 2)
 582     {
 583       TString execDir(gSystem->pwd());
 584       TSystemDirectory* baseDir = new TSystemDirectory(".", aDataDir);
 585       TList* dirList            = baseDir->GetListOfFiles();
 586       Int_t nDirs               = dirList->GetEntries();
 587       gSystem->cd(execDir);
 588
 589       Int_t count = 0;
 590
 591       for (Int_t iDir=0; iDir<nDirs; ++iDir)
 592         {
 593           TSystemFile* presentDir = (TSystemFile*) dirList->At(iDir);
 594           if (!presentDir || !presentDir->IsDirectory() || strcmp(presentDir->GetName(), ".") == 0 || strcmp(presentDir->GetName(), "..") == 0)
 595             continue;
 596
 597           if (offset > 0)
 598             {
 599               --offset;
 600               continue;
 601             }
 602
 603           if (count++ == aRuns)
 604             break;
 605
 606           TString presentDirName(aDataDir);
 607           presentDirName += "/";
 608           presentDirName += presentDir->GetName();
 609           chain->Add(presentDirName + "/AliAOD.root/aodTree");
 610           cerr<<presentDirName<<endl;
 611         }
 612
 613     }
 614   else
 615     {
 616       // Open the input stream
 617       ifstream in;
 618       in.open(aDataDir);
 619
 620       Int_t count = 0;
 621
 622       // Read the input list of files and add them to the chain
 623       TString aodfile;
 624       while(in.good()) {
 625         in >> aodfile;
 626         if (!aodfile.Contains("root")) continue; // protection
 627
 628         if (offset > 0)
 629           {
 630             --offset;
 631             continue;
 632           }
 633
 634         if (count++ == aRuns)
 635           break;
 636
 637         // add aod file
 638         chain->Add(aodfile);
 639       }
 640
 641       in.close();
 642     }
 643
 644   return chain;
 645 }
 646
 647
 648
 649 void LookupWrite(TChain* chain, const char* target)
 650 {
 651   // looks up the chain and writes the remaining files to the text file target
 652
 653   chain->Lookup();
 654
 655   TObjArray* list = chain->GetListOfFiles();
 656   TIterator* iter = list->MakeIterator();
 657   TObject* obj = 0;
 658
 659   ofstream outfile;
 660   outfile.open(target);
 661
 662   while ((obj = iter->Next()))
 663     outfile << obj->GetTitle() << "#AliESDs.root" << endl;
 664
 665   outfile.close();
 666
 667   delete iter;
 668 }