[u/mrichter/AliRoot.git] / STEER / STEERBase / AliEventplane.cxx

/**************************************************************************
 * Copyright(c) 1998-2008, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/

//*****************************************************
//   Class AliEventplane
//   author: Alberica Toia, Johanna Gramling
//*****************************************************
/// A container for the event plane stored in AOD in ESD
 
#include "AliLog.h"
#include "AliEventplane.h"
#include "TVector2.h"
#include "AliVTrack.h"
#include "TObjArray.h"
#include "TArrayF.h"
#include "AliVEvent.h"
#include "AliVVZERO.h"

ClassImp(AliEventplane)

AliEventplane::AliEventplane() : TNamed("Eventplane", "Eventplane"),
  fQVector(0),
  fQContributionX(0),
  fQContributionY(0),
  fQContributionXsub1(0),
  fQContributionYsub1(0),
  fQContributionXsub2(0),
  fQContributionYsub2(0),
  fEventplaneQ(-1),
  fQsub1(0),
  fQsub2(0),
  fQsubRes(0)
{
  /// constructor
  fQContributionX = new TArrayF(0);
  fQContributionY = new TArrayF(0);
  fQContributionXsub1 = new TArrayF(0);
  fQContributionYsub1 = new TArrayF(0);
  fQContributionXsub2 = new TArrayF(0);
  fQContributionYsub2 = new TArrayF(0);
  for(Int_t i = 0; i < 11; ++i) {
    fMeanX2[i]=fMeanY2[i]=0.;
    fAPlus[i]=fAMinus[i]=1.;
    fLambdaPlus[i]=fLambdaMinus[i]=0.;
    fCos8Psi[i]=0.;
  }
}

AliEventplane::AliEventplane(const AliEventplane& ep) : 
  TNamed(ep),
  fQVector(0),
  fQContributionX(0),
  fQContributionY(0),
  fQContributionXsub1(0),
  fQContributionYsub1(0),
  fQContributionXsub2(0),
  fQContributionYsub2(0),
  fEventplaneQ(0),
  fQsub1(0),
  fQsub2(0),
  fQsubRes(0)
{
  /// Copy constructor
  ((AliEventplane &) ep).CopyEP(*this);
}

AliEventplane& AliEventplane::operator=(const AliEventplane& ep)
{
  /// Assignment operator
  if (this!=&ep){
    TNamed::operator=(ep);
    ((AliEventplane &) ep).CopyEP(*this);
  }
  return *this;
}

void AliEventplane::CopyEP(AliEventplane& ep) const 
{ // copy function

  AliEventplane& target = (AliEventplane &) ep;
  
  if(fQContributionX)
  {
    if(ep.fQContributionX)
    { 
      *(ep.fQContributionX) = *fQContributionX;
    }
    else 
    {
      ep.fQContributionX = new TArrayF(*fQContributionX);
    }
  }
  else
  {
    if(ep.fQContributionX) delete ep.fQContributionX;
    ep.fQContributionX = 0;
  }
  
  if(fQContributionY)
  {
    if(ep.fQContributionY)
    { 
      *(ep.fQContributionY) = *fQContributionY;
    }
    else 
    {
      ep.fQContributionY = new TArrayF(*fQContributionY);
    }
  }
  else
  {
    if(ep.fQContributionY) delete ep.fQContributionY;
    ep.fQContributionY = 0;
  }

  
  if(fQContributionXsub1)
  {
    if(ep.fQContributionXsub1)
    { 
      *(ep.fQContributionXsub1) = *fQContributionXsub1;
    }
    else 
    {
      ep.fQContributionXsub1 = new TArrayF(*fQContributionXsub1);
    }
  }
  else
  {
    if(ep.fQContributionXsub1) delete ep.fQContributionXsub1;
    ep.fQContributionXsub1 = 0;
  }
  
  if(fQContributionYsub1)
  {
    if(ep.fQContributionYsub1)
    { 
      *(ep.fQContributionYsub1) = *fQContributionYsub1;
    }
    else 
    {
      ep.fQContributionYsub1 = new TArrayF(*fQContributionYsub1);
    }
  }
  else
  {
    if(ep.fQContributionYsub1) delete ep.fQContributionYsub1;
    ep.fQContributionYsub1 = 0;
  }
  
  
  if(fQContributionXsub2)
  {
    if(ep.fQContributionXsub2)
    { 
      *(ep.fQContributionXsub2) = *fQContributionXsub2;
    }
    else 
    {
      ep.fQContributionXsub2 = new TArrayF(*fQContributionXsub2);
    }
  }
  else
  {
    if(ep.fQContributionXsub2) delete ep.fQContributionXsub2;
    ep.fQContributionXsub2 = 0;
  }
  
  if(fQContributionYsub2)
  {
    if(ep.fQContributionYsub2)
    { 
      *(ep.fQContributionYsub2) = *fQContributionYsub2;
    }
    else 
    {
      ep.fQContributionYsub2 = new TArrayF(*fQContributionYsub2);
    }
  }
  else
  {
    if(ep.fQContributionYsub2) delete ep.fQContributionYsub2;
    ep.fQContributionYsub2 = 0;
  }
  
  if (fEventplaneQ)
      target.fEventplaneQ = fEventplaneQ;
  
  if (fQVector)
      target.fQVector = dynamic_cast<TVector2*> (fQVector->Clone());
  if (fQsub1)
      target.fQsub1 = dynamic_cast<TVector2*> (fQsub1->Clone());
  if (fQsub2)
      target.fQsub2 = dynamic_cast<TVector2*> (fQsub2->Clone());
  if (fQsubRes)
      target.fQsubRes = fQsubRes;

  for(Int_t i = 0; i < 11; ++i) {
    target.fMeanX2[i]=fMeanX2[i];
    target.fMeanY2[i]=fMeanY2[i];
    target.fAPlus[i]=fAPlus[i];
    target.fAMinus[i]=fAMinus[i];
    target.fLambdaPlus[i]=fLambdaPlus[i];
    target.fLambdaMinus[i]=fLambdaMinus[i];
    target.fCos8Psi[i]=fCos8Psi[i];
  }
}

AliEventplane::~AliEventplane()
{
  /// destructor
  if (fQContributionX){
      delete fQContributionX;
      fQContributionX = 0;
  }
  if (fQContributionY){
      delete fQContributionY;
      fQContributionY = 0;
  }
  if (fQContributionXsub1){
      delete fQContributionXsub1;
      fQContributionXsub1 = 0;
  }
  if (fQContributionYsub1){
      delete fQContributionYsub1;
      fQContributionYsub1 = 0;
  }
  if (fQContributionXsub2){
      delete fQContributionXsub2;
      fQContributionXsub2 = 0;
  }
  if (fQContributionYsub2){
      delete fQContributionYsub2;
      fQContributionYsub2 = 0;
  }
  if (fQVector){
      delete fQVector;
      fQVector = 0;
  }
  if (fQsub1){
      delete fQsub1;
      fQsub1 = 0;
  }
    if (fQsub2){
      delete fQsub2;
      fQsub2 = 0;
  }
}

TVector2* AliEventplane::GetQVector()
{
  return fQVector;
}

Double_t AliEventplane::GetEventplane(const char *x, const AliVEvent *event, Int_t harmonic) const
{
  TString method = x;
  Double_t qx = 0, qy = 0;
  if(method.CompareTo("Q")==0)      return fEventplaneQ;
  else if(method.CompareTo("V0A")==0) return CalculateVZEROEventPlane(event, 8, harmonic, qx, qy);
  else if(method.CompareTo("V0C")==0) return CalculateVZEROEventPlane(event, 9, harmonic, qx, qy);
  else if(method.CompareTo("V0")==0)  return CalculateVZEROEventPlane(event,10, harmonic, qx, qy);

  return -1000.;
}

Double_t AliEventplane::CalculateVZEROEventPlane(const AliVEvent *  event, Int_t firstRing, Int_t lastRing, Int_t harmonic, Double_t &qxTot, Double_t &qyTot) const
{
  // Calculate the VZERO event-plane by summing rings
  // The flatenning is done on every ring separately
  if(!event) {
    AliError("No Event received");
    return -1000.;
  }
  AliVVZERO *vzeroData = event->GetVZEROData();
  if(!vzeroData) {
    AliError("Enable to get VZERO Data");
    return -1000.;
  }
  if(harmonic <= 0) {
    AliError("Required harmonic is less or equal to 0");
    return -1000.;
  }

  qxTot=qyTot=0.;
  Double_t qx, qy;
  Double_t totMult = 0.;
  for(Int_t ring = firstRing; ring <=lastRing; ++ring) {
    qx=qy=0.;
    Double_t multRing = 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(harmonic*phi);
      qy += mult*TMath::Sin(harmonic*phi);
      multRing += mult;
    }

    if (multRing < 1e-6) continue;
    totMult += multRing;

    if (harmonic == 2) {
      // Recentering
      Double_t qxPrime = qx - fMeanX2[ring];
      Double_t qyPrime = qy - fMeanY2[ring];
      // Twist
      Double_t trans[2][2];
      trans[0][0] = 1./(1.-fLambdaPlus[ring]*fLambdaMinus[ring]);
      trans[0][1] = -fLambdaMinus[ring]/(1.-fLambdaPlus[ring]*fLambdaMinus[ring]);
      trans[1][0] = -fLambdaPlus[ring]/(1.-fLambdaPlus[ring]*fLambdaMinus[ring]);
      trans[1][1] = 1./(1.-fLambdaPlus[ring]*fLambdaMinus[ring]);
      Double_t qxSeconde = trans[0][0]*qxPrime + trans[0][1]*qyPrime;
      Double_t qySeconde = trans[1][0]*qxPrime + trans[1][1]*qyPrime;
      // Rescaling
      Double_t qxTierce = qxSeconde/fAPlus[ring];
      Double_t qyTierce = qySeconde/fAMinus[ring];
      // 4th Fourier momenta in order to flatten the EP within a sector
      Double_t psi = TMath::ATan2(qyTierce,qxTierce)/harmonic;
      Double_t deltaPsi = (fCos8Psi[ring]*TMath::Sin(2.*4.*psi))/2.;
      Double_t qxQuarte = qxTierce*TMath::Cos(2.*deltaPsi) - qyTierce*TMath::Sin(2.*deltaPsi);
      Double_t qyQuarte = qxTierce*TMath::Sin(2.*deltaPsi) + qyTierce*TMath::Cos(2.*deltaPsi);
      qxTot += qxQuarte;
      qyTot += qyQuarte;
    }
    else {
      qxTot += qx;
      qyTot += qy;
    }
  }

  // In case of no hits return an invalid event-plane angle
  if (totMult<1e-6) return -999;

  return (TMath::ATan2(qyTot,qxTot)/harmonic);
}

Double_t AliEventplane::CalculateVZEROEventPlane(const AliVEvent *  event, Int_t ring, Int_t harmonic, Double_t &qx, Double_t &qy) const
{
  // Calculate the VZERO event-plane from an individual ring
  // Ring 8 - total V0A (rings 4-7), Ring 9 - total V0C (rings 0-3)
  // Ring 10 - total V0 (rings 0-7)
  if(!event) {
    AliError("No Event received");
    return -1000.;
  }
  AliVVZERO *vzeroData = event->GetVZEROData();
  if(!vzeroData) {
    AliError("Enable to get VZERO Data");
    return -1000.;
  }
  if(harmonic <= 0) {
    AliError("Required harmonic is less or equal to 0");
    return -1000.;
  }

  Int_t firstCh=-1,lastCh=-1;
  if (ring < 8) {
    firstCh = ring*8;
    lastCh = (ring+1)*8;
  }
  else if (ring == 8) {
    firstCh = 32;
    lastCh = 64;
  }
  else if (ring == 9) {
    firstCh = 0;
    lastCh = 32;
  }
  else if (ring == 10) {
    firstCh = 0;
    lastCh = 64;
  }
  qx=qy=0.;
  Double_t multRing = 0.;
  for(Int_t iCh = firstCh; iCh < lastCh; ++iCh) {
    Double_t phi = TMath::Pi()/8. + TMath::Pi()/4.*(iCh%8);
    Double_t mult = event->GetVZEROEqMultiplicity(iCh);
    qx += mult*TMath::Cos(harmonic*phi);
    qy += mult*TMath::Sin(harmonic*phi);
    multRing += mult;
  }

  // In case of no hits return an invalid event-plane angle
  if (multRing < 1e-6) return -999;

  if (harmonic == 2) {
    // Recentering
    Double_t qxPrime = qx - fMeanX2[ring];
    Double_t qyPrime = qy - fMeanY2[ring];
    // Twist
    Double_t trans[2][2];
    trans[0][0] = 1./(1.-fLambdaPlus[ring]*fLambdaMinus[ring]);
    trans[0][1] = -fLambdaMinus[ring]/(1.-fLambdaPlus[ring]*fLambdaMinus[ring]);
    trans[1][0] = -fLambdaPlus[ring]/(1.-fLambdaPlus[ring]*fLambdaMinus[ring]);
    trans[1][1] = 1./(1.-fLambdaPlus[ring]*fLambdaMinus[ring]);
    Double_t qxSeconde = trans[0][0]*qxPrime + trans[0][1]*qyPrime;
    Double_t qySeconde = trans[1][0]*qxPrime + trans[1][1]*qyPrime;
    // Rescaling
    Double_t qxTierce = qxSeconde/fAPlus[ring];
    Double_t qyTierce = qySeconde/fAMinus[ring];
    // 4th Fourier momenta in order to flatten the EP within a sector
    Double_t psi = TMath::ATan2(qyTierce,qxTierce)/harmonic;
    Double_t deltaPsi = (fCos8Psi[ring]*TMath::Sin(2.*4.*psi))/2.;
    Double_t qxQuarte = qxTierce*TMath::Cos(2.*deltaPsi) - qyTierce*TMath::Sin(2.*deltaPsi);
    Double_t qyQuarte = qxTierce*TMath::Sin(2.*deltaPsi) + qyTierce*TMath::Cos(2.*deltaPsi);
    qx = qxQuarte;
    qy = qyQuarte;
  }

  return (TMath::ATan2(qy,qx)/harmonic);
}

void AliEventplane::SetVZEROEPParams(Int_t ring,
				     Double_t meanX2, Double_t meanY2,
				     Double_t aPlus, Double_t aMinus,
				     Double_t lambdaPlus, Double_t lambdaMinus,
				     Double_t cos8Psi)
{
  // Set the VZERO event-plane
  // flatenning parameters
  if (ring < 0 || ring >= 11) AliFatal(Form("Bad ring index (%d)",ring));

  fMeanX2[ring] = meanX2;
  fMeanY2[ring] = meanY2;
  fAPlus[ring] = aPlus;
  fAMinus[ring] = aMinus;
  fLambdaPlus[ring] = lambdaPlus;
  fLambdaMinus[ring] = lambdaMinus;
  fCos8Psi[ring] = cos8Psi;
}

TVector2* AliEventplane::GetQsub1()
{
  return fQsub1;
}

TVector2* AliEventplane::GetQsub2()
{
  return fQsub2;
}

Double_t AliEventplane::GetQsubRes()
{
  return fQsubRes;
}

Bool_t AliEventplane::IsEventInEventplaneClass(Double_t a, Double_t b, const char *x)
{
  TString method = x;
  if ((method.CompareTo("Q")==0) && (fEventplaneQ >=a && fEventplaneQ < b)) return kTRUE;
  else return kFALSE;
}

Double_t AliEventplane::GetQContributionX(AliVTrack* track)
{ 
  return fQContributionX->GetAt(track->GetID());
}

Double_t AliEventplane::GetQContributionY(AliVTrack* track)
{ 
  return fQContributionY->GetAt(track->GetID());
}

Double_t AliEventplane::GetQContributionXsub1(AliVTrack* track)
{ 
  return fQContributionXsub1->GetAt(track->GetID());
}

Double_t AliEventplane::GetQContributionYsub1(AliVTrack* track)
{ 
  return fQContributionYsub1->GetAt(track->GetID());
}

Double_t AliEventplane::GetQContributionXsub2(AliVTrack* track)
{ 
  return fQContributionXsub2->GetAt(track->GetID());
}

Double_t AliEventplane::GetQContributionYsub2(AliVTrack* track)
{ 
  return fQContributionYsub2->GetAt(track->GetID());
}

void AliEventplane::Reset()
{ 
  delete fQVector; fQVector=0;
  fQContributionX->Reset();
  fQContributionY->Reset();
  fQContributionXsub1->Reset();
  fQContributionYsub1->Reset();
  fQContributionXsub2->Reset();
  fQContributionYsub2->Reset();
  fEventplaneQ = -1;
  delete fQsub1; fQsub1=0;
  delete fQsub2; fQsub2=0;
  fQsubRes = 0;
}
Commit	Line	Data
ce7adfe9	1	/**************************************************************************
	2	* Copyright(c) 1998-2008, ALICE Experiment at CERN, All rights reserved. *
	3	* *
	4	* Author: The ALICE Off-line Project. *
	5	* Contributors are mentioned in the code where appropriate. *
	6	* *
	7	* Permission to use, copy, modify and distribute this software and its *
	8	* documentation strictly for non-commercial purposes is hereby granted *
	9	* without fee, provided that the above copyright notice appears in all *
	10	* copies and that both the copyright notice and this permission notice *
	11	* appear in the supporting documentation. The authors make no claims *
	12	* about the suitability of this software for any purpose. It is *
	13	* provided "as is" without express or implied warranty. *
	14	**************************************************************************/
	15
	16	//*****************************************************
	17	// Class AliEventplane
	18	// author: Alberica Toia, Johanna Gramling
	19	//*****************************************************
	20	/// A container for the event plane stored in AOD in ESD
	21
