trimmubarid.C
Go to the documentation of this file.
1 ///............................................
2 /// trimvar.C
3 /// \author Connor Johnson
4 /// \brief Rewritten and cleaned up version of Bishu's code to
5 /// make trimfiles for muonid training.
6 /// This one should be much more streamlined in general
7 ///...........................................
9 #include "CAFAna/Core/Progress.h"
12 
14 
15 #include <algorithm>
16 #include <iostream>
17 #include <string>
18 
19 //From Root
20 #include "TStyle.h"
21 #include "TFile.h"
22 #include "TTree.h"
23 #include "TVector3.h"
24 #include "TH1.h"
25 
27 
28 void trimmubarid(const std::string outputFilename = "fMubarIDTrain.root", const std::string setname = "fhctrain")
29 {
30 
31  // Nominal RHC MC defnames
32  std::map<string, string> training_datasets =
33  {
34  {"fhctrain", "muonid_prod5_fhc_training_sample"},
35  {"rhctrain", "muonid_prod5_rhc_training_sample"}
36  };
37  std::map<string, string> testing_datasets =
38  {
39  {"fhctest", "muonid_prod5_fhc_testing_sample"},
40  {"rhctest", "muonid_prod5_rhc_testing_sample"}
41  };
42 
43  std::string samQuery = "";
44  if (training_datasets.find(setname) != training_datasets.end()){
45  cout << "Making Training Trim file from dataset: " << training_datasets[setname] << endl;
46  samQuery = training_datasets[setname];
47  }
48  else if (testing_datasets.find(setname) != testing_datasets.end()){
49  cout << "Making Testing Trim file from dataset: " << testing_datasets[setname] << endl;
50  samQuery = testing_datasets[setname];
51  }
52  else{
53  cerr << "Set not found in training or testing definitions. Exiting." << endl;
54  return;
55  }
56 
57  TFile outputFile(outputFilename.c_str(), "RECREATE");
58 
59  TTree *sigTree = new TTree("sigTree", "Signal Events");
60  TTree *bkgTree = new TTree("bkgTree", "Background Events");
61 
62  //...............................................................................
63 
64  // Training variables for MuonID
65  float fScatLL = -999.0;
66  float fDedxLL = -999.0;
67  float fAvededxlast10cm = -999.0;
68  float fAvededxlast40cm = -999.0;
69 
70  //-------------------------------------------------------------------
71 
72  // Set the branch addresses for both signal / background trees
73  sigTree->Branch("DedxLL", &fDedxLL, "DedxLL/F");
74  sigTree->Branch("ScatLL", &fScatLL, "ScatLL/F");
75  sigTree->Branch("Avededxlast10cm", &fAvededxlast10cm, "Avededxlast10cm/F");
76  sigTree->Branch("Avededxlast40cm", &fAvededxlast40cm, "Avededxlast40cm/F");
77 
78  bkgTree->Branch("DedxLL", &fDedxLL, "DedxLL/F");
79  bkgTree->Branch("ScatLL", &fScatLL, "ScatLL/F");
80  bkgTree->Branch("Avededxlast10cm", &fAvededxlast10cm, "Avededxlast10cm/F");
81  bkgTree->Branch("Avededxlast40cm", &fAvededxlast40cm, "Avededxlast40cm/F");
82 
83  //---------------------------------------------------------------------
84 
85  // Begin looping over the selected files from the SAM Database.
86 
87  int stride = -1;
88  int offset = -1;
89  int limit = -1;
90 
91  if(getenv("CAFANA_STRIDE")){
92  stride = atoi(getenv("CAFANA_STRIDE"));
93  if(stride > 1 && getenv("CAFANA_OFFSET")){
94  offset = atoi(getenv("CAFANA_OFFSET"));
95  }
96  }
97  if(getenv("CAFANA_LIMIT"))
98  limit = atoi(getenv("CAFANA_LIMIT"));
99 
100  ana::SAMQuerySource* filesrc = new ana::SAMQuerySource(samQuery, stride, offset, limit);
101 
102  int Nfiles = filesrc->NFiles();
103 
104  ana::Progress* prog = new ana::Progress("Filling signal and background trees from " + to_string(Nfiles) + " files.");
105 
106  unsigned int fileIdx = 0;
107  TFile* file1;
108 
109  // Swapped from while loop to do-while, so we don't lose the first file
110  while((file1 = (TFile*)filesrc->GetNextFile())) {
111  prog->SetProgress((fileIdx+1.)/Nfiles);
112 
113  // Get the recTree object
114  TTree *recTree = (TTree*)file1->Get("recTree");
115 
116  // Define the sr branch
117  caf::StandardRecord* sr = 0;
118  recTree->SetBranchAddress("rec", &sr);
119 
120  //--------------------------------------------------------
121 
122  //Begin filling the branches with variables and cuts applied.
123  for (int event = 0; event < recTree->GetEntries(); ++event) {
124  recTree->GetEntry(event);
125 
126  // Basic quality cuts:
127  // 1. Valid elastic vertex
128  // 2. At least one true neutrino
129  // 3. At least one kalman track
130  // 4. At least 20 hits
131  // 5. At leats 4 continuous planes
132  if(!sr->vtx.elastic.IsValid || sr->mc.nnu < 1 || sr->trk.kalman.ntracks < 1 || sr->slc.nhit < 20 || sr->slc.ncontplanes < 4)
133  continue;
134 
135  // Containment Cut
136  bool iscontprong = true;
137  for(const caf::SRFuzzyKProng& prong : sr->vtx.elastic.fuzzyk.png) {
138  TVector3 start = prong.shwlid.start;
139  TVector3 stop = prong.shwlid.stop;
140  if( std::min( start.X(), stop.X() ) < -180.0 ||
141  std::max( start.X(), stop.X() ) > 180.0 ||
142  std::min( start.Y(), stop.Y() ) < -180.0 ||
143  std::max( start.Y(), stop.Y() ) > 180.0 ||
144  std::min( start.Z(), stop.Z() ) < 20.0 ||
145  std::max( start.Z(), stop.Z() ) > 1525.0 )
146  {
147  iscontprong = false;
148  break;
149  }
150  }
151  if (!iscontprong)
152  continue;
153 
154  // loop over track
155  for (const caf::SRKalmanTrack& track : sr->trk.kalman.tracks){
156 
157  // Check the remid calculation was successful (necessary for de/dx LL)
158  if(track.rempid < 0 ) continue;
159 
160  fDedxLL = track.dedxllh;
161  fScatLL = track.scatllh;
162 
163  if (track.avedEdxlast10cm < 0.) continue;
164  if (track.avedEdxlast40cm < 0.) continue;
165 
166  fAvededxlast10cm = min(track.avedEdxlast10cm, float(30.));
167  fAvededxlast40cm = min(track.avedEdxlast40cm, float(30.));
168 
169  // If true muon or mubar, fill the signal tree. Otherwise fill in into the bkgd tree.
170  if(abs(track.truth.pdg) == 13)
171  sigTree->Fill();
172  else
173  bkgTree->Fill();
174 
175  } // end of trk loop
176  } // end of event loop
177  ++fileIdx;
178  } // end of file loop
179 
180  prog->Done();
181 
182  outputFile.Write();
183  outputFile.Close();
184 
185  cout << "Finished making Trim for " << samQuery << endl;
186 
187 } // end of void
188 
189 
SRShowerLID shwlid
Shower information.
Definition: SRFuzzyKProng.h:19
T max(const caf::Proxy< T > &a, T b)
void trimmubarid(const std::string outputFilename="fMubarIDTrain.root", const std::string setname="fhctrain")
Rewritten and cleaned up version of Bishu&#39;s code to make trimfiles for muonid training. This one should be much more streamlined in general ...........................................
Definition: trimmubarid.C:28
size_t ntracks
Definition: SRKalman.h:23
std::vector< SRFuzzyKProng > png
Vector of 3D prong objects.
Definition: SRFuzzyK.h:19
Definition: event.h:19
OStream cerr
Definition: OStream.cxx:7
void abs(TH1 *hist)
float avedEdxlast40cm
Average dE/dx in the last 40 cm approximately.
Definition: SRTrack.h:75
SRKalman kalman
Tracks produced by KalmanTrack.
Definition: SRTrackBranch.h:24
SRParticleTruth truth
Truth information for the track.
Definition: SRTrack.h:48
int pdg
PDG Code of the best matched truth particle.
int NFiles() const override
May return -1 indicating the number of files is not known.
std::string getenv(std::string const &name)
caf::StandardRecord * sr
File source based on a SAM query or dataset (definition)
short nnu
Number of neutrinos in nu vector (0 or 1)
Definition: SRTruthBranch.h:37
SRVector3D start
Shower start point in detector coordinates. [cm].
Definition: SRShower.h:32
unsigned int nhit
number of hits
Definition: SRSlice.h:22
OStream cout
Definition: OStream.cxx:6
static float min(const float a, const float b, const float c)
Definition: absgeo.cxx:45
The StandardRecord is the primary top-level object in the Common Analysis File trees.
void SetProgress(double frac)
Update the progress fraction between zero and one.
Definition: Progress.cxx:41
SRElastic elastic
Single vertex found by Elastic Arms.
virtual TFile * GetNextFile() override
Returns the next file in sequence, ready for reading.
A simple ascii-art progress bar.
Definition: Progress.h:9
SRSlice slc
Slice branch: nhit, extents, time, etc.
SRTruthBranch mc
Truth branch for MC: energy, flavor, etc.
SRFuzzyK fuzzyk
Primary 3D prong object.
Definition: SRElastic.h:44
SRVector3D stop
shower stop point
Definition: SRShower.h:41
std::string to_string(ModuleType mt)
Definition: ModuleType.h:32
float avedEdxlast10cm
Average dE/dx in the last 10 cm approximately.
Definition: SRTrack.h:72
T min(const caf::Proxy< T > &a, T b)
unsigned int ncontplanes
number of continuous planes
Definition: SRSlice.h:25
SRTrackBranch trk
Track branch: nhit, len, etc.
void Done()
Call this when action is completed.
Definition: Progress.cxx:90
std::vector< SRKalmanTrack > tracks
3D Tracks produced by KalmanTrack
Definition: SRKalman.h:16
SRVertexBranch vtx
Vertex branch: location, time, etc.
enum BeamMode string