[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),
  fQ(NULL),
  fQsum(NULL),
  fQ2sum(0),
  fQtheta(0),
  fTrack(NULL)
//  fLYZEP(NULL)
{

  // Constructor.
  //  fQ.Set(0.,0.);           // flow vector
  //  fQsum.Set(0.,0.);        // flow vector sum
  fQ = new AliFlowVector();           // flow vector
  fQsum = new TVector2();        // flow vector sum
  //  fLYZEP = new AliFlowLYZEventPlane();
}

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


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

//-----------------------------------------------------------------------
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 fNtheta = AliFlowLYZConstants::kTheta;
  Int_t fNbinsPt = AliFlowCommonConstants::GetNbinsPt();
  Double_t  fPtMin = AliFlowCommonConstants::GetPtMin();	     
  Double_t  fPtMax = AliFlowCommonConstants::GetPtMax();


  fHistProFlow = new TProfile("FlowPro_VPt_LYZEP","FlowPro_VPt_LYZEP",fNbinsPt,fPtMin,fPtMax);
  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",fNtheta,-0.5,fNtheta-0.5);
  fHistProR0theta->SetXTitle("#theta");
  fHistProR0theta->SetYTitle("r_{0}^{#theta}");

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

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

  fEventNumber = 0;  //set number of events to zero

      
} 
 
//-----------------------------------------------------------------------
 
void AliFlowAnalysisWithLYZEventPlane::Make(AliFlowEventSimple* anEvent, AliFlowLYZEventPlane* fLYZEP) {
//void AliFlowAnalysisWithLYZEventPlane::Make(AliFlowEventSimple* anEvent) {

  //Get the event plane and weight for each event
  if (anEvent) {
         
    //fill control histograms     
    fCommonHists->FillControlHistograms(anEvent);

    //get the Q vector from the FlowEvent
    *fQ = anEvent->GetQ(); 
    //Weight with the multiplicity
    Double_t fQX = 0.;
    Double_t fQY = 0.;
    if (fQ->GetMult()!=0.) {
      fQX = fQ->X()/fQ->GetMult();
      fQY = fQ->Y()/fQ->GetMult();
    } else {cerr<<"fQ->GetMult() is zero!"<<endl; }
    fQ->Set(fQX,fQY);
    //cout<<"fQ->Mod() = " << fQ->Mod() << endl;
    //for chi calculation:
    *fQsum += *fQ;
    //cout<<"fQsum.Mod() = "<<fQsum.Mod()<<endl;
    fQ2sum += fQ->Mod2();
    cout<<"fQ2sum = "<<fQ2sum<<endl;

    //call AliFlowLYZEventPlane::CalculateRPandW() here!
    fLYZEP->CalculateRPandW(*fQ);

    Double_t fWR = fLYZEP->GetWR();     
    Double_t fRP = fLYZEP->GetPsi();

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

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

    //Flip sign of WR
    Double_t low = TMath::Pi()/4.;
    Double_t high = 3.*(TMath::Pi()/4.);
    if ((fDeltaPhi > low) && (fDeltaPhi < high)){
      fRP -= (TMath::Pi()/2);
      fWR = -fWR;
      cerr<<"*** fRP modified ***"<<endl;
    }
    fHistProWr->Fill(fQ->Mod(),fWR); //corrected weight
       
    //calculate flow
    //loop over the tracks of the event
    Int_t fNumberOfTracks = anEvent->NumberOfTracks(); 
    for (Int_t i=0;i<fNumberOfTracks;i++) 
      {
	fTrack = anEvent->GetTrack(i) ; 
	if (fTrack){
	  if (fTrack->UseForDifferentialFlow()) {
	    Double_t fPhi = fTrack->Phi();
	    //if (fPhi<0.) fPhi+=2*TMath::Pi();
	    //calculate flow v2:
	    Double_t fv2 = fWR * TMath::Cos(2*(fPhi-fRP));
	    Double_t fPt = fTrack->Pt();
	    //fill histogram
	    fHistProFlow->Fill(fPt,100*fv2);  
	  }  
	}//track selected
      }//loop over tracks
	  
    fEventNumber++;
    cout<<"@@@@@ "<<fEventNumber<<" events processed"<<endl;
  }
}

  //--------------------------------------------------------------------    
