[u/mrichter/AliRoot.git] / PWG2 / FLOW / AliFlowAnalysisWithLYZEventPlane.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.                  * 
**************************************************************************/

/*
$Log$
*/ 

//#define AliFlowAnalysisWithLYZEventPlane_cxx
 
#include "Riostream.h"  //needed as include

#include "TFile.h"
#include "TList.h"
#include "TComplex.h"   //needed as include
#include "TCanvas.h"   //needed as include
#include "TLegend.h"   //needed as include
#include "TProfile.h"  //needed as include
#include "TVector2.h"

class TH1F;

#include "AliFlowLYZConstants.h"    //needed as include
#include "AliFlowCommonConstants.h" //needed as include
#include "AliFlowEventSimple.h"
#include "AliFlowTrackSimple.h"
#include "AliFlowCommonHist.h"
#include "AliFlowCommonHistResults.h"
#include "AliFlowLYZEventPlane.h"
#include "AliFlowAnalysisWithLYZEventPlane.h"

class AliFlowVector;

// AliFlowAnalysisWithLYZEventPlane:
//
// Class to do flow analysis with the event plane from the LYZ method
//
// author: N. van der Kolk (kolk@nikhef.nl)


ClassImp(AliFlowAnalysisWithLYZEventPlane)

  //-----------------------------------------------------------------------
 
 AliFlowAnalysisWithLYZEventPlane::AliFlowAnalysisWithLYZEventPlane():
   fHistList(NULL),
   fSecondRunList(NULL),
   fSecondReDtheta(NULL),
   fSecondImDtheta(NULL),
   fFirstr0theta(NULL),
   fHistProFlow(NULL),
   fHistProFlow2(NULL),
   fHistProWr(NULL),
   fHistProWrCorr(NULL),
   fHistQsumforChi(NULL),
   fHistDeltaPhi(NULL),
   fHistDeltaPhi2(NULL),
   fHistDeltaPhihere(NULL),
   fHistPhiEP(NULL),
   fHistPhiEPhere(NULL),
   fHistPhiLYZ(NULL),
   fHistPhiLYZ2(NULL),
   fCommonHists(NULL),
   fCommonHistsRes(NULL),
   fEventNumber(0),
   fQsum(NULL),
   fQ2sum(0)
{

  // Constructor.
  fQsum = new TVector2();        // flow vector sum

  fHistList = new TList();
  fSecondRunList = new TList();
}

 
 //-----------------------------------------------------------------------


 AliFlowAnalysisWithLYZEventPlane::~AliFlowAnalysisWithLYZEventPlane() 
 {
   //destructor
   delete fQsum;
   delete fHistList;
   delete fSecondRunList;
 }
 

//-----------------------------------------------------------------------

void AliFlowAnalysisWithLYZEventPlane::WriteHistograms(TString* outputFileName)
{
 //store the final results in output .root file

  TFile *output = new TFile(outputFileName->Data(),"RECREATE");
  output->WriteObject(fHistList, "cobjLYZEP","SingleKey");
  delete output;
}

//-----------------------------------------------------------------------

void AliFlowAnalysisWithLYZEventPlane::WriteHistograms(TString outputFileName)
{
 //store the final results in output .root file

  TFile *output = new TFile(outputFileName.Data(),"RECREATE");
  output->WriteObject(fHistList, "cobjLYZEP","SingleKey");
  delete output;
}

