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/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("libEMCALUtils");
 218     //gSystem->Load("libPWG4PartCorrBase");
 219     //gSystem->Load("libPWG4PartCorrDep");
 220
 221     //--------------------------------------------------------
 222     //If you want to use root and par files from aliroot
 223     //--------------------------------------------------------
 224     SetupPar("STEERBase");
 225     SetupPar("ESD");
 226     SetupPar("AOD");
 227     SetupPar("ANALYSIS");
 228     SetupPar("ANALYSISalice");
 229     SetupPar("JETAN");
 230         SetupPar("PHOSUtils");
 231         SetupPar("EMCALUtils");
 232 //      //Create Geometry
 233 //    TGeoManager::Import("geometry.root") ; //need file "geometry.root" in local dir!!!!
 234     SetupPar("PWG4PartCorrBase");
 235     SetupPar("PWG4PartCorrDep");
 236
 237   }
 238
 239   //---------------------------------------------------------
 240   // <<<<<<<<<< PROOF mode >>>>>>>>>>>>
 241   //---------------------------------------------------------
 242   else if (mode==mPROOF) {
 243     //
 244     // Connect to proof
 245     // Put appropriate username here
 246     // TProof::Reset("proof://mgheata@lxb6046.cern.ch");
 247     TProof::Open("proof://mgheata@lxb6046.cern.ch");
 248
 249     //    gProof->ClearPackages();
 250     //    gProof->ClearPackage("ESD");
 251     //    gProof->ClearPackage("AOD");
 252     //    gProof->ClearPackage("ANALYSIS");
 253     //    gProof->ClearPackage("JETAN");
 254     //    gProof->ClearPackage("PWG4PartCorrDep");
 255     //    gProof->ClearPackage("PWG4PartCorrBase");
 256
 257     // Enable the STEERBase Package
 258     gProof->UploadPackage("STEERBase.par");
 259     gProof->EnablePackage("STEERBase");
 260     // Enable the ESD Package
 261     gProof->UploadPackage("ESD.par");
 262     gProof->EnablePackage("ESD");
 263     // Enable the AOD Package
 264     gProof->UploadPackage("AOD.par");
 265     gProof->EnablePackage("AOD");
 266     // Enable the Analysis Package
 267     gProof->UploadPackage("ANALYSIS.par");
 268     gProof->EnablePackage("ANALYSIS");
 269     // Enable JETAN
 270     gProof->UploadPackage("JETAN.par");
 271     gProof->EnablePackage("JETAN");
 272     // Enable PHOSUtils
 273     gProof->UploadPackage("PHOSUtils.par");
 274     gProof->EnablePackage("PHOSUtils");
 275         // Enable PHOSUtils
 276         gProof->UploadPackage("EMCALUtils.par");
 277         gProof->EnablePackage("EMCALUtils");
 278     // Enable PartCorr analysis
 279     gProof->UploadPackage("PWG4PartCorrBase.par");
 280     gProof->EnablePackage("PWG4PartCorrBase");
 281     gProof->UploadPackage("PWG4PartCorrDep.par");
 282     gProof->EnablePackage("PWG4PartCorrDep");
 283     //
 284     gProof->ShowEnabledPackages();
 285   }
 286
 287 }
 288
 289 void SetupPar(char* pararchivename)
 290 {
 291   //Load par files, create analysis libraries
 292   //For testing, if par file already decompressed and modified
 293   //classes then do not decompress.
 294
 295   TString cdir(Form("%s", gSystem->WorkingDirectory() )) ;
 296   TString parpar(Form("%s.par", pararchivename)) ;
 297   if ( gSystem->AccessPathName(parpar.Data()) ) {
 298     gSystem->ChangeDirectory(gSystem->Getenv("ALICE_ROOT")) ;
 299     TString processline(Form(".! make %s", parpar.Data())) ;
 300     gROOT->ProcessLine(processline.Data()) ;
 301     gSystem->ChangeDirectory(cdir) ;
 302     processline = Form(".! mv $ALICE_ROOT/%s .", parpar.Data()) ;
 303     gROOT->ProcessLine(processline.Data()) ;
 304   }
 305   if ( gSystem->AccessPathName(pararchivename) ) {
 306     TString processline = Form(".! tar xvzf %s",parpar.Data()) ;
 307     gROOT->ProcessLine(processline.Data());
 308   }
 309
 310   TString ocwd = gSystem->WorkingDirectory();
 311   gSystem->ChangeDirectory(pararchivename);
 312
 313   // check for BUILD.sh and execute
 314   if (!gSystem->AccessPathName("PROOF-INF/BUILD.sh")) {
 315     printf("*******************************\n");
 316     printf("*** Building PAR archive    ***\n");
 317     cout<<pararchivename<<endl;
 318     printf("*******************************\n");
 319
 320     if (gSystem->Exec("PROOF-INF/BUILD.sh")) {
 321       Error("runProcess","Cannot Build the PAR Archive! - Abort!");
 322       return -1;
 323     }
 324   }
 325   // check for SETUP.C and execute
 326   if (!gSystem->AccessPathName("PROOF-INF/SETUP.C")) {
 327     printf("*******************************\n");
 328     printf("*** Setup PAR archive       ***\n");
 329     cout<<pararchivename<<endl;
 330     printf("*******************************\n");
 331     gROOT->Macro("PROOF-INF/SETUP.C");
 332   }
 333
 334   gSystem->ChangeDirectory(ocwd.Data());
 335   printf("Current dir: %s\n", ocwd.Data());
 336 }
 337
 338
 339
 340 void CreateChain(const anaModes mode, TChain * chain, TChain * chainxs){
 341   //Fills chain with data
 342   TString ocwd = gSystem->WorkingDirectory();
 343
 344   //-----------------------------------------------------------
 345   //Analysis of CAF data locally and with PROOF
 346   //-----------------------------------------------------------
 347   if(mode ==mPROOF || mode ==mLocalCAF){
 348     // Chain from CAF
 349     gROOT->LoadMacro("$ALICE_ROOT/PWG0/CreateESDChain.C");
 350     // The second parameter is the number of input files in the chain
 351     chain = CreateESDChain("ESD12001.txt", 5);
 352   }
 353
 354   //---------------------------------------
 355   //Local files analysis
 356   //---------------------------------------
 357   else if(mode == mLocal){
 358     //If you want to add several ESD files sitting in a common directory INDIR
 359     //Specify as environmental variables the directory (INDIR), the number of files
 360     //to analyze (NEVENT) and the pattern name of the directories with files (PATTERN)
 361
 362     if(gSystem->Getenv("INDIR"))
 363       kInDir = gSystem->Getenv("INDIR") ;
 364     else cout<<"INDIR not set, use default: "<<kInDir<<endl;
 365
 366     if(gSystem->Getenv("PATTERN"))
 367       kPattern = gSystem->Getenv("PATTERN") ;
 368     else  cout<<"PATTERN not set, use default: "<<kPattern<<endl;
 369
 370     if(gSystem->Getenv("NEVENT"))
 371       kEvent = atoi(gSystem->Getenv("NEVENT")) ;
 372     else cout<<"NEVENT not set, use default: "<<kEvent<<endl;
 373
 374     //Check if env variables are set and are correct
 375     if ( kInDir  && kEvent) {
 376       printf("Get %d files from directory %s\n",kEvent,kInDir);
 377       if ( ! gSystem->cd(kInDir) ) {//check if ESDs directory exist
 378         printf("%s does not exist\n", kInDir) ;
 379         return ;
 380       }
 381       cout<<"INDIR : "<<kInDir<<endl;
 382       cout<<"NEVENT : "<<kEvent<<endl;
 383       cout<<"PATTERN: " <<kPattern<<endl;
 384
 385           TString datafile="";
 386       if(kInputData == "ESD") datafile = "AliESDs.root" ;
 387       else if(kInputData == "AOD") datafile = "aod.root" ;
 388       else if(kInputData == "MC")  datafile = "galice.root" ;
 389
 390       //Loop on ESD files, add them to chain
 391       Int_t event =0;
 392       Int_t skipped=0 ;
 393       char file[120] ;
 394       char filexs[120] ;
 395
 396       for (event = 0 ; event < kEvent ; event++) {
 397         sprintf(file, "%s/%s%d/%s", kInDir,kPattern,event,datafile.Data()) ;
 398         sprintf(filexs, "%s/%s%d/%s", kInDir,kPattern,event,kXSFileName) ;
 399         TFile * fESD = 0 ;
 400         //Check if file exists and add it, if not skip it
 401         if ( fESD = TFile::Open(file)) {
 402           if ( fESD->Get(kTreeName) ) {
 403             printf("++++ Adding %s\n", file) ;
 404             chain->AddFile(file);
 405             chainxs->Add(filexs) ;
 406           }
 407         }
 408         else {
 409           printf("---- Skipping %s\n", file) ;
 410           skipped++ ;
 411         }
 412       }
 413       printf("number of entries # %lld, skipped %d\n", chain->GetEntries(), skipped) ;
 414     }
 415     else {
 416       TString input = "AliESDs.root" ;
 417       cout<<">>>>>> No list added, take a single file <<<<<<<<< "<<input<<endl;
 418       chain->AddFile(input);
 419     }
 420
 421   }// local files analysis
 422
 423   //------------------------------
 424   //GRID xml files
 425   //-----------------------------
 426   else if(mode == mGRID){
 427     //Get colection file. It is specified by the environmental
 428     //variable XML
 429
 430     if(gSystem->Getenv("XML") )
 431       kXML = gSystem->Getenv("XML");
 432     else
 433       sprintf(kXML, "collection.xml") ;
 434
 435     if (!TFile::Open(kXML)) {
 436       printf("No collection file with name -- %s -- was found\n",kXML);
 437       return ;
 438     }
 439     else cout<<"XML file "<<kXML<<endl;
 440
 441     //Load necessary libraries and connect to the GRID
 442     gSystem->Load("libNetx.so") ;
 443     gSystem->Load("libRAliEn.so");
 444     TGrid::Connect("alien://") ;
 445
 446     //Feed Grid with collection file
 447     //TGridCollection * collection =  (TGridCollection*)gROOT->ProcessLine(Form("TAlienCollection::Open(\"%s\", 0)", kXML));
 448     TGridCollection * collection = (TGridCollection*) TAlienCollection::Open(kXML);
 449     if (! collection) {
 450       AliError(Form("%s not found", kXML)) ;
 451       return kFALSE ;
 452     }
 453     TGridResult* result = collection->GetGridResult("",0 ,0);
 454
 455     // Makes the ESD chain
 456     printf("*** Getting the Chain       ***\n");
 457     for (Int_t index = 0; index < result->GetEntries(); index++) {
 458       TString alienURL = result->GetKey(index, "turl") ;
 459       cout << "================== " << alienURL << endl ;
 460       chain->Add(alienURL) ;
 461       alienURL.ReplaceAll("AliESDs.root",kXSFileName);
 462       chainxs->Add(alienURL) ;
 463     }
 464   }// xml analysis
 465
 466   gSystem->ChangeDirectory(ocwd.Data());
 467 }
 468
 469 //________________________________________________
 470 void GetAverageXsection(TTree * tree, Double_t & xs, Float_t & ntr)
 471 {
 472   // Read the PYTHIA statistics from the file pyxsec.root created by
 473   // the function WriteXsection():
 474   // integrated cross section (xsection) and
 475   // the  number of Pyevent() calls (ntrials)
 476   // and calculate the weight per one event xsection/ntrials
 477   // The spectrum calculated by a user should be
 478   // multiplied by this weight, something like this:
 479   // TH1F *userSpectrum ... // book and fill the spectrum
 480   // userSpectrum->Scale(weight)
 481   //
 482   // Yuri Kharlov 19 June 2007
 483   // Gustavo Conesa 15 April 2008
 484   Double_t xsection = 0;
 485   UInt_t    ntrials = 0;
 486   xs = 0;
 487   ntr = 0;
 488
 489   Int_t nfiles =  tree->GetEntries()  ;
 490   if (tree && nfiles > 0) {
 491     tree->SetBranchAddress("xsection",&xsection);
 492     tree->SetBranchAddress("ntrials",&ntrials);
 493     for(Int_t i = 0; i < nfiles; i++){
 494       tree->GetEntry(i);
 495       xs += xsection ;
 496       ntr += ntrials ;
 497       cout << "xsection " <<xsection<<" ntrials "<<ntrials<<endl;
 498     }
 499
 500     xs =   xs /  nfiles;
 501     ntr =  ntr / nfiles;
 502     cout << "-----------------------------------------------------------------"<<endl;
 503     cout << "Average of "<< nfiles<<" files: xsection " <<xs<<" ntrials "<<ntr<<endl;
 504     cout << "-----------------------------------------------------------------"<<endl;
 505   }
 506   else cout << " >>>> Empty tree !!!! <<<<< "<<endl;
 507
 508 }
 509
 510
 511