/**************************************************************************
 * Copyright(c) 1998-1999, 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$
*/ 

#include "Riostream.h"
#include "AliFlowLeeYangZerosMaker.h"
#include "AliFlowEvent.h"
#include "AliFlowSelection.h"
#include "AliFlowConstants.h"     //??
#include "AliFlowLYZHist1.h"
#include "AliFlowLYZHist2.h"
#include "AliFlowLYZConstants.h"  //??
//#include "AliFlowLYZSummary.h"

#include "TMath.h"
#include "TComplex.h" 
#include "TProfile.h"
#include "TObjArray.h"
#include "TFile.h"
#include "TVector.h"
#include "TVector2.h"
#include "TGraphErrors.h"
#include "TCanvas.h"

class TTree;
class TH1F;
class TH1D;
class TVector3;
class TProfile2D;
class TObject;

//class Riostream; //does not compile
//class TMath;     //does not compile
//class TVector;   //does not compile

//Description: Maker to analyze Flow using the LeeYangZeros method  
//             Equation numbers are from Big Paper (BP): Nucl. Phys. A 727, 373 (2003)
//             Practical Guide (PG):    J. Phys. G: Nucl. Part. Phys. 30, S1213 (2004)  
//             Adapted from StFlowLeeYangZerosMaker.cxx           
//             by Markus Oldenberg and Art Poskanzer, LBNL        
//             with advice from Jean-Yves Ollitrault and Nicolas Borghini   
//
//Author: Naomi van der Kolk (kolk@nikhef.nl)



ClassImp(AliFlowLeeYangZerosMaker)

 //-----------------------------------------------------------------------
 
  AliFlowLeeYangZerosMaker::AliFlowLeeYangZerosMaker():
    fFirstRun(kTRUE),
    fUseSum(kTRUE),
    fDebug(kFALSE),
    fHistFileName(0),
    fHistFile(0),
    fSummaryFile(0),
    firstRunFile(0)

{
  //default constructor
  if (fDebug) cout<<"****AliFlowLeeYangZerosMaker::AliFlowLeeYangZerosMaker default constructor****"<<endl;

  fFlowSelect = new AliFlowSelection();
  if (fDebug) { cerr<<"****fFlowSelect in constructor AliFlowLeeYangZerosMaker ("<<fFlowSelect<<")****"<<endl;}
  // output file (histograms)
  TString fHistFileName = "flowLYZAnalysPlot.root" ;
}
 

AliFlowLeeYangZerosMaker::AliFlowLeeYangZerosMaker(const AliFlowSelection* flowSelect):
  fFirstRun(kTRUE),
  fUseSum(kTRUE),
  fDebug(kFALSE),
  fHistFileName(0),
  fHistFile(0),
  fSummaryFile(0),
  firstRunFile(0)
{
  //custum constructor
  if (fDebug) cout<<"****AliFlowLeeYangZerosMaker::AliFlowLeeYangZerosMaker custum constructor****"<<endl;

  if(flowSelect) { fFlowSelect = new AliFlowSelection(*flowSelect); }
  else 	  { 
    fFlowSelect = new AliFlowSelection() ; 
    if (fDebug) cerr<<"****fFlowSelect in constructor AliFlowLeeYangZerosMaker ("<<fFlowSelect<<")****"<<endl;
  }
  // output file (histograms)
  TString fHistFileName = "flowLYZAnalysPlot.root" ;
}
 
 //-----------------------------------------------------------------------


 AliFlowLeeYangZerosMaker::~AliFlowLeeYangZerosMaker() 
 {
   //default destructor
   if (fDebug) cout<<"****~AliFlowLeeYangZerosMaker****"<<endl;
   delete fHistFile;
 }
 
 //-----------------------------------------------------------------------

