[u/mrichter/AliRoot.git] / PWGLF / FORWARD / analysis2 / AliFMDEventPlaneFinder.cxx

// This task finds the eventplane
// using the FMD
// 
#include "AliFMDEventPlaneFinder.h"
#include <TList.h>
#include <TMath.h>
#include "AliLog.h"
#include <TH2D.h>
#include <TH3D.h>
#include "TROOT.h"
#include <iostream>
#include <iomanip>
#include <TString.h>
#include <TFile.h>
#include "AliOADBContainer.h"
#include "AliAnalysisManager.h"
#include "AliVEvent.h"
#include "AliEventplane.h"
#include "AliCentrality.h"
#include "AliAODForwardEP.h"
#include "AliAODForwardMult.h"
#include "AliAODEvent.h"

ClassImp(AliFMDEventPlaneFinder)
#if 0
; // For Emacs
#endif 

//____________________________________________________________________
AliFMDEventPlaneFinder::AliFMDEventPlaneFinder()
  : TNamed(),
    fList(0),
    fEvent(0),
    fQt(),
    fQa(),
    fQc(),
    fQ1(),
    fQ2(),
    fQeta(),
    fHistFMDEventplane(0),
    fHistFMDEventplaneA(0),
    fHistFMDEventplaneC(0),
    fHistFMDEventplane1(0),
    fHistFMDEventplane2(0),
    fHistPhiDist(0),
    fHistFMDmTPCep(0),
    fHistFMDvsTPCep(0),
    fHistFMDmVZEROep(0),
    fHistFMDvsVZEROep(0),
    fDebug(0),
    fOADBFileName(0),
    fOADBContainer(0),
    fPhiDist(0),
    fRunNumber(0),
    fUsePhiWeights(0)
{
  // 
  // Constructor 
  //
   
}

//____________________________________________________________________
AliFMDEventPlaneFinder::AliFMDEventPlaneFinder(const char* title)
  : TNamed("fmdEventPlaneFinder", title), 
    fList(0),
    fEvent(0),
    fQt(),
    fQa(),
    fQc(),
    fQ1(),
    fQ2(),
    fQeta(),
    fHistFMDEventplane(0),
    fHistFMDEventplaneA(0),
    fHistFMDEventplaneC(0),
    fHistFMDEventplane1(0),
    fHistFMDEventplane2(0),
    fHistPhiDist(0),
    fHistFMDmTPCep(0),
    fHistFMDvsTPCep(0),
    fHistFMDmVZEROep(0),
    fHistFMDvsVZEROep(0),
    fDebug(0),
    fOADBFileName(0),
    fOADBContainer(0),
    fPhiDist(0),
    fRunNumber(0),
    fUsePhiWeights(1)
{
  // 
  // Constructor 
  // 
  // Parameters:
  //    name Name of object
  //
}

//____________________________________________________________________
AliFMDEventPlaneFinder::AliFMDEventPlaneFinder(const 
						 AliFMDEventPlaneFinder& o)
  : TNamed(o),
    fList(o.fList),
    fEvent(o.fEvent),
    fQt(o.fQt),
    fQa(o.fQa),
    fQc(o.fQc),
    fQ1(o.fQ1),
    fQ2(o.fQ2),
    fQeta(o.fQeta),
    fHistFMDEventplane(o.fHistFMDEventplane),
    fHistFMDEventplaneA(o.fHistFMDEventplaneA),
    fHistFMDEventplaneC(o.fHistFMDEventplaneC),
    fHistFMDEventplane1(o.fHistFMDEventplane1),
    fHistFMDEventplane2(o.fHistFMDEventplane2),
    fHistPhiDist(o.fHistPhiDist),
    fHistFMDmTPCep(o.fHistFMDmTPCep),
    fHistFMDvsTPCep(o.fHistFMDvsTPCep),
    fHistFMDmVZEROep(o.fHistFMDmVZEROep),
    fHistFMDvsVZEROep(o.fHistFMDvsVZEROep),
    fDebug(o.fDebug),
    fOADBFileName(o.fOADBFileName),
    fOADBContainer(o.fOADBContainer),
    fPhiDist(o.fPhiDist),
    fRunNumber(o.fRunNumber),
    fUsePhiWeights(o.fUsePhiWeights)
{
  // 
  // Copy constructor 
  // 
  // Parameters:
  //    o Object to copy from 
  //
}

//____________________________________________________________________
AliFMDEventPlaneFinder::~AliFMDEventPlaneFinder()
{
  // 
  // Destructor 
  //
}

