StRoot  1
StNegativeBinomial.cxx
1 /******************************************************************************
2  * \$Id: StNegativeBinomial.cxx,v 1.2 2012/04/25 05:13:24 hmasui Exp \$
3  * \$Log: StNegativeBinomial.cxx,v \$
4  * Revision 1.2 2012/04/25 05:13:24 hmasui
5  * Use STAR logger. Take into account additional constant inefficiency by using the trigger bias parameter
6  *
7 ******************************************************************************/
8
9 #include <assert.h>
10
11 #include "TH1.h"
12 #include "TMath.h"
13 #include "TString.h"
14
15 #include "StMessMgr.h"
16
17 #include "StGlauberUtilities/StGlauberUtilities.h"
18 #include "StNegativeBinomial.h"
19
20 ClassImp(StNegativeBinomial)
21
22  UInt_t StNegativeBinomial::mCounter = 0 ;
23
24 //____________________________________________________________________________________________________
25 StNegativeBinomial::StNegativeBinomial(const Double_t npp, const Double_t k, const Double_t x,
26  const Double_t efficiency, const Double_t triggerbias, const Bool_t isConstEfficiency)
27  : mEfficiency(efficiency), mTriggerBias(triggerbias), mNpp(npp), mK(k), mX(x), mIsConstEfficiency(isConstEfficiency)
28 {
29  LOG_INFO << Form("StNegativeBinomial (npp, k, x, eff, trig. bias) = (%1.3f, %1.3f, %1.3f, %1.3f, %1.3f)",
30  mNpp, mK, mX, mEfficiency, mTriggerBias)
31  << endm;
32
33  if(mIsConstEfficiency){
34  LOG_INFO << "StNegativeBinomial Use constant efficiency" << endm;
35  }
36  else{
37  LOG_INFO << "StNegativeBinomial Use multiplicity dependent efficiency" << endm;
38  }
39
41  mhNbd = 0 ;
42  InitHistogram() ;
43 }
44
45 //____________________________________________________________________________________________________
47 {
48  delete mhNbd ;
49 }
50
51 //____________________________________________________________________________________________________
52 void StNegativeBinomial::InitHistogram()
53 {
54  const Int_t nbin = 100 ;
55
56  if( mhNbd ){
57  // Reset current stuffs
58  mhNbd->Reset();
59  }
60  else{
61  mhNbd = new TH1D(Form("mhNbd_%d", mCounter++), "", nbin, 0, nbin);
62  }
63
64  for(Int_t i=0;i<nbin;i++){
65  mhNbd->SetBinContent(i+1, GetNegativeBinomial(i));
66  }
67 }
68
69 //____________________________________________________________________________________________________
70 void StNegativeBinomial::SetParameters(const Double_t npp, const Double_t k, const Double_t x)
71 {
72  mNpp = npp ;
73  mK = k ;
74
75  if( x > 0.0 ){
76  mX = x ;
77  }
78
80  InitHistogram() ;
81 }
82
83 //____________________________________________________________________________________________________
84 Double_t StNegativeBinomial::GetTwoComponentMultiplicity(const Double_t npart, const Double_t ncoll) const
85 {
87
88  return (mX==0.0) ? npart : 0.5*(1.0-mX)*npart + mX*ncoll ;
89 }
90
91 //____________________________________________________________________________________________________
92 Int_t StNegativeBinomial::GetMultiplicity(const Double_t npart, const Double_t ncoll) const
93 {
95  // Take into account trigger efficiency, multiplicity efficiency
96
97  const Double_t nchPP = GetTwoComponentMultiplicity(npart, ncoll) ;
98  const Double_t nchSampled = nchPP ;
99
100  // Sample multiplicity (including trigger bias)
101  const Int_t nchInt = TMath::Nint(nchSampled);
102  Int_t multIdeal = 0;
103  for(Int_t ich=0; ich<nchInt; ich++){
104  multIdeal += (Int_t)mhNbd->GetRandom();
105  }
106
107  // Multiplicity dependent efficiency for TPC
108  const Double_t efficiency = (mIsConstEfficiency) ? mEfficiency : GetEfficiency(multIdeal) ;
109  //cout<<efficiency<<" ;"<<npart<<" ;"<<ncoll<<" TEST "<<endl;
110
111  Int_t mult = 0;
112  for(Int_t im=0; im<2*multIdeal; im++){
113  if( StGlauberUtilities::instance()->GetUniform2() <= efficiency ){
114  mult++;
115  }
116  }
117
118  if ( mTriggerBias == 1.0 ) return mult ;
119
120  const Int_t nmult = mult ;
121  mult = 0 ;
122  for(Int_t im=0; im<nmult; im++){
123  if( StGlauberUtilities::instance()->GetUniform() <= mTriggerBias ){
124  mult++;
125  }
126  }
127
128  return mult ;
129 }
130
131 //____________________________________________________________________________________________________
132 TH1* StNegativeBinomial::GetMultiplicity(const Double_t npart, const Double_t ncoll,
133  const Double_t weight) const
134 {
135  // Do not forget to delete histogram
136  const Double_t nchPP = GetTwoComponentMultiplicity(npart, ncoll) ;
137  const Int_t nch = TMath::Nint(nchPP * mTriggerBias) ; // with trigger bias
138
139  // Multiplicity dependent efficiency for TPC
140 // const Double_t efficiency = (mUseTpc) ? GetEfficiency(nch) : mEfficiency ;
141 // const Double_t efficiency = mEfficiency ;
142  const Double_t efficiency = (mIsConstEfficiency) ? mEfficiency : GetEfficiency(nch) ;
143  const Int_t nchSampled = TMath::Nint(nch * efficiency) ;
144
145  const Int_t nbin = 1000 ;
146  TH1* h = new TH1D("hmultTmp", "", nbin, 0, static_cast<Double_t>(nbin));
147
148  for(Int_t ix=0; ix<nbin; ix++){
149  // Probability to obtain 'nchSampled' in 'ix'-th bin
150  const Double_t probability = GetNegativeBinomial(ix, nchSampled) ;
151
152  // To avoid nan
153  if(probability>0.0 && probability<DBL_MAX){
154  h->Fill(ix+0.5, probability*weight);
155  }
156  }
157
158  return h ;
159 }
160
161 //____________________________________________________________________________________________________
162 Double_t StNegativeBinomial::GetNegativeBinomial(const Int_t n) const
163 {
164  const Double_t Const = TMath::Exp( TMath::LnGamma(n+mK) - TMath::LnGamma(n+1) - TMath::LnGamma(mK) );
165  const Double_t term = n * TMath::Log(mNpp/mK) - (n+mK) * TMath::Log(mNpp/mK+1.0) ;
166
167  return Const * TMath::Exp(term);
168 }
169
170 //____________________________________________________________________________________________________
171 Double_t StNegativeBinomial::GetNegativeBinomial(const Int_t n, const Double_t m) const
172 {
173  // npp and k are scaled by m
174
175  const Double_t Const = TMath::Exp( TMath::LnGamma(n+mK*m) - TMath::LnGamma(n+1) - TMath::LnGamma(mK*m) );
176  const Double_t term = n * TMath::Log(mNpp/mK) - (n+mK*m) * TMath::Log(mNpp/mK+1.0) ;
177
178  return Const * TMath::Exp(term);
179 }
180
181 //______________________________________________________________________________________________________
182 Double_t StNegativeBinomial::GetEfficiency(const Int_t mult) const
183 {
185  if( mIsConstEfficiency ){
186  Error("StNegativeBinomial::GetEfficiency", "Something is wrong, supposed to be constant efficiency. abort");
187  assert(0) ;
188  }
189
191  if ( mult < 0 ) return mEfficiency ;
192
193  // Multiplicity dependent efficiency correction (valid only for main TPC refmult)
194
195  // const Double_t eff0 = 1.05 ; // 98 % efficiency at p+p
196  const Double_t eff0 = 0.98 ; // 98 % efficiency at p+p
197 // const Double_t d = 0.14 ; // 14 % drop at most central Au + Au
198  const Double_t d = mEfficiency ;
199
200  return eff0 * (1.0 - mult * d/540.0) ;
201  //return eff0 * d * (1.0 - mult * 0.14/560.0) ;
202 }
203
204 //______________________________________________________________________________________________________
206 {
207  if(mhNbd) mhNbd->Draw() ;
208 }
209
void SetParameters(const Double_t npp, const Double_t k, const Double_t x=-1.0)
Get NBD(npp*m, k*m; n)
Definition: FJcore.h:367
void DrawNbd() const
Get flag for efficiency.
Double_t GetUniform() const
Get uniform distribution in 0 &lt; x &lt; 1.
Double_t GetUniform2() const
Double_t GetEfficiency() const
Get x parameter.
virtual ~StNegativeBinomial()
Default destructor.
Int_t GetMultiplicity(const Double_t npart, const Double_t ncoll) const
Get multiplcity by convoluting NBD.
Double_t GetTwoComponentMultiplicity(const Double_t npart, const Double_t ncoll) const
(1-x)*npart/2 + x*ncoll