adding the option for hybrid tracks global+global constrained as well as global+TPC...
[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=0, mGRID=3};
14 //mLocal    = 0: Analyze locally files in your computer
15 //mGRID     = 3: Analyze files on GRID
16
17 //---------------------------------------------------------------------------
18 //Settings to read locally several files, only for "mLocal" mode
19 //The different values are default, they can be set with environmental 
20 //variables: INDIR, PATTERN, NFILES, respectivelly
21
22 char * kInDir = "/user/data/files/"; 
23 char * kPattern = ""; // Data are in files kInDir/kPattern+i 
24 Int_t kFile = 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
36 //---------------------------------------------------------------------------
37
38 const Bool_t kMC = kFALSE; //With real data kMC = kFALSE
39 const TString kInputData = "ESD"; //ESD, AOD, MC, deltaAOD
40 const Bool_t  outAOD = kFALSE; //Some tasks doesnt need it.
41 TString kTreeName;
42 const Bool_t kUsePAR = kFALSE; //set to kFALSE for libraries
43 const Int_t kFilter = kFALSE; //Use ESD filter
44
45 void ana(Int_t mode=mLocal)
46 {
47   // Main
48   
49   //--------------------------------------------------------------------
50   // Load analysis libraries
51   // Look at the method below, 
52   // change whatever you need for your analysis case
53   // ------------------------------------------------------------------
54   LoadLibraries() ;
55   // TGeoManager::Import("geometry.root") ; //need file "geometry.root" in local dir!!!!
56   
57   //-------------------------------------------------------------------------------------------------
58   //Create chain from ESD and from cross sections files, look below for options.
59   //-------------------------------------------------------------------------------------------------
60   if(kInputData == "ESD") kTreeName = "esdTree" ;
61   else if(kInputData.Contains("AOD")) kTreeName = "aodTree" ;
62   else if (kInputData == "MC") kTreeName = "TE" ;
63   else {
64     cout<<"Wrong  data type "<<kInputData<<endl;
65     break;
66   }
67   
68   if(kFilter) outAOD = kTRUE;
69   
70   TChain *chain       = new TChain(kTreeName) ;
71   TChain * chainxs = new TChain("Xsection") ;
72   CreateChain(mode, chain, chainxs);  
73   
74   if(chain){
75     AliLog::SetGlobalLogLevel(AliLog::kError);//Minimum prints on screen
76     
77     //--------------------------------------
78     // Make the analysis manager
79     //-------------------------------------
80     AliAnalysisManager *mgr  = new AliAnalysisManager("Manager", "Manager");
81     //AliAnalysisManager::SetUseProgressBar(kTRUE);
82     //mgr->SetSkipTerminate(kTRUE);
83     //mgr->SetNSysInfo(1);
84     
85     // MC handler
86     if((kMC || kInputData == "MC") && !kInputData.Contains("AOD")){
87       AliMCEventHandler* mcHandler = new AliMCEventHandler();
88       mcHandler->SetReadTR(kFALSE);//Do not search TrackRef file
89       mgr->SetMCtruthEventHandler(mcHandler);
90       if( kInputData == "MC") {
91         cout<<"MC INPUT EVENT HANDLER"<<endl;
92         mgr->SetInputEventHandler(NULL);
93       }
94     }
95     
96     // AOD output handler
97     if(kInputData!="deltaAOD" && outAOD){
98       cout<<"Init output handler"<<endl;
99       AliAODHandler* aodoutHandler   = new AliAODHandler();
100       aodoutHandler->SetOutputFileName("aod.root");
101       ////aodoutHandler->SetCreateNonStandardAOD();
102       mgr->SetOutputEventHandler(aodoutHandler);
103     }
104     
105     //input
106     
107     if(kInputData == "ESD"){
108       // ESD handler
109       AliESDInputHandler *esdHandler = new AliESDInputHandler();
110       esdHandler->SetReadFriends(kFALSE);
111       mgr->SetInputEventHandler(esdHandler);
112       cout<<"ESD handler "<<mgr->GetInputEventHandler()<<endl;
113     }
114     else if(kInputData.Contains("AOD")){
115       // AOD handler
116       AliAODInputHandler *aodHandler = new AliAODInputHandler();
117       mgr->SetInputEventHandler(aodHandler);
118       if(kInputData == "deltaAOD") aodHandler->AddFriend("deltaAODPartCorr.root");
119       cout<<"AOD handler "<<mgr->GetInputEventHandler()<<endl;
120     }
121     //mgr->RegisterExternalFile("deltaAODPartCorr.root");
122     //mgr->SetDebugLevel(-1); // For debugging, do not uncomment if you want no messages.
123     
124     //-------------------------------------------------------------------------
125     //Define task, put here any other task that you want to use.
126     //-------------------------------------------------------------------------
127     AliAnalysisDataContainer *cinput1 = mgr->GetCommonInputContainer();
128     AliAnalysisDataContainer *coutput1;
129     if(outAOD){
130       
131       coutput1 = mgr->GetCommonOutputContainer();
132       
133       if(kInputData=="ESD"){
134         
135         gROOT->LoadMacro("$ALICE_ROOT/ANALYSIS/macros/AddTaskPhysicsSelection.C");
136         AliPhysicsSelectionTask* physSelTask = AddTaskPhysicsSelection();
137         
138         if(kFilter){
139
140           // Set of cuts
141           //for standard global track cuts
142           AliESDtrackCuts* esdTrackCutsGlobal =  AliESDtrackCuts::GetStandardITSTPCTrackCuts2009(kTRUE);
143           esdTrackCutsGlobal->SetName("StandardFromAliESDTrackCuts");
144           
145           //for TPC tracks only
146           //    AliESDtrackCuts* esdTrackCutsTPC =  AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();
147           //      esdTrackCutsTPC->SetRequireTPCRefit(kTRUE); 
148           //      esdTrackCutsTPC->SetMinNClustersTPC(70);
149           
150           AliAnalysisFilter* trackFilter = new AliAnalysisFilter("trackFilter");
151           trackFilter->AddCuts(esdTrackCutsGlobal); 
152           //trackFilter->AddCuts(esdTrackCutsTPC); 
153           
154           gROOT->LoadMacro("$ALICE_ROOT/ANALYSIS/macros/AddTaskESDFilter.C");
155           AliAnalysisTaskESDfilter *taskesdfilter = AddTaskESDFilter(kFALSE);
156           esdfilter->SetTrackFilter(trackFilter);
157           kInputData = "AOD" ;
158           kTreeName = "aodTree" ;
159         }
160       }
161     }
162     
163     gROOT->LoadMacro("$ALICE_ROOT/PWG4/macros/AddTaskPartCorr.C");
164    
165     AliAnalysisTaskParticleCorrelation *taskEMCAL = AddTaskPartCorr(kInputData,"EMCAL", kFALSE,kFALSE, kFALSE); 
166     //mgr->ProfileTask("PartCorrEMCAL");
167     
168     //AliAnalysisTaskParticleCorrelation *taskPHOS  = AddTaskPartCorr(kInputData,"PHOS", kFALSE,kFALSE,kFALSE);
169     //mgr->ProfileTask("PartCorrPHOS");
170     
171     //gROOT->LoadMacro("$ALICE_ROOT/PWG4/macros/QA/AddTaskCalorimeterQA.C");
172     //AliAnalysisTaskParticleCorrelation *taskQA = AddTaskCalorimeterQA(kInputData,kFALSE,kFALSE);
173     //mgr->ProfileTask("CalorimeterPerformance");      
174     
175     //-----------------------
176     // Run the analysis
177     //-----------------------    
178     mgr->InitAnalysis();
179     mgr->PrintStatus();
180     mgr->StartAnalysis("local",chain);
181
182     cout <<" Analysis ended sucessfully "<< endl ;
183     
184   }
185   else cout << "Chain was not produced ! "<<endl;
186   //printf("*** DELETE MANAGER ***\n");
187   // delete mgr;
188   
189 }
190
191 void  LoadLibraries() {
192   
193   //--------------------------------------
194   // Load the needed libraries most of them already loaded by aliroot
195   //--------------------------------------
196   gSystem->Load("libTree.so");
197   gSystem->Load("libGeom.so");
198   gSystem->Load("libVMC.so");
199   gSystem->Load("libXMLIO.so");
200   gSystem->Load("libMatrix.so");
201   gSystem->Load("libPhysics.so");
202
203   if(kUsePAR){
204     //--------------------------------------------------------
205     //If you want to use root and par files from aliroot
206     //--------------------------------------------------------  
207     SetupPar("STEERBase");
208     SetupPar("ESD");
209     SetupPar("AOD");
210     SetupPar("ANALYSIS");
211     SetupPar("ANALYSISalice");
212     SetupPar("PHOSUtils");
213     SetupPar("EMCALUtils");
214     SetupPar("PWG4PartCorrBase");
215     SetupPar("PWG4PartCorrDep");
216   }
217   else{
218     //--------------------------------------------------------
219     // If you want to use already compiled libraries 
220     // in the aliroot distribution
221     //--------------------------------------------------------
222     gSystem->Load("libSTEERBase.so");
223     gSystem->Load("libESD.so");
224     gSystem->Load("libAOD.so");
225     gSystem->Load("libANALYSIS.so");
226     gSystem->Load("libANALYSISalice.so");
227     gSystem->Load("libPHOSUtils");
228     gSystem->Load("libEMCALUtils");
229     gSystem->Load("libPWG4PartCorrBase.so");
230     gSystem->Load("libPWG4PartCorrDep.so");
231     if(kFilter){
232       gSystem->Load("libCORRFW.so");
233       gSystem->Load("libPWG3base.so");
234       gSystem->Load("libPWG3muon.so");
235     }
236   }
237 }
238
239 void SetupPar(char* pararchivename)
240 {
241   //Load par files, create analysis libraries
242   //For testing, if par file already decompressed and modified
243   //classes then do not decompress.
244  
245   TString cdir(Form("%s", gSystem->WorkingDirectory() )) ; 
246   TString parpar(Form("%s.par", pararchivename)) ; 
247 //  if ( gSystem->AccessPathName(parpar.Data()) ) {
248 //    gSystem->ChangeDirectory(gSystem->Getenv("ALICE_ROOT")) ;
249 //     TString processline(Form(".! make %s", parpar.Data())) ; 
250 //     gROOT->ProcessLine(processline.Data()) ;
251 //     gSystem->ChangeDirectory(cdir) ; 
252 //     processline = Form(".! mv $ALICE_ROOT/%s .", parpar.Data()) ;
253 //     gROOT->ProcessLine(processline.Data()) ;
254 //   } 
255   if ( gSystem->AccessPathName(pararchivename) ) {  
256     TString processline = Form(".! tar xvzf %s",parpar.Data()) ;
257     gROOT->ProcessLine(processline.Data());
258   }
259   
260   TString ocwd = gSystem->WorkingDirectory();
261   gSystem->ChangeDirectory(pararchivename);
262   
263   // check for BUILD.sh and execute
264   if (!gSystem->AccessPathName("PROOF-INF/BUILD.sh")) {
265     printf("*******************************\n");
266     printf("*** Building PAR archive    ***\n");
267     cout<<pararchivename<<endl;
268     printf("*******************************\n");
269     
270     if (gSystem->Exec("PROOF-INF/BUILD.sh")) {
271       Error("runProcess","Cannot Build the PAR Archive! - Abort!");
272       return -1;
273     }
274   }
275   // check for SETUP.C and execute
276   if (!gSystem->AccessPathName("PROOF-INF/SETUP.C")) {
277     printf("*******************************\n");
278     printf("*** Setup PAR archive       ***\n");
279     cout<<pararchivename<<endl;
280     printf("*******************************\n");
281     gROOT->Macro("PROOF-INF/SETUP.C");
282   }
283   
284   gSystem->ChangeDirectory(ocwd.Data());
285   printf("Current dir: %s\n", ocwd.Data());
286 }
287
288
289
290 void CreateChain(const anaModes mode, TChain * chain, TChain * chainxs){
291   //Fills chain with data
292   TString ocwd = gSystem->WorkingDirectory();
293   
294   //---------------------------------------
295   //Local files analysis
296   //---------------------------------------
297   if(mode == mLocal){    
298     //If you want to add several ESD files sitting in a common directory INDIR
299     //Specify as environmental variables the directory (INDIR), the number of files 
300     //to analyze (NFILES) and the pattern name of the directories with files (PATTERN)
301
302     if(gSystem->Getenv("INDIR"))  
303       kInDir = gSystem->Getenv("INDIR") ; 
304     else cout<<"INDIR not set, use default: "<<kInDir<<endl;    
305     
306     if(gSystem->Getenv("PATTERN"))   
307       kPattern = gSystem->Getenv("PATTERN") ; 
308     else  cout<<"PATTERN not set, use default: "<<kPattern<<endl;
309     
310     if(gSystem->Getenv("NFILES"))
311       kFile = atoi(gSystem->Getenv("NFILES")) ;
312     else cout<<"NFILES not set, use default: "<<kFile<<endl;
313     
314     //Check if env variables are set and are correct
315     if ( kInDir  && kFile) {
316       printf("Get %d files from directory %s\n",kFile,kInDir);
317       if ( ! gSystem->cd(kInDir) ) {//check if ESDs directory exist
318         printf("%s does not exist\n", kInDir) ;
319         return ;
320       }
321
322       //if(gSystem->Getenv("XSFILE"))  
323       //kXSFileName = gSystem->Getenv("XSFILE") ; 
324       //else cout<<" XS file name not set, use default: "<<kXSFileName<<endl;   
325       char * kGener = gSystem->Getenv("GENER");
326       if(kGener) {
327         cout<<"GENER "<<kGener<<endl;
328         if(!strcmp(kGener,"PYTHIA")) kXSFileName = "pyxsec.root";
329         else if(!strcmp(kGener,"HERWIG")) kXSFileName = "hexsec.root";
330         else cout<<" UNKNOWN GENER, use default: "<<kXSFileName<<endl;
331       }
332       else cout<<" GENER not set, use default xs file name: "<<kXSFileName<<endl;
333
334       cout<<"INDIR   : "<<kInDir<<endl;
335       cout<<"NFILES  : "<<kFile<<endl;
336       cout<<"PATTERN : " <<kPattern<<endl;
337       cout<<"XSFILE  : "<<kXSFileName<<endl;
338
339       TString datafile="";
340       if(kInputData == "ESD") datafile = "AliESDs.root" ;
341       else if(kInputData.Contains("AOD")) datafile = "AliAOD.root" ;
342       else if(kInputData == "MC")  datafile = "galice.root" ;
343       
344       //Loop on ESD files, add them to chain
345       Int_t event =0;
346       Int_t skipped=0 ; 
347       char file[120] ;
348       char filexs[120] ;
349       
350       for (event = 0 ; event < kFile ; event++) {
351         sprintf(file, "%s/%s%d/%s", kInDir,kPattern,event,datafile.Data()) ; 
352         sprintf(filexs, "%s/%s%d/%s", kInDir,kPattern,event,kXSFileName) ; 
353         TFile * fESD = 0 ; 
354         //Check if file exists and add it, if not skip it
355         if ( fESD = TFile::Open(file)) {
356           if ( fESD->Get(kTreeName) ) { 
357             printf("++++ Adding %s\n", file) ;
358             chain->AddFile(file);
359             chainxs->Add(filexs) ; 
360           }
361         }
362         else { 
363           printf("---- Skipping %s\n", file) ;
364           skipped++ ;
365         }
366       }
367       printf("number of entries # %lld, skipped %d\n", chain->GetEntries(), skipped*100) ;      
368     }
369     else {
370       TString input = "AliESDs.root" ;
371       cout<<">>>>>> No list added, take a single file <<<<<<<<< "<<input<<endl;
372       chain->AddFile(input);
373     }
374     
375   }// local files analysis
376   
377   //------------------------------
378   //GRID xml files
379   //-----------------------------
380   else if(mode == mGRID){
381     //Get colection file. It is specified by the environmental
382     //variable XML
383
384     if(gSystem->Getenv("XML") )
385       kXML = gSystem->Getenv("XML");
386     else
387       sprintf(kXML, "collection.xml") ; 
388     
389     if (!TFile::Open(kXML)) {
390       printf("No collection file with name -- %s -- was found\n",kXML);
391       return ;
392     }
393     else cout<<"XML file "<<kXML<<endl;
394
395     //Load necessary libraries and connect to the GRID
396     gSystem->Load("libNetx.so") ; 
397     gSystem->Load("libRAliEn.so"); 
398     TGrid::Connect("alien://") ;
399
400     //Feed Grid with collection file
401     //TGridCollection * collection =  (TGridCollection*)gROOT->ProcessLine(Form("TAlienCollection::Open(\"%s\", 0)", kXML));
402     TGridCollection * collection = (TGridCollection*) TAlienCollection::Open(kXML);
403     if (! collection) {
404       AliError(Form("%s not found", kXML)) ; 
405       return kFALSE ; 
406     }
407     TGridResult* result = collection->GetGridResult("",0 ,0);
408    
409     // Makes the ESD chain 
410     printf("*** Getting the Chain       ***\n");
411     for (Int_t index = 0; index < result->GetEntries(); index++) {
412       TString alienURL = result->GetKey(index, "turl") ; 
413       cout << "================== " << alienURL << endl ; 
414       chain->Add(alienURL) ; 
415       alienURL.ReplaceAll("AliESDs.root",kXSFileName);
416       chainxs->Add(alienURL) ; 
417     }
418   }// xml analysis
419   
420   gSystem->ChangeDirectory(ocwd.Data());
421 }
422
423 //________________________________________________
424 void GetAverageXsection(TTree * tree, Double_t & xs, Float_t & ntr)
425 {
426   // Read the PYTHIA statistics from the file pyxsec.root created by
427   // the function WriteXsection():
428   // integrated cross section (xsection) and
429   // the  number of Pyevent() calls (ntrials)
430   // and calculate the weight per one event xsection/ntrials
431   // The spectrum calculated by a user should be
432   // multiplied by this weight, something like this:
433   // TH1F *userSpectrum ... // book and fill the spectrum
434   // userSpectrum->Scale(weight)
435   //
436   // Yuri Kharlov 19 June 2007
437   // Gustavo Conesa 15 April 2008
438   Double_t xsection = 0;
439   UInt_t    ntrials = 0;
440   xs = 0;
441   ntr = 0;
442   
443   Int_t nfiles =  tree->GetEntries()  ;
444   if (tree && nfiles > 0) {
445     tree->SetBranchAddress("xsection",&xsection);
446     tree->SetBranchAddress("ntrials",&ntrials);
447     for(Int_t i = 0; i < nfiles; i++){
448       tree->GetEntry(i);
449       xs += xsection ;
450       ntr += ntrials ;
451       cout << "xsection " <<xsection<<" ntrials "<<ntrials<<endl; 
452     }
453     
454     xs =   xs /  nfiles;
455     ntr =  ntr / nfiles;
456     cout << "-----------------------------------------------------------------"<<endl;
457     cout << "Average of "<< nfiles<<" files: xsection " <<xs<<" ntrials "<<ntr<<endl; 
458     cout << "-----------------------------------------------------------------"<<endl;
459   } 
460   else cout << " >>>> Empty tree !!!! <<<<< "<<endl;
461   
462 }
463
464
465