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

//.........这里部分代码省略.........
	reader->AddVariable( "qtb1",                &qtb1 );
	reader->AddVariable( "nSV",                &nSV );
	reader->AddVariable( "mu0pt",                &mu0pt );
	reader->AddVariable( "mu1pt",                &mu1pt );
	reader->AddVariable( "PtbalZH",                &PtbalZH );
	reader->AddVariable( "Angle",                &Angle );
	reader->AddVariable( "Centrality",                &Centrality );
	reader->AddVariable( "MET",                &MET );
	reader->AddVariable( "EvntShpAplanarity",                &EvntShpAplanarity );
	
   // Spectator variables declared in the training have to be added to the reader, too
  Float_t UnweightedEta, mu0pt;
//   reader->AddSpectator( "UnweightedEta",   &UnweightedEta );
  // reader->AddSpectator( "mu0pt",   &mu0pt );

/*   Float_t Category_cat1, Category_cat2, Category_cat3;
   if (Use["Category"]){
      // Add artificial spectators for distinguishing categories
      reader->AddSpectator( "Category_cat1 := var3<=0",             &Category_cat1 );
      reader->AddSpectator( "Category_cat2 := (var3>0)&&(var4<0)",  &Category_cat2 );
      reader->AddSpectator( "Category_cat3 := (var3>0)&&(var4>=0)", &Category_cat3 );
   }
*/
   // --- Book the MVA methods

   TString dir    = "weights/";
   TString prefix = "TMVAClassification";

   // Book method(s)
   for (std::map<std::string,int>::iterator it = Use.begin(); it != Use.end(); it++) {
      if (it->second) {
         TString methodName = TString(it->first) + TString(" method");
         TString weightfile = dir + prefix + TString("_") + TString(it->first) + TString(".weights.xml");
         reader->BookMVA( methodName, weightfile ); 
      }
   }
   
	// Book output histograms
   	TH1F* hCHFb0_OpenSelection= new TH1F		("hCHFb0_OpenSelection", "charged Hadron Energy Fraction b1", 40, 0.0, 1.2);
   	TH1F* hCHFb1_OpenSelection= new TH1F		("hCHFb1_OpenSelection", "charged Hadron Energy Fraction b2", 40, 0.0, 1.2);
 	TH1F* hPtjj_OpenSelection= new TH1F		("hPtjj_OpenSelection","Pt of two b jets with highest CSV ", 50, 0.0, 400);
 	TH1F* hPtmumu_OpenSelection= new TH1F	("hPtmumu_OpenSelection","Pt of two muons with highest pt ", 50, 0.0, 400);
	TH1F* hPtbalZH_OpenSelection= new TH1F	("hPtbalZH_OpenSelection", "Pt balance of Z and H", 40, -80, 80);
	TH1F* hPtmu0_OpenSelection= new TH1F		("hPtmu0_OpenSelection","Pt of muon with highest pt ", 30, 0.0, 300);
	TH1F* hPtmu1_OpenSelection= new TH1F		("hPtmu1_OpenSelection","Pt of muon with second highest pt ", 30, 0.0, 300);
	TH1F* hPFRelIsomu0_OpenSelection= new TH1F		("hPFRelIsomu0_OpenSelection", "PF Rel Iso of muon with highest Pt", 40, 0., 0.2);
	TH1F* hPFRelIsomu1_OpenSelection= new TH1F		("hPFRelIsomu1_OpenSelection", "PF Rel Iso of muon with second highest Pt", 40, 0., 0.2);
	TH1F* hCSV0_OpenSelection= new TH1F		("hCSV0_OpenSelection","Jet with highest CSV ",			40, 0, 1.5);
	TH1F* hCSV1_OpenSelection= new TH1F		("hCSV1_OpenSelection","Jet with second highest CSV ",		40, 0, 1.5);
	TH1F* hdphiVH_OpenSelection= new TH1F	("hdphiVH_OpenSelection","Delta phi between Z and Higgs ", 50, -0.1, 4.5);
	TH1F* hdetaJJ_OpenSelection= new TH1F	("hdetaJJ_OpenSelection","Delta eta between two jets ", 60, -4, 4);
	TH1F* hNjets_OpenSelection= new TH1F	("hNjets_OpenSelection", "Number of Jets",		13, -2.5, 10.5);
	TH1F* hMjj_OpenSelection	= new TH1F	("hMjj_OpenSelection",  "Invariant Mass of two Jets ",		50, 0, 300);
	TH1F* hMmumu_OpenSelection	= new TH1F	("hMmumu_OpenSelection",  "Invariant Mass of two muons ",	75, 0, 200);
    TH1F* hRMSeta_OpenSelection= new TH1F	("hRMSeta_OpenSelection", "RMS Eta",		30, 0, 3);
    TH1F* hStaDeveta_OpenSelection= new TH1F	("hStaDeveta_OpenSelection", "Standard Deviation Eta",		30, 0, 3);
	TH1F* hUnweightedEta_OpenSelection= new TH1F	("hUnweightedEta_OpenSelection",  "Unweighted Eta ",		50, 0, 15);		
	TH1F* hdphiJJ_vect_OpenSelection= new TH1F	("hdphiJJ_vect_OpenSelection", "Delta phi between two jets",  30, -3.5, 4);
	TH1F* hCircularity_OpenSelection= new TH1F("hCircularity_OpenSelection","EventShapeVariables circularity", 30, 0.0, 1.2);
	TH1F* hHt_OpenSelection= new TH1F("hHt_OpenSelection","scalar sum of pt of four particles", 50, 0.0, 500);
    TH1F* hCentrality_OpenSelection= new TH1F	("hCentrality_OpenSelection", "Centrality", 40, 0.0, 0.8);
    TH1F* hEventPt_OpenSelection= new TH1F	("hEventPt_OpenSelection", "Pt of HV system", 50, 0.0, 100);
    TH1F* hAngle_OpenSelection= new TH1F	("hAngle_OpenSelection", "Angle between H and Z", 45, 0, 4.5);
    TH1F* hSphericity_OpenSelection= new TH1F	("hSphericity_OpenSelection", "EventShapeVariables sphericity", 50, 0.0, 1);
    TH1F* hAplanarity_OpenSelection= new TH1F	("hAplanarity_OpenSelection", "EventShapeVariables Aplanarity", 50, -0.1, .4);
    TH1F* hIsotropy_OpenSelection= new TH1F	("hIsotropy_OpenSelection",  "EventShapeVariables isotropy", 30, 0.0, 1.3);

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;
}

