Adding comments (Laurent)
[u/mrichter/AliRoot.git] / PWG2 / FLOW / AliFlowWeighter.cxx
CommitLineData
327288af 1//////////////////////////////////////////////////////////////////////
9777bfcb 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
327288af 24#ifndef ALIFLOWWEIGHTER_CXX
25#define ALIFLOWWEIGHTER_CXX
9777bfcb 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
51a9b7d8 48#include <math.h>
9777bfcb 49#include <stdio.h>
50#include <stdlib.h>
51#include <iostream>
52using namespace std; //required for resolving the 'cout' symbol
53
da4d7f17 54ClassImp(AliFlowWeighter)
9777bfcb 55//-----------------------------------------------------------------------
4e566f2f 56AliFlowWeighter::AliFlowWeighter(const AliFlowSelection* flowSelect):
c1a88675 57 fEventNumber(0), fTrackNumber(0), fNumberOfV0s(0), fNumberOfTracks(0), fPhiBins(0), fPhiMin(0.), fPhiMax(0.),
4e566f2f 58 fFlowEvent(0x0), fFlowTrack(0x0), fFlowSelect(0x0), fFlowTracks(0x0),
c1a88675 59 fWgtFile(0x0), fWgtFileName("flowPhiWgt.hist.root"), fPhiWgtHistList(0x0)
9777bfcb 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)
c1a88675 67 //fWgtFile = 0 ;
9777bfcb 68}
69//-----------------------------------------------------------------------
70AliFlowWeighter::~AliFlowWeighter()
71{
72 // default distructor (no actions)
73}
74//-----------------------------------------------------------------------
75Bool_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
327288af 96 fPhiBins = AliFlowConstants::kPhiBins ;
9777bfcb 97 fPhiMin = 0.;
98 fPhiMax = 2*TMath::Pi() ;
99
100 TString* histTitle ;
327288af 101 for(int k = 0; k < AliFlowConstants::kSels; k++)
9777bfcb 102 {
103 histTitle = new TString("Flow_BayPidMult_Sel");
104 *histTitle += k+1;
327288af 105 fHistFull[k].fHistBayPidMult = new TH1F(histTitle->Data(), histTitle->Data(),AliFlowConstants::kPid,-0.5,((Float_t)AliFlowConstants::kPid-0.5));
9777bfcb 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
327288af 111 for(int j = 0; j < AliFlowConstants::kHars; j++)
9777bfcb 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//-----------------------------------------------------------------------
194Bool_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.
327288af 209 for(int k=0;k<AliFlowConstants::kSels;k++) { fHistFull[k].fHistBayPidMult->Write() ; }
9777bfcb 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//-----------------------------------------------------------------------
224Bool_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//-----------------------------------------------------------------------
255void 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
327288af 284 for (int k = 0; k < AliFlowConstants::kSels; k++)
9777bfcb 285 {
286 fFlowSelect->SetSelection(k) ;
327288af 287 for (int j = 0; j < AliFlowConstants::kHars; j++)
9777bfcb 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//-----------------------------------------------------------------------
321Bool_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
327288af 329 fPhiWgtHistList = new TOrdCollection(4*AliFlowConstants::kSels*AliFlowConstants::kHars) ;
9777bfcb 330
331 // Creates PhiWgt Histograms
332 TString* histTitle ;
327288af 333 for(Int_t k = 0; k < AliFlowConstants::kSels; k++)
9777bfcb 334 {
327288af 335 for(Int_t j = 0; j < AliFlowConstants::kHars; j++)
9777bfcb 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//-----------------------------------------------------------------------
327288af 416void AliFlowWeighter::PrintBayesian(Int_t selN) const
9777bfcb 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
327288af 421 if(selN>AliFlowConstants::kSels) { selN = 0 ; }
422 Char_t* names[AliFlowConstants::kPid] = {"e","mu","pi","k","p","d"} ;
9777bfcb 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) : " ;
327288af 428 for(Int_t ii=0;ii<AliFlowConstants::kPid;ii++)
9777bfcb 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