void AliFlowAnalysisWithLYZEventPlane::Finish() {
   
  //plot histograms etc. 
  cout<<"AliFlowAnalysisWithLYZEventPlane::Terminate()"<<endl;
  //constands:
  Double_t  fJ01 = 2.405; 
  Int_t fNtheta = AliFlowLYZConstants::kTheta;
  Int_t fNbinsPt = AliFlowCommonConstants::GetNbinsPt();
  //Double_t fErr2difcomb = 0.;
  //Double_t fErrdifcomb = 0.;
  
  //calculate fV the mean of fVtheta
  Double_t  fVtheta = 0; 
  Double_t  fV = 0; 
  for (Int_t theta=0;theta<fNtheta;theta++)	{
    Double_t fR0 = fFirstr0theta->GetBinContent(theta+1); 
    if (fR0!=0.) { fVtheta = fJ01/fR0 ;}
    fV += fVtheta;
  }
  fV /= fNtheta;

  //calculate fChi 
  Double_t  fSigma2 = 0;
  Double_t  fChi= 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;
    fSigma2 = fQ2sum - TMath::Power(fQsum->X(),2.) - TMath::Power(fQsum->Y(),2.) - TMath::Power(fV,2.);  //BP eq. 62
    if (fSigma2>0) fChi = fV/TMath::Sqrt(fSigma2);
    else fChi = -1.;
    fCommonHistsRes->FillChi(fChi);
    cerr<<"fV = "<<fV<<" and chi = "<<fChi<<endl;
  }
  
  for(Int_t b=0;b<fNbinsPt;b++){
    Double_t fv2pro = 0.;
    Double_t fErr2difcomb = 0.;   
    Double_t fErrdifcomb = 0.;
    if(fHistProFlow) {
      fv2pro = fHistProFlow->GetBinContent(b);
      //calculate error
      for (Int_t theta=0;theta<fNtheta;theta++) {
	Double_t fTheta = ((double)theta/fNtheta)*TMath::Pi(); 
	Int_t fNprime = TMath::Nint(fHistProFlow->GetBinEntries(b));
	//cerr<<"fNprime = "<<fNprime<<endl;
	if (fNprime!=0.) { 
	  Double_t fApluscomb = TMath::Exp((fJ01*fJ01)/(2*fChi*fChi)*
					   TMath::Cos(fTheta));
	  Double_t fAmincomb = TMath::Exp(-(fJ01*fJ01)/(2*fChi*fChi)*
					  TMath::Cos(fTheta));
	  fErr2difcomb += (TMath::Cos(fTheta)/(4*fNprime*TMath::BesselJ1(fJ01)*
						 TMath::BesselJ1(fJ01)))*
	    ((fApluscomb*TMath::BesselJ0(2*fJ01*TMath::Sin(fTheta/2))) - 
	     (fAmincomb*TMath::BesselJ0(2*fJ01*TMath::Cos(fTheta/2))));
	} //if !=0
	//else { cout<<"fNprime = 0."<<endl; }
      } //loop over theta
      
      if (fErr2difcomb!=0.) {
	fErr2difcomb /= fNtheta;
	fErrdifcomb = TMath::Sqrt(fErr2difcomb)*100;
	//cerr<<"fErrdifcomb = "<<fErrdifcomb<<endl;
      }
      else {fErrdifcomb = 0.; }

      //fill TH1D
      fCommonHistsRes->FillDifferentialFlow(b, fv2pro, fErrdifcomb); 
    
    } //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->GetfHistDiffFlow();
  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),
	83	fQ(NULL),
	84	fQsum(NULL),
	85	fQ2sum(0),
	86	fQtheta(0),
