simple-tof-analysis
 All Classes Namespaces Functions Variables Groups Pages
FitSpectra.h
1 using namespace std;
2 
3 #include <TFitResult.h>
4 #include "TVirtualFitter.h"
5 
6 
7 
8 
9 //__________________________________________________________________________________________________________________________
10 TF1 * GetDeuteronFitFunction(double fitMinPt = 0.0, double fitMaxPt = 0.6, const string &particle_name = "deuteron")
11 {
12  double mass = 0;
13  if (particle_name=="proton")
14  mass = 0.93827208;
15  else if (particle_name=="deuteron")
16  mass = 1.8756129;
17 
18  // TF1 * FitD = new TF1("FitD", "([0]*x)/([1]*[1] + 1.8756*[1])*exp(-(sqrt(x*x + 1.8756*1.8756) - 1.8756)/[1])", 0, 0.6); ///massD = 1.8756;
19  // TF1 * FitD = new TF1("FitD", "[0]*x*exp(-(sqrt(x*x + 1.8756*1.8756) - 1.8756)/[1])", 0, 0.6); ///massD = 1.8756;
20 
21  TF1 * FitD = new TF1("FitD", Form("[0]*x*exp(-(sqrt(x*x + %f*%f) - %f)/[1])", mass, mass, mass), fitMinPt, fitMaxPt);
22 
23  // TF1 * FitD = new TF1("FitD", Form("[0]*x/([1]*[1] + [1]*%f) * exp(-(sqrt(x*x + %f*%f) - %f)/[1])", mass, mass, mass, mass), fitMinPt, fitMaxPt);
24 
25  return FitD;
26 }
27 
28 
29 
30 
31 //__________________________________________________________________________________________________________________________
32 vector<TF1*> GetPtSpectraFit (const vector<TGraphErrors*> vDataGraphs,
33  vector<TGraphErrors*> &vGConfidenceBands,
34  const string &particle_name,
35  bool fixT = true,
36  double T_initial_guess = 0.0373538,
37  double fitMinPt = 0.0,
38  double fitMaxPt = 0.6,
39  bool LimitParameters = true)
40 {
41  vector<TF1*> vFits;
42 
43  for (int i = 0; i < int(vDataGraphs.size()); i++)
44  // for (int i = 0; i < 1; i++)
45  {
47  TF1 * FitD = GetDeuteronFitFunction(fitMinPt, fitMaxPt, particle_name);
48 
49  FitD->SetParameters(1e-4, T_initial_guess);
50  if (LimitParameters == true)
51  {
52  FitD->SetParLimits(0, 1e-4, 6e-3);
53  FitD->SetParLimits(1, 0.010, 0.600);
54  }
55  if (fixT == true)
56  FitD->FixParameter (1, T_initial_guess); //0.0373538 for deuterons.
57 
58 
60  TFitResultPtr fit_ptr = vDataGraphs[i]->Fit(FitD, "NQSWEM", "", fitMinPt, fitMaxPt);
61 
62 
64  int ngr = vDataGraphs[i]->GetN();
65  TGraphErrors * grint = new TGraphErrors(ngr);
66  grint->SetName(Form("graph_fit_confdinterval_%d", i));
67  for (int k = 0; k < ngr; k++)
68  grint->SetPoint(k, vDataGraphs[i]->GetX()[k], 0);
69 
71  (TVirtualFitter::GetFitter())->GetConfidenceIntervals(grint, 0.95);
72 
85  vGConfidenceBands.push_back(grint);
86  vFits.push_back(FitD);
87 
88  double Param[2];
89  FitD->GetParameters(Param);
90 
91  double ParEr[2];
92  ParEr[0] = FitD->GetParError(0);
93  ParEr[1] = FitD->GetParError(1);
94  // cerr<<ParEr[0]<<", "<<ParEr[1]<<endl;
95 
96  cerr<<Form("Fit completed. Normalization = %e +/-%6.2f%%, T = %.8f +/-%6.2f%%, Chi-squared/NDF = %e", Param[0], 100*ParEr[0]/Param[0], Param[1], 100*ParEr[1]/Param[1], fit_ptr->Chi2()/fit_ptr->Ndf())<<endl;
97  }
98 
99  if (vGConfidenceBands.size() != vFits.size())
100  cerr<<"\n\n[WARNING] vGConfidenceBands size != vFits size"<<endl;
101  return vFits;
102 }
103 
104 
105 
106 
107 //__________________________________________________________________________________________________________________________
108 vector<TF1*> GetPtSpectraFitAndEval (vector<TGraphErrors*> vDataGraphs,
109  vector<TGraphErrors*> &vGConfidenceBands,
110  vector<int> &vMissingPointIndex,
111  const string &particle_name,
112  bool fixT = true,
113  double T_initial_guess = 0.0373538,
114  double fitMinPt = 0.0,
115  double fitMaxPt = 0.6,
116  bool add_missing_points_from_fit = false)
117 {
119  fixT = false; if (particle_name == "deuteron") fixT = true;
120  T_initial_guess = 0.14; if (particle_name == "deuteron") T_initial_guess = 0.0373538;
121  fitMinPt = 0.2; if (particle_name == "deuteron") fitMinPt = 0.0;
122  fitMaxPt = 1.0; if (particle_name == "deuteron") fitMaxPt = 0.6;
123  add_missing_points_from_fit = false; if (particle_name == "deuteron") add_missing_points_from_fit = true;
124 
125 
126  vector<TF1*> vFits;
127 
128  for (int i = 0; i < int(vDataGraphs.size()); i++)
129  // for (int i = 0; i < 1; i++)
130  {
132  TF1 * FitD = GetDeuteronFitFunction(fitMinPt, fitMaxPt);
133 
134  FitD->SetParameters(1e-4, T_initial_guess);
135  if (fixT == true)
136  FitD->FixParameter (1, T_initial_guess); //0.0373538 for deuterons.
137 
138 
139  int ngr = vDataGraphs[i]->GetN();
140  TGraphErrors * grint = new TGraphErrors(ngr);
141  grint->SetName(Form("graph_fit_confdinterval_%d", i));
142  if (ngr == 0)
143  {
144  cerr<<"\n\n[WARNING] TGraphError is empty for i="<<i<<". Setting fit to zero!"<<endl;
145  FitD->FixParameter (0, 0);
146  vFits.push_back(FitD);
147  vGConfidenceBands.push_back(grint);
148  continue;
149  }
150 
151  TFitResultPtr fit_ptr = vDataGraphs[i]->Fit(FitD, "NQS", "", fitMinPt, fitMaxPt);
152 
153 
155  if (add_missing_points_from_fit == true)
156  AddMissingPointAndErrorToGraphFromFit(vDataGraphs[i], fit_ptr, FitD, vMissingPointIndex, true);
157 
158 
160  for (int k = 0; k < ngr; k++)
161  grint->SetPoint(k, vDataGraphs[i]->GetX()[k], 0);
162 
164  (TVirtualFitter::GetFitter())->GetConfidenceIntervals(grint, 0.95);
165 
178  vGConfidenceBands.push_back(grint);
179  vFits.push_back(FitD);
180 
181  double Param[2];
182  FitD->GetParameters(Param);
183  // cerr<<"Spectrum fit completed. Normalization = "<<Param[0]<<", T = "<<Param[1]<<endl;
184 
185  double ParEr[2];
186  ParEr[0] = FitD->GetParError(0);
187  ParEr[1] = FitD->GetParError(1);
188  // cerr<<ParEr[0]<<", "<<ParEr[1]<<endl;
189 
190  cerr<<Form("Spectrum fit completed. Normalization = %e +/-%6.2f%%, T = %.8f +/-%6.2f%%, Chi-squared/NDF = %e", Param[0], 100*ParEr[0]/Param[0], Param[1], 100*ParEr[1]/Param[1], fit_ptr->Chi2()/fit_ptr->Ndf())<<endl;
191  }
192 
193  if (vGConfidenceBands.size() != vFits.size())
194  cerr<<"\n\n[WARNING] Size of vGConfidenceBands = "<<vGConfidenceBands.size()<<". Size of vFits = "<<vFits.size()<<endl<<endl;
195  return vFits;
196 }
197 
198 
199 
200 
201 //__________________________________________________________________________________________________________________________
202 vector<TGraphErrors*> GetFitEvaluatedIntegratedYSpectraGraphs (const vector<TGraphErrors*> vDataGraphs,
203  vector<TGraphErrors*> &vGConfidenceBands,
204  bool fixT = true,
205  double T_initial_guess = 0.0373538,
206  double fitMinPt = 0.0,
207  double fitMaxPt = 0.6
208  )
209 {
210  vector<TGraphErrors*> vIntegratedSpectraGraphs;
211 
212  for (int i = 0; i < int(vDataGraphs.size()); i++)
213  {
215  TF1 * FitD = GetDeuteronFitFunction(fitMinPt, fitMaxPt);
216 
217  FitD->SetParameters(1e-4, T_initial_guess);
218  if (fixT == true)
219  FitD->FixParameter (1, T_initial_guess); //0.0373538 for deuterons.
220 
221 
222  int ngr = vDataGraphs[i]->GetN();
223  if (ngr == 0)
224  {
225  cerr<<"\n\n[WARNING] TGraphError is empty for i="<<i<<". Moving to next..."<<endl;
226 
227  TGraphErrors * GraphYIntegrated = new TGraphErrors(ngr);
228  TGraphErrors * GConfidenceBand = new TGraphErrors(ngr);
229  vIntegratedSpectraGraphs.push_back(GraphYIntegrated);
230  vGConfidenceBands.push_back(GConfidenceBand);
231  continue;
232  }
233 
234 
235 
236 
237  TFitResultPtr fit_ptr = vDataGraphs[i]->Fit(FitD, "NQS", "", fitMinPt, fitMaxPt);
238  double Param[2];
239  FitD->GetParameters(Param);
240  // cerr<<"Spectrum fit completed. Normalization = "<<Param[0]<<", T = "<<Param[1]<<". Chi-squared/NDF = "<<fit_ptr->Chi2()/fit_ptr->Ndf()<<endl;
241 
242  double ParEr[2];
243  ParEr[0] = FitD->GetParError(0);
244  ParEr[1] = FitD->GetParError(1);
245  // cerr<<ParEr[0]<<", "<<ParEr[1]<<endl;
246 
247  cerr<<Form("Spectrum fit completed. Normalization = %e +/-%6.2f%%, T = %.8f +/-%6.2f%%, Chi-squared/NDF = %e", Param[0], 100*ParEr[0]/Param[0], Param[1], 100*ParEr[1]/Param[1], fit_ptr->Chi2()/fit_ptr->Ndf())<<endl;
248 
249 
251  auto GConfidenceBand = new TGraphErrors(1);
252 
254  int nPoints = 100;
255  TGraphErrors * GraphYIntegrated = GetdndyGraphFromFit(fit_ptr, FitD, vDataGraphs[i], GConfidenceBand, nPoints, i);
256 
257  vGConfidenceBands.push_back(GConfidenceBand);
258  vIntegratedSpectraGraphs.push_back(GraphYIntegrated);
259  }
260 
261  if (vGConfidenceBands.size() != vIntegratedSpectraGraphs.size())
262  cerr<<"\n\n[WARNING] Size of vGConfidenceBands = "<<vGConfidenceBands.size()<<". Size of vIntegratedSpectraGraphs = "<<vIntegratedSpectraGraphs.size()<<endl<<endl;
263  return vIntegratedSpectraGraphs;
264 }
Generated by   doxygen 1.8.5