PWG4/macros/ana.C

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