[u/mrichter/AliRoot.git] / PWGPP / VZERO / AliAnaVZEROEPFlatenning.cxx

#include <stdio.h>
#include <stdlib.h>

#include "TH2F.h"
#include "TFile.h"
#include "TProfile.h"
#include "TList.h"
#include "TChain.h"

#include "AliAnalysisTask.h"
#include "AliAnalysisManager.h"

#include "AliVEvent.h"
#include "AliESDEvent.h"
#include "AliAODEvent.h"
#include "AliInputEventHandler.h"
#include "AliAnaVZEROEPFlatenning.h"
#include "AliCentrality.h"
#include "AliEventplane.h"

// VZERO includes
#include "AliVVZERO.h"

ClassImp(AliAnaVZEROEPFlatenning)

AliAnaVZEROEPFlatenning::AliAnaVZEROEPFlatenning() 
  : AliAnalysisTaskSE("AliAnaVZEROEPFlatenning"), fEvent(0), fOutputList(0),
  fMBTrigName("CPBI"),
  fUsePhysSel(kFALSE),
  fPsiAC(NULL),
  fPsiACOrg(NULL)
{
  // Default constructor
  // Init pointers
  for(Int_t i = 0; i < 8; ++i) fX2[i] = fY2[i] = fX2Y2[i] = fCos8Psi[i] = NULL;
  for(Int_t i = 0; i < 8; ++i) fX2In[i] = fY2In[i] = fX2Y2In[i] = fCos8PsiIn[i] = NULL;
  for(Int_t i = 0; i < 8; ++i) fPsiRingRawCentr[i] = fPsiRingFlatCentr[i] = fPsiRingFlatFinalCentr[i] = NULL;
  for(Int_t i = 0; i < 8; ++i) fC2[i] = fS2[i] = fC4[i] = fS4[i] = NULL;
  // Define input and output slots here
  // Input slot #0 works with a TChain
  DefineInput(0, TChain::Class());
  // Output slot #0 id reserved by the base class for AOD
  // Output slot #1 writes into a TH1 container
  DefineOutput(1, TList::Class());
}

//________________________________________________________________________
AliAnaVZEROEPFlatenning::AliAnaVZEROEPFlatenning(const char *name) 
  : AliAnalysisTaskSE(name), fEvent(0), fOutputList(0),
  fMBTrigName("CPBI"),
  fUsePhysSel(kFALSE),
  fPsiAC(NULL),
  fPsiACOrg(NULL)
{
  // Constructor
  // Init pointers
  for(Int_t i = 0; i < 8; ++i) fX2[i] = fY2[i] = fX2Y2[i] = fCos8Psi[i] = NULL;
  for(Int_t i = 0; i < 8; ++i) fX2In[i] = fY2In[i] = fX2Y2In[i] = fCos8PsiIn[i] = NULL;
  for(Int_t i = 0; i < 8; ++i) fPsiRingRawCentr[i] = fPsiRingFlatCentr[i] = fPsiRingFlatFinalCentr[i] = NULL;
  for(Int_t i = 0; i < 8; ++i) fC2[i] = fS2[i] = fC4[i] = fS4[i] = NULL;
  // Define input and output slots here
  // Input slot #0 works with a TChain
  DefineInput(0, TChain::Class());
  // Output slot #0 id reserved by the base class for AOD
  // Output slot #1 writes into a TH1 container
  DefineOutput(1, TList::Class());
}

//________________________________________________________________________
void AliAnaVZEROEPFlatenning::SetInput(const char *filename)
{
  // Read recentering
  // parameters from an input file

  AliInfo(Form("Reading input histograms from %s",filename));
  TFile *f = TFile::Open(filename);
  TList *list = (TList*)f->Get("coutput");

  for(Int_t i = 0; i < 8; ++i) {
    fX2In[i] = (TProfile*)list->FindObject(Form("fX2_%d",i))->Clone(Form("fX2In_%d",i));
    fX2In[i]->SetDirectory(0);
    fY2In[i] = (TProfile*)list->FindObject(Form("fY2_%d",i))->Clone(Form("fY2In_%d",i));
    fY2In[i]->SetDirectory(0);
    fX2Y2In[i] = (TProfile*)list->FindObject(Form("fX2Y2_%d",i))->Clone(Form("fX2Y2In_%d",i));
    fX2Y2In[i]->SetDirectory(0);
    fCos8PsiIn[i] = (TProfile*)list->FindObject(Form("fCos8Psi_%d",i))->Clone(Form("fCos8PsiIn_%d",i));
    fCos8PsiIn[i]->SetDirectory(0);
  }
  f->Close();
}

