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