levy_func_problem.C
Go to the documentation of this file.
1 //
2 // testing problems with Levy function used at AGKY/KNO
3 //
4 // INPUTS:
5 // c -> Levy function parameter
6 // method -> 0 : compute P(n) from <n>P(n) = KNO(n/<n>)
7 // 1 : compute P(n) from a Poisson transformation of the
8 // asymptotic scaling function (KNO/Levy) as in
9 // hep-ph/9608479 v1 29-Aug-1996, Eq.5
10 // -> 2 : like 1 but psi in Eq.5 is not the Levy function
11 // but a Gamma function
12 //
13 // Costas Andreopoulos <costas.andreopoulos \at stfc.ac.uk>
14 // University of Liverpool & STFC Rutherford Appleton Lab
15 //
16 
17 double integrand (double * x, double * par);
18 double prob_method0 (double n, double nav, double c);
19 double prob_method1 (double n, double nav, double c);
20 double prob_method2 (double n, double nav, double c);
21 double kno_func (double z, double c);
22 void plot_multiplicity_prob (double c);
23 
24 //.................................................................
25 void levy_func_problem(double c, int method)
26 {
27  cout << "Using method = " << method << ", Levy(c) = " << c << endl;
28 
29  const int nnav = 40;
30 
31  double nav_min = 0.25;
32  double nav_max = 3.00;
33  double dnav = (nav_max-nav_min)/(nnav-1);
34 
35  // canvas for plotting <n>_out = f(<n>_in)
36  TCanvas * canvas = new TCanvas("canvas","",20,20,700,700);
37  canvas->SetFillColor(0);
38  canvas->SetBorderMode(0);
39  TH1F * frame = (TH1F*) canvas->DrawFrame(nav_min,nav_min,nav_max,nav_max);
40  frame->Draw();
41 
42  double nav_in [nnav];
43  double nav_out[nnav];
44 
45  for(int i=0; i<nnav; i++) {
46 
47  double nav_input = nav_min + i*dnav;
48 
49  cout << "** nav(in) = " << nav_input << endl;
50 
51  double minmult = 0;
52  double maxmult = 10;
53  int nbins = TMath::Nint(maxmult-minmult+1);
54 
55  TH1D * mult_prob = new TH1D("mult_prob",
56  "hadronic multiplicity distribution", nbins, minmult-0.5, maxmult+0.5);
57  mult_prob->SetDirectory(0);
58 
59  int nbins = mult_prob->FindBin(maxmult);
60  for(int j=1; j<=nbins; j++) {
61  double n = mult_prob->GetBinCenter(j); // bin centre
62 
63  double P = 0;
64  if (method==0) P = prob_method0(n,nav_input,c);
65  else if (method==1) P = prob_method1(n,nav_input,c);
66  else if (method==2) P = prob_method2(n,nav_input,c);
67 
68  cout << " n = " << n << " -> P(n) = " << P << endl;
69  mult_prob->Fill(n,P);
70  }
71  double nav_output = mult_prob->GetMean();
72  delete mult_prob;
73 
74  cout << " ** nav(out) = " << nav_output << endl;
75 
76  nav_in[i] = nav_input;
77  nav_out[i] = nav_output;
78  }
79 
80  TGraph * gr = new TGraph(nnav,nav_in,nav_out);
81  gr->SetMarkerStyle(20);
82  gr->SetMarkerSize(1);
83 
84  TGraph * id = new TGraph(nnav,nav_in,nav_in);
85  id->SetLineStyle(1);
86  id->SetLineColor(2);
87 
88  gr->Draw("P");
89  id->Draw("L");
90 
91  canvas->Update();
92 }
93 //.................................................................
95 {
96  const int nnav = 2;
97  double nav[nnav] = { 0.5, 4. };
98  TH1D * mult_prob[nnav][2];
99 
100  double minmult = 0;
101  double maxmult = 10;
102  int nbins = TMath::Nint(maxmult-minmult+1);
103 
104  // canvas for plotting sample multiplicity probability distributions
105  TCanvas * canvas = new TCanvas("canvas","",120,120,800,800);
106  canvas->SetFillColor(0);
107  canvas->SetBorderMode(0);
108 
109  for(int i=0; i<nnav; i++) {
110 
111  double nav_input = nav[i];
112 
113  cout << "** nav = " << nav_input << endl;
114 
115  for(int method=0; method<=1; method++) {
116  mult_prob[i][method] = new TH1D("mult_prob",
117  "hadronic multiplicity distribution", nbins, minmult-0.5, maxmult+0.5);
118  mult_prob[i][method]->SetDirectory(0);
119 
120  int nbins = mult_prob[i][method]->FindBin(maxmult);
121  for(int j=1; j<=nbins; j++) {
122  double n = mult_prob[i][method]->GetBinCenter(j); // bin centre
123 
124  double P = 0;
125  if (method==0) P = prob_method0(n,nav_input,c);
126  else if (method==1) P = prob_method1(n,nav_input,c);
127 
128  cout << " n = " << n << " -> P(n) = " << P << endl;
129  mult_prob[i][method]->Fill(n,P);
130  }//bins
131  mult_prob[i][method]->Scale(1/mult_prob[i][method]->Integral("width"));
132  }//meth
133  }//nnav
134 
135  for(int i=0; i<nnav; i++) {
136  for(int method=0; method<=1; method++) {
137  mult_prob[i][method]->SetLineStyle(1);
138  mult_prob[i][method]->SetLineWidth(1);
139  mult_prob[i][method]->SetLineColor(method+1);
140  }
141  }
142 
143  canvas->cd();
144  canvas->Divide(1,nnav);
145 
146  for(int i=0; i<nnav; i++) {
147  canvas->cd(1+i);
148  mult_prob[i][0]->Draw();
149  mult_prob[i][1]->Draw("SAME");
150  }
151  canvas->Update();
152 }
153 //.................................................................
154 double prob_method2(double n, double nav, double c)
155 {
156 // Compute P(n) as in hep-ph/9608479 v1 29-Aug-1996, eq.5
157 //
158  TF1 * intg = new TF1("integrand",integrand,0,10,3);
159 
160  intg->SetParameter(0,n);
161  intg->SetParameter(1,nav);
162  intg->SetParameter(2,-999999);
163 
164  double P = intg->Integral(0,10);
165 
166  delete intg;
167 
168  return P;
169 }
170 //.................................................................
171 double prob_method1(double n, double nav, double c)
172 {
173 // Compute P(n) as in hep-ph/9608479 v1 29-Aug-1996, eq.5
174 //
175  TF1 * intg = new TF1("integrand",integrand,0,10,3);
176 
177  intg->SetParameter(0,n);
178  intg->SetParameter(1,nav);
179  intg->SetParameter(2,c);
180 
181  double P = intg->Integral(0,10);
182 
183  delete intg;
184 
185  return P;
186 }
187 //.................................................................
188 double prob_method0(double n, double nav, double c)
189 {
190 // compute P(n) from <n>P(n) = KNO(n/<n>)
191 //
192  double z = n/nav;
193  double kno = kno_func(z,c);
194  double P = kno/nav;
195  return P;
196 }
197 //.................................................................
198 double integrand(double * x, double * par)
199 {
200 // preasymptotic multiplicity distribution reconstructed from the
201 // asymptotic scaling function via Poisson transform
202 
203  // inputs
204  double xx = x[0];
205 
206  // parameters
207  double n = par[0];
208  double nav = par[1];
209  double c = par[2];
210 
211  // compute integrand
212 
213  double psi = 0;
214  if(c<0) {
215  if(n==0) psi = 0;
216  else psi = TMath::Gamma(xx);
217  }
218  else psi = kno_func(xx,c);
219 
220  double fct = TMath::Factorial((int)n);
221  double pwr = TMath::Power(nav*xx,n);
222  double expn = TMath::Exp(-nav*xx);
223 
224  double f = psi * (pwr/fct) * expn;
225  return f;
226 }
227 //.................................................................
228 double kno_func(double z, double c)
229 {
230 // z reduced multiplicity:
231 // c KNO param
232 
233  double x = c*z+1;
234  double psi = 2*TMath::Exp(-c)*TMath::Power(c,x)/TMath::Gamma(x);
235 
236  return psi;
237 }
238 //.................................................................
Double_t xx
Definition: macro.C:12
void levy_func_problem(double c, int method)
double prob_method2(double n, double nav, double c)
double Integral(const Spectrum &s, const double pot, int cvnregion)
Int_t par
Definition: SimpleIterate.C:24
void plot_multiplicity_prob(double c)
float Gamma() const
double kno_func(double z, double c)
#define P(a, b, c, d, e, x)
const int nbins
Definition: cellShifts.C:15
double prob_method1(double n, double nav, double c)
const double j
Definition: BetheBloch.cxx:29
z
Definition: test.py:28
OStream cout
Definition: OStream.cxx:6
double integrand(double *x, double *par)
double prob_method0(double n, double nav, double c)