//________________________________________________________________________
void AliAnaVZEROEPFlatenning::UserCreateOutputObjects()
{
  // Create histograms
  // Called once

  fOutputList = new TList();
  fOutputList->SetOwner(kTRUE);

  for(Int_t i = 0; i < 8; ++i) {
    fX2[i] = new TProfile(Form("fX2_%d",i),"",21,0,105,"s");
    fOutputList->Add(fX2[i]);
    fY2[i] = new TProfile(Form("fY2_%d",i),"",21,0,105,"s");
    fOutputList->Add(fY2[i]);
    fX2Y2[i] = new TProfile(Form("fX2Y2_%d",i),"",21,0,105,"s");
    fOutputList->Add(fX2Y2[i]);
    fCos8Psi[i] = new TProfile(Form("fCos8Psi_%d",i),"",21,0,105,"s");
    fOutputList->Add(fCos8Psi[i]);
  }

  for(Int_t i = 0; i < 8; ++i) {
    fPsiRingRawCentr[i] = new TH2F(Form("fPsiRingRawCentr_%d",i),"",105,0,105,100,-TMath::Pi()/2,TMath::Pi()/2);
    fOutputList->Add(fPsiRingRawCentr[i]);
    fPsiRingFlatCentr[i] = new TH2F(Form("fPsiRingFlatCentr_%d",i),"",105,0,105,100,-TMath::Pi()/2,TMath::Pi()/2);
    fOutputList->Add(fPsiRingFlatCentr[i]);
    fPsiRingFlatFinalCentr[i] = new TH2F(Form("fPsiRingFlatFinalCentr_%d",i),"",105,0,105,100,-TMath::Pi()/2,TMath::Pi()/2);
    fOutputList->Add(fPsiRingFlatFinalCentr[i]);
  }

  for(Int_t i = 0; i < 8; ++i) {
    fC2[i] = new TProfile(Form("fC2_%d",i),"",21,0,105,"s");
    fOutputList->Add(fC2[i]);
    fS2[i] = new TProfile(Form("fS2_%d",i),"",21,0,105,"s");
    fOutputList->Add(fS2[i]);
    fC4[i] = new TProfile(Form("fC4_%d",i),"",21,0,105,"s");
    fOutputList->Add(fC4[i]);
    fS4[i] = new TProfile(Form("fS4_%d",i),"",21,0,105,"s");
    fOutputList->Add(fS4[i]);
  }

  fPsiAC = new TH2F("fPsiAC","",100,-TMath::Pi()/2,TMath::Pi()/2,100,-TMath::Pi()/2,TMath::Pi()/2);
  fOutputList->Add(fPsiAC);
  fPsiACOrg = new TH2F("fPsiACOrg","",100,-TMath::Pi()/2,TMath::Pi()/2,100,-TMath::Pi()/2,TMath::Pi()/2);
  fOutputList->Add(fPsiACOrg);

  PostData(1, fOutputList);
}

//________________________________________________________________________
void AliAnaVZEROEPFlatenning::Init()
{
  // Nothing here
  // ....
}

