[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(),
    fHepFMD(0),
    fHepFMDA(0),
    fHepFMDC(0),
    fHepFMDQC1(0),
    fHepFMDQC2(0),
    fHdiffFMDAC(0),
    fHdiffFMDTPC(0),
    fHdiffFMDVZERO(0),
    fHcorrFMDAC(0),
    fHcorrFMDTPC(0),
    fHcorrFMDVZERO(0),
    fHPhi(0),
    fDebug(0),
    fOADBFileName(0),
    fOADBContainer(0),
    fPhiDist(0),
    fRunNumber(0),
    fUsePhiWeights(0)
{
  // 
  // Constructor 
  //
  DGUARD(fDebug,0,"Default CTOR of AliFMDEventPlaneFinder");
   
}

//____________________________________________________________________
AliFMDEventPlaneFinder::AliFMDEventPlaneFinder(const char* title)
  : TNamed("fmdEventPlaneFinder", title), 
    fList(0),
    fEvent(0),
    fQt(),
    fQa(),
    fQc(),
    fQ1(),
    fQ2(),
    fQeta(),
    fHepFMD(0),
    fHepFMDA(0),
    fHepFMDC(0),
    fHepFMDQC1(0),
    fHepFMDQC2(0),
    fHdiffFMDAC(0),
    fHdiffFMDTPC(0),
    fHdiffFMDVZERO(0),
    fHcorrFMDAC(0),
    fHcorrFMDTPC(0),
    fHcorrFMDVZERO(0),
    fHPhi(0),
    fDebug(0),
    fOADBFileName(0),
    fOADBContainer(0),
    fPhiDist(0),
    fRunNumber(0),
    fUsePhiWeights(1)
{
  // 
  // Constructor 
  // 
  // Parameters:
  //    name Name of object
  //
  DGUARD(fDebug,0,"Named CTOR of AliFMDEventPlaneFinder: %s", title);
}

//____________________________________________________________________
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),
    fHepFMD(o.fHepFMD),
    fHepFMDA(o.fHepFMDA),
    fHepFMDC(o.fHepFMDC),
    fHepFMDQC1(o.fHepFMDQC1),
    fHepFMDQC2(o.fHepFMDQC2),
    fHdiffFMDAC(o.fHdiffFMDAC),
    fHdiffFMDTPC(o.fHdiffFMDTPC),
    fHdiffFMDVZERO(o.fHdiffFMDVZERO),
    fHcorrFMDAC(o.fHcorrFMDAC),
    fHcorrFMDTPC(o.fHcorrFMDTPC),
    fHcorrFMDVZERO(o.fHcorrFMDVZERO),
    fHPhi(o.fHPhi),
    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 
  //
  DGUARD(fDebug,0,"Copy CTOR of AliFMDEventPlaneFinder");
}

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

//____________________________________________________________________
AliFMDEventPlaneFinder&
AliFMDEventPlaneFinder::operator=(const AliFMDEventPlaneFinder& o)
{
  // 
  // Assignement operator
  // 
  // Parameters:
  //    o Object to assign from 
  // 
  // Return:
  //    Reference to this object
  //
  DGUARD(fDebug,3,"Assignment of AliFMDEventPlaneFinder");
  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;
  fHepFMD              = o.fHepFMD;
  fHepFMDA             = o.fHepFMDA;
  fHepFMDC             = o.fHepFMDC;
  fHepFMDQC1           = o.fHepFMDQC1;
  fHepFMDQC2           = o.fHepFMDQC2;
  fHdiffFMDAC          = o.fHdiffFMDAC;
  fHdiffFMDTPC         = o.fHdiffFMDTPC;
  fHdiffFMDVZERO       = o.fHdiffFMDVZERO;
  fHcorrFMDAC          = o.fHcorrFMDAC; 
  fHcorrFMDTPC         = o.fHcorrFMDTPC;
  fHcorrFMDVZERO       = o.fHcorrFMDVZERO;
  fHPhi                = o.fHPhi;
  fDebug               = o.fDebug;
  fOADBFileName        = o.fOADBFileName;
  fOADBContainer       = o.fOADBContainer;
  fPhiDist             = o.fPhiDist;
  fRunNumber           = o.fRunNumber;
  fUsePhiWeights       = o.fUsePhiWeights;

  return *this;
}

//____________________________________________________________________
TH1D*
AliFMDEventPlaneFinder::MakePsiRHist(const char* name, 
				     const char* title,
				     Int_t color)
{
  // Generate a 1D histogram of Psi_R. 
  TH1D* ret = new TH1D(name, Form("#Psi_{R} - %s",title), 100, 0, TMath::Pi());
  ret->SetDirectory(0);
  ret->SetXTitle("#Psi_{R} [radians]");
  ret->SetYTitle("Events");
  ret->SetLineColor(color);
  ret->SetFillColor(color);
  ret->SetFillStyle(3001);
  ret->Sumw2();
  fList->Add(ret);
  return ret;
}
//____________________________________________________________________
TH1D*
AliFMDEventPlaneFinder::MakeDiffHist(const char* name, 
				     const char* first,
				     const char* second,
				     Int_t color)
{
  TH1D* ret = new TH1D(name, Form("#Delta#Psi_{R} - %s minus %s",first,second), 
		       100, -TMath::Pi()/2, +TMath::Pi()/2);
  ret->SetDirectory(0);
  ret->SetXTitle(Form("#Psi_{R,%s}-#Psi_{R,%s} [radians]", first, second));
  ret->SetYTitle("Events");
  ret->SetLineColor(color);
  ret->SetFillColor(color);
  ret->SetFillStyle(3001);
  ret->Sumw2();
  fList->Add(ret);
  return ret;

}
//____________________________________________________________________
TH2F*
AliFMDEventPlaneFinder::MakeCorrHist(const char* name, 
				     const char* first,
				     const char* second)
{
  TH2F* ret = new TH2F(name, Form("#Psi_{R} - %s vs %s", first, second), 
		       100, 0, TMath::Pi(), 100, 0, TMath::Pi());
  ret->SetDirectory(0);
  ret->Sumw2();
  ret->SetXTitle(Form("#Psi_{R,%s} [radians]", first));
  ret->SetYTitle(Form("#Psi_{R,%s} [radians]", second));
  ret->SetZTitle("Events");
  fList->Add(ret);
  return ret;
}
   