Bool_t AliFlowLeeYangZerosMaker::Init() 
{
  //init method 
  if (fDebug) cout<<"****AliFlowLeeYangZerosMaker::Init()****"<<endl;

  // Open output files (->plots)
  fHistFile = new TFile(fHistFileName.Data(), "RECREATE");
  //fHistFile->cd() ;  //all histograms will be saved in this file

  //for each harmonic ???
  fQsum.Set(0.,0.);
  fQ2sum = 0.;
  
  // Book histograms
  Int_t fNtheta = AliFlowLYZConstants::kTheta;
  
  
  //for control histograms
  fHistOrigMult = new TH1F("Control_FlowLYZ_OrigMult", "Control_FlowLYZ_OrigMult",1000, 0., 10000.);
  fHistOrigMult->SetXTitle("Original Multiplicity");
  fHistOrigMult->SetYTitle("Counts");

  fHistMult = new TH1F("Control_FlowLYZ_Mult", "Control_FlowLYZ_Mult",1000, 0., 10000.);
  fHistMult->SetXTitle("Multiplicity from selection");
  fHistMult->SetYTitle("Counts");

  fHistQ = new TH1F("Control_FlowLYZ_Q","Control_FlowLYZ_Q",500, 0., 10.);
  fHistQ->SetXTitle("Qvector");
  fHistQ->SetYTitle("Counts");

  fHistPt = new TH1F("Control_FlowLYZ_Pt","Control_FlowLYZ_Pt",200, 0., 10.);
  fHistPt->SetXTitle("Pt (GeV/c)");
  fHistPt->SetYTitle("Counts");

  fHistPhi = new TH1F("Control_FlowLYZ_Phi","Control_FlowLYZ_Phi",70, 0., 7.);
  fHistPhi->SetXTitle("Phi");
  fHistPhi->SetYTitle("Counts");

  fHistEta = new TH1F("Control_FlowLYZ_Eta","Control_FlowLYZ_Eta",40, 0., 2.);
  fHistEta->SetXTitle("Eta");
  fHistEta->SetYTitle("Counts");

  fHistQtheta = new TH1F("Control_FlowLYZ_Qtheta","Control_FlowLYZ_Qtheta",50,-1000.,1000.);
  fHistQtheta->SetXTitle("Qtheta");
  fHistQtheta->SetYTitle("Counts");
  
  if (fFirstRun){
    //for first loop over events
    fHistProR0thetaHar1  = new TProfile("First_FlowProLYZ_r0theta_Har1","First_FlowProLYZ_r0theta_Har1",fNtheta,-0.5,fNtheta-0.5);
    fHistProR0thetaHar1->SetXTitle("#theta");
    fHistProR0thetaHar1->SetYTitle("r_{0}^{#theta}");
  
    fHistProR0thetaHar2  = new TProfile("First_FlowProLYZ_r0theta_Har2","First_FlowProLYZ_r0theta_Har2",fNtheta,-0.5,fNtheta-0.5);
    fHistProR0thetaHar2->SetXTitle("#theta");
    fHistProR0thetaHar2->SetYTitle("r_{0}^{#theta}");

    fHistProVthetaHar1  = new TProfile("First_FlowProLYZ_Vtheta_Har1","First_FlowProLYZ_Vtheta_Har1",fNtheta,-0.5,fNtheta-0.5);
    fHistProVthetaHar1->SetXTitle("#theta");
    fHistProVthetaHar1->SetYTitle("V_{n}^{#theta}");

    fHistProVthetaHar2  = new TProfile("First_FlowProLYZ_Vtheta_Har2","First_FlowProLYZ_Vtheta_Har2",fNtheta,-0.5,fNtheta-0.5);
    fHistProVthetaHar2->SetXTitle("#theta");
    fHistProVthetaHar2->SetYTitle("V_{n}^{#theta}");

    fHistProVR0 = new TProfile("First_FlowProLYZ_vR0","First_FlowProLYZ_vR0",2,0.5,2.5,-100.,100.);
    fHistProVR0->SetXTitle("Harmonic");
    fHistProVR0->SetYTitle("v integrated from r0 (%)");

    fHistVR0 = new TH1D("First_FlowLYZ_vR0","First_FlowLYZ_vR0",2,0.5,2.5);
    fHistVR0->SetXTitle("Harmonic");
    fHistVR0->SetYTitle("v integrated from r0 (%)");

    fHistProV = new TProfile("First_FlowProLYZ_V","First_FlowProLYZ_V",2,0.5,2.5,-1000.,1000.);
    fHistProV->SetXTitle("Harmonic");
    fHistProV->SetYTitle("v integrated");
  
    //class AliFlowLYZHist1 defines the histograms: fHistProGtheta, fHistProReGtheta, fHistProImGtheta, fHistProR0theta
    for (Int_t j=0;j<AliFlowConstants::kHars;j++)  
      {
	for (Int_t theta=0;theta<fNtheta;theta++)
	  {  
	    fHist1[j][theta]=new AliFlowLYZHist1(theta,j+1);
	  }
      }
  }
  else {
    //for second loop over events
    fHistProReDenomHar1 = new TProfile("Second_FlowProLYZ_ReDenom_Har1","Second_FlowProLYZ_ReDenom_Har1" , fNtheta, -0.5, fNtheta-0.5);
    fHistProReDenomHar1->SetXTitle("#theta");
    fHistProReDenomHar1->SetYTitle("Re(Q^{#theta}e^{ir_{0}^{#theta}Q^{#theta}})");

    fHistProReDenomHar2 = new TProfile("Second_FlowProLYZ_ReDenom_Har2","Second_FlowProLYZ_ReDenom_Har2" , fNtheta, -0.5, fNtheta-0.5);
    fHistProReDenomHar2->SetXTitle("#theta");
    fHistProReDenomHar2->SetYTitle("Re(Q^{#theta}e^{ir_{0}^{#theta}Q^{#theta}})");

    fHistProImDenomHar1 = new TProfile("Second_FlowProLYZ_ImDenom_Har1","Second_FlowProLYZ_ImDenom_Har1" , fNtheta, -0.5, fNtheta-0.5);
    fHistProImDenomHar1->SetXTitle("#theta");
    fHistProImDenomHar1->SetYTitle("Im(Q^{#theta}e^{ir_{0}^{#theta}Q^{#theta}})");

    fHistProImDenomHar2 = new TProfile("Second_FlowProLYZ_ImDenom_Har2","Second_FlowProLYZ_ImDenom_Har2" , fNtheta, -0.5, fNtheta-0.5);
    fHistProImDenomHar2->SetXTitle("#theta");
    fHistProImDenomHar2->SetYTitle("Im(Q^{#theta}e^{ir_{0}^{#theta}Q^{#theta}})");

    fHistProVetaHar1 = new TProfile("Second_FlowProLYZ_Veta_Har1","Second_FlowProLYZ_Veta_Har1",40,-2.,2.);
    fHistProVetaHar1->SetXTitle("rapidity");
    fHistProVetaHar1->SetYTitle("v (%)");

    fHistProVetaHar2 = new TProfile("Second_FlowProLYZ_Veta_Har2","Second_FlowProLYZ_Veta_Har2",40,-2.,2.);
    fHistProVetaHar2->SetXTitle("rapidity");
    fHistProVetaHar2->SetYTitle("v (%)");

    fHistProVPtHar1 = new TProfile("Second_FlowProLYZ_VPt_Har1","Second_FlowProLYZ_VPt_Har1",100,0.,10.);
    fHistProVPtHar1->SetXTitle("Pt");
    fHistProVPtHar1->SetYTitle("v (%)");

    fHistProVPtHar2 = new TProfile("Second_FlowProLYZ_VPt_Har2","Second_FlowProLYZ_VPt_Har2",100,0.,10.);
    fHistProVPtHar2->SetXTitle("Pt");
    fHistProVPtHar2->SetYTitle("v (%)");

    fHistVPtHar2 = new TH1D("Second_FlowLYZ_VPt_Har2","Second_FlowLYZ_VPt_Har2",100,0.,10.);
    fHistVPtHar2->SetXTitle("Pt");
    fHistVPtHar2->SetYTitle("v (%)");

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

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


    //class AliFlowLYZHist2 defines the histograms: 
    for (Int_t j=0;j<AliFlowConstants::kHars;j++) 
      {
	for (Int_t theta=0;theta<fNtheta;theta++)
	  {  
	    fHist2[j][theta]=new AliFlowLYZHist2(theta,j+1);
	  }
      }

    //read hists from first run file
    //firstRunFile = new TFile("fof_flowLYZAnal_firstrun.root","READ");  //default is read
    if (firstRunFile->IsZombie()){ //check if file exists
      cout << "Error opening file, run first with fFirstrun = kTRUE" << endl;
      exit(-1);
    } else if (firstRunFile->IsOpen()){
      cout<<"----firstRunFile is open----"<<endl<<endl;
      fHistProVthetaHar1  = (TProfile*)firstRunFile->Get("First_FlowProLYZ_Vtheta_Har1");
      fHistProVthetaHar2  = (TProfile*)firstRunFile->Get("First_FlowProLYZ_Vtheta_Har2");
      fHistProR0thetaHar2  = (TProfile*)firstRunFile->Get("First_FlowProLYZ_r0theta_Har2");
      fHistProV = (TProfile*)firstRunFile->Get("First_FlowProLYZ_V");
    }    
  }
   

  if (fDebug) cout<<"****Histograms initialised****"<<endl;
  if (fDebug) cout<<"****fFlowSelect in Init() "<<fFlowSelect<<"****"<<endl;
   
  fEventNumber = 0; //set event counter to zero
  /*
  if (fUseSum)
    {
      //initialize LYZ summary class
      fLYZSummary = new AliFlowLYZSummary(); 
      fSummaryFile = new TFile("fFlowSummary.root","RECREATE","Flow LYZ summary file");
      fSummaryFile->SetCompressionLevel(1);
      fFlowTree = new TTree("FlowTree", "Flow Summary Tree");
      fFlowTree->SetAutoSave(1000000);  // autosave when 1 Mbyte written
      fFlowTree->Branch("fLYZSummary","AliFlowLYZSummary",&fLYZSummary,25000,99);
    }
  */
  return kTRUE; 
}
 
 //-----------------------------------------------------------------------
 
