1 //////////////////////////////////////////////////////////////////////
5 // Author: Emanuele Simili
7 //////////////////////////////////////////////////////////////////////
8 //_____________________________________________________________
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.
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.
22 //////////////////////////////////////////////////////////////////////
24 #ifndef ALIFLOWWEIGHTER_CXX
25 #define ALIFLOWWEIGHTER_CXX
34 #include <TObjArray.h>
35 #include <TOrdCollection.h>
40 #include "AliFlowEvent.h"
41 #include "AliFlowTrack.h"
42 #include "AliFlowV0.h"
43 #include "AliFlowConstants.h"
44 #include "AliFlowSelection.h"
45 #include "AliFlowWeighter.h"
52 using namespace std; //required for resolving the 'cout' symbol
54 ClassImp(AliFlowWeighter)
55 //-----------------------------------------------------------------------
56 AliFlowWeighter::AliFlowWeighter(const AliFlowSelection* flowSelect):
57 fFlowEvent(0x0), fFlowTrack(0x0), fFlowSelect(0x0), fFlowTracks(0x0),
58 fWgtFile(0x0), fPhiWgtHistList(0x0)
60 // default constructor (selection given or default selection)
62 if(flowSelect) { fFlowSelect = new AliFlowSelection(*flowSelect) ; }
63 else { fFlowSelect = new AliFlowSelection() ; }
65 // output file (histograms)
66 fWgtFileName = "flowPhiWgt.hist.root" ;
69 //-----------------------------------------------------------------------
70 AliFlowWeighter::~AliFlowWeighter()
72 // default distructor (no actions)
74 //-----------------------------------------------------------------------
75 Bool_t AliFlowWeighter::Init()
77 // sets some defaults for the analysis
79 cout << "* FlowWeighter * - Init()" << endl ; cout << endl ;
81 // Open output files (->plots)
82 fWgtFile = new TFile(fWgtFileName.Data(), "RECREATE");
85 // counters and pointers to 0
93 //for(Int_t ii=0;ii<3;ii++) { fVertex[ii] = 0 ; }
96 fPhiBins = AliFlowConstants::kPhiBins ;
98 fPhiMax = 2*TMath::Pi() ;
101 for(int k = 0; k < AliFlowConstants::kSels; k++)
103 histTitle = new TString("Flow_BayPidMult_Sel");
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");
111 for(int j = 0; j < AliFlowConstants::kHars; j++)
114 histTitle = new TString("Flow_Phi_TPC_Sel");
116 histTitle->Append("_Har");
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");
123 histTitle = new TString("Flow_Phi_TPCplus_Sel");
125 histTitle->Append("_Har");
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");
132 histTitle = new TString("Flow_Phi_TPCminus_Sel");
134 histTitle->Append("_Har");
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");
141 histTitle = new TString("Flow_Phi_TPCcross_Sel");
143 histTitle->Append("_Har");
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");
150 // Tpc - Phi lab flattened
151 histTitle = new TString("Flow_Phi_Flat_TPC_Sel");
153 histTitle->Append("_Har");
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");
160 histTitle = new TString("Flow_Phi_Flat_TPCplus_Sel");
162 histTitle->Append("_Har");
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");
169 histTitle = new TString("Flow_Phi_Flat_TPCminus_Sel");
171 histTitle->Append("_Har");
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");
178 histTitle = new TString("Flow_Phi_Flat_TPCcross_Sel");
180 histTitle->Append("_Har");
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");
191 //-----------------------------------------------------------------------
193 //-----------------------------------------------------------------------
194 Bool_t AliFlowWeighter::Finish()
196 // Calls the method to fill wgt histograms, then saves them
197 // on the fWgtFile and closes the file .
199 cout << "* FlowWeighter * - Finish()" << endl ; cout << endl ;
204 // Write PhiWgt histograms
205 fPhiWgtHistList->Write();
206 delete fPhiWgtHistList ;
208 // Write Bayesian Weights for P.Id.
209 for(int k=0;k<AliFlowConstants::kSels;k++) { fHistFull[k].fHistBayPidMult->Write() ; }
211 // Write the AliFlowSelection object
212 fFlowSelect->Write();
217 cout << " Finish() - Wgt file closed : " << fWgtFileName.Data() << endl ; cout << endl ;
221 //-----------------------------------------------------------------------
223 //-----------------------------------------------------------------------
224 Bool_t AliFlowWeighter::WeightEvent(AliFlowEvent* fFlowEvent)
226 // Reads the AliFlowEvent (* fFlowEvent) and loops over the
227 // AliFlowTraks to fill phi histograms .
229 cout << " AliFlowWeighter::WeightEvent(" << fFlowEvent << " ) - " << fEventNumber << endl ;
230 if(!fFlowEvent) { return kFALSE ; }
232 if(fFlowSelect->Select(fFlowEvent)) // event selected - here below the ANALYSIS FLAGS are setted -
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)
240 TracksLoop(fFlowTracks) ;
244 cout << " * " << fEventNumber << " (event ID = " << fFlowEvent->EventID() << ") discarded . " << endl ;
245 delete fFlowEvent ; fFlowEvent = 0 ;
252 //-----------------------------------------------------------------------
254 //-----------------------------------------------------------------------
255 void AliFlowWeighter::TracksLoop(TObjArray* fFlowTracks)
257 // fills phi and PId histograms
259 cout << " Tracks Loop . " << endl ;
261 Float_t phi, eta, zFirstPoint ; // , zLastPoint ;
263 Char_t pid[10] = "0" ;
264 for(fTrackNumber=0;fTrackNumber<fNumberOfTracks;fTrackNumber++)
266 fFlowTrack = (AliFlowTrack*)fFlowTracks->At(fTrackNumber) ;
267 // cout << "Track n. " << fTrackNumber << endl ; fFlowTrack->Dump() ;
269 phi = fFlowTrack->Phi() ;
270 eta = fFlowTrack->Eta() ;
271 zFirstPoint = fFlowTrack->ZFirstPoint() ;
272 //zLastPoint = fFlowTrack->ZLastPoint() ;
273 strcpy(pid,fFlowTrack->Pid()) ;
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 ; }
283 // Looping over Selections and Harmonics
284 for (int k = 0; k < AliFlowConstants::kSels; k++)
286 fFlowSelect->SetSelection(k) ;
287 for (int j = 0; j < AliFlowConstants::kHars; j++)
289 fFlowSelect->SetHarmonic(j);
290 if(fFlowSelect->Select(fFlowTrack))
292 Bool_t kTpcPlus = kFALSE ;
293 Bool_t kTpcMinus = kFALSE ;
294 Bool_t kTpcAll = kFALSE ;
297 if(fFlowTrack->FitPtsTPC())
299 if(zFirstPoint >= 0. && eta > 0.) { kTpcPlus = kTRUE ; }
300 else if(zFirstPoint <= 0. && eta < 0.) { kTpcMinus = kTRUE ; }
301 else { kTpcAll = kTRUE ; }
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) ;
311 // PID Multiplicities (particle for R.P.) - just once for each selection
312 if(j==0) { fHistFull[k].fHistBayPidMult->Fill(fPidId) ; }
320 //-----------------------------------------------------------------------
321 Bool_t AliFlowWeighter::Weightening()
323 // Calculates weights, and fills PhiWgt histograms .
324 // This is called at the end of the event loop .
326 cout << " AliFlowWeighter::Weightening() " << endl ; cout << endl ;
328 // PhiWgt histogram collection
329 fPhiWgtHistList = new TOrdCollection(4*AliFlowConstants::kSels*AliFlowConstants::kHars) ;
331 // Creates PhiWgt Histograms
333 for(Int_t k = 0; k < AliFlowConstants::kSels; k++)
335 for(Int_t j = 0; j < AliFlowConstants::kHars; j++)
338 histTitle = new TString("Flow_Phi_Weight_TPCplus_Sel");
340 histTitle->Append("_Har");
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");
348 histTitle = new TString("Flow_Phi_Weight_TPCminus_Sel");
350 histTitle->Append("_Har");
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");
358 histTitle = new TString("Flow_Phi_Weight_TPCcross_Sel");
360 histTitle->Append("_Har");
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");
368 histTitle = new TString("Flow_Phi_Weight_TPC_Sel");
370 histTitle->Append("_Har");
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");
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 ;
385 for(Int_t i=0;i<fPhiBins;i++)
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.);
397 if(meanTPC==0) { cout << " Sel." << k << " , Har." << j << " : empty histogram ! " << endl ; }
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);
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);
415 //-----------------------------------------------------------------------
416 void AliFlowWeighter::PrintBayesian(Int_t selN) const
418 // Prints the normalized particle abundance of all events (selection selN).
419 // Call this at the end of the loop, just before Finish() .
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() ;
427 cout << " Sel." << selN << " particles normalized abundance (tot. " << totCount << " tracks) : " ;
428 for(Int_t ii=0;ii<AliFlowConstants::kPid;ii++)
430 bayes = (fHistFull[selN].fHistBayPidMult->GetBinContent(ii+1) / totCount) ;
431 cout << bayes << "_" << names[ii] << " ; " ;
434 else { cout << " Sel." << selN << " : empty P.Id. histogram ! " << endl ; }
439 //-----------------------------------------------------------------------