StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
bichselG.C
1 /*
2  root.exe lBichsel.C bichselG.C+
3  bichselG("N"); // dN/dx
4  bichselG("I70"; // I70
5  bichselG("Bz"); // Ifit
6  TH1D *pB70 = (TH1D *) ((TF1 *) gROOT->GetListOfFunctions()->FindObject("B70p1"))->GetHistogram();
7  TH1D *piB70 = (TH1D *) ((TF1 *) gROOT->GetListOfFunctions()->FindObject("B70#pi1"))->GetHistogram();
8  TH1D *diffB70 = new TH1D(*pB70);
9  diffB70->SetName("diffB70");
10  diffB70->Add(pB70,piB70,1,-1);
11  diffB70->SetLineColor(1);
12  TH1D *pBz = (TH1D *) ((TF1 *) gROOT->GetListOfFunctions()->FindObject("Bzp1"))->GetHistogram();
13  TH1D *piBz = (TH1D *) ((TF1 *) gROOT->GetListOfFunctions()->FindObject("Bz#pi1"))->GetHistogram();
14  TH1D *diffBz = new TH1D(*pBz);
15  diffBz->SetName("diffBz");
16  diffBz->Add(pBz,piBz,1,-1);
17  diffBz->SetLineColor(2);
18  TH1D *pdNdx = (TH1D *) ((TF1 *) gROOT->GetListOfFunctions()->FindObject("dNdxp1"))->GetHistogram();
19  TH1D *pidNdx = (TH1D *) ((TF1 *) gROOT->GetListOfFunctions()->FindObject("dNdx#pi1"))->GetHistogram();
20  TH1D *diffdNdx = new TH1D(*pdNdx);
21  diffdNdx->SetName("diffdNdx");
22  diffdNdx->Add(pdNdx,pidNdx,1,-1);
23  diffdNdx->SetLineColor(3);
24  cppi = new TCanvas("cppi","cppi");
25  diffB70->SetXTitle("log_{10}p")
26  diffB70->SetTitle("z_{p} - z_{#pi}");
27  diffB70->Draw("l");
28  diffBz->Draw("samel");
29  diffdNdx->Draw("samel");
30  TLegend *l = new TLegend(0.6,0.6,0.8,0.8);
31  l->AddEntry(diff70,"I70");
32  l->AddEntry(diffBz,"Ifit");
33  l->AddEntry(diffdNdx,"dNdx");
34 
35  TH1D *eB70 = (TH1D *) ((TF1 *) gROOT->GetListOfFunctions()->FindObject("B70e1"))->GetHistogram();
36  TH1D *piB70 = (TH1D *) ((TF1 *) gROOT->GetListOfFunctions()->FindObject("B70#pi1"))->GetHistogram();
37  TH1D *diffB70 = new TH1D(*eB70);
38  diffB70->SetName("diffB70");
39  diffB70->Add(eB70,piB70,1,-1);
40  diffB70->SetLineColor(1);
41  TH1D *eBz = (TH1D *) ((TF1 *) gROOT->GetListOfFunctions()->FindObject("Bze1"))->GetHistogram();
42  TH1D *piBz = (TH1D *) ((TF1 *) gROOT->GetListOfFunctions()->FindObject("Bz#pi1"))->GetHistogram();
43  TH1D *diffBz = new TH1D(*eBz);
44  diffBz->SetName("diffBz");
45  diffBz->Add(eBz,piBz,1,-1);
46  diffBz->SetLineColor(2);
47  TH1D *edNdx = (TH1D *) ((TF1 *) gROOT->GetListOfFunctions()->FindObject("dNdxe1"))->GetHistogram();
48  TH1D *pidNdx = (TH1D *) ((TF1 *) gROOT->GetListOfFunctions()->FindObject("dNdx#pi1"))->GetHistogram();
49  TH1D *diffdNdx = new TH1D(*edNdx);
50  diffdNdx->SetName("diffdNdx");
51  diffdNdx->Add(edNdx,pidNdx,1,-1);
52  diffdNdx->SetLineColor(3);
53  cepi = new TCanvas("cepi","cepi");
54  diffB70->SetXTitle("log_{10}p")
55  diffB70->SetTitle("z_{e} - z_{#pi}");
56  diffB70->Draw("l");
57  diffBz->Draw("samel");
58  diffdNdx->Draw("samel");
59  TLegend *l = new TLegend(0.6,0.6,0.8,0.8);
60  l->AddEntry(diff70,"I70");
61  l->AddEntry(diffBz,"Ifit");
62  l->AddEntry(diffdNdx,"dNdx");
63 
64 */
65 #if !defined(__CINT__)
66 // code that should be seen ONLY by the compiler
67 #else
68 #if !defined(__CINT__) || defined(__MAKECINT__)
69 // code that should be seen by the compiler AND rootcint
70 #else
71 // code that should always be seen
72 #endif
73 #endif
74 //________________________________________________________________________________
75 #if !defined(__CINT__) || defined(__MAKECINT__)
76 #include "Riostream.h"
77 #include <stdio.h>
78 #include "TF1.h"
79 #include "TMath.h"
80 #include "TSystem.h"
81 #include "TCanvas.h"
82 #include "TClassTable.h"
83 #include "StBichsel/Bichsel.h"
84 #include "StBichsel/StdEdxModel.h"
85 #include "TLegend.h"
86 #include "TROOT.h"
87 #else
88 class Bichsel;
89 #endif
90 Bichsel *m_Bichsel = 0;
91 const Int_t NMasses = 11;
92 const Double_t Masses[NMasses] = {0.93827231,
93  0.493677,
94  0.13956995,
95  0.51099907e-3,
96  1.87561339,
97  0.1056584,
98  2.80925,
99  2.80923, //GEANT3
100  3.727417, //GEANT3
101  0.13956995,
102  0.93827231*6.94/1.008
103 };
104 const Int_t Index[NMasses] = { 4, 3, 2, 0, 5, 1, 6, 7, 8, -2, 9};
105 const Char_t *Names[NMasses] = {"p", "K","#pi","e", "d","#mu","t","He3","#alpha","2#pi", "Li"};
106 const Int_t NF = 4; // 0 1 2 3 4 5 6 7 8 -2 9
107 const Char_t *FNames[8] = {"Girrf","Sirrf","Bz","B70","B60","B70M","dNdx","BzM"};
108 const Int_t Nlog2dx = 3;
109 const Double_t log2dx[Nlog2dx] = {0,1,2};
110 //________________________________________________________________________________
111 Double_t bichselZ(Double_t *x,Double_t *par) {
112  Double_t pove = x[0];
113  Double_t scale = 1;
114  Double_t mass = par[0];
115  if (mass < 0) {mass = - mass; scale = 2;}
116  Double_t poverm = pove/mass;
117  Double_t charge = 1.;
118  Double_t dx2 = 1;
119  if (par[1] > 1.0) {
120  charge = par[1];
121  poverm *= charge;
122  dx2 = TMath::Log2(5.);
123  }
124  return (scale*charge*charge*TMath::Exp(m_Bichsel->GetMostProbableZ(TMath::Log10(poverm),dx2)));//TMath::Exp(7.81779499999999961e-01));
125  //return charge*charge*TMath::Log10(m_Bichsel->GetI70(TMath::Log10(poverm),1.));
126 }
127 //________________________________________________________________________________
128 Double_t bichselZM(Double_t *x,Double_t *par) {
129  Double_t pove = x[0];
130  Double_t scale = 1;
131  Double_t mass = par[0];
132  if (mass < 0) {mass = - mass; scale = 2;}
133  Double_t poverm = pove/mass;
134  Double_t charge = 1.;
135  Double_t dx2 = 1;
136  if (par[1] > 1.0) {
137  charge = par[1];
138  poverm *= charge;
139  dx2 = TMath::Log2(5.);
140  }
141  return 1e6*(scale*charge*charge*TMath::Exp(m_Bichsel->GetMostProbableZM(TMath::Log10(poverm),dx2)));//TMath::Exp(7.81779499999999961e-01));
142  //return charge*charge*TMath::Log10(m_Bichsel->GetI70(TMath::Log10(poverm),1.));
143 }
144 //________________________________________________________________________________
145 Double_t bichsel70(Double_t *x,Double_t *par) {
146  Double_t pove = x[0];
147  Double_t scale = 1;
148  Double_t mass = par[0];
149  if (mass < 0) {mass = - mass; scale = 2;}
150  Double_t poverm = pove/mass;
151  Double_t charge = 1.;
152  Double_t dx2 = 1;
153  if (par[1] > 1.0) {
154  charge = par[1];
155  poverm *= charge;
156  dx2 = TMath::Log2(5.);
157  }
158  // return 1e6*(scale*charge*charge*m_Bichsel->GetI70M(TMath::Log10(poverm),dx2));//TMath::Exp(7.81779499999999961e-01));
159  return charge*charge*(m_Bichsel->GetI70(TMath::Log10(poverm),1.));
160 }
161 //________________________________________________________________________________
162 Double_t bichsel70M(Double_t *x,Double_t *par) {
163  Double_t pove = x[0];
164  Double_t scale = 1;
165  Double_t mass = par[0];
166  if (mass < 0) {mass = - mass; scale = 2;}
167  Double_t poverm = pove/mass;
168  Double_t charge = 1.;
169  Double_t dx2 = 1;
170  if (par[1] > 1.0) {
171  charge = par[1];
172  poverm *= charge;
173  dx2 = TMath::Log2(5.);
174  }
175  return (scale*charge*charge*m_Bichsel->GetI70M(TMath::Log10(poverm),dx2));//TMath::Exp(7.81779499999999961e-01));
176 }
177 //________________________________________________________________________________
178 Double_t dNdx(Double_t *x,Double_t *par) {
179  Double_t pove = x[0];
180  Double_t scale = 1;
181  Double_t mass = par[0];
182  if (mass < 0) {mass = - mass; scale = 2;}
183  Double_t poverm = pove/mass;
184  Double_t charge = 1.;
185  Double_t dx2 = 1;
186  if (par[1] > 1.0) charge = par[1];
187  poverm *= charge;
188  return (scale*StdEdxModel::instance()->dNdx(poverm,charge));//TMath::Exp(7.81779499999999961e-01));
189 }
190 #if !defined(__CINT__) && !defined(__CLING__)
191 //________________________________________________________________________________
192 Double_t aleph70P(Double_t *x,Double_t *par) {
193  /*
194  W.Blum, L. Rolandi "Particle Detection with Drift Chambers", page 246, eq. (9.5)
195  F_g(v) = p[0]/beta**p[3]*(p[1] - beta**p[3] - log(p[2] + (beta*gamma)**-p[4]);
196  F_g(v) = p[0]*(1/beta**p[3]*(p[1] - log(p[2] + 1/(beta*gamma)**p[4])) - 1)
197  */
198  Double_t bg = x[0];
199  Double_t b2inv = 1. + 1./(bg*bg);
200  Double_t beta = 1./TMath::Sqrt(b2inv);
201  Double_t dEdx = par[0]*(-1 + TMath::Power(beta,-par[3])*(par[1] - TMath::Log(TMath::Max(1e-10,par[2] + TMath::Power(bg,-par[4])))));
202  return dEdx;
203 };
204 //________________________________________________________________________________
205 Double_t aleph70(Double_t *x,Double_t *par) {
206  static const Double_t dEdxMIP = 2.39761562607903311;
207  static Double_t MIPBetaGamma = 4.;
208 #if 0
209  struct Par_t {Int_t h, N; Double_t xmin, xmax, pars[10];};
210  const Par_t Par[9] = {
211  /* name h n+1 xmin xmax pars[n+1] */
212  /* electron */{ 0, 4, 3.0, 6.0,{ 0.14105, -0.09078, 0.01901, -0.00128, 0, 0, 0, 0, 0, 0}},
213  /* muon */{ 1, 2, 0.0, 4.5,{ -0.00689, 0.00256, 0, 0, 0, 0, 0, 0, 0, 0}},
214  /* pion */{ 2, -1, 0.0, 4.5,{ 0.00000, 0, 0, 0, 0, 0, 0, 0, 0, 0}},
215  /* kaon */{ 3, 9, -0.1, 3.7,{ 0.00869, 0.21918, -0.88919, 1.30023, -0.97075, 0.41298, -0.10214, 0.01379, -0.00079, 0}},
216  /* proton */{ 4, 9, -0.6, 3.3,{ 0.03052, -0.02423, -0.05636, -0.11585, 0.41292, -0.38956, 0.16837, -0.03494, 0.00283, 0}},
217  /* deuteron */{ 5, 8, -1.0, 2.9,{ 0.03523, -0.10625, 0.04182, 0.07820, -0.03816, -0.02735, 0.01940, -0.00304, 0, 0}},
218  /* triton */{ 6, 10, -1.0, 2.8,{ 0.03092, -0.07846, 0.01668, 0.00331, 0.12771, -0.05480, -0.13897, 0.13928, -0.04715, 0.00555}},
219  /* He3 */{ 7, 9, -0.8, 2.9,{ 0.09158, -0.07816, 0.07039, 0.00578, -0.16160, 0.26547, -0.18728, 0.05988, -0.00710, 0}},
220  /* alpha */{ 8, 10, -0.8, 2.9,{ 0.09366, -0.08276, 0.06191, 0.02631, -0.17044, 0.30536, -0.26867, 0.11847, -0.02505, 0.00201}}
221  };
222  // const Double_t ppar[7] = { 0.0857988, 9.46138, 0.000206655, 2.12222, 0.974, -1, 0.13957}; /* pion */
223  static Double_t ppar[7] = { 0.0857988, 9.46138, 0.000206655, 2.12222, 0.974, -1, 0.13957}; /* pion */
224 #else
225  // static Double_t ppar[7] = { 0.0762, 10.632, 0.134e-4, 1.863, 1.948, -1, -1}; /* Aleph parameters from Alice TPC TDR */
226  // static Double_t ppar[7] = { 0.08942, 8.91971, 0.00024, 2.27383, 1.54174, -1.00000, 0}; /* pion */
227  static Double_t ppar[7] = {0.12337, 6.61371, 0.00201, 2.27381, 1.54174, -1.00000, 0 }; /* g All */
228 #endif
229  static Double_t Norm = dEdxMIP/aleph70P(&MIPBetaGamma,ppar);
230  Int_t hyp = (Int_t ) par[0];
231  Int_t h = Index[hyp];
232  Double_t ScaleL10 = 0;
233  if (h < 0) {
234  h = -h;
235  ScaleL10 = TMath::Log10(2.);
236  }
237  Double_t pove = x[0];
238  Double_t mass = Masses[hyp];
239  Double_t poverm = pove/mass;
240  Double_t charge = 1.;
241  if (h > 6 && h > 9) charge = 2;
242  else if (h == 10) charge = 3;
243  poverm *= charge;
244  Double_t bg = poverm;
245  /*
246  W.Blum, L. Rolandi "Particle Detection with Drift Chambers", page 246, eq. (9.5)
247  F_g(v) = p[0]/beta**p[3]*(p[1] - beta**p[3] - log(p[2] + (beta*gamma)**-p[4]));
248  F_g(v) = p[0]*(1/beta**p[3]*(p[1] - log(p[2] + 1/(beta*gamma)**p[4])) - 1)
249  */
250  Double_t dEdxL10 = TMath::Log10(Norm*aleph70P(&bg,ppar));
251 #if 0
252  if (Par[h].N > 0) {
253  TString fName(Form("pol%i",Par[h].N-1));
254  TF1 *f = (TF1 *) gROOT->GetListOfFunctions()->FindObject(fName);
255  if (! f) {
256  f = new TF1(fName,fName,0,1);
257  }
258  f->SetParameters(&Par[h].pars[0]);
259  Double_t bgL10 = TMath::Log10(bg);
260  dEdxL10 += f->Eval(bgL10);
261  }
262 #endif
263  return 2*TMath::Log10(charge) + dEdxL10 + ScaleL10;
264 }
265 #endif /* __CINT__ */
266 //________________________________________________________________________________
267 void bichselG(const Char_t *type="I70m") {
268  if (gClassTable->GetID("StBichsel") < 0 || !m_Bichsel) {
269  gSystem->Load("libTable");
270  gSystem->Load("St_base");
271  gSystem->Load("StarClassLibrary");
272  gSystem->Load("StBichsel");
273  m_Bichsel = Bichsel::Instance();
274  }
275  TString Type(type);
276  TLegend *leg = new TLegend(0.65,0.45,0.75,0.9,"");
277  Double_t xmax = 10;
278  for (int h = 0; h < NMasses; h++) { // Masses
279  // for (int h = 0; h < 7; h++) { // Masses
280  Int_t f = 3;
281  if (Type.Contains("BzM",TString::kIgnoreCase)) f = 7;
282  else if (Type.Contains("Bz",TString::kIgnoreCase)) f = 2;
283  else if (Type.Contains("I70M",TString::kIgnoreCase)) f = 5;
284  else if (Type.Contains("I70",TString::kIgnoreCase)) f = 3;
285  else if (Type.Contains("I60",TString::kIgnoreCase)) f = 4;
286  else if (Type.Contains("N",TString::kIgnoreCase)) f = 6;
287  Int_t dx = 1;
288  Char_t *FunName = Form("%s%s%i",FNames[f],Names[h],(int)log2dx[dx]);
289  cout << "Make " << h << "\t" << FunName << endl;
290  Double_t xmin = 0.1;
291  TF1 *func = 0;
292  if (f == 3) func = new TF1(FunName,bichsel70,xmin, xmax,2);
293  else if (f == 2) func = new TF1(FunName,bichselZ ,xmin, xmax,2);
294  else if (f == 5) func = new TF1(FunName,bichsel70M ,xmin, xmax,2);
295  else if (f == 6) func = new TF1(FunName,dNdx ,xmin, xmax,2);
296  else if (f == 7) func = new TF1(FunName,bichselZM,xmin, xmax,2);
297  else {
298  return;
299  }
300  if (h == 9) func->SetParameter(0,-Masses[h]);
301  else func->SetParameter(0,Masses[h]);
302  func->SetParameter(1,1.);
303  if (h >= 7 && h < 9) func->SetParameter(1,2.);
304  if (h == 10) func->SetParameter(1,3.);
305  Int_t color = h+1;
306  if (color > 8) color -= 8;
307  // if (color > 7) color++;
308 #if 1
309  func->SetLineColor(color);
310  func->SetMarkerColor(color);
311 #endif
312  func->Draw("same");
313  leg->AddEntry(func,Names[h]);
314 #if !defined( __CINT__) && defined(__Aleph__)
315  TF1 *fA = new TF1(Form("Aleph%s",Names[h]),aleph70,xmin,xmax, 1);
316  fA->SetParameter(0,h);
317  fA->SetLineColor(color);
318  fA->SetMarkerColor(color);
319  fA->SetLineStyle(2);
320  fA->Draw("same");
321  leg->AddEntry(fA,Names[h]);
322 #endif
323  }
324  leg->Draw();
325 }
Generated by   doxygen 1.8.5