Fixed some coding violations
[u/mrichter/AliRoot.git] / PWG2 / FLOW / AliSelectorFoF.cxx
1 /* AliSelectorFoF.cxx, v1.1 01/02/2007 esimili Exp */
2 /* derived from AliSelector.cxx,v 1.17 2006/08/31 jgrosseo Exp */
3
4 // The class definition in esdV0.h has been generated automatically
5 // by the ROOT utility TTree::MakeSelector(). This class is derived
6 // from the ROOT class TSelector. For more information on the TSelector
7 // framework see $ROOTSYS/README/README.SELECTOR or the ROOT User Manual.
8
9 // The following methods are defined in this file:
10 //    Begin():        called everytime a loop on the tree starts,
11 //                    a convenient place to create your histograms.
12 //    SlaveBegin():   called after Begin(), when on PROOF called only on the
13 //                    slave servers.
14 //    Process():      called for each event, in this function you decide what
15 //                    to read and fill your histograms.
16 //    SlaveTerminate: called at the end of the loop on the tree, when on PROOF
17 //                    called only on the slave servers.
18 //    Terminate():    called at the end of the loop on the tree,
19 //                    a convenient place to draw/fit your histograms.
20 //
21 // To use this file, try the following session on your Tree T:
22 //
23 // Root > T->Process("AliSelector.C")
24 // Root > T->Process("AliSelector.C","some options")
25 // Root > T->Process("AliSelector.C+")
26 //
27
28 #include "AliSelectorFoF.h"
29
30 #include <TStyle.h>
31 #include <TSystem.h>
32 #include <TCanvas.h>
33 #include <TRegexp.h>
34 #include <TTime.h>
35 #include <TFriendElement.h>
36 #include <TTree.h>
37 #include <TChain.h>
38 #include <TFile.h>
39 #include <TTimeStamp.h>
40 #include <TMath.h>
41 #include <TObject.h>
42 #include <TObjArray.h>
43 #include <TSelector.h>
44
45 #include "AliLog.h"      
46 #include "AliESD.h"      
47 #include "AliESDtrack.h"  
48 #include "AliESDv0.h"     
49
50 #include <stdio.h>
51 #include <stdlib.h>
52 #include <iostream>
53 using namespace std; //required for resolving the 'cout' symbol
54
55 #include "AliFlowEvent.h"
56 #include "AliFlowTrack.h"
57 #include "AliFlowV0.h"
58 #include "AliFlowConstants.h"
59 #include "AliFlowSelection.h"
60 #include "AliFlowMaker.h"
61 #include "AliFlowAnalyser.h"
62 #include "AliFlowWeighter.h"
63
64 ClassImp(AliSelectorFoF)
65
66 //-----------------------------------------------------------------------
67
68 AliSelectorFoF::AliSelectorFoF() :
69   TSelector(),
70   fTree(0),
71   fESD(0),
72   fCountFiles(0),
73   fKineFile(0)
74 {
75   //
76   // Constructor. Initialization of pointers
77   //
78   
79   fFlowEventFileName    = "fof_flowEvts.root" ;    
80   fFlowWgtFileName      = "fof_flowWgts.root" ;    
81   fFlowAnalysisFileName = "fof_flowAnal.root" ;    
82 }
83
84 //-----------------------------------------------------------------------
85
86 AliSelectorFoF::~AliSelectorFoF()
87 {
88   //
89   // Destructor
90   //
91
92  if (fTree) { fTree->ResetBranchAddresses() ; }
93
94  if (fESD)
95  {
96    delete fESD;
97    fESD = 0;
98  }
99 }
100
101 //-----------------------------------------------------------------------
102
103 void AliSelectorFoF::CheckOptions()
104 {
105   // checks the option string for the debug flag
106
107   AliLog::SetClassDebugLevel(ClassName(), AliLog::kInfo);
108
109   TString option = GetOption();
110
111   if (option.Contains("moredebug"))
112   {
113     printf("Enabling verbose debug mode for %s\n", ClassName());
114     AliLog::SetClassDebugLevel(ClassName(), AliLog::kDebug+1);
115     AliInfo(Form("Called with option %s.", option.Data()));
116   }
117   else if (option.Contains("debug"))
118   {
119     printf("Enabling debug mode for %s\n", ClassName());
120     AliLog::SetClassDebugLevel(ClassName(), AliLog::kDebug);
121     AliInfo(Form("Called with option %s.", option.Data()));
122   }
123 }
124
125 //-----------------------------------------------------------------------
126
127 void AliSelectorFoF::Begin(TTree*)
128 {
129   // The Begin() function is called at the start of the query.
130   // When running with PROOF Begin() is only called on the client.
131   // The tree argument is deprecated (on PROOF 0 is passed).
132   
133   cout << " HERE I begin !!! " << endl ; cout << endl ;
134
135   CheckOptions();
136
137   AliDebug(AliLog::kDebug, "============BEGIN===========");
138
139  // flags
140  fDoNothing     = kFALSE ;  // kFALSE ;
141  fOnFlyAnalysis = kTRUE ;  // kTRUE ;    
142  fOnFlyWeight   = kFALSE ;  // kTRUE ;
143  fSaveFlowEvents= kTRUE ;   // kTRUE ;
144
145 // Maker part :
146  fFlowMaker = new AliFlowMaker() ;
147  // ESD Cuts
148  fFlowMaker->SetNHitsCut(1) ;
149  fFlowMaker->SetECut(0.01,100.) ; 
150  fFlowMaker->PrintCutList() ;
151
152 // Opens flow event (output) file
153  if(fSaveFlowEvents)  
154  { 
155  // flow events file (output)
156   fFlowfile = new TFile(fFlowEventFileName.Data(),"RECREATE") ;
157   fFlowfile->cd() ; 
158   cout << " . Writing AliFlowEvents on  : " << fFlowEventFileName.Data()   << "  . " << endl ;
159  }
160  
161 // Analysis part :
162  if(fOnFlyWeight || fOnFlyAnalysis)  
163  { 
164   cout << " . Here the flow selection ... " << endl ;
165   cout << endl ;
166
167  // AliFlowSelection...
168   fFlowSelect = new AliFlowSelection() ;
169   // Event Cuts
170   fFlowSelect->SetCentralityCut(-1) ;
171   fFlowSelect->SetRunIdCut(-1) ;
172   // R.P. calculation cuts
173   for(int j=0;j<AliFlowConstants::kHars;j++)
174   {
175    fFlowSelect->SetEtaCut(0., 1.1, j, 1) ;
176    fFlowSelect->SetPtCut(0.1, 10. , j, 1);
177   }
178   fFlowSelect->SetConstrainCut(kTRUE) ;
179   fFlowSelect->SetDcaGlobalCut(0.,0.1);
180   // Correlation analysis cuts (not all of them)
181   fFlowSelect->SetEtaPart(-1.1,1.1);
182   fFlowSelect->SetPtPart(0.1,10.);
183   fFlowSelect->SetConstrainablePart(kTRUE);
184   fFlowSelect->SetDcaGlobalPart(0.,0.1);
185   // V0 analysis cuts (not all of them ... they are useless anyway)
186   fFlowSelect->SetV0Mass(0.4875,0.5078) ;        // Mk0 = 0.49765
187   fFlowSelect->SetV0SideBands(0.08) ;
188   fFlowSelect->SetV0Pt(0.1,10.) ;
189   fFlowSelect->SetV0Eta(-2.1,2.1) ;
190   // print list :
191   cout << " . Selection for R.P. calculation: " << endl ;
192   fFlowSelect->PrintSelectionList() ;
193   cout << " . Selection for correlation analysis: " << endl ;
194   fFlowSelect->PrintList() ;
195   cout << " . Selection for V0 analysis: " << endl ;
196   fFlowSelect->PrintV0List() ;
197
198 // Weight part :
199   if(fOnFlyWeight)  
200   { 
201    cout << " . Here the flow weighter ... " << endl ;
202    cout << endl ;
203
204   // AliFlowWeighter...
205    fFlowWeighter = new AliFlowWeighter(fFlowSelect) ;
206
207   // Wgt file (create)
208    fFlowWgtFile = new TFile(fFlowWgtFileName.Data(),"RECREATE") ;
209    fFlowWeighter->SetWgtFile(fFlowWgtFile) ;
210    cout << " . Writing Wgt Histograms on  : " << fFlowWeighter->GetWgtFileName()  << "  . " << endl ;
211    fFlowWeighter->Init() ;
212   }
213
214   if(fOnFlyAnalysis)  
215   { 
216    cout << " . Here the flow analysis ... " << endl ;
217    cout << endl ;
218
219   // AliFlowAnalyser...
220    fFlowAnal = new AliFlowAnalyser(fFlowSelect) ;
221
222   // Wgt file (read)
223    //TString wgtFileName = "fFlowPhiWgt.root" ;
224    //TFile* wgtFile = new TFile(wgtFileName.Data(),"READ");
225    //fFlowAnal->FillWgtArrays(wgtFile) ; // fix this !!!
226    //cout << " . Weights from  : " << fFlowAnal->GetWgtFileName() << "  . " << endl ;
227
228   // analysis file (output)
229    fFlowAnal->SetHistFileName(fFlowAnalysisFileName.Data()) ;
230    cout << " . Writing Analysis Histograms on  : " << fFlowAnal->GetHistFileName()   << "  . " << endl ;
231
232   // Analysis settings
233    fFlowAnal->SetFlowForV0(kFALSE) ;         // default kTRUE.
234    fFlowAnal->SetEtaSub() ;            // default kFALSE
235    //fFlowAnal->SetV1Ep1Ep2() ;        // default kFALSE.
236    fFlowAnal->SetShuffle() ;           // default kFALSE. shuffles track array
237    //fFlowAnal->SetRedoWgt();             // default kFALSE. recalculates phiWgt (even if phiWgt file is already there)
238    fFlowAnal->SetUsePhiWgt(kFALSE) ;   // default kTRUE if phiWgt file is there, kFALSE if not (& phiWgt file is created)
239    fFlowAnal->SetUseOnePhiWgt() ; // or // fFlowAnal->SetUseFirstLastPhiWgt() ; // uses 1 or 3 wgt histograms (default is 1)
240    fFlowAnal->SetUseBayWgt(kFALSE) ;   // default kFALSE. uses bayesian weights in P.id.
241    //fFlowAnal->SetUsePtWgt();  // default kFALSE. uses pT as a weight for RP determination
242    //fFlowAnal->SetUseEtaWgt(); // default kFALSE. uses eta as a weight for RP determination
243    fFlowAnal->Init() ;
244   }
245  }
246
247  cout << "                ... init is done " << endl ; 
248 }
249
250 //-----------------------------------------------------------------------
251
252 void AliSelectorFoF::SlaveBegin(TTree* tree)
253 {
254   // The SlaveBegin() function is called after the Begin() function.
255   // When running with PROOF SlaveBegin() is called on each slave server.
256   // The tree argument is deprecated (on PROOF 0 is passed).
257
258   CheckOptions();
259
260   AliDebug(AliLog::kDebug, "=======SLAVEBEGIN========");
261   AliDebug(AliLog::kDebug, Form("Hostname: %s", gSystem->HostName()));
262   AliDebug(AliLog::kDebug, Form("Time: %s", gSystem->Now().AsString()));
263
264   if (tree != 0) { Init(tree) ; }
265 }
266
267 //-----------------------------------------------------------------------
268
269 void AliSelectorFoF::Init(TTree *tree)
270 {
271   // The Init() function is called when the selector needs to initialize
272   // a new tree or chain. Typically here the branch addresses of the tree
273   // will be set. It is normaly not necessary to make changes to the
274   // generated code, but the routine can be extended by the user if needed.
275   // Init() will be called many times when running with PROOF.
276
277   AliDebug(AliLog::kDebug, "=========Init==========");
278
279   fTree = tree;
280
281   if (fTree == 0)
282   {
283    AliDebug(AliLog::kError, "ERROR: tree argument is 0.");
284    return;
285   }
286
287   // Set branch address
288   fTree->SetBranchAddress("ESD", &fESD);
289   if (fESD != 0) { AliDebug(AliLog::kInfo, "INFO: Found ESD branch in chain.") ; }
290 }
291
292 //-----------------------------------------------------------------------
293
294 Bool_t AliSelectorFoF::Notify()
295 {
296   // The Notify() function is called when a new file is opened. This
297   // can be either for a new TTree in a TChain or when when a new TTree
298   // is started when using PROOF. Typically here the branch pointers
299   // will be retrieved. It is normaly not necessary to make changes
300   // to the generated code, but the routine can be extended by the
301   // user if needed.
302
303   AliDebug(AliLog::kDebug, "=========NOTIFY==========");
304   AliDebug(AliLog::kDebug, Form("Hostname: %s", gSystem->HostName()));
305   AliDebug(AliLog::kDebug, Form("Time: %s", TTimeStamp(time(0)).AsString()));
306
307   ++fCountFiles;
308   if (fTree)
309   {
310     TFile *f = fTree->GetCurrentFile();
311     AliDebug(AliLog::kInfo, Form("Processing %d. file %s", fCountFiles, f->GetName()));
312   }
313   else
314   {
315     AliDebug(AliLog::kError, "fTree not available");
316   }
317
318   DeleteKinematicsFile();
319
320   return kTRUE;
321 }
322
323 //-----------------------------------------------------------------------
324
325 Bool_t AliSelectorFoF::Process(Long64_t entry)
326 {
327   // The Process() function is called for each entry in the tree (or possibly
328   // keyed object in the case of PROOF) to be processed. The entry argument
329   // specifies which entry in the currently loaded tree is to be processed.
330   // It can be passed to either TTree::GetEntry() or TBranch::GetEntry()
331   // to read either all or the required parts of the data. When processing
332   // keyed objects with PROOF, the object is already loaded and is available
333   // via the fObject pointer.
334   //
335   // This function should contain the "body" of the analysis. It can contain
336   // simple or elaborate selection criteria, run algorithms on the data
337   // of the event and typically fill histograms.
338
339   // WARNING when a selector is used with a TChain, you must use
340   //  the pointer to the current TTree to call GetEntry(entry).
341   //  The entry is always the local entry number in the current tree.
342   //  Assuming that fTree is the pointer to the TChain being processed,
343   //  use fTree->GetTree()->GetEntry(entry).
344
345   AliDebug(AliLog::kDebug, Form("=========PROCESS========== Entry %lld", entry));
346
347   if(!fTree)
348   {
349    AliDebug(AliLog::kError, "ERROR: fTree is 0.") ;
350    return kFALSE ;
351   }
352
353   fEventNumber = entry ;
354   fTree->GetTree()->GetEntry(fEventNumber) ;
355
356   if(fESD) { AliDebug(AliLog::kDebug, Form("ESD: We have %d tracks.", fESD->GetNumberOfTracks())); }
357   cout << " event !!! " << entry << endl ;
358
359   fRunID = fESD->GetRunNumber() ;
360   fEventNumber = fESD->GetEventNumber() ;
361   fNumberOfTracks = fESD->GetNumberOfTracks() ;
362   fNumberOfV0s = fESD->GetNumberOfV0s() ;
363
364   cout << " *evt n. " << fEventNumber << " (run " << fRunID << ") " << endl ;
365   cout << "  tracks: " << fNumberOfTracks << " ,   v0s " << fNumberOfV0s << endl ;
366
367  // Dummy Loop
368   if(fDoNothing)  { cout << " ... doing nothing ... " << endl ; return kTRUE; } 
369
370  // Instantiate a new AliFlowEvent
371   cout << " filling the flow event :| " << endl ;
372   fFlowEvent = fFlowMaker->FillFlowEvent(fESD) ;
373   if(!fFlowEvent) { cout << "!something bad occurred" << endl ; return kFALSE ; }
374   else            { cout << " event filled :) " << entry << "                # t.G. ! #       "<< endl ; }
375
376  // Saves the AliFlowEvent
377   if(fSaveFlowEvents)
378   { 
379    cout << " saving flow event :| " << endl ;
380    TString strID = "" ; strID += entry ; 
381    fFlowfile->cd() ; fFlowEvent->Write(strID.Data()) ;
382    cout << " saved :) " << strID.Data() << "                # t.g. ! #       " << endl ; cout << endl ;  
383   }
384
385  // Weights
386   if(fOnFlyWeight)
387   { 
388    cout << " weightening :| " << endl ;
389    bool done = fFlowWeighter->WeightEvent(fFlowEvent) ;  
390    if(done) { cout << " weights done :) " << entry << "                # t.h.v. ! #       " << endl ; } 
391   }
392   
393  // On-Fly Analysis
394   if(fOnFlyAnalysis)  
395   {
396    cout << " doing analysis :| " << endl ;
397    bool done = fFlowAnal->Analyse(fFlowEvent) ;
398    if(done) { cout << " analysis done :) " << entry << "                # t.h.v. ! #       " << endl ; } 
399   }
400    
401   return kTRUE;
402 }
403
404 //-----------------------------------------------------------------------
405
406 void AliSelectorFoF::SlaveTerminate()
407 {
408   // The SlaveTerminate() function is called after all entries or objects
409   // have been processed. When running with PROOF SlaveTerminate() is called
410   // on each slave server.
411
412   AliDebug(AliLog::kDebug, "=======SLAVETERMINATE=======");
413
414   DeleteKinematicsFile();
415 }
416
417 //-----------------------------------------------------------------------
418
419 void AliSelectorFoF::Terminate()
420 {
421   // The Terminate() function is the last function to be called during
422   // a query. It always runs on the client, it can be used to present
423   // the results graphically or save the results to file.
424
425   AliDebug(AliLog::kDebug, "=========TERMINATE==========");
426
427   cout << " Finished ... " << endl ;
428
429   if(fDoNothing) { cout << "             ... & nothing was done ! " << endl ; cout << endl ; return ; }
430
431   cout << "  nTracks:  " << fFlowMaker->GetNgoodTracks() << endl ;   
432   cout << "  nV0s:  " << fFlowMaker->GetNgoodV0s()  << endl ;        
433   cout << "  nTracks (|eta|<0.5):  " << fFlowMaker->GetNgoodTracksEta() << endl ; 
434   cout << "  nTracks+:  " << fFlowMaker->GetNposiTracks() << endl ;          
435   cout << "  nTracks-:  " << fFlowMaker->GetNnegaTracks() << endl ;          
436   cout << "  nTracks unconstrained:  " << fFlowMaker->GetNunconstrained() << endl ;      
437   cout << "  Bayesian :  " ; 
438   for(int ii=0;ii<5;ii++) { cout << fFlowMaker->GetBayesianNorm(ii) << "   " ; } 
439   cout << " . " << endl ; 
440
441   if(fOnFlyWeight) { fFlowWeighter->Finish() ; delete fFlowWeighter ; }
442   
443   if(fOnFlyAnalysis) 
444   { 
445    cout << endl ; cout << "Particles normalized abundance:" << endl;  // shows the bayesian vector
446    fFlowAnal->Resolution() ; cout << endl;
447    fFlowAnal->PrintRunBayesian() ; cout << endl;
448    fFlowAnal->Finish() ;  cout << endl;
449    delete fFlowAnal ;
450   }
451      
452   if(fSaveFlowEvents) { fFlowfile->Close() ; cout << " file closed . " << endl ; }
453   delete fFlowMaker ;
454
455   delete fFlowSelect ;
456   
457   cout << endl ; 
458   return ;
459
460
461 //-----------------------------------------------------------------------
462
463 TTree* AliSelectorFoF::GetKinematics()
464 {
465   // Returns kinematics tree corresponding to current ESD active in fTree
466   // Loads the kinematics from the kinematics file, the file is identified by replacing "AliESDs" to
467   // "Kinematics" in the file path of the ESD file. This is a hack, to be changed!
468
469   if (!fKineFile)
470   {
471     if(!fTree->GetCurrentFile()) { return 0 ; }
472
473     TString fileName(fTree->GetCurrentFile()->GetName());
474     fileName.ReplaceAll("AliESDs", "Kinematics");
475
476     // temporary workaround for PROOF bug #18505
477     fileName.ReplaceAll("#Kinematics.root#Kinematics.root", "#Kinematics.root");
478
479     AliDebug(AliLog::kInfo, Form("Opening %s", fileName.Data()));
480
481     fKineFile = TFile::Open(fileName);
482     if(!fKineFile) { return 0 ; }
483   }
484
485   return dynamic_cast<TTree*> (fKineFile->Get(Form("Event%d/TreeK", fTree->GetTree()->GetReadEntry())));
486 }
487
488 //-----------------------------------------------------------------------
489
490 void AliSelectorFoF::DeleteKinematicsFile()
491 {
492   //
493   // Closes the kinematics file and deletes the pointer.
494   //
495
496   if (fKineFile)
497   {
498     fKineFile->Close();
499     delete fKineFile;
500     fKineFile = 0;
501   }
502 }
503
504 //-----------------------------------------------------------------------
505
506