fd5c821a	22	#include "AliLog.h"
ce7adfe9	23	#include "AliEventplane.h"
ce7adfe9	24	#include "TVector2.h"
aa260dd1	25	#include "AliVTrack.h"
ce7adfe9	26	#include "TObjArray.h"
ce7adfe9	27	#include "TArrayF.h"
fd5c821a	28	#include "AliVEvent.h"
fd5c821a	29	#include "AliVVZERO.h"
ce7adfe9	30
	31	ClassImp(AliEventplane)
	32
	33	AliEventplane::AliEventplane() : TNamed("Eventplane", "Eventplane"),
	34	fQVector(0),
	35	fQContributionX(0),
	36	fQContributionY(0),
e51055a0	37	fQContributionXsub1(0),
	38	fQContributionYsub1(0),
	39	fQContributionXsub2(0),
	40	fQContributionYsub2(0),
3e64f60c	41	fEventplaneQ(-1),
ce7adfe9	42	fQsub1(0),
	43	fQsub2(0),
	44	fQsubRes(0)
	45	{
	46	/// constructor
3e64f60c	47	fQContributionX = new TArrayF(0);
3e64f60c	48	fQContributionY = new TArrayF(0);
e51055a0	49	fQContributionXsub1 = new TArrayF(0);
	50	fQContributionYsub1 = new TArrayF(0);
	51	fQContributionXsub2 = new TArrayF(0);
	52	fQContributionYsub2 = new TArrayF(0);
cc4d0e9c	53	for(Int_t i = 0; i < 11; ++i) {
0643f5a7	54	fMeanX2[i]=fMeanY2[i]=0.;
	55	fAPlus[i]=fAMinus[i]=1.;
	56	fLambdaPlus[i]=fLambdaMinus[i]=0.;
	57	fCos8Psi[i]=0.;
	58	}
ce7adfe9	59	}
	60
	61	AliEventplane::AliEventplane(const AliEventplane& ep) :
