[u/mrichter/AliRoot.git] / PWG2 / FLOW / AliFlowCommon / AliFlowAnalysisWithMCEventPlane.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.                  * 
**************************************************************************/

#define AliFlowAnalysisWithMCEventPlane_cxx
 
#include "Riostream.h"  //needed as include
#include "TFile.h"      //needed as include
#include "TProfile.h"   //needed as include
#include "TProfile2D.h"   
#include "TComplex.h"   //needed as include
#include "TList.h"

class TH1F;

#include "AliFlowCommonConstants.h"    //needed as include
#include "AliFlowEventSimple.h"
#include "AliFlowTrackSimple.h"
#include "AliFlowCommonHist.h"
#include "AliFlowCommonHistResults.h"
#include "AliFlowAnalysisWithMCEventPlane.h"

class AliFlowVector;

// AliFlowAnalysisWithMCEventPlane:
// Description: Maker to analyze Flow from the generated MC reaction plane.
//              This class is used to get the real value of the flow 
//              to compare the other methods to when analysing simulated events
// author: N. van der Kolk (kolk@nikhef.nl)

ClassImp(AliFlowAnalysisWithMCEventPlane)

  //-----------------------------------------------------------------------
 
 AliFlowAnalysisWithMCEventPlane::AliFlowAnalysisWithMCEventPlane():
   fQsum(NULL),
   fQ2sum(0),
   fEventNumber(0),
   fDebug(kFALSE),
   fHistList(NULL),
   fCommonHists(NULL),
   fCommonHistsRes(NULL),
   fHistRP(NULL),
   fHistProIntFlow(NULL),
   fHistProDiffFlowPtEtaRP(NULL),
   fHistProDiffFlowPtRP(NULL),
   fHistProDiffFlowEtaRP(NULL),
   fHistProDiffFlowPtEtaPOI(NULL),
   fHistProDiffFlowPtPOI(NULL),
   fHistProDiffFlowEtaPOI(NULL)
{

  // Constructor.
  fHistList = new TList();

  fQsum = new TVector2;        // flow vector sum
}

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


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

void AliFlowAnalysisWithMCEventPlane::WriteHistograms(TString* outputFileName)
{
 //store the final results in output .root file
 TFile *output = new TFile(outputFileName->Data(),"RECREATE");
 output->WriteObject(fHistList, "cobjMCEP","SingleKey");
 delete output;
}

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

void AliFlowAnalysisWithMCEventPlane::WriteHistograms(TString outputFileName)
{
 //store the final results in output .root file
 TFile *output = new TFile(outputFileName.Data(),"RECREATE");
 output->WriteObject(fHistList, "cobjMCEP","SingleKey");
 delete output;
}

