WrongSignBDTVars.cxx

Go to the documentation of this file.

#include "3FlavorAna/Vars/WrongSignBDTVars.h"

#include "CAFAna/Vars/CVNFinalStates.h"
#include "3FlavorAna/Vars/NueVarsExtra.h"

#include "CAFAna/Core/Utilities.h"

#include "StandardRecord/Proxy/SRProxy.h"

#include "TMVA/Reader.h"

namespace ana
{

  float BDTNueVars[5] = {0};
  float BDTNumuVars[5] = {0};
  float BDTCVNNue_p1Vars[4] = {0};
  float BDTCVNNue_p2Vars[6] = {0};
  float BDTCVNNumu_p1Vars[4] = {0};
  float BDTCVNNumu_p2Vars[6] = {0};
  static TMVA::Reader* fReaderBDTNue  = 0;
  static TMVA::Reader* fReaderBDTNumu  = 0;
  static TMVA::Reader* fReaderBDTCVNNue_p1  = 0;
  static TMVA::Reader* fReaderBDTCVNNue_p2  = 0;
  static TMVA::Reader* fReaderBDTCVNNumu_p1  = 0;
  static TMVA::Reader* fReaderBDTCVNNumu_p2  = 0;

  float GetAntiNueBDT::operator() (const caf::SRProxy* sr) const
  {
    float antinuebdt = -5.;

    BDTNueVars[0] = kStretch(sr);
    BDTNueVars[1] = kInElasticity(sr);
    BDTNueVars[2] = kPrimDedx0(sr);
    BDTNueVars[3] = kSecShwGap(sr);
    BDTNueVars[4] = kSecAvgDedx(sr);

    if(!fReaderBDTNue)  InitTMVA();
    antinuebdt = fReaderBDTNue->EvaluateMVA("BDTG");
    return antinuebdt;
  }

  void GetAntiNueBDT::InitTMVA() const
  {
    if(!fReaderBDTNue) fReaderBDTNue = new TMVA::Reader;
    std::string pidlib = FindCAFAnaDir()+"/data/wrong_sign/weights/TMVAClassification_BDTG2.weights.xml";
    fReaderBDTNue->AddVariable("stretch", &BDTNueVars[0]);
    fReaderBDTNue->AddVariable("inelasticity", &BDTNueVars[1]);
    fReaderBDTNue->AddVariable("prim_dedx0", &BDTNueVars[2]);
    fReaderBDTNue->AddVariable("sec_shw_gap", &BDTNueVars[3]);
    fReaderBDTNue->AddVariable("sec_avg_dedx", &BDTNueVars[4]);
    fReaderBDTNue->BookMVA("BDTG", pidlib);
  }

  float GetAntiNumuBDT::operator() (const caf::SRProxy* sr) const
  {
    float antinumubdt = -5.;

    BDTNumuVars[0] = kStretch(sr);
    BDTNumuVars[1] = kNumuHadVisEFrac(sr);
    BDTNumuVars[2] = kPrimDedx0(sr);
    BDTNumuVars[3] = kSecShwGap(sr);
    BDTNumuVars[4] = kSecAvgDedx(sr);

    if(!fReaderBDTNumu)  InitTMVA();
    antinumubdt = fReaderBDTNumu->EvaluateMVA("BDTG");
    return antinumubdt;
  }

  void GetAntiNumuBDT::InitTMVA() const
  {
    if(!fReaderBDTNumu) fReaderBDTNumu = new TMVA::Reader;
    std::string pidlib = FindCAFAnaDir()+"/data/wrong_sign/weights/TMVAClassification_BDTG2.weights.xml";
    fReaderBDTNumu->AddVariable("stretch", &BDTNumuVars[0]);
    fReaderBDTNumu->AddVariable("inelasticity", &BDTNumuVars[1]);
    fReaderBDTNumu->AddVariable("prim_dedx0", &BDTNumuVars[2]);
    fReaderBDTNumu->AddVariable("sec_shw_gap", &BDTNumuVars[3]);
    fReaderBDTNumu->AddVariable("sec_avg_dedx", &BDTNumuVars[4]);
    fReaderBDTNumu->BookMVA("BDTG", pidlib);
  }


