PWG4/macros/anaExampleTask.C

   1 /* $Id: anaGammaAnalysis.C 25095 2008-04-11 12:54:47Z schutz $ */
   2 /* $Log$
   3 /* Revision 1.2  2007/12/13 09:45:45  gustavo
   4 /* Scaling option and more comentaries added
   5 /*
   6 /* Revision 1.1  2007/12/07 14:13:02  gustavo
   7 /* Example macros for execution and configuration of the analysis
   8 /* */
   9
  10 //---------------------------------------------------
  11 // Example macro to do analysis with the
  12 // AliAnalysisTaksSE
  13 // Can be executed with Root and AliRoot
  14 //
  15 // Pay attention to the options and definitions
  16 // set in the lines below
  17 //
  18 //  Author : Gustavo Conesa Balbastre (INFN-LNF)
  19 //
  20 //-------------------------------------------------
  21 enum anaModes {mLocal, mLocalCAF,mPROOF,mGRID};
  22 //mLocal: Analyze locally files in your computer
  23 //mLocalCAF: Analyze locally CAF files
  24 //mPROOF: Analyze CAF files with PROOF
  25
  26 //---------------------------------------------------------------------------
  27 //Settings to read locally several files, only for "mLocal" mode
  28 //The different values are default, they can be set with environmental
  29 //variables: INDIR, PATTERN, NEVENT, respectivelly
  30 char * kInDir = "/user/data";
  31 char * kPattern = ""; // Data are in diles /data/Run0,
  32 // /Data/Run1 ...
  33 Int_t kEvent = 1; // Number of files
  34 //---------------------------------------------------------------------------
  35 //Collection file for grid analysis
  36 char * kXML = "collection.xml";
  37 //---------------------------------------------------------------------------
  38 //Scale histograms from file. Change to kTRUE when xsection file exists
  39 //Put name of file containing xsection
  40 //Put number of events per ESD file
  41 //This is an specific case for normalization of Pythia files.
  42 const Bool_t kGetXSectionFromFileAndScale = kFALSE ;
  43 const char * kXSFileName = "pyxsec.root";
  44 const Int_t kNumberOfEventsPerFile = 100;
  45 //---------------------------------------------------------------------------
  46
  47 const Bool_t kMC = kTRUE; //With real data kMC = kFALSE
  48 const TString kInputData = "ESD";
  49
  50 void anaExampleTask(Int_t mode=mLocal)
  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   LoadLibraries(mode) ;
  60
  61   //-------------------------------------------------------------------------------------------------
  62   //Create chain from ESD and from cross sections files, look below for options.
  63   //-------------------------------------------------------------------------------------------------
  64   TChain *chain   = new TChain("esdTree") ;
  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){
  77       AliMCEventHandler* mcHandler = new AliMCEventHandler();
  78       mcHandler->SetReadTR(kFALSE);//Do not search TrackRef file
  79       mgr->SetMCtruthEventHandler(mcHandler);
  80     }
  81
  82     // AOD output handler
  83     AliAODHandler* aodoutHandler   = new AliAODHandler();
  84     aodoutHandler->SetOutputFileName("aod.root");
  85     //aodoutHandler->SetCreateNonStandardAOD();
  86     mgr->SetOutputEventHandler(aodoutHandler);
  87
  88     //input
  89     if(kInputData == "ESD"){
  90       // ESD handler
  91       AliESDInputHandler *esdHandler = new AliESDInputHandler();
  92       mgr->SetInputEventHandler(esdHandler);
  93     }
  94     if(kInputData == "AOD"){
  95       // AOD handler
  96       AliAODInputHandler *aodHandler = new AliAODInputHandler();
  97       mgr->SetInputEventHandler(aodHandler);
  98     }
  99
 100     //mgr->SetDebugLevel(10); // For debugging
 101
 102     //-------------------------------------------------------------------------
 103     //Define task, put here any other task that you want to use.
 104     //-------------------------------------------------------------------------
 105     AliAnalysisTaskPHOSExample * taskex = new AliAnalysisTaskPHOSExample ("Gamma");
 106     mgr->AddTask(taskex);
 107
 108     // Create containers for input/output
 109     AliAnalysisDataContainer *cinput1 = mgr->GetCommonInputContainer();
 110     AliAnalysisDataContainer *coutput1 = mgr->GetCommonOutputContainer();
 111     AliAnalysisDataContainer *coutput2 = mgr->CreateContainer("histos", TList::Class(),
 112                                                               AliAnalysisManager::kOutputContainer, "histos.root");
 113
 114     mgr->ConnectInput  (taskex,     0, cinput1);
 115     mgr->ConnectOutput (taskex,     0, coutput1 );
 116     mgr->ConnectOutput (taskex,     1, coutput2 );
 117
 118     //------------------------
 119     //Scaling task
 120     //-----------------------
 121     Int_t nfiles = chainxs->GetEntries();
 122     if(kGetXSectionFromFileAndScale && nfiles > 0){
 123       //Get the cross section
 124       Double_t xsection=0;
 125       Float_t ntrials = 0;
 126       GetAverageXsection(chainxs, xsection, ntrials);
 127
 128       AliAnaScale * scale = new AliAnaScale("scale") ;
 129       scale->Set(xsection/ntrials/kNumberOfEventsPerFile/nfiles) ;
 130       mgr->AddTask(scale);
 131
 132       AliAnalysisDataContainer *coutput3 = mgr->CreateContainer("histosscaled", TList::Class(),
 133                                                                 AliAnalysisManager::kOutputContainer, "gammahistosscaled.root");
 134       mgr->ConnectInput  (scale,     0, coutput2);
 135       mgr->ConnectOutput (scale,     0, coutput3 );
 136     }
 137
 138     //-----------------------
 139     // Run the analysis
 140     //-----------------------
 141     TString smode = "";
 142     if (mode==mLocal || mode == mLocalCAF)
 143       smode = "local";
 144     else if (mode==mPROOF)
 145       smode = "proof";
 146     else if (mode==mGRID)
 147       smode = "local";
 148
 149     mgr->InitAnalysis();
 150     mgr->PrintStatus();
 151     mgr->StartAnalysis(smode.Data(),chain);
 152   }
 153   else cout << "Chain was not produced ! "<<endl;
 154
 155 }
 156
 157 void  LoadLibraries(const anaModes mode) {
 158
 159   //--------------------------------------
 160   // Load the needed libraries most of them already loaded by aliroot
 161   //--------------------------------------
 162   gSystem->Load("libTree.so");
 163   gSystem->Load("libGeom.so");
 164   gSystem->Load("libVMC.so");
 165   gSystem->Load("libXMLIO.so");
 166
 167   //----------------------------------------------------------
 168   // >>>>>>>>>>> Local mode <<<<<<<<<<<<<<
 169   //----------------------------------------------------------
 170   if (mode==mLocal || mode == mLocalCAF || mode == mGRID) {
 171     //--------------------------------------------------------
 172     // If you want to use already compiled libraries
 173     // in the aliroot distribution
 174     //--------------------------------------------------------
 175     //gSystem->Load("libSTEERBase");
 176     //gSystem->Load("libESD");
 177     //gSystem->Load("libAOD");
 178     //gSystem->Load("libANALYSIS");
 179     //gSystem->Load("libANALYSISalice");
 180     //gSystem->Load("libPWG4PartCorrBase");
 181      //gSystem->Load("libPWG4PartCorrDep");
 182
 183     //--------------------------------------------------------
 184     //If you want to use root and par files from aliroot
 185     //--------------------------------------------------------
 186     SetupPar("STEERBase");
 187     SetupPar("ESD");
 188     SetupPar("AOD");
 189     SetupPar("ANALYSIS");
 190     SetupPar("ANALYSISalice");
 191     SetupPar("PWG4PartCorrBase");
 192     SetupPar("PWG4PartCorrDep");
 193   }
 194
 195   //---------------------------------------------------------
 196   // <<<<<<<<<< PROOF mode >>>>>>>>>>>>
 197   //---------------------------------------------------------
 198   else if (mode==mPROOF) {
 199     //
 200     // Connect to proof
 201     // Put appropriate username here
 202     // TProof::Reset("proof://mgheata@lxb6046.cern.ch");
 203     TProof::Open("proof://mgheata@lxb6046.cern.ch");
 204
 205     // Enable the STEERBase Package
 206     gProof->UploadPackage("STEERBase.par");
 207     gProof->EnablePackage("STEERBase");
 208     // Enable the ESD Package
 209     gProof->UploadPackage("ESD.par");
 210     gProof->EnablePackage("ESD");
 211     // Enable the AOD Package
 212     gProof->UploadPackage("AOD.par");
 213     gProof->EnablePackage("AOD");
 214     // Enable the Analysis Package
 215     gProof->UploadPackage("ANALYSIS.par");
 216     gProof->EnablePackage("ANALYSIS");
 217     // Enable PartCorr analysis
 218     gProof->UploadPackage("PWG4PartCorrBase.par");
 219     gProof->EnablePackage("PWG4PartCorrBase");
 220     gProof->UploadPackage("PWG4PartCorrDep.par");
 221     gProof->EnablePackage("PWG4PartCorrDep");
 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