PWGHF/hfe/macros/DrawContamination.C

   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 }