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