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