fe536aed	62	TNamed(ep),
ce7adfe9	63	fQVector(0),
	64	fQContributionX(0),
	65	fQContributionY(0),
e51055a0	66	fQContributionXsub1(0),
	67	fQContributionYsub1(0),
	68	fQContributionXsub2(0),
	69	fQContributionYsub2(0),
ce7adfe9	70	fEventplaneQ(0),
	71	fQsub1(0),
	72	fQsub2(0),
	73	fQsubRes(0)
	74	{
	75	/// Copy constructor
	76	((AliEventplane &) ep).CopyEP(*this);
	77	}
	78
	79	AliEventplane& AliEventplane::operator=(const AliEventplane& ep)
	80	{
	81	/// Assignment operator
fe536aed	82	if (this!=&ep){
fe536aed	83	TNamed::operator=(ep);
ce7adfe9	84	((AliEventplane &) ep).CopyEP(*this);
fe536aed	85	}
ce7adfe9	86	return *this;
	87	}
	88
fe536aed	89	void AliEventplane::CopyEP(AliEventplane& ep) const
ce7adfe9	90	{ // copy function
	91
	92	AliEventplane& target = (AliEventplane &) ep;
fe536aed	93
	94	if(fQContributionX)
	95	{
	96	if(ep.fQContributionX)
	97	{
	98	(ep.fQContributionX) = fQContributionX;
	99	}
	100	else
	101	{
	102	ep.fQContributionX = new TArrayF(*fQContributionX);
	103	}
	104	}
	105	else
	106	{
	107	if(ep.fQContributionX) delete ep.fQContributionX;
	108	ep.fQContributionX = 0;
	109	}
	110
	111	if(fQContributionY)
	112	{
	113	if(ep.fQContributionY)
	114	{
	115	(ep.fQContributionY) = fQContributionY;
	116	}
	117	else
	118	{
	119	ep.fQContributionY = new TArrayF(*fQContributionY);
	120	}
	121	}
	122	else
	123	{
	124	if(ep.fQContributionY) delete ep.fQContributionY;
	125	ep.fQContributionY = 0;
	126	}
	127
	128
	129	if(fQContributionXsub1)
	130	{
	131	if(ep.fQContributionXsub1)
	132	{
	133	(ep.fQContributionXsub1) = fQContributionXsub1;
	134	}
	135	else
	136	{
	137	ep.fQContributionXsub1 = new TArrayF(*fQContributionXsub1);
	138	}
	139	}
	140	else
	141	{
	142	if(ep.fQContributionXsub1) delete ep.fQContributionXsub1;
	143	ep.fQContributionXsub1 = 0;
	144	}
	145
	146	if(fQContributionYsub1)
	147	{
	148	if(ep.fQContributionYsub1)
	149	{
	150	(ep.fQContributionYsub1) = fQContributionYsub1;
	151	}
	152	else
	153	{
	154	ep.fQContributionYsub1 = new TArrayF(*fQContributionYsub1);
	155	}
	156	}
157	else
158	{
159	if(ep.fQContributionYsub1) delete ep.fQContributionYsub1;
160	ep.fQContributionYsub1 = 0;
161	}
162
163
164	if(fQContributionXsub2)
165	{
166	if(ep.fQContributionXsub2)
167	{
168	(ep.fQContributionXsub2) = fQContributionXsub2;
169	}
170	else
171	{
172	ep.fQContributionXsub2 = new TArrayF(*fQContributionXsub2);
173	}
174	}
175	else
176	{
177	if(ep.fQContributionXsub2) delete ep.fQContributionXsub2;
178	ep.fQContributionXsub2 = 0;
179	}
180
181	if(fQContributionYsub2)
182	{
183	if(ep.fQContributionYsub2)
184	{
185	(ep.fQContributionYsub2) = fQContributionYsub2;
186	}
187	else
188	{
189	ep.fQContributionYsub2 = new TArrayF(*fQContributionYsub2);
190	}
191	}
192	else
193	{
194	if(ep.fQContributionYsub2) delete ep.fQContributionYsub2;
195	ep.fQContributionYsub2 = 0;
196	}
197
ce7adfe9	198	if (fEventplaneQ)
ce7adfe9	199	target.fEventplaneQ = fEventplaneQ;
fe536aed	200
ce7adfe9	201	if (fQVector)
	202	target.fQVector = dynamic_cast<TVector2*> (fQVector->Clone());
	203	if (fQsub1)
	204	target.fQsub1 = dynamic_cast<TVector2*> (fQsub1->Clone());
	205	if (fQsub2)
	206	target.fQsub2 = dynamic_cast<TVector2*> (fQsub2->Clone());
	207	if (fQsubRes)
	208	target.fQsubRes = fQsubRes;
0643f5a7	209
cc4d0e9c	210	for(Int_t i = 0; i < 11; ++i) {
0643f5a7	211	target.fMeanX2[i]=fMeanX2[i];
	212	target.fMeanY2[i]=fMeanY2[i];
	213	target.fAPlus[i]=fAPlus[i];
	214	target.fAMinus[i]=fAMinus[i];
	215	target.fLambdaPlus[i]=fLambdaPlus[i];
	216	target.fLambdaMinus[i]=fLambdaMinus[i];
	217	target.fCos8Psi[i]=fCos8Psi[i];
	218	}
ce7adfe9	219	}
	220
	221	AliEventplane::~AliEventplane()
	222	{
	223	/// destructor
	224	if (fQContributionX){
	225	delete fQContributionX;
	226	fQContributionX = 0;
	227	}
	228	if (fQContributionY){
	229	delete fQContributionY;
	230	fQContributionY = 0;
	231	}
e51055a0	232	if (fQContributionXsub1){
	233	delete fQContributionXsub1;
	234	fQContributionXsub1 = 0;
	235	}
	236	if (fQContributionYsub1){
	237	delete fQContributionYsub1;
	238	fQContributionYsub1 = 0;
	239	}
	240	if (fQContributionXsub2){
	241	delete fQContributionXsub2;
	242	fQContributionXsub2 = 0;
	243	}
	244	if (fQContributionYsub2){
	245	delete fQContributionYsub2;
	246	fQContributionYsub2 = 0;
	247	}
ce7adfe9	248	if (fQVector){
	249	delete fQVector;
	250	fQVector = 0;
	251	}
	252	if (fQsub1){
	253	delete fQsub1;
	254	fQsub1 = 0;
	255	}
	256	if (fQsub2){
	257	delete fQsub2;
	258	fQsub2 = 0;
	259	}
	260	}
	261
	262	TVector2* AliEventplane::GetQVector()
	263	{
	264	return fQVector;
	265	}
	266