//________________________________________________________________________
void AliAnaVZEROEPFlatenning::UserExec(Option_t *) 
{
  // Main loop
  // Called for each event


  fEvent = InputEvent();
  if (!fEvent) {
    printf("ERROR: fEvent not available\n");
    return;
  }

  AliVVZERO* esdV0 = fEvent->GetVZEROData();
  if (!esdV0) {
    Printf("ERROR: esd/aod V0  not available");
    return;
  }

  // Phys sel
  Bool_t goodEvent = kTRUE;
  Bool_t isSelected;
  if (fUsePhysSel)
    isSelected = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & (AliVEvent::kMB | AliVEvent::kSemiCentral));
  else
    isSelected = ((esdV0->GetV0ADecision()==1) && (esdV0->GetV0CDecision()==1));

  if (!isSelected) goodEvent = kFALSE;

  AliCentrality *centrality = fEvent->GetCentrality();
  //  Float_t percentile = centrality->GetCentralityPercentile("V0M");
  Float_t spdPercentile = centrality->GetCentralityPercentile("CL1");

  TString inputDataType = AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()->GetDataType(); // can be "ESD" or "AOD"
  if (inputDataType == "ESD") {
    AliESDEvent* esdEvent = static_cast<AliESDEvent*>(fEvent);
    // Trigger
    TString trigStr(esdEvent->GetFiredTriggerClasses());
    if (!trigStr.Contains("-B-")) return;
    if (!trigStr.Contains(fMBTrigName.Data())) return;

    const AliESDVertex *primaryVtx = esdEvent->GetPrimaryVertexSPD();
    if (!primaryVtx) goodEvent = kFALSE;
    if (!primaryVtx->GetStatus()) goodEvent = kFALSE;
    Double_t tPrimaryVtxPosition[3];
    primaryVtx->GetXYZ(tPrimaryVtxPosition);
    if (TMath::Abs(tPrimaryVtxPosition[2]) > 10.0) goodEvent = kFALSE;
  }
  else if (inputDataType == "AOD") {
    AliAODEvent* aodEvent = static_cast<AliAODEvent*>(fEvent);
    // Trigger
    TString trigStr(aodEvent->GetHeader()->GetFiredTriggerClasses());
    if (!trigStr.Contains("-B-")) return;
    if (!trigStr.Contains(fMBTrigName.Data())) return;

    const AliAODVertex *primaryVtx = aodEvent->GetPrimaryVertexSPD();
    if (!primaryVtx) goodEvent = kFALSE;
    if (primaryVtx->GetNContributors()==0) goodEvent = kFALSE;
    Double_t tPrimaryVtxPosition[3];
    primaryVtx->GetXYZ(tPrimaryVtxPosition);
    if (TMath::Abs(tPrimaryVtxPosition[2]) > 10.0) goodEvent = kFALSE;
  }
  else {
    AliFatal("Trying to get the vertex from neither ESD nor AOD event!");
    return;
  }

  if (goodEvent) {
    Double_t qxTierce[8],qyTierce[8];
    for(Int_t iring = 0; iring < 8; ++iring) {
      qxTierce[iring] = qyTierce[iring] = 0.;
      Double_t c2 = 0;
      Double_t s2 = 0;
      Double_t totMult = 0;
      for(Int_t iCh = iring*8; iCh < (iring+1)*8; ++iCh) {
	Double_t phi = TMath::Pi()/8. + TMath::Pi()/4.*(iCh%8);
	c2 += fEvent->GetVZEROEqMultiplicity(iCh)*TMath::Cos(2.*phi);
	s2 += fEvent->GetVZEROEqMultiplicity(iCh)*TMath::Sin(2.*phi);
	if (fEvent->GetVZEROEqMultiplicity(iCh) > 1e-6) {
	  fC2[iring]->Fill(spdPercentile,TMath::Cos(2.*phi),fEvent->GetVZEROEqMultiplicity(iCh));
	  fS2[iring]->Fill(spdPercentile,TMath::Sin(2.*phi),fEvent->GetVZEROEqMultiplicity(iCh));
	  fC4[iring]->Fill(spdPercentile,TMath::Cos(4.*phi),fEvent->GetVZEROEqMultiplicity(iCh));
	  fS4[iring]->Fill(spdPercentile,TMath::Sin(4.*phi),fEvent->GetVZEROEqMultiplicity(iCh));
	}
	totMult += fEvent->GetVZEROEqMultiplicity(iCh);
      }
      if (totMult < 1e-6) continue;

      fX2[iring]->Fill(spdPercentile,c2);
      fY2[iring]->Fill(spdPercentile,s2);
      fX2Y2[iring]->Fill(spdPercentile,c2*s2);

      Double_t qxOut = 0, qyOut = 0;
      Double_t psiRingRaw = fEvent->GetEventplane()->CalculateVZEROEventPlane(fEvent,iring,iring,2,qxOut,qyOut);
      fPsiRingRawCentr[iring]->Fill(spdPercentile,psiRingRaw);
      Double_t psiRingFlat2 = CalculateVZEROEventPlane(fEvent,iring,spdPercentile,qxTierce[iring],qyTierce[iring]);
      fPsiRingFlatCentr[iring]->Fill(spdPercentile,psiRingFlat2);
      fCos8Psi[iring]->Fill(spdPercentile, 2./4.*TMath::Cos(2.*4.*psiRingFlat2));
      Int_t ibin = fCos8PsiIn[iring]->FindBin(spdPercentile);
      Double_t psiRingFlatFinal = psiRingFlat2 + (fCos8PsiIn[iring]->GetBinContent(ibin)*TMath::Sin(2.*4.*psiRingFlat2))/2.;
      if (psiRingFlatFinal > TMath::Pi()/2) psiRingFlatFinal -= TMath::Pi();
      if (psiRingFlatFinal <-TMath::Pi()/2) psiRingFlatFinal += TMath::Pi();
      fPsiRingFlatFinalCentr[iring]->Fill(spdPercentile,psiRingFlatFinal);
    }
    Double_t qxA = qxTierce[4] + qxTierce[5] + qxTierce[6] + qxTierce[7];
    Double_t qyA = qyTierce[4] + qyTierce[5] + qyTierce[6] + qyTierce[7];
    Double_t qxC = qxTierce[0] + qxTierce[1] + qxTierce[2] + qxTierce[3];
    Double_t qyC = qyTierce[0] + qyTierce[1] + qyTierce[2] + qyTierce[3];
    Double_t psiA = TMath::ATan2(qyA,qxA)/2.;
    Double_t psiC = TMath::ATan2(qyC,qxC)/2.;

    fPsiAC->Fill(psiA,psiC);
    Double_t psiAOrg = fEvent->GetEventplane()->GetEventplane("V0A",fEvent,2);
    Double_t psiCOrg = fEvent->GetEventplane()->GetEventplane("V0C",fEvent,2);
    fPsiACOrg->Fill(psiAOrg,psiCOrg);
  }

  PostData(1, fOutputList);
}      

