Improved version, kinetree included
[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
23 using namespace std; //required for resolving the 'cout' symbol
24
25 void testAnal(TString output = "flowEvtsAnal.root")
26 {
27  cout << " . Here the new flow analysis (2007) ... " << endl ;
28  cout << endl ;
29
30  bool kOne = kFALSE ;
31
32  TStopwatch timer;
33  timer.Start();
34
35  gSystem->Load("libPhysics.so");
36  //gSystem->Load("libPWG2.so");
37  gSystem->Load("libPWG2flow.so");
38
39  // output file //
40
41  TFile * fFlowfile = new TFile(output.Data(),"RECREATE") ;
42  //fFlowfile->cd() ; 
43
44  // selection object (cuts) //
45
46  AliFlowSelection* select = new AliFlowSelection() ;
47  
48  // Event Cuts
49  select->SetCentralityCut(-1) ;
50  select->SetRunIdCut(-1) ;
51
52  // R.P. calculation cuts
53  for(int j=0;j<AliFlowConstants::kHars;j++)
54  {
55   select->SetEtaCut(0., 1.1, j, 1) ;
56   select->SetPtCut(0.1, 10. , j, 1);
57  }
58  select->SetConstrainCut(kTRUE) ;
59  select->SetDcaGlobalCut(0.,0.1);
60  //select->SetNhitsCut(3., 1) ;
61  //select->SetPidCut("pi");
62
63  // Tracks Correlation analysis cuts
64  select->SetEtaPart(-1.1,1.1);
65  select->SetPtPart(0.1,10.);
66  select->SetConstrainablePart(kTRUE);
67  select->SetDcaGlobalPart(0.,0.1);
68  //select->SetPPart(0.1,10.);
69  //select->SetEtaAbsPart(0.5,1.1);
70  //select->SetPidPart("pi");
71  //select->SetYPart(..,..);
72  //select->SetPidProbPart(0.5,1.);
73  //select->SetFitPtsPart(3, 480);
74  //select->SetFitOverMaxPtsPart(0.1,3.);
75  //select->SetDedxPtsPart(1.,0.);
76  //select->SetChiSqPart(0.,100.);
77
78  // V0 Correlation analysis cuts
79  //select->SetV0DcaCross(0.,0.1) ;
80  select->SetV0Mass(0.4875,0.5078) ;     // Mk0 = 0.49765
81  select->SetV0SideBands(0.08) ;
82  //select->SetV0P(0.,10.) ;
83  select->SetV0Pt(0.1,10.) ;
84  select->SetV0Eta(-2.1,2.1) ;
85  //select->SetV0EtaAbs(0.,0.9) ;
86  //select->SetV0Pid("0") ;
87  //select->SetV0Y(0.,10.) ;
88  //select->SetV0Lenght(0.,100.) ;
89  //select->SetV0LenghtOverSigma(0.,5.) ;
90  //select->SetV0ChiSqPart(0.,100.) ;
91
92  // couts
93  cout << " . Selection for R.P. calculation: " << endl ;
94  select->PrintSelectionList() ;
95  cout << " . Selection for correlation analysis: " << endl ;
96  select->PrintList() ;
97  cout << " . Selection for V0 analysis: " << endl ;
98  select->PrintV0List() ;
99
100  // flow analyser //
101
102  AliFlowAnalyser* flow = new AliFlowAnalyser(select) ;
103  AliFlowEvent* flowEvt = new AliFlowEvent() ;
104
105  // output file
106  flow->SetHistFileName(output.Data()) ;
107  cout << " . Writing Histograms on  : " << flow->GetHistFileName()   << "  . " << endl ;
108
109  // Wgt file
110  TString wgtFileName = "flowPhiWgt.root" ;
111  TFile* wgtFile = new TFile(wgtFileName.Data(),"READ");
112  if(!wgtFile || wgtFile->IsZombie()) { cout << " . NO phi Weights . " << endl ;}
113  else { flow->FillWgtArrays(wgtFile) ; cout << " . Weights from  : " << flow->GetWgtFileName() << "  . " << endl ; }
114
115  // Analysis settings
116  flow->SetFlowForV0() ;         // default kTRUE.
117  flow->SetEtaSub() ;            // default kFALSE
118  //flow->SetV1Ep1Ep2() ;        // default kFALSE.
119  flow->SetShuffle(kFALSE) ;     // default kFALSE. shuffles track array
120  //flow->SetRedoWgt();          // default kFALSE. recalculates phiWgt (even if phiWgt file is already there)
121  flow->SetUsePhiWgt(kFALSE) ;   // default kTRUE if phiWgt file is there, kFALSE if not (& phiWgt file is created)
122  flow->SetUseOnePhiWgt() ; // or // flow->SetUseFirstLastPhiWgt() ; // uses 1 or 3 wgt histograms (default is 1)
123  flow->SetUseBayWgt(kFALSE) ;   // default kFALSE. uses bayesian weights in P.id.
124  //flow->SetUsePtWgt();         // default kFALSE. uses pT as a weight for RP determination
125  //flow->SetUseEtaWgt();        // default kFALSE. uses eta as a weight for RP determination
126  //flow->SetCustomRespFunc();   // default kFALSE. uses the combined response function from the ESD
127
128  // init (make histograms, start the analysis)
129  cout << endl ;
130  kOne = flow->Init() ;  if(kOne) { cout << "   ok! " << endl ;}
131  
132  // flowEvents chain (imput) //
133
134  TString fFlowFileName = "flowEvts.root" ;
135
136  // organizing event loop
137  TFile* flowEventsFile = new TFile(fFlowFileName.Data(), "READ") ; // flowEventsFile->ls() ;
138  Int_t nEvts = flowEventsFile->GetNkeys() ; 
139  cout << " . Found  " << nEvts << " AliFlowEvents in file " << fFlowFileName.Data() << endl ;
140  TList* flowEventsList = (TList*)flowEventsFile->GetListOfKeys() ; 
141
142  // event loop
143  cout << endl ;
144  for(Int_t ie=0;ie<nEvts;ie++)
145  {
146   TString evtName = flowEventsList->At(ie)->GetName() ;
147   flowEventsFile->GetObject(evtName.Data(),flowEvt) ;
148   // cout << "dumping event " << ie << " : " << endl ; flowEvt->Dump() ; cout << endl ; 
149   Bool_t succ = flow->Analyze(flowEvt) ;  
150   flow->PrintEventQuantities() ;
151   if(succ) { cout << ie << " done ... " << endl ; } 
152  }
153
154  // p.Id
155  cout << endl ;
156  cout << "Particles normalized abundance:" << endl;  // shows the bayesian vector
157  flow->PrintRunBayesian() ; cout << endl;
158  
159  // resolution
160  cout << endl ;
161  flow->Resolution() ;
162
163  // saves the wgt file
164  cout << endl ;
165  flow->Weightening() ;
166
167  // finish (save hists, close file)
168  cout << endl ;
169  kOne = flow->Finish() ; if(kOne) { cout << "   ok! " << endl ;}
170
171  timer.Stop();
172  cout << endl ;
173  cout << " . " ; timer.Print();
174  cout << " . here it was (analyser) ... " << endl ;  //juice!
175  cout << endl ;
176
177  // break ;
178  
179 }