C++ Reader::EvaluateMVA方法代码示例

Float_t hhMVA::GetBDTValue(Float_t mTT, Float_t ptTT, Float_t mBB, Float_t ptBB,
			   Float_t mHH, Float_t ptHH, Float_t mt2,
			   Float_t dRbb, Float_t dRtt, Float_t dRhh) {
  
  if (!fTMVAReader) {
    cout << "TMVA reader not initialized properly" << endl;
    return -999;
  }

  fMVAVar_mTT   =mTT;
  fMVAVar_ptTT  =ptTT;
  fMVAVar_mBB1  =mBB;
  fMVAVar_ptBB1 =ptBB;
  fMVAVar_mHH   =mHH;
  fMVAVar_ptHH  =ptHH;
  fMVAVar_mt2   =mt2;
  fMVAVar_dRbb  =dRbb;
  fMVAVar_dRtt  =dRtt;
  fMVAVar_dRhh  =dRhh;

  TMVA::Reader *reader = 0;
  reader = fTMVAReader;
  
  return reader->EvaluateMVA("BDT method");
  
}

int reader_wrapper::GetEntry(Long64_t e) {
  /// don't care about spectators here
  for (auto b: m_branches) {
    b->GetEntry(e);
  }
  for (auto& v : m_variables) {
    v.value = v.ttreeformula->EvalInstance();
  }
  m_response = m_reader->EvaluateMVA(m_methodName.Data());
  m_responseBranch->Fill();
  return 0;
}

void TMVAPredict()
{
  std::ofstream outfile ("baseline_c.csv");
  outfile << "id,prediction\n";

  TMVA::Tools::Instance();

  std::cout << "==> Start TMVAPredict" << std::endl;
  TMVA::Reader *reader = new TMVA::Reader( "!Color:!Silent" );  
  string variables_name[3] = {"LifeTime",
                           "FlightDistance",
                           "pt"}
  Float_t variables[3];
  for (int i=0; i < 3; i++){
    reader->AddVariable(variables_name[i].c_str(), &variables[i]);
    variables[i] = 0.0;
  }

  TString dir    = "weights/";
  TString prefix = "TMVAClassification";
  TString method_name = "GBDT";
  TString weightfile = dir + prefix + TString("_") + method_name + TString(".weights.xml");
  reader->BookMVA( method_name, weightfile ); 

  TFile *input(0);
  input = TFile::Open("../tau_data/test.root");
  TTree* tree = (TTree*)input->Get("data");
  
  Int_t ids;
  Float_t prediction;
  tree->SetBranchAddress("id", &ids);

  for (int i=0; i < 3; i++){
    tree->SetBranchAddress(variables_name[i].c_str(), &variables[i]);
  }
 
  for (Long64_t ievt=0; ievt < tree->GetEntries(); ievt++) {
    tree->GetEntry(ievt);
    prediction = reader->EvaluateMVA(method_name);
    outfile << ids << "," << (prediction + 1.) / 2. << "\n";
  }

  outfile.close();
  input->Close();
  delete reader;
}

//.........这里部分代码省略.........
    intree->SetBranchAddress("SubJet2_CSV"         ,&SubJet2_CSV         );
    intree->SetBranchAddress("SubJet2_CSVIVF"      ,&SubJet2_CSVIVF      );
    //--------------------------------------
    // CSV TaggingVariables
    //--------------------------------------
    intree->SetBranchAddress("TagVarCSV2_jetNTracks"               ,&TagVarCSV2_jetNTracks              );
    intree->SetBranchAddress("TagVarCSV2_jetNTracksEtaRel"         ,&TagVarCSV2_jetNTracksEtaRel        );
    intree->SetBranchAddress("TagVarCSV2_trackSumJetEtRatio"       ,&TagVarCSV2_trackSumJetEtRatio      );
    intree->SetBranchAddress("TagVarCSV2_trackSumJetDeltaR"        ,&TagVarCSV2_trackSumJetDeltaR       );
    intree->SetBranchAddress("TagVarCSV2_trackSip2dValAboveCharm"  ,&TagVarCSV2_trackSip2dValAboveCharm );
    intree->SetBranchAddress("TagVarCSV2_trackSip2dSigAboveCharm"  ,&TagVarCSV2_trackSip2dSigAboveCharm );
    intree->SetBranchAddress("TagVarCSV2_trackSip3dValAboveCharm"  ,&TagVarCSV2_trackSip3dValAboveCharm );
    intree->SetBranchAddress("TagVarCSV2_trackSip3dSigAboveCharm"  ,&TagVarCSV2_trackSip3dSigAboveCharm );
    intree->SetBranchAddress("TagVarCSV2_vertexCategory"           ,&TagVarCSV2_vertexCategory          );
    intree->SetBranchAddress("TagVarCSV2_jetNSecondaryVertices"    ,&TagVarCSV2_jetNSecondaryVertices   );
    intree->SetBranchAddress("TagVarCSV2_vertexMass"               ,&TagVarCSV2_vertexMass              );
    intree->SetBranchAddress("TagVarCSV2_vertexNTracks"            ,&TagVarCSV2_vertexNTracks           );
    intree->SetBranchAddress("TagVarCSV2_vertexEnergyRatio"        ,&TagVarCSV2_vertexEnergyRatio       );
    intree->SetBranchAddress("TagVarCSV2_vertexJetDeltaR"          ,&TagVarCSV2_vertexJetDeltaR         );
    intree->SetBranchAddress("TagVarCSV2_flightDistance2dVal"      ,&TagVarCSV2_flightDistance2dVal     );
    intree->SetBranchAddress("TagVarCSV2_flightDistance2dSig"      ,&TagVarCSV2_flightDistance2dSig     );
    intree->SetBranchAddress("TagVarCSV2_flightDistance3dVal"      ,&TagVarCSV2_flightDistance3dVal     );
    intree->SetBranchAddress("TagVarCSV2_flightDistance3dSig"      ,&TagVarCSV2_flightDistance3dSig     );
    intree->SetBranchAddress("TagVarCSV2_trackEtaRel_0"            ,&TagVarCSV2_trackEtaRel_0           );
    intree->SetBranchAddress("TagVarCSV2_trackEtaRel_1"            ,&TagVarCSV2_trackEtaRel_1           );
    intree->SetBranchAddress("TagVarCSV2_trackEtaRel_2"            ,&TagVarCSV2_trackEtaRel_2           );
    
    std::cout << "Now looping over " << intree->GetEntries() << " entries..." << std::endl;
    for(Long64_t iEntry = 0; iEntry < intree->GetEntries(); iEntry++){
	    if (iEntry % 1000 == 0) std::cout << "Processing Entry #" << iEntry << std::endl;
	    intree->GetEntry(iEntry); // all variables now filled!

	    bool isBkg = ( Jet_massGroomed>80 && Jet_massGroomed<150 );
	    float BDTG_Disc = reader->EvaluateMVA("BDTG_T1000D3_fat_BBvsQCD method");

	    if (isBkg) {
		    hBDTGDiscBkg->Fill(BDTG_Disc);
		    hFatCSVIVFDiscBkg->Fill(Jet_CSVIVF);
		    hSubCSVIVFDiscBkg->Fill(std::min(SubJet1_CSVIVF,SubJet2_CSVIVF));
	    }
    }

    // signal input tree
    TString infilenameSig="RadionToHH_4b_M-800_TuneZ2star_8TeV-Madgraph_pythia6_JetTaggingVariables_evaluation.root";
    TFile inSig(infilenameSig);
    intree = (TTree*)inSig.Get("tagVars/ttree");

    // set the branches to point to address of the variables declared above
    //######################################
    // Fat jet variables
    //######################################
    intree->SetBranchAddress("Jet_pt"          ,&Jet_pt          );
    intree->SetBranchAddress("Jet_eta"         ,&Jet_eta         );
    intree->SetBranchAddress("Jet_phi"         ,&Jet_phi         );
    intree->SetBranchAddress("Jet_mass"        ,&Jet_mass        );
    intree->SetBranchAddress("Jet_massGroomed" ,&Jet_massGroomed );
    intree->SetBranchAddress("Jet_flavour"     ,&Jet_flavour     );
    intree->SetBranchAddress("Jet_nbHadrons"   ,&Jet_nbHadrons   );
    intree->SetBranchAddress("Jet_JP"          ,&Jet_JP          );
    intree->SetBranchAddress("Jet_JBP"         ,&Jet_JBP         );
    intree->SetBranchAddress("Jet_CSV"         ,&Jet_CSV         );
    intree->SetBranchAddress("Jet_CSVIVF"      ,&Jet_CSVIVF      );
    intree->SetBranchAddress("Jet_tau1"        ,&Jet_tau1        );
    intree->SetBranchAddress("Jet_tau2"        ,&Jet_tau2        );
    //--------------------------------------
    // CSV TaggingVariables

