namespace added for constants
[u/mrichter/AliRoot.git] / PWG2 / FLOW / macros / testAnal.C
1 /////////////////////////////////////////////////////////////
2 //
3 // $Id$
4 //
5 // Author: Emanuele Simili
6 //
7 /////////////////////////////////////////////////////////////
8 //
9 // Description: ROOT macro to perform the Flow analysis on AliFlowEvents (new way)
10 //
11 /////////////////////////////////////////////////////////////
12
13 #include <vector>
14 #include <iostream>
15 #include <fstream>
16 #include "TVector.h"
17 #include "TVector2.h"
18 #include "TVector3.h"
19 #include "TMath.h"
20 #include "TObjArray"
21 #include "TStopwatch.h"
22 #include "TKey.h"
23 #include "TList.h"
24
25 using namespace std; //required for resolving the 'cout' symbol
26
27 int testAnal(int cen = -1)
28 {
29  cout << " . Here the new flow analysis (2007 rs) ... " << endl ;
30  cout << endl ;
31
32  bool kOne = kFALSE ;
33  int limit = -1 ;
34
35  if(limit>0)
36  {
37   cout << " . limited to " << limit << "events . " << endl ;
38   cout << endl ;
39  }
40  
41  // histogram file (output) //
42
43  TString output = "flowEvtsAnal" ; 
44  if(cen >= 0) { output += cen ; }
45  output += ".root" ;
46  
47  // flowEvents chain (imput) //
48
49  TString fFlowFileName = "flowEvts" ;
50  if(cen >= 0) { fFlowFileName += cen ; }
51  fFlowFileName += ".root" ;
52
53  // Wgt file (optional) //
54
55  TString wgtFileName = "flowPhiWgt" ;
56  if(cen >= 0) { wgtFileName += cen ; }
57  wgtFileName += ".root" ;
58
59  // start, load libs //
60
61  TStopwatch timer;
62  timer.Start();
63
64  //gSystem->Load("libPhysics.so");
65  gSystem->Load("libPWG2flow.so");
66
67  // open output file //
68
69  TFile * fFlowfile = new TFile(output.Data(),"RECREATE") ;
70  //fFlowfile->cd() ; 
71
72  // selection object (cuts) // 
73
74  AliFlowSelection* select = new AliFlowSelection() ;
75  
76  // Event Cuts
77  select->SetCentralityCut(-1) ;  // cen
78  select->SetRunIdCut(-1) ;
79
80  // R.P. calculation cuts
81  for(int j=0;j<AliFlowConstants::kHars;j++)
82  {
83   select->SetEtaCut(0., 0.9, j, 1) ;
84   select->SetPtCut(0.1, 10. , j, 1);
85  }
86  select->SetConstrainCut(kTRUE) ;
87  select->SetDcaGlobalCut(0.,0.1);
88  //select->SetNhitsCut(0, 1) ;
89  //select->SetPidCut("pi");
90
91  // Tracks Correlation analysis cuts
92  select->SetEtaPart(-0.9,0.9);
93  select->SetPtPart(0.1,10.);
94  select->SetConstrainablePart(kTRUE);
95  select->SetDcaGlobalPart(0.,0.1);
96  //select->SetPPart(0.1,10.);
97  //select->SetEtaAbsPart(0.5,1.1);
98  //select->SetPidPart("pi");
99  //select->SetYPart(..,..);
100  //select->SetPidProbPart(0.5,1.);
101  //select->SetFitPtsPart(3, 480);
102  //select->SetFitOverMaxPtsPart(0.1,3.);
103  //select->SetDedxPtsPart(1.,0.);
104  //select->SetChiSqPart(0.,100.);
105
106  // V0 Correlation analysis cuts
107  //select->SetV0DcaCross(0.,0.1) ;
108  select->SetV0Mass(0.48,0.52) ;     // Mk0 = 0.49765
109  select->SetV0SideBands(0.10) ;
110  //select->SetV0P(0.,10.) ;
111  select->SetV0Pt(0.1,10.) ;
112  select->SetV0Eta(-2.1,2.1) ;
113  //select->SetV0EtaAbs(0.,0.9) ;
114  //select->SetV0Pid("0") ;
115  //select->SetV0Y(0.,10.) ;
116  //select->SetV0Lenght(0.,100.) ;
117  //select->SetV0LenghtOverSigma(0.,5.) ;
118  //select->SetV0ChiSqPart(0.,100.) ;
119
120  // couts
121  cout << " . Selection for R.P. calculation: " << endl ;
122  select->PrintSelectionList() ;
123  cout << " . Selection for correlation analysis: " << endl ;
124  select->PrintList() ;
125  cout << " . Selection for V0 analysis: " << endl ;
126  select->PrintV0List() ;
127
128  // flow analyser //
129
130  AliFlowAnalyser* flow = new AliFlowAnalyser(select) ;
131
132  // output file
133  flow->SetHistFileName(output.Data()) ;
134  cout << " . Writing Histograms on  : " << flow->GetHistFileName()   << "  . " << endl ;
135
136  // Wgt s
137  TFile* wgtFile = new TFile(wgtFileName.Data(),"READ");
138  if(!wgtFile || wgtFile->IsZombie()) { cout << " . NO phi Weights . " << endl ;}
139  else { flow->FillWgtArrays(wgtFile) ; cout << " . Weights from  : " << flow->GetWgtFileName() << "  . " << endl ; }
140
141  // Analysis settings
142  flow->SetFlowForV0(kFALSE) ;   // default kTRUE.
143  flow->SetSub(1) ; //eta // flow->SetSub(0) ; //rnd // flow->SetSub(-1) ; //charged
144  //flow->SetV1Ep1Ep2() ;        // default kFALSE.
145  //flow->SetRedoWgt();          // default kFALSE. recalculates phiWgt (even if phiWgt file is already there)
146  flow->SetUsePhiWgt(kFALSE) ;   // default kTRUE if phiWgt file is there, kFALSE if not (& phiWgt file is created)
147  flow->SetUseOnePhiWgt() ; // or // flow->SetUseFirstLastPhiWgt() ; // uses 1 or 3 wgt histograms (default is 1)
148  flow->SetUseBayWgt(kFALSE) ;   // default kFALSE. uses bayesian weights in P.id.
149  //flow->SetUsePtWgt();         // default kFALSE. uses pT as a weight for RP determination
150  //flow->SetUseEtaWgt();        // default kFALSE. uses eta as a weight for RP determination
151  //flow->SetCustomRespFunc();   // default kFALSE. uses the combined response function from the ESD
152
153  // init (make histograms, start the analysis)
154  cout << endl ;
155  kOne = flow->Init() ;  if(kOne) { cout << "   ok! " << endl ;}
156
157  // organizing event loop
158  TFile* flowEventsFile = new TFile(fFlowFileName.Data(), "READ") ; // flowEventsFile->ls() ;
159  Int_t nEvts = flowEventsFile->GetNkeys() ; 
160  cout << " . Found  " << nEvts << " AliFlowEvents in file " << fFlowFileName.Data() << endl ;
161  TList* flowEventsList = (TList*)flowEventsFile->GetListOfKeys() ; 
162  AliFlowEvent* flowEvt = new AliFlowEvent() ;
163  TIter next(flowEventsList); 
164  TKey* key ;
165
166  cout << endl ;
167
168  // Loop over the events
169  Int_t count = 0 ;
170  while( key=(TKey *)next() ) 
171  {
172   // TString evtName = flowEventsList->At(ie)->GetName() ;
173   // flowEventsFile->GetObject(evtName.Data(),flowEvt) ;
174   // cout << "dumping event " << ie << " : " << endl ; flowEvt->Dump() ; cout << endl ; 
175
176   flowEvt = (AliFlowEvent *)key->ReadObj();
177   if(!flowEvt) break;
178
179   // Process event .......
180   Bool_t succ = flow->Analyze(flowEvt) ;  
181   flow->PrintEventQuantities() ;
182   if(succ) { cout << count << " done ... " << endl ; }  
183   
184   delete flowEvt ;
185
186   if(count == limit) { break ; }  
187   count++ ;
188  }
189
190  // p.Id
191  cout << endl ;
192  cout << "Particles normalized abundance:" << endl;  // shows the bayesian vector
193  flow->PrintRunBayesian() ; cout << endl;
194  
195  // resolution
196  cout << endl ;
197  flow->Resolution() ;
198
199  // saves the wgt file
200  cout << endl ;
201  flow->Weightening() ;
202
203  // finish (save hists, close file)
204  cout << endl ;
205  kOne = flow->Finish() ; if(kOne) { cout << "   ok! " << endl ;}
206
207  timer.Stop();
208  cout << endl ;
209  cout << " . " ; timer.Print();
210  cout << " . here it was (analyser rs) ... " << endl ;  //juice!
211  cout << endl ;
212  
213  // cout << endl ; cout << " Memory Check (from Paul)" << endl ; 
214  // gObjectTable->Print();
215  // cout << endl ; cout << endl ;
216
217  return cen ;
218 }