AliPhysics  ff07904 (ff07904)
AliFlowAnalysisWithLYZEventPlane.cxx
Go to the documentation of this file.
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 // AliFlowAnalysisWithLYZEventPlane:
17 //
18 // Class to do flow analysis with the event plane from the LYZ method
19 //
20 // author: N. van der Kolk (kolk@nikhef.nl)
21 
22 /*
23 $Log$
24 */
25 
26 //#define AliFlowAnalysisWithLYZEventPlane_cxx
27 
28 #include "Riostream.h" //needed as include
29 #include "TMath.h" //needed as include
30 #include "TProfile.h" //needed as include
31 #include "TH1F.h"
32 #include "TFile.h"
33 #include "TList.h"
34 #include "TVector2.h"
35 
36 #include "AliFlowLYZConstants.h" //needed as include
37 #include "AliFlowCommonConstants.h" //needed as include
38 #include "AliFlowEventSimple.h"
39 #include "AliFlowTrackSimple.h"
40 #include "AliFlowCommonHist.h"
42 #include "AliFlowLYZEventPlane.h"
44 #include "AliFlowVector.h"
45 
46 using std::endl;
47 using std::cout;
48 using std::cerr;
50 
51  //-----------------------------------------------------------------------
52 
54  fHistList(NULL),
55  fSecondRunList(NULL),
56  fSecondReDtheta(NULL),
57  fSecondImDtheta(NULL),
58  fFirstr0theta(NULL),
59  fHistProVetaRP(NULL),
60  fHistProVetaPOI(NULL),
61  fHistProVPtRP(NULL),
62  fHistProVPtPOI(NULL),
63  fHistProWr(NULL),
64  fHistProWrCorr(NULL),
65  fHistQsumforChi(NULL),
66  fHistDeltaPhi(NULL),
67  fHistDeltaPhi2(NULL),
68  fHistDeltaPhihere(NULL),
69  fHistPhiEP(NULL),
70  fHistPhiEPhere(NULL),
71  fHistPhiLYZ(NULL),
72  fHistPhiLYZ2(NULL),
73  fCommonHists(NULL),
74  fCommonHistsRes(NULL),
75  fEventNumber(0),
76  fQsum(NULL),
77  fQ2sum(0)
78 {
79 
80  // Constructor.
81  fQsum = new TVector2(); // flow vector sum
82 
83  fHistList = new TList();
84  fSecondRunList = new TList();
85 }
86 
87 
88 
89  //-----------------------------------------------------------------------
90 
91 
93  {
94  //destructor
95  delete fQsum;
96  delete fHistList;
97  delete fSecondRunList;
98  }
99 
100 
101 //-----------------------------------------------------------------------
102 
104 {
105  //store the final results in output .root file
106 
107  TFile *output = new TFile(outputFileName->Data(),"RECREATE");
108  //output->WriteObject(fHistList, "cobjLYZEP","SingleKey");
109  fHistList->SetName("cobjLYZEP");
110  fHistList->SetOwner(kTRUE);
111  fHistList->Write(fHistList->GetName(), TObject::kSingleKey);
112  delete output;
113 }
114 
115 //-----------------------------------------------------------------------
116 
118 {
119  //store the final results in output .root file
120 
121  TFile *output = new TFile(outputFileName.Data(),"RECREATE");
122  //output->WriteObject(fHistList, "cobjLYZEP","SingleKey");
123  fHistList->SetName("cobjLYZEP");
124  fHistList->SetOwner(kTRUE);
125  fHistList->Write(fHistList->GetName(), TObject::kSingleKey);
126  delete output;
127 }
128 
129 //-----------------------------------------------------------------------
130 
131 void AliFlowAnalysisWithLYZEventPlane::WriteHistograms(TDirectoryFile *outputFileName)
132 {
133  //store the final results in output .root file
134  fHistList->SetName("cobjLYZEP");
135  fHistList->SetOwner(kTRUE);
136  outputFileName->Add(fHistList);
137  outputFileName->Write(outputFileName->GetName(), TObject::kSingleKey);
138 }
139 
140 //-----------------------------------------------------------------------
141 
143 
144  //Initialise all histograms
145  cout<<"---Analysis with Lee Yang Zeros Event Plane Method---"<<endl;
146 
147  //save old value and prevent histograms from being added to directory
148  //to avoid name clashes in case multiple analaysis objects are used
149  //in an analysis
150  Bool_t oldHistAddStatus = TH1::AddDirectoryStatus();
151  TH1::AddDirectory(kFALSE);
152 
153  //input histograms
154  if (fSecondRunList) {
155  fSecondReDtheta = (TProfile*)fSecondRunList->FindObject("Second_FlowPro_ReDtheta_LYZSUM");
157 
158  fSecondImDtheta = (TProfile*)fSecondRunList->FindObject("Second_FlowPro_ImDtheta_LYZSUM");
160 
161  fFirstr0theta = (TProfile*)fSecondRunList->FindObject("First_FlowPro_r0theta_LYZSUM");
162  fHistList->Add(fFirstr0theta);
163 
164  //warnings
165  if (!fSecondReDtheta) {cout<<"fSecondReDtheta is NULL!"<<endl; }
166  if (!fSecondImDtheta) {cout<<"fSecondImDtheta is NULL!"<<endl; }
167  if (!fFirstr0theta) {cout<<"fFirstr0theta is NULL!"<<endl; }
168 
169  }
170 
171  fCommonHists = new AliFlowCommonHist("AliFlowCommonHistLYZEP");
172  fHistList->Add(fCommonHists);
173 
174  fCommonHistsRes = new AliFlowCommonHistResults("AliFlowCommonHistResultsLYZEP");
175  fHistList->Add(fCommonHistsRes);
176 
183 
184  fHistProVetaRP = new TProfile("FlowPro_VetaRP_LYZEP","FlowPro_VetaRP_LYZEP",iNbinsEta,dEtaMin,dEtaMax);
185  fHistProVetaRP->SetXTitle("rapidity");
186  fHistProVetaRP->SetYTitle("v_{2}(#eta) for RP selection");
188 
189  fHistProVetaPOI = new TProfile("FlowPro_VetaPOI_LYZEP","FlowPro_VetaPOI_LYZEP",iNbinsEta,dEtaMin,dEtaMax);
190  fHistProVetaPOI->SetXTitle("rapidity");
191  fHistProVetaPOI->SetYTitle("v_{2}(#eta) for POI selection");
193 
194  fHistProVPtRP = new TProfile("FlowPro_VPtRP_LYZEP","FlowPro_VPtRP_LYZEP",iNbinsPt,dPtMin,dPtMax);
195  fHistProVPtRP->SetXTitle("Pt");
196  fHistProVPtRP->SetYTitle("v_{2}(p_{T}) for RP selection");
197  fHistList->Add(fHistProVPtRP);
198 
199  fHistProVPtPOI = new TProfile("FlowPro_VPtPOI_LYZEP","FlowPro_VPtPOI_LYZEP",iNbinsPt,dPtMin,dPtMax);
200  fHistProVPtPOI->SetXTitle("p_{T}");
201  fHistProVPtPOI->SetYTitle("v_{2}(p_{T}) for POI selection");
203 
204  fHistProWr = new TProfile("FlowPro_Wr_LYZEP","FlowPro_Wr_LYZEP",100,0.,0.25);
205  fHistProWr->SetXTitle("Q");
206  fHistProWr->SetYTitle("Wr");
207  fHistList->Add(fHistProWr);
208 
209  fHistQsumforChi = new TH1F("Flow_QsumforChi_LYZEP","Flow_QsumforChi_LYZEP",3,-1.,2.);
210  fHistQsumforChi->SetXTitle("Qsum.X , Qsum.Y, Q2sum");
211  fHistQsumforChi->SetYTitle("value");
213 
214  fHistDeltaPhi = new TH1F("Flow_DeltaPhi_LYZEP","Flow_DeltaPhi_LYZEP",100,0.,3.14);
215  fHistDeltaPhi->SetXTitle("DeltaPhi");
216  fHistDeltaPhi->SetYTitle("Counts");
217  fHistList->Add(fHistDeltaPhi);
218 
219  fHistPhiLYZ = new TH1F("Flow_PhiLYZ_LYZEP","Flow_PhiLYZ_LYZEP",100,0.,3.14);
220  fHistPhiLYZ->SetXTitle("Phi from LYZ");
221  fHistPhiLYZ->SetYTitle("Counts");
222  fHistList->Add(fHistPhiLYZ);
223 
224  fHistPhiEP = new TH1F("Flow_PhiEP_LYZEP","Flow_PhiEP_LYZEP",100,0.,3.14);
225  fHistPhiEP->SetXTitle("Phi from EP");
226  fHistPhiEP->SetYTitle("Counts");
227  fHistList->Add(fHistPhiEP);
228 
229  fEventNumber = 0; //set number of events to zero
230 
231  //restore old status
232  TH1::AddDirectory(oldHistAddStatus);
233 }
234 
235 //-----------------------------------------------------------------------
236 
238 
239  //Get the event plane and weight for each event
240  if (anEvent) {
241 
242  //fill control histograms
244 
245  //get the Q vector from the FlowEvent
246  AliFlowVector vQ = anEvent->GetQ();
247  //if (vQ.X()== 0. && vQ.Y()== 0. ) { cout<<"Q vector is NULL!"<<endl; } //coding violation
248  //Weight with the multiplicity
249  Double_t dQX = 0.;
250  Double_t dQY = 0.;
251  if (TMath::AreEqualAbs(vQ.GetMult(),0.0,1e-10)) {
252  dQX = vQ.X()/vQ.GetMult();
253  dQY = vQ.Y()/vQ.GetMult();
254  } else {cerr<<"vQ.GetMult() is zero!"<<endl; }
255  vQ.Set(dQX,dQY);
256  //cout<<"vQ("<<dQX<<","<<dQY<<")"<<endl;
257 
258  //for chi calculation:
259  *fQsum += vQ;
260  fHistQsumforChi->SetBinContent(1,fQsum->X());
261  fHistQsumforChi->SetBinContent(2,fQsum->Y());
262  fQ2sum += vQ.Mod2();
263  fHistQsumforChi->SetBinContent(3,fQ2sum);
264  //cout<<"fQ2sum = "<<fQ2sum<<endl;
265 
266  //call AliFlowLYZEventPlane::CalculateRPandW() here!
267  aLYZEP->CalculateRPandW(vQ);
268 
269  Double_t dWR = aLYZEP->GetWR();
270  Double_t dRP = aLYZEP->GetPsi();
271 
272  //fHistProWr->Fill(vQ.Mod(),dWR); //this weight is always positive
273  fHistPhiLYZ->Fill(dRP);
274 
275  //plot difference between event plane from EP-method and LYZ-method
276  Double_t dRPEP = vQ.Phi()/2; //gives distribution from (0 to pi)
277  //Double_t dRPEP = 0.5*TMath::ATan2(vQ.Y(),vQ.X()); //gives distribution from (-pi/2 to pi/2)
278  //cout<<"dRPEP = "<< dRPEP <<endl;
279  fHistPhiEP->Fill(dRPEP);
280 
281  Double_t dDeltaPhi = dRPEP - dRP;
282  if (dDeltaPhi < 0.) { dDeltaPhi += TMath::Pi(); } //to shift distribution from (-pi/2 to pi/2) to (0 to pi)
283  //cout<<"dDeltaPhi = "<<dDeltaPhi<<endl;
284  fHistDeltaPhi->Fill(dDeltaPhi);
285 
286  //Flip sign of WR
287  Double_t dLow = TMath::Pi()/4.;
288  Double_t dHigh = 3.*(TMath::Pi()/4.);
289  if ((dDeltaPhi > dLow) && (dDeltaPhi < dHigh)){
290  dRP -= (TMath::Pi()/2);
291  dWR = -dWR;
292  cerr<<"*** dRP modified ***"<<endl;
293  }
294  fHistProWr->Fill(vQ.Mod(),dWR); //corrected weight
295 
296  //calculate flow for RP and POI selections
297  //loop over the tracks of the event
298  Int_t iNumberOfTracks = anEvent->NumberOfTracks();
299  for (Int_t i=0;i<iNumberOfTracks;i++)
300  {
301  AliFlowTrackSimple* pTrack = anEvent->GetTrack(i) ;
302  if (pTrack){
303  Double_t dPhi = pTrack->Phi();
304  //if (dPhi<0.) fPhi+=2*TMath::Pi();
305  Double_t dPt = pTrack->Pt();
306  Double_t dEta = pTrack->Eta();
307  //calculate flow v2:
308  Double_t dv2 = dWR * TMath::Cos(2*(dPhi-dRP));
309  if (pTrack->InRPSelection()) {
310  //fill histograms for RP selection
311  fHistProVetaRP -> Fill(dEta,dv2);
312  fHistProVPtRP -> Fill(dPt,dv2);
313  }
314  if (pTrack->InPOISelection()) {
315  //fill histograms for POI selection
316  fHistProVetaPOI -> Fill(dEta,dv2);
317  fHistProVPtPOI -> Fill(dPt,dv2);
318  }
319  }//track
320  }//loop over tracks
321 
322  fEventNumber++;
323  // cout<<"@@@@@ "<<fEventNumber<<" events processed"<<endl;
324  }
325 }
326 
327  //--------------------------------------------------------------------
329  //get pointers to all output histograms (called before Finish())
330  if (outputListHistos) {
331  //Get the common histograms from the output list
332  AliFlowCommonHist *pCommonHist = dynamic_cast<AliFlowCommonHist*>
333  (outputListHistos->FindObject("AliFlowCommonHistLYZEP"));
334  AliFlowCommonHistResults *pCommonHistResults = dynamic_cast<AliFlowCommonHistResults*>
335  (outputListHistos->FindObject("AliFlowCommonHistResultsLYZEP"));
336 
337  TProfile* pHistProR0theta = dynamic_cast<TProfile*>
338  (outputListHistos->FindObject("First_FlowPro_r0theta_LYZSUM"));
339 
340  TProfile* pHistProVetaRP = dynamic_cast<TProfile*>
341  (outputListHistos->FindObject("FlowPro_VetaRP_LYZEP"));
342  TProfile* pHistProVetaPOI = dynamic_cast<TProfile*>
343  (outputListHistos->FindObject("FlowPro_VetaPOI_LYZEP"));
344  TProfile* pHistProVPtRP = dynamic_cast<TProfile*>
345  (outputListHistos->FindObject("FlowPro_VPtRP_LYZEP"));
346  TProfile* pHistProVPtPOI = dynamic_cast<TProfile*>
347  (outputListHistos->FindObject("FlowPro_VPtPOI_LYZEP"));
348 
349  TH1F* pHistQsumforChi = dynamic_cast<TH1F*>
350  (outputListHistos->FindObject("Flow_QsumforChi_LYZEP"));
351 
352  if (pCommonHist && pCommonHistResults && pHistProR0theta &&
353  pHistProVetaRP && pHistProVetaPOI && pHistProVPtRP &&
354  pHistProVPtPOI && pHistQsumforChi ) {
355  this -> SetCommonHists(pCommonHist);
356  this -> SetCommonHistsRes(pCommonHistResults);
357  this -> SetFirstr0theta(pHistProR0theta);
358  this -> SetHistProVetaRP(pHistProVetaRP);
359  this -> SetHistProVetaPOI(pHistProVetaPOI);
360  this -> SetHistProVPtRP(pHistProVPtRP);
361  this -> SetHistProVPtPOI(pHistProVPtPOI);
362  this -> SetHistQsumforChi(pHistQsumforChi);
363  }
364  } else {
365  cout<<"WARNING: Histograms needed to run Finish() are not accessable!"<<endl;
366  }
367 }
368 
369  //--------------------------------------------------------------------
371 
372  //plot histograms etc.
373  cout<<"AliFlowAnalysisWithLYZEventPlane::Finish()"<<endl;
374 
375  //constants:
376  Double_t dJ01 = 2.405;
380  //set the event number
381  if (fCommonHists) {
382  SetEventNumber((int)fCommonHists->GetHistQ()->GetEntries());
383  //cout<<"number of events processed is "<<fEventNumber<<endl;
384  } else {
385  cout<<"Commonhist pointer is NULL."<<endl;
386  cout<<"Leaving LYZ Event plane analysis!"<<endl;
387  return;
388  }
389 
390  //set the sum of Q vectors
391  fQsum->Set(fHistQsumforChi->GetBinContent(1),fHistQsumforChi->GetBinContent(2));
392  SetQ2sum(fHistQsumforChi->GetBinContent(3));
393 
394  //calculate dV the mean of dVtheta
395  Double_t dVtheta = 0;
396  Double_t dV = 0;
397  for (Int_t theta=0;theta<iNtheta;theta++) {
398  Double_t dR0 = fFirstr0theta->GetBinContent(theta+1);
399  if (TMath::AreEqualAbs(dR0,0.0,1e-10)) { dVtheta = dJ01/dR0 ;}
400  dV += dVtheta;
401  }
402  dV /= iNtheta;
403 
404  //calculate dChi
405  Double_t dSigma2 = 0;
406  Double_t dChi= 0;
407  if (fEventNumber!=0) {
408  *fQsum /= fEventNumber;
409  //cerr<<"fQsum->X() = "<<fQsum->X()<<endl;
410  //cerr<<"fQsum->Y() = "<<fQsum->Y()<<endl;
411  fQ2sum /= fEventNumber;
412  //cerr<<"fEventNumber = "<<fEventNumber<<endl;
413  //cerr<<"fQ2sum = "<<fQ2sum<<endl;
414  dSigma2 = fQ2sum - TMath::Power(fQsum->X(),2.) - TMath::Power(fQsum->Y(),2.) - TMath::Power(dV,2.); //BP eq. 62
415  //cerr<<"dSigma2"<<dSigma2<<endl;
416  if (dSigma2>0) dChi = dV/TMath::Sqrt(dSigma2);
417  else dChi = -1.;
418  fCommonHistsRes->FillChi(dChi);
419 
420  // recalculate statistical errors on integrated flow
421  //combining 5 theta angles to 1 relative error BP eq. 89
422  Double_t dRelErr2comb = 0.;
423  Int_t iEvts = fEventNumber;
424  if (iEvts!=0) {
425  for (Int_t theta=0;theta<iNtheta;theta++){
426  Double_t dTheta = ((double)theta/iNtheta)*TMath::Pi();
427  Double_t dApluscomb = TMath::Exp((dJ01*dJ01)/(2*dChi*dChi)*
428  TMath::Cos(dTheta));
429  Double_t dAmincomb = TMath::Exp(-(dJ01*dJ01)/(2*dChi*dChi)*
430  TMath::Cos(dTheta));
431  dRelErr2comb += (1/(2*iEvts*(dJ01*dJ01)*TMath::BesselJ1(dJ01)*
432  TMath::BesselJ1(dJ01)))*
433  (dApluscomb*TMath::BesselJ0(2*dJ01*TMath::Sin(dTheta/2)) +
434  dAmincomb*TMath::BesselJ0(2*dJ01*TMath::Cos(dTheta/2)));
435  }
436  dRelErr2comb /= iNtheta;
437  }
438  Double_t dRelErrcomb = TMath::Sqrt(dRelErr2comb);
439  Double_t dVErr = dV*dRelErrcomb ;
440  fCommonHistsRes->FillIntegratedFlow(dV, dVErr);
441 
442  cout<<"*************************************"<<endl;
443  cout<<"*************************************"<<endl;
444  cout<<" Integrated flow from "<<endl;
445  cout<<" Lee-Yang Zeroes Event Plane "<<endl;
446  cout<<endl;
447  cout<<"dChi = "<<dChi<<endl;
448  cout<<"dV = "<<dV<<" +- "<<dVErr<<endl;
449  cout<<endl;
450 
451  }
452 
453  //copy content of profile into TH1D, add error and fill the AliFlowCommonHistResults
454 
455  //v as a function of eta for RP selection
456  for(Int_t b=0;b<iNbinsEta;b++) {
457  Double_t dv2pro = fHistProVetaRP->GetBinContent(b);
459  Double_t dErrdifcomb = 0.; //set error to zero
460  Double_t dErr2difcomb = 0.; //set error to zero
461  //calculate error
462  if (TMath::AreEqualAbs(dNprime,0.0,1e-10)) {
463  for (Int_t theta=0;theta<iNtheta;theta++) {
464  Double_t dTheta = ((double)theta/iNtheta)*TMath::Pi();
465  Double_t dApluscomb = TMath::Exp((dJ01*dJ01)/(2*dChi*dChi)*
466  TMath::Cos(dTheta));
467  Double_t dAmincomb = TMath::Exp(-(dJ01*dJ01)/(2*dChi*dChi)*
468  TMath::Cos(dTheta));
469  dErr2difcomb += (TMath::Cos(dTheta)/(4*dNprime*TMath::BesselJ1(dJ01)*
470  TMath::BesselJ1(dJ01)))*
471  ((dApluscomb*TMath::BesselJ0(2*dJ01*TMath::Sin(dTheta/2))) -
472  (dAmincomb*TMath::BesselJ0(2*dJ01*TMath::Cos(dTheta/2))));
473  } //loop over theta
474  }
475 
476  if (TMath::AreEqualAbs(dErr2difcomb, 0.0, 1e-10)) {
477  dErr2difcomb /= iNtheta;
478  dErrdifcomb = TMath::Sqrt(dErr2difcomb);
479  }
480  else {dErrdifcomb = 0.;}
481  //fill TH1D
482  fCommonHistsRes->FillDifferentialFlowEtaRP(b, dv2pro, dErrdifcomb);
483  } //loop over bins b
484 
485 
486  //v as a function of eta for POI selection
487  for(Int_t b=0;b<iNbinsEta;b++) {
488  Double_t dv2pro = fHistProVetaPOI->GetBinContent(b);
490  Double_t dErrdifcomb = 0.; //set error to zero
491  Double_t dErr2difcomb = 0.; //set error to zero
492  //calculate error
493  if (dNprime!=0.) {
494  for (Int_t theta=0;theta<iNtheta;theta++) {
495  Double_t dTheta = ((double)theta/iNtheta)*TMath::Pi();
496  Double_t dApluscomb = TMath::Exp((dJ01*dJ01)/(2*dChi*dChi)*
497  TMath::Cos(dTheta));
498  Double_t dAmincomb = TMath::Exp(-(dJ01*dJ01)/(2*dChi*dChi)*
499  TMath::Cos(dTheta));
500  dErr2difcomb += (TMath::Cos(dTheta)/(4*dNprime*TMath::BesselJ1(dJ01)*
501  TMath::BesselJ1(dJ01)))*
502  ((dApluscomb*TMath::BesselJ0(2*dJ01*TMath::Sin(dTheta/2))) -
503  (dAmincomb*TMath::BesselJ0(2*dJ01*TMath::Cos(dTheta/2))));
504  } //loop over theta
505  }
506 
507  if (dErr2difcomb!=0.) {
508  dErr2difcomb /= iNtheta;
509  dErrdifcomb = TMath::Sqrt(dErr2difcomb);
510  }
511  else {dErrdifcomb = 0.;}
512  //fill TH1D
513  fCommonHistsRes->FillDifferentialFlowEtaPOI(b, dv2pro, dErrdifcomb);
514  } //loop over bins b
515 
516  //v as a function of Pt for RP selection
517  TH1F* fHistPtRP = fCommonHists->GetHistPtRP(); //for calculating integrated flow
518  Double_t dVRP = 0.;
519  Double_t dSum = 0.;
520  Double_t dErrV =0.;
521 
522  for(Int_t b=0;b<iNbinsPt;b++) {
523  Double_t dv2pro = fHistProVPtRP->GetBinContent(b);
525  Double_t dErrdifcomb = 0.; //set error to zero
526  Double_t dErr2difcomb = 0.; //set error to zero
527  //calculate error
528  if (dNprime!=0.) {
529  for (Int_t theta=0;theta<iNtheta;theta++) {
530  Double_t dTheta = ((double)theta/iNtheta)*TMath::Pi();
531  Double_t dApluscomb = TMath::Exp((dJ01*dJ01)/(2*dChi*dChi)*
532  TMath::Cos(dTheta));
533  Double_t dAmincomb = TMath::Exp(-(dJ01*dJ01)/(2*dChi*dChi)*
534  TMath::Cos(dTheta));
535  dErr2difcomb += (TMath::Cos(dTheta)/(4*dNprime*TMath::BesselJ1(dJ01)*
536  TMath::BesselJ1(dJ01)))*
537  ((dApluscomb*TMath::BesselJ0(2*dJ01*TMath::Sin(dTheta/2))) -
538  (dAmincomb*TMath::BesselJ0(2*dJ01*TMath::Cos(dTheta/2))));
539  } //loop over theta
540  }
541 
542  if (dErr2difcomb!=0.) {
543  dErr2difcomb /= iNtheta;
544  dErrdifcomb = TMath::Sqrt(dErr2difcomb);
545  //cerr<<"dErrdifcomb = "<<dErrdifcomb<<endl;
546  }
547  else {dErrdifcomb = 0.;}
548 
549  //fill TH1D
550  fCommonHistsRes->FillDifferentialFlowPtRP(b, dv2pro, dErrdifcomb);
551  //calculate integrated flow for RP selection
552  if (fHistPtRP){
553  Double_t dYieldPt = fHistPtRP->GetBinContent(b);
554  dVRP += dv2pro*dYieldPt;
555  dSum +=dYieldPt;
556  dErrV += dYieldPt*dYieldPt*dErrdifcomb*dErrdifcomb;
557  } else { cout<<"fHistPtRP is NULL"<<endl; }
558 
559  } //loop over bins b
560 
561  if (TMath::AreEqualAbs(dSum, 0.0, 1e-10)) {
562  dVRP /= dSum; //the pt distribution should be normalised
563  dErrV /= (dSum*dSum);
564  dErrV = TMath::Sqrt(dErrV);
565  }
567 
568  cout<<"dV(RP) = "<<dVRP<<" +- "<<dErrV<<endl;
569  //cout<<endl;
570 
571 
572  //v as a function of Pt for POI selection
573  TH1F* fHistPtPOI = fCommonHists->GetHistPtPOI(); //for calculating integrated flow
574  Double_t dVPOI = 0.;
575  dSum = 0.;
576  dErrV =0.;
577 
578  for(Int_t b=0;b<iNbinsPt;b++) {
579  Double_t dv2pro = fHistProVPtPOI->GetBinContent(b);
581 
582  //cerr<<"dNprime = "<<dNprime<<endl;
583  //Int_t iNprime = TMath::Nint(fHistProVPtPOI->GetBinEntries(b));
584  //cerr<<"iNprime = "<<iNprime<<endl;
585 
586  Double_t dErrdifcomb = 0.; //set error to zero
587  Double_t dErr2difcomb = 0.; //set error to zero
588  //calculate error
589  if (dNprime!=0.) {
590  for (Int_t theta=0;theta<iNtheta;theta++) {
591  Double_t dTheta = ((double)theta/iNtheta)*TMath::Pi();
592  Double_t dApluscomb = TMath::Exp((dJ01*dJ01)/(2*dChi*dChi)*
593  TMath::Cos(dTheta));
594  Double_t dAmincomb = TMath::Exp(-(dJ01*dJ01)/(2*dChi*dChi)*
595  TMath::Cos(dTheta));
596  dErr2difcomb += (TMath::Cos(dTheta)/(4*dNprime*TMath::BesselJ1(dJ01)*
597  TMath::BesselJ1(dJ01)))*
598  ((dApluscomb*TMath::BesselJ0(2*dJ01*TMath::Sin(dTheta/2))) -
599  (dAmincomb*TMath::BesselJ0(2*dJ01*TMath::Cos(dTheta/2))));
600  } //loop over theta
601  }
602 
603  if (dErr2difcomb!=0.) {
604  dErr2difcomb /= iNtheta;
605  dErrdifcomb = TMath::Sqrt(dErr2difcomb);
606  //cerr<<"dErrdifcomb = "<<dErrdifcomb<<endl;
607  }
608  else {dErrdifcomb = 0.;}
609 
610  //fill TH1D
611  fCommonHistsRes->FillDifferentialFlowPtPOI(b, dv2pro, dErrdifcomb);
612 
613  //calculate integrated flow for POI selection
614  if (fHistPtPOI){
615  Double_t dYieldPt = fHistPtPOI->GetBinContent(b);
616  dVPOI += dv2pro*dYieldPt;
617  dSum +=dYieldPt;
618  dErrV += dYieldPt*dYieldPt*dErrdifcomb*dErrdifcomb;
619  } else { cout<<"fHistPtPOI is NULL"<<endl; }
620  } //loop over bins b
621 
622  if (dSum != 0.) {
623  dVPOI /= dSum; //the pt distribution should be normalised
624  dErrV /= (dSum*dSum);
625  dErrV = TMath::Sqrt(dErrV);
626  }
628 
629  cout<<"dV(POI) = "<<dVPOI<<" +- "<<dErrV<<endl;
630  cout<<endl;
631  cout<<"*************************************"<<endl;
632  cout<<"*************************************"<<endl;
633 
634  //cout<<".....finished"<<endl;
635  }
636 
void SetCommonHistsRes(AliFlowCommonHistResults *const aCommonHistResult)
double Double_t
Definition: External.C:58
virtual AliFlowVector GetQ(Int_t n=2, TList *weightsList=NULL, Bool_t usePhiWeights=kFALSE, Bool_t usePtWeights=kFALSE, Bool_t useEtaWeights=kFALSE)
Double_t GetEntriesInEtaBinPOI(Int_t iBin)
Int_t GetNtheta() const
AliFlowTrackSimple * GetTrack(Int_t i)
virtual void Make(AliFlowEventSimple *fEvent, AliFlowLYZEventPlane *fLYZEP)
Double_t GetEntriesInPtBinPOI(Int_t iBin)
void SetCommonHists(AliFlowCommonHist *const aCommonHist)
Double_t GetMult() const
Definition: AliFlowVector.h:46
Bool_t FillDifferentialFlowPtRP(Int_t aBin, Double_t av, Double_t anError)
virtual void GetOutputHistograms(TList *outputListHistos)
void SetHistProVPtPOI(TProfile *const aHistProVPtPOI)
Bool_t InRPSelection() const
Bool_t FillControlHistograms(AliFlowEventSimple *anEvent, TList *weightsList=NULL, Bool_t usePhiWeights=kFALSE, Bool_t usePtWeights=kFALSE, Bool_t useEtaWeights=kFALSE)
void SetHistProVPtRP(TProfile *const aHistProVPtRP)
static AliFlowLYZConstants * GetMaster()
Double_t Phi() const
Bool_t FillIntegratedFlowRP(Double_t aV, Double_t anError)
int Int_t
Definition: External.C:63
void SetHistQsumforChi(TH1F *const aHistQsumforChi)
void SetFirstr0theta(TProfile *const aFirstr0theta)
Double_t GetEntriesInPtBinRP(Int_t iBin)
static AliFlowCommonConstants * GetMaster()
Double_t GetEntriesInEtaBinRP(Int_t iBin)
Double_t GetWR() const
void SetHistProVetaRP(TProfile *const aHistProVetaRP)
Bool_t FillIntegratedFlow(Double_t aV, Double_t anError)
Bool_t FillDifferentialFlowEtaPOI(Int_t aBin, Double_t av, Double_t anError)
Bool_t FillIntegratedFlowPOI(Double_t aV, Double_t anError)
Bool_t FillDifferentialFlowPtPOI(Int_t aBin, Double_t av, Double_t anError)
Double_t Pt() const
Double_t GetPsi() const
void SetHistProVetaPOI(TProfile *const aHistProVetaPOI)
Bool_t InPOISelection(Int_t poiType=1) const
Double_t Eta() const
bool Bool_t
Definition: External.C:53
Bool_t FillDifferentialFlowEtaRP(Int_t aBin, Double_t av, Double_t anError)
void CalculateRPandW(AliFlowVector aQ)
Int_t NumberOfTracks() const