simple-tof-analysis
 All Classes Namespaces Functions Variables Groups Pages
MassSquaredFits.h
1 using namespace std;
2 
3 
4 const double massSqP = 0.88; //0.8803544;
5 const double massSqK = 0.24; //0.2437169;
6 const double massSqD = 3.52;
7 
8 
9 //__________________________________________________________________________________________________________________________
10 vector<TF1*> GetProtonKaonDeuteronCombinedFits(TH1D * hMassSq_ProtonKaonDeuteron)
11 {
12  /*
13  massPi = 0.13957;
14  massK = 0.493677;
15  massP = 0.93827208;
16  massD = 1.8756129;
17  */
18 
19  double nData = hMassSq_ProtonKaonDeuteron->GetEntries();
20 
21  double massDelta = 0.05;
22  double massSqDDelta = 0.4;
23  double m2_sig1 = 0.15;
24  double m2_sig1_min = 0.10;
25  double m2_sig1_max = 0.20;
26  double m2_sig2 = 0.50;
27  double m2_sig2_min = 0.20;
28  double m2_sig2_max = 5.50; //1.50
29  gErrorIgnoreLevel = kFatal;
30 
31 
32  TF1 * FitK = new TF1("FitK", "[0]*([1]*TMath::Gaus(x, [2], [3]) + (1-[1])*TMath::Gaus(x, [2], [4]))", -5, 5);
33  FitK->SetParameters(2000, 0.5, massSqK, m2_sig1, m2_sig2);
34  FitK->SetParLimits(0, 1, nData);
35  FitK->SetParLimits(1, 0, 1);
36  FitK->SetParLimits(2, massSqK-massDelta, massSqK+massDelta);
37  FitK->SetParLimits(3, m2_sig1_min, m2_sig1_max);
38  FitK->SetParLimits(4, m2_sig2_min, m2_sig2_max);
39  hMassSq_ProtonKaonDeuteron->Fit(FitK, "NQ", "", -5, 5); //(FitK, "NQ", "", -5, 5)
40 
41  TF1 * FitP = new TF1("FitP", "[0]*([1]*TMath::Gaus(x, [2], [3]) + (1-[1])*TMath::Gaus(x, [2], [4]))", -5, 5);
42  FitP->SetParameters(2000, 0.5, massSqP, m2_sig1, m2_sig2);
43  FitP->SetParLimits(0, 0, nData);
44  FitP->SetParLimits(1, 0, 1);
45  FitP->SetParLimits(2, massSqP-massDelta, massSqP+massDelta);
46  FitP->SetParLimits(3, m2_sig1_min, m2_sig1_max);
47  FitP->SetParLimits(4, m2_sig2_min, m2_sig2_max);
48  hMassSq_ProtonKaonDeuteron->Fit(FitP, "NQ", "", -5, 5);
49 
50  TF1 * FitD = new TF1("FitD", "[0]*([1]*TMath::Gaus(x, [2], [3]) + (1-[1])*TMath::Gaus(x, [2], [4]))", -5, 5);
51  FitD->SetParameters(10, 0.5, massSqD, 2*m2_sig1, 2*m2_sig2);
52  FitD->SetParLimits(0, 0, nData/10);
53  FitD->SetParLimits(1, 0, 1);
54  FitD->SetParLimits(2, massSqD-massSqDDelta, massSqD+massSqDDelta);
55  FitD->SetParLimits(3, m2_sig1_min/2, 2*m2_sig1_max);
56  FitD->SetParLimits(4, m2_sig2_min/2, 2*m2_sig2_max);
57  hMassSq_ProtonKaonDeuteron->Fit(FitD, "NQ", "", -5, 5);
58 
59 
60  TF1 * FitT = new TF1("FitT", "[0]*([1]*TMath::Gaus(x, [2], [3]) + (1-[1])*TMath::Gaus(x, [2], [4])) + [5]*([6]*TMath::Gaus(x, [7], [8]) + (1-[6])*TMath::Gaus(x, [7], [9])) + [10]*([11]*TMath::Gaus(x, [12], [13]) + (1-[11])*TMath::Gaus(x, [12], [14]))", -5, 5);
61  // cout<<FitT->GetNumberFreeParameters()<<endl;
62  // cout<<FitT->GetNumberFreeParameters()<<endl;
63 
64 
65  //
66  // FitT->SetParameter(0, FitK->GetParameter(0));
67  // FitT->SetParameter(1, FitK->GetParameter(1));
68  // FitT->SetParameter(2, FitK->GetParameter(2));
69  // FitT->SetParameter(3, FitK->GetParameter(3));
70  // FitT->SetParameter(4, FitK->GetParameter(4));
71  // FitT->SetParameter(5, FitP->GetParameter(0));
72  // FitT->SetParameter(6, FitP->GetParameter(1));
73  // FitT->SetParameter(7, FitP->GetParameter(2));
74  // FitT->SetParameter(8, FitP->GetParameter(3));
75  // FitT->SetParameter(9, FitP->GetParameter(4));
76  // FitT->SetParameter(10, FitD->GetParameter(0));
77  // FitT->SetParameter(11, FitD->GetParameter(1));
78  // FitT->SetParameter(12, FitD->GetParameter(2));
79  // FitT->SetParameter(13, FitD->GetParameter(3));
80  // FitT->SetParameter(14, FitD->GetParameter(4));
81  double guessParam[15];
82  FitK->GetParameters(&guessParam[0]);
83  FitP->GetParameters(&guessParam[5]);
84  FitD->GetParameters(&guessParam[10]);
85  FitT->SetParameters(guessParam);
86 
87 
88  // fit kaon
89  FitT->SetParLimits(0, 1, nData);
90  FitT->SetParLimits(1, 0, 1);
91  FitT->SetParLimits(2, massSqK-massDelta, massSqK+massDelta);
92  FitT->SetParLimits(3, m2_sig1_min, m2_sig1_max);
93  FitT->SetParLimits(4, m2_sig2_min, m2_sig2_max);
94  // fit proton
95  FitT->SetParLimits(5, 1, nData);
96  FitT->SetParLimits(6, 0, 1);
97  FitT->SetParLimits(7, massSqP-massDelta, massSqP+massDelta);
98  FitT->SetParLimits(8, m2_sig1_min, m2_sig1_max);
99  FitT->SetParLimits(9, m2_sig2_min, m2_sig2_max);
100  // fit deuteron
101  FitT->SetParLimits(10, 1, nData/100);
102  FitT->SetParLimits(11, 0, 1);
103  FitT->SetParLimits(12, massSqD-massDelta, massSqD+massDelta);
104  FitT->SetParLimits(13, m2_sig1_min/2, 2*m2_sig1_max);
105  FitT->SetParLimits(14, m2_sig2_min/2, 3*m2_sig2_max);
106 
107  // execute total fit
108  hMassSq_ProtonKaonDeuteron->Fit(FitT, "NQ", "", -5, 5); //"NQ" "NV"
109 
110 
111 
112  //
113  if (FitT->GetParameter(14) >= 3*m2_sig2_max)
114  cout<<"======================> Deuteron sigma-2: "<<FitT->GetParameter(14)<<endl;
115  else
116  cout<<"Deuteron sigma-2: "<<FitT->GetParameter(14)<<endl;
117  cout<<"Ratio sigma-1 proton/kaon: "<<FitT->GetParameter(8)/FitT->GetParameter(3)<<endl;
118 
119 
120 
121  //Set all individual particle fits to fitted parameters
122  double fitParam[15];
123  FitT->GetParameters(fitParam);
124  FitK->SetParameters(&fitParam[0]);
125  FitP->SetParameters(&fitParam[5]);
126  FitD->SetParameters(&fitParam[10]);
127  gErrorIgnoreLevel = kInfo;
128 
129 
130  vector<TF1*> vFits = {FitT, FitK, FitP, FitD};
131  return vFits;
132 }
133 
134 
135 
136 
137 
138 //__________________________________________________________________________________________________________________________
139 TH1D * GetParticleHistogramFromPionData(const TH1D * hMassSq_Pion, double particle_mass)
140 {
141  //shift pion mass-sq distribution, to get it to peak at a given particle's mass
142  TH1D * h = (TH1D*)hMassSq_Pion->Clone();
143  h->Reset("M");
144 
145  const int ibinPionPeak = hMassSq_Pion->GetMaximumBin();
146  const int iParticlePeak = hMassSq_Pion->GetXaxis()->FindBin(particle_mass);
147 
148  int shift = iParticlePeak - ibinPionPeak + 1;
149  if (fabs(particle_mass - massSqK) < 0.01) shift = shift - 1;
150 
151  for (int i = 1; i <= hMassSq_Pion->GetNbinsX(); i++)
152  {
153  if (hMassSq_Pion->GetBinContent(i - shift) >= 0 && (i - shift) >= 0)
154  h->SetBinContent(i, hMassSq_Pion->GetBinContent(i - shift));
155  else
156  h->SetBinContent(i, 0);
157  }
158 
159  h->SetBinContent(h->GetNbinsX()+1, 0); //remove anything stored in overflow bin
160  return h;
161 }
162 
163 
164 
165 
166 //__________________________________________________________________________________________________________________________
167 TH1D * GetBroadenedParticleHistogramFromPionData(const TH1D * hMassSq_Pion, double particle_mass, double width)
168 {
170  TH1D * hZeroed = (TH1D*)hMassSq_Pion->Clone();
171  hZeroed->Reset("M");
172 
173  const int ibinPionPeak = hMassSq_Pion->GetMaximumBin();
174  const int iBinZeroMass = hMassSq_Pion->GetXaxis()->FindBin(0);
175  int shift_left = ibinPionPeak - iBinZeroMass;
176 
177  for (int i = 1; i <= hMassSq_Pion->GetNbinsX(); i++)
178  {
179  if (hMassSq_Pion->GetBinContent(i + shift_left) >= 0 && (i + shift_left) >= 0)
180  hZeroed->SetBinContent(i, hMassSq_Pion->GetBinContent(i + shift_left));
181  else
182  hZeroed->SetBinContent(i, 0);
183  }
184  hZeroed->SetBinContent(hZeroed->GetNbinsX(), hZeroed->GetBinContent(hZeroed->GetNbinsX()-1)); //set last bin content to same as second last bin
185  TGraph * gZeroed = new TGraph(hZeroed);
186  gZeroed->SetBit(TGraph::kIsSortedX);
187 
189  TH1D * hZeroedBroad = (TH1D*)hMassSq_Pion->Clone();
190  hZeroedBroad->Reset("M");
191  for (int j = 1; j <= hZeroedBroad->GetNbinsX(); j++)
192  {
193  double bin_content = gZeroed->Eval(hZeroedBroad->GetBinCenter(j)/width, 0, "S");
194  double bin_error = GetErrorN(bin_content)[0];
195  hZeroedBroad->SetBinContent(j, bin_content);
196  hZeroedBroad->SetBinError(j, bin_error);
197  }
198 
199 
201  TH1D * hBroadShifted = (TH1D*)hMassSq_Pion->Clone();
202  hBroadShifted->Reset("M");
203 
204  const int iBinZero = hZeroedBroad->GetMaximumBin();
205  const int iParticlePeak = hZeroedBroad->GetXaxis()->FindBin(particle_mass);
206  int shift = iParticlePeak - iBinZero + 1;
207  if (fabs(particle_mass - massSqK) < 0.01) shift = shift - 1;
208 
209  for (int i = 1; i <= hBroadShifted->GetNbinsX(); i++)
210  {
211  if (hZeroedBroad->GetBinContent(i - shift) >= 0 && (i - shift) >= 0)
212  hBroadShifted->SetBinContent(i, hZeroedBroad->GetBinContent(i - shift));
213  else
214  hBroadShifted->SetBinContent(i, 0);
215  }
216 
217 
218  hBroadShifted->SetBinContent(hBroadShifted->GetNbinsX()+1, 0); //remove anything stored in overflow bin
219  hBroadShifted->SetBinContent(0, 0); //remove anything stored in underflow bin
220  return hBroadShifted;
221 }
222 
223 
224 
225 
226 
227 
228 //__________________________________________________________________________________________________________________________
230 vector<TH1D*> GetParticleFitsFromPionData(TH1D * hMassSq_ProtonKaonDeuteron, TH1D * hMassSq_Pion, int ivol)
231 {
232  TH1D * hMassSq_PionShiftedToProton = GetParticleHistogramFromPionData(hMassSq_Pion, massSqP);
233  TGraph * GPionShiftedToProton = new TGraph(hMassSq_PionShiftedToProton);
234 
235  TH1D * hMassSq_PionShiftedToKaon = GetParticleHistogramFromPionData(hMassSq_Pion, massSqK);
236  TGraph * GPionShiftedToKaon = new TGraph(hMassSq_PionShiftedToKaon);
237 
238  TH1D * hMassSq_PionShiftedToDeuteron = GetParticleHistogramFromPionData(hMassSq_Pion, massSqD);
239  hMassSq_PionShiftedToDeuteron->Scale(0.9);
240  TGraph * GPionShiftedToDeuteron = new TGraph(hMassSq_PionShiftedToDeuteron);
241 
242 
243  //
244  Print1DHist(hMassSq_Pion, "plots_MassSq_shifted/", "HighStatistics", 1, true);
245  hMassSq_PionShiftedToProton->SetName(Form("Vol-%d_dEdx-massSq_PionShiftedToProton", ivol));
246  Print1DHist(hMassSq_PionShiftedToProton, "plots_MassSq_shifted/", "HighStatistics", 1, true);
247  hMassSq_PionShiftedToKaon->SetName(Form("Vol-%d_dEdx-massSq_PionShiftedToKaon", ivol));
248  Print1DHist(hMassSq_PionShiftedToKaon, "plots_MassSq_shifted/", "HighStatistics", 1, true);
249  hMassSq_PionShiftedToDeuteron->SetName(Form("Vol-%d_dEdx-massSq_PionShiftedToDeuteron", ivol));
250  Print1DHist(hMassSq_PionShiftedToDeuteron, "plots_MassSq_shifted/", "HighStatistics", 1, true);
251 
252 
253  //
254  GPionShiftedToProton->SetBit(TGraph::kIsSortedX);
255  GPionShiftedToKaon->SetBit(TGraph::kIsSortedX);
256  GPionShiftedToDeuteron->SetBit(TGraph::kIsSortedX);
257  TF1 * FitT = new TF1("FitT",[&](double * x, double * p){return p[0]*GPionShiftedToProton->Eval(x[0], 0, "S") + p[1]*GPionShiftedToKaon->Eval(x[0], 0, "S") + p[2]*GPionShiftedToDeuteron->Eval(x[0], 0, "S");}, -5, 5, 3);
258 
259  // execute total fit
260  hMassSq_ProtonKaonDeuteron->Fit(FitT, "NQ", "", -2.5, 5); //"NQ" "NV"
261 
262 
263  //
264  double Param[3];
265  FitT->GetParameters(Param);
266 
267  // Set deuteron to zero if fit is weird
268  if (Param[2] > Param[0] || Param[2] < 0)
269  Param[2] = 0.00000001;
270 
271 
272  vector<TH1D*> vFits;
273  if (Param[0] >= 0 && Param[1] >= 0 && Param[2] >= 0)
274  {
275  cout<<"Fit params: "<<Param[0]<<", "<<Param[1]<<", "<<Param[2]<<endl;
276 
277  //______________________
278  hMassSq_PionShiftedToProton->Scale(Param[0]);
279  TH1D * hhMassSq_TotalFit = (TH1D*)hMassSq_PionShiftedToProton->Clone();
280  hMassSq_PionShiftedToKaon->Scale(Param[1]);
281  hhMassSq_TotalFit->Add(hMassSq_PionShiftedToKaon);
282  hMassSq_PionShiftedToDeuteron->Scale(Param[2]);
283  hhMassSq_TotalFit->Add(hMassSq_PionShiftedToDeuteron);
284  /*
285  hMassSq_PionShiftedToProton = GetScaledTH1D(Param[0], hMassSq_PionShiftedToProton);
286  TH1D * hhMassSq_TotalFit = (TH1D*)hMassSq_PionShiftedToProton->Clone();
287  hMassSq_PionShiftedToKaon = GetScaledTH1D(Param[1], hMassSq_PionShiftedToKaon);
288  hhMassSq_TotalFit->Add(hMassSq_PionShiftedToKaon);
289  hMassSq_PionShiftedToDeuteron = GetScaledTH1D(Param[2], hMassSq_PionShiftedToDeuteron);
290  hhMassSq_TotalFit->Add(hMassSq_PionShiftedToDeuteron);
291  */
292 
293 
294  //______________________
295  hhMassSq_TotalFit->SetName(Form("Vol-%d_dEdx-massSq_FinalFit", ivol));
296  Print1DHist(hhMassSq_TotalFit, "plots_MassSq_shifted/", "HighStatistics", 1, true);
297  vFits = {hhMassSq_TotalFit, hMassSq_PionShiftedToKaon, hMassSq_PionShiftedToProton, hMassSq_PionShiftedToDeuteron};
298  }
299 
300  else
301  {
302  cout<<"[WARNING]: Fit params: "<<Param[0]<<", "<<Param[1]<<", "<<Param[2]<<". Fitting unsuccessful for Volume "<<ivol<<endl;
303  vFits = {};
304  }
305 
306  return vFits;
307 }
308 
309 
310 
311 
312 
313 
314 //__________________________________________________________________________________________________________________________
316 vector<TH1D*> GetParticleFitsFromBroadenedPionData(TH1D * hMassSq_ProtonKaonDeuteron,
317  TH1D * hMassSq_Pion,
318  int ivol,
319  double pBinCenter,
320  double &d_norm_PercentErr,
321  double &chiSquaredperNDF)
322 {
323  TFile * fPionTemplateWidths = new TFile("MassSquarePeakBroadening_Philip/PionTemplateWidths_backup.root");
324  if (pBinCenter > 11) pBinCenter = 11;
325  TGraphErrors * GParticlePeakWidth = (TGraphErrors*)fPionTemplateWidths->Get(Form("%2.1f-GeV_MSquaredResolutionNormalizedPion", pBinCenter));
326  GParticlePeakWidth->SetBit(TGraph::kIsSortedX);
327  double width_k = GParticlePeakWidth->Eval(massSqK, 0, "S");
328  double width_p = GParticlePeakWidth->Eval(massSqP, 0, "S");
329  double width_d = GParticlePeakWidth->Eval(massSqD, 0, "S");
330 
331 
332  TH1D * hMassSq_PionShiftedToProton = GetBroadenedParticleHistogramFromPionData(hMassSq_Pion, massSqP, width_p);
333  TGraph * GPionShiftedToProton = new TGraph(hMassSq_PionShiftedToProton);
334 
335  TH1D * hMassSq_PionShiftedToKaon = GetBroadenedParticleHistogramFromPionData(hMassSq_Pion, massSqK, width_k);
336  TGraph * GPionShiftedToKaon = new TGraph(hMassSq_PionShiftedToKaon);
337 
338  TH1D * hMassSq_PionShiftedToDeuteron = GetBroadenedParticleHistogramFromPionData(hMassSq_Pion, massSqD, width_d);
339  hMassSq_PionShiftedToDeuteron->Scale(0.9);
340  TGraph * GPionShiftedToDeuteron = new TGraph(hMassSq_PionShiftedToDeuteron);
341 
342 
343  //
344  Print1DHist(hMassSq_Pion, "plots_MassSq_shifted_broad/", "HighStatistics", 1, true);
345  cout<<"RMS of hMassSq_Pion: "<<hMassSq_Pion->GetRMS()<<endl;
346  hMassSq_PionShiftedToProton->SetName(Form("Vol-%d_dEdx-massSq_PionShiftedToProton", ivol));
347  Print1DHist(hMassSq_PionShiftedToProton, "plots_MassSq_shifted_broad/", "HighStatistics", 1, true);
348  cout<<"RMS of hMassSq_PionShiftedToProton: "<<hMassSq_PionShiftedToProton->GetRMS()<<endl;
349  hMassSq_PionShiftedToKaon->SetName(Form("Vol-%d_dEdx-massSq_PionShiftedToKaon", ivol));
350  Print1DHist(hMassSq_PionShiftedToKaon, "plots_MassSq_shifted_broad/", "HighStatistics", 1, true);
351  cout<<"RMS of hMassSq_PionShiftedToKaon: "<<hMassSq_PionShiftedToKaon->GetRMS()<<endl;
352  hMassSq_PionShiftedToDeuteron->SetName(Form("Vol-%d_dEdx-massSq_PionShiftedToDeuteron", ivol));
353  Print1DHist(hMassSq_PionShiftedToDeuteron, "plots_MassSq_shifted_broad/", "HighStatistics", 1, true);
354  cout<<"RMS of hMassSq_PionShiftedToDeuteron: "<<hMassSq_PionShiftedToDeuteron->GetRMS()<<endl;
355 
356 
357  //
358  GPionShiftedToProton->SetBit(TGraph::kIsSortedX);
359  GPionShiftedToKaon->SetBit(TGraph::kIsSortedX);
360  GPionShiftedToDeuteron->SetBit(TGraph::kIsSortedX);
361  TF1 * FitT = new TF1("FitT",[&](double * x, double * p){return p[0]*GPionShiftedToProton->Eval(x[0], 0, "S") + p[1]*GPionShiftedToKaon->Eval(x[0], 0, "S") + p[2]*GPionShiftedToDeuteron->Eval(x[0], 0, "S");}, -5, 5, 3);
362 
363  // execute total fit
364  auto fit_result_ptr = hMassSq_ProtonKaonDeuteron->Fit(FitT, "SNQ", "", -2.5, 5); //"NQ" "NV"
365 
366 
367  //
368  cout<<"\nChi-squared/NDF = "<<fit_result_ptr->Chi2()/fit_result_ptr->Ndf()<<endl<<endl;
369 
370 
371  //
372  double Param[3];
373  FitT->GetParameters(Param);
374 
375  double ParEr[3];
376  ParEr[0] = FitT->GetParError(0);
377  ParEr[1] = FitT->GetParError(1);
378  ParEr[2] = FitT->GetParError(2);
379 
380  // Set deuteron to zero if fit is weird
381  if (Param[2] > Param[0] || Param[2] < 0)
382  Param[2] = 0.00000001;
383 
384 
385  vector<TH1D*> vFits;
386  if (Param[0] >= 0 && Param[1] >= 0 && Param[2] >= 0)
387  {
388  cout<<endl;
389  cout<<"Fit params : "<<Param[0]<<", "<<Param[1]<<", "<<Param[2]<<endl;
390  cout<<"Fit errors : "<<ParEr[0]<<", "<<ParEr[1]<<", "<<ParEr[2]<<endl;
391  cout<<"Percent err : "<<100*ParEr[0]/Param[0]<<", "<<100*ParEr[1]/Param[1]<<", "<<100*ParEr[2]/Param[2]<<endl;
392  cout<<endl;
393 
394  //______________________
395  hMassSq_PionShiftedToProton->Scale(Param[0]);
396  TH1D * hhMassSq_TotalFit = (TH1D*)hMassSq_PionShiftedToProton->Clone();
397  hMassSq_PionShiftedToKaon->Scale(Param[1]);
398  hhMassSq_TotalFit->Add(hMassSq_PionShiftedToKaon);
399  hMassSq_PionShiftedToDeuteron->Scale(Param[2]);
400  hhMassSq_TotalFit->Add(hMassSq_PionShiftedToDeuteron);
401  /*
402  hMassSq_PionShiftedToProton = GetScaledTH1D(Param[0], hMassSq_PionShiftedToProton);
403  TH1D * hhMassSq_TotalFit = (TH1D*)hMassSq_PionShiftedToProton->Clone();
404  hMassSq_PionShiftedToKaon = GetScaledTH1D(Param[1], hMassSq_PionShiftedToKaon);
405  hhMassSq_TotalFit->Add(hMassSq_PionShiftedToKaon);
406  hMassSq_PionShiftedToDeuteron = GetScaledTH1D(Param[2], hMassSq_PionShiftedToDeuteron);
407  hhMassSq_TotalFit->Add(hMassSq_PionShiftedToDeuteron);
408  */
409 
410 
411  //______________________
412  hhMassSq_TotalFit->SetName(Form("Vol-%d_massSq_Fit_total", ivol));
413  Print1DHist(hhMassSq_TotalFit, "plots_MassSq_shifted_broad/", "HighStatistics", 1, true);
414 
415  hMassSq_PionShiftedToKaon->SetName(Form("Vol-%d_massSq_Fit_K", ivol));
416  hMassSq_PionShiftedToProton->SetName(Form("Vol-%d_massSq_Fit_P", ivol));
417  hMassSq_PionShiftedToDeuteron->SetName(Form("Vol-%d_massSq_Fit_D", ivol));
418 
419 
420  d_norm_PercentErr = 100*ParEr[2]/Param[2];
421  chiSquaredperNDF = fit_result_ptr->Chi2()/fit_result_ptr->Ndf();
422  vFits = {hhMassSq_TotalFit, hMassSq_PionShiftedToKaon, hMassSq_PionShiftedToProton, hMassSq_PionShiftedToDeuteron};
423  }
424 
425  else
426  {
427  cout<<"[WARNING]: Fit params: "<<Param[0]<<", "<<Param[1]<<", "<<Param[2]<<". Fitting unsuccessful for Volume "<<ivol<<endl;
428  vFits = {};
429  }
430 
431  return vFits;
432 }
433 
434 
435 
436 
437 
438 
439 //__________________________________________________________________________________________________________________________
440 void RunFitDiagnostics(TH1D * hMassSq_FittedProton, TH1D * hMassSq_FittedKaon, TH1D * hMassSq_FittedDeuteron, TH1D * hMassSq_TotalFit, int ivol)
441 {
442  /*
443  TH1D * hMassSq_FitResidual = (TH1D*)hMassSq_TotalFit->Clone();
444  hMassSq_FitResidual->SetName(Form("Vol-%d_dEdx-massSq_FinalFitResidual", ivol));
445  hMassSq_FitResidual->Add(hMassSq_FittedProton, -1);
446  hMassSq_FitResidual->Add(hMassSq_FittedKaon, -1);
447  hMassSq_FitResidual->Add(hMassSq_FittedDeuteron, -1);
448  Print1DHist(hMassSq_FitResidual, "plots_1DMassSqFitsTFractionFitter/FitResiduals/", "HighStatistics", 1, true);
449  */
450 
451  TH1D * hMassSq_FitSum = (TH1D*)hMassSq_FittedProton->Clone();
452  hMassSq_FitSum->Add(hMassSq_FittedKaon);
453  hMassSq_FitSum->Add(hMassSq_FittedDeuteron);
454  TH1D * hRatio = GetRatioTH1D(hMassSq_TotalFit, hMassSq_FitSum);
455  hRatio->SetName(Form("Vol-%d_dEdx-massSq_RatioTotalFitOverSumFit", ivol));
456  hRatio->SetAxisRange(0, 2, "Y");
457  for (int i = 0; i <= hRatio->GetNbinsX(); i++)
458  {
459  hRatio->SetBinError(i, 0);
460  if (hRatio->GetBinCenter(i) < -4)
461  hRatio->SetBinContent(i, 1);
462  }
463  Print1DHist(hRatio, "plots_1DMassSqFitsTFractionFitter/FitResiduals/", "HighStatistics", 1, true, 0);
464 }
465 
466 
467 
468 
469 
470 //__________________________________________________________________________________________________________________________
472 vector<TH1D*> GetParticleFitsUsingTFractionFitter(TH1D * hMassSq_ProtonKaonDeuteron, TH1D * hMassSq_Pion, int ivol)
473 {
474  const double massSqP = 0.8803544;
475  const double massSqK = 0.2437169;
476  const double massSqD = 3.6; //3.52;
477 
478  //
479  TH1D * hMassSq_PionShiftedToProton = GetParticleHistogramFromPionData(hMassSq_Pion, massSqP);
480  TH1D * hMassSq_PionShiftedToKaon = GetParticleHistogramFromPionData(hMassSq_Pion, massSqK);
481  TH1D * hMassSq_PionShiftedToDeuteron = GetParticleHistogramFromPionData(hMassSq_Pion, massSqD);
482  hMassSq_PionShiftedToDeuteron->Scale(0.1);
483 
484 
485  //
486  TObjArray * vHistMC = new TObjArray(3); // MC histograms are put in this array
487  vHistMC->Add(hMassSq_PionShiftedToProton);
488  vHistMC->Add(hMassSq_PionShiftedToKaon);
489  vHistMC->Add(hMassSq_PionShiftedToDeuteron);
490  TFractionFitter * fit = new TFractionFitter(hMassSq_ProtonKaonDeuteron, vHistMC, "Q"); // initialise, option "Q" or "V".
491  fit->Constrain(0, 0.0, 1.0); // constrain fraction 1 to be between 0 and 1
492  fit->Constrain(1, 0.0, 1.0); // constrain fraction 1 to be between 0 and 1
493  fit->Constrain(2, 0.0, 1.0); // constrain fraction 1 to be between 0 and 0.1
494 
495 
496  vector<TH1D*> vFits;
497  Int_t status = fit->Fit(); // perform the fit
498  if (status == 0) // check on fit status
499  {
500  // cout<<"\n\n\nFit completed!\n\n"<<endl;
501  TH1D * hMassSq_TotalFit = (TH1D*) fit->GetPlot();
502 
503  //
504  vector<double> Param;
505  double param_values;
506  double param_errors;
507  fit->GetResult(0, param_values, param_errors); Param.push_back(param_values); //proton
508  fit->GetResult(1, param_values, param_errors); Param.push_back(param_values); //kaon
509  fit->GetResult(2, param_values, param_errors); Param.push_back(param_values); //deuteron
510 
511  TH1D * hMassSq_FittedProton = (TH1D*)fit->GetMCPrediction(0);
512  TH1D * hMassSq_FittedKaon = (TH1D*)fit->GetMCPrediction(1);
513  TH1D * hMassSq_FittedDeuteron = (TH1D*)fit->GetMCPrediction(2);
514 
515  //scale correctly to get the final fitted fractional compoments
516  hMassSq_FittedProton->Scale(1./hMassSq_FittedProton->Integral()*Param[0]*hMassSq_ProtonKaonDeuteron->Integral());
517  hMassSq_FittedKaon->Scale(1./hMassSq_FittedKaon->Integral()*Param[1]*hMassSq_ProtonKaonDeuteron->Integral());
518  hMassSq_FittedDeuteron->Scale(1./hMassSq_FittedDeuteron->Integral()*Param[2]*hMassSq_ProtonKaonDeuteron->Integral());
519 
520  //diagnostics
521  RunFitDiagnostics(hMassSq_FittedProton, hMassSq_FittedKaon, hMassSq_FittedDeuteron, hMassSq_TotalFit, ivol);
522 
523  //
524  cout<<"Fit params: "<<Param[0]<<", "<<Param[1]<<", "<<Param[2]<<endl;
525  vFits = {hMassSq_TotalFit, hMassSq_FittedKaon, hMassSq_FittedProton, hMassSq_FittedDeuteron};
526  }
527 
528  else
529  {
530  cout<<"[WARNING]: Fit status: "<<status<<" for Volume "<<ivol<<endl;
531  vFits = {};
532  }
533 
534  return vFits;
535 }
536 
537 
538 
539 
540 
541 
542 
543 
Generated by   doxygen 1.8.5