[u/mrichter/AliRoot.git] / PWGCF / EBYE / BalanceFunctions / AliBalance.h

#ifndef ALIBALANCE_H
#define ALIBALANCE_H
/*  See cxx source for full Copyright notice */


/* $Id$ */

//-------------------------------------------------------------------------
//                          Class AliBalance
//   This is the class for the Balance Function analysis
//
//    Origin: Panos Christakoglou, UOA-CERN, Panos.Christakoglou@cern.ch
//-------------------------------------------------------------------------

#include <vector>
#include <TObject.h>
#include "TString.h"

using std::vector;

#define ANALYSIS_TYPES	7
#define MAXIMUM_NUMBER_OF_STEPS	1024

class TGraphErrors;
class TH1D;
class TH2D;

const TString kBFAnalysisType[ANALYSIS_TYPES] = {"y","eta","qlong","qout","qside","qinv","phi"};

class AliBalance : public TObject {
 public:
  enum EAnalysisType {
    kRapidity = 0, 
    kEta      = 1, 
    kQlong    = 2, 
    kQout     = 3, 
    kQside    = 4, 
    kQinv     = 5, 
    kPhi      = 6 
  };

  AliBalance();
  AliBalance(const AliBalance& balance);
  ~AliBalance();

  void SetCentralityIdentifier(const char* centralityId) {
    fCentralityId = centralityId;}
  
  void SetAnalysisLevel(const char* analysisLevel) {
    fAnalysisLevel = analysisLevel;}
  void SetShuffle(Bool_t shuffle) {fShuffle = shuffle;}
  void SetHBTcut(Bool_t HBTcut) {fHBTcut = HBTcut;}
  void SetConversionCut(Bool_t ConversionCut) {fConversionCut = ConversionCut;}
  void SetInterval(Int_t iAnalysisType, Double_t p1Start, Double_t p1Stop,
		   Int_t ibins, Double_t p2Start, Double_t p2Stop);
  void SetCentralityInterval(Double_t cStart, Double_t cStop)  { fCentStart = cStart; fCentStop = cStop;};
  void SetNp(Int_t analysisType, Double_t NpSet)   { fNp[analysisType] = NpSet; }
  void SetNn(Int_t analysisType, Double_t NnSet)   { fNn[analysisType] = NnSet; }
  void SetNpp(Int_t analysisType, Int_t ibin, Double_t NppSet) { if(ibin > -1 && ibin < MAXIMUM_NUMBER_OF_STEPS) fNpp[analysisType][ibin] = NppSet; }
  void SetNpn(Int_t analysisType, Int_t ibin, Double_t NpnSet) { if(ibin > -1 && ibin < MAXIMUM_NUMBER_OF_STEPS) fNpn[analysisType][ibin] = NpnSet; }
  void SetNnp(Int_t analysisType, Int_t ibin, Double_t NnpSet) { if(ibin > -1 && ibin < MAXIMUM_NUMBER_OF_STEPS) fNnp[analysisType][ibin] = NnpSet; }
  void SetNnn(Int_t analysisType, Int_t ibin, Double_t NnnSet) { if(ibin > -1 && ibin < MAXIMUM_NUMBER_OF_STEPS) fNnn[analysisType][ibin] = NnnSet; }

  void InitHistograms(void);

  const char* GetAnalysisLevel() {return fAnalysisLevel.Data();}
  Int_t GetNumberOfAnalyzedEvent()  const {return fAnalyzedEvents;}

  Int_t GetNumberOfBins(Int_t ibin) const {return fNumberOfBins[ibin];}
  Double_t GetP1Start(Int_t ibin)   const {return fP1Start[ibin];}
  Double_t GetP1Stop(Int_t ibin)    const {return fP1Stop[ibin];}   
  Double_t GetP2Start(Int_t ibin)   const {return fP2Start[ibin];}
  Double_t GetP2Stop(Int_t ibin)    const {return fP2Stop[ibin];}    
 
