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