//-----------------------------------------------------------------------
void AliFlowAnalysisWithLYZEventPlane::Init() {

  //Initialise all histograms
  cout<<"---Analysis with Lee Yang Zeros Event Plane Method---"<<endl;

  //input histograms
  if (fSecondRunList) {
    fSecondReDtheta = (TProfile*)fSecondRunList->FindObject("Second_FlowPro_ReDtheta_LYZ");
    fHistList->Add(fSecondReDtheta);

    fSecondImDtheta = (TProfile*)fSecondRunList->FindObject("Second_FlowPro_ImDtheta_LYZ");
    fHistList->Add(fSecondImDtheta);
    
    fFirstr0theta = (TProfile*)fSecondRunList->FindObject("First_FlowPro_r0theta_LYZ");
    fHistList->Add(fFirstr0theta);

    //warnings
    if (!fSecondReDtheta) {cout<<"fSecondReDtheta is NULL!"<<endl; }
    if (!fSecondImDtheta) {cout<<"fSecondImDtheta is NULL!"<<endl; }
    if (!fFirstr0theta)   {cout<<"fFirstr0theta is NULL!"<<endl; }

  }

  fCommonHists = new AliFlowCommonHist("AliFlowCommonHistLYZEP");
  fHistList->Add(fCommonHists);
  
  fCommonHistsRes = new AliFlowCommonHistResults("AliFlowCommonHistResultsLYZEP"); 
  fHistList->Add(fCommonHistsRes); 
    
  Int_t iNbinsPt = AliFlowCommonConstants::GetNbinsPt();
  Double_t  dPtMin = AliFlowCommonConstants::GetPtMin();	     
  Double_t  dPtMax = AliFlowCommonConstants::GetPtMax();

  fHistProFlow = new TProfile("FlowPro_VPt_LYZEP","FlowPro_VPt_LYZEP",iNbinsPt,dPtMin,dPtMax);
  fHistProFlow->SetXTitle("Pt");
  fHistProFlow->SetYTitle("v2 (%)");
  fHistList->Add(fHistProFlow);
  
  fHistProWr = new TProfile("FlowPro_Wr_LYZEP","FlowPro_Wr_LYZEP",100,0.,0.25);
  fHistProWr->SetXTitle("Q");
  fHistProWr->SetYTitle("Wr");
  fHistList->Add(fHistProWr);

  fHistQsumforChi = new TH1F("Flow_QsumforChi_LYZEP","Flow_QsumforChi_LYZEP",3,-1.,2.);
  fHistQsumforChi->SetXTitle("Qsum.X , Qsum.Y, Q2sum");
  fHistQsumforChi->SetYTitle("value");
  fHistList->Add(fHistQsumforChi);

  fHistDeltaPhi = new TH1F("Flow_DeltaPhi_LYZEP","Flow_DeltaPhi_LYZEP",100,0.,3.14);
  fHistDeltaPhi->SetXTitle("DeltaPhi");
  fHistDeltaPhi->SetYTitle("Counts");
  fHistList->Add(fHistDeltaPhi);

  fHistPhiLYZ = new TH1F("Flow_PhiLYZ_LYZEP","Flow_PhiLYZ_LYZEP",100,0.,3.14);
  fHistPhiLYZ->SetXTitle("Phi from LYZ");
  fHistPhiLYZ->SetYTitle("Counts");
  fHistList->Add(fHistPhiLYZ);

  fHistPhiEP = new TH1F("Flow_PhiEP_LYZEP","Flow_PhiEP_LYZEP",100,0.,3.14);
  fHistPhiEP->SetXTitle("Phi from EP");
  fHistPhiEP->SetYTitle("Counts");
  fHistList->Add(fHistPhiEP);

  fEventNumber = 0;  //set number of events to zero
      
} 
 
//-----------------------------------------------------------------------
 
