Train_CosRej_BDTMLP.C
Go to the documentation of this file.
1 #include <cstdlib>
2 #include <iostream>
3 #include <map>
4 #include <string>
5 
6 #include "TChain.h"
7 #include "TFile.h"
8 #include "TTree.h"
9 #include "TString.h"
10 #include "TObjString.h"
11 #include "TSystem.h"
12 #include "TROOT.h"
13 
14 //#include "TRInterface.h"
15 //#include "TMVA/MethodC50.h"
16 //#include "TMVA/MethodRSNNS.h"
17 //#include "TMVA/MethodRXGB.h"
18 #include "TMVA/Types.h"
19 #include "TMVA/Tools.h"
20 #include "TMVA/Factory.h"
21 #include "TMVA/TMVAGui.h"
22 
23 using namespace TMVA;
24 
25 void Train_CosRej_BDTMLP( std::string sIdentifier, bool IsFHC, bool IsBPF ) {
26 
27 
28  std::cout << "I passed the following; sIdentifier is " << sIdentifier << "..... IsFHC is " << IsFHC << "..... IsBPF is " << IsBPF << "......." << std::endl;
29 
30  TMVA::Tools::Instance();
31 
32  // --- Am I doing RHC or not?
33  std::string sFHC_RHC = "fhc";
34  if (!IsFHC) sFHC_RHC = "rhc";
35 
36  std::string sBPF = "BPF";
37  if (!IsBPF) sBPF = "Kal";
38 
39  // --- Open my files
40  // Make some strings
41  std::string sSigTrain = "/nova/ana/users/karlwarb/PostAna18/CosRej/TrimCAFs-190215/Merged-"+sFHC_RHC+"-"+sIdentifier+"-BeamMC-Train.root";
42  std::string sSigTest = "/nova/ana/users/karlwarb/PostAna18/CosRej/TrimCAFs-190215/Merged-"+sFHC_RHC+"-"+sIdentifier+"-BeamMC-Test.root";
43 
44  std::string sBkgTrain = "/nova/ana/users/karlwarb/PostAna18/CosRej/TrimCAFs-190326/Merged-"+sFHC_RHC+"-"+sIdentifier+"-Cosmics-Train.root";
45  std::string sBkgTest = "/nova/ana/users/karlwarb/PostAna18/CosRej/TrimCAFs-190326/Merged-"+sFHC_RHC+"-"+sIdentifier+"-Cosmics-Test.root";
46  // Some output.
47  std::cout << "\n You have set IsFHC to " << IsFHC << " therefore;"
48  << "\n\t sSigTrain is " << sSigTrain
49  << "\n\t sSigTest is " << sSigTest
50  << "\n\t sBkgTrain is " << sBkgTrain
51  << "\n\t sBkgTest is " << sBkgTest
52  << std::endl;
53 
54  // --- Declare and get my Trees
55  // Open the TFiles
56  TFile *SigTrainFile = new TFile( sSigTrain.c_str() );
57  TFile *SigTestFile = new TFile( sSigTest .c_str() );
58  TFile *BkgTrainFile = new TFile( sBkgTrain.c_str() );
59  TFile *BkgTestFile = new TFile( sBkgTest .c_str() );
60  // Define and grab the trees.
61  TTree *SigTrainTree, *SigTestTree;
62  TTree *BkgTrainTree, *BkgTestTree;
63  SigTrainFile -> GetObject("TrimTree", SigTrainTree );
64  SigTestFile -> GetObject("TrimTree", SigTestTree );
65  BkgTrainFile -> GetObject("TrimTree", BkgTrainTree );
66  BkgTestFile -> GetObject("TrimTree", BkgTestTree );
67 
68  // --- Create the output file, and the Factory.
69  std::string OutName = "cosrej_tmva_output-"+sFHC_RHC+"-"+sBPF+"-"+sIdentifier+".root";
70  TFile* outputFile = TFile::Open( OutName.c_str(), "RECREATE" );
71  TMVA::Factory factory("TMVAClassification", outputFile,
72  "!V:ROC:!Correlations:!Silent:Color:DrawProgressBar:AnalysisType=Classification" );
73 
74  // --- Make the loader and add my variables.
75  TMVA::DataLoader loader("dataset");
76  // Set a bunch of variables that I want to train on.
77  loader.AddVariable("min(CosFwdCell+CosBakCell,KalTrkFwdCell+KalTrkBakCell)");
78  // --- What are my BPF Variables?
79  if (IsBPF) {
80  loader.AddVariable("BPFCosNumi");
81  loader.AddVariable("BPFLen");
82  loader.AddVariable("max(BPFStY,BPFEndY)");
83  loader.AddVariable("cos(BPFDirY)");
84  loader.AddVariable("BPFNHit/SliceNHit");
85  loader.AddVariable("BPFPtOverP");
86  }
87  // --- What are my Kalman Track Variables?
88  else {
89  loader.AddVariable("KalTrkCosNumi");
90  loader.AddVariable("KalTrkLen");
91  loader.AddVariable("max(KalTrkStY,KalTrkEndY)");
92  loader.AddVariable("cos(KalTrkDirY)");
93  loader.AddVariable("KalTrkNHit/SliceNHit");
94  loader.AddVariable("KalTrkPtOverP");
95  }
96  // --- Set a bunch of variables which just look over the process.
97  loader.AddSpectator("IsNu");
98  loader.AddSpectator("IsNuMu");
99  loader.AddSpectator("CVNMuon17");
100  loader.AddSpectator("CVNMuon18");
101  loader.AddSpectator("KalTrkRemID");
102  //loader.AddSpectator("IsCC");
103  //loader.AddSpectator("TrueNuE");
104 
105  // --- Define my Signal Cut
106  // Signal Cut
107  TCut signal = "CVNMuon17>0.1&&CVNMuon18>0.4&&KalTrkRemID>0.4&&IsNu==1&&IsNuMu==1&&IsCC==1";
108  // Background Cut
109  TCut bkgrnd = "CVNMuon17>0.1&&CVNMuon18>0.4&&KalTrkRemID>0.4&&IsNu==0";
110 
111  // --- Now want to add my Signal and Background trees.
112  // It is useful to figure out the ratio of signal / background events that we expect.
113  // These numbers can be taken out of the Cut Flow Table for the previous analysis.
114  // Unfortunately if any "Loose PID" cut is applied these will likely not be in the table.
115  // Therefore to calculate the expected number you have to calculate a ratio from
116  // [N_Cont] the number of expected events after containment (Cut Flow Table)
117  // [TrCont] the number of training events after containment (In training file)
118  // [TrPID ] the number of training events after "Loose PID" (In training file)
119  // The number of expected events is then;
120  // N_Cont * TrPID / TrCont.
121 
122  double SigTreeWeight = 1.;
123  double BkgTreeWeight = 1.;
124  //if ( IsFHC && sIdentifier == "HighGain" ) { SigTreeWeight = 1.690e-3; BkgTreeWeight = 1; }
125  /*
126  else if ( IsFHC && sIdentifier == "Period1" ) { SigTreeWeight = 1.852e-6; BkgTreeWeight = 7.143e-3; }
127  else if ( IsFHC && sIdentifier == "Period2" ) { SigTreeWeight = 1.852e-6; BkgTreeWeight = 7.143e-3; }
128  else if (!IsFHC && sIdentifier == "HighGain" ) { SigTreeWeight = 5.455e-5; BkgTreeWeight = 3.571e-3; }
129  */
130 
131  std::cout << "\n==========================================\n"
132  << "\n\t The signal and background weights are;"
133  << "\n\t Signal: " << SigTreeWeight
134  << "\n\t Backgr: " << BkgTreeWeight
135  << "\n==========================================\n"
136  << std::endl;
137 
138  // --- Add the tree, with the appropriate weights.
139  std::cerr << "Adding Sig" << std::endl;
140  loader.AddSignalTree ( SigTrainTree, SigTreeWeight, TMVA::Types::kTraining );
141  std::cerr << "Adding Sig 2" << std::endl;
142  loader.AddSignalTree ( SigTestTree , SigTreeWeight, TMVA::Types::kTesting );
143  std::cerr << "Adding Bkg" << std::endl;
144  loader.AddBackgroundTree( BkgTrainTree, BkgTreeWeight, TMVA::Types::kTraining );
145  std::cerr << "Adding Bkg 2" << std::endl;
146  loader.AddBackgroundTree( BkgTestTree , BkgTreeWeight, TMVA::Types::kTesting );
147  std::cerr << "Added Bkg" << std::endl;
148 
149  // ---- Add the trees.
150  loader.PrepareTrainingAndTestTree( signal, bkgrnd, "SplitMode=Random:NormMode=NumEvents:!V" );
151 
152  //factory.SetWeightExpression("osc"); // oscillate MC; could apply ppfx/x-sec weights here too if we had them
153 
154  // --- TMVA options: http://tmva.sourceforge.net/old_site/
155  // --- TMVA User Guide: https://root.cern.ch/download/doc/tmva/TMVAUsersGuide.pdf
156 
157  // Multi-Layer Perceptron (Neural Network)
158  factory.BookMethod(&loader, TMVA::Types::kMLP, "MLP",
159  "!H:!V:NeuronType=tanh:VarTransform=N:NCycles=100:HiddenLayers=N+5:TestRate=5!UseRegulator" );
160 
161  // The BDT which Kirk settled on.
162  factory.BookMethod(&loader, TMVA::Types::kBDT, "BDT",
163  //"!H:!V:NTrees=500:InverseBoostNegWeights::BoostType=RealAdaBoost:Shrinkage=0.1:nCuts=40:MaxDepth=5:MinNodeSize=2%:NodePurityLimit=0.9");
164  "!H:!V:NTrees=500:InverseBoostNegWeights::BoostType=uBoost:Shrinkage=0.1:nCuts=40:MaxDepth=5:MinNodeSize=2%:NodePurityLimit=0.9");
165 
166  factory.TrainAllMethods();
167  factory.TestAllMethods();
168  factory.EvaluateAllMethods();
169 
170  outputFile->Close();
171 
172  TMVA::TMVAGui( OutName.c_str() );
173 }
std::string sBPF
bool IsFHC
OStream cerr
Definition: OStream.cxx:7
f GetObject("mytree", mytree)
loader
Definition: demo0.py:10
OStream cout
Definition: OStream.cxx:6
void Train_CosRej_BDTMLP(std::string sIdentifier, bool IsFHC, bool IsBPF)
Definition: tmvaglob.h:28
enum BeamMode string