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