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