//____________________________________________________________________
AliFMDEventPlaneFinder&
AliFMDEventPlaneFinder::operator=(const AliFMDEventPlaneFinder& o)
{
  // 
  // Assignement operator
  // 
  // Parameters:
  //    o Object to assign from 
  // 
  // Return:
  //    Reference to this object
  //
  if (&o == this) return *this; 
  TNamed::operator=(o);

  fList               = o.fList;
  fEvent              = o.fEvent;
  fQt                 = o.fQt;
  fQa                 = o.fQa;
  fQc                 = o.fQc;
  fQ1                 = o.fQ1;
  fQ2                 = o.fQ2;
  fQeta               = o.fQeta;
  fHistFMDEventplane  = o.fHistFMDEventplane;
  fHistFMDEventplaneA = o.fHistFMDEventplaneA;
  fHistFMDEventplaneC = o.fHistFMDEventplaneC;
  fHistFMDEventplane1 = o.fHistFMDEventplane1;
  fHistFMDEventplane2 = o.fHistFMDEventplane2;
  fHistPhiDist        = o.fHistPhiDist;
  fHistFMDmTPCep      = o.fHistFMDmTPCep;
  fHistFMDvsTPCep     = o.fHistFMDvsTPCep;
  fHistFMDmVZEROep    = o.fHistFMDmVZEROep;
  fHistFMDvsVZEROep   = o.fHistFMDvsVZEROep;
  fDebug              = o.fDebug;
  fOADBFileName       = o.fOADBFileName;
  fOADBContainer      = o.fOADBContainer;
  fPhiDist            = o.fPhiDist;
  fRunNumber          = o.fRunNumber;
  fUsePhiWeights      = o.fUsePhiWeights;

  return *this;
}

//____________________________________________________________________
void
AliFMDEventPlaneFinder::Init(const TAxis& etaAxis)
{
  // Intialize this sub-algorithm 
  //
  // Parameters:
  //   etaAxis   fmd eta axis binning
  //
  fHistFMDEventplane = new TH1D("hFMDEventplane", "FMD eventplane", 
                                100, 0., TMath::Pi());
  fHistFMDEventplane->Sumw2();
  fHistFMDEventplane->SetDirectory(0);
  fList->Add(fHistFMDEventplane);

  //
  fHistFMDEventplaneA = new TH1D("hFMDEventplaneA", "FMD eventplaneA", 
                                100, 0., TMath::Pi());
  fHistFMDEventplaneA->Sumw2();
  fHistFMDEventplaneA->SetDirectory(0);
  fList->Add(fHistFMDEventplaneA);

  //
  fHistFMDEventplaneC = new TH1D("hFMDEventplaneC", "FMD eventplaneC", 
                                100, 0., TMath::Pi());
  fHistFMDEventplaneC->Sumw2();
  fHistFMDEventplaneC->SetDirectory(0);
  fList->Add(fHistFMDEventplaneC);

  //
  fHistFMDEventplane1 = new TH1D("hFMDEventplane1", "FMD eventplane1", 
                                100, 0., TMath::Pi());
  fHistFMDEventplane1->Sumw2();
  fHistFMDEventplane1->SetDirectory(0);
  fList->Add(fHistFMDEventplane1);

  //
  fHistFMDEventplane2 = new TH1D("hFMDEventplane2", "FMD eventplane2", 
                                100, 0., TMath::Pi());
  fHistFMDEventplane2->Sumw2();
  fHistFMDEventplane2->SetDirectory(0);
  fList->Add(fHistFMDEventplane2);

  //
  fHistPhiDist = new TH2D("hPhiDist", "Phi distribution in FMD",
	etaAxis.GetNbins(), etaAxis.GetXmin(), etaAxis.GetXmax(), 
	20, 0., TMath::TwoPi());
  fHistPhiDist->Sumw2();
  fHistPhiDist->SetDirectory(0);
  fList->Add(fHistPhiDist);

  //
  fHistFMDmTPCep = new TH1D("hFMDmTPCep", 
                             "Eventplane from FMD - TPC",
                             100, -TMath::Pi()/2., TMath::Pi()/2.);
  fHistFMDmTPCep->Sumw2();
  fHistFMDmTPCep->SetDirectory(0);
  fList->Add(fHistFMDmTPCep);

  //
  fHistFMDvsTPCep = new TH2F("hFMDvsTPCep", 
                             "Eventplane from FMD vs. TPC",
                             100, 0., TMath::Pi(),
                             100, 0., TMath::Pi());
  fHistFMDvsTPCep->Sumw2();
  fHistFMDvsTPCep->SetDirectory(0);
  fList->Add(fHistFMDvsTPCep);

  //
  fHistFMDmVZEROep = new TH1D("hFMDmVZEROep", 
                             "Eventplane from FMD - VZERO",
                             100, -TMath::Pi()/2., TMath::Pi()/2.);
  fHistFMDmVZEROep->Sumw2();
  fHistFMDmVZEROep->SetDirectory(0);
  fList->Add(fHistFMDmVZEROep);

  //
  fHistFMDvsVZEROep = new TH2F("hFMDvsVZEROep", 
                             "Eventplane from FMD vs. VZERO",
                             100, 0., TMath::Pi(),
                             100, 0., TMath::Pi());
  fHistFMDvsVZEROep->Sumw2();
  fHistFMDvsVZEROep->SetDirectory(0);
  fList->Add(fHistFMDvsVZEROep);

  if (!fUsePhiWeights) return;
  
  if (fOADBFileName.Length()==0)
    fOADBFileName = Form("%s/PWGLF/FORWARD/FMDEVENTPLANE/data/fmdEPoadb.root", 
     			AliAnalysisManager::GetOADBPath());

  TFile foadb(fOADBFileName);
  if(!foadb.IsOpen()) {
    AliError(Form("Cannot open OADB file %s", fOADBFileName.Data()));
    return;
  }
  fOADBContainer = (AliOADBContainer*)foadb.Get("FMDphidist");
  if (!fOADBContainer) AliError(Form("No OADB container found in %s", fOADBFileName.Data()));
  foadb.Close();

  return;
}

