[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

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():
  fSecondRunFile(0),
  fSecondRunFileName(0),
  fHistList(NULL),
  fSecondReDtheta(NULL),
  fSecondImDtheta(NULL),
  fFirstr0theta(NULL),
  fSecondVPt(NULL),
  fHistProFlow(NULL),
  fHistProFlow2(NULL),
  fHistProWr(NULL),
  fHistProWrCorr(NULL),
  fHistFlow(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();
}

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


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

//-----------------------------------------------------------------------
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, no input file from LYZ analysis" << endl;
    exit(-1);
  } else if (fSecondRunFile->IsOpen()){
    cout<<"----secondRunFile is open----"<<endl;
    //get List
    TList* list   = (TList*)fSecondRunFile->Get("cobj1");
    fSecondVPt    = (TProfile*)list->FindObject("Flow_Differential_Pt_LYZ"); //to compare to
    fSecondReDtheta = (TProfile*)list->FindObject("Second_FlowPro_ReDtheta_LYZ");
    fHistList->Add(fSecondReDtheta);
    fSecondImDtheta = (TProfile*)list->FindObject("Second_FlowPro_ImDtheta_LYZ");
    fHistList->Add(fSecondImDtheta);
    fFirstr0theta = (TProfile*)list->FindObject("First_FlowPro_r0theta_LYZ");
    fHistList->Add(fFirstr0theta);
  }


  fCommonHists = new AliFlowCommonHist("LYZEP");
  fHistList->Add(fCommonHists->GetHistMultOrig());
  fHistList->Add(fCommonHists->GetHistMultInt());
  fHistList->Add(fCommonHists->GetHistMultDiff());
  fHistList->Add(fCommonHists->GetHistPtInt());
  fHistList->Add(fCommonHists->GetHistPtDiff());
  fHistList->Add(fCommonHists->GetHistPhiInt());
  fHistList->Add(fCommonHists->GetHistPhiDiff());
  fHistList->Add(fCommonHists->GetHistEtaInt());
  fHistList->Add(fCommonHists->GetHistEtaDiff());
  fHistList->Add(fCommonHists->GetHistProMeanPtperBin());
  fHistList->Add(fCommonHists->GetHistQ());

  fCommonHistsRes = new AliFlowCommonHistResults("LYZEP"); 
  fHistList->Add(fCommonHistsRes->GetHistDiffFlow()); 
  fHistList->Add(fCommonHistsRes->GetHistChi()); 
  fHistList->Add(fCommonHistsRes->GetHistIntFlow()); 

  
  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);

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