namespace added for constants
[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., 0.9, 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(-0.9,0.9);
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.1) ;
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 = "fof_flowWgts.root" ;
224    TFile* wgtFile = new TFile(wgtFileName.Data(),"READ");
225    if(!wgtFile || wgtFile->IsZombie()) { cout << " . NO phi Weights . " << endl ;}
226    else { fFlowAnal->FillWgtArrays(wgtFile) ; 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() ;         // default kTRUE.
234    fFlowAnal->SetSub(1) ; //eta  // ->SetChrSub() ; //charge  // ->SetRndSub() ; //random (default)
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->SetCustomRespFunc();    // default kFALSE. uses the combined response function from the ESD
244    fFlowAnal->Init() ;
245   }
246  }
247
248  cout << "                ... init is done " << endl ; 
249 }
250
251 //-----------------------------------------------------------------------
252
253 void AliSelectorFoF::SlaveBegin(TTree* tree)
254 {
255   // The SlaveBegin() function is called after the Begin() function.
256   // When running with PROOF SlaveBegin() is called on each slave server.
257   // The tree argument is deprecated (on PROOF 0 is passed).
258
259   CheckOptions();
260
261   AliDebug(AliLog::kDebug, "=======SLAVEBEGIN========");
262   AliDebug(AliLog::kDebug, Form("Hostname: %s", gSystem->HostName()));
263   AliDebug(AliLog::kDebug, Form("Time: %s", gSystem->Now().AsString()));
264
265   if (tree != 0) { Init(tree) ; }
266 }
267
268 //-----------------------------------------------------------------------
269
270 void AliSelectorFoF::Init(TTree *tree)
271 {
272   // The Init() function is called when the selector needs to initialize
273   // a new tree or chain. Typically here the branch addresses of the tree
274   // will be set. It is normaly not necessary to make changes to the
275   // generated code, but the routine can be extended by the user if needed.
276   // Init() will be called many times when running with PROOF.
277
278   AliDebug(AliLog::kDebug, "=========Init==========");
279
280   fTree = tree;
281
282   if (fTree == 0)
283   {
284    AliDebug(AliLog::kError, "ERROR: tree argument is 0.");
285    return;
286   }
287
288   // Set branch address
289   fTree->SetBranchAddress("ESD", &fESD);
290   if (fESD != 0) { AliDebug(AliLog::kInfo, "INFO: Found ESD branch in chain.") ; }
291 }
292
293 //-----------------------------------------------------------------------
294
295 Bool_t AliSelectorFoF::Notify()
296 {
297   // The Notify() function is called when a new file is opened. This
298   // can be either for a new TTree in a TChain or when when a new TTree
299   // is started when using PROOF. Typically here the branch pointers
300   // will be retrieved. It is normaly not necessary to make changes
301   // to the generated code, but the routine can be extended by the
302   // user if needed.
303
304   AliDebug(AliLog::kDebug, "=========NOTIFY==========");
305   AliDebug(AliLog::kDebug, Form("Hostname: %s", gSystem->HostName()));
306   AliDebug(AliLog::kDebug, Form("Time: %s", TTimeStamp(time(0)).AsString()));
307
308   ++fCountFiles;
309   if (fTree)
310   {
311     TFile *f = fTree->GetCurrentFile();
312     AliDebug(AliLog::kInfo, Form("Processing %d. file %s", fCountFiles, f->GetName()));
313   }
314   else
315   {
316     AliDebug(AliLog::kError, "fTree not available");
317   }
318
319   DeleteKinematicsFile();
320
321   return kTRUE;
322 }
323
324 //-----------------------------------------------------------------------
325
326 Bool_t AliSelectorFoF::Process(Long64_t entry)
327 {
328   // The Process() function is called for each entry in the tree (or possibly
329   // keyed object in the case of PROOF) to be processed. The entry argument
330   // specifies which entry in the currently loaded tree is to be processed.
331   // It can be passed to either TTree::GetEntry() or TBranch::GetEntry()
332   // to read either all or the required parts of the data. When processing
333   // keyed objects with PROOF, the object is already loaded and is available
334   // via the fObject pointer.
335   //
336   // This function should contain the "body" of the analysis. It can contain
337   // simple or elaborate selection criteria, run algorithms on the data
338   // of the event and typically fill histograms.
339
340   // WARNING when a selector is used with a TChain, you must use
341   //  the pointer to the current TTree to call GetEntry(entry).
342   //  The entry is always the local entry number in the current tree.
343   //  Assuming that fTree is the pointer to the TChain being processed,
344   //  use fTree->GetTree()->GetEntry(entry).
345
346   AliDebug(AliLog::kDebug, Form("=========PROCESS========== Entry %lld", entry));
347
348   if(!fTree)
349   {
350    AliDebug(AliLog::kError, "ERROR: fTree is 0.") ;
351    return kFALSE ;
352   }
353
354   fEventNumber = entry ;
355   fTree->GetTree()->GetEntry(fEventNumber) ;
356
357   if(fESD) { AliDebug(AliLog::kDebug, Form("ESD: We have %d tracks.", fESD->GetNumberOfTracks())); }
358   cout << " event !!! " << entry << endl ;
359
360   fRunID = fESD->GetRunNumber() ;
361   //  fEventNumber = fESD->GetEventNumber() ;
362   fEventNumber = -1 ;
363   fNumberOfTracks = fESD->GetNumberOfTracks() ;
364   fNumberOfV0s = fESD->GetNumberOfV0s() ;
365
366   cout << " *evt n. " << fEventNumber << " (run " << fRunID << ") " << endl ;
367   cout << "  tracks: " << fNumberOfTracks << " ,   v0s " << fNumberOfV0s << endl ;
368
369  // Dummy Loop
370   if(fDoNothing)  { cout << " ... doing nothing ... " << endl ; return kTRUE; } 
371
372  // Instantiate a new AliFlowEvent
373   cout << " filling the flow event :| " << endl ;
374   fFlowEvent = fFlowMaker->FillFlowEvent(fESD) ;
375   if(!fFlowEvent) { cout << "! something bad occurred !" << endl ; return kFALSE ; }
376   else            { cout << "# event done :) " << entry << "                # ok ! #       " << endl ; }
377
378  // Saves the AliFlowEvent
379   if(fSaveFlowEvents)
380   { 
381    cout << " saving flow event :| " << endl ;
382    TString strID = "" ; strID += entry ; 
383    fFlowfile->cd() ; fFlowEvent->Write(strID.Data()) ;
384    cout << "# event saved :) " << strID.Data() << "                # ok ! #       " << endl ; cout << endl ;  
385   }
386
387  // Weights
388   if(fOnFlyWeight)
389   { 
390    cout << " weightening :| " << endl ;
391    bool done = fFlowWeighter->WeightEvent(fFlowEvent) ;  
392    if(done) { cout << "# weights done :) " << entry << "                # ok ! #       " << endl ; } 
393   }
394   
395  // On-Fly Analysis
396   if(fOnFlyAnalysis)  
397   {
398    cout << " doing analysis :| " << endl ;
399    bool done = fFlowAnal->Analyze(fFlowEvent) ;
400    fFlowAnal->PrintEventQuantities() ;
401    if(done) { cout << "# analysis done :) " << entry << "                # ok ! #       " << endl ; } 
402   }
403    
404   return kTRUE;
405 }
406
407 //-----------------------------------------------------------------------
408
409 void AliSelectorFoF::SlaveTerminate()
410 {
411   // The SlaveTerminate() function is called after all entries or objects
412   // have been processed. When running with PROOF SlaveTerminate() is called
413   // on each slave server.
414
415   AliDebug(AliLog::kDebug, "=======SLAVETERMINATE=======");
416
417   DeleteKinematicsFile();
418 }
419
420 //-----------------------------------------------------------------------
421
422 void AliSelectorFoF::Terminate()
423 {
424   // The Terminate() function is the last function to be called during
425   // a query. It always runs on the client, it can be used to present
426   // the results graphically or save the results to file.
427
428   AliDebug(AliLog::kDebug, "=========TERMINATE==========");
429
430   cout << " Finished ... " << endl ;
431
432   if(fDoNothing) { cout << "             ... & nothing was done ! " << endl ; cout << endl ; return ; }
433
434   cout << endl ;
435   cout << "  nTracks:  " << fFlowMaker->GetNgoodTracks() << endl ;   
436   cout << "  nV0s:  " << fFlowMaker->GetNgoodV0s()  << endl ;        
437   cout << "  nTracks (|eta|<0.5):  " << fFlowMaker->GetNgoodTracksEta() << endl ; 
438   cout << "  nTracks+:  " << fFlowMaker->GetNposiTracks() << endl ;          
439   cout << "  nTracks-:  " << fFlowMaker->GetNnegaTracks() << endl ;          
440   cout << "  nTracks unconstrained:  " << fFlowMaker->GetNunconstrained() << endl ;      
441   cout << "  Bayesian (e,mu,pi,k,p) :  " ; 
442   for(int ii=0;ii<5;ii++) { cout << fFlowMaker->GetBayesianNorm(ii) << "   " ; } 
443   cout << " . " << endl ; cout << endl ;
444
445   if(fOnFlyWeight) { fFlowWeighter->Finish() ; delete fFlowWeighter ; }
446   
447   if(fOnFlyAnalysis) 
448   { 
449    fFlowAnal->Resolution() ; cout << endl;
450    fFlowAnal->Finish() ;  cout << endl;
451    //cout << endl ; cout << "Particles normalized abundance:" << endl;  // shows the bayesian vector
452    //fFlowAnal->PrintRunBayesian() ; cout << endl;
453    delete fFlowAnal ;
454   }
455      
456   if(fSaveFlowEvents) { fFlowfile->Close() ; cout << " file closed . " << endl ; }
457   delete fFlowMaker ;
458
459   delete fFlowSelect ;
460   
461   cout << endl ; 
462   return ;
463
464
465 //-----------------------------------------------------------------------
466
467 TTree* AliSelectorFoF::GetKinematics()
468 {
469   // Returns kinematics tree corresponding to current ESD active in fTree
470   // Loads the kinematics from the kinematics file, the file is identified by replacing "AliESDs" to
471   // "Kinematics" in the file path of the ESD file. This is a hack, to be changed!
472
473   if (!fKineFile)
474   {
475     if(!fTree->GetCurrentFile()) { return 0 ; }
476
477     TString fileName(fTree->GetCurrentFile()->GetName());
478     fileName.ReplaceAll("AliESDs", "Kinematics");
479
480     // temporary workaround for PROOF bug #18505
481     fileName.ReplaceAll("#Kinematics.root#Kinematics.root", "#Kinematics.root");
482
483     AliDebug(AliLog::kInfo, Form("Opening %s", fileName.Data()));
484
485     fKineFile = TFile::Open(fileName);
486     if(!fKineFile) { return 0 ; }
487   }
488
489   return dynamic_cast<TTree*> (fKineFile->Get(Form("Event%d/TreeK", fTree->GetTree()->GetReadEntry())));
490 }
491
492 //-----------------------------------------------------------------------
493
494 void AliSelectorFoF::DeleteKinematicsFile()
495 {
496   //
497   // Closes the kinematics file and deletes the pointer.
498   //
499
500   if (fKineFile)
501   {
502     fKineFile->Close();
503     delete fKineFile;
504     fKineFile = 0;
505   }
506 }
507
508 //-----------------------------------------------------------------------
509
510