//________________________________________________________________________
void AliAnaVZEROEPFlatenning::Terminate(Option_t *) 
{
  // Check the output list and store output config file
  // Called once at the end of the query

  fOutputList = dynamic_cast<TList*> (GetOutputData(1));
  if (!fOutputList) {
    printf("ERROR: Output list not available\n");
    return;
  }
}

Double_t AliAnaVZEROEPFlatenning::CalculateVZEROEventPlane(const AliVEvent *  event, Int_t ring, Float_t centrality, Double_t &qxTierce, Double_t &qyTierce) const
{
  // Calculate the VZERO event plane
  // taking into account the recentering/twist and rescaling
  // of the cumulants
  qxTierce = qyTierce = 0.;
  if(!event) {
    AliError("No Event received");
    return -1000.;
  }
  AliVVZERO *vzeroData = event->GetVZEROData();
  if(!vzeroData) {
    AliError("Enable to get VZERO Data");
    return -1000.;
  }

  Int_t ibin = fX2In[ring]->FindBin(centrality);
  Double_t meanX2 = fX2In[ring]->GetBinContent(ibin);
  Double_t meanY2 = fY2In[ring]->GetBinContent(ibin);
  Double_t sigmaX2 = fX2In[ring]->GetBinError(ibin);
  Double_t sigmaY2 = fY2In[ring]->GetBinError(ibin);
  Double_t rho = (fX2Y2In[ring]->GetBinContent(ibin)-meanX2*meanY2)/sigmaX2/sigmaY2;

  Double_t b = rho*sigmaX2*sigmaY2*
    TMath::Sqrt(2.*(sigmaX2*sigmaX2+sigmaY2*sigmaY2-2.*sigmaX2*sigmaY2*TMath::Sqrt(1.-rho*rho))/
		((sigmaX2*sigmaX2-sigmaY2*sigmaY2)*(sigmaX2*sigmaX2-sigmaY2*sigmaY2)+
		 4.*sigmaX2*sigmaX2*sigmaY2*sigmaY2*rho*rho));
  Double_t aPlus = TMath::Sqrt(2.*sigmaX2*sigmaX2-b*b);
  Double_t aMinus= TMath::Sqrt(2.*sigmaY2*sigmaY2-b*b);

  Double_t lambdaPlus = b/aPlus;
  Double_t lambdaMinus = b/aMinus;

  Double_t trans[2][2];
  trans[0][0] = 1./(1.-lambdaPlus*lambdaMinus);
  trans[0][1] = -lambdaMinus/(1.-lambdaPlus*lambdaMinus);
  trans[1][0] = -lambdaPlus/(1.-lambdaPlus*lambdaMinus);
  trans[1][1] = 1./(1.-lambdaPlus*lambdaMinus);

  Double_t qx=0., qy=0.;
  for(Int_t iCh = ring*8; iCh < (ring+1)*8; ++iCh) {
    Double_t phi = TMath::Pi()/8. + TMath::Pi()/4.*(iCh%8);
    Double_t mult = event->GetVZEROEqMultiplicity(iCh);
    qx += mult*TMath::Cos(2.*phi);
    qy += mult*TMath::Sin(2.*phi);
  }
  // Recentering
  Double_t qxPrime = qx - meanX2;
  Double_t qyPrime = qy - meanY2;
  // Twist
  Double_t qxSeconde = trans[0][0]*qxPrime + trans[0][1]*qyPrime;
  Double_t qySeconde = trans[1][0]*qxPrime + trans[1][1]*qyPrime;
  // Rescaling
  qxTierce = qxSeconde/aPlus;
  qyTierce = qySeconde/aMinus;

  return (TMath::ATan2(qyTierce,qxTierce)/2.);
}
Commit	Line	Data
345e2385	1	#include <stdio.h>
	2	#include <stdlib.h>
	3
	4	#include "TH2F.h"
	5	#include "TFile.h"
	6	#include "TProfile.h"
	7	#include "TList.h"
	8	#include "TChain.h"
	9
	10	#include "AliAnalysisTask.h"
	11	#include "AliAnalysisManager.h"
	12
6b4456fd	13	#include "AliVEvent.h"
345e2385	14	#include "AliESDEvent.h"
6b4456fd	15	#include "AliAODEvent.h"
6b4456fd	16	#include "AliInputEventHandler.h"
345e2385	17	#include "AliAnaVZEROEPFlatenning.h"
	18	#include "AliCentrality.h"
	19	#include "AliEventplane.h"
	20
	21	// VZERO includes
6b4456fd	22	#include "AliVVZERO.h"
345e2385	23
	24	ClassImp(AliAnaVZEROEPFlatenning)
	25
	26	AliAnaVZEROEPFlatenning::AliAnaVZEROEPFlatenning()
