[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),
  fPsiARawCentr(NULL),
  fPsiAFlatCentr(NULL),
  fPsiCRawCentr(NULL),
  fPsiCFlatCentr(NULL),
  fPsiACRawCentr(NULL),
  fPsiACFlatCentr(NULL)
{
  // Default constructor
  // Init pointers
  for(Int_t i = 0; i < 11; ++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;
  for(Int_t i = 0; i < 11; ++i) fX2Corr[i] = fY2Corr[i] = fX2Y2Corr[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),
  fPsiARawCentr(NULL),
  fPsiAFlatCentr(NULL),
  fPsiCRawCentr(NULL),
  fPsiCFlatCentr(NULL),
  fPsiACRawCentr(NULL),
  fPsiACFlatCentr(NULL)
{
  // Constructor
  // Init pointers
  for(Int_t i = 0; i < 11; ++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;
  for(Int_t i = 0; i < 11; ++i) fX2Corr[i] = fY2Corr[i] = fX2Y2Corr[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 < 11; ++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 < 11; ++i) {
    fX2Corr[i] = new TProfile(Form("fX2Corr_%d",i),"",21,0,105,"s");
    fOutputList->Add(fX2Corr[i]);
    fY2Corr[i] = new TProfile(Form("fY2Corr_%d",i),"",21,0,105,"s");
    fOutputList->Add(fY2Corr[i]);
    fX2Y2Corr[i] = new TProfile(Form("fX2Y2Corr_%d",i),"",21,0,105,"s");
    fOutputList->Add(fX2Y2Corr[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]);
  }
  fPsiARawCentr = new TH2F("fPsiARawCentr","",105,0,105,100,-TMath::Pi()/2,TMath::Pi()/2);
  fOutputList->Add(fPsiARawCentr);
  fPsiAFlatCentr = new TH2F("fPsiAFlatCentr","",105,0,105,100,-TMath::Pi()/2,TMath::Pi()/2);
  fOutputList->Add(fPsiAFlatCentr);
  fPsiCRawCentr = new TH2F("fPsiCRawCentr","",105,0,105,100,-TMath::Pi()/2,TMath::Pi()/2);
  fOutputList->Add(fPsiCRawCentr);
  fPsiCFlatCentr = new TH2F("fPsiCFlatCentr","",105,0,105,100,-TMath::Pi()/2,TMath::Pi()/2);
  fOutputList->Add(fPsiCFlatCentr);
  fPsiACRawCentr = new TH2F("fPsiACRawCentr","",105,0,105,100,-TMath::Pi()/2,TMath::Pi()/2);
  fOutputList->Add(fPsiACRawCentr);
  fPsiACFlatCentr = new TH2F("fPsiACFlatCentr","",105,0,105,100,-TMath::Pi()/2,TMath::Pi()/2);
  fOutputList->Add(fPsiACFlatCentr);

  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]);
  }

  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 | AliVEvent::kCentral));
  else
    isSelected = ((esdV0->GetV0ADecision()==1) && (esdV0->GetV0CDecision()==1));

  if (!isSelected) goodEvent = kFALSE;

  AliCentrality *centrality = fEvent->GetCentrality();
  //  Float_t spdPercentile = 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 totMult = 0, totMultA = 0, totMultC = 0;
    Double_t c2A = 0, c2C = 0, s2A = 0, s2C = 0;
    Double_t qxA = 0, qyA = 0;
    Double_t qxC = 0, qyC = 0;
    for(Int_t iring = 0; iring < 8; ++iring) {
      Double_t ringMult = 0;
      Double_t c2 = 0, s2 = 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));
	}
	ringMult += fEvent->GetVZEROEqMultiplicity(iCh);
      }
      if (ringMult < 1e-6) continue;
      totMult += ringMult;
      if (iring >= 4) {
	totMultA += ringMult;
	c2A += c2;
	s2A += s2;
      }
      else {
	totMultC += ringMult;
	c2C += c2;
	s2C += s2;
      }
      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);
      fCos8Psi[iring]->Fill(spdPercentile, 2./4.*TMath::Cos(2.*4.*psiRingRaw));
      Double_t qxTierce = 0, qyTierce = 0;
      Double_t psiRingFlat = CalculateVZEROEventPlane(fEvent,iring,spdPercentile,qxTierce,qyTierce);
      fPsiRingFlatCentr[iring]->Fill(spdPercentile,psiRingFlat);
      Int_t ibin = fCos8PsiIn[iring]->FindBin(spdPercentile);
      Double_t psiRingFlatFinal = psiRingFlat + (fCos8PsiIn[iring]->GetBinContent(ibin)*TMath::Sin(2.*4.*psiRingFlat))/2.;
      if (psiRingFlatFinal > TMath::Pi()/2) psiRingFlatFinal -= TMath::Pi();
      if (psiRingFlatFinal <-TMath::Pi()/2) psiRingFlatFinal += TMath::Pi();
      fPsiRingFlatFinalCentr[iring]->Fill(spdPercentile,psiRingFlatFinal);
      fX2Corr[iring]->Fill(spdPercentile,qxTierce);
      fY2Corr[iring]->Fill(spdPercentile,qyTierce);
      fX2Y2Corr[iring]->Fill(spdPercentile,qxTierce*qyTierce);
      if (iring >= 4) {
	qxA += qxTierce;
	qyA += qyTierce;
      }
      else {
 	qxC += qxTierce;
	qyC += qyTierce;
     }
    }

    if (totMultA > 1e-6) {
      fX2[8]->Fill(spdPercentile,c2A);
      fY2[8]->Fill(spdPercentile,s2A);
      fX2Y2[8]->Fill(spdPercentile,c2A*s2A);
      Double_t psiA = TMath::ATan2(qyA,qxA)/2.;
      fPsiAFlatCentr->Fill(spdPercentile,psiA);
      //      Double_t psiAOrg = fEvent->GetEventplane()->GetEventplane("V0A",fEvent,2);
      Double_t qxEP = 0, qyEP = 0;
      Double_t psiAOrg = fEvent->GetEventplane()->CalculateVZEROEventPlane(fEvent,8,2,qxEP,qyEP);
      fX2Corr[8]->Fill(spdPercentile,qxEP);
      fY2Corr[8]->Fill(spdPercentile,qyEP);
      fX2Y2Corr[8]->Fill(spdPercentile,qxEP*qyEP);
      fPsiARawCentr->Fill(spdPercentile,psiAOrg);
      fCos8Psi[8]->Fill(spdPercentile, 2./4.*TMath::Cos(2.*4.*psiAOrg));
    }
    if (totMultC > 1e-6) {
      fX2[9]->Fill(spdPercentile,c2C);
      fY2[9]->Fill(spdPercentile,s2C);
      fX2Y2[9]->Fill(spdPercentile,c2C*s2C);
      Double_t psiC = TMath::ATan2(qyC,qxC)/2.;
      fPsiCFlatCentr->Fill(spdPercentile,psiC);
      //      Double_t psiCOrg = fEvent->GetEventplane()->GetEventplane("V0C",fEvent,2);
      Double_t qxEP = 0, qyEP = 0;
      Double_t psiCOrg = fEvent->GetEventplane()->CalculateVZEROEventPlane(fEvent,9,2,qxEP,qyEP);
      fX2Corr[9]->Fill(spdPercentile,qxEP);
      fY2Corr[9]->Fill(spdPercentile,qyEP);
      fX2Y2Corr[9]->Fill(spdPercentile,qxEP*qyEP);
      fPsiCRawCentr->Fill(spdPercentile,psiCOrg);
      fCos8Psi[9]->Fill(spdPercentile, 2./4.*TMath::Cos(2.*4.*psiCOrg));
   }
    if (totMult > 1e-6) {
      fX2[10]->Fill(spdPercentile,c2A+c2C);
      fY2[10]->Fill(spdPercentile,s2A+s2C);
      fX2Y2[10]->Fill(spdPercentile,(c2A+c2C)*(s2A+s2C));
      Double_t psiAC = TMath::ATan2(qyA+qyC,qxA+qxC)/2.;
      fPsiACFlatCentr->Fill(spdPercentile,psiAC);
      //      Double_t psiACOrg = fEvent->GetEventplane()->GetEventplane("V0",fEvent,2);
      Double_t qxEP = 0, qyEP = 0;
      Double_t psiACOrg = fEvent->GetEventplane()->CalculateVZEROEventPlane(fEvent,10,2,qxEP,qyEP);
      fX2Corr[10]->Fill(spdPercentile,qxEP);
      fY2Corr[10]->Fill(spdPercentile,qyEP);
      fX2Y2Corr[10]->Fill(spdPercentile,qxEP*qyEP);
      fPsiACRawCentr->Fill(spdPercentile,psiACOrg);
      fCos8Psi[10]->Fill(spdPercentile, 2./4.*TMath::Cos(2.*4.*psiACOrg));
    }
  }

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