PWG4/macros/ana.C

   1 /* $Id:  $ */
   2 //--------------------------------------------------
   3 // Example macro to do analysis with the
   4 // analysis classes in CaloTrackCorrelations
   5 // Can be executed with Root and AliRoot
   6 //
   7 // Pay attention to the options and definitions
   8 // set in the lines below
   9 //
  10 //  Author : Gustavo Conesa Balbastre (INFN-LNF)
  11 //
  12 //-------------------------------------------------
  13 enum anaModes {mLocal=0, mPROOF=1, mPlugin=2, mGRID=3};
  14 //mLocal    = 0: Analyze locally files in your computer
  15 //mPROOF    = 1: Analyze files on GRID with Plugin
  16 //mPlugin   = 2: Analyze files on GRID with Plugin
  17 //mGRID     = 3: Analyze files on GRID, jobs launched from aliensh
  18
  19 //---------------------------------------------------------------------------
  20 // Settings to read locally several files, only for "mLocal" mode
  21 // The different values are default, they can be set with environmental
  22 // variables: INDIR, PATTERN, NFILES, respectivelly
  23
  24 char * kInDir   = "/user/data/files/";
  25 char * kPattern = ""; // Data are in files kInDir/kPattern+i
  26 Int_t  kFile    = 1;
  27
  28 //---------------------------------------------------------------------------
  29 // Dataset for proof analysis, mode=mPROOF
  30 // char * kDataset = "/alice/vernet/PbPb_LHC10h_ESD";
  31
  32 char *  kDatasetPROOF     = "/alice/vernet/LHC11b_149646";
  33 //char *  kDatasetPROOF     = "/alice/vernet/LHC11b10a_AOD046";//LHC11d_AOD076
  34 Int_t   kDatasetNMaxFiles = 20;
  35 TString ccin2p3UserName   = "arbor" ;
  36 TString alienUserName     = "narbor" ;
  37
  38 //---------------------------------------------------------------------------
  39 // Collection file for grid analysis
  40
  41 char * kXML = "collection.xml";
  42
  43 //---------------------------------------------------------------------------
  44 //Scale histograms from file. Change to kTRUE when xsection file exists
  45 //Put name of file containing xsection
  46 //Put number of events per ESD file
  47 //This is an specific case for normalization of Pythia files.
  48 const Bool_t kGetXSectionFromFileAndScale = kFALSE ;
  49 const char * kXSFileName = "pyxsec.root";
  50
  51 //---------------------------------------------------------------------------
  52
  53 //Set some default values, but used values are set in the code!
  54
  55 Bool_t  kMC        = kFALSE; //With real data kMC = kFALSE
  56 TString kInputData = "ESD"; //ESD, AOD, MC, deltaAOD
  57 Int_t   kYear      = 2011;
  58 TString kCollision = "pp";
  59 Bool_t  outAOD     = kTRUE; //Some tasks doesnt need it.
  60 TString kTreeName;
  61 TString kPass      = "";
  62 char    kTrigger[1024];
  63 Int_t   kRun       = 0;
  64
  65 //________________________
  66 void ana(Int_t mode=mGRID)
  67 {
  68   // Main
  69
  70   //--------------------------------------------------------------------
  71   // Load analysis libraries
  72   // Look at the method below,
  73   // change whatever you need for your analysis case
  74   // ------------------------------------------------------------------
  75
  76   LoadLibraries(mode) ;
  77   //gSystem->ListLibraries();
  78
  79   //-------------------------------------------------------------------------------------------------
  80   //Create chain from ESD and from cross sections files, look below for options.
  81   //-------------------------------------------------------------------------------------------------
  82
  83   // Set kInputData and kTreeName looking to the kINDIR
  84
  85   CheckInputData(mode);
  86
  87   // Check global analysis settings
  88
  89   CheckEnvironmentVariables();
  90
  91   printf("*********************************************\n");
  92   printf("*** Input data < %s >, pass %s, tree < %s >, MC?  < %d > ***\n",kInputData.Data(),kPass.Data(),kTreeName.Data(),kMC);
  93   printf("*********************************************\n");
  94
  95   TChain * chain   = new TChain(kTreeName) ;
  96   TChain * chainxs = new TChain("Xsection") ;
  97   CreateChain(mode, chain, chainxs);
  98
  99   printf("*********************************************\n");
 100   printf("number of entries # %lld, skipped %d\n", chain->GetEntries()) ;
 101   printf("*********************************************\n");
 102
 103   if(!chain)
 104   {
 105     printf("STOP, no chain available\n");
 106     return;
 107   }
 108
 109   AliLog::SetGlobalLogLevel(AliLog::kError);//Minimum prints on screen
 110
 111   //------------------------------------------
 112   //  Alien handler part
 113   //------------------------------------------
 114   AliAnalysisGrid *alienHandler=0x0;
 115   if(mode==mPlugin)
 116   {
 117     // Create and configure the alien handler plugin
 118     gROOT->LoadMacro("CreateAlienHandler.C");
 119     alienHandler = CreateAlienHandler();
 120     if (!alienHandler) return;
 121   }
 122
 123   //--------------------------------------
 124   // Make the analysis manager
 125   //-------------------------------------
 126   AliAnalysisManager *mgr  = new AliAnalysisManager("Manager", "Manager");
 127   //AliAnalysisManager::SetUseProgressBar(kTRUE);
 128   //mgr->SetSkipTerminate(kTRUE);
 129   //mgr->SetNSysInfo(1);
 130
 131   if(mode==mPlugin)
 132   {
 133     // Connect plugin to the analysis manager
 134     mgr->SetGridHandler(alienHandler);
 135   }
 136
 137   // MC handler
 138   if((kMC || kInputData == "MC") && !kInputData.Contains("AOD"))
 139   {
 140     AliMCEventHandler* mcHandler = new AliMCEventHandler();
 141     mcHandler->SetReadTR(kFALSE);//Do not search TrackRef file
 142     mgr->SetMCtruthEventHandler(mcHandler);
 143     if( kInputData == "MC")
 144     {
 145       cout<<"MC INPUT EVENT HANDLER"<<endl;
 146       mgr->SetInputEventHandler(NULL);
 147     }
 148   }
 149
 150
 151   // AOD output handler
 152   if(kInputData!="deltaAOD" && outAOD)
 153   {
 154     cout<<"Init output handler"<<endl;
 155     AliAODHandler* aodoutHandler   = new AliAODHandler();
 156     aodoutHandler->SetOutputFileName("aod.root");
 157     ////aodoutHandler->SetCreateNonStandardAOD();
 158     mgr->SetOutputEventHandler(aodoutHandler);
 159   }
 160
 161   //input
 162
 163   if(kInputData == "ESD")
 164   {
 165     // ESD handler
 166     AliESDInputHandler *esdHandler = new AliESDInputHandler();
 167     esdHandler->SetReadFriends(kFALSE);
 168     mgr->SetInputEventHandler(esdHandler);
 169     cout<<"ESD handler "<<mgr->GetInputEventHandler()<<endl;
 170   }
 171   else if(kInputData.Contains("AOD"))
 172   {
 173     // AOD handler
 174     AliAODInputHandler *aodHandler = new AliAODInputHandler();
 175     mgr->SetInputEventHandler(aodHandler);
 176     if(kInputData == "deltaAOD") aodHandler->AddFriend("deltaAODCaloTrackCorr.root");
 177     cout<<"AOD handler "<<mgr->GetInputEventHandler()<<endl;
 178   }
 179   //mgr->RegisterExternalFile("deltaAODCaloTrackCorr.root");
 180   //mgr->SetDebugLevel(1); // For debugging, do not uncomment if you want no messages.
 181
 182   TString outputFile = AliAnalysisManager::GetCommonFileName();
 183
 184   //-------------------------------------------------------------------------
 185   //Define task, put here any other task that you want to use.
 186   //-------------------------------------------------------------------------
 187
 188   // Physics selection
 189   gROOT->LoadMacro("$ALICE_ROOT/ANALYSIS/macros/AddTaskPhysicsSelection.C");
 190   AliPhysicsSelectionTask* physSelTask = AddTaskPhysicsSelection(kMC);
 191
 192   // Centrality
 193   if(kCollision=="PbPb")
 194   {
 195     gROOT->LoadMacro("$ALICE_ROOT/ANALYSIS/macros/AddTaskCentrality.C");
 196     AliCentralitySelectionTask *taskCentrality = AddTaskCentrality();
 197   }
 198
 199   // Simple event counting tasks
 200   AddTaskCounter("");   // All
 201   //AddTaskCounter("MB"); // Min Bias
 202   if(!kMC)
 203   {
 204     AddTaskCounter("INT7"); // Min Bias
 205     //AddTaskCounter("EMC1"); // Trig Th > 1.5 GeV approx
 206     AddTaskCounter("EMC7"); // Trig Th > 4-5 GeV
 207     AddTaskCounter("PHOS"); //
 208   }
 209
 210   // -----------------
 211   // Photon conversion
 212   // -----------------
 213 /*
 214   if(kInputData=="ESD"){
 215     printf("* Configure photon conversion analysis in macro \n");
 216     TString arguments = "-run-on-train -use-own-xyz  -force-aod -mc-off ";
 217     gROOT->LoadMacro("$ALICE_ROOT/PWG4/macros/ConfigGammaConversion.C");
 218     AliAnalysisTaskGammaConversion * taskGammaConversion =
 219     ConfigGammaConversion(arguments,mgr->GetCommonInputContainer());
 220     taskGammaConversion->SelectCollisionCandidates();
 221
 222     // Gamma Conversion AOD to AODPWG4Particle
 223     AliAnalysisTaskGCPartToPWG4Part * taskGCToPC = new AliAnalysisTaskGCPartToPWG4Part("GCPartToPWG4Part");
 224     taskGCToPC->SetGammaCutId("90035620401003321022000000090");
 225     mgr->AddTask(taskGCToPC);
 226     mgr->ConnectInput  (taskGCToPC, 0, mgr->GetCommonInputContainer() );
 227     mgr->ConnectOutput (taskGCToPC, 0, mgr->GetCommonOutputContainer());
 228   }
 229 */
 230
 231   Bool_t kPrint   = kFALSE;
 232   Bool_t deltaAOD = kFALSE;
 233   gROOT->LoadMacro("AddTaskCaloTrackCorr.C");   // $ALICE_ROOT/PWG4/macros
 234   gROOT->LoadMacro("AddTaskEMCALClusterize.C"); // $ALICE_ROOT/PWG4/CaloCalib/macros
 235
 236
 237   // ------
 238   // Tracks
 239   // ------
 240
 241   // Track isolation-correlation analysis and EMCAL QA analysis
 242   AliAnalysisTaskCaloTrackCorrelation *anamb  = AddTaskCaloTrackCorr(kInputData, "EMCAL",   kPrint,kMC, deltaAOD,  outputFile.Data(),
 243                                                                kYear,kRun,kCollision,"INT7","");   // PHOS trigger
 244
 245   AliAnalysisTaskCaloTrackCorrelation *anatr  = AddTaskCaloTrackCorr(kInputData, "EMCAL",   kPrint,kMC, deltaAOD,  outputFile.Data(),
 246                                                                kYear,kRun,kCollision,"EMC7","");
 247
 248
 249   // -----
 250   // EMCAL
 251   // -----
 252
 253   Bool_t  bTrackMatch = kTRUE;
 254   Int_t   minEcell    = 50;  // 50  MeV (10 MeV used in reconstruction)
 255   Int_t   minEseed    = 100; // 100 MeV
 256   Int_t   dTime       = 0;   // default, 250 ns
 257   Int_t   wTime       = 0;   // default 425 < T < 825 ns
 258   TString clTrigger   = "";  // Do not select, do Min Bias and triggered
 259
 260   //Analysis with clusterizer V1
 261   AliAnalysisTaskEMCALClusterize * clv1 = AddTaskEMCALClusterize(kMC,"V1",clTrigger,kRun,kPass, bTrackMatch,
 262                                                                  minEcell,minEseed,dTime,wTime);
 263
 264   TString arrayNameV1(Form("V1_Ecell%d_Eseed%d_DT%d_WT%d",minEcell,minEseed, dTime,wTime));
 265   printf("Name of clusterizer array: %s\n",arrayNameV1.Data());
 266
 267   if(!kMC)
 268   {
 269
 270     AliAnalysisTaskCaloTrackCorrelation *anav1tr  = AddTaskCaloTrackCorr(kInputData, "EMCAL",   kPrint,kMC, deltaAOD,  outputFile.Data(),
 271                                                                    kYear,kRun,kCollision,"EMC7",arrayNameV1);
 272
 273     AliAnalysisTaskCaloTrackCorrelation *anav1mb  = AddTaskCaloTrackCorr(kInputData, "EMCAL",   kPrint,kMC, deltaAOD,  outputFile.Data(),
 274                                                                    kYear,kRun,kCollision,"INT7",arrayNameV1);
 275   }
 276   else
 277   {// No trigger (should be MB, but for single particle productions it does not work)
 278
 279     AliAnalysisTaskCaloTrackCorrelation *anav1  = AddTaskCaloTrackCorr(kInputData, "EMCAL",   kPrint,kMC, deltaAOD,  outputFile.Data(),
 280                                                                  kYear,kRun,kCollision,"",arrayNameV1);
 281   }
 282
 283
 284
 285   //Analysis with clusterizer V2
 286   AliAnalysisTaskEMCALClusterize * clv2 = AddTaskEMCALClusterize(kMC,"V2",clTrigger,kRun,kPass, bTrackMatch,
 287                                                                  minEcell,minEseed,dTime,wTime);
 288
 289   TString arrayNameV2(Form("V2_Ecell%d_Eseed%d_DT%d_WT%d",minEcell,minEseed, dTime,wTime));
 290   printf("Name of clusterizer array: %s\n",arrayNameV2.Data());
 291
 292   if(!kMC)
 293   {
 294
 295     AliAnalysisTaskCaloTrackCorrelation *anav2tr  = AddTaskCaloTrackCorr(kInputData, "EMCAL",   kPrint,kMC, deltaAOD,  outputFile.Data(),
 296                                                                    kYear,kRun,kCollision,"EMC7",arrayNameV2);
 297
 298     AliAnalysisTaskCaloTrackCorrelation *anav2mb  = AddTaskCaloTrackCorr(kInputData, "EMCAL",   kPrint,kMC, deltaAOD,  outputFile.Data(),
 299                                                                    kYear,kRun,kCollision,"INT7",arrayNameV2);
 300   }
 301   else
 302   {// No trigger (should be MB, but for single particle productions it does not work)
 303     AliAnalysisTaskCaloTrackCorrelation *anav2  = AddTaskCaloTrackCorr(kInputData, "EMCAL",   kPrint,kMC, deltaAOD,  outputFile.Data(),
 304                                                                  kYear,kRun,kCollision,"",arrayNameV2);
 305   }
 306
 307   /*
 308   // -----
 309   // PHOS
 310   // -----
 311
 312   //Add here PHOS tender or whatever is needed
 313
 314   if(!kMC)
 315   {
 316
 317     AliAnalysisTaskCaloTrackCorrelation *anav1tr = AddTaskCaloTrackCorr(kInputData, "PHOS", kPrint,kMC, deltaAOD,  outputFile.Data(),
 318                                                                   kYear,kRun,kCollision,"PHOS",""); // PHOS trigger
 319
 320     AliAnalysisTaskCaloTrackCorrelation *anav1mb = AddTaskCaloTrackCorr(kInputData, "PHOS",   kPrint,kMC, deltaAOD,  outputFile.Data(),
 321                                                                   kYear,kRun,kCollision,"INT7","");
 322
 323   }
 324   else
 325   {// No trigger
 326
 327     AliAnalysisTaskCaloTrackCorrelation *anav1mb = AddTaskCaloTrackCorr(kInputData, "PHOS",   kPrint,kMC, deltaAOD,  outputFile.Data(),
 328                                                                   kYear,kRun,kCollision,"","");
 329
 330   }
 331
 332   */
 333   //-----------------------
 334   // Run the analysis
 335   //-----------------------
 336   mgr->InitAnalysis();
 337   mgr->PrintStatus();
 338
 339   if      (mode == mPlugin) mgr->StartAnalysis("grid");
 340   else if (mode == mPROOF ) mgr->StartAnalysis("proof",chain);
 341   else                      mgr->StartAnalysis("local",chain);
 342
 343   cout <<" Analysis ended sucessfully "<< endl ;
 344
 345 }
 346
 347 //_____________________________
 348 void  LoadLibraries(Int_t mode)
 349 {
 350
 351   if (mode == mPROOF) {
 352     //TProof::Mgr("ccalpmaster")->SetROOTVersion("ALICE_v5-27-06b");
 353     gROOT->LoadMacro("/afs/in2p3.fr/group/alice/laf/EnableAliRootForLAF.C");
 354     TProof* proof = EnableAliRootForLAF("ccaplmaster",nPROOFWorkers.Data(),ccin2p3UserName.Data(),alienUserName.Data(),"",kFALSE,kTRUE,kTRUE,"OADB:ANALYSIS:ANALYSISalice:AOD:ESD:CORRFW:STEERBase:EMCALUtils:PHOSUtils:PWG4PartCorrBase:PWG4PartCorrDep:PWG4CaloCalib");
 355
 356     //  TProof* proof = TProof::Open("ccaplmaster",Form("workers=%s",nPROOFWorkers.Data()));
 357
 358     //     //proof->ClearPackages();
 359     //     proof->UploadPackage("STEERBase");
 360     //     proof->UploadPackage("ESD");
 361     //     proof->UploadPackage("AOD");
 362     //     proof->UploadPackage("ANALYSIS");
 363     //     proof->UploadPackage("OADB");
 364     //     proof->UploadPackage("ANALYSISalice");
 365     //     proof->UploadPackage("CORRFW");
 366     //     //proof->UploadPackage("JETAN");
 367     //     proof->UploadPackage("PHOSUtils");
 368     //     proof->UploadPackage("EMCALUtils");
 369     //     proof->UploadPackage("PWG4PartCorrBase");
 370     //     proof->UploadPackage("PWG4PartCorrDep");
 371     //     proof->UploadPackage("PWG4CaloCalib");
 372
 373     //     proof->EnablePackage("STEERBase");
 374     //     proof->EnablePackage("ESD");
 375     //     proof->EnablePackage("AOD");
 376     //     proof->EnablePackage("ANALYSIS");
 377     //     proof->EnablePackage("OADB");
 378     //     proof->EnablePackage("ANALYSISalice");
 379     //     proof->EnablePackage("CORRFW");
 380     //     //proof->EnablePackage("JETAN");
 381     //     proof->EnablePackage("PHOSUtils");
 382     //     proof->EnablePackage("EMCALUtils");
 383     //     proof->EnablePackage("PWG4PartCorrBase");
 384     //     proof->EnablePackage("PWG4PartCorrDep");
 385     //     proof->EnablePackage("PWG4CaloCalib");
 386     return;
 387   }
 388
 389   //--------------------------------------
 390   // Load the needed libraries most of them already loaded by aliroot
 391   //--------------------------------------
 392   gSystem->Load("libTree.so");
 393   gSystem->Load("libGeom.so");
 394   gSystem->Load("libVMC.so");
 395   gSystem->Load("libXMLIO.so");
 396   gSystem->Load("libMatrix.so");
 397   gSystem->Load("libPhysics.so");
 398   gSystem->Load("libMinuit.so"); // Root + libraries to if reclusterization is done
 399
 400   gSystem->Load("libSTEERBase.so");
 401   gSystem->Load("libGui.so"); // Root + libraries to if reclusterization is done
 402   gSystem->Load("libCDB.so"); // Root + libraries to if reclusterization is done
 403   gSystem->Load("libESD.so"); // Root + libraries to if reclusterization is done
 404   gSystem->Load("libAOD.so");
 405   gSystem->Load("libRAWDatabase.so"); // Root + libraries to if reclusterization is done
 406   gSystem->Load("libProof.so");
 407   gSystem->Load("libANALYSIS.so");
 408   gSystem->Load("libSTEER.so"); // Root + libraries to if reclusterization is done
 409
 410   gSystem->Load("libRAWDatarec.so"); // Root + libraries to if reclusterization is done
 411   gSystem->Load("libRAWDatasim.so"); // Root + libraries to if reclusterization is done
 412   gSystem->Load("libVZERObase.so");  // Root + libraries to if reclusterization is done
 413   gSystem->Load("libVZEROrec.so");   // Root + libraries to if reclusterization is done
 414
 415   gSystem->Load("libEMCALUtils");
 416   //SetupPar("EMCALUtils");
 417   gSystem->Load("libEMCALraw");  // Root + libraries to if reclusterization is done
 418   gSystem->Load("libEMCALbase"); // Root + libraries to if reclusterization is done
 419   gSystem->Load("libEMCALsim");  // Root + libraries to if reclusterization is done
 420   gSystem->Load("libEMCALrec");  // Root + libraries to if reclusterization is done
 421   //SetupPar("EMCALraw");
 422   //SetupPar("EMCALbase");
 423   //SetupPar("EMCALsim");
 424   //SetupPar("EMCALrec");
 425
 426   gSystem->Load("libANALYSISalice.so");
 427   //gSystem->Load("libTENDER.so");
 428   //gSystem->Load("libTENDERSupplies.so");
 429
 430   gSystem->Load("libPHOSUtils");
 431   gSystem->Load("libEMCALUtils");
 432   gSystem->Load("libPWG4PartCorrBase");
 433   gSystem->Load("libPWG4PartCorrDep");
 434   gSystem->Load("libPWG4CaloCalib");
 435   //SetupPar("PWG4PartCorrBase");
 436   //SetupPar("PWG4PartCorrDep");
 437   //SetupPar("PWG4CaloCalib");
 438
 439   //gSystem->Load("libJETAN");
 440   //gSystem->Load("FASTJETAN");
 441   //gSystem->Load("PWG4JetTasks");
 442
 443   gSystem->Load("libCORRFW.so");
 444   gSystem->Load("libPWG4GammaConv.so");
 445   //SetupPar("PWG4GammaConv");
 446
 447   // needed for plugin?
 448   gSystem->AddIncludePath("-I$ALICE_ROOT");
 449   gSystem->AddIncludePath("-I./");
 450
 451 }
 452
 453 //_________________________________
 454 void SetupPar(char* pararchivename)
 455 {
 456   //Load par files, create analysis libraries
 457   //For testing, if par file already decompressed and modified
 458   //classes then do not decompress.
 459
 460   TString cdir(Form("%s", gSystem->WorkingDirectory() )) ;
 461   TString parpar(Form("%s.par", pararchivename)) ;
 462
 463   if ( gSystem->AccessPathName(pararchivename) ) {
 464     TString processline = Form(".! tar xvzf %s",parpar.Data()) ;
 465     gROOT->ProcessLine(processline.Data());
 466   }
 467
 468   TString ocwd = gSystem->WorkingDirectory();
 469   gSystem->ChangeDirectory(pararchivename);
 470
 471   // check for BUILD.sh and execute
 472   if (!gSystem->AccessPathName("PROOF-INF/BUILD.sh")) {
 473     printf("*******************************\n");
 474     printf("*** Building PAR archive    ***\n");
 475     cout<<pararchivename<<endl;
 476     printf("*******************************\n");
 477
 478     if (gSystem->Exec("PROOF-INF/BUILD.sh")) {
 479       Error("runProcess","Cannot Build the PAR Archive! - Abort!");
 480       return -1;
 481     }
 482   }
 483   // check for SETUP.C and execute
 484   if (!gSystem->AccessPathName("PROOF-INF/SETUP.C")) {
 485     printf("*******************************\n");
 486     printf("*** Setup PAR archive       ***\n");
 487     cout<<pararchivename<<endl;
 488     printf("*******************************\n");
 489     gROOT->Macro("PROOF-INF/SETUP.C");
 490   }
 491
 492   gSystem->ChangeDirectory(ocwd.Data());
 493   printf("Current dir: %s\n", ocwd.Data());
 494 }
 495
 496 //______________________________________
 497 void CheckInputData(const anaModes mode)
 498 {
 499   //Sets input data and tree
 500
 501   TString ocwd = gSystem->WorkingDirectory();
 502
 503   //---------------------------------------
 504   //Local files analysis
 505   //---------------------------------------
 506   if(mode == mLocal){
 507     //If you want to add several ESD files sitting in a common directory INDIR
 508     //Specify as environmental variables the directory (INDIR), the number of files
 509     //to analyze (NFILES) and the pattern name of the directories with files (PATTERN)
 510
 511     if(gSystem->Getenv("INDIR"))
 512       kInDir = gSystem->Getenv("INDIR") ;
 513     else cout<<"INDIR not set, use default: "<<kInDir<<endl;
 514
 515     TString sindir(kInDir);
 516     if     (sindir.Contains("pass1")) kPass = "pass1";
 517     else if(sindir.Contains("pass2")) kPass = "pass2";
 518     else if(sindir.Contains("pass3")) kPass = "pass3";
 519
 520     if(gSystem->Getenv("PATTERN"))
 521       kPattern = gSystem->Getenv("PATTERN") ;
 522     else  cout<<"PATTERN not set, use default: "<<kPattern<<endl;
 523
 524     cout<<"INDIR   : "<<kInDir<<endl;
 525     cout<<"NFILES  : "<<kFile<<endl;
 526
 527     char fileE[120] ;
 528     char fileA[120] ;
 529     char fileG[120] ;
 530     char fileEm[120] ;
 531     for (Int_t event = 0 ; event < kFile ; event++) {
 532       sprintf(fileE,  "%s/%s%d/AliESDs.root",    kInDir,kPattern,event) ;
 533       sprintf(fileA,  "%s/%s%d/AliAOD.root",     kInDir,kPattern,event) ;
 534       sprintf(fileG,  "%s/%s%d/galice.root",     kInDir,kPattern,event) ;
 535       sprintf(fileEm, "%s/%s%d/embededAOD.root", kInDir,kPattern,event) ;
 536
 537       TFile * fESD = TFile::Open(fileE) ;
 538       TFile * fAOD = TFile::Open(fileA) ;
 539
 540       //Check if file exists and add it, if not skip it
 541       if (fESD)
 542       {
 543         kTreeName  = "esdTree";
 544         kInputData = "ESD";
 545         TFile * fG = TFile::Open(fileG);
 546         if(fG) { kMC = kTRUE; fG->Close();}
 547         else     kMC = kFALSE;
 548
 549         // Get run number
 550         TTree* esdTree = (TTree*)fESD->Get("esdTree");
 551         AliESDEvent* esd = new AliESDEvent();
 552         esd->ReadFromTree(esdTree);
 553         esdTree->GetEvent(0);
 554         kRun = esd->GetRunNumber();
 555
 556         return;
 557       }
 558       else if(fAOD)
 559       {
 560         kTreeName  = "aodTree";
 561         kInputData = "AOD";
 562         if(((TTree*) fAOD->Get("aodTree"))->GetBranch("mcparticles")) kMC=kTRUE;
 563         else kMC = kFALSE;
 564
 565         // Get run number
 566         TTree* aodTree = (TTree*)fAOD->Get("aodTree");
 567         AliAODEvent* aod = new AliAODEvent();
 568         aod->ReadFromTree(aodTree);
 569         aodTree->GetEvent(0);
 570         kRun = aod->GetRunNumber();
 571         return;
 572       }
 573       else if(TFile::Open(fileEm))
 574       {
 575         kTreeName  = "aodTree";
 576         kInputData = "AOD";
 577         kMC        = kTRUE;
 578
 579         return;
 580       }
 581       else if(TFile::Open(fileG))
 582       {
 583         kTreeName  = "TE";
 584         kInputData = "MC";
 585         kMC        = kTRUE;
 586         return;
 587       }
 588     }
 589
 590     if(fESD) fESD->Close();
 591     if(fAOD) fAOD->Close();
 592
 593   }// local files analysis
 594
 595   //------------------------------
 596   //GRID xml files
 597   //-----------------------------
 598   else if(mode == mGRID){
 599     //Get colection file. It is specified by the environmental
 600     //variable XML
 601
 602     if(gSystem->Getenv("XML") )
 603       kXML = gSystem->Getenv("XML");
 604     else
 605       sprintf(kXML, "collection.xml") ;
 606
 607     if (!TFile::Open(kXML)) {
 608       printf("No collection file with name -- %s -- was found\n",kXML);
 609       return ;
 610     }
 611     else cout<<"XML file "<<kXML<<endl;
 612
 613     //Load necessary libraries and connect to the GRID
 614     gSystem->Load("libNetx.so") ;
 615     gSystem->Load("libRAliEn.so");
 616     TGrid::Connect("alien://") ;
 617
 618     //Feed Grid with collection file
 619     TGridCollection * collection = (TGridCollection*) TAlienCollection::Open(kXML);
 620     if (! collection) {
 621       AliError(Form("%s not found", kXML)) ;
 622       return kFALSE ;
 623     }
 624     TGridResult* result = collection->GetGridResult("",0 ,0);
 625
 626     for (Int_t index = 0; index < result->GetEntries(); index++) {
 627       TString alienURL = result->GetKey(index, "turl") ;
 628       cout << "================== " << alienURL << endl ;
 629
 630       if     (alienURL.Contains("pass1")) kPass = "pass1";
 631       else if(alienURL.Contains("pass2")) kPass = "pass2";
 632       else if(alienURL.Contains("pass3")) kPass = "pass3";
 633
 634       kRun = AliAnalysisManager::GetRunFromAlienPath(alienURL.Data());
 635       printf("Run number from alien path = %d\n",kRun);
 636
 637       TFile * fAOD = 0 ;
 638       //Check if file exists and add it, if not skip it
 639       if (alienURL.Contains("AliESDs.root"))
 640       {
 641         kTreeName  = "esdTree";
 642         kInputData = "ESD";
 643         alienURL.ReplaceAll("AliESDs.root","galice.root");
 644         if(TFile::Open(alienURL)) kMC=kTRUE;
 645         else kMC = kFALSE;
 646         return;
 647       }
 648       else if(alienURL.Contains("AliAOD.root"))
 649       {
 650         kTreeName  = "aodTree";
 651         kInputData = "AOD";
 652         fAOD = TFile::Open(alienURL);
 653         if(((TTree*) fAOD->Get("aodTree"))->GetBranch("mcparticles")) kMC=kTRUE;
 654         else kMC = kFALSE;
 655         return;
 656       }
 657       else if(alienURL.Contains("embededAOD.root"))
 658       {
 659         kTreeName  = "aodTree";
 660         kInputData = "AOD";
 661         kMC=kTRUE;
 662         return;
 663       }
 664       else if(alienURL.Contains("galice.root"))
 665       {
 666         kTreeName  = "TE";
 667         kInputData = "MC";
 668         kMC=kTRUE;
 669         return;
 670       }
 671     }
 672   }// xml analysis
 673   //------------------------------
 674   //PROOF files
 675   //-----------------------------
 676   else if(mode == mPROOF){
 677
 678     TFileCollection* coll  = gProof->GetDataSet(kDatasetPROOF)->GetStagedSubset();
 679
 680     TIter iter(coll->GetList());
 681
 682     TFileInfo* fileInfo = 0;
 683     while ((fileInfo = dynamic_cast<TFileInfo*> (iter())))
 684     {
 685       if (fileInfo->GetFirstUrl()) {
 686         TString ProofURL = fileInfo->GetFirstUrl()->GetUrl();
 687         cout << "================== " << ProofURL << endl ;
 688
 689         if     (ProofURL.Contains("pass1")) kPass = "pass1";
 690         else if(ProofURL.Contains("pass2")) kPass = "pass2";
 691         else if(ProofURL.Contains("pass3")) kPass = "pass3";
 692
 693         kRun = AliAnalysisManager::GetRunFromAlienPath(ProofURL.Data());
 694         printf("Run number from alien path = %d\n",kRun);
 695
 696         TFile * fAOD = 0 ;
 697         //Check if file exists and add it, if not skip it
 698         if (ProofURL.Contains("AliESDs.root"))
 699         {
 700           kTreeName  = "esdTree";
 701           kInputData = "ESD";
 702           alienURL.ReplaceAll("AliESDs.root","galice.root");
 703           if(TFile::Open(ProofURL)) kMC=kTRUE;
 704           else kMC = kFALSE;
 705
 706           return;
 707         }
 708         else if(ProofURL.Contains("AliAOD.root"))
 709         {
 710           kTreeName  = "aodTree";
 711           kInputData = "AOD";
 712           fAOD = TFile::Open(ProofURL);
 713           if(((TTree*) fAOD->Get("aodTree"))->GetBranch("mcparticles")) kMC=kTRUE;
 714           else kMC = kFALSE;
 715           return;
 716         }
 717         else if(ProofURL.Contains("embededAOD.root"))
 718         {
 719           kTreeName  = "aodTree";
 720           kInputData = "AOD";
 721           kMC=kTRUE;
 722           return;
 723         }
 724         else if(ProofURL.Contains("galice.root"))
 725         {
 726           kTreeName  = "TE";
 727           kInputData = "MC";
 728           kMC=kTRUE;
 729           return;
 730         }
 731       }
 732     }
 733   }// proof analysis
 734
 735   gSystem->ChangeDirectory(ocwd.Data());
 736
 737 }
 738
 739 //_____________________________________________________________________
 740 void CreateChain(const anaModes mode, TChain * chain, TChain * chainxs)
 741 {
 742   //Fills chain with data
 743   TString ocwd = gSystem->WorkingDirectory();
 744
 745   //---------------------------------------
 746   // Local files analysis
 747   //---------------------------------------
 748   if(mode == mLocal){
 749     //If you want to add several ESD files sitting in a common directory INDIR
 750     //Specify as environmental variables the directory (INDIR), the number of files
 751     //to analyze (NFILES) and the pattern name of the directories with files (PATTERN)
 752
 753     if(gSystem->Getenv("INDIR"))
 754       kInDir = gSystem->Getenv("INDIR") ;
 755     else cout<<"INDIR not set, use default: "<<kInDir<<endl;
 756
 757     if(gSystem->Getenv("PATTERN"))
 758       kPattern = gSystem->Getenv("PATTERN") ;
 759     else  cout<<"PATTERN not set, use default: "<<kPattern<<endl;
 760
 761     if(gSystem->Getenv("NFILES"))
 762       kFile = atoi(gSystem->Getenv("NFILES")) ;
 763     else cout<<"NFILES not set, use default: "<<kFile<<endl;
 764
 765     //Check if env variables are set and are correct
 766     if ( kInDir  && kFile) {
 767       printf("Get %d files from directory %s\n",kFile,kInDir);
 768       if ( ! gSystem->cd(kInDir) ) {//check if ESDs directory exist
 769         printf("%s does not exist\n", kInDir) ;
 770         return ;
 771       }
 772
 773       //if(gSystem->Getenv("XSFILE"))
 774       //kXSFileName = gSystem->Getenv("XSFILE") ;
 775       //else cout<<" XS file name not set, use default: "<<kXSFileName<<endl;
 776       char * kGener = gSystem->Getenv("GENER");
 777       if(kGener) {
 778         cout<<"GENER "<<kGener<<endl;
 779         if     (!strcmp(kGener,"PYTHIA")) kXSFileName = "pyxsec.root";
 780         else if(!strcmp(kGener,"HERWIG")) kXSFileName = "hexsec.root";
 781         else cout<<" UNKNOWN GENER, use default: "<<kXSFileName<<endl;
 782       }
 783       else cout<<" GENER not set, use default xs file name: "<<kXSFileName<<endl;
 784
 785       cout<<"INDIR   : "<<kInDir     <<endl;
 786       cout<<"NFILES  : "<<kFile      <<endl;
 787       cout<<"PATTERN : "<<kPattern   <<endl;
 788       cout<<"XSFILE  : "<<kXSFileName<<endl;
 789
 790       TString datafile="";
 791       if     (kInputData == "ESD")        datafile = "AliESDs.root" ;
 792       else if(kInputData.Contains("AOD")) datafile = "AliAOD.root"  ;
 793       else if(kInputData == "MC")         datafile = "galice.root"  ;
 794
 795       //Loop on ESD/AOD/MC files, add them to chain
 796       Int_t event =0;
 797       Int_t skipped=0 ;
 798       char file[120] ;
 799       char filexs[120] ;
 800
 801       for (event = 0 ; event < kFile ; event++) {
 802         sprintf(file,   "%s/%s%d/%s", kInDir,kPattern,event,datafile.Data()) ;
 803         sprintf(filexs, "%s/%s%d/%s", kInDir,kPattern,event,kXSFileName) ;
 804         TFile * fData = 0 ;
 805         //Check if file exists and add it, if not skip it
 806         if ( fData = TFile::Open(file)) {
 807           if ( fData->Get(kTreeName) ) {
 808             printf("++++ Adding %s\n", file) ;
 809             chain->AddFile(file);
 810             chainxs->Add(filexs) ;
 811           }
 812         }
 813         else {
 814           printf("---- Skipping %s\n", file) ;
 815           skipped++ ;
 816         }
 817       }
 818     }
 819     else {
 820       TString input = "AliESDs.root" ;
 821       cout<<">>>>>> No list added, take a single file <<<<<<<<< "<<input<<endl;
 822       chain->AddFile(input);
 823     }
 824
 825   }// local files analysis
 826
 827   //------------------------------
 828   // GRID xml files
 829   //------------------------------
 830   else if(mode == mGRID){
 831     //Get colection file. It is specified by the environmental
 832     //variable XML
 833
 834     //Feed Grid with collection file
 835     TGridCollection * collection = (TGridCollection*) TAlienCollection::Open(kXML);
 836     if (! collection) {
 837       AliError(Form("%s not found", kXML)) ;
 838       return kFALSE ;
 839     }
 840
 841     TGridResult* result = collection->GetGridResult("",0 ,0);
 842
 843     // Makes the ESD chain
 844     printf("*** Getting the Chain       ***\n");
 845     for (Int_t index = 0; index < result->GetEntries(); index++) {
 846       TString alienURL = result->GetKey(index, "turl") ;
 847       cout << "================== " << alienURL << endl ;
 848       chain->Add(alienURL) ;
 849       alienURL.ReplaceAll("AliESDs.root",kXSFileName);
 850       chainxs->Add(alienURL) ;
 851     }
 852   }// xml analysis
 853
 854   //------------------------------
 855   // PROOF
 856   //------------------------------
 857   else if (mode == mPROOF) {
 858
 859     TFileCollection* ds= gProof->GetDataSet(kDatasetPROOF)->GetStagedSubset();
 860
 861     gROOT->LoadMacro("/afs/in2p3.fr/group/alice/laf/dataset_management/CreateChainFromDataSet.C");
 862     chain = CreateChainFromDataSet(ds, kTreeName , kDatasetNMaxFiles);
 863     printf("chain has %d entries\n",chain->GetEntries());
 864   }
 865
 866   gSystem->ChangeDirectory(ocwd.Data());
 867
 868 }
 869
 870 //_________________________________________________________________
 871 void GetAverageXsection(TTree * tree, Double_t & xs, Float_t & ntr)
 872 {
 873   // Read the PYTHIA statistics from the file pyxsec.root created by
 874   // the function WriteXsection():
 875   // integrated cross section (xsection) and
 876   // the  number of Pyevent() calls (ntrials)
 877   // and calculate the weight per one event xsection/ntrials
 878   // The spectrum calculated by a user should be
 879   // multiplied by this weight, something like this:
 880   // TH1F *userSpectrum ... // book and fill the spectrum
 881   // userSpectrum->Scale(weight)
 882   //
 883   // Yuri Kharlov 19 June 2007
 884   // Gustavo Conesa 15 April 2008
 885   Double_t xsection = 0;
 886   UInt_t    ntrials = 0;
 887   xs = 0;
 888   ntr = 0;
 889
 890   Int_t nfiles =  tree->GetEntries()  ;
 891   if (tree && nfiles > 0)
 892   {
 893     tree->SetBranchAddress("xsection",&xsection);
 894     tree->SetBranchAddress("ntrials",&ntrials);
 895     for(Int_t i = 0; i < nfiles; i++){
 896       tree->GetEntry(i);
 897       xs += xsection ;
 898       ntr += ntrials ;
 899       cout << "xsection " <<xsection<<" ntrials "<<ntrials<<endl;
 900     }
 901
 902     xs =   xs /  nfiles;
 903     ntr =  ntr / nfiles;
 904     cout << "-----------------------------------------------------------------"<<endl;
 905     cout << "Average of "<< nfiles<<" files: xsection " <<xs<<" ntrials "<<ntr<<endl;
 906     cout << "-----------------------------------------------------------------"<<endl;
 907   }
 908   else cout << " >>>> Empty tree !!!! <<<<< "<<endl;
 909
 910 }
 911
 912 //______________________________
 913 void CheckEnvironmentVariables()
 914 {
 915
 916   sprintf(kTrigger,"");
 917
 918   Bool_t bRecalibrate = kFALSE;
 919   Bool_t bBadChannel = kFALSE;
 920
 921   for (int i=0; i< gApplication->Argc();i++){
 922
 923 #ifdef VERBOSEARGS
 924
 925     printf("Arg  %d:  %s\n",i,gApplication->Argv(i));
 926
 927 #endif
 928
 929     TString sRun = "";
 930
 931     if (!(strcmp(gApplication->Argv(i),"--trigger")))
 932       sprintf(trigger,gApplication->Argv(i+1));
 933
 934     if (!(strcmp(gApplication->Argv(i),"--recalibrate")))
 935       bRecalibrate = atoi(gApplication->Argv(i+1));
 936
 937     if (!(strcmp(gApplication->Argv(i),"--badchannel")))
 938       bBadChannel = atoi(gApplication->Argv(i+1));
 939
 940     if (!(strcmp(gApplication->Argv(i),"--run"))){
 941       sRun = gApplication->Argv(i+1);
 942       if(sRun.Contains("LHC10")) {
 943         kYear = 2010;
 944       }
 945       else {
 946         if(kRun <=0){
 947           kRun = atoi(gApplication->Argv(i+1));
 948         }
 949         else printf("** Run number already set  to %d, do not set to %d\n",kRun,atoi(gApplication->Argv(i+1)));
 950       }//numeric run
 951     }//--run available
 952
 953   }// args loop
 954
 955   if(!sRun.Contains("LHC10")){
 956     if ( kRun < 140000) {
 957       kYear = 2010;
 958       if( kRun >= 136851 ) kCollision = "PbPb";
 959     }
 960     else{
 961       kYear = 2011;
 962     }
 963   }
 964
 965   if(kMC) sprintf(kTrigger,"");
 966
 967   printf("*********************************************\n");
 968   //printf("*** Settings trigger %s, recalibrate %d, remove bad channels %d, year %d, collision %s, run %d ***\n",
 969   //       kTrigger,bRecalibrate,bBadChannel, kYear,kCollision.Data(), kRun);
 970   printf("*** Settings year %d, collision %s, run %d ***\n",kYear,kCollision.Data(), kRun);
 971   printf("*********************************************\n");
 972
 973 }
 974
 975 //_______________________________________________
 976 void AddTaskCounter(const TString trigger = "MB")
 977 {
 978
 979   AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
 980
 981   AliAnalysisTaskCounter * counter =  new AliAnalysisTaskCounter(Form("Counter%s",trigger.Data()));
 982   if(kRun > 140000 && kRun < 146900) counter ->RejectFastCluster();
 983   if     (kCollision=="pp"  )   counter->SetZVertexCut(50.);  //Open cut
 984   else if(kCollision=="PbPb")   counter->SetZVertexCut(10.);  //Centrality defined in this range.
 985
 986   if(trigger=="EMC7")
 987   {
 988     printf("counter trigger EMC7\n");
 989     counter->SelectCollisionCandidates(AliVEvent::kEMC7);
 990   }
 991   else if (trigger=="INT7")
 992   {
 993     printf("counter trigger INT7\n");
 994     counter->SelectCollisionCandidates(AliVEvent::kINT7);
 995   }
 996   if(trigger=="EMC1")
 997   {
 998     printf("counter trigger EMC1\n");
 999     counter->SelectCollisionCandidates(AliVEvent::kEMC1);
1000   }
1001   else if(trigger=="MB")
1002   {
1003     printf("counter trigger MB\n");
1004     counter->SelectCollisionCandidates(AliVEvent::kMB);
1005   }
1006   else if(trigger=="PHOS")
1007   {
1008     printf("counter trigger PHOS\n");
1009     counter->SelectCollisionCandidates(AliVEvent::kPHI7);
1010   }
1011
1012   TString outputFile = AliAnalysisManager::GetCommonFileName();
1013   AliAnalysisDataContainer *cinput1 = mgr->GetCommonInputContainer();
1014
1015   AliAnalysisDataContainer *coutput =
1016   mgr->CreateContainer(Form("Counter%s",trigger.Data()), TList::Class(), AliAnalysisManager::kOutputContainer,  outputFile.Data());
1017   mgr->AddTask(counter);
1018   mgr->ConnectInput  (counter, 0, cinput1);
1019   mgr->ConnectOutput (counter, 1, coutput);
1020
1021 }
1022
1023
1024
1025
1026