Bool_t AliFlowLeeYangZerosMaker::Make(AliFlowEvent* fFlowEvent) 
{
  //make method
  if (fDebug) cout<<"****AliFlowLeeYangZerosMaker::Make()****"<<endl;
        
  //get tracks from event
  if (fFlowEvent) {
    fFlowTracks = fFlowEvent->TrackCollection();
    if (fDebug) cout<<"****fFlowSelect in Make() "<<fFlowSelect<<"****"<<endl;
    if (fDebug) fFlowSelect->PrintSelectionList() ;
    if (fDebug) fFlowSelect->PrintList() ;
     
    if (fFlowSelect && fFlowSelect->Select(fFlowEvent))   // check if event is selected
      {  
	fFlowTracks = fFlowEvent->TrackCollection();     //get tracks from event
	fFlowEvent->SetSelections(fFlowSelect) ;		// does the selection of tracks for r.p. calculation (sets flags in AliFlowTrack)
	if (fFirstRun){
	  MakeControlHistograms(fFlowEvent);
	  FillFromFlowEvent(fFlowEvent);
	}
	else {
	  MakeControlHistograms(fFlowEvent);
	  SecondFillFromFlowEvent(fFlowEvent);
	}
      }
  }
  else {
    cout<<"##### FlowLeeYangZero: FlowEvent pointer null"<<endl;
    return kFALSE;
  }

  fEventNumber++;
     
  return kTRUE; 
}


  //-----------------------------------------------------------------------     
 Bool_t AliFlowLeeYangZerosMaker::Finish() 
{
  //finish method
  if (fDebug) cout<<"****AliFlowLeeYangZerosMaker::Finish()****"<<endl; 
  
 //define variables for both runs
  Float_t  fR0 = 0;
  Float_t  fv = 0;
  Float_t  fVtheta = 0; 
  Float_t  fSigma2 = 0;
  Float_t  fChi= 0;
  Float_t  fJ01 = 2.405; 
  Int_t fNtheta = AliFlowLYZConstants::kTheta;
  Float_t  fV = 0; 
  
  if (fFirstRun){
    for (Int_t j=0;j<AliFlowConstants::kHars;j++)    
      //loop over harmonics j=0 ->first harmonic, j=1 ->second harmonic
      {
	Float_t  fMeanMult = fHistMult->GetMean();
	//cerr<<"fMeanMult = "<<fMeanMult<<endl;

	for (Int_t theta=0;theta<fNtheta;theta++)
	  {
	    //get the first minimum r0
	    fHist1[j][theta]->FillGtheta();
	    fR0 = fHist1[j][theta]->GetR0();
	    //cerr<<"fR0 = "<<fR0<<endl;
	   	   
	    //calculate integrated flow
	    if (fR0!=0) fVtheta = fJ01/fR0;
	    if (fMeanMult!=0.) 
	      {
		fv = fVtheta/fMeanMult;
		cerr<<"fv = "<<fv<<endl;
	      }
	   

	    //fill the histograms
	    if (j==0) 
	      {
		fHistProR0thetaHar1->Fill(theta,fR0);
		fHistProVthetaHar1->Fill(theta,fVtheta);
		fHistProV->Fill(j+1,fVtheta);    //profile takes care of the averaging over theta.
	      }
	    if (j==1) 
	      {
		fHistProR0thetaHar2->Fill(theta,fR0);
		fHistProVthetaHar2->Fill(theta,fVtheta);
		fHistProV->Fill(j+1,fVtheta);   //profile takes care of the averaging over theta.
	      }
	     
	    fHistProVR0->Fill(j+1,fv*100);    //*100 to get %     //profile takes care of the averaging over theta.              
	  
	  } //end of loop over theta
       
	//sigma2 and chi 
	if (j==1)   //second harmonic only temporarily
	  { 
	    if (fEventNumber!=0) {
	      fQsum /= fEventNumber;
	      //cerr<<"fQsum.X() = "<<fQsum.X()<<endl;
	      //cerr<<"fQsum.Y() = "<<fQsum.Y()<<endl;
	      fQ2sum /= fEventNumber;
	      //cerr<<"fQ2sum = "<<fQ2sum<<endl;
	      fV = fHistProV->GetBinContent(j+1);
	      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.;
	      cerr<<"fV = "<<fV<<" and chi = "<<fChi<<endl;
	    }
	   
	  } //j==1

      } //end of loop over harmonics

    // recalculate statistical errors on integrated flow
    //combining 5 theta angles to 1 relative error BP eq. 89
    Double_t fRelErr2comb = 0.;
    Int_t nEvts = fEventNumber;
    for (Int_t theta=0;theta<fNtheta;theta++){
      Double_t fTheta = ((double)theta/fNtheta)*TMath::Pi(); 
      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));
      fRelErr2comb += (1/(2*nEvts*(fJ01*fJ01)*TMath::BesselJ1(fJ01)*
			  TMath::BesselJ1(fJ01)))*
	(fApluscomb*TMath::BesselJ0(2*fJ01*TMath::Sin(fTheta/2)) + 
	 fAmincomb*TMath::BesselJ0(2*fJ01)*TMath::Cos(fTheta/2));
    }
    fRelErr2comb /= fNtheta;
    Double_t fRelErrcomb = TMath::Sqrt(fRelErr2comb);
    cerr<<"fRelErrcomb = "<<fRelErrcomb<<endl;         

    //copy content of profile into TH1D and add error
    for(Int_t b=0;b<2;b++){
      Double_t fv2pro = fHistProVR0->GetBinContent(b+1);
      fHistVR0->SetBinContent(b+1,fv2pro);
      Double_t fv2Err = fv2pro*fRelErrcomb ; 
      cerr<<"fv2pro +- fv2Err = "<<fv2pro<<" +- "<<fv2Err<<" for bin "<<b+1<<endl;
      fHistVR0->SetBinError(b+1,fv2Err);
    }
   

    if (fDebug) cout<<"****histograms filled****"<<endl;  
    /*
      if (fUseSum)
      {
      if (fSummaryFile->IsOpen())
      {
      fSummaryFile->Write(0,TObject::kOverwrite);
      fSummaryFile->Close();
      }
      }
    */
    
    //save histograms in file //temp for testing selector
    fHistFile->cd();
    fHistFile->Write();

    return kTRUE;
    fEventNumber =0; //set to zero for second round over events
  }  //firstrun
   
 
  else {  //second run
   
    //calculate differential flow
    //declare variables
    Float_t fEta, fPt, fReRatio, fVeta, fVPt;
    Float_t fReDenom = 0;
    Float_t fImDenom = 0; 
    Double_t fR0 = 0;
    TComplex fDenom, fNumer, fDtheta;
    Int_t m = 1;
    TComplex i = TComplex::I();
    Float_t fBesselRatio[3] = {1., 1.202, 2.69};
    Double_t fErrdifcomb = 0.;  //set error to zero
    Double_t fErr2difcomb = 0.; //set error to zero
 
    /*
    //v1 integrated
    for (Int_t j=0;j<AliFlowConstants::kHars;j++)   //loop over harmonics j=0 ->first harmonic, j=1 ->second harmonic
    {
    for (Int_t theta=0;theta<fNtheta;theta++)
    {
    //get the first minimum r0
    //fR0 = GetR0(fHist1[j][theta]->FillGtheta()); 

    fReDenom = fHistProReDenomHar2->GetBinContent(theta+1);
    fImDenom = fHistProImDenomHar2->GetBinContent(theta+1);
    TComplex fDenom(fReDenom,fImDenom);
	  
    //complete!!

    }//end of loop over theta
    }//end of loop over harmonics
    */

    //differential flow
    //for (Int_t j=0;j<AliFlowConstants::kHars;j++)    //loop over harmonics j=0 ->first harmonic, j=1 ->second harmonic
    //{
    Int_t j=1; //temp only harm 2
    //fFlowSelect->SetHarmonic(j);  //not needed here?
     
    for (Int_t theta=0;theta<fNtheta;theta++)  { //loop over theta
      if (j==0) {
	fR0 = fHistProR0thetaHar1->GetBinContent(theta+1);
	fVtheta = 2.405/fR0;
	//fVtheta =  fHistProVthetaHar1->GetBinContent(theta+1);  // BP Eq. 9  -> Vn^theta = j01/r0^theta
	fReDenom = fHistProReDenomHar1->GetBinContent(theta+1);
	fImDenom = fHistProImDenomHar1->GetBinContent(theta+1);
      }
      if (j==1) {
	fR0 = fHistProR0thetaHar2->GetBinContent(theta+1);
	if (fDebug) cerr<<"fR0 = "<<fR0<<endl;
	fVtheta = 2.405/fR0;                                    // BP Eq. 9  -> Vn^theta = j01/r0^theta
	fReDenom = fHistProReDenomHar2->GetBinContent(theta+1);
	fImDenom = fHistProImDenomHar2->GetBinContent(theta+1);
      }
	 
      fDenom(fReDenom,fImDenom);
	 

      //for new method and use by others (only with the sum generating function):
      if (fUseSum) {
	fR0 = fHistProR0thetaHar2->GetBinContent(theta+1); 
	fDtheta = fR0*fDenom;
	fHistProReDtheta->Fill(theta,fDtheta.Re());
	fHistProImDtheta->Fill(theta,fDtheta.Im());
      }
	 
      fDenom *= TComplex::Power(i, m-1);
      //cerr<<"TComplex::Power(i, m-1) = "<<TComplex::Power(i, m-1).Rho()<<endl; //checked ok
	 
      //v as a function of eta	  
      Int_t fEtaBins = AliFlowLYZConstants::kEtaBins;
      for (Int_t be=1;be<=fEtaBins;be++)  {
	fEta = fHist2[j][theta]->GetBinCenter(be);
	fNumer = fHist2[j][theta]->GetfNumer(be);
	if (fNumer.Rho()==0) {
	  if (fDebug) cerr<<"WARNING: modulus of fNumer is zero in Finish()"<<endl;
	  fReRatio = 0;
	}
	else if (fDenom.Rho()==0) {
	  if (fDebug) cerr<<"WARNING: modulus of fDenom is zero"<<endl;
	  fReRatio = 0;
	}
	else {
	  //if ( j==1 && theta==0) cerr<<"modulus of fNumer = "<<fNumer.Rho() <<endl; //always a number smaller than 1, or 0.
	  fReRatio = (fNumer/fDenom).Re();
	}

	fVeta = fBesselRatio[m-1]*fReRatio*fVtheta*100.; //BP eq. 12
	//if ( j==1 && theta==0) cerr<<"eta = "<<fEta<<" cerr::fReRatio for eta = "<<fReRatio<<" cerr::fVeta for eta = "<<fVeta<<endl;
	      	      
	if (j==0) fHistProVetaHar1->Fill(fEta,fVeta);
	if (j==1) fHistProVetaHar2->Fill(fEta,fVeta);
      } //loop over bins be
	 
	//v as a function of Pt
      Int_t fPtBins = AliFlowLYZConstants::kPtBins;
      for (Int_t bp=1;bp<=fPtBins;bp++)  {
	fPt = fHist2[j][theta]->GetBinCenterPt(bp);
	fNumer = fHist2[j][theta]->GetfNumerPt(bp);
	if (fNumer.Rho()==0) {
	  if (fDebug) cerr<<"modulus of fNumer is zero"<<endl;
	  fReRatio = 0;
	}
	else if (fDenom.Rho()==0) {
	  if (fDebug) cerr<<"modulus of fDenom is zero"<<endl;
	  fReRatio = 0;
	}
	else {
	  //if ( j==1 && theta==0) cerr<<"modulus of fNumer = "<<fNumer.Rho() <<endl; //always a number smaller than 1, or 0.
	  fReRatio = (fNumer/fDenom).Re();
	}
   
	fVPt = fBesselRatio[m-1]*fReRatio*fVtheta*100.; //BP eq. 12
	//cerr<<"fBesselRatio[m-1] = "<<fBesselRatio[m-1]<<endl;   //checked ok
	//if ( j==1 && theta==0) cerr<<"pt = "<<fPt<<" cerr::fReRatio for pt = "<<fReRatio<<" cerr::fVPt for pt = "<<fVeta<<endl;
	      
	if (j==0) fHistProVPtHar1->Fill(fPt,fVPt);
	if (j==1) fHistProVPtHar2->Fill(fPt,fVPt);
      } //loop over bins bp

    }//end of loop over theta


    //sigma2 and chi (for statistical error calculations)
    if (j==1) {  //second harmonic only temporarily
      
      if (fEventNumber!=0) {
	fQsum /= fEventNumber;
	//cerr<<"fQsum.X() = "<<fQsum.X()<<endl;
	//cerr<<"fQsum.Y() = "<<fQsum.Y()<<endl;
	fQ2sum /= fEventNumber;
	cerr<<"fQ2sum = "<<fQ2sum<<endl;
	fV = fHistProV->GetBinContent(j+1);
	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.;
	cerr<<"fV = "<<fV<<" and chi = "<<fChi<<endl;
      }
	   
    } //j==1

    
    //copy content of profile into TH1D and add error
    for(Int_t b=0;b<100;b++){
      Double_t fv2pro = fHistProVPtHar2->GetBinContent(b);
      //calculate error
      for (Int_t theta=0;theta<fNtheta;theta++) {
	Double_t fTheta = ((double)theta/fNtheta)*TMath::Pi(); 
	Int_t Nprime = fHist2[j][theta]->GetNprimePt(b);
	//cerr<<"Nprime = "<<Nprime<<endl;
	if (Nprime!=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*Nprime*TMath::BesselJ1(fJ01)*
						 TMath::BesselJ1(fJ01)))*
	    ((fApluscomb*TMath::BesselJ0(2*fJ01*TMath::Sin(fTheta/2))) - 
	     (fAmincomb*TMath::BesselJ0(2*fJ01*TMath::Cos(fTheta/2))));
	}
      } //loop over theta
      
      if (fErr2difcomb!=0.) {
	fErr2difcomb /= fNtheta;
	fErrdifcomb = TMath::Sqrt(fErr2difcomb)*100;
	//cerr<<"fErrdifcomb = "<<fErrdifcomb<<endl;
      }
      else {fErrdifcomb = 0.;}
	  
      //fill TH1D
      if (j==1) {
	fHistVPtHar2->SetBinContent(b,fv2pro);
	fHistVPtHar2->SetBinError(b,fErrdifcomb);
      }
      //check that error is set
      //if (j==1) { cout<<"difference between calculated error and error in hostogram: "<< fErrdifcomb - fHistVPtHar2->GetBinError(b)<<endl; }
	    
    } //loop over bins b

	    
    //} //end of loop over harmonics         //temporarily out

    //save histograms in file
    fHistFile->cd();
    fHistFile->Write();
    fHistFile->Close();
    //Note that when the file is closed, all histograms and Trees in memory associated with this file are deleted
    if (fDebug) cout<<"****Histograms saved and fHistFile closed, all histograms deleted****"<<endl;
     
    //close the first run file 
    firstRunFile->Close();

     
  } //secondrun
   
  delete fFlowSelect;
  cout<<"----LYZ analysis finished....----"<<endl<<endl;

  return kTRUE;
}


