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