New base class for trigger modules handling generic events.
[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
18#include "Riostream.h" //needed as include
1dfa3c16 19#include "TFile.h" //needed as include
d7eb18ec 20#include "TList.h"
8d312f00 21#include "TMath.h"
22#include "TProfile.h"
23#include "TVector2.h"
24
25class TH1F;
26
27#include "AliFlowCommonConstants.h" //needed as include
28#include "AliFlowEventSimple.h"
29#include "AliFlowTrackSimple.h"
30#include "AliFlowCommonHist.h"
395fadba 31#include "AliFlowCommonHistResults.h"
8d312f00 32#include "AliFlowAnalysisWithScalarProduct.h"
33
34class AliFlowVector;
35
36// AliFlowAnalysisWithScalarProduct:
37// Description:
38// Maker to analyze Flow with the Scalar product method.
39//
40// author: N. van der Kolk (kolk@nikhef.nl)
41
42ClassImp(AliFlowAnalysisWithScalarProduct)
43
44 //-----------------------------------------------------------------------
45
46 AliFlowAnalysisWithScalarProduct::AliFlowAnalysisWithScalarProduct():
e35ddff0 47 fEventNumber(0),
48 fDebug(kFALSE),
29195b69 49 fWeightsList(NULL),
50 fUsePhiWeights(kFALSE),
51 fPhiWeights(NULL),
e35ddff0 52 fHistList(NULL),
395fadba 53 fHistProUQetaRP(NULL),
54 fHistProUQetaPOI(NULL),
55 fHistProUQPtRP(NULL),
56 fHistProUQPtPOI(NULL),
e4cecfa0 57 fHistProQaQb(NULL),
92041674 58 fHistProM(NULL),
29195b69 59 fHistM(NULL),
395fadba 60 fCommonHists(NULL),
61 fCommonHistsRes(NULL)
8d312f00 62{
8d312f00 63 // Constructor.
29195b69 64 fWeightsList = new TList();
e35ddff0 65 fHistList = new TList();
8d312f00 66}
67 //-----------------------------------------------------------------------
68
69
70 AliFlowAnalysisWithScalarProduct::~AliFlowAnalysisWithScalarProduct()
71 {
72 //destructor
29195b69 73 delete fWeightsList;
e35ddff0 74 delete fHistList;
8d312f00 75 }
76
77
1dfa3c16 78//-----------------------------------------------------------------------
79
80void AliFlowAnalysisWithScalarProduct::WriteHistograms(TString* outputFileName)
81{
82 //store the final results in output .root file
83
84 TFile *output = new TFile(outputFileName->Data(),"RECREATE");
0fe80f88 85 //output->WriteObject(fHistList, "cobjSP","SingleKey");
86 fHistList->SetName("cobjSP");
9455e15e 87 fHistList->SetOwner(kTRUE);
0fe80f88 88 fHistList->Write(fHistList->GetName(), TObject::kSingleKey);
1dfa3c16 89 delete output;
90}
91
b0fda271 92//-----------------------------------------------------------------------
93
94void AliFlowAnalysisWithScalarProduct::WriteHistograms(TString outputFileName)
95{
96 //store the final results in output .root file
97
98 TFile *output = new TFile(outputFileName.Data(),"RECREATE");
0fe80f88 99 //output->WriteObject(fHistList, "cobjSP","SingleKey");
100 fHistList->SetName("cobjSP");
9455e15e 101 fHistList->SetOwner(kTRUE);
0fe80f88 102 fHistList->Write(fHistList->GetName(), TObject::kSingleKey);
b0fda271 103 delete output;
104}
105
ad87ae62 106//-----------------------------------------------------------------------
107
108void AliFlowAnalysisWithScalarProduct::WriteHistograms(TDirectoryFile *outputFileName)
109{
110 //store the final results in output .root file
111 fHistList->SetName("cobjSP");
112 fHistList->SetOwner(kTRUE);
113 outputFileName->Add(fHistList);
114 outputFileName->Write(outputFileName->GetName(), TObject::kSingleKey);
115}
116
8d312f00 117//-----------------------------------------------------------------------
118void AliFlowAnalysisWithScalarProduct::Init() {
119
120 //Define all histograms
d7eb18ec 121 cout<<"---Analysis with the Scalar Product Method--- Init"<<endl;
8d312f00 122
bee2efdc 123 Int_t iNbinsPt = AliFlowCommonConstants::GetMaster()->GetNbinsPt();
124 Double_t dPtMin = AliFlowCommonConstants::GetMaster()->GetPtMin();
125 Double_t dPtMax = AliFlowCommonConstants::GetMaster()->GetPtMax();
126 Int_t iNbinsEta = AliFlowCommonConstants::GetMaster()->GetNbinsEta();
127 Double_t dEtaMin = AliFlowCommonConstants::GetMaster()->GetEtaMin();
128 Double_t dEtaMax = AliFlowCommonConstants::GetMaster()->GetEtaMax();
8d312f00 129
395fadba 130 fHistProUQetaRP = new TProfile("Flow_UQetaRP_SP","Flow_UQetaRP_SP",iNbinsEta,dEtaMin,dEtaMax);
131 fHistProUQetaRP->SetXTitle("{eta}");
132 fHistProUQetaRP->SetYTitle("<uQ>");
133 fHistList->Add(fHistProUQetaRP);
134
135 fHistProUQetaPOI = new TProfile("Flow_UQetaPOI_SP","Flow_UQetaPOI_SP",iNbinsEta,dEtaMin,dEtaMax);
136 fHistProUQetaPOI->SetXTitle("{eta}");
137 fHistProUQetaPOI->SetYTitle("<uQ>");
138 fHistList->Add(fHistProUQetaPOI);
139
140 fHistProUQPtRP = new TProfile("Flow_UQPtRP_SP","Flow_UQPtRP_SP",iNbinsPt,dPtMin,dPtMax);
141 fHistProUQPtRP->SetXTitle("p_t (GeV)");
142 fHistProUQPtRP->SetYTitle("<uQ>");
143 fHistList->Add(fHistProUQPtRP);
144
145 fHistProUQPtPOI = new TProfile("Flow_UQPtPOI_SP","Flow_UQPtPOI_SP",iNbinsPt,dPtMin,dPtMax);
146 fHistProUQPtPOI->SetXTitle("p_t (GeV)");
147 fHistProUQPtPOI->SetYTitle("<uQ>");
148 fHistList->Add(fHistProUQPtPOI);
149
e4cecfa0 150 fHistProQaQb = new TProfile("Flow_QaQb_SP","Flow_QaQb_SP", 1, -0.5, 0.5);
151 fHistProQaQb->SetYTitle("<QaQb>");
152 fHistList->Add(fHistProQaQb);
8d312f00 153
92041674 154 fHistProM = new TProfile("Flow_M_SP","Flow_M_SP",2,0.5, 2.5);
155 fHistProM -> SetYTitle("<*>");
156 fHistProM -> SetXTitle("<M-1>, <Ma*Mb>");
157 fHistList->Add(fHistProM);
158
29195b69 159 fHistM = new TH1D("Flow_Msum_SP","Flow_Msum_SP",2,0.5, 2.5);
160 fHistM -> SetYTitle("sum (*)");
161 fHistM -> SetXTitle("sum (M-1), sum (Ma*Mb)");
162 fHistList->Add(fHistM);
163
04f6283b 164 fCommonHists = new AliFlowCommonHist("AliFlowCommonHistSP");
099e1753 165 fHistList->Add(fCommonHists);
395fadba 166 fCommonHistsRes = new AliFlowCommonHistResults("AliFlowCommonHistResultsSP");
167 fHistList->Add(fCommonHistsRes);
168
29195b69 169 //weights
170 if(fUsePhiWeights) {
171 if(!fWeightsList) {
172 cout<<"WARNING: fWeightsList is NULL in the Scalar Product method."<<endl;
173 exit(0);
174 }
175 if(fWeightsList->FindObject("phi_weights")) {
176 fPhiWeights = dynamic_cast<TH1F*>(fWeightsList->FindObject("phi_weights"));
177 fHistList->Add(fPhiWeights);
178 } else {
179 cout<<"WARNING: histogram with phi weights is not accessible in Scalar Product"<<endl;
180 exit(0);
181 }
182 } // end of if(fUsePhiWeights)
183
e2d51347 184 fEventNumber = 0; //set number of events to zero
185}
186
8d312f00 187//-----------------------------------------------------------------------
188
8232a5ec 189void AliFlowAnalysisWithScalarProduct::Make(AliFlowEventSimple* anEvent) {
8d312f00 190
191 //Fill histogram
8232a5ec 192 if (anEvent) {
8d312f00 193
194 //fill control histograms
8232a5ec 195 fCommonHists->FillControlHistograms(anEvent);
cbbaf54a 196
197 //get Q vectors for the eta-subevents
b125a454 198 AliFlowVector* vQarray = new AliFlowVector[2];
29195b69 199 if (fUsePhiWeights) {
200 anEvent->Get2Qsub(vQarray,2,fWeightsList,kTRUE);
201 } else {
202 anEvent->Get2Qsub(vQarray);
203 }
b125a454 204 AliFlowVector vQa = vQarray[0];
205 AliFlowVector vQb = vQarray[1];
cbbaf54a 206 //get total Q vector
207 AliFlowVector vQ = vQa + vQb;
208
209 //fill the multiplicity histograms for the prefactor
29195b69 210 Double_t dMmin1 = vQ.GetMult()-1; //M-1
211 Double_t dMaMb = vQa.GetMult()*vQb.GetMult(); //Ma*Mb
212 fHistM -> Fill(1,dMmin1);//sum of (M-1)
213 fHistM -> Fill(2,vQa.GetMult()*vQb.GetMult()); //sum of (Ma*Mb)
214 fHistProM -> Fill(1,dMmin1); //<M-1>
215 fHistProM -> Fill(2,vQa.GetMult()*vQb.GetMult()); //<Ma*Mb>
cbbaf54a 216 //scalar product of the two subevents
29195b69 217 Double_t dQaQb = (vQa*vQb)*dMaMb;
e4cecfa0 218 fHistProQaQb -> Fill(0.,dQaQb);
e35ddff0 219
8d312f00 220 //loop over the tracks of the event
e35ddff0 221 AliFlowTrackSimple* pTrack = NULL;
222 Int_t iNumberOfTracks = anEvent->NumberOfTracks();
223 for (Int_t i=0;i<iNumberOfTracks;i++)
8d312f00 224 {
e35ddff0 225 pTrack = anEvent->GetTrack(i) ;
226 if (pTrack){
e35ddff0 227 Double_t dPhi = pTrack->Phi();
29195b69 228 //set default weight to 1
229 Double_t dW = 1.;
230 //phi weight of pTrack
231 if(fUsePhiWeights && fPhiWeights) {
232 Int_t iNBinsPhi = fPhiWeights->GetNbinsX();
233 dW = fPhiWeights->GetBinContent(1+(Int_t)(TMath::Floor(dPhi*iNBinsPhi/TMath::TwoPi()))); //bin = 1 + value*nbins/range
234 } //TMath::Floor rounds to the lower integer
e35ddff0 235 //calculate vU
236 TVector2 vU;
237 Double_t dUX = TMath::Cos(2*dPhi);
238 Double_t dUY = TMath::Sin(2*dPhi);
239 vU.Set(dUX,dUY);
240 Double_t dModulus = vU.Mod();
241 if (dModulus!=0.) vU.Set(dUX/dModulus,dUY/dModulus); // make length 1
242 else cerr<<"dModulus is zero!"<<endl;
243
244 TVector2 vQm = vQ;
cbbaf54a 245 //subtract particle from the flowvector if used to define it
6a1a854d 246 if (pTrack->InRPSelection()) {
b125a454 247 if (pTrack->InSubevent(0) || pTrack->InSubevent(1)) {
29195b69 248 Double_t dQmX = vQm.X() - dW*dUX;
249 Double_t dQmY = vQm.Y() - dW*dUY;
b125a454 250 vQm.Set(dQmX,dQmY);
251 }
8d312f00 252 }
253
e35ddff0 254 //dUQ = scalar product of vU and vQm
29195b69 255 Double_t dUQ = (vU * vQm)*dMmin1;
e35ddff0 256 Double_t dPt = pTrack->Pt();
395fadba 257 Double_t dEta = pTrack->Eta();
258 //fill the profile histograms
b7cb54d5 259 if (pTrack->InRPSelection()) {
395fadba 260 fHistProUQetaRP -> Fill(dEta,dUQ);
261 fHistProUQPtRP -> Fill(dPt,dUQ);
395fadba 262 }
b7cb54d5 263 if (pTrack->InPOISelection()) {
395fadba 264 fHistProUQetaPOI -> Fill(dEta,dUQ);
265 fHistProUQPtPOI -> Fill(dPt,dUQ);
8d312f00 266 }
267 }//track selected
268 }//loop over tracks
d7eb18ec 269
8d312f00 270 fEventNumber++;
deebd72f 271 // cout<<"@@@@@ "<<fEventNumber<<" events processed"<<endl;
b125a454 272 delete [] vQarray;
8d312f00 273 }
274}
275
fd46c3dd 276 //--------------------------------------------------------------------
277void AliFlowAnalysisWithScalarProduct::GetOutputHistograms(TList *outputListHistos){
278
279 //get pointers to all output histograms (called before Finish())
280 if (outputListHistos) {
281 //Get the common histograms from the output list
282 AliFlowCommonHist *pCommonHist = dynamic_cast<AliFlowCommonHist*>
283 (outputListHistos->FindObject("AliFlowCommonHistSP"));
284 AliFlowCommonHistResults *pCommonHistResults = dynamic_cast<AliFlowCommonHistResults*>
285 (outputListHistos->FindObject("AliFlowCommonHistResultsSP"));
286 TProfile* pHistProQaQb = dynamic_cast<TProfile*>(outputListHistos->FindObject("Flow_QaQb_SP"));
287 TProfile* pHistProM = dynamic_cast<TProfile*>(outputListHistos->FindObject("Flow_M_SP"));
288 TProfile* pHistProUQetaRP = dynamic_cast<TProfile*>(outputListHistos->FindObject("Flow_UQetaRP_SP"));
289 TProfile* pHistProUQetaPOI = dynamic_cast<TProfile*>(outputListHistos->FindObject("Flow_UQetaPOI_SP"));
290 TProfile* pHistProUQPtRP = dynamic_cast<TProfile*>(outputListHistos->FindObject("Flow_UQPtRP_SP"));
291 TProfile* pHistProUQPtPOI = dynamic_cast<TProfile*>(outputListHistos->FindObject("Flow_UQPtPOI_SP"));
292 if (pCommonHist && pCommonHistResults && pHistProQaQb && pHistProM &&
293 pHistProUQetaRP && pHistProUQetaPOI && pHistProUQPtRP && pHistProUQPtPOI) {
294 this -> SetCommonHists(pCommonHist);
295 this -> SetCommonHistsRes(pCommonHistResults);
296 this -> SetHistProQaQb(pHistProQaQb);
297 this -> SetHistProM(pHistProM);
298 this -> SetHistProUQetaRP(pHistProUQetaRP);
299 this -> SetHistProUQetaPOI(pHistProUQetaPOI);
300 this -> SetHistProUQPtRP(pHistProUQPtRP);
301 this -> SetHistProUQPtPOI(pHistProUQPtPOI);
302 }
303 }
304}
305
306//--------------------------------------------------------------------
8d312f00 307void AliFlowAnalysisWithScalarProduct::Finish() {
308
e4cecfa0 309 //calculate flow and fill the AliFlowCommonHistResults
395fadba 310 if (fDebug) cout<<"AliFlowAnalysisWithScalarProduct::Finish()"<<endl;
311
312 cout<<"*************************************"<<endl;
313 cout<<"*************************************"<<endl;
314 cout<<" Integrated flow from "<<endl;
315 cout<<" Scalar product "<<endl;
316 cout<<endl;
e4cecfa0 317
bee2efdc 318 Int_t iNbinsPt = AliFlowCommonConstants::GetMaster()->GetNbinsPt();
319 Int_t iNbinsEta = AliFlowCommonConstants::GetMaster()->GetNbinsEta();
80f72ed6 320
321 Double_t dMmin1 = fHistProM->GetBinContent(1); //average over M-1
322 Double_t dMmin1Err = fHistProM->GetBinError(1); //error on average over M-1
c09d90fd 323 Double_t dMaMb = fHistProM->GetBinContent(2); //average over Ma*Mb
324 Double_t dMaMbErr = fHistProM->GetBinError(2); //error on average over Ma*Mb
80f72ed6 325
29195b69 326 //Double_t dSumMmin1 = fHistM->GetBinContent(1); //sum over (M-1)
327 //Double_t dSumMaMb = fHistM->GetBinContent(2); //sum over (Ma*Mb)
328
80f72ed6 329 Double_t dMcorrection = 0.; //correction factor for Ma != Mb
330 Double_t dMcorrectionErr = 0.;
331 Double_t dMcorrectionErrRel = 0.;
332 Double_t dMcorrectionErrRel2 = 0.;
333
334 if (dMaMb != 0. && dMmin1 != 0.) {
335 dMcorrection = dMmin1/(TMath::Sqrt(dMaMb));
336 dMcorrectionErr = dMcorrection*(dMmin1Err/dMmin1 + dMaMbErr/(2*dMaMb));
337 dMcorrectionErrRel = dMcorrectionErr/dMcorrection;
338 dMcorrectionErrRel2 = dMcorrectionErrRel*dMcorrectionErrRel;
339 }
340
c09d90fd 341 Double_t dQaQbAv = TMath::Abs(fHistProQaQb->GetBinContent(1)); //average over events //TEST TAKE ABS
29195b69 342 dQaQbAv /= dMaMb; //normalize for weighted average
80f72ed6 343 Double_t dQaQbErr = fHistProQaQb->GetBinError(1);
29195b69 344 dQaQbErr /= dMaMb; //normalize for weighted average
80f72ed6 345 Double_t dQaQbErrRel = 0.;
346 if (dQaQbAv != 0.) {
347 dQaQbErrRel = dQaQbErr/dQaQbAv; }
348 Double_t dQaQbErrRel2 = dQaQbErrRel*dQaQbErrRel;
349
e4cecfa0 350 if (dQaQbAv <= 0.){
351 //set v to -0
352 fCommonHistsRes->FillIntegratedFlowRP(-0.,0.);
353 fCommonHistsRes->FillIntegratedFlow(-0.,0.);
354 cout<<"dV(RP) = -0. +- 0."<<endl;
355 fCommonHistsRes->FillIntegratedFlowPOI(-0.,0.);
356 cout<<"dV(POI) = -0. +- 0."<<endl;
357 } else {
c09d90fd 358 Double_t dQaQbSqrt = TMath::Sqrt(dQaQbAv); //DOES NOT WORK IF dQaQbAv IS NEGATIVE
80f72ed6 359 if (dMaMb>0.) { dQaQbSqrt *= dMcorrection; }
92041674 360 else { dQaQbSqrt = 0.; }
80f72ed6 361 Double_t dQaQbSqrtErrRel2 = dMcorrectionErrRel2 + (1/4)*dQaQbErrRel2;
362
e4cecfa0 363 //v as a function of eta for RP selection
364 for(Int_t b=0;b<iNbinsEta;b++) {
365 Double_t duQpro = fHistProUQetaRP->GetBinContent(b);
29195b69 366 duQpro /= dMmin1; //normalize for weighted average
e4cecfa0 367 Double_t duQerr = fHistProUQetaRP->GetBinError(b); //copy error for now
29195b69 368 duQerr /= dMmin1; //normalize for weighted average
e4cecfa0 369 Double_t duQerrRel = 0.;
370 if (duQpro != 0.) {duQerrRel = duQerr/duQpro;}
371 Double_t duQerrRel2 = duQerrRel*duQerrRel;
372
92041674 373 Double_t dv2pro = 0.;
374 if (dQaQbSqrt!=0.) { dv2pro = duQpro/dQaQbSqrt; }
e4cecfa0 375 Double_t dv2errRel2 = duQerrRel2 + dQaQbSqrtErrRel2;
92041674 376 Double_t dv2errRel = 0.;
377 if (dv2errRel2>0.) { dv2errRel = TMath::Sqrt(dv2errRel2); }
e4cecfa0 378 Double_t dv2err = dv2pro*dv2errRel;
379 //fill TH1D
380 fCommonHistsRes->FillDifferentialFlowEtaRP(b, dv2pro, dv2err);
381 } //loop over bins b
395fadba 382
e4cecfa0 383 //v as a function of eta for POI selection
384 for(Int_t b=0;b<iNbinsEta;b++) {
385 Double_t duQpro = fHistProUQetaPOI->GetBinContent(b);
29195b69 386 duQpro /= dMmin1; //normalize for weighted average
e4cecfa0 387 Double_t duQerr = fHistProUQetaPOI->GetBinError(b); //copy error for now
29195b69 388 duQerr /= dMmin1; //normalize for weighted average
e4cecfa0 389 Double_t duQerrRel = 0.;
390 if (duQpro != 0.) {duQerrRel = duQerr/duQpro;}
391 Double_t duQerrRel2 = duQerrRel*duQerrRel;
392
92041674 393 Double_t dv2pro = 0.;
394 if (dQaQbSqrt!=0.) { dv2pro = duQpro/dQaQbSqrt; }
e4cecfa0 395 Double_t dv2errRel2 = duQerrRel2 + dQaQbSqrtErrRel2;
92041674 396 Double_t dv2errRel = 0.;
397 if (dv2errRel2>0.) { dv2errRel = TMath::Sqrt(dv2errRel2); }
e4cecfa0 398 Double_t dv2err = dv2pro*dv2errRel;
399 //fill TH1D
400 fCommonHistsRes->FillDifferentialFlowEtaPOI(b, dv2pro, dv2err);
401 } //loop over bins b
395fadba 402
e4cecfa0 403 //v as a function of Pt for RP selection
b7cb54d5 404 TH1F* fHistPtRP = fCommonHists->GetHistPtRP(); //for calculating integrated flow
e4cecfa0 405 Double_t dVRP = 0.;
406 Double_t dSum = 0.;
407 Double_t dErrV =0.;
408
409 for(Int_t b=0;b<iNbinsPt;b++) {
410 Double_t duQpro = fHistProUQPtRP->GetBinContent(b);
29195b69 411 duQpro /= dMmin1; //normalize for weighted average
e4cecfa0 412 Double_t duQerr = fHistProUQPtRP->GetBinError(b); //copy error for now
29195b69 413 duQerr /= dMmin1; //normalize for weighted average
e4cecfa0 414 Double_t duQerrRel = 0.;
415 if (duQpro != 0.) {duQerrRel = duQerr/duQpro;}
416 Double_t duQerrRel2 = duQerrRel*duQerrRel;
417
92041674 418 Double_t dv2pro = 0.;
419 if (dQaQbSqrt!=0.) { dv2pro = duQpro/dQaQbSqrt; }
e4cecfa0 420 Double_t dv2errRel2 = duQerrRel2 + dQaQbSqrtErrRel2;
92041674 421 Double_t dv2errRel = 0.;
422 if (dv2errRel2>0.) { dv2errRel = TMath::Sqrt(dv2errRel2); }
e4cecfa0 423 Double_t dv2err = dv2pro*dv2errRel;
424 //fill TH1D
425 fCommonHistsRes->FillDifferentialFlowPtRP(b, dv2pro, dv2err);
426 //calculate integrated flow for RP selection
b7cb54d5 427 if (fHistPtRP){
428 Double_t dYieldPt = fHistPtRP->GetBinContent(b);
e4cecfa0 429 dVRP += dv2pro*dYieldPt;
430 dSum +=dYieldPt;
431 dErrV += dYieldPt*dYieldPt*dv2err*dv2err;
b7cb54d5 432 } else { cout<<"fHistPtRP is NULL"<<endl; }
e4cecfa0 433 } //loop over bins b
434
435 if (dSum != 0.) {
436 dVRP /= dSum; //the pt distribution should be normalised
437 dErrV /= (dSum*dSum);
438 dErrV = TMath::Sqrt(dErrV);
439 }
440 fCommonHistsRes->FillIntegratedFlowRP(dVRP,dErrV);
441 fCommonHistsRes->FillIntegratedFlow(dVRP,dErrV);
442
443 cout<<"dV(RP) = "<<dVRP<<" +- "<<dErrV<<endl;
395fadba 444
e4cecfa0 445 //v as a function of Pt for POI selection
b7cb54d5 446 TH1F* fHistPtPOI = fCommonHists->GetHistPtPOI(); //for calculating integrated flow
e4cecfa0 447 Double_t dVPOI = 0.;
448 dSum = 0.;
449 dErrV =0.;
395fadba 450
e4cecfa0 451 for(Int_t b=0;b<iNbinsPt;b++) {
452 Double_t duQpro = fHistProUQPtPOI->GetBinContent(b);
29195b69 453 duQpro /= dMmin1; //normalize for weighted average
e4cecfa0 454 Double_t duQerr = fHistProUQPtPOI->GetBinError(b); //copy error for now
29195b69 455 duQerr /= dMmin1; //normalize for weighted average
e4cecfa0 456 Double_t duQerrRel = 0.;
457 if (duQpro != 0.) {duQerrRel = duQerr/duQpro;}
458 Double_t duQerrRel2 = duQerrRel*duQerrRel;
459
92041674 460 Double_t dv2pro = 0.;
461 if (dQaQbSqrt!=0.) { dv2pro = duQpro/dQaQbSqrt; }
e4cecfa0 462 Double_t dv2errRel2 = duQerrRel2 + dQaQbSqrtErrRel2;
92041674 463 Double_t dv2errRel = 0.;
464 if (dv2errRel2>0.) { dv2errRel = TMath::Sqrt(dv2errRel2); }
e4cecfa0 465 Double_t dv2err = dv2pro*dv2errRel;
466 //fill TH1D
467 fCommonHistsRes->FillDifferentialFlowPtPOI(b, dv2pro, dv2err);
468
469 //calculate integrated flow for POI selection
b7cb54d5 470 if (fHistPtPOI){
471 Double_t dYieldPt = fHistPtPOI->GetBinContent(b);
e4cecfa0 472 dVPOI += dv2pro*dYieldPt;
473 dSum +=dYieldPt;
474 dErrV += dYieldPt*dYieldPt*dv2err*dv2err;
b7cb54d5 475 } else { cout<<"fHistPtPOI is NULL"<<endl; }
e4cecfa0 476 } //loop over bins b
477
478 if (dSum != 0.) {
479 dVPOI /= dSum; //the pt distribution should be normalised
480 dErrV /= (dSum*dSum);
481 dErrV = TMath::Sqrt(dErrV);
482 }
483 fCommonHistsRes->FillIntegratedFlowPOI(dVPOI,dErrV);
484
485 cout<<"dV(POI) = "<<dVPOI<<" +- "<<dErrV<<endl;
c09d90fd 486 }
395fadba 487 cout<<endl;
488 cout<<"*************************************"<<endl;
489 cout<<"*************************************"<<endl;
8d312f00 490
395fadba 491 //cout<<".....finished"<<endl;
8d312f00 492 }
8232a5ec 493
494