PWGHF/hfe/macros/DrawSlice.C

   1 void DrawSlice(Double_t p, Bool_t singlePiGauss){
   2   //TVirtualFitter::Set
   3   TFile *in = TFile::Open("HFEtask.root");
   4   TList *qa = (TList *)in->Get("HFE_QA");
   5
   6   // Make Plots for TPC
   7   TList *pidqa = (TList *)qa->FindObject("HFEpidQA");
   8   AliHFEtpcPIDqa *tpcqa = (AliHFEtpcPIDqa *)pidqa->FindObject("TPCQA");
   9   TH2 *hTPCsig = tpcqa->MakeSpectrumNSigma(AliHFEdetPIDqa::kBeforePID);
  10
  11   // define cut model
  12   TF1 cutmodel("cutmodel", "[0] * TMath::Exp([1]*x) + [2]", 0, 20);
  13   cutmodel.SetParameter(0, -2.75);
  14   cutmodel.SetParameter(1, -0.8757);
  15   cutmodel.SetParameter(2, -0.9);
  16
  17   Int_t nBinsP = 0;
  18   Int_t pbin = hTPCsig->GetXaxis()->FindBin(p);
  19   TH1 *hslice = hTPCsig->ProjectionY("hslice", pbin - nBinsP, pbin + nBinsP);
  20
  21   // Functions needed
  22   TF1 *pi1 = new TF1("pi1", "[0] * TMath::Gaus(x, [1], [2]) + [3]*TMath::Gaus(x, [4], [5])", -20, 20),
  23       *pi2 = new TF1("pi2", "[0] * TMath::Gaus(x, [1], [2]) + [3]*TMath::Gaus(x, [4], [5])", -20, 20),
  24       *el1 = new TF1("el1", "[0] * TMath::Gaus(x, [1], [2])", -20, 20),
  25       *el2 = new TF1("el1", "[0] * TMath::Gaus(x, [1], [2])", -20, 20),
  26       *combined = new TF1("combined",  "[0] * TMath::Gaus(x, [1], [2]) + [3]*TMath::Gaus(x, [4], [5]) +[6] * TMath::Gaus(x, [7], [8])" , -20, 20);
  27
  28   // Constraints
  29   Double_t minele = 1e0;
  30   Double_t maxele = 1e4;
  31   Double_t minpi = 1e2;
  32   Double_t maxpi = 1e5;
  33   pi1->SetParLimits(0, minpi, maxpi); pi1->SetParLimits(3, minpi, maxpi);
  34   pi1->SetParLimits(2, 0.4, 0.9); pi1->SetParLimits(5, 0.4, 0.9);
  35   pi1->SetParLimits(1, -6.5,-3.5); pi1->SetParLimits(4, -6.5,-3.);
  36   if(singlePiGauss){
  37     pi1->FixParameter(3,0);
  38     pi1->FixParameter(4,0);
  39     pi1->FixParameter(5,0);
  40   }
  41   el1->SetParLimits(0, minele, maxele);
  42   el1->SetParLimits(1, -0.4, 0.);
  43   el1->SetParLimits(2, 0.9, 1.3);
  44   el1->SetParameter(1, -0.1);
  45
  46   // test
  47   //el1->FixParameter(1, -3.78551e-01);
  48   hslice->Fit(pi1, "N", "", -8., -3.);
  49   hslice->Fit(el1, "LN", "", -1.1, 5);
  50
  51   // Use single fits to constrain combined fit
  52   combined->FixParameter(1, pi1->GetParameter(1));
  53   combined->FixParameter(2, pi1->GetParameter(2));
  54   combined->FixParameter(4, pi1->GetParameter(4));
  55   combined->FixParameter(5, pi1->GetParameter(5));
  56   combined->FixParameter(7, el1->GetParameter(1));
  57   combined->FixParameter(8, el1->GetParameter(2));
  58   if(singlePiGauss)
  59     combined->FixParameter(3, 0);
  60   else
  61     combined->SetParLimits(3, minpi, maxpi);
  62   combined->SetParLimits(0, minpi, maxpi);
  63   combined->SetParLimits(6, minele, maxele);
  64
  65   hslice->Fit(combined, "", "N", -8., 3.);
  66
  67   for(Int_t ipar = 0; ipar < 6; ipar++) pi2->SetParameter(ipar, combined->GetParameter(ipar));
  68   for(Int_t ipar = 0; ipar < 3; ipar++) el2->SetParameter(ipar, combined->GetParameter(ipar+6));
  69
  70   combined->SetLineColor(kBlack);
  71   el2->SetLineColor(kGreen);
  72   pi2->SetLineColor(kBlue);
  73   combined->SetLineWidth(2);
  74   el2->SetLineWidth(2);
  75   pi2->SetLineWidth(2);
  76   el2->SetLineStyle(2);
  77   pi2->SetLineStyle(2);
  78
  79   hslice->SetStats(kFALSE);
  80   hslice->SetTitle("TPC Signal");
  81   hslice->GetXaxis()->SetTitle("TPC dE/dx - <dE/dx>|_{electrons} [#sigma]");
  82
  83   TCanvas *output = new TCanvas("output", "Slice output");
  84   output->SetLogy();
  85   output->cd();
  86   hslice->Draw();
  87   combined->Draw("same");
  88   el2->Draw("same");
  89   pi2->Draw("same");
  90
  91   TLegend *leg = new TLegend(0.5, 0.7, 0.89, 0.89);
  92   leg->SetFillStyle(0);
  93   leg->SetBorderSize(0);
  94   leg->AddEntry(hslice, "Measurement", "l");
  95   leg->AddEntry(combined, "Electrons + Pions", "l");
  96   leg->AddEntry(el2, "Electrons", "l");
  97   leg->AddEntry(pi2, "Pions", "l");
  98   leg->Draw();
  99
 100   TPaveText *momentum = new TPaveText(0.6, 0.45, 0.8, 0.55, "NDC");
 101   momentum->SetBorderSize(0);
 102   momentum->SetFillStyle(0);
 103   momentum->AddText(Form("p = %.2fGeV/c", p));
 104   momentum->Draw();
 105
 106   // Draw also Selection Bands
 107   TF1 *fBandPions = new TF1(*pi2), *fBandElectrons = new TF1(*el2);
 108   fBandPions->SetRange(cutmodel.Eval(p), 5);
 109   fBandElectrons->SetRange(cutmodel.Eval(p), 5);
 110
 111   fBandPions->SetFillColor(kBlue);
 112   fBandPions->SetFillStyle(3004);
 113   fBandElectrons->SetFillColor(kGreen);
 114   fBandElectrons->SetFillStyle(3005);
 115   fBandPions->SetLineWidth(0);
 116   fBandElectrons->SetLineWidth(0);
 117
 118   fBandElectrons->Draw("same");
 119   fBandPions->Draw("same");
 120
 121   // Now calculate the contamination
 122   Double_t nPions = pi2->Integral(cutmodel.Eval(p), 5);
 123   Double_t nCandidates = combined->Integral(cutmodel.Eval(p), 5);
 124   Double_t contamination = nPions/nCandidates;
 125   Double_t nCandidatesHisto = hslice->Integral(hslice->GetXaxis()->FindBin(cutmodel.Eval(p)), hslice->GetXaxis()->FindBin(5));
 126   Double_t nPionsBefore = pi1->Integral(cutmodel.Eval(p), 5);
 127   Double_t contaminationHisto = nPionsBefore/nCandidatesHisto;
 128   printf("Contamination @ %.3fGeV/c: %f\n", p, contamination);
 129   printf("Contamination @ %.3fGeV/c: %f (determined from double gauss and histogram)\n", p, contaminationHisto);
 130
 131   TPaveText *tcont = new TPaveText(0.45, 0.35, 0.9, 0.45, "NDC");
 132   tcont->SetBorderSize(0);
 133   tcont->SetFillStyle(0);
 134   tcont->AddText(Form("contamination: %f", contamination));
 135   tcont->Draw();
 136
 137   // Save Canvas
 138   TString filename = Form("slice%d.root", ((Int_t)(100*p)));
 139   if(singlePiGauss)
 140     filename = Form("slice%dsingle.root", ((Int_t)(100*p)));
 141
 142   TFile *out = new TFile(filename.Data(), "RECREATE");
 143   out->cd();
 144   output->Write();
 145   out->Close(); delete out;
 146 }