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