1 #include "AliMassFitControl.h"
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]){
9 return par[0] + par[1]*x[0] + par[2]*x[0]*x[0]; //a+bx+cx**2
12 TH1F* PtMassAna2(TH2F *PtMass, Int_t mode,Int_t ihist, const Int_t NControl, TObjArray *ControlArray,Char_t* fulllabel){
14 gROOT->SetStyle("Plain");
15 gStyle->SetOptFit(1111);
16 //TString *tLabel = new TString(fulllabel); //not reqd
18 //Arguments TH2F *PtMass - 2D histogram of some variable versus mass
19 // Int_t mode - mode for switching things depending on particle
20 // const Int_t NControl - number of bins to project into which must
21 // be the same as the number of controller objects (see next arg).
22 // TObjArray *ContolArray - holds the controllers objects, 1 per bin
24 //const Int_t NDefault=27; // Used in default array (change as required)
26 if (ControlArray->GetEntries() != NControl){
27 cout << "PtMassAna2: Problem: Specified number of bins and number of bins in supplied TObjArray do not match" << endl;
30 //Set up mode mapping to particle and names for histograms
32 const Char_t* part = "K0";
33 const Char_t* partName = "K^{0}_{s}";
36 const Char_t* part = "LAM";
37 const Char_t* partName = "#Lambda";
39 else if (mode == 2 ){//since lambda and anti-Lambda have same limits etc.
40 const Char_t* part = "LAM";
41 const Char_t* partName = "#bar{#Lambda}";
43 else if (mode ==3) {// This is for combined Lambda + anti-Lambda
44 const Char_t* part = "LAM";
45 const Char_t* partName = "#Lambda + #bar{#Lambda}";
48 const Char_t* part = "XI";
49 const Char_t* partName = "#Xi";
52 cout << "Mode not recognized, defaulting to K0" << endl;
53 const Char_t* part = "K0";
55 cout << "Debug info. Particle: " << part << endl;
57 Float_t xxlo, xxhi; //Exclusion limits for fitting
58 // FUTURE - these can also be part of AliMassFitControl to allow different regions
60 Float_t NSigmaEx=4.0; // Number of sigma to exclude from bkgd fit
61 Float_t NSigmaInt=3.5; // Number of sigma to integrate histo over
63 // Fitting function for initial combined peak + background fit
64 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);
65 //>>>>>>>>> NB Need to set proper range when it is used <<<<<<<<<<<<<<<<<<<<<
66 //>>>>>>>>> AND use 'R' in fit option <<<<<<<<<<<<<<<<<<<<<
67 gausQuad->SetLineWidth(1);
68 gausQuad->SetLineColor(kRed);
70 // Fitting function for background
71 TF1* quad = new TF1("quad",quad,0.3,1.8,5);
72 quad->SetLineWidth(2);
73 quad->SetLineColor(kBlue);
74 quad->SetLineStyle(2);
76 // Function for background under peak - to be integrated and drawn
77 TF1* qback = new TF1("qback","[0]+[1]*x+[2]*x*x",0.3,1.8);
78 qback->SetLineWidth(2);
79 qback->SetLineColor(kRed);
80 qback->SetLineStyle(2);
81 qback->SetFillStyle(3013);
82 qback->SetFillColor(kRed);
84 // Set ranges and limits according to particle
86 Float_t defMean,defArea,defWidth,defa0,defa1,defa2,defBW; //defaults
87 Float_t defMnLo, defMnHi;
89 // FUTURE: xmin, xmax define the range considered in fit so could go into
93 // xxlo=0.46;xxhi=0.53;
94 // gausQuad->SetParameters(0.4977,2000,0.003,1,1,0,0.001*RBFact);
95 //gausQuad->SetParameters(0.4977,2000,0.003,1,1,0,0.001);
96 defMean=0.4977; defArea=1000; defWidth=0.01;
97 defa0=2; defa1=2; defa2=0;
98 //gausQuad->SetParLimits(0,0.483,0.523); // Constrain mean
99 defMnLo = 0.483; defMnHi=0.523;
105 // xxlo=1.10;xxhi=1.135;
106 // gausQuad->SetParameters(1.115,2000,0.0028,10,100,-100,0.001*RBFact);
107 // gausQuad->SetParameters(1.115,2000,0.0028,10,100,-100,0.001);
108 defMean=1.115; defArea=1000; defWidth=0.08;
109 defa0=100; defa1=-10; defa2=0;
110 // gausQuad->SetParLimits(0,1.113,1.123); // Constrain mean
111 defMnLo=1.113; defMnHi=1.123;
114 defMean=1.32171; defArea=500; defWidth=0.08;
115 defa0=100; defa1=-10; defa2=0;
116 defMnLo=1.301; defMnHi=1.341;
118 //gausQuad->SetRange(xmin,xmax);
119 gausQuad->SetParLimits(1,0.0,1E7); // Constrain area 0 to 10 million
120 gausQuad->SetParLimits(2,0.0001,0.04); // Constrain width 0.1 MeV to 40 MeV
121 gausQuad->SetParNames("Mean","Area1","Sigma1","p0","p1","p2","BinWid");
123 // Deciding how to divide up the canvas
125 c1->SetWindowSize(1024,768);
126 Int_t NDraw = NControl;
127 if(NDraw==16){c1->Divide(4,4);}
128 elseif(NDraw<=36 && NDraw>30){ c1->Divide(6,6);}
129 elseif(NDraw<=30 && NDraw>25) { c1->Divide(6,5);}
130 elseif(NDraw<=25 && NDraw>20){ c1->Divide(5,5);}
131 elseif(NDraw<=20 && NDraw>16){ c1->Divide(5,4);}
132 elseif(NDraw<=15 && NDraw>12){ c1->Divide(5,3);}
133 elseif(NDraw<=12 && NDraw>8){ c1->Divide(4,3);}
134 elseif(NDraw<=8 && NDraw>6){ c1->Divide(4,2);}
135 elseif(NDraw<=6 && NDraw>4){c1->Divide(3,2);}
136 elseif(NDraw<=4 && NDraw>2){c1->Divide(2,2);}
137 elseif(NDraw==2){c1->Divide(2,1);}
138 elseif(NDraw==1){/*do nothing*/;}
139 else{c1->Divide(7,6);}
142 // Project out the histograms from 2D into 1D 'slices'
144 Char_t id[20], title[80];
145 //cout << "NControl: " << NControl << endl;
146 TH1D* hMassSlice[NControl];
147 //Arrays to store various quantities which can later be histogrammed.
148 // 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
149 const Int_t NBinsArrays = NControl+1;
150 Stat_t sigBkgd[NBinsArrays], errSigBkgd[NBinsArrays]; // Signal+background and error on this
151 Stat_t bkgd[NBinsArrays], errBkgd[NBinsArrays]; // Background and error on this
152 Stat_t signal[NBinsArrays], errSignal[NBinsArrays]; // Net signal and error
154 Stat_t chi2PerDOF1[NBinsArrays], errChi2PerDOF1[NBinsArrays]; // Chi-squared per defree of freedom from the initial fit
155 Stat_t chi2PerDOF2[NBinsArrays], errChi2PerDOF2[NBinsArrays]; // Chi-squared per defree of freedom from the 2nd (background) fit
157 // REDO won't need this?
158 // Zero the signal - avoids problems with plotting
159 for(Int_t j=0;j<NBinsArrays;j++){
162 } // <<< Not actually in use a the moment <<<
164 Stat_t mean[NBinsArrays], errMean[NBinsArrays]; // Mean and error
165 Stat_t sigma[NBinsArrays], errSigma[NBinsArrays]; // Gaussian width and error
166 Stat_t area[NBinsArrays], errArea[NBinsArrays]; // Peak area from fit and error
167 Stat_t integLow[NBinsArrays], integUp[NBinsArrays]; // Lower/upper limits for integration
169 // 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};
170 Float_t BinPtEdges[NBinsArrays];
172 // 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};
174 Int_t BinLo, BinHi; //For pt bins
175 // **** Main loop over the AliMassFitControllers ****
177 TIter controlIter(ControlArray);
178 AliMassFitControl *controller;
179 while (controller=(AliMassFitControl*)controlIter.Next()) {
180 if(ihist == 0)controller->CalcBinLimits(20); //This BinsPerGeV argument should be calculated from the data
181 if(ihist == 1)controller->CalcBinLimits(50); //This BinsPerGeV argument should be calculated from the data
182 if(ihist == 2)controller->CalcBinLimits(50); //This BinsPerGeV argument should be calculated from the data
183 if(ihist == 3)controller->CalcBinLimits(1); //This BinsPerGeV argument should be calculated from the data
184 if(ihist == 4)controller->CalcBinLimits(1); //This BinsPerGeV argument should be calculated from the data
185 if(ihist == 6)controller->CalcBinLimits(20000); //This BinsPerGeV argument should be calculated from the data
186 if(ihist == 5)controller->CalcBinLimits(100); //This BinsPerGeV argument should be calculated from the data
187 //Had to introduce this Nint fn otherwise got inconsistencies after type implicit conversion
188 // BinLo=TMath::Nint(1.+BinPtEdges[N]*20.); //cout << "BinLo: " << BinLo << ", " << 1+BinPtEdges[N]*20. << endl;
189 // BinHi=TMath::Nint(BinPtEdges[N+1]*20.); //cout << "BinHi: " << BinHi << ", " << BinPtEdges[N+1]*20. << endl;
190 BinLo=controller->BinLower();
191 BinHi=controller->BinUpper();
192 Int_t N=ControlArray->IndexOf(controller) - 1;
193 sprintf(id,"Mass%d",N);
194 cout << "Mass histo:" << N << " Firstbin: " << BinLo << " Last bin:" << BinHi << " " << controller->PtLower() << "-" << controller->PtUpper() << " GeV" << endl;
195 //cout << "About to create mass projection " << N << endl;
196 if(ihist == 0) hMassSlice[N] = PtMass->ProjectionX(id,BinLo,BinHi);
197 if(ihist != 0) hMassSlice[N] = PtMass->ProjectionY(id,BinLo,BinHi);
198 //cout << "Mass projection " << N << " created." << endl;
199 sprintf(title,"%s Mass, %.2f < p_{t} < %.2f",partName,controller->PtLower(),controller->PtUpper());
200 hMassSlice[N]->SetTitle(title);
201 hMassSlice[N]->SetXTitle("GeV/c^{2}");
202 // cout << "Mass projection " << N << " titles set." << endl;
203 hMassSlice[N]->Rebin(controller->RebinFactor());
204 // } // End of creating the projections
205 Int_t massBinLo,massBinHi; //For mass histogram bins
207 // Do all the fitting
208 //REDO this is all one loop (iteration over TObjArray) now
209 //for(Int_t N=0;N<NBins-NFirst;N++){
210 NArray=N+1; // Arrays numbered from 1 not 0..
211 c1->cd(N+1); //Canvas subdivisions numbered from 1 and NOT 0
213 // Set range - can be different each time
214 xmin=controller->MinMass();
215 xmax=controller->MaxMass();
216 gausQuad->SetRange(xmin,xmax);
217 // Everytime need to set some sensible parameters for both fits
218 gausQuad->SetParameters(defMean,defArea,defWidth,defa0,defa1,defa2,0.1);
219 quad->SetParameters(defa0,defa1,defa2,0.49,0.51);// Last 2 parameters are
220 // the range to exclude but these are fixed later by FixParameter(3,.. etc.
221 gausQuad->FixParameter(6,hMassSlice[N]->GetBinWidth(1));
222 // Release Linear parameter
223 gausQuad->ReleaseParameter(4);
224 quad->ReleaseParameter(1);
225 // Release Quadratic parameter
226 gausQuad->ReleaseParameter(5);
227 quad->ReleaseParameter(2);
228 gausQuad->SetParLimits(0,defMnLo,defMnHi);
229 gausQuad->SetParLimits(1,0.0,1E7); // Constrain area 0 to 10 million
230 gausQuad->SetParLimits(2,0.0001,0.04); // Constrain width 0.1 MeV to 40 MeV
231 // gausQuad->SetParLimits(4,0,0); // Releases linear parameter
232 // quad->SetParLimits(1,0,0); // Releases linear parameter
233 // gausQuad->SetParLimits(5,0,0); // Releases quadratic parameter
234 // quad->SetParLimits(2,0,0); // Releases quadratic parameter
235 if(controller->FixedLin()){
236 quad->FixParameter(1,0.0);
237 gausQuad->FixParameter(4,0.0);
238 cout << "Info: linear part fixed to zero." << endl;
239 gausQuad->SetParLimits(3,0.0,1000); // Ensure constant is +ve
240 // NB documentation indicates that limits only used when option B
241 // is specified. Also how are they can they be freed later?
243 if(controller->FixedQuad()){
244 quad->FixParameter(2,0.0);
245 gausQuad->FixParameter(5,0.0);
246 cout << "Info: quadratic part fixed to zero." << endl;
248 hMassSlice[N]->Fit(gausQuad,"LR");
249 hMassSlice[N]->Draw("E,SAME");
250 gausQuad->DrawCopy("SAME");
251 hMassSlice[N]->Clone("hMySlice");
252 //Fit the background:
253 xxlo=gausQuad->GetParameter(0)-NSigmaEx*gausQuad->GetParameter(2);
254 xxhi=gausQuad->GetParameter(0)+NSigmaEx*gausQuad->GetParameter(2);
255 // Limits have to be adjusted to match bins
256 massBinLo=hMassSlice[N]->FindBin(xxlo);
257 xxlo = hMassSlice[N]->GetBinLowEdge(massBinLo);
258 massBinHi=hMassSlice[N]->FindBin(xxhi)+1;
259 xxhi = hMassSlice[N]->GetBinLowEdge(massBinHi);
260 quad->SetParameters(gausQuad->GetParameter(3),gausQuad->GetParameter(4),gausQuad->GetParameter(5),xxlo,xxhi);
261 quad->FixParameter(3,xxlo);
262 quad->FixParameter(4,xxhi);
263 quad->SetRange(xmin,xmax);
264 hMassSlice[N]->Fit(quad,"LRN+");
265 quad->SetFillColor(30);
266 quad->SetFillStyle(3013);
267 quad->SetRange(xmin,xxlo);
268 quad->DrawCopy("SAME");
269 quad->SetRange(xxhi,xmax);
270 quad->DrawCopy("SAME");
272 // Draw the background - use a new function since other one not defined everywhere
273 xxlo=gausQuad->GetParameter(0)-NSigmaInt*gausQuad->GetParameter(2);
274 xxhi=gausQuad->GetParameter(0)+NSigmaInt*gausQuad->GetParameter(2);
275 massBinLo=hMassSlice[N]->FindBin(xxlo);
276 xxlo = hMassSlice[N]->GetBinLowEdge(massBinLo);
277 massBinHi=hMassSlice[N]->FindBin(xxhi)+1;
278 xxhi = hMassSlice[N]->GetBinLowEdge(massBinHi);
279 qback->SetRange(xxlo,xxhi);
280 qback->SetParameter(0,quad->GetParameter(0));
281 qback->SetParameter(1,quad->GetParameter(1));
282 qback->SetParameter(2,quad->GetParameter(2));
283 qback->DrawCopy("SAME");
285 TCanvas *myCan = new TCanvas();
286 TPad *mypad = new TPad();
289 hMassSlice[N]->Draw();
290 quad->SetRange(xmin,xxlo);
291 quad->DrawCopy("SAME");
292 quad->SetRange(xxhi,xmax);
293 quad->DrawCopy("SAME");
294 qback->DrawCopy("SAME");
299 //Integrate the signal+background (i.e. original histo)
300 sigBkgd[NArray]=hMassSlice[N]->Integral(massBinLo,massBinHi);
301 errSigBkgd[NArray]=TMath::Sqrt(sigBkgd[NArray]);
302 //Integrate background - better to do this analytically ?
303 bkgd[NArray]=qback->Integral(xxlo,xxhi)/gausQuad->GetParameter(6);//Divide by bin width
304 errBkgd[NArray]=TMath::Sqrt(bkgd[NArray]); //Get errors from fit parameters?
306 // Fill arrays for diagnostic histograms
307 mean[NArray]=gausQuad->GetParameter(0);
308 errMean[NArray]=gausQuad->GetParError(0);
309 area[NArray]=gausQuad->GetParameter(1);
310 errArea[NArray]=gausQuad->GetParError(1);
311 sigma[NArray]=gausQuad->GetParameter(2);
312 errSigma[NArray]=gausQuad->GetParError(2);
313 chi2PerDOF1[NArray]=gausQuad->GetChisquare()/(Double_t)gausQuad->GetNDF();
314 errChi2PerDOF1[NArray]= 0.0; // Don't know how to calc. error for this?
315 chi2PerDOF2[NArray]=quad->GetChisquare()/(Double_t)quad->GetNDF();
316 errChi2PerDOF2[NArray]= 0.0; // Don't know how to calc. error for this?
318 BinPtEdges[N]=controller->PtLower();
319 // BinPtEdges(N+1)=controller->PtUpper();
321 BinPtEdges[NBinsArrays-1]=((AliMassFitControl*)ControlArray->Last())->PtUpper();
322 // for (Int_t jj=0;jj<NBinsArrays;jj++){
323 // cout << "BinPtEdges " << jj << " = " << BinPtEdges[jj] << endl;
324 // cout << "Mean " << jj << " = " << mean[jj] << endl;
325 // cout << "Sigma " << jj << " = " << sigma[jj] << endl;
326 // cout << "Signal + bkgd " << jj << " = " << sigBkgd[jj] << endl;
327 // cout << "Background " << jj << " = " << bkgd[jj] << endl;
329 // Yields in each bin returned in a histogram
330 TH1F* hYields= new TH1F("Yield","Raw Particle Yield",NBinsArrays-1,BinPtEdges);
331 //hYields->Reset("ICE");
332 hYields->SetXTitle("p_{t} / (GeV/c)");
333 //hYields->SetContent(sigBkgd);
334 //hYields->SetError(errSigBkgd);
338 // Fill the diagnostic histograms: clone, set title, contents, errors.
339 TH1F* hMeans = hYields->Clone("Means");
340 hMeans->SetTitle("Mean from Gaussian Fit");
341 hMeans->SetContent(mean);
342 hMeans->SetError(errMean);
344 TH1F* hSigmas = hYields->Clone("Sigmas");
345 hSigmas->SetTitle("Sigma from Gaussian Fit");
346 hSigmas->SetContent(sigma);
347 hSigmas->SetError(errSigma);
348 TH1F* hAreas = hYields->Clone("Areas");
349 hAreas->SetTitle("Area from Gaussian Fit");
350 hAreas->SetContent(area);
351 hAreas->SetError(errArea);
352 TH1F* hSigBkgd = hYields->Clone("SigBkgd");
353 hSigBkgd->SetTitle("Signal + Background (Histo. Integral)");
354 hSigBkgd->SetContent(sigBkgd);
355 hSigBkgd->SetError(errSigBkgd);
356 TH1F* hBkgd = hYields->Clone("Bkgd");
357 hBkgd->SetTitle("Background (Fn. integral)");
358 hBkgd->SetContent(bkgd);
359 hBkgd->SetError(errBkgd);
362 hYields->Add(hSigBkgd,hBkgd,1.0,-1.0);
364 TH1F* hSoverB = hYields->Clone("SoverB");
365 hSoverB->SetTitle("Signal to Background ratio");
367 hSoverB->Divide(hBkgd);
369 TH1F* hSoverA = hYields->Clone("SoverA"); // Signal over area
370 hSoverA->SetTitle("Ratio of Signal to Area from Fit");
372 hSoverA->Divide(hAreas);
374 TH1F* hChi2PerDOF1 = hYields->Clone("Chi2PerDOF1");
375 hChi2PerDOF1->SetTitle("Chi-squared per d.o.f. - initial fit");
376 hChi2PerDOF1->SetContent(chi2PerDOF1);
377 hChi2PerDOF1->SetError(errChi2PerDOF1);
379 TH1F* hChi2PerDOF2 = hYields->Clone("Chi2PerDOF2");
380 hChi2PerDOF2->SetTitle("Chi-squared per d.o.f. - background fit");
381 hChi2PerDOF2->SetContent(chi2PerDOF2);
382 hChi2PerDOF2->SetError(errChi2PerDOF2);
384 // Draw the diagnostic histograms on their own canvas
385 TCanvas* cDiag = new TCanvas("Diag","Diagnostics",600,600);
392 hMeans->SetMinimum(bookmass-0.01);
393 hMeans->SetMaximum(bookmass+0.01);
394 } else if (part=="K0"){
396 hMeans->SetMinimum(bookmass-0.01);
397 hMeans->SetMaximum(bookmass+0.01);
398 } else if (part=="XI") {
399 hMeans->SetMaximum(1.34);
400 hMeans->SetMinimum(1.30);
403 Float_t maxPt = hMeans->GetBinLowEdge(hMeans->GetNbinsX()+1);
405 PDGmass.SetLineStyle(2);
407 PDGmass.DrawLine(0,bookmass,maxPt,bookmass);
411 hSigmas->SetMaximum(0.06);
413 hSigmas->SetMaximum(0.006);}
414 hSigmas->SetMinimum(0.0);
418 hSoverB->SetMinimum(0.0);
422 hSoverA->SetMinimum(0.5);
423 hSoverA->SetMaximum(1.5);
427 hChi2PerDOF1->SetMinimum(0);
428 hChi2PerDOF1->SetMaximum(6);
429 hChi2PerDOF1->Draw();
432 hChi2PerDOF2->SetMinimum(0);
433 hChi2PerDOF2->SetMaximum(6);
434 hChi2PerDOF2->Draw();
436 Char_t fileNameBase[80];
437 sprintf(fileNameBase,"Diagnostics%s_%s",part,fulllabel);
438 Char_t fileNamePng[80];
439 sprintf(fileNamePng,"%s.png",fileNameBase);
440 Char_t fileNameEps[80];
441 sprintf(fileNameEps,"%s.eps",fileNameBase);
442 Char_t fileNamePdf[80];
443 sprintf(fileNamePdf,"%s.pdf",fileNameBase);
445 cDiag->SaveAs(fileNamePng);
446 cDiag->SaveAs(fileNameEps);
447 cDiag->SaveAs(fileNamePdf);
449 // Scale result by the bin size to get dN/dpT and by bin centre to get 1/pt...
451 for(Int_t j=0;j<=hYields->GetNbinsX();j++){
452 y = hYields->GetBinContent(j)/hYields->GetBinWidth(j);
453 ey = hYields->GetBinError(j)/hYields->GetBinWidth(j);
454 hYields->SetBinContent(j,y);
455 hYields->SetBinError(j,ey);
456 // y = hYields->GetBinContent(j)/hYields->GetBinCenter(j);
457 // ey = hYields->GetBinError(j)/hYields->GetBinCenter(j);
458 hYields->SetBinContent(j,y);
459 hYields->SetBinError(j,ey);