[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 "AliESDEvent.h"
#include "AliESDInputHandler.h"
#include "AliAnaVZEROEPFlatenning.h"
#include "AliCentrality.h"
#include "AliEventplane.h"

// VZERO includes
#include "AliESDVZERO.h"

ClassImp(AliAnaVZEROEPFlatenning)

AliAnaVZEROEPFlatenning::AliAnaVZEROEPFlatenning() 
  : AliAnalysisTaskSE("AliAnaVZEROEPFlatenning"), fESD(0), fOutputList(0),
  fMBTrigName("CPBI"),
  fUsePhysSel(kFALSE)
{
  // 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;

  // 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), fESD(0), fOutputList(0),
  fMBTrigName("CPBI"),
  fUsePhysSel(kFALSE)
{
  // 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;

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

  PostData(1, fOutputList);
}

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

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


  fESD = dynamic_cast<AliESDEvent*>(InputEvent());
  if (!fESD) {
    printf("ERROR: fESD not available\n");
    return;
  }

  AliESDVZERO* esdV0 = fESD->GetVZEROData();
  if (!esdV0) {
    Printf("ERROR: esd V0  not available");
    return;
  }

  // Trigger
  TString trigStr(fESD->GetFiredTriggerClasses());
  if (!trigStr.Contains("-B-")) return;
  if (!trigStr.Contains(fMBTrigName.Data())) return;

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

  if (!isSelected) goodEvent = kFALSE;

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

  const AliESDVertex *primaryVtx = fESD->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;

  if (goodEvent) {
    for(Int_t iring = 0; iring < 8; ++iring) {
      Double_t c2 = 0;
      Double_t s2 = 0;
      Double_t c4 = 0;
      Double_t s4 = 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 += fESD->GetVZEROEqMultiplicity(iCh)*TMath::Cos(2.*phi);
	s2 += fESD->GetVZEROEqMultiplicity(iCh)*TMath::Sin(2.*phi);
	c4 += fESD->GetVZEROEqMultiplicity(iCh)*TMath::Cos(4.*phi);
	s4 += fESD->GetVZEROEqMultiplicity(iCh);
	totMult += fESD->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 psiRingRaw = fESD->GetEventplane()->CalculateVZEROEventPlane(fESD,iring,iring,2);
      fPsiRingRawCentr[iring]->Fill(spdPercentile,psiRingRaw);
      Double_t psiRingFlat2 = CalculateVZEROEventPlane(fESD,iring,spdPercentile);
      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);
    }
  }

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