void AliFlowAnalysisWithLYZEventPlane::Make(AliFlowEventSimple* anEvent, AliFlowLYZEventPlane* aLYZEP) {
  
  //Get the event plane and weight for each event
  if (anEvent) {
         
    //fill control histograms     
    fCommonHists->FillControlHistograms(anEvent);

    //get the Q vector from the FlowEvent
    AliFlowVector vQ = anEvent->GetQ(); 
    if (vQ.X()== 0. && vQ.Y()== 0. ) { cout<<"Q vector is NULL!"<<endl; }
    //Weight with the multiplicity
    Double_t dQX = 0.;
    Double_t dQY = 0.;
    if (vQ.GetMult()!=0.) {
      dQX = vQ.X()/vQ.GetMult();
      dQY = vQ.Y()/vQ.GetMult();
    } else {cerr<<"vQ.GetMult() is zero!"<<endl; }
    vQ.Set(dQX,dQY);
    //cout<<"vQ("<<dQX<<","<<dQY<<")"<<endl;

    //for chi calculation:
    *fQsum += vQ;
    fHistQsumforChi->SetBinContent(1,fQsum->X());
    fHistQsumforChi->SetBinContent(2,fQsum->Y());
    fQ2sum += vQ.Mod2();
    fHistQsumforChi->SetBinContent(3,fQ2sum);
    //cout<<"fQ2sum = "<<fQ2sum<<endl;

    //call AliFlowLYZEventPlane::CalculateRPandW() here!
    aLYZEP->CalculateRPandW(vQ);

    Double_t dWR = aLYZEP->GetWR();     
    Double_t dRP = aLYZEP->GetPsi();

    //fHistProWr->Fill(vQ.Mod(),dWR); //this weight is always positive
    fHistPhiLYZ->Fill(dRP);   
    
    //plot difference between event plane from EP-method and LYZ-method
    Double_t dRPEP = vQ.Phi()/2;                              //gives distribution from (0 to pi)
    //Double_t dRPEP = 0.5*TMath::ATan2(vQ.Y(),vQ.X());       //gives distribution from (-pi/2 to pi/2)
    //cout<<"dRPEP = "<< dRPEP <<endl;
    fHistPhiEP->Fill(dRPEP);

    Double_t dDeltaPhi = dRPEP - dRP;
    if (dDeltaPhi < 0.) { dDeltaPhi += TMath::Pi(); }        //to shift distribution from (-pi/2 to pi/2) to (0 to pi)
    //cout<<"dDeltaPhi = "<<dDeltaPhi<<endl;
    fHistDeltaPhi->Fill(dDeltaPhi); 

    //Flip sign of WR
    Double_t dLow = TMath::Pi()/4.;
    Double_t dHigh = 3.*(TMath::Pi()/4.);
    if ((dDeltaPhi > dLow) && (dDeltaPhi < dHigh)){
      dRP -= (TMath::Pi()/2);
      dWR = -dWR;
      cerr<<"*** dRP modified ***"<<endl;
    }
    fHistProWr->Fill(vQ.Mod(),dWR); //corrected weight
       
    //calculate flow
    //loop over the tracks of the event
    Int_t iNumberOfTracks = anEvent->NumberOfTracks(); 
    for (Int_t i=0;i<iNumberOfTracks;i++) 
      {
	AliFlowTrackSimple* pTrack = anEvent->GetTrack(i) ; 
	if (pTrack){
	  if (pTrack->UseForDifferentialFlow()) {
	    Double_t dPhi = pTrack->Phi();
	    //if (dPhi<0.) fPhi+=2*TMath::Pi();
	    //calculate flow v2:
	    Double_t dv2 = dWR * TMath::Cos(2*(dPhi-dRP));
	    Double_t dPt = pTrack->Pt();
	    //fill histogram
	    fHistProFlow->Fill(dPt,100*dv2);  
	  }  
	}//track selected
      }//loop over tracks
	  
    fEventNumber++;
    cout<<"@@@@@ "<<fEventNumber<<" events processed"<<endl;
  }
}

  //--------------------------------------------------------------------    