  float GetAntiNueBDTCVN::operator() (const caf::SRProxy* sr) const
  {
    float antinuebdt = -5.;

    BDTCVNNue_p1Vars[0] = kStretch(sr);
    BDTCVNNue_p1Vars[1] = kInElasticity(sr);
    BDTCVNNue_p1Vars[2] = kPrimDedx0(sr);
    BDTCVNNue_p1Vars[3] = kCVNFSProtonScore2018(sr);

    BDTCVNNue_p2Vars[0] = kStretch(sr);
    BDTCVNNue_p2Vars[1] = kInElasticity(sr);
    BDTCVNNue_p2Vars[2] = kPrimDedx0(sr);
    BDTCVNNue_p2Vars[3] = kCVNFSProtonScore2018(sr);
    BDTCVNNue_p2Vars[4] = kProngCVNMaxProtonDist(sr);
    BDTCVNNue_p2Vars[5] = kProngCVNMaxProton(sr);

    if(!fReaderBDTCVNNue_p1 || !fReaderBDTCVNNue_p2)  InitTMVA();
    if(kNProng(sr) == 1)
      antinuebdt = fReaderBDTCVNNue_p1->EvaluateMVA("BDTG");
    if(kNProng(sr) > 1)
      antinuebdt = fReaderBDTCVNNue_p2->EvaluateMVA("BDTG");

    return antinuebdt;
  }

  void GetAntiNueBDTCVN::InitTMVA() const
  {
    if(!fReaderBDTCVNNue_p1) fReaderBDTCVNNue_p1 = new TMVA::Reader;
    if(!fReaderBDTCVNNue_p2) fReaderBDTCVNNue_p2 = new TMVA::Reader;
    std::string pidlib_p1 = FindCAFAnaDir()+"/data/wrong_sign/weights/TMVAClassification_BDTGCVN2_p1.weights.xml";
    std::string pidlib_p2 = FindCAFAnaDir()+"/data/wrong_sign/weights/TMVAClassification_BDTGCVN2_p2.weights.xml";
    fReaderBDTCVNNue_p1->AddVariable("stretch", &BDTCVNNue_p1Vars[0]);
    fReaderBDTCVNNue_p1->AddVariable("inelasticity", &BDTCVNNue_p1Vars[1]);
    fReaderBDTCVNNue_p1->AddVariable("prim_dedx0", &BDTCVNNue_p1Vars[2]);
    fReaderBDTCVNNue_p1->AddVariable("cvnfsp", &BDTCVNNue_p1Vars[3]);
    fReaderBDTCVNNue_p1->BookMVA("BDTG", pidlib_p1);

    fReaderBDTCVNNue_p2->AddVariable("stretch", &BDTCVNNue_p2Vars[0]);
    fReaderBDTCVNNue_p2->AddVariable("inelasticity", &BDTCVNNue_p2Vars[1]);
    fReaderBDTCVNNue_p2->AddVariable("prim_dedx0", &BDTCVNNue_p2Vars[2]);
    fReaderBDTCVNNue_p2->AddVariable("cvnfsp", &BDTCVNNue_p2Vars[3]);
    fReaderBDTCVNNue_p2->AddVariable("maxprotondist", &BDTCVNNue_p2Vars[4]);
    fReaderBDTCVNNue_p2->AddVariable("maxproton", &BDTCVNNue_p2Vars[5]);
    fReaderBDTCVNNue_p2->BookMVA("BDTG", pidlib_p2);
  }

