Smalle changes required by HLT (Jochen)
[u/mrichter/AliRoot.git] / PWG2 / FLOW / AliFlowCommon / AliFlowAnalysisWithScalarProduct.cxx
CommitLineData
8d312f00 1/*************************************************************************
2* Copyright(c) 1998-2008, ALICE Experiment at CERN, All rights reserved. *
3* *
4* Author: The ALICE Off-line Project. *
5* Contributors are mentioned in the code where appropriate. *
6* *
7* Permission to use, copy, modify and distribute this software and its *
8* documentation strictly for non-commercial purposes is hereby granted *
9* without fee, provided that the above copyright notice appears in all *
10* copies and that both the copyright notice and this permission notice *
11* appear in the supporting documentation. The authors make no claims *
12* about the suitability of this software for any purpose. It is *
13* provided "as is" without express or implied warranty. *
14**************************************************************************/
15
16#define AliFlowAnalysisWithScalarProduct_cxx
17
929098e4 18#include "Riostream.h"
19#include "TFile.h"
d7eb18ec 20#include "TList.h"
8d312f00 21#include "TMath.h"
22#include "TProfile.h"
23#include "TVector2.h"
a554b203 24#include "TH1D.h"
25#include "TH2D.h"
8d312f00 26
929098e4 27#include "AliFlowCommonConstants.h"
8d312f00 28#include "AliFlowEventSimple.h"
a93de0f0 29#include "AliFlowVector.h"
8d312f00 30#include "AliFlowTrackSimple.h"
31#include "AliFlowCommonHist.h"
395fadba 32#include "AliFlowCommonHistResults.h"
8d312f00 33#include "AliFlowAnalysisWithScalarProduct.h"
8d312f00 34
12dc476e 35//////////////////////////////////////////////////////////////////////////////
8d312f00 36// AliFlowAnalysisWithScalarProduct:
37// Description:
38// Maker to analyze Flow with the Scalar product method.
39//
12dc476e 40// authors: N. van der Kolk (kolk@nikhef.nl), A. Bilandzic (anteb@nikhef.nl)
41//////////////////////////////////////////////////////////////////////////////
8d312f00 42
43ClassImp(AliFlowAnalysisWithScalarProduct)
44
45 //-----------------------------------------------------------------------
12dc476e 46 AliFlowAnalysisWithScalarProduct::AliFlowAnalysisWithScalarProduct():
e35ddff0 47 fEventNumber(0),
48 fDebug(kFALSE),
a93de0f0 49 fApplyCorrectionForNUA(kFALSE),
a554b203 50 fRelDiffMsub(1.),
29195b69 51 fWeightsList(NULL),
52 fUsePhiWeights(kFALSE),
04c6b875 53 fPhiWeightsSub0(NULL),
54 fPhiWeightsSub1(NULL),
29ecccee 55 fHistProFlags(NULL),
395fadba 56 fHistProUQetaRP(NULL),
57 fHistProUQetaPOI(NULL),
58 fHistProUQPtRP(NULL),
59 fHistProUQPtPOI(NULL),
e4cecfa0 60 fHistProQaQb(NULL),
a554b203 61 fHistProQaQbNorm(NULL),
29ecccee 62 fHistProQaQbReImNorm(NULL),
63 fHistProNonIsotropicTermsQ(NULL),
12dc476e 64 fHistSumOfLinearWeights(NULL),
65 fHistSumOfQuadraticWeights(NULL),
66 fHistProUQQaQbPtRP(NULL),
67 fHistProUQQaQbEtaRP(NULL),
68 fHistProUQQaQbPtPOI(NULL),
69 fHistProUQQaQbEtaPOI(NULL),
395fadba 70 fCommonHists(NULL),
a554b203 71 fCommonHistsRes(NULL),
72 fHistQaQb(NULL),
73 fHistQaQbNorm(NULL),
74 fHistQaQbCos(NULL),
75 fHistResolution(NULL),
76 fHistQaNorm(NULL),
77 fHistQaNormvsMa(NULL),
78 fHistQbNorm(NULL),
79 fHistQbNormvsMb(NULL),
a93de0f0 80 fHistMavsMb(NULL),
81 fHistList(NULL)
8d312f00 82{
8d312f00 83 // Constructor.
29195b69 84 fWeightsList = new TList();
e35ddff0 85 fHistList = new TList();
12dc476e 86
87 // Initialize arrays:
88 for(Int_t i=0;i<3;i++)
89 {
90 fHistSumOfWeightsPtRP[i] = NULL;
91 fHistSumOfWeightsEtaRP[i] = NULL;
92 fHistSumOfWeightsPtPOI[i] = NULL;
93 fHistSumOfWeightsEtaPOI[i] = NULL;
94 }
29ecccee 95 for(Int_t rp=0;rp<2;rp++)
96 {
97 for(Int_t pe=0;pe<2;pe++)
98 {
99 for(Int_t sc=0;sc<2;sc++)
100 {
101 fHistProNonIsotropicTermsU[rp][pe][sc] = NULL;
102 }
103 }
104 }
8d312f00 105}
106 //-----------------------------------------------------------------------
8d312f00 107 AliFlowAnalysisWithScalarProduct::~AliFlowAnalysisWithScalarProduct()
108 {
109 //destructor
29195b69 110 delete fWeightsList;
e35ddff0 111 delete fHistList;
8d312f00 112 }
113
1dfa3c16 114//-----------------------------------------------------------------------
1dfa3c16 115void AliFlowAnalysisWithScalarProduct::WriteHistograms(TString* outputFileName)
116{
117 //store the final results in output .root file
118
119 TFile *output = new TFile(outputFileName->Data(),"RECREATE");
0fe80f88 120 //output->WriteObject(fHistList, "cobjSP","SingleKey");
121 fHistList->SetName("cobjSP");
9455e15e 122 fHistList->SetOwner(kTRUE);
0fe80f88 123 fHistList->Write(fHistList->GetName(), TObject::kSingleKey);
1dfa3c16 124 delete output;
125}
126
b0fda271 127//-----------------------------------------------------------------------
b0fda271 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);
b0fda271 137 delete output;
138}
139
ad87ae62 140//-----------------------------------------------------------------------
ad87ae62 141void AliFlowAnalysisWithScalarProduct::WriteHistograms(TDirectoryFile *outputFileName)
142{
143 //store the final results in output .root file
144 fHistList->SetName("cobjSP");
145 fHistList->SetOwner(kTRUE);
146 outputFileName->Add(fHistList);
147 outputFileName->Write(outputFileName->GetName(), TObject::kSingleKey);
148}
149
8d312f00 150//-----------------------------------------------------------------------
151void AliFlowAnalysisWithScalarProduct::Init() {
152
153 //Define all histograms
d7eb18ec 154 cout<<"---Analysis with the Scalar Product Method--- Init"<<endl;
8d312f00 155
489d5531 156 //save old value and prevent histograms from being added to directory
157 //to avoid name clashes in case multiple analaysis objects are used
158 //in an analysis
a93de0f0 159
489d5531 160 Bool_t oldHistAddStatus = TH1::AddDirectoryStatus();
161 TH1::AddDirectory(kFALSE);
162
bee2efdc 163 Int_t iNbinsPt = AliFlowCommonConstants::GetMaster()->GetNbinsPt();
164 Double_t dPtMin = AliFlowCommonConstants::GetMaster()->GetPtMin();
165 Double_t dPtMax = AliFlowCommonConstants::GetMaster()->GetPtMax();
166 Int_t iNbinsEta = AliFlowCommonConstants::GetMaster()->GetNbinsEta();
167 Double_t dEtaMin = AliFlowCommonConstants::GetMaster()->GetEtaMin();
168 Double_t dEtaMax = AliFlowCommonConstants::GetMaster()->GetEtaMax();
8d312f00 169
29ecccee 170 fHistProFlags = new TProfile("Flow_Flags_SP","Flow_Flags_SP",1,0,1,"s");
171 fHistProFlags->GetXaxis()->SetBinLabel(1,"fApplyCorrectionForNUA");
172 fHistList->Add(fHistProFlags);
173
12dc476e 174 fHistProUQetaRP = new TProfile("Flow_UQetaRP_SP","Flow_UQetaRP_SP",iNbinsEta,dEtaMin,dEtaMax,"s");
395fadba 175 fHistProUQetaRP->SetXTitle("{eta}");
176 fHistProUQetaRP->SetYTitle("<uQ>");
177 fHistList->Add(fHistProUQetaRP);
178
12dc476e 179 fHistProUQetaPOI = new TProfile("Flow_UQetaPOI_SP","Flow_UQetaPOI_SP",iNbinsEta,dEtaMin,dEtaMax,"s");
395fadba 180 fHistProUQetaPOI->SetXTitle("{eta}");
181 fHistProUQetaPOI->SetYTitle("<uQ>");
182 fHistList->Add(fHistProUQetaPOI);
183
12dc476e 184 fHistProUQPtRP = new TProfile("Flow_UQPtRP_SP","Flow_UQPtRP_SP",iNbinsPt,dPtMin,dPtMax,"s");
395fadba 185 fHistProUQPtRP->SetXTitle("p_t (GeV)");
186 fHistProUQPtRP->SetYTitle("<uQ>");
187 fHistList->Add(fHistProUQPtRP);
188
12dc476e 189 fHistProUQPtPOI = new TProfile("Flow_UQPtPOI_SP","Flow_UQPtPOI_SP",iNbinsPt,dPtMin,dPtMax,"s");
395fadba 190 fHistProUQPtPOI->SetXTitle("p_t (GeV)");
191 fHistProUQPtPOI->SetYTitle("<uQ>");
192 fHistList->Add(fHistProUQPtPOI);
193
12dc476e 194 fHistProQaQb = new TProfile("Flow_QaQb_SP","Flow_QaQb_SP", 1, 0.5, 1.5,"s");
e4cecfa0 195 fHistProQaQb->SetYTitle("<QaQb>");
29ecccee 196 fHistList->Add(fHistProQaQb);
8d312f00 197
a554b203 198 fHistProQaQbNorm = new TProfile("FlowPro_QaQbNorm_SP","FlowPro_QaQbNorm_SP", 1, 0.5, 1.5,"s");
199 fHistProQaQbNorm->SetYTitle("<QaQb/MaMb>");
200 fHistList->Add(fHistProQaQbNorm);
29ecccee 201
202 fHistProQaQbReImNorm = new TProfile("FlowPro_QaQbReImNorm_SP","FlowPro_QaQbReImNorm_SP", 4, 0.5, 4.5,"s");
203 fHistProQaQbReImNorm->GetXaxis()->SetBinLabel(1,"#LT#LTsin(#phi_{a})#GT#GT");
204 fHistProQaQbReImNorm->GetXaxis()->SetBinLabel(2,"#LT#LTcos(#phi_{a})#GT#GT");
205 fHistProQaQbReImNorm->GetXaxis()->SetBinLabel(3,"#LT#LTsin(#phi_{b})#GT#GT");
206 fHistProQaQbReImNorm->GetXaxis()->SetBinLabel(4,"#LT#LTcos(#phi_{b})#GT#GT");
207 fHistList->Add(fHistProQaQbReImNorm);
208
209 fHistProNonIsotropicTermsQ = new TProfile("FlowPro_NonIsotropicTermsQ_SP","FlowPro_NonIsotropicTermsQ_SP", 2, 0.5, 2.5,"s");
210 fHistProNonIsotropicTermsQ->GetXaxis()->SetBinLabel(1,"#LT#LTsin(#phi_{a+b})#GT#GT");
211 fHistProNonIsotropicTermsQ->GetXaxis()->SetBinLabel(2,"#LT#LTcos(#phi_{a+b})#GT#GT");
212 fHistList->Add(fHistProNonIsotropicTermsQ);
213
214 TString rpPoi[2] = {"RP","POI"};
215 TString ptEta[2] = {"Pt","Eta"};
216 TString sinCos[2] = {"sin","cos"};
217 Int_t nBinsPtEta[2] = {iNbinsPt,iNbinsEta};
218 Double_t minPtEta[2] = {dPtMin,dEtaMin};
219 Double_t maxPtEta[2] = {dPtMax,dEtaMax};
220 for(Int_t rp=0;rp<2;rp++)
221 {
222 for(Int_t pe=0;pe<2;pe++)
223 {
224 for(Int_t sc=0;sc<2;sc++)
225 {
226 fHistProNonIsotropicTermsU[rp][pe][sc] = new TProfile(Form("Flow_NonIsotropicTerms_%s_%s_%s_SP",rpPoi[rp].Data(),ptEta[pe].Data(),sinCos[sc].Data()),Form("Flow_NonIsotropicTerms_%s_%s_%s_SP",rpPoi[rp].Data(),ptEta[pe].Data(),sinCos[sc].Data()),nBinsPtEta[pe],minPtEta[pe],maxPtEta[pe]);
227 fHistList->Add(fHistProNonIsotropicTermsU[rp][pe][sc]);
228 }
229 }
230 }
231
12dc476e 232 fHistSumOfLinearWeights = new TH1D("Flow_SumOfLinearWeights_SP","Flow_SumOfLinearWeights_SP",1,-0.5, 0.5);
233 fHistSumOfLinearWeights -> SetYTitle("sum (*)");
234 fHistSumOfLinearWeights -> SetXTitle("sum (Ma*Mb)");
235 fHistList->Add(fHistSumOfLinearWeights);
236
237 fHistSumOfQuadraticWeights = new TH1D("Flow_SumOfQuadraticWeights_SP","Flow_SumOfQuadraticWeights_SP",1,-0.5, 0.5);
238 fHistSumOfQuadraticWeights -> SetYTitle("sum (*)");
239 fHistSumOfQuadraticWeights -> SetXTitle("sum (Ma*Mb)^2");
240 fHistList->Add(fHistSumOfQuadraticWeights);
241
242 fHistProUQQaQbPtRP = new TProfile("Flow_UQQaQbPtRP_SP","Flow_UQQaQbPtRP_SP",iNbinsPt,dPtMin,dPtMax);
243 fHistProUQQaQbPtRP -> SetYTitle("<*>");
244 fHistProUQQaQbPtRP -> SetXTitle("<Qu QaQb>");
245 fHistList->Add(fHistProUQQaQbPtRP);
246
247 fHistProUQQaQbEtaRP = new TProfile("Flow_UQQaQbEtaRP_SP","Flow_UQQaQbEtaRP_SP",iNbinsEta,dEtaMin,dEtaMax);
248 fHistProUQQaQbEtaRP -> SetYTitle("<*>");
249 fHistProUQQaQbEtaRP -> SetXTitle("<Qu QaQb>");
250 fHistList->Add(fHistProUQQaQbEtaRP);
251
252 fHistProUQQaQbPtPOI = new TProfile("Flow_UQQaQbPtPOI_SP","Flow_UQQaQbPtPOI_SP",iNbinsPt,dPtMin,dPtMax);
253 fHistProUQQaQbPtPOI -> SetYTitle("<*>");
254 fHistProUQQaQbPtPOI -> SetXTitle("<Qu QaQb>");
255 fHistList->Add(fHistProUQQaQbPtPOI);
256
257 fHistProUQQaQbEtaPOI = new TProfile("Flow_UQQaQbEtaPOI_SP","Flow_UQQaQbEtaPOI_SP",iNbinsEta,dEtaMin,dEtaMax);
258 fHistProUQQaQbEtaPOI -> SetYTitle("<*>");
259 fHistProUQQaQbEtaPOI -> SetXTitle("<Qu QaQb>");
260 fHistList->Add(fHistProUQQaQbEtaPOI);
261
262 TString weightFlag[3] = {"w_Qu_","w_Qu^2_","w_QuQaQb_"};
263 for(Int_t i=0;i<3;i++)
264 {
cff1bec4 265 fHistSumOfWeightsPtRP[i] = new TH1D(Form("Flow_SumOfWeights%sPtRP_SP",weightFlag[i].Data()),
266 Form("Flow_SumOfWeights%sPtRP_SP",weightFlag[i].Data()),iNbinsPt,dPtMin,dPtMax);
12dc476e 267 fHistSumOfWeightsPtRP[i] -> SetYTitle("sum (*)");
268 fHistSumOfWeightsPtRP[i] -> SetXTitle("p_{T}");
269 fHistList->Add(fHistSumOfWeightsPtRP[i]);
270
cff1bec4 271 fHistSumOfWeightsEtaRP[i] = new TH1D(Form("Flow_SumOfWeights%sEtaRP_SP",weightFlag[i].Data()),
272 Form("Flow_SumOfWeights%sEtaRP_SP",weightFlag[i].Data()),iNbinsEta,dEtaMin,dEtaMax);
12dc476e 273 fHistSumOfWeightsEtaRP[i] -> SetYTitle("sum (*)");
274 fHistSumOfWeightsEtaRP[i] -> SetXTitle("#eta");
275 fHistList->Add(fHistSumOfWeightsEtaRP[i]);
276
cff1bec4 277 fHistSumOfWeightsPtPOI[i] = new TH1D(Form("Flow_SumOfWeights%sPtPOI_SP",weightFlag[i].Data()),
278 Form("Flow_SumOfWeights%sPtPOI_SP",weightFlag[i].Data()),iNbinsPt,dPtMin,dPtMax);
12dc476e 279 fHistSumOfWeightsPtPOI[i] -> SetYTitle("sum (*)");
280 fHistSumOfWeightsPtPOI[i] -> SetXTitle("p_{T}");
281 fHistList->Add(fHistSumOfWeightsPtPOI[i]);
282
cff1bec4 283 fHistSumOfWeightsEtaPOI[i] = new TH1D(Form("Flow_SumOfWeights%sEtaPOI_SP",weightFlag[i].Data()),
284 Form("Flow_SumOfWeights%sEtaPOI_SP",weightFlag[i].Data()),iNbinsEta,dEtaMin,dEtaMax);
12dc476e 285 fHistSumOfWeightsEtaPOI[i] -> SetYTitle("sum (*)");
286 fHistSumOfWeightsEtaPOI[i] -> SetXTitle("#eta");
287 fHistList->Add(fHistSumOfWeightsEtaPOI[i]);
288 }
29ecccee 289
04f6283b 290 fCommonHists = new AliFlowCommonHist("AliFlowCommonHistSP");
099e1753 291 fHistList->Add(fCommonHists);
395fadba 292 fCommonHistsRes = new AliFlowCommonHistResults("AliFlowCommonHistResultsSP");
293 fHistList->Add(fCommonHistsRes);
294
a554b203 295 fHistQaQb = new TH1D("Flow_QaQb_SP","Flow_QaQb_SP",200,-100.,100.);
296 fHistQaQb -> SetYTitle("dN/dQaQb");
297 fHistQaQb -> SetXTitle("QaQb");
298 fHistList->Add(fHistQaQb);
299
300 fHistQaQbNorm = new TH1D("Flow_QaQbNorm_SP","Flow_QaQbNorm_SP",44,-1.1,1.1);
301 fHistQaQbNorm -> SetYTitle("dN/d(QaQb/MaMb)");
302 fHistQaQbNorm -> SetXTitle("QaQb/MaMb");
303 fHistList->Add(fHistQaQbNorm);
304
305 fHistQaQbCos = new TH1D("Flow_QaQbCos_SP","Flow_QaQbCos_SP",63,0.,TMath::Pi());
306 fHistQaQbCos -> SetYTitle("dN/d(#phi)");
307 fHistQaQbCos -> SetXTitle("#phi");
308 fHistList->Add(fHistQaQbCos);
309
310 fHistResolution = new TH1D("Flow_resolution_SP","Flow_resolution_SP",100,-1.0,1.0);
311 fHistResolution -> SetYTitle("dN/d(cos(2(#phi_a - #phi_b))");
312 fHistResolution -> SetXTitle("cos(2*(#phi_a - #phi_b))");
313 fHistList->Add(fHistResolution);
314
315 fHistQaNorm = new TH1D("Flow_QaNorm_SP","Flow_QaNorm_SP",22,0.,1.1);
316 fHistQaNorm -> SetYTitle("dN/d(|Qa/Ma|)");
317 fHistQaNorm -> SetXTitle("|Qa/Ma|");
318 fHistList->Add(fHistQaNorm);
319
320 fHistQaNormvsMa = new TH2D("Flow_QaNormvsMa_SP","Flow_QaNormvsMa_SP",100,0.,100.,22,0.,1.1);
321 fHistQaNormvsMa -> SetYTitle("|Qa/Ma|");
322 fHistQaNormvsMa -> SetXTitle("Ma");
323 fHistList->Add(fHistQaNormvsMa);
324
325 fHistQbNorm = new TH1D("Flow_QbNorm_SP","Flow_QbNorm_SP",22,0.,1.1);
326 fHistQbNorm -> SetYTitle("dN/d(|Qb/Mb|)");
327 fHistQbNorm -> SetXTitle("|Qb/Mb|");
328 fHistList->Add(fHistQbNorm);
329
330 fHistQbNormvsMb = new TH2D("Flow_QbNormvsMb_SP","Flow_QbNormvsMb_SP",100,0.,100.,22,0.,1.1);
331 fHistQbNormvsMb -> SetYTitle("|Qb/Mb|");
332 fHistQbNormvsMb -> SetXTitle("|Mb|");
333 fHistList->Add(fHistQbNormvsMb);
334
335 fHistMavsMb = new TH2D("Flow_MavsMb_SP","Flow_MavsMb_SP",100,0.,100.,100,0.,100.);
336 fHistMavsMb -> SetYTitle("Ma");
337 fHistMavsMb -> SetXTitle("Mb");
338 fHistList->Add(fHistMavsMb);
339
340
29195b69 341 //weights
342 if(fUsePhiWeights) {
343 if(!fWeightsList) {
344 cout<<"WARNING: fWeightsList is NULL in the Scalar Product method."<<endl;
345 exit(0);
346 }
04c6b875 347 if(fWeightsList->FindObject("phi_weights_sub0")) {
348 fPhiWeightsSub0 = dynamic_cast<TH1F*>
349 (fWeightsList->FindObject("phi_weights_sub0"));
350 fHistList->Add(fPhiWeightsSub0);
29195b69 351 } else {
352 cout<<"WARNING: histogram with phi weights is not accessible in Scalar Product"<<endl;
353 exit(0);
354 }
04c6b875 355 if(fWeightsList->FindObject("phi_weights_sub1")) {
356 fPhiWeightsSub1 = dynamic_cast<TH1F*>
357 (fWeightsList->FindObject("phi_weights_sub1"));
358 fHistList->Add(fPhiWeightsSub1);
359 } else {
360 cout<<"WARNING: histogram with phi weights is not accessible in Scalar Product"<<endl;
361 exit(0);
362 }
363
29195b69 364 } // end of if(fUsePhiWeights)
365
29ecccee 366 fEventNumber = 0; //set number of events to zero
367
368 //store all boolean flags needed in Finish():
369 this->StoreFlags();
489d5531 370
371 TH1::AddDirectory(oldHistAddStatus);
e2d51347 372}
373
8d312f00 374//-----------------------------------------------------------------------
8232a5ec 375void AliFlowAnalysisWithScalarProduct::Make(AliFlowEventSimple* anEvent) {
8d312f00 376
489d5531 377 //Calculate flow based on <QaQb/MaMb> = <v^2>
12dc476e 378
379 //Fill histograms
8232a5ec 380 if (anEvent) {
8d312f00 381
cbbaf54a 382 //get Q vectors for the eta-subevents
b125a454 383 AliFlowVector* vQarray = new AliFlowVector[2];
29195b69 384 if (fUsePhiWeights) {
385 anEvent->Get2Qsub(vQarray,2,fWeightsList,kTRUE);
386 } else {
387 anEvent->Get2Qsub(vQarray);
388 }
12dc476e 389 //subevent a
b125a454 390 AliFlowVector vQa = vQarray[0];
12dc476e 391 //subevent b
b125a454 392 AliFlowVector vQb = vQarray[1];
12dc476e 393
7a01f4a7 394 //check that the subevents are not empty:
12dc476e 395 Double_t dMa = vQa.GetMult();
12dc476e 396 Double_t dMb = vQb.GetMult();
a93de0f0 397 if (dMa > 0. && dMb > 0.) {
7a01f4a7 398
a554b203 399 //request that the subevent multiplicities are not too different
400 Double_t dRelDiff = TMath::Abs((dMa - dMb)/(dMa + dMb));
401 if (dRelDiff < fRelDiffMsub) {
402
3c5d5752 403 //fill control histograms
a554b203 404 fCommonHists->FillControlHistograms(anEvent);
405
406 //fill some SP control histograms
407 fHistProQaQb -> Fill(1.,vQa*vQb,1.); //Fill with weight 1 -> Weight with MaMb????
04c6b875 408 fHistQaQbCos ->Fill(TMath::ACos((vQa/vQa.Mod())*(vQb/vQb.Mod()))); //Acos(Qa*Qb) = angle
a554b203 409 fHistResolution -> Fill(TMath::Cos( vQa.Phi()- vQb.Phi() )); //vQa.Phi() returns 2*phi
410 fHistQaQb -> Fill(vQa*vQb);
411 fHistMavsMb -> Fill(dMb,dMa);
412
413 //get total Q vector = the sum of subevent a and subevent b
414 AliFlowVector vQ = vQa + vQb;
a93de0f0 415
29ecccee 416 //needed to correct for non-uniform acceptance:
a93de0f0 417 if(dMa+dMb > 0.) {
418 fHistProNonIsotropicTermsQ->Fill(1.,vQ.Y()/(dMa+dMb),dMa+dMb);
419 fHistProNonIsotropicTermsQ->Fill(2.,vQ.X()/(dMa+dMb),dMa+dMb);
29ecccee 420 }
a93de0f0 421
a554b203 422 //weight the Q vectors for the subevents by the multiplicity
a93de0f0 423 //Note: Weight Q only in the particle loop when it is clear
424 //if it should be (m-1) or M
a554b203 425 Double_t dQXa = vQa.X()/dMa;
426 Double_t dQYa = vQa.Y()/dMa;
427 vQa.Set(dQXa,dQYa);
428
429 Double_t dQXb = vQb.X()/dMb;
430 Double_t dQYb = vQb.Y()/dMb;
431 vQb.Set(dQXb,dQYb);
12dc476e 432
a554b203 433 //scalar product of the two subevents
434 Double_t dQaQb = (vQa*vQb);
435 fHistProQaQbNorm -> Fill(1.,dQaQb,dMa*dMb); //Fill (QaQb/MaMb) with weight (MaMb).
436 //needed for the error calculation:
437 fHistSumOfLinearWeights -> Fill(0.,dMa*dMb);
438 fHistSumOfQuadraticWeights -> Fill(0.,pow(dMa*dMb,2.));
a93de0f0 439 //needed for correcting non-uniform acceptance:
440 fHistProQaQbReImNorm->Fill(1.,dQYa,dMa); // to get <<sin(phi_a)>>
441 fHistProQaQbReImNorm->Fill(2.,dQXa,dMa); // to get <<cos(phi_a)>>
442 fHistProQaQbReImNorm->Fill(3.,dQYb,dMb); // to get <<sin(phi_b)>>
443 fHistProQaQbReImNorm->Fill(4.,dQXb,dMb); // to get <<cos(phi_b)>>
444
a554b203 445 //fill some SP control histograms
446 fHistQaQbNorm ->Fill(vQa*vQb);
447 fHistQaNorm ->Fill(vQa.Mod());
448 fHistQaNormvsMa->Fill(dMa,vQa.Mod());
449 fHistQbNorm ->Fill(vQb.Mod());
450 fHistQbNormvsMb->Fill(dMb,vQb.Mod());
451
452 //loop over the tracks of the event
453 AliFlowTrackSimple* pTrack = NULL;
454 Int_t iNumberOfTracks = anEvent->NumberOfTracks();
455 Double_t dMq = vQ.GetMult();
a93de0f0 456
a554b203 457 for (Int_t i=0;i<iNumberOfTracks;i++)
458 {
459 pTrack = anEvent->GetTrack(i) ;
460 if (pTrack){
461 Double_t dPhi = pTrack->Phi();
462 //set default phi weight to 1
463 Double_t dW = 1.;
464 //phi weight of pTrack
04c6b875 465 if(fUsePhiWeights && fPhiWeightsSub0 && fPhiWeightsSub1) {
466 if (pTrack->InSubevent(0) ) {
467 Int_t iNBinsPhiSub0 = fPhiWeightsSub0->GetNbinsX();
468 dW = fPhiWeightsSub0->GetBinContent(1+(Int_t)(TMath::Floor(dPhi*iNBinsPhiSub0/TMath::TwoPi())));
469 }
470 else if ( pTrack->InSubevent(1)) {
471 Int_t iNBinsPhiSub1 = fPhiWeightsSub1->GetNbinsX();
472 dW = fPhiWeightsSub1->GetBinContent(1+(Int_t)(TMath::Floor(dPhi*iNBinsPhiSub1/TMath::TwoPi())));
473 }
474 //bin = 1 + value*nbins/range
a554b203 475 //TMath::Floor rounds to the lower integer
476 }
12dc476e 477
a554b203 478 //calculate vU
479 TVector2 vU;
480 Double_t dUX = TMath::Cos(2*dPhi);
481 Double_t dUY = TMath::Sin(2*dPhi);
482 vU.Set(dUX,dUY);
483 Double_t dModulus = vU.Mod();
a93de0f0 484 if (dModulus > 0.) vU.Set(dUX/dModulus,dUY/dModulus); // make length 1
a554b203 485 else cerr<<"dModulus is zero!"<<endl;
12dc476e 486
a554b203 487 //redefine the Q vector and devide by its multiplicity
488 TVector2 vQm;
489 Double_t dQmX = 0.;
490 Double_t dQmY = 0.;
491 //subtract particle from the flowvector if used to define it
492 if (pTrack->InSubevent(0) || pTrack->InSubevent(1)) {
04c6b875 493 dQmX = (vQ.X() - dW*(pTrack->Weight())*dUX)/(dMq-1);
494 dQmY = (vQ.Y() - dW*(pTrack->Weight())*dUY)/(dMq-1);
a554b203 495 vQm.Set(dQmX,dQmY);
7a01f4a7 496
a554b203 497 //dUQ = scalar product of vU and vQm
498 Double_t dUQ = (vU * vQm);
499 Double_t dPt = pTrack->Pt();
500 Double_t dEta = pTrack->Eta();
a93de0f0 501
a554b203 502 //fill the profile histograms
503 if (pTrack->InRPSelection()) {
504 fHistProUQetaRP -> Fill(dEta,dUQ,(dMq-1)); //Fill (Qu/(Mq-1)) with weight (Mq-1)
505 //needed for the error calculation:
506 fHistProUQQaQbEtaRP -> Fill(dEta,dUQ*dQaQb,(dMq-1)*dMa*dMb); //Fill [Qu/(Mq-1)]*[QaQb/MaMb] with weight (Mq-1)MaMb
507 fHistProUQPtRP -> Fill(dPt,dUQ,(dMq-1)); //Fill (Qu/(Mq-1)) with weight (Mq-1)
508 fHistProUQQaQbPtRP -> Fill(dPt,dUQ*dQaQb,(dMq-1)*dMa*dMb); //Fill [Qu/(Mq-1)]*[QaQb/MaMb] with weight (Mq-1)MaMb
509
510 fHistSumOfWeightsEtaRP[0]->Fill(dEta,(dMq-1)); // sum of Mq-1
511 fHistSumOfWeightsEtaRP[1]->Fill(dEta,pow((dMq-1),2.));// sum of (Mq-1)^2
512 fHistSumOfWeightsEtaRP[2]->Fill(dEta,(dMq-1)*dMa*dMb);// sum of (Mq-1)*MaMb
513 fHistSumOfWeightsPtRP[0]->Fill(dPt,(dMq-1)); // sum of Mq-1
514 fHistSumOfWeightsPtRP[1]->Fill(dPt,pow((dMq-1),2.)); // sum of (Mq-1)^2
29ecccee 515 fHistSumOfWeightsPtRP[2]->Fill(dPt,(dMq-1)*dMa*dMb); // sum of (Mq-1)*MaMb
516 //nonisotropic terms:
517 fHistProNonIsotropicTermsU[0][0][0]->Fill(dPt,dUY,1.);
518 fHistProNonIsotropicTermsU[0][0][1]->Fill(dPt,dUX,1.);
519 fHistProNonIsotropicTermsU[0][1][0]->Fill(dEta,dUY,1.);
520 fHistProNonIsotropicTermsU[0][1][1]->Fill(dEta,dUX,1.);
a554b203 521 }
522 if (pTrack->InPOISelection()) {
523 fHistProUQetaPOI -> Fill(dEta,dUQ,(dMq-1));//Fill (Qu/(Mq-1)) with weight (Mq-1)
524 //needed for the error calculation:
525 fHistProUQQaQbEtaPOI -> Fill(dEta,dUQ*dQaQb,(dMq-1)*dMa*dMb); //Fill [Qu/(Mq-1)]*[QaQb/MaMb] with weight (Mq-1)MaMb
526 fHistProUQPtPOI -> Fill(dPt,dUQ,(dMq-1)); //Fill (Qu/(Mq-1)) with weight (Mq-1)
527 fHistProUQQaQbPtPOI -> Fill(dPt,dUQ*dQaQb,(dMq-1)*dMa*dMb); //Fill [Qu/(Mq-1)]*[QaQb/MaMb] with weight (Mq-1)MaMb
528
529 fHistSumOfWeightsEtaPOI[0]->Fill(dEta,(dMq-1)); // sum of Mq-1
530 fHistSumOfWeightsEtaPOI[1]->Fill(dEta,pow((dMq-1),2.));// sum of (Mq-1)^2
531 fHistSumOfWeightsEtaPOI[2]->Fill(dEta,(dMq-1)*dMa*dMb);// sum of (Mq-1)*MaMb
532 fHistSumOfWeightsPtPOI[0]->Fill(dPt,(dMq-1)); // sum of Mq-1
533 fHistSumOfWeightsPtPOI[1]->Fill(dPt,pow((dMq-1),2.)); // sum of (Mq-1)^2
29ecccee 534 fHistSumOfWeightsPtPOI[2]->Fill(dPt,(dMq-1)*dMa*dMb); // sum of (Mq-1)*MaMb
535 //nonisotropic terms:
536 fHistProNonIsotropicTermsU[1][0][0]->Fill(dPt,dUY,1.);
537 fHistProNonIsotropicTermsU[1][0][1]->Fill(dPt,dUX,1.);
538 fHistProNonIsotropicTermsU[1][1][0]->Fill(dEta,dUY,1.);
539 fHistProNonIsotropicTermsU[1][1][1]->Fill(dEta,dUX,1.);
a554b203 540 }
7a01f4a7 541
a554b203 542 } else { //do not subtract the particle from the flowvector
543 dQmX = vQ.X()/dMq;
544 dQmY = vQ.Y()/dMq;
545 vQm.Set(dQmX,dQmY);
7a01f4a7 546
a554b203 547 //dUQ = scalar product of vU and vQm
548 Double_t dUQ = (vU * vQm);
549 Double_t dPt = pTrack->Pt();
550 Double_t dEta = pTrack->Eta();
a93de0f0 551
a554b203 552 //fill the profile histograms
553 if (pTrack->InRPSelection()) {
554 fHistProUQetaRP -> Fill(dEta,dUQ,dMq); //Fill (Qu/Mq) with weight Mq
555 //needed for the error calculation:
556 fHistProUQQaQbEtaRP -> Fill(dEta,dUQ*dQaQb,dMq*dMa*dMb); //Fill [Qu/Mq]*[QaQb/MaMb] with weight Mq*MaMb
557 fHistProUQPtRP -> Fill(dPt,dUQ,dMq); //Fill (Qu/Mq) with weight Mq
558 fHistProUQQaQbPtRP -> Fill(dPt,dUQ*dQaQb,dMq*dMa*dMb); //Fill [Qu/Mq]*[QaQb/MaMb] with weight Mq*MaMb
559
560 fHistSumOfWeightsEtaRP[0]->Fill(dEta,dMq); // sum of Mq
561 fHistSumOfWeightsEtaRP[1]->Fill(dEta,pow(dMq,2.));// sum of Mq^2
562 fHistSumOfWeightsEtaRP[2]->Fill(dEta,dMq*dMa*dMb);// sum of Mq*MaMb
563 fHistSumOfWeightsPtRP[0]->Fill(dPt,dMq); // sum of Mq
564 fHistSumOfWeightsPtRP[1]->Fill(dPt,pow(dMq,2.)); // sum of Mq^2
29ecccee 565 fHistSumOfWeightsPtRP[2]->Fill(dPt,dMq*dMa*dMb); // sum of Mq*MaMb
566 //nonisotropic terms:
567 fHistProNonIsotropicTermsU[0][0][0]->Fill(dPt,dUY,1.);
568 fHistProNonIsotropicTermsU[0][0][1]->Fill(dPt,dUX,1.);
569 fHistProNonIsotropicTermsU[0][1][0]->Fill(dEta,dUY,1.);
570 fHistProNonIsotropicTermsU[0][1][1]->Fill(dEta,dUX,1.);
a554b203 571 }
572 if (pTrack->InPOISelection()) {
573 fHistProUQetaPOI -> Fill(dEta,dUQ,dMq); //Fill (Qu/Mq) with weight Mq
574 //needed for the error calculation:
575 fHistProUQQaQbEtaPOI -> Fill(dEta,dUQ*dQaQb,dMq*dMa*dMb); //Fill [Qu/Mq]*[QaQb/MaMb] with weight Mq*MaMb
576 fHistProUQPtPOI -> Fill(dPt,dUQ,dMq); //Fill (Qu/Mq) with weight Mq
577 fHistProUQQaQbPtPOI -> Fill(dPt,dUQ*dQaQb,dMq*dMa*dMb); //Fill [Qu/Mq]*[QaQb/MaMb] with weight Mq*MaMb
578
579 fHistSumOfWeightsEtaPOI[0]->Fill(dEta,dMq); // sum of Mq
580 fHistSumOfWeightsEtaPOI[1]->Fill(dEta,pow(dMq,2.));// sum of Mq^2
581 fHistSumOfWeightsEtaPOI[2]->Fill(dEta,dMq*dMa*dMb);// sum of Mq*MaMb
582 fHistSumOfWeightsPtPOI[0]->Fill(dPt,dMq); // sum of Mq
583 fHistSumOfWeightsPtPOI[1]->Fill(dPt,pow(dMq,2.)); // sum of Mq^2
584 fHistSumOfWeightsPtPOI[2]->Fill(dPt,dMq*dMa*dMb); // sum of Mq*MaMb
29ecccee 585 //nonisotropic terms:
586 fHistProNonIsotropicTermsU[1][0][0]->Fill(dPt,dUY,1.);
587 fHistProNonIsotropicTermsU[1][0][1]->Fill(dPt,dUX,1.);
588 fHistProNonIsotropicTermsU[1][1][0]->Fill(dEta,dUY,1.);
589 fHistProNonIsotropicTermsU[1][1][1]->Fill(dEta,dUX,1.);
a554b203 590 }
591 }//track not in subevents
7a01f4a7 592
a554b203 593 }//track
7a01f4a7 594
a554b203 595 }//loop over tracks
596
597 fEventNumber++;
598
599 } //difference Ma and Mb
12dc476e 600
7a01f4a7 601 }// subevents not empty
b125a454 602 delete [] vQarray;
a554b203 603
7a01f4a7 604 } //event
a554b203 605
7a01f4a7 606}//end of Make()
8d312f00 607
12dc476e 608//--------------------------------------------------------------------
fd46c3dd 609void AliFlowAnalysisWithScalarProduct::GetOutputHistograms(TList *outputListHistos){
610
611 //get pointers to all output histograms (called before Finish())
12dc476e 612
fd46c3dd 613 if (outputListHistos) {
614 //Get the common histograms from the output list
615 AliFlowCommonHist *pCommonHist = dynamic_cast<AliFlowCommonHist*>
616 (outputListHistos->FindObject("AliFlowCommonHistSP"));
617 AliFlowCommonHistResults *pCommonHistResults = dynamic_cast<AliFlowCommonHistResults*>
618 (outputListHistos->FindObject("AliFlowCommonHistResultsSP"));
619 TProfile* pHistProQaQb = dynamic_cast<TProfile*>(outputListHistos->FindObject("Flow_QaQb_SP"));
a554b203 620 TProfile* pHistProQaQbNorm = dynamic_cast<TProfile*>(outputListHistos->FindObject("FlowPro_QaQbNorm_SP"));
29ecccee 621 TProfile* pHistProQaQbReImNorm = dynamic_cast<TProfile*>(outputListHistos->FindObject("FlowPro_QaQbReImNorm_SP"));
622 TProfile* pHistProNonIsotropicTermsQ = dynamic_cast<TProfile*>(outputListHistos->FindObject("FlowPro_NonIsotropicTermsQ_SP"));
a554b203 623 TH1D* pHistSumOfLinearWeights = dynamic_cast<TH1D*>(outputListHistos->FindObject("Flow_SumOfLinearWeights_SP"));
624 TH1D* pHistSumOfQuadraticWeights = dynamic_cast<TH1D*>(outputListHistos->FindObject("Flow_SumOfQuadraticWeights_SP"));
625
29ecccee 626 TProfile* pHistProFlags = dynamic_cast<TProfile*>(outputListHistos->FindObject("Flow_Flags_SP"));
fd46c3dd 627 TProfile* pHistProUQetaRP = dynamic_cast<TProfile*>(outputListHistos->FindObject("Flow_UQetaRP_SP"));
628 TProfile* pHistProUQetaPOI = dynamic_cast<TProfile*>(outputListHistos->FindObject("Flow_UQetaPOI_SP"));
629 TProfile* pHistProUQPtRP = dynamic_cast<TProfile*>(outputListHistos->FindObject("Flow_UQPtRP_SP"));
630 TProfile* pHistProUQPtPOI = dynamic_cast<TProfile*>(outputListHistos->FindObject("Flow_UQPtPOI_SP"));
12dc476e 631 TProfile* pHistProUQQaQbPtRP = dynamic_cast<TProfile*>(outputListHistos->FindObject("Flow_UQQaQbPtRP_SP"));
632 TProfile* pHistProUQQaQbEtaRP = dynamic_cast<TProfile*>(outputListHistos->FindObject("Flow_UQQaQbEtaRP_SP"));
633 TProfile* pHistProUQQaQbPtPOI = dynamic_cast<TProfile*>(outputListHistos->FindObject("Flow_UQQaQbPtPOI_SP"));
634 TProfile* pHistProUQQaQbEtaPOI = dynamic_cast<TProfile*>(outputListHistos->FindObject("Flow_UQQaQbEtaPOI_SP"));
12dc476e 635 TString weightFlag[3] = {"w_Qu_","w_Qu^2_","w_QuQaQb_"};
636 TH1D* pHistSumOfWeightsPtRP[3] = {NULL};
637 TH1D* pHistSumOfWeightsEtaRP[3] = {NULL};
638 TH1D* pHistSumOfWeightsPtPOI[3] = {NULL};
a93de0f0 639 TH1D* pHistSumOfWeightsEtaPOI[3] = {NULL};
640
641 for(Int_t i=0;i<3;i++) {
642 pHistSumOfWeightsPtRP[i] = dynamic_cast<TH1D*>(outputListHistos->FindObject(Form("Flow_SumOfWeights%sPtRP_SP",weightFlag[i].Data())));
643 pHistSumOfWeightsEtaRP[i] = dynamic_cast<TH1D*>(outputListHistos->FindObject(Form("Flow_SumOfWeights%sEtaRP_SP",weightFlag[i].Data())));
644 pHistSumOfWeightsPtPOI[i] = dynamic_cast<TH1D*>(outputListHistos->FindObject(Form("Flow_SumOfWeights%sPtPOI_SP",weightFlag[i].Data())));
645 pHistSumOfWeightsEtaPOI[i] = dynamic_cast<TH1D*>(outputListHistos->FindObject(Form("Flow_SumOfWeights%sEtaPOI_SP",weightFlag[i].Data())));
646 }
647
29ecccee 648 TString rpPoi[2] = {"RP","POI"};
649 TString ptEta[2] = {"Pt","Eta"};
650 TString sinCos[2] = {"sin","cos"};
651 TProfile *pHistProNonIsotropicTermsU[2][2][2] = {{{NULL}}};
a93de0f0 652 for(Int_t rp=0;rp<2;rp++) {
653 for(Int_t pe=0;pe<2;pe++) {
654 for(Int_t sc=0;sc<2;sc++) {
655 pHistProNonIsotropicTermsU[rp][pe][sc] = dynamic_cast<TProfile*>(outputListHistos->FindObject(Form("Flow_NonIsotropicTerms_%s_%s_%s_SP",rpPoi[rp].Data(),ptEta[pe].Data(),sinCos[sc].Data())));
656 }
657 }
29ecccee 658 }
659
a554b203 660 TH1D* pHistQaQb = dynamic_cast<TH1D*>(outputListHistos->FindObject("Flow_QaQb_SP"));
661 TH1D* pHistQaQbNorm = dynamic_cast<TH1D*>(outputListHistos->FindObject("Flow_QaQbNorm_SP"));
662 TH1D* pHistQaQbCos = dynamic_cast<TH1D*>(outputListHistos->FindObject("Flow_QaQbCos_SP"));
663 TH1D* pHistResolution = dynamic_cast<TH1D*>(outputListHistos->FindObject("Flow_resolution_SP"));
664 TH1D* pHistQaNorm = dynamic_cast<TH1D*>(outputListHistos->FindObject("Flow_QaNorm_SP"));
665 TH2D* pHistQaNormvsMa = dynamic_cast<TH2D*>(outputListHistos->FindObject("Flow_QaNormvsMa_SP"));
666 TH1D* pHistQbNorm = dynamic_cast<TH1D*>(outputListHistos->FindObject("Flow_QbNorm_SP"));
667 TH2D* pHistQbNormvsMb = dynamic_cast<TH2D*>(outputListHistos->FindObject("Flow_QbNormvsMb_SP"));
668 TH2D* pHistMavsMb = dynamic_cast<TH2D*>(outputListHistos->FindObject("Flow_MavsMb_SP"));
669
12dc476e 670 //pass the pointers to the task
a554b203 671 if (pCommonHist && pCommonHistResults && pHistProQaQb && pHistProQaQbNorm &&
a93de0f0 672 pHistProQaQbReImNorm && pHistProNonIsotropicTermsQ &&
a554b203 673 pHistSumOfLinearWeights && pHistSumOfQuadraticWeights &&
674 pHistQaQb && pHistQaQbNorm && pHistQaQbCos && pHistResolution &&
675 pHistQaNorm && pHistQaNormvsMa && pHistQbNorm && pHistQbNormvsMb &&
29ecccee 676 pHistMavsMb && pHistProFlags &&
a554b203 677 pHistProUQetaRP && pHistProUQetaPOI &&
678 pHistProUQPtRP && pHistProUQPtPOI &&
679 pHistProUQQaQbPtRP && pHistProUQQaQbEtaRP &&
680 pHistProUQQaQbPtPOI && pHistProUQQaQbEtaPOI) {
fd46c3dd 681 this -> SetCommonHists(pCommonHist);
682 this -> SetCommonHistsRes(pCommonHistResults);
683 this -> SetHistProQaQb(pHistProQaQb);
a554b203 684 this -> SetHistProQaQbNorm(pHistProQaQbNorm);
29ecccee 685 this -> SetHistProQaQbReImNorm(pHistProQaQbReImNorm);
686 this -> SetHistProNonIsotropicTermsQ(pHistProNonIsotropicTermsQ);
a554b203 687 this -> SetHistSumOfLinearWeights(pHistSumOfLinearWeights);
688 this -> SetHistSumOfQuadraticWeights(pHistSumOfQuadraticWeights);
689 this -> SetHistQaQb(pHistQaQb);
690 this -> SetHistQaQbNorm(pHistQaQbNorm);
691 this -> SetHistQaQbCos(pHistQaQbCos);
692 this -> SetHistResolution(pHistResolution);
693 this -> SetHistQaNorm(pHistQaNorm);
694 this -> SetHistQaNormvsMa(pHistQaNormvsMa);
695 this -> SetHistQbNorm(pHistQbNorm);
696 this -> SetHistQbNormvsMb(pHistQbNormvsMb);
697 this -> SetHistMavsMb(pHistMavsMb);
29ecccee 698 this -> SetHistProFlags(pHistProFlags);
fd46c3dd 699 this -> SetHistProUQetaRP(pHistProUQetaRP);
700 this -> SetHistProUQetaPOI(pHistProUQetaPOI);
701 this -> SetHistProUQPtRP(pHistProUQPtRP);
702 this -> SetHistProUQPtPOI(pHistProUQPtPOI);
12dc476e 703 this -> SetHistProUQQaQbPtRP(pHistProUQQaQbPtRP);
704 this -> SetHistProUQQaQbEtaRP(pHistProUQQaQbEtaRP);
705 this -> SetHistProUQQaQbPtPOI(pHistProUQQaQbPtPOI);
706 this -> SetHistProUQQaQbEtaPOI(pHistProUQQaQbEtaPOI);
a554b203 707
708 for(Int_t i=0;i<3;i++) {
709 if(pHistSumOfWeightsPtRP[i]) this -> SetHistSumOfWeightsPtRP(pHistSumOfWeightsPtRP[i],i);
710 if(pHistSumOfWeightsEtaRP[i]) this -> SetHistSumOfWeightsEtaRP(pHistSumOfWeightsEtaRP[i],i);
711 if(pHistSumOfWeightsPtPOI[i]) this -> SetHistSumOfWeightsPtPOI(pHistSumOfWeightsPtPOI[i],i);
712 if(pHistSumOfWeightsEtaPOI[i]) this -> SetHistSumOfWeightsEtaPOI(pHistSumOfWeightsEtaPOI[i],i);
29ecccee 713 }
a93de0f0 714 for(Int_t rp=0;rp<2;rp++) {
715 for(Int_t pe=0;pe<2;pe++) {
716 for(Int_t sc=0;sc<2;sc++) {
717 if(pHistProNonIsotropicTermsU[rp][pe][sc]) {
718 this->SetHistProNonIsotropicTermsU(pHistProNonIsotropicTermsU[rp][pe][sc],rp,pe,sc);
719 }
720 }
721 }
29ecccee 722 }
a554b203 723
724 } else {
725 cout<<"WARNING: Histograms needed to run Finish() in SP are not accessable!"<<endl; }
726
12dc476e 727 } // end of if(outputListHistos)
fd46c3dd 728}
729
730//--------------------------------------------------------------------
8d312f00 731void AliFlowAnalysisWithScalarProduct::Finish() {
732
489d5531 733 //calculate flow and fill the AliFlowCommonHistResults
395fadba 734 if (fDebug) cout<<"AliFlowAnalysisWithScalarProduct::Finish()"<<endl;
29ecccee 735
736 //access all boolean flags needed in Finish():
737 this->AccessFlags();
738
395fadba 739 cout<<"*************************************"<<endl;
740 cout<<"*************************************"<<endl;
741 cout<<" Integrated flow from "<<endl;
742 cout<<" Scalar product "<<endl;
743 cout<<endl;
e4cecfa0 744
bee2efdc 745 Int_t iNbinsPt = AliFlowCommonConstants::GetMaster()->GetNbinsPt();
746 Int_t iNbinsEta = AliFlowCommonConstants::GetMaster()->GetNbinsEta();
12dc476e 747
a554b203 748 //Calculate the event plane resolution
749 //----------------------------------
750 Double_t dCos2phi = fHistResolution->GetMean();
7f9d91f9 751 if (dCos2phi > 0.0){
a554b203 752 Double_t dResolution = TMath::Sqrt(2*dCos2phi);
753 cout<<"An estimate of the event plane resolution is: "<<dResolution<<endl;
754 cout<<endl;
755 }
756
12dc476e 757 //Calculate reference flow (noname)
758 //----------------------------------
759 //weighted average over (QaQb/MaMb) with weight (MaMb)
a554b203 760 Double_t dQaQb = fHistProQaQbNorm->GetBinContent(1);
a93de0f0 761 Double_t dSpreadQaQb = fHistProQaQbNorm->GetBinError(1);
762 Double_t dEntriesQaQb = fHistProQaQbNorm->GetEntries();
763
29ecccee 764 //non-isotropic terms:
765 Double_t dImQa = fHistProQaQbReImNorm->GetBinContent(1);
766 Double_t dReQa = fHistProQaQbReImNorm->GetBinContent(2);
767 Double_t dImQb = fHistProQaQbReImNorm->GetBinContent(3);
768 Double_t dReQb = fHistProQaQbReImNorm->GetBinContent(4);
3c5d5752 769
29ecccee 770 if(fApplyCorrectionForNUA)
771 {
772 dQaQb = dQaQb - dImQa*dImQb - dReQa*dReQb;
773 }
774
a93de0f0 775 if (dEntriesQaQb > 0.) {
776 cout<<"QaQb = "<<dQaQb<<" +- "<<(dSpreadQaQb/TMath::Sqrt(dEntriesQaQb))<<endl;
777 cout<<endl;
778 }
779
a554b203 780 Double_t dV = -999.;
524798bf 781 if(dQaQb>=0.)
782 {
783 dV = TMath::Sqrt(dQaQb);
784 }
12dc476e 785 //statistical error of dQaQb:
786 // statistical error = term1 * spread * term2:
787 // term1 = sqrt{sum_{i=1}^{N} w^2}/(sum_{i=1}^{N} w)
788 // term2 = 1/sqrt(1-term1^2)
12dc476e 789 Double_t dSumOfLinearWeights = fHistSumOfLinearWeights->GetBinContent(1);
790 Double_t dSumOfQuadraticWeights = fHistSumOfQuadraticWeights->GetBinContent(1);
791 Double_t dTerm1 = 0.;
792 Double_t dTerm2 = 0.;
793 if(dSumOfLinearWeights) {
794 dTerm1 = pow(dSumOfQuadraticWeights,0.5)/dSumOfLinearWeights;
795 }
796 if(1.-pow(dTerm1,2.) > 0.) {
797 dTerm2 = 1./pow(1-pow(dTerm1,2.),0.5);
798 }
799 Double_t dStatErrorQaQb = dTerm1 * dSpreadQaQb * dTerm2;
800 //calculate the statistical error
801 Double_t dVerr = 0.;
802 if(dQaQb > 0.) {
803 dVerr = (1./(2.*pow(dQaQb,0.5)))*dStatErrorQaQb;
804 }
805 fCommonHistsRes->FillIntegratedFlow(dV,dVerr);
a93de0f0 806 cout<<"v2(subevents) = "<<dV<<" +- "<<dVerr<<endl;
12dc476e 807
808 //Calculate differential flow and integrated flow (RP, POI)
809 //---------------------------------------------------------
12dc476e 810 //v as a function of eta for RP selection
a9de3704 811 for(Int_t b=1;b<iNbinsEta+1;b++) {
12dc476e 812 Double_t duQpro = fHistProUQetaRP->GetBinContent(b);
a93de0f0 813 if(fApplyCorrectionForNUA) {
814 duQpro = duQpro
815 - fHistProNonIsotropicTermsU[0][1][1]->GetBinContent(b)*fHistProNonIsotropicTermsQ->GetBinContent(2)
816 - fHistProNonIsotropicTermsU[0][1][0]->GetBinContent(b)*fHistProNonIsotropicTermsQ->GetBinContent(1);
29ecccee 817 }
a554b203 818 Double_t dv2pro = -999.;
12dc476e 819 if (dV!=0.) { dv2pro = duQpro/dV; }
820 //calculate the statistical error
821 Double_t dv2ProErr = CalculateStatisticalError(b, dStatErrorQaQb, fHistProUQetaRP, fHistProUQQaQbEtaRP, fHistSumOfWeightsEtaRP);
12dc476e 822 //fill TH1D
cff1bec4 823 fCommonHistsRes->FillDifferentialFlowEtaRP(b, dv2pro, dv2ProErr);
12dc476e 824 } //loop over bins b
29ecccee 825
80f72ed6 826
12dc476e 827 //v as a function of eta for POI selection
a9de3704 828 for(Int_t b=1;b<iNbinsEta+1;b++) {
12dc476e 829 Double_t duQpro = fHistProUQetaPOI->GetBinContent(b);
a93de0f0 830 if(fApplyCorrectionForNUA) {
831 duQpro = duQpro
832 - fHistProNonIsotropicTermsU[1][1][1]->GetBinContent(b)*fHistProNonIsotropicTermsQ->GetBinContent(2)
833 - fHistProNonIsotropicTermsU[1][1][0]->GetBinContent(b)*fHistProNonIsotropicTermsQ->GetBinContent(1);
29ecccee 834 }
a554b203 835 Double_t dv2pro = -999.;
12dc476e 836 if (dV!=0.) { dv2pro = duQpro/dV; }
837 //calculate the statistical error
838 Double_t dv2ProErr = CalculateStatisticalError(b, dStatErrorQaQb, fHistProUQetaPOI, fHistProUQQaQbEtaPOI, fHistSumOfWeightsEtaPOI);
839
840 //fill TH1D
841 fCommonHistsRes->FillDifferentialFlowEtaPOI(b, dv2pro, dv2ProErr);
842 } //loop over bins b
843
80f72ed6 844
12dc476e 845 //v as a function of Pt for RP selection
846 TH1F* fHistPtRP = fCommonHists->GetHistPtRP(); //for calculating integrated flow
847 Double_t dVRP = 0.;
848 Double_t dSumRP = 0.;
849 Double_t dErrVRP =0.;
850
a9de3704 851 for(Int_t b=1;b<iNbinsPt+1;b++) {
12dc476e 852 Double_t duQpro = fHistProUQPtRP->GetBinContent(b);
a93de0f0 853 if(fApplyCorrectionForNUA) {
854 duQpro = duQpro
855 - fHistProNonIsotropicTermsU[0][0][1]->GetBinContent(b)*fHistProNonIsotropicTermsQ->GetBinContent(2)
856 - fHistProNonIsotropicTermsU[0][0][0]->GetBinContent(b)*fHistProNonIsotropicTermsQ->GetBinContent(1);
29ecccee 857 }
a554b203 858 Double_t dv2pro = -999.;
12dc476e 859 if (dV!=0.) { dv2pro = duQpro/dV; }
860 //calculate the statistical error
861 Double_t dv2ProErr = CalculateStatisticalError(b, dStatErrorQaQb, fHistProUQPtRP, fHistProUQQaQbPtRP, fHistSumOfWeightsPtRP);
862
863 //fill TH1D
864 fCommonHistsRes->FillDifferentialFlowPtRP(b, dv2pro, dv2ProErr);
865
866 //calculate integrated flow for RP selection
867 if (fHistPtRP){
868 Double_t dYieldPt = fHistPtRP->GetBinContent(b);
869 dVRP += dv2pro*dYieldPt;
870 dSumRP +=dYieldPt;
871 dErrVRP += dYieldPt*dYieldPt*dv2ProErr*dv2ProErr;
872 } else { cout<<"fHistPtRP is NULL"<<endl; }
873 } //loop over bins b
874
875 if (dSumRP != 0.) {
876 dVRP /= dSumRP; //the pt distribution should be normalised
877 dErrVRP /= (dSumRP*dSumRP);
878 dErrVRP = TMath::Sqrt(dErrVRP);
80f72ed6 879 }
12dc476e 880 fCommonHistsRes->FillIntegratedFlowRP(dVRP,dErrVRP);
a93de0f0 881 cout<<"v2(RP) = "<<dVRP<<" +- "<<dErrVRP<<endl;
12dc476e 882
80f72ed6 883
12dc476e 884 //v as a function of Pt for POI selection
885 TH1F* fHistPtPOI = fCommonHists->GetHistPtPOI(); //for calculating integrated flow
886 Double_t dVPOI = 0.;
887 Double_t dSumPOI = 0.;
888 Double_t dErrVPOI =0.;
889
a9de3704 890 for(Int_t b=1;b<iNbinsPt+1;b++) {
12dc476e 891 Double_t duQpro = fHistProUQPtPOI->GetBinContent(b);
a93de0f0 892 if(fApplyCorrectionForNUA) {
893 duQpro = duQpro
894 - fHistProNonIsotropicTermsU[1][0][1]->GetBinContent(b)*fHistProNonIsotropicTermsQ->GetBinContent(2)
895 - fHistProNonIsotropicTermsU[1][0][0]->GetBinContent(b)*fHistProNonIsotropicTermsQ->GetBinContent(1);
29ecccee 896 }
a554b203 897 Double_t dv2pro = -999.;
12dc476e 898 if (dV!=0.) { dv2pro = duQpro/dV; }
899 //calculate the statistical error
900 Double_t dv2ProErr = CalculateStatisticalError(b, dStatErrorQaQb, fHistProUQPtPOI, fHistProUQQaQbPtPOI, fHistSumOfWeightsPtPOI);
901
902 //fill TH1D
903 fCommonHistsRes->FillDifferentialFlowPtPOI(b, dv2pro, dv2ProErr);
395fadba 904
12dc476e 905 //calculate integrated flow for POI selection
906 if (fHistPtPOI){
907 Double_t dYieldPt = fHistPtPOI->GetBinContent(b);
908 dVPOI += dv2pro*dYieldPt;
909 dSumPOI +=dYieldPt;
910 dErrVPOI += dYieldPt*dYieldPt*dv2ProErr*dv2ProErr;
911 } else { cout<<"fHistPtPOI is NULL"<<endl; }
912 } //loop over bins b
395fadba 913
a93de0f0 914 if (dSumPOI > 0.) {
489d5531 915 dVPOI /= dSumPOI; //the pt distribution should be normalised
12dc476e 916 dErrVPOI /= (dSumPOI*dSumPOI);
917 dErrVPOI = TMath::Sqrt(dErrVPOI);
918 }
919 fCommonHistsRes->FillIntegratedFlowPOI(dVPOI,dErrVPOI);
a93de0f0 920 cout<<"v2(POI) = "<<dVPOI<<" +- "<<dErrVPOI<<endl;
e4cecfa0 921
395fadba 922 cout<<endl;
923 cout<<"*************************************"<<endl;
924 cout<<"*************************************"<<endl;
12dc476e 925
395fadba 926 //cout<<".....finished"<<endl;
12dc476e 927}
8232a5ec 928
929
12dc476e 930//--------------------------------------------------------------------
931Double_t AliFlowAnalysisWithScalarProduct::CalculateStatisticalError(Int_t b, Double_t aStatErrorQaQb, TProfile* pHistProUQ, TProfile* pHistProUQQaQb, TH1D** pHistSumOfWeights) {
12dc476e 932 //calculate the statistical error for differential flow for bin b
933 Double_t duQproSpread = pHistProUQ->GetBinError(b);
934 Double_t sumOfMq = pHistSumOfWeights[0]->GetBinContent(b);
935 Double_t sumOfMqSquared = pHistSumOfWeights[1]->GetBinContent(b);
a554b203 936 Double_t dQaQb = fHistProQaQbNorm->GetBinContent(1);
12dc476e 937 Double_t dTerm1 = 0.;
938 Double_t dTerm2 = 0.;
939 if(sumOfMq) {
940 dTerm1 = (pow(sumOfMqSquared,0.5)/sumOfMq);
941 }
942 if(1.-pow(dTerm1,2.)>0.) {
943 dTerm2 = 1./pow(1.-pow(dTerm1,2.),0.5);
944 }
945 Double_t duQproErr = dTerm1*duQproSpread*dTerm2;
12dc476e 946 // covariances:
947 Double_t dTerm1Cov = pHistSumOfWeights[2]->GetBinContent(b);
948 Double_t dTerm2Cov = fHistSumOfLinearWeights->GetBinContent(1);
949 Double_t dTerm3Cov = sumOfMq;
950 Double_t dWeightedCovariance = 0.;
951 if(dTerm2Cov*dTerm3Cov>0.) {
952 Double_t dDenominator = 1.-dTerm1Cov/(dTerm2Cov*dTerm3Cov);
953 Double_t dPrefactor = dTerm1Cov/(dTerm2Cov*dTerm3Cov);
954 if(dDenominator!=0) {
a554b203 955 Double_t dCovariance = (pHistProUQQaQb->GetBinContent(b)-dQaQb*pHistProUQ->GetBinContent(b))/dDenominator;
12dc476e 956 dWeightedCovariance = dCovariance*dPrefactor;
957 }
958 }
cff1bec4 959
12dc476e 960 Double_t dv2ProErr = 0.; // final statitical error:
a554b203 961 if(dQaQb>0.) {
962 Double_t dv2ProErrorSquared = (1./4.)*pow(dQaQb,-3.)*
12dc476e 963 (pow(pHistProUQ->GetBinContent(b),2.)*pow(aStatErrorQaQb,2.)
a554b203 964 + 4.*pow(dQaQb,2.)*pow(duQproErr,2.)
965 - 4.*dQaQb*pHistProUQ->GetBinContent(b)*dWeightedCovariance);
12dc476e 966 if(dv2ProErrorSquared>0.) dv2ProErr = pow(dv2ProErrorSquared,0.5);
967 }
968
969 return dv2ProErr;
970}
29ecccee 971
972
973//--------------------------------------------------------------------
974
975void AliFlowAnalysisWithScalarProduct::StoreFlags()
976{
977 // Store all boolean flags needed in Finish() in profile fHistProFlags.
978
979 // Apply correction for non-uniform acceptance or not:
980 fHistProFlags->Fill(0.5,fApplyCorrectionForNUA);
981
982}
983
984//--------------------------------------------------------------------
985
986void AliFlowAnalysisWithScalarProduct::AccessFlags()
987{
988 // Access all boolean flags needed in Finish() from profile fHistProFlags.
989
990 // Apply correction for non-uniform acceptance or not:
991 fApplyCorrectionForNUA = (Bool_t) fHistProFlags->GetBinContent(1);
992
993}
994
12dc476e 995//--------------------------------------------------------------------