[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 "TComplex.h"   //needed as include
#include "TCanvas.h"   //needed as include
#include "TLegend.h"   //needed as include
#include "TProfile.h"  //needed as include

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():
  fOutFile(0), 
  fFirstRunFile(0),
  fSecondRunFile(0),
  fFirstRunFileName(0),
  fSecondRunFileName(0),
  fOutFileName(0),
  fSecondReDtheta(0),
  fSecondImDtheta(0),
  fFirstr0theta(0),
  fSecondVPt(0),
  fHistProFlow(0),
  fHistProFlow2(0),
  fHistProWr(0),
  fHistProWrCorr(0),
  fHistFlow(0),
  fHistDeltaPhi(0),
  fHistDeltaPhi2(0),
  fHistDeltaPhihere(0),
  fHistPhiEP(0),
  fHistPhiEPhere(0),
  fHistPhiLYZ(0),
  fHistPhiLYZ2(0),
  fHistProR0theta(0),
  fHistProReDtheta(0),
  fHistProImDtheta(0),
  fCommonHists(0),
  fCommonHistsRes(0),
  fEventNumber(0),
  fQsum(NULL),
  fQ2sum(0)
{

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

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


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

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

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

  //input histograms
  if (fSecondRunFile->IsZombie()){ //check if file exists
    cout << "Error opening file, run first with fFirstrun = kTRUE" << endl;
    exit(-1);
  } else if (fSecondRunFile->IsOpen()){
    cout<<"----secondRunFile is open----"<<endl;
    fSecondVPt = (TProfile*)fSecondRunFile->Get("Flow_Differential_Pt_LYZ"); //to compare to
  }

  if (fFirstRunFile->IsZombie()){ //check if file exists
    cout << "Error opening file, run first with fFirstrun = kTRUE" << endl;
    exit(-1);
  } else if (fFirstRunFile->IsOpen()){
    cout<<"----firstRunFile is open----"<<endl<<endl;
    fFirstr0theta = (TProfile*)fFirstRunFile->Get("First_FlowPro_r0theta_LYZ");
  }


  //output file
  // ********make output file 
  // analysis file (output)
  fOutFile = new TFile(fOutFileName.Data(),"RECREATE") ;

  fCommonHists = new AliFlowCommonHist("LYZEP");
  fCommonHistsRes = new AliFlowCommonHistResults("LYZEP");
  
  // output histograms
  Int_t iNtheta = AliFlowLYZConstants::kTheta;
  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 (%)");
  
  fHistProWr = new TProfile("FlowPro_Wr_LYZEP","FlowPro_Wr_LYZEP",100,0.,0.25);
  fHistProWr->SetXTitle("Q");
  fHistProWr->SetYTitle("Wr");

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

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

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

  
  fHistProR0theta  = new TProfile("FlowPro_r0theta_LYZEP","FlowPro_r0theta_LYZEP",iNtheta,-0.5,iNtheta-0.5);
  fHistProR0theta->SetXTitle("#theta");
  fHistProR0theta->SetYTitle("r_{0}^{#theta}");

  fHistProReDtheta = new TProfile("FlowPro_ReDtheta_LYZEP","FlowPro_ReDtheta_LYZEP",iNtheta, -0.5, iNtheta-0.5);
  fHistProReDtheta->SetXTitle("#theta");
  fHistProReDtheta->SetYTitle("Re(D^{#theta})");

  fHistProImDtheta = new TProfile("FlowPro_ImDtheta_LYZEP","FlowPro_ImDtheta_LYZEP",iNtheta, -0.5, iNtheta-0.5);
  fHistProImDtheta->SetXTitle("#theta");
  fHistProImDtheta->SetYTitle("Im(D^{#theta})");

  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(); 
    //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);
    //for chi calculation:
    *fQsum += vQ;
    fQ2sum += vQ.Mod2();
    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::Terminate()"<<endl;
  //constands:
  Double_t  dJ01 = 2.405; 
  Int_t iNtheta = AliFlowLYZConstants::kTheta;
  Int_t iNbinsPt = AliFlowCommonConstants::GetNbinsPt();
    
  //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
    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

  // write to file
  fOutFile->Write();

  cout<<"Making some plots to check the results:"<<endl<<endl;

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