]> 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
fix to prevent sqrt of a negative number
[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),
52 fPhiWeights(NULL),
e35ddff0 53 fHistList(NULL),
395fadba 54 fHistProUQetaRP(NULL),
55 fHistProUQetaPOI(NULL),
56 fHistProUQPtRP(NULL),
57 fHistProUQPtPOI(NULL),
e4cecfa0 58 fHistProQaQb(NULL),
a554b203 59 fHistProQaQbNorm(NULL),
12dc476e 60 fHistSumOfLinearWeights(NULL),
61 fHistSumOfQuadraticWeights(NULL),
62 fHistProUQQaQbPtRP(NULL),
63 fHistProUQQaQbEtaRP(NULL),
64 fHistProUQQaQbPtPOI(NULL),
65 fHistProUQQaQbEtaPOI(NULL),
395fadba 66 fCommonHists(NULL),
a554b203 67 fCommonHistsRes(NULL),
68 fHistQaQb(NULL),
69 fHistQaQbNorm(NULL),
70 fHistQaQbCos(NULL),
71 fHistResolution(NULL),
72 fHistQaNorm(NULL),
73 fHistQaNormvsMa(NULL),
74 fHistQbNorm(NULL),
75 fHistQbNormvsMb(NULL),
76 fHistMavsMb(NULL)
77
8d312f00 78{
8d312f00 79 // Constructor.
29195b69 80 fWeightsList = new TList();
e35ddff0 81 fHistList = new TList();
12dc476e 82
83 // Initialize arrays:
84 for(Int_t i=0;i<3;i++)
85 {
86 fHistSumOfWeightsPtRP[i] = NULL;
87 fHistSumOfWeightsEtaRP[i] = NULL;
88 fHistSumOfWeightsPtPOI[i] = NULL;
89 fHistSumOfWeightsEtaPOI[i] = NULL;
90 }
8d312f00 91}
92 //-----------------------------------------------------------------------
8d312f00 93 AliFlowAnalysisWithScalarProduct::~AliFlowAnalysisWithScalarProduct()
94 {
95 //destructor
29195b69 96 delete fWeightsList;
e35ddff0 97 delete fHistList;
8d312f00 98 }
99
1dfa3c16 100//-----------------------------------------------------------------------
1dfa3c16 101void AliFlowAnalysisWithScalarProduct::WriteHistograms(TString* outputFileName)
102{
103 //store the final results in output .root file
104
105 TFile *output = new TFile(outputFileName->Data(),"RECREATE");
0fe80f88 106 //output->WriteObject(fHistList, "cobjSP","SingleKey");
107 fHistList->SetName("cobjSP");
9455e15e 108 fHistList->SetOwner(kTRUE);
0fe80f88 109 fHistList->Write(fHistList->GetName(), TObject::kSingleKey);
1dfa3c16 110 delete output;
111}
112
b0fda271 113//-----------------------------------------------------------------------
b0fda271 114void AliFlowAnalysisWithScalarProduct::WriteHistograms(TString outputFileName)
115{
116 //store the final results in output .root file
117
118 TFile *output = new TFile(outputFileName.Data(),"RECREATE");
0fe80f88 119 //output->WriteObject(fHistList, "cobjSP","SingleKey");
120 fHistList->SetName("cobjSP");
9455e15e 121 fHistList->SetOwner(kTRUE);
0fe80f88 122 fHistList->Write(fHistList->GetName(), TObject::kSingleKey);
b0fda271 123 delete output;
124}
125
ad87ae62 126//-----------------------------------------------------------------------
ad87ae62 127void AliFlowAnalysisWithScalarProduct::WriteHistograms(TDirectoryFile *outputFileName)
128{
129 //store the final results in output .root file
130 fHistList->SetName("cobjSP");
131 fHistList->SetOwner(kTRUE);
132 outputFileName->Add(fHistList);
133 outputFileName->Write(outputFileName->GetName(), TObject::kSingleKey);
134}
135
8d312f00 136//-----------------------------------------------------------------------
137void AliFlowAnalysisWithScalarProduct::Init() {
138
139 //Define all histograms
d7eb18ec 140 cout<<"---Analysis with the Scalar Product Method--- Init"<<endl;
8d312f00 141
489d5531 142 //save old value and prevent histograms from being added to directory
143 //to avoid name clashes in case multiple analaysis objects are used
144 //in an analysis
145 Bool_t oldHistAddStatus = TH1::AddDirectoryStatus();
146 TH1::AddDirectory(kFALSE);
147
bee2efdc 148 Int_t iNbinsPt = AliFlowCommonConstants::GetMaster()->GetNbinsPt();
149 Double_t dPtMin = AliFlowCommonConstants::GetMaster()->GetPtMin();
150 Double_t dPtMax = AliFlowCommonConstants::GetMaster()->GetPtMax();
151 Int_t iNbinsEta = AliFlowCommonConstants::GetMaster()->GetNbinsEta();
152 Double_t dEtaMin = AliFlowCommonConstants::GetMaster()->GetEtaMin();
153 Double_t dEtaMax = AliFlowCommonConstants::GetMaster()->GetEtaMax();
8d312f00 154
12dc476e 155 fHistProUQetaRP = new TProfile("Flow_UQetaRP_SP","Flow_UQetaRP_SP",iNbinsEta,dEtaMin,dEtaMax,"s");
395fadba 156 fHistProUQetaRP->SetXTitle("{eta}");
157 fHistProUQetaRP->SetYTitle("<uQ>");
158 fHistList->Add(fHistProUQetaRP);
159
12dc476e 160 fHistProUQetaPOI = new TProfile("Flow_UQetaPOI_SP","Flow_UQetaPOI_SP",iNbinsEta,dEtaMin,dEtaMax,"s");
395fadba 161 fHistProUQetaPOI->SetXTitle("{eta}");
162 fHistProUQetaPOI->SetYTitle("<uQ>");
163 fHistList->Add(fHistProUQetaPOI);
164
12dc476e 165 fHistProUQPtRP = new TProfile("Flow_UQPtRP_SP","Flow_UQPtRP_SP",iNbinsPt,dPtMin,dPtMax,"s");
395fadba 166 fHistProUQPtRP->SetXTitle("p_t (GeV)");
167 fHistProUQPtRP->SetYTitle("<uQ>");
168 fHistList->Add(fHistProUQPtRP);
169
12dc476e 170 fHistProUQPtPOI = new TProfile("Flow_UQPtPOI_SP","Flow_UQPtPOI_SP",iNbinsPt,dPtMin,dPtMax,"s");
395fadba 171 fHistProUQPtPOI->SetXTitle("p_t (GeV)");
172 fHistProUQPtPOI->SetYTitle("<uQ>");
173 fHistList->Add(fHistProUQPtPOI);
174
12dc476e 175 fHistProQaQb = new TProfile("Flow_QaQb_SP","Flow_QaQb_SP", 1, 0.5, 1.5,"s");
e4cecfa0 176 fHistProQaQb->SetYTitle("<QaQb>");
177 fHistList->Add(fHistProQaQb);
8d312f00 178
a554b203 179 fHistProQaQbNorm = new TProfile("FlowPro_QaQbNorm_SP","FlowPro_QaQbNorm_SP", 1, 0.5, 1.5,"s");
180 fHistProQaQbNorm->SetYTitle("<QaQb/MaMb>");
181 fHistList->Add(fHistProQaQbNorm);
182
12dc476e 183 fHistSumOfLinearWeights = new TH1D("Flow_SumOfLinearWeights_SP","Flow_SumOfLinearWeights_SP",1,-0.5, 0.5);
184 fHistSumOfLinearWeights -> SetYTitle("sum (*)");
185 fHistSumOfLinearWeights -> SetXTitle("sum (Ma*Mb)");
186 fHistList->Add(fHistSumOfLinearWeights);
187
188 fHistSumOfQuadraticWeights = new TH1D("Flow_SumOfQuadraticWeights_SP","Flow_SumOfQuadraticWeights_SP",1,-0.5, 0.5);
189 fHistSumOfQuadraticWeights -> SetYTitle("sum (*)");
190 fHistSumOfQuadraticWeights -> SetXTitle("sum (Ma*Mb)^2");
191 fHistList->Add(fHistSumOfQuadraticWeights);
192
193 fHistProUQQaQbPtRP = new TProfile("Flow_UQQaQbPtRP_SP","Flow_UQQaQbPtRP_SP",iNbinsPt,dPtMin,dPtMax);
194 fHistProUQQaQbPtRP -> SetYTitle("<*>");
195 fHistProUQQaQbPtRP -> SetXTitle("<Qu QaQb>");
196 fHistList->Add(fHistProUQQaQbPtRP);
197
198 fHistProUQQaQbEtaRP = new TProfile("Flow_UQQaQbEtaRP_SP","Flow_UQQaQbEtaRP_SP",iNbinsEta,dEtaMin,dEtaMax);
199 fHistProUQQaQbEtaRP -> SetYTitle("<*>");
200 fHistProUQQaQbEtaRP -> SetXTitle("<Qu QaQb>");
201 fHistList->Add(fHistProUQQaQbEtaRP);
202
203 fHistProUQQaQbPtPOI = new TProfile("Flow_UQQaQbPtPOI_SP","Flow_UQQaQbPtPOI_SP",iNbinsPt,dPtMin,dPtMax);
204 fHistProUQQaQbPtPOI -> SetYTitle("<*>");
205 fHistProUQQaQbPtPOI -> SetXTitle("<Qu QaQb>");
206 fHistList->Add(fHistProUQQaQbPtPOI);
207
208 fHistProUQQaQbEtaPOI = new TProfile("Flow_UQQaQbEtaPOI_SP","Flow_UQQaQbEtaPOI_SP",iNbinsEta,dEtaMin,dEtaMax);
209 fHistProUQQaQbEtaPOI -> SetYTitle("<*>");
210 fHistProUQQaQbEtaPOI -> SetXTitle("<Qu QaQb>");
211 fHistList->Add(fHistProUQQaQbEtaPOI);
212
213 TString weightFlag[3] = {"w_Qu_","w_Qu^2_","w_QuQaQb_"};
214 for(Int_t i=0;i<3;i++)
215 {
cff1bec4 216 fHistSumOfWeightsPtRP[i] = new TH1D(Form("Flow_SumOfWeights%sPtRP_SP",weightFlag[i].Data()),
217 Form("Flow_SumOfWeights%sPtRP_SP",weightFlag[i].Data()),iNbinsPt,dPtMin,dPtMax);
12dc476e 218 fHistSumOfWeightsPtRP[i] -> SetYTitle("sum (*)");
219 fHistSumOfWeightsPtRP[i] -> SetXTitle("p_{T}");
220 fHistList->Add(fHistSumOfWeightsPtRP[i]);
221
cff1bec4 222 fHistSumOfWeightsEtaRP[i] = new TH1D(Form("Flow_SumOfWeights%sEtaRP_SP",weightFlag[i].Data()),
223 Form("Flow_SumOfWeights%sEtaRP_SP",weightFlag[i].Data()),iNbinsEta,dEtaMin,dEtaMax);
12dc476e 224 fHistSumOfWeightsEtaRP[i] -> SetYTitle("sum (*)");
225 fHistSumOfWeightsEtaRP[i] -> SetXTitle("#eta");
226 fHistList->Add(fHistSumOfWeightsEtaRP[i]);
227
cff1bec4 228 fHistSumOfWeightsPtPOI[i] = new TH1D(Form("Flow_SumOfWeights%sPtPOI_SP",weightFlag[i].Data()),
229 Form("Flow_SumOfWeights%sPtPOI_SP",weightFlag[i].Data()),iNbinsPt,dPtMin,dPtMax);
12dc476e 230 fHistSumOfWeightsPtPOI[i] -> SetYTitle("sum (*)");
231 fHistSumOfWeightsPtPOI[i] -> SetXTitle("p_{T}");
232 fHistList->Add(fHistSumOfWeightsPtPOI[i]);
233
cff1bec4 234 fHistSumOfWeightsEtaPOI[i] = new TH1D(Form("Flow_SumOfWeights%sEtaPOI_SP",weightFlag[i].Data()),
235 Form("Flow_SumOfWeights%sEtaPOI_SP",weightFlag[i].Data()),iNbinsEta,dEtaMin,dEtaMax);
12dc476e 236 fHistSumOfWeightsEtaPOI[i] -> SetYTitle("sum (*)");
237 fHistSumOfWeightsEtaPOI[i] -> SetXTitle("#eta");
238 fHistList->Add(fHistSumOfWeightsEtaPOI[i]);
239 }
240
04f6283b 241 fCommonHists = new AliFlowCommonHist("AliFlowCommonHistSP");
099e1753 242 fHistList->Add(fCommonHists);
395fadba 243 fCommonHistsRes = new AliFlowCommonHistResults("AliFlowCommonHistResultsSP");
244 fHistList->Add(fCommonHistsRes);
245
a554b203 246 fHistQaQb = new TH1D("Flow_QaQb_SP","Flow_QaQb_SP",200,-100.,100.);
247 fHistQaQb -> SetYTitle("dN/dQaQb");
248 fHistQaQb -> SetXTitle("QaQb");
249 fHistList->Add(fHistQaQb);
250
251 fHistQaQbNorm = new TH1D("Flow_QaQbNorm_SP","Flow_QaQbNorm_SP",44,-1.1,1.1);
252 fHistQaQbNorm -> SetYTitle("dN/d(QaQb/MaMb)");
253 fHistQaQbNorm -> SetXTitle("QaQb/MaMb");
254 fHistList->Add(fHistQaQbNorm);
255
256 fHistQaQbCos = new TH1D("Flow_QaQbCos_SP","Flow_QaQbCos_SP",63,0.,TMath::Pi());
257 fHistQaQbCos -> SetYTitle("dN/d(#phi)");
258 fHistQaQbCos -> SetXTitle("#phi");
259 fHistList->Add(fHistQaQbCos);
260
261 fHistResolution = new TH1D("Flow_resolution_SP","Flow_resolution_SP",100,-1.0,1.0);
262 fHistResolution -> SetYTitle("dN/d(cos(2(#phi_a - #phi_b))");
263 fHistResolution -> SetXTitle("cos(2*(#phi_a - #phi_b))");
264 fHistList->Add(fHistResolution);
265
266 fHistQaNorm = new TH1D("Flow_QaNorm_SP","Flow_QaNorm_SP",22,0.,1.1);
267 fHistQaNorm -> SetYTitle("dN/d(|Qa/Ma|)");
268 fHistQaNorm -> SetXTitle("|Qa/Ma|");
269 fHistList->Add(fHistQaNorm);
270
271 fHistQaNormvsMa = new TH2D("Flow_QaNormvsMa_SP","Flow_QaNormvsMa_SP",100,0.,100.,22,0.,1.1);
272 fHistQaNormvsMa -> SetYTitle("|Qa/Ma|");
273 fHistQaNormvsMa -> SetXTitle("Ma");
274 fHistList->Add(fHistQaNormvsMa);
275
276 fHistQbNorm = new TH1D("Flow_QbNorm_SP","Flow_QbNorm_SP",22,0.,1.1);
277 fHistQbNorm -> SetYTitle("dN/d(|Qb/Mb|)");
278 fHistQbNorm -> SetXTitle("|Qb/Mb|");
279 fHistList->Add(fHistQbNorm);
280
281 fHistQbNormvsMb = new TH2D("Flow_QbNormvsMb_SP","Flow_QbNormvsMb_SP",100,0.,100.,22,0.,1.1);
282 fHistQbNormvsMb -> SetYTitle("|Qb/Mb|");
283 fHistQbNormvsMb -> SetXTitle("|Mb|");
284 fHistList->Add(fHistQbNormvsMb);
285
286 fHistMavsMb = new TH2D("Flow_MavsMb_SP","Flow_MavsMb_SP",100,0.,100.,100,0.,100.);
287 fHistMavsMb -> SetYTitle("Ma");
288 fHistMavsMb -> SetXTitle("Mb");
289 fHistList->Add(fHistMavsMb);
290
291
29195b69 292 //weights
293 if(fUsePhiWeights) {
294 if(!fWeightsList) {
295 cout<<"WARNING: fWeightsList is NULL in the Scalar Product method."<<endl;
296 exit(0);
297 }
298 if(fWeightsList->FindObject("phi_weights")) {
299 fPhiWeights = dynamic_cast<TH1F*>(fWeightsList->FindObject("phi_weights"));
300 fHistList->Add(fPhiWeights);
301 } else {
302 cout<<"WARNING: histogram with phi weights is not accessible in Scalar Product"<<endl;
303 exit(0);
304 }
305 } // end of if(fUsePhiWeights)
306
e2d51347 307 fEventNumber = 0; //set number of events to zero
489d5531 308
309 TH1::AddDirectory(oldHistAddStatus);
e2d51347 310}
311
8d312f00 312//-----------------------------------------------------------------------
8232a5ec 313void AliFlowAnalysisWithScalarProduct::Make(AliFlowEventSimple* anEvent) {
8d312f00 314
489d5531 315 //Calculate flow based on <QaQb/MaMb> = <v^2>
12dc476e 316
317 //Fill histograms
8232a5ec 318 if (anEvent) {
8d312f00 319
cbbaf54a 320 //get Q vectors for the eta-subevents
b125a454 321 AliFlowVector* vQarray = new AliFlowVector[2];
29195b69 322 if (fUsePhiWeights) {
323 anEvent->Get2Qsub(vQarray,2,fWeightsList,kTRUE);
324 } else {
325 anEvent->Get2Qsub(vQarray);
326 }
12dc476e 327 //subevent a
b125a454 328 AliFlowVector vQa = vQarray[0];
12dc476e 329 //subevent b
b125a454 330 AliFlowVector vQb = vQarray[1];
12dc476e 331
7a01f4a7 332 //check that the subevents are not empty:
12dc476e 333 Double_t dMa = vQa.GetMult();
12dc476e 334 Double_t dMb = vQb.GetMult();
7a01f4a7 335 if (dMa != 0 && dMb != 0) {
336
a554b203 337
338 //request that the subevent multiplicities are not too different
339 Double_t dRelDiff = TMath::Abs((dMa - dMb)/(dMa + dMb));
340 if (dRelDiff < fRelDiffMsub) {
341
342 //fill control histograms
343 fCommonHists->FillControlHistograms(anEvent);
344
345 //fill some SP control histograms
346 fHistProQaQb -> Fill(1.,vQa*vQb,1.); //Fill with weight 1 -> Weight with MaMb????
347 fHistQaQbCos ->Fill(0.5*TMath::ACos((vQa/vQa.Mod())*(vQb/vQb.Mod()))); //Acos(Qa*Qb) = 2*phi
348 fHistResolution -> Fill(TMath::Cos( vQa.Phi()- vQb.Phi() )); //vQa.Phi() returns 2*phi
349 fHistQaQb -> Fill(vQa*vQb);
350 fHistMavsMb -> Fill(dMb,dMa);
351
352 //get total Q vector = the sum of subevent a and subevent b
353 AliFlowVector vQ = vQa + vQb;
354 //weight the Q vectors for the subevents by the multiplicity
355 //Note: Weight Q only in the particle loop when it is clear if it should be (m-1) or M
356 Double_t dQXa = vQa.X()/dMa;
357 Double_t dQYa = vQa.Y()/dMa;
358 vQa.Set(dQXa,dQYa);
359
360 Double_t dQXb = vQb.X()/dMb;
361 Double_t dQYb = vQb.Y()/dMb;
362 vQb.Set(dQXb,dQYb);
12dc476e 363
a554b203 364 //scalar product of the two subevents
365 Double_t dQaQb = (vQa*vQb);
366 fHistProQaQbNorm -> Fill(1.,dQaQb,dMa*dMb); //Fill (QaQb/MaMb) with weight (MaMb).
367 //needed for the error calculation:
368 fHistSumOfLinearWeights -> Fill(0.,dMa*dMb);
369 fHistSumOfQuadraticWeights -> Fill(0.,pow(dMa*dMb,2.));
370
371 //fill some SP control histograms
372 fHistQaQbNorm ->Fill(vQa*vQb);
373 fHistQaNorm ->Fill(vQa.Mod());
374 fHistQaNormvsMa->Fill(dMa,vQa.Mod());
375 fHistQbNorm ->Fill(vQb.Mod());
376 fHistQbNormvsMb->Fill(dMb,vQb.Mod());
377
378 //loop over the tracks of the event
379 AliFlowTrackSimple* pTrack = NULL;
380 Int_t iNumberOfTracks = anEvent->NumberOfTracks();
381 Double_t dMq = vQ.GetMult();
7a01f4a7 382
a554b203 383 for (Int_t i=0;i<iNumberOfTracks;i++)
384 {
385 pTrack = anEvent->GetTrack(i) ;
386 if (pTrack){
387 Double_t dPhi = pTrack->Phi();
388 //set default phi weight to 1
389 Double_t dW = 1.;
390 //phi weight of pTrack
391 if(fUsePhiWeights && fPhiWeights) {
392 Int_t iNBinsPhi = fPhiWeights->GetNbinsX();
393 dW = fPhiWeights->GetBinContent(1+(Int_t)(TMath::Floor(dPhi*iNBinsPhi/TMath::TwoPi())));
394 //bin = 1 + value*nbins/range
395 //TMath::Floor rounds to the lower integer
396 }
12dc476e 397
a554b203 398 //calculate vU
399 TVector2 vU;
400 Double_t dUX = TMath::Cos(2*dPhi);
401 Double_t dUY = TMath::Sin(2*dPhi);
402 vU.Set(dUX,dUY);
403 Double_t dModulus = vU.Mod();
404 if (dModulus!=0.) vU.Set(dUX/dModulus,dUY/dModulus); // make length 1
405 else cerr<<"dModulus is zero!"<<endl;
12dc476e 406
a554b203 407 //redefine the Q vector and devide by its multiplicity
408 TVector2 vQm;
409 Double_t dQmX = 0.;
410 Double_t dQmY = 0.;
411 //subtract particle from the flowvector if used to define it
412 if (pTrack->InSubevent(0) || pTrack->InSubevent(1)) {
413 dQmX = (vQ.X() - dW*dUX)/(dMq-1);
414 dQmY = (vQ.Y() - dW*dUY)/(dMq-1);
415 vQm.Set(dQmX,dQmY);
7a01f4a7 416
a554b203 417 //dUQ = scalar product of vU and vQm
418 Double_t dUQ = (vU * vQm);
419 Double_t dPt = pTrack->Pt();
420 Double_t dEta = pTrack->Eta();
421 //fill the profile histograms
422 if (pTrack->InRPSelection()) {
423 fHistProUQetaRP -> Fill(dEta,dUQ,(dMq-1)); //Fill (Qu/(Mq-1)) with weight (Mq-1)
424 //needed for the error calculation:
425 fHistProUQQaQbEtaRP -> Fill(dEta,dUQ*dQaQb,(dMq-1)*dMa*dMb); //Fill [Qu/(Mq-1)]*[QaQb/MaMb] with weight (Mq-1)MaMb
426 fHistProUQPtRP -> Fill(dPt,dUQ,(dMq-1)); //Fill (Qu/(Mq-1)) with weight (Mq-1)
427 fHistProUQQaQbPtRP -> Fill(dPt,dUQ*dQaQb,(dMq-1)*dMa*dMb); //Fill [Qu/(Mq-1)]*[QaQb/MaMb] with weight (Mq-1)MaMb
428
429 fHistSumOfWeightsEtaRP[0]->Fill(dEta,(dMq-1)); // sum of Mq-1
430 fHistSumOfWeightsEtaRP[1]->Fill(dEta,pow((dMq-1),2.));// sum of (Mq-1)^2
431 fHistSumOfWeightsEtaRP[2]->Fill(dEta,(dMq-1)*dMa*dMb);// sum of (Mq-1)*MaMb
432 fHistSumOfWeightsPtRP[0]->Fill(dPt,(dMq-1)); // sum of Mq-1
433 fHistSumOfWeightsPtRP[1]->Fill(dPt,pow((dMq-1),2.)); // sum of (Mq-1)^2
434 fHistSumOfWeightsPtRP[2]->Fill(dPt,(dMq-1)*dMa*dMb); // sum of (Mq-1)*MaMb
435 }
436 if (pTrack->InPOISelection()) {
437 fHistProUQetaPOI -> Fill(dEta,dUQ,(dMq-1));//Fill (Qu/(Mq-1)) with weight (Mq-1)
438 //needed for the error calculation:
439 fHistProUQQaQbEtaPOI -> Fill(dEta,dUQ*dQaQb,(dMq-1)*dMa*dMb); //Fill [Qu/(Mq-1)]*[QaQb/MaMb] with weight (Mq-1)MaMb
440 fHistProUQPtPOI -> Fill(dPt,dUQ,(dMq-1)); //Fill (Qu/(Mq-1)) with weight (Mq-1)
441 fHistProUQQaQbPtPOI -> Fill(dPt,dUQ*dQaQb,(dMq-1)*dMa*dMb); //Fill [Qu/(Mq-1)]*[QaQb/MaMb] with weight (Mq-1)MaMb
442
443 fHistSumOfWeightsEtaPOI[0]->Fill(dEta,(dMq-1)); // sum of Mq-1
444 fHistSumOfWeightsEtaPOI[1]->Fill(dEta,pow((dMq-1),2.));// sum of (Mq-1)^2
445 fHistSumOfWeightsEtaPOI[2]->Fill(dEta,(dMq-1)*dMa*dMb);// sum of (Mq-1)*MaMb
446 fHistSumOfWeightsPtPOI[0]->Fill(dPt,(dMq-1)); // sum of Mq-1
447 fHistSumOfWeightsPtPOI[1]->Fill(dPt,pow((dMq-1),2.)); // sum of (Mq-1)^2
448 fHistSumOfWeightsPtPOI[2]->Fill(dPt,(dMq-1)*dMa*dMb); // sum of (Mq-1)*MaMb
449 }
7a01f4a7 450
a554b203 451 } else { //do not subtract the particle from the flowvector
452 dQmX = vQ.X()/dMq;
453 dQmY = vQ.Y()/dMq;
454 vQm.Set(dQmX,dQmY);
7a01f4a7 455
a554b203 456 //dUQ = scalar product of vU and vQm
457 Double_t dUQ = (vU * vQm);
458 Double_t dPt = pTrack->Pt();
459 Double_t dEta = pTrack->Eta();
460 //fill the profile histograms
461 if (pTrack->InRPSelection()) {
462 fHistProUQetaRP -> Fill(dEta,dUQ,dMq); //Fill (Qu/Mq) with weight Mq
463 //needed for the error calculation:
464 fHistProUQQaQbEtaRP -> Fill(dEta,dUQ*dQaQb,dMq*dMa*dMb); //Fill [Qu/Mq]*[QaQb/MaMb] with weight Mq*MaMb
465 fHistProUQPtRP -> Fill(dPt,dUQ,dMq); //Fill (Qu/Mq) with weight Mq
466 fHistProUQQaQbPtRP -> Fill(dPt,dUQ*dQaQb,dMq*dMa*dMb); //Fill [Qu/Mq]*[QaQb/MaMb] with weight Mq*MaMb
467
468 fHistSumOfWeightsEtaRP[0]->Fill(dEta,dMq); // sum of Mq
469 fHistSumOfWeightsEtaRP[1]->Fill(dEta,pow(dMq,2.));// sum of Mq^2
470 fHistSumOfWeightsEtaRP[2]->Fill(dEta,dMq*dMa*dMb);// sum of Mq*MaMb
471 fHistSumOfWeightsPtRP[0]->Fill(dPt,dMq); // sum of Mq
472 fHistSumOfWeightsPtRP[1]->Fill(dPt,pow(dMq,2.)); // sum of Mq^2
473 fHistSumOfWeightsPtRP[2]->Fill(dPt,dMq*dMa*dMb); // sum of Mq*MaMb
474 }
475 if (pTrack->InPOISelection()) {
476 fHistProUQetaPOI -> Fill(dEta,dUQ,dMq); //Fill (Qu/Mq) with weight Mq
477 //needed for the error calculation:
478 fHistProUQQaQbEtaPOI -> Fill(dEta,dUQ*dQaQb,dMq*dMa*dMb); //Fill [Qu/Mq]*[QaQb/MaMb] with weight Mq*MaMb
479 fHistProUQPtPOI -> Fill(dPt,dUQ,dMq); //Fill (Qu/Mq) with weight Mq
480 fHistProUQQaQbPtPOI -> Fill(dPt,dUQ*dQaQb,dMq*dMa*dMb); //Fill [Qu/Mq]*[QaQb/MaMb] with weight Mq*MaMb
481
482 fHistSumOfWeightsEtaPOI[0]->Fill(dEta,dMq); // sum of Mq
483 fHistSumOfWeightsEtaPOI[1]->Fill(dEta,pow(dMq,2.));// sum of Mq^2
484 fHistSumOfWeightsEtaPOI[2]->Fill(dEta,dMq*dMa*dMb);// sum of Mq*MaMb
485 fHistSumOfWeightsPtPOI[0]->Fill(dPt,dMq); // sum of Mq
486 fHistSumOfWeightsPtPOI[1]->Fill(dPt,pow(dMq,2.)); // sum of Mq^2
487 fHistSumOfWeightsPtPOI[2]->Fill(dPt,dMq*dMa*dMb); // sum of Mq*MaMb
488 }
489 }//track not in subevents
7a01f4a7 490
a554b203 491 }//track
7a01f4a7 492
a554b203 493 }//loop over tracks
494
495 fEventNumber++;
496
497 } //difference Ma and Mb
12dc476e 498
7a01f4a7 499 }// subevents not empty
b125a454 500 delete [] vQarray;
a554b203 501
7a01f4a7 502 } //event
a554b203 503
7a01f4a7 504}//end of Make()
8d312f00 505
12dc476e 506//--------------------------------------------------------------------
fd46c3dd 507void AliFlowAnalysisWithScalarProduct::GetOutputHistograms(TList *outputListHistos){
508
509 //get pointers to all output histograms (called before Finish())
12dc476e 510
fd46c3dd 511 if (outputListHistos) {
512 //Get the common histograms from the output list
513 AliFlowCommonHist *pCommonHist = dynamic_cast<AliFlowCommonHist*>
514 (outputListHistos->FindObject("AliFlowCommonHistSP"));
515 AliFlowCommonHistResults *pCommonHistResults = dynamic_cast<AliFlowCommonHistResults*>
516 (outputListHistos->FindObject("AliFlowCommonHistResultsSP"));
517 TProfile* pHistProQaQb = dynamic_cast<TProfile*>(outputListHistos->FindObject("Flow_QaQb_SP"));
a554b203 518 TProfile* pHistProQaQbNorm = dynamic_cast<TProfile*>(outputListHistos->FindObject("FlowPro_QaQbNorm_SP"));
519 TH1D* pHistSumOfLinearWeights = dynamic_cast<TH1D*>(outputListHistos->FindObject("Flow_SumOfLinearWeights_SP"));
520 TH1D* pHistSumOfQuadraticWeights = dynamic_cast<TH1D*>(outputListHistos->FindObject("Flow_SumOfQuadraticWeights_SP"));
521
fd46c3dd 522 TProfile* pHistProUQetaRP = dynamic_cast<TProfile*>(outputListHistos->FindObject("Flow_UQetaRP_SP"));
523 TProfile* pHistProUQetaPOI = dynamic_cast<TProfile*>(outputListHistos->FindObject("Flow_UQetaPOI_SP"));
524 TProfile* pHistProUQPtRP = dynamic_cast<TProfile*>(outputListHistos->FindObject("Flow_UQPtRP_SP"));
525 TProfile* pHistProUQPtPOI = dynamic_cast<TProfile*>(outputListHistos->FindObject("Flow_UQPtPOI_SP"));
12dc476e 526 TProfile* pHistProUQQaQbPtRP = dynamic_cast<TProfile*>(outputListHistos->FindObject("Flow_UQQaQbPtRP_SP"));
527 TProfile* pHistProUQQaQbEtaRP = dynamic_cast<TProfile*>(outputListHistos->FindObject("Flow_UQQaQbEtaRP_SP"));
528 TProfile* pHistProUQQaQbPtPOI = dynamic_cast<TProfile*>(outputListHistos->FindObject("Flow_UQQaQbPtPOI_SP"));
529 TProfile* pHistProUQQaQbEtaPOI = dynamic_cast<TProfile*>(outputListHistos->FindObject("Flow_UQQaQbEtaPOI_SP"));
12dc476e 530 TString weightFlag[3] = {"w_Qu_","w_Qu^2_","w_QuQaQb_"};
531 TH1D* pHistSumOfWeightsPtRP[3] = {NULL};
532 TH1D* pHistSumOfWeightsEtaRP[3] = {NULL};
533 TH1D* pHistSumOfWeightsPtPOI[3] = {NULL};
534 TH1D* pHistSumOfWeightsEtaPOI[3] = {NULL};
535 for(Int_t i=0;i<3;i++)
536 {
cff1bec4 537 pHistSumOfWeightsPtRP[i] = dynamic_cast<TH1D*>(outputListHistos->FindObject(Form("Flow_SumOfWeights%sPtRP_SP",weightFlag[i].Data())));
538 pHistSumOfWeightsEtaRP[i] = dynamic_cast<TH1D*>(outputListHistos->FindObject(Form("Flow_SumOfWeights%sEtaRP_SP",weightFlag[i].Data())));
539 pHistSumOfWeightsPtPOI[i] = dynamic_cast<TH1D*>(outputListHistos->FindObject(Form("Flow_SumOfWeights%sPtPOI_SP",weightFlag[i].Data())));
540 pHistSumOfWeightsEtaPOI[i] = dynamic_cast<TH1D*>(outputListHistos->FindObject(Form("Flow_SumOfWeights%sEtaPOI_SP",weightFlag[i].Data())));
12dc476e 541 }
542
a554b203 543 TH1D* pHistQaQb = dynamic_cast<TH1D*>(outputListHistos->FindObject("Flow_QaQb_SP"));
544 TH1D* pHistQaQbNorm = dynamic_cast<TH1D*>(outputListHistos->FindObject("Flow_QaQbNorm_SP"));
545 TH1D* pHistQaQbCos = dynamic_cast<TH1D*>(outputListHistos->FindObject("Flow_QaQbCos_SP"));
546 TH1D* pHistResolution = dynamic_cast<TH1D*>(outputListHistos->FindObject("Flow_resolution_SP"));
547 TH1D* pHistQaNorm = dynamic_cast<TH1D*>(outputListHistos->FindObject("Flow_QaNorm_SP"));
548 TH2D* pHistQaNormvsMa = dynamic_cast<TH2D*>(outputListHistos->FindObject("Flow_QaNormvsMa_SP"));
549 TH1D* pHistQbNorm = dynamic_cast<TH1D*>(outputListHistos->FindObject("Flow_QbNorm_SP"));
550 TH2D* pHistQbNormvsMb = dynamic_cast<TH2D*>(outputListHistos->FindObject("Flow_QbNormvsMb_SP"));
551 TH2D* pHistMavsMb = dynamic_cast<TH2D*>(outputListHistos->FindObject("Flow_MavsMb_SP"));
552
12dc476e 553 //pass the pointers to the task
a554b203 554 if (pCommonHist && pCommonHistResults && pHistProQaQb && pHistProQaQbNorm &&
555 pHistSumOfLinearWeights && pHistSumOfQuadraticWeights &&
556 pHistQaQb && pHistQaQbNorm && pHistQaQbCos && pHistResolution &&
557 pHistQaNorm && pHistQaNormvsMa && pHistQbNorm && pHistQbNormvsMb &&
558 pHistMavsMb &&
559 pHistProUQetaRP && pHistProUQetaPOI &&
560 pHistProUQPtRP && pHistProUQPtPOI &&
561 pHistProUQQaQbPtRP && pHistProUQQaQbEtaRP &&
562 pHistProUQQaQbPtPOI && pHistProUQQaQbEtaPOI) {
fd46c3dd 563 this -> SetCommonHists(pCommonHist);
564 this -> SetCommonHistsRes(pCommonHistResults);
565 this -> SetHistProQaQb(pHistProQaQb);
a554b203 566 this -> SetHistProQaQbNorm(pHistProQaQbNorm);
567 this -> SetHistSumOfLinearWeights(pHistSumOfLinearWeights);
568 this -> SetHistSumOfQuadraticWeights(pHistSumOfQuadraticWeights);
569 this -> SetHistQaQb(pHistQaQb);
570 this -> SetHistQaQbNorm(pHistQaQbNorm);
571 this -> SetHistQaQbCos(pHistQaQbCos);
572 this -> SetHistResolution(pHistResolution);
573 this -> SetHistQaNorm(pHistQaNorm);
574 this -> SetHistQaNormvsMa(pHistQaNormvsMa);
575 this -> SetHistQbNorm(pHistQbNorm);
576 this -> SetHistQbNormvsMb(pHistQbNormvsMb);
577 this -> SetHistMavsMb(pHistMavsMb);
fd46c3dd 578 this -> SetHistProUQetaRP(pHistProUQetaRP);
579 this -> SetHistProUQetaPOI(pHistProUQetaPOI);
580 this -> SetHistProUQPtRP(pHistProUQPtRP);
581 this -> SetHistProUQPtPOI(pHistProUQPtPOI);
12dc476e 582 this -> SetHistProUQQaQbPtRP(pHistProUQQaQbPtRP);
583 this -> SetHistProUQQaQbEtaRP(pHistProUQQaQbEtaRP);
584 this -> SetHistProUQQaQbPtPOI(pHistProUQQaQbPtPOI);
585 this -> SetHistProUQQaQbEtaPOI(pHistProUQQaQbEtaPOI);
a554b203 586
587 for(Int_t i=0;i<3;i++) {
588 if(pHistSumOfWeightsPtRP[i]) this -> SetHistSumOfWeightsPtRP(pHistSumOfWeightsPtRP[i],i);
589 if(pHistSumOfWeightsEtaRP[i]) this -> SetHistSumOfWeightsEtaRP(pHistSumOfWeightsEtaRP[i],i);
590 if(pHistSumOfWeightsPtPOI[i]) this -> SetHistSumOfWeightsPtPOI(pHistSumOfWeightsPtPOI[i],i);
591 if(pHistSumOfWeightsEtaPOI[i]) this -> SetHistSumOfWeightsEtaPOI(pHistSumOfWeightsEtaPOI[i],i);
592 }
593
594 } else {
595 cout<<"WARNING: Histograms needed to run Finish() in SP are not accessable!"<<endl; }
596
12dc476e 597 } // end of if(outputListHistos)
fd46c3dd 598}
599
600//--------------------------------------------------------------------
8d312f00 601void AliFlowAnalysisWithScalarProduct::Finish() {
602
489d5531 603 //calculate flow and fill the AliFlowCommonHistResults
395fadba 604 if (fDebug) cout<<"AliFlowAnalysisWithScalarProduct::Finish()"<<endl;
605
606 cout<<"*************************************"<<endl;
607 cout<<"*************************************"<<endl;
608 cout<<" Integrated flow from "<<endl;
609 cout<<" Scalar product "<<endl;
610 cout<<endl;
e4cecfa0 611
bee2efdc 612 Int_t iNbinsPt = AliFlowCommonConstants::GetMaster()->GetNbinsPt();
613 Int_t iNbinsEta = AliFlowCommonConstants::GetMaster()->GetNbinsEta();
12dc476e 614
a554b203 615 //Calculate the event plane resolution
616 //----------------------------------
617 Double_t dCos2phi = fHistResolution->GetMean();
7f9d91f9 618 if (dCos2phi > 0.0){
a554b203 619 Double_t dResolution = TMath::Sqrt(2*dCos2phi);
620 cout<<"An estimate of the event plane resolution is: "<<dResolution<<endl;
621 cout<<endl;
622 }
623
12dc476e 624 //Calculate reference flow (noname)
625 //----------------------------------
626 //weighted average over (QaQb/MaMb) with weight (MaMb)
a554b203 627 Double_t dQaQb = fHistProQaQbNorm->GetBinContent(1);
628 Double_t dV = -999.;
524798bf 629 if(dQaQb>=0.)
630 {
631 dV = TMath::Sqrt(dQaQb);
632 }
12dc476e 633 //statistical error of dQaQb:
634 // statistical error = term1 * spread * term2:
635 // term1 = sqrt{sum_{i=1}^{N} w^2}/(sum_{i=1}^{N} w)
636 // term2 = 1/sqrt(1-term1^2)
a554b203 637 Double_t dSpreadQaQb = fHistProQaQbNorm->GetBinError(1);
12dc476e 638 Double_t dSumOfLinearWeights = fHistSumOfLinearWeights->GetBinContent(1);
639 Double_t dSumOfQuadraticWeights = fHistSumOfQuadraticWeights->GetBinContent(1);
640 Double_t dTerm1 = 0.;
641 Double_t dTerm2 = 0.;
642 if(dSumOfLinearWeights) {
643 dTerm1 = pow(dSumOfQuadraticWeights,0.5)/dSumOfLinearWeights;
644 }
645 if(1.-pow(dTerm1,2.) > 0.) {
646 dTerm2 = 1./pow(1-pow(dTerm1,2.),0.5);
647 }
648 Double_t dStatErrorQaQb = dTerm1 * dSpreadQaQb * dTerm2;
649 //calculate the statistical error
650 Double_t dVerr = 0.;
651 if(dQaQb > 0.) {
652 dVerr = (1./(2.*pow(dQaQb,0.5)))*dStatErrorQaQb;
653 }
654 fCommonHistsRes->FillIntegratedFlow(dV,dVerr);
655 cout<<"dV = "<<dV<<" +- "<<dVerr<<endl;
656
657 //Calculate differential flow and integrated flow (RP, POI)
658 //---------------------------------------------------------
12dc476e 659 //v as a function of eta for RP selection
a9de3704 660 for(Int_t b=1;b<iNbinsEta+1;b++) {
12dc476e 661 Double_t duQpro = fHistProUQetaRP->GetBinContent(b);
a554b203 662 Double_t dv2pro = -999.;
12dc476e 663 if (dV!=0.) { dv2pro = duQpro/dV; }
664 //calculate the statistical error
665 Double_t dv2ProErr = CalculateStatisticalError(b, dStatErrorQaQb, fHistProUQetaRP, fHistProUQQaQbEtaRP, fHistSumOfWeightsEtaRP);
12dc476e 666 //fill TH1D
cff1bec4 667 fCommonHistsRes->FillDifferentialFlowEtaRP(b, dv2pro, dv2ProErr);
12dc476e 668 } //loop over bins b
669
80f72ed6 670
12dc476e 671 //v as a function of eta for POI selection
a9de3704 672 for(Int_t b=1;b<iNbinsEta+1;b++) {
12dc476e 673 Double_t duQpro = fHistProUQetaPOI->GetBinContent(b);
a554b203 674 Double_t dv2pro = -999.;
12dc476e 675 if (dV!=0.) { dv2pro = duQpro/dV; }
676 //calculate the statistical error
677 Double_t dv2ProErr = CalculateStatisticalError(b, dStatErrorQaQb, fHistProUQetaPOI, fHistProUQQaQbEtaPOI, fHistSumOfWeightsEtaPOI);
678
679 //fill TH1D
680 fCommonHistsRes->FillDifferentialFlowEtaPOI(b, dv2pro, dv2ProErr);
681 } //loop over bins b
682
80f72ed6 683
12dc476e 684 //v as a function of Pt for RP selection
685 TH1F* fHistPtRP = fCommonHists->GetHistPtRP(); //for calculating integrated flow
686 Double_t dVRP = 0.;
687 Double_t dSumRP = 0.;
688 Double_t dErrVRP =0.;
689
a9de3704 690 for(Int_t b=1;b<iNbinsPt+1;b++) {
12dc476e 691 Double_t duQpro = fHistProUQPtRP->GetBinContent(b);
a554b203 692 Double_t dv2pro = -999.;
12dc476e 693 if (dV!=0.) { dv2pro = duQpro/dV; }
694 //calculate the statistical error
695 Double_t dv2ProErr = CalculateStatisticalError(b, dStatErrorQaQb, fHistProUQPtRP, fHistProUQQaQbPtRP, fHistSumOfWeightsPtRP);
696
697 //fill TH1D
698 fCommonHistsRes->FillDifferentialFlowPtRP(b, dv2pro, dv2ProErr);
699
700 //calculate integrated flow for RP selection
701 if (fHistPtRP){
702 Double_t dYieldPt = fHistPtRP->GetBinContent(b);
703 dVRP += dv2pro*dYieldPt;
704 dSumRP +=dYieldPt;
705 dErrVRP += dYieldPt*dYieldPt*dv2ProErr*dv2ProErr;
706 } else { cout<<"fHistPtRP is NULL"<<endl; }
707 } //loop over bins b
708
709 if (dSumRP != 0.) {
710 dVRP /= dSumRP; //the pt distribution should be normalised
711 dErrVRP /= (dSumRP*dSumRP);
712 dErrVRP = TMath::Sqrt(dErrVRP);
80f72ed6 713 }
12dc476e 714 fCommonHistsRes->FillIntegratedFlowRP(dVRP,dErrVRP);
715 cout<<"dV(RP) = "<<dVRP<<" +- "<<dErrVRP<<endl;
716
80f72ed6 717
12dc476e 718 //v as a function of Pt for POI selection
719 TH1F* fHistPtPOI = fCommonHists->GetHistPtPOI(); //for calculating integrated flow
720 Double_t dVPOI = 0.;
721 Double_t dSumPOI = 0.;
722 Double_t dErrVPOI =0.;
723
a9de3704 724 for(Int_t b=1;b<iNbinsPt+1;b++) {
12dc476e 725 Double_t duQpro = fHistProUQPtPOI->GetBinContent(b);
a554b203 726 Double_t dv2pro = -999.;
12dc476e 727 if (dV!=0.) { dv2pro = duQpro/dV; }
728 //calculate the statistical error
729 Double_t dv2ProErr = CalculateStatisticalError(b, dStatErrorQaQb, fHistProUQPtPOI, fHistProUQQaQbPtPOI, fHistSumOfWeightsPtPOI);
730
731 //fill TH1D
732 fCommonHistsRes->FillDifferentialFlowPtPOI(b, dv2pro, dv2ProErr);
395fadba 733
12dc476e 734 //calculate integrated flow for POI selection
735 if (fHistPtPOI){
736 Double_t dYieldPt = fHistPtPOI->GetBinContent(b);
737 dVPOI += dv2pro*dYieldPt;
738 dSumPOI +=dYieldPt;
739 dErrVPOI += dYieldPt*dYieldPt*dv2ProErr*dv2ProErr;
740 } else { cout<<"fHistPtPOI is NULL"<<endl; }
741 } //loop over bins b
395fadba 742
12dc476e 743 if (dSumPOI != 0.) {
489d5531 744 dVPOI /= dSumPOI; //the pt distribution should be normalised
12dc476e 745 dErrVPOI /= (dSumPOI*dSumPOI);
746 dErrVPOI = TMath::Sqrt(dErrVPOI);
747 }
748 fCommonHistsRes->FillIntegratedFlowPOI(dVPOI,dErrVPOI);
749 cout<<"dV(POI) = "<<dVPOI<<" +- "<<dErrVPOI<<endl;
e4cecfa0 750
395fadba 751 cout<<endl;
752 cout<<"*************************************"<<endl;
753 cout<<"*************************************"<<endl;
12dc476e 754
395fadba 755 //cout<<".....finished"<<endl;
12dc476e 756}
8232a5ec 757
758
12dc476e 759//--------------------------------------------------------------------
760Double_t AliFlowAnalysisWithScalarProduct::CalculateStatisticalError(Int_t b, Double_t aStatErrorQaQb, TProfile* pHistProUQ, TProfile* pHistProUQQaQb, TH1D** pHistSumOfWeights) {
12dc476e 761 //calculate the statistical error for differential flow for bin b
762 Double_t duQproSpread = pHistProUQ->GetBinError(b);
763 Double_t sumOfMq = pHistSumOfWeights[0]->GetBinContent(b);
764 Double_t sumOfMqSquared = pHistSumOfWeights[1]->GetBinContent(b);
a554b203 765 Double_t dQaQb = fHistProQaQbNorm->GetBinContent(1);
12dc476e 766 Double_t dTerm1 = 0.;
767 Double_t dTerm2 = 0.;
768 if(sumOfMq) {
769 dTerm1 = (pow(sumOfMqSquared,0.5)/sumOfMq);
770 }
771 if(1.-pow(dTerm1,2.)>0.) {
772 dTerm2 = 1./pow(1.-pow(dTerm1,2.),0.5);
773 }
774 Double_t duQproErr = dTerm1*duQproSpread*dTerm2;
12dc476e 775 // covariances:
776 Double_t dTerm1Cov = pHistSumOfWeights[2]->GetBinContent(b);
777 Double_t dTerm2Cov = fHistSumOfLinearWeights->GetBinContent(1);
778 Double_t dTerm3Cov = sumOfMq;
779 Double_t dWeightedCovariance = 0.;
780 if(dTerm2Cov*dTerm3Cov>0.) {
781 Double_t dDenominator = 1.-dTerm1Cov/(dTerm2Cov*dTerm3Cov);
782 Double_t dPrefactor = dTerm1Cov/(dTerm2Cov*dTerm3Cov);
783 if(dDenominator!=0) {
a554b203 784 Double_t dCovariance = (pHistProUQQaQb->GetBinContent(b)-dQaQb*pHistProUQ->GetBinContent(b))/dDenominator;
12dc476e 785 dWeightedCovariance = dCovariance*dPrefactor;
786 }
787 }
cff1bec4 788
12dc476e 789 Double_t dv2ProErr = 0.; // final statitical error:
a554b203 790 if(dQaQb>0.) {
791 Double_t dv2ProErrorSquared = (1./4.)*pow(dQaQb,-3.)*
12dc476e 792 (pow(pHistProUQ->GetBinContent(b),2.)*pow(aStatErrorQaQb,2.)
a554b203 793 + 4.*pow(dQaQb,2.)*pow(duQproErr,2.)
794 - 4.*dQaQb*pHistProUQ->GetBinContent(b)*dWeightedCovariance);
12dc476e 795 if(dv2ProErrorSquared>0.) dv2ProErr = pow(dv2ProErrorSquared,0.5);
796 }
797
798 return dv2ProErr;
799}
800//--------------------------------------------------------------------