  Double_t GetNp(Int_t analysisType) const { return 1.0*fNp[analysisType]; }
  Double_t GetNn(Int_t analysisType) const { return 1.0*fNn[analysisType]; }
  Double_t GetNnn(Int_t analysisType, Int_t p2) const { 
    return 1.0*fNnn[analysisType][p2]; }
  Double_t GetNpp(Int_t analysisType, Int_t p2) const { 
    return 1.0*fNpp[analysisType][p2]; }
  Double_t GetNpn(Int_t analysisType, Int_t p2) const { 
    return 1.0*fNpn[analysisType][p2]; }  
  Double_t GetNnp(Int_t analysisType, Int_t p2) const { 
    return 1.0*fNnp[analysisType][p2]; }

  void CalculateBalance(Float_t fCentrality, vector<Double_t> **chargeVector,Float_t bSign = 0.);
  
  Double_t GetBalance(Int_t a, Int_t p2);
  Double_t GetError(Int_t a, Int_t p2);

  TH2D *GetHistNp(Int_t iAnalysisType)  const { return fHistP[iAnalysisType];}
  TH2D *GetHistNn(Int_t iAnalysisType)  const { return fHistN[iAnalysisType];}
  TH2D *GetHistNpn(Int_t iAnalysisType) const { return fHistPN[iAnalysisType];}
  TH2D *GetHistNnp(Int_t iAnalysisType) const { return fHistNP[iAnalysisType];}
  TH2D *GetHistNpp(Int_t iAnalysisType) const { return fHistPP[iAnalysisType];}
  TH2D *GetHistNnn(Int_t iAnalysisType) const { return fHistNN[iAnalysisType];}

  TH2D *GetQAHistHBTbefore()         {return fHistHBTbefore;};
  TH2D *GetQAHistHBTafter()          {return fHistHBTafter;};
  TH2D *GetQAHistConversionbefore()  {return fHistConversionbefore;};
  TH2D *GetQAHistConversionafter()   {return fHistConversionafter;};

  void PrintAnalysisSettings();
  TGraphErrors *DrawBalance(Int_t fAnalysisType);

  void SetHistNp(Int_t iAnalysisType, TH2D *gHist) { 
    fHistP[iAnalysisType] = gHist;}
  void SetHistNn(Int_t iAnalysisType, TH2D *gHist) { 
    fHistN[iAnalysisType] = gHist;}
  void SetHistNpn(Int_t iAnalysisType, TH2D *gHist) { 
    fHistPN[iAnalysisType] = gHist;}
  void SetHistNnp(Int_t iAnalysisType, TH2D *gHist) { 
    fHistNP[iAnalysisType] = gHist;}
  void SetHistNpp(Int_t iAnalysisType, TH2D *gHist) { 
    fHistPP[iAnalysisType] = gHist;}
  void SetHistNnn(Int_t iAnalysisType, TH2D *gHist) { 
    fHistNN[iAnalysisType] = gHist;}

  TH1D *GetBalanceFunctionHistogram(Int_t iAnalysisType,Double_t centrMin, Double_t centrMax, Double_t etaWindow = -1, Bool_t correctWithEfficiency = kFALSE, Bool_t correctWithAcceptanceOnly = kFALSE, Bool_t correctWithMixed = kFALSE, TH1D *hMixed[4]=NULL);
  void PrintResults(Int_t iAnalysisType, TH1D *gHist);

 private:
  inline Float_t GetDPhiStar(Float_t phi1, Float_t pt1, Float_t charge1, Float_t phi2, Float_t pt2, Float_t charge2, Float_t radius, Float_t bSign); 

  Bool_t fShuffle; // shuffled balance function object
  Bool_t fHBTcut;  // apply HBT like cuts
  Bool_t fConversionCut;  // apply conversion cuts

  TString fAnalysisLevel; //ESD, AOD or MC
  Int_t fAnalyzedEvents; //number of events that have been analyzed

  TString fCentralityId;//Centrality identifier to be used for the histo naming

