PWGCF/EBYE/BalanceFunctions/AliBalance.h

   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
  15 #include <vector>
  16 #include <TObject.h>
  17 #include "TString.h"
  18
  19 using std::vector;
  20
  21 #define ANALYSIS_TYPES  7
  22 #define MAXIMUM_NUMBER_OF_STEPS 1024
  23
  24 class TGraphErrors;
  25 class TH1D;
  26 class TH2D;
  27
  28 const TString kBFAnalysisType[ANALYSIS_TYPES] = {"y","eta","qlong","qout","qside","qinv","phi"};
  29
  30 class AliBalance : public TObject {
  31  public:
  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
  42   AliBalance();
  43   AliBalance(const AliBalance& balance);
  44   ~AliBalance();
  45
  46   void SetCentralityIdentifier(const char* centralityId) {
  47     fCentralityId = centralityId;}
  48
  49   void SetAnalysisLevel(const char* analysisLevel) {
  50     fAnalysisLevel = analysisLevel;}
  51   void SetShuffle(Bool_t shuffle) {fShuffle = shuffle;}
  52   void SetHBTcut(Bool_t HBTcut) {fHBTcut = HBTcut;}
  53   void SetConversionCut(Bool_t ConversionCut) {fConversionCut = ConversionCut;}
  54   void SetInterval(Int_t iAnalysisType, Double_t p1Start, Double_t p1Stop,
  55                    Int_t ibins, Double_t p2Start, Double_t p2Stop);
  56   void SetCentralityInterval(Double_t cStart, Double_t cStop)  { fCentStart = cStart; fCentStop = cStop;};
  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();}
  67   Int_t GetNumberOfAnalyzedEvent()  const {return fAnalyzedEvents;}
  68
  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];}
  74
  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 {
  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]; }
  85
  86   void CalculateBalance(Float_t fCentrality, vector<Double_t> **chargeVector,Float_t bSign = 0.);
  87
  88   Double_t GetBalance(Int_t a, Int_t p2);
  89   Double_t GetError(Int_t a, Int_t p2);
  90
  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];}
  97
  98   TH2D *GetQAHistHBTbefore()         {return fHistHBTbefore;};
  99   TH2D *GetQAHistHBTafter()          {return fHistHBTafter;};
 100   TH2D *GetQAHistConversionbefore()  {return fHistConversionbefore;};
 101   TH2D *GetQAHistConversionafter()   {return fHistConversionafter;};
 102
 103   void PrintAnalysisSettings();
 104   TGraphErrors *DrawBalance(Int_t fAnalysisType);
 105
 106   void SetHistNp(Int_t iAnalysisType, TH2D *gHist) {
 107     fHistP[iAnalysisType] = gHist;}
 108   void SetHistNn(Int_t iAnalysisType, TH2D *gHist) {
 109     fHistN[iAnalysisType] = gHist;}
 110   void SetHistNpn(Int_t iAnalysisType, TH2D *gHist) {
 111     fHistPN[iAnalysisType] = gHist;}
 112   void SetHistNnp(Int_t iAnalysisType, TH2D *gHist) {
 113     fHistNP[iAnalysisType] = gHist;}
 114   void SetHistNpp(Int_t iAnalysisType, TH2D *gHist) {
 115     fHistPP[iAnalysisType] = gHist;}
 116   void SetHistNnn(Int_t iAnalysisType, TH2D *gHist) {
 117     fHistNN[iAnalysisType] = gHist;}
 118
 119   TH1D *GetBalanceFunctionHistogram(Int_t iAnalysisType,Double_t centrMin, Double_t centrMax, Double_t etaWindow = -1, Bool_t correctWithEfficiency = kFALSE, Bool_t correctWithAcceptanceOnly = kFALSE);
 120   void PrintResults(Int_t iAnalysisType, TH1D *gHist);
 121
 122  private:
 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);
 124
 125   Bool_t fShuffle; // shuffled balance function object
 126   Bool_t fHBTcut;  // apply HBT like cuts
 127   Bool_t fConversionCut;  // apply conversion cuts
 128
 129   TString fAnalysisLevel; //ESD, AOD or MC
 130   Int_t fAnalyzedEvents; //number of events that have been analyzed
 131
 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
 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
 142
 143
 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(+-)
 147   Double_t fNnp[ANALYSIS_TYPES][MAXIMUM_NUMBER_OF_STEPS]; //N(-+)
 148   Double_t fNn[ANALYSIS_TYPES], fNp[ANALYSIS_TYPES]; //number of pos./neg. inside the analyzed interval
 149
 150   Double_t fB[ANALYSIS_TYPES][MAXIMUM_NUMBER_OF_STEPS]; //BF matrix
 151   Double_t ferror[ANALYSIS_TYPES][MAXIMUM_NUMBER_OF_STEPS]; //error of the BF
 152
 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--
 159
 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
 166
 167
 168   AliBalance & operator=(const AliBalance & ) {return *this;}
 169
 170   ClassDef(AliBalance, 3)
 171 };
 172
 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
 199 #endif