CalcPropLikeRatio.C

Go to the documentation of this file.

void FormatCanvas(TCanvas* canvas);

void FormatHistLong(TH1D* hist, int colour, int style);

void FormatHist(TH1D* hist, int colour, int style);

void CalcPropLikeRatio(const int nProbBins = 25)
{

  const int nParticles = 5;

  TString testFileName =
    "/data/t2k/phsxvg/DUNEPid/data/4APA/EventSamples/train_muon_realflux_split_2.root";

  TString trainFileName[nParticles] = {
    "/data/t2k/phsxvg/DUNEPid/data/4APA/TrainMVA/mvaPlots_muon_all.root",
    "/data/t2k/phsxvg/DUNEPid/data/4APA/TrainMVA/mvaPlots_electron_all.root",
    "/data/t2k/phsxvg/DUNEPid/data/4APA/TrainMVA/mvaPlots_proton_all.root",
    "/data/t2k/phsxvg/DUNEPid/data/4APA/TrainMVA/mvaPlots_pich_all.root",
    "/data/t2k/phsxvg/DUNEPid/data/4APA/TrainMVA/mvaPlots_photon_all.root"};

  TString weightFileName[nParticles] = {
    "/data/t2k/phsxvg/DUNEPid/data/4APA/weights/muon_all_BDT.weights.xml",
    "/data/t2k/phsxvg/DUNEPid/data/4APA/weights/electron_all_BDT.weights.xml",
    "/data/t2k/phsxvg/DUNEPid/data/4APA/weights/proton_all_BDT.weights.xml",
    "/data/t2k/phsxvg/DUNEPid/data/4APA/weights/pich_all_BDT.weights.xml",
    "/data/t2k/phsxvg/DUNEPid/data/4APA/weights/photon_all_BDT.weights.xml"};

  TFile* fileForBinning = TFile::Open(trainFileName[0]);
  fileForBinning->cd("Method_BDT/BDT");

  int numbins = MVA_BDT_Train_S->GetNbinsX();
  const int nbins = numbins;

  fileForBinning->Close();

  double likeRatioS[nParticles][nbins] = {0};

  double binWidth[nParticles], MVAMin[nParticles], MVAMax[nParticles];

  for (int p = 0; p < nParticles; p++) {
    cout << endl;
    cout << "Calculating likelihood ratios from file " << trainFileName[p] << endl;
    cout << endl;

    TFile* trainFile = TFile::Open(trainFileName[p]);
    trainFile->cd("Method_BDT/BDT");

    if (MVA_BDT_Train_S->GetNbinsX() != nbins) {
      cerr << "MVA_BDT_Train_S->GetNbinsX() not equal to nbins in file " << trainFileName[p]
           << endl;
      exit(1);
    }
    if (MVA_BDT_Train_B->GetNbinsX() != nbins) {
      cerr << "MVA_BDT_Train_B->GetNbinsX() not equal to nbins in file " << trainFileName[p]
           << endl;
      exit(1);
    }

    binWidth[p] = MVA_BDT_Train_S->GetBinWidth(1);
    MVAMin[p] = MVA_BDT_Train_S->GetBinLowEdge(1);
    MVAMax[p] = MVA_BDT_Train_S->GetBinLowEdge(nbins) + MVA_BDT_Train_S->GetBinWidth(nbins);

    if (MVA_BDT_Train_B->GetBinWidth(1) != binWidth[p]) {
      cerr << "MVA_BDT_Train_B->GetBinWidth(1) not equal to binWidth in file " << trainFileName[p]
           << endl;
      exit(1);
    }
    if (MVA_BDT_Train_B->GetBinLowEdge(1) != MVAMin[p]) {
      cerr << "MVA_BDT_Train_B->GetBinLowEdge(1) not equal to MVAMin in file " << trainFileName[p]
           << endl;
      exit(1);
    }
    if (MVA_BDT_Train_B->GetBinLowEdge(nbins) + MVA_BDT_Train_B->GetBinWidth(nbins) != MVAMax[p]) {
      cerr << "MVA_BDT_Train_B->GetBinLowEdge(nbins) + MVA_BDT_Train_B->GetBinWidth(nbins) not "
              "equal to MVAMax in file "
           << trainFileName[p] << endl;
      exit(1);
    }

    for (int i = 1; i <= nbins; i++) {
      likeRatioS[p][i] =
        MVA_BDT_Train_S->GetBinContent(i) / TMath::Max(0.01, MVA_BDT_Train_B->GetBinContent(i));
    }
  }

  TFile* testfile = TFile::Open(testFileName);
  testfile->cd();
  TTree* PIDVar = (TTree*)testfile->Get("mvaTree");

  float evalRatio, coreHaloRatio, concentration, conicalness, dEdxStart, dEdxEnd, dEdxEndRatio;

  PIDVar->SetBranchAddress("evalRatio", &evalRatio);
  PIDVar->SetBranchAddress("coreHaloRatio", &coreHaloRatio);
  PIDVar->SetBranchAddress("concentration", &concentration);
  PIDVar->SetBranchAddress("conicalness", &conicalness);
  PIDVar->SetBranchAddress("dEdxStart", &dEdxStart);
  PIDVar->SetBranchAddress("dEdxEnd", &dEdxEnd);
  PIDVar->SetBranchAddress("dEdxEndRatio", &dEdxEndRatio);

  double MVAValue;
  int MVABin;

  float evalRatioR, coreHaloRatioR, concentrationR, conicalnessR, dEdxStartR, dEdxEndR,
    dEdxEndRatioR;

  vector<double> likeRatio[nParticles];
  vector<double>::iterator likeIter;

  TH1D* hMVAValue[nParticles] = {0};
  for (int p = 0; p < nParticles; p++)
    hMVAValue[p] = new TH1D(Form("hMVAValue%d", p), "", nbins, -0.5, 0.5);

  for (int p = 0; p < nParticles; p++) {
    cout << endl;
    cout << "Reading weights from file " << weightFileName[p] << endl;

    TMVA::Reader* reader = new TMVA::Reader();

    reader->AddVariable("evalRatio", &evalRatioR);
    reader->AddVariable("coreHaloRatio", &coreHaloRatioR);
    reader->AddVariable("concentration", &concentrationR);
    reader->AddVariable("conicalness", &conicalnessR);
    reader->AddVariable("dEdxStart", &dEdxStartR);
    reader->AddVariable("dEdxEnd", &dEdxEndR);
    reader->AddVariable("dEdxEndRatio", &dEdxEndRatioR);

    reader->BookMVA("BDT", weightFileName[p]);

    for (int i = 0; i < PIDVar->GetEntries(); ++i) {
      PIDVar->GetEntry(i);

      evalRatioR = evalRatio;
      coreHaloRatioR = coreHaloRatio;
      concentrationR = concentration;
      conicalnessR = conicalness;
      dEdxStartR = dEdxStart;
      dEdxEndR = dEdxEnd;
      dEdxEndRatioR = dEdxEndRatio;

      MVAValue = reader->EvaluateMVA("BDT");

      MVABin = 1 + int((MVAValue - MVAMin[p]) / binWidth[p]);

      hMVAValue[p]->Fill(MVAValue);

      likeRatio[p].push_back(likeRatioS[p][MVABin]);
    }

    delete reader;
  }

  TCanvas* MVAVal = new TCanvas("MVAVal", "MVAVal", 800, 600);
  FormatCanvas(MVAVal);
  hMVAValue[0]->GetXaxis()->SetTitle("BDT variable");
  hMVAValue[0]->GetYaxis()->SetRangeUser(0, 250);
  FormatHistLong(hMVAValue[0], 2, 1);
  hMVAValue[0]->Draw();
  FormatHist(hMVAValue[1], 4, 1);
  hMVAValue[1]->Draw("same");
  FormatHist(hMVAValue[2], 6, 1);
  hMVAValue[2]->Draw("same");
  FormatHist(hMVAValue[3], 8, 1);
  hMVAValue[3]->Draw("same");
  FormatHist(hMVAValue[4], 1, 1);
  hMVAValue[4]->Draw("same");

  TLegend* lM = new TLegend(0.54, 0.57, 0.84, 0.87);
  lM->AddEntry(hMVAValue[0], "Muon weights", "L");
  lM->AddEntry(hMVAValue[1], "Electron weights", "L");
  lM->AddEntry(hMVAValue[2], "Proton weights", "L");
  lM->AddEntry(hMVAValue[3], "#pi^{#pm} weights", "L");
  lM->AddEntry(hMVAValue[4], "Photon weights", "L");
  lM->SetBorderSize(0);
  lM->Draw();

  int likeVectorSize = likeRatio[0].size();
  for (int p = 1; p < nParticles; p++) {
    if (likeRatio[p].size() != likeVectorSize) {
      cerr << "Likelihood ratio vector sizes not equal" << endl;
      exit(1);
    }
  }

  int vectorIndex = 0;

  double particleLikeRatio[nParticles] = {0};
  double totalLikeRatio = 0.0;

  TH1D* hParticlePLR[nParticles] = {0};

  for (int p = 0; p < nParticles; p++)
    hParticlePLR[p] = new TH1D(Form("hParticlePLR%d", p), "", nProbBins, 0.0, 1.0);

  TH1D* hLikeRatio[nParticles] = {0};

  for (int p = 0; p < nParticles; p++)
    hLikeRatio[p] = new TH1D(Form("hLikeRatio%d", p), "", 40, -10.0, 10.0);

  while (vectorIndex < likeRatio[0].size()) {
    totalLikeRatio = 0.0;

    for (int p = 0; p < nParticles; p++) {
      likeIter = likeRatio[p].begin() + vectorIndex;
      particleLikeRatio[p] = (*likeIter);
      totalLikeRatio += (*likeIter);
    }

    for (int p = 0; p < nParticles; p++) {
      //if particleLikeRatio[p] / totalLikeRatio = 1.0, this goes to overflow bin of hParticlePLR[p]
      //for this reason, set it to be no more than 0.9999
      hParticlePLR[p]->Fill(TMath::Min(0.9999, particleLikeRatio[p] / totalLikeRatio));
    }

    ++vectorIndex;
  }

  TCanvas* ParticlePLR = new TCanvas("ParticlePLR", "ParticlePLR", 800, 600);
  FormatCanvas(ParticlePLR);
  hParticlePLR[0]->GetXaxis()->SetTitle("Proportion of likelihood ratio");
  hParticlePLR[0]->GetYaxis()->SetRangeUser(0, 800);
  FormatHistLong(hParticlePLR[0], 2, 1);
  hParticlePLR[0]->Draw();
  FormatHist(hParticlePLR[1], 4, 1);
  hParticlePLR[1]->Draw("same");
  FormatHist(hParticlePLR[2], 6, 1);
  hParticlePLR[2]->Draw("same");
  FormatHist(hParticlePLR[3], 8, 1);
  hParticlePLR[3]->Draw("same");
  FormatHist(hParticlePLR[4], 1, 1);
  hParticlePLR[4]->Draw("same");

  TLegend* lP = new TLegend(0.6, 0.55, 0.83, 0.85);
  lP->AddEntry(hParticlePLR[0], "Muon", "L");
  lP->AddEntry(hParticlePLR[1], "Electron", "L");
  lP->AddEntry(hParticlePLR[2], "Proton", "L");
  lP->AddEntry(hParticlePLR[3], "#pi^{#pm}", "L");
  lP->AddEntry(hParticlePLR[4], "Photon", "L");
  lP->SetBorderSize(0);
  lP->Draw();

} //end of CalcPropLikeRatio()

void FormatCanvas(TCanvas* canvas)
{
  canvas->Divide(1);
  canvas->cd(1);
  gPad->SetLeftMargin(0.18);
  gPad->SetBottomMargin(0.15);
  gPad->SetRightMargin(0.15);
}

void FormatHistLong(TH1D* hist, int colour, int style)
{
  hist->SetTitle("");
  hist->SetLineColor(colour);
  hist->SetLineWidth(2);
  hist->SetLineStyle(style);
  hist->SetStats(0);
  hist->GetXaxis()->SetTitleOffset(1.2);
  hist->GetYaxis()->SetTitleOffset(1.4);
  hist->GetXaxis()->CenterTitle();
  hist->GetYaxis()->CenterTitle();
}

void FormatHist(TH1D* hist, int colour, int style)
{
  hist->SetTitle("");
  hist->SetLineColor(colour);
  hist->SetLineWidth(2);
  hist->SetLineStyle(style);
  hist->SetStats(0);
}