//_____________________________________________________________________
void
AliFMDEventPlaneFinder::DefineOutput(TList* dir)
{
  // 
  // Output diagnostic histograms to directory 
  // 
  // Parameters:
  //    dir List to write in
  //  
  fList = new TList();
  fList->SetOwner();
  fList->SetName(GetName());
  dir->Add(fList);

  return;
}

//____________________________________________________________________
Bool_t
AliFMDEventPlaneFinder::FindEventplane(AliVEvent* ev,
				       AliAODForwardEP& aodep,
                                       TH2D* h,
                                       AliForwardUtil::Histos* hists)
{
  // 
  // Do the calculations 
  // 
  // Parameters:
  //    hists    Histogram cache
  //    ep       calculated eventplane 
  // 
  // Return:
  //    true on successs 

  fQt.Set(0., 0.);
  fQa.Set(0., 0.);
  fQc.Set(0., 0.);
  fQ1.Set(0., 0.);
  fQ2.Set(0., 0.);

  TH1D epEtaHist = aodep.GetHistogram();

  fEvent = ev;
  if (hists) {
    for (UShort_t d=1; d<=3; d++) { 
      UShort_t nr = (d == 1 ? 1 : 2);
      for (UShort_t q=0; q<nr; q++) { 
        Char_t r = (q == 0 ? 'I' : 'O');
        CalcQVectors(hists->Get(d,r), &epEtaHist);
      } // for q
    } // for d
  }
  else if (h) {
    CalcQVectors(h, &epEtaHist);
  }

  aodep.SetEventplane(CalcEventplane(fQt));
  aodep.SetEventplaneA(CalcEventplane(fQa));
  aodep.SetEventplaneC(CalcEventplane(fQc));
  // aodep.SetEventplane1(CalcEventplane(fQ1));
  // aodep.SetEventplane2(CalcEventplane(fQ2));
  
  FillHists(&aodep);

  return kTRUE;
}

//_____________________________________________________________________
void
AliFMDEventPlaneFinder::CalcQVectors(TH2D* h, TH1D* eHist)
{
  //
  // Calculate the Q vectors
  //
  Double_t phi = 0, eta = 0, weight = 0;
  for (Int_t e = 1; e <= h->GetNbinsX(); e++) {
    Double_t qx = 0, qy = 0;
    eta = h->GetXaxis()->GetBinCenter(e);
    for (Int_t p = 1; p <= h->GetNbinsY(); p++) {
      phi = h->GetYaxis()->GetBinCenter(p);
      weight = h->GetBinContent(e, p);
      if (fUsePhiWeights) weight *= GetPhiWeight(e, p);
      // fix for FMD1 hole
      if (e > 168 && p == 14) {
	h->GetBinContent(e, 4);
        if (fUsePhiWeights) weight *= GetPhiWeight(e, 4);
      }
      fHistPhiDist->Fill(eta, phi, weight);
      // for underflowbin total Nch/eta
      fHistPhiDist->Fill(eta, -1., weight);
      
      // increment Q vectors
      qx += weight*TMath::Cos(2.*phi);
      qy += weight*TMath::Sin(2.*phi);
    }
    TVector2 qVec = TVector2(qx, qy);
    fQt += qVec;
    if (eta < 0) fQa += qVec;
    if (eta > 0) fQc += qVec;
    // TODO: Add random ep increments
    fQeta += qVec;
    if (e % 10 == 0) {
      eHist->Fill(eta, CalcEventplane(fQeta));
      fQeta.Set(0., 0.);
    }
  }

  return;
}

//_____________________________________________________________________
Double_t
AliFMDEventPlaneFinder::CalcEventplane(TVector2 v) const
{
  //
  // Calculate the eventplane
  //
  Double_t ep = -1;
 
  if (v.Mod() == 0.) return ep;

  ep = v.Phi();
  ep /= 2.;

  if (fDebug > 0)
    AliInfo(Form("Eventplane found to be: %f", ep));
 
  return ep;
}