6b4456fd	27	: AliAnalysisTaskSE("AliAnaVZEROEPFlatenning"), fEvent(0), fOutputList(0),
345e2385	28	fMBTrigName("CPBI"),
fe509094	29	fUsePhysSel(kFALSE),
	30	fPsiAC(NULL),
	31	fPsiACOrg(NULL)
345e2385	32	{
	33	// Default constructor
	34	// Init pointers
	35	for(Int_t i = 0; i < 8; ++i) fX2[i] = fY2[i] = fX2Y2[i] = fCos8Psi[i] = NULL;
	36	for(Int_t i = 0; i < 8; ++i) fX2In[i] = fY2In[i] = fX2Y2In[i] = fCos8PsiIn[i] = NULL;
	37	for(Int_t i = 0; i < 8; ++i) fPsiRingRawCentr[i] = fPsiRingFlatCentr[i] = fPsiRingFlatFinalCentr[i] = NULL;
fe509094	38	for(Int_t i = 0; i < 8; ++i) fC2[i] = fS2[i] = fC4[i] = fS4[i] = NULL;
345e2385	39	// Define input and output slots here
	40	// Input slot #0 works with a TChain
	41	DefineInput(0, TChain::Class());
	42	// Output slot #0 id reserved by the base class for AOD
	43	// Output slot #1 writes into a TH1 container
	44	DefineOutput(1, TList::Class());
	45	}
	46
	47	//________________________________________________________________________
	48	AliAnaVZEROEPFlatenning::AliAnaVZEROEPFlatenning(const char *name)
6b4456fd	49	: AliAnalysisTaskSE(name), fEvent(0), fOutputList(0),
345e2385	50	fMBTrigName("CPBI"),
fe509094	51	fUsePhysSel(kFALSE),
	52	fPsiAC(NULL),
	53	fPsiACOrg(NULL)