//-----------------------------------------------------------------------
  Bool_t AliFlowLeeYangZerosMaker::MakeControlHistograms(AliFlowEvent* fFlowEvent) 
{ 
  //contol histograms of pt, eta, phi, Q, mult
  if (fDebug) cout<<"****AliFlowLeeYangZerosMaker::MakeControlHistograms()****"<<endl; 
   
  //define variables
  TVector2 fQ;
  Float_t fPt, fPhi, fEta, fQmult;

  if (!fFlowEvent){
    cout<<"##### FlowLeeYangZero: FlowEvent pointer null"<<endl;
    return kFALSE;
  }

  if (!fFlowSelect){
    cout<<"##### FlowLeeYangZero: FlowSelect pointer null"<<endl;
    return kFALSE;
  }

  //set selection and harmonic
  fFlowSelect->SetSelection(1); 
  fFlowSelect->SetHarmonic(1); //second harmonic
     
  //cerr<<"selection in MakeControlHistograms()"<<endl;
  //fFlowSelect->PrintSelectionList() ;
  //fFlowSelect->PrintList() ;

  fFlowTracks = fFlowEvent->TrackCollection();
  Int_t fNumberOfTracks = fFlowTracks->GetEntries();
  fHistOrigMult->Fill(fNumberOfTracks);
  //cerr<<"fNumberOfTracks = "<<fNumberOfTracks<<endl;
  Int_t fMult = fFlowEvent->Mult(fFlowSelect);  // Multiplicity of tracks in the specified Selection
  fHistMult->Fill(fMult);   
  //cerr<<"Mult = "<<fMult<<endl;
  
  fQ = fFlowEvent ->Q(fFlowSelect);
  fQmult = fQ.Mod()/TMath::Sqrt(fNumberOfTracks);
  fHistQ->Fill(fQmult);
   
  Int_t tempmult = 0; //for testing
  for (Int_t i=0;i<fNumberOfTracks;i++) 
    {
      fFlowTrack = (AliFlowTrack*)fFlowTracks->At(i) ;   //get object at array position i
      //if (fFlowSelect->SelectPart(fFlowTrack))           //if track is selected
      if (fFlowSelect->Select(fFlowTrack))  
	{
	  fPt = fFlowTrack->Pt();
	  fHistPt->Fill(fPt);
	  tempmult++;
	  fPhi = fFlowTrack->Phi();
	  if (fPhi<0.) fPhi+=2*TMath::Pi();
	  fHistPhi->Fill(fPhi);
	  fEta = fFlowTrack->Eta();
	  fHistEta->Fill(fEta);
	}
    }
  if (fMult!=tempmult){cerr<<"ERROR: Mult() is not tempmult! "<<fMult<<" :: "<<tempmult<<endl<<endl;}
  //else {cerr<<"Mult()= tempmult "<<fMult<<" :: "<<tempmult<<endl<<endl;}

  return kTRUE; 


}



