simple-tof-analysis
 All Classes Namespaces Functions Variables Groups Pages
BetheBloch.h
1 using namespace std;
2 
3 
4 
5 //__________________________________________________________________________________________________________________________
6 double BetheBlochAntoniMod(const double betaGamma)
7 {
8  const double lnbg = std::log(betaGamma);
9 
10  static const double table[] = {
11  6.4124, 6.0301, 5.6736, 5.3411, 5.0312, 4.7422, 4.4729, 4.2219,
12  3.9880, 3.7701, 3.5671, 3.3781, 3.2020, 3.0382, 2.8857, 2.7437,
13  2.6117, 2.4888, 2.3746, 2.2684, 2.1697, 2.0779, 1.9927, 1.9135,
14  1.8400, 1.7718, 1.7086, 1.6499, 1.5954, 1.5450, 1.4983, 1.4550,
15  1.4150, 1.3780, 1.3437, 1.3121, 1.2829, 1.2560, 1.2312, 1.2083,
16  1.1873, 1.1680, 1.1502, 1.1339, 1.1190, 1.1053, 1.0929, 1.0815,
17  1.0712, 1.0618, 1.0533, 1.0456, 1.0387, 1.0324, 1.0269, 1.0219,
18  1.0176, 1.0137, 1.0103, 1.0074, 1.0049, 1.0028, 1.0011, 1.0000,
19  1.0012, 1.0028, 1.0046, 1.0066, 1.0089, 1.0114, 1.0141, 1.0170,
20  1.0201, 1.0234, 1.0268, 1.0303, 1.0340, 1.0378, 1.0418, 1.0458,
21  1.0500, 1.0542, 1.0585, 1.0629, 1.0674, 1.0720, 1.0766, 1.0813,
22  1.0860, 1.0908, 1.0957, 1.1006, 1.1055, 1.1105, 1.1155, 1.1206,
23  1.1257, 1.1308, 1.1359, 1.1411, 1.1463, 1.1515, 1.1568, 1.1620,
24  1.1673, 1.1726, 1.1779, 1.1832, 1.1886, 1.1939, 1.1993, 1.2047,
25  1.2101, 1.2155, 1.2209, 1.2263, 1.2317, 1.2372, 1.2426, 1.2481,
26  1.2535, 1.2590, 1.2644, 1.2699, 1.2753, 1.2807, 1.2860, 1.2912,
27  1.2963, 1.3013, 1.3063, 1.3112, 1.3160, 1.3208, 1.3255, 1.3301,
28  1.3346, 1.3391, 1.3435, 1.3478, 1.3521, 1.3563, 1.3604, 1.3644,
29  1.3684, 1.3723, 1.3762, 1.3800, 1.3841, 1.3884, 1.3927, 1.3969,
30  1.4011, 1.4052, 1.4092, 1.4129, 1.4167, 1.4203, 1.4239, 1.4275,
31  1.4310, 1.4347, 1.4383, 1.4419, 1.4454, 1.4488, 1.4523, 1.4559,
32  1.4595, 1.4630, 1.4664, 1.4698, 1.4731, 1.4762, 1.4792, 1.4821,
33  1.4851, 1.4880, 1.4908, 1.4936, 1.4963, 1.4990, 1.5017, 1.5043,
34  1.5069, 1.5093, 1.5116, 1.5140, 1.5163, 1.5185, 1.5207, 1.5229,
35  1.5250, 1.5271, 1.5292, 1.5312, 1.5332, 1.5352, 1.5371, 1.5390,
36  1.5408, 1.5427, 1.5445, 1.5462, 1.5479, 1.5496, 1.5513, 1.5529,
37  1.5545, 1.5555, 1.5564, 1.5573, 1.5581, 1.5590, 1.5598, 1.5606,
38  1.5613, 1.5620, 1.5627, 1.5634, 1.5641, 1.5647, 1.5653, 1.5659,
39  1.5665, 1.5670, 1.5675, 1.5680, 1.5685, 1.5690, 1.5694, 1.5699,
40  1.5703, 1.5706, 1.5710, 1.5714, 1.5717, 1.5720, 1.5723, 1.5726,
41  1.5729, 1.5732, 1.5734, 1.5737, 1.5739, 1.5741, 1.5743, 1.5745,
42  1.5746, 1.5748, 1.5750, 1.5751, 1.5752, 1.5754, 1.5755, 1.5756
43  };
44 
45  const int N = 256;
46  const double lnbgMin = -1.011739;
47  const double lnbgStep = 0.038083;
48 
49  int j = (lnbg - lnbgMin) / lnbgStep;
50 
51  if (j > 0)
52  {
53  if (j >= N-1) j = N-2;
54  const double dx = (lnbg - (j * lnbgStep + lnbgMin));
55  return dx * (table[j+1] - table[j]) / lnbgStep + table[j];
56  }
57 
58  else
59  {
60  const double lnDedx = -1.8*(lnbg - lnbgMin) + std::log(table[0]);
61  return std::exp(lnDedx);
62  }
63 }
64 
65 
66 
67 //__________________________________________________________________________________________________________________________
68 TGraph * GetBBGraph(const TH2D * h, const double mass, const double xMin, const double xMax, const Color_t color)
69 {
70  TGraph* g = new TGraphErrors();
71  for (int i = 0; i <= h->GetNbinsX(); i++)
72  {
73  g->SetPoint(i, h->GetXaxis()->GetBinCenter(i+1), BetheBlochAntoniMod(h->GetXaxis()->GetBinCenter(i+1)/mass));
74  }
75 
76  g->GetXaxis()->SetRangeUser(xMin, xMax);
77  g->SetLineColor(color);
78  g->SetLineWidth(5);
79  return g;
80 }
81 
82 
83 
84 
85 
86 //__________________________________________________________________________________________________________________________
87 void MakeBetheBlochPlot (TH2D * h, int charge, bool targetIn, const string &year)
88 {
89  if (h == NULL)
90  {
92  cerr<<"[WARNING]: Specified histogram not found!"<<endl;
93  return;
94  }
95 
96  TCanvas * c1 = new TCanvas("c1", "c1", 1800, 1200);
97 
98  gStyle->SetOptStat(0); // set stat box to show number of entries
99  gPad->SetGridx(1);
100  gPad->SetGridy(1);
101 
102  c1->SetLogz(1);
103  c1->SetTopMargin(0.04);
104  c1->SetBottomMargin(0.09);
105  c1->SetRightMargin(0.12);
106  c1->SetLeftMargin(0.07);
107 
108  h->GetZaxis()->SetTitle("Counts");
109  h->GetXaxis()->SetTitleOffset(1.0);
110  h->GetYaxis()->SetTitleOffset(1.0);
111  h->GetZaxis()->SetTitleOffset(0.9);
112 
113  // h->RebinY(4);
114  // h->SetAxisRange(4, 15, "X");
115  h->SetMinimum(1);
116  //h->SetMaximum(zmax);
117 
118  h->Draw("colz");
119 
120  TGraph * G_bbDeuteron = GetBBGraph(h, 1.875612, 4, 14, kTeal);
121  TGraph * G_bbProton = GetBBGraph(h, 0.938272, 4, 14, kRed);
122  TGraph * G_bbKaon = GetBBGraph(h, 0.493677, 4, 14, kBlue);
123  TGraph * G_bbPion = GetBBGraph(h, 0.139570, 4, 14, kGreen+3);
124  TGraph * G_bbElectron = GetBBGraph(h, 0.000001, 4, 14, kMagenta);
125 
126  G_bbDeuteron->Draw("same");
127  G_bbProton->Draw("same");
128  G_bbKaon->Draw("same");
129  G_bbPion->Draw("same");
130  G_bbElectron->Draw("same");
131 
132  TGraph * G_bbPionCut = GetBBGraph(h, 0.335, 4, 14, kBlack);
133  G_bbPionCut->Draw("same");
134 
135  TGraph * G_bbElectronCut = GetBBGraph(h, 0.02, 4, 14, kBlack);
136  G_bbElectronCut->Draw("same");
137 
139  // TF1 * bb_lin_fit = new TF1("bb_lin_fit", "[0]*x+[1]", 4, 14);
140  // G_bbProton->Fit(bb_lin_fit, "Q"); //"V"
141  // bb_lin_fit->Draw("same");
142 
143 
144 
145  string name = "";
146  if (charge == 1)
147  {
148  if (targetIn) name = "plots/plots_BethBloch/momentum_vs_dEdx_positiveCharge_targetIN_" + year + ".png";
149  else name = "plots/plots_BethBloch/momentum_vs_dEdx_positiveCharge_targetOUT_" + year + ".png";
150  }
151  else
152  {
153  if (targetIn) name = "plots/plots_BethBloch/momentum_vs_dEdx_negativeCharge_targetIN_" + year + ".png";
154  else name = "plots/plots_BethBloch/momentum_vs_dEdx_negativeCharge_targetOUT_" + year + ".png";
155  }
156 
157  c1->Update();
158  c1->SaveAs(name.c_str());
159  delete c1;
160 }
161 
162 
163 
164 
165 
Generated by   doxygen 1.8.5