345e2385	54	{
	55	// Constructor
	56	// Init pointers
	57	for(Int_t i = 0; i < 8; ++i) fX2[i] = fY2[i] = fX2Y2[i] = fCos8Psi[i] = NULL;
	58	for(Int_t i = 0; i < 8; ++i) fX2In[i] = fY2In[i] = fX2Y2In[i] = fCos8PsiIn[i] = NULL;
	59	for(Int_t i = 0; i < 8; ++i) fPsiRingRawCentr[i] = fPsiRingFlatCentr[i] = fPsiRingFlatFinalCentr[i] = NULL;
fe509094	60	for(Int_t i = 0; i < 8; ++i) fC2[i] = fS2[i] = fC4[i] = fS4[i] = NULL;
345e2385	61	// Define input and output slots here
	62	// Input slot #0 works with a TChain
	63	DefineInput(0, TChain::Class());
	64	// Output slot #0 id reserved by the base class for AOD
	65	// Output slot #1 writes into a TH1 container
	66	DefineOutput(1, TList::Class());
	67	}
	68
	69	//________________________________________________________________________
	70	void AliAnaVZEROEPFlatenning::SetInput(const char *filename)
	71	{
	72	// Read recentering
	73	// parameters from an input file
	74
	75	AliInfo(Form("Reading input histograms from %s",filename));
	76	TFile *f = TFile::Open(filename);
	77	TList list = (TList)f->Get("coutput");
	78
	79	for(Int_t i = 0; i < 8; ++i) {
	80	fX2In[i] = (TProfile*)list->FindObject(Form("fX2_%d",i))->Clone(Form("fX2In_%d",i));
	81	fX2In[i]->SetDirectory(0);
	82	fY2In[i] = (TProfile*)list->FindObject(Form("fY2_%d",i))->Clone(Form("fY2In_%d",i));
	83	fY2In[i]->SetDirectory(0);
	84	fX2Y2In[i] = (TProfile*)list->FindObject(Form("fX2Y2_%d",i))->Clone(Form("fX2Y2In_%d",i));
	85	fX2Y2In[i]->SetDirectory(0);
	86	fCos8PsiIn[i] = (TProfile*)list->FindObject(Form("fCos8Psi_%d",i))->Clone(Form("fCos8PsiIn_%d",i));
	87	fCos8PsiIn[i]->SetDirectory(0);
	88	}
	89	f->Close();
	90	}
	91
	92	//________________________________________________________________________
	93	void AliAnaVZEROEPFlatenning::UserCreateOutputObjects()
	94	{
	95	// Create histograms
	96	// Called once
	97
	98	fOutputList = new TList();
	99	fOutputList->SetOwner(kTRUE);
	100
	101	for(Int_t i = 0; i < 8; ++i) {
	102	fX2[i] = new TProfile(Form("fX2_%d",i),"",21,0,105,"s");
	103	fOutputList->Add(fX2[i]);
	104	fY2[i] = new TProfile(Form("fY2_%d",i),"",21,0,105,"s");
	105	fOutputList->Add(fY2[i]);
	106	fX2Y2[i] = new TProfile(Form("fX2Y2_%d",i),"",21,0,105,"s");
	107	fOutputList->Add(fX2Y2[i]);
	108	fCos8Psi[i] = new TProfile(Form("fCos8Psi_%d",i),"",21,0,105,"s");
	109	fOutputList->Add(fCos8Psi[i]);
	110	}
	111
	112	for(Int_t i = 0; i < 8; ++i) {
	113	fPsiRingRawCentr[i] = new TH2F(Form("fPsiRingRawCentr_%d",i),"",105,0,105,100,-TMath::Pi()/2,TMath::Pi()/2);
	114	fOutputList->Add(fPsiRingRawCentr[i]);
	115	fPsiRingFlatCentr[i] = new TH2F(Form("fPsiRingFlatCentr_%d",i),"",105,0,105,100,-TMath::Pi()/2,TMath::Pi()/2);
	116	fOutputList->Add(fPsiRingFlatCentr[i]);
	117	fPsiRingFlatFinalCentr[i] = new TH2F(Form("fPsiRingFlatFinalCentr_%d",i),"",105,0,105,100,-TMath::Pi()/2,TMath::Pi()/2);
	118	fOutputList->Add(fPsiRingFlatFinalCentr[i]);
	119	}
	120
fe509094	121	for(Int_t i = 0; i < 8; ++i) {
	122	fC2[i] = new TProfile(Form("fC2_%d",i),"",21,0,105,"s");
	123	fOutputList->Add(fC2[i]);
	124	fS2[i] = new TProfile(Form("fS2_%d",i),"",21,0,105,"s");
	125	fOutputList->Add(fS2[i]);
	126	fC4[i] = new TProfile(Form("fC4_%d",i),"",21,0,105,"s");
	127	fOutputList->Add(fC4[i]);
	128	fS4[i] = new TProfile(Form("fS4_%d",i),"",21,0,105,"s");
	129	fOutputList->Add(fS4[i]);
	130	}
	131
	132	fPsiAC = new TH2F("fPsiAC","",100,-TMath::Pi()/2,TMath::Pi()/2,100,-TMath::Pi()/2,TMath::Pi()/2);
	133	fOutputList->Add(fPsiAC);
	134	fPsiACOrg = new TH2F("fPsiACOrg","",100,-TMath::Pi()/2,TMath::Pi()/2,100,-TMath::Pi()/2,TMath::Pi()/2);
	135	fOutputList->Add(fPsiACOrg);
	136
345e2385	137	PostData(1, fOutputList);
	138	}
	139
	140	//________________________________________________________________________
	141	void AliAnaVZEROEPFlatenning::Init()
	142	{
	143	// Nothing here
	144	// ....
	145	}
	146
	147	//________________________________________________________________________
	148	void AliAnaVZEROEPFlatenning::UserExec(Option_t *)
	149	{
	150	// Main loop
	151	// Called for each event
	152
	153
6b4456fd	154	fEvent = InputEvent();
	155	if (!fEvent) {
	156	printf("ERROR: fEvent not available\n");
345e2385	157	return;
	158	}
	159
6b4456fd	160	AliVVZERO* esdV0 = fEvent->GetVZEROData();
345e2385	161	if (!esdV0) {
6b4456fd	162	Printf("ERROR: esd/aod V0 not available");
345e2385	163	return;
	164	}
	165
345e2385	166	// Phys sel
	167	Bool_t goodEvent = kTRUE;
	168	Bool_t isSelected;
	169	if (fUsePhysSel)
6b4456fd	170	isSelected = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & (AliVEvent::kMB \| AliVEvent::kSemiCentral));
345e2385	171	else
	172	isSelected = ((esdV0->GetV0ADecision()==1) && (esdV0->GetV0CDecision()==1));
	173
	174	if (!isSelected) goodEvent = kFALSE;
	175
6b4456fd	176	AliCentrality *centrality = fEvent->GetCentrality();
345e2385	177	// Float_t percentile = centrality->GetCentralityPercentile("V0M");
	178	Float_t spdPercentile = centrality->GetCentralityPercentile("CL1");
	179