  Int_t fNumberOfBins[ANALYSIS_TYPES];//number of bins of the analyzed interval
  Double_t fP1Start[ANALYSIS_TYPES];//lower boundaries for single particle histograms 
  Double_t fP1Stop[ANALYSIS_TYPES];//upper boundaries for single particle histograms 
  Double_t fP2Start[ANALYSIS_TYPES];//lower boundaries for pair histograms 
  Double_t fP2Stop[ANALYSIS_TYPES];//upper boundaries for pair histograms 
  Double_t fP2Step[ANALYSIS_TYPES];//bin size for pair histograms 
  Double_t fCentStart;//lower boundary for centrality
  Double_t fCentStop;//upper boundary for centrality

 	
  Double_t fNnn[ANALYSIS_TYPES][MAXIMUM_NUMBER_OF_STEPS]; //N(--)
  Double_t fNpp[ANALYSIS_TYPES][MAXIMUM_NUMBER_OF_STEPS]; //N(++)
  Double_t fNpn[ANALYSIS_TYPES][MAXIMUM_NUMBER_OF_STEPS]; //N(+-)
  Double_t fNnp[ANALYSIS_TYPES][MAXIMUM_NUMBER_OF_STEPS]; //N(-+)
  Double_t fNn[ANALYSIS_TYPES], fNp[ANALYSIS_TYPES]; //number of pos./neg. inside the analyzed interval
  
  Double_t fB[ANALYSIS_TYPES][MAXIMUM_NUMBER_OF_STEPS]; //BF matrix
  Double_t ferror[ANALYSIS_TYPES][MAXIMUM_NUMBER_OF_STEPS]; //error of the BF
  
  TH2D *fHistP[ANALYSIS_TYPES]; //N+
  TH2D *fHistN[ANALYSIS_TYPES]; //N-
  TH2D *fHistPN[ANALYSIS_TYPES]; //N+-
  TH2D *fHistNP[ANALYSIS_TYPES]; //N-+
  TH2D *fHistPP[ANALYSIS_TYPES]; //N++
  TH2D *fHistNN[ANALYSIS_TYPES]; //N--

  //QA histograms
  TH2D *fHistHBTbefore; // Delta Eta vs. Delta Phi before HBT inspired cuts
  TH2D *fHistHBTafter; // Delta Eta vs. Delta Phi after HBT inspired cuts
  TH2D *fHistConversionbefore; // Delta Eta vs. Delta Phi before Conversion cuts
  TH2D *fHistConversionafter; // Delta Eta vs. Delta Phi before Conversion cuts


  AliBalance & operator=(const AliBalance & ) {return *this;}

  ClassDef(AliBalance, 3)
};

Float_t AliBalance::GetDPhiStar(Float_t phi1, Float_t pt1, Float_t charge1, Float_t phi2, Float_t pt2, Float_t charge2, Float_t radius, Float_t bSign)
{ 
  //
  // calculates dphistar
  //
  
  Float_t dphistar = phi1 - phi2 - charge1 * bSign * TMath::ASin(0.075 * radius / pt1) + charge2 * bSign * TMath::ASin(0.075 * radius / pt2);
  
  static const Double_t kPi = TMath::Pi();
  
  // circularity
//   if (dphistar > 2 * kPi)
//     dphistar -= 2 * kPi;
//   if (dphistar < -2 * kPi)
//     dphistar += 2 * kPi;
  
  if (dphistar > kPi)
    dphistar = kPi * 2 - dphistar;
  if (dphistar < -kPi)
    dphistar = -kPi * 2 - dphistar;
  if (dphistar > kPi) // might look funny but is needed
    dphistar = kPi * 2 - dphistar;
  
  return dphistar;
}

#endif
Commit	Line	Data
6c178944	1	#ifndef ALIBALANCE_H
	2	#define ALIBALANCE_H
	3	/* See cxx source for full Copyright notice */
	4
	5
	6	/* $Id$ */
	7
	8	//-------------------------------------------------------------------------
	9	// Class AliBalance
	10	// This is the class for the Balance Function analysis
	11	//
	12	// Origin: Panos Christakoglou, UOA-CERN, Panos.Christakoglou@cern.ch
	13	//-------------------------------------------------------------------------
	14