fd5c821a	267	Double_t AliEventplane::GetEventplane(const char x, const AliVEvent event, Int_t harmonic) const
ce7adfe9	268	{
ce7adfe9	269	TString method = x;
fe509094	270	Double_t qx = 0, qy = 0;
ce7adfe9	271	if(method.CompareTo("Q")==0) return fEventplaneQ;
cc4d0e9c	272	else if(method.CompareTo("V0A")==0) return CalculateVZEROEventPlane(event, 8, harmonic, qx, qy);
	273	else if(method.CompareTo("V0C")==0) return CalculateVZEROEventPlane(event, 9, harmonic, qx, qy);
	274	else if(method.CompareTo("V0")==0) return CalculateVZEROEventPlane(event,10, harmonic, qx, qy);
fd5c821a	275
	276	return -1000.;
	277	}
	278
fe509094	279	Double_t AliEventplane::CalculateVZEROEventPlane(const AliVEvent * event, Int_t firstRing, Int_t lastRing, Int_t harmonic, Double_t &qxTot, Double_t &qyTot) const
fd5c821a	280	{
cc4d0e9c	281	// Calculate the VZERO event-plane by summing rings
cc4d0e9c	282	// The flatenning is done on every ring separately
fd5c821a	283	if(!event) {
	284	AliError("No Event received");
	285	return -1000.;
	286	}
	287	AliVVZERO *vzeroData = event->GetVZEROData();
	288	if(!vzeroData) {
	289	AliError("Enable to get VZERO Data");
	290	return -1000.;
	291	}
	292	if(harmonic <= 0) {
	293	AliError("Required harmonic is less or equal to 0");
	294	return -1000.;
	295	}
	296
fe509094	297	qxTot=qyTot=0.;
0643f5a7	298	Double_t qx, qy;
	299	Double_t totMult = 0.;
	300	for(Int_t ring = firstRing; ring <=lastRing; ++ring) {
	301	qx=qy=0.;
0643f5a7	302	Double_t multRing = 0.;
	303	for(Int_t iCh = ring8; iCh < (ring+1)8; ++iCh) {
	304	Double_t phi = TMath::Pi()/8. + TMath::Pi()/4.*(iCh%8);
	305	Double_t mult = event->GetVZEROEqMultiplicity(iCh);
	306	qx += multTMath::Cos(harmonicphi);
	307	qy += multTMath::Sin(harmonicphi);
	308	multRing += mult;
	309	}
	310
	311	if (multRing < 1e-6) continue;
	312	totMult += multRing;
	313
	314	if (harmonic == 2) {
	315	// Recentering
	316	Double_t qxPrime = qx - fMeanX2[ring];
	317	Double_t qyPrime = qy - fMeanY2[ring];
	318	// Twist
cc4d0e9c	319	Double_t trans[2][2];
	320	trans[0][0] = 1./(1.-fLambdaPlus[ring]*fLambdaMinus[ring]);
	321	trans[0][1] = -fLambdaMinus[ring]/(1.-fLambdaPlus[ring]*fLambdaMinus[ring]);
	322	trans[1][0] = -fLambdaPlus[ring]/(1.-fLambdaPlus[ring]*fLambdaMinus[ring]);
	323	trans[1][1] = 1./(1.-fLambdaPlus[ring]*fLambdaMinus[ring]);
0643f5a7	324	Double_t qxSeconde = trans[0][0]qxPrime + trans[0][1]qyPrime;
	325	Double_t qySeconde = trans[1][0]qxPrime + trans[1][1]qyPrime;
	326	// Rescaling
	327	Double_t qxTierce = qxSeconde/fAPlus[ring];
	328	Double_t qyTierce = qySeconde/fAMinus[ring];
93736b53	329	// 4th Fourier momenta in order to flatten the EP within a sector
	330	Double_t psi = TMath::ATan2(qyTierce,qxTierce)/harmonic;
	331	Double_t deltaPsi = (fCos8Psi[ring]TMath::Sin(2.4.*psi))/2.;
	332	Double_t qxQuarte = qxTierceTMath::Cos(2.deltaPsi) - qyTierceTMath::Sin(2.deltaPsi);
	333	Double_t qyQuarte = qxTierceTMath::Sin(2.deltaPsi) + qyTierceTMath::Cos(2.deltaPsi);
	334	qxTot += qxQuarte;
	335	qyTot += qyQuarte;
0643f5a7	336	}
	337	else {
	338	qxTot += qx;
	339	qyTot += qy;
	340	}
	341	}
fd5c821a	342
0643f5a7	343	// In case of no hits return an invalid event-plane angle
0643f5a7	344	if (totMult<1e-6) return -999;
fd5c821a	345
93736b53	346	return (TMath::ATan2(qyTot,qxTot)/harmonic);
ce7adfe9	347	}
ce7adfe9	348
cc4d0e9c	349	Double_t AliEventplane::CalculateVZEROEventPlane(const AliVEvent * event, Int_t ring, Int_t harmonic, Double_t &qx, Double_t &qy) const
	350	{
	351	// Calculate the VZERO event-plane from an individual ring
	352	// Ring 8 - total V0A (rings 4-7), Ring 9 - total V0C (rings 0-3)
	353	// Ring 10 - total V0 (rings 0-7)
	354	if(!event) {
	355	AliError("No Event received");
	356	return -1000.;
	357	}
	358	AliVVZERO *vzeroData = event->GetVZEROData();
	359	if(!vzeroData) {
	360	AliError("Enable to get VZERO Data");
	361	return -1000.;
	362	}
	363	if(harmonic <= 0) {
	364	AliError("Required harmonic is less or equal to 0");
	365	return -1000.;
	366	}
	367
6e4f0e6b	368	Int_t firstCh=-1,lastCh=-1;
	369	if (ring < 8) {
	370	firstCh = ring*8;
	371	lastCh = (ring+1)*8;
	372	}
	373	else if (ring == 8) {
	374	firstCh = 32;
	375	lastCh = 64;
	376	}
	377	else if (ring == 9) {
	378	firstCh = 0;
	379	lastCh = 32;
	380	}
	381	else if (ring == 10) {
	382	firstCh = 0;
	383	lastCh = 64;
	384	}
cc4d0e9c	385	qx=qy=0.;
cc4d0e9c	386	Double_t multRing = 0.;
6e4f0e6b	387	for(Int_t iCh = firstCh; iCh < lastCh; ++iCh) {
cc4d0e9c	388	Double_t phi = TMath::Pi()/8. + TMath::Pi()/4.*(iCh%8);
	389	Double_t mult = event->GetVZEROEqMultiplicity(iCh);
	390	qx += multTMath::Cos(harmonicphi);
	391	qy += multTMath::Sin(harmonicphi);
	392	multRing += mult;
	393	}
	394
	395	// In case of no hits return an invalid event-plane angle
	396	if (multRing < 1e-6) return -999;
	397
	398	if (harmonic == 2) {
	399	// Recentering
	400	Double_t qxPrime = qx - fMeanX2[ring];
	401	Double_t qyPrime = qy - fMeanY2[ring];
	402	// Twist
	403	Double_t trans[2][2];
	404	trans[0][0] = 1./(1.-fLambdaPlus[ring]*fLambdaMinus[ring]);
	405	trans[0][1] = -fLambdaMinus[ring]/(1.-fLambdaPlus[ring]*fLambdaMinus[ring]);
	406	trans[1][0] = -fLambdaPlus[ring]/(1.-fLambdaPlus[ring]*fLambdaMinus[ring]);
	407	trans[1][1] = 1./(1.-fLambdaPlus[ring]*fLambdaMinus[ring]);
	408	Double_t qxSeconde = trans[0][0]qxPrime + trans[0][1]qyPrime;
	409	Double_t qySeconde = trans[1][0]qxPrime + trans[1][1]qyPrime;
	410	// Rescaling
	411	Double_t qxTierce = qxSeconde/fAPlus[ring];
	412	Double_t qyTierce = qySeconde/fAMinus[ring];
	413	// 4th Fourier momenta in order to flatten the EP within a sector
	414	Double_t psi = TMath::ATan2(qyTierce,qxTierce)/harmonic;
	415	Double_t deltaPsi = (fCos8Psi[ring]TMath::Sin(2.4.*psi))/2.;
	416	Double_t qxQuarte = qxTierceTMath::Cos(2.deltaPsi) - qyTierceTMath::Sin(2.deltaPsi);
	417	Double_t qyQuarte = qxTierceTMath::Sin(2.deltaPsi) + qyTierceTMath::Cos(2.deltaPsi);
	418	qx = qxQuarte;
	419	qy = qyQuarte;
	420	}
	421
	422	return (TMath::ATan2(qy,qx)/harmonic);
	423	}
	424