//____________________________________________________________________
void
AliFMDEventPlaneFinder::Init(const TAxis& etaAxis)
{
  // Intialize this sub-algorithm 
  //
  // Parameters:
  //   etaAxis   fmd eta axis binning
  //
  DGUARD(fDebug,1,"Initalization of AliFMDEventPlaneFinder");

  fHepFMD    = MakePsiRHist("epFMD", "FMD", kRed+1);
  fHepFMDA   = MakePsiRHist("epFMDA","FMD - A side", kGreen+1);
  fHepFMDC   = MakePsiRHist("epFMDC","FMD - C side", kBlue+1);
  fHepFMDQC1 = MakePsiRHist("epFMDQC1", "FMD QC{1}", kMagenta+1);
  fHepFMDQC2 = MakePsiRHist("epEMDQC2", "FMD QC{2}", kCyan+1);

  fHdiffFMDAC     = MakeDiffHist("diffFMDAC",    "FMD-A", "FMD-C", kRed+1);
  fHdiffFMDTPC    = MakeDiffHist("diffFMDTPC",   "FMD",   "TPC",   kGreen+1);
  fHdiffFMDVZERO  = MakeDiffHist("diffFMDVZERO", "FMD",   "VZERO", kBlue+1);
  fHcorrFMDAC     = MakeCorrHist("corrFMDAC",    "FMD-A", "FMD-C");
  fHcorrFMDTPC    = MakeCorrHist("corrFMDTPC",   "FMD",   "TPC");
  fHcorrFMDVZERO  = MakeCorrHist("corrFMDVZERO", "FMD",   "VZERO");

  //
  fHPhi = new TH2D("hPhi", "Phi distribution in FMD",
		   etaAxis.GetNbins(), etaAxis.GetXmin(), etaAxis.GetXmax(), 
		   20, 0., TMath::TwoPi());
  fHPhi->Sumw2();
  fHPhi->SetXTitle("#eta");
  fHPhi->SetYTitle("#phi [radians]");
  fHPhi->SetZTitle("Events");
  fHPhi->SetDirectory(0);
  fList->Add(fHPhi);

  //

  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
  //  
  DGUARD(fDebug,1,"Define output of AliFMDEventPlaneFinder");
  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 
  DGUARD(fDebug,1,"Find the event plane in AliFMDEventPlaneFinder");

  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
  //
  DGUARD(fDebug,2,"Calculate Q-vectors in AliFMDEventPlaneFinder");
  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) {
	weight = h->GetBinContent(e, 4);
        if (fUsePhiWeights) weight *= GetPhiWeight(e, 4);
      }
      fHPhi->Fill(eta, phi, weight);
      // for underflowbin total Nch/eta
      fHPhi->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(const TVector2& v) const
{
  //
  // Calculate the eventplane
  //
  DGUARD(fDebug,2,"Calculate Event plane in AliFMDEventPlaneFinder");
  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;
}

//_____________________________________________________________________
Double_t 
AliFMDEventPlaneFinder::CalcDifference(Double_t a1, Double_t a2) const
{
  Double_t diff = a1 - a2; 
  if (diff <  TMath::Pi()/2) diff = TMath::Pi() - diff;
  if (diff >= TMath::Pi()/2) diff = TMath::Pi() - diff;
  return diff;
}

//_____________________________________________________________________
void 
AliFMDEventPlaneFinder::FillHists(AliAODForwardEP* fmdEP)
{
  //
  // Fill diagnostics histograms
  //
  DGUARD(fDebug,2,"Fill histograms in AliFMDEventPlaneFinder");
  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
  fHepFMD->Fill(fmdEPt);
  fHepFMDA->Fill(fmdEPa);
  fHepFMDC->Fill(fmdEPc);

  fHdiffFMDAC->Fill(CalcDifference(fmdEPa,fmdEPc));
  fHcorrFMDAC->Fill(fmdEPa, fmdEPc);
  // fHepFMDQC1->Fill(fmdEP1);
  // fHepFMDQC2->Fill(fmdEP2);

  // TPC comparison
  AliEventplane* ep = fEvent->GetEventplane();
  Double_t tpcEP   = (ep ? ep->GetEventplane("Q") : -1);
  fHdiffFMDTPC->Fill(CalcDifference(fmdEPt, tpcEP));
  fHcorrFMDTPC->Fill(fmdEPt, tpcEP);

  // VZERO comparison
  Double_t vzeroEP = ep ? ep->GetEventplane("V0", fEvent, 2) : -1;
  fHdiffFMDVZERO->Fill(CalcDifference(fmdEPt, vzeroEP));
  fHcorrFMDVZERO->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) return phiWeight;

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