[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 "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),
   fHistProFlow(NULL),
   fHistRP(NULL),
   fHistProIntFlow(NULL),
   fHistProDiffFlowPtRP(NULL),
   fHistProDiffFlowEtaRP(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);
  
  fHistProFlow = new TProfile("FlowPro_VPt_MCEP","FlowPro_VPt_MCEP",iNbinsPt,dPtMin,dPtMax);
  fHistProFlow->SetXTitle("P_{t}");
  fHistProFlow->SetYTitle("v_{2}");
  fHistList->Add(fHistProFlow);

  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("fHistProIntFlow","fHistProIntFlow",1,0.,1.);
  fHistProIntFlow->SetLabelSize(0.06);
  (fHistProIntFlow->GetXaxis())->SetBinLabel(1,"v_{n}{2}");
  fHistProIntFlow->SetYTitle("");
  fHistList->Add(fHistProIntFlow);
 
  fHistProDiffFlowPtRP = new TProfile("fHistProDiffFlowPtRP","fHistProDiffFlowPtRP",iNbinsPt,dPtMin,dPtMax);
  fHistProDiffFlowPtRP->SetXTitle("P_{t}");
  fHistProDiffFlowPtRP->SetYTitle("");
  fHistList->Add(fHistProDiffFlowPtRP);  
  
  fHistProDiffFlowEtaRP = new TProfile("fHistProDiffFlowEtaRP","fHistProDiffFlowEtaRP",iNbinsEta,dEtaMin,dEtaMax);
  fHistProDiffFlowEtaRP->SetXTitle("#eta");
  fHistProDiffFlowEtaRP->SetYTitle("");
  fHistList->Add(fHistProDiffFlowEtaRP);
  
  fHistProDiffFlowPtPOI = new TProfile("fHistProDiffFlowPtPOI","fHistProDiffFlowPtPOI",iNbinsPt,dPtMin,dPtMax);
  fHistProDiffFlowPtPOI->SetXTitle("P_{t}");
  fHistProDiffFlowPtPOI->SetYTitle("");
  fHistList->Add(fHistProDiffFlowPtPOI);  
  
  fHistProDiffFlowEtaPOI = new TProfile("fHistProDiffFlowEtaPOI","fHistProDiffFlowEtaPOI",iNbinsEta,dEtaMin,dEtaMax);
  fHistProDiffFlowEtaPOI->SetXTitle("#eta");
  fHistProDiffFlowEtaPOI->SetYTitle("");
  fHistList->Add(fHistProDiffFlowEtaPOI);          
 
  fEventNumber = 0;  //set number of events to zero
      
} 
 
//-----------------------------------------------------------------------
 
void AliFlowAnalysisWithMCEventPlane::Make(AliFlowEventSimple* anEvent, Double_t aRP) {

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