void AliFlowAnalysisWithLYZEventPlane::Finish() {
   
  //plot histograms etc. 
  cout<<"AliFlowAnalysisWithLYZEventPlane::Finish()"<<endl;
  
  //constants:
  Double_t  dJ01 = 2.405; 
  Int_t iNtheta = AliFlowLYZConstants::kTheta;
  Int_t iNbinsPt = AliFlowCommonConstants::GetNbinsPt();
  //set the event number
  SetEventNumber((int)fCommonHists->GetHistMultOrig()->GetEntries());
  //cout<<"number of events processed is "<<fEventNumber<<endl;
  
  //set the sum of Q vectors
  fQsum->Set(fHistQsumforChi->GetBinContent(1),fHistQsumforChi->GetBinContent(2));
  SetQ2sum(fHistQsumforChi->GetBinContent(3));  

  //calculate dV the mean of dVtheta
  Double_t  dVtheta = 0; 
  Double_t  dV = 0; 
  for (Int_t theta=0;theta<iNtheta;theta++)	{
    Double_t dR0 = fFirstr0theta->GetBinContent(theta+1); 
    if (dR0!=0.) { dVtheta = dJ01/dR0 ;}
    dV += dVtheta;
  }
  dV /= iNtheta;

  //calculate dChi 
  Double_t  dSigma2 = 0;
  Double_t  dChi= 0;
  if (fEventNumber!=0) {
    *fQsum /= fEventNumber;
    //cerr<<"fQsum->X() = "<<fQsum->X()<<endl;
    //cerr<<"fQsum->Y() = "<<fQsum->Y()<<endl;
    fQ2sum /= fEventNumber;
    //cerr<<"fEventNumber = "<<fEventNumber<<endl;
    //cerr<<"fQ2sum = "<<fQ2sum<<endl;
    dSigma2 = fQ2sum - TMath::Power(fQsum->X(),2.) - TMath::Power(fQsum->Y(),2.) - TMath::Power(dV,2.);  //BP eq. 62
    //cerr<<"dSigma2"<<dSigma2<<endl;
    if (dSigma2>0) dChi = dV/TMath::Sqrt(dSigma2);
    else dChi = -1.;
    fCommonHistsRes->FillChi(dChi);
    cerr<<"dV = "<<dV<<" and chi = "<<dChi<<endl;
  }
  
  for(Int_t b=0;b<iNbinsPt;b++){
    Double_t dv2pro = 0.;
    Double_t dErr2difcomb = 0.;   
    Double_t dErrdifcomb = 0.;
    if(fHistProFlow) {
      dv2pro = fHistProFlow->GetBinContent(b);
      //calculate error
      for (Int_t theta=0;theta<iNtheta;theta++) {
	Double_t dTheta = ((double)theta/iNtheta)*TMath::Pi(); 
	Int_t iNprime = TMath::Nint(fHistProFlow->GetBinEntries(b));
	//cerr<<"iNprime = "<<iNprime<<endl;
	if (iNprime!=0) { 
	  Double_t dApluscomb = TMath::Exp((dJ01*dJ01)/(2*dChi*dChi)*
					   TMath::Cos(dTheta));
	  Double_t dAmincomb = TMath::Exp(-(dJ01*dJ01)/(2*dChi*dChi)*
					  TMath::Cos(dTheta));
	  dErr2difcomb += (TMath::Cos(dTheta)/(4*iNprime*TMath::BesselJ1(dJ01)*
						 TMath::BesselJ1(dJ01)))*
	    ((dApluscomb*TMath::BesselJ0(2*dJ01*TMath::Sin(dTheta/2))) - 
	     (dAmincomb*TMath::BesselJ0(2*dJ01*TMath::Cos(dTheta/2))));
	} //if !=0
	//else { cout<<"iNprime = 0"<<endl; }
      } //loop over theta
      
      if (dErr2difcomb!=0.) {
	dErr2difcomb /= iNtheta;
	dErrdifcomb = TMath::Sqrt(dErr2difcomb)*100;
	//cerr<<"dErrdifcomb = "<<dErrdifcomb<<endl;
      }
      else {dErrdifcomb = 0.; }

      //fill TH1D
      fCommonHistsRes->FillDifferentialFlow(b, dv2pro, dErrdifcomb); 
    
    } //if fHistProFLow
    else  {
      cout << "Profile Hist missing" << endl;
      break;
    }
    
  } //loop over b

  cout<<".....finished"<<endl;
 }
Commit	Line	Data
f1d945a1	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	$Log$
	18	*/
	19
	20	//#define AliFlowAnalysisWithLYZEventPlane_cxx
	21
	22	#include "Riostream.h" //needed as include
	23
	24	#include "TFile.h"
28ca24ad	25	#include "TList.h"
f1d945a1	26	#include "TComplex.h" //needed as include
	27	#include "TCanvas.h" //needed as include
	28	#include "TLegend.h" //needed as include
	29	#include "TProfile.h" //needed as include
9d062fe3	30	#include "TVector2.h"
f1d945a1	31
	32	class TH1F;
	33
	34	#include "AliFlowLYZConstants.h" //needed as include
	35	#include "AliFlowCommonConstants.h" //needed as include
	36	#include "AliFlowEventSimple.h"
	37	#include "AliFlowTrackSimple.h"
	38	#include "AliFlowCommonHist.h"
	39	#include "AliFlowCommonHistResults.h"
	40	#include "AliFlowLYZEventPlane.h"
	41	#include "AliFlowAnalysisWithLYZEventPlane.h"
	42
	43	class AliFlowVector;
	44
	45	// AliFlowAnalysisWithLYZEventPlane:
	46	//
	47	// Class to do flow analysis with the event plane from the LYZ method
	48	//
	49	// author: N. van der Kolk (kolk@nikhef.nl)
	50
	51
	52	ClassImp(AliFlowAnalysisWithLYZEventPlane)
	53
	54	//-----------------------------------------------------------------------
	55
	56	AliFlowAnalysisWithLYZEventPlane::AliFlowAnalysisWithLYZEventPlane():
