PWG3/hfe/macros/DrawSliceProton.C

   1 void DrawSliceProton(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 = 2;
  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]) + [6]*TMath::Gaus(x, [7], [8])", -20, 20),
  23       *pi2 = new TF1("pi2", "[0] * TMath::Gaus(x, [1], [2]) + [3]*TMath::Gaus(x, [4], [5]) + [6]*TMath::Gaus(x, [7], [8])", -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]) + [9]*TMath::Gaus(x, [10], [11])" , -20, 20);
  27
  28   // Constraints
  29   Double_t minele = 1e0;
  30   Double_t maxele = 1e3;
  31   Double_t minpr = 1e1;
  32   Double_t maxpr = 1e4;
  33   Double_t minpi = 1e2;
  34   Double_t maxpi = 1e5;
  35   pi1->SetParLimits(0, minpi, maxpi); pi1->SetParLimits(3, minpi, maxpi);
  36   pi1->SetParLimits(2, 0.4, 1.5); pi1->SetParLimits(5, 0.4, 1.5);
  37   pi1->SetParLimits(1, -8,-4); pi1->SetParLimits(4, -5.5,-3);
  38   if(singlePiGauss){
  39     pi1->FixParameter(3,0);
  40     pi1->FixParameter(4,0);
  41     pi1->FixParameter(5,0);
  42   }
  43   // Set Protons
  44   pi1->SetParLimits(6, minpr, maxpr);
  45   pi1->SetParLimits(7, -7., -5.8);
  46   pi1->SetParLimits(8, 0.1, 1.3);
  47   // Set Electrons
  48   el1->SetParLimits(0, minele, maxele);
  49   el1->SetParLimits(1, -0.9., 0.);
  50   el1->SetParLimits(2, 0.9, 1.5);
  51   el1->SetParameter(1, -0.1);
  52   //el1->FixParameter(1, -1.34726e-01);
  53   //el1->FixParameter(2, 9.77243e-01);
  54   hslice->Fit(pi1, "N", "", -8, -1);
  55   hslice->Fit(el1, "NL", "", -.8, 3);
  56
  57   // Use single fits to constrain combined fit
  58   combined->FixParameter(1, pi1->GetParameter(1));
  59   combined->FixParameter(2, pi1->GetParameter(2));
  60   combined->FixParameter(4, pi1->GetParameter(4));
  61   combined->FixParameter(5, pi1->GetParameter(5));
  62   if(singlePiGauss)
  63     combined->FixParameter(3, 0);
  64   else
  65     combined->SetParLimits(3, minpi, maxpi);
  66   // Protons
  67   combined->FixParameter(7, pi1->GetParameter(7));
  68   combined->FixParameter(8, pi1->GetParameter(8));
  69   // Electrons
  70   combined->FixParameter(10, el1->GetParameter(1));
  71   combined->FixParameter(11, el1->GetParameter(2));
  72   combined->SetParLimits(0, minpi, maxpi);
  73   combined->SetParLimits(6, minpr, maxpr);
  74   combined->SetParLimits(9, minele, maxele);
  75
  76   hslice->Fit(combined, "N", "", -8., 3.);
  77
  78   for(Int_t ipar = 0; ipar < 9; ipar++) pi2->SetParameter(ipar, combined->GetParameter(ipar));
  79   for(Int_t ipar = 0; ipar < 3; ipar++) el2->SetParameter(ipar, combined->GetParameter(ipar+9));
  80
  81   combined->SetLineColor(kBlack);
  82   el2->SetLineColor(kGreen);
  83   pi2->SetLineColor(kBlue);
  84   combined->SetLineWidth(2);
  85   el2->SetLineWidth(2);
  86   pi2->SetLineWidth(2);
  87   el2->SetLineStyle(2);
  88   pi2->SetLineStyle(2);
  89
  90   hslice->SetStats(kFALSE);
  91   hslice->SetTitle("TPC Signal");
  92   hslice->GetXaxis()->SetTitle("TPC dE/dx - <dE/dx>|_{electrons} [#sigma]");
  93
  94   TCanvas *output = new TCanvas("output", "Slice output");
  95   output->SetLogy();
  96   output->cd();
  97   hslice->Draw();
  98   combined->Draw("same");
  99   el2->Draw("same");
 100   pi2->Draw("same");
 101
 102   TLegend *leg = new TLegend(0.5, 0.7, 0.89, 0.89);
 103   leg->SetFillStyle(0);
 104   leg->SetBorderSize(0);
 105   leg->AddEntry(hslice, "Measurement", "l");
 106   leg->AddEntry(combined, "Electrons + Pions", "l");
 107   leg->AddEntry(el2, "Electrons", "l");
 108   leg->AddEntry(pi2, "Pions", "l");
 109   leg->Draw();
 110
 111   TPaveText *momentum = new TPaveText(0.6, 0.45, 0.8, 0.55, "NDC");
 112   momentum->SetBorderSize(0);
 113   momentum->SetFillStyle(0);
 114   momentum->AddText(Form("p = %.2fGeV/c", p));
 115   momentum->Draw();
 116
 117   // Draw also Selection Bands
 118   TF1 *fBandPions = new TF1(*pi2), *fBandElectrons = new TF1(*el2);
 119   fBandPions->SetRange(cutmodel.Eval(p), 5);
 120   fBandElectrons->SetRange(cutmodel.Eval(p), 5);
 121
 122   fBandPions->SetFillColor(kBlue);
 123   fBandPions->SetFillStyle(3004);
 124   fBandElectrons->SetFillColor(kGreen);
 125   fBandElectrons->SetFillStyle(3005);
 126   fBandPions->SetLineWidth(0);
 127   fBandElectrons->SetLineWidth(0);
 128
 129   fBandElectrons->Draw("same");
 130   fBandPions->Draw("same");
 131
 132   // Now calculate the contamination
 133   Double_t nPions = pi2->Integral(cutmodel.Eval(p), 5);
 134   Double_t nCandidates = combined->Integral(cutmodel.Eval(p), 5);
 135   Double_t contamination = nPions/nCandidates;
 136   Double_t nCandidatesHisto = hslice->Integral(hslice->GetXaxis()->FindBin(cutmodel.Eval(p)), hslice->GetXaxis()->FindBin(5));
 137   Double_t nPionsBefore = pi1->Integral(cutmodel.Eval(p), 5);
 138   Double_t contaminationHisto = nPionsBefore/nCandidatesHisto;
 139   printf("Contamination @ %.3fGeV/c: %f\n", p, contamination);
 140   printf("Contamination @ %.3fGeV/c: %f (determined from double gauss and histogram)\n", p, contaminationHisto);
 141
 142   TPaveText *tcont = new TPaveText(0.45, 0.35, 0.9, 0.45, "NDC");
 143   tcont->SetBorderSize(0);
 144   tcont->SetFillStyle(0);
 145   tcont->AddText(Form("contamination: %f", contamination));
 146   tcont->Draw();
 147
 148   // Save Canvas
 149   TFile *out = new TFile(Form("slice%dpr.root", ((Int_t)(100*p))), "RECREATE");
 150   out->cd();
 151   output->Write();
 152   out->Close(); delete out;
 153 }