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