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
  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 const Bool_t kMC = kFALSE; //With real data kMC = kFALSE
  39 const TString kInputData = "ESD"; //ESD, AOD, MC, deltaAOD
  40 const Bool_t  outAOD = kFALSE; //Some tasks doesnt need it.
  41 TString kTreeName;
  42 const Bool_t kUsePAR = kFALSE; //set to kFALSE for libraries
  43 const Int_t kFilter = kFALSE; //Use ESD filter
  44
  45 void ana(Int_t mode=mLocal)
  46 {
  47   // Main
  48
  49   //--------------------------------------------------------------------
  50   // Load analysis libraries
  51   // Look at the method below,
  52   // change whatever you need for your analysis case
  53   // ------------------------------------------------------------------
  54   LoadLibraries() ;
  55   // TGeoManager::Import("geometry.root") ; //need file "geometry.root" in local dir!!!!
  56
  57   //-------------------------------------------------------------------------------------------------
  58   //Create chain from ESD and from cross sections files, look below for options.
  59   //-------------------------------------------------------------------------------------------------
  60   if(kInputData == "ESD") kTreeName = "esdTree" ;
  61   else if(kInputData.Contains("AOD")) kTreeName = "aodTree" ;
  62   else if (kInputData == "MC") kTreeName = "TE" ;
  63   else {
  64     cout<<"Wrong  data type "<<kInputData<<endl;
  65     break;
  66   }
  67
  68   if(kFilter) outAOD = kTRUE;
  69
  70   TChain *chain       = new TChain(kTreeName) ;
  71   TChain * chainxs = new TChain("Xsection") ;
  72   CreateChain(mode, chain, chainxs);
  73
  74   if(chain){
  75     AliLog::SetGlobalLogLevel(AliLog::kError);//Minimum prints on screen
  76
  77     //--------------------------------------
  78     // Make the analysis manager
  79     //-------------------------------------
  80     AliAnalysisManager *mgr  = new AliAnalysisManager("Manager", "Manager");
  81     //AliAnalysisManager::SetUseProgressBar(kTRUE);
  82     //mgr->SetSkipTerminate(kTRUE);
  83     //mgr->SetNSysInfo(1);
  84
  85     // MC handler
  86     if((kMC || kInputData == "MC") && !kInputData.Contains("AOD")){
  87       AliMCEventHandler* mcHandler = new AliMCEventHandler();
  88       mcHandler->SetReadTR(kFALSE);//Do not search TrackRef file
  89       mgr->SetMCtruthEventHandler(mcHandler);
  90       if( kInputData == "MC") {
  91         cout<<"MC INPUT EVENT HANDLER"<<endl;
  92         mgr->SetInputEventHandler(NULL);
  93       }
  94     }
  95
  96     // AOD output handler
  97     if(kInputData!="deltaAOD" && outAOD){
  98       cout<<"Init output handler"<<endl;
  99       AliAODHandler* aodoutHandler   = new AliAODHandler();
 100       aodoutHandler->SetOutputFileName("aod.root");
 101       ////aodoutHandler->SetCreateNonStandardAOD();
 102       mgr->SetOutputEventHandler(aodoutHandler);
 103     }
 104
 105     //input
 106
 107     if(kInputData == "ESD"){
 108       // ESD handler
 109       AliESDInputHandler *esdHandler = new AliESDInputHandler();
 110       esdHandler->SetReadFriends(kFALSE);
 111       mgr->SetInputEventHandler(esdHandler);
 112       cout<<"ESD handler "<<mgr->GetInputEventHandler()<<endl;
 113     }
 114     else if(kInputData.Contains("AOD")){
 115       // AOD handler
 116       AliAODInputHandler *aodHandler = new AliAODInputHandler();
 117       mgr->SetInputEventHandler(aodHandler);
 118       if(kInputData == "deltaAOD") aodHandler->AddFriend("deltaAODPartCorr.root");
 119       cout<<"AOD handler "<<mgr->GetInputEventHandler()<<endl;
 120     }
 121     //mgr->RegisterExternalFile("deltaAODPartCorr.root");
 122     //mgr->SetDebugLevel(-1); // For debugging, do not uncomment if you want no messages.
 123
 124     //-------------------------------------------------------------------------
 125     //Define task, put here any other task that you want to use.
 126     //-------------------------------------------------------------------------
 127     AliAnalysisDataContainer *cinput1 = mgr->GetCommonInputContainer();
 128     AliAnalysisDataContainer *coutput1;
 129     if(outAOD){
 130
 131       coutput1 = mgr->GetCommonOutputContainer();
 132
 133       if(kInputData=="ESD"){
 134
 135         gROOT->LoadMacro("$ALICE_ROOT/ANALYSIS/macros/AddTaskPhysicsSelection.C");
 136         AliPhysicsSelectionTask* physSelTask = AddTaskPhysicsSelection();
 137
 138         if(kFilter){
 139
 140           // Set of cuts
 141           //for standard global track cuts
 142           AliESDtrackCuts* esdTrackCutsGlobal =  AliESDtrackCuts::GetStandardITSTPCTrackCuts2009(kTRUE);
 143           esdTrackCutsGlobal->SetName("StandardFromAliESDTrackCuts");
 144
 145           //for TPC tracks only
 146           //    AliESDtrackCuts* esdTrackCutsTPC =  AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();
 147           //      esdTrackCutsTPC->SetRequireTPCRefit(kTRUE);
 148           //      esdTrackCutsTPC->SetMinNClustersTPC(70);
 149
 150           AliAnalysisFilter* trackFilter = new AliAnalysisFilter("trackFilter");
 151           trackFilter->AddCuts(esdTrackCutsGlobal);
 152           //trackFilter->AddCuts(esdTrackCutsTPC);
 153
 154           gROOT->LoadMacro("$ALICE_ROOT/ANALYSIS/macros/AddTaskESDFilter.C");
 155           AliAnalysisTaskESDfilter *taskesdfilter = AddTaskESDFilter(kFALSE);
 156           esdfilter->SetTrackFilter(trackFilter);
 157           kInputData = "AOD" ;
 158           kTreeName = "aodTree" ;
 159         }
 160       }
 161     }
 162
 163     gROOT->LoadMacro("$ALICE_ROOT/PWG4/macros/AddTaskPartCorr.C");
 164
 165     AliAnalysisTaskParticleCorrelation *taskEMCAL = AddTaskPartCorr(kInputData,"EMCAL", kFALSE,kFALSE, kFALSE);
 166     //mgr->ProfileTask("PartCorrEMCAL");
 167
 168     //AliAnalysisTaskParticleCorrelation *taskPHOS  = AddTaskPartCorr(kInputData,"PHOS", kFALSE,kFALSE,kFALSE);
 169     //mgr->ProfileTask("PartCorrPHOS");
 170
 171     //gROOT->LoadMacro("$ALICE_ROOT/PWG4/macros/QA/AddTaskCalorimeterQA.C");
 172     //AliAnalysisTaskParticleCorrelation *taskQA = AddTaskCalorimeterQA(kInputData,kFALSE,kFALSE);
 173     //mgr->ProfileTask("CalorimeterPerformance");
 174
 175     //-----------------------
 176     // Run the analysis
 177     //-----------------------
 178     mgr->InitAnalysis();
 179     mgr->PrintStatus();
 180     mgr->StartAnalysis("local",chain);
 181
 182     cout <<" Analysis ended sucessfully "<< endl ;
 183
 184   }
 185   else cout << "Chain was not produced ! "<<endl;
 186   //printf("*** DELETE MANAGER ***\n");
 187   // delete mgr;
 188
 189 }
 190
 191 void  LoadLibraries() {
 192
 193   //--------------------------------------
 194   // Load the needed libraries most of them already loaded by aliroot
 195   //--------------------------------------
 196   gSystem->Load("libTree.so");
 197   gSystem->Load("libGeom.so");
 198   gSystem->Load("libVMC.so");
 199   gSystem->Load("libXMLIO.so");
 200   gSystem->Load("libMatrix.so");
 201   gSystem->Load("libPhysics.so");
 202
 203   if(kUsePAR){
 204     //--------------------------------------------------------
 205     //If you want to use root and par files from aliroot
 206     //--------------------------------------------------------
 207     SetupPar("STEERBase");
 208     SetupPar("ESD");
 209     SetupPar("AOD");
 210     SetupPar("ANALYSIS");
 211     SetupPar("ANALYSISalice");
 212     SetupPar("PHOSUtils");
 213     SetupPar("EMCALUtils");
 214     SetupPar("PWG4PartCorrBase");
 215     SetupPar("PWG4PartCorrDep");
 216   }
 217   else{
 218     //--------------------------------------------------------
 219     // If you want to use already compiled libraries
 220     // in the aliroot distribution
 221     //--------------------------------------------------------
 222     gSystem->Load("libSTEERBase.so");
 223     gSystem->Load("libESD.so");
 224     gSystem->Load("libAOD.so");
 225     gSystem->Load("libANALYSIS.so");
 226     gSystem->Load("libANALYSISalice.so");
 227     gSystem->Load("libPHOSUtils");
 228     gSystem->Load("libEMCALUtils");
 229     gSystem->Load("libPWG4PartCorrBase.so");
 230     gSystem->Load("libPWG4PartCorrDep.so");
 231     if(kFilter){
 232       gSystem->Load("libCORRFW.so");
 233       gSystem->Load("libPWG3base.so");
 234       gSystem->Load("libPWG3muon.so");
 235     }
 236   }
 237 }
 238
 239 void SetupPar(char* pararchivename)
 240 {
 241   //Load par files, create analysis libraries
 242   //For testing, if par file already decompressed and modified
 243   //classes then do not decompress.
 244
 245   TString cdir(Form("%s", gSystem->WorkingDirectory() )) ;
 246   TString parpar(Form("%s.par", pararchivename)) ;
 247 //  if ( gSystem->AccessPathName(parpar.Data()) ) {
 248 //    gSystem->ChangeDirectory(gSystem->Getenv("ALICE_ROOT")) ;
 249 //     TString processline(Form(".! make %s", parpar.Data())) ;
 250 //     gROOT->ProcessLine(processline.Data()) ;
 251 //     gSystem->ChangeDirectory(cdir) ;
 252 //     processline = Form(".! mv $ALICE_ROOT/%s .", parpar.Data()) ;
 253 //     gROOT->ProcessLine(processline.Data()) ;
 254 //   }
 255   if ( gSystem->AccessPathName(pararchivename) ) {
 256     TString processline = Form(".! tar xvzf %s",parpar.Data()) ;
 257     gROOT->ProcessLine(processline.Data());
 258   }
 259
 260   TString ocwd = gSystem->WorkingDirectory();
 261   gSystem->ChangeDirectory(pararchivename);
 262
 263   // check for BUILD.sh and execute
 264   if (!gSystem->AccessPathName("PROOF-INF/BUILD.sh")) {
 265     printf("*******************************\n");
 266     printf("*** Building PAR archive    ***\n");
 267     cout<<pararchivename<<endl;
 268     printf("*******************************\n");
 269
 270     if (gSystem->Exec("PROOF-INF/BUILD.sh")) {
 271       Error("runProcess","Cannot Build the PAR Archive! - Abort!");
 272       return -1;
 273     }
 274   }
 275   // check for SETUP.C and execute
 276   if (!gSystem->AccessPathName("PROOF-INF/SETUP.C")) {
 277     printf("*******************************\n");
 278     printf("*** Setup PAR archive       ***\n");
 279     cout<<pararchivename<<endl;
 280     printf("*******************************\n");
 281     gROOT->Macro("PROOF-INF/SETUP.C");
 282   }
 283
 284   gSystem->ChangeDirectory(ocwd.Data());
 285   printf("Current dir: %s\n", ocwd.Data());
 286 }
 287
 288
 289
 290 void CreateChain(const anaModes mode, TChain * chain, TChain * chainxs){
 291   //Fills chain with data
 292   TString ocwd = gSystem->WorkingDirectory();
 293
 294   //---------------------------------------
 295   //Local files analysis
 296   //---------------------------------------
 297   if(mode == mLocal){
 298     //If you want to add several ESD files sitting in a common directory INDIR
 299     //Specify as environmental variables the directory (INDIR), the number of files
 300     //to analyze (NFILES) and the pattern name of the directories with files (PATTERN)
 301
 302     if(gSystem->Getenv("INDIR"))
 303       kInDir = gSystem->Getenv("INDIR") ;
 304     else cout<<"INDIR not set, use default: "<<kInDir<<endl;
 305
 306     if(gSystem->Getenv("PATTERN"))
 307       kPattern = gSystem->Getenv("PATTERN") ;
 308     else  cout<<"PATTERN not set, use default: "<<kPattern<<endl;
 309
 310     if(gSystem->Getenv("NFILES"))
 311       kFile = atoi(gSystem->Getenv("NFILES")) ;
 312     else cout<<"NFILES not set, use default: "<<kFile<<endl;
 313
 314     //Check if env variables are set and are correct
 315     if ( kInDir  && kFile) {
 316       printf("Get %d files from directory %s\n",kFile,kInDir);
 317       if ( ! gSystem->cd(kInDir) ) {//check if ESDs directory exist
 318         printf("%s does not exist\n", kInDir) ;
 319         return ;
 320       }
 321
 322       //if(gSystem->Getenv("XSFILE"))
 323       //kXSFileName = gSystem->Getenv("XSFILE") ;
 324       //else cout<<" XS file name not set, use default: "<<kXSFileName<<endl;
 325       char * kGener = gSystem->Getenv("GENER");
 326       if(kGener) {
 327         cout<<"GENER "<<kGener<<endl;
 328         if(!strcmp(kGener,"PYTHIA")) kXSFileName = "pyxsec.root";
 329         else if(!strcmp(kGener,"HERWIG")) kXSFileName = "hexsec.root";
 330         else cout<<" UNKNOWN GENER, use default: "<<kXSFileName<<endl;
 331       }
 332       else cout<<" GENER not set, use default xs file name: "<<kXSFileName<<endl;
 333
 334       cout<<"INDIR   : "<<kInDir<<endl;
 335       cout<<"NFILES  : "<<kFile<<endl;
 336       cout<<"PATTERN : " <<kPattern<<endl;
 337       cout<<"XSFILE  : "<<kXSFileName<<endl;
 338
 339       TString datafile="";
 340       if(kInputData == "ESD") datafile = "AliESDs.root" ;
 341       else if(kInputData.Contains("AOD")) datafile = "AliAOD.root" ;
 342       else if(kInputData == "MC")  datafile = "galice.root" ;
 343
 344       //Loop on ESD files, add them to chain
 345       Int_t event =0;
 346       Int_t skipped=0 ;
 347       char file[120] ;
 348       char filexs[120] ;
 349
 350       for (event = 0 ; event < kFile ; event++) {
 351         sprintf(file, "%s/%s%d/%s", kInDir,kPattern,event,datafile.Data()) ;
 352         sprintf(filexs, "%s/%s%d/%s", kInDir,kPattern,event,kXSFileName) ;
 353         TFile * fESD = 0 ;
 354         //Check if file exists and add it, if not skip it
 355         if ( fESD = TFile::Open(file)) {
 356           if ( fESD->Get(kTreeName) ) {
 357             printf("++++ Adding %s\n", file) ;
 358             chain->AddFile(file);
 359             chainxs->Add(filexs) ;
 360           }
 361         }
 362         else {
 363           printf("---- Skipping %s\n", file) ;
 364           skipped++ ;
 365         }
 366       }
 367       printf("number of entries # %lld, skipped %d\n", chain->GetEntries(), skipped*100) ;
 368     }
 369     else {
 370       TString input = "AliESDs.root" ;
 371       cout<<">>>>>> No list added, take a single file <<<<<<<<< "<<input<<endl;
 372       chain->AddFile(input);
 373     }
 374
 375   }// local files analysis
 376
 377   //------------------------------
 378   //GRID xml files
 379   //-----------------------------
 380   else if(mode == mGRID){
 381     //Get colection file. It is specified by the environmental
 382     //variable XML
 383
 384     if(gSystem->Getenv("XML") )
 385       kXML = gSystem->Getenv("XML");
 386     else
 387       sprintf(kXML, "collection.xml") ;
 388
 389     if (!TFile::Open(kXML)) {
 390       printf("No collection file with name -- %s -- was found\n",kXML);
 391       return ;
 392     }
 393     else cout<<"XML file "<<kXML<<endl;
 394
 395     //Load necessary libraries and connect to the GRID
 396     gSystem->Load("libNetx.so") ;
 397     gSystem->Load("libRAliEn.so");
 398     TGrid::Connect("alien://") ;
 399
 400     //Feed Grid with collection file
 401     //TGridCollection * collection =  (TGridCollection*)gROOT->ProcessLine(Form("TAlienCollection::Open(\"%s\", 0)", kXML));
 402     TGridCollection * collection = (TGridCollection*) TAlienCollection::Open(kXML);
 403     if (! collection) {
 404       AliError(Form("%s not found", kXML)) ;
 405       return kFALSE ;
 406     }
 407     TGridResult* result = collection->GetGridResult("",0 ,0);
 408
 409     // Makes the ESD chain
 410     printf("*** Getting the Chain       ***\n");
 411     for (Int_t index = 0; index < result->GetEntries(); index++) {
 412       TString alienURL = result->GetKey(index, "turl") ;
 413       cout << "================== " << alienURL << endl ;
 414       chain->Add(alienURL) ;
 415       alienURL.ReplaceAll("AliESDs.root",kXSFileName);
 416       chainxs->Add(alienURL) ;
 417     }
 418   }// xml analysis
 419
 420   gSystem->ChangeDirectory(ocwd.Data());
 421 }
 422
 423 //________________________________________________
 424 void GetAverageXsection(TTree * tree, Double_t & xs, Float_t & ntr)
 425 {
 426   // Read the PYTHIA statistics from the file pyxsec.root created by
 427   // the function WriteXsection():
 428   // integrated cross section (xsection) and
 429   // the  number of Pyevent() calls (ntrials)
 430   // and calculate the weight per one event xsection/ntrials
 431   // The spectrum calculated by a user should be
 432   // multiplied by this weight, something like this:
 433   // TH1F *userSpectrum ... // book and fill the spectrum
 434   // userSpectrum->Scale(weight)
 435   //
 436   // Yuri Kharlov 19 June 2007
 437   // Gustavo Conesa 15 April 2008
 438   Double_t xsection = 0;
 439   UInt_t    ntrials = 0;
 440   xs = 0;
 441   ntr = 0;
 442
 443   Int_t nfiles =  tree->GetEntries()  ;
 444   if (tree && nfiles > 0) {
 445     tree->SetBranchAddress("xsection",&xsection);
 446     tree->SetBranchAddress("ntrials",&ntrials);
 447     for(Int_t i = 0; i < nfiles; i++){
 448       tree->GetEntry(i);
 449       xs += xsection ;
 450       ntr += ntrials ;
 451       cout << "xsection " <<xsection<<" ntrials "<<ntrials<<endl;
 452     }
 453
 454     xs =   xs /  nfiles;
 455     ntr =  ntr / nfiles;
 456     cout << "-----------------------------------------------------------------"<<endl;
 457     cout << "Average of "<< nfiles<<" files: xsection " <<xs<<" ntrials "<<ntr<<endl;
 458     cout << "-----------------------------------------------------------------"<<endl;
 459   }
 460   else cout << " >>>> Empty tree !!!! <<<<< "<<endl;
 461
 462 }
 463
 464
 465