]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG/FLOW/Base/AliFlowAnalysisWithMSP.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
removed some compiler warnings
[u/mrichter/AliRoot.git] / PWG / FLOW / Base / AliFlowAnalysisWithMSP.cxx
CommitLineData
dc8ba4dd 1/*************************************************************************
2* Copyright(c) 1998-2008, ALICE Experiment at CERN, All rights reserved. *
3* *
4* Author: The ALICE Off-line Project. *
5* Contributors are mentioned in the code where appropriate. *
6* *
7* Permission to use, copy, modify and distribute this software and its *
8* documentation strictly for non-commercial purposes is hereby granted *
9* without fee, provided that the above copyright notice appears in all *
10* copies and that both the copyright notice and this permission notice *
11* appear in the supporting documentation. The authors make no claims *
12* about the suitability of this software for any purpose. It is *
13* provided "as is" without express or implied warranty. *
14**************************************************************************/
15
16#include "AliFlowAnalysisWithMSP.h"
17
18#include "AliFlowVector.h"
19#include "AliFlowMSPHistograms.h"
20#include "AliFlowEventSimple.h"
21#include "AliFlowTrackSimple.h"
22#include "AliFlowCommonConstants.h"
23#include "AliFlowCommonHist.h"
24#include "AliFlowCommonHistResults.h"
25
26#include <TList.h>
27#include <TH1D.h>
28#include <TProfile.h>
29#include <TVector2.h>
30#include <TDirectoryFile.h>
31#include <TMath.h>
32
33#include <iostream>
34#include <iomanip>
35
36AliFlowAnalysisWithMSP::AliFlowAnalysisWithMSP()
37 : TNamed(), fHarmonic(2), fNUA(kFALSE), fBookCommonHistograms(kFALSE), fCommonHist(0), fQaComponents(0),
38 fQbComponents(0), fQaQb(0), fPtUComponents(0), fEtaUComponents(0), fAllStatistics(0),
39 fPtStatistics(0), fEtaStatistics(0),
40 fIntegratedFlow(0), fDiffFlowPt(0), fDiffFlowEta(0), fFlags(0)
41{
42 // Default constructor. Intended for root IO purposes only
43}
44
45AliFlowAnalysisWithMSP::AliFlowAnalysisWithMSP(TDirectoryFile *file)
46 : TNamed(), fHarmonic(2), fNUA(kFALSE), fUseCommonConstants(kFALSE), fBookCommonHistograms(kFALSE), fCommonHist(0), fQaComponents(0),
47 fQbComponents(0), fQaQb(0), fPtUComponents(0), fEtaUComponents(0), fAllStatistics(0),
48 fPtStatistics(0), fEtaStatistics(0),
49 fIntegratedFlow(0), fDiffFlowPt(0), fDiffFlowEta(0), fFlags(0)
50{
51 ReadHistograms(file);
52}
53
54AliFlowAnalysisWithMSP::AliFlowAnalysisWithMSP(const unsigned int harmonic, const bool commonConst, const bool commonHist)
55 : TNamed(), fHarmonic(harmonic), fNUA(kFALSE), fUseCommonConstants(commonConst), fBookCommonHistograms(commonHist), fCommonHist(0), fQaComponents(0),
56 fQbComponents(0), fQaQb(0), fPtUComponents(0), fEtaUComponents(0), fAllStatistics(0),
57 fPtStatistics(0), fEtaStatistics(0),
58 fIntegratedFlow(0), fDiffFlowPt(0), fDiffFlowEta(0), fFlags(0)
59{
60 // Constructor defining the harmonic, usage of non uniform acceptance corrections and the AliFlowCommonHist() histograms
61 // This is the constructor intended for the user
62 SetNameTitle("MSP","Flow analysis with the Modified Scalar Product method");
63}
64
65AliFlowAnalysisWithMSP::~AliFlowAnalysisWithMSP()
66{
67 delete fCommonHist;
68 delete fQaComponents;
69 delete fQbComponents;
70 delete fQaQb;
71 delete fPtUComponents;
72 delete fEtaUComponents;
73 delete fAllStatistics;
74 delete fPtStatistics;
75 delete fEtaStatistics;
76 delete fIntegratedFlow;
77 delete fDiffFlowPt;
78 delete fDiffFlowEta;
79 delete fFlags;
80}
81
82void AliFlowAnalysisWithMSP::Init()
83{
84 // Create all output objects. Memory consumption can be reduced by switching off some of the control histograms.
85 //
86 delete fCommonHist; fCommonHist=0; // Delete existing histograms
87 delete fQaComponents; fQaComponents=0;
88 delete fQbComponents; fQbComponents=0;
89 delete fQaQb; fQaQb=0;
90 delete fPtUComponents; fPtUComponents=0;
91 delete fEtaUComponents; fEtaUComponents=0;
92 delete fAllStatistics; fAllStatistics=0;
93 delete fPtStatistics; fPtStatistics=0;
94 delete fEtaStatistics; fEtaStatistics=0;
95 delete fIntegratedFlow; fIntegratedFlow=0;
96 delete fDiffFlowPt; fDiffFlowPt=0;
97 delete fDiffFlowEta; fDiffFlowEta=0;
98 delete fFlags; fFlags=0;
99
100 // Default binning for histograms
101 // TODO: allow variable binning
102
103 // Defaults
104 int nBinsPt=100;
105 double ptMin=0;
106 double ptMax=10.0;
107 int nBinsEta=100;
108 double etaMin=-0.8;
109 double etaMax=+0.8;
110
111 if( fUseCommonConstants || fBookCommonHistograms ) {
112 AliFlowCommonConstants *c=AliFlowCommonConstants::GetMaster();
113 nBinsPt=c->GetNbinsPt();
114 ptMin=c->GetPtMin();
115 ptMax=c->GetPtMax();
116 nBinsEta=c->GetNbinsEta();
117 etaMin=c->GetEtaMin();
118 etaMax=c->GetEtaMax();
119 }
120
121 if( fBookCommonHistograms ) { // Use the common constants from the flow package for histogram definitions
122 std::cerr << "AliFlowAnalysisWithMSP::Init() Creating common histograms" << std::endl;
123 fCommonHist=new AliFlowCommonHist("AliFlowCommonHist_MSP","AliFlowCommonHist",kTRUE);
124 }
125
126 std::cerr << "AliFlowAnalysisWithMSP::Init() creating MSPHistograms" << std::endl;
127 // Averages for NUA corrections:
128 fQaComponents=new AliFlowMSPHistograms(2,"QaComponents",1,0.5,1.5); // QaX, QaY
129 fQaComponents->SetXName("all");
130 fQaComponents->SetVarName("QaX",0);
131 fQaComponents->SetVarName("QaY",1);
132 fQbComponents=new AliFlowMSPHistograms(2,"QbComponents",1,0.5,1.5); // QbX, QbY
133 fQbComponents->SetXName("all");
134 fQbComponents->SetVarName("QbX",0);
135 fQbComponents->SetVarName("QbY",1);
136 fQaQb=new AliFlowMSPHistograms(1,"QaQb",1,0.5,1.5); // QaQb
137 fQaQb->SetXName("all");
138 fQaQb->SetVarName("QaQb",0);
139 fPtUComponents=new AliFlowMSPHistograms(2,"PtUComponents",nBinsPt,ptMin,ptMax); // ux(pt), uy(pt)
140 fPtUComponents->SetXName("pt");
141 fPtUComponents->SetVarName("ux",0);
142 fPtUComponents->SetVarName("uy",1);
143 fEtaUComponents=new AliFlowMSPHistograms(2,"EtaUComponents",nBinsEta,etaMin,etaMax); // ux(eta), uy(eta)
144 fEtaUComponents->SetXName("eta");
145 fEtaUComponents->SetVarName("ux",0);
146 fEtaUComponents->SetVarName("uy",1);
147
148 // Correlation terms
149 fAllStatistics=new AliFlowMSPHistograms(3,"AllStatistics",1,0.5,1.5); // terms integrated over pt and eta
150 fAllStatistics->SetXName("all");
151 fAllStatistics->SetVarName("uQa",0);
152 fAllStatistics->SetVarName("uQb",1);
153 fAllStatistics->SetVarName("QaQb",2);
154 fPtStatistics=new AliFlowMSPHistograms(3,"PtStatistics",nBinsPt,ptMin,ptMax); // terms per pt bin
155 fPtStatistics->SetXName("pt");
156 fPtStatistics->SetVarName("uQa",0);
157 fPtStatistics->SetVarName("uQb",1);
158 fPtStatistics->SetVarName("QaQb",2);
159 fEtaStatistics=new AliFlowMSPHistograms(3,"EtaStatistics",nBinsEta,etaMin,etaMax); // terms per eta bin
160 fEtaStatistics->SetXName("eta");
161 fEtaStatistics->SetVarName("uQa",0);
162 fEtaStatistics->SetVarName("uQb",1);
163 fEtaStatistics->SetVarName("QaQb",2);
164
165 fIntegratedFlow=0; // Created in Finish()
166 fDiffFlowPt=0; // Created in Finish
167 fDiffFlowEta=0; // Created in Finish
168
169 fFlags=new TProfile("Flags","Flags for AliFlowAnalysisWithMSP",10,0.5,10.5,"s");
170 fFlags->Fill("Harmonic",fHarmonic); // bin 1
171}
172
173void AliFlowAnalysisWithMSP::Make(AliFlowEventSimple *event)
174{
175 // Analyze one event. The modified scalar product method estimates flow using the formula:
176 // BEGIN_LATEX v_2(MSP) = \sqrt( uQ^{a} uQ^{b} / (Q^{a}Q^{b}) ) END_LATEX
177 // The Q vectors are calculated for the harmonic set previously: fHarmonic.
178 // Make(event) does not use the new operator, thus avoiding memory leaks in a simple way
179 // Depending on the compiler about 200 bytes may be pushed/popped on the stack for local variables.
180 //
181
182 if( !event ) return; // Protect against running without events. Can this happen???
183
184 AliFlowVector flowVectors[2]; // Calculate the two subevent Q vectors:
185 event->Get2Qsub(flowVectors, fHarmonic); // No phi, pt or eta weights implemented yet
186
187 AliFlowVector &Qa=flowVectors[0]; // Define some mnemonics for the subevent flow vectors:
188 AliFlowVector &Qb=flowVectors[1];
189
190 const double QaW=Qa.GetMult(); // Weight for Qa and combinations of Qa
191 const double QbW=Qb.GetMult(); // Weight for Qb and combinations of Qb
192
193 if( QaW<2 || QbW<2 ) return; // Require at least 2 particles in each subevent
194
195 if(fCommonHist)fCommonHist->FillControlHistograms(event); // Standard for all flow analysis
196
197
198 const double qaxy[]={Qa.X()/QaW,Qa.Y()/QaW}; // Two variables expected
199 const double wqaxy[]={QaW,QaW};
200 fQaComponents->Fill(1, qaxy, wqaxy); // only one bin (all pt)
201
202 const double qbxy[]={Qb.X()/QbW,Qb.Y()/QbW}; // Two variables expected
203 const double wqbxy[]={QbW,QbW};
204 fQbComponents->Fill(1, qbxy, wqbxy); // only one bin (all pt)
205
206 const double QaQbW=QaW*QbW;
207 const double weightedQaQb = (Qa*Qb)/QaQbW; // Scalar product of subevent Q vectors with weight
208 fQaQb->Fill(1,&weightedQaQb,&QaQbW); // Average of QaQb per event
209
210
211 int iTrack=0;
212 while( AliFlowTrackSimple *track=event->GetTrack(iTrack++) ) {// Loop over the tracks in the event
213 // Get the track vector
214 if(! track->InPOISelection() ) continue;
215 const double trackWeight=track->Weight();
216 const double phi=track->Phi();
217 const double pt=track->Pt();
218 const double eta=track->Eta();
219
220 AliFlowVector u;
221 u.SetMagPhi(1, fHarmonic*phi, trackWeight);
222 u.SetMult(1);
223
224 // Remove track from subevent a
225 AliFlowVector mQa(Qa); // Initialize with Qa flow vector
226 if( track->InSubevent(0) ) {
227 mQa-=u; // Should introduce phi weights here
228 }
229
230 // Remove track from subevent b
231 AliFlowVector mQb(Qb); // Initialize with Qb flow vector
232 if( track->InSubevent(1) ) {
233 mQb-=u; // Should introduce phi weights here
234 }
235
236 const double uQaW = mQa.GetMult(); // Weight is multiplicity of Q vector
237 const double uQbW = mQb.GetMult();
238 const double uQa=u*mQa; // Correlate POI with subevent
239 const double uQb=u*mQb;
240
241 const double uxy[]={u.X(),u.Y()};
242 const double wxy[]={1.0,1.0};
243 fPtUComponents->Fill(pt, uxy, wxy); // vs pt
244 fEtaUComponents->Fill(eta, uxy, wxy); // vs eta
245
246 const double par[]={uQa/uQaW, uQb/uQbW, weightedQaQb};
247 const double wgt[]={uQaW, uQbW, QaQbW};
248 fAllStatistics->Fill(1, par, wgt );
249 fPtStatistics->Fill(pt, par, wgt );
250 fEtaStatistics->Fill(eta, par, wgt );
251 }
252}
253
254void AliFlowAnalysisWithMSP::Finish()
255{
256 // Calculate the final result from the stored correlations.
257 // The NUA corrections are applied if the flag fNUA was set before the call to Finish()
258 // If the output histograms already exist then they are replaced by the newly calculated result
259
260 Print(); // Print a summary of the NUA terms and integrated flow
261 // TODO: Create result histograms for integrated flow and store the flags to be able to restore the initial state from histograms only
262
263 // Create result histograms
264 fFlags->Fill("NUA",fNUA);
265 delete fIntegratedFlow; fIntegratedFlow=0; // First delete existing results (if any)
266 delete fDiffFlowPt; fDiffFlowPt=0;
267 delete fDiffFlowEta; fDiffFlowEta=0;
268
269 fIntegratedFlow=new TH1D("IntegratedFlow","Integrated flow results",10,0.5,10.5);
270 double vn, vnerror;
271 Calculate(vn, vnerror, fAllStatistics, fPtUComponents, 0, 0);
272 fIntegratedFlow->SetBinContent(1,vn);
273 fIntegratedFlow->SetBinError(1,vnerror);
274 fIntegratedFlow->GetXaxis()->SetBinLabel(1,"POI");
275 Calculate(vn, vnerror, fAllStatistics, fPtUComponents, 0, 1);
276 fIntegratedFlow->SetBinContent(2,vn);
277 fIntegratedFlow->SetBinError(2,vnerror);
278 fIntegratedFlow->GetXaxis()->SetBinLabel(2,"A");
279 Calculate(vn, vnerror, fAllStatistics, fPtUComponents, 0, 2);
280 fIntegratedFlow->SetBinContent(3,vn);
281 fIntegratedFlow->SetBinError(3,vnerror);
282 fIntegratedFlow->GetXaxis()->SetBinLabel(3,"B");
283
284 int nbinspt=fPtStatistics->Nbins();
285 double ptlow=fPtStatistics->XLow();
286 double pthigh=fPtStatistics->XHigh();
287 fDiffFlowPt = new TH1D("DiffFlowPt","flow of POI vs pt",nbinspt,ptlow,pthigh);
288 fDiffFlowPt->SetDirectory(0); // Do not automatically store in file
289 fDiffFlowPt->SetStats(kFALSE);
290 fDiffFlowPt->SetXTitle("p_{t}");
291 fDiffFlowPt->SetYTitle(Form("v_{%d}",fHarmonic));
292
293 for(int i=1; i<=nbinspt; ++i) {
294 double vn=0;
295 double evn=0;
296 if( Calculate(vn, evn, fPtStatistics, fPtUComponents, i) && evn<1 ) {
297 fDiffFlowPt->SetBinContent(i,vn);
298 fDiffFlowPt->SetBinError(i,evn);
299 }
300 }
301
302 int nbinseta=fEtaStatistics->Nbins();
303 double etalow=fEtaStatistics->XLow();
304 double etahigh=fEtaStatistics->XHigh();
305 fDiffFlowEta= new TH1D("DiffFlowEta","flow of POI vs #eta",nbinseta,etalow,etahigh);
306 fDiffFlowEta->SetDirectory(0); // Do not automatically store in file
307 fDiffFlowEta->SetStats(kFALSE);
308 fDiffFlowEta->SetXTitle("p_{t}");
309 fDiffFlowEta->SetYTitle(Form("v_{%d}",fHarmonic));
310
311 for(int i=1; i<=nbinseta; ++i) {
312 double vn=0;
313 double evn=0;
314 if( Calculate(vn, evn, fEtaStatistics, fEtaUComponents, i) && evn<1 ) {
315 fDiffFlowEta->SetBinContent(i,vn);
316 fDiffFlowEta->SetBinError(i,evn);
317 }
318 }
319}
320
321void AliFlowAnalysisWithMSP::WriteHistograms(TDirectoryFile *file) const
322{
323 //file->Write(file->GetName(), TObject::kSingleKey); // Make sure the directory itself is written
324
325 if(fCommonHist) file->WriteTObject(fCommonHist);
326
327 file->WriteTObject(fQaComponents,0,"Overwrite"); // Averages of Qa components per event
328 file->WriteTObject(fQbComponents,0,"Overwrite"); // Averages of Qb components per event
329 file->WriteTObject(fQaQb,0,"Overwrite"); // Average of QaQb per event
330 file->WriteTObject(fPtUComponents,0,"Overwrite"); // u components vs pt
331 file->WriteTObject(fEtaUComponents,0,"Overwrite"); // u components vs eta
332 file->WriteTObject(fAllStatistics,0,"Overwrite"); // Integrated uQa, uQb and QaQa
333 file->WriteTObject(fPtStatistics,0,"Overwrite"); // uQa, uQb and QaQb vs pt
334 file->WriteTObject(fEtaStatistics,0,"Overwrite"); // uQa, uQb and QaQb vs eta
335
336 if( fIntegratedFlow ) file->WriteTObject(fIntegratedFlow,0,"Overwrite"); // Integrated flow for POI and subevents
337 if( fDiffFlowPt ) file->WriteTObject(fDiffFlowPt,0,"Overwrite"); // Differential flow vs pt if calculated
338 if( fDiffFlowEta ) file->WriteTObject(fDiffFlowEta,0,"Overwrite"); // Differential flow vs eta if calculated
339
340 file->WriteTObject(fFlags,0,"Overwrite"); // fHarmonic, fNUA
341
342 file->WriteKeys(); // Make sure it happens now
343}
344
345void AliFlowAnalysisWithMSP::WriteCommonResults(TDirectoryFile *file) const
346{
347 // Copy the results to a AliFlowCommonHistResults() class and write to file
348 // If the results were not calculated again then Finish() is called to generate them
349 // The global AliFlowCommonConstants() is modified to adapt to the binning used in this analysis
350
351 AliFlowCommonConstants *c=AliFlowCommonConstants::GetMaster();
352 int nBinsPt=fPtStatistics->Nbins();
353 double ptMin=fPtStatistics->XLow();
354 double ptMax=fPtStatistics->XHigh();
355
356 int nBinsEta=fEtaStatistics->Nbins();
357 double etaMin=fEtaStatistics->XLow();
358 double etaMax=fEtaStatistics->XHigh();
359
360 c->SetNbinsPt(nBinsPt);
361 c->SetPtMin(ptMin);
362 c->SetPtMax(ptMax);
363 c->SetNbinsEta(nBinsEta);
364 c->SetEtaMin(etaMin);
365 c->SetEtaMax(etaMax);
366
367 bool oldAddStatus=TH1::AddDirectoryStatus();
368 TH1::AddDirectory(kFALSE);
369 AliFlowCommonHistResults *h=new AliFlowCommonHistResults("AliFlowCommonHistResults_MSP","AliFlowCommonHistResults from the MSP method",fHarmonic);
370
371 double ivn, ivnerror;
372 Calculate(ivn, ivnerror, fAllStatistics, fPtUComponents, 0, 0);
373 h->FillIntegratedFlowPOI(ivn, ivnerror);
374
375 for(int bin=1; bin<=nBinsPt; ++bin) {
376 double vn=0;
377 double evn=0;
378 if( Calculate(vn, evn, fPtStatistics, fPtUComponents, bin) && evn>0 ) {
379 h->FillDifferentialFlowPtPOI(bin, vn, evn);
380 }
381 }
382
383 for(int bin=1; bin<=nBinsEta; ++bin) {
384 double vn=0;
385 double evn=0;
386 if( Calculate(vn, evn, fEtaStatistics, fEtaUComponents, bin) && evn>0 ) {
387 h->FillDifferentialFlowEtaPOI(bin, vn, evn);
388 }
389 }
390
391 file->WriteTObject(h,0,"Overwrite");
392
393 TH1::AddDirectory(oldAddStatus);
394}
395
396
397
398void AliFlowAnalysisWithMSP::Print(const Option_t *opt)const
399{
400 std::cout << setprecision(3);
401
402 std::cout << "****************************************************" << std::endl;
403 std::cout << " Integrated flow from Modified Scalar Product " << std::endl;
404 std::cout << " " << std::endl;
405
406 double vn=0;
407 double vnerror=0;
408
409 Calculate(vn, vnerror, fAllStatistics, fPtUComponents, 0, 0);
410 std::cout << "v" << fHarmonic << " for POI : " << setw(10) << vn << " +- " << setw(8) << vnerror << std::endl;
411 Calculate(vn, vnerror, fAllStatistics, fPtUComponents, 0, 1);
412 std::cout << "v" << fHarmonic << " for subevent A: " << setw(10) << vn << " +- " << setw(8) << vnerror << std::endl;
413 Calculate(vn, vnerror, fAllStatistics, fPtUComponents, 0, 2);
414 std::cout << "v" << fHarmonic << " for subevent B: " << setw(10) << vn << " +- " << setw(8) << vnerror << std::endl;
415 std::cout << std::endl;
416
417 std::cout << "NUA terms: " << (fNUA?"(applied)":"(NOT applied)") << std::endl;
418
419 const double ux=fPtUComponents->Average(0); // Average over all bins
420 const double eux=TMath::Sqrt(fPtUComponents->Variance(0));
421 std::cout << "<ux> " << setw(12) << ux << " +- " << setw(12) << eux << (TMath::Abs(ux)<2*eux?" NOT significant ":" ") << std::endl;
422 const double ux0eta=fEtaUComponents->Average(fEtaUComponents->FindBin(0.0),0);
423 const double eux0eta=TMath::Sqrt(fEtaUComponents->Variance(fEtaUComponents->FindBin(0.0),0));
424 std::cout << "<ux> eta=0 " << setw(12) << ux0eta << " +- " << setw(12) << eux0eta << (TMath::Abs(ux0eta)<2*eux0eta?" NOT significant":" ") << std::endl;
425 const double ux0pt=fPtUComponents->Average(fPtUComponents->FindBin(1.0),0);
426 const double eux0pt=TMath::Sqrt(fPtUComponents->Variance(fPtUComponents->FindBin(1.0),0));
427 std::cout << "<ux> pt=1 " << setw(12) << ux0pt << " +- " << setw(12) << eux0pt << (TMath::Abs(ux0pt)<2*eux0pt?" NOT significant":" ") << std::endl;
428
429 const double uy=fPtUComponents->Average(0); // Average over all bins
430 const double euy=TMath::Sqrt(fPtUComponents->Variance(0));
431 std::cout << "<uy> " << setw(12) << uy << " +- " << setw(12) << euy << (TMath::Abs(uy)<2*euy?" NOT significant ":" ") << std::endl;;
432 const double uy0eta=fEtaUComponents->Average(fEtaUComponents->FindBin(0.0),1);
433 const double euy0eta=TMath::Sqrt(fEtaUComponents->Variance(fEtaUComponents->FindBin(0.0),1));
434 std::cout << "<uy> eta=0 " << setw(12) << uy0eta << " +- " << setw(12) << euy0eta << (TMath::Abs(uy0eta)<2*euy0eta?" NOT significant ":" ") << std::endl;
435 const double uy0pt=fPtUComponents->Average(fPtUComponents->FindBin(1.0),1);
436 const double euy0pt=TMath::Sqrt(fPtUComponents->Variance(fPtUComponents->FindBin(1.0),1));
437 std::cout << "<uy> pt=1 " << setw(12) << uy0pt << " +- " << setw(12) << euy0pt << (TMath::Abs(uy0pt)<2*euy0pt?" NOT significant ":" ") << std::endl;
438
439 std::cout << "<QaX> " << setw(12) << fQaComponents->Average(0) << " +- " << setw(12) << TMath::Sqrt(fQaComponents->Variance(0)) << std::endl;
440 std::cout << "<QaY> " << setw(12) << fQaComponents->Average(1) << " +- " << setw(12) << TMath::Sqrt(fQaComponents->Variance(1)) << std::endl;
441 std::cout << "<QbX> " << setw(12) << fQbComponents->Average(0) << " +- " << setw(12) << TMath::Sqrt(fQbComponents->Variance(0)) << std::endl;
442 std::cout << "<QbY> " << setw(12) << fQbComponents->Average(1) << " +- " << setw(12) << TMath::Sqrt(fQbComponents->Variance(1)) << std::endl;
443 std::cout << "<QaQb> " << setw(12) << fQaQb->Average(0) << " +- " << setw(12) << TMath::Sqrt(fQbComponents->Variance(0)) << std::endl;
444 std::cout << std::endl;
445
446 std::cout << "Covariance matrix: " << std::endl;
447 std::cout << " " << setw(12) << "uQa" << setw(12) << "uQb" << setw(12) << "QaQb" << std::endl;
448 std::cout << "uQa " << setw(12) << fAllStatistics->Covariance(0,0) << std::endl;
449 std::cout << "uQb " << setw(12) << fAllStatistics->Covariance(1,0) << setw(12) << fAllStatistics->Covariance(1,1) << std::endl;
450 std::cout << "QaQb " << setw(12) << fAllStatistics->Covariance(2,0) << setw(12) << fAllStatistics->Covariance(2,1) << setw(12) << fQaQb->Variance(0) << std::endl;
451 std::cout << "****************************************************" << std::endl;
452 std::cout << std::endl;
453}
454
455
456// private functions --------------------------------------------------------------------------------------
457bool AliFlowAnalysisWithMSP::Calculate(double &vn, double &vnerror, const AliFlowMSPHistograms *hist, const AliFlowMSPHistograms *components, const int bin, const int poi) const
458{
459 // Get all averages and correlations need for the flow calculation
460 double uQa=hist->Average(bin,0); // <<uQa>>
461 double VuQa=hist->Variance(bin,0); // Var(<<uQa>>)
462 double uQb=hist->Average(bin,1); // <<uQb>>
463 double VuQb=hist->Variance(bin,1); // Var(<<uQb>>)
464 double QaQb=fQaQb->Average(1,0); // <QaQb> Should not be taken from hist(bin) buit from fQaQa because there QaQb is entered multiple times: <<QaQb>>!!
465 double VQaQb=fQaQb->Variance(1,0); // V(<QaQb>)
466 double CuQauQb=hist->Covariance(bin,0,1); // Cov(<<uQa>>,<<uQb>>)
467 double CuQaQaQb=hist->Covariance(bin,0,2); // Cov(<<uQa>>,<QaQb>)
468 double CuQbQaQb=hist->Covariance(bin,1,2); // Cov(<<uQb>>,<QaQb>)
469
470 if( fNUA && components ) {
471 // Apply NUA correction to QaQb.
472 double QaX=fQaComponents->Average(0);
473 double QaY=fQaComponents->Average(1);
474 double QbX=fQbComponents->Average(0);
475 double QbY=fQbComponents->Average(1);
476
477 QaQb=QaQb-QaX*QbX-QaY*QbY;
478
479 // Apply NUA to uQa and uQb (per bin)
480 double uX=components->Average(bin,0); // bin 0 is integrated over all bins
481 double uY=components->Average(bin,1);
482
483 uQa = uQa - uX*QaX - uY*QaY;
484 uQb = uQb - uX*QbX - uY*QbY;
485 // Error calculation not fully modified: only above terms, the spread in <<u>> and <Qa> and <Qb> are not used
486 // therefore this should only be applied if significant!
487 // TODO: Check if not fully NUA correcting the error calculation is justified (but this is compatible with the original SP method)!
488 }
489
490 // Some sanity checks:
491 if( uQa*uQb*QaQb <= 0 ) { // Catch imaginary results
492 vn=-99;
493 vnerror=-99;
494 return false;
495 }
496
497 // Decent results, calculate, print and store
498 // TODO: The three cases are cyclic permutations of u, Qa, Qb so there should be a simpler way to do this.
499 // However the difficulty is that we deal here with uQa, uQb and QaQb
500 switch (poi) {
501 case 1: // Subevent A reference flow
502 {
503 if( TMath::Abs(uQb) < 1e-30*TMath::Abs(uQa*QaQb) ) { // Protect against infinity
504 vn=0;
505 vnerror=-1;
506 return false;
507 }
508 double vnB = TMath::Sqrt( uQa*QaQb / uQb ); // vn
509 double VvnB = QaQb*VuQa/(4*uQa*uQb) // Variance of vn
510 + uQa*VQaQb/(4*QaQb*uQb)
511 + uQa*QaQb*VuQb/(4*TMath::Power(uQb,3))
512 + CuQaQaQb/(2*uQb)
513 - QaQb*CuQauQb/(2*TMath::Power(uQb,2))
514 - uQa*CuQbQaQb/(2*TMath::Power(uQb,2));
515 vn=vnB;
516 if( VvnB<0 ) {
517 vnerror=VvnB;
518 return false;
519 }
520 vnerror=TMath::Sqrt(VvnB);
521 }
522 break;
523 case 2: // Subevent B reference flow
524 {
525 if( TMath::Abs(uQa) < 1e-30*TMath::Abs(uQb*QaQb) ) { // Protect against infinity
526 vn=0;
527 vnerror=-1;
528 return false;
529 }
530 double vnA = TMath::Sqrt( uQb*QaQb / uQa ); // vn
531 double VvnA = uQb*VQaQb/(4*QaQb*uQa) // Variance of vn
532 + QaQb*VuQb/(4*uQb*uQa)
533 + QaQb*uQb*VuQa/(4*TMath::Power(uQa,3))
534 + CuQbQaQb/(2*uQa)
535 - uQb*CuQaQaQb/(2*TMath::Power(uQa,2))
536 - uQa*CuQauQb/(2*TMath::Power(uQa,2));
537 vn=vnA;
538 if( VvnA<0 ) {
539 vnerror=VvnA;
540 return false;
541 }
542 vnerror=TMath::Sqrt(VvnA);
543 }
544 break;
545 default: // POI flow
546 {
547 if( TMath::Abs(QaQb) < 1e-30*TMath::Abs(uQa*uQb) ) { // Catch infinity
548 vn=0;
549 vnerror=-1;
550 return false;
551 }
552 double vnP = TMath::Sqrt( uQa*uQb / QaQb ); // vn
553 if( TMath::Abs(uQa*QaQb) < 1e-30*TMath::Abs(uQb*VuQa)
554 || TMath::Abs(uQb*QaQb) < 1e-30*TMath::Abs(uQa*VuQb)
555 || TMath::Abs(QaQb) < 1e-30*TMath::Abs(uQa*uQb*VQaQb)
556 || TMath::Abs(QaQb) < 1e-30*TMath::Abs(CuQauQb)
557 || TMath::Abs(QaQb) < 1e-30*TMath::Abs(uQb*CuQaQaQb)
558 || TMath::Abs(QaQb) < 1e-30*TMath::Abs(uQa*CuQbQaQb)
559 ){
560 vnerror=-98;
561 return false;
562 }
563 double VvnP = uQb*VuQa/(4*uQa*QaQb) // Variance of vn
564 + uQa*VuQb/(4*uQb*QaQb)
565 + uQa*uQb*VQaQb/(4*TMath::Power(QaQb,3))
566 + CuQauQb/(2*QaQb)
567 - uQb*CuQaQaQb/(2*TMath::Power(QaQb,2))
568 - uQa*CuQbQaQb/(2*TMath::Power(QaQb,2));
569 vn=vnP;
570 if( VvnP<0 ) {
571 vnerror=VvnP;
572 return false;
573 }
574 vnerror=TMath::Sqrt(VvnP);
575 }
576 } // Switch between POI and subevents
577
578 return (vnerror>=0);
579}
580
581void AliFlowAnalysisWithMSP::ReadHistograms(TDirectoryFile *file)
582{
583 delete fCommonHist; fCommonHist=0; // Delete existing histograms
584 delete fQaComponents; fQaComponents=0;
585 delete fQbComponents; fQbComponents=0;
586 delete fQaQb; fQaQb=0;
587 delete fPtUComponents; fPtUComponents=0;
588 delete fEtaUComponents; fEtaUComponents=0;
589 delete fAllStatistics; fAllStatistics=0;
590 delete fPtStatistics; fPtStatistics=0;
591 delete fEtaStatistics; fEtaStatistics=0;
592 delete fIntegratedFlow; fIntegratedFlow=0;
593 delete fDiffFlowPt; fDiffFlowPt=0;
594 delete fDiffFlowEta; fDiffFlowEta=0;
595 delete fFlags; fFlags=0;
596
597 file->GetObject("QaComponents",fQaComponents);
598 file->GetObject("QbComponents",fQbComponents);
599 file->GetObject("QaQb",fQaQb);
600 file->GetObject("PtUComponents",fPtUComponents);
601 file->GetObject("EtaUComponents",fEtaUComponents);
602 file->GetObject("AllStatistics",fAllStatistics);
603 file->GetObject("PtStatistics",fPtStatistics);
604 file->GetObject("EtaStatistics",fEtaStatistics);
605 if( !fQaComponents || !fQbComponents || !fQaQb || !fPtUComponents || !fEtaUComponents || !fAllStatistics || !fPtStatistics || !fEtaStatistics ) {
606 std::cerr << "AliFlowAnalysisWithMSP::ReadHistograms() : One or more histograms were not read correctly from " << file->GetPath() << std::endl;
607 }
608
609 // Optional histograms, may return a zero pointer
610 file->GetObject("AliFlowCommonHist_MSP",fCommonHist); // The AliFlowCommonHist is optional
611 file->GetObject("IntegratedFlow",fIntegratedFlow);
612 file->GetObject("DiffFlowPt",fDiffFlowPt);
613 file->GetObject("DiffFlowEta",fDiffFlowEta);
614
615 file->GetObject("Flags",fFlags);
616
617 if( !fFlags ){
618 std::cerr << "AliFlowAnalysisWithMSP::ReadHistograms() : Flags histogram not found, using defaults" << std::endl;
619 fHarmonic=2;
620 fNUA=false;
621 }else{
622 fHarmonic=fFlags->GetBinContent(1);
623 double harmonicSpread=fFlags->GetBinError(1);
624 if( harmonicSpread!=0 ) {
625 std::cerr << "These histograms seem to be merged from analysis with different Harmonics. Results are invalid!" << std::endl;
626 delete fIntegratedFlow; fIntegratedFlow=0;
627 delete fDiffFlowPt; fDiffFlowPt=0;
628 delete fDiffFlowEta; fDiffFlowEta=0;
629 }
630 fNUA=fFlags->GetBinContent(2); // Mixing NUA does not matter since it needs to be recalculated anyway
631 double nuaSpread=fFlags->GetBinError(2);
632 if( nuaSpread!=0 ) {
633 std::cerr << "These histograms seem to be merged from analysis with different NUA corrections. Results removed" << std::endl;
634 delete fIntegratedFlow; fIntegratedFlow=0;
635 delete fDiffFlowPt; fDiffFlowPt=0;
636 delete fDiffFlowEta; fDiffFlowEta=0;
637 }
638 }
639 fBookCommonHistograms=(fCommonHist!=0);
640}