//.........这里部分代码省略.........

TFile *tmp  = new TFile( "tmp.root","RECREATE" );

TTree* theTree = BigTree->CopyTree("((cat == 1) + (cat == 2))*(ln==0)*(Cosmic==0)*(fabs(Mass_Z - 91.18)<10)*(Pt_Z>30)*(DeltaPhi_metjet>0.5)*(Pt_J1 < 30)*(pfMEToverPt_Z > 0.4)*(pfMEToverPt_Z < 1.8)*((Pt_Jet_btag_CSV_max > 20)*(btag_CSV_max < 0.244) + (1-(Pt_Jet_btag_CSV_max > 20)))*(sqrt(pow(dilepPROJLong + 1.25*recoilPROJLong + 0.0*uncertPROJLong,2)*(dilepPROJLong + 1.25*recoilPROJLong + 0.0*uncertPROJLong > 0) + 1.0*pow(dilepPROJPerp + 1.25*recoilPROJPerp + 0.0*uncertPROJPerp,2)*(dilepPROJPerp + 1.25*recoilPROJPerp + 0.0*uncertPROJPerp > 0)) > 45.0)");
   std::cout << "--- Select signal sample" << std::endl;
   Float_t userVar1, userVar2;
//   theTree->SetBranchAddress( "var1", &userVar1 );
//   theTree->SetBranchAddress( "var2", &userVar2 );
//   theTree->SetBranchAddress( "var3", &var3 );
//   theTree->SetBranchAddress( "var4", &var4 );

theTree->SetBranchAddress( " Z_rapidity_z", &Z_rapidity_z);

theTree->SetBranchAddress( " THRUST_2D", &THRUST_2D);

theTree->SetBranchAddress( " L1_L2_cosangle", &L1_L2_cosangle);

theTree->SetBranchAddress( " TransMass_ZH150_uncl", &TransMass_ZH150_uncl);

theTree->SetBranchAddress( " TransMass_ZH150", &TransMass_ZH150);

theTree->SetBranchAddress( " DeltaPhi_ZH", &DeltaPhi_ZH);

theTree->SetBranchAddress( " DeltaPhi_ZH_uncl", &DeltaPhi_ZH_uncl);

theTree->SetBranchAddress( " CMAngle", &CMAngle);

