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