Adding comments (Laurent)
[u/mrichter/AliRoot.git] / PWG2 / FLOW / AliFlowWeighter.cxx
1 //////////////////////////////////////////////////////////////////////
2 //
3 // $Id$
4 //
5 // Author: Emanuele Simili
6 //
7 //////////////////////////////////////////////////////////////////////
8 //_____________________________________________________________
9 //
10 // Description: 
11 //         the AliFlowWeighter class generates the phi-weights which 
12 // are later used to correct for azimuthal anisotropy in the reconstruction 
13 // efficiency of the ALICE TPC. It also fills an histogram of normalised 
14 // particle abundancies, which can be used as bayesian weights for particle Id. 
15 //
16 // - The method Init() generates the histograms and opens a new fPhiWgt file.
17 // - The method WeightEvent(AliFlowEvent*) fills phi and PId histograms. 
18 //   It must be inserted in a loop over the event sample.
19 // - The method Finish() calculates the weights, saves the histograms and 
20 //   closes the file. The AliFlowSelection object is saved as well.
21 //
22 //////////////////////////////////////////////////////////////////////
23
24 #ifndef ALIFLOWWEIGHTER_CXX
25 #define ALIFLOWWEIGHTER_CXX
26
27 // ROOT things
28 #include <TROOT.h>
29 #include <TMath.h>
30 #include <TFile.h>
31 #include <TTree.h>
32 #include <TString.h>
33 #include <TObject.h>
34 #include <TObjArray.h>
35 #include <TOrdCollection.h>
36 #include <TH1.h>
37 #include <TVector2.h>
38
39 // Flow things
40 #include "AliFlowEvent.h"
41 #include "AliFlowTrack.h"
42 #include "AliFlowV0.h"
43 #include "AliFlowConstants.h"
44 #include "AliFlowSelection.h"
45 #include "AliFlowWeighter.h"
46
47 // ANSI things
48 #include <math.h>
49 #include <stdio.h>
50 #include <stdlib.h>
51 #include <iostream>
52 using namespace std; //required for resolving the 'cout' symbol
53
54 ClassImp(AliFlowWeighter) 
55 //-----------------------------------------------------------------------
56 AliFlowWeighter::AliFlowWeighter(const AliFlowSelection* flowSelect):
57   fEventNumber(0), fTrackNumber(0), fNumberOfV0s(0), fNumberOfTracks(0), fPhiBins(0), fPhiMin(0.), fPhiMax(0.), 
58   fFlowEvent(0x0), fFlowTrack(0x0), fFlowSelect(0x0), fFlowTracks(0x0),
59   fWgtFile(0x0), fWgtFileName("flowPhiWgt.hist.root"), fPhiWgtHistList(0x0)
60 {
61  // default constructor (selection given or default selection)
62
63  if(flowSelect) { fFlowSelect = new AliFlowSelection(*flowSelect) ; }
64  else           { fFlowSelect = new AliFlowSelection() ; }
65  
66  // output file (histograms)
67  //fWgtFile     = 0 ;
68 }
69 //-----------------------------------------------------------------------
70 AliFlowWeighter::~AliFlowWeighter()
71 {
72  // default distructor (no actions) 
73 }
74 //-----------------------------------------------------------------------
75 Bool_t AliFlowWeighter::Init() 
76 {
77 // sets some defaults for the analysis
78  
79  cout << "* FlowWeighter *  -  Init()" << endl ; cout << endl ; 
80  
81  // Open output files (->plots)
82  fWgtFile = new TFile(fWgtFileName.Data(), "RECREATE");
83  fWgtFile->cd() ;
84
85  // counters and pointers to 0
86  fEventNumber = 0 ;
87  fTrackNumber = 0 ;     
88  fNumberOfTracks = 0 ;  
89  fNumberOfV0s    = 0 ;  
90  fFlowEvent  = 0 ;
91  fFlowTrack  = 0 ;
92  fFlowTracks = 0 ;
93  //for(Int_t ii=0;ii<3;ii++) { fVertex[ii] = 0 ; }
94
95  // Histogram settings
96  fPhiBins = AliFlowConstants::kPhiBins ;     
97  fPhiMin  = 0.;
98  fPhiMax  = 2*TMath::Pi() ; 
99
100  TString* histTitle ;
101  for(int k = 0; k < AliFlowConstants::kSels; k++)
102  {
103   histTitle = new TString("Flow_BayPidMult_Sel");
104   *histTitle += k+1;
105   fHistFull[k].fHistBayPidMult = new TH1F(histTitle->Data(), histTitle->Data(),AliFlowConstants::kPid,-0.5,((Float_t)AliFlowConstants::kPid-0.5));
106   fHistFull[k].fHistBayPidMult->Sumw2() ;
107   fHistFull[k].fHistBayPidMult->SetXTitle("e+/-  ,  mu+/-  ,  pi+/-  ,  K+/-  ,  p+/-  ,  d+/- ");
108   fHistFull[k].fHistBayPidMult->SetYTitle("Counts");
109   delete histTitle;
110
111   for(int j = 0; j < AliFlowConstants::kHars; j++) 
112   {
113    // Tpc - Phi lab
114    histTitle = new TString("Flow_Phi_TPC_Sel");
115    *histTitle += k+1;
116    histTitle->Append("_Har");
117    *histTitle += j+1;
118    fHistFull[k].fHistFullHar[j].fHistPhi = new TH1D(histTitle->Data(), histTitle->Data(), fPhiBins, fPhiMin, fPhiMax);
119    fHistFull[k].fHistFullHar[j].fHistPhi->SetXTitle("Azimuthal Angles (rad)");
120    fHistFull[k].fHistFullHar[j].fHistPhi->SetYTitle("Counts");
121    delete histTitle;
122    // Tpc (plus)
123    histTitle = new TString("Flow_Phi_TPCplus_Sel");
124    *histTitle += k+1;
125    histTitle->Append("_Har");
126    *histTitle += j+1;
127    fHistFull[k].fHistFullHar[j].fHistPhiPlus = new TH1D(histTitle->Data(), histTitle->Data(), fPhiBins, fPhiMin, fPhiMax);
128    fHistFull[k].fHistFullHar[j].fHistPhiPlus->SetXTitle("Azimuthal Angles (rad)");
129    fHistFull[k].fHistFullHar[j].fHistPhiPlus->SetYTitle("Counts");
130    delete histTitle;
131    // Tpc (minus)
132    histTitle = new TString("Flow_Phi_TPCminus_Sel");
133    *histTitle += k+1;
134    histTitle->Append("_Har");
135    *histTitle += j+1;
136    fHistFull[k].fHistFullHar[j].fHistPhiMinus = new TH1D(histTitle->Data(), histTitle->Data(), fPhiBins, fPhiMin, fPhiMax);
137    fHistFull[k].fHistFullHar[j].fHistPhiMinus->SetXTitle("Azimuthal Angles (rad)");
138    fHistFull[k].fHistFullHar[j].fHistPhiMinus->SetYTitle("Counts");
139    delete histTitle;
140    // Tpc (cross)
141    histTitle = new TString("Flow_Phi_TPCcross_Sel");
142    *histTitle += k+1;
143    histTitle->Append("_Har");
144    *histTitle += j+1;
145    fHistFull[k].fHistFullHar[j].fHistPhiAll = new TH1D(histTitle->Data(), histTitle->Data(), fPhiBins, fPhiMin, fPhiMax);
146    fHistFull[k].fHistFullHar[j].fHistPhiAll->SetXTitle("Azimuthal Angles (rad)");
147    fHistFull[k].fHistFullHar[j].fHistPhiAll->SetYTitle("Counts");
148    delete histTitle;
149
150    // Tpc - Phi lab flattened
151    histTitle = new TString("Flow_Phi_Flat_TPC_Sel");
152    *histTitle += k+1;
153    histTitle->Append("_Har");
154    *histTitle += j+1;
155    fHistFull[k].fHistFullHar[j].fHistPhiFlat = new TH1D(histTitle->Data(), histTitle->Data(), fPhiBins, fPhiMin, fPhiMax);
156    fHistFull[k].fHistFullHar[j].fHistPhiFlat->SetXTitle("Azimuthal Angles (rad)");
157    fHistFull[k].fHistFullHar[j].fHistPhiFlat->SetYTitle("Counts");
158    delete histTitle;
159    // Tpc Plus
160    histTitle = new TString("Flow_Phi_Flat_TPCplus_Sel");
161    *histTitle += k+1;
162    histTitle->Append("_Har");
163    *histTitle += j+1;
164    fHistFull[k].fHistFullHar[j].fHistPhiFlatPlus = new TH1D(histTitle->Data(), histTitle->Data(), fPhiBins, fPhiMin, fPhiMax);
165    fHistFull[k].fHistFullHar[j].fHistPhiFlatPlus->SetXTitle("Azimuthal Angles (rad)");
166    fHistFull[k].fHistFullHar[j].fHistPhiFlatPlus->SetYTitle("Counts");
167    delete histTitle;
168    // Tpc Minus
169    histTitle = new TString("Flow_Phi_Flat_TPCminus_Sel");
170    *histTitle += k+1;
171    histTitle->Append("_Har");
172    *histTitle += j+1;
173    fHistFull[k].fHistFullHar[j].fHistPhiFlatMinus = new TH1D(histTitle->Data(), histTitle->Data(), fPhiBins, fPhiMin, fPhiMax);
174    fHistFull[k].fHistFullHar[j].fHistPhiFlatMinus->SetXTitle("Azimuthal Angles (rad)");
175    fHistFull[k].fHistFullHar[j].fHistPhiFlatMinus->SetYTitle("Counts");
176    delete histTitle;
177    // Tpc cross
178    histTitle = new TString("Flow_Phi_Flat_TPCcross_Sel");
179    *histTitle += k+1;
180    histTitle->Append("_Har");
181    *histTitle += j+1;
182    fHistFull[k].fHistFullHar[j].fHistPhiFlatAll = new TH1D(histTitle->Data(), histTitle->Data(), fPhiBins, fPhiMin, fPhiMax);
183    fHistFull[k].fHistFullHar[j].fHistPhiFlatAll->SetXTitle("Azimuthal Angles (rad)");
184    fHistFull[k].fHistFullHar[j].fHistPhiFlatAll->SetYTitle("Counts");
185    delete histTitle;
186   }
187  }
188  
189  return kTRUE ;
190
191 //-----------------------------------------------------------------------
192
193 //-----------------------------------------------------------------------
194 Bool_t AliFlowWeighter::Finish() 
195 {
196  // Calls the method to fill wgt histograms, then saves them 
197  // on the fWgtFile and closes the file .
198  
199  cout << "* FlowWeighter *  -  Finish()" << endl ; cout << endl ;
200
201  Weightening() ;
202  fWgtFile->cd() ; 
203
204  // Write PhiWgt histograms
205  fPhiWgtHistList->Write();
206  delete fPhiWgtHistList ;
207  
208  // Write Bayesian Weights for P.Id.
209  for(int k=0;k<AliFlowConstants::kSels;k++) { fHistFull[k].fHistBayPidMult->Write() ; }
210
211  // Write the AliFlowSelection object
212  fFlowSelect->Write();
213  delete fFlowSelect ;
214
215  fWgtFile->Close();
216
217  cout << "    Finish()  -  Wgt file closed : " << fWgtFileName.Data() << endl ; cout << endl ; 
218
219  return kTRUE ;
220 }
221 //-----------------------------------------------------------------------
222 // ###
223 //-----------------------------------------------------------------------
224 Bool_t AliFlowWeighter::WeightEvent(AliFlowEvent* fFlowEvent)         
225 {
226  // Reads the AliFlowEvent (* fFlowEvent) and loops over the 
227  // AliFlowTraks to fill phi histograms .
228  
229  cout << " AliFlowWeighter::WeightEvent(" << fFlowEvent << " )   -   " << fEventNumber << endl ;
230  if(!fFlowEvent) { return kFALSE ; }
231
232  if(fFlowSelect->Select(fFlowEvent))     // event selected - here below the ANALYSIS FLAGS are setted -
233  {
234   fFlowTracks = fFlowEvent->TrackCollection() ; 
235   fNumberOfTracks = fFlowTracks->GetEntries() ;
236   cout << "       event ID = " << fFlowEvent->EventID() << " :  found " << fNumberOfTracks << " AliFlowTracks . " << endl ;  
237   fFlowEvent->SetNoWgt() ;
238   fFlowEvent->SetSelections(fFlowSelect) ; // does the selection of tracks for r.p. calculation (sets flags in AliFlowTrack)
239   
240   TracksLoop(fFlowTracks) ;  
241  }
242  else 
243  {
244   cout << "       * " << fEventNumber << " (event ID = " << fFlowEvent->EventID() << ") discarded . " << endl ; 
245   delete fFlowEvent  ; fFlowEvent = 0 ; 
246   return kFALSE ;
247  }
248  fEventNumber++ ;
249  
250  return kTRUE ;
251 }
252 //-----------------------------------------------------------------------
253 // ###
254 //-----------------------------------------------------------------------
255 void AliFlowWeighter::TracksLoop(TObjArray* fFlowTracks) 
256 {
257  // fills phi and PId histograms
258
259  cout << " Tracks Loop . " << endl ; 
260   
261  Float_t phi, eta, zFirstPoint ; // , zLastPoint ;
262  Int_t  fPidId ;
263  Char_t pid[10] = "0" ;
264  for(fTrackNumber=0;fTrackNumber<fNumberOfTracks;fTrackNumber++) 
265  {
266   fFlowTrack = (AliFlowTrack*)fFlowTracks->At(fTrackNumber) ;
267   // cout << "Track n. " << fTrackNumber << endl ; fFlowTrack->Dump() ; 
268  
269   phi = fFlowTrack->Phi() ;             
270   eta = fFlowTrack->Eta() ;             
271   zFirstPoint = fFlowTrack->ZFirstPoint() ; 
272   //zLastPoint = fFlowTrack->ZLastPoint() ;
273   strcpy(pid,fFlowTrack->Pid()) ; 
274
275   fPidId = -1 ;
276   if(strstr(pid,"e"))       { fPidId = 0 ; }
277   else if(strstr(pid,"mu")) { fPidId = 1 ; }
278   else if(strstr(pid,"pi")) { fPidId = 2 ; }
279   else if(strstr(pid,"k"))  { fPidId = 3 ; }
280   else if(strstr(pid,"pr")) { fPidId = 4 ; }
281   else if(strstr(pid,"d"))  { fPidId = 5 ; }
282
283   // Looping over Selections and Harmonics
284   for (int k = 0; k < AliFlowConstants::kSels; k++) 
285   {
286    fFlowSelect->SetSelection(k) ;
287    for (int j = 0; j < AliFlowConstants::kHars; j++) 
288    {
289     fFlowSelect->SetHarmonic(j);
290     if(fFlowSelect->Select(fFlowTrack))
291     {
292      Bool_t kTpcPlus  = kFALSE ;
293      Bool_t kTpcMinus = kFALSE ;
294      Bool_t kTpcAll   = kFALSE ;
295
296      // Set Tpc (+ and -)
297      if(fFlowTrack->FitPtsTPC())
298      {
299       if(zFirstPoint >= 0. && eta > 0.)      { kTpcPlus  = kTRUE ; } 
300       else if(zFirstPoint <= 0. && eta < 0.) { kTpcMinus = kTRUE ; }
301       else                                   { kTpcAll   = kTRUE ; }      
302      }
303     
304      // Phi distribution (particle for R.P.)
305      Float_t wt = 1. ; // TMath::Abs(fFlowEvent->Weight(k, j, fFlowTrack)) ; 
306      if(kTpcPlus)       { fHistFull[k].fHistFullHar[j].fHistPhiPlus->Fill(phi,wt) ;  } 
307      else if(kTpcMinus) { fHistFull[k].fHistFullHar[j].fHistPhiMinus->Fill(phi,wt) ; } 
308      else if(kTpcAll)   { fHistFull[k].fHistFullHar[j].fHistPhiAll->Fill(phi,wt) ;   } 
309      fHistFull[k].fHistFullHar[j].fHistPhi->Fill(phi,wt) ;
310
311      // PID Multiplicities (particle for R.P.) - just once for each selection
312      if(j==0)           { fHistFull[k].fHistBayPidMult->Fill(fPidId) ; }
313     }
314    }
315   }
316  }
317   
318  return ;
319 }
320 //-----------------------------------------------------------------------
321 Bool_t AliFlowWeighter::Weightening()
322 {
323  // Calculates weights, and fills PhiWgt histograms .
324  // This is called at the end of the event loop .
325  
326  cout << " AliFlowWeighter::Weightening() " << endl ; cout << endl ;
327  
328  // PhiWgt histogram collection
329  fPhiWgtHistList = new TOrdCollection(4*AliFlowConstants::kSels*AliFlowConstants::kHars) ;
330  
331  // Creates PhiWgt Histograms
332  TString* histTitle ;
333  for(Int_t k = 0; k < AliFlowConstants::kSels; k++)
334  {
335   for(Int_t j = 0; j < AliFlowConstants::kHars; j++) 
336   {
337    // Tpc plus
338    histTitle = new TString("Flow_Phi_Weight_TPCplus_Sel");
339    *histTitle += k+1;
340    histTitle->Append("_Har");
341    *histTitle += j+1;
342    fHistFull[k].fHistFullHar[j].fHistPhiWgtPlus = new TH1D(histTitle->Data(),histTitle->Data(), fPhiBins, fPhiMin, fPhiMax);
343    fHistFull[k].fHistFullHar[j].fHistPhiWgtPlus->Sumw2();
344    fHistFull[k].fHistFullHar[j].fHistPhiWgtPlus->SetXTitle("Azimuthal Angles (rad)");
345    fHistFull[k].fHistFullHar[j].fHistPhiWgtPlus->SetYTitle("PhiWgt");
346    delete histTitle;
347    // Tpc minus
348    histTitle = new TString("Flow_Phi_Weight_TPCminus_Sel");
349    *histTitle += k+1;
350    histTitle->Append("_Har");
351    *histTitle += j+1;
352    fHistFull[k].fHistFullHar[j].fHistPhiWgtMinus = new TH1D(histTitle->Data(),histTitle->Data(), fPhiBins, fPhiMin, fPhiMax);
353    fHistFull[k].fHistFullHar[j].fHistPhiWgtMinus->Sumw2();
354    fHistFull[k].fHistFullHar[j].fHistPhiWgtMinus->SetXTitle("Azimuthal Angles (rad)");
355    fHistFull[k].fHistFullHar[j].fHistPhiWgtMinus->SetYTitle("PhiWgt");
356    delete histTitle;
357    // Tpc cross
358    histTitle = new TString("Flow_Phi_Weight_TPCcross_Sel");
359    *histTitle += k+1;
360    histTitle->Append("_Har");
361    *histTitle += j+1;
362    fHistFull[k].fHistFullHar[j].fHistPhiWgtAll = new TH1D(histTitle->Data(),histTitle->Data(), fPhiBins, fPhiMin, fPhiMax);
363    fHistFull[k].fHistFullHar[j].fHistPhiWgtAll->Sumw2();
364    fHistFull[k].fHistFullHar[j].fHistPhiWgtAll->SetXTitle("Azimuthal Angles (rad)");
365    fHistFull[k].fHistFullHar[j].fHistPhiWgtAll->SetYTitle("PhiWgt");
366    delete histTitle;
367    // Tpc
368    histTitle = new TString("Flow_Phi_Weight_TPC_Sel");
369    *histTitle += k+1;
370    histTitle->Append("_Har");
371    *histTitle += j+1;
372    fHistFull[k].fHistFullHar[j].fHistPhiWgt = new TH1D(histTitle->Data(),histTitle->Data(), fPhiBins, fPhiMin, fPhiMax);
373    fHistFull[k].fHistFullHar[j].fHistPhiWgt->Sumw2();
374    fHistFull[k].fHistFullHar[j].fHistPhiWgt->SetXTitle("Azimuthal Angles (rad)");
375    fHistFull[k].fHistFullHar[j].fHistPhiWgt->SetYTitle("PhiWgt");
376    delete histTitle;
377
378    // Calculate PhiWgt
379    Double_t meanPlus  = fHistFull[k].fHistFullHar[j].fHistPhiPlus->Integral() / (Double_t)fPhiBins ;
380    Double_t meanMinus = fHistFull[k].fHistFullHar[j].fHistPhiMinus->Integral() / (Double_t)fPhiBins ;
381    Double_t meanCross = fHistFull[k].fHistFullHar[j].fHistPhiAll->Integral() / (Double_t)fPhiBins ;
382    Double_t meanTPC = fHistFull[k].fHistFullHar[j].fHistPhi->Integral() / (Double_t)fPhiBins ;
383
384    // Tpc
385    for(Int_t i=0;i<fPhiBins;i++) 
386    {
387     fHistFull[k].fHistFullHar[j].fHistPhiWgtPlus->SetBinContent(i+1,meanPlus);
388     fHistFull[k].fHistFullHar[j].fHistPhiWgtPlus->SetBinError(i+1, 0.);
389     fHistFull[k].fHistFullHar[j].fHistPhiWgtMinus->SetBinContent(i+1,meanMinus);
390     fHistFull[k].fHistFullHar[j].fHistPhiWgtMinus->SetBinError(i+1, 0.);
391     fHistFull[k].fHistFullHar[j].fHistPhiWgtAll->SetBinContent(i+1,meanCross);
392     fHistFull[k].fHistFullHar[j].fHistPhiWgtAll->SetBinError(i+1, 0.);
393     fHistFull[k].fHistFullHar[j].fHistPhiWgt->SetBinContent(i+1,meanTPC);
394     fHistFull[k].fHistFullHar[j].fHistPhiWgt->SetBinError(i+1, 0.);
395    }
396    
397    if(meanTPC==0) { cout << " Sel." << k << " , Har." << j << " :  empty histogram ! " << endl ; }
398    else 
399    {
400     fHistFull[k].fHistFullHar[j].fHistPhiWgtPlus->Divide(fHistFull[k].fHistFullHar[j].fHistPhiPlus);
401     fHistFull[k].fHistFullHar[j].fHistPhiWgtMinus->Divide(fHistFull[k].fHistFullHar[j].fHistPhiMinus);
402     fHistFull[k].fHistFullHar[j].fHistPhiWgtAll->Divide(fHistFull[k].fHistFullHar[j].fHistPhiAll);
403     fHistFull[k].fHistFullHar[j].fHistPhiWgt->Divide(fHistFull[k].fHistFullHar[j].fHistPhi);
404    }
405
406    fPhiWgtHistList->AddLast(fHistFull[k].fHistFullHar[j].fHistPhiWgtPlus);
407    fPhiWgtHistList->AddLast(fHistFull[k].fHistFullHar[j].fHistPhiWgtMinus);
408    fPhiWgtHistList->AddLast(fHistFull[k].fHistFullHar[j].fHistPhiWgtAll);
409    fPhiWgtHistList->AddLast(fHistFull[k].fHistFullHar[j].fHistPhiWgt);
410   }
411  }
412
413  return kTRUE ;
414 }
415 //-----------------------------------------------------------------------
416 void AliFlowWeighter::PrintBayesian(Int_t selN) const
417 {
418  // Prints the normalized particle abundance of all events (selection selN).
419  // Call this at the end of the loop, just before Finish() .
420
421  if(selN>AliFlowConstants::kSels) { selN = 0 ; } 
422  Char_t* names[AliFlowConstants::kPid] = {"e","mu","pi","k","p","d"} ;
423  Double_t bayes = 0. ;
424  Double_t totCount = (fHistFull[selN].fHistBayPidMult)->GetSumOfWeights() ;
425  if(totCount) 
426  { 
427   cout << " Sel." << selN << " particles normalized abundance (tot. " << totCount << " tracks) : " ;
428   for(Int_t ii=0;ii<AliFlowConstants::kPid;ii++)
429   {
430    bayes = (fHistFull[selN].fHistBayPidMult->GetBinContent(ii+1) / totCount) ; 
431    cout << bayes << "_" << names[ii] << " ; " ;
432   }
433  }
434  else {  cout << " Sel." << selN << " :  empty P.Id. histogram ! " << endl ; }
435  cout << endl ;
436  
437  return ;
438 }
439 //-----------------------------------------------------------------------
440
441
442 #endif