theTree->SetBranchAddress( " CS_cosangle", &CS_cosangle);

   // efficiency calculator for cut method
   Int_t    nSelCutsGA = 0;
   Double_t effS       = 0.7;

   std::vector<Float_t> vecVar(9); // vector for EvaluateMVA tests

   std::cout << "--- Processing: " << theTree->GetEntries() << " events" << std::endl;
   TStopwatch sw;
   sw.Start();
   for (Long64_t ievt=0; ievt<theTree->GetEntries();ievt++) {

      if (ievt%1000 == 0){
         std::cout << "--- ... Processing event: " << ievt << std::endl;
      }

      theTree->GetEntry(ievt);

      var1 = userVar1 + userVar2;
      var2 = userVar1 - userVar2;

      if (ievt <20){
         // test the twodifferent Reader::EvaluateMVA functions 
         // access via registered variables compared to access via vector<float>
//         vecVar[0]=var1;
//         vecVar[1]=var2;
//         vecVar[2]=var3;
//         vecVar[3]=var4;      

vecVar[0]=Z_rapidity_z;

vecVar[1]=THRUST_2D;

vecVar[2]=L1_L2_cosangle;

vecVar[3]=TransMass_ZH150_uncl;

vecVar[4]=TransMass_ZH150;

//.........这里部分代码省略.........
    Float_t lepsoft_fbrem_;
    Float_t lepsoft_eOverPIn_;
    Float_t lepsoft_q_;
    Float_t lepsoft_dPhiIn_;

    theTree->SetBranchAddress( "lephard_pt_"             ,   &lephard_pt_              ); 
    theTree->SetBranchAddress( "lepsoft_pt_"             ,   &lepsoft_pt_              ); 
    theTree->SetBranchAddress( "lepsoft_fr_"             ,   &lepsoft_fr_              ); 
    theTree->SetBranchAddress( "dil_dphi_"               ,   &dil_dphi_                ); 
    theTree->SetBranchAddress( "dil_mass_"               ,   &dil_mass_                ); 
    theTree->SetBranchAddress( "event_type_"             ,   &event_type_              ); 
    theTree->SetBranchAddress( "met_projpt_"             ,   &met_projpt_              ); 
    theTree->SetBranchAddress( "jets_num_"               ,   &jets_num_                ); 
    theTree->SetBranchAddress( "extralep_num_"           ,   &extralep_num_            ); 
    theTree->SetBranchAddress( "lowptbtags_num_"         ,   &lowptbtags_num_          ); 
    theTree->SetBranchAddress( "softmu_num_"             ,   &softmu_num_              ); 
    theTree->SetBranchAddress( "event_scale1fb_"         ,   &event_scale1fb_          ); 
    theTree->SetBranchAddress( "lepsoft_passTighterId_"  ,   &lepsoft_passTighterId_   );
    theTree->SetBranchAddress( "met_pt_"                 ,   &met_pt_                  );
    theTree->SetBranchAddress( "mt_lephardmet_"          ,   &mt_lephardmet_           );
    theTree->SetBranchAddress( "mt_lepsoftmet_"          ,   &mt_lepsoftmet_           );
    theTree->SetBranchAddress( "mthiggs_"                ,   &mthiggs_                 );
    theTree->SetBranchAddress( "dphi_lephardmet_"        ,   &dphi_lephardmet_         );
    theTree->SetBranchAddress( "dphi_lepsoftmet_"        ,   &dphi_lepsoftmet_         );
    theTree->SetBranchAddress( "lepsoft_fbrem_"          ,   &lepsoft_fbrem_           );
    theTree->SetBranchAddress( "lepsoft_eOverPIn_"       ,   &lepsoft_eOverPIn_        );
    theTree->SetBranchAddress( "lepsoft_q_"              ,   &lepsoft_q_               );
    theTree->SetBranchAddress( "lepsoft_dPhiIn_"         ,   &lepsoft_dPhiIn_          );

    // Efficiency calculator for cut method
    Int_t    nSelCutsGA = 0;
    Double_t effS       = 0.7;

    std::vector<Float_t> vecVar(4); // vector for EvaluateMVA tests

    std::cout << "--- Processing: " << theTree->GetEntries() << " events" << std::endl;
    TStopwatch sw;
    sw.Start();

    int npass   = 0;
    float yield = 0.;
    
    for (Long64_t ievt=0; ievt<theTree->GetEntries();ievt++) {

      if (ievt%1000 == 0) std::cout << "--- ... Processing event: " << ievt << std::endl;

      theTree->GetEntry(ievt);

      //-------------------------------------------------------
      // event selection
      //-------------------------------------------------------

      if( dil_dphi_ > 1. ) continue;

      //em
      if( event_type_ > 0.5 && event_type_ < 2.5 ){
        if( met_projpt_ < 20. )   continue;
      }
      //ee/mm
      if( event_type_ < 0.5 || event_type_ > 2.5 ){
        if( met_projpt_ < 35. )   continue;
      }
      if( lephard_pt_ < 20.           )             continue;
      if( jets_num_ > 0               )             continue;
      if( extralep_num_ > 0           )             continue;
      if( lowptbtags_num_ > 0         )             continue;

void apply(std::string iName="train/OutputTmp.root") { 
  TMVA::Tools::Instance();
  TMVA::Reader *reader = new TMVA::Reader( "!Color:!Silent" );    

  float lPt        = 0; reader->AddVariable("pt"                 , &lPt);
  //float lEta       = 0; reader->AddVariable("eta"                , &lEta);
  //float lDR        = 0; reader->AddVariable("dR"                 , &lDR);
  //float lPtc       = 0; reader->AddVariable("ptc"                , &lPtc);
  // float lPtdR      = 0; reader->AddVariable("ptdR"               , &lPtdR);
  //float lPuppi     = 0; reader->AddVariable("puppi"              , &lPuppi);
  float lPtODR     = 0; reader->AddVariable("ptodR"              , &lPtODR);
  //float lPtODRS    = 0; reader->AddVariable("ptodRS"             , &lPtODRS);
  float lPtODRSO   = 0; reader->AddVariable("ptodRSO"            , &lPtODRSO);
  //float lDRLV      = 0; reader->AddVariable("dR_lv"              , &lDRLV);
  //float lPtcLV     = 0; reader->AddVariable("ptc_lv"             , &lPtcLV);
  //float lPtdRLV    = 0; reader->AddVariable("ptdR_lv"            , &lPtdRLV);
  //float lPuppiLV   = 0; reader->AddVariable("puppi_lv"           , &lPuppiLV);
  float lPtODRLV   = 0; reader->AddVariable("ptodR_lv"           , &lPtODRLV);
  //float lPtODRSLV  = 0; reader->AddVariable("ptodRS_lv"          , &lPtODRSLV);
  float lPtODRSOLV = 0; reader->AddVariable("ptodRSO_lv"         , &lPtODRSOLV);
  //float lDRPU      = 0; reader->AddVariable("dR_pu"              , &lDRPU);
  //float lPtcPU     = 0; reader->AddVariable("pt_pu"              , &lPtcPU);
  //float lPtdRPU    = 0; reader->AddVariable("ptdR_pu"            , &lPtdRPU);
  //float lPuppiPU   = 0; reader->AddVariable("puppi_pu"           , &lPuppiPU);
  //float lPtODRPU   = 0; reader->AddVariable("ptodR_pu"           , &lPtODRPU);
  //float lPtODRSPU  = 0; reader->AddVariable("ptodRS_pu"          , &lPtODRSPU);
  //float lPtODRSOPU = 0; reader->AddVariable("ptodRSO_pu"         , &lPtODRSOPU);
  
  std::string lJetName = "BDT";
  reader->BookMVA(lJetName .c_str(),(std::string("weights/TMVAClassificationCategory_PUDisc_v1")+std::string(".weights.xml")).c_str());
  
  TFile *lFile = new TFile(iName.c_str());
  TTree *lTree = (TTree*) lFile->Get("tree");
   lTree->SetBranchAddress("pt"                 , &lPt);
  //lTree->SetBranchAddress("eta"                , &lEta);
   //lTree->SetBranchAddress("dR"                 , &lDR);
   //lTree->SetBranchAddress("ptc"                , &lPtc);
   //lTree->SetBranchAddress("ptdR"               , &lPtdR);
  //lTree->SetBranchAddress("puppi"              , &lPuppi);
  lTree->SetBranchAddress("ptodR"              , &lPtODR);
  //lTree->SetBranchAddress("ptodRS"             , &lPtODRS);
  lTree->SetBranchAddress("ptodRSO"            , &lPtODRSO);
  //lTree->SetBranchAddress("dR_lv"              , &lDRLV);
  // lTree->SetBranchAddress("ptc_lv"             , &lPtcLV);
  //lTree->SetBranchAddress("ptdR_lv"            , &lPtdRLV);
  //lTree->SetBranchAddress("puppi_lv"           , &lPuppiLV);
  lTree->SetBranchAddress("ptodR_lv"           , &lPtODRLV);
  //lTree->SetBranchAddress("ptodRS_lv"          , &lPtODRSLV);
  lTree->SetBranchAddress("ptodRSO_lv"         , &lPtODRSOLV);
  //lTree->SetBranchAddress("dR_pu"              , &lDRPU);
  //lTree->SetBranchAddress("pt_pu"              , &lPtcPU);
  //lTree->SetBranchAddress("ptdR_pu"            , &lPtdRPU);
  //lTree->SetBranchAddress("puppi_pu"           , &lPuppiPU);
  //lTree->SetBranchAddress("ptodR_pu"           , &lPtODRPU);
  //lTree->SetBranchAddress("ptodRS_pu"          , &lPtODRSPU);
  //lTree->SetBranchAddress("ptodRSO_pu"         , &lPtODRSOPU);
    
  int lNEvents = lTree->GetEntries();
  TFile *lOFile = new TFile("Output.root","RECREATE");
  TTree *lOTree = lTree->CloneTree(0);
  float lMVA    = 0; lOTree->Branch("bdt"     ,&lMVA ,"lMVA/F");
  for (Long64_t i0=0; i0<lNEvents;i0++) {
    if (i0 % 10000 == 0) std::cout << "--- ... Processing event: " << double(i0)/double(lNEvents) << std::endl;
    lTree->GetEntry(i0);
    lMVA      = float(reader->EvaluateMVA(lJetName.c_str()));
    lOTree->Fill();
  }
  lOTree->Write();
  lOFile->Close();
  delete reader;
}

//.........这里部分代码省略.........
  TString pExpress0 = "2*fjet1QGtagSub2+fjet1QGtagSub1";
  TString pExpr0(pExpress0);
  TTreeFormula* lFVars0 = new TTreeFormula(pExpr0,pExpr0,lTree);
  TString pExpress1 = "fjet1QGtagSub1";
  TString pExpr1(pExpress1);
  TTreeFormula* lFVars1 = new TTreeFormula(pExpr1,pExpr1,lTree);
  TString pExpress2 = "fjet1QGtagSub2";
  TString pExpr2(pExpress2);
  TTreeFormula* lFVars2 = new TTreeFormula(pExpr2,pExpr2,lTree);
  TString pExpress3 = "fjet1QGtag";
  TString pExpr3(pExpress3);
  TTreeFormula* lFVars3 = new TTreeFormula(pExpr3,pExpr3,lTree);
  TString pExpress4 = "fjet1PullAngle";
  TString pExpr4(pExpress4);
  TTreeFormula* lFVars4 = new TTreeFormula(pExpr4,pExpr4,lTree);
  TString pExpress5 = "fjet1Pull";
  TString pExpr5(pExpress5);
  TTreeFormula* lFVars5 = new TTreeFormula(pExpr5,pExpr5,lTree);
  TString pExpress6 = "fjet1MassTrimmed";
  TString pExpr6(pExpress6);
  TTreeFormula* lFVars6 = new TTreeFormula(pExpr6,pExpr6,lTree);
  TString pExpress7 = "fjet1MassPruned";
  TString pExpr7(pExpress7);
  TTreeFormula* lFVars7 = new TTreeFormula(pExpr7,pExpr7,lTree);
  TString pExpress8 = "fjet1MassSDbm1";
  TString pExpr8(pExpress8);
  TTreeFormula* lFVars8 = new TTreeFormula(pExpr8,pExpr8,lTree);
  TString pExpress9 = "fjet1MassSDb2";
  TString pExpr9(pExpress9);
  TTreeFormula* lFVars9 = new TTreeFormula(pExpr9,pExpr9,lTree);
  TString pExpress10 = "fjet1MassSDb0";
  TString pExpr10(pExpress10);
  TTreeFormula* lFVars10 = new TTreeFormula(pExpr10,pExpr10,lTree);
  TString pExpress11 = "fjet1QJetVol";
  TString pExpr11(pExpress11);
  TTreeFormula* lFVars11 = new TTreeFormula(pExpr11,pExpr11,lTree);
  TString pExpress12 = "fjet1C2b2";
  TString pExpr12(pExpress12);
  TTreeFormula* lFVars12 = new TTreeFormula(pExpr12,pExpr12,lTree);
  TString pExpress13 = "fjet1C2b1";
  TString pExpr13(pExpress13);
  TTreeFormula* lFVars13 = new TTreeFormula(pExpr13,pExpr13,lTree);
  TString pExpress14 = "fjet1C2b0p5";
  TString pExpr14(pExpress14);
  TTreeFormula* lFVars14 = new TTreeFormula(pExpr14,pExpr14,lTree);
  TString pExpress15 = "fjet1C2b0p2";
  TString pExpr15(pExpress15);
  TTreeFormula* lFVars15 = new TTreeFormula(pExpr15,pExpr15,lTree);
  TString pExpress16 = "fjet1C2b0";
  TString pExpr16(pExpress16);
  TTreeFormula* lFVars16 = new TTreeFormula(pExpr16,pExpr16,lTree);
  TString pExpress17 = "fjet1Tau2";
  TString pExpr17(pExpress17);
  TTreeFormula* lFVars17 = new TTreeFormula(pExpr17,pExpr17,lTree);
  TString pExpress18 = "fjet1Tau1";
  TString pExpr18(pExpress18);
  TTreeFormula* lFVars18 = new TTreeFormula(pExpr18,pExpr18,lTree);
  TString pExpress19 = "fjet1Tau2/fjet1Tau1";
  TString pExpr19(pExpress19);
  TTreeFormula* lFVars19 = new TTreeFormula(pExpr19,pExpr19,lTree);

  //lTree->SetBranchAddress( "jet1mprune"         , &lJP);
  //lTree->SetBranchAddress( iVar1.c_str()           , &lJT1);
  //if(iVar1 != iVar2) lTree->SetBranchAddress( iVar2.c_str()           , &lJT2); 

  int lNEvents = lTree->GetEntries();
  TFile *lOFile = new TFile("Output.root","RECREATE");
  TTree *lOTree = new TTree("DMSTree","DMSTree");
  float lMVA    = 0; lOTree->Branch("bdt_all",&lMVA ,"bdt_all/F");
  for (Long64_t i0=0; i0<lNEvents;i0++) {
    if (i0 % 10000 == 0) std::cout << "--- ... Processing event: " << double(i0)/double(lNEvents) << std::endl;
    lTree->GetEntry(i0);
    jet1QGtagComb = lFVars0->EvalInstance();
    fjet1QGtagSub1 = lFVars1->EvalInstance();
    fjet1QGtagSub2 = lFVars2->EvalInstance();
    fjet1QGtag = lFVars3->EvalInstance();
    fjet1PullAngle = lFVars4->EvalInstance();
    fjet1Pull = lFVars5->EvalInstance();
    fjet1MassTrimmed = lFVars6->EvalInstance();
    fjet1MassPruned = lFVars7->EvalInstance();
    fjet1MassSDbm1 = lFVars8->EvalInstance();
    fjet1MassSDb2 = lFVars9->EvalInstance();
    fjet1MassSDb0 = lFVars10->EvalInstance();
    fjet1QJetVol = lFVars11->EvalInstance();
    fjet1C2b2 = lFVars12->EvalInstance();
    fjet1C2b1 = lFVars13->EvalInstance();
    fjet1C2b0p5 = lFVars14->EvalInstance();
    fjet1C2b0p2 = lFVars15->EvalInstance();
    fjet1C2b0 = lFVars16->EvalInstance();
    fjet1Tau2 = lFVars17->EvalInstance();
    fjet1Tau1 = lFVars18->EvalInstance();
    tau2tau1 = lFVars19->EvalInstance();

    lMVA      = float(reader->EvaluateMVA(lJetName.c_str()));
    lOTree->Fill();
  }
  lOTree->Write();
  lOFile->Close();
  delete reader;
}

void PlotDecisionBoundary( TString weightFile = "weights/TMVAClassification_BDT.weights.xml",TString v0="var0", TString v1="var1", TString dataFileName = "/home/hvoss/TMVA/TMVA_data/data/data_circ.root") 
{   
   //---------------------------------------------------------------
   // default MVA methods to be trained + tested

   // this loads the library
   TMVA::Tools::Instance();

   std::cout << std::endl;
   std::cout << "==> Start TMVAClassificationApplication" << std::endl;


   //
   // create the Reader object
   //
   TMVA::Reader *reader = new TMVA::Reader( "!Color:!Silent" );    

   // create a set of variables and declare them to the reader
   // - the variable names must corresponds in name and type to 
   // those given in the weight file(s) that you use
   Float_t var0, var1;
   reader->AddVariable( v0,                &var0 );
   reader->AddVariable( v1,                &var1 );

   //
   // book the MVA method
   //
   reader->BookMVA( "M1", weightFile ); 
   
   TFile *f = new TFile(dataFileName);
   TTree *signal     = (TTree*)f->Get("TreeS");
   TTree *background = (TTree*)f->Get("TreeB");


//Declaration of leaves types
   Float_t         svar0;
   Float_t         svar1;
   Float_t         bvar0;
   Float_t         bvar1;
   Float_t         sWeight=1.0; // just in case you have weight defined, also set these branchaddresses
   Float_t         bWeight=1.0*signal->GetEntries()/background->GetEntries(); // just in case you have weight defined, also set these branchaddresses

   // Set branch addresses.
   signal->SetBranchAddress(v0,&svar0);
   signal->SetBranchAddress(v1,&svar1);
   background->SetBranchAddress(v0,&bvar0);
   background->SetBranchAddress(v1,&bvar1);


   UInt_t nbin = 50;
   Float_t xmax = signal->GetMaximum(v0.Data());
   Float_t xmin = signal->GetMinimum(v0.Data());
   Float_t ymax = signal->GetMaximum(v1.Data());
   Float_t ymin = signal->GetMinimum(v1.Data());
 
   xmax = TMath::Max(xmax,background->GetMaximum(v0.Data()));  
   xmin = TMath::Min(xmin,background->GetMinimum(v0.Data()));
   ymax = TMath::Max(ymax,background->GetMaximum(v1.Data()));
   ymin = TMath::Min(ymin,background->GetMinimum(v1.Data()));


   TH2D *hs=new TH2D("hs","",nbin,xmin,xmax,nbin,ymin,ymax);   
   TH2D *hb=new TH2D("hb","",nbin,xmin,xmax,nbin,ymin,ymax);   
   hs->SetXTitle(v0);
   hs->SetYTitle(v1);
   hb->SetXTitle(v0);
   hb->SetYTitle(v1);
   hs->SetMarkerColor(4);
   hb->SetMarkerColor(2);


   TH2F * hist = new TH2F( "MVA",    "MVA",    nbin,xmin,xmax,nbin,ymin,ymax);

   // Prepare input tree (this must be replaced by your data source)
   // in this example, there is a toy tree with signal and one with background events
   // we'll later on use only the "signal" events for the test in this example.

   Float_t MinMVA=10000, MaxMVA=-100000;
   for (Int_t ibin=1; ibin<nbin+1; ibin++){
      for (Int_t jbin=1; jbin<nbin+1; jbin++){
         var0 = hs->GetXaxis()->GetBinCenter(ibin);
         var1 = hs->GetYaxis()->GetBinCenter(jbin);
         Float_t mvaVal=reader->EvaluateMVA( "M1" ) ;
         if (MinMVA>mvaVal) MinMVA=mvaVal;
         if (MaxMVA<mvaVal) MaxMVA=mvaVal;
         hist->SetBinContent(ibin,jbin, mvaVal);
      }
   }

   // creating a fine histograms containing the error rate
   const Int_t nValBins=100;
   Double_t    sum = 0.;

   TH1F *mvaS= new TH1F("mvaS","",nValBins,MinMVA,MaxMVA);
   TH1F *mvaB= new TH1F("mvaB","",nValBins,MinMVA,MaxMVA);
   TH1F *mvaSC= new TH1F("mvaSC","",nValBins,MinMVA,MaxMVA);
   TH1F *mvaBC= new TH1F("mvaBC","",nValBins,MinMVA,MaxMVA);

   Long64_t nentries;
   nentries = TreeS->GetEntries();
//.........这里部分代码省略.........

//.........这里部分代码省略.........
		    //cout<<FATnJet<<endl;
		    for (FATi=0;FATi<FATnJet;FATi++){
		      //cout<<"FATi="<<FATi<<endl;
		      if(FATjetCISVV2[FATi]<0 ||FATjetCISVV2[FATi]>1 )continue;
		      FATjetP4_1 = (TLorentzVector*)FATjetP4->At(FATi);
		      bool isOverlap=0;
		      for(int i=0;i<mus.size();i++){
			TLorentzVector* thisMu =(TLorentzVector*)muP4->At(mus[i]) ;
			if(FATjetP4_1->DeltaR(*thisMu)<0.8){
			  isOverlap=1;
			  break;
			}
		      }
		      if(!isOverlap){
			for(int i=0;i<eles.size();i++){
			  TLorentzVector* thisEle =(TLorentzVector*)eleP4->At(eles[i]) ;
			  if(FATjetP4_1->DeltaR(*thisEle)<0.8){
			    isOverlap=1;
			    break;
			  }
			}
			
		      }
		      if(isOverlap)continue;
		      
		      isFAT=1;
		    break;
		    }
		    if(!isFAT)continue;
		    
		    if(FATjetP4_1->Pt()<200)continue;
		    
		    Int_t* FATnSubSDJet=data.GetPtrInt("FATnSubSDJet");
		    
		    if(FATnSubSDJet[FATi]<2)continue;
		    
		    int nsub=FATnSubSDJet[FATi];
		    if(nsub!=2)cout<<"nsub="<<nsub<<endl;if(nsub!=2)cout<<"nsub="<<nsub<<endl;
		    //
		    
		    float* FATjetTau1=data.GetPtrFloat("FATjetTau1");
		    float* FATjetTau2=data.GetPtrFloat("FATjetTau2");
		    vector<float>   *FATsubjetSDCSV =  data.GetPtrVectorFloat("FATsubjetSDCSV");
		    vector<float>   *FATsubjetSDPx =  data.GetPtrVectorFloat("FATsubjetSDPx");
		    vector<float>   *FATsubjetSDPy =  data.GetPtrVectorFloat("FATsubjetSDPy");
		    vector<float>   *FATsubjetSDPz =  data.GetPtrVectorFloat("FATsubjetSDPz");
		    vector<float>   *FATsubjetSDCE =  data.GetPtrVectorFloat("FATsubjetSDCE");
		    
		 	

		    int subi=0,subj=0;
		    bool isSubi=0,isSubj=0;
		    
	         
		    for(int subij=0;subij<FATnSubSDJet[FATi];subij++){
		      if(FATsubjetSDCSV[FATi][subij]>1 ||FATsubjetSDCSV[FATi][subij]<0 )continue;
		      if(!isSubi && !isSubj){
			subi=subij;isSubi=1;continue;
		      }
		      if(isSubi && !isSubj){
			subj=subij;isSubj=1;break;
		      }
		    }
		    //
		    if (!isSubi || !isSubj)continue;
		    
		    TLorentzVector  FATsubjet_1,FATsubjet_2;
		    
		    
		    FATsubjet_1.SetPxPyPzE(FATsubjetSDPx[FATi][subi],FATsubjetSDPy[FATi][subi],FATsubjetSDPz[FATi][subi],FATsubjetSDCE[FATi][subi]);
	            FATsubjet_2.SetPxPyPzE(FATsubjetSDPx[FATi][subj],FATsubjetSDPy[FATi][subj],FATsubjetSDPz[FATi][subj],FATsubjetSDCE[FATi][subj]);
	           
		    
		    fatPt=FATjetP4_1->Pt();//BfatPt->Fill();
		    fatCSV=FATjetCISVV2[FATi];//BfatCSV->Fill();
		    sub1Pt=FATsubjet_1.Pt();//Bsub1Pt->Fill();
		    //sub1Eta=FATsubjet_1.Eta();//Bsub1Eta->Fill();
		    sub1CSV=FATsubjetSDCSV[FATi][subi];//Bsub1CSV->Fill();
		    sub2Pt=FATsubjet_2.Pt();//Bsub2Pt->Fill();
		    //sub2Eta=FATsubjet_2.Eta();//Bsub2Eta->Fill();
		    sub2CSV=FATsubjetSDCSV[FATi][subj];//Bsub2CSV->Fill();
		    delta_R=FATsubjet_1.DeltaR(FATsubjet_2);//BdeltaR->Fill();
		    tau21=FATjetTau2[FATi]/FATjetTau1[FATi];
		    tau1=FATjetTau1[FATi];
		    tau2=FATjetTau2[FATi];
		    
		    
		    th1 ->Fill( reader->EvaluateMVA("BDT method"));
		   
		  }
		  th1->Write();
		  
		  fw->Close();
		  
		}	
		
		
    }	    
    
}

