PWG4/macros/anaGammaAnalysis.C

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