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