//.........这里部分代码省略.........
    currentTree->SetBranchAddress( "pt2", &pt2_ );
    currentTree->SetBranchAddress( "Deta",&Deta_ );
    currentTree->SetBranchAddress( "Mjj", &Mjj_ );
    currentTree->SetBranchAddress( "diTauSVFitPt",&diTauSVFitPt);
    //currentTree->SetBranchAddress( "diTauSVFitEta",&diTauSVFitEta);
    currentTree->SetBranchAddress( "diTauSVFitMass",&diTauSVFitMass);
    currentTree->SetBranchAddress( "diTauVisMass",&diTauVisMass);
    currentTree->SetBranchAddress( "ptL1", &ptL1 );
    currentTree->SetBranchAddress( "ptL2",  &ptL2 );
    currentTree->SetBranchAddress( "etaL1", &etaL1 );
    currentTree->SetBranchAddress( "etaL2", &etaL2 );
    currentTree->SetBranchAddress( "combRelIsoLeg1",&combRelIsoLeg1);
    currentTree->SetBranchAddress( "tightestHPSWP",&tightestHPSWP);
    currentTree->SetBranchAddress( "diTauCharge",&diTauCharge);
    currentTree->SetBranchAddress( "MtLeg1",&MtLeg1);
    currentTree->SetBranchAddress( "numPV",&numPV);
    currentTree->SetBranchAddress( "sampleWeight",&sampleWeight);
    currentTree->SetBranchAddress( "ptVeto",&ptVeto);
    currentTree->SetBranchAddress( "HLT",&HLT);

    for (Long64_t ievt=0; ievt<currentTree->GetEntries();ievt++) {

      currentTree->GetEntry(ievt);

      if (ievt%10000 == 0){
	std::cout << (jt->first) << " ---> processing event: " << ievt << " ..." <<std::endl;
      }

      pt1  = pt1_;
      pt2  = pt2_;
      Deta = Deta_;
      Mjj  = Mjj_;

      bool pass = effS_>0 ? reader->EvaluateMVA( "Cuts", effS_ ) : (pt1>0);

      if(pass){
	tot+=sampleWeight;
	counter++;
	if(ptVeto<20){
	  tot2+=sampleWeight;
	  counter2++;
	  if(HLT>0.5 && HLT<1.5){
	    tot3+=sampleWeight;
	    counter3++;
	  }
	}
      }
    
    }// end loop   

    cutMap_VBF[jt->first] = tot;
    cutMap_VBFE[jt->first] = counter>0 ? sqrt(counter)*tot/counter : 0;
    cutMap_JetVeto[jt->first] = tot2;
    cutMap_JetVetoE[jt->first] = counter2>0 ? sqrt(counter2)*tot2/counter2 : 0;
    cutMap_HLT[jt->first] = tot3;
    cutMap_HLTE[jt->first] = counter3>0 ? sqrt(counter3)*tot3/counter3 : 0;
  }


  
 

  std::vector< std::map<std::string,float> > allFilters;
  allFilters.push_back(cutMap_allEventsFilter);
  allFilters.push_back(cutMap_vertexScrapingFilter);
  allFilters.push_back(cutMap_oneElectronFilter);