6b4456fd	180	TString inputDataType = AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()->GetDataType(); // can be "ESD" or "AOD"
	181	if (inputDataType == "ESD") {
	182	AliESDEvent* esdEvent = static_cast<AliESDEvent*>(fEvent);
	183	// Trigger
	184	TString trigStr(esdEvent->GetFiredTriggerClasses());
	185	if (!trigStr.Contains("-B-")) return;
	186	if (!trigStr.Contains(fMBTrigName.Data())) return;
	187
	188	const AliESDVertex *primaryVtx = esdEvent->GetPrimaryVertexSPD();
	189	if (!primaryVtx) goodEvent = kFALSE;
	190	if (!primaryVtx->GetStatus()) goodEvent = kFALSE;
	191	Double_t tPrimaryVtxPosition[3];
	192	primaryVtx->GetXYZ(tPrimaryVtxPosition);
	193	if (TMath::Abs(tPrimaryVtxPosition[2]) > 10.0) goodEvent = kFALSE;
	194	}
	195	else if (inputDataType == "AOD") {
	196	AliAODEvent* aodEvent = static_cast<AliAODEvent*>(fEvent);
	197	// Trigger
	198	TString trigStr(aodEvent->GetHeader()->GetFiredTriggerClasses());
	199	if (!trigStr.Contains("-B-")) return;
	200	if (!trigStr.Contains(fMBTrigName.Data())) return;
	201
	202	const AliAODVertex *primaryVtx = aodEvent->GetPrimaryVertexSPD();
	203	if (!primaryVtx) goodEvent = kFALSE;
	204	if (primaryVtx->GetNContributors()==0) goodEvent = kFALSE;
	205	Double_t tPrimaryVtxPosition[3];
	206	primaryVtx->GetXYZ(tPrimaryVtxPosition);
	207	if (TMath::Abs(tPrimaryVtxPosition[2]) > 10.0) goodEvent = kFALSE;
	208	}
	209	else {
	210	AliFatal("Trying to get the vertex from neither ESD nor AOD event!");
	211	return;
	212	}
345e2385	213
345e2385	214	if (goodEvent) {
fe509094	215	Double_t qxTierce[8],qyTierce[8];
345e2385	216	for(Int_t iring = 0; iring < 8; ++iring) {
fe509094	217	qxTierce[iring] = qyTierce[iring] = 0.;
345e2385	218	Double_t c2 = 0;
345e2385	219	Double_t s2 = 0;
345e2385	220	Double_t totMult = 0;
	221	for(Int_t iCh = iring8; iCh < (iring+1)8; ++iCh) {
	222	Double_t phi = TMath::Pi()/8. + TMath::Pi()/4.*(iCh%8);
6b4456fd	223	c2 += fEvent->GetVZEROEqMultiplicity(iCh)TMath::Cos(2.phi);
	224	s2 += fEvent->GetVZEROEqMultiplicity(iCh)TMath::Sin(2.phi);
	225	if (fEvent->GetVZEROEqMultiplicity(iCh) > 1e-6) {
	226	fC2[iring]->Fill(spdPercentile,TMath::Cos(2.*phi),fEvent->GetVZEROEqMultiplicity(iCh));
	227	fS2[iring]->Fill(spdPercentile,TMath::Sin(2.*phi),fEvent->GetVZEROEqMultiplicity(iCh));
	228	fC4[iring]->Fill(spdPercentile,TMath::Cos(4.*phi),fEvent->GetVZEROEqMultiplicity(iCh));
	229	fS4[iring]->Fill(spdPercentile,TMath::Sin(4.*phi),fEvent->GetVZEROEqMultiplicity(iCh));
fe509094	230	}
6b4456fd	231	totMult += fEvent->GetVZEROEqMultiplicity(iCh);
345e2385	232	}
	233	if (totMult < 1e-6) continue;
	234
	235	fX2[iring]->Fill(spdPercentile,c2);
	236	fY2[iring]->Fill(spdPercentile,s2);
	237	fX2Y2[iring]->Fill(spdPercentile,c2*s2);
	238
fe509094	239	Double_t qxOut = 0, qyOut = 0;
6b4456fd	240	Double_t psiRingRaw = fEvent->GetEventplane()->CalculateVZEROEventPlane(fEvent,iring,iring,2,qxOut,qyOut);
345e2385	241	fPsiRingRawCentr[iring]->Fill(spdPercentile,psiRingRaw);
6b4456fd	242	Double_t psiRingFlat2 = CalculateVZEROEventPlane(fEvent,iring,spdPercentile,qxTierce[iring],qyTierce[iring]);
345e2385	243	fPsiRingFlatCentr[iring]->Fill(spdPercentile,psiRingFlat2);
	244	fCos8Psi[iring]->Fill(spdPercentile, 2./4.TMath::Cos(2.4.*psiRingFlat2));
	245	Int_t ibin = fCos8PsiIn[iring]->FindBin(spdPercentile);
	246	Double_t psiRingFlatFinal = psiRingFlat2 + (fCos8PsiIn[iring]->GetBinContent(ibin)TMath::Sin(2.4.*psiRingFlat2))/2.;
	247	if (psiRingFlatFinal > TMath::Pi()/2) psiRingFlatFinal -= TMath::Pi();
	248	if (psiRingFlatFinal <-TMath::Pi()/2) psiRingFlatFinal += TMath::Pi();
	249	fPsiRingFlatFinalCentr[iring]->Fill(spdPercentile,psiRingFlatFinal);
	250	}
fe509094	251	Double_t qxA = qxTierce[4] + qxTierce[5] + qxTierce[6] + qxTierce[7];
	252	Double_t qyA = qyTierce[4] + qyTierce[5] + qyTierce[6] + qyTierce[7];
	253	Double_t qxC = qxTierce[0] + qxTierce[1] + qxTierce[2] + qxTierce[3];
	254	Double_t qyC = qyTierce[0] + qyTierce[1] + qyTierce[2] + qyTierce[3];
	255	Double_t psiA = TMath::ATan2(qyA,qxA)/2.;
	256	Double_t psiC = TMath::ATan2(qyC,qxC)/2.;
	257
	258	fPsiAC->Fill(psiA,psiC);
6b4456fd	259	Double_t psiAOrg = fEvent->GetEventplane()->GetEventplane("V0A",fEvent,2);
6b4456fd	260	Double_t psiCOrg = fEvent->GetEventplane()->GetEventplane("V0C",fEvent,2);
fe509094	261	fPsiACOrg->Fill(psiAOrg,psiCOrg);
345e2385	262	}
	263
	264	PostData(1, fOutputList);
	265	}
	266
	267	//________________________________________________________________________
	268	void AliAnaVZEROEPFlatenning::Terminate(Option_t *)
	269	{
	270	// Check the output list and store output config file
	271	// Called once at the end of the query
	272
	273	fOutputList = dynamic_cast<TList*> (GetOutputData(1));
	274	if (!fOutputList) {
	275	printf("ERROR: Output list not available\n");
	276	return;
	277	}
	278	}
	279