0643f5a7	425	void AliEventplane::SetVZEROEPParams(Int_t ring,
	426	Double_t meanX2, Double_t meanY2,
	427	Double_t aPlus, Double_t aMinus,
	428	Double_t lambdaPlus, Double_t lambdaMinus,
	429	Double_t cos8Psi)
	430	{
	431	// Set the VZERO event-plane
	432	// flatenning parameters
cc4d0e9c	433	if (ring < 0 \|\| ring >= 11) AliFatal(Form("Bad ring index (%d)",ring));
0643f5a7	434
	435	fMeanX2[ring] = meanX2;
	436	fMeanY2[ring] = meanY2;
	437	fAPlus[ring] = aPlus;
	438	fAMinus[ring] = aMinus;
	439	fLambdaPlus[ring] = lambdaPlus;
	440	fLambdaMinus[ring] = lambdaMinus;
	441	fCos8Psi[ring] = cos8Psi;
	442	}
fd5c821a	443
ce7adfe9	444	TVector2* AliEventplane::GetQsub1()
	445	{
	446	return fQsub1;
	447	}
	448
	449	TVector2* AliEventplane::GetQsub2()
	450	{
	451	return fQsub2;
	452	}
	453
	454	Double_t AliEventplane::GetQsubRes()
	455	{
	456	return fQsubRes;
	457	}
	458
	459	Bool_t AliEventplane::IsEventInEventplaneClass(Double_t a, Double_t b, const char *x)
	460	{
	461	TString method = x;
	462	if ((method.CompareTo("Q")==0) && (fEventplaneQ >=a && fEventplaneQ < b)) return kTRUE;
	463	else return kFALSE;
	464	}
	465
