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