remove analysis classes new implemention
[u/mrichter/AliRoot.git] / PWG4 / macros / anaGammaAnalysis.C
CommitLineData
34ce8d31 1/* $Id$ */
37ff46b0 2/* $Log$
030af48c 3/* Revision 1.2 2007/12/13 09:45:45 gustavo
4/* Scaling option and more comentaries added
5/*
37ff46b0 6/* Revision 1.1 2007/12/07 14:13:02 gustavo
7/* Example macros for execution and configuration of the analysis
8/* */
34ce8d31 9
10//---------------------------------------------------
11// Example macro to do analysis with the
12// analysis classes in PWG4Gamma
13// Can be executed with Root and AliRoot
37ff46b0 14//
ce8dba2f 15// Pay attention to the options and definitions
16// set in the lines below
34ce8d31 17//
18// Author : Gustavo Conesa Balbastre (INFN-LNF)
19//
20//-------------------------------------------------
34ce8d31 21enum anaModes {mLocal, mLocalCAF,mPROOF,mGRID};
22//mLocal: Analyze locally files in your computer
23//mLocalCAF: Analyze locally CAF files
24//mPROOF: Analyze CAF files with PROOF
37ff46b0 25
ce8dba2f 26//---------------------------------------------------------------------------
27//Settings to read locally several files, only for "mLocal" mode
28//The different values are default, they can be set with environmental
d6072244 29//variables: INDIR, PATTERN, NEVENT, respectivelly
f78666e7 30char * kInDir = "/home/group/alice/schutz/analysis/PWG4/data";
31char * kPattern = ""; // Data are in diles /data/Run0,
030af48c 32// /Data/Run1 ...
f78666e7 33Int_t kEvent = 1; // Number of files
ce8dba2f 34//---------------------------------------------------------------------------
35//Collection file for grid analysis
36char * kXML = "collection.xml";
37//---------------------------------------------------------------------------
37ff46b0 38//Scale histograms from file. Change to kTRUE when xsection file exists
39//Put name of file containing xsection
40//Put number of events per ESD file
41//This is an specific case for normalization of Pythia files.
d6072244 42const Bool_t kGetXSectionFromFileAndScale = kTRUE ;
37ff46b0 43const char * kXSFileName = "pyxsec.root";
44const Int_t kNumberOfEventsPerFile = 100;
ce8dba2f 45//---------------------------------------------------------------------------
37ff46b0 46
ce8dba2f 47const Bool_t kMC = kTRUE; //With real data kMC = kFALSE
34ce8d31 48
d6072244 49void anaGammaAnalysis(Int_t mode=mLocal, TString configName = "ConfigKineGammaDirect")
34ce8d31 50{
51 // Main
ce8dba2f 52
53 //--------------------------------------------------------------------
34ce8d31 54 // Load analysis libraries
55 // Look at the method below,
56 // change whatever you need for your analysis case
ce8dba2f 57 // ------------------------------------------------------------------
34ce8d31 58 LoadLibraries(mode) ;
59
ce8dba2f 60 //-------------------------------------------------------------------------------------------------
61 //Create chain from ESD and from cross sections files, look below for options.
62 //-------------------------------------------------------------------------------------------------
63 TChain *chain = new TChain("esdTree") ;
64 TChain * chainxs = new TChain("Xsection") ;
65 CreateChain(mode, chain, chainxs);
37ff46b0 66
34ce8d31 67 if(chain){
ce8dba2f 68 AliLog::SetGlobalLogLevel(AliLog::kError);//Minimum prints on screen
37ff46b0 69
ce8dba2f 70 //--------------------------------------
71 // Make the analysis manager
72 //-------------------------------------
73 AliAnalysisManager *mgr = new AliAnalysisManager("Manager", "Manager");
74 // MC handler
75 if(kMC){
76 AliMCEventHandler* mcHandler = new AliMCEventHandler();
77 mgr->SetMCtruthEventHandler(mcHandler);
78 }
79 // AOD output handler
80 AliAODHandler* aodHandler = new AliAODHandler();
81 aodHandler->SetOutputFileName("aod.root");
82 mgr->SetOutputEventHandler(aodHandler);
83 // ESD handler
84 AliESDInputHandler *esdHandler = new AliESDInputHandler();
85 mgr->SetInputEventHandler(esdHandler);
86
87 //mgr->SetDebugLevel(10); // For debugging
88
89 //-------------------------------------------------------------------------
90 //Define task, put here any other task that you want to use.
91 //-------------------------------------------------------------------------
92 AliAnalysisTaskGamma * taskgamma = new AliAnalysisTaskGamma ("Gamma");
93 taskgamma->SetConfigFileName(configName); //Default name is ConfigGammaAnalysis
94 mgr->AddTask(taskgamma);
95
96 // Create containers for input/output
97 AliAnalysisDataContainer *cinput1 = mgr->CreateContainer("cchain",TChain::Class(),
98 AliAnalysisManager::kInputContainer);
99 AliAnalysisDataContainer *coutput1 = mgr->CreateContainer("tree", TTree::Class(),
100 AliAnalysisManager::kOutputContainer, "default");
101 AliAnalysisDataContainer *coutput2 = mgr->CreateContainer("gammahistos", TList::Class(),
102 AliAnalysisManager::kOutputContainer, "gammahistos.root");
103
104 mgr->ConnectInput (taskgamma, 0, cinput1);
105 mgr->ConnectOutput (taskgamma, 0, coutput1 );
106 mgr->ConnectOutput (taskgamma, 1, coutput2 );
34ce8d31 107
ce8dba2f 108 //------------------------
109 //Scaling task
110 //-----------------------
111 Int_t nfiles = chainxs->GetEntries();
112 if(kGetXSectionFromFileAndScale && nfiles > 0){
113 //Get the cross section
114 Double_t xsection=0;
115 Float_t ntrials = 0;
116 GetAverageXsection(chainxs, xsection, ntrials);
117
118 AliAnaScale * scale = new AliAnaScale("scale") ;
119 scale->Set(xsection/ntrials/kNumberOfEventsPerFile/nfiles) ;
120 mgr->AddTask(scale);
121
122 AliAnalysisDataContainer *coutput3 = mgr->CreateContainer("gammahistosscaled", TList::Class(),
123 AliAnalysisManager::kOutputContainer, "gammahistosscaled.root");
124 mgr->ConnectInput (scale, 0, coutput2);
125 mgr->ConnectOutput (scale, 0, coutput3 );
126 }
127
128 //-----------------------
129 // Run the analysis
130 //-----------------------
131 TString smode = "";
132 if (mode==mLocal || mode == mLocalCAF)
133 smode = "local";
134 else if (mode==mPROOF)
135 smode = "proof";
136 else if (mode==mGRID)
137 smode = "grid";
138
139 mgr->InitAnalysis();
140 mgr->PrintStatus();
141 mgr->StartAnalysis(smode.Data(),chain);
34ce8d31 142 }
143 else cout << "Chain was not produced ! "<<endl;
ce8dba2f 144
34ce8d31 145}
146
147void LoadLibraries(const anaModes mode) {
148
149 //--------------------------------------
37ff46b0 150 // Load the needed libraries most of them already loaded by aliroot
34ce8d31 151 //--------------------------------------
152 gSystem->Load("libTree.so");
153 gSystem->Load("libGeom.so");
154 gSystem->Load("libVMC.so");
155 gSystem->Load("libXMLIO.so");
156
ce8dba2f 157 //----------------------------------------------------------
34ce8d31 158 // >>>>>>>>>>> Local mode <<<<<<<<<<<<<<
ce8dba2f 159 //----------------------------------------------------------
34ce8d31 160 if (mode==mLocal || mode == mLocalCAF || mode == mGRID) {
37ff46b0 161 //--------------------------------------------------------
162 // If you want to use already compiled libraries
163 // in the aliroot distribution
164 //--------------------------------------------------------
165 //gSystem->Load("libSTEERBase");
166 //gSystem->Load("libESD");
167 //gSystem->Load("libAOD");
168 //gSystem->Load("libANALYSIS");
d6072244 169 //gSystem->Load("libANALYSISalice");
37ff46b0 170 //gSystem->Load("libPWG4Gamma");
ce8dba2f 171
37ff46b0 172 //--------------------------------------------------------
d6072244 173 //If you want to use root and par files from aliroot
37ff46b0 174 //--------------------------------------------------------
175 SetupPar("STEERBase");
176 SetupPar("ESD");
177 SetupPar("AOD");
178 SetupPar("ANALYSIS");
d6072244 179 SetupPar("ANALYSISalice");
37ff46b0 180 SetupPar("PWG4Gamma");
181
34ce8d31 182 }
ce8dba2f 183
184 //---------------------------------------------------------
34ce8d31 185 // <<<<<<<<<< PROOF mode >>>>>>>>>>>>
ce8dba2f 186 //---------------------------------------------------------
34ce8d31 187 else if (mode==mPROOF) {
188 //
189 // Connect to proof
190 // Put appropriate username here
191 // TProof::Reset("proof://mgheata@lxb6046.cern.ch");
192 TProof::Open("proof://mgheata@lxb6046.cern.ch");
193
194 // gProof->ClearPackages();
195 // gProof->ClearPackage("ESD");
196 // gProof->ClearPackage("AOD");
197 // gProof->ClearPackage("ANALYSIS");
198 // gProof->ClearPackage("PWG4Gamma");
199
200 // Enable the STEERBase Package
201 gProof->UploadPackage("STEERBase.par");
202 gProof->EnablePackage("STEERBase");
203 // Enable the ESD Package
204 gProof->UploadPackage("ESD.par");
205 gProof->EnablePackage("ESD");
206 // Enable the AOD Package
207 gProof->UploadPackage("AOD.par");
208 gProof->EnablePackage("AOD");
209 // Enable the Analysis Package
210 gProof->UploadPackage("ANALYSIS.par");
211 gProof->EnablePackage("ANALYSIS");
212 // Enable gamma jet analysis
213 gProof->UploadPackage("PWG4Gamma.par");
214 gProof->EnablePackage("PWG4Gamma");
215 //
216 gProof->ShowEnabledPackages();
217 }
218
219}
220
d6072244 221void SetupPar(char* pararchivename)
34ce8d31 222{
34ce8d31 223 //Load par files, create analysis libraries
224 //For testing, if par file already decompressed and modified
225 //classes then do not decompress.
ce8dba2f 226
f78666e7 227 TString cdir(Form("%s", gSystem->WorkingDirectory() )) ;
228 TString parpar(Form("%s.par", pararchivename)) ;
229 if ( gSystem->AccessPathName(parpar.Data()) ) {
230 gSystem->ChangeDirectory(gSystem->Getenv("ALICE_ROOT")) ;
231 TString processline(Form(".! make %s", parpar.Data())) ;
232 gROOT->ProcessLine(processline.Data()) ;
233 gSystem->ChangeDirectory(cdir) ;
234 processline = Form(".! mv /tmp/%s .", parpar.Data()) ;
235 gROOT->ProcessLine(processline.Data()) ;
236 }
237 if ( gSystem->AccessPathName(pararchivename) ) {
238 TString processline = Form(".! tar xvzf %s",parpar.Data()) ;
239 gROOT->ProcessLine(processline.Data());
240 }
ce8dba2f 241
f78666e7 242 TString ocwd = gSystem->WorkingDirectory();
243 gSystem->ChangeDirectory(pararchivename);
244
245 // check for BUILD.sh and execute
246 if (!gSystem->AccessPathName("PROOF-INF/BUILD.sh")) {
247 printf("*******************************\n");
248 printf("*** Building PAR archive ***\n");
d6072244 249 cout<<pararchivename<<endl;
f78666e7 250 printf("*******************************\n");
34ce8d31 251
f78666e7 252 if (gSystem->Exec("PROOF-INF/BUILD.sh")) {
253 Error("runProcess","Cannot Build the PAR Archive! - Abort!");
254 return -1;
34ce8d31 255 }
34ce8d31 256 }
f78666e7 257 // check for SETUP.C and execute
258 if (!gSystem->AccessPathName("PROOF-INF/SETUP.C")) {
259 printf("*******************************\n");
260 printf("*** Setup PAR archive ***\n");
d6072244 261 cout<<pararchivename<<endl;
f78666e7 262 printf("*******************************\n");
263 gROOT->Macro("PROOF-INF/SETUP.C");
264 }
265
266 gSystem->ChangeDirectory(ocwd.Data());
267 printf("Current dir: %s\n", ocwd.Data());
34ce8d31 268}
269
270
f78666e7 271
ce8dba2f 272void CreateChain(const anaModes mode, TChain * chain, TChain * chainxs){
273 //Fills chain with data
34ce8d31 274 TString ocwd = gSystem->WorkingDirectory();
ce8dba2f 275
276 //-----------------------------------------------------------
34ce8d31 277 //Analysis of CAF data locally and with PROOF
ce8dba2f 278 //-----------------------------------------------------------
34ce8d31 279 if(mode ==mPROOF || mode ==mLocalCAF){
280 // Chain from CAF
281 gROOT->LoadMacro("$ALICE_ROOT/PWG0/CreateESDChain.C");
282 // The second parameter is the number of input files in the chain
283 chain = CreateESDChain("ESD12001.txt", 5);
284 }
ce8dba2f 285
286 //---------------------------------------
34ce8d31 287 //Local files analysis
ce8dba2f 288 //---------------------------------------
289 else if(mode == mLocal){
d6072244 290 //If you want to add several ESD files sitting in a common directory INDIR
291 //Specify as environmental variables the directory (INDIR), the number of files
34ce8d31 292 //to analyze (NEVENT) and the pattern name of the directories with files (PATTERN)
ce8dba2f 293
d6072244 294 if(gSystem->Getenv("INDIR"))
295 kInDir = gSystem->Getenv("INDIR") ;
296 else cout<<"INDIR not set, use default: "<<kInDir<<endl;
ce8dba2f 297
030af48c 298 if(gSystem->Getenv("PATTERN"))
299 kPattern = gSystem->Getenv("PATTERN") ;
300 else cout<<"PATTERN not set, use default: "<<kPattern<<endl;
ce8dba2f 301
37ff46b0 302 if(gSystem->Getenv("NEVENT"))
303 kEvent = atoi(gSystem->Getenv("NEVENT")) ;
030af48c 304 else cout<<"NEVENT not set, use default: "<<kEvent<<endl;
37ff46b0 305
34ce8d31 306 //Check if env variables are set and are correct
307 if ( kInDir && kEvent) {
308 printf("Get %d files from directory %s\n",kEvent,kInDir);
309 if ( ! gSystem->cd(kInDir) ) {//check if ESDs directory exist
310 printf("%s does not exist\n", kInDir) ;
311 return ;
312 }
d6072244 313 cout<<"INDIR : "<<kInDir<<endl;
34ce8d31 314 cout<<"NEVENT : "<<kEvent<<endl;
315 cout<<"PATTERN: " <<kPattern<<endl;
37ff46b0 316
34ce8d31 317 //Loop on ESD files, add them to chain
37ff46b0 318 Int_t event =0;
319 Int_t skipped=0 ;
34ce8d31 320 char file[120] ;
ce8dba2f 321 char filexs[120] ;
322
34ce8d31 323 for (event = 0 ; event < kEvent ; event++) {
ce8dba2f 324 sprintf(file, "%s/%s%d/AliESDs.root", kInDir,kPattern,event) ;
325 sprintf(filexs, "%s/%s%d/%s", kInDir,kPattern,event,kXSFileName) ;
34ce8d31 326 TFile * fESD = 0 ;
327 //Check if file exists and add it, if not skip it
030af48c 328 if ( fESD = TFile::Open(file)) {
34ce8d31 329 if ( fESD->Get("esdTree") ) {
37ff46b0 330 printf("++++ Adding %s\n", file) ;
331 chain->AddFile(file);
ce8dba2f 332 chainxs->Add(filexs) ;
34ce8d31 333 }
ce8dba2f 334 }
335 else {
336 printf("---- Skipping %s\n", file) ;
337 skipped++ ;
338 }
34ce8d31 339 }
340 printf("number of entries # %lld, skipped %d\n", chain->GetEntries(), skipped*100) ;
341 }
342 else {
ce8dba2f 343 TString input = "AliESDs.root" ;
34ce8d31 344 cout<<">>>>>> No list added, take a single file <<<<<<<<< "<<input<<endl;
345 chain->AddFile(input);
346 }
ce8dba2f 347
34ce8d31 348 }// local files analysis
ce8dba2f 349
350 //------------------------------
34ce8d31 351 //GRID xml files
ce8dba2f 352 //-----------------------------
34ce8d31 353 else if(mode == mGRID){
ce8dba2f 354 //Get colection file. It is specified by the environmental
355 //variable XML
356
34ce8d31 357 if(gSystem->Getenv("XML") )
358 kXML = gSystem->Getenv("XML");
359 else
ce8dba2f 360 sprintf(kXML, "collection.xml") ;
34ce8d31 361
362 if (!TFile::Open(kXML)) {
363 printf("No collection file with name -- %s -- was found\n",kXML);
364 return ;
365 }
ce8dba2f 366 else cout<<"XML file "<<kXML<<endl;
34ce8d31 367
ce8dba2f 368 //Load necessary libraries and connect to the GRID
f78666e7 369 gSystem->Load("libNetx.so") ;
370 gSystem->Load("libRAliEn.so");
371 TGrid::Connect("alien://") ;
ce8dba2f 372
373 //Feed Grid with collection file
374 //TGridCollection * collection = (TGridCollection*)gROOT->ProcessLine(Form("TAlienCollection::Open(\"%s\", 0)", kXML));
375 TGridCollection * collection = (TGridCollection*) TAlienCollection::Open(kXML);
34ce8d31 376 if (! collection) {
377 AliError(Form("%s not found", kXML)) ;
378 return kFALSE ;
379 }
34ce8d31 380 TGridResult* result = collection->GetGridResult("",0 ,0);
ce8dba2f 381
34ce8d31 382 // Makes the ESD chain
383 printf("*** Getting the Chain ***\n");
f78666e7 384 for (Int_t index = 0; index < result->GetEntries(); index++) {
385 TString alienURL = result->GetKey(index, "turl") ;
386 cout << "================== " << alienURL << endl ;
387 chain->Add(alienURL) ;
ce8dba2f 388 alienURL.ReplaceAll("AliESDs.root",kXSFileName);
389 chainxs->Add(alienURL) ;
f78666e7 390 }
34ce8d31 391 }// xml analysis
392
393 gSystem->ChangeDirectory(ocwd.Data());
34ce8d31 394}
37ff46b0 395
396//________________________________________________
ce8dba2f 397void GetAverageXsection(TTree * tree, Double_t & xs, Float_t & ntr)
37ff46b0 398{
399 // Read the PYTHIA statistics from the file pyxsec.root created by
400 // the function WriteXsection():
401 // integrated cross section (xsection) and
402 // the number of Pyevent() calls (ntrials)
403 // and calculate the weight per one event xsection/ntrials
404 // The spectrum calculated by a user should be
405 // multiplied by this weight, something like this:
406 // TH1F *userSpectrum ... // book and fill the spectrum
407 // userSpectrum->Scale(weight)
408 //
409 // Yuri Kharlov 19 June 2007
ce8dba2f 410 // Gustavo Conesa 15 April 2008
411 Double_t xsection = 0;
412 UInt_t ntrials = 0;
413 xs = 0;
414 ntr = 0;
415
416 Int_t nfiles = tree->GetEntries() ;
417 if (tree && nfiles > 0) {
418 tree->SetBranchAddress("xsection",&xsection);
419 tree->SetBranchAddress("ntrials",&ntrials);
420 for(Int_t i = 0; i < nfiles; i++){
421 tree->GetEntry(i);
422 xs += xsection ;
423 ntr += ntrials ;
424 cout << "xsection " <<xsection<<" ntrials "<<ntrials<<endl;
425 }
426
427 xs = xs / nfiles;
428 ntr = ntr / nfiles;
429 cout << "-----------------------------------------------------------------"<<endl;
430 cout << "Average of "<< nfiles<<" files: xsection " <<xs<<" ntrials "<<ntr<<endl;
431 cout << "-----------------------------------------------------------------"<<endl;
432 }
433 else cout << " >>>> Empty tree !!!! <<<<< "<<endl;
434
435}
37ff46b0 436
37ff46b0 437
37ff46b0 438