aa260dd1	466	Double_t AliEventplane::GetQContributionX(AliVTrack* track)
ce7adfe9	467	{
	468	return fQContributionX->GetAt(track->GetID());
	469	}
	470
aa260dd1	471	Double_t AliEventplane::GetQContributionY(AliVTrack* track)
ce7adfe9	472	{
	473	return fQContributionY->GetAt(track->GetID());
	474	}
3e64f60c	475
e51055a0	476	Double_t AliEventplane::GetQContributionXsub1(AliVTrack* track)
	477	{
	478	return fQContributionXsub1->GetAt(track->GetID());
	479	}
	480
	481	Double_t AliEventplane::GetQContributionYsub1(AliVTrack* track)
	482	{
	483	return fQContributionYsub1->GetAt(track->GetID());
	484	}
	485
	486	Double_t AliEventplane::GetQContributionXsub2(AliVTrack* track)
	487	{
	488	return fQContributionXsub2->GetAt(track->GetID());
	489	}
	490
	491	Double_t AliEventplane::GetQContributionYsub2(AliVTrack* track)
	492	{
	493	return fQContributionYsub2->GetAt(track->GetID());
	494	}
	495
3e64f60c	496	void AliEventplane::Reset()
	497	{
	498	delete fQVector; fQVector=0;
	499	fQContributionX->Reset();
	500	fQContributionY->Reset();
e51055a0	501	fQContributionXsub1->Reset();
	502	fQContributionYsub1->Reset();
	503	fQContributionXsub2->Reset();
	504	fQContributionYsub2->Reset();
3e64f60c	505	fEventplaneQ = -1;
	506	delete fQsub1; fQsub1=0;
	507	delete fQsub2; fQsub2=0;
	508	fQsubRes = 0;
	509	}