compiler warnings (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"
29#include "AliFlowTrackSimple.h"
30#include "AliFlowCommonHist.h"
395fadba 31#include "AliFlowCommonHistResults.h"
8d312f00 32#include "AliFlowAnalysisWithScalarProduct.h"
929098e4 33#include "AliFlowVector.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),
a554b203 49 fRelDiffMsub(1.),
29195b69 50 fWeightsList(NULL),
51 fUsePhiWeights(kFALSE),
04c6b875 52 fPhiWeightsSub0(NULL),
53 fPhiWeightsSub1(NULL),
e35ddff0 54 fHistList(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 fApplyCorrectionForNUA(kFALSE),
63 fHistProQaQbReImNorm(NULL),
64 fHistProNonIsotropicTermsQ(NULL),
12dc476e 65 fHistSumOfLinearWeights(NULL),
66 fHistSumOfQuadraticWeights(NULL),
67 fHistProUQQaQbPtRP(NULL),
68 fHistProUQQaQbEtaRP(NULL),
69 fHistProUQQaQbPtPOI(NULL),
70 fHistProUQQaQbEtaPOI(NULL),
395fadba 71 fCommonHists(NULL),
a554b203 72 fCommonHistsRes(NULL),
73 fHistQaQb(NULL),
74 fHistQaQbNorm(NULL),
75 fHistQaQbCos(NULL),
76 fHistResolution(NULL),
77 fHistQaNorm(NULL),
78 fHistQaNormvsMa(NULL),
79 fHistQbNorm(NULL),
80 fHistQbNormvsMb(NULL),
81 fHistMavsMb(NULL)
82
8d312f00 83{
8d312f00 84 // Constructor.
29195b69 85 fWeightsList = new TList();
e35ddff0 86 fHistList = new TList();
12dc476e 87
88 // Initialize arrays:
89 for(Int_t i=0;i<3;i++)
90 {
91 fHistSumOfWeightsPtRP[i] = NULL;
92 fHistSumOfWeightsEtaRP[i] = NULL;
93 fHistSumOfWeightsPtPOI[i] = NULL;
94 fHistSumOfWeightsEtaPOI[i] = NULL;
95 }
29ecccee 96 for(Int_t rp=0;rp<2;rp++)
97 {
98 for(Int_t pe=0;pe<2;pe++)
99 {
100 for(Int_t sc=0;sc<2;sc++)
101 {
102 fHistProNonIsotropicTermsU[rp][pe][sc] = NULL;
103 }
104 }
105 }
8d312f00 106}
107 //-----------------------------------------------------------------------
8d312f00 108 AliFlowAnalysisWithScalarProduct::~AliFlowAnalysisWithScalarProduct()
109 {
110 //destructor
29195b69 111 delete fWeightsList;
e35ddff0 112 delete fHistList;
8d312f00 113 }
114
1dfa3c16 115//-----------------------------------------------------------------------
1dfa3c16 116void AliFlowAnalysisWithScalarProduct::WriteHistograms(TString* outputFileName)
117{
118 //store the final results in output .root file
119
120 TFile *output = new TFile(outputFileName->Data(),"RECREATE");
0fe80f88 121 //output->WriteObject(fHistList, "cobjSP","SingleKey");
122 fHistList->SetName("cobjSP");
9455e15e 123 fHistList->SetOwner(kTRUE);
0fe80f88 124 fHistList->Write(fHistList->GetName(), TObject::kSingleKey);
1dfa3c16 125 delete output;
126}
127
b0fda271 128//-----------------------------------------------------------------------
b0fda271 129void AliFlowAnalysisWithScalarProduct::WriteHistograms(TString outputFileName)
130{
131 //store the final results in output .root file
132
133 TFile *output = new TFile(outputFileName.Data(),"RECREATE");
0fe80f88 134 //output->WriteObject(fHistList, "cobjSP","SingleKey");
135 fHistList->SetName("cobjSP");
9455e15e 136 fHistList->SetOwner(kTRUE);
0fe80f88 137 fHistList->Write(fHistList->GetName(), TObject::kSingleKey);
b0fda271 138 delete output;
139}
140
ad87ae62 141//-----------------------------------------------------------------------
ad87ae62 142void AliFlowAnalysisWithScalarProduct::WriteHistograms(TDirectoryFile *outputFileName)
143{
144 //store the final results in output .root file
145 fHistList->SetName("cobjSP");
146 fHistList->SetOwner(kTRUE);
147 outputFileName->Add(fHistList);
148 outputFileName->Write(outputFileName->GetName(), TObject::kSingleKey);
149}
150
8d312f00 151//-----------------------------------------------------------------------
152void AliFlowAnalysisWithScalarProduct::Init() {
153
154 //Define all histograms
d7eb18ec 155 cout<<"---Analysis with the Scalar Product Method--- Init"<<endl;
8d312f00 156
489d5531 157 //save old value and prevent histograms from being added to directory
158 //to avoid name clashes in case multiple analaysis objects are used
159 //in an analysis
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();
7a01f4a7 397 if (dMa != 0 && dMb != 0) {
398
a554b203 399
400 //request that the subevent multiplicities are not too different
401 Double_t dRelDiff = TMath::Abs((dMa - dMb)/(dMa + dMb));
402 if (dRelDiff < fRelDiffMsub) {
403
3c5d5752 404 //fill control histograms
a554b203 405 fCommonHists->FillControlHistograms(anEvent);
406
407 //fill some SP control histograms
408 fHistProQaQb -> Fill(1.,vQa*vQb,1.); //Fill with weight 1 -> Weight with MaMb????
04c6b875 409 fHistQaQbCos ->Fill(TMath::ACos((vQa/vQa.Mod())*(vQb/vQb.Mod()))); //Acos(Qa*Qb) = angle
a554b203 410 fHistResolution -> Fill(TMath::Cos( vQa.Phi()- vQb.Phi() )); //vQa.Phi() returns 2*phi
411 fHistQaQb -> Fill(vQa*vQb);
412 fHistMavsMb -> Fill(dMb,dMa);
413
414 //get total Q vector = the sum of subevent a and subevent b
415 AliFlowVector vQ = vQa + vQb;
29ecccee 416 //needed to correct for non-uniform acceptance:
417 if(dMa+dMb>0.)
418 {
419 fHistProNonIsotropicTermsQ->Fill(1.,vQ.Y()/(dMa+dMb),dMa+dMb);
420 fHistProNonIsotropicTermsQ->Fill(2.,vQ.X()/(dMa+dMb),dMa+dMb);
421 }
a554b203 422 //weight the Q vectors for the subevents by the multiplicity
423 //Note: Weight Q only in the particle loop when it is clear if it should be (m-1) or M
424 Double_t dQXa = vQa.X()/dMa;
425 Double_t dQYa = vQa.Y()/dMa;
426 vQa.Set(dQXa,dQYa);
427
428 Double_t dQXb = vQb.X()/dMb;
429 Double_t dQYb = vQb.Y()/dMb;
430 vQb.Set(dQXb,dQYb);
12dc476e 431
a554b203 432 //scalar product of the two subevents
433 Double_t dQaQb = (vQa*vQb);
434 fHistProQaQbNorm -> Fill(1.,dQaQb,dMa*dMb); //Fill (QaQb/MaMb) with weight (MaMb).
435 //needed for the error calculation:
436 fHistSumOfLinearWeights -> Fill(0.,dMa*dMb);
437 fHistSumOfQuadraticWeights -> Fill(0.,pow(dMa*dMb,2.));
29ecccee 438 //needed for correcting non-uniform acceptance:
439 fHistProQaQbReImNorm->Fill(1.,dQYa,dMa); // to get <<sin(phi_a)>>
440 fHistProQaQbReImNorm->Fill(2.,dQXa,dMa); // to get <<cos(phi_a)>>
441 fHistProQaQbReImNorm->Fill(3.,dQYb,dMb); // to get <<sin(phi_b)>>
442 fHistProQaQbReImNorm->Fill(4.,dQXb,dMb); // to get <<cos(phi_b)>>
a554b203 443
444 //fill some SP control histograms
445 fHistQaQbNorm ->Fill(vQa*vQb);
446 fHistQaNorm ->Fill(vQa.Mod());
447 fHistQaNormvsMa->Fill(dMa,vQa.Mod());
448 fHistQbNorm ->Fill(vQb.Mod());
449 fHistQbNormvsMb->Fill(dMb,vQb.Mod());
450
451 //loop over the tracks of the event
452 AliFlowTrackSimple* pTrack = NULL;
453 Int_t iNumberOfTracks = anEvent->NumberOfTracks();
454 Double_t dMq = vQ.GetMult();
7a01f4a7 455
a554b203 456 for (Int_t i=0;i<iNumberOfTracks;i++)
457 {
458 pTrack = anEvent->GetTrack(i) ;
459 if (pTrack){
460 Double_t dPhi = pTrack->Phi();
461 //set default phi weight to 1
462 Double_t dW = 1.;
463 //phi weight of pTrack
04c6b875 464 if(fUsePhiWeights && fPhiWeightsSub0 && fPhiWeightsSub1) {
465 if (pTrack->InSubevent(0) ) {
466 Int_t iNBinsPhiSub0 = fPhiWeightsSub0->GetNbinsX();
467 dW = fPhiWeightsSub0->GetBinContent(1+(Int_t)(TMath::Floor(dPhi*iNBinsPhiSub0/TMath::TwoPi())));
468 }
469 else if ( pTrack->InSubevent(1)) {
470 Int_t iNBinsPhiSub1 = fPhiWeightsSub1->GetNbinsX();
471 dW = fPhiWeightsSub1->GetBinContent(1+(Int_t)(TMath::Floor(dPhi*iNBinsPhiSub1/TMath::TwoPi())));
472 }
473 //bin = 1 + value*nbins/range
a554b203 474 //TMath::Floor rounds to the lower integer
475 }
12dc476e 476
a554b203 477 //calculate vU
478 TVector2 vU;
479 Double_t dUX = TMath::Cos(2*dPhi);
480 Double_t dUY = TMath::Sin(2*dPhi);
481 vU.Set(dUX,dUY);
482 Double_t dModulus = vU.Mod();
483 if (dModulus!=0.) vU.Set(dUX/dModulus,dUY/dModulus); // make length 1
484 else cerr<<"dModulus is zero!"<<endl;
12dc476e 485
a554b203 486 //redefine the Q vector and devide by its multiplicity
487 TVector2 vQm;
488 Double_t dQmX = 0.;
489 Double_t dQmY = 0.;
490 //subtract particle from the flowvector if used to define it
491 if (pTrack->InSubevent(0) || pTrack->InSubevent(1)) {
04c6b875 492 dQmX = (vQ.X() - dW*(pTrack->Weight())*dUX)/(dMq-1);
493 dQmY = (vQ.Y() - dW*(pTrack->Weight())*dUY)/(dMq-1);
a554b203 494 vQm.Set(dQmX,dQmY);
7a01f4a7 495
a554b203 496 //dUQ = scalar product of vU and vQm
497 Double_t dUQ = (vU * vQm);
498 Double_t dPt = pTrack->Pt();
499 Double_t dEta = pTrack->Eta();
500 //fill the profile histograms
501 if (pTrack->InRPSelection()) {
502 fHistProUQetaRP -> Fill(dEta,dUQ,(dMq-1)); //Fill (Qu/(Mq-1)) with weight (Mq-1)
503 //needed for the error calculation:
504 fHistProUQQaQbEtaRP -> Fill(dEta,dUQ*dQaQb,(dMq-1)*dMa*dMb); //Fill [Qu/(Mq-1)]*[QaQb/MaMb] with weight (Mq-1)MaMb
505 fHistProUQPtRP -> Fill(dPt,dUQ,(dMq-1)); //Fill (Qu/(Mq-1)) with weight (Mq-1)
506 fHistProUQQaQbPtRP -> Fill(dPt,dUQ*dQaQb,(dMq-1)*dMa*dMb); //Fill [Qu/(Mq-1)]*[QaQb/MaMb] with weight (Mq-1)MaMb
507
508 fHistSumOfWeightsEtaRP[0]->Fill(dEta,(dMq-1)); // sum of Mq-1
509 fHistSumOfWeightsEtaRP[1]->Fill(dEta,pow((dMq-1),2.));// sum of (Mq-1)^2
510 fHistSumOfWeightsEtaRP[2]->Fill(dEta,(dMq-1)*dMa*dMb);// sum of (Mq-1)*MaMb
511 fHistSumOfWeightsPtRP[0]->Fill(dPt,(dMq-1)); // sum of Mq-1
512 fHistSumOfWeightsPtRP[1]->Fill(dPt,pow((dMq-1),2.)); // sum of (Mq-1)^2
29ecccee 513 fHistSumOfWeightsPtRP[2]->Fill(dPt,(dMq-1)*dMa*dMb); // sum of (Mq-1)*MaMb
514 //nonisotropic terms:
515 fHistProNonIsotropicTermsU[0][0][0]->Fill(dPt,dUY,1.);
516 fHistProNonIsotropicTermsU[0][0][1]->Fill(dPt,dUX,1.);
517 fHistProNonIsotropicTermsU[0][1][0]->Fill(dEta,dUY,1.);
518 fHistProNonIsotropicTermsU[0][1][1]->Fill(dEta,dUX,1.);
a554b203 519 }
520 if (pTrack->InPOISelection()) {
521 fHistProUQetaPOI -> Fill(dEta,dUQ,(dMq-1));//Fill (Qu/(Mq-1)) with weight (Mq-1)
522 //needed for the error calculation:
523 fHistProUQQaQbEtaPOI -> Fill(dEta,dUQ*dQaQb,(dMq-1)*dMa*dMb); //Fill [Qu/(Mq-1)]*[QaQb/MaMb] with weight (Mq-1)MaMb
524 fHistProUQPtPOI -> Fill(dPt,dUQ,(dMq-1)); //Fill (Qu/(Mq-1)) with weight (Mq-1)
525 fHistProUQQaQbPtPOI -> Fill(dPt,dUQ*dQaQb,(dMq-1)*dMa*dMb); //Fill [Qu/(Mq-1)]*[QaQb/MaMb] with weight (Mq-1)MaMb
526
527 fHistSumOfWeightsEtaPOI[0]->Fill(dEta,(dMq-1)); // sum of Mq-1
528 fHistSumOfWeightsEtaPOI[1]->Fill(dEta,pow((dMq-1),2.));// sum of (Mq-1)^2
529 fHistSumOfWeightsEtaPOI[2]->Fill(dEta,(dMq-1)*dMa*dMb);// sum of (Mq-1)*MaMb
530 fHistSumOfWeightsPtPOI[0]->Fill(dPt,(dMq-1)); // sum of Mq-1
531 fHistSumOfWeightsPtPOI[1]->Fill(dPt,pow((dMq-1),2.)); // sum of (Mq-1)^2
29ecccee 532 fHistSumOfWeightsPtPOI[2]->Fill(dPt,(dMq-1)*dMa*dMb); // sum of (Mq-1)*MaMb
533 //nonisotropic terms:
534 fHistProNonIsotropicTermsU[1][0][0]->Fill(dPt,dUY,1.);
535 fHistProNonIsotropicTermsU[1][0][1]->Fill(dPt,dUX,1.);
536 fHistProNonIsotropicTermsU[1][1][0]->Fill(dEta,dUY,1.);
537 fHistProNonIsotropicTermsU[1][1][1]->Fill(dEta,dUX,1.);
a554b203 538 }
7a01f4a7 539
a554b203 540 } else { //do not subtract the particle from the flowvector
541 dQmX = vQ.X()/dMq;
542 dQmY = vQ.Y()/dMq;
543 vQm.Set(dQmX,dQmY);
7a01f4a7 544
a554b203 545 //dUQ = scalar product of vU and vQm
546 Double_t dUQ = (vU * vQm);
547 Double_t dPt = pTrack->Pt();
548 Double_t dEta = pTrack->Eta();
549 //fill the profile histograms
550 if (pTrack->InRPSelection()) {
551 fHistProUQetaRP -> Fill(dEta,dUQ,dMq); //Fill (Qu/Mq) with weight Mq
552 //needed for the error calculation:
553 fHistProUQQaQbEtaRP -> Fill(dEta,dUQ*dQaQb,dMq*dMa*dMb); //Fill [Qu/Mq]*[QaQb/MaMb] with weight Mq*MaMb
554 fHistProUQPtRP -> Fill(dPt,dUQ,dMq); //Fill (Qu/Mq) with weight Mq
555 fHistProUQQaQbPtRP -> Fill(dPt,dUQ*dQaQb,dMq*dMa*dMb); //Fill [Qu/Mq]*[QaQb/MaMb] with weight Mq*MaMb
556
557 fHistSumOfWeightsEtaRP[0]->Fill(dEta,dMq); // sum of Mq
558 fHistSumOfWeightsEtaRP[1]->Fill(dEta,pow(dMq,2.));// sum of Mq^2
559 fHistSumOfWeightsEtaRP[2]->Fill(dEta,dMq*dMa*dMb);// sum of Mq*MaMb
560 fHistSumOfWeightsPtRP[0]->Fill(dPt,dMq); // sum of Mq
561 fHistSumOfWeightsPtRP[1]->Fill(dPt,pow(dMq,2.)); // sum of Mq^2
29ecccee 562 fHistSumOfWeightsPtRP[2]->Fill(dPt,dMq*dMa*dMb); // sum of Mq*MaMb
563 //nonisotropic terms:
564 fHistProNonIsotropicTermsU[0][0][0]->Fill(dPt,dUY,1.);
565 fHistProNonIsotropicTermsU[0][0][1]->Fill(dPt,dUX,1.);
566 fHistProNonIsotropicTermsU[0][1][0]->Fill(dEta,dUY,1.);
567 fHistProNonIsotropicTermsU[0][1][1]->Fill(dEta,dUX,1.);
a554b203 568 }
569 if (pTrack->InPOISelection()) {
570 fHistProUQetaPOI -> Fill(dEta,dUQ,dMq); //Fill (Qu/Mq) with weight Mq
571 //needed for the error calculation:
572 fHistProUQQaQbEtaPOI -> Fill(dEta,dUQ*dQaQb,dMq*dMa*dMb); //Fill [Qu/Mq]*[QaQb/MaMb] with weight Mq*MaMb
573 fHistProUQPtPOI -> Fill(dPt,dUQ,dMq); //Fill (Qu/Mq) with weight Mq
574 fHistProUQQaQbPtPOI -> Fill(dPt,dUQ*dQaQb,dMq*dMa*dMb); //Fill [Qu/Mq]*[QaQb/MaMb] with weight Mq*MaMb
575
576 fHistSumOfWeightsEtaPOI[0]->Fill(dEta,dMq); // sum of Mq
577 fHistSumOfWeightsEtaPOI[1]->Fill(dEta,pow(dMq,2.));// sum of Mq^2
578 fHistSumOfWeightsEtaPOI[2]->Fill(dEta,dMq*dMa*dMb);// sum of Mq*MaMb
579 fHistSumOfWeightsPtPOI[0]->Fill(dPt,dMq); // sum of Mq
580 fHistSumOfWeightsPtPOI[1]->Fill(dPt,pow(dMq,2.)); // sum of Mq^2
581 fHistSumOfWeightsPtPOI[2]->Fill(dPt,dMq*dMa*dMb); // sum of Mq*MaMb
29ecccee 582 //nonisotropic terms:
583 fHistProNonIsotropicTermsU[1][0][0]->Fill(dPt,dUY,1.);
584 fHistProNonIsotropicTermsU[1][0][1]->Fill(dPt,dUX,1.);
585 fHistProNonIsotropicTermsU[1][1][0]->Fill(dEta,dUY,1.);
586 fHistProNonIsotropicTermsU[1][1][1]->Fill(dEta,dUX,1.);
a554b203 587 }
588 }//track not in subevents
7a01f4a7 589
a554b203 590 }//track
7a01f4a7 591
a554b203 592 }//loop over tracks
593
594 fEventNumber++;
595
596 } //difference Ma and Mb
12dc476e 597
7a01f4a7 598 }// subevents not empty
b125a454 599 delete [] vQarray;
a554b203 600
7a01f4a7 601 } //event
a554b203 602
7a01f4a7 603}//end of Make()
8d312f00 604
12dc476e 605//--------------------------------------------------------------------
fd46c3dd 606void AliFlowAnalysisWithScalarProduct::GetOutputHistograms(TList *outputListHistos){
607
608 //get pointers to all output histograms (called before Finish())
12dc476e 609
fd46c3dd 610 if (outputListHistos) {
611 //Get the common histograms from the output list
612 AliFlowCommonHist *pCommonHist = dynamic_cast<AliFlowCommonHist*>
613 (outputListHistos->FindObject("AliFlowCommonHistSP"));
614 AliFlowCommonHistResults *pCommonHistResults = dynamic_cast<AliFlowCommonHistResults*>
615 (outputListHistos->FindObject("AliFlowCommonHistResultsSP"));
616 TProfile* pHistProQaQb = dynamic_cast<TProfile*>(outputListHistos->FindObject("Flow_QaQb_SP"));
a554b203 617 TProfile* pHistProQaQbNorm = dynamic_cast<TProfile*>(outputListHistos->FindObject("FlowPro_QaQbNorm_SP"));
29ecccee 618 TProfile* pHistProQaQbReImNorm = dynamic_cast<TProfile*>(outputListHistos->FindObject("FlowPro_QaQbReImNorm_SP"));
619 TProfile* pHistProNonIsotropicTermsQ = dynamic_cast<TProfile*>(outputListHistos->FindObject("FlowPro_NonIsotropicTermsQ_SP"));
a554b203 620 TH1D* pHistSumOfLinearWeights = dynamic_cast<TH1D*>(outputListHistos->FindObject("Flow_SumOfLinearWeights_SP"));
621 TH1D* pHistSumOfQuadraticWeights = dynamic_cast<TH1D*>(outputListHistos->FindObject("Flow_SumOfQuadraticWeights_SP"));
622
29ecccee 623 TProfile* pHistProFlags = dynamic_cast<TProfile*>(outputListHistos->FindObject("Flow_Flags_SP"));
fd46c3dd 624 TProfile* pHistProUQetaRP = dynamic_cast<TProfile*>(outputListHistos->FindObject("Flow_UQetaRP_SP"));
625 TProfile* pHistProUQetaPOI = dynamic_cast<TProfile*>(outputListHistos->FindObject("Flow_UQetaPOI_SP"));
626 TProfile* pHistProUQPtRP = dynamic_cast<TProfile*>(outputListHistos->FindObject("Flow_UQPtRP_SP"));
627 TProfile* pHistProUQPtPOI = dynamic_cast<TProfile*>(outputListHistos->FindObject("Flow_UQPtPOI_SP"));
12dc476e 628 TProfile* pHistProUQQaQbPtRP = dynamic_cast<TProfile*>(outputListHistos->FindObject("Flow_UQQaQbPtRP_SP"));
629 TProfile* pHistProUQQaQbEtaRP = dynamic_cast<TProfile*>(outputListHistos->FindObject("Flow_UQQaQbEtaRP_SP"));
630 TProfile* pHistProUQQaQbPtPOI = dynamic_cast<TProfile*>(outputListHistos->FindObject("Flow_UQQaQbPtPOI_SP"));
631 TProfile* pHistProUQQaQbEtaPOI = dynamic_cast<TProfile*>(outputListHistos->FindObject("Flow_UQQaQbEtaPOI_SP"));
12dc476e 632 TString weightFlag[3] = {"w_Qu_","w_Qu^2_","w_QuQaQb_"};
633 TH1D* pHistSumOfWeightsPtRP[3] = {NULL};
634 TH1D* pHistSumOfWeightsEtaRP[3] = {NULL};
635 TH1D* pHistSumOfWeightsPtPOI[3] = {NULL};
636 TH1D* pHistSumOfWeightsEtaPOI[3] = {NULL};
637 for(Int_t i=0;i<3;i++)
638 {
cff1bec4 639 pHistSumOfWeightsPtRP[i] = dynamic_cast<TH1D*>(outputListHistos->FindObject(Form("Flow_SumOfWeights%sPtRP_SP",weightFlag[i].Data())));
640 pHistSumOfWeightsEtaRP[i] = dynamic_cast<TH1D*>(outputListHistos->FindObject(Form("Flow_SumOfWeights%sEtaRP_SP",weightFlag[i].Data())));
641 pHistSumOfWeightsPtPOI[i] = dynamic_cast<TH1D*>(outputListHistos->FindObject(Form("Flow_SumOfWeights%sPtPOI_SP",weightFlag[i].Data())));
642 pHistSumOfWeightsEtaPOI[i] = dynamic_cast<TH1D*>(outputListHistos->FindObject(Form("Flow_SumOfWeights%sEtaPOI_SP",weightFlag[i].Data())));
29ecccee 643 }
644 TString rpPoi[2] = {"RP","POI"};
645 TString ptEta[2] = {"Pt","Eta"};
646 TString sinCos[2] = {"sin","cos"};
647 TProfile *pHistProNonIsotropicTermsU[2][2][2] = {{{NULL}}};
648 for(Int_t rp=0;rp<2;rp++)
649 {
650 for(Int_t pe=0;pe<2;pe++)
651 {
652 for(Int_t sc=0;sc<2;sc++)
653 {
654 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())));
655 }
656 }
657 }
658
a554b203 659 TH1D* pHistQaQb = dynamic_cast<TH1D*>(outputListHistos->FindObject("Flow_QaQb_SP"));
660 TH1D* pHistQaQbNorm = dynamic_cast<TH1D*>(outputListHistos->FindObject("Flow_QaQbNorm_SP"));
661 TH1D* pHistQaQbCos = dynamic_cast<TH1D*>(outputListHistos->FindObject("Flow_QaQbCos_SP"));
662 TH1D* pHistResolution = dynamic_cast<TH1D*>(outputListHistos->FindObject("Flow_resolution_SP"));
663 TH1D* pHistQaNorm = dynamic_cast<TH1D*>(outputListHistos->FindObject("Flow_QaNorm_SP"));
664 TH2D* pHistQaNormvsMa = dynamic_cast<TH2D*>(outputListHistos->FindObject("Flow_QaNormvsMa_SP"));
665 TH1D* pHistQbNorm = dynamic_cast<TH1D*>(outputListHistos->FindObject("Flow_QbNorm_SP"));
666 TH2D* pHistQbNormvsMb = dynamic_cast<TH2D*>(outputListHistos->FindObject("Flow_QbNormvsMb_SP"));
667 TH2D* pHistMavsMb = dynamic_cast<TH2D*>(outputListHistos->FindObject("Flow_MavsMb_SP"));
668
12dc476e 669 //pass the pointers to the task
a554b203 670 if (pCommonHist && pCommonHistResults && pHistProQaQb && pHistProQaQbNorm &&
29ecccee 671 pHistProQaQbReImNorm && pHistProNonIsotropicTermsQ &&
a554b203 672 pHistSumOfLinearWeights && pHistSumOfQuadraticWeights &&
673 pHistQaQb && pHistQaQbNorm && pHistQaQbCos && pHistResolution &&
674 pHistQaNorm && pHistQaNormvsMa && pHistQbNorm && pHistQbNormvsMb &&
29ecccee 675 pHistMavsMb && pHistProFlags &&
a554b203 676 pHistProUQetaRP && pHistProUQetaPOI &&
677 pHistProUQPtRP && pHistProUQPtPOI &&
678 pHistProUQQaQbPtRP && pHistProUQQaQbEtaRP &&
679 pHistProUQQaQbPtPOI && pHistProUQQaQbEtaPOI) {
fd46c3dd 680 this -> SetCommonHists(pCommonHist);
681 this -> SetCommonHistsRes(pCommonHistResults);
682 this -> SetHistProQaQb(pHistProQaQb);
a554b203 683 this -> SetHistProQaQbNorm(pHistProQaQbNorm);
29ecccee 684 this -> SetHistProQaQbReImNorm(pHistProQaQbReImNorm);
685 this -> SetHistProNonIsotropicTermsQ(pHistProNonIsotropicTermsQ);
a554b203 686 this -> SetHistSumOfLinearWeights(pHistSumOfLinearWeights);
687 this -> SetHistSumOfQuadraticWeights(pHistSumOfQuadraticWeights);
688 this -> SetHistQaQb(pHistQaQb);
689 this -> SetHistQaQbNorm(pHistQaQbNorm);
690 this -> SetHistQaQbCos(pHistQaQbCos);
691 this -> SetHistResolution(pHistResolution);
692 this -> SetHistQaNorm(pHistQaNorm);
693 this -> SetHistQaNormvsMa(pHistQaNormvsMa);
694 this -> SetHistQbNorm(pHistQbNorm);
695 this -> SetHistQbNormvsMb(pHistQbNormvsMb);
696 this -> SetHistMavsMb(pHistMavsMb);
29ecccee 697 this -> SetHistProFlags(pHistProFlags);
fd46c3dd 698 this -> SetHistProUQetaRP(pHistProUQetaRP);
699 this -> SetHistProUQetaPOI(pHistProUQetaPOI);
700 this -> SetHistProUQPtRP(pHistProUQPtRP);
701 this -> SetHistProUQPtPOI(pHistProUQPtPOI);
12dc476e 702 this -> SetHistProUQQaQbPtRP(pHistProUQQaQbPtRP);
703 this -> SetHistProUQQaQbEtaRP(pHistProUQQaQbEtaRP);
704 this -> SetHistProUQQaQbPtPOI(pHistProUQQaQbPtPOI);
705 this -> SetHistProUQQaQbEtaPOI(pHistProUQQaQbEtaPOI);
a554b203 706
707 for(Int_t i=0;i<3;i++) {
708 if(pHistSumOfWeightsPtRP[i]) this -> SetHistSumOfWeightsPtRP(pHistSumOfWeightsPtRP[i],i);
709 if(pHistSumOfWeightsEtaRP[i]) this -> SetHistSumOfWeightsEtaRP(pHistSumOfWeightsEtaRP[i],i);
710 if(pHistSumOfWeightsPtPOI[i]) this -> SetHistSumOfWeightsPtPOI(pHistSumOfWeightsPtPOI[i],i);
711 if(pHistSumOfWeightsEtaPOI[i]) this -> SetHistSumOfWeightsEtaPOI(pHistSumOfWeightsEtaPOI[i],i);
29ecccee 712 }
713 for(Int_t rp=0;rp<2;rp++)
714 {
715 for(Int_t pe=0;pe<2;pe++)
716 {
717 for(Int_t sc=0;sc<2;sc++)
718 {
719 if(pHistProNonIsotropicTermsU[rp][pe][sc]) this->SetHistProNonIsotropicTermsU(pHistProNonIsotropicTermsU[rp][pe][sc],rp,pe,sc);
720 }
721 }
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);
29ecccee 761 //non-isotropic terms:
762 Double_t dImQa = fHistProQaQbReImNorm->GetBinContent(1);
763 Double_t dReQa = fHistProQaQbReImNorm->GetBinContent(2);
764 Double_t dImQb = fHistProQaQbReImNorm->GetBinContent(3);
765 Double_t dReQb = fHistProQaQbReImNorm->GetBinContent(4);
3c5d5752 766
29ecccee 767 if(fApplyCorrectionForNUA)
768 {
769 dQaQb = dQaQb - dImQa*dImQb - dReQa*dReQb;
770 }
771
a554b203 772 Double_t dV = -999.;
524798bf 773 if(dQaQb>=0.)
774 {
775 dV = TMath::Sqrt(dQaQb);
776 }
12dc476e 777 //statistical error of dQaQb:
778 // statistical error = term1 * spread * term2:
779 // term1 = sqrt{sum_{i=1}^{N} w^2}/(sum_{i=1}^{N} w)
780 // term2 = 1/sqrt(1-term1^2)
a554b203 781 Double_t dSpreadQaQb = fHistProQaQbNorm->GetBinError(1);
12dc476e 782 Double_t dSumOfLinearWeights = fHistSumOfLinearWeights->GetBinContent(1);
783 Double_t dSumOfQuadraticWeights = fHistSumOfQuadraticWeights->GetBinContent(1);
784 Double_t dTerm1 = 0.;
785 Double_t dTerm2 = 0.;
786 if(dSumOfLinearWeights) {
787 dTerm1 = pow(dSumOfQuadraticWeights,0.5)/dSumOfLinearWeights;
788 }
789 if(1.-pow(dTerm1,2.) > 0.) {
790 dTerm2 = 1./pow(1-pow(dTerm1,2.),0.5);
791 }
792 Double_t dStatErrorQaQb = dTerm1 * dSpreadQaQb * dTerm2;
793 //calculate the statistical error
794 Double_t dVerr = 0.;
795 if(dQaQb > 0.) {
796 dVerr = (1./(2.*pow(dQaQb,0.5)))*dStatErrorQaQb;
797 }
798 fCommonHistsRes->FillIntegratedFlow(dV,dVerr);
799 cout<<"dV = "<<dV<<" +- "<<dVerr<<endl;
800
801 //Calculate differential flow and integrated flow (RP, POI)
802 //---------------------------------------------------------
12dc476e 803 //v as a function of eta for RP selection
a9de3704 804 for(Int_t b=1;b<iNbinsEta+1;b++) {
12dc476e 805 Double_t duQpro = fHistProUQetaRP->GetBinContent(b);
29ecccee 806 if(fApplyCorrectionForNUA)
807 {
808 duQpro = duQpro
809 - fHistProNonIsotropicTermsU[0][1][1]->GetBinContent(b)*fHistProNonIsotropicTermsQ->GetBinContent(2)
810 - fHistProNonIsotropicTermsU[0][1][0]->GetBinContent(b)*fHistProNonIsotropicTermsQ->GetBinContent(1);
811 }
a554b203 812 Double_t dv2pro = -999.;
12dc476e 813 if (dV!=0.) { dv2pro = duQpro/dV; }
814 //calculate the statistical error
815 Double_t dv2ProErr = CalculateStatisticalError(b, dStatErrorQaQb, fHistProUQetaRP, fHistProUQQaQbEtaRP, fHistSumOfWeightsEtaRP);
12dc476e 816 //fill TH1D
cff1bec4 817 fCommonHistsRes->FillDifferentialFlowEtaRP(b, dv2pro, dv2ProErr);
12dc476e 818 } //loop over bins b
29ecccee 819
80f72ed6 820
12dc476e 821 //v as a function of eta for POI selection
a9de3704 822 for(Int_t b=1;b<iNbinsEta+1;b++) {
12dc476e 823 Double_t duQpro = fHistProUQetaPOI->GetBinContent(b);
29ecccee 824 if(fApplyCorrectionForNUA)
825 {
826 duQpro = duQpro
827 - fHistProNonIsotropicTermsU[1][1][1]->GetBinContent(b)*fHistProNonIsotropicTermsQ->GetBinContent(2)
828 - fHistProNonIsotropicTermsU[1][1][0]->GetBinContent(b)*fHistProNonIsotropicTermsQ->GetBinContent(1);
829 }
a554b203 830 Double_t dv2pro = -999.;
12dc476e 831 if (dV!=0.) { dv2pro = duQpro/dV; }
832 //calculate the statistical error
833 Double_t dv2ProErr = CalculateStatisticalError(b, dStatErrorQaQb, fHistProUQetaPOI, fHistProUQQaQbEtaPOI, fHistSumOfWeightsEtaPOI);
834
835 //fill TH1D
836 fCommonHistsRes->FillDifferentialFlowEtaPOI(b, dv2pro, dv2ProErr);
837 } //loop over bins b
838
80f72ed6 839
12dc476e 840 //v as a function of Pt for RP selection
841 TH1F* fHistPtRP = fCommonHists->GetHistPtRP(); //for calculating integrated flow
842 Double_t dVRP = 0.;
843 Double_t dSumRP = 0.;
844 Double_t dErrVRP =0.;
845
a9de3704 846 for(Int_t b=1;b<iNbinsPt+1;b++) {
12dc476e 847 Double_t duQpro = fHistProUQPtRP->GetBinContent(b);
29ecccee 848 if(fApplyCorrectionForNUA)
849 {
850 duQpro = duQpro
851 - fHistProNonIsotropicTermsU[0][0][1]->GetBinContent(b)*fHistProNonIsotropicTermsQ->GetBinContent(2)
852 - fHistProNonIsotropicTermsU[0][0][0]->GetBinContent(b)*fHistProNonIsotropicTermsQ->GetBinContent(1);
853 }
a554b203 854 Double_t dv2pro = -999.;
12dc476e 855 if (dV!=0.) { dv2pro = duQpro/dV; }
856 //calculate the statistical error
857 Double_t dv2ProErr = CalculateStatisticalError(b, dStatErrorQaQb, fHistProUQPtRP, fHistProUQQaQbPtRP, fHistSumOfWeightsPtRP);
858
859 //fill TH1D
860 fCommonHistsRes->FillDifferentialFlowPtRP(b, dv2pro, dv2ProErr);
861
862 //calculate integrated flow for RP selection
863 if (fHistPtRP){
864 Double_t dYieldPt = fHistPtRP->GetBinContent(b);
865 dVRP += dv2pro*dYieldPt;
866 dSumRP +=dYieldPt;
867 dErrVRP += dYieldPt*dYieldPt*dv2ProErr*dv2ProErr;
868 } else { cout<<"fHistPtRP is NULL"<<endl; }
869 } //loop over bins b
870
871 if (dSumRP != 0.) {
872 dVRP /= dSumRP; //the pt distribution should be normalised
873 dErrVRP /= (dSumRP*dSumRP);
874 dErrVRP = TMath::Sqrt(dErrVRP);
80f72ed6 875 }
12dc476e 876 fCommonHistsRes->FillIntegratedFlowRP(dVRP,dErrVRP);
877 cout<<"dV(RP) = "<<dVRP<<" +- "<<dErrVRP<<endl;
878
80f72ed6 879
12dc476e 880 //v as a function of Pt for POI selection
881 TH1F* fHistPtPOI = fCommonHists->GetHistPtPOI(); //for calculating integrated flow
882 Double_t dVPOI = 0.;
883 Double_t dSumPOI = 0.;
884 Double_t dErrVPOI =0.;
885
a9de3704 886 for(Int_t b=1;b<iNbinsPt+1;b++) {
12dc476e 887 Double_t duQpro = fHistProUQPtPOI->GetBinContent(b);
29ecccee 888 if(fApplyCorrectionForNUA)
889 {
890 duQpro = duQpro
891 - fHistProNonIsotropicTermsU[1][0][1]->GetBinContent(b)*fHistProNonIsotropicTermsQ->GetBinContent(2)
892 - fHistProNonIsotropicTermsU[1][0][0]->GetBinContent(b)*fHistProNonIsotropicTermsQ->GetBinContent(1);
893 }
a554b203 894 Double_t dv2pro = -999.;
12dc476e 895 if (dV!=0.) { dv2pro = duQpro/dV; }
896 //calculate the statistical error
897 Double_t dv2ProErr = CalculateStatisticalError(b, dStatErrorQaQb, fHistProUQPtPOI, fHistProUQQaQbPtPOI, fHistSumOfWeightsPtPOI);
898
899 //fill TH1D
900 fCommonHistsRes->FillDifferentialFlowPtPOI(b, dv2pro, dv2ProErr);
395fadba 901
12dc476e 902 //calculate integrated flow for POI selection
903 if (fHistPtPOI){
904 Double_t dYieldPt = fHistPtPOI->GetBinContent(b);
905 dVPOI += dv2pro*dYieldPt;
906 dSumPOI +=dYieldPt;
907 dErrVPOI += dYieldPt*dYieldPt*dv2ProErr*dv2ProErr;
908 } else { cout<<"fHistPtPOI is NULL"<<endl; }
909 } //loop over bins b
395fadba 910
12dc476e 911 if (dSumPOI != 0.) {
489d5531 912 dVPOI /= dSumPOI; //the pt distribution should be normalised
12dc476e 913 dErrVPOI /= (dSumPOI*dSumPOI);
914 dErrVPOI = TMath::Sqrt(dErrVPOI);
915 }
916 fCommonHistsRes->FillIntegratedFlowPOI(dVPOI,dErrVPOI);
917 cout<<"dV(POI) = "<<dVPOI<<" +- "<<dErrVPOI<<endl;
e4cecfa0 918
395fadba 919 cout<<endl;
920 cout<<"*************************************"<<endl;
921 cout<<"*************************************"<<endl;
12dc476e 922
395fadba 923 //cout<<".....finished"<<endl;
12dc476e 924}
8232a5ec 925
926
12dc476e 927//--------------------------------------------------------------------
928Double_t AliFlowAnalysisWithScalarProduct::CalculateStatisticalError(Int_t b, Double_t aStatErrorQaQb, TProfile* pHistProUQ, TProfile* pHistProUQQaQb, TH1D** pHistSumOfWeights) {
12dc476e 929 //calculate the statistical error for differential flow for bin b
930 Double_t duQproSpread = pHistProUQ->GetBinError(b);
931 Double_t sumOfMq = pHistSumOfWeights[0]->GetBinContent(b);
932 Double_t sumOfMqSquared = pHistSumOfWeights[1]->GetBinContent(b);
a554b203 933 Double_t dQaQb = fHistProQaQbNorm->GetBinContent(1);
12dc476e 934 Double_t dTerm1 = 0.;
935 Double_t dTerm2 = 0.;
936 if(sumOfMq) {
937 dTerm1 = (pow(sumOfMqSquared,0.5)/sumOfMq);
938 }
939 if(1.-pow(dTerm1,2.)>0.) {
940 dTerm2 = 1./pow(1.-pow(dTerm1,2.),0.5);
941 }
942 Double_t duQproErr = dTerm1*duQproSpread*dTerm2;
12dc476e 943 // covariances:
944 Double_t dTerm1Cov = pHistSumOfWeights[2]->GetBinContent(b);
945 Double_t dTerm2Cov = fHistSumOfLinearWeights->GetBinContent(1);
946 Double_t dTerm3Cov = sumOfMq;
947 Double_t dWeightedCovariance = 0.;
948 if(dTerm2Cov*dTerm3Cov>0.) {
949 Double_t dDenominator = 1.-dTerm1Cov/(dTerm2Cov*dTerm3Cov);
950 Double_t dPrefactor = dTerm1Cov/(dTerm2Cov*dTerm3Cov);
951 if(dDenominator!=0) {
a554b203 952 Double_t dCovariance = (pHistProUQQaQb->GetBinContent(b)-dQaQb*pHistProUQ->GetBinContent(b))/dDenominator;
12dc476e 953 dWeightedCovariance = dCovariance*dPrefactor;
954 }
955 }
cff1bec4 956
12dc476e 957 Double_t dv2ProErr = 0.; // final statitical error:
a554b203 958 if(dQaQb>0.) {
959 Double_t dv2ProErrorSquared = (1./4.)*pow(dQaQb,-3.)*
12dc476e 960 (pow(pHistProUQ->GetBinContent(b),2.)*pow(aStatErrorQaQb,2.)
a554b203 961 + 4.*pow(dQaQb,2.)*pow(duQproErr,2.)
962 - 4.*dQaQb*pHistProUQ->GetBinContent(b)*dWeightedCovariance);
12dc476e 963 if(dv2ProErrorSquared>0.) dv2ProErr = pow(dv2ProErrorSquared,0.5);
964 }
965
966 return dv2ProErr;
967}
29ecccee 968
969
970//--------------------------------------------------------------------
971
972void AliFlowAnalysisWithScalarProduct::StoreFlags()
973{
974 // Store all boolean flags needed in Finish() in profile fHistProFlags.
975
976 // Apply correction for non-uniform acceptance or not:
977 fHistProFlags->Fill(0.5,fApplyCorrectionForNUA);
978
979}
980
981//--------------------------------------------------------------------
982
983void AliFlowAnalysisWithScalarProduct::AccessFlags()
984{
985 // Access all boolean flags needed in Finish() from profile fHistProFlags.
986
987 // Apply correction for non-uniform acceptance or not:
988 fApplyCorrectionForNUA = (Bool_t) fHistProFlags->GetBinContent(1);
989
990}
991
12dc476e 992//--------------------------------------------------------------------