//-----------------------------------------------------------------------
void AliFlowAnalysisWithMCEventPlane::Init() {

  //Define all histograms
  cout<<"---Analysis with the real MC Event Plane---"<<endl;

  Int_t iNbinsPt = AliFlowCommonConstants::GetNbinsPt();
  Double_t dPtMin = AliFlowCommonConstants::GetPtMin();	     
  Double_t dPtMax = AliFlowCommonConstants::GetPtMax();
  
  Int_t iNbinsEta = AliFlowCommonConstants::GetNbinsEta();
  Double_t dEtaMin = AliFlowCommonConstants::GetEtaMin();	     
  Double_t dEtaMax = AliFlowCommonConstants::GetEtaMax();  

  fCommonHists = new AliFlowCommonHist("AliFlowCommonHistMCEP");
  fHistList->Add(fCommonHists);
  fCommonHistsRes = new AliFlowCommonHistResults("AliFlowCommonHistResultsMCEP");
  fHistList->Add(fCommonHistsRes);
  
  fHistRP = new TH1F("Flow_RP_MCEP","Flow_RP_MCEP",100,0.,3.14);
  fHistRP->SetXTitle("Reaction Plane Angle");
  fHistRP->SetYTitle("Counts");
  fHistList->Add(fHistRP);
  
  fHistProIntFlow = new TProfile("FlowPro_V_MCEP","FlowPro_V_MCEP",1,0.,1.);
  fHistProIntFlow->SetLabelSize(0.06);
  (fHistProIntFlow->GetXaxis())->SetBinLabel(1,"v_{n}{2}");
  fHistProIntFlow->SetYTitle("");
  fHistList->Add(fHistProIntFlow);
  
  fHistProDiffFlowPtEtaRP = new TProfile2D("FlowPro_VPtEtaRP_MCEP","FlowPro_VPtEtaRP_MCEP",iNbinsPt,dPtMin,dPtMax,iNbinsEta,dEtaMin,dEtaMax);
  fHistProDiffFlowPtEtaRP->SetXTitle("P_{t}");
  fHistProDiffFlowPtEtaRP->SetYTitle("#eta");
  fHistList->Add(fHistProDiffFlowPtEtaRP);
 
  fHistProDiffFlowPtRP = new TProfile("FlowPro_VPtRP_MCEP","FlowPro_VPtRP_MCEP",iNbinsPt,dPtMin,dPtMax);
  fHistProDiffFlowPtRP->SetXTitle("P_{t}");
  fHistProDiffFlowPtRP->SetYTitle("");
  fHistList->Add(fHistProDiffFlowPtRP);  
  
  fHistProDiffFlowEtaRP = new TProfile("FlowPro_VetaRP_MCEP","FlowPro_VetaRP_MCEP",iNbinsEta,dEtaMin,dEtaMax);
  fHistProDiffFlowEtaRP->SetXTitle("#eta");
  fHistProDiffFlowEtaRP->SetYTitle("");
  fHistList->Add(fHistProDiffFlowEtaRP);
  
  fHistProDiffFlowPtEtaPOI = new TProfile2D("FlowPro_VPtEtaPOI_MCEP","FlowPro_VPtEtaPOI_MCEP",iNbinsPt,dPtMin,dPtMax,iNbinsEta,dEtaMin,dEtaMax);
  fHistProDiffFlowPtEtaPOI->SetXTitle("P_{t}");
  fHistProDiffFlowPtEtaPOI->SetYTitle("#eta");
  fHistList->Add(fHistProDiffFlowPtEtaPOI);
  
  fHistProDiffFlowPtPOI = new TProfile("FlowPro_VPtPOI_MCEP","FlowPro_VPtPOI_MCEP",iNbinsPt,dPtMin,dPtMax);
  fHistProDiffFlowPtPOI->SetXTitle("P_{t}");
  fHistProDiffFlowPtPOI->SetYTitle("");
  fHistList->Add(fHistProDiffFlowPtPOI);  
  
  fHistProDiffFlowEtaPOI = new TProfile("FlowPro_VetaPOI_MCEP","FlowPro_VetaPOI_MCEP",iNbinsEta,dEtaMin,dEtaMax);
  fHistProDiffFlowEtaPOI->SetXTitle("#eta");
  fHistProDiffFlowEtaPOI->SetYTitle("");
  fHistList->Add(fHistProDiffFlowEtaPOI);          
 
  fEventNumber = 0;  //set number of events to zero
      
} 
 
//-----------------------------------------------------------------------
 
void AliFlowAnalysisWithMCEventPlane::Make(AliFlowEventSimple* anEvent) {

  //Calculate v2 from the MC reaction plane
  if (anEvent) {
  
    // get the MC reaction plane angle
    Double_t aRP = anEvent->GetMCReactionPlaneAngle();  
    //fill control histograms     
    fCommonHists->FillControlHistograms(anEvent);

    //get the Q vector from the FlowEvent
    AliFlowVector vQ = anEvent->GetQ(); 
    //cout<<"vQ.Mod() = " << vQ.Mod() << endl;
    //for chi calculation:
    *fQsum += vQ;
    //cout<<"fQsum.Mod() = "<<fQsum.Mod()<<endl;
    fQ2sum += vQ.Mod2();
    //cout<<"fQ2sum = "<<fQ2sum<<endl;
        
    fHistRP->Fill(aRP);   
    
    Double_t dPhi = 0.;
    Double_t dv2  = 0.;
    Double_t dPt  = 0.;
    Double_t dEta = 0.;
    //Double_t dPi = TMath::Pi();                   
                                                             
    //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->InRPSelection()){
            dPhi = pTrack->Phi();
            dv2  = TMath::Cos(2*(dPhi-aRP));
	    dPt  = pTrack->Pt();
	    dEta = pTrack->Eta();
            //no-name int. flow (to be improved):
            fHistProIntFlow->Fill(0.,dv2);
            //differential flow (Pt, Eta, RP):
            fHistProDiffFlowPtEtaRP->Fill(dPt,dEta,dv2,1.);
            //differential flow (Pt, RP):
            fHistProDiffFlowPtRP->Fill(dPt,dv2,1.);
            //differential flow (Eta, RP):
            fHistProDiffFlowEtaRP->Fill(dEta,dv2,1.);
          }
	  if (pTrack->InPOISelection()) {
	    dPhi = pTrack->Phi();
	    //if (dPhi<0.) dPhi+=2*TMath::Pi();
	    //calculate flow v2:
	    dv2  = TMath::Cos(2*(dPhi-aRP));
	    dPt  = pTrack->Pt();
	    dEta = pTrack->Eta();
	    //differential flow (Pt, Eta, POI):
            fHistProDiffFlowPtEtaPOI->Fill(dPt,dEta,dv2,1.);
	    //differential flow (Pt, POI):
            fHistProDiffFlowPtPOI->Fill(dPt,dv2,1.);
            //differential flow (Eta, POI):
            fHistProDiffFlowEtaPOI->Fill(dEta,dv2,1.); 
	  }	      
	}//track selected
      }//loop over tracks
	  
    fEventNumber++;
    //    cout<<"@@@@@ "<<fEventNumber<<" events processed"<<endl;
  }
}

  //--------------------------------------------------------------------    