  float GetAntiNumuBDTCVN::operator() (const caf::SRProxy* sr) const
  {
    float antinumubdt = -5.;

    BDTCVNNumu_p1Vars[0] = kStretch(sr);
    BDTCVNNumu_p1Vars[1] = kNumuHadVisEFrac(sr);
    BDTCVNNumu_p1Vars[2] = kPrimDedx0(sr);
    BDTCVNNumu_p1Vars[3] = kCVNFSProtonScore2018(sr);

    BDTCVNNumu_p2Vars[0] = kStretch(sr);
    BDTCVNNumu_p2Vars[1] = kNumuHadVisEFrac(sr);
    BDTCVNNumu_p2Vars[2] = kPrimDedx0(sr);
    BDTCVNNumu_p2Vars[3] = kCVNFSProtonScore2018(sr);
    BDTCVNNumu_p2Vars[4] = kProngCVNMaxProtonDist(sr);
    BDTCVNNumu_p2Vars[5] = kProngCVNMaxProton(sr);

    if(!fReaderBDTCVNNumu_p1 || !fReaderBDTCVNNumu_p2)  InitTMVA();
    if(kNProng(sr) == 1)
      antinumubdt = fReaderBDTCVNNumu_p1->EvaluateMVA("BDTG");
    if(kNProng(sr) > 1)
      antinumubdt = fReaderBDTCVNNumu_p2->EvaluateMVA("BDTG");

    return antinumubdt;
  }

  void GetAntiNumuBDTCVN::InitTMVA() const
  {
    if(!fReaderBDTCVNNumu_p1) fReaderBDTCVNNumu_p1 = new TMVA::Reader;
    if(!fReaderBDTCVNNumu_p2) fReaderBDTCVNNumu_p2 = new TMVA::Reader;
    std::string pidlib_p1 = FindCAFAnaDir()+"/data/wrong_sign/weights/TMVAClassification_BDTGCVN2_p1.weights.xml";
    std::string pidlib_p2 = FindCAFAnaDir()+"/data/wrong_sign/weights/TMVAClassification_BDTGCVN2_p2.weights.xml";
    fReaderBDTCVNNumu_p1->AddVariable("stretch", &BDTCVNNumu_p1Vars[0]);
    fReaderBDTCVNNumu_p1->AddVariable("inelasticity", &BDTCVNNumu_p1Vars[1]);
    fReaderBDTCVNNumu_p1->AddVariable("prim_dedx0", &BDTCVNNumu_p1Vars[2]);
    fReaderBDTCVNNumu_p1->AddVariable("cvnfsp", &BDTCVNNumu_p1Vars[3]);
    fReaderBDTCVNNumu_p1->BookMVA("BDTG", pidlib_p1);

    fReaderBDTCVNNumu_p2->AddVariable("stretch", &BDTCVNNumu_p2Vars[0]);
    fReaderBDTCVNNumu_p2->AddVariable("inelasticity", &BDTCVNNumu_p2Vars[1]);
    fReaderBDTCVNNumu_p2->AddVariable("prim_dedx0", &BDTCVNNumu_p2Vars[2]);
    fReaderBDTCVNNumu_p2->AddVariable("cvnfsp", &BDTCVNNumu_p2Vars[3]);
    fReaderBDTCVNNumu_p2->AddVariable("maxprotondist", &BDTCVNNumu_p2Vars[4]);
    fReaderBDTCVNNumu_p2->AddVariable("maxproton", &BDTCVNNumu_p2Vars[5]);
    fReaderBDTCVNNumu_p2->BookMVA("BDTG", pidlib_p2);
  }

  const Var kStretch([](const caf::SRProxy* sr)
                    {
                      if(sr->slc.meanpos.z == sr->slc.boxmin.z)
                        return -5.;
                      return (double)((sr->slc.boxmax.z-sr->slc.meanpos.z)/(sr->slc.meanpos.z-sr->slc.boxmin.z));
                    });

  const Var kEthetaSqr([](const caf::SRProxy* sr)
                       {
                        if(!sr->vtx.elastic.IsValid)
                          return -5.;
                        if(sr->vtx.elastic.fuzzyk.png.size() == 0)
                          return -5.;
                        return (double)(sr->vtx.elastic.fuzzyk.png[0].calE*std::pow(acos(sr->vtx.elastic.fuzzyk.png[0].dir.z),2));
                        });