//-----------------------------------------------------------------------
 
 Bool_t AliFlowLeeYangZerosMaker::FillFromFlowEvent(AliFlowEvent* fFlowEvent) 
{ 
  // Get event quantities from AliFlowEvent for all particles

  if (fDebug) cout<<"****AliFlowLeeYangZerosMaker::FillFromFlowEvent()****"<<endl;
   
  if (!fFlowEvent){
    cout<<"##### FlowLeeYangZero: FlowEvent pointer null"<<endl;
    return kFALSE;
  }
   
  if (!fFlowSelect){
    cout<<"##### FlowLeeYangZero: FlowSelect pointer null"<<endl;
    return kFALSE;
  }


  //define variables
  TComplex fExpo, fGtheta, fGthetaNew, fZ;
  //Int_t m;   
  Int_t fNtheta = AliFlowLYZConstants::kTheta;
  Int_t fNbins = AliFlowLYZConstants::kNbins;
   

  //calculate flow
  fFlowSelect->SetSelection(1); 
      
  for (Int_t j=0;j<AliFlowConstants::kHars;j++)   //loop over harmonics j=0 ->first harmonic, j=1 ->second harmonic
    {
      Int_t m=1;
      fFlowSelect->SetHarmonic(j);
      Float_t fOrder = (double)((j+1)/m);
      //cerr<<"fOrder = "<<fOrder<<endl;

      //get the Q vector from the FlowEvent
      TVector2 fQ = fFlowEvent->Q(fFlowSelect); 
      //for chi calculation:
      if (j==1)  //second harmonic only temporarily
	{
	  fQsum += fQ;
	  fQ2sum += fQ.Mod2();
	  //cerr<<"fQ2sum = "<<fQ2sum<<endl;
	}

      fFlowTracks = fFlowEvent->TrackCollection();

      for (Int_t theta=0;theta<fNtheta;theta++)
	{
	  fTheta = ((float)theta/fNtheta)*TMath::Pi()/fOrder; 
	  //cerr<<"fTheta = "<<fTheta<<endl;
	   
	  //calculate fQtheta = cos(fOrder*(fPhi-fTheta);the projection of the Q vector on the reference direction fTheta
	  fQtheta = GetQtheta(fFlowSelect,fFlowTracks,fTheta);
	  //save fQtheta in AliFlowLYZSummary class

	  //something
	  //AliFlowLYZSummary::SetQtheta(theta,fQtheta);

	  if (j==1 && theta==0) fHistQtheta->Fill(fQtheta);
	  //cerr<<"fQtheta = "<<fQtheta<<endl;
	   
	  for (Int_t bin=1;bin<=fNbins;bin++)
	    {
	      Float_t fR = fHist1[j][theta]->GetBinCenter(bin); //bincentre of bins in histogram
	      //if (theta == 0) cerr<<"cerr::fR = "<<fR<<endl;
	      if (fUseSum)
		{
		  //calculate the sum generating function
		  fExpo(0.,fR*fQtheta);                           //Re=0 ; Im=fR*fQtheta
		  fGtheta = TComplex::Exp(fExpo);
		}
	      else
		{
		  //calculate the product generating function
		  fGtheta = GetGrtheta(fFlowSelect,fR,fTheta);  //make this function
		  if (fGtheta.Rho2() > 100.) break;
		}

	      fHist1[j][theta]->Fill(fR,fGtheta);              //fill real and imaginary part of fGtheta

	    } //loop over bins
	} //loop over theta 
    } //loop over harmonics
   
  return kTRUE;

          
}

 //-----------------------------------------------------------------------   
 Bool_t AliFlowLeeYangZerosMaker::SecondFillFromFlowEvent(AliFlowEvent* fFlowEvent) 
{ 
  //for differential flow

  if (fDebug) cout<<"****AliFlowLeeYangZerosMaker::SecondFillFromFlowEvent()****"<<endl;
    
  if (!fFlowEvent){
    cout<<"##### FlowLeeYangZero: FlowEvent pointer null"<<endl;
    return kFALSE;
  }
   
  if (!fFlowSelect){
    cout<<"##### FlowLeeYangZero: FlowSelect pointer null"<<endl;
    return kFALSE;
  }

  //define variables
  //TVector2 fQ;
  TComplex fExpo, fDenom, fNumer,fCosTermComplex;
  Float_t  fOrder, fR0, fPhi, fEta, fPt;
  Double_t fCosTerm;
  Double_t m = 1.;
  Int_t fNtheta = AliFlowLYZConstants::kTheta;
   
  //calculate flow
  fFlowSelect->SetSelection(1);

  //cerr<<"selection in SecondFillFromFlowEvent()"<<endl;
  //fFlowSelect->PrintSelectionList() ;
  //fFlowSelect->PrintList() ;

      
  for (Int_t j=0;j<AliFlowConstants::kHars;j++)   //loop over harmonics j=0 ->first harmonic, j=1 ->second harmonic
    {
      m=1;
      fFlowSelect->SetHarmonic(j);
      fOrder = (double)((j+1)/m);

      //get the Q vector from the FlowEvent
      TVector2 fQ = fFlowEvent->Q(fFlowSelect); 
      //for chi calculation:
      if (j==1)  //second harmonic only temporarily
	{
	  fQsum += fQ;
	  fQ2sum += fQ.Mod2();
	  //cerr<<"fQ2sum = "<<fQ2sum<<endl;
	}

      fFlowTracks = fFlowEvent->TrackCollection();
       
      for (Int_t theta=0;theta<fNtheta;theta++)
	{
	  fTheta = ((float)theta/fNtheta)*TMath::Pi()/fOrder;   

	  //calculate fQtheta = cos(fOrder*(fPhi-fTheta);the projection of the Q vector on the reference direction fTheta
	  fQtheta = GetQtheta(fFlowSelect,fFlowTracks,fTheta);
	  //cerr<<"fQtheta for fdenom = "<<fQtheta<<endl;
	   
	  /*
	    if (j==0)  //first harmonic
	    {
	    //denominator for differential v
	    fR0 = fHistProR0thetaHar1->GetBinContent(theta+1);
	    fExpo(0.,fR0*fQtheta);
	    fDenom = fQtheta*(TComplex::Exp(fExpo)); //BP eq 12
	    //cerr<<"fDenom.Re() = "<<fDenom.Re()<<endl;
	    //cerr<<"fDenom.Im() = "<<fDenom.Im()<<endl;

	    //denominator for differential v
	    // ****put in product generating function!!

	       
	    fHistProReDenomHar1->Fill(theta,fDenom.Re());               //fill the real part of fDenom
	    fHistProImDenomHar1->Fill(theta,fDenom.Im());               //fill the imaginary part of fDenom
	    }
	  */

	  if (j==1)  //second harmonic
	    {
	      //denominator for differential v
	      fR0 = fHistProR0thetaHar2->GetBinContent(theta+1);
	      //cerr<<"fR0 = "<<fR0 <<endl;

	      if (fUseSum)                                                    //sum generating function
		{
		  fExpo(0.,fR0*fQtheta);
		  fDenom = fQtheta*(TComplex::Exp(fExpo)); //BP eq 12
		  //loop over tracks in event
		  fFlowTracks = fFlowEvent->TrackCollection();
		  Int_t fNumberOfTracks = fFlowTracks->GetEntries();
		  for (Int_t i=0;i<fNumberOfTracks;i++) 
		    {
		      fFlowTrack = (AliFlowTrack*)fFlowTracks->At(i) ;                //get object at array position i
		      if (fFlowSelect->SelectPart(fFlowTrack))                        //if track is selected
			{
			  fPhi = fFlowTrack->Phi();
			  fCosTerm = cos(m*fOrder*(fPhi-fTheta));
			  //cerr<<"fCosTerm = "<<fCosTerm <<endl;
			  fNumer = fCosTerm*(TComplex::Exp(fExpo));
			  if (fNumer.Rho()==0) {cerr<<"WARNING: modulus of fNumer is zero in SecondFillFromFlowEvent"<<endl;}
			  if (fDebug) cerr<<"modulus of fNumer is "<<fNumer.Rho()<<endl;
			  fEta = fFlowTrack->Eta();  //rapidity
			  fPt = fFlowTrack->Pt();
			  fHist2[j][theta]->Fill(fEta,fPt,fNumer);
			} //if
		    } //loop over tracks
		   
		} //sum
	      else                                                        //product generating function
		{
		  fDenom = GetDiffFlow(fFlowSelect,fR0,theta); 
		   
		}//product
	      
	      fHistProReDenomHar2->Fill(theta,fDenom.Re());               //fill the real part of fDenom
	      fHistProImDenomHar2->Fill(theta,fDenom.Im());               //fill the imaginary part of fDenom
	    }//j==1
   	   	   
	   	       
	}//end of loop over theta

    }//loop over harmonics

  
  return kTRUE;
 
  
 
}
 //-----------------------------------------------------------------------   
 Double_t AliFlowLeeYangZerosMaker::GetQtheta(AliFlowSelection* fFlowSelect, TObjArray* fFlowTracks, Float_t fTheta) 
{
  //calculate Qtheta. Qtheta is the sum over all particles of cos(fOrder*(fPhi-fTheta)) BP eq. 3
  if (fDebug) cout<<"****AliFlowLeeYangZerosMaker::GetQtheta()****"<<endl;

  if (!fFlowSelect){
     cout<<"##### FlowLeeYangZero: FlowSelect pointer null"<<endl;
     return kFALSE;
   }


  Double_t fQtheta = 0.;
  Int_t fHarmonic = fFlowSelect->Har(); 
  Double_t fOrder = (double)(fHarmonic+1);
       
  Int_t fNumberOfTracks = fFlowTracks->GetEntries();
  //cerr<<"GetQtheta::fNumberOfTracks = "<<fNumberOfTracks<<endl;
  
  for (Int_t i=0;i<fNumberOfTracks;i++)                  //loop over tracks in event
    {
      fFlowTrack = (AliFlowTrack*)fFlowTracks->At(i) ;   //get object at array position i
      if(fFlowSelect->Select(fFlowTrack))                //if track is selected  //gives the same number of particles as Mult(fFlowSelect) method
	{
     	  Float_t fPhi = fFlowTrack->Phi();
	  fQtheta += cos(fOrder*(fPhi-fTheta));
	}
     
    }//loop over tracks
  
  return fQtheta;
 
}
 
