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