//.........这里部分代码省略.........
  theTree->SetBranchAddress("phigg",&phigg);
  theTree->SetBranchAddress("jetleadNoIDpt",&jetleadNoIDpt);
  theTree->SetBranchAddress("jetleadNoIDphi",&jetleadNoIDphi);
  theTree->SetBranchAddress("jetleadNoIDeta",&jetleadNoIDeta);

  theTree->SetBranchAddress("ph1.sceta",&ph1sceta);
  theTree->SetBranchAddress("ph1.scphi",&ph1scphi);

  theTree->SetBranchAddress("ph2.sceta",&ph2sceta);
  theTree->SetBranchAddress("ph2.scphi",&ph2scphi);
 

  // Setup the diphoton BDT
  Float_t rVtxSigmaMoM, wVtxSigmaMoM, cosDPhi;
  Float_t pho1_ptOverM;
  Float_t pho2_ptOverM;
  Float_t diphoMVA;
  
  TMVA::Reader* reader = new TMVA::Reader("Silent");
  reader->AddVariable("masserrsmeared/mass"        , &rVtxSigmaMoM);
  reader->AddVariable("masserrsmearedwrongvtx/mass", &wVtxSigmaMoM);
  reader->AddVariable("vtxprob"                    , &vtxprob     );
  reader->AddVariable("ph1.pt/mass"                , &pho1_ptOverM);
  reader->AddVariable("ph2.pt/mass"                , &pho2_ptOverM);
  reader->AddVariable("ph1.eta"                    , &teta1       );
  reader->AddVariable("ph2.eta"                    , &teta2       );
  reader->AddVariable("TMath::Cos(ph1.phi-ph2.phi)", &cosDPhi     );
  reader->AddVariable("ph1.idmva"                  , &idmva_1     );
  reader->AddVariable("ph2.idmva"                  , &idmva_2     );
  const char *diphotonWeights = (
    "/home/veverka/cms/cmssw/031/CMSSW_5_3_10_patch1/src/MitPhysics/data/"
    "HggBambu_SMDipho_Oct01_redqcdweightallsigevenbkg_BDTG.weights.xml"
    );
  reader->BookMVA("BDTG", diphotonWeights);

  TRandom3 rng(0);

  int eventCounter=0;

  // Loop over the entries.
  std::cout << "Looping over " << theTree->GetEntries() << " entries." << std::endl;
  for (int i=0; i < theTree->GetEntries(); ++i) {
   
    if (eventCounter > 9 && debug ) break;
    if (debug) {
      cout << "Processing entry " << i << " :" << endl
           << "    mass:   " << mass << endl
           << "    ph1pt:  " << ph1pt << endl
           << "    ph2pt:  " << ph2pt << endl
           << "    idmva_1:" << idmva_1 << endl
           << "    idmva_2:" << idmva_2 << endl;
    }
    
    theTree->GetEntry(i);

    // MET category
    vhMet = 0;
    double dEtaJPh1 = ph1sceta - jetleadNoIDeta;
    double dPhiJPh1 = TMath::ACos(TMath::Cos(ph1scphi - jetleadNoIDphi));
    double dRJPh1 = TMath::Sqrt(TMath::Power(dEtaJPh1, 2) +
                                TMath::Power(dPhiJPh1, 2));
    double dEtaJPh2 = ph2sceta - jetleadNoIDeta;
    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));

