]>
Commit | Line | Data |
---|---|---|
8c1c76e9 | 1 | void DrawContamination(const char *fname = "HFEtask.root"){ |
2 | //TVirtualFitter::Set | |
3 | TFile *in = TFile::Open(fname); | |
4 | TList *qa = (TList *)in->Get("HFE_QA"); | |
5 | ||
6 | // Make Plots for TPC | |
7 | TList *pidqa = (TList *)qa->FindObject("PIDqa"); | |
8 | TList *tpc = (TList *)pidqa->FindObject("list_fQAhistosTPC"); | |
9 | // Create histograms by projecting the THnSparse | |
10 | THnSparseF *hTPC = (THnSparseF *)tpc->FindObject("fHistTPCsignal"); | |
11 | hTPC->GetAxis(4)->SetRange(1,1); | |
12 | TH2 *hTPCsig = hTPC->Projection(3,1); | |
13 | hTPCsig->SetName("hTPCsig"); | |
14 | ||
15 | // define cut model | |
16 | TF1 cutmodel("cutmodel", "[0] * TMath::Exp([1]*x) + [2]", 0, 20); | |
17 | cutmodel.SetParameter(0, -2.75); | |
18 | cutmodel.SetParameter(1, -0.8757); | |
19 | cutmodel.SetParameter(2, -0.9); | |
20 | ||
21 | Int_t offset = hTPCsig->GetXaxis()->FindBin(0.35); | |
22 | Int_t nPoints = hTPCsig->GetXaxis()->FindBin(4.) - offset; | |
23 | Int_t points = 0; | |
24 | TArrayD ap(nPoints), apErr(nPoints), acont(nPoints), aeff(nPoints), aemean(nPoints), aesigma(nPoints), aemeanErr(nPoints), aeSigmaErr(nPoints); | |
25 | Double_t p = 0., pmin = 0, pmax = 0, lowerbound = 0.; | |
26 | TH1 *hde = NULL; | |
27 | TCanvas outfit("outfit","Output of the fit", 640, 480); | |
28 | TCanvas singlefit("singlefit","Output of the fit", 640, 480); | |
29 | for(Int_t ipbin = hTPCsig->GetXaxis()->FindBin(0.5); ipbin < hTPCsig->GetXaxis()->FindBin(4.);ipbin+=12){ | |
30 | p = hTPCsig->GetXaxis()->GetBinCenter(ipbin); | |
31 | lowerbound = cutmodel.Eval(p); | |
32 | printf("ipbin = %d, p = %f, Lower limit %f\n", ipbin, p, lowerbound); | |
33 | hde = hTPCsig->ProjectionY("py", ipbin-6, ipbin+6); | |
34 | pmin = hTPCsig->GetXaxis()->GetBinLowEdge(ipbin-6); | |
35 | pmax = hTPCsig->GetXaxis()->GetBinUpEdge(ipbin+6); | |
36 | ||
37 | // First fit gaus for electrons and pions | |
38 | TF1 gpsingle("gpsingle", "[0]*TMath::Gaus(x, [1], [2])", -7, -1.5), | |
39 | gesingle("gpsingle", "[0]*TMath::Gaus(x, [1], [2])", -1.5, 3); | |
40 | ||
41 | gesingle.SetParLimits(0, 1, 1e9); gpsingle.SetParLimits(0, 1, 1e9); | |
42 | gesingle.SetParLimits(1, -0.7, 1.2); gpsingle.SetParLimits(1, -6, -4); | |
43 | gesingle.SetParLimits(2, 0.95, 1.35); gpsingle.SetParLimits(2, 0.75, 1.35); | |
44 | ||
45 | Double_t lowerlimitpi = -6, upperlimitpi = -2; | |
46 | if(p > 2.7) {lowerlimitpi = -5; upperlimitpi = -1; gpsingle.SetParLimits(-5, -2)} | |
47 | hde->Fit(&gpsingle, "", "N", lowerlimitpi, upperlimitpi); | |
48 | hde->Fit(&gesingle, "", "N", -1, 2); | |
49 | ||
50 | singlefit.cd(); | |
51 | singlefit.SetLogy(); | |
52 | hde->Draw(); | |
53 | gesingle.SetLineColor(kBlue); | |
54 | gesingle.Draw("lsame"); | |
55 | gpsingle.SetLineColor(kGreen); | |
56 | gpsingle.Draw("lsame"); | |
57 | singlefit.SaveAs(Form("singlegaus_input%dGeV.png", static_cast<Int_t>(p*100))); | |
58 | ||
59 | printf("Model restricted, fit with double gauss\n"); | |
60 | // Define 2 gaussian model: 1st gauss electron, 2nd gauss pion | |
61 | TF1 fitfun("2gausmodel", "[0]*TMath::Gaus(x, [1], [2]) + [3]*TMath::Gaus(x,[4],[5])", -20, 20); | |
62 | Double_t errmp = 0.1 * TMath::Abs(gpsingle.GetParameter(1)), | |
63 | errme = 0.2 * TMath::Abs(gesingle.GetParameter(1)), | |
64 | errsp = 0.1 * gpsingle.GetParameter(2), | |
65 | errse = 0.1 * gesingle.GetParameter(2); | |
66 | printf("Limits:\n=============================================\n"); | |
67 | printf("Electrons: Mean:\t\t%.3f, Sigma:\t\t%.3f\n", errme, errse); | |
68 | printf("Mean Lower:\t\t%.3f, Upper:\t\t%.3f\n", gesingle.GetParameter(1) - errme, gesingle.GetParameter(1) + errme); | |
69 | printf("Sigma Lower:\t\t%.3f, Upper:\t\t%.3f\n", gesingle.GetParameter(2) - errse, gesingle.GetParameter(2) + errse); | |
70 | printf("Pions: Mean:\t\t%.3f, Sigma:\t\t%.3f\n", errmp, errsp); | |
71 | printf("Mean Lower:\t\t%.3f, Upper:\t\t%.3f\n", gpsingle.GetParameter(1) - errmp, gpsingle.GetParameter(1) + errmp); | |
72 | printf("Mean Lower:\t\t%.3f, Upper:\t\t%.3f\n", gpsingle.GetParameter(2) - errsp, gpsingle.GetParameter(2) + errsp); | |
73 | ||
74 | fitfun.SetParLimits(0, 0.1, 1e9); | |
75 | fitfun.SetParLimits(3, 0.1, 1e9); | |
76 | fitfun.SetParLimits(1, gesingle.GetParameter(1) - errme, gesingle.GetParameter(1) + errme); | |
77 | fitfun.SetParLimits(4, gpsingle.GetParameter(1) - errmp, gpsingle.GetParameter(1) + errmp); | |
78 | fitfun.SetParLimits(2, gesingle.GetParameter(2) - errse, gesingle.GetParameter(2) + errse); | |
79 | fitfun.SetParLimits(5, gpsingle.GetParameter(2) - errsp, gpsingle.GetParameter(2) + errsp); | |
80 | fitfun.SetParameter(1, gesingle.GetParameter(1)); | |
81 | fitfun.SetParameter(4, gpsingle.GetParameter(1)); | |
82 | fitfun.SetParameter(2, gesingle.GetParameter(2)); | |
83 | fitfun.SetParameter(5, gpsingle.GetParameter(2)); | |
84 | // Fit the histogram | |
85 | hde->Fit(&fitfun, "","", lowerlimitpi, 3); | |
86 | // Save monitoring plot | |
87 | outfit.cd(); | |
88 | outfit.SetLogy(); | |
89 | outfit.Clear(); | |
90 | hde->SetTitle(); | |
91 | hde->GetYaxis()->SetTitle("Yield"); | |
92 | hde->SetStats(kFALSE); | |
93 | hde->Draw(); | |
94 | TPaveText plabel(0.1, 0.7, 0.3, 0.89, "NDC"); | |
95 | plabel.AddText(Form("p = %0.2fGeV/c", p)); | |
96 | plabel.SetBorderSize(0); | |
97 | plabel.SetFillStyle(0); | |
98 | plabel.Draw(); | |
99 | TPaveText fitpar(0.6, 0.5, 0.89, 0.89, "NDC"); fitpar.AddText("Fit Patameters"); | |
100 | fitpar.SetBorderSize(0); fitpar.SetFillStyle(0); | |
101 | for(Int_t ipar = 0; ipar < 6; ipar++) fitpar.AddText(Form("p%d: %0.3e #pm %0.3e", ipar, fitfun.GetParameter(ipar), fitfun.GetParError(ipar))); | |
102 | fitpar.Draw(); | |
103 | outfit.SaveAs(Form("contaminationFit%dGeV.png", static_cast<Int_t>(p*100))); | |
104 | ||
105 | // Make single gaus functions, integrate them within the electron selection band | |
106 | TF1 gausele("gausele", "[0]*TMath::Gaus(x, [1], [2])", -20, 20), | |
107 | gauspi("gausele", "[0]*TMath::Gaus(x, [1], [2])", -20, 20); | |
108 | for(Int_t ipar = 0; ipar < 3; ipar++){ | |
109 | gausele.SetParameter(ipar, fitfun.GetParameter(ipar)); | |
110 | gauspi.SetParameter(ipar, fitfun.GetParameter(ipar+3)); | |
111 | } | |
112 | Double_t cele = gausele.Integral(lowerbound, 2.); | |
113 | Double_t celt = gausele.Integral(-10, 10); | |
114 | Double_t cpi = gauspi.Integral(lowerbound, 2.); | |
115 | Double_t ct = fitfun.Integral(lowerbound, 2); | |
116 | Double_t cont = cpi / ct; | |
117 | Double_t eff = celt ? cele/celt : 0; | |
118 | printf("Contamination level at p = %0.2f: %f\n", p, cont); | |
119 | printf("Efficiency level at p = %0.2f: %f\n", p, eff); | |
120 | ||
121 | ap[points] = p; | |
122 | apErr[points] = (pmax - pmin)/2; | |
123 | acont[points] = cont; | |
124 | aeff[points] = eff; | |
125 | aemean[points] = fitfun.GetParameter(1); | |
126 | aesigma[points] = fitfun.GetParameter(2); | |
127 | aemeanErr[points] = fitfun.GetParError(1); | |
128 | aeSigmaErr[points] = fitfun.GetParError(2); | |
129 | points++; | |
130 | delete hde; | |
131 | } | |
132 | TGraph *contamination = new TGraph(points); | |
133 | TGraph *efficiency = new TGraph(points); | |
134 | TGraphErrors *electronMean = new TGraphErrors(points); | |
135 | TGraphErrors *electronSigma = new TGraphErrors(points); | |
136 | for(Int_t ip = 0; ip < points; ip++){ | |
137 | contamination->SetPoint(ip, ap[ip], acont[ip]); | |
138 | efficiency->SetPoint(ip, ap[ip], aeff[ip]); | |
139 | electronMean->SetPoint(ip, ap[ip], aemean[ip]); | |
140 | electronMean->SetPointError(ip, apErr[ip], aemeanErr[ip]); | |
141 | electronSigma->SetPoint(ip, ap[ip], aesigma[ip]); | |
142 | electronSigma->SetPointError(ip, apErr[ip], aeSigmaErr[ip]); | |
143 | } | |
144 | ||
145 | contamination->SetName("contamination"); | |
146 | efficiency->SetName("efficiency"); | |
147 | ||
148 | TFile *outfile = new TFile("Fitresults.root", "recreate"); | |
149 | outfile->cd(); | |
150 | efficiency->Write(); | |
151 | contamination->Write(); | |
152 | gROOT->cd(); | |
153 | ||
154 | // Draw Plots | |
155 | TCanvas *cCont = new TCanvas("cContamination", "Contamination", 800, 600); | |
156 | cCont->cd(); | |
157 | contamination->SetMarkerColor(kBlue); | |
158 | contamination->SetMarkerStyle(22); | |
159 | contamination->GetXaxis()->SetTitle("p_{T} / GeV/c"); | |
160 | contamination->GetYaxis()->SetTitle("contamination"); | |
161 | contamination->Draw("ap"); | |
162 | TCanvas *cEff = new TCanvas("cEfficiency", "Efficiency", 800, 600); | |
163 | cEff->cd(); | |
164 | efficiency->SetMarkerColor(kBlue); | |
165 | efficiency->SetMarkerStyle(22); | |
166 | efficiency->GetXaxis()->SetTitle("p_{T} / GeV/c"); | |
167 | efficiency->GetYaxis()->SetTitle("efficiency"); | |
168 | efficiency->Draw("ap"); | |
169 | outfile->Close(); | |
170 | TCanvas *cElectrons = new TCanvas("cElectrons", "Electron Mean and Sigma"); | |
171 | cElectrons->cd(); | |
172 | electronMean->SetMarkerStyle(22); | |
173 | electronMean->SetMarkerColor(kBlue); | |
174 | electronMean->SetLineColor(kBlue); | |
175 | electronMean->SetLineWidth(1); | |
176 | electronMean->GetXaxis()->SetTitle("p [GeV/c]"); | |
177 | electronMean->GetYaxis()->SetTitle("Electron Mean/Sigma [n#sigma]"); | |
178 | electronMean->GetYaxis()->SetRangeUser(-1., 2); | |
179 | electronSigma->SetMarkerStyle(22); | |
180 | electronSigma->SetMarkerColor(kRed); | |
181 | electronSigma->SetLineColor(kRed); | |
182 | electronSigma->SetLineWidth(1); | |
183 | electronSigma->GetXaxis()->SetTitle("p [GeV/c]"); | |
184 | electronSigma->GetYaxis()->SetTitle("Electron Mean/Sigma [n#sigma]"); | |
185 | electronSigma->GetYaxis()->SetRangeUser(-1., 2.); | |
186 | electronMean->Draw("ape"); | |
187 | electronSigma->Draw("epsame"); | |
188 | TLegend *leg = new TLegend(0.7, 0.75, 0.89, 0.89); | |
189 | leg->SetBorderSize(0); | |
190 | leg->SetFillStyle(0); | |
191 | leg->AddEntry(electronMean, "Mean", "p"); | |
192 | leg->AddEntry(electronSigma, "Sigma", "p"); | |
193 | leg->Draw(); | |
194 | delete outfile; | |
195 | } |