PWG4/macros/anaPartCorrJetAn.C

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