9d062fe3	57	fHistList(NULL),
	58	fSecondRunList(NULL),
	59	fSecondReDtheta(NULL),
	60	fSecondImDtheta(NULL),
	61	fFirstr0theta(NULL),
	62	fHistProFlow(NULL),
	63	fHistProFlow2(NULL),
	64	fHistProWr(NULL),
	65	fHistProWrCorr(NULL),
	66	fHistQsumforChi(NULL),
	67	fHistDeltaPhi(NULL),
	68	fHistDeltaPhi2(NULL),
	69	fHistDeltaPhihere(NULL),
	70	fHistPhiEP(NULL),
	71	fHistPhiEPhere(NULL),
	72	fHistPhiLYZ(NULL),
	73	fHistPhiLYZ2(NULL),
	74	fCommonHists(NULL),
	75	fCommonHistsRes(NULL),
	76	fEventNumber(0),
	77	fQsum(NULL),
	78	fQ2sum(0)
f1d945a1	79	{
	80
	81	// Constructor.
0092f3c2	82	fQsum = new TVector2(); // flow vector sum
28ca24ad	83
28ca24ad	84	fHistList = new TList();
9d062fe3	85	fSecondRunList = new TList();
f1d945a1	86	}
	87
	88
	89
	90	//-----------------------------------------------------------------------
	91
	92
	93	AliFlowAnalysisWithLYZEventPlane::~AliFlowAnalysisWithLYZEventPlane()
	94	{
	95	//destructor
0092f3c2	96	delete fQsum;
28ca24ad	97	delete fHistList;
9d062fe3	98	delete fSecondRunList;
f1d945a1	99	}
	100
	101
1dfa3c16	102	//-----------------------------------------------------------------------
	103
	104	void AliFlowAnalysisWithLYZEventPlane::WriteHistograms(TString* outputFileName)
	105	{
	106	//store the final results in output .root file
	107
	108	TFile *output = new TFile(outputFileName->Data(),"RECREATE");
	109	output->WriteObject(fHistList, "cobjLYZEP","SingleKey");
	110	delete output;
	111	}
	112
b0fda271	113	//-----------------------------------------------------------------------
	114
	115	void AliFlowAnalysisWithLYZEventPlane::WriteHistograms(TString outputFileName)
	116	{
	117	//store the final results in output .root file
	118
	119	TFile *output = new TFile(outputFileName.Data(),"RECREATE");
	120	output->WriteObject(fHistList, "cobjLYZEP","SingleKey");
	121	delete output;
	122	}
	123
