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