[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::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
f1d945a1	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
9d062fe3	120	if (fSecondRunList) {
9d062fe3	121	fSecondReDtheta = (TProfile*)fSecondRunList->FindObject("Second_FlowPro_ReDtheta_LYZ");
28ca24ad	122	fHistList->Add(fSecondReDtheta);
9d062fe3	123
9d062fe3	124	fSecondImDtheta = (TProfile*)fSecondRunList->FindObject("Second_FlowPro_ImDtheta_LYZ");
28ca24ad	125	fHistList->Add(fSecondImDtheta);
9d062fe3	126
9d062fe3	127	fFirstr0theta = (TProfile*)fSecondRunList->FindObject("First_FlowPro_r0theta_LYZ");
28ca24ad	128	fHistList->Add(fFirstr0theta);
f1d945a1	129
9d062fe3	130	//warnings
	131	if (!fSecondReDtheta) {cout<<"fSecondReDtheta is NULL!"<<endl; }
	132	if (!fSecondImDtheta) {cout<<"fSecondImDtheta is NULL!"<<endl; }
	133	if (!fFirstr0theta) {cout<<"fFirstr0theta is NULL!"<<endl; }
f1d945a1	134
9d062fe3	135	}
28ca24ad	136
9d062fe3	137	fCommonHists = new AliFlowCommonHist("AliFlowCommonHistLYZEP");
9d062fe3	138	fHistList->Add(fCommonHists);
f1d945a1	139
9d062fe3	140	fCommonHistsRes = new AliFlowCommonHistResults("AliFlowCommonHistResultsLYZEP");
	141	fHistList->Add(fCommonHistsRes);
	142
882ffd6a	143	Int_t iNbinsPt = AliFlowCommonConstants::GetNbinsPt();
	144	Double_t dPtMin = AliFlowCommonConstants::GetPtMin();
	145	Double_t dPtMax = AliFlowCommonConstants::GetPtMax();
f1d945a1	146
882ffd6a	147	fHistProFlow = new TProfile("FlowPro_VPt_LYZEP","FlowPro_VPt_LYZEP",iNbinsPt,dPtMin,dPtMax);
f1d945a1	148	fHistProFlow->SetXTitle("Pt");
f1d945a1	149	fHistProFlow->SetYTitle("v2 (%)");
28ca24ad	150	fHistList->Add(fHistProFlow);
f1d945a1	151
	152	fHistProWr = new TProfile("FlowPro_Wr_LYZEP","FlowPro_Wr_LYZEP",100,0.,0.25);
	153	fHistProWr->SetXTitle("Q");
	154	fHistProWr->SetYTitle("Wr");
28ca24ad	155	fHistList->Add(fHistProWr);
f1d945a1	156
9d062fe3	157	fHistQsumforChi = new TH1F("Flow_QsumforChi_LYZEP","Flow_QsumforChi_LYZEP",3,-1.,2.);
	158	fHistQsumforChi->SetXTitle("Qsum.X , Qsum.Y, Q2sum");
	159	fHistQsumforChi->SetYTitle("value");
	160	fHistList->Add(fHistQsumforChi);
	161
f1d945a1	162	fHistDeltaPhi = new TH1F("Flow_DeltaPhi_LYZEP","Flow_DeltaPhi_LYZEP",100,0.,3.14);
	163	fHistDeltaPhi->SetXTitle("DeltaPhi");
	164	fHistDeltaPhi->SetYTitle("Counts");
28ca24ad	165	fHistList->Add(fHistDeltaPhi);
f1d945a1	166
	167	fHistPhiLYZ = new TH1F("Flow_PhiLYZ_LYZEP","Flow_PhiLYZ_LYZEP",100,0.,3.14);
	168	fHistPhiLYZ->SetXTitle("Phi from LYZ");
	169	fHistPhiLYZ->SetYTitle("Counts");
28ca24ad	170	fHistList->Add(fHistPhiLYZ);
f1d945a1	171
	172	fHistPhiEP = new TH1F("Flow_PhiEP_LYZEP","Flow_PhiEP_LYZEP",100,0.,3.14);
	173	fHistPhiEP->SetXTitle("Phi from EP");
	174	fHistPhiEP->SetYTitle("Counts");
28ca24ad	175	fHistList->Add(fHistPhiEP);
f1d945a1	176
f1d945a1	177	fEventNumber = 0; //set number of events to zero
f1d945a1	178
	179	}
	180
	181	//-----------------------------------------------------------------------
	182