//.........这里部分代码省略.........
   theTree->SetBranchAddress( "deltaRphotonmuon", &mydeltaRphotonmuon );
//   theTree->SetBranchAddress( "ht", &myht );
   theTree->SetBranchAddress( "costopphoton", &mycostopphoton );
   theTree->SetBranchAddress( "jetmultiplicity", &myjetmultiplicity );
//   theTree->SetBranchAddress( "bjetmultiplicity", &mybjetmultiplicity );
   theTree->SetBranchAddress( "deltaphiphotonmet", &mydeltaphiphotonmet );
   theTree->SetBranchAddress( "cvsdiscriminant", &mycvsdiscriminant );
//   theTree->SetBranchAddress( "leptoncharge", &myleptoncharge );
   theTree->SetBranchAddress( "weight", &myweight);
   theTree->SetBranchAddress( "btagSF", &mybtagSF);
   theTree->SetBranchAddress( "btagSFup", &mybtagSFup);
   theTree->SetBranchAddress( "btagSFdown", &mybtagSFdown);
   theTree->SetBranchAddress( "mistagSFup", &mymistagSFup);
   theTree->SetBranchAddress( "mistagSFdown", &mymistagSFdown);
   theTree->SetBranchAddress( "triggerSF", &mytriggerSF);
   theTree->SetBranchAddress( "triggerSFup", &mytriggerSFup);
   theTree->SetBranchAddress( "triggerSFdown", &mytriggerSFdown);
   theTree->SetBranchAddress( "photonSF", &myphotonSF);
   theTree->SetBranchAddress( "photonSFup", &myphotonSFup);
   theTree->SetBranchAddress( "photonSFdown", &myphotonSFdown);
   theTree->SetBranchAddress( "muonSF", &mymuonSF);
   theTree->SetBranchAddress( "muonSFup", &mymuonSFup);
   theTree->SetBranchAddress( "muonSFdown", &mymuonSFdown);
   theTree->SetBranchAddress( "pileupSF", &mypileupSF);
   theTree->SetBranchAddress( "pileupSFup", &mypileupSFup);
   theTree->SetBranchAddress( "pileupSFdown", &mypileupSFdown);



 // Efficiency calculator for cut method
   Int_t    nSelCutsGA = 0;
   Double_t effS       = 0.7;

   std::vector<Float_t> vecVar(4); // vector for EvaluateMVA tests