//should include inttype 11, 12, 13 as the signal
void TMVAClassificationApplication_cc1presv2_bdt_ver3noveract( TString myMethodList = "", TString fname )
{
#ifdef __CINT__
    gROOT->ProcessLine( ".O0" ); // turn off optimization in CINT
#endif
    
    //---------------------------------------------------------------
    
    // This loads the library
    TMVA::Tools::Instance();
    
    // Default MVA methods to be trained + tested
    std::map<std::string,int> Use;
    //
    // --- Boosted Decision Trees
    Use["BDT"]             = 1; // uses Adaptive Boost
    
    
    std::cout << std::endl;
    std::cout << "==> Start TMVAClassificationApplication" << std::endl;
    
    // Select methods (don't look at this code - not of interest)
    if (myMethodList != "") {
        for (std::map<std::string,int>::iterator it = Use.begin(); it != Use.end(); it++) it->second = 0;
        
        std::vector<TString> mlist = gTools().SplitString( myMethodList, ',' );
        for (UInt_t i=0; i<mlist.size(); i++) {
            std::string regMethod(mlist[i]);
            
            if (Use.find(regMethod) == Use.end()) {
                std::cout << "Method \"" << regMethod
                << "\" not known in TMVA under this name. Choose among the following:" << std::endl;
                for (std::map<std::string,int>::iterator it = Use.begin(); it != Use.end(); it++) {
                    std::cout << it->first << " ";
                }
                std::cout << std::endl;
                return;
            }
            Use[regMethod] = 1;
        }
    }
    
    // --------------------------------------------------------------------------------------------------
    
    // --- 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) used
    Float_t mumucl, pmucl;
    Float_t pang_t, muang_t;
    Float_t veract;
    Float_t ppe, mupe;
    Float_t range, coplanarity;
    Float_t opening;//newadd
    
    reader->AddVariable( "mumucl", &mumucl );
    reader->AddVariable( "pmucl", &pmucl );
    reader->AddVariable( "pang_t", &pang_t );
    reader->AddVariable( "muang_t", &muang_t );
    //reader->AddVariable( "veract", &veract );
    reader->AddVariable( "ppe", &ppe);
    reader->AddVariable( "mupe", &mupe);
    reader->AddVariable( "range", &range);
    reader->AddVariable( "coplanarity", &coplanarity);
    reader->AddVariable( "opening", &opening);//newadd
    
    // Spectator variables declared in the training have to be added to the reader, too
    Int_t fileIndex, inttype;
    Float_t nuE, norm, totcrsne;
    reader->AddSpectator( "fileIndex", &fileIndex );
    reader->AddSpectator( "nuE", &nuE );
    reader->AddSpectator( "inttype", &inttype );
    reader->AddSpectator( "norm", &norm );
    reader->AddSpectator( "totcrsne", &totcrsne );
    reader->AddSpectator( "veract", &veract );
    
    // --- Book the MVA methods
    
    TString dir    = "weights/";
    TString prefix = "TMVAClassification_cc1presv2_ver3noveract";//newchange
    
    // Book method(s)
    for (std::map<std::string,int>::iterator it = Use.begin(); it != Use.end(); it++) {
        if (it->second) {
            TString methodName = TString(it->first) + TString(" method");
            TString weightfile = dir + prefix + TString("_") + TString(it->first) + TString(".weights.xml");
            reader->BookMVA( methodName, weightfile );
        }
    }
    
    // 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.
    //
    
    /*//TString fname = "/home/cvson/cc1picoh/frkikawa/meAna/ip4tmva/data_merged_ccqe_addpid_ver3noveract_pid1pres.root";
    //TString fname = "/home/cvson/cc1picoh/frkikawa/meAna/ip4tmva/pm_merged_ccqe_tot_addpid_ver3noveract_pid1pres.root";
//.........这里部分代码省略.........

void applyBDT(std::string iName="/mnt/hscratch/dabercro/skims2/BDT_Signal.root",
              TString inputVariables = "trainingVars.txt",
              TString inputTree = "DMSTree",
              std::string iWeightFile="weights/TMVAClassificationCategory_BDT_simple_alpha.weights.xml",
              TString outName = "Output.root",
              TString fOutputName = "TMVA.root",
              TString fMethodName = "BDT",
              TString fUniformVariable = "fjet1MassTrimmed",
              TString fWeight = "fuck",
              Int_t NumBins = 10, Double_t VarMin = 10, Double_t VarMax = 190, Int_t NumMapPoints = 501) {

  Double_t binWidth = (VarMax - VarMin)/NumBins;
  Double_t VarVals[NumBins+1];
  for (Int_t i0 = 0; i0 < NumBins + 1; i0++)
    VarVals[i0] = VarMin + i0 * binWidth;

  std::vector<TGraph*> transformGraphs;
  TGraph *tempGraph;

  if (fUniformVariable != "") {
    // First scale the BDT to be uniform
    // Make the background shape

    TFile *TMVAFile = TFile::Open(fOutputName);
    TTree *BackgroundTree = (TTree*) TMVAFile->Get("TrainTree");

    mithep::PlotHists *HistPlotter = new mithep::PlotHists();
    HistPlotter->SetDefaultTree(BackgroundTree);
    HistPlotter->SetDefaultExpr(fMethodName);

    Double_t binWidth = 2.0/(NumMapPoints - 1);
    Double_t BDTBins[NumMapPoints];
    for (Int_t i0 = 0; i0 < NumMapPoints; i0++)
      BDTBins[i0] = i0 * binWidth - 1;

    for (Int_t iVarBin = 0; iVarBin < NumBins; iVarBin++) {
      TString BinCut = TString::Format("%s*(%s>=%f&&%s<%f)", fWeight.Data(),
                                       fUniformVariable.Data(), VarVals[iVarBin],
                                       fUniformVariable.Data(), VarVals[iVarBin+1]);
      HistPlotter->AddWeight(fWeight + TString("*(classID == 1 && ") + BinCut + ")");
    }

    std::vector<TH1D*> BDTHists = HistPlotter->MakeHists(NumMapPoints,-1,1);

    for (Int_t iVarBin = 0; iVarBin < NumBins; iVarBin++) {
      tempGraph = new TGraph(NumMapPoints);
      transformGraphs.push_back(tempGraph);
      Double_t FullIntegral = BDTHists[iVarBin]->Integral();
      for (Int_t iMapPoint = 0; iMapPoint < NumMapPoints; iMapPoint++) {
        transformGraphs[iVarBin]->SetPoint(iMapPoint, BDTBins[iMapPoint],
                                           BDTHists[iVarBin]->Integral(0,iMapPoint)/FullIntegral);
      }
    }

    for (UInt_t iHist = 0; iHist < BDTHists.size(); iHist++)
      delete BDTHists[iHist];

    TMVAFile->Close();
  }

  TMVA::Tools::Instance();
  TMVA::Reader *reader = new TMVA::Reader( "!Color:!Silent" );    

  TString BDTName;
  ifstream configFile;
  configFile.open(inputVariables.Data());
  TString tempFormula;

  std::vector<float> Evaluated;
  std::vector<TString> Strings;
  TTreeFormula *Formulas[40];
  
  configFile >> BDTName;  // Is the name of the BDT

  while(!configFile.eof()){
    configFile >> tempFormula;
    if(tempFormula != ""){
      Evaluated.push_back(0.);
      Strings.push_back(tempFormula);
    }
  }
  
  TFile *lFile = new TFile(iName.c_str());
  TTree *lTree = (TTree*) lFile->FindObjectAny(inputTree);

  if(lTree->GetBranch(BDTName) == NULL){

    for(unsigned int i0 = 0;i0 < Strings.size();i0++){
      if (i0 == 0)
        reader->AddSpectator(Strings[i0],&Evaluated[i0]);
      else
        reader->AddVariable(Strings[i0],&Evaluated[i0]);
      Formulas[i0] = new TTreeFormula(Strings[i0],Strings[i0],lTree);
    }
    
    std::string lJetName = "BDT";
    reader->BookMVA(lJetName .c_str(),iWeightFile.c_str());
    
    int lNEvents = lTree->GetEntries();
    TFile *lOFile = new TFile(outName,"RECREATE");
//.........这里部分代码省略.........

//.........这里部分代码省略.........
//  reader->AddVariable ("topphotonmass", &topphotonmass);
//reader->AddVariable ("pttop", &pttop);
//reader->AddVariable ("etatop", &etatop);
  reader->AddVariable ("jetmultiplicity", &jetmultiplicity);
//  reader->AddVariable ("bjetmultiplicity", &bjetmultiplicity);
  reader->AddVariable ("deltaphiphotonmet", &deltaphiphotonmet);
  reader->AddVariable ("cvsdiscriminant", &cvsdiscriminant);
//  reader->AddVariable ("leptoncharge", &leptoncharge);


/*
   // Spectator variables declared in the training have to be added to the reader, too
   reader->AddSpectator( "spec1 := var1*2",   &spec1 );
   reader->AddSpectator( "spec2 := var1*3",   &spec2 );

   Float_t Category_cat1, Category_cat2, Category_cat3;
   if (Use["Category"]){
      // Add artificial spectators for distinguishing categories
      reader->AddSpectator( "Category_cat1 := var3<=0",             &Category_cat1 );
      reader->AddSpectator( "Category_cat2 := (var3>0)&&(var4<0)",  &Category_cat2 );
      reader->AddSpectator( "Category_cat3 := (var3>0)&&(var4>=0)", &Category_cat3 );
   }
*/
   // --- Book the MVA methods

   TString dir    = "weights/";
   TString prefix = "TMVA";

   // Book method(s)
   for (std::map<std::string,int>::iterator it = Use.begin(); it != Use.end(); it++) {
      if (it->second) {
         TString methodName = TString(it->first) + TString(" method");
         TString weightfile = dir + prefix + TString("_") + TString(it->first) + TString(".weights.xml");
         reader->BookMVA( methodName, weightfile ); 
      }
   }
   
   // Book output histograms
   UInt_t nbin = 40;
   TH1F   *histLk(0), *histLkD(0), *histLkPCA(0), *histLkKDE(0), *histLkMIX(0), *histPD(0), *histPDD(0);
   TH1F   *histPDPCA(0), *histPDEFoam(0), *histPDEFoamErr(0), *histPDEFoamSig(0), *histKNN(0), *histHm(0);
   TH1F   *histFi(0), *histFiG(0), *histFiB(0), *histLD(0), *histNn(0),*histNnbfgs(0),*histNnbnn(0);
   TH1F   *histNnC(0), *histNnT(0), *histBdt(0), *histBdtG(0), *histBdtD(0), *histRf(0), *histSVMG(0);
   TH1F   *histSVMP(0), *histSVML(0), *histFDAMT(0), *histFDAGA(0), *histCat(0), *histPBdt(0);

   if (Use["Likelihood"])    histLk      = new TH1F( "MVA_Likelihood",    "MVA_Likelihood",    nbin, -1, 1 );
   if (Use["LikelihoodD"])   histLkD     = new TH1F( "MVA_LikelihoodD",   "MVA_LikelihoodD",   nbin, -1, 0.9999 );
   if (Use["LikelihoodPCA"]) histLkPCA   = new TH1F( "MVA_LikelihoodPCA", "MVA_LikelihoodPCA", nbin, -1, 1 );
   if (Use["LikelihoodKDE"]) histLkKDE   = new TH1F( "MVA_LikelihoodKDE", "MVA_LikelihoodKDE", nbin,  -0.00001, 0.99999 );
   if (Use["LikelihoodMIX"]) histLkMIX   = new TH1F( "MVA_LikelihoodMIX", "MVA_LikelihoodMIX", nbin,  0, 1 );
   if (Use["PDERS"])         histPD      = new TH1F( "MVA_PDERS",         "MVA_PDERS",         nbin,  0, 1 );
   if (Use["PDERSD"])        histPDD     = new TH1F( "MVA_PDERSD",        "MVA_PDERSD",        nbin,  0, 1 );
   if (Use["PDERSPCA"])      histPDPCA   = new TH1F( "MVA_PDERSPCA",      "MVA_PDERSPCA",      nbin,  0, 1 );
   if (Use["KNN"])           histKNN     = new TH1F( "MVA_KNN",           "MVA_KNN",           nbin,  0, 1 );
   if (Use["HMatrix"])       histHm      = new TH1F( "MVA_HMatrix",       "MVA_HMatrix",       nbin, -0.95, 1.55 );
   if (Use["Fisher"])        histFi      = new TH1F( "MVA_Fisher",        "MVA_Fisher",        nbin, -4, 4 );
   if (Use["FisherG"])       histFiG     = new TH1F( "MVA_FisherG",       "MVA_FisherG",       nbin, -1, 1 );
   if (Use["BoostedFisher"]) histFiB     = new TH1F( "MVA_BoostedFisher", "MVA_BoostedFisher", nbin, -2, 2 );
   if (Use["LD"])            histLD      = new TH1F( "MVA_LD",            "MVA_LD",            nbin, -2, 2 );
   if (Use["MLP"])           histNn      = new TH1F( "MVA_MLP",           "MVA_MLP",           nbin, -1.25, 1.5 );
   if (Use["MLPBFGS"])       histNnbfgs  = new TH1F( "MVA_MLPBFGS",       "MVA_MLPBFGS",       nbin, -1.25, 1.5 );
   if (Use["MLPBNN"])        histNnbnn   = new TH1F( "MVA_MLPBNN",        "MVA_MLPBNN",        nbin, -1.25, 1.5 );
   if (Use["CFMlpANN"])      histNnC     = new TH1F( "MVA_CFMlpANN",      "MVA_CFMlpANN",      nbin,  0, 1 );
   if (Use["TMlpANN"])       histNnT     = new TH1F( "MVA_TMlpANN",       "MVA_TMlpANN",       nbin, -1.3, 1.3 );
   if (Use["BDT"])           histBdt     = new TH1F( "MVA_BDT",           "MVA_BDT",           nbin, -0.8, 0.8 );
   if (Use["BDTD"])          histBdtD    = new TH1F( "MVA_BDTD",          "MVA_BDTD",          nbin, -0.8, 0.8 );

void cutFlowStudyMu( TString weightFile = "TMVAClassificationPtOrd_qqH115vsWZttQCD_Cuts.weights.xml",
		     Double_t effS_ = 0.3) 
{
  
  ofstream out("cutFlow-MuTauStream.txt");
  out.precision(4);

  TMVA::Tools::Instance();
  TMVA::Reader *reader = new TMVA::Reader( "!Color:!Silent" );   

  Float_t pt1, pt2;
  Float_t Deta, Mjj;
  Float_t eta1,eta2;

  reader->AddVariable( "pt1", &pt1);
  reader->AddVariable( "pt2", &pt2);
  reader->AddVariable( "Deta",&Deta);
  reader->AddVariable( "Mjj", &Mjj);
  reader->AddSpectator("eta1",&eta1);
  reader->AddSpectator("eta2",&eta2);
  reader->BookMVA( "Cuts", TString("/home/llr/cms/lbianchini/CMSSW_3_9_9/src/Bianchi/TauTauStudies/test/Macro/weights/")+weightFile ); 
 
  TFile *fFullSignalVBF           = new TFile("/data_CMS/cms/lbianchini//MuTauStream2011/treeMuTauStream_VBFH115-Mu-powheg-PUS1.root","READ"); 
  TFile *fFullSignalGGH           = new TFile("/data_CMS/cms/lbianchini//MuTauStream2011/treeMuTauStream_GGFH115-Mu-powheg-PUS1.root","READ");  
  TFile *fFullBackgroundDYTauTau  = new TFile("/data_CMS/cms/lbianchini//MuTauStream2011/treeMuTauStream_DYToTauTau-Mu-20-PUS1.root","READ"); 
  TFile *fFullBackgroundDYMuMu    = new TFile("/data_CMS/cms/lbianchini//MuTauStream2011/treeMuTauStream_DYToMuMu-20-PUS1.root","READ"); 
  TFile *fFullBackgroundWJets     = new TFile("/data_CMS/cms/lbianchini//MuTauStream2011/treeMuTauStream_WJets-Mu-madgraph-PUS1.root","READ"); 
  TFile *fFullBackgroundQCD       = new TFile("/data_CMS/cms/lbianchini//MuTauStream2011/treeMuTauStream_QCDmu.root","READ"); 
  TFile *fFullBackgroundTTbar     = new TFile("/data_CMS/cms/lbianchini//MuTauStream2011/treeMuTauStream_TTJets-Mu-madgraph-PUS1.root","READ"); 
  TFile *fFullBackgroundDiBoson   = new TFile("/data_CMS/cms/lbianchini//MuTauStream2011/treeMuTauStream_DiBoson-Mu.root","READ"); 

  // OpenNTuples
  TString fSignalNameVBF           = "/data_CMS/cms/lbianchini/VbfJetsStudy/OpenNtuples/MuTauStream2011/v2/nTupleVBFH115-Mu-powheg-PUS1_Open_MuTauStream.root";
  TString fSignalNameGGH           = "/data_CMS/cms/lbianchini/VbfJetsStudy/OpenNtuples/MuTauStream2011/v2/nTupleGGFH115-Mu-powheg-PUS1_Open_MuTauStream.root";
  TString fBackgroundNameDYTauTau  = "/data_CMS/cms/lbianchini/VbfJetsStudy/OpenNtuples/MuTauStream2011/v2/nTupleDYToTauTau-Mu-20-PUS1_Open_MuTauStream.root";
  TString fBackgroundNameDYMuMu  = "/data_CMS/cms/lbianchini/VbfJetsStudy/OpenNtuples/MuTauStream2011/v2/nTupleDYToMuMu-20-PUS1_Open_MuTauStream.root";
  TString fBackgroundNameWJets     = "/data_CMS/cms/lbianchini/VbfJetsStudy/OpenNtuples/MuTauStream2011/v2/nTupleWJets-Mu-madgraph-PUS1_Open_MuTauStream.root";
  TString fBackgroundNameQCD       = "/data_CMS/cms/lbianchini/VbfJetsStudy/OpenNtuples/MuTauStream2011/v2/nTupleQCDmu_Open_MuTauStream.root";
  TString fBackgroundNameTTbar     = "/data_CMS/cms/lbianchini/VbfJetsStudy/OpenNtuples/MuTauStream2011/v2/nTupleTTJets-Mu-madgraph-PUS1_Open_MuTauStream.root";
  TString fBackgroundNameDiBoson   = "/data_CMS/cms/lbianchini/VbfJetsStudy/OpenNtuples/MuTauStream2011/v2/nTupleDiBoson-Mu_Open_MuTauStream.root";

  TFile *fSignalVBF(0); 
  TFile *fSignalGGH(0); 
  TFile *fBackgroundDYTauTau(0);
  TFile *fBackgroundDYMuMu(0);
  TFile *fBackgroundWJets(0);
  TFile *fBackgroundQCD(0);
  TFile *fBackgroundTTbar(0);
  TFile *fBackgroundDiBoson(0);
 
  fSignalVBF          = TFile::Open( fSignalNameVBF ); 
  fSignalGGH          = TFile::Open( fSignalNameGGH ); 
  fBackgroundDYTauTau = TFile::Open( fBackgroundNameDYTauTau ); 
  fBackgroundDYMuMu   = TFile::Open( fBackgroundNameDYMuMu ); 
  fBackgroundWJets    = TFile::Open( fBackgroundNameWJets ); 
  fBackgroundQCD      = TFile::Open( fBackgroundNameQCD ); 
  fBackgroundTTbar    = TFile::Open( fBackgroundNameTTbar ); 
  fBackgroundDiBoson  = TFile::Open( fBackgroundNameDiBoson ); 

  if(!fSignalVBF || !fBackgroundDYTauTau || !fBackgroundWJets || !fBackgroundQCD || !fBackgroundTTbar ||
     !fSignalGGH || !fBackgroundDYMuMu || !fBackgroundDiBoson ){
    std::cout << "ERROR: could not open files" << std::endl;
    exit(1);
  }

  TString tree = "outTreePtOrd";

  TTree *signalVBF           = (TTree*)fSignalVBF->Get(tree);
  TTree *signalGGH           = (TTree*)fSignalGGH->Get(tree);
  TTree *backgroundDYTauTau  = (TTree*)fBackgroundDYTauTau->Get(tree);
  TTree *backgroundDYMuMu    = (TTree*)fBackgroundDYMuMu->Get(tree);
  TTree *backgroundWJets     = (TTree*)fBackgroundWJets->Get(tree);
  TTree *backgroundQCD       = (TTree*)fBackgroundQCD->Get(tree);
  TTree *backgroundTTbar     = (TTree*)fBackgroundTTbar->Get(tree);
  TTree *backgroundDiBoson   = (TTree*)fBackgroundDiBoson->Get(tree);

  // here I define the map between a sample name and its tree
  std::map<std::string,TTree*> tMap;
  tMap["ggH115"]=signalGGH;
  tMap["qqH115"]=signalVBF;
  tMap["Ztautau"]=backgroundDYTauTau;
  tMap["Zmumu"]=backgroundDYMuMu;
  tMap["Wjets"]=backgroundWJets;
  tMap["QCD"]=backgroundQCD;
  tMap["TTbar"]=backgroundTTbar;
  tMap["DiBoson"]=backgroundDiBoson;

  std::map<std::string,TTree*>::iterator jt;

  Float_t pt1_, pt2_;
  Float_t Deta_, Mjj_;
  Float_t Dphi,diTauSVFitPt,diTauSVFitEta,diTauVisMass,diTauSVFitMass,ptL1,ptL2,etaL1,etaL2,diTauCharge,MtLeg1,numPV,combRelIsoLeg1,sampleWeight,ptVeto,HLT;
  Int_t tightestHPSWP;

  /////////////////////////////////////////////////////////////////////////////////////////////////////////////////



  // here I choose the order in the stack
  std::vector<string> samples;
//.........这里部分代码省略.........

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
   Double_t var0, var1;
   reader->AddVariable( v0,                &var0 );
   reader->AddVariable( v1,                &var1 );

   //
   // book the MVA method
   //
   reader->BookMVA( "MyMVAMethod", 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,(Float_t)background->GetMaximum(v0.Data()));  
   xmin = TMath::Min(xmin,(Float_t)background->GetMinimum(v0.Data()));
   ymax = TMath::Max(ymax,(Float_t)background->GetMaximum(v1.Data()));
   ymin = TMath::Min(ymin,(Float_t)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 (UInt_t ibin=1; ibin<nbin+1; ibin++){
      for (UInt_t jbin=1; jbin<nbin+1; jbin++){
         var0 = hs->GetXaxis()->GetBinCenter(ibin);
         var1 = hs->GetYaxis()->GetBinCenter(jbin);
         Float_t mvaVal=reader->EvaluateMVA( "MyMVAMethod" ) ;
         if (MinMVA>mvaVal) MinMVA=mvaVal;
         if (MaxMVA<mvaVal) MaxMVA=mvaVal;
         hist->SetBinContent(ibin,jbin, mvaVal);
      }
   }


   // now you need to try to find the MVA-value at which you cut for the plotting of the decision boundary
   // (Use the smallest number of misclassifications as criterion)
   const Int_t nValBins=100;
   Double_t    sum = 0.;
   TH1F *mvaS= new TH1F("mvaS","",nValBins,MinMVA,MaxMVA); mvaS->SetXTitle("MVA-ouput"); mvaS->SetYTitle("#entries");
   TH1F *mvaB= new TH1F("mvaB","",nValBins,MinMVA,MaxMVA); mvaB->SetXTitle("MVA-ouput"); mvaB->SetYTitle("#entries");
   TH1F *mvaSC= new TH1F("mvaSC","",nValBins,MinMVA,MaxMVA); mvaSC->SetXTitle("MVA-ouput"); mvaSC->SetYTitle("cummulation");
   TH1F *mvaBC= new TH1F("mvaBC","",nValBins,MinMVA,MaxMVA); mvaBC->SetXTitle("MVA-ouput"); mvaBC->SetYTitle("cummulation");

   Long64_t nentries;
//.........这里部分代码省略.........

void TMVARegressionApplication( int wMs,int wM, string st,string st2,string option="",TString myMethodList = "" )
{
    //---------------------------------------------------------------
    // This loads the library
    TMVA::Tools::Instance();

    // Default MVA methods to be trained + tested
    std::map<std::string,int> Use;

    // --- Mutidimensional likelihood and Nearest-Neighbour methods
    Use["PDERS"]           = 0;
    Use["PDEFoam"]         = 0;
    Use["KNN"]             = 0;
    //
    // --- Linear Discriminant Analysis
    Use["LD"]		        = 0;
    //
    // --- Function Discriminant analysis
    Use["FDA_GA"]          = 0;
    Use["FDA_MC"]          = 0;
    Use["FDA_MT"]          = 0;
    Use["FDA_GAMT"]        = 0;
    //
    // --- Neural Network
    Use["MLP"]             = 0;
    //
    // --- Support Vector Machine
    Use["SVM"]             = 0;
    //
    // --- Boosted Decision Trees
    Use["BDT"]             = 0;
    Use["BDTG"]            = 1;
    // ---------------------------------------------------------------

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

    // Select methods (don't look at this code - not of interest)
    if (myMethodList != "") {
        for (std::map<std::string,int>::iterator it = Use.begin(); it != Use.end(); it++) it->second = 0;

        std::vector<TString> mlist = gTools().SplitString( myMethodList, ',' );
        for (UInt_t i=0; i<mlist.size(); i++) {
            std::string regMethod(mlist[i]);

            if (Use.find(regMethod) == Use.end()) {
                std::cout << "Method \"" << regMethod << "\" not known in TMVA under this name. Choose among the following:" << std::endl;
                for (std::map<std::string,int>::iterator it = Use.begin(); it != Use.end(); it++) std::cout << it->first << " ";
                std::cout << std::endl;
                return;
            }
            Use[regMethod] = 1;
        }
    }

    // --------------------------------------------------------------------------------------------------

    // --- 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) used
    //Float_t var1, var2;
    //reader->AddVariable( "var1", &var1 );
    //reader->AddVariable( "var2", &var2 );
    Float_t pt_AK8MatchedToHbb,eta_AK8MatchedToHbb,nsv_AK8MatchedToHbb,sv0mass_AK8MatchedToHbb,sv1mass_AK8MatchedToHbb,
            nch_AK8MatchedToHbb,nmu_AK8MatchedToHbb,nel_AK8MatchedToHbb,muenfr_AK8MatchedToHbb,emenfr_AK8MatchedToHbb;
    reader->AddVariable( "pt_AK8MatchedToHbb", &pt_AK8MatchedToHbb );
    reader->AddVariable( "eta_AK8MatchedToHbb", &eta_AK8MatchedToHbb );
    reader->AddVariable( "nsv_AK8MatchedToHbb", &nsv_AK8MatchedToHbb );
    reader->AddVariable( "sv0mass_AK8MatchedToHbb", &sv0mass_AK8MatchedToHbb );
    reader->AddVariable( "sv1mass_AK8MatchedToHbb", &sv1mass_AK8MatchedToHbb );
    reader->AddVariable( "nch_AK8MatchedToHbb", &nch_AK8MatchedToHbb );
    reader->AddVariable( "nmu_AK8MatchedToHbb", &nmu_AK8MatchedToHbb );
    reader->AddVariable( "nel_AK8MatchedToHbb", &nel_AK8MatchedToHbb );
    reader->AddVariable( "muenfr_AK8MatchedToHbb", &muenfr_AK8MatchedToHbb );
    reader->AddVariable( "emenfr_AK8MatchedToHbb", &emenfr_AK8MatchedToHbb );


    // Spectator variables declared in the training have to be added to the reader, too
    Float_t spec1,spec2;
    reader->AddSpectator( "spec1:=n_pv",  &spec1 );
    reader->AddSpectator( "spec2:=msoftdrop_AK8MatchedToHbb",  &spec2 );

    // --- Book the MVA methods

    TString dir    = "weights/";
    TString prefix = "TMVARegression";

    // Book method(s)
    for (std::map<std::string,int>::iterator it = Use.begin(); it != Use.end(); it++) {
        if (it->second) {
            TString methodName = it->first + " method";
            TString weightfile = dir + prefix + "_" + TString(it->first) + ".weights.xml";
            reader->BookMVA( methodName, weightfile );
        }
    }

    TH1* hists[100];
//.........这里部分代码省略.........

void TMVAClassificationCategoryApplication()
{
   // ---------------------------------------------------------------
   // default MVA methods to be trained + tested
   std::map<std::string,int> Use;
   // ---
   Use["LikelihoodCat"] = 1;
   Use["FisherCat"]     = 1;
   // ---------------------------------------------------------------

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

   // --- Create the Reader object

   TMVA::Reader *reader = new TMVA::Reader( "!Color:!Silent" );    

   // Create a set of variables and spectators and declare them to the reader
   // - the variable names MUST corresponds in name and type to those given in the weight file(s) used
   Float_t var1, var2, var3, var4, eta;
   reader->AddVariable( "var1", &var1 );
   reader->AddVariable( "var2", &var2 );
   reader->AddVariable( "var3", &var3 );
   reader->AddVariable( "var4", &var4 );

   reader->AddSpectator( "eta", &eta );

   // --- Book the MVA methods

   for (std::map<std::string,int>::iterator it = Use.begin(); it != Use.end(); it++) {
      if (it->second) {
         TString methodName = it->first + " method";
         TString weightfile = "dataset/weights/TMVAClassificationCategory_" + TString(it->first) + ".weights.xml";
         reader->BookMVA( methodName, weightfile ); 
      }
   }

   // Book output histograms
   UInt_t nbin = 100;
   std::map<std::string,TH1*> hist;
   hist["LikelihoodCat"] = new TH1F( "MVA_LikelihoodCat",   "MVA_LikelihoodCat", nbin, -1, 0.9999 );
   hist["FisherCat"]     = new TH1F( "MVA_FisherCat",       "MVA_FisherCat",     nbin, -4, 4 );

   // 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.
   //
   TString fname = TString(gSystem->DirName(__FILE__) ) + "/data/";
   // if directory data not found try using tutorials dir
   if (gSystem->AccessPathName( fname )) {
      fname = TString(gROOT->GetTutorialsDir()) + "/tmva/data/";
   }
   if (UseOffsetMethod) fname += "toy_sigbkg_categ_offset.root";
   else                 fname += "toy_sigbkg_categ_varoff.root";
   std::cout << "--- TMVAClassificationApp    : Accessing " << fname << "!" << std::endl;
   TFile *input = TFile::Open(fname);
   if (!input) {
      std::cout << "ERROR: could not open data file: " << fname << std::endl;
      exit(1);
   }

   // --- Event loop

   // Prepare the tree
   // - here the variable names have to corresponds to your tree
   // - you can use the same variables as above which is slightly faster,
   //   but of course you can use different ones and copy the values inside the event loop
   //
   TTree* theTree = (TTree*)input->Get("TreeS");
   std::cout << "--- Use signal sample for evalution" << std::endl;
   theTree->SetBranchAddress( "var1", &var1 );
   theTree->SetBranchAddress( "var2", &var2 );
   theTree->SetBranchAddress( "var3", &var3 );
   theTree->SetBranchAddress( "var4", &var4 );

   theTree->SetBranchAddress( "eta",  &eta ); // spectator

   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);

      // --- Return the MVA outputs and fill into histograms

      for (std::map<std::string,int>::iterator it = Use.begin(); it != Use.end(); it++) {
         if (!it->second) continue;
         TString methodName = it->first + " method";
         hist[it->first]->Fill( reader->EvaluateMVA( methodName ) );         
      }

   }
   sw.Stop();
   std::cout << "--- End of event loop: "; sw.Print();

   // --- Write histograms

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

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

   // CMS STATS:
   // cut definitions. Define categories and overall selection.
#include "TMVA_tprime_cuts.C"



   // Add artificial spectators for distinguishing categories
   Float_t Category_mycat1, Category_mycat2, Category_mycat3, Category_mycat4;
   TString sCat1("Category_cat1:=");
   TString sCat2("Category_cat2:=");
   TString sCat3("Category_cat3:=");
   TString sCat4("Category_cat4:=");
   sCat1.Append(cut_os_cat1);
   sCat2.Append(cut_os_cat2);
   sCat3.Append(cut_ss);
   sCat4.Append(cut_tri);
   reader->AddSpectator( sCat1, &Category_mycat1 );
   reader->AddSpectator( sCat2, &Category_mycat2 );
   reader->AddSpectator( sCat3, &Category_mycat3 );
   reader->AddSpectator( sCat4, &Category_mycat4 );


   // --- Book the MVA methods

   TString dir    = "weights/";
   TString prefix = "TMVAClassification";

   // Book method(s)
   for (std::map<std::string,int>::iterator it = Use.begin(); it != Use.end(); it++) {
      if (it->second) {
         TString methodName = TString(it->first) + TString(" method");
         TString weightfile = dir + prefix + TString("_") + TString(it->first) + TString(".weights.xml");
         reader->BookMVA( methodName, weightfile ); 
      }
   }

   
   // Book output histograms
   UInt_t nbin = 100;
   TH1F * histBdt(0);
   TH1F * histCat(0);

   if (Use["BDT"])           histBdt     = new TH1F( "MVA_BDT",           "MVA_BDT",           nbin, -0.8, 0.8 );
   if (Use["Category"])      histCat     = new TH1F( "MVA_Category",      "MVA_Category",      nbin, -2., 2. );


   // 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.
   //   
   TFile *input(0);


   // CMS STATS:
   // Specify files with data, for which you want to compute the classifier values
   TString fname = "./data/pass7_OS_test1/TTbar_MuMu/all.root";   
   //TString fname = "./data/pass7_TRI_test1/TTbar_MuMu/all.root";   


   if (!gSystem->AccessPathName( fname )) 
      input = TFile::Open( fname ); // check if file in local directory exists
   
   if (!input) {
      std::cout << "ERROR: could not open data file" << std::endl;
      exit(1);

void TMVAMulticlassApplication( TString myMethodList = "" )
{
#ifdef __CINT__
   gROOT->ProcessLine( ".O0" ); // turn off optimization in CINT
#endif

   TMVA::Tools::Instance();
   
   //---------------------------------------------------------------
   // default MVA methods to be trained + tested
   std::map<std::string,int> Use;
   Use["MLP"]             = 1;
   Use["BDTG"]            = 1;
   Use["FDA_GA"]          = 0;
   //---------------------------------------------------------------
  
   std::cout << std::endl;
   std::cout << "==> Start TMVAMulticlassApplication" << std::endl; 
   if (myMethodList != "") {
      for (std::map<std::string,int>::iterator it = Use.begin(); it != Use.end(); it++) it->second = 0;

      std::vector<TString> mlist = gTools().SplitString( myMethodList, ',' );
      for (UInt_t i=0; i<mlist.size(); i++) {
         std::string regMethod(mlist[i]);

         if (Use.find(regMethod) == Use.end()) {
            std::cout << "Method \"" << regMethod << "\" not known in TMVA under this name. Choose among the following:" << std::endl;
            for (std::map<std::string,int>::iterator it = Use.begin(); it != Use.end(); it++) std::cout << it->first << " " << std::endl;
            std::cout << std::endl;
            return;
         }
         Use[regMethod] = 1;
      }
   }

   
   // 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 var1, var2, var3, var4;
   reader->AddVariable( "var1", &var1 );
   reader->AddVariable( "var2", &var2 );
   reader->AddVariable( "var3", &var3 );
   reader->AddVariable( "var4", &var4 );

   // book the MVA methods
   TString dir    = "weights/";
   TString prefix = "TMVAMulticlass";
   
   for (std::map<std::string,int>::iterator it = Use.begin(); it != Use.end(); it++) {
      if (it->second) {
        TString methodName = TString(it->first) + TString(" method");
        TString weightfile = dir + prefix + TString("_") + TString(it->first) + TString(".weights.xml"); 
        reader->BookMVA( methodName, weightfile ); 
      }
   }

   // book output histograms
   UInt_t nbin = 100;
   TH1F *histMLP_signal(0), *histBDTG_signal(0), *histFDAGA_signal(0);
   if (Use["MLP"])    
      histMLP_signal    = new TH1F( "MVA_MLP_signal",    "MVA_MLP_signal",    nbin, 0., 1.1 );
   if (Use["BDTG"])
      histBDTG_signal  = new TH1F( "MVA_BDTG_signal",   "MVA_BDTG_signal",   nbin, 0., 1.1 );
   if (Use["FDA_GA"])
      histFDAGA_signal = new TH1F( "MVA_FDA_GA_signal", "MVA_FDA_GA_signal", nbin, 0., 1.1 );

   TFile *input(0); 
   TString fname = "./tmva_example_multiple_background.root";
   if (!gSystem->AccessPathName( fname )) {
      input = TFile::Open( fname ); // check if file in local directory exists
   }
   if (!input) {
      std::cout << "ERROR: could not open data file, please generate example data first!" << std::endl;
      exit(1);
   }
   std::cout << "--- TMVAMulticlassApp : Using input file: " << input->GetName() << std::endl;
   
   // prepare the tree
   // - here the variable names have to corresponds to your tree
   // - you can use the same variables as above which is slightly faster,
   //   but of course you can use different ones and copy the values inside the event loop
  
   TTree* theTree = (TTree*)input->Get("TreeS");
   std::cout << "--- Select signal sample" << std::endl;
   theTree->SetBranchAddress( "var1", &var1 );
   theTree->SetBranchAddress( "var2", &var2 );
   theTree->SetBranchAddress( "var3", &var3 );
   theTree->SetBranchAddress( "var4", &var4 );

   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;
//.........这里部分代码省略.........

TString TMVAPredict(TString method_name, EnumPredictMode predictMode = EnumPredictMode::FINAL)
{
    std::cout << "------------ predict with : " << method_name << " ------ " << std::endl;
    std::vector<std::string> inputNames = {"training","test","check_correlation","check_agreement"};
    std::map<std::string,std::vector<std::string>> varsForInput;

    std::vector<std::string> variableOrder = {"id", "signal", "mass", "min_ANNmuon", "prediction"};
    
    varsForInput["training"].emplace_back ("prediction");
    if (predictMode != EnumPredictMode::INTERMEDIATE)
    {
        varsForInput["training"].emplace_back ("id");
        varsForInput["training"].emplace_back ("signal");
        varsForInput["training"].emplace_back ("mass");
        varsForInput["training"].emplace_back ("min_ANNmuon");

        varsForInput["test"].emplace_back ("prediction");
        varsForInput["test"].emplace_back ("id");

        varsForInput["check_agreement"].emplace_back ("signal");
        varsForInput["check_agreement"].emplace_back ("weight");
        varsForInput["check_agreement"].emplace_back ("prediction");

        varsForInput["check_correlation"].emplace_back ("mass");
        varsForInput["check_correlation"].emplace_back ("prediction");
    }

    
    std::map<std::string,std::vector<std::string>> createForInput;
    createForInput["training"].emplace_back ("root");

    if (predictMode != EnumPredictMode::INTERMEDIATE)
    {
        createForInput["training"].emplace_back ("csv");
        createForInput["test"].emplace_back ("csv");
        createForInput["check_agreement"].emplace_back ("csv");
        createForInput["check_correlation"].emplace_back ("csv");
    }


    // -------- prepare the Reader ------
    TMVA::Tools::Instance();

    std::cout << "==> Start TMVAPredict" << std::endl;
    TMVA::Reader *reader = new TMVA::Reader( "!Color:!Silent" );  


    std::vector<Float_t> variables (variableNames.size ());
    auto itVar = begin (variables);
    for (auto varName : variableNames)
    {
	Float_t* pVar = &(*itVar);
        auto localVarName = varName;
        localVarName.substr(0,localVarName.find(":="));
	reader->AddVariable(varName.c_str(), pVar);
	(*itVar) = 0.0;
	++itVar;
    }

    // spectators not known for the reader (in test.csv)
    for (auto varName : spectatorNames)
    {
	Float_t spectator (0.0);
	reader->AddSpectator (varName.c_str(), &spectator);
	++itVar;
    }


    TString dir    = "weights/";
    TString prefix = "TMVAClassification";
    TString weightfile = dir + prefix + TString("_") + method_name + TString(".weights.xml");
    std::cout << "weightfile name : " << weightfile.Data () << std::endl;
    reader->BookMVA( method_name, weightfile ); 
  


    // --------- for each of the input files
    for (auto inputName : inputNames)
    {
        // --- define variables
	Int_t id;
	Float_t prediction;
	Float_t weight;
	Float_t min_ANNmuon;
	Float_t mass;
	Float_t signal;

       
        // --- open input file
	TFile *input(0);
        std::stringstream infilename;
        infilename << pathToData.Data () << inputName << ".root";
	std::cout << "infilename = " << infilename.str ().c_str () << std::endl;
	input = TFile::Open (infilename.str ().c_str ());
	TTree* tree = (TTree*)input->Get("data");


        // --- prepare branches on input file
	// id field if needed
	if (contains (varsForInput, inputName, "id"))
//.........这里部分代码省略.........

TString useAutoencoder (TString method_name)
{
    TMVA::Tools::Instance();

    std::cout << "==> Start useAutoencoder" << std::endl;
    TMVA::Reader *reader = new TMVA::Reader( "!Color:!Silent" );

    Float_t signal = 0.0;
    Float_t outSignal = 0.0;
    Float_t inSignal = 0.0;

    std::vector<std::string> localVariableNames (variableNames+additionalVariableNames);
  
    std::vector<Float_t> variables (localVariableNames.size ());
    auto itVar = begin (variables);
    for (auto varName : localVariableNames)
    {
	Float_t* pVar = &(*itVar);
	reader->AddVariable(varName.c_str(), pVar);
	(*itVar) = 0.0;
	++itVar;
    }
    int idxSignal = std::distance (localVariableNames.begin (),
				   std::find (localVariableNames.begin (), localVariableNames.end (),std::string ("signal")));

  
    TString dir    = "weights/";
    TString prefix = "TMVAAutoencoder";
    TString weightfile = dir + prefix + TString("_") + method_name + TString(".weights.xml");
    TString outPrefix = "transformed";
    TString outfilename = pathToData + outPrefix + TString("_") + method_name + TString(".root");
    reader->BookMVA( method_name, weightfile );

  
    TFile* outFile = new TFile (outfilename.Data (), "RECREATE");

  
  
    std::vector<std::string> inputNames = {"training"};
    std::map<std::string,std::vector<std::string>> varsForInput;
    varsForInput["training"].emplace_back ("id");
    varsForInput["training"].emplace_back ("signal");

  
    for (auto inputName : inputNames)
    {
	std::stringstream outfilename;
	outfilename << inputName << "_transformed__" << method_name.Data () << ".root";
	std::cout << outfilename.str () << std::endl;
	/* return; */
      
	std::stringstream infilename;
	infilename << pathToData.Data () << inputName << ".root";

	TTree* outTree = new TTree("transformed","transformed");
      
	std::vector<Float_t> outVariables (localVariableNames.size ());
	itVar = begin (variables);
	auto itOutVar = begin (outVariables);
	for (auto varName : localVariableNames)
        {
	    Float_t* pOutVar = &(*itOutVar);
	    outTree->Branch (varName.c_str (), pOutVar, "F");
	    (*itOutVar) = 0.0;
	    ++itOutVar;

	    Float_t* pVar = &(*itVar);
	    std::stringstream svar;
	    svar << varName << "_in";
	    outTree->Branch (svar.str ().c_str (), pVar, "F");
	    (*itVar) = 0.0;
	    ++itVar;
        }
	Float_t signal_original = 0.0;
	outTree->Branch ("signal_original", &signal_original, "F");

	TFile *input(0);
	std::cout << "infilename = " << infilename.str ().c_str () << std::endl;
	input = TFile::Open (infilename.str ().c_str ());
	TTree* tree = (TTree*)input->Get("data");
  
	Int_t ids;

	// id field if needed
	if (std::find (varsForInput[inputName].begin (), varsForInput[inputName].end (), "id") != varsForInput[inputName].end ())
	    tree->SetBranchAddress("id", &ids);

      
	// variables for prediction
	itVar = begin (variables);
	for (auto inputName : localVariableNames)
        {
	    Float_t* pVar = &(*itVar);
	    tree->SetBranchAddress (inputName.c_str(), pVar);
	    ++itVar;
        }
 
	for (Long64_t ievt=0; ievt < tree->GetEntries(); ievt++)
        {
	    tree->GetEntry(ievt);
//.........这里部分代码省略.........

/*
void allBranches(TTree* inTree){
  
  inTree->SetBranchAddress("fullpmet",&fullpmet);
  inTree->SetBranchAddress("trkpmet",&trkpmet);
  inTree->SetBranchAddress("ratioMet",&ratioMet);
  inTree->SetBranchAddress("ptll",&ptll);
  inTree->SetBranchAddress("mth",&mth);
  inTree->SetBranchAddress("jetpt1",&jetpt1);
  inTree->SetBranchAddress("ptWW",&ptWW);
  inTree->SetBranchAddress("dphilljet",&dphilljet);
  inTree->SetBranchAddress("dphillmet",&dphillmet);
  inTree->SetBranchAddress("dphijet1met",&dphijet1met);
  inTree->SetBranchAddress("nvtx",&nvtx);

  inTree->SetBranchAddress("baseW",&baseW);        
}

void createOutput(TTree* outTree){

  outTree->Branch("fullpmet",&fullpmet);
  outTree->Branch("trkpmet",&trkpmet);
  outTree->Branch("ratioMet",&ratioMet);
  outTree->Branch("ptll",&ptll);
  outTree->Branch("mth",&mth);
  outTree->Branch("jetpt1",&jetpt1);
  outTree->Branch("ptWW",&ptWW);
  outTree->Branch("dphilljet",&dphilljet);
  outTree->Branch("dphillmet",&dphillmet);
  outTree->Branch("dphijet1met",&dphijet1met);
  outTree->Branch("nvtx",&nvtx);

  outTree->Branch("baseW",&baseW);        

}
*/
void read(TString sampleName = "WW")
{
  //calling the reader of the MVA analysis
  TMVA::Reader* reader = new TMVA::Reader("");
  
  reader->AddVariable("fullpmet",&fullpmet);
  reader->AddVariable("trkpmet",&trkpmet);
  reader->AddVariable("ratioMet",&ratioMet);
  reader->AddVariable("ptll",&ptll);
  reader->AddVariable("mth",&mth);
  reader->AddVariable("jetpt1",&jetpt1);
  reader->AddVariable("ptWW",&ptWW);
  reader->AddVariable("dphilljet",&dphilljet);
  reader->AddVariable("dphillmet",&dphillmet);
  reader->AddVariable("dphijet1met",&dphijet1met);
  reader->AddVariable("nvtx",&nvtx);
  
  //reader->BookMVA("Fisher", "weights/MVAnalysis_Fisher.weights.xml");
  reader->BookMVA("BDT", "weights/" + sampleName + "_BDT.weights.xml");
  
  /*
  //Calling WW Signal File and Tree and Creating Output Signal Trees 
  TFile *MySigFile = new TFile("../rootFiles/AllJet/OF/WW.root","READ");
  TTree* MySigTree = (TTree*)MySigFile->Get("nt");

  MySigTree->SetBranchAddress("fullpmet",&fullpmet);
  MySigTree->SetBranchAddress("trkpmet",&trkpmet);
  MySigTree->SetBranchAddress("ratioMet",&ratioMet);
  MySigTree->SetBranchAddress("ptll",&ptll);
  MySigTree->SetBranchAddress("mth",&mth);
  MySigTree->SetBranchAddress("jetpt1",&jetpt1);
  MySigTree->SetBranchAddress("ptWW",&ptWW);
  MySigTree->SetBranchAddress("dphilljet",&dphilljet);
  MySigTree->SetBranchAddress("dphillmet",&dphillmet);
  MySigTree->SetBranchAddress("dphijet1met",&dphijet1met);
  MySigTree->SetBranchAddress("nvtx",&nvtx);

  TTree *outSigTree = new TTree ("Out","Out");
  outSigTree -> SetDirectory(0);
  createOutput(outSigTree);
  */
  const int nProcesses = 2;

  enum{iWW, iDY};

  TFile *MyFile[nProcesses];
  TTree *MyTree[nProcesses];
  TTree *outTree[nProcesses];

  TString MyName[nProcesses];

  MyName[iWW]    = "WW";
  MyName[iDY]    = "DY";

  Float_t pt1;
  Float_t pt2;
  Float_t ptll;
  Float_t mll;
  Float_t mth;
  Float_t pfType1Met;
  Float_t drll;
  Float_t dphill;
  Float_t dphilljet;
  Float_t dphillmet;
//.........这里部分代码省略.........