PWG4/macros/ana.C

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