StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
bichsel.C
1 #include "Riostream.h"
2 class Bichsel;
3 Bichsel *m_Bichsel = 0;
4 const Int_t NMasses = 5;
5 const Double_t Masses[NMasses] = {0.93827231,
6  0.493677,
7  0.13956995,
8  0.51099907e-3,
9  1.87561339};
10 const Char_t *Names[NMasses] = {"p", "K","pi","e","d"};
11 const Int_t NF = 4;
12 const Char_t *FNames[5] = {"Girrf","Sirrf","Bz","B70","B60"};
13 const Int_t Nlog2dx = 1;//3;
14 const Double_t log2dx[Nlog2dx] = {1};//0,1,2};
15 //________________________________________________________________________________
16 Double_t gfunc(Double_t *x,Double_t *par) {
17  Double_t pove = x[0];
18  Double_t poverm = pove/par[0];
19  Int_t type = par[1];
20  Int_t k = 0;
21  if (type == 3) k = 1;
22  Int_t type = par[1];
23  Int_t k = 0;
24  if (type == 3) k = 1;
25  // Double_t Scale = BetheBloch::Sirrf(4.,par[2],k)/
26  // BetheBloch::Girrf(4.,par[4],k);
27  return BetheBloch::Girrf(poverm,par[4],k);
28 }
29 //________________________________________________________________________________
30 Double_t sifunc(Double_t *x,Double_t *par) {
31  Double_t pove = x[0];
32  Double_t poverm = pove/par[0];
33  Int_t type = par[1];
34  Int_t k = 0;
35  if (type == 3) k = 1;
36  return BetheBloch::Sirrf(poverm,par[2],k);
37 }
38 // Double_t bbfunc(Double_t *x,Double_t *par) {
39 // Double_t pove = x[0];
40 // Double_t poverm = pove/par[0];
41 // BetheBloch BB;
42 // Double_t value = 1.e6*BB(poverm);
43 // // printf("x : %f p: %f val : %f \n",x[0],poverm,value);
44 // return value;
45 // }
46 Double_t bichselZ(Double_t *x,Double_t *par) {
47  Double_t pove = x[0];
48  Double_t poverm = pove/par[0];
49  Int_t type = par[1];
50  Int_t k = 0;
51  if (type == 3) k = 1;
52  // Double_t Scale = BetheBloch::Sirrf(4.,par[2],k)/
53  // TMath::Exp(Bichsel::Instance()->GetMostProbableZ(TMath::Log10(4),par[3]));
54  return TMath::Exp(Bichsel::Instance()->GetMostProbableZ(TMath::Log10(poverm),par[3]));
55 }
56 Double_t bichsel70(Double_t *x,Double_t *par) {
57  Double_t pove = x[0];
58  Double_t poverm = pove/par[0];
59  Int_t type = par[1];
60  Int_t k = 0;
61  if (type == 3) k = 1;
62  // Double_t Scale = BetheBloch::Sirrf(4.,par[2],k)/
63  // TMath::Exp(Bichsel::Instance()->GetI70(TMath::Log10(4),par[3]));
64  return Bichsel::Instance()->GetI70(TMath::Log10(poverm),par[3]);
65 }
66 Double_t bichsel60(Double_t *x,Double_t *par) {
67  Double_t pove = x[0];
68  Double_t poverm = pove/par[0];
69  Int_t type = par[1];
70  Int_t k = 0;
71  if (type == 3) k = 1;
72  // Double_t Scale = BetheBloch::Sirrf(4.,par[2],k)/
73  // TMath::Exp(Bichsel::Instance()->GetI60(TMath::Log10(4),par[3]));
74  return Bichsel::Instance()->GetI60(TMath::Log10(poverm),par[3]);
75 }
76 void bichsel() {
77  if (gClassTable->GetID("StBichsel") < 0) {
78  gSystem->Load("libTable");
79  gSystem->Load("St_base");
80  gSystem->Load("StarClassLibrary");
81  gSystem->Load("StBichsel");
82  }
83  if (!m_Bichsel) m_Bichsel = Bichsel::Instance();
84  TCanvas *c1 = new TCanvas("c1");
85  c1->SetLogx();
86  c1->SetLogy();
87  c1->SetGrid();
88  // TH1F *hr = c1->DrawFrame(2.e-2,1,1.e3,1.e2);
89  // TH1F *hr = c1->DrawFrame(1.e-2,1,1.e3,1.e2);
90  TH1F *hr = c1->DrawFrame(1.e-1,1,1.e2,2.e2);
91  // hr->SetXTitle("Momentum (GeV/c)");
92  hr->SetTitle("dE/dx predictions");
93  hr->SetXTitle("Momentum (GeV/c) ");
94  hr->SetYTitle("dE/dx (keV/cm)");
95  hr->Draw();
96  // Mass Type Length log2(dx)
97  Double_t params[5] = { 1.0, 0., 60., 1., 1e-3};
98  TLegend *leg = new TLegend(0.4,0.7,0.9,0.9,"");//TLegend(0.79,0.91,0.89,0.89,"");
99  for (Int_t h = 0; h < 4; h++) { // Masses
100  params[0] = Masses[h];
101  params[1] = h;
102  for (Int_t f = 0; f < NF; f++) { // Functions
103  for (Int_t dx = 0; dx < Nlog2dx; dx++) {
104  params[3] = log2dx[dx];
105  Char_t *FunName = Form("%s%s%i",FNames[f],Names[h],(int)log2dx[dx]);
106  // cout << "Make " << FunName << endl;
107  TF1 *func = 0;
108  if (TString(FNames[f]) == "Sirrf") {
109  func = new TF1(FunName,sifunc,1.e-2,1.e3,5);
110  func->SetLineColor(2);
111  }
112  if (TString(FNames[f]) == "Girrf") {
113  func = new TF1(FunName,gfunc,1.e-2,1.e3,5);
114  params[4] = 1.e-3;
115  if (dx == 1) params[4] = 1.e-2;
116  if (dx == 2) params[4] = 1.e-3;
117  if (dx == 3) params[4] = 1.e-4;
118  func->SetLineColor(3);
119  }
120  else {if (TString(FNames[f]) == "Bz") {
121  func = new TF1(FunName,bichselZ,1.e-2,1.e3,5);
122  func->SetLineColor(4);
123  }
124  else {if (TString(FNames[f]) == "B70") {
125  func = new TF1(FunName,bichsel70,1.e-2,1.e3,5);
126  func->SetLineColor(6);
127  }
128  else {if (TString(FNames[f]) == "B60") {
129  func = new TF1(FunName,bichsel60,1.e-2,1.e3,5);
130  func->SetLineColor(7);
131  }}}}
132  if (! func) continue;
133  func->SetParameters(params);
134  func->Draw("same");
135  if (dx == 0 && h == 0) {
136  TString name(FNames[f]);
137  if (name == "Sirrf") name += ": Fitted from Year 1 data";
138  else if (name == "Girrf") name += ": Geant3 prediction";
139  else if (name == "Bz") name += ": Bichsel most probable (dX = 2 cm)";
140  else if (name == "B70") name += ": Bichsel, 30 % truncation (dX = 2 cm)";
141  else if (name == "B60") name += ": Bichsel, 40 % truncation (dX = 2 cm)";
142  leg->AddEntry(func,name.Data(),"L");
143  }
144  }
145  }
146  }
147  leg->Draw();
148 }