c64cb1f6	15	#include <vector>
e559055a	16	#include <TObject.h>
9d1f0df5	17	#include "TString.h"
6c178944	18
c64cb1f6	19	using std::vector;
c64cb1f6	20
b966652f	21	#define ANALYSIS_TYPES 7
6c178944	22	#define MAXIMUM_NUMBER_OF_STEPS 1024
6c178944	23
6c178944	24	class TGraphErrors;
b966652f	25	class TH1D;
74b04dd4	26	class TH2D;
6c178944	27
9fd4b54e	28	const TString kBFAnalysisType[ANALYSIS_TYPES] = {"y","eta","qlong","qout","qside","qinv","phi"};
6432ac6a	29
5c33329d	30	class AliBalance : public TObject {
6c178944	31	public:
3feee083	32	enum EAnalysisType {
	33	kRapidity = 0,
	34	kEta = 1,
	35	kQlong = 2,
	36	kQout = 3,
	37	kQside = 4,
	38	kQinv = 5,
	39	kPhi = 6
	40	};
	41
e559055a	42	AliBalance();
7f0257ea	43	AliBalance(const AliBalance& balance);
6c178944	44	~AliBalance();
6a3f819e	45
	46	void SetCentralityIdentifier(const char* centralityId) {
	47	fCentralityId = centralityId;}
6c178944	48
6432ac6a	49	void SetAnalysisLevel(const char* analysisLevel) {
6432ac6a	50	fAnalysisLevel = analysisLevel;}
9fd4b54e	51	void SetShuffle(Bool_t shuffle) {fShuffle = shuffle;}
	52	void SetHBTcut(Bool_t HBTcut) {fHBTcut = HBTcut;}
	53	void SetConversionCut(Bool_t ConversionCut) {fConversionCut = ConversionCut;}
6432ac6a	54	void SetInterval(Int_t iAnalysisType, Double_t p1Start, Double_t p1Stop,
6432ac6a	55	Int_t ibins, Double_t p2Start, Double_t p2Stop);
74b04dd4	56	void SetCentralityInterval(Double_t cStart, Double_t cStop) { fCentStart = cStart; fCentStop = cStop;};
6432ac6a	57	void SetNp(Int_t analysisType, Double_t NpSet) { fNp[analysisType] = NpSet; }
	58	void SetNn(Int_t analysisType, Double_t NnSet) { fNn[analysisType] = NnSet; }
	59	void SetNpp(Int_t analysisType, Int_t ibin, Double_t NppSet) { if(ibin > -1 && ibin < MAXIMUM_NUMBER_OF_STEPS) fNpp[analysisType][ibin] = NppSet; }
	60	void SetNpn(Int_t analysisType, Int_t ibin, Double_t NpnSet) { if(ibin > -1 && ibin < MAXIMUM_NUMBER_OF_STEPS) fNpn[analysisType][ibin] = NpnSet; }
	61	void SetNnp(Int_t analysisType, Int_t ibin, Double_t NnpSet) { if(ibin > -1 && ibin < MAXIMUM_NUMBER_OF_STEPS) fNnp[analysisType][ibin] = NnpSet; }
	62	void SetNnn(Int_t analysisType, Int_t ibin, Double_t NnnSet) { if(ibin > -1 && ibin < MAXIMUM_NUMBER_OF_STEPS) fNnn[analysisType][ibin] = NnnSet; }
	63
	64	void InitHistograms(void);
	65
	66	const char* GetAnalysisLevel() {return fAnalysisLevel.Data();}
9fd4b54e	67	Int_t GetNumberOfAnalyzedEvent() const {return fAnalyzedEvents;}
6c178944	68
9fd4b54e	69	Int_t GetNumberOfBins(Int_t ibin) const {return fNumberOfBins[ibin];}
	70	Double_t GetP1Start(Int_t ibin) const {return fP1Start[ibin];}
	71	Double_t GetP1Stop(Int_t ibin) const {return fP1Stop[ibin];}
	72	Double_t GetP2Start(Int_t ibin) const {return fP2Start[ibin];}
	73	Double_t GetP2Stop(Int_t ibin) const {return fP2Stop[ibin];}
7f0257ea	74
b966652f	75	Double_t GetNp(Int_t analysisType) const { return 1.0*fNp[analysisType]; }
	76	Double_t GetNn(Int_t analysisType) const { return 1.0*fNn[analysisType]; }
	77	Double_t GetNnn(Int_t analysisType, Int_t p2) const {
	78	return 1.0*fNnn[analysisType][p2]; }
	79	Double_t GetNpp(Int_t analysisType, Int_t p2) const {
	80	return 1.0*fNpp[analysisType][p2]; }
	81	Double_t GetNpn(Int_t analysisType, Int_t p2) const {
6432ac6a	82	return 1.0*fNpn[analysisType][p2]; }
	83	Double_t GetNnp(Int_t analysisType, Int_t p2) const {
	84	return 1.0*fNnp[analysisType][p2]; }
09bb7bf4	85
845fdeca	86	void CalculateBalance(Float_t fCentrality, vector<Double_t> **chargeVector,Float_t bSign = 0.);
7f0257ea	87
3feee083	88	Double_t GetBalance(Int_t a, Int_t p2);
3feee083	89	Double_t GetError(Int_t a, Int_t p2);
6c178944	90
9fd4b54e	91	TH2D *GetHistNp(Int_t iAnalysisType) const { return fHistP[iAnalysisType];}
	92	TH2D *GetHistNn(Int_t iAnalysisType) const { return fHistN[iAnalysisType];}
	93	TH2D *GetHistNpn(Int_t iAnalysisType) const { return fHistPN[iAnalysisType];}
	94	TH2D *GetHistNnp(Int_t iAnalysisType) const { return fHistNP[iAnalysisType];}
	95	TH2D *GetHistNpp(Int_t iAnalysisType) const { return fHistPP[iAnalysisType];}
	96	TH2D *GetHistNnn(Int_t iAnalysisType) const { return fHistNN[iAnalysisType];}
b966652f	97
6de53600	98	TH2D *GetQAHistHBTbefore() {return fHistHBTbefore;};
	99	TH2D *GetQAHistHBTafter() {return fHistHBTafter;};
	100	TH2D *GetQAHistConversionbefore() {return fHistConversionbefore;};
	101	TH2D *GetQAHistConversionafter() {return fHistConversionafter;};
	102
09bb7bf4	103	void PrintAnalysisSettings();
33a86282	104	TGraphErrors *DrawBalance(Int_t fAnalysisType);
6432ac6a	105
74b04dd4	106	void SetHistNp(Int_t iAnalysisType, TH2D *gHist) {
6432ac6a	107	fHistP[iAnalysisType] = gHist;}
74b04dd4	108	void SetHistNn(Int_t iAnalysisType, TH2D *gHist) {
6432ac6a	109	fHistN[iAnalysisType] = gHist;}
74b04dd4	110	void SetHistNpn(Int_t iAnalysisType, TH2D *gHist) {
6432ac6a	111	fHistPN[iAnalysisType] = gHist;}
74b04dd4	112	void SetHistNnp(Int_t iAnalysisType, TH2D *gHist) {
6432ac6a	113	fHistNP[iAnalysisType] = gHist;}
74b04dd4	114	void SetHistNpp(Int_t iAnalysisType, TH2D *gHist) {
6432ac6a	115	fHistPP[iAnalysisType] = gHist;}
74b04dd4	116	void SetHistNnn(Int_t iAnalysisType, TH2D *gHist) {
6432ac6a	117	fHistNN[iAnalysisType] = gHist;}
6432ac6a	118
47e040b5	119	TH1D GetBalanceFunctionHistogram(Int_t iAnalysisType,Double_t centrMin, Double_t centrMax, Double_t etaWindow = -1, Bool_t correctWithEfficiency = kFALSE, Bool_t correctWithAcceptanceOnly = kFALSE, Bool_t correctWithMixed = kFALSE, TH1D hMixed[4]=NULL);
6432ac6a	120	void PrintResults(Int_t iAnalysisType, TH1D *gHist);
6c178944	121
6c178944	122	private:
9fd4b54e	123	inline Float_t GetDPhiStar(Float_t phi1, Float_t pt1, Float_t charge1, Float_t phi2, Float_t pt2, Float_t charge2, Float_t radius, Float_t bSign);
845fdeca	124
9fd4b54e	125	Bool_t fShuffle; // shuffled balance function object
	126	Bool_t fHBTcut; // apply HBT like cuts
	127	Bool_t fConversionCut; // apply conversion cuts
49a3750f	128
9d1f0df5	129	TString fAnalysisLevel; //ESD, AOD or MC
6c178944	130	Int_t fAnalyzedEvents; //number of events that have been analyzed
6432ac6a	131
6a3f819e	132	TString fCentralityId;//Centrality identifier to be used for the histo naming
	133
	134	Int_t fNumberOfBins[ANALYSIS_TYPES];//number of bins of the analyzed interval
9fd4b54e	135	Double_t fP1Start[ANALYSIS_TYPES];//lower boundaries for single particle histograms
	136	Double_t fP1Stop[ANALYSIS_TYPES];//upper boundaries for single particle histograms
	137	Double_t fP2Start[ANALYSIS_TYPES];//lower boundaries for pair histograms
	138	Double_t fP2Stop[ANALYSIS_TYPES];//upper boundaries for pair histograms
	139	Double_t fP2Step[ANALYSIS_TYPES];//bin size for pair histograms
	140	Double_t fCentStart;//lower boundary for centrality
	141	Double_t fCentStop;//upper boundary for centrality
74b04dd4	142
6c178944	143
b966652f	144	Double_t fNnn[ANALYSIS_TYPES][MAXIMUM_NUMBER_OF_STEPS]; //N(--)
	145	Double_t fNpp[ANALYSIS_TYPES][MAXIMUM_NUMBER_OF_STEPS]; //N(++)
	146	Double_t fNpn[ANALYSIS_TYPES][MAXIMUM_NUMBER_OF_STEPS]; //N(+-)
6432ac6a	147	Double_t fNnp[ANALYSIS_TYPES][MAXIMUM_NUMBER_OF_STEPS]; //N(-+)
b966652f	148	Double_t fNn[ANALYSIS_TYPES], fNp[ANALYSIS_TYPES]; //number of pos./neg. inside the analyzed interval
6c178944	149
b966652f	150	Double_t fB[ANALYSIS_TYPES][MAXIMUM_NUMBER_OF_STEPS]; //BF matrix
b966652f	151	Double_t ferror[ANALYSIS_TYPES][MAXIMUM_NUMBER_OF_STEPS]; //error of the BF
09bb7bf4	152
74b04dd4	153	TH2D *fHistP[ANALYSIS_TYPES]; //N+
	154	TH2D *fHistN[ANALYSIS_TYPES]; //N-
	155	TH2D *fHistPN[ANALYSIS_TYPES]; //N+-
	156	TH2D *fHistNP[ANALYSIS_TYPES]; //N-+
	157	TH2D *fHistPP[ANALYSIS_TYPES]; //N++
	158	TH2D *fHistNN[ANALYSIS_TYPES]; //N--
b966652f	159
6de53600	160	//QA histograms
	161	TH2D *fHistHBTbefore; // Delta Eta vs. Delta Phi before HBT inspired cuts
	162	TH2D *fHistHBTafter; // Delta Eta vs. Delta Phi after HBT inspired cuts
	163	TH2D *fHistConversionbefore; // Delta Eta vs. Delta Phi before Conversion cuts
	164	TH2D *fHistConversionafter; // Delta Eta vs. Delta Phi before Conversion cuts
	165
845fdeca	166
845fdeca	167
7f0257ea	168	AliBalance & operator=(const AliBalance & ) {return *this;}
7f0257ea	169
3feee083	170	ClassDef(AliBalance, 3)
7f0257ea	171	};
6c178944	172
845fdeca	173	Float_t AliBalance::GetDPhiStar(Float_t phi1, Float_t pt1, Float_t charge1, Float_t phi2, Float_t pt2, Float_t charge2, Float_t radius, Float_t bSign)
	174	{
	175	//
	176	// calculates dphistar
	177	//
	178
	179	Float_t dphistar = phi1 - phi2 - charge1 * bSign * TMath::ASin(0.075 * radius / pt1) + charge2 * bSign * TMath::ASin(0.075 * radius / pt2);
	180
	181	static const Double_t kPi = TMath::Pi();
	182
	183	// circularity
	184	// if (dphistar > 2 * kPi)
	185	// dphistar -= 2 * kPi;
	186	// if (dphistar < -2 * kPi)
	187	// dphistar += 2 * kPi;
	188
	189	if (dphistar > kPi)
	190	dphistar = kPi * 2 - dphistar;
	191	if (dphistar < -kPi)
	192	dphistar = -kPi * 2 - dphistar;
	193	if (dphistar > kPi) // might look funny but is needed
	194	dphistar = kPi * 2 - dphistar;
	195
	196	return dphistar;
	197	}
	198
6c178944	199	#endif