]>
Commit | Line | Data |
---|---|---|
9d9613cc | 1 | #include "AliMassFitControl.h" |
2 | ||
3 | // Function for fitting back ground away from peak | |
4 | Float_t quad(Double_t *x, Double_t *par){ | |
5 | if (x[0] > par[3] && x[0] < par[4]){ | |
6 | TF1::RejectPoint(); | |
7 | return 0; | |
8 | } | |
9 | return par[0] + par[1]*x[0] + par[2]*x[0]*x[0]; //a+bx+cx**2 | |
10 | } | |
11 | ||
12 | TH1F* PtMassAna2(TH2F *PtMass, Int_t mode,Int_t ihist, const Int_t NControl, TObjArray *ControlArray,Char_t* fulllabel){ | |
13 | ||
14 | ||
15 | // gStyle->SetOptTitle(0); | |
16 | // gStyle->SetOptStat(0); | |
17 | gStyle->SetFrameBorderSize(0); | |
18 | gStyle->SetFrameBorderMode(0); | |
19 | gStyle->SetCanvasBorderMode(0); | |
20 | gStyle->SetCanvasColor(0); | |
21 | ||
22 | //TString *tLabel = new TString(fulllabel); //not reqd | |
23 | ||
24 | //Arguments TH2F *PtMass - 2D histogram of some variable versus mass | |
25 | // Int_t mode - mode for switching things depending on particle | |
26 | // const Int_t NControl - number of bins to project into which must | |
27 | // be the same as the number of controller objects (see next arg). | |
28 | // TObjArray *ContolArray - holds the controllers objects, 1 per bin | |
29 | ||
30 | //const Int_t NDefault=27; // Used in default array (change as required) | |
31 | // Trap user errors | |
32 | if (ControlArray->GetEntries() != NControl){ | |
33 | cout << "PtMassAna2: Problem: Specified number of bins and number of bins in supplied TObjArray do not match" << endl; | |
34 | return 0; | |
35 | } | |
36 | //Set up mode mapping to particle and names for histograms | |
37 | if(mode == 0){ | |
38 | const Char_t* part = "K0"; | |
39 | const Char_t* partName = "K^{0}_{s}"; | |
40 | } | |
41 | else if (mode == 1 ){ | |
42 | const Char_t* part = "LAM"; | |
43 | const Char_t* partName = "#Lambda"; | |
44 | } | |
45 | else if (mode == 2 ){//since lambda and anti-Lambda have same limits etc. | |
46 | const Char_t* part = "LAM"; | |
47 | const Char_t* partName = "#bar{#Lambda}"; | |
48 | } | |
49 | else if (mode ==3) {// This is for combined Lambda + anti-Lambda | |
50 | const Char_t* part = "LAM"; | |
51 | const Char_t* partName = "#Lambda + #bar{#Lambda}"; | |
52 | } | |
53 | else if (mode ==4) { | |
54 | const Char_t* part = "XI"; | |
55 | const Char_t* partName = "#Xi"; | |
56 | } | |
57 | else{ | |
58 | cout << "Mode not recognized, defaulting to K0" << endl; | |
59 | const Char_t* part = "K0"; | |
60 | } | |
61 | cout << "Debug info. Particle: " << part << endl; | |
62 | ||
63 | Float_t xxlo, xxhi; //Exclusion limits for fitting | |
64 | // FUTURE - these can also be part of AliMassFitControl to allow different regions | |
65 | // in different bins | |
66 | Float_t NSigmaEx=4.0; // Number of sigma to exclude from bkgd fit | |
67 | Float_t NSigmaInt=3.5; // Number of sigma to integrate histo over | |
68 | ||
69 | // Fitting function for initial combined peak + background fit | |
70 | TF1 *gausQuad = new TF1("gausQuad","(([1]/([2]*pow(6.2831,0.5)))*[6]*exp(-0.5*pow((x-[0])/[2],2)))+ [3] + [4]*x + [5]*x*x",0.3,1.8); | |
71 | //>>>>>>>>> NB Need to set proper range when it is used <<<<<<<<<<<<<<<<<<<<< | |
72 | //>>>>>>>>> AND use 'R' in fit option <<<<<<<<<<<<<<<<<<<<< | |
73 | gausQuad->SetLineWidth(1); | |
74 | gausQuad->SetLineColor(kRed); | |
75 | ||
76 | // Fitting function for background | |
77 | TF1* quad = new TF1("quad",quad,0.3,1.8,5); | |
78 | quad->SetLineWidth(2); | |
79 | quad->SetLineColor(kBlue); | |
80 | quad->SetLineStyle(2); | |
81 | ||
82 | // Function for background under peak - to be integrated and drawn | |
83 | TF1* qback = new TF1("qback","[0]+[1]*x+[2]*x*x",0.3,1.8); | |
84 | qback->SetLineWidth(2); | |
85 | qback->SetLineColor(kRed); | |
86 | qback->SetLineStyle(2); | |
87 | qback->SetFillStyle(3013); | |
88 | qback->SetFillColor(kRed); | |
89 | ||
90 | // Set ranges and limits according to particle | |
91 | Float_t xmin,xmax; | |
92 | Float_t defMean,defArea,defWidth,defa0,defa1,defa2,defBW; //defaults | |
93 | Float_t defMnLo, defMnHi; | |
94 | if (part=="K0"){ | |
95 | // FUTURE: xmin, xmax define the range considered in fit so could go into | |
96 | // AliMassFitControl | |
97 | //xmin=0.418; | |
98 | //xmax=0.578; | |
99 | // xxlo=0.46;xxhi=0.53; | |
100 | // gausQuad->SetParameters(0.4977,2000,0.003,1,1,0,0.001*RBFact); | |
101 | //gausQuad->SetParameters(0.4977,2000,0.003,1,1,0,0.001); | |
102 | defMean=0.4977; defArea=1000; defWidth=0.01; | |
103 | defa0=2; defa1=2; defa2=0; | |
104 | //gausQuad->SetParLimits(0,0.483,0.523); // Constrain mean | |
105 | defMnLo = 0.483; defMnHi=0.523; | |
106 | } | |
107 | if (part=="LAM"){ | |
108 | // FUTURE: See above | |
109 | //xmin=1.085; | |
110 | //xmax=1.17; | |
111 | // xxlo=1.10;xxhi=1.135; | |
112 | // gausQuad->SetParameters(1.115,2000,0.0028,10,100,-100,0.001*RBFact); | |
113 | // gausQuad->SetParameters(1.115,2000,0.0028,10,100,-100,0.001); | |
114 | defMean=1.115; defArea=1000; defWidth=0.08; | |
115 | defa0=100; defa1=-10; defa2=0; | |
116 | // gausQuad->SetParLimits(0,1.113,1.123); // Constrain mean | |
117 | defMnLo=1.113; defMnHi=1.123; | |
118 | } | |
119 | if (part=="XI"){ | |
120 | defMean=1.32171; defArea=500; defWidth=0.08; | |
121 | defa0=100; defa1=-10; defa2=0; | |
122 | defMnLo=1.301; defMnHi=1.341; | |
123 | } | |
124 | //gausQuad->SetRange(xmin,xmax); | |
125 | gausQuad->SetParLimits(1,0.0,1E7); // Constrain area 0 to 10 million | |
126 | gausQuad->SetParLimits(2,0.0001,0.04); // Constrain width 0.1 MeV to 40 MeV | |
127 | gausQuad->SetParNames("Mean","Area1","Sigma1","p0","p1","p2","BinWid"); | |
128 | ||
129 | // Deciding how to divide up the canvas | |
130 | c1->Clear(); | |
131 | c1->SetWindowSize(1024,768); | |
132 | Int_t NDraw = NControl; | |
133 | if(NDraw==16){c1->Divide(4,4);} | |
134 | elseif(NDraw<=36 && NDraw>25){ c1->Divide(6,6);} | |
135 | elseif(NDraw<=25 && NDraw>20){ c1->Divide(5,5);} | |
136 | elseif(NDraw<=20 && NDraw>16){ c1->Divide(5,4);} | |
137 | elseif(NDraw<=15 && NDraw>12){ c1->Divide(5,3);} | |
138 | elseif(NDraw<=12 && NDraw>8){ c1->Divide(4,3);} | |
139 | elseif(NDraw<=8 && NDraw>6){ c1->Divide(4,2);} | |
140 | elseif(NDraw<=6 && NDraw>4){c1->Divide(3,2);} | |
141 | elseif(NDraw<=4 && NDraw>2){c1->Divide(2,2);} | |
142 | elseif(NDraw==2){c1->Divide(2,1);} | |
143 | elseif(NDraw==1){/*do nothing*/;} | |
144 | else{c1->Divide(7,6);} | |
145 | ||
146 | // | |
147 | // Project out the histograms from 2D into 1D 'slices' | |
148 | // | |
149 | Char_t id[20], title[80]; | |
150 | //cout << "NControl: " << NControl << endl; | |
151 | TH1D* hMassSlice[NControl]; | |
152 | //Arrays to store various quantities which can later be histogrammed. | |
153 | // const Int_t NBinsArrays = NBins+1-NFirst; // array of 1D projection histograms is counted from 0 but other arrays are used for histograms contents and therefore N+1 are need because element zero is left empty | |
154 | const Int_t NBinsArrays = NControl+1; | |
155 | Stat_t sigBkgd[NBinsArrays], errSigBkgd[NBinsArrays]; // Signal+background and error on this | |
156 | Stat_t bkgd[NBinsArrays], errBkgd[NBinsArrays]; // Background and error on this | |
157 | Stat_t signal[NBinsArrays], errSignal[NBinsArrays]; // Net signal and error | |
158 | ||
159 | Stat_t chi2PerDOF1[NBinsArrays], errChi2PerDOF1[NBinsArrays]; // Chi-squared per defree of freedom from the initial fit | |
160 | Stat_t chi2PerDOF2[NBinsArrays], errChi2PerDOF2[NBinsArrays]; // Chi-squared per defree of freedom from the 2nd (background) fit | |
161 | ||
162 | Double_t BinCont = 0; | |
163 | Double_t BinErr = 0; | |
164 | ||
165 | ||
166 | ||
167 | // REDO won't need this? | |
168 | // Zero the signal - avoids problems with plotting | |
169 | for(Int_t j=0;j<NBinsArrays;j++){ | |
170 | signal[j]=0.0; | |
171 | errSignal[j]=1.0; | |
172 | } // <<< Not actually in use a the moment <<< | |
173 | ||
174 | Stat_t mean[NBinsArrays], errMean[NBinsArrays]; // Mean and error | |
175 | Stat_t sigma[NBinsArrays], errSigma[NBinsArrays]; // Gaussian width and error | |
176 | Stat_t area[NBinsArrays], errArea[NBinsArrays]; // Peak area from fit and error | |
177 | Stat_t integLow[NBinsArrays], integUp[NBinsArrays]; // Lower/upper limits for integration | |
178 | ||
179 | // Float_t DefaultBinPtEdges[NDefault]={0.0,0.1,0.2,0.3,0.4,0.5,0.6,0.8,1.0,1.2,1.4,1.6,1.8,2.0,2.2,2.4,2.6,2.8,3.0,3.5,4.0,4.5,5.0,5.5,6.0,6.5,7.0}; | |
180 | Float_t BinPtEdges[NBinsArrays]; | |
181 | ||
182 | // Float_t BinPtEdges[NBins+1]={0.0,0.2,0.4,0.6,0.8,1.0,1.2,1.4,1.6,1.8,2.0,2.2,2.4,2.6,2.8,3.0,3.5,4.0,4.5,5.0,5.5,6.0,6.5,7.0}; | |
183 | ||
184 | Int_t BinLo, BinHi; //For pt bins | |
185 | // **** Main loop over the AliMassFitControllers **** | |
186 | ||
187 | TIter controlIter(ControlArray); | |
188 | AliMassFitControl *controller; | |
189 | while (controller=(AliMassFitControl*)controlIter.Next()) { | |
190 | if(ihist == 0)controller->CalcBinLimits(20); //This BinsPerGeV argument should be calculated from the data | |
191 | if(ihist == 1)controller->CalcBinLimits(50); //This BinsPerGeV argument should be calculated from the data | |
192 | if(ihist == 2)controller->CalcBinLimits(50); //This BinsPerGeV argument should be calculated from the data | |
193 | if(ihist == 3)controller->CalcBinLimits(1); //This BinsPerGeV argument should be calculated from the data | |
194 | if(ihist == 4)controller->CalcBinLimits(1); //This BinsPerGeV argument should be calculated from the data | |
195 | if(ihist == 6)controller->CalcBinLimits(20000); //This BinsPerGeV argument should be calculated from the data | |
196 | if(ihist == 5)controller->CalcBinLimits(100); //This BinsPerGeV argument should be calculated from the data | |
197 | //Had to introduce this Nint fn otherwise got inconsistencies after type implicit conversion | |
198 | // BinLo=TMath::Nint(1.+BinPtEdges[N]*20.); //cout << "BinLo: " << BinLo << ", " << 1+BinPtEdges[N]*20. << endl; | |
199 | // BinHi=TMath::Nint(BinPtEdges[N+1]*20.); //cout << "BinHi: " << BinHi << ", " << BinPtEdges[N+1]*20. << endl; | |
200 | BinLo=controller->BinLower(); | |
201 | BinHi=controller->BinUpper(); | |
202 | Int_t N=ControlArray->IndexOf(controller) - 1; | |
203 | sprintf(id,"Mass%d",N); | |
204 | cout << "Mass histo:" << N << " Firstbin: " << BinLo << " Last bin:" << BinHi << " " << controller->PtLower() << "-" << controller->PtUpper() << " GeV" << endl; | |
205 | //cout << "About to create mass projection " << N << endl; | |
206 | if(ihist == 0) hMassSlice[N] = PtMass->ProjectionX(id,BinLo,BinHi); | |
207 | if(ihist != 0) hMassSlice[N] = PtMass->ProjectionY(id,BinLo,BinHi); | |
208 | //cout << "Mass projection " << N << " created." << endl; | |
209 | sprintf(title,"%s Mass, %.2f < p_{t} < %.2f",partName,controller->PtLower(),controller->PtUpper()); | |
210 | hMassSlice[N]->SetTitle(title); | |
211 | hMassSlice[N]->SetXTitle("GeV/c^{2}"); | |
212 | // cout << "Mass projection " << N << " titles set." << endl; | |
213 | hMassSlice[N]->Rebin(controller->RebinFactor()); | |
214 | // } // End of creating the projections | |
215 | Int_t massBinLo,massBinHi; //For mass histogram bins | |
216 | Int_t NArray; | |
217 | // Do all the fitting | |
218 | //REDO this is all one loop (iteration over TObjArray) now | |
219 | //for(Int_t N=0;N<NBins-NFirst;N++){ | |
220 | NArray=N+1; // Arrays numbered from 1 not 0.. | |
221 | c1->cd(N+1); //Canvas subdivisions numbered from 1 and NOT 0 | |
222 | ||
223 | // Set range - can be different each time | |
224 | xmin=controller->MinMass(); | |
225 | xmax=controller->MaxMass(); | |
226 | gausQuad->SetRange(xmin,xmax); | |
227 | // Everytime need to set some sensible parameters for both fits | |
228 | gausQuad->SetParameters(defMean,defArea,defWidth,defa0,defa1,defa2,0.1); | |
229 | quad->SetParameters(defa0,defa1,defa2,0.49,0.51);// Last 2 parameters are | |
230 | // the range to exclude but these are fixed later by FixParameter(3,.. etc. | |
231 | gausQuad->FixParameter(6,hMassSlice[N]->GetBinWidth(1)); | |
232 | // Release Linear parameter | |
233 | gausQuad->ReleaseParameter(4); | |
234 | quad->ReleaseParameter(1); | |
235 | // Release Quadratic parameter | |
236 | gausQuad->ReleaseParameter(5); | |
237 | quad->ReleaseParameter(2); | |
238 | gausQuad->SetParLimits(0,defMnLo,defMnHi); | |
239 | gausQuad->SetParLimits(1,0.0,1E7); // Constrain area 0 to 10 million | |
240 | gausQuad->SetParLimits(2,0.0001,0.04); // Constrain width 0.1 MeV to 40 MeV | |
241 | // gausQuad->SetParLimits(4,0,0); // Releases linear parameter | |
242 | // quad->SetParLimits(1,0,0); // Releases linear parameter | |
243 | // gausQuad->SetParLimits(5,0,0); // Releases quadratic parameter | |
244 | // quad->SetParLimits(2,0,0); // Releases quadratic parameter | |
245 | if(controller->FixedLin()){ | |
246 | quad->FixParameter(1,0.0); | |
247 | gausQuad->FixParameter(4,0.0); | |
248 | cout << "Info: linear part fixed to zero." << endl; | |
249 | gausQuad->SetParLimits(3,0.0,1000); // Ensure constant is +ve | |
250 | // NB documentation indicates that limits only used when option B | |
251 | // is specified. Also how are they can they be freed later? | |
252 | } | |
253 | if(controller->FixedQuad()){ | |
254 | quad->FixParameter(2,0.0); | |
255 | gausQuad->FixParameter(5,0.0); | |
256 | cout << "Info: quadratic part fixed to zero." << endl; | |
257 | } | |
258 | ||
259 | //setting the errors for histogram | |
260 | for(int i =1; i<= hMassSlice[N]->GetNbinsX();i++){ | |
261 | BinCont = hMassSlice[N]->GetBinContent(i); | |
262 | BinErr = TMath::Sqrt(BinCont); | |
263 | hMassSlice[N]->SetBinError(i,BinErr); | |
264 | } | |
265 | ||
266 | hMassSlice[N]->Fit(gausQuad,"LR"); | |
267 | hMassSlice[N]->Draw("E,SAME"); | |
268 | gausQuad->DrawCopy("SAME"); | |
269 | hMassSlice[N]->Clone("hMySlice"); | |
270 | //Fit the background: | |
271 | xxlo=gausQuad->GetParameter(0)-NSigmaEx*gausQuad->GetParameter(2); | |
272 | xxhi=gausQuad->GetParameter(0)+NSigmaEx*gausQuad->GetParameter(2); | |
273 | // Limits have to be adjusted to match bins | |
274 | massBinLo=hMassSlice[N]->FindBin(xxlo); | |
275 | xxlo = hMassSlice[N]->GetBinLowEdge(massBinLo); | |
276 | massBinHi=hMassSlice[N]->FindBin(xxhi)+1; | |
277 | xxhi = hMassSlice[N]->GetBinLowEdge(massBinHi); | |
278 | quad->SetParameters(gausQuad->GetParameter(3),gausQuad->GetParameter(4),gausQuad->GetParameter(5),xxlo,xxhi); | |
279 | quad->FixParameter(3,xxlo); | |
280 | quad->FixParameter(4,xxhi); | |
281 | quad->SetRange(xmin,xmax); | |
282 | hMassSlice[N]->Fit(quad,"LRN+"); | |
283 | quad->SetFillColor(30); | |
284 | quad->SetFillStyle(3013); | |
285 | quad->SetRange(xmin,xxlo); | |
286 | quad->DrawCopy("SAME"); | |
287 | quad->SetRange(xxhi,xmax); | |
288 | quad->DrawCopy("SAME"); | |
289 | ||
290 | // Draw the background - use a new function since other one not defined everywhere | |
291 | xxlo=gausQuad->GetParameter(0)-NSigmaInt*gausQuad->GetParameter(2); | |
292 | xxhi=gausQuad->GetParameter(0)+NSigmaInt*gausQuad->GetParameter(2); | |
293 | massBinLo=hMassSlice[N]->FindBin(xxlo); | |
294 | xxlo = hMassSlice[N]->GetBinLowEdge(massBinLo); | |
295 | massBinHi=hMassSlice[N]->FindBin(xxhi)+1; | |
296 | xxhi = hMassSlice[N]->GetBinLowEdge(massBinHi); | |
297 | qback->SetRange(xxlo,xxhi); | |
298 | qback->SetParameter(0,quad->GetParameter(0)); | |
299 | qback->SetParameter(1,quad->GetParameter(1)); | |
300 | qback->SetParameter(2,quad->GetParameter(2)); | |
301 | qback->DrawCopy("SAME"); | |
302 | if(N==3){ | |
303 | TCanvas *myCan = new TCanvas(); | |
304 | TPad *mypad = new TPad(); | |
305 | myCan->cd(); | |
306 | ||
307 | // TLatex *myTitle = new TLatex(1.04753,20527.5,"ALICE PbPb at #sqrt{s_{NN}} = 2.76TeV, |y|<0.75"); | |
308 | TLatex *myTitle = new TLatex(0.4753,2527.5,"ALICE PbPb at #sqrt{s_{NN}} = 2.76TeV, |y|<0.75"); | |
309 | myTitle->SetTextSize(0.04); | |
310 | myTitle->SetTextColor(1); | |
311 | ||
312 | // TLatex *myTitle1 = new TLatex(1.14634,17940.8, title); | |
313 | TLatex *myTitle1 = new TLatex(0.4634,2940.8, title); | |
314 | myTitle1->SetTextSize(0.04); | |
315 | myTitle1->SetTextColor(1); | |
316 | ||
317 | mypad->Draw(); | |
318 | ||
319 | ||
320 | hMassSlice[N]->Draw("E"); | |
321 | hMassSlice[N]->SetXTitle("Inv. Mass [GeV/c^{2}]"); | |
322 | hMassSlice[N]->SetYTitle("N"); | |
323 | myTitle->Draw("same"); | |
324 | myTitle1->Draw("same"); | |
325 | ||
326 | quad->SetRange(xmin,xxlo); | |
327 | quad->DrawCopy("SAME"); | |
328 | quad->SetRange(xxhi,xmax); | |
329 | quad->DrawCopy("SAME"); | |
330 | qback->SetRange(xxlo,xxhi); | |
331 | qback->DrawCopy("SAME"); | |
332 | TFile *hfile = new TFile("06_07MassFitPID.root","RECREATE"); | |
333 | hMassSlice[N]->Write(); | |
334 | hfile->Close(); | |
335 | ||
336 | } | |
337 | ||
338 | ||
339 | //Integrate the signal+background (i.e. original histo) | |
340 | sigBkgd[NArray]=hMassSlice[N]->Integral(massBinLo,massBinHi); | |
341 | errSigBkgd[NArray]=TMath::Sqrt(sigBkgd[NArray]); | |
342 | //Integrate background - better to do this analytically ? | |
343 | bkgd[NArray]=qback->Integral(xxlo,xxhi)/gausQuad->GetParameter(6);//Divide by bin width | |
344 | errBkgd[NArray]=TMath::Sqrt(bkgd[NArray]); //Get errors from fit parameters? | |
345 | ||
346 | // Fill arrays for diagnostic histograms | |
347 | mean[NArray]=gausQuad->GetParameter(0); | |
348 | errMean[NArray]=gausQuad->GetParError(0); | |
349 | area[NArray]=gausQuad->GetParameter(1); | |
350 | errArea[NArray]=gausQuad->GetParError(1); | |
351 | sigma[NArray]=gausQuad->GetParameter(2); | |
352 | errSigma[NArray]=gausQuad->GetParError(2); | |
353 | chi2PerDOF1[NArray]=gausQuad->GetChisquare()/(Double_t)gausQuad->GetNDF(); | |
354 | errChi2PerDOF1[NArray]= 0.0; // Don't know how to calc. error for this? | |
355 | chi2PerDOF2[NArray]=quad->GetChisquare()/(Double_t)quad->GetNDF(); | |
356 | errChi2PerDOF2[NArray]= 0.0; // Don't know how to calc. error for this? | |
357 | ||
358 | BinPtEdges[N]=controller->PtLower(); | |
359 | // BinPtEdges(N+1)=controller->PtUpper(); | |
360 | } | |
361 | BinPtEdges[NBinsArrays-1]=((AliMassFitControl*)ControlArray->Last())->PtUpper(); | |
362 | // for (Int_t jj=0;jj<NBinsArrays;jj++){ | |
363 | // cout << "BinPtEdges " << jj << " = " << BinPtEdges[jj] << endl; | |
364 | // cout << "Mean " << jj << " = " << mean[jj] << endl; | |
365 | // cout << "Sigma " << jj << " = " << sigma[jj] << endl; | |
366 | // cout << "Signal + bkgd " << jj << " = " << sigBkgd[jj] << endl; | |
367 | // cout << "Background " << jj << " = " << bkgd[jj] << endl; | |
368 | // } | |
369 | // Yields in each bin returned in a histogram | |
370 | TH1F* hYields= new TH1F("Yield","Raw Particle Yield",NBinsArrays-1,BinPtEdges); | |
371 | //hYields->Reset("ICE"); | |
372 | hYields->SetXTitle("p_{t} / (GeV/c)"); | |
373 | //hYields->SetContent(sigBkgd); | |
374 | //hYields->SetError(errSigBkgd); | |
375 | //c1->cd(4); | |
376 | //hYields->Draw(); | |
377 | ||
378 | // Fill the diagnostic histograms: clone, set title, contents, errors. | |
379 | TH1F* hMeans = hYields->Clone("Means"); | |
380 | hMeans->SetTitle("Mean from Gaussian Fit"); | |
381 | hMeans->SetContent(mean); | |
382 | hMeans->SetError(errMean); | |
383 | //hMeans->Sumw2(); | |
384 | TH1F* hSigmas = hYields->Clone("Sigmas"); | |
385 | hSigmas->SetTitle("Sigma from Gaussian Fit"); | |
386 | hSigmas->SetContent(sigma); | |
387 | hSigmas->SetError(errSigma); | |
388 | TH1F* hAreas = hYields->Clone("Areas"); | |
389 | hAreas->SetTitle("Area from Gaussian Fit"); | |
390 | hAreas->SetContent(area); | |
391 | hAreas->SetError(errArea); | |
392 | TH1F* hSigBkgd = hYields->Clone("SigBkgd"); | |
393 | hSigBkgd->SetTitle("Signal + Background (Histo. Integral)"); | |
394 | hSigBkgd->SetContent(sigBkgd); | |
395 | hSigBkgd->SetError(errSigBkgd); | |
396 | TH1F* hBkgd = hYields->Clone("Bkgd"); | |
397 | hBkgd->SetTitle("Background (Fn. integral)"); | |
398 | hBkgd->SetContent(bkgd); | |
399 | hBkgd->SetError(errBkgd); | |
400 | ||
401 | hYields->Sumw2(); | |
402 | hYields->Add(hSigBkgd,hBkgd,1.0,-1.0); | |
403 | ||
404 | TH1F* hSoverB = hYields->Clone("SoverB"); | |
405 | hSoverB->SetTitle("Signal to Background ratio"); | |
406 | hSoverB->Sumw2(); | |
407 | hSoverB->Divide(hBkgd); | |
408 | ||
409 | TH1F* hSoverA = hYields->Clone("SoverA"); // Signal over area | |
410 | hSoverA->SetTitle("Ratio of Signal to Area from Fit"); | |
411 | hSoverA->Sumw2(); | |
412 | hSoverA->Divide(hAreas); | |
413 | ||
414 | TH1F* hChi2PerDOF1 = hYields->Clone("Chi2PerDOF1"); | |
415 | hChi2PerDOF1->SetTitle("Chi-squared per d.o.f. - initial fit"); | |
416 | hChi2PerDOF1->SetContent(chi2PerDOF1); | |
417 | hChi2PerDOF1->SetError(errChi2PerDOF1); | |
418 | ||
419 | TH1F* hChi2PerDOF2 = hYields->Clone("Chi2PerDOF2"); | |
420 | hChi2PerDOF2->SetTitle("Chi-squared per d.o.f. - background fit"); | |
421 | hChi2PerDOF2->SetContent(chi2PerDOF2); | |
422 | hChi2PerDOF2->SetError(errChi2PerDOF2); | |
423 | ||
424 | // Draw the diagnostic histograms on their own canvas | |
425 | TCanvas* cDiag = new TCanvas("Diag","Diagnostics",600,600); | |
426 | cDiag->Divide(2,3); | |
427 | ||
428 | cDiag->cd(1); | |
429 | Float_t bookmass; | |
430 | if(part=="LAM"){ | |
431 | hMeans->SetMinimum(1.11); | |
432 | hMeans->SetMaximum(1.13); | |
433 | bookmass=1.11563; | |
434 | } else if (part=="K0"){ | |
435 | hMeans->SetMaximum(0.51); | |
436 | hMeans->SetMinimum(0.475); | |
437 | bookmass=0.4976; | |
438 | } else if (part=="XI") { | |
439 | hMeans->SetMaximum(1.34); | |
440 | hMeans->SetMinimum(1.30); | |
441 | bookmass=1.32171; | |
442 | } | |
443 | TLine PDGmass; | |
444 | PDGmass.SetLineStyle(2); | |
445 | hMeans->Draw(); | |
446 | PDGmass.DrawLine(0,bookmass,5,bookmass); | |
447 | ||
448 | cDiag->cd(2); | |
449 | if (part=="K0"){ | |
450 | //hSigmas->SetMaximum(0.012); | |
451 | hSigmas->SetMaximum(0.112); | |
452 | }else{ | |
453 | hSigmas->SetMaximum(0.006);} | |
454 | hSigmas->SetMinimum(0.0); | |
455 | hSigmas->Draw(); | |
456 | ||
457 | cDiag->cd(3); | |
458 | hSoverB->SetMinimum(0.0); | |
459 | hSoverB->Draw(); | |
460 | ||
461 | cDiag->cd(4); | |
462 | hSoverA->SetMinimum(0.5); | |
463 | hSoverA->SetMaximum(1.5); | |
464 | hSoverA->Draw(); | |
465 | ||
466 | cDiag->cd(5); | |
467 | hChi2PerDOF1->SetMinimum(0); | |
468 | hChi2PerDOF1->SetMaximum(6); | |
469 | hChi2PerDOF1->Draw(); | |
470 | ||
471 | cDiag->cd(6); | |
472 | hChi2PerDOF2->SetMinimum(0); | |
473 | hChi2PerDOF2->SetMaximum(60); | |
474 | hChi2PerDOF2->Draw(); | |
475 | ||
476 | Char_t fileNameBase[80]; | |
477 | sprintf(fileNameBase,"Diagnostics%s_%s",part,fulllabel); | |
478 | Char_t fileNamePng[80]; | |
479 | sprintf(fileNamePng,"%s.png",fileNameBase); | |
480 | Char_t fileNameEps[80]; | |
481 | sprintf(fileNameEps,"%s.eps",fileNameBase); | |
482 | Char_t fileNamePdf[80]; | |
483 | sprintf(fileNamePdf,"%s.pdf",fileNameBase); | |
484 | ||
485 | cDiag->SaveAs(fileNamePng); | |
486 | cDiag->SaveAs(fileNameEps); | |
487 | cDiag->SaveAs(fileNamePdf); | |
488 | ||
489 | // Scale result by the bin size to get dN/dpT and by bin centre to get 1/pt... | |
490 | Float_t y, ey; | |
491 | for(Int_t j=0;j<=hYields->GetNbinsX();j++){ | |
492 | y = hYields->GetBinContent(j)/hYields->GetBinWidth(j); | |
493 | ey = hYields->GetBinError(j)/hYields->GetBinWidth(j); | |
494 | hYields->SetBinContent(j,y); | |
495 | hYields->SetBinError(j,ey); | |
496 | // y = hYields->GetBinContent(j)/hYields->GetBinCenter(j); | |
497 | // ey = hYields->GetBinError(j)/hYields->GetBinCenter(j); | |
498 | hYields->SetBinContent(j,y); | |
499 | hYields->SetBinError(j,ey); | |
500 | } | |
501 | return hYields; | |
502 | ||
503 | } |