AliPhysics  6bc8652 (6bc8652)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AliFlowAnalysisWithMCEventPlane.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 // AliFlowAnalysisWithMCEventPlane:
16 // Description: Maker to analyze Flow from the generated MC reaction plane.
17 // This class is used to get the real value of the flow
18 // to compare the other methods to when analysing simulated events
19 // author: N. van der Kolk (kolk@nikhef.nl)
20 
21 #define AliFlowAnalysisWithMCEventPlane_cxx
22 
23 #include "Riostream.h"
24 #include "TFile.h"
25 #include "TProfile.h"
26 #include "TProfile2D.h"
27 #include "TList.h"
28 #include "TH1F.h"
29 #include "TMath.h"
30 #include "TVector2.h"
31 
32 #include "AliFlowCommonConstants.h"
33 #include "AliFlowEventSimple.h"
34 #include "AliFlowTrackSimple.h"
35 #include "AliFlowCommonHist.h"
38 #include "AliFlowVector.h"
39 
40 class AliFlowVector;
41 
42 using std::endl;
43 using std::cout;
45 
46  //-----------------------------------------------------------------------
47 
49  fQsum(NULL),
50  fQ2sum(0),
51  fEventNumber(0),
52  fDebug(kFALSE),
53  fHistList(NULL),
54  fCommonHists(NULL),
55  fCommonHistsRes(NULL),
56  fHistRP(NULL),
57  fHistProIntFlow(NULL),
58  fHistProIntFlowVsM(NULL),
59  fHistProDiffFlowPtEtaRP(NULL),
60  fHistProDiffFlowPtRP(NULL),
61  fHistProDiffFlowEtaRP(NULL),
62  fHistProDiffFlowPtEtaPOI(NULL),
63  fHistProDiffFlowPtPOI(NULL),
64  fHistProDiffFlowEtaPOI(NULL),
65  fHistSpreadOfFlow(NULL),
66  fHarmonic(2),
67  fMixedHarmonicsList(NULL),
68  fEvaluateMixedHarmonics(kFALSE),
69  fMixedHarmonicsSettings(NULL),
70  fnBinsMult(10000),
71  fMinMult(0.),
72  fMaxMult(10000.),
73  fNinCorrelator(2),
74  fMinCorrelator(2),
75  fXinPairAngle(0.5)
76 {
77 
78  // Constructor.
79  fHistList = new TList();
80 
81  fQsum = new TVector2; // flow vector sum
82 
83  fMixedHarmonicsList = new TList();
84  this->InitalizeArraysForMixedHarmonics();
85 
86 }
87 
88 //-----------------------------------------------------------------------
89 
91 {
92  //destructor
93  delete fHistList;
94  delete fQsum;
95 }
96 
97 //-----------------------------------------------------------------------
98 
100 {
101  //store the final results in output .root file
102  TFile *output = new TFile(outputFileName->Data(),"RECREATE");
103  //output->WriteObject(fHistList, "cobjMCEP","SingleKey");
104  fHistList->SetName("cobjMCEP");
105  fHistList->SetOwner(kTRUE);
106  fHistList->Write(fHistList->GetName(), TObject::kSingleKey);
107  delete output;
108 }
109 
110 //-----------------------------------------------------------------------
111 
113 {
114  //store the final results in output .root file
115  TFile *output = new TFile(outputFileName.Data(),"RECREATE");
116  //output->WriteObject(fHistList, "cobjMCEP","SingleKey");
117  fHistList->SetName("cobjMCEP");
118  fHistList->SetOwner(kTRUE);
119  fHistList->Write(fHistList->GetName(), TObject::kSingleKey);
120  delete output;
121 }
122 
123 //-----------------------------------------------------------------------
124 
125 void AliFlowAnalysisWithMCEventPlane::WriteHistograms(TDirectoryFile *outputFileName)
126 {
127  //store the final results in output .root file
128  fHistList->SetName("cobjMCEP");
129  fHistList->SetOwner(kTRUE);
130  outputFileName->Add(fHistList);
131  outputFileName->Write(outputFileName->GetName(), TObject::kSingleKey);
132 }
133 
134 //-----------------------------------------------------------------------
136 
137  //Define all histograms
138  cout<<"---Analysis with the real MC Event Plane---"<<endl;
139 
140  //save old value and prevent histograms from being added to directory
141  //to avoid name clashes in case multiple analaysis objects are used
142  //in an analysis
143  Bool_t oldHistAddStatus = TH1::AddDirectoryStatus();
144  TH1::AddDirectory(kFALSE);
145 
149 
153 
154  fCommonHists = new AliFlowCommonHist("AliFlowCommonHistMCEP");
155  fHistList->Add(fCommonHists);
156  fCommonHistsRes = new AliFlowCommonHistResults("AliFlowCommonHistResultsMCEP","",fHarmonic);
158 
159  // store harmonic in common control histogram:
160  (fCommonHists->GetHarmonic())->Fill(0.5,fHarmonic);
161 
162  fHistRP = new TH1F("Flow_RP_MCEP","Flow_RP_MCEP",100,0.,3.14);
163  fHistRP->SetXTitle("Reaction Plane Angle");
164  fHistRP->SetYTitle("Counts");
165  fHistList->Add(fHistRP);
166 
167  fHistProIntFlow = new TProfile("FlowPro_V_MCEP","FlowPro_V_MCEP",1,0.,1.);
168  fHistProIntFlow->SetLabelSize(0.06);
169  (fHistProIntFlow->GetXaxis())->SetBinLabel(1,"v_{n}{MCEP}");
170  fHistProIntFlow->SetYTitle("");
172 
173  fHistProIntFlowVsM = new TProfile("FlowPro_VsM_MCEP","FlowPro_VsM_MCEP",10000,0.,10000.); // to be improved - hardwired 10000
174  //fHistProIntFlowVsM->SetLabelSize(0.06);
175  (fHistProIntFlowVsM->GetXaxis())->SetTitle("M");
176  fHistProIntFlowVsM->SetYTitle("");
178 
179  fHistProDiffFlowPtEtaRP = new TProfile2D("FlowPro_VPtEtaRP_MCEP","FlowPro_VPtEtaRP_MCEP",iNbinsPt,dPtMin,dPtMax,iNbinsEta,dEtaMin,dEtaMax);
180  fHistProDiffFlowPtEtaRP->SetXTitle("P_{t}");
181  fHistProDiffFlowPtEtaRP->SetYTitle("#eta");
183 
184  fHistProDiffFlowPtRP = new TProfile("FlowPro_VPtRP_MCEP","FlowPro_VPtRP_MCEP",iNbinsPt,dPtMin,dPtMax);
185  fHistProDiffFlowPtRP->SetXTitle("P_{t}");
186  fHistProDiffFlowPtRP->SetYTitle("");
188 
189  fHistProDiffFlowEtaRP = new TProfile("FlowPro_VetaRP_MCEP","FlowPro_VetaRP_MCEP",iNbinsEta,dEtaMin,dEtaMax);
190  fHistProDiffFlowEtaRP->SetXTitle("#eta");
191  fHistProDiffFlowEtaRP->SetYTitle("");
193 
194  fHistProDiffFlowPtEtaPOI = new TProfile2D("FlowPro_VPtEtaPOI_MCEP","FlowPro_VPtEtaPOI_MCEP",iNbinsPt,dPtMin,dPtMax,iNbinsEta,dEtaMin,dEtaMax);
195  fHistProDiffFlowPtEtaPOI->SetXTitle("P_{t}");
196  fHistProDiffFlowPtEtaPOI->SetYTitle("#eta");
198 
199  fHistProDiffFlowPtPOI = new TProfile("FlowPro_VPtPOI_MCEP","FlowPro_VPtPOI_MCEP",iNbinsPt,dPtMin,dPtMax);
200  fHistProDiffFlowPtPOI->SetXTitle("P_{t}");
201  fHistProDiffFlowPtPOI->SetYTitle("");
203 
204  fHistProDiffFlowEtaPOI = new TProfile("FlowPro_VetaPOI_MCEP","FlowPro_VetaPOI_MCEP",iNbinsEta,dEtaMin,dEtaMax);
205  fHistProDiffFlowEtaPOI->SetXTitle("#eta");
206  fHistProDiffFlowEtaPOI->SetYTitle("");
208 
209  fHistSpreadOfFlow = new TH1D("fHistSpreadOfFlow","fHistSpreadOfFlow",1000,-1,1);
210  fHistSpreadOfFlow->SetXTitle("v_{2}");
211  fHistSpreadOfFlow->SetYTitle("counts");
213 
214  fEventNumber = 0; //set number of events to zero
215 
217 
218  TH1::AddDirectory(oldHistAddStatus);
219 }
220 
221 //-----------------------------------------------------------------------
222 
224 
225  //Calculate v2 from the MC reaction plane
226  if (anEvent) {
227 
228  // get the MC reaction plane angle
229  Double_t aRP = anEvent->GetMCReactionPlaneAngle();
230  //fill control histograms
232 
233  //get the Q vector from the FlowEvent
234  AliFlowVector vQ = anEvent->GetQ(fHarmonic);
235  //cout<<"vQ.Mod() = " << vQ.Mod() << endl;
236  //for chi calculation:
237  *fQsum += vQ;
238  //cout<<"fQsum.Mod() = "<<fQsum.Mod()<<endl;
239  fQ2sum += vQ.Mod2();
240  //cout<<"fQ2sum = "<<fQ2sum<<endl;
241 
242  fHistRP->Fill(aRP);
243 
244  Double_t dPhi = 0.;
245  Double_t dv = 0.;
246  Double_t dPt = 0.;
247  Double_t dEta = 0.;
248  //Double_t dPi = TMath::Pi();
249 
250  // profile to calculate flow e-b-y:
251  TProfile *flowEBE = new TProfile("flowEBE","flowEBE",1,0,1);
252 
253  //calculate flow
254  //loop over the tracks of the event
255  Int_t iNumberOfTracks = anEvent->NumberOfTracks();
256  Int_t iNumberOfRPs = anEvent->GetEventNSelTracksRP();
257  for (Int_t i=0;i<iNumberOfTracks;i++)
258  {
259  AliFlowTrackSimple* pTrack = anEvent->GetTrack(i) ;
260  if (pTrack){
261  if (pTrack->InRPSelection()){
262  dPhi = pTrack->Phi();
263  dv = TMath::Cos(fHarmonic*(dPhi-aRP));
264  dPt = pTrack->Pt();
265  dEta = pTrack->Eta();
266  //reference flow:
267  fHistProIntFlow->Fill(0.,dv);
268  //reference flow versus multiplicity:
269  fHistProIntFlowVsM->Fill(iNumberOfRPs+0.5,dv);
270  //reference flow e-b-e:
271  flowEBE->Fill(0.,dv);
272  //differential flow (Pt, Eta, RP):
273  fHistProDiffFlowPtEtaRP->Fill(dPt,dEta,dv,1.);
274  //differential flow (Pt, RP):
275  fHistProDiffFlowPtRP->Fill(dPt,dv,1.);
276  //differential flow (Eta, RP):
277  fHistProDiffFlowEtaRP->Fill(dEta,dv,1.);
278  }
279  if (pTrack->InPOISelection()) {
280  dPhi = pTrack->Phi();
281  //if (dPhi<0.) dPhi+=2*TMath::Pi();
282  //calculate flow v2:
283  dv = TMath::Cos(fHarmonic*(dPhi-aRP));
284  dPt = pTrack->Pt();
285  dEta = pTrack->Eta();
286  //differential flow (Pt, Eta, POI):
287  fHistProDiffFlowPtEtaPOI->Fill(dPt,dEta,dv,1.);
288  //differential flow (Pt, POI):
289  fHistProDiffFlowPtPOI->Fill(dPt,dv,1.);
290  //differential flow (Eta, POI):
291  fHistProDiffFlowEtaPOI->Fill(dEta,dv,1.);
292  }
293  }//track selected
294  }//loop over tracks
295 
296  fEventNumber++;
297  // cout<<"@@@@@ "<<fEventNumber<<" events processed"<<endl;
298 
299  // store flow value for this event:
300  fHistSpreadOfFlow->Fill(flowEBE->GetBinContent(1),flowEBE->GetBinEntries(1));
301  delete flowEBE;
302 
304  }
305 }
306  //--------------------------------------------------------------------
307 
309  // get the pointers to all output histograms before calling Finish()
310  if (outputListHistos) {
311  //Get the common histograms from the output list
312  AliFlowCommonHist *pCommonHists = dynamic_cast<AliFlowCommonHist*>
313  (outputListHistos->FindObject("AliFlowCommonHistMCEP"));
314  AliFlowCommonHistResults *pCommonHistResults =
315  dynamic_cast<AliFlowCommonHistResults*>
316  (outputListHistos->FindObject("AliFlowCommonHistResultsMCEP"));
317 
318  TProfile *pHistProIntFlow = dynamic_cast<TProfile*>
319  (outputListHistos->FindObject("FlowPro_V_MCEP"));
320 
321  TProfile2D *pHistProDiffFlowPtEtaRP = dynamic_cast<TProfile2D*>
322  (outputListHistos->FindObject("FlowPro_VPtEtaRP_MCEP"));
323 
324  TProfile *pHistProDiffFlowPtRP = dynamic_cast<TProfile*>
325  (outputListHistos->FindObject("FlowPro_VPtRP_MCEP"));
326 
327  TProfile *pHistProDiffFlowEtaRP = dynamic_cast<TProfile*>
328  (outputListHistos->FindObject("FlowPro_VetaRP_MCEP"));
329 
330  TProfile2D *pHistProDiffFlowPtEtaPOI = dynamic_cast<TProfile2D*>
331  (outputListHistos->FindObject("FlowPro_VPtEtaPOI_MCEP"));
332 
333  TProfile *pHistProDiffFlowPtPOI = dynamic_cast<TProfile*>
334  (outputListHistos->FindObject("FlowPro_VPtPOI_MCEP"));
335 
336  TProfile *pHistProDiffFlowEtaPOI = dynamic_cast<TProfile*>
337  (outputListHistos->FindObject("FlowPro_VetaPOI_MCEP"));
338 
339  if (pCommonHists && pCommonHistResults && pHistProIntFlow &&
340  pHistProDiffFlowPtRP && pHistProDiffFlowEtaRP &&
341  pHistProDiffFlowPtPOI && pHistProDiffFlowEtaPOI) {
342  this->SetCommonHists(pCommonHists);
343  this->SetCommonHistsRes(pCommonHistResults);
344  this->SetHistProIntFlow(pHistProIntFlow);
345  this->SetHistProDiffFlowPtEtaRP(pHistProDiffFlowPtEtaRP);
346  this->SetHistProDiffFlowPtRP(pHistProDiffFlowPtRP);
347  this->SetHistProDiffFlowEtaRP(pHistProDiffFlowEtaRP);
348  this->SetHistProDiffFlowPtEtaPOI(pHistProDiffFlowPtEtaPOI);
349  this->SetHistProDiffFlowPtPOI(pHistProDiffFlowPtPOI);
350  this->SetHistProDiffFlowEtaPOI(pHistProDiffFlowEtaPOI);
351  } else {
352  cout<<"WARNING: Histograms needed to run Finish() are not accessible!"<<endl; }
353 
354  //fListHistos->Print();
355 
356  TList *pMixedHarmonicsList = dynamic_cast<TList*>
357  (outputListHistos->FindObject("Mixed Harmonics"));
358  if(pMixedHarmonicsList) {this->GetOutputHistoramsForMixedHarmonics(pMixedHarmonicsList);}
359 
360  } else { cout << "histogram list pointer is empty" << endl;}
361 
362 }
363 
364 //--------------------------------------------------------------------
365 
367 {
368  // Get pointers to all objects relevant for mixed harmonics study.
369 
370  if(mixedHarmonicsList)
371  {
372  // ...
373  } else
374  {
375  cout<<endl;
376  cout<<"WARNING (MCEP): mixedHarmonicsList in NULL in MCEP::GetOutputHistoramsForMixedHarmonics() !!!! "<<endl;
377  cout<<endl;
378  }
379 
380 } // end of void AliFlowAnalysisWithMCEventPlane::GetOutputHistoramsForMixedHarmonics(TList *mixedHarmonicsList)
381 
382 //--------------------------------------------------------------------
383 
385 
386  //*************make histograms etc.
387  if (fDebug) cout<<"AliFlowAnalysisWithMCEventPlane::Terminate()"<<endl;
388 
391 
392  // access harmonic:
394  {
395  fHarmonic = (Int_t)(fCommonHists->GetHarmonic())->GetBinContent(1); // to be improved (moved somewhere else?)
396  }
397 
398  //reference flow :
399  Double_t dV = fHistProIntFlow->GetBinContent(1);
400  Double_t dErrV = fHistProIntFlow->GetBinError(1); // to be improved (treatment of errors for non-Gaussian distribution needed!)
401  //fill reference flow:
403  cout<<"dV"<<fHarmonic<<"{MC} is "<<dV<<" +- "<<dErrV<<endl;
404 
405  //RP:
406  TH1F* fHistPtRP = NULL;
408  {
409  fHistPtRP = fCommonHists->GetHistPtRP();
410  }
411  Double_t dYieldPtRP = 0.;
412  Double_t dVRP = 0.;
413  Double_t dErrVRP = 0.;
414  Double_t dSumRP = 0.;
415  //differential flow (RP, Pt):
416  Double_t dvPtRP = 0.;
417  Double_t dErrvPtRP = 0.;
418  for(Int_t b=1;b<=iNbinsPt;b++)
419  {
420  dvPtRP = fHistProDiffFlowPtRP->GetBinContent(b);
421  dErrvPtRP = fHistProDiffFlowPtRP->GetBinError(b);//to be improved (treatment of errors for non-Gaussian distribution needed!)
422  fCommonHistsRes->FillDifferentialFlowPtRP(b, dvPtRP, dErrvPtRP);
423  if(fHistPtRP){
424  //integrated flow (RP)
425  dYieldPtRP = fHistPtRP->GetBinContent(b);
426  dVRP += dvPtRP*dYieldPtRP;
427  dSumRP += dYieldPtRP;
428  //error on integrated flow
429  dErrVRP += dYieldPtRP*dYieldPtRP*dErrvPtRP*dErrvPtRP;
430  }
431  }
432 
433  if (!(TMath::AreEqualAbs(dSumRP, 0.0, 1e-10)) ) {
434  dVRP /= dSumRP; //because pt distribution should be normalised
435  dErrVRP /= (dSumRP*dSumRP);
436  dErrVRP = TMath::Sqrt(dErrVRP);
437  }
438  // fill integrated flow (RP):
439  fCommonHistsRes->FillIntegratedFlowRP(dVRP,dErrVRP);
440  cout<<"dV"<<fHarmonic<<"{MC} (RP) is "<<dVRP<<" +- "<<dErrVRP<<endl;
441 
442  //differential flow (RP, Eta):
443  Double_t dvEtaRP = 0.;
444  Double_t dErrvEtaRP = 0.;
445  for(Int_t b=1;b<=iNbinsEta;b++)
446  {
447  dvEtaRP = fHistProDiffFlowEtaRP->GetBinContent(b);
448  dErrvEtaRP = fHistProDiffFlowEtaRP->GetBinError(b);//to be improved (treatment of errors for non-Gaussian distribution needed!)
449  fCommonHistsRes->FillDifferentialFlowEtaRP(b, dvEtaRP, dErrvEtaRP);
450  }
451 
452  //POI:
453  TH1F* fHistPtPOI = NULL;
455  {
456  fHistPtPOI = fCommonHists->GetHistPtPOI();
457  }
458  Double_t dYieldPtPOI = 0.;
459  Double_t dVPOI = 0.;
460  Double_t dErrVPOI = 0.;
461  Double_t dSumPOI = 0.;
462  Double_t dvproPtPOI = 0.;
463  Double_t dErrdifcombPtPOI = 0.;
464  Double_t dvproEtaPOI = 0.;
465  Double_t dErrdifcombEtaPOI = 0.;
466  //Pt:
468  for(Int_t b=1;b<=iNbinsPt;b++){
469  dvproPtPOI = fHistProDiffFlowPtPOI->GetBinContent(b);
470  dErrdifcombPtPOI = fHistProDiffFlowPtPOI->GetBinError(b);//to be improved (treatment of errors for non-Gaussian distribution needed!)
471  //fill TH1D
472  fCommonHistsRes->FillDifferentialFlowPtPOI(b, dvproPtPOI, dErrdifcombPtPOI);
473  if (fHistPtPOI){
474  //integrated flow (POI)
475  dYieldPtPOI = fHistPtPOI->GetBinContent(b);
476  dVPOI += dvproPtPOI*dYieldPtPOI;
477  dSumPOI += dYieldPtPOI;
478  //error on integrated flow
479  dErrVPOI += dYieldPtPOI*dYieldPtPOI*dErrdifcombPtPOI*dErrdifcombPtPOI;
480  }
481  }//end of for(Int_t b=0;b<iNbinsPt;b++)
482  } else { cout<<"fHistProFlow is NULL"<<endl; }
483 
484  if (!(TMath::AreEqualAbs(dSumPOI, 0.0, 1e-10)) ) {
485  dVPOI /= dSumPOI; //because pt distribution should be normalised
486  dErrVPOI /= (dSumPOI*dSumPOI);
487  dErrVPOI = TMath::Sqrt(dErrVPOI);
488  }
489  cout<<"dV"<<fHarmonic<<"{MC} (POI) is "<<dVPOI<<" +- "<<dErrVPOI<<endl;
490 
491  fCommonHistsRes->FillIntegratedFlowPOI(dVPOI,dErrVPOI);
492 
493  //Eta:
495  {
496  for(Int_t b=1;b<=iNbinsEta;b++)
497  {
498  dvproEtaPOI = fHistProDiffFlowEtaPOI->GetBinContent(b);
499  dErrdifcombEtaPOI = fHistProDiffFlowEtaPOI->GetBinError(b);//to be improved (treatment of errors for non-Gaussian distribution needed!)
500  //fill common hist results:
501  fCommonHistsRes->FillDifferentialFlowEtaPOI(b, dvproEtaPOI, dErrdifcombEtaPOI);
502  }
503  }
504 
505  cout<<endl;
506  //cout<<".....finished"<<endl;
507 }
508 
509 //-----------------------------------------------------------------------
510 
512 {
513  // Iinitialize all arrays for mixed harmonics.
514 
515  for(Int_t cs=0;cs<2;cs++) // cos/sin
516  {
517  fPairCorrelator[cs] = NULL;
518  fPairCorrelatorVsM[cs] = NULL;
519  for(Int_t sd=0;sd<2;sd++) // pt sum/difference
520  {
521  fPairCorrelatorVsPtSumDiff[cs][sd] = NULL;
522  }
523  }
524 
525 } // end of void InitalizeArraysForMixedHarmonics()
526 
527 //-----------------------------------------------------------------------
528 
530 {
531  // Book all objects needed for mixed harmonics.
532 
533  // List holding all objects relevant for mixed harmonics:
534  fMixedHarmonicsList->SetName("Mixed Harmonics");
535  fMixedHarmonicsList->SetOwner(kTRUE);
537 
538  // Profile holding settings relevant for mixed harmonics:
539  TString mixedHarmonicsSettingsName = "fMixedHarmonicsSettings";
540  fMixedHarmonicsSettings = new TProfile(mixedHarmonicsSettingsName.Data(),"Settings for Mixed Harmonics",4,0,4);
541  //fMixedHarmonicsSettings->GetXaxis()->SetLabelSize(0.025);
542  fMixedHarmonicsSettings->GetXaxis()->SetBinLabel(1,"fEvaluateMixedHarmonics");
544  fMixedHarmonicsSettings->GetXaxis()->SetBinLabel(2,"fNinCorrelator");
546  fMixedHarmonicsSettings->GetXaxis()->SetBinLabel(3,"fMinCorrelator");
548  fMixedHarmonicsSettings->GetXaxis()->SetBinLabel(4,"x in #phi_{pair}"); // phi_{pair} = x*phi1+(1-x)*phi2
551 
552  // Profiles used to calculate <cos[m*phi_{pair}-n*RP]> and <sin[m*phi_{pair}-n*RP]>, where phi_{pair} = x*phi1+(1-x)*phi2:
553  TString cosSinFlag[2] = {"Cos","Sin"};
554  TString cosSinTitleFlag[2] = {"cos[m#phi_{pair}-n#psi_{RP}]","sin[m#phi_{pair}-n#psi_{RP}]"};
555  if(fNinCorrelator == 2 && fMinCorrelator == 2 && TMath::Abs(fXinPairAngle-0.5)<1.e-44) // default values
556  {
557  cosSinTitleFlag[0] = "cos[#phi_{1}+#phi_{2}-2#psi_{RP})]";
558  cosSinTitleFlag[1] = "sin[#phi_{1}+#phi_{2}-2#psi_{RP})]";
559  }
560  TString psdFlag[2] = {"Sum","Diff"};
561  TString psdTitleFlag[2] = {"(p_{t,1}+p_{t,2})/2","#left|p_{t,1}-p_{t,2}#right|"};
562  TString pairCorrelatorName = "fPairCorrelator";
563  TString pairCorrelatorVsMName = "fPairCorrelatorVsM";
564  TString pairCorrelatorVsPtSumDiffName = "fPairCorrelatorVsPt";
568  for(Int_t cs=0;cs<2;cs++)
569  {
570  fPairCorrelator[cs] = new TProfile(Form("%s, %s",pairCorrelatorName.Data(),cosSinFlag[cs].Data()),cosSinTitleFlag[cs].Data(),1,0.,1.);
571  fPairCorrelator[cs]->GetXaxis()->SetBinLabel(1,cosSinTitleFlag[cs].Data());
573 
574  fPairCorrelatorVsM[cs] = new TProfile(Form("%s, %s",pairCorrelatorVsMName.Data(),cosSinFlag[cs].Data()),cosSinTitleFlag[cs].Data(),fnBinsMult,fMinMult,fMaxMult);
575  fPairCorrelatorVsM[cs]->GetXaxis()->SetTitle("# of RPs");
577 
578  for(Int_t sd=0;sd<2;sd++)
579  {
580  fPairCorrelatorVsPtSumDiff[cs][sd] = new TProfile(Form("%s%s, %s",pairCorrelatorVsPtSumDiffName.Data(),psdFlag[sd].Data(),cosSinFlag[cs].Data()),cosSinTitleFlag[cs].Data(),iNbinsPt,dPtMin,dPtMax);
581  fPairCorrelatorVsPtSumDiff[cs][sd]->GetXaxis()->SetTitle(psdTitleFlag[sd].Data());
583  } // end of for(Int_t sd=0;sd<2;sd++)
584  } // end of for(Int_t cs=0;cs<2;cs++)
585 
586 } // end of void AliFlowAnalysisWithMCEventPlane::BookObjectsForMixedHarmonics()
587 
588 //-----------------------------------------------------------------------
589 
591 {
592  // Evaluate correlators relevant for the mixed harmonics.
593 
594  // Get the MC reaction plane angle:
595  Double_t dReactionPlane = anEvent->GetMCReactionPlaneAngle();
596  // Get the number of tracks:
597  Int_t iNumberOfTracks = anEvent->NumberOfTracks();
598  Int_t nRP = anEvent->GetEventNSelTracksRP(); // number of Reference Particles
599  AliFlowTrackSimple *pTrack = NULL;
600  Double_t dPhi1 = 0.;
601  Double_t dPhi2 = 0.;
602  Double_t dPt1 = 0.;
603  Double_t dPt2 = 0.;
604  Double_t n = fNinCorrelator; // shortcut
605  Double_t m = fMinCorrelator; // shortcut
606  Double_t x = fXinPairAngle; // shortcut
607  for(Int_t i=0;i<iNumberOfTracks;i++)
608  {
609  pTrack = anEvent->GetTrack(i);
610  if(pTrack && pTrack->InRPSelection())
611  {
612  dPhi1 = pTrack->Phi();
613  dPt1 = pTrack->Pt();
614  }
615  for(Int_t j=0;j<iNumberOfTracks;j++)
616  {
617  if(j==i) continue;
618  pTrack = anEvent->GetTrack(j);
619  if(pTrack && pTrack->InPOISelection())
620  {
621  dPhi2 = pTrack->Phi();
622  dPt2 = pTrack->Pt();
623  }
624  Double_t dPhiPair = x*dPhi1+(1.-x)*dPhi2;
625  Double_t dPtSum = 0.5*(dPt1+dPt2);
626  Double_t dPtDiff = TMath::Abs(dPt1-dPt2);
627  fPairCorrelator[0]->Fill(0.5,TMath::Cos(m*dPhiPair-n*dReactionPlane),1.);
628  fPairCorrelator[1]->Fill(0.5,TMath::Sin(m*dPhiPair-n*dReactionPlane),1.);
629  fPairCorrelatorVsM[0]->Fill(nRP+0.5,TMath::Cos(m*dPhiPair-n*dReactionPlane),1.);
630  fPairCorrelatorVsM[1]->Fill(nRP+0.5,TMath::Sin(m*dPhiPair-n*dReactionPlane),1.);
631  fPairCorrelatorVsPtSumDiff[0][0]->Fill(dPtSum,TMath::Cos(m*dPhiPair-n*dReactionPlane),1.);
632  fPairCorrelatorVsPtSumDiff[1][0]->Fill(dPtSum,TMath::Sin(m*dPhiPair-n*dReactionPlane),1.);
633  fPairCorrelatorVsPtSumDiff[0][1]->Fill(dPtDiff,TMath::Cos(m*dPhiPair-n*dReactionPlane),1.);
634  fPairCorrelatorVsPtSumDiff[1][1]->Fill(dPtDiff,TMath::Sin(m*dPhiPair-n*dReactionPlane),1.);
635  } // end of for(Int_t j=i+1;j<iNumberOfTracks;j++)
636  } // end of for(Int_t i=0;i<iNumberOfTracks;i++)
637 
638 } // end of void AliFlowAnalysisWithMCEventPlane::EvaluateMixedHarmonics()
639 
640 
641 
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)
void SetHistProDiffFlowPtPOI(TProfile *const aHistProDiffFlowPtPOI)
void SetHistProDiffFlowPtRP(TProfile *const aHistProDiffFlowPtRP)
AliFlowTrackSimple * GetTrack(Int_t i)
Bool_t FillDifferentialFlowPtRP(Int_t aBin, Double_t av, Double_t anError)
Int_t GetEventNSelTracksRP() const
void SetHistProDiffFlowEtaPOI(TProfile *const aHistProDiffFlowEtaPOI)
void GetOutputHistograms(TList *outputListHistos)
Bool_t InRPSelection() const
Bool_t FillControlHistograms(AliFlowEventSimple *anEvent, TList *weightsList=NULL, Bool_t usePhiWeights=kFALSE, Bool_t usePtWeights=kFALSE, Bool_t useEtaWeights=kFALSE)
ClassImp(AliFlowAnalysisWithMCEventPlane) AliFlowAnalysisWithMCEventPlane
Double_t Phi() const
Bool_t FillIntegratedFlowRP(Double_t aV, Double_t anError)
int Int_t
Definition: External.C:63
Definition: External.C:212
static AliFlowCommonConstants * GetMaster()
virtual void EvaluateMixedHarmonics(AliFlowEventSimple *anEvent)
void SetHistProDiffFlowEtaRP(TProfile *const aHistProDiffFlowEtaRP)
void SetHistProDiffFlowPtEtaPOI(TProfile2D *const aHistProDiffFlowPtEtaPOI)
TList * fHistList
flag for lyz analysis: more print statements
Bool_t Data(TH1F *h, Double_t *rangefit, Bool_t writefit, Double_t &sgn, Double_t &errsgn, Double_t &bkg, Double_t &errbkg, Double_t &sgnf, Double_t &errsgnf, Double_t &sigmafit, Int_t &status)
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)
void SetHistProDiffFlowPtEtaRP(TProfile2D *const aHistProDiffFlowPtEtaRP)
Bool_t FillDifferentialFlowPtPOI(Int_t aBin, Double_t av, Double_t anError)
Double_t Pt() const
void SetCommonHistsRes(AliFlowCommonHistResults *const aCommonHistResult)
Double_t GetMCReactionPlaneAngle() const
void SetCommonHists(AliFlowCommonHist *const aCommonHist)
Bool_t InPOISelection(Int_t poiType=1) const
Double_t Eta() const
virtual void GetOutputHistoramsForMixedHarmonics(TList *mixedHarmonicsList)
bool Bool_t
Definition: External.C:53
TProfile * GetHarmonic()
Bool_t FillDifferentialFlowEtaRP(Int_t aBin, Double_t av, Double_t anError)
Int_t NumberOfTracks() const
void SetHistProIntFlow(TProfile *const aHistProIntFlow)