//_____________________________________________________________________
void 
AliFMDEventPlaneFinder::FillHists(AliAODForwardEP* fmdEP)
{
  //
  // Fill diagnostics histograms
  //
  Double_t fmdEPt = fmdEP->GetEventplane();
  Double_t fmdEPa = fmdEP->GetEventplaneA();
  Double_t fmdEPc = fmdEP->GetEventplaneC();
  // Double_t fmdEP1 = fmdEP->GetEventplane1();
  // Double_t fmdEP2 = fmdEP->GetEventplane2();

  // FMD hists
  fHistFMDEventplane->Fill(fmdEPt);
  fHistFMDEventplaneA->Fill(fmdEPa);
  fHistFMDEventplaneC->Fill(fmdEPc);
//  fHistFMDEventplane1->Fill(fmdEP1);
//  fHistFMDEventplane2->Fill(fmdEP2);

  // TPC comparison
  AliEventplane* ep = fEvent->GetEventplane();
  Double_t tpcEP = -1;
  if (ep) ep->GetEventplane("Q");

  Double_t diffTPC = fmdEPt - tpcEP;

  if (diffTPC <  TMath::Pi()/2.) diffTPC = TMath::Pi() - diffTPC;
  if (diffTPC >= TMath::Pi()/2.) diffTPC = TMath::Pi() - diffTPC;

  fHistFMDmTPCep->Fill(diffTPC);

//  if (fmdEPt >= TMath::Pi()/2. && fmdEPt - tpcEP >= TMath::Pi()/2.) tpcEP = TMath::Pi()-tpcEP;
//  if (fmdEPt <  TMath::Pi()/2. && tpcEP - fmdEPt >= TMath::Pi()/2.) tpcEP = TMath::Pi()-tpcEP;

  fHistFMDvsTPCep->Fill(fmdEPt, tpcEP);

  // VZERO comparison
  Double_t vzeroEP = -1;
  if (ep) ep->GetEventplane("V0");

  Double_t diffVZERO = fmdEPt - vzeroEP;

  if (diffVZERO <  TMath::Pi()/2.) diffVZERO = TMath::Pi() - diffVZERO;
  if (diffVZERO >= TMath::Pi()/2.) diffVZERO = TMath::Pi() - diffVZERO;

  fHistFMDmVZEROep->Fill(diffVZERO);

//  if (fmdEPt >= TMath::Pi()/2. && fmdEPt - vzeroEP >= TMath::Pi()/2.) vzeroEP = TMath::Pi()-vzeroEP;
//  if (fmdEPt <  TMath::Pi()/2. && vzeroEP - fmdEPt >= TMath::Pi()/2.) vzeroEP = TMath::Pi()-vzeroEP;

  fHistFMDvsVZEROep->Fill(fmdEPt, vzeroEP);
  return;
}

//_____________________________________________________________________
Double_t 
AliFMDEventPlaneFinder::GetPhiWeight(Int_t etaBin, Int_t phiBin) const
{
  //
  // Get phi weight for flattening
  //
  Double_t phiWeight = 1;
  
  if (fPhiDist) {
    Double_t nParticles = fPhiDist->GetBinContent(etaBin, 0);
    Double_t nPhiBins = fPhiDist->GetNbinsY();
    Double_t phiDistValue = fPhiDist->GetBinContent(etaBin, phiBin);
    
    if (phiDistValue > 0) phiWeight = nParticles/nPhiBins/phiDistValue;
  }

  return phiWeight;
}

//____________________________________________________________________
void
AliFMDEventPlaneFinder::SetRunNumber(Int_t run)
{
  // 
  // Set run number
  //
  if (fRunNumber == run) return;
  
  fRunNumber = run;
  if (fUsePhiWeights) GetPhiDist();

  return;
}

//____________________________________________________________________
void
AliFMDEventPlaneFinder::GetPhiDist()
{
  //
  // Get phi dist from OADB
  //
  if (!fOADBContainer) return;

  fPhiDist = (TH2D*)fOADBContainer->GetObject(fRunNumber, "Default")->Clone(Form("fPhiDist_%d", fRunNumber));
  fList->Add(fPhiDist);

  return;
}

//____________________________________________________________________
void
AliFMDEventPlaneFinder::Print(Option_t* /*option*/) const
{
  // 
  // Print information 
  // 
  // Parameters:
  //    option Not used
  //
  char ind[gROOT->GetDirLevel()+3];
  for (Int_t i = 0; i < gROOT->GetDirLevel(); i++) ind[i] = ' ';
  ind[gROOT->GetDirLevel()] = '\0';
  std::cout << ind << ClassName() << ": " << GetName() << '\n'
	    << std::boolalpha 
	    << ind << "Eventplane finder active!" << '\n'
	    << ind << "Loading OADB object for run number: " << fRunNumber << '\n'
	    << ind << "Using phi weights: " << (fUsePhiWeights ? "true" : "false") << '\n'
	    << std::noboolalpha
	    << std::flush;
  std::cout << std::endl;
}
//____________________________________________________________________
//
// EOF
//
Commit	Line	Data
2b556440	1	// This task finds the eventplane
	2	// using the FMD
	3	//
	4	#include "AliFMDEventPlaneFinder.h"
	5	#include <TList.h>
	6	#include <TMath.h>
	7	#include "AliLog.h"
	8	#include <TH2D.h>
	9	#include <TH3D.h>
	10	#include "TROOT.h"
	11	#include <iostream>
	12	#include <iomanip>
	13	#include <TString.h>
	14	#include <TFile.h>
	15	#include "AliOADBContainer.h"
	16	#include "AliAnalysisManager.h"
	17	#include "AliVEvent.h"
	18	#include "AliEventplane.h"
	19	#include "AliCentrality.h"
	20	#include "AliAODForwardEP.h"
	21	#include "AliAODForwardMult.h"
	22	#include "AliAODEvent.h"
	23
	24	ClassImp(AliFMDEventPlaneFinder)
	25	#if 0
	26	; // For Emacs
	27	#endif
	28
	29	//____________________________________________________________________
	30	AliFMDEventPlaneFinder::AliFMDEventPlaneFinder()
	31	: TNamed(),
	32	fList(0),
	33	fEvent(0),
	34	fQt(),
	35	fQa(),
	36	fQc(),
	37	fQ1(),
	38	fQ2(),
	39	fQeta(),
	40	fHistFMDEventplane(0),
	41	fHistFMDEventplaneA(0),
	42	fHistFMDEventplaneC(0),
	43	fHistFMDEventplane1(0),
	44	fHistFMDEventplane2(0),
	45	fHistPhiDist(0),
	46	fHistFMDmTPCep(0),
	47	fHistFMDvsTPCep(0),
	48	fHistFMDmVZEROep(0),
	49	fHistFMDvsVZEROep(0),
	50	fDebug(0),
	51	fOADBFileName(0),
	52	fOADBContainer(0),
	53	fPhiDist(0),
	54	fRunNumber(0),
	55	fUsePhiWeights(0)
	56	{
	57	//
	58	// Constructor
	59	//
	60
	61	}
	62
	63	//____________________________________________________________________
	64	AliFMDEventPlaneFinder::AliFMDEventPlaneFinder(const char* title)