d97a8fc5	87	fTrack(NULL)
d97a8fc5	88	// fLYZEP(NULL)
f1d945a1	89	{
	90
	91	// Constructor.
0092f3c2	92	// fQ.Set(0.,0.); // flow vector
	93	// fQsum.Set(0.,0.); // flow vector sum
	94	fQ = new AliFlowVector(); // flow vector
	95	fQsum = new TVector2(); // flow vector sum
d97a8fc5	96	// fLYZEP = new AliFlowLYZEventPlane();
f1d945a1	97	}
	98
	99
	100
	101	//-----------------------------------------------------------------------
	102
	103
	104	AliFlowAnalysisWithLYZEventPlane::~AliFlowAnalysisWithLYZEventPlane()
	105	{
	106	//destructor
0092f3c2	107	delete fQ;
0092f3c2	108	delete fQsum;
d97a8fc5	109	// delete fLYZEP;
f1d945a1	110	}
	111
	112
	113	//-----------------------------------------------------------------------
	114	void AliFlowAnalysisWithLYZEventPlane::Init() {
	115
	116	//Initialise all histograms
	117	cout<<"---Analysis with Lee Yang Zeros Event Plane Method---"<<endl;
	118
	119	//input histograms
	120	if (fSecondRunFile->IsZombie()){ //check if file exists
	121	cout << "Error opening file, run first with fFirstrun = kTRUE" << endl;
	122	exit(-1);
	123	} else if (fSecondRunFile->IsOpen()){
	124	cout<<"----secondRunFile is open----"<<endl;
	125	fSecondVPt = (TProfile*)fSecondRunFile->Get("Flow_Differential_Pt_LYZ"); //to compare to
	126	}
	127
	128	if (fFirstRunFile->IsZombie()){ //check if file exists
	129	cout << "Error opening file, run first with fFirstrun = kTRUE" << endl;
	130	exit(-1);
	131	} else if (fFirstRunFile->IsOpen()){
	132	cout<<"----firstRunFile is open----"<<endl<<endl;
	133	fFirstr0theta = (TProfile*)fFirstRunFile->Get("First_FlowPro_r0theta_LYZ");
	134	}
	135
	136
	137	//output file
	138	// ********make output file
	139	// analysis file (output)
	140	fOutFile = new TFile(fOutFileName.Data(),"RECREATE") ;
	141
	142	fCommonHists = new AliFlowCommonHist("LYZEP");
	143	fCommonHistsRes = new AliFlowCommonHistResults("LYZEP");
	144
	145	// output histograms
	146	Int_t fNtheta = AliFlowLYZConstants::kTheta;
	147	Int_t fNbinsPt = AliFlowCommonConstants::GetNbinsPt();
	148	Double_t fPtMin = AliFlowCommonConstants::GetPtMin();
	149	Double_t fPtMax = AliFlowCommonConstants::GetPtMax();
	150
	151
	152	fHistProFlow = new TProfile("FlowPro_VPt_LYZEP","FlowPro_VPt_LYZEP",fNbinsPt,fPtMin,fPtMax);
	153	fHistProFlow->SetXTitle("Pt");
	154	fHistProFlow->SetYTitle("v2 (%)");
	155
	156	fHistProWr = new TProfile("FlowPro_Wr_LYZEP","FlowPro_Wr_LYZEP",100,0.,0.25);
	157	fHistProWr->SetXTitle("Q");
	158	fHistProWr->SetYTitle("Wr");
	159
	160	fHistDeltaPhi = new TH1F("Flow_DeltaPhi_LYZEP","Flow_DeltaPhi_LYZEP",100,0.,3.14);
	161	fHistDeltaPhi->SetXTitle("DeltaPhi");
	162	fHistDeltaPhi->SetYTitle("Counts");
	163
	164	fHistPhiLYZ = new TH1F("Flow_PhiLYZ_LYZEP","Flow_PhiLYZ_LYZEP",100,0.,3.14);
	165	fHistPhiLYZ->SetXTitle("Phi from LYZ");
	166	fHistPhiLYZ->SetYTitle("Counts");
	167
	168	fHistPhiEP = new TH1F("Flow_PhiEP_LYZEP","Flow_PhiEP_LYZEP",100,0.,3.14);
	169	fHistPhiEP->SetXTitle("Phi from EP");
	170	fHistPhiEP->SetYTitle("Counts");
	171
	172
	173	fHistProR0theta = new TProfile("FlowPro_r0theta_LYZEP","FlowPro_r0theta_LYZEP",fNtheta,-0.5,fNtheta-0.5);
174	fHistProR0theta->SetXTitle("#theta");
175	fHistProR0theta->SetYTitle("r_{0}^{#theta}");
176
177	fHistProReDtheta = new TProfile("FlowPro_ReDtheta_LYZEP","FlowPro_ReDtheta_LYZEP",fNtheta, -0.5, fNtheta-0.5);
178	fHistProReDtheta->SetXTitle("#theta");
179	fHistProReDtheta->SetYTitle("Re(D^{#theta})");
180
181	fHistProImDtheta = new TProfile("FlowPro_ImDtheta_LYZEP","FlowPro_ImDtheta_LYZEP",fNtheta, -0.5, fNtheta-0.5);
182	fHistProImDtheta->SetXTitle("#theta");
183	fHistProImDtheta->SetYTitle("Im(D^{#theta})");
184
185	fEventNumber = 0; //set number of events to zero
186
187
188	}
189
190	//-----------------------------------------------------------------------
191
8232a5ec	192	void AliFlowAnalysisWithLYZEventPlane::Make(AliFlowEventSimple* anEvent, AliFlowLYZEventPlane* fLYZEP) {
d97a8fc5	193	//void AliFlowAnalysisWithLYZEventPlane::Make(AliFlowEventSimple* anEvent) {
f1d945a1	194
f1d945a1	195	//Get the event plane and weight for each event
8232a5ec	196	if (anEvent) {
f1d945a1	197
f1d945a1	198	//fill control histograms
8232a5ec	199	fCommonHists->FillControlHistograms(anEvent);
f1d945a1	200
f1d945a1	201	//get the Q vector from the FlowEvent
0092f3c2	202	*fQ = anEvent->GetQ();
f1d945a1	203	//Weight with the multiplicity
7d1e46b8	204	Double_t fQX = 0.;
7d1e46b8	205	Double_t fQY = 0.;
0092f3c2	206	if (fQ->GetMult()!=0.) {
	207	fQX = fQ->X()/fQ->GetMult();
	208	fQY = fQ->Y()/fQ->GetMult();
	209	} else {cerr<<"fQ->GetMult() is zero!"<<endl; }
	210	fQ->Set(fQX,fQY);
	211	//cout<<"fQ->Mod() = " << fQ->Mod() << endl;
f1d945a1	212	//for chi calculation:
0092f3c2	213	fQsum += fQ;
f1d945a1	214	//cout<<"fQsum.Mod() = "<<fQsum.Mod()<<endl;
0092f3c2	215	fQ2sum += fQ->Mod2();
f1d945a1	216	cout<<"fQ2sum = "<<fQ2sum<<endl;
	217
	218	//call AliFlowLYZEventPlane::CalculateRPandW() here!
0092f3c2	219	fLYZEP->CalculateRPandW(*fQ);
f1d945a1	220
	221	Double_t fWR = fLYZEP->GetWR();
	222	Double_t fRP = fLYZEP->GetPsi();
	223
	224	//fHistProWr->Fill(fQ.Mod(),fWR); //this weight is always positive
	225	fHistPhiLYZ->Fill(fRP);
	226
	227	//plot difference between event plane from EP-method and LYZ-method
0092f3c2	228	Double_t fRPEP = fQ->Phi()/2; //gives distribution from (0 to pi)
f1d945a1	229	//Float_t fRPEP = 0.5*TMath::ATan2(fQ.Y(),fQ.X()); //gives distribution from (-pi/2 to pi/2)
	230	//cout<<"fRPEP = "<< fRPEP <<endl;
	231	fHistPhiEP->Fill(fRPEP);
	232
	233	Double_t fDeltaPhi = fRPEP - fRP;
	234	if (fDeltaPhi < 0.) { fDeltaPhi += TMath::Pi(); } //to shift distribution from (-pi/2 to pi/2) to (0 to pi)
	235	//cout<<"fDeltaPhi = "<<fDeltaPhi<<endl;
	236	fHistDeltaPhi->Fill(fDeltaPhi);
	237
	238	//Flip sign of WR
	239	Double_t low = TMath::Pi()/4.;
	240	Double_t high = 3.*(TMath::Pi()/4.);
	241	if ((fDeltaPhi > low) && (fDeltaPhi < high)){
	242	fRP -= (TMath::Pi()/2);
	243	fWR = -fWR;
	244	cerr<<"* fRP modified *"<<endl;
	245	}
0092f3c2	246	fHistProWr->Fill(fQ->Mod(),fWR); //corrected weight
f1d945a1	247
	248	//calculate flow
	249	//loop over the tracks of the event
8232a5ec	250	Int_t fNumberOfTracks = anEvent->NumberOfTracks();
f1d945a1	251	for (Int_t i=0;i<fNumberOfTracks;i++)
f1d945a1	252	{
8232a5ec	253	fTrack = anEvent->GetTrack(i) ;
f1d945a1	254	if (fTrack){
	255	if (fTrack->UseForDifferentialFlow()) {
	256	Double_t fPhi = fTrack->Phi();
	257	//if (fPhi<0.) fPhi+=2*TMath::Pi();
	258	//calculate flow v2:
	259	Double_t fv2 = fWR * TMath::Cos(2*(fPhi-fRP));
	260	Double_t fPt = fTrack->Pt();
	261	//fill histogram
	262	fHistProFlow->Fill(fPt,100*fv2);
	263	}
	264	}//track selected
	265	}//loop over tracks
	266
	267	fEventNumber++;
	268	cout<<"@@@@@ "<<fEventNumber<<" events processed"<<endl;
	269	}
	270	}
	271
	272	//--------------------------------------------------------------------
	273	void AliFlowAnalysisWithLYZEventPlane::Finish() {
	274
	275	//plot histograms etc.
	276	cout<<"AliFlowAnalysisWithLYZEventPlane::Terminate()"<<endl;
	277	//constands:
	278	Double_t fJ01 = 2.405;
	279	Int_t fNtheta = AliFlowLYZConstants::kTheta;
	280	Int_t fNbinsPt = AliFlowCommonConstants::GetNbinsPt();
	281	//Double_t fErr2difcomb = 0.;
	282	//Double_t fErrdifcomb = 0.;
	283
	284	//calculate fV the mean of fVtheta
	285	Double_t fVtheta = 0;
	286	Double_t fV = 0;
	287	for (Int_t theta=0;theta<fNtheta;theta++) {
	288	Double_t fR0 = fFirstr0theta->GetBinContent(theta+1);
	289	if (fR0!=0.) { fVtheta = fJ01/fR0 ;}
	290	fV += fVtheta;
	291	}
	292	fV /= fNtheta;
	293
	294	//calculate fChi
	295	Double_t fSigma2 = 0;
	296	Double_t fChi= 0;
	297	if (fEventNumber!=0) {
0092f3c2	298	*fQsum /= fEventNumber;
f1d945a1	299	//cerr<<"fQsum.X() = "<<fQsum.X()<<endl;
	300	//cerr<<"fQsum.Y() = "<<fQsum.Y()<<endl;
	301	fQ2sum /= fEventNumber;
	302	cerr<<"fEventNumber = "<<fEventNumber<<endl;
	303	cerr<<"fQ2sum = "<<fQ2sum<<endl;
0092f3c2	304	fSigma2 = fQ2sum - TMath::Power(fQsum->X(),2.) - TMath::Power(fQsum->Y(),2.) - TMath::Power(fV,2.); //BP eq. 62
f1d945a1	305	if (fSigma2>0) fChi = fV/TMath::Sqrt(fSigma2);
	306	else fChi = -1.;
	307	fCommonHistsRes->FillChi(fChi);
	308	cerr<<"fV = "<<fV<<" and chi = "<<fChi<<endl;
	309	}
	310
	311	for(Int_t b=0;b<fNbinsPt;b++){
	312	Double_t fv2pro = 0.;
	313	Double_t fErr2difcomb = 0.;
	314	Double_t fErrdifcomb = 0.;
	315	if(fHistProFlow) {
	316	fv2pro = fHistProFlow->GetBinContent(b);
	317	//calculate error
	318	for (Int_t theta=0;theta<fNtheta;theta++) {
	319	Double_t fTheta = ((double)theta/fNtheta)*TMath::Pi();
	320	Int_t fNprime = TMath::Nint(fHistProFlow->GetBinEntries(b));
	321	//cerr<<"fNprime = "<<fNprime<<endl;
	322	if (fNprime!=0.) {
	323	Double_t fApluscomb = TMath::Exp((fJ01fJ01)/(2fChifChi)
	324	TMath::Cos(fTheta));
	325	Double_t fAmincomb = TMath::Exp(-(fJ01fJ01)/(2fChifChi)
	326	TMath::Cos(fTheta));
	327	fErr2difcomb += (TMath::Cos(fTheta)/(4fNprimeTMath::BesselJ1(fJ01)*
	328	TMath::BesselJ1(fJ01)))*
	329	((fApluscombTMath::BesselJ0(2fJ01*TMath::Sin(fTheta/2))) -
	330	(fAmincombTMath::BesselJ0(2fJ01*TMath::Cos(fTheta/2))));
	331	} //if !=0
	332	//else { cout<<"fNprime = 0."<<endl; }
	333	} //loop over theta
	334
	335	if (fErr2difcomb!=0.) {
	336	fErr2difcomb /= fNtheta;
	337	fErrdifcomb = TMath::Sqrt(fErr2difcomb)*100;
	338	//cerr<<"fErrdifcomb = "<<fErrdifcomb<<endl;
	339	}
	340	else {fErrdifcomb = 0.; }
	341
	342	//fill TH1D
	343	fCommonHistsRes->FillDifferentialFlow(b, fv2pro, fErrdifcomb);
	344
	345	} //if fHistProFLow
	346	else {
	347	cout << "Profile Hist missing" << endl;
	348	break;
	349	}
	350
	351	} //loop over b
	352
	353	// write to file
	354	fOutFile->Write();
	355
	356	cout<<"Making some plots to check the results:"<<endl<<endl;
	357
	358	TCanvas *canvas = new TCanvas("canvas","compare v2 vs pt",800,800);
	359	canvas->cd();
	360	fSecondVPt->SetLineColor(3);
	361	fSecondVPt->SetLineWidth(2);
	362	fSecondVPt->Draw();
	363	fHistFlow = fCommonHistsRes->GetfHistDiffFlow();
	364	fHistFlow->Draw("SAME");
	365	// draw the legend
	366	TLegend *legend2 = new TLegend(0.6,0.65,0.88,0.85);
	367	legend2->SetTextFont(72);
	368	legend2->SetTextSize(0.04);
369	legend2->AddEntry(fSecondVPt,"stand. LYZ","lpe");
370	legend2->AddEntry(fHistFlow,"new LYZ with calculated errors","lpe");
371	legend2->Draw();
372
373	TCanvas *canvas5 = new TCanvas("canvas5","phi and Delta Phi",800,600);
374	canvas5->Divide(3,1);
375	canvas5->cd(1);
376	fHistDeltaPhi->Draw();
377	canvas5->cd(2);
378	fHistPhiLYZ->Draw();
379	canvas5->cd(3);
380	fHistPhiEP->Draw();
381
382	cout<<".....finished"<<endl;
383	}