void AliFlowAnalysisWithMCEventPlane::Finish() {
   
  //*************make histograms etc. 
  if (fDebug) cout<<"AliFlowAnalysisWithMCEventPlane::Terminate()"<<endl;
   
  Int_t iNbinsPt  = AliFlowCommonConstants::GetNbinsPt();  
  Int_t iNbinsEta = AliFlowCommonConstants::GetNbinsEta(); 
         
  // no-name int. flow (to be improved):
  Double_t dV = fHistProIntFlow->GetBinContent(1);  
  Double_t dErrV = fHistProIntFlow->GetBinError(1); // to be improved (treatment of errors for non-Gaussian distribution needed!)  
  // fill no-name int. flow (to be improved):
  fCommonHistsRes->FillIntegratedFlow(dV,dErrV);
  cout<<"dV{MC} is       "<<dV<<" +- "<<dErrV<<endl;
  
  //RP:
  TH1F* fHistPtRP = fCommonHists->GetHistPtRP(); 
  Double_t dYieldPtRP = 0.;
  Double_t dVRP = 0.;
  Double_t dErrVRP = 0.;
  Double_t dSumRP = 0.;
  //differential flow (RP, Pt): 
  Double_t dvPtRP = 0.;           
  Double_t dErrvPtRP = 0.;
  for(Int_t b=1;b<iNbinsPt;b++)
  {
   dvPtRP    = fHistProDiffFlowPtRP->GetBinContent(b);
   dErrvPtRP = fHistProDiffFlowPtRP->GetBinError(b);//to be improved (treatment of errors for non-Gaussian distribution needed!)
   fCommonHistsRes->FillDifferentialFlowPtRP(b, dvPtRP, dErrvPtRP);
   if(fHistPtRP){
	//integrated flow (RP)
	dYieldPtRP = fHistPtRP->GetBinContent(b);
	dVRP += dvPtRP*dYieldPtRP;
	dSumRP += dYieldPtRP;
	//error on integrated flow
	dErrVRP += dYieldPtRP*dYieldPtRP*dErrvPtRP*dErrvPtRP;
      }
  }
  if (dSumRP != 0. ) {
    dVRP /= dSumRP;  //because pt distribution should be normalised
    dErrVRP /= (dSumRP*dSumRP);
    dErrVRP = TMath::Sqrt(dErrVRP); 
  }
  // fill integrated flow (RP):
  fCommonHistsRes->FillIntegratedFlowRP(dVRP,dErrVRP);
  cout<<"dV{MC} (RP) is  "<<dVRP<<" +- "<<dErrVRP<<endl;
  
  //differential flow (RP, Eta): 
  Double_t dvEtaRP = 0.;           
  Double_t dErrvEtaRP = 0.;
  for(Int_t b=1;b<iNbinsEta;b++)
  {
   dvEtaRP    = fHistProDiffFlowEtaRP->GetBinContent(b);
   dErrvEtaRP = fHistProDiffFlowEtaRP->GetBinError(b);//to be improved (treatment of errors for non-Gaussian distribution needed!)
   fCommonHistsRes->FillDifferentialFlowEtaRP(b, dvEtaRP, dErrvEtaRP);
  }
                                                                                                                                   
  //POI:
  TH1F* fHistPtPOI = fCommonHists->GetHistPtPOI(); 
  Double_t dYieldPtPOI = 0.;
  Double_t dVPOI = 0.;
  Double_t dErrVPOI = 0.;
  Double_t dSumPOI = 0.;
  Double_t dv2proPtPOI = 0.;
  Double_t dErrdifcombPtPOI = 0.; 
  Double_t dv2proEtaPOI = 0.;
  Double_t dErrdifcombEtaPOI = 0.;   
  //Pt:
  if(fHistProDiffFlowPtPOI) {
    for(Int_t b=1;b<iNbinsPt;b++){
      dv2proPtPOI = fHistProDiffFlowPtPOI->GetBinContent(b);
      dErrdifcombPtPOI = fHistProDiffFlowPtPOI->GetBinError(b);//to be improved (treatment of errors for non-Gaussian distribution needed!)
      //fill TH1D
      fCommonHistsRes->FillDifferentialFlowPtPOI(b, dv2proPtPOI, dErrdifcombPtPOI); 
      if (fHistPtPOI){
	//integrated flow (POI)
	dYieldPtPOI = fHistPtPOI->GetBinContent(b);
	dVPOI += dv2proPtPOI*dYieldPtPOI;
	dSumPOI += dYieldPtPOI;
	//error on integrated flow
	dErrVPOI += dYieldPtPOI*dYieldPtPOI*dErrdifcombPtPOI*dErrdifcombPtPOI;
      }
    }//end of for(Int_t b=0;b<iNbinsPt;b++)  
  } else { cout<<"fHistProFlow is NULL"<<endl; }
  if (dSumPOI != 0. ) {
    dVPOI /= dSumPOI;  //because pt distribution should be normalised
    dErrVPOI /= (dSumPOI*dSumPOI);
    dErrVPOI = TMath::Sqrt(dErrVPOI); 
  }
  cout<<"dV{MC} (POI) is "<<dVPOI<<" +- "<<dErrVPOI<<endl;

  fCommonHistsRes->FillIntegratedFlowPOI(dVPOI,dErrVPOI);
  
  //Eta:
  if(fHistProDiffFlowEtaPOI)
  {
   for(Int_t b=1;b<iNbinsEta;b++)
   {
    dv2proEtaPOI = fHistProDiffFlowEtaPOI->GetBinContent(b);
    dErrdifcombEtaPOI = fHistProDiffFlowEtaPOI->GetBinError(b);//to be improved (treatment of errors for non-Gaussian distribution needed!)
    //fill common hist results:
    fCommonHistsRes->FillDifferentialFlowEtaPOI(b, dv2proEtaPOI, dErrdifcombEtaPOI); 
   }
  }   
  
  cout<<endl;     	      	  
  //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	#define AliFlowAnalysisWithMCEventPlane_cxx
	17
	18	#include "Riostream.h" //needed as include
	19	#include "TFile.h" //needed as include
	20	#include "TProfile.h" //needed as include
9cc9e6cc	21	#include "TProfile2D.h"
f1d945a1	22	#include "TComplex.h" //needed as include
28ca24ad	23	#include "TList.h"
f1d945a1	24
	25	class TH1F;
	26
	27	#include "AliFlowCommonConstants.h" //needed as include
	28	#include "AliFlowEventSimple.h"
	29	#include "AliFlowTrackSimple.h"
	30	#include "AliFlowCommonHist.h"
	31	#include "AliFlowCommonHistResults.h"
	32	#include "AliFlowAnalysisWithMCEventPlane.h"
	33
	34	class AliFlowVector;
	35
	36	// AliFlowAnalysisWithMCEventPlane:
	37	// Description: Maker to analyze Flow from the generated MC reaction plane.
	38	// This class is used to get the real value of the flow
	39	// to compare the other methods to when analysing simulated events
	40	// author: N. van der Kolk (kolk@nikhef.nl)
	41
	42	ClassImp(AliFlowAnalysisWithMCEventPlane)
	43
	44	//-----------------------------------------------------------------------
	45
	46	AliFlowAnalysisWithMCEventPlane::AliFlowAnalysisWithMCEventPlane():
0092f3c2	47	fQsum(NULL),
f1d945a1	48	fQ2sum(0),
f1d945a1	49	fEventNumber(0),
f1d945a1	50	fDebug(kFALSE),
28ca24ad	51	fHistList(NULL),
8232a5ec	52	fCommonHists(NULL),
8232a5ec	53	fCommonHistsRes(NULL),
45cdbef2	54	fHistRP(NULL),
3963f25c	55	fHistProIntFlow(NULL),
9cc9e6cc	56	fHistProDiffFlowPtEtaRP(NULL),
45cdbef2	57	fHistProDiffFlowPtRP(NULL),
45cdbef2	58	fHistProDiffFlowEtaRP(NULL),
9cc9e6cc	59	fHistProDiffFlowPtEtaPOI(NULL),
45cdbef2	60	fHistProDiffFlowPtPOI(NULL),
45cdbef2	61	fHistProDiffFlowEtaPOI(NULL)
f1d945a1	62	{
	63
	64	// Constructor.
28ca24ad	65	fHistList = new TList();
28ca24ad	66
0092f3c2	67	fQsum = new TVector2; // flow vector sum
f1d945a1	68	}
	69
	70
	71	//-----------------------------------------------------------------------
	72
	73
	74	AliFlowAnalysisWithMCEventPlane::~AliFlowAnalysisWithMCEventPlane()
	75	{
	76	//destructor
28ca24ad	77	delete fHistList;
0092f3c2	78	delete fQsum;
f1d945a1	79	}
f1d945a1	80
7ebf0358	81	//-----------------------------------------------------------------------
	82
	83	void AliFlowAnalysisWithMCEventPlane::WriteHistograms(TString* outputFileName)
	84	{
	85	//store the final results in output .root file
	86	TFile *output = new TFile(outputFileName->Data(),"RECREATE");
7a2c2652	87	output->WriteObject(fHistList, "cobjMCEP","SingleKey");
7ebf0358	88	delete output;
7ebf0358	89	}
f1d945a1	90
f1d945a1	91	//-----------------------------------------------------------------------
b0fda271	92
	93	void AliFlowAnalysisWithMCEventPlane::WriteHistograms(TString outputFileName)
	94	{
	95	//store the final results in output .root file
	96	TFile *output = new TFile(outputFileName.Data(),"RECREATE");
	97	output->WriteObject(fHistList, "cobjMCEP","SingleKey");
	98	delete output;
	99	}
	100
	101	//-----------------------------------------------------------------------
f1d945a1	102	void AliFlowAnalysisWithMCEventPlane::Init() {
	103
	104	//Define all histograms
	105	cout<<"---Analysis with the real MC Event Plane---"<<endl;
	106
00731146	107	Int_t iNbinsPt = AliFlowCommonConstants::GetNbinsPt();
45cdbef2	108	Double_t dPtMin = AliFlowCommonConstants::GetPtMin();
	109	Double_t dPtMax = AliFlowCommonConstants::GetPtMax();
	110
	111	Int_t iNbinsEta = AliFlowCommonConstants::GetNbinsEta();
	112	Double_t dEtaMin = AliFlowCommonConstants::GetEtaMin();
	113	Double_t dEtaMax = AliFlowCommonConstants::GetEtaMax();
f1d945a1	114
9d062fe3	115	fCommonHists = new AliFlowCommonHist("AliFlowCommonHistMCEP");
	116	fHistList->Add(fCommonHists);
	117	fCommonHistsRes = new AliFlowCommonHistResults("AliFlowCommonHistResultsMCEP");
	118	fHistList->Add(fCommonHistsRes);
	119
9d062fe3	120	fHistRP = new TH1F("Flow_RP_MCEP","Flow_RP_MCEP",100,0.,3.14);
f1d945a1	121	fHistRP->SetXTitle("Reaction Plane Angle");
f1d945a1	122	fHistRP->SetYTitle("Counts");
28ca24ad	123	fHistList->Add(fHistRP);
45cdbef2	124
b7cb54d5	125	fHistProIntFlow = new TProfile("FlowPro_V_MCEP","FlowPro_V_MCEP",1,0.,1.);
3963f25c	126	fHistProIntFlow->SetLabelSize(0.06);
	127	(fHistProIntFlow->GetXaxis())->SetBinLabel(1,"v_{n}{2}");
	128	fHistProIntFlow->SetYTitle("");
	129	fHistList->Add(fHistProIntFlow);
9cc9e6cc	130
	131	fHistProDiffFlowPtEtaRP = new TProfile2D("FlowPro_VPtEtaRP_MCEP","FlowPro_VPtEtaRP_MCEP",iNbinsPt,dPtMin,dPtMax,iNbinsEta,dEtaMin,dEtaMax);
	132	fHistProDiffFlowPtEtaRP->SetXTitle("P_{t}");
	133	fHistProDiffFlowPtEtaRP->SetYTitle("#eta");
	134	fHistList->Add(fHistProDiffFlowPtEtaRP);
45cdbef2	135
b7cb54d5	136	fHistProDiffFlowPtRP = new TProfile("FlowPro_VPtRP_MCEP","FlowPro_VPtRP_MCEP",iNbinsPt,dPtMin,dPtMax);
45cdbef2	137	fHistProDiffFlowPtRP->SetXTitle("P_{t}");
	138	fHistProDiffFlowPtRP->SetYTitle("");
	139	fHistList->Add(fHistProDiffFlowPtRP);
	140
b7cb54d5	141	fHistProDiffFlowEtaRP = new TProfile("FlowPro_VetaRP_MCEP","FlowPro_VetaRP_MCEP",iNbinsEta,dEtaMin,dEtaMax);
45cdbef2	142	fHistProDiffFlowEtaRP->SetXTitle("#eta");
	143	fHistProDiffFlowEtaRP->SetYTitle("");
	144	fHistList->Add(fHistProDiffFlowEtaRP);
	145
9cc9e6cc	146	fHistProDiffFlowPtEtaPOI = new TProfile2D("FlowPro_VPtEtaPOI_MCEP","FlowPro_VPtEtaPOI_MCEP",iNbinsPt,dPtMin,dPtMax,iNbinsEta,dEtaMin,dEtaMax);
	147	fHistProDiffFlowPtEtaPOI->SetXTitle("P_{t}");
	148	fHistProDiffFlowPtEtaPOI->SetYTitle("#eta");
	149	fHistList->Add(fHistProDiffFlowPtEtaPOI);
	150
b7cb54d5	151	fHistProDiffFlowPtPOI = new TProfile("FlowPro_VPtPOI_MCEP","FlowPro_VPtPOI_MCEP",iNbinsPt,dPtMin,dPtMax);
45cdbef2	152	fHistProDiffFlowPtPOI->SetXTitle("P_{t}");
	153	fHistProDiffFlowPtPOI->SetYTitle("");
	154	fHistList->Add(fHistProDiffFlowPtPOI);
	155
b7cb54d5	156	fHistProDiffFlowEtaPOI = new TProfile("FlowPro_VetaPOI_MCEP","FlowPro_VetaPOI_MCEP",iNbinsEta,dEtaMin,dEtaMax);
45cdbef2	157	fHistProDiffFlowEtaPOI->SetXTitle("#eta");
	158	fHistProDiffFlowEtaPOI->SetYTitle("");
	159	fHistList->Add(fHistProDiffFlowEtaPOI);
f1d945a1	160
	161	fEventNumber = 0; //set number of events to zero
	162
	163	}
	164
	165	//-----------------------------------------------------------------------
	166
d7671632	167	void AliFlowAnalysisWithMCEventPlane::Make(AliFlowEventSimple* anEvent) {
f1d945a1	168
f1d945a1	169	//Calculate v2 from the MC reaction plane
8232a5ec	170	if (anEvent) {
d7671632	171
	172	// get the MC reaction plane angle
	173	Double_t aRP = anEvent->GetMCReactionPlaneAngle();
f1d945a1	174	//fill control histograms
8232a5ec	175	fCommonHists->FillControlHistograms(anEvent);
f1d945a1	176
f1d945a1	177	//get the Q vector from the FlowEvent
00731146	178	AliFlowVector vQ = anEvent->GetQ();
00731146	179	//cout<<"vQ.Mod() = " << vQ.Mod() << endl;
f1d945a1	180	//for chi calculation:
00731146	181	*fQsum += vQ;
f1d945a1	182	//cout<<"fQsum.Mod() = "<<fQsum.Mod()<<endl;
00731146	183	fQ2sum += vQ.Mod2();
9d062fe3	184	//cout<<"fQ2sum = "<<fQ2sum<<endl;
f1d945a1	185
00731146	186	fHistRP->Fill(aRP);
45cdbef2	187
	188	Double_t dPhi = 0.;
	189	Double_t dv2 = 0.;
	190	Double_t dPt = 0.;
	191	Double_t dEta = 0.;
	192	//Double_t dPi = TMath::Pi();
	193
f1d945a1	194	//calculate flow
f1d945a1	195	//loop over the tracks of the event
00731146	196	Int_t iNumberOfTracks = anEvent->NumberOfTracks();
00731146	197	for (Int_t i=0;i<iNumberOfTracks;i++)
f1d945a1	198	{
00731146	199	AliFlowTrackSimple* pTrack = anEvent->GetTrack(i) ;
00731146	200	if (pTrack){
b7cb54d5	201	if (pTrack->InRPSelection()){
45cdbef2	202	dPhi = pTrack->Phi();
b7cb54d5	203	dv2 = TMath::Cos(2*(dPhi-aRP));
b7cb54d5	204	dPt = pTrack->Pt();
45cdbef2	205	dEta = pTrack->Eta();
3963f25c	206	//no-name int. flow (to be improved):
3963f25c	207	fHistProIntFlow->Fill(0.,dv2);
9cc9e6cc	208	//differential flow (Pt, Eta, RP):
9cc9e6cc	209	fHistProDiffFlowPtEtaRP->Fill(dPt,dEta,dv2,1.);
45cdbef2	210	//differential flow (Pt, RP):
	211	fHistProDiffFlowPtRP->Fill(dPt,dv2,1.);
	212	//differential flow (Eta, RP):
	213	fHistProDiffFlowEtaRP->Fill(dEta,dv2,1.);
	214	}
b7cb54d5	215	if (pTrack->InPOISelection()) {
45cdbef2	216	dPhi = pTrack->Phi();
00731146	217	//if (dPhi<0.) dPhi+=2*TMath::Pi();
f1d945a1	218	//calculate flow v2:
b7cb54d5	219	dv2 = TMath::Cos(2*(dPhi-aRP));
b7cb54d5	220	dPt = pTrack->Pt();
45cdbef2	221	dEta = pTrack->Eta();
9cc9e6cc	222	//differential flow (Pt, Eta, POI):
9cc9e6cc	223	fHistProDiffFlowPtEtaPOI->Fill(dPt,dEta,dv2,1.);
45cdbef2	224	//differential flow (Pt, POI):
	225	fHistProDiffFlowPtPOI->Fill(dPt,dv2,1.);
	226	//differential flow (Eta, POI):
	227	fHistProDiffFlowEtaPOI->Fill(dEta,dv2,1.);
	228	}
f1d945a1	229	}//track selected
	230	}//loop over tracks
	231
	232	fEventNumber++;
deebd72f	233	// cout<<"@@@@@ "<<fEventNumber<<" events processed"<<endl;
f1d945a1	234	}
	235	}
	236
	237	//--------------------------------------------------------------------
	238	void AliFlowAnalysisWithMCEventPlane::Finish() {
	239
	240	//*************make histograms etc.
	241	if (fDebug) cout<<"AliFlowAnalysisWithMCEventPlane::Terminate()"<<endl;
9d062fe3	242
45cdbef2	243	Int_t iNbinsPt = AliFlowCommonConstants::GetNbinsPt();
	244	Int_t iNbinsEta = AliFlowCommonConstants::GetNbinsEta();
	245
3963f25c	246	// no-name int. flow (to be improved):
	247	Double_t dV = fHistProIntFlow->GetBinContent(1);
	248	Double_t dErrV = fHistProIntFlow->GetBinError(1); // to be improved (treatment of errors for non-Gaussian distribution needed!)
	249	// fill no-name int. flow (to be improved):
	250	fCommonHistsRes->FillIntegratedFlow(dV,dErrV);
	251	cout<<"dV{MC} is "<<dV<<" +- "<<dErrV<<endl;
65af4339	252
3963f25c	253	//RP:
b7cb54d5	254	TH1F* fHistPtRP = fCommonHists->GetHistPtRP();
3963f25c	255	Double_t dYieldPtRP = 0.;
	256	Double_t dVRP = 0.;
	257	Double_t dErrVRP = 0.;
	258	Double_t dSumRP = 0.;
	259	//differential flow (RP, Pt):
45cdbef2	260	Double_t dvPtRP = 0.;
45cdbef2	261	Double_t dErrvPtRP = 0.;
3963f25c	262	for(Int_t b=1;b<iNbinsPt;b++)
45cdbef2	263	{
3963f25c	264	dvPtRP = fHistProDiffFlowPtRP->GetBinContent(b);
	265	dErrvPtRP = fHistProDiffFlowPtRP->GetBinError(b);//to be improved (treatment of errors for non-Gaussian distribution needed!)
	266	fCommonHistsRes->FillDifferentialFlowPtRP(b, dvPtRP, dErrvPtRP);
	267	if(fHistPtRP){
	268	//integrated flow (RP)
	269	dYieldPtRP = fHistPtRP->GetBinContent(b);
	270	dVRP += dvPtRP*dYieldPtRP;
	271	dSumRP += dYieldPtRP;
	272	//error on integrated flow
	273	dErrVRP += dYieldPtRPdYieldPtRPdErrvPtRP*dErrvPtRP;
	274	}
45cdbef2	275	}
3963f25c	276	if (dSumRP != 0. ) {
	277	dVRP /= dSumRP; //because pt distribution should be normalised
	278	dErrVRP /= (dSumRP*dSumRP);
	279	dErrVRP = TMath::Sqrt(dErrVRP);
	280	}
	281	// fill integrated flow (RP):
	282	fCommonHistsRes->FillIntegratedFlowRP(dVRP,dErrVRP);
	283	cout<<"dV{MC} (RP) is "<<dVRP<<" +- "<<dErrVRP<<endl;
65af4339	284
3963f25c	285	//differential flow (RP, Eta):
45cdbef2	286	Double_t dvEtaRP = 0.;
45cdbef2	287	Double_t dErrvEtaRP = 0.;
3963f25c	288	for(Int_t b=1;b<iNbinsEta;b++)
45cdbef2	289	{
3963f25c	290	dvEtaRP = fHistProDiffFlowEtaRP->GetBinContent(b);
	291	dErrvEtaRP = fHistProDiffFlowEtaRP->GetBinError(b);//to be improved (treatment of errors for non-Gaussian distribution needed!)
	292	fCommonHistsRes->FillDifferentialFlowEtaRP(b, dvEtaRP, dErrvEtaRP);
45cdbef2	293	}
	294
	295	//POI:
b7cb54d5	296	TH1F* fHistPtPOI = fCommonHists->GetHistPtPOI();
3963f25c	297	Double_t dYieldPtPOI = 0.;
	298	Double_t dVPOI = 0.;
	299	Double_t dErrVPOI = 0.;
	300	Double_t dSumPOI = 0.;
	301	Double_t dv2proPtPOI = 0.;
	302	Double_t dErrdifcombPtPOI = 0.;
	303	Double_t dv2proEtaPOI = 0.;
	304	Double_t dErrdifcombEtaPOI = 0.;
45cdbef2	305	//Pt:
b7cb54d5	306	if(fHistProDiffFlowPtPOI) {
3963f25c	307	for(Int_t b=1;b<iNbinsPt;b++){
3963f25c	308	dv2proPtPOI = fHistProDiffFlowPtPOI->GetBinContent(b);
3963f25c	309	dErrdifcombPtPOI = fHistProDiffFlowPtPOI->GetBinError(b);//to be improved (treatment of errors for non-Gaussian distribution needed!)
f1d945a1	310	//fill TH1D
3963f25c	311	fCommonHistsRes->FillDifferentialFlowPtPOI(b, dv2proPtPOI, dErrdifcombPtPOI);
3963f25c	312	if (fHistPtPOI){
65af4339	313	//integrated flow (POI)
3963f25c	314	dYieldPtPOI = fHistPtPOI->GetBinContent(b);
	315	dVPOI += dv2proPtPOI*dYieldPtPOI;
	316	dSumPOI += dYieldPtPOI;
f1d945a1	317	//error on integrated flow
3963f25c	318	dErrVPOI += dYieldPtPOIdYieldPtPOIdErrdifcombPtPOI*dErrdifcombPtPOI;
f1d945a1	319	}
65af4339	320	}//end of for(Int_t b=0;b<iNbinsPt;b++)
45cdbef2	321	} else { cout<<"fHistProFlow is NULL"<<endl; }
3963f25c	322	if (dSumPOI != 0. ) {
	323	dVPOI /= dSumPOI; //because pt distribution should be normalised
	324	dErrVPOI /= (dSumPOI*dSumPOI);
	325	dErrVPOI = TMath::Sqrt(dErrVPOI);
65af4339	326	}
3963f25c	327	cout<<"dV{MC} (POI) is "<<dVPOI<<" +- "<<dErrVPOI<<endl;
	328
	329	fCommonHistsRes->FillIntegratedFlowPOI(dVPOI,dErrVPOI);
45cdbef2	330
45cdbef2	331	//Eta:
	332	if(fHistProDiffFlowEtaPOI)
	333	{
3963f25c	334	for(Int_t b=1;b<iNbinsEta;b++)
45cdbef2	335	{
3963f25c	336	dv2proEtaPOI = fHistProDiffFlowEtaPOI->GetBinContent(b);
3963f25c	337	dErrdifcombEtaPOI = fHistProDiffFlowEtaPOI->GetBinError(b);//to be improved (treatment of errors for non-Gaussian distribution needed!)
45cdbef2	338	//fill common hist results:
3963f25c	339	fCommonHistsRes->FillDifferentialFlowEtaPOI(b, dv2proEtaPOI, dErrdifcombEtaPOI);
45cdbef2	340	}
65af4339	341	}
3963f25c	342
3963f25c	343	cout<<endl;
b7cb54d5	344	//cout<<".....finished"<<endl;
9d062fe3	345	}
f1d945a1	346
f1d945a1	347
9d062fe3	348