65	: TNamed("fmdEventPlaneFinder", title),
66	fList(0),
67	fEvent(0),
68	fQt(),
69	fQa(),
70	fQc(),
71	fQ1(),
72	fQ2(),
73	fQeta(),
74	fHistFMDEventplane(0),
75	fHistFMDEventplaneA(0),
76	fHistFMDEventplaneC(0),
77	fHistFMDEventplane1(0),
78	fHistFMDEventplane2(0),
79	fHistPhiDist(0),
80	fHistFMDmTPCep(0),
81	fHistFMDvsTPCep(0),
82	fHistFMDmVZEROep(0),
83	fHistFMDvsVZEROep(0),
84	fDebug(0),
85	fOADBFileName(0),
86	fOADBContainer(0),
87	fPhiDist(0),
88	fRunNumber(0),
89	fUsePhiWeights(1)
90	{
91	//
92	// Constructor
93	//
94	// Parameters:
95	// name Name of object
96	//
97	}
98
99	//____________________________________________________________________
100	AliFMDEventPlaneFinder::AliFMDEventPlaneFinder(const
101	AliFMDEventPlaneFinder& o)
102	: TNamed(o),
103	fList(o.fList),
104	fEvent(o.fEvent),
105	fQt(o.fQt),
106	fQa(o.fQa),
107	fQc(o.fQc),
108	fQ1(o.fQ1),
109	fQ2(o.fQ2),
110	fQeta(o.fQeta),
111	fHistFMDEventplane(o.fHistFMDEventplane),
112	fHistFMDEventplaneA(o.fHistFMDEventplaneA),
113	fHistFMDEventplaneC(o.fHistFMDEventplaneC),
114	fHistFMDEventplane1(o.fHistFMDEventplane1),
115	fHistFMDEventplane2(o.fHistFMDEventplane2),
116	fHistPhiDist(o.fHistPhiDist),
117	fHistFMDmTPCep(o.fHistFMDmTPCep),
118	fHistFMDvsTPCep(o.fHistFMDvsTPCep),
119	fHistFMDmVZEROep(o.fHistFMDmVZEROep),
120	fHistFMDvsVZEROep(o.fHistFMDvsVZEROep),
121	fDebug(o.fDebug),
122	fOADBFileName(o.fOADBFileName),
123	fOADBContainer(o.fOADBContainer),
124	fPhiDist(o.fPhiDist),
125	fRunNumber(o.fRunNumber),
126	fUsePhiWeights(o.fUsePhiWeights)
127	{
128	//
129	// Copy constructor
130	//
131	// Parameters:
132	// o Object to copy from
133	//
134	}
135
136	//____________________________________________________________________
137	AliFMDEventPlaneFinder::~AliFMDEventPlaneFinder()
138	{
139	//
140	// Destructor
141	//
142	}
143
144	//____________________________________________________________________
145	AliFMDEventPlaneFinder&
146	AliFMDEventPlaneFinder::operator=(const AliFMDEventPlaneFinder& o)
147	{
148	//
149	// Assignement operator
150	//
151	// Parameters:
152	// o Object to assign from
153	//
154	// Return:
155	// Reference to this object
156	//
157	if (&o == this) return *this;
158	TNamed::operator=(o);
159
160	fList = o.fList;
161	fEvent = o.fEvent;
162	fQt = o.fQt;
163	fQa = o.fQa;
164	fQc = o.fQc;
165	fQ1 = o.fQ1;
166	fQ2 = o.fQ2;
167	fQeta = o.fQeta;
168	fHistFMDEventplane = o.fHistFMDEventplane;
169	fHistFMDEventplaneA = o.fHistFMDEventplaneA;
170	fHistFMDEventplaneC = o.fHistFMDEventplaneC;
171	fHistFMDEventplane1 = o.fHistFMDEventplane1;
172	fHistFMDEventplane2 = o.fHistFMDEventplane2;
173	fHistPhiDist = o.fHistPhiDist;
174	fHistFMDmTPCep = o.fHistFMDmTPCep;
175	fHistFMDvsTPCep = o.fHistFMDvsTPCep;
176	fHistFMDmVZEROep = o.fHistFMDmVZEROep;
177	fHistFMDvsVZEROep = o.fHistFMDvsVZEROep;
178	fDebug = o.fDebug;
179	fOADBFileName = o.fOADBFileName;
180	fOADBContainer = o.fOADBContainer;
181	fPhiDist = o.fPhiDist;
182	fRunNumber = o.fRunNumber;
183	fUsePhiWeights = o.fUsePhiWeights;
184
185	return *this;
186	}
187
188	//____________________________________________________________________
189	void
190	AliFMDEventPlaneFinder::Init(const TAxis& etaAxis)
191	{
192	// Intialize this sub-algorithm
193	//
194	// Parameters:
195	// etaAxis fmd eta axis binning
196	//
197	fHistFMDEventplane = new TH1D("hFMDEventplane", "FMD eventplane",
198	100, 0., TMath::Pi());
199	fHistFMDEventplane->Sumw2();
200	fHistFMDEventplane->SetDirectory(0);
201	fList->Add(fHistFMDEventplane);
202
203	//
204	fHistFMDEventplaneA = new TH1D("hFMDEventplaneA", "FMD eventplaneA",
205	100, 0., TMath::Pi());
206	fHistFMDEventplaneA->Sumw2();
207	fHistFMDEventplaneA->SetDirectory(0);
208	fList->Add(fHistFMDEventplaneA);
209
210	//
211	fHistFMDEventplaneC = new TH1D("hFMDEventplaneC", "FMD eventplaneC",
212	100, 0., TMath::Pi());
213	fHistFMDEventplaneC->Sumw2();
214	fHistFMDEventplaneC->SetDirectory(0);
215	fList->Add(fHistFMDEventplaneC);
216
217	//
218	fHistFMDEventplane1 = new TH1D("hFMDEventplane1", "FMD eventplane1",
219	100, 0., TMath::Pi());
220	fHistFMDEventplane1->Sumw2();
221	fHistFMDEventplane1->SetDirectory(0);
222	fList->Add(fHistFMDEventplane1);
223
224	//
225	fHistFMDEventplane2 = new TH1D("hFMDEventplane2", "FMD eventplane2",
226	100, 0., TMath::Pi());
227	fHistFMDEventplane2->Sumw2();
228	fHistFMDEventplane2->SetDirectory(0);
229	fList->Add(fHistFMDEventplane2);
230
231	//
232	fHistPhiDist = new TH2D("hPhiDist", "Phi distribution in FMD",
233	etaAxis.GetNbins(), etaAxis.GetXmin(), etaAxis.GetXmax(),
234	20, 0., TMath::TwoPi());
235	fHistPhiDist->Sumw2();
236	fHistPhiDist->SetDirectory(0);
237	fList->Add(fHistPhiDist);
238
239	//
240	fHistFMDmTPCep = new TH1D("hFMDmTPCep",
241	"Eventplane from FMD - TPC",
242	100, -TMath::Pi()/2., TMath::Pi()/2.);
243	fHistFMDmTPCep->Sumw2();
244	fHistFMDmTPCep->SetDirectory(0);
245	fList->Add(fHistFMDmTPCep);
246
247	//
248	fHistFMDvsTPCep = new TH2F("hFMDvsTPCep",
249	"Eventplane from FMD vs. TPC",
250	100, 0., TMath::Pi(),
251	100, 0., TMath::Pi());
252	fHistFMDvsTPCep->Sumw2();
253	fHistFMDvsTPCep->SetDirectory(0);
254	fList->Add(fHistFMDvsTPCep);
255
256	//
257	fHistFMDmVZEROep = new TH1D("hFMDmVZEROep",
258	"Eventplane from FMD - VZERO",
259	100, -TMath::Pi()/2., TMath::Pi()/2.);
260	fHistFMDmVZEROep->Sumw2();
261	fHistFMDmVZEROep->SetDirectory(0);
262	fList->Add(fHistFMDmVZEROep);
263
264	//
265	fHistFMDvsVZEROep = new TH2F("hFMDvsVZEROep",
266	"Eventplane from FMD vs. VZERO",
267	100, 0., TMath::Pi(),
268	100, 0., TMath::Pi());
269	fHistFMDvsVZEROep->Sumw2();
270	fHistFMDvsVZEROep->SetDirectory(0);
271	fList->Add(fHistFMDvsVZEROep);
272
273	if (!fUsePhiWeights) return;
274
275	if (fOADBFileName.Length()==0)
276	fOADBFileName = Form("%s/PWGLF/FORWARD/FMDEVENTPLANE/data/fmdEPoadb.root",
277	AliAnalysisManager::GetOADBPath());
278
279	TFile foadb(fOADBFileName);
280	if(!foadb.IsOpen()) {
281	AliError(Form("Cannot open OADB file %s", fOADBFileName.Data()));
282	return;
283	}
284	fOADBContainer = (AliOADBContainer*)foadb.Get("FMDphidist");
285	if (!fOADBContainer) AliError(Form("No OADB container found in %s", fOADBFileName.Data()));
286	foadb.Close();
287
288	return;
289	}
290
291	//_____________________________________________________________________
292	void
293	AliFMDEventPlaneFinder::DefineOutput(TList* dir)
294	{
295	//
296	// Output diagnostic histograms to directory
297	//
298	// Parameters:
299	// dir List to write in
300	//
301	fList = new TList();
302	fList->SetOwner();
303	fList->SetName(GetName());
304	dir->Add(fList);
305
306	return;
307	}
308
309	//____________________________________________________________________
310	Bool_t
311	AliFMDEventPlaneFinder::FindEventplane(AliVEvent* ev,
312	AliAODForwardEP& aodep,
313	TH2D* h,
314	AliForwardUtil::Histos* hists)
315	{
316	//
317	// Do the calculations
318	//
319	// Parameters:
320	// hists Histogram cache
321	// ep calculated eventplane
322	//
323	// Return:
324	// true on successs
325
326	fQt.Set(0., 0.);
327	fQa.Set(0., 0.);
328	fQc.Set(0., 0.);
329	fQ1.Set(0., 0.);
330	fQ2.Set(0., 0.);
331
332	TH1D epEtaHist = aodep.GetHistogram();
333
334	fEvent = ev;
335	if (hists) {
336	for (UShort_t d=1; d<=3; d++) {
337	UShort_t nr = (d == 1 ? 1 : 2);
338	for (UShort_t q=0; q<nr; q++) {
339	Char_t r = (q == 0 ? 'I' : 'O');
340	CalcQVectors(hists->Get(d,r), &epEtaHist);
341	} // for q
342	} // for d
343	}
344	else if (h) {
345	CalcQVectors(h, &epEtaHist);
346	}
347
348	aodep.SetEventplane(CalcEventplane(fQt));
349	aodep.SetEventplaneA(CalcEventplane(fQa));
350	aodep.SetEventplaneC(CalcEventplane(fQc));
351	// aodep.SetEventplane1(CalcEventplane(fQ1));
352	// aodep.SetEventplane2(CalcEventplane(fQ2));
353
354	FillHists(&aodep);
355
356	return kTRUE;
357	}
358
359	//_____________________________________________________________________
360	void
361	AliFMDEventPlaneFinder::CalcQVectors(TH2D* h, TH1D* eHist)
362	{
363	//
364	// Calculate the Q vectors
365	//
366	Double_t phi = 0, eta = 0, weight = 0;
367	for (Int_t e = 1; e <= h->GetNbinsX(); e++) {
368	Double_t qx = 0, qy = 0;
369	eta = h->GetXaxis()->GetBinCenter(e);
370	for (Int_t p = 1; p <= h->GetNbinsY(); p++) {
371	phi = h->GetYaxis()->GetBinCenter(p);
372	weight = h->GetBinContent(e, p);
373	if (fUsePhiWeights) weight *= GetPhiWeight(e, p);
374	// fix for FMD1 hole
375	if (e > 168 && p == 14) {
376	h->GetBinContent(e, 4);
377	if (fUsePhiWeights) weight *= GetPhiWeight(e, 4);
378	}
379	fHistPhiDist->Fill(eta, phi, weight);
380	// for underflowbin total Nch/eta
381	fHistPhiDist->Fill(eta, -1., weight);
382
383	// increment Q vectors
384	qx += weightTMath::Cos(2.phi);
385	qy += weightTMath::Sin(2.phi);
386	}
387	TVector2 qVec = TVector2(qx, qy);
388	fQt += qVec;
389	if (eta < 0) fQa += qVec;
390	if (eta > 0) fQc += qVec;
391	// TODO: Add random ep increments
392	fQeta += qVec;
393	if (e % 10 == 0) {
394	eHist->Fill(eta, CalcEventplane(fQeta));
395	fQeta.Set(0., 0.);
396	}
397	}
398
399	return;
400	}
401
402	//_____________________________________________________________________
403	Double_t
404	AliFMDEventPlaneFinder::CalcEventplane(TVector2 v) const
405	{
406	//
407	// Calculate the eventplane
408	//
409	Double_t ep = -1;
410
411	if (v.Mod() == 0.) return ep;
412
413	ep = v.Phi();
414	ep /= 2.;
415
416	if (fDebug > 0)
417	AliInfo(Form("Eventplane found to be: %f", ep));
418
419	return ep;
420	}
421
422	//_____________________________________________________________________
423	void
424	AliFMDEventPlaneFinder::FillHists(AliAODForwardEP* fmdEP)
425	{
426	//
427	// Fill diagnostics histograms
428	//
429	Double_t fmdEPt = fmdEP->GetEventplane();
430	Double_t fmdEPa = fmdEP->GetEventplaneA();
431	Double_t fmdEPc = fmdEP->GetEventplaneC();
432	// Double_t fmdEP1 = fmdEP->GetEventplane1();
433	// Double_t fmdEP2 = fmdEP->GetEventplane2();
434
435	// FMD hists
436	fHistFMDEventplane->Fill(fmdEPt);
437	fHistFMDEventplaneA->Fill(fmdEPa);
438	fHistFMDEventplaneC->Fill(fmdEPc);
439	// fHistFMDEventplane1->Fill(fmdEP1);
440	// fHistFMDEventplane2->Fill(fmdEP2);
441
442	// TPC comparison
443	AliEventplane* ep = fEvent->GetEventplane();
444	Double_t tpcEP = -1;
445	if (ep) ep->GetEventplane("Q");
446
447	Double_t diffTPC = fmdEPt - tpcEP;
448
449	if (diffTPC < TMath::Pi()/2.) diffTPC = TMath::Pi() - diffTPC;
450	if (diffTPC >= TMath::Pi()/2.) diffTPC = TMath::Pi() - diffTPC;
451
452	fHistFMDmTPCep->Fill(diffTPC);
453
454	// if (fmdEPt >= TMath::Pi()/2. && fmdEPt - tpcEP >= TMath::Pi()/2.) tpcEP = TMath::Pi()-tpcEP;
455	// if (fmdEPt < TMath::Pi()/2. && tpcEP - fmdEPt >= TMath::Pi()/2.) tpcEP = TMath::Pi()-tpcEP;
456
457	fHistFMDvsTPCep->Fill(fmdEPt, tpcEP);
458
459	// VZERO comparison
460	Double_t vzeroEP = -1;
461	if (ep) ep->GetEventplane("V0");
462
463	Double_t diffVZERO = fmdEPt - vzeroEP;
464
465	if (diffVZERO < TMath::Pi()/2.) diffVZERO = TMath::Pi() - diffVZERO;
466	if (diffVZERO >= TMath::Pi()/2.) diffVZERO = TMath::Pi() - diffVZERO;
467
468	fHistFMDmVZEROep->Fill(diffVZERO);
469
470	// if (fmdEPt >= TMath::Pi()/2. && fmdEPt - vzeroEP >= TMath::Pi()/2.) vzeroEP = TMath::Pi()-vzeroEP;
471	// if (fmdEPt < TMath::Pi()/2. && vzeroEP - fmdEPt >= TMath::Pi()/2.) vzeroEP = TMath::Pi()-vzeroEP;
472
473	fHistFMDvsVZEROep->Fill(fmdEPt, vzeroEP);
474	return;
475	}
476
477	//_____________________________________________________________________
478	Double_t
479	AliFMDEventPlaneFinder::GetPhiWeight(Int_t etaBin, Int_t phiBin) const
480	{
481	//
482	// Get phi weight for flattening
483	//
484	Double_t phiWeight = 1;
485
486	if (fPhiDist) {
487	Double_t nParticles = fPhiDist->GetBinContent(etaBin, 0);
488	Double_t nPhiBins = fPhiDist->GetNbinsY();
489	Double_t phiDistValue = fPhiDist->GetBinContent(etaBin, phiBin);
490
491	if (phiDistValue > 0) phiWeight = nParticles/nPhiBins/phiDistValue;
492	}
493
494	return phiWeight;
495	}
496
497	//____________________________________________________________________
498	void
499	AliFMDEventPlaneFinder::SetRunNumber(Int_t run)
500	{
501	//
502	// Set run number
503	//
504	if (fRunNumber == run) return;
505
506	fRunNumber = run;
507	if (fUsePhiWeights) GetPhiDist();
508
509	return;
510	}
511
512	//____________________________________________________________________
513	void
514	AliFMDEventPlaneFinder::GetPhiDist()
515	{
516	//
517	// Get phi dist from OADB
518	//
519	if (!fOADBContainer) return;
520
521	fPhiDist = (TH2D*)fOADBContainer->GetObject(fRunNumber, "Default")->Clone(Form("fPhiDist_%d", fRunNumber));
522	fList->Add(fPhiDist);
523
524	return;
525	}
526
527	//____________________________________________________________________
528	void
529	AliFMDEventPlaneFinder::Print(Option_t* /option/) const
530	{
531	//
532	// Print information
533	//
534	// Parameters:
535	// option Not used
536	//
537	char ind[gROOT->GetDirLevel()+3];
538	for (Int_t i = 0; i < gROOT->GetDirLevel(); i++) ind[i] = ' ';
539	ind[gROOT->GetDirLevel()] = '\0';
540	std::cout << ind << ClassName() << ": " << GetName() << '\n'
541	<< std::boolalpha
542	<< ind << "Eventplane finder active!" << '\n'
543	<< ind << "Loading OADB object for run number: " << fRunNumber << '\n'
544	<< ind << "Using phi weights: " << (fUsePhiWeights ? "true" : "false") << '\n'
545	<< std::noboolalpha
546	<< std::flush;
547	std::cout << std::endl;
548	}
549	//____________________________________________________________________
550	//
551	// EOF
552	//
553
554
555