train.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 "TMVARegGui.C"
15 
16 //#if not defined(__CINT__) || defined(__MAKECINT__)
17 #include "TMVA/Tools.h"
18 #include "TMVA/Factory.h"
19 #include "TMVA/DataLoader.h"
20 //#endif
21 
22 using namespace TMVA;
23 
24 void train()
25 {
26  TFile *file = new TFile("files/hadded_all_files.root");
27 
28  // Need to balance the beam MC / cosmic weights using livetimes.
29  // Use recent estimates for each
30  // Note, overall exposure scale doesn't affect training, but lt/POT does
31  TFile *filepot = new TFile("files/hadded_fluxswap_files.root");
32 
33  double pot = 9.47964e20;
34  double lt = 438.1584;
35 
36  TH1 *hPOT = (TH1F*)filepot->Get("pot");
37  double totPOT = hPOT->Integral();
38  cout<<"pot from fluxswap is "<<totPOT<<endl;
39  delete hPOT;
40  TH1 *hLt = (TH1F*)file->Get("livetime");
41  double totLt = hLt->Integral();
42  delete hLt;
43 
44  double nuScale = pot/totPOT;
45  double cosScale= lt/totLt;
46 
47  cout<<" scale for mc "<<nuScale;
48  cout<<" scale for cosmic "<<cosScale;
49 
50  TTree *sigTree = (TTree*)file->Get("sigTree");
51  TTree *bkgTree = (TTree*)file->Get("bkgTree");
52 
53  TFile* outputFile = TFile::Open("files/test_current_bdt.root","recreate");
54  TMVA::Factory *factory = new TMVA::Factory("MVAnalysis",outputFile,"V");
55  TMVA::DataLoader *dataloader=new TMVA::DataLoader("dataset");
56 
57  dataloader->AddVariable("nhit", 'F');
58  dataloader->AddVariable("cvnsse", 'F');
59  //dataloader->AddVariable("cvne", 'F');
60  dataloader->AddVariable("sparseassym",'F');
61  //dataloader->AddVariable("pxp", 'F');
62  //dataloader->AddVariable("pyp", 'F');
63  dataloader->AddVariable("ptp", 'F');
64  dataloader->AddVariable("disttop", 'F');
65  dataloader->AddVariable("distnottop", 'F');
66  //dataloader->AddVariable("vtxx", 'F');
67  //dataloader->AddVariable("vtxy", 'F');
68  //dataloader->AddVariable("vtxz", 'F');
69  //dataloader->AddVariable("prongl", 'F');
70  //dataloader->AddVariable("inelast", 'F');
71  //dataloader->AddVariable("hitsperpl", 'F');
72  //dataloader->AddVariable("numofshow", 'I');
73  //dataloader->AddVariable("widthofshow",'F');
74  //dataloader->AddVariable("costheta", 'F');
75 
76  // These are questionable
77  //dataloader->AddVariable("calE", 'F');
78  //dataloader->AddVariable("costhetaZ", 'F');
79  //dataloader->AddVariable("prongl", 'F');
80  //dataloader->AddVariable("remid", 'F');
81 
82  //dataloader->AddVariable("distfront", 'F');
83  //dataloader->AddVariable("distback", 'F');
84  //dataloader->AddVariable("disteast", 'F');
85  //dataloader->AddVariable("cosang", 'F');
86  //dataloader->AddVariable("vtxz", 'F');
87  // dataloader->AddVariable("vtxx", 'F');
88  // dataloader->AddVariable("vtxy", 'F');
89  // dataloader->AddVariable("vtxz", 'F');
90  //dataloader->AddVariable("efrac", 'F');
91  //dataloader->AddVariable("nplfr", 'F');
92  //dataloader->AddVariable("maxy", 'F');
93  //dataloader->AddVariable("distwest", 'F');
94  //dataloader->AddVariable("distbottom", 'F');
95  //dataloader->AddVariable("lid", 'F');
96  //dataloader->AddVariable("starty", 'F');
97  //dataloader->AddVariable("stopy", 'F');
98 
99  // Dumb oscillations, for ok 1st order prediction estimate
100  dataloader->AddSpectator("woscdumb", 'F');
101  // All necessary for determining signal / cosmic bkg
102  dataloader->AddSpectator("isGENIE", 'F');
103  dataloader->AddSpectator("pdg", 'F');
104  dataloader->AddSpectator("pdgorig", 'F');
105  dataloader->AddSpectator("isCC", 'F');
106 
107 
108  // Define a preselection here, may well influence results
109  // Also add in a truth cut to reject beam bkg from training
110  TCut presel = "cvnsse>0.5 && prongl > 0";
111  //dataloader->AddCut(presel);
112 
113  //TCut sigCut = "pdg==-12 && pdgorig==-14";
114  //TCut bkgCut = "pdg!=-12 || pdgorig!=-14"
115 
116  dataloader->AddSignalTree (sigTree);
117  dataloader->AddBackgroundTree(bkgTree);
118 
119  dataloader->AddCut(presel);
120 
121  // Apply osc weight to beam events, scale cosmics by lt
122  std::string weightExp = "woscdumb*" + std::to_string(nuScale) + "*isGENIE*" + "(pdg==12&&pdgorig==14)+" + std::to_string(cosScale) + "*!isGENIE";
123  std::cout << weightExp << std::endl;
124  dataloader->SetWeightExpression(weightExp.c_str());
125 
126  dataloader->PrepareTrainingAndTestTree("", "SplitMode=random");
127 
128  // A few methods here, to easily swap out and try
129  //factory->BookMethod(TMVA::Types::kCuts, "Cuts","H:!V:FitMethod=MC:EffSel:SampleSize=200000:VarProp=FSmart" );
130  //factory->BookMethod(TMVA::Types::kLikelihood, "Likelihood", "V");
131  //factory->BookMethod(TMVA::Types::kLD, "LD", "H:!V:VarTransform=None:CreateMVAPdfs:PDFInterpolMVAPdf=Spline2:NbinsMVAPdf=50:NsmoothMVAPdf=10");
132  // This was the method we settled on in 2016, and retrained with in 2017
133  factory->BookMethod(dataloader, TMVA::Types::kBDT, "test_current_bdt","H:V:NTrees=500:InverseBoostNegWeights::BoostType=RealAdaBoost:Shrinkage=0.1:nCuts=40:MaxDepth=5:MinNodeSize=2%:NodePurityLimit=0.9");
134  //factory->BookMethod( TMVA::Types::kMLP, "MLP", "H:!V:NeuronType=tanh:VarTransform=N:NCycles=600:HiddenLayers=N+5:TestRate=5:!UseRegulator" );
135 
136 
137  // factory->BookMethod( TMVA::Types::kBDT, "BDTG",
138  // "H:V:NTrees=4000:InverseBoostNegWeights::BoostType=Grad:Shrinkage=0.15:nCuts=40:MaxDepth=30:MinNodeSize=2%:GradBaggingFraction=0.5:NodePurityLimit=0.9");
139 
140  factory->TrainAllMethods();
141  factory->TestAllMethods();
142  factory->EvaluateAllMethods();
143 
144  outputFile->Close();
145  delete factory;
146  delete dataloader;
147 }
148 
Preselection Object.
Definition: FillPIDs.h:20
double lt
#define pot
OStream cout
Definition: OStream.cxx:6
void train()
Definition: train.C:24
TFile * file
Definition: cellShifts.C:17
Definition: tmvaglob.h:28
std::string to_string(ModuleType mt)
Definition: ModuleType.h:32
enum BeamMode string