f1d945a1	124	//-----------------------------------------------------------------------
	125	void AliFlowAnalysisWithLYZEventPlane::Init() {
	126
	127	//Initialise all histograms
	128	cout<<"---Analysis with Lee Yang Zeros Event Plane Method---"<<endl;
	129
	130	//input histograms
9d062fe3	131	if (fSecondRunList) {
9d062fe3	132	fSecondReDtheta = (TProfile*)fSecondRunList->FindObject("Second_FlowPro_ReDtheta_LYZ");
28ca24ad	133	fHistList->Add(fSecondReDtheta);
9d062fe3	134
9d062fe3	135	fSecondImDtheta = (TProfile*)fSecondRunList->FindObject("Second_FlowPro_ImDtheta_LYZ");
28ca24ad	136	fHistList->Add(fSecondImDtheta);
9d062fe3	137
9d062fe3	138	fFirstr0theta = (TProfile*)fSecondRunList->FindObject("First_FlowPro_r0theta_LYZ");
28ca24ad	139	fHistList->Add(fFirstr0theta);
f1d945a1	140
9d062fe3	141	//warnings
	142	if (!fSecondReDtheta) {cout<<"fSecondReDtheta is NULL!"<<endl; }
	143	if (!fSecondImDtheta) {cout<<"fSecondImDtheta is NULL!"<<endl; }
	144	if (!fFirstr0theta) {cout<<"fFirstr0theta is NULL!"<<endl; }
f1d945a1	145
9d062fe3	146	}
28ca24ad	147
9d062fe3	148	fCommonHists = new AliFlowCommonHist("AliFlowCommonHistLYZEP");
9d062fe3	149	fHistList->Add(fCommonHists);
f1d945a1	150
9d062fe3	151	fCommonHistsRes = new AliFlowCommonHistResults("AliFlowCommonHistResultsLYZEP");
	152	fHistList->Add(fCommonHistsRes);
	153
882ffd6a	154	Int_t iNbinsPt = AliFlowCommonConstants::GetNbinsPt();
	155	Double_t dPtMin = AliFlowCommonConstants::GetPtMin();
	156	Double_t dPtMax = AliFlowCommonConstants::GetPtMax();
f1d945a1	157
882ffd6a	158	fHistProFlow = new TProfile("FlowPro_VPt_LYZEP","FlowPro_VPt_LYZEP",iNbinsPt,dPtMin,dPtMax);
f1d945a1	159	fHistProFlow->SetXTitle("Pt");
f1d945a1	160	fHistProFlow->SetYTitle("v2 (%)");
28ca24ad	161	fHistList->Add(fHistProFlow);
f1d945a1	162
	163	fHistProWr = new TProfile("FlowPro_Wr_LYZEP","FlowPro_Wr_LYZEP",100,0.,0.25);
	164	fHistProWr->SetXTitle("Q");
	165	fHistProWr->SetYTitle("Wr");
28ca24ad	166	fHistList->Add(fHistProWr);
f1d945a1	167
9d062fe3	168	fHistQsumforChi = new TH1F("Flow_QsumforChi_LYZEP","Flow_QsumforChi_LYZEP",3,-1.,2.);
	169	fHistQsumforChi->SetXTitle("Qsum.X , Qsum.Y, Q2sum");
	170	fHistQsumforChi->SetYTitle("value");
	171	fHistList->Add(fHistQsumforChi);
	172
f1d945a1	173	fHistDeltaPhi = new TH1F("Flow_DeltaPhi_LYZEP","Flow_DeltaPhi_LYZEP",100,0.,3.14);
	174	fHistDeltaPhi->SetXTitle("DeltaPhi");
	175	fHistDeltaPhi->SetYTitle("Counts");
28ca24ad	176	fHistList->Add(fHistDeltaPhi);
f1d945a1	177
	178	fHistPhiLYZ = new TH1F("Flow_PhiLYZ_LYZEP","Flow_PhiLYZ_LYZEP",100,0.,3.14);
	179	fHistPhiLYZ->SetXTitle("Phi from LYZ");
	180	fHistPhiLYZ->SetYTitle("Counts");
28ca24ad	181	fHistList->Add(fHistPhiLYZ);
f1d945a1	182
	183	fHistPhiEP = new TH1F("Flow_PhiEP_LYZEP","Flow_PhiEP_LYZEP",100,0.,3.14);
	184	fHistPhiEP->SetXTitle("Phi from EP");
	185	fHistPhiEP->SetYTitle("Counts");
28ca24ad	186	fHistList->Add(fHistPhiEP);
f1d945a1	187
f1d945a1	188	fEventNumber = 0; //set number of events to zero
f1d945a1	189
	190	}
	191
	192	//-----------------------------------------------------------------------
	193
