StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
gausSlice.C
1 Double_t gausBack (Double_t *x, Double_t *par)
2 {
3  Double_t sigma = par[3];
4  Double_t mean = par[2];
5  Double_t sigfrac = par[1];
6  Double_t backfrac = 1.-sigfrac;
7  Double_t amp = par[0];
8 
9 // Both have integrals of 1 between +- 2 Sigma
10  Double_t gaus = 1./(sqrt(2.*3.14159265)*sigma*TMath::Erf(sqrt(2)))*
11  exp(-0.5*(x[0]-mean)*(x[0]-mean)/(sigma*sigma));
12  Double_t back = 1./(4.*sigma);
13 
14  Double_t fitval = amp*(backfrac*back + sigfrac*gaus);
15  return fitval;
16 }
17 Double_t erfBack (Double_t *x, Double_t *par)
18 {
19  Double_t sigma = par[3];
20  Double_t mean = par[2];
21  Double_t sigfrac = par[1];
22  Double_t backfrac = 1.-sigfrac;
23  Double_t amp = par[0];
24  Double_t blah1 = TMath::Erf(-(x[0]-0.4-mean)/sigma);
25  Double_t blah2 = TMath::Erf((x[0]+0.4-mean)/sigma);
26  Double_t gaus = 0.5*(blah1+blah2);
27 
28  Double_t fitval =amp*(backfrac + sigfrac*gaus);
29  return fitval;
30 }
31 
32 TGraphErrors**
33 fitSlice (void* func,TH2D* hist, Int_t nhists)
34 {
35  TGraphErrors** retGraphs = new TGraphErrors*[5];
36  TH1D** hists = new TH1D*[nhists];
37  Int_t nbins = hist->GetNbinsX();
38  Double_t* x = new Double_t[nhists];
39  Double_t* ex = new Double_t[nhists];
40 
41  for (Int_t i=0; i< nhists; ++i) {
42  Int_t low = i*(nbins/nhists)+1;//The 0th bin is underflow
43  Int_t high = (i+1)*(nbins/nhists);
44  Double_t xlow = hist->GetXaxis()->GetBinLowEdge(low);
45  Double_t xhigh = hist->GetXaxis()->GetBinUpEdge(high);
46  x[i] = (xlow + xhigh)/2.;
47  ex[i] = (xhigh-xlow)/2.;
48 
49  Char_t* blah = new Char_t[100];
50  Char_t* blah2 = hist->GetName();
51 
52  sprintf(blah,"%s_proj%d",blah2,i);
53 
54  hists[i] = hist->ProjectionY(blah,low,high,"E");
55  }
56 
57  TF1* gausBack = new TF1("gausBack",func,-3,3,4);
58  gausBack->SetParNames("2 Sigma Sum","2 Sigma Signal Fraction","Mean","Sigma");
59  Double_t* amp = new Double_t[nhists];
60  Double_t* sigfrac = new Double_t[nhists];
61  Double_t* mean = new Double_t[nhists];
62  Double_t* sigma = new Double_t[nhists];
63  Double_t* chisq = new Double_t[nhists];
64 
65  Double_t* eamp = new Double_t[nhists];
66  Double_t* esigfrac = new Double_t[nhists];
67  Double_t* emean = new Double_t[nhists];
68  Double_t* esigma = new Double_t[nhists];
69  Double_t* echisq = new Double_t[nhists];
70 
71  for (Int_t i=0; i< nhists; ++i) {
72  TH1D* histslice = hists[i];
73  gausBack->SetParameters(histslice->GetMaximum()*histslice->GetRMS(),1,
74  histslice->GetMean(),histslice->GetRMS());
75  printf("Original parameters: %f %f %f %f",histslice->GetMaximum()*histslice->GetRMS(),
76  1,histslice->GetMean(),histslice->GetRMS());
77 
78  if (histslice->GetMaximum()>4.) {
79 
80  histslice->Fit("gausBack","EL");
81  amp[i] = gausBack->GetParameter(0);
82  sigfrac[i] = gausBack->GetParameter(1);
83  mean[i] = gausBack->GetParameter(2);
84  sigma[i] = gausBack->GetParameter(3);
85  chisq[i] = gausBack->GetChisquare();
86 
87  eamp[i] = gausBack->GetParError(0);
88  esigfrac[i] = gausBack->GetParError(1);
89  emean[i] = gausBack->GetParError(2);
90  esigma[i] = gausBack->GetParError(3);
91  echisq[i] = 0;
92  }
93  else {
94  amp[i] = 0;
95  sigfrac[i] = 0;
96  mean[i] = 0;
97  sigma[i] = 0;
98  chisq[i] = 0;
99 
100  eamp[i] = 0;
101  esigfrac[i] = 0;
102  emean[i] = 0;
103  esigma[i] = 0;
104  echisq[i] = 0;
105  }
106 
107 
108  }
109  retGraphs[0] = new TGraphErrors(nhists,x,amp,ex,eamp);
110  retGraphs[1] = new TGraphErrors(nhists,x,sigfrac,ex,esigfrac);
111  retGraphs[2] = new TGraphErrors(nhists,x,mean,ex,emean);
112  retGraphs[3] = new TGraphErrors(nhists,x,sigma,ex,esigma);
113  retGraphs[4] = new TGraphErrors(nhists,x,chisq,ex,echisq);
114  retGraphs[0]->SetTitle("2 Sigma Sum");
115  retGraphs[1]->SetTitle("2 Sigma Signal Fraction");
116  retGraphs[2]->SetTitle("Mean");
117  retGraphs[3]->SetTitle("Sigma");
118  retGraphs[4]->SetTitle("Chi-squared");
119 
120  return retGraphs;
121 
122 
123 }
124 TGraphErrors**
125 gausSlice (TH2D* hist, Int_t nhists)
126 {
127 
128  TGraphErrors** retval = fitSlice(gausBack,hist,nhists);
129  return retval;
130 }
131 TGraphErrors**
132 erfSlice (TH2D* hist, Int_t nhists)
133 {
134 
135  TGraphErrors** retval = fitSlice(erfBack,hist,nhists);
136  return retval;
137 }
138 
139 TGraphErrors**
140 funcFit (void* func,TH1D* histslice)
141 {
142  TGraphErrors** retGraphs = new TGraphErrors*[5];
143  Int_t nhists=1;
144  Int_t i=1;
145 
146  TF1* gausBack = new TF1("gausBack",func,-3,3,4);
147  gausBack->SetParNames("2 Sigma Sum","2 Sigma Signal Fraction","Mean","Sigma");
148  Double_t* amp = new Double_t[nhists];
149  Double_t* sigfrac = new Double_t[nhists];
150  Double_t* mean = new Double_t[nhists];
151  Double_t* sigma = new Double_t[nhists];
152  Double_t* chisq = new Double_t[nhists];
153 
154  Double_t* eamp = new Double_t[nhists];
155  Double_t* esigfrac = new Double_t[nhists];
156  Double_t* emean = new Double_t[nhists];
157  Double_t* esigma = new Double_t[nhists];
158  Double_t* echisq = new Double_t[nhists];
159  gausBack->SetParameters(histslice->GetMaximum()*histslice->GetRMS(),1,
160  histslice->GetMean(),histslice->GetRMS());
161 
162  printf("Original parameters: %f %f %f %f",histslice->GetMaximum(),
163  1,histslice->GetMean(),histslice->GetRMS());
164 
165  if (histslice->GetMaximum()) {
166 
167  histslice->Fit("gausBack");
168  amp[i] = gausBack->GetParameter(0);
169  sigfrac[i] = gausBack->GetParameter(1);
170  mean[i] = gausBack->GetParameter(2);
171  sigma[i] = gausBack->GetParameter(3);
172  chisq[i] = gausBack->GetChisquare();
173 
174  eamp[i] = gausBack->GetParError(0);
175  esigfrac[i] = gausBack->GetParError(1);
176  emean[i] = gausBack->GetParError(2);
177  esigma[i] = gausBack->GetParError(3);
178  echisq[i] = 0;
179  }
180  else {
181  amp[i] = 0;
182  sigfrac[i] = 0;
183  mean[i] = 0;
184  sigma[i] = 0;
185  chisq[i] = 0;
186 
187  eamp[i] = 0;
188  esigfrac[i] = 0;
189  emean[i] = 0;
190  esigma[i] = 0;
191  echisq[i] = 0;
192  }
193 
194  return 0;
195 
196 
197 
198 }
199 
200 TGraphErrors**
201 erfFit(TH1D* hist)
202 {
203  TGraphErrors** retval = funcFit(erfBack,hist);
204  return retval;
205 
206 }
207 TGraphErrors**
208 gausFit(TH1D* hist)
209 {
210  TGraphErrors** retval = funcFit(gausBack,hist);
211  return retval;
212 
213 }
214 
215 
216 
217 
218 
Generated by   doxygen 1.8.5