PWG4/macros/electrons/anaElectron.C

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