882ffd6a	194	void AliFlowAnalysisWithLYZEventPlane::Make(AliFlowEventSimple* anEvent, AliFlowLYZEventPlane* aLYZEP) {
882ffd6a	195
f1d945a1	196	//Get the event plane and weight for each event
8232a5ec	197	if (anEvent) {
f1d945a1	198
f1d945a1	199	//fill control histograms
8232a5ec	200	fCommonHists->FillControlHistograms(anEvent);
f1d945a1	201
f1d945a1	202	//get the Q vector from the FlowEvent
882ffd6a	203	AliFlowVector vQ = anEvent->GetQ();
9d062fe3	204	if (vQ.X()== 0. && vQ.Y()== 0. ) { cout<<"Q vector is NULL!"<<endl; }
f1d945a1	205	//Weight with the multiplicity
882ffd6a	206	Double_t dQX = 0.;
	207	Double_t dQY = 0.;
	208	if (vQ.GetMult()!=0.) {
	209	dQX = vQ.X()/vQ.GetMult();
	210	dQY = vQ.Y()/vQ.GetMult();
	211	} else {cerr<<"vQ.GetMult() is zero!"<<endl; }
	212	vQ.Set(dQX,dQY);
9d062fe3	213	//cout<<"vQ("<<dQX<<","<<dQY<<")"<<endl;
28ca24ad	214
f1d945a1	215	//for chi calculation:
882ffd6a	216	*fQsum += vQ;
9d062fe3	217	fHistQsumforChi->SetBinContent(1,fQsum->X());
9d062fe3	218	fHistQsumforChi->SetBinContent(2,fQsum->Y());
882ffd6a	219	fQ2sum += vQ.Mod2();
9d062fe3	220	fHistQsumforChi->SetBinContent(3,fQ2sum);
9d062fe3	221	//cout<<"fQ2sum = "<<fQ2sum<<endl;
f1d945a1	222
f1d945a1	223	//call AliFlowLYZEventPlane::CalculateRPandW() here!
882ffd6a	224	aLYZEP->CalculateRPandW(vQ);
f1d945a1	225
882ffd6a	226	Double_t dWR = aLYZEP->GetWR();
882ffd6a	227	Double_t dRP = aLYZEP->GetPsi();
f1d945a1	228
882ffd6a	229	//fHistProWr->Fill(vQ.Mod(),dWR); //this weight is always positive
882ffd6a	230	fHistPhiLYZ->Fill(dRP);
f1d945a1	231
f1d945a1	232	//plot difference between event plane from EP-method and LYZ-method
882ffd6a	233	Double_t dRPEP = vQ.Phi()/2; //gives distribution from (0 to pi)
	234	//Double_t dRPEP = 0.5*TMath::ATan2(vQ.Y(),vQ.X()); //gives distribution from (-pi/2 to pi/2)
	235	//cout<<"dRPEP = "<< dRPEP <<endl;
	236	fHistPhiEP->Fill(dRPEP);
f1d945a1	237
882ffd6a	238	Double_t dDeltaPhi = dRPEP - dRP;
	239	if (dDeltaPhi < 0.) { dDeltaPhi += TMath::Pi(); } //to shift distribution from (-pi/2 to pi/2) to (0 to pi)
	240	//cout<<"dDeltaPhi = "<<dDeltaPhi<<endl;
	241	fHistDeltaPhi->Fill(dDeltaPhi);
f1d945a1	242
f1d945a1	243	//Flip sign of WR
882ffd6a	244	Double_t dLow = TMath::Pi()/4.;
	245	Double_t dHigh = 3.*(TMath::Pi()/4.);
	246	if ((dDeltaPhi > dLow) && (dDeltaPhi < dHigh)){
	247	dRP -= (TMath::Pi()/2);
	248	dWR = -dWR;
	249	cerr<<"* dRP modified *"<<endl;
f1d945a1	250	}
882ffd6a	251	fHistProWr->Fill(vQ.Mod(),dWR); //corrected weight
f1d945a1	252
	253	//calculate flow
	254	//loop over the tracks of the event
882ffd6a	255	Int_t iNumberOfTracks = anEvent->NumberOfTracks();
882ffd6a	256	for (Int_t i=0;i<iNumberOfTracks;i++)
f1d945a1	257	{
882ffd6a	258	AliFlowTrackSimple* pTrack = anEvent->GetTrack(i) ;
	259	if (pTrack){
	260	if (pTrack->UseForDifferentialFlow()) {
	261	Double_t dPhi = pTrack->Phi();
	262	//if (dPhi<0.) fPhi+=2*TMath::Pi();
f1d945a1	263	//calculate flow v2:
882ffd6a	264	Double_t dv2 = dWR * TMath::Cos(2*(dPhi-dRP));
882ffd6a	265	Double_t dPt = pTrack->Pt();
f1d945a1	266	//fill histogram
882ffd6a	267	fHistProFlow->Fill(dPt,100*dv2);
f1d945a1	268	}
	269	}//track selected
	270	}//loop over tracks
	271
	272	fEventNumber++;
	273	cout<<"@@@@@ "<<fEventNumber<<" events processed"<<endl;
	274	}
	275	}
	276
	277	//--------------------------------------------------------------------
	278	void AliFlowAnalysisWithLYZEventPlane::Finish() {
	279
	280	//plot histograms etc.
bc92c0cb	281	cout<<"AliFlowAnalysisWithLYZEventPlane::Finish()"<<endl;
9d062fe3	282
9d062fe3	283	//constants:
882ffd6a	284	Double_t dJ01 = 2.405;
	285	Int_t iNtheta = AliFlowLYZConstants::kTheta;
	286	Int_t iNbinsPt = AliFlowCommonConstants::GetNbinsPt();
9d062fe3	287	//set the event number
	288	SetEventNumber((int)fCommonHists->GetHistMultOrig()->GetEntries());
	289	//cout<<"number of events processed is "<<fEventNumber<<endl;
	290
	291	//set the sum of Q vectors
	292	fQsum->Set(fHistQsumforChi->GetBinContent(1),fHistQsumforChi->GetBinContent(2));
	293	SetQ2sum(fHistQsumforChi->GetBinContent(3));
	294
882ffd6a	295	//calculate dV the mean of dVtheta
	296	Double_t dVtheta = 0;
	297	Double_t dV = 0;
	298	for (Int_t theta=0;theta<iNtheta;theta++) {
	299	Double_t dR0 = fFirstr0theta->GetBinContent(theta+1);
	300	if (dR0!=0.) { dVtheta = dJ01/dR0 ;}
	301	dV += dVtheta;
f1d945a1	302	}
882ffd6a	303	dV /= iNtheta;
f1d945a1	304
882ffd6a	305	//calculate dChi
	306	Double_t dSigma2 = 0;
	307	Double_t dChi= 0;
f1d945a1	308	if (fEventNumber!=0) {
0092f3c2	309	*fQsum /= fEventNumber;
bc92c0cb	310	//cerr<<"fQsum->X() = "<<fQsum->X()<<endl;
bc92c0cb	311	//cerr<<"fQsum->Y() = "<<fQsum->Y()<<endl;
f1d945a1	312	fQ2sum /= fEventNumber;
bc92c0cb	313	//cerr<<"fEventNumber = "<<fEventNumber<<endl;
bc92c0cb	314	//cerr<<"fQ2sum = "<<fQ2sum<<endl;
882ffd6a	315	dSigma2 = fQ2sum - TMath::Power(fQsum->X(),2.) - TMath::Power(fQsum->Y(),2.) - TMath::Power(dV,2.); //BP eq. 62
bc92c0cb	316	//cerr<<"dSigma2"<<dSigma2<<endl;
882ffd6a	317	if (dSigma2>0) dChi = dV/TMath::Sqrt(dSigma2);
	318	else dChi = -1.;
	319	fCommonHistsRes->FillChi(dChi);
	320	cerr<<"dV = "<<dV<<" and chi = "<<dChi<<endl;
f1d945a1	321	}
f1d945a1	322
882ffd6a	323	for(Int_t b=0;b<iNbinsPt;b++){
	324	Double_t dv2pro = 0.;
	325	Double_t dErr2difcomb = 0.;
	326	Double_t dErrdifcomb = 0.;
f1d945a1	327	if(fHistProFlow) {
882ffd6a	328	dv2pro = fHistProFlow->GetBinContent(b);
f1d945a1	329	//calculate error
882ffd6a	330	for (Int_t theta=0;theta<iNtheta;theta++) {
	331	Double_t dTheta = ((double)theta/iNtheta)*TMath::Pi();
	332	Int_t iNprime = TMath::Nint(fHistProFlow->GetBinEntries(b));
	333	//cerr<<"iNprime = "<<iNprime<<endl;
	334	if (iNprime!=0) {
	335	Double_t dApluscomb = TMath::Exp((dJ01dJ01)/(2dChidChi)
	336	TMath::Cos(dTheta));
	337	Double_t dAmincomb = TMath::Exp(-(dJ01dJ01)/(2dChidChi)
	338	TMath::Cos(dTheta));
	339	dErr2difcomb += (TMath::Cos(dTheta)/(4iNprimeTMath::BesselJ1(dJ01)*
	340	TMath::BesselJ1(dJ01)))*
	341	((dApluscombTMath::BesselJ0(2dJ01*TMath::Sin(dTheta/2))) -
	342	(dAmincombTMath::BesselJ0(2dJ01*TMath::Cos(dTheta/2))));
f1d945a1	343	} //if !=0
882ffd6a	344	//else { cout<<"iNprime = 0"<<endl; }
f1d945a1	345	} //loop over theta
f1d945a1	346
882ffd6a	347	if (dErr2difcomb!=0.) {
	348	dErr2difcomb /= iNtheta;
	349	dErrdifcomb = TMath::Sqrt(dErr2difcomb)*100;
	350	//cerr<<"dErrdifcomb = "<<dErrdifcomb<<endl;
f1d945a1	351	}
882ffd6a	352	else {dErrdifcomb = 0.; }
f1d945a1	353
f1d945a1	354	//fill TH1D
882ffd6a	355	fCommonHistsRes->FillDifferentialFlow(b, dv2pro, dErrdifcomb);
f1d945a1	356
	357	} //if fHistProFLow
	358	else {
	359	cout << "Profile Hist missing" << endl;
	360	break;
	361	}
	362
	363	} //loop over b
	364
f1d945a1	365	cout<<".....finished"<<endl;
f1d945a1	366	}
882ffd6a	367