#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)
{
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)
{
((AliEventplane &) ep).CopyEP(*this);
}
AliEventplane& AliEventplane::operator=(const AliEventplane& ep)
{
if (this!=&ep){
TNamed::operator=(ep);
((AliEventplane &) ep).CopyEP(*this);
}
return *this;
}
void AliEventplane::CopyEP(AliEventplane& ep) const
{
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()
{
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
{
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) {
Double_t qxPrime = qx - fMeanX2[ring];
Double_t qyPrime = qy - fMeanY2[ring];
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;
Double_t qxTierce = qxSeconde/fAPlus[ring];
Double_t qyTierce = qySeconde/fAMinus[ring];
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;
}
}
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
{
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;
}
if (multRing < 1e-6) return -999;
if (harmonic == 2) {
Double_t qxPrime = qx - fMeanX2[ring];
Double_t qyPrime = qy - fMeanY2[ring];
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;
Double_t qxTierce = qxSeconde/fAPlus[ring];
Double_t qyTierce = qySeconde/fAMinus[ring];
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)
{
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;
}