//   std::cout << "--- Processing: " << theTree->GetEntries() << " events" << std::endl;
   TStopwatch sw;
   sw.Start();
   for (Long64_t ievt=0; ievt<theTree->GetEntries();ievt++) {
//   std::cout << "--- ... Processing event: " << ievt << std::endl;
double finalweight;

      if (ievt%1000 == 0) std::cout << "--- ... Processing event: " << ievt << std::endl;

      theTree->GetEntry(ievt);
//for (int l=0;l<sizeof(myptphoton);l++){
//std::cout << "--- ... reza: " << myptphoton[l] <<std::endl;
//}
//std::cout << "--- ......................."<< (*mycvsdiscriminant)[0]<<std::endl;
      // --- Return the MVA outputs and fill into histograms
ptphoton=(float)(*myptphoton)[0];
etaphoton=(float)(*myetaphoton )[0];
ptmuon=(float)(*myptmuon )[0];
etamuon=(float)(*myetamuon )[0];
ptjet=(float)(*myptjet )[0];
etajet=(float)(*myetajet )[0];
masstop=(float)(*mymasstop )[0];
//mtw=(float)(*mymtw )[0];
deltaRphotonjet=(float)(*mydeltaRphotonjet )[0];
deltaRphotonmuon=(float)(*mydeltaRphotonmuon )[0];
//ht=(float)(*myht )[0];
costopphoton=(float)(*mycostopphoton )[0];
jetmultiplicity=(float)(*myjetmultiplicity )[0];
//bjetmultiplicity=(float)(*mybjetmultiplicity )[0];
deltaphiphotonmet=(float)(*mydeltaphiphotonmet )[0];
cvsdiscriminant=(float)(*mycvsdiscriminant)[0];

//.........这里部分代码省略.........
    outTree[i] -> Branch("dphillmet",&dphillmet);  
    outTree[i] -> Branch("trkMet",&trkMet);     
    //outTree[i] -> Branch("Mt1",&Mt1);        
    //outTree[i] -> Branch("Mt2",&Mt2);        
    outTree[i] -> Branch("mpmet",&mpmet);      
    //outTree[i] -> Branch("Mc",&Mc);         
    outTree[i] -> Branch("ptWW",&ptWW);       
    outTree[i] -> Branch("Ht",&Ht);         

    //Applying Selections
    Int_t cont = 0;
    for (int j = 0; j < MyTree[i]->GetEntries(); ++j){
      if (j == 0) cout<<MyName[i]<<": "<<MyTree[i]->GetEntries()<<endl;
      MyTree[i]->GetEntry(j);
      
      if (pt1        < cutpt1)        continue;
      if (pt2        < cutpt2)        continue;
      if (ptll       < cutptll)       continue;
      if (mll        < cutmll)        continue;
      if (mth        < cutmth)        continue;
      if (pfType1Met < cutpfType1Met) continue;
      if (drll       < cutdrll)       continue;
      if (dphill     < cutdphill)     continue;
      if (dphilljet  < cutdphilljet)  continue;
      if (dphillmet  < cutdphillmet)  continue;
      if (trkMet     < cuttrkMet)     continue;
      //      if (Mt1        < cutMt1)        continue;
      //if (Mt2        < cutMt2)        continue;
      if (mpmet      < cutmpmet)      continue;
      //if (Mc         < cutMc)         continue;
      if (ptWW       < cutptWW)       continue;
      if (Ht         < cutHt)         continue;
      
      value = reader->EvaluateMVA("BDT");
      if(value       < cutValue)      continue;
      
      ++cont;
      outTree[i]->Fill();
    }
    cout<<MyName[i]<<" survived: "<<cont<<" ("<<100 * cont / MyTree[i]->GetEntries()<<"%)"<<endl;  
  }

  /*
  //applying selections on ZH sample
  Int_t contZH = 0;
  for (int i = 0; i < ZHTree->GetEntries(); ++i){
    if (i == 0) cout<<"ZH Entries : "<<ZHTree->GetEntries()<<endl;
    ZHTree->GetEntry(i);

    if (pt1        < cutpt1)        continue;
    if (pt2        < cutpt2)        continue;
    if (ptll       < cutptll)       continue;
    if (mll        < cutmll)        continue;
    if (mth        < cutmth)        continue;
    if (pfType1Met < cutpfType1Met) continue;
    if (drll       < cutdrll)       continue;
    if (dphill     < cutdphill)     continue;
    if (dphilljet  < cutdphilljet)  continue;
    if (dphillmet  < cutdphillmet)  continue;
    if (trkMet     < cuttrkMet)     continue;
    if (Mt1        < cutMt1)        continue;
    if (Mt2        < cutMt2)        continue;
    if (mpmet      < cutmpmet)      continue;
    if (Mc         < cutMc)         continue;
    if (ptWW       < cutptWW)       continue;
    if (Ht         < cutHt)         continue;

//.........这里部分代码省略.........
    double dPhiJPh2 = TMath::ACos(TMath::Cos(ph2scphi - jetleadNoIDphi));
    double dRJPh2 = TMath::Sqrt(TMath::Power(dEtaJPh2, 2) +
                                TMath::Power(dPhiJPh2, 2));
    double dPhiMetGG = TMath::ACos(TMath::Cos(phigg - corrpfmetphi));
    double dPhiMetJet = TMath::ACos(
             TMath::Cos(TMath::Abs(jetleadNoIDphi - corrpfmetphi))
             );
    if (TMath::Abs(ph1sceta) < 1.4442 &&
        TMath::Abs(ph2sceta) < 1.4442 &&
        corrpfmet > 70. &&
        ph1pt/mass > 45./120. &&
        dPhiMetGG > 2.1 &&
        (
          jetleadNoIDpt < 50. ||
          dRJPh1 < 0.5 ||
          dRJPh2 < 0.5 ||
          dPhiMetJet < 2.7
        ) &&
        ph2pt > mass/4) {
       vhMet = 1;
    }

    // Calculate needed variables for the diphoMVA
    if (smearMassError) {
      rVtxSigmaMoM = masserr / mass;          // with smearing
      wVtxSigmaMoM = masserrwvtx / mass;      // with smearing
    } else {
      rVtxSigmaMoM = masserr_ns / mass;       // no smearing
      wVtxSigmaMoM = masserrwvtx_ns / mass;   // no smearing
    }
    cosDPhi = TMath::Cos(phi1 - phi2);
    pho1_ptOverM = ph1pt / mass;
    pho2_ptOverM = ph2pt / mass;
    diphoMVA = reader->EvaluateMVA("BDTG");

    bool passPreselection = (mass > 100 &&
                             mass < 180 &&
                             ph1pt > mass/3 &&
                             ph2pt > mass/4 &&
                             idmva_1 > -0.2 &&
                             idmva_2 > -0.2);

    if (passPreselection == false) {
      if (debug) {
        cout << "    passPreselection: " << passPreselection << endl;
      }
      continue;
    }

    if (debug) {
      cout << "    ... passed preselection." << endl;
    }

    eventCounter++;
    if      (tth   == 1) tth = 2;
    else if (tth   == 2) tth = 1;

    if      (vhHad == 2) vhHad = 1;

    cat = kIncl0;
    if      (tth   == 2) cat = kTTHLep;
    else if (vhLep == 2) cat = kVHLepTight;
    else if (vhLep == 1) cat = kVHLepLoose;
    else if (vbf   >  0) cat = kDijet0;
    else if (vhMet == 1) cat = kVHMet;
    else if (tth   == 1) cat = kTTHHad;

int testPyKerasRegression(){
   // Get data file
   std::cout << "Get test data..." << std::endl;
   TString fname = "./tmva_reg_example.root";
   if (gSystem->AccessPathName(fname))  // file does not exist in local directory
      gSystem->Exec("curl -O http://root.cern.ch/files/tmva_reg_example.root");
   TFile *input = TFile::Open(fname);

   // Build model from python file
   std::cout << "Generate keras model..." << std::endl;
   UInt_t ret;
   ret = gSystem->Exec("echo '"+pythonSrc+"' > generateKerasModelRegression.py");
   if(ret!=0){
       std::cout << "[ERROR] Failed to write python code to file" << std::endl;
       return 1;
   }
   ret = gSystem->Exec("python generateKerasModelRegression.py");
   if(ret!=0){
       std::cout << "[ERROR] Failed to generate model using python" << std::endl;
       return 1;
   }

   // Setup PyMVA and factory
   std::cout << "Setup TMVA..." << std::endl;
   TMVA::PyMethodBase::PyInitialize();
   TFile* outputFile = TFile::Open("ResultsTestPyKerasRegression.root", "RECREATE");
   TMVA::Factory *factory = new TMVA::Factory("testPyKerasRegression", outputFile,
      "!V:Silent:Color:!DrawProgressBar:AnalysisType=Regression");

   // Load data
   TMVA::DataLoader *dataloader = new TMVA::DataLoader("datasetTestPyKerasRegression");

   TTree *tree = (TTree*)input->Get("TreeR");
   dataloader->AddRegressionTree(tree);

   dataloader->AddVariable("var1");
   dataloader->AddVariable("var2");
   dataloader->AddTarget("fvalue");

   dataloader->PrepareTrainingAndTestTree("",
      "SplitMode=Random:NormMode=NumEvents:!V");

   // Book and train method
   factory->BookMethod(dataloader, TMVA::Types::kPyKeras, "PyKeras",
      "!H:!V:VarTransform=D,G:FilenameModel=kerasModelRegression.h5:FilenameTrainedModel=trainedKerasModelRegression.h5:NumEpochs=10:BatchSize=32:SaveBestOnly=false:Verbose=0");
   std::cout << "Train model..." << std::endl;
   factory->TrainAllMethods();

   // Clean-up
   delete factory;
   delete dataloader;
   delete outputFile;

   // Setup reader
   UInt_t numEvents = 100;
   std::cout << "Run reader and estimate target of " << numEvents << " events..." << std::endl;
   TMVA::Reader *reader = new TMVA::Reader("!Color:Silent");
   Float_t vars[3];
   reader->AddVariable("var1", vars+0);
   reader->AddVariable("var2", vars+1);
   reader->BookMVA("PyKeras", "datasetTestPyKerasRegression/weights/testPyKerasRegression_PyKeras.weights.xml");

   // Get mean squared error on events
   tree->SetBranchAddress("var1", vars+0);
   tree->SetBranchAddress("var2", vars+1);
   tree->SetBranchAddress("fvalue", vars+2);

   Float_t meanMvaError = 0;
   for(UInt_t i=0; i<numEvents; i++){
      tree->GetEntry(i);
      meanMvaError += std::pow(vars[2]-reader->EvaluateMVA("PyKeras"),2);
   }
   meanMvaError = meanMvaError/float(numEvents);

   // Check whether the response is obviously better than guessing
   std::cout << "Mean squared error: " << meanMvaError << std::endl;
   if(meanMvaError > 30.0){
      std::cout << "[ERROR] Mean squared error is " << meanMvaError << " (>30.0)" << std::endl;
      return 1;
   }

   return 0;
}