fe509094	280	Double_t AliAnaVZEROEPFlatenning::CalculateVZEROEventPlane(const AliVEvent * event, Int_t ring, Float_t centrality, Double_t &qxTierce, Double_t &qyTierce) const
345e2385	281	{
	282	// Calculate the VZERO event plane
	283	// taking into account the recentering/twist and rescaling
	284	// of the cumulants
fe509094	285	qxTierce = qyTierce = 0.;
345e2385	286	if(!event) {
	287	AliError("No Event received");
	288	return -1000.;
	289	}
	290	AliVVZERO *vzeroData = event->GetVZEROData();
	291	if(!vzeroData) {
	292	AliError("Enable to get VZERO Data");
	293	return -1000.;
	294	}
	295
	296	Int_t ibin = fX2In[ring]->FindBin(centrality);
	297	Double_t meanX2 = fX2In[ring]->GetBinContent(ibin);
	298	Double_t meanY2 = fY2In[ring]->GetBinContent(ibin);
	299	Double_t sigmaX2 = fX2In[ring]->GetBinError(ibin);
	300	Double_t sigmaY2 = fY2In[ring]->GetBinError(ibin);
	301	Double_t rho = (fX2Y2In[ring]->GetBinContent(ibin)-meanX2*meanY2)/sigmaX2/sigmaY2;
	302
	303	Double_t b = rhosigmaX2sigmaY2*
	304	TMath::Sqrt(2.(sigmaX2sigmaX2+sigmaY2sigmaY2-2.sigmaX2sigmaY2TMath::Sqrt(1.-rho*rho))/
	305	((sigmaX2sigmaX2-sigmaY2sigmaY2)(sigmaX2sigmaX2-sigmaY2*sigmaY2)+
	306	4.sigmaX2sigmaX2sigmaY2sigmaY2rhorho));
	307	Double_t aPlus = TMath::Sqrt(2.sigmaX2sigmaX2-b*b);
	308	Double_t aMinus= TMath::Sqrt(2.sigmaY2sigmaY2-b*b);
	309
	310	Double_t lambdaPlus = b/aPlus;
	311	Double_t lambdaMinus = b/aMinus;
	312
	313	Double_t trans[2][2];
	314	trans[0][0] = 1./(1.-lambdaPlus*lambdaMinus);
	315	trans[0][1] = -lambdaMinus/(1.-lambdaPlus*lambdaMinus);
	316	trans[1][0] = -lambdaPlus/(1.-lambdaPlus*lambdaMinus);
	317	trans[1][1] = 1./(1.-lambdaPlus*lambdaMinus);
	318
	319	Double_t qx=0., qy=0.;
	320	for(Int_t iCh = ring8; iCh < (ring+1)8; ++iCh) {
	321	Double_t phi = TMath::Pi()/8. + TMath::Pi()/4.*(iCh%8);
	322	Double_t mult = event->GetVZEROEqMultiplicity(iCh);
	323	qx += multTMath::Cos(2.phi);
	324	qy += multTMath::Sin(2.phi);
	325	}
	326	// Recentering
	327	Double_t qxPrime = qx - meanX2;
	328	Double_t qyPrime = qy - meanY2;
	329	// Twist
	330	Double_t qxSeconde = trans[0][0]qxPrime + trans[0][1]qyPrime;
	331	Double_t qySeconde = trans[1][0]qxPrime + trans[1][1]qyPrime;
	332	// Rescaling
fe509094	333	qxTierce = qxSeconde/aPlus;
fe509094	334	qyTierce = qySeconde/aMinus;
345e2385	335
	336	return (TMath::ATan2(qyTierce,qxTierce)/2.);
	337	}