882ffd6a	183	void AliFlowAnalysisWithLYZEventPlane::Make(AliFlowEventSimple* anEvent, AliFlowLYZEventPlane* aLYZEP) {
882ffd6a	184
f1d945a1	185	//Get the event plane and weight for each event
8232a5ec	186	if (anEvent) {
f1d945a1	187
f1d945a1	188	//fill control histograms
8232a5ec	189	fCommonHists->FillControlHistograms(anEvent);
f1d945a1	190
f1d945a1	191	//get the Q vector from the FlowEvent
882ffd6a	192	AliFlowVector vQ = anEvent->GetQ();
9d062fe3	193	if (vQ.X()== 0. && vQ.Y()== 0. ) { cout<<"Q vector is NULL!"<<endl; }
f1d945a1	194	//Weight with the multiplicity
882ffd6a	195	Double_t dQX = 0.;
	196	Double_t dQY = 0.;
	197	if (vQ.GetMult()!=0.) {
	198	dQX = vQ.X()/vQ.GetMult();
	199	dQY = vQ.Y()/vQ.GetMult();
	200	} else {cerr<<"vQ.GetMult() is zero!"<<endl; }
	201	vQ.Set(dQX,dQY);
9d062fe3	202	//cout<<"vQ("<<dQX<<","<<dQY<<")"<<endl;
28ca24ad	203
f1d945a1	204	//for chi calculation:
882ffd6a	205	*fQsum += vQ;
9d062fe3	206	fHistQsumforChi->SetBinContent(1,fQsum->X());
9d062fe3	207	fHistQsumforChi->SetBinContent(2,fQsum->Y());
882ffd6a	208	fQ2sum += vQ.Mod2();
9d062fe3	209	fHistQsumforChi->SetBinContent(3,fQ2sum);
9d062fe3	210	//cout<<"fQ2sum = "<<fQ2sum<<endl;
f1d945a1	211
f1d945a1	212	//call AliFlowLYZEventPlane::CalculateRPandW() here!
882ffd6a	213	aLYZEP->CalculateRPandW(vQ);
f1d945a1	214
882ffd6a	215	Double_t dWR = aLYZEP->GetWR();
882ffd6a	216	Double_t dRP = aLYZEP->GetPsi();
f1d945a1	217
882ffd6a	218	//fHistProWr->Fill(vQ.Mod(),dWR); //this weight is always positive
882ffd6a	219	fHistPhiLYZ->Fill(dRP);
f1d945a1	220
f1d945a1	221	//plot difference between event plane from EP-method and LYZ-method
882ffd6a	222	Double_t dRPEP = vQ.Phi()/2; //gives distribution from (0 to pi)
	223	//Double_t dRPEP = 0.5*TMath::ATan2(vQ.Y(),vQ.X()); //gives distribution from (-pi/2 to pi/2)
	224	//cout<<"dRPEP = "<< dRPEP <<endl;
	225	fHistPhiEP->Fill(dRPEP);
f1d945a1	226
882ffd6a	227	Double_t dDeltaPhi = dRPEP - dRP;
	228	if (dDeltaPhi < 0.) { dDeltaPhi += TMath::Pi(); } //to shift distribution from (-pi/2 to pi/2) to (0 to pi)
	229	//cout<<"dDeltaPhi = "<<dDeltaPhi<<endl;
	230	fHistDeltaPhi->Fill(dDeltaPhi);
f1d945a1	231
f1d945a1	232	//Flip sign of WR
882ffd6a	233	Double_t dLow = TMath::Pi()/4.;
	234	Double_t dHigh = 3.*(TMath::Pi()/4.);
	235	if ((dDeltaPhi > dLow) && (dDeltaPhi < dHigh)){
	236	dRP -= (TMath::Pi()/2);
	237	dWR = -dWR;
	238	cerr<<"* dRP modified *"<<endl;
f1d945a1	239	}
882ffd6a	240	fHistProWr->Fill(vQ.Mod(),dWR); //corrected weight
f1d945a1	241
	242	//calculate flow
	243	//loop over the tracks of the event
882ffd6a	244	Int_t iNumberOfTracks = anEvent->NumberOfTracks();
882ffd6a	245	for (Int_t i=0;i<iNumberOfTracks;i++)
f1d945a1	246	{
882ffd6a	247	AliFlowTrackSimple* pTrack = anEvent->GetTrack(i) ;
	248	if (pTrack){
	249	if (pTrack->UseForDifferentialFlow()) {
	250	Double_t dPhi = pTrack->Phi();
	251	//if (dPhi<0.) fPhi+=2*TMath::Pi();
f1d945a1	252	//calculate flow v2:
882ffd6a	253	Double_t dv2 = dWR * TMath::Cos(2*(dPhi-dRP));
882ffd6a	254	Double_t dPt = pTrack->Pt();
f1d945a1	255	//fill histogram
882ffd6a	256	fHistProFlow->Fill(dPt,100*dv2);
f1d945a1	257	}
	258	}//track selected
	259	}//loop over tracks
	260
	261	fEventNumber++;
	262	cout<<"@@@@@ "<<fEventNumber<<" events processed"<<endl;
	263	}
	264	}
	265
	266	//--------------------------------------------------------------------
	267	void AliFlowAnalysisWithLYZEventPlane::Finish() {
	268
	269	//plot histograms etc.
bc92c0cb	270	cout<<"AliFlowAnalysisWithLYZEventPlane::Finish()"<<endl;
9d062fe3	271
9d062fe3	272	//constants:
882ffd6a	273	Double_t dJ01 = 2.405;
	274	Int_t iNtheta = AliFlowLYZConstants::kTheta;
	275	Int_t iNbinsPt = AliFlowCommonConstants::GetNbinsPt();
9d062fe3	276	//set the event number
	277	SetEventNumber((int)fCommonHists->GetHistMultOrig()->GetEntries());
	278	//cout<<"number of events processed is "<<fEventNumber<<endl;
	279
	280	//set the sum of Q vectors
	281	fQsum->Set(fHistQsumforChi->GetBinContent(1),fHistQsumforChi->GetBinContent(2));
	282	SetQ2sum(fHistQsumforChi->GetBinContent(3));
	283
882ffd6a	284	//calculate dV the mean of dVtheta
	285	Double_t dVtheta = 0;
	286	Double_t dV = 0;
	287	for (Int_t theta=0;theta<iNtheta;theta++) {
	288	Double_t dR0 = fFirstr0theta->GetBinContent(theta+1);
	289	if (dR0!=0.) { dVtheta = dJ01/dR0 ;}
	290	dV += dVtheta;
f1d945a1	291	}
882ffd6a	292	dV /= iNtheta;
f1d945a1	293
882ffd6a	294	//calculate dChi
	295	Double_t dSigma2 = 0;
	296	Double_t dChi= 0;
f1d945a1	297	if (fEventNumber!=0) {
0092f3c2	298	*fQsum /= fEventNumber;
bc92c0cb	299	//cerr<<"fQsum->X() = "<<fQsum->X()<<endl;
bc92c0cb	300	//cerr<<"fQsum->Y() = "<<fQsum->Y()<<endl;
f1d945a1	301	fQ2sum /= fEventNumber;
bc92c0cb	302	//cerr<<"fEventNumber = "<<fEventNumber<<endl;
bc92c0cb	303	//cerr<<"fQ2sum = "<<fQ2sum<<endl;
882ffd6a	304	dSigma2 = fQ2sum - TMath::Power(fQsum->X(),2.) - TMath::Power(fQsum->Y(),2.) - TMath::Power(dV,2.); //BP eq. 62
bc92c0cb	305	//cerr<<"dSigma2"<<dSigma2<<endl;
882ffd6a	306	if (dSigma2>0) dChi = dV/TMath::Sqrt(dSigma2);
	307	else dChi = -1.;
	308	fCommonHistsRes->FillChi(dChi);
	309	cerr<<"dV = "<<dV<<" and chi = "<<dChi<<endl;
f1d945a1	310	}
f1d945a1	311
882ffd6a	312	for(Int_t b=0;b<iNbinsPt;b++){
	313	Double_t dv2pro = 0.;
	314	Double_t dErr2difcomb = 0.;
	315	Double_t dErrdifcomb = 0.;
f1d945a1	316	if(fHistProFlow) {
882ffd6a	317	dv2pro = fHistProFlow->GetBinContent(b);
f1d945a1	318	//calculate error
882ffd6a	319	for (Int_t theta=0;theta<iNtheta;theta++) {
	320	Double_t dTheta = ((double)theta/iNtheta)*TMath::Pi();
	321	Int_t iNprime = TMath::Nint(fHistProFlow->GetBinEntries(b));
	322	//cerr<<"iNprime = "<<iNprime<<endl;
	323	if (iNprime!=0) {
	324	Double_t dApluscomb = TMath::Exp((dJ01dJ01)/(2dChidChi)
	325	TMath::Cos(dTheta));
	326	Double_t dAmincomb = TMath::Exp(-(dJ01dJ01)/(2dChidChi)
	327	TMath::Cos(dTheta));
	328	dErr2difcomb += (TMath::Cos(dTheta)/(4iNprimeTMath::BesselJ1(dJ01)*
	329	TMath::BesselJ1(dJ01)))*
	330	((dApluscombTMath::BesselJ0(2dJ01*TMath::Sin(dTheta/2))) -
	331	(dAmincombTMath::BesselJ0(2dJ01*TMath::Cos(dTheta/2))));
f1d945a1	332	} //if !=0
882ffd6a	333	//else { cout<<"iNprime = 0"<<endl; }
f1d945a1	334	} //loop over theta
f1d945a1	335
882ffd6a	336	if (dErr2difcomb!=0.) {
	337	dErr2difcomb /= iNtheta;
	338	dErrdifcomb = TMath::Sqrt(dErr2difcomb)*100;
	339	//cerr<<"dErrdifcomb = "<<dErrdifcomb<<endl;
f1d945a1	340	}
882ffd6a	341	else {dErrdifcomb = 0.; }
f1d945a1	342
f1d945a1	343	//fill TH1D
882ffd6a	344	fCommonHistsRes->FillDifferentialFlow(b, dv2pro, dErrdifcomb);
f1d945a1	345
	346	} //if fHistProFLow
	347	else {
	348	cout << "Profile Hist missing" << endl;
	349	break;
	350	}
	351
	352	} //loop over b
	353
f1d945a1	354	cout<<".....finished"<<endl;
f1d945a1	355	}
882ffd6a	356