  const Var kInElasticity([](const caf::SRProxy* sr)
                        {
                          if(!sr->vtx.elastic.IsValid) return 1.f;
                          if(sr->vtx.elastic.fuzzyk.nshwlid < 1) return 1.f;
                          return ((sr->slc.calE -
                                   sr->vtx.elastic.fuzzyk.png[0].shwlid.calE)/(sr->slc.calE));
                        });

   const Var kNumuHadVisEFrac(
              [](const caf::SRProxy* sr)
              {
                if(sr->slc.calE == 0) return -5.;
                return (double)(sr->energy.numu.hadcalE+sr->energy.numu.hadtrkE)/(sr->slc.calE);
              });

   const Var kNumuHadEFrac(
              [](const caf::SRProxy* sr)
              {
                if(sr->slc.calE == 0) return -5.;
                return (double)(sr->energy.numu.hadcalE)/(sr->slc.calE);
              });

  const Var kLeptonAngle([](const caf::SRProxy* sr)
                        {
                          if(!sr->vtx.elastic.IsValid)
                            return -5.;
                          if(sr->vtx.elastic.fuzzyk.png.size() == 0)
                            return -5.;
                          return (double)acos(sr->vtx.elastic.fuzzyk.png[0].dir.z)*180./M_PI;
                        });

  #define PrimVAR(VAR)                                  \
      Var(                                                      \
          [](const caf::SRProxy* sr)                    \
          {                                                     \
              if(!sr->vtx.elastic.IsValid) return -5.;        \
              if(sr->vtx.elastic.fuzzyk.nshwlid < 1) return -5.;        \
              return double(sr->vtx.elastic.fuzzyk.png[0].shwlid.VAR);  \
          });

  #define SecVAR(VAR)                                   \
      Var(                                                      \
          [](const caf::SRProxy* sr)                    \
          {                                                     \
              if(!sr->vtx.elastic.IsValid) return -5.;        \
              if(sr->vtx.elastic.fuzzyk.nshwlid< 2) return 0.;        \
              return double(sr->vtx.elastic.fuzzyk.png[1].shwlid.VAR);  \
          });

  const Var kPrimDedx0 = PrimVAR(lid.dedx0)
  const Var kPrimDedx1 = PrimVAR(lid.dedx1)

  const Var kPrimAvgDedx([](const caf::SRProxy* sr)
                     {
                      if(!sr->vtx.elastic.IsValid)
                        return -5.;
                      if(sr->vtx.elastic.fuzzyk.nshwlid < 1) return -5.;
                      if(sr->vtx.elastic.fuzzyk.npng < 1) return -5.;
                      if(sr->vtx.elastic.fuzzyk.png[0].len < 1e-10) return -5.;
                      return (double)(sr->vtx.elastic.fuzzyk.png[0].calE/sr->vtx.elastic.fuzzyk.png[0].len);
                      });
  const Var kSecAvgDedx([](const caf::SRProxy* sr)
                     {
                      if(!sr->vtx.elastic.IsValid)
                        return -5.;
                      if(sr->vtx.elastic.fuzzyk.nshwlid < 2) return 0.;
                      if(sr->vtx.elastic.fuzzyk.npng < 2) return 0.;
                      if(sr->vtx.elastic.fuzzyk.png[1].len < 1e-10) return -5.;
                      return (double)(sr->vtx.elastic.fuzzyk.png[1].calE/sr->vtx.elastic.fuzzyk.png[1].len);
                      });

  const Var kSecShwGap(
          [](const caf::SRProxy* sr)
          {
              if(!sr->vtx.elastic.IsValid) return -5.;
              if(sr->vtx.elastic.fuzzyk.nshwlid< 1) return -5.;
              if(sr->vtx.elastic.fuzzyk.nshwlid < 2) return (double)sr->vtx.elastic.fuzzyk.png[0].shwlid.gap;
              return double(sr->vtx.elastic.fuzzyk.png[1].shwlid.gap);
          });
  const Var kSecShwLen = SecVAR(len)

  const Var kNProng(
      [](const caf::SRProxy* sr)
      {
        if(!sr->vtx.elastic.IsValid) return -5;
        return (int)sr->vtx.elastic.fuzzyk.npng;
      });

}