//-----------------------------------------------------------------------   
TComplex AliFlowLeeYangZerosMaker::GetGrtheta(AliFlowSelection* fFlowSelect, Float_t fR, Float_t fTheta) 
{
  // Product Generating Function for LeeYangZeros method
  // PG Eq. 3 (J. Phys. G Nucl. Part. Phys 30 S1213 (2004))

  if (fDebug) cout<<"****AliFlowLeeYangZerosMaker::GetGrtheta()****"<<endl;
  
  if (!fFlowSelect){
    cout<<"##### FlowLeeYangZero: FlowSelect pointer null"<<endl;
    return kFALSE;
  }
  
  TComplex fG = TComplex::One();
  Int_t fHarmonic = fFlowSelect->Har(); 
  Double_t fOrder = (double)(fHarmonic+1);
  Double_t fWgt = 1.;

  Int_t fNumberOfTracks = fFlowTracks->GetEntries();
  //cerr<<"GetGrtheta::fNumberOfTracks = "<<fNumberOfTracks<<endl; 

  for (Int_t i=0;i<fNumberOfTracks;i++) //loop over tracks in event
    {
      fFlowTrack = (AliFlowTrack*)fFlowTracks->At(i) ;                //get object at array position i
      if (fFlowSelect->SelectPart(fFlowTrack))                        //if track is selected
	{
	  Float_t fPhi = fFlowTrack->Phi();
	  Double_t fGIm = fR * fWgt*cos(fOrder*(fPhi - fTheta));
	  TComplex fGi(1., fGIm);
	  fG *= fGi;     //product over all tracks
	}//if
    }//loop over tracks

  return fG;

  } 


//-----------------------------------------------------------------------   
TComplex AliFlowLeeYangZerosMaker::GetDiffFlow(AliFlowSelection* fFlowSelect, Float_t fR0, Int_t theta) 
{
  // Sum for the denominator for diff. flow for the Product Generating Function for LeeYangZeros method
  // PG Eq. 9 (J. Phys. G Nucl. Part. Phys 30 S1213 (2004))
  // Also for v1 mixed harmonics: DF Eq. 5
  // It is the deriverative of Grtheta at r0 divided by Grtheta at r0

  if (fDebug) cout<<"****AliFlowLeeYangZerosMaker::GetGrtheta()****"<<endl;

  if (!fFlowSelect){
    cout<<"##### FlowLeeYangZero: FlowSelect pointer null"<<endl;
    return kFALSE;
  }


  TComplex fG = TComplex::One();
  TComplex fdGr0(0.,0.);
  Int_t fHarmonic = fFlowSelect->Har(); 
  Double_t fOrder = (double)(fHarmonic+1);
  Double_t fWgt = 1.;

  Int_t fNumberOfTracks = fFlowTracks->GetEntries();
  //cerr<<"GetGrtheta::fNumberOfTracks = "<<fNumberOfTracks<<endl; 

  Int_t fNtheta = AliFlowLYZConstants::kTheta;
  Float_t fTheta = ((float)theta/fNtheta)*TMath::Pi()/fOrder;

  for (Int_t i=0;i<fNumberOfTracks;i++) //loop over tracks in event
    {
      fFlowTrack = (AliFlowTrack*)fFlowTracks->At(i) ;                //get object at array position i
      if (fFlowSelect->SelectPart(fFlowTrack))                        //if track is selected
	{
	  Float_t fPhi = fFlowTrack->Phi();
	  Double_t fCosTerm = fWgt*cos(fOrder*(fPhi - fTheta));
	  //GetGr0theta
	  Double_t fGIm = fR0 * fCosTerm;
	  TComplex fGi(1., fGIm);
	  fG *= fGi;     //product over all tracks
	  //GetdGr0theta
	  TComplex fCosTermComplex(1., fR0*fCosTerm);
	  fdGr0 +=(fCosTerm / fCosTermComplex);  //sum over all tracks
	}//if
    }//loop over tracks

  for (Int_t i=0;i<fNumberOfTracks;i++) 
    {
      fFlowTrack = (AliFlowTrack*)fFlowTracks->At(i) ;                //get object at array position i
      if (fFlowSelect->SelectPart(fFlowTrack))                        //if track is selected
	{
	  Float_t fPhi = fFlowTrack->Phi();
	  Double_t fCosTerm = cos(fOrder*(fPhi-fTheta));
	  TComplex fCosTermComplex(1.,fR0*fCosTerm);
	  
	  TComplex fNumer = fG*fCosTerm/fCosTermComplex;  //PG Eq. 9
	  Float_t fEta = fFlowTrack->Eta();  //rapidity
	  Float_t fPt = fFlowTrack->Pt();
	  fHist2[1][theta]->Fill(fEta,fPt,fNumer);
	}//if
    }//loop over tracks

  TComplex fDenom = fG*fdGr0;
  return fDenom;

  } 

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