]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG2/FLOW/AliFlowCommon/AliFlowAnalysisWithScalarProduct.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
coverty fixes
[u/mrichter/AliRoot.git] / PWG2 / FLOW / AliFlowCommon / AliFlowAnalysisWithScalarProduct.cxx
CommitLineData
8d312f00 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#define AliFlowAnalysisWithScalarProduct_cxx
17
929098e4 18#include "Riostream.h"
19#include "TFile.h"
d7eb18ec 20#include "TList.h"
8d312f00 21#include "TMath.h"
22#include "TProfile.h"
23#include "TVector2.h"
a554b203 24#include "TH1D.h"
25#include "TH2D.h"
8d312f00 26
929098e4 27#include "AliFlowCommonConstants.h"
8d312f00 28#include "AliFlowEventSimple.h"
a93de0f0 29#include "AliFlowVector.h"
8d312f00 30#include "AliFlowTrackSimple.h"
31#include "AliFlowCommonHist.h"
395fadba 32#include "AliFlowCommonHistResults.h"
8d312f00 33#include "AliFlowAnalysisWithScalarProduct.h"
8d312f00 34
12dc476e 35//////////////////////////////////////////////////////////////////////////////
8d312f00 36// AliFlowAnalysisWithScalarProduct:
37// Description:
38// Maker to analyze Flow with the Scalar product method.
39//
12dc476e 40// authors: N. van der Kolk (kolk@nikhef.nl), A. Bilandzic (anteb@nikhef.nl)
41//////////////////////////////////////////////////////////////////////////////
8d312f00 42
43ClassImp(AliFlowAnalysisWithScalarProduct)
44
45 //-----------------------------------------------------------------------
12dc476e 46 AliFlowAnalysisWithScalarProduct::AliFlowAnalysisWithScalarProduct():
e35ddff0 47 fEventNumber(0),
48 fDebug(kFALSE),
a93de0f0 49 fApplyCorrectionForNUA(kFALSE),
50fe33aa 50 fHarmonic(2),
a554b203 51 fRelDiffMsub(1.),
29195b69 52 fWeightsList(NULL),
53 fUsePhiWeights(kFALSE),
04c6b875 54 fPhiWeightsSub0(NULL),
55 fPhiWeightsSub1(NULL),
29ecccee 56 fHistProFlags(NULL),
395fadba 57 fHistProUQetaRP(NULL),
58 fHistProUQetaPOI(NULL),
1b302085 59 fHistProUQetaAllEventsPOI(NULL),
395fadba 60 fHistProUQPtRP(NULL),
61 fHistProUQPtPOI(NULL),
1b302085 62 fHistProUQPtAllEventsPOI(NULL),
63 fHistProQNorm(NULL),
e4cecfa0 64 fHistProQaQb(NULL),
a554b203 65 fHistProQaQbNorm(NULL),
29ecccee 66 fHistProQaQbReImNorm(NULL),
67 fHistProNonIsotropicTermsQ(NULL),
12dc476e 68 fHistSumOfLinearWeights(NULL),
69 fHistSumOfQuadraticWeights(NULL),
70 fHistProUQQaQbPtRP(NULL),
71 fHistProUQQaQbEtaRP(NULL),
72 fHistProUQQaQbPtPOI(NULL),
73 fHistProUQQaQbEtaPOI(NULL),
1b302085 74 fCommonHistsSP(NULL),
75 fCommonHistsResSP(NULL),
76 fCommonHistsmuQ(NULL),
77 fHistQNorm(NULL),
a554b203 78 fHistQaQb(NULL),
79 fHistQaQbNorm(NULL),
1b302085 80 fHistQNormvsQaQbNorm(NULL),
a554b203 81 fHistQaQbCos(NULL),
82 fHistResolution(NULL),
83 fHistQaNorm(NULL),
84 fHistQaNormvsMa(NULL),
85 fHistQbNorm(NULL),
86 fHistQbNormvsMb(NULL),
a93de0f0 87 fHistMavsMb(NULL),
88 fHistList(NULL)
8d312f00 89{
8d312f00 90 // Constructor.
29195b69 91 fWeightsList = new TList();
e35ddff0 92 fHistList = new TList();
12dc476e 93
94 // Initialize arrays:
95 for(Int_t i=0;i<3;i++)
96 {
97 fHistSumOfWeightsPtRP[i] = NULL;
98 fHistSumOfWeightsEtaRP[i] = NULL;
99 fHistSumOfWeightsPtPOI[i] = NULL;
100 fHistSumOfWeightsEtaPOI[i] = NULL;
101 }
29ecccee 102 for(Int_t rp=0;rp<2;rp++)
103 {
104 for(Int_t pe=0;pe<2;pe++)
105 {
106 for(Int_t sc=0;sc<2;sc++)
107 {
108 fHistProNonIsotropicTermsU[rp][pe][sc] = NULL;
109 }
110 }
111 }
8d312f00 112}
113 //-----------------------------------------------------------------------
8d312f00 114 AliFlowAnalysisWithScalarProduct::~AliFlowAnalysisWithScalarProduct()
115 {
116 //destructor
29195b69 117 delete fWeightsList;
e35ddff0 118 delete fHistList;
8d312f00 119 }
120
1dfa3c16 121//-----------------------------------------------------------------------
1dfa3c16 122void AliFlowAnalysisWithScalarProduct::WriteHistograms(TString* outputFileName)
123{
124 //store the final results in output .root file
125
126 TFile *output = new TFile(outputFileName->Data(),"RECREATE");
0fe80f88 127 //output->WriteObject(fHistList, "cobjSP","SingleKey");
128 fHistList->SetName("cobjSP");
9455e15e 129 fHistList->SetOwner(kTRUE);
0fe80f88 130 fHistList->Write(fHistList->GetName(), TObject::kSingleKey);
1dfa3c16 131 delete output;
132}
133
b0fda271 134//-----------------------------------------------------------------------
b0fda271 135void AliFlowAnalysisWithScalarProduct::WriteHistograms(TString outputFileName)
136{
137 //store the final results in output .root file
138
139 TFile *output = new TFile(outputFileName.Data(),"RECREATE");
0fe80f88 140 //output->WriteObject(fHistList, "cobjSP","SingleKey");
141 fHistList->SetName("cobjSP");
9455e15e 142 fHistList->SetOwner(kTRUE);
0fe80f88 143 fHistList->Write(fHistList->GetName(), TObject::kSingleKey);
b0fda271 144 delete output;
145}
146
ad87ae62 147//-----------------------------------------------------------------------
ad87ae62 148void AliFlowAnalysisWithScalarProduct::WriteHistograms(TDirectoryFile *outputFileName)
149{
150 //store the final results in output .root file
151 fHistList->SetName("cobjSP");
152 fHistList->SetOwner(kTRUE);
153 outputFileName->Add(fHistList);
154 outputFileName->Write(outputFileName->GetName(), TObject::kSingleKey);
155}
156
8d312f00 157//-----------------------------------------------------------------------
158void AliFlowAnalysisWithScalarProduct::Init() {
159
160 //Define all histograms
d7eb18ec 161 cout<<"---Analysis with the Scalar Product Method--- Init"<<endl;
8d312f00 162
489d5531 163 //save old value and prevent histograms from being added to directory
164 //to avoid name clashes in case multiple analaysis objects are used
165 //in an analysis
a93de0f0 166
489d5531 167 Bool_t oldHistAddStatus = TH1::AddDirectoryStatus();
168 TH1::AddDirectory(kFALSE);
169
bee2efdc 170 Int_t iNbinsPt = AliFlowCommonConstants::GetMaster()->GetNbinsPt();
171 Double_t dPtMin = AliFlowCommonConstants::GetMaster()->GetPtMin();
172 Double_t dPtMax = AliFlowCommonConstants::GetMaster()->GetPtMax();
173 Int_t iNbinsEta = AliFlowCommonConstants::GetMaster()->GetNbinsEta();
174 Double_t dEtaMin = AliFlowCommonConstants::GetMaster()->GetEtaMin();
175 Double_t dEtaMax = AliFlowCommonConstants::GetMaster()->GetEtaMax();
8d312f00 176
1b302085 177 fHistProFlags = new TProfile("FlowPro_Flags_SP","Flow_Flags_SP",1,0,1,"s");
29ecccee 178 fHistProFlags->GetXaxis()->SetBinLabel(1,"fApplyCorrectionForNUA");
179 fHistList->Add(fHistProFlags);
180
1b302085 181 fHistProUQetaRP = new TProfile("FlowPro_UQetaRP_SP","Flow_UQetaRP_SP",iNbinsEta,dEtaMin,dEtaMax,"s");
395fadba 182 fHistProUQetaRP->SetXTitle("{eta}");
183 fHistProUQetaRP->SetYTitle("<uQ>");
184 fHistList->Add(fHistProUQetaRP);
185
1b302085 186 fHistProUQetaPOI = new TProfile("FlowPro_UQetaPOI_SP","Flow_UQetaPOI_SP",iNbinsEta,dEtaMin,dEtaMax,"s");
395fadba 187 fHistProUQetaPOI->SetXTitle("{eta}");
188 fHistProUQetaPOI->SetYTitle("<uQ>");
189 fHistList->Add(fHistProUQetaPOI);
190
1b302085 191 fHistProUQetaAllEventsPOI = new TProfile("FlowPro_UQetaAllEventsPOI_SP","FlowPro_UQetaAllEventsPOI_SP",iNbinsEta,dEtaMin,dEtaMax);
192 fHistProUQetaAllEventsPOI->SetXTitle("{eta}");
193 fHistProUQetaAllEventsPOI->SetYTitle("<uQ>");
194 fHistList->Add(fHistProUQetaAllEventsPOI);
195
196 fHistProUQPtRP = new TProfile("FlowPro_UQPtRP_SP","Flow_UQPtRP_SP",iNbinsPt,dPtMin,dPtMax,"s");
395fadba 197 fHistProUQPtRP->SetXTitle("p_t (GeV)");
198 fHistProUQPtRP->SetYTitle("<uQ>");
199 fHistList->Add(fHistProUQPtRP);
200
1b302085 201 fHistProUQPtPOI = new TProfile("FlowPro_UQPtPOI_SP","Flow_UQPtPOI_SP",iNbinsPt,dPtMin,dPtMax,"s");
395fadba 202 fHistProUQPtPOI->SetXTitle("p_t (GeV)");
203 fHistProUQPtPOI->SetYTitle("<uQ>");
204 fHistList->Add(fHistProUQPtPOI);
205
1b302085 206 fHistProUQPtAllEventsPOI = new TProfile("FlowPro_UQPtAllEventsPOI_SP","FlowPro_UQPtAllEventsPOI_SP",iNbinsPt,dPtMin,dPtMax);
207 fHistProUQPtAllEventsPOI->SetXTitle("p_t (GeV)");
208 fHistProUQPtAllEventsPOI->SetYTitle("<uQ>");
209 fHistList->Add(fHistProUQPtAllEventsPOI);
210
211 fHistProQNorm = new TProfile("FlowPro_QNorm_SP","FlowPro_QNorm_SP", 1, 0.5, 1.5,"s");
212 fHistProQNorm ->SetYTitle("<|Qa+Qb|>");
213 fHistList->Add(fHistProQNorm);
214
215 fHistProQaQb = new TProfile("FlowPro_QaQb_SP","FlowPro_QaQb_SP", 1, 0.5, 1.5,"s");
e4cecfa0 216 fHistProQaQb->SetYTitle("<QaQb>");
29ecccee 217 fHistList->Add(fHistProQaQb);
8d312f00 218
a554b203 219 fHistProQaQbNorm = new TProfile("FlowPro_QaQbNorm_SP","FlowPro_QaQbNorm_SP", 1, 0.5, 1.5,"s");
220 fHistProQaQbNorm->SetYTitle("<QaQb/MaMb>");
221 fHistList->Add(fHistProQaQbNorm);
29ecccee 222
223 fHistProQaQbReImNorm = new TProfile("FlowPro_QaQbReImNorm_SP","FlowPro_QaQbReImNorm_SP", 4, 0.5, 4.5,"s");
224 fHistProQaQbReImNorm->GetXaxis()->SetBinLabel(1,"#LT#LTsin(#phi_{a})#GT#GT");
225 fHistProQaQbReImNorm->GetXaxis()->SetBinLabel(2,"#LT#LTcos(#phi_{a})#GT#GT");
226 fHistProQaQbReImNorm->GetXaxis()->SetBinLabel(3,"#LT#LTsin(#phi_{b})#GT#GT");
227 fHistProQaQbReImNorm->GetXaxis()->SetBinLabel(4,"#LT#LTcos(#phi_{b})#GT#GT");
228 fHistList->Add(fHistProQaQbReImNorm);
229
230 fHistProNonIsotropicTermsQ = new TProfile("FlowPro_NonIsotropicTermsQ_SP","FlowPro_NonIsotropicTermsQ_SP", 2, 0.5, 2.5,"s");
231 fHistProNonIsotropicTermsQ->GetXaxis()->SetBinLabel(1,"#LT#LTsin(#phi_{a+b})#GT#GT");
232 fHistProNonIsotropicTermsQ->GetXaxis()->SetBinLabel(2,"#LT#LTcos(#phi_{a+b})#GT#GT");
233 fHistList->Add(fHistProNonIsotropicTermsQ);
234
235 TString rpPoi[2] = {"RP","POI"};
236 TString ptEta[2] = {"Pt","Eta"};
237 TString sinCos[2] = {"sin","cos"};
238 Int_t nBinsPtEta[2] = {iNbinsPt,iNbinsEta};
239 Double_t minPtEta[2] = {dPtMin,dEtaMin};
240 Double_t maxPtEta[2] = {dPtMax,dEtaMax};
241 for(Int_t rp=0;rp<2;rp++)
242 {
243 for(Int_t pe=0;pe<2;pe++)
244 {
245 for(Int_t sc=0;sc<2;sc++)
246 {
1b302085 247 fHistProNonIsotropicTermsU[rp][pe][sc] = new TProfile(Form("FlowPro_NonIsotropicTerms_%s_%s_%s_SP",rpPoi[rp].Data(),ptEta[pe].Data(),sinCos[sc].Data()),Form("FlowPro_NonIsotropicTerms_%s_%s_%s_SP",rpPoi[rp].Data(),ptEta[pe].Data(),sinCos[sc].Data()),nBinsPtEta[pe],minPtEta[pe],maxPtEta[pe]);
29ecccee 248 fHistList->Add(fHistProNonIsotropicTermsU[rp][pe][sc]);
249 }
250 }
251 }
252
12dc476e 253 fHistSumOfLinearWeights = new TH1D("Flow_SumOfLinearWeights_SP","Flow_SumOfLinearWeights_SP",1,-0.5, 0.5);
254 fHistSumOfLinearWeights -> SetYTitle("sum (*)");
255 fHistSumOfLinearWeights -> SetXTitle("sum (Ma*Mb)");
256 fHistList->Add(fHistSumOfLinearWeights);
257
258 fHistSumOfQuadraticWeights = new TH1D("Flow_SumOfQuadraticWeights_SP","Flow_SumOfQuadraticWeights_SP",1,-0.5, 0.5);
259 fHistSumOfQuadraticWeights -> SetYTitle("sum (*)");
260 fHistSumOfQuadraticWeights -> SetXTitle("sum (Ma*Mb)^2");
261 fHistList->Add(fHistSumOfQuadraticWeights);
262
1b302085 263 fHistProUQQaQbPtRP = new TProfile("FlowPro_UQQaQbPtRP_SP","FlowPro_UQQaQbPtRP_SP",iNbinsPt,dPtMin,dPtMax);
12dc476e 264 fHistProUQQaQbPtRP -> SetYTitle("<*>");
265 fHistProUQQaQbPtRP -> SetXTitle("<Qu QaQb>");
266 fHistList->Add(fHistProUQQaQbPtRP);
267
1b302085 268 fHistProUQQaQbEtaRP = new TProfile("FlowPro_UQQaQbEtaRP_SP","FlowPro_UQQaQbEtaRP_SP",iNbinsEta,dEtaMin,dEtaMax);
12dc476e 269 fHistProUQQaQbEtaRP -> SetYTitle("<*>");
270 fHistProUQQaQbEtaRP -> SetXTitle("<Qu QaQb>");
271 fHistList->Add(fHistProUQQaQbEtaRP);
272
1b302085 273 fHistProUQQaQbPtPOI = new TProfile("FlowPro_UQQaQbPtPOI_SP","FlowPro_UQQaQbPtPOI_SP",iNbinsPt,dPtMin,dPtMax);
12dc476e 274 fHistProUQQaQbPtPOI -> SetYTitle("<*>");
275 fHistProUQQaQbPtPOI -> SetXTitle("<Qu QaQb>");
276 fHistList->Add(fHistProUQQaQbPtPOI);
277
1b302085 278 fHistProUQQaQbEtaPOI = new TProfile("FlowPro_UQQaQbEtaPOI_SP","FlowPro_UQQaQbEtaPOI_SP",iNbinsEta,dEtaMin,dEtaMax);
12dc476e 279 fHistProUQQaQbEtaPOI -> SetYTitle("<*>");
280 fHistProUQQaQbEtaPOI -> SetXTitle("<Qu QaQb>");
281 fHistList->Add(fHistProUQQaQbEtaPOI);
282
283 TString weightFlag[3] = {"w_Qu_","w_Qu^2_","w_QuQaQb_"};
284 for(Int_t i=0;i<3;i++)
285 {
cff1bec4 286 fHistSumOfWeightsPtRP[i] = new TH1D(Form("Flow_SumOfWeights%sPtRP_SP",weightFlag[i].Data()),
287 Form("Flow_SumOfWeights%sPtRP_SP",weightFlag[i].Data()),iNbinsPt,dPtMin,dPtMax);
12dc476e 288 fHistSumOfWeightsPtRP[i] -> SetYTitle("sum (*)");
289 fHistSumOfWeightsPtRP[i] -> SetXTitle("p_{T}");
290 fHistList->Add(fHistSumOfWeightsPtRP[i]);
291
cff1bec4 292 fHistSumOfWeightsEtaRP[i] = new TH1D(Form("Flow_SumOfWeights%sEtaRP_SP",weightFlag[i].Data()),
293 Form("Flow_SumOfWeights%sEtaRP_SP",weightFlag[i].Data()),iNbinsEta,dEtaMin,dEtaMax);
12dc476e 294 fHistSumOfWeightsEtaRP[i] -> SetYTitle("sum (*)");
295 fHistSumOfWeightsEtaRP[i] -> SetXTitle("#eta");
296 fHistList->Add(fHistSumOfWeightsEtaRP[i]);
297
cff1bec4 298 fHistSumOfWeightsPtPOI[i] = new TH1D(Form("Flow_SumOfWeights%sPtPOI_SP",weightFlag[i].Data()),
299 Form("Flow_SumOfWeights%sPtPOI_SP",weightFlag[i].Data()),iNbinsPt,dPtMin,dPtMax);
12dc476e 300 fHistSumOfWeightsPtPOI[i] -> SetYTitle("sum (*)");
301 fHistSumOfWeightsPtPOI[i] -> SetXTitle("p_{T}");
302 fHistList->Add(fHistSumOfWeightsPtPOI[i]);
303
cff1bec4 304 fHistSumOfWeightsEtaPOI[i] = new TH1D(Form("Flow_SumOfWeights%sEtaPOI_SP",weightFlag[i].Data()),
305 Form("Flow_SumOfWeights%sEtaPOI_SP",weightFlag[i].Data()),iNbinsEta,dEtaMin,dEtaMax);
12dc476e 306 fHistSumOfWeightsEtaPOI[i] -> SetYTitle("sum (*)");
307 fHistSumOfWeightsEtaPOI[i] -> SetXTitle("#eta");
308 fHistList->Add(fHistSumOfWeightsEtaPOI[i]);
309 }
29ecccee 310
1b302085 311 fCommonHistsSP = new AliFlowCommonHist("AliFlowCommonHistSP");
312 fHistList->Add(fCommonHistsSP);
313 fCommonHistsResSP = new AliFlowCommonHistResults("AliFlowCommonHistResultsSP");
314 fHistList->Add(fCommonHistsResSP);
315 fCommonHistsmuQ = new AliFlowCommonHist("AliFlowCommonHistmuQ");
316 fHistList->Add(fCommonHistsmuQ);
317
50fe33aa 318 (fCommonHistsSP->GetHarmonic())->Fill(0.5,fHarmonic); // store harmonic
319 (fCommonHistsmuQ->GetHarmonic())->Fill(0.5,fHarmonic); // store harmonic
320
1b302085 321 fHistQNorm = new TH1D("Flow_QNorm_SP","Flow_QNorm_SP",110,0.,1.1);
322 fHistQNorm -> SetYTitle("dN/d(|(Qa+Qb)/(Ma+Mb)|)");
323 fHistQNorm -> SetXTitle("|(Qa+Qb)/(Ma+Mb)|");
324 fHistList->Add(fHistQNorm);
395fadba 325
a554b203 326 fHistQaQb = new TH1D("Flow_QaQb_SP","Flow_QaQb_SP",200,-100.,100.);
327 fHistQaQb -> SetYTitle("dN/dQaQb");
328 fHistQaQb -> SetXTitle("QaQb");
329 fHistList->Add(fHistQaQb);
330
331 fHistQaQbNorm = new TH1D("Flow_QaQbNorm_SP","Flow_QaQbNorm_SP",44,-1.1,1.1);
332 fHistQaQbNorm -> SetYTitle("dN/d(QaQb/MaMb)");
333 fHistQaQbNorm -> SetXTitle("QaQb/MaMb");
334 fHistList->Add(fHistQaQbNorm);
335
1b302085 336 fHistQNormvsQaQbNorm = new TH2D("Flow_QNormvsQaQbNorm_SP","Flow_QNormvsQaQbNorm_SP",88,-1.1,1.1,22,0.,1.1);
337 fHistQNormvsQaQbNorm -> SetYTitle("|Q/Mq|");
338 fHistQNormvsQaQbNorm -> SetXTitle("QaQb/MaMb");
339 fHistList->Add(fHistQNormvsQaQbNorm);
340
a554b203 341 fHistQaQbCos = new TH1D("Flow_QaQbCos_SP","Flow_QaQbCos_SP",63,0.,TMath::Pi());
342 fHistQaQbCos -> SetYTitle("dN/d(#phi)");
343 fHistQaQbCos -> SetXTitle("#phi");
344 fHistList->Add(fHistQaQbCos);
345
346 fHistResolution = new TH1D("Flow_resolution_SP","Flow_resolution_SP",100,-1.0,1.0);
347 fHistResolution -> SetYTitle("dN/d(cos(2(#phi_a - #phi_b))");
348 fHistResolution -> SetXTitle("cos(2*(#phi_a - #phi_b))");
349 fHistList->Add(fHistResolution);
350
351 fHistQaNorm = new TH1D("Flow_QaNorm_SP","Flow_QaNorm_SP",22,0.,1.1);
352 fHistQaNorm -> SetYTitle("dN/d(|Qa/Ma|)");
353 fHistQaNorm -> SetXTitle("|Qa/Ma|");
354 fHistList->Add(fHistQaNorm);
355
356 fHistQaNormvsMa = new TH2D("Flow_QaNormvsMa_SP","Flow_QaNormvsMa_SP",100,0.,100.,22,0.,1.1);
357 fHistQaNormvsMa -> SetYTitle("|Qa/Ma|");
358 fHistQaNormvsMa -> SetXTitle("Ma");
359 fHistList->Add(fHistQaNormvsMa);
360
361 fHistQbNorm = new TH1D("Flow_QbNorm_SP","Flow_QbNorm_SP",22,0.,1.1);
362 fHistQbNorm -> SetYTitle("dN/d(|Qb/Mb|)");
363 fHistQbNorm -> SetXTitle("|Qb/Mb|");
364 fHistList->Add(fHistQbNorm);
365
366 fHistQbNormvsMb = new TH2D("Flow_QbNormvsMb_SP","Flow_QbNormvsMb_SP",100,0.,100.,22,0.,1.1);
367 fHistQbNormvsMb -> SetYTitle("|Qb/Mb|");
368 fHistQbNormvsMb -> SetXTitle("|Mb|");
369 fHistList->Add(fHistQbNormvsMb);
370
371 fHistMavsMb = new TH2D("Flow_MavsMb_SP","Flow_MavsMb_SP",100,0.,100.,100,0.,100.);
372 fHistMavsMb -> SetYTitle("Ma");
373 fHistMavsMb -> SetXTitle("Mb");
374 fHistList->Add(fHistMavsMb);
375
376
29195b69 377 //weights
378 if(fUsePhiWeights) {
379 if(!fWeightsList) {
380 cout<<"WARNING: fWeightsList is NULL in the Scalar Product method."<<endl;
381 exit(0);
382 }
04c6b875 383 if(fWeightsList->FindObject("phi_weights_sub0")) {
384 fPhiWeightsSub0 = dynamic_cast<TH1F*>
385 (fWeightsList->FindObject("phi_weights_sub0"));
386 fHistList->Add(fPhiWeightsSub0);
29195b69 387 } else {
388 cout<<"WARNING: histogram with phi weights is not accessible in Scalar Product"<<endl;
389 exit(0);
390 }
04c6b875 391 if(fWeightsList->FindObject("phi_weights_sub1")) {
392 fPhiWeightsSub1 = dynamic_cast<TH1F*>
393 (fWeightsList->FindObject("phi_weights_sub1"));
394 fHistList->Add(fPhiWeightsSub1);
395 } else {
396 cout<<"WARNING: histogram with phi weights is not accessible in Scalar Product"<<endl;
397 exit(0);
398 }
399
29195b69 400 } // end of if(fUsePhiWeights)
401
29ecccee 402 fEventNumber = 0; //set number of events to zero
403
404 //store all boolean flags needed in Finish():
405 this->StoreFlags();
489d5531 406
407 TH1::AddDirectory(oldHistAddStatus);
e2d51347 408}
409
8d312f00 410//-----------------------------------------------------------------------
8232a5ec 411void AliFlowAnalysisWithScalarProduct::Make(AliFlowEventSimple* anEvent) {
8d312f00 412
1b302085 413
414 if (anEvent) {
415
416 //Calculate muQ (for comparing pp and PbPb)
417 FillmuQ(anEvent);
418
419 //Calculate flow based on <QaQb/MaMb> = <v^2>
420 FillSP(anEvent);
421
422 }
423}
424
425//-----------------------------------------------------------------------
426void AliFlowAnalysisWithScalarProduct::FillSP(AliFlowEventSimple* anEvent) {
427
489d5531 428 //Calculate flow based on <QaQb/MaMb> = <v^2>
12dc476e 429
430 //Fill histograms
8232a5ec 431 if (anEvent) {
8d312f00 432
cbbaf54a 433 //get Q vectors for the eta-subevents
b125a454 434 AliFlowVector* vQarray = new AliFlowVector[2];
29195b69 435 if (fUsePhiWeights) {
50fe33aa 436 anEvent->Get2Qsub(vQarray,fHarmonic,fWeightsList,kTRUE);
29195b69 437 } else {
50fe33aa 438 anEvent->Get2Qsub(vQarray,fHarmonic);
29195b69 439 }
12dc476e 440 //subevent a
b125a454 441 AliFlowVector vQa = vQarray[0];
12dc476e 442 //subevent b
b125a454 443 AliFlowVector vQb = vQarray[1];
12dc476e 444
1b302085 445 //For calculating v2 only events should be taken where both subevents are not empty
7a01f4a7 446 //check that the subevents are not empty:
12dc476e 447 Double_t dMa = vQa.GetMult();
12dc476e 448 Double_t dMb = vQb.GetMult();
a93de0f0 449 if (dMa > 0. && dMb > 0.) {
7a01f4a7 450
a554b203 451 //request that the subevent multiplicities are not too different
1b302085 452 //fRelDiffMsub can be set from the configuration macro
a554b203 453 Double_t dRelDiff = TMath::Abs((dMa - dMb)/(dMa + dMb));
454 if (dRelDiff < fRelDiffMsub) {
455
3c5d5752 456 //fill control histograms
1b302085 457 if (fUsePhiWeights) {
458 fCommonHistsSP->FillControlHistograms(anEvent,fWeightsList,kTRUE);
459 } else {
460 fCommonHistsSP->FillControlHistograms(anEvent);
461 }
a554b203 462
463 //fill some SP control histograms
464 fHistProQaQb -> Fill(1.,vQa*vQb,1.); //Fill with weight 1 -> Weight with MaMb????
04c6b875 465 fHistQaQbCos ->Fill(TMath::ACos((vQa/vQa.Mod())*(vQb/vQb.Mod()))); //Acos(Qa*Qb) = angle
a554b203 466 fHistResolution -> Fill(TMath::Cos( vQa.Phi()- vQb.Phi() )); //vQa.Phi() returns 2*phi
467 fHistQaQb -> Fill(vQa*vQb);
468 fHistMavsMb -> Fill(dMb,dMa);
469
470 //get total Q vector = the sum of subevent a and subevent b
471 AliFlowVector vQ = vQa + vQb;
a93de0f0 472
29ecccee 473 //needed to correct for non-uniform acceptance:
1b302085 474 fHistProNonIsotropicTermsQ->Fill(1.,vQ.Y()/(dMa+dMb),dMa+dMb);
475 fHistProNonIsotropicTermsQ->Fill(2.,vQ.X()/(dMa+dMb),dMa+dMb);
a93de0f0 476
a554b203 477 //weight the Q vectors for the subevents by the multiplicity
a93de0f0 478 //Note: Weight Q only in the particle loop when it is clear
479 //if it should be (m-1) or M
a554b203 480 Double_t dQXa = vQa.X()/dMa;
481 Double_t dQYa = vQa.Y()/dMa;
482 vQa.Set(dQXa,dQYa);
483
484 Double_t dQXb = vQb.X()/dMb;
485 Double_t dQYb = vQb.Y()/dMb;
486 vQb.Set(dQXb,dQYb);
12dc476e 487
a554b203 488 //scalar product of the two subevents
489 Double_t dQaQb = (vQa*vQb);
490 fHistProQaQbNorm -> Fill(1.,dQaQb,dMa*dMb); //Fill (QaQb/MaMb) with weight (MaMb).
491 //needed for the error calculation:
492 fHistSumOfLinearWeights -> Fill(0.,dMa*dMb);
493 fHistSumOfQuadraticWeights -> Fill(0.,pow(dMa*dMb,2.));
a93de0f0 494 //needed for correcting non-uniform acceptance:
495 fHistProQaQbReImNorm->Fill(1.,dQYa,dMa); // to get <<sin(phi_a)>>
496 fHistProQaQbReImNorm->Fill(2.,dQXa,dMa); // to get <<cos(phi_a)>>
497 fHistProQaQbReImNorm->Fill(3.,dQYb,dMb); // to get <<sin(phi_b)>>
498 fHistProQaQbReImNorm->Fill(4.,dQXb,dMb); // to get <<cos(phi_b)>>
499
a554b203 500 //fill some SP control histograms
501 fHistQaQbNorm ->Fill(vQa*vQb);
502 fHistQaNorm ->Fill(vQa.Mod());
503 fHistQaNormvsMa->Fill(dMa,vQa.Mod());
504 fHistQbNorm ->Fill(vQb.Mod());
505 fHistQbNormvsMb->Fill(dMb,vQb.Mod());
506
507 //loop over the tracks of the event
508 AliFlowTrackSimple* pTrack = NULL;
509 Int_t iNumberOfTracks = anEvent->NumberOfTracks();
510 Double_t dMq = vQ.GetMult();
a93de0f0 511
a554b203 512 for (Int_t i=0;i<iNumberOfTracks;i++)
513 {
514 pTrack = anEvent->GetTrack(i) ;
515 if (pTrack){
516 Double_t dPhi = pTrack->Phi();
1b302085 517
a554b203 518 //calculate vU
519 TVector2 vU;
1b302085 520 //do not need to use weight for v as the length will be made 1
50fe33aa 521 Double_t dUX = TMath::Cos(fHarmonic*dPhi);
522 Double_t dUY = TMath::Sin(fHarmonic*dPhi);
a554b203 523 vU.Set(dUX,dUY);
524 Double_t dModulus = vU.Mod();
a93de0f0 525 if (dModulus > 0.) vU.Set(dUX/dModulus,dUY/dModulus); // make length 1
a554b203 526 else cerr<<"dModulus is zero!"<<endl;
12dc476e 527
a554b203 528 //redefine the Q vector and devide by its multiplicity
529 TVector2 vQm;
530 Double_t dQmX = 0.;
531 Double_t dQmY = 0.;
532 //subtract particle from the flowvector if used to define it
533 if (pTrack->InSubevent(0) || pTrack->InSubevent(1)) {
1b302085 534 //set default phi weight to 1
535 Double_t dW = 1.;
536 //if phi weights are used
537 if(fUsePhiWeights && fPhiWeightsSub0 && fPhiWeightsSub1) {
538 //value of the center of the phi bin
539 Double_t dPhiCenter = 0.;
540 if (pTrack->InSubevent(0) ) {
541 Int_t iNBinsPhiSub0 = fPhiWeightsSub0->GetNbinsX();
542 Int_t phiBin = 1+(Int_t)(TMath::Floor(dPhi*iNBinsPhiSub0/TMath::TwoPi()));
543 dW = fPhiWeightsSub0->GetBinContent(phiBin);
544 dPhiCenter = fPhiWeightsSub0->GetBinCenter(phiBin);
50fe33aa 545 dQmX = (vQ.X() - dW*(pTrack->Weight())* TMath::Cos(fHarmonic*dPhiCenter) )/(dMq-1);
546 dQmY = (vQ.Y() - dW*(pTrack->Weight())* TMath::Sin(fHarmonic*dPhiCenter) )/(dMq-1);
1b302085 547
548 vQm.Set(dQmX,dQmY);
549 }
550
551 else if ( pTrack->InSubevent(1)) {
552 Int_t iNBinsPhiSub1 = fPhiWeightsSub1->GetNbinsX();
553 Int_t phiBin = 1+(Int_t)(TMath::Floor(dPhi*iNBinsPhiSub1/TMath::TwoPi()));
554 dW = fPhiWeightsSub1->GetBinContent(phiBin);
555 dPhiCenter = fPhiWeightsSub1->GetBinCenter(phiBin);
50fe33aa 556 dQmX = (vQ.X() - dW*(pTrack->Weight())* TMath::Cos(fHarmonic*dPhiCenter) )/(dMq-1);
557 dQmY = (vQ.Y() - dW*(pTrack->Weight())* TMath::Sin(fHarmonic*dPhiCenter) )/(dMq-1);
1b302085 558
559 vQm.Set(dQmX,dQmY);
560 }
561 //bin = 1 + value*nbins/range
562 //TMath::Floor rounds to the lower integer
563 }
564 // if no phi weights are used
565 else {
566 dQmX = (vQ.X() - (pTrack->Weight())*dUX)/(dMq-1);
567 dQmY = (vQ.Y() - (pTrack->Weight())*dUY)/(dMq-1);
568 vQm.Set(dQmX,dQmY);
569 }
570
a554b203 571 //dUQ = scalar product of vU and vQm
572 Double_t dUQ = (vU * vQm);
573 Double_t dPt = pTrack->Pt();
574 Double_t dEta = pTrack->Eta();
a93de0f0 575
a554b203 576 //fill the profile histograms
577 if (pTrack->InRPSelection()) {
578 fHistProUQetaRP -> Fill(dEta,dUQ,(dMq-1)); //Fill (Qu/(Mq-1)) with weight (Mq-1)
579 //needed for the error calculation:
580 fHistProUQQaQbEtaRP -> Fill(dEta,dUQ*dQaQb,(dMq-1)*dMa*dMb); //Fill [Qu/(Mq-1)]*[QaQb/MaMb] with weight (Mq-1)MaMb
581 fHistProUQPtRP -> Fill(dPt,dUQ,(dMq-1)); //Fill (Qu/(Mq-1)) with weight (Mq-1)
582 fHistProUQQaQbPtRP -> Fill(dPt,dUQ*dQaQb,(dMq-1)*dMa*dMb); //Fill [Qu/(Mq-1)]*[QaQb/MaMb] with weight (Mq-1)MaMb
583
584 fHistSumOfWeightsEtaRP[0]->Fill(dEta,(dMq-1)); // sum of Mq-1
585 fHistSumOfWeightsEtaRP[1]->Fill(dEta,pow((dMq-1),2.));// sum of (Mq-1)^2
586 fHistSumOfWeightsEtaRP[2]->Fill(dEta,(dMq-1)*dMa*dMb);// sum of (Mq-1)*MaMb
587 fHistSumOfWeightsPtRP[0]->Fill(dPt,(dMq-1)); // sum of Mq-1
588 fHistSumOfWeightsPtRP[1]->Fill(dPt,pow((dMq-1),2.)); // sum of (Mq-1)^2
29ecccee 589 fHistSumOfWeightsPtRP[2]->Fill(dPt,(dMq-1)*dMa*dMb); // sum of (Mq-1)*MaMb
590 //nonisotropic terms:
591 fHistProNonIsotropicTermsU[0][0][0]->Fill(dPt,dUY,1.);
592 fHistProNonIsotropicTermsU[0][0][1]->Fill(dPt,dUX,1.);
593 fHistProNonIsotropicTermsU[0][1][0]->Fill(dEta,dUY,1.);
594 fHistProNonIsotropicTermsU[0][1][1]->Fill(dEta,dUX,1.);
a554b203 595 }
596 if (pTrack->InPOISelection()) {
597 fHistProUQetaPOI -> Fill(dEta,dUQ,(dMq-1));//Fill (Qu/(Mq-1)) with weight (Mq-1)
598 //needed for the error calculation:
599 fHistProUQQaQbEtaPOI -> Fill(dEta,dUQ*dQaQb,(dMq-1)*dMa*dMb); //Fill [Qu/(Mq-1)]*[QaQb/MaMb] with weight (Mq-1)MaMb
600 fHistProUQPtPOI -> Fill(dPt,dUQ,(dMq-1)); //Fill (Qu/(Mq-1)) with weight (Mq-1)
601 fHistProUQQaQbPtPOI -> Fill(dPt,dUQ*dQaQb,(dMq-1)*dMa*dMb); //Fill [Qu/(Mq-1)]*[QaQb/MaMb] with weight (Mq-1)MaMb
602
603 fHistSumOfWeightsEtaPOI[0]->Fill(dEta,(dMq-1)); // sum of Mq-1
604 fHistSumOfWeightsEtaPOI[1]->Fill(dEta,pow((dMq-1),2.));// sum of (Mq-1)^2
605 fHistSumOfWeightsEtaPOI[2]->Fill(dEta,(dMq-1)*dMa*dMb);// sum of (Mq-1)*MaMb
606 fHistSumOfWeightsPtPOI[0]->Fill(dPt,(dMq-1)); // sum of Mq-1
607 fHistSumOfWeightsPtPOI[1]->Fill(dPt,pow((dMq-1),2.)); // sum of (Mq-1)^2
29ecccee 608 fHistSumOfWeightsPtPOI[2]->Fill(dPt,(dMq-1)*dMa*dMb); // sum of (Mq-1)*MaMb
609 //nonisotropic terms:
610 fHistProNonIsotropicTermsU[1][0][0]->Fill(dPt,dUY,1.);
611 fHistProNonIsotropicTermsU[1][0][1]->Fill(dPt,dUX,1.);
612 fHistProNonIsotropicTermsU[1][1][0]->Fill(dEta,dUY,1.);
613 fHistProNonIsotropicTermsU[1][1][1]->Fill(dEta,dUX,1.);
a554b203 614 }
7a01f4a7 615
a554b203 616 } else { //do not subtract the particle from the flowvector
617 dQmX = vQ.X()/dMq;
618 dQmY = vQ.Y()/dMq;
619 vQm.Set(dQmX,dQmY);
1b302085 620
621 //fill histograms with vQm
622 fHistProQNorm->Fill(1.,vQm.Mod(),dMq);
623 fHistQNorm->Fill(vQm.Mod());
624 fHistQNormvsQaQbNorm->Fill(vQa*vQb ,vQm.Mod());
7a01f4a7 625
a554b203 626 //dUQ = scalar product of vU and vQm
627 Double_t dUQ = (vU * vQm);
628 Double_t dPt = pTrack->Pt();
629 Double_t dEta = pTrack->Eta();
a93de0f0 630
a554b203 631 //fill the profile histograms
632 if (pTrack->InRPSelection()) {
633 fHistProUQetaRP -> Fill(dEta,dUQ,dMq); //Fill (Qu/Mq) with weight Mq
634 //needed for the error calculation:
635 fHistProUQQaQbEtaRP -> Fill(dEta,dUQ*dQaQb,dMq*dMa*dMb); //Fill [Qu/Mq]*[QaQb/MaMb] with weight Mq*MaMb
636 fHistProUQPtRP -> Fill(dPt,dUQ,dMq); //Fill (Qu/Mq) with weight Mq
637 fHistProUQQaQbPtRP -> Fill(dPt,dUQ*dQaQb,dMq*dMa*dMb); //Fill [Qu/Mq]*[QaQb/MaMb] with weight Mq*MaMb
638
639 fHistSumOfWeightsEtaRP[0]->Fill(dEta,dMq); // sum of Mq
640 fHistSumOfWeightsEtaRP[1]->Fill(dEta,pow(dMq,2.));// sum of Mq^2
641 fHistSumOfWeightsEtaRP[2]->Fill(dEta,dMq*dMa*dMb);// sum of Mq*MaMb
642 fHistSumOfWeightsPtRP[0]->Fill(dPt,dMq); // sum of Mq
643 fHistSumOfWeightsPtRP[1]->Fill(dPt,pow(dMq,2.)); // sum of Mq^2
29ecccee 644 fHistSumOfWeightsPtRP[2]->Fill(dPt,dMq*dMa*dMb); // sum of Mq*MaMb
645 //nonisotropic terms:
646 fHistProNonIsotropicTermsU[0][0][0]->Fill(dPt,dUY,1.);
647 fHistProNonIsotropicTermsU[0][0][1]->Fill(dPt,dUX,1.);
648 fHistProNonIsotropicTermsU[0][1][0]->Fill(dEta,dUY,1.);
649 fHistProNonIsotropicTermsU[0][1][1]->Fill(dEta,dUX,1.);
a554b203 650 }
651 if (pTrack->InPOISelection()) {
652 fHistProUQetaPOI -> Fill(dEta,dUQ,dMq); //Fill (Qu/Mq) with weight Mq
653 //needed for the error calculation:
654 fHistProUQQaQbEtaPOI -> Fill(dEta,dUQ*dQaQb,dMq*dMa*dMb); //Fill [Qu/Mq]*[QaQb/MaMb] with weight Mq*MaMb
655 fHistProUQPtPOI -> Fill(dPt,dUQ,dMq); //Fill (Qu/Mq) with weight Mq
656 fHistProUQQaQbPtPOI -> Fill(dPt,dUQ*dQaQb,dMq*dMa*dMb); //Fill [Qu/Mq]*[QaQb/MaMb] with weight Mq*MaMb
657
658 fHistSumOfWeightsEtaPOI[0]->Fill(dEta,dMq); // sum of Mq
659 fHistSumOfWeightsEtaPOI[1]->Fill(dEta,pow(dMq,2.));// sum of Mq^2
660 fHistSumOfWeightsEtaPOI[2]->Fill(dEta,dMq*dMa*dMb);// sum of Mq*MaMb
661 fHistSumOfWeightsPtPOI[0]->Fill(dPt,dMq); // sum of Mq
662 fHistSumOfWeightsPtPOI[1]->Fill(dPt,pow(dMq,2.)); // sum of Mq^2
663 fHistSumOfWeightsPtPOI[2]->Fill(dPt,dMq*dMa*dMb); // sum of Mq*MaMb
29ecccee 664 //nonisotropic terms:
665 fHistProNonIsotropicTermsU[1][0][0]->Fill(dPt,dUY,1.);
666 fHistProNonIsotropicTermsU[1][0][1]->Fill(dPt,dUX,1.);
667 fHistProNonIsotropicTermsU[1][1][0]->Fill(dEta,dUY,1.);
668 fHistProNonIsotropicTermsU[1][1][1]->Fill(dEta,dUX,1.);
a554b203 669 }
670 }//track not in subevents
7a01f4a7 671
a554b203 672 }//track
7a01f4a7 673
a554b203 674 }//loop over tracks
675
676 fEventNumber++;
677
678 } //difference Ma and Mb
12dc476e 679
7a01f4a7 680 }// subevents not empty
b125a454 681 delete [] vQarray;
a554b203 682
7a01f4a7 683 } //event
a554b203 684
1b302085 685}//end of FillSP()
686
687//-----------------------------------------------------------------------
688void AliFlowAnalysisWithScalarProduct::FillmuQ(AliFlowEventSimple* anEvent) {
689
690 if (anEvent) {
691
692 //get Q vectors for the eta-subevents
693 AliFlowVector* vQarray = new AliFlowVector[2];
694 if (fUsePhiWeights) {
50fe33aa 695 anEvent->Get2Qsub(vQarray,fHarmonic,fWeightsList,kTRUE);
1b302085 696 } else {
50fe33aa 697 anEvent->Get2Qsub(vQarray,fHarmonic);
1b302085 698 }
699 //subevent a
700 AliFlowVector vQa = vQarray[0];
701 //subevent b
702 AliFlowVector vQb = vQarray[1];
703
704 //get total Q vector = the sum of subevent a and subevent b
705 AliFlowVector vQ;
706 if (vQa.GetMult() > 0 && vQb.GetMult() > 0) {
707 vQ = vQa + vQb;
708 } else if ( vQa.GetMult() > 0 && !(vQb.GetMult() > 0) ){
709 vQ = vQa;
710 } else if ( !(vQa.GetMult() > 0) && vQb.GetMult() > 0 ){
711 vQ = vQb;
712 }
713
714 //For calculating uQ for comparison all events should be taken also if one of the subevents is empty
715 //check if the total Q vector is not empty
716 Double_t dMq = vQ.GetMult();
717 if (dMq > 0.) {
718
719 //Fill control histograms
720 if (fUsePhiWeights) {
721 fCommonHistsmuQ->FillControlHistograms(anEvent,fWeightsList,kTRUE);
722 } else {
723 fCommonHistsmuQ->FillControlHistograms(anEvent);
724 }
725
726 //loop over all POI tracks and fill uQ
727 AliFlowTrackSimple* pTrack = NULL;
728 for (Int_t i=0;i<anEvent->NumberOfTracks();i++) {
729 pTrack = anEvent->GetTrack(i) ;
730 if (pTrack){
731
732 if (pTrack->InPOISelection()) {
733
734 Double_t dPhi = pTrack->Phi();
735 //weights do not need to be used as the length of vU will be set to 1
736
737 //calculate vU
738 TVector2 vU;
50fe33aa 739 Double_t dUX = TMath::Cos(fHarmonic*dPhi);
740 Double_t dUY = TMath::Sin(fHarmonic*dPhi);
1b302085 741 vU.Set(dUX,dUY);
742 Double_t dModulus = vU.Mod();
743 // make length 1
744 if (dModulus!=0.) vU.Set(dUX/dModulus,dUY/dModulus);
745 else cerr<<"dModulus is zero!"<<endl;
746
747 //redefine the Q vector
748 TVector2 vQm;
749 Double_t dQmX = 0.;
750 Double_t dQmY = 0.;
751 //subtract particle from the flowvector if used to define it
752 if (pTrack->InSubevent(0) || pTrack->InSubevent(1)) {
753 //the number of tracks contributing to vQ must be more than 1
754 if (dMq > 1) {
755 //set default phi weight to 1
756 Double_t dW = 1.;
757 //if phi weights are used
758 if(fUsePhiWeights && fPhiWeightsSub0 && fPhiWeightsSub1) {
759 //value of the center of the phi bin
760 Double_t dPhiCenter = 0.;
761 if (pTrack->InSubevent(0) ) {
762 Int_t iNBinsPhiSub0 = fPhiWeightsSub0->GetNbinsX();
763 Int_t phiBin = 1+(Int_t)(TMath::Floor(dPhi*iNBinsPhiSub0/TMath::TwoPi()));
764 dW = fPhiWeightsSub0->GetBinContent(phiBin);
765 dPhiCenter = fPhiWeightsSub0->GetBinCenter(phiBin);
50fe33aa 766 dQmX = (vQ.X() - dW*(pTrack->Weight())* TMath::Cos(fHarmonic*dPhiCenter) );
767 dQmY = (vQ.Y() - dW*(pTrack->Weight())* TMath::Sin(fHarmonic*dPhiCenter) );
1b302085 768
769 vQm.Set(dQmX,dQmY);
770 }
771
772 else if ( pTrack->InSubevent(1)) {
773 Int_t iNBinsPhiSub1 = fPhiWeightsSub1->GetNbinsX();
774 Int_t phiBin = 1+(Int_t)(TMath::Floor(dPhi*iNBinsPhiSub1/TMath::TwoPi()));
775 dW = fPhiWeightsSub1->GetBinContent(phiBin);
776 dPhiCenter = fPhiWeightsSub1->GetBinCenter(phiBin);
50fe33aa 777 dQmX = (vQ.X() - dW*(pTrack->Weight())* TMath::Cos(fHarmonic*dPhiCenter) );
778 dQmY = (vQ.Y() - dW*(pTrack->Weight())* TMath::Sin(fHarmonic*dPhiCenter) );
1b302085 779
780 vQm.Set(dQmX,dQmY);
781 }
782 //bin = 1 + value*nbins/range
783 //TMath::Floor rounds to the lower integer
784 }
785 // if no phi weights are used
786 else {
787 dQmX = (vQ.X() - (pTrack->Weight())*dUX);
788 dQmY = (vQ.Y() - (pTrack->Weight())*dUY);
789 vQm.Set(dQmX,dQmY);
790 }
791
792 //dUQ = scalar product of vU and vQm
793 Double_t dUQ = (vU * vQm);
794 Double_t dPt = pTrack->Pt();
795 Double_t dEta = pTrack->Eta();
796 //fill the profile histograms
797 fHistProUQetaAllEventsPOI -> Fill(dEta,dUQ); //Fill (Qu)
798 fHistProUQPtAllEventsPOI -> Fill(dPt,dUQ); //Fill (Qu)
799
800 } //dMq > 1
801 }
802 else { //do not subtract the particle from the flowvector
803
804 dQmX = vQ.X();
805 dQmY = vQ.Y();
806 vQm.Set(dQmX,dQmY);
807
808 //dUQ = scalar product of vU and vQm
809 Double_t dUQ = (vU * vQm);
810 Double_t dPt = pTrack->Pt();
811 Double_t dEta = pTrack->Eta();
812 //fill the profile histograms
813 fHistProUQetaAllEventsPOI -> Fill(dEta,dUQ); //Fill (Qu)
814 fHistProUQPtAllEventsPOI -> Fill(dPt,dUQ); //Fill (Qu)
815
816 }
817
818 } //in POI selection
819 } //track valid
820 } //end of loop over tracks
821 } //Q vector is not empty
822
823 } //anEvent valid
824
825} //end of FillmuQ
8d312f00 826
12dc476e 827//--------------------------------------------------------------------
fd46c3dd 828void AliFlowAnalysisWithScalarProduct::GetOutputHistograms(TList *outputListHistos){
829
830 //get pointers to all output histograms (called before Finish())
12dc476e 831
fd46c3dd 832 if (outputListHistos) {
833 //Get the common histograms from the output list
1b302085 834 AliFlowCommonHist *pCommonHistSP = dynamic_cast<AliFlowCommonHist*>
fd46c3dd 835 (outputListHistos->FindObject("AliFlowCommonHistSP"));
1b302085 836 AliFlowCommonHistResults *pCommonHistResultsSP = dynamic_cast<AliFlowCommonHistResults*>
fd46c3dd 837 (outputListHistos->FindObject("AliFlowCommonHistResultsSP"));
1b302085 838 AliFlowCommonHist *pCommonHistmuQ = dynamic_cast<AliFlowCommonHist*>
839 (outputListHistos->FindObject("AliFlowCommonHistmuQ"));
840
841 TProfile* pHistProQNorm = dynamic_cast<TProfile*>(outputListHistos->FindObject("FlowPro_QNorm_SP"));
842 TProfile* pHistProQaQb = dynamic_cast<TProfile*>(outputListHistos->FindObject("FlowPro_QaQb_SP"));
a554b203 843 TProfile* pHistProQaQbNorm = dynamic_cast<TProfile*>(outputListHistos->FindObject("FlowPro_QaQbNorm_SP"));
29ecccee 844 TProfile* pHistProQaQbReImNorm = dynamic_cast<TProfile*>(outputListHistos->FindObject("FlowPro_QaQbReImNorm_SP"));
845 TProfile* pHistProNonIsotropicTermsQ = dynamic_cast<TProfile*>(outputListHistos->FindObject("FlowPro_NonIsotropicTermsQ_SP"));
a554b203 846 TH1D* pHistSumOfLinearWeights = dynamic_cast<TH1D*>(outputListHistos->FindObject("Flow_SumOfLinearWeights_SP"));
847 TH1D* pHistSumOfQuadraticWeights = dynamic_cast<TH1D*>(outputListHistos->FindObject("Flow_SumOfQuadraticWeights_SP"));
848
1b302085 849 TProfile* pHistProFlags = dynamic_cast<TProfile*>(outputListHistos->FindObject("FlowPro_Flags_SP"));
850 TProfile* pHistProUQetaRP = dynamic_cast<TProfile*>(outputListHistos->FindObject("FlowPro_UQetaRP_SP"));
851 TProfile* pHistProUQetaPOI = dynamic_cast<TProfile*>(outputListHistos->FindObject("FlowPro_UQetaPOI_SP"));
852 TProfile* pHistProUQPtRP = dynamic_cast<TProfile*>(outputListHistos->FindObject("FlowPro_UQPtRP_SP"));
853 TProfile* pHistProUQPtPOI = dynamic_cast<TProfile*>(outputListHistos->FindObject("FlowPro_UQPtPOI_SP"));
854 TProfile* pHistProUQQaQbPtRP = dynamic_cast<TProfile*>(outputListHistos->FindObject("FlowPro_UQQaQbPtRP_SP"));
855 TProfile* pHistProUQQaQbEtaRP = dynamic_cast<TProfile*>(outputListHistos->FindObject("FlowPro_UQQaQbEtaRP_SP"));
856 TProfile* pHistProUQQaQbPtPOI = dynamic_cast<TProfile*>(outputListHistos->FindObject("FlowPro_UQQaQbPtPOI_SP"));
857 TProfile* pHistProUQQaQbEtaPOI = dynamic_cast<TProfile*>(outputListHistos->FindObject("FlowPro_UQQaQbEtaPOI_SP"));
12dc476e 858 TString weightFlag[3] = {"w_Qu_","w_Qu^2_","w_QuQaQb_"};
1b302085 859
860
12dc476e 861 TH1D* pHistSumOfWeightsPtRP[3] = {NULL};
862 TH1D* pHistSumOfWeightsEtaRP[3] = {NULL};
863 TH1D* pHistSumOfWeightsPtPOI[3] = {NULL};
a93de0f0 864 TH1D* pHistSumOfWeightsEtaPOI[3] = {NULL};
865
866 for(Int_t i=0;i<3;i++) {
867 pHistSumOfWeightsPtRP[i] = dynamic_cast<TH1D*>(outputListHistos->FindObject(Form("Flow_SumOfWeights%sPtRP_SP",weightFlag[i].Data())));
868 pHistSumOfWeightsEtaRP[i] = dynamic_cast<TH1D*>(outputListHistos->FindObject(Form("Flow_SumOfWeights%sEtaRP_SP",weightFlag[i].Data())));
869 pHistSumOfWeightsPtPOI[i] = dynamic_cast<TH1D*>(outputListHistos->FindObject(Form("Flow_SumOfWeights%sPtPOI_SP",weightFlag[i].Data())));
870 pHistSumOfWeightsEtaPOI[i] = dynamic_cast<TH1D*>(outputListHistos->FindObject(Form("Flow_SumOfWeights%sEtaPOI_SP",weightFlag[i].Data())));
871 }
872
29ecccee 873 TString rpPoi[2] = {"RP","POI"};
874 TString ptEta[2] = {"Pt","Eta"};
875 TString sinCos[2] = {"sin","cos"};
876 TProfile *pHistProNonIsotropicTermsU[2][2][2] = {{{NULL}}};
a93de0f0 877 for(Int_t rp=0;rp<2;rp++) {
878 for(Int_t pe=0;pe<2;pe++) {
879 for(Int_t sc=0;sc<2;sc++) {
67d3cb61 880 pHistProNonIsotropicTermsU[rp][pe][sc] = dynamic_cast<TProfile*>(outputListHistos->FindObject(Form("FlowPro_NonIsotropicTerms_%s_%s_%s_SP",rpPoi[rp].Data(),ptEta[pe].Data(),sinCos[sc].Data())));
a93de0f0 881 }
882 }
29ecccee 883 }
884
1b302085 885 TH1D* pHistQNorm = dynamic_cast<TH1D*>(outputListHistos->FindObject("Flow_QNorm_SP"));
a554b203 886 TH1D* pHistQaQb = dynamic_cast<TH1D*>(outputListHistos->FindObject("Flow_QaQb_SP"));
887 TH1D* pHistQaQbNorm = dynamic_cast<TH1D*>(outputListHistos->FindObject("Flow_QaQbNorm_SP"));
1b302085 888 TH2D* pHistQNormvsQaQbNorm = dynamic_cast<TH2D*>(outputListHistos->FindObject("Flow_QNormvsQaQbNorm_SP"));
a554b203 889 TH1D* pHistQaQbCos = dynamic_cast<TH1D*>(outputListHistos->FindObject("Flow_QaQbCos_SP"));
890 TH1D* pHistResolution = dynamic_cast<TH1D*>(outputListHistos->FindObject("Flow_resolution_SP"));
891 TH1D* pHistQaNorm = dynamic_cast<TH1D*>(outputListHistos->FindObject("Flow_QaNorm_SP"));
892 TH2D* pHistQaNormvsMa = dynamic_cast<TH2D*>(outputListHistos->FindObject("Flow_QaNormvsMa_SP"));
893 TH1D* pHistQbNorm = dynamic_cast<TH1D*>(outputListHistos->FindObject("Flow_QbNorm_SP"));
894 TH2D* pHistQbNormvsMb = dynamic_cast<TH2D*>(outputListHistos->FindObject("Flow_QbNormvsMb_SP"));
895 TH2D* pHistMavsMb = dynamic_cast<TH2D*>(outputListHistos->FindObject("Flow_MavsMb_SP"));
896
12dc476e 897 //pass the pointers to the task
67d3cb61 898 if (pCommonHistSP &&
899 pCommonHistResultsSP &&
900 pCommonHistmuQ &&
901 pHistProQNorm &&
902 pHistProQaQb &&
903 pHistProQaQbNorm &&
904 pHistProQaQbReImNorm &&
1b302085 905 pHistProNonIsotropicTermsQ &&
906 pHistSumOfLinearWeights &&
907 pHistSumOfQuadraticWeights &&
67d3cb61 908 pHistProFlags &&
909 pHistProUQetaRP &&
910 pHistProUQetaPOI &&
911 pHistProUQPtRP &&
912 pHistProUQPtPOI &&
913 pHistProUQQaQbPtRP &&
914 pHistProUQQaQbEtaRP &&
915 pHistProUQQaQbPtPOI &&
916 pHistProUQQaQbEtaPOI &&
917 pHistSumOfWeightsPtRP[0] && pHistSumOfWeightsPtRP[1] && pHistSumOfWeightsPtRP[2] &&
918 pHistSumOfWeightsEtaRP[0] && pHistSumOfWeightsEtaRP[1] && pHistSumOfWeightsEtaRP[2] &&
919 pHistSumOfWeightsPtPOI[0] && pHistSumOfWeightsPtPOI[1] && pHistSumOfWeightsPtPOI[2] &&
920 pHistSumOfWeightsEtaPOI[0] && pHistSumOfWeightsEtaPOI[1] && pHistSumOfWeightsEtaPOI[2] &&
921 pHistProNonIsotropicTermsU[0][0][0] && pHistProNonIsotropicTermsU[1][0][0] && pHistProNonIsotropicTermsU[0][1][0] && pHistProNonIsotropicTermsU[0][0][1] &&
922 pHistProNonIsotropicTermsU[1][1][0] && pHistProNonIsotropicTermsU[1][0][1] && pHistProNonIsotropicTermsU[0][1][1] && pHistProNonIsotropicTermsU[1][1][1] &&
923 pHistQNorm &&
924 pHistQaQb &&
925 pHistQaQbNorm &&
926 pHistQNormvsQaQbNorm &&
927 pHistQaQbCos &&
928 pHistResolution &&
929 pHistQaNorm &&
930 pHistQaNormvsMa &&
931 pHistQbNorm &&
932 pHistQbNormvsMb &&
933 pHistMavsMb
934 ) {
1b302085 935
936 this -> SetCommonHistsSP(pCommonHistSP);
937 this -> SetCommonHistsResSP(pCommonHistResultsSP);
938 this -> SetCommonHistsmuQ(pCommonHistmuQ);
939 this -> SetHistProQNorm(pHistProQNorm);
fd46c3dd 940 this -> SetHistProQaQb(pHistProQaQb);
a554b203 941 this -> SetHistProQaQbNorm(pHistProQaQbNorm);
29ecccee 942 this -> SetHistProQaQbReImNorm(pHistProQaQbReImNorm);
943 this -> SetHistProNonIsotropicTermsQ(pHistProNonIsotropicTermsQ);
a554b203 944 this -> SetHistSumOfLinearWeights(pHistSumOfLinearWeights);
945 this -> SetHistSumOfQuadraticWeights(pHistSumOfQuadraticWeights);
29ecccee 946 this -> SetHistProFlags(pHistProFlags);
fd46c3dd 947 this -> SetHistProUQetaRP(pHistProUQetaRP);
948 this -> SetHistProUQetaPOI(pHistProUQetaPOI);
949 this -> SetHistProUQPtRP(pHistProUQPtRP);
950 this -> SetHistProUQPtPOI(pHistProUQPtPOI);
12dc476e 951 this -> SetHistProUQQaQbPtRP(pHistProUQQaQbPtRP);
952 this -> SetHistProUQQaQbEtaRP(pHistProUQQaQbEtaRP);
953 this -> SetHistProUQQaQbPtPOI(pHistProUQQaQbPtPOI);
67d3cb61 954 this -> SetHistProUQQaQbEtaPOI(pHistProUQQaQbEtaPOI);
a554b203 955 for(Int_t i=0;i<3;i++) {
956 if(pHistSumOfWeightsPtRP[i]) this -> SetHistSumOfWeightsPtRP(pHistSumOfWeightsPtRP[i],i);
957 if(pHistSumOfWeightsEtaRP[i]) this -> SetHistSumOfWeightsEtaRP(pHistSumOfWeightsEtaRP[i],i);
958 if(pHistSumOfWeightsPtPOI[i]) this -> SetHistSumOfWeightsPtPOI(pHistSumOfWeightsPtPOI[i],i);
959 if(pHistSumOfWeightsEtaPOI[i]) this -> SetHistSumOfWeightsEtaPOI(pHistSumOfWeightsEtaPOI[i],i);
67d3cb61 960 }
a93de0f0 961 for(Int_t rp=0;rp<2;rp++) {
962 for(Int_t pe=0;pe<2;pe++) {
963 for(Int_t sc=0;sc<2;sc++) {
964 if(pHistProNonIsotropicTermsU[rp][pe][sc]) {
965 this->SetHistProNonIsotropicTermsU(pHistProNonIsotropicTermsU[rp][pe][sc],rp,pe,sc);
966 }
967 }
968 }
29ecccee 969 }
67d3cb61 970 this -> SetHistQNorm(pHistQNorm);
971 this -> SetHistQaQb(pHistQaQb);
972 this -> SetHistQaQbNorm(pHistQaQbNorm);
973 this -> SetHistQNormvsQaQbNorm(pHistQNormvsQaQbNorm);
974 this -> SetHistQaQbCos(pHistQaQbCos);
975 this -> SetHistResolution(pHistResolution);
976 this -> SetHistQaNorm(pHistQaNorm);
977 this -> SetHistQaNormvsMa(pHistQaNormvsMa);
978 this -> SetHistQbNorm(pHistQbNorm);
979 this -> SetHistQbNormvsMb(pHistQbNormvsMb);
980 this -> SetHistMavsMb(pHistMavsMb);
a554b203 981
982 } else {
983 cout<<"WARNING: Histograms needed to run Finish() in SP are not accessable!"<<endl; }
984
12dc476e 985 } // end of if(outputListHistos)
fd46c3dd 986}
987
988//--------------------------------------------------------------------
8d312f00 989void AliFlowAnalysisWithScalarProduct::Finish() {
990
489d5531 991 //calculate flow and fill the AliFlowCommonHistResults
395fadba 992 if (fDebug) cout<<"AliFlowAnalysisWithScalarProduct::Finish()"<<endl;
29ecccee 993
50fe33aa 994 // access harmonic:
995 if(fCommonHistsSP && fCommonHistsSP->GetHarmonic())
996 {
997 fHarmonic = (Int_t)(fCommonHistsSP->GetHarmonic())->GetBinContent(1);
998 }
999
29ecccee 1000 //access all boolean flags needed in Finish():
1001 this->AccessFlags();
1002
395fadba 1003 cout<<"*************************************"<<endl;
1004 cout<<"*************************************"<<endl;
1005 cout<<" Integrated flow from "<<endl;
1006 cout<<" Scalar product "<<endl;
1007 cout<<endl;
e4cecfa0 1008
bee2efdc 1009 Int_t iNbinsPt = AliFlowCommonConstants::GetMaster()->GetNbinsPt();
1010 Int_t iNbinsEta = AliFlowCommonConstants::GetMaster()->GetNbinsEta();
12dc476e 1011
a554b203 1012 //Calculate the event plane resolution
1013 //----------------------------------
1014 Double_t dCos2phi = fHistResolution->GetMean();
7f9d91f9 1015 if (dCos2phi > 0.0){
a554b203 1016 Double_t dResolution = TMath::Sqrt(2*dCos2phi);
1017 cout<<"An estimate of the event plane resolution is: "<<dResolution<<endl;
1018 cout<<endl;
1019 }
1020
12dc476e 1021 //Calculate reference flow (noname)
1022 //----------------------------------
1023 //weighted average over (QaQb/MaMb) with weight (MaMb)
a554b203 1024 Double_t dQaQb = fHistProQaQbNorm->GetBinContent(1);
a93de0f0 1025 Double_t dSpreadQaQb = fHistProQaQbNorm->GetBinError(1);
1026 Double_t dEntriesQaQb = fHistProQaQbNorm->GetEntries();
1027
29ecccee 1028 //non-isotropic terms:
1029 Double_t dImQa = fHistProQaQbReImNorm->GetBinContent(1);
1030 Double_t dReQa = fHistProQaQbReImNorm->GetBinContent(2);
1031 Double_t dImQb = fHistProQaQbReImNorm->GetBinContent(3);
1032 Double_t dReQb = fHistProQaQbReImNorm->GetBinContent(4);
3c5d5752 1033
29ecccee 1034 if(fApplyCorrectionForNUA)
1035 {
1036 dQaQb = dQaQb - dImQa*dImQb - dReQa*dReQb;
1037 }
1038
a93de0f0 1039 if (dEntriesQaQb > 0.) {
1040 cout<<"QaQb = "<<dQaQb<<" +- "<<(dSpreadQaQb/TMath::Sqrt(dEntriesQaQb))<<endl;
1041 cout<<endl;
1042 }
1043
a554b203 1044 Double_t dV = -999.;
524798bf 1045 if(dQaQb>=0.)
1046 {
1047 dV = TMath::Sqrt(dQaQb);
1048 }
12dc476e 1049 //statistical error of dQaQb:
1050 // statistical error = term1 * spread * term2:
1051 // term1 = sqrt{sum_{i=1}^{N} w^2}/(sum_{i=1}^{N} w)
1052 // term2 = 1/sqrt(1-term1^2)
12dc476e 1053 Double_t dSumOfLinearWeights = fHistSumOfLinearWeights->GetBinContent(1);
1054 Double_t dSumOfQuadraticWeights = fHistSumOfQuadraticWeights->GetBinContent(1);
1055 Double_t dTerm1 = 0.;
1056 Double_t dTerm2 = 0.;
1057 if(dSumOfLinearWeights) {
1058 dTerm1 = pow(dSumOfQuadraticWeights,0.5)/dSumOfLinearWeights;
1059 }
1060 if(1.-pow(dTerm1,2.) > 0.) {
1061 dTerm2 = 1./pow(1-pow(dTerm1,2.),0.5);
1062 }
1063 Double_t dStatErrorQaQb = dTerm1 * dSpreadQaQb * dTerm2;
1064 //calculate the statistical error
1065 Double_t dVerr = 0.;
1066 if(dQaQb > 0.) {
1067 dVerr = (1./(2.*pow(dQaQb,0.5)))*dStatErrorQaQb;
1068 }
1b302085 1069 fCommonHistsResSP->FillIntegratedFlow(dV,dVerr);
50fe33aa 1070 cout<<Form("v%i(subevents) = ",fHarmonic)<<dV<<" +- "<<dVerr<<endl;
12dc476e 1071
1072 //Calculate differential flow and integrated flow (RP, POI)
1073 //---------------------------------------------------------
12dc476e 1074 //v as a function of eta for RP selection
a9de3704 1075 for(Int_t b=1;b<iNbinsEta+1;b++) {
12dc476e 1076 Double_t duQpro = fHistProUQetaRP->GetBinContent(b);
a93de0f0 1077 if(fApplyCorrectionForNUA) {
1078 duQpro = duQpro
1079 - fHistProNonIsotropicTermsU[0][1][1]->GetBinContent(b)*fHistProNonIsotropicTermsQ->GetBinContent(2)
1080 - fHistProNonIsotropicTermsU[0][1][0]->GetBinContent(b)*fHistProNonIsotropicTermsQ->GetBinContent(1);
29ecccee 1081 }
a554b203 1082 Double_t dv2pro = -999.;
12dc476e 1083 if (dV!=0.) { dv2pro = duQpro/dV; }
1084 //calculate the statistical error
1085 Double_t dv2ProErr = CalculateStatisticalError(b, dStatErrorQaQb, fHistProUQetaRP, fHistProUQQaQbEtaRP, fHistSumOfWeightsEtaRP);
12dc476e 1086 //fill TH1D
1b302085 1087 fCommonHistsResSP->FillDifferentialFlowEtaRP(b, dv2pro, dv2ProErr);
12dc476e 1088 } //loop over bins b
29ecccee 1089
80f72ed6 1090
12dc476e 1091 //v as a function of eta for POI selection
a9de3704 1092 for(Int_t b=1;b<iNbinsEta+1;b++) {
12dc476e 1093 Double_t duQpro = fHistProUQetaPOI->GetBinContent(b);
a93de0f0 1094 if(fApplyCorrectionForNUA) {
1095 duQpro = duQpro
1096 - fHistProNonIsotropicTermsU[1][1][1]->GetBinContent(b)*fHistProNonIsotropicTermsQ->GetBinContent(2)
1097 - fHistProNonIsotropicTermsU[1][1][0]->GetBinContent(b)*fHistProNonIsotropicTermsQ->GetBinContent(1);
29ecccee 1098 }
a554b203 1099 Double_t dv2pro = -999.;
12dc476e 1100 if (dV!=0.) { dv2pro = duQpro/dV; }
1101 //calculate the statistical error
1102 Double_t dv2ProErr = CalculateStatisticalError(b, dStatErrorQaQb, fHistProUQetaPOI, fHistProUQQaQbEtaPOI, fHistSumOfWeightsEtaPOI);
1103
1104 //fill TH1D
1b302085 1105 fCommonHistsResSP->FillDifferentialFlowEtaPOI(b, dv2pro, dv2ProErr);
12dc476e 1106 } //loop over bins b
1107
80f72ed6 1108
12dc476e 1109 //v as a function of Pt for RP selection
1b302085 1110 TH1F* fHistPtRP = fCommonHistsSP->GetHistPtRP(); //for calculating integrated flow
12dc476e 1111 Double_t dVRP = 0.;
1112 Double_t dSumRP = 0.;
1113 Double_t dErrVRP =0.;
1114
a9de3704 1115 for(Int_t b=1;b<iNbinsPt+1;b++) {
12dc476e 1116 Double_t duQpro = fHistProUQPtRP->GetBinContent(b);
a93de0f0 1117 if(fApplyCorrectionForNUA) {
1118 duQpro = duQpro
1119 - fHistProNonIsotropicTermsU[0][0][1]->GetBinContent(b)*fHistProNonIsotropicTermsQ->GetBinContent(2)
1120 - fHistProNonIsotropicTermsU[0][0][0]->GetBinContent(b)*fHistProNonIsotropicTermsQ->GetBinContent(1);
29ecccee 1121 }
a554b203 1122 Double_t dv2pro = -999.;
12dc476e 1123 if (dV!=0.) { dv2pro = duQpro/dV; }
1124 //calculate the statistical error
1125 Double_t dv2ProErr = CalculateStatisticalError(b, dStatErrorQaQb, fHistProUQPtRP, fHistProUQQaQbPtRP, fHistSumOfWeightsPtRP);
1126
1127 //fill TH1D
1b302085 1128 fCommonHistsResSP->FillDifferentialFlowPtRP(b, dv2pro, dv2ProErr);
12dc476e 1129
1130 //calculate integrated flow for RP selection
1131 if (fHistPtRP){
1132 Double_t dYieldPt = fHistPtRP->GetBinContent(b);
1133 dVRP += dv2pro*dYieldPt;
1134 dSumRP +=dYieldPt;
1135 dErrVRP += dYieldPt*dYieldPt*dv2ProErr*dv2ProErr;
1136 } else { cout<<"fHistPtRP is NULL"<<endl; }
1137 } //loop over bins b
1138
1139 if (dSumRP != 0.) {
1140 dVRP /= dSumRP; //the pt distribution should be normalised
1141 dErrVRP /= (dSumRP*dSumRP);
1142 dErrVRP = TMath::Sqrt(dErrVRP);
80f72ed6 1143 }
1b302085 1144 fCommonHistsResSP->FillIntegratedFlowRP(dVRP,dErrVRP);
50fe33aa 1145 cout<<Form("v%i(RP) = ",fHarmonic)<<dVRP<<" +- "<<dErrVRP<<endl;
12dc476e 1146
80f72ed6 1147
12dc476e 1148 //v as a function of Pt for POI selection
1b302085 1149 TH1F* fHistPtPOI = fCommonHistsSP->GetHistPtPOI(); //for calculating integrated flow
12dc476e 1150 Double_t dVPOI = 0.;
1151 Double_t dSumPOI = 0.;
1152 Double_t dErrVPOI =0.;
1153
a9de3704 1154 for(Int_t b=1;b<iNbinsPt+1;b++) {
12dc476e 1155 Double_t duQpro = fHistProUQPtPOI->GetBinContent(b);
a93de0f0 1156 if(fApplyCorrectionForNUA) {
1157 duQpro = duQpro
1158 - fHistProNonIsotropicTermsU[1][0][1]->GetBinContent(b)*fHistProNonIsotropicTermsQ->GetBinContent(2)
1159 - fHistProNonIsotropicTermsU[1][0][0]->GetBinContent(b)*fHistProNonIsotropicTermsQ->GetBinContent(1);
29ecccee 1160 }
a554b203 1161 Double_t dv2pro = -999.;
12dc476e 1162 if (dV!=0.) { dv2pro = duQpro/dV; }
1163 //calculate the statistical error
1164 Double_t dv2ProErr = CalculateStatisticalError(b, dStatErrorQaQb, fHistProUQPtPOI, fHistProUQQaQbPtPOI, fHistSumOfWeightsPtPOI);
1165
1166 //fill TH1D
1b302085 1167 fCommonHistsResSP->FillDifferentialFlowPtPOI(b, dv2pro, dv2ProErr);
395fadba 1168
12dc476e 1169 //calculate integrated flow for POI selection
1170 if (fHistPtPOI){
1171 Double_t dYieldPt = fHistPtPOI->GetBinContent(b);
1172 dVPOI += dv2pro*dYieldPt;
1173 dSumPOI +=dYieldPt;
1174 dErrVPOI += dYieldPt*dYieldPt*dv2ProErr*dv2ProErr;
1175 } else { cout<<"fHistPtPOI is NULL"<<endl; }
1176 } //loop over bins b
395fadba 1177
a93de0f0 1178 if (dSumPOI > 0.) {
489d5531 1179 dVPOI /= dSumPOI; //the pt distribution should be normalised
12dc476e 1180 dErrVPOI /= (dSumPOI*dSumPOI);
1181 dErrVPOI = TMath::Sqrt(dErrVPOI);
1182 }
1b302085 1183 fCommonHistsResSP->FillIntegratedFlowPOI(dVPOI,dErrVPOI);
50fe33aa 1184 cout<<Form("v%i(POI) = ",fHarmonic)<<dVPOI<<" +- "<<dErrVPOI<<endl;
e4cecfa0 1185
395fadba 1186 cout<<endl;
1187 cout<<"*************************************"<<endl;
1188 cout<<"*************************************"<<endl;
12dc476e 1189
395fadba 1190 //cout<<".....finished"<<endl;
12dc476e 1191}
8232a5ec 1192
1193
12dc476e 1194//--------------------------------------------------------------------
1195Double_t AliFlowAnalysisWithScalarProduct::CalculateStatisticalError(Int_t b, Double_t aStatErrorQaQb, TProfile* pHistProUQ, TProfile* pHistProUQQaQb, TH1D** pHistSumOfWeights) {
12dc476e 1196 //calculate the statistical error for differential flow for bin b
1197 Double_t duQproSpread = pHistProUQ->GetBinError(b);
1198 Double_t sumOfMq = pHistSumOfWeights[0]->GetBinContent(b);
1199 Double_t sumOfMqSquared = pHistSumOfWeights[1]->GetBinContent(b);
a554b203 1200 Double_t dQaQb = fHistProQaQbNorm->GetBinContent(1);
12dc476e 1201 Double_t dTerm1 = 0.;
1202 Double_t dTerm2 = 0.;
1203 if(sumOfMq) {
1204 dTerm1 = (pow(sumOfMqSquared,0.5)/sumOfMq);
1205 }
1206 if(1.-pow(dTerm1,2.)>0.) {
1207 dTerm2 = 1./pow(1.-pow(dTerm1,2.),0.5);
1208 }
1209 Double_t duQproErr = dTerm1*duQproSpread*dTerm2;
12dc476e 1210 // covariances:
1211 Double_t dTerm1Cov = pHistSumOfWeights[2]->GetBinContent(b);
1212 Double_t dTerm2Cov = fHistSumOfLinearWeights->GetBinContent(1);
1213 Double_t dTerm3Cov = sumOfMq;
1214 Double_t dWeightedCovariance = 0.;
1215 if(dTerm2Cov*dTerm3Cov>0.) {
1216 Double_t dDenominator = 1.-dTerm1Cov/(dTerm2Cov*dTerm3Cov);
1217 Double_t dPrefactor = dTerm1Cov/(dTerm2Cov*dTerm3Cov);
1218 if(dDenominator!=0) {
a554b203 1219 Double_t dCovariance = (pHistProUQQaQb->GetBinContent(b)-dQaQb*pHistProUQ->GetBinContent(b))/dDenominator;
12dc476e 1220 dWeightedCovariance = dCovariance*dPrefactor;
1221 }
1222 }
cff1bec4 1223
12dc476e 1224 Double_t dv2ProErr = 0.; // final statitical error:
a554b203 1225 if(dQaQb>0.) {
1226 Double_t dv2ProErrorSquared = (1./4.)*pow(dQaQb,-3.)*
12dc476e 1227 (pow(pHistProUQ->GetBinContent(b),2.)*pow(aStatErrorQaQb,2.)
a554b203 1228 + 4.*pow(dQaQb,2.)*pow(duQproErr,2.)
1229 - 4.*dQaQb*pHistProUQ->GetBinContent(b)*dWeightedCovariance);
12dc476e 1230 if(dv2ProErrorSquared>0.) dv2ProErr = pow(dv2ProErrorSquared,0.5);
1231 }
1232
1233 return dv2ProErr;
1234}
29ecccee 1235
1236
1237//--------------------------------------------------------------------
1238
1239void AliFlowAnalysisWithScalarProduct::StoreFlags()
1240{
1241 // Store all boolean flags needed in Finish() in profile fHistProFlags.
1242
1243 // Apply correction for non-uniform acceptance or not:
1244 fHistProFlags->Fill(0.5,fApplyCorrectionForNUA);
1245
1246}
1247
1248//--------------------------------------------------------------------
1249
1250void AliFlowAnalysisWithScalarProduct::AccessFlags()
1251{
1252 // Access all boolean flags needed in Finish() from profile fHistProFlags.
1253
1254 // Apply correction for non-uniform acceptance or not:
1255 fApplyCorrectionForNUA = (Bool_t) fHistProFlags->GetBinContent(1);
1256
1257}
1258
12dc476e 1259//--------------------------------------------------------------------