NumubarCCIncCuts.cxx
Go to the documentation of this file.
2 
3 namespace ana { namespace xsec { namespace numubarcc
4 {
5 
6  const bool testo = true;
7 
8  ////////////////////////////////////
9  /// Main analysis truth cuts
10  ////////////////////////////////////
11 
12  bool VtxInBounds(const caf::SRVector3DProxy * vec, const TVector3 * vmin, const TVector3 * vmax){
13  return (vec->X() >= vmin->X() && vec->Y() >= vmin->Y() && vec->Z() >= vmin->Z() &&
14  vec->X() <= vmax->X() && vec->Y() <= vmax->Y() && vec->Z() <= vmax->Z());
15  }
16 
18  return (nu->iscc && nu->pdg == -14);
19  });
21 
22  const NuTruthCut kIsNumuCC_NT([](const caf::SRNeutrinoProxy* nu)
23  {
24  return (nu->iscc && nu->pdg == +14);
25  });
27 
29  const Cut kIsNumuNumubarCC = CutFromNuTruthCut(kIsNumuNumubarCC_NT);
30 
31  ////////////////////////////////////
32  /// Signal cuts
33  ////////////////////////////////////
34  const NuTruthCut kIsNueCC_NT([](const caf::SRNeutrinoProxy* nu){
35  return (nu->iscc && (nu->pdg == 12));
36  });
38  return (nu->iscc && (nu->pdg == -12));
39  });
41  return (nu->iscc && (std::abs(nu->pdg) == 12));
42  });
43  const NuTruthCut kIsCC_NT([](const caf::SRNeutrinoProxy* nu){
44  return (nu->iscc);
45  });
47 
48  const NuTruthCut kTrueMuon_NT([](const caf::SRNeutrinoProxy* nu){
49  for(int iprim = 0; iprim < (int) nu->prim.size(); ++iprim)
50  if(std::abs(nu->prim[iprim].pdg) == 13)
51  return true;
52  return false;
53  });
54 
59  const Cut kIsNC = CutFromNuTruthCut(kIsNC_NT);
60 
62 
63  const Cut kTrueKalmanMuon([](const caf::SRProxy* sr)
64  {
65  // for(int itrack = 0; itrack < (int) sr->trk.kalman.ntracks; ++itrack)
66  // if (std::abs(sr->trk.kalman.tracks[itrack].truth.pdg) == 13)
67  // return true;
68 
69  // return false;
70 
71  int ibesttrk = kBestMuonIDIndex(sr);
72  if(ibesttrk < 0 || ibesttrk >= int(sr->trk.kalman.ntracks))
73  return false;
74 
75  return std::abs(sr->trk.kalman.tracks[ibesttrk].truth.pdg) == 13;
76  });
77 
78  ////////////////////////////////////
79  /// Geometry
80  ////////////////////////////////////
81  const TVector3 * detector_vtx_min = new TVector3(-191, -187, 0);
82  const TVector3 * detector_vtx_max = new TVector3( 192, 194, 1270);
83 
84  const TVector3 * loose_vtx_min = new TVector3(-160, -160, 25);
85  const TVector3 * loose_vtx_max = new TVector3(160, 160, 1150);
86 
87  const TVector3 * numucc_fiducial_min = new TVector3(-130, -130, 100);
88  const TVector3 * numucc_fiducial_max = new TVector3( 140, 140, 1000);
89 
90  // All prongs and shower containment
91  const TVector3 * containLow = new TVector3(-180, -180, 20);
92  const TVector3 * containHigh = new TVector3(180, 180, 1525);
93 
95  {
96  return VtxInBounds(&nu->vtx, detector_vtx_min, detector_vtx_max);
97  });
99 
101  {
102  return VtxInBounds(&nu->vtx, loose_vtx_min, loose_vtx_max);
103  });
106  const Cut kTrueVtxCut = ana::CutFromNuTruthCut(kTrueVtxCut_NT);
107 
109  {
110  return VtxInBounds(&nu->vtx, numucc_fiducial_min, numucc_fiducial_max);
111  });
113 
114  ////////////////////////////////////
115  /// Preselection
116  ////////////////////////////////////
117  const Cut kQualityCut([](const caf::SRProxy* sr)
118  {
119  return (sr->trk.kalman.ntracks > 0 && sr->slc.nhit > 20 && sr->slc.ncontplanes > 4);
120  });
121  const Cut kContainmentCut([](const caf::SRProxy* sr)
122  {
123  if(!sr->vtx.elastic.IsValid) return false;
124  int ibesttrk = kBestMuonIDIndex(sr);
125  if(ibesttrk < 0 || ibesttrk >= int(sr->trk.kalman.ntracks))
126  return false;
127 
128  for(const caf::SRFuzzyKProngProxy& prong : sr->vtx.elastic.fuzzyk.png)
129  if(!VtxInBounds(&prong.shwlid.start, containLow, containHigh) || !VtxInBounds(&prong.shwlid.stop, containLow, containHigh))
130  return false;
131 
132  // Only most muon-like track present in muon catcher
133  if(sr->trk.kalman.ntracks < 1)
134  return false;
135  const unsigned short muon_catcher_edge = 1275;
136  for(unsigned int i = 0; i < sr->trk.kalman.ntracks; ++i)
137  {
138  if(int(i) == ibesttrk)
139  continue;
140  if( sr->trk.kalman.tracks[i].start.Z() > muon_catcher_edge ||
141  sr->trk.kalman.tracks[i].stop.Z() > muon_catcher_edge )
142  return false;
143  }
144 
145  // Muon catcher has less modules, require muon
146  const caf::SRKalmanTrackProxy& besttrack = sr->trk.kalman.tracks[ibesttrk];
147  return ((besttrack.stop.Z() < muon_catcher_edge || besttrack.trkyposattrans < 55 ) // air gap
148  && besttrack.trkfwdcellnd > 5 && besttrack.trkbakcellnd > 10);
149  });
150 
152 
153  ////////////////////////////////////
154  /// Main analysis reconstruction cuts
155  ////////////////////////////////////
156  // const Cut kTrueMubarTrk([](const caf::SRProxy* sr)
157  // {
158  // unsigned int ibesttrk = 0;
159  // if (sr->trk.kalman.ntracks < 1)
160  // return false;
161  // return (sr->trk.kalman.tracks[ibesttrk].truth.pdg == -13);
162  // });
163 
164  // const Cut kTrueMuonTrk([](const caf::SRProxy* sr)
165  // {
166  // unsigned int ibesttrk = 0;
167  // if (sr->trk.kalman.ntracks < 1)
168  // return false;
169  // return (sr->trk.kalman.tracks[ibesttrk].truth.pdg == +13);
170  // });
171 
172  // const Cut kTrueMuOrMubarTrk([](const caf::SRProxy* sr)
173  // {
174  // unsigned int ibesttrk = 0;
175  // if (sr->trk.kalman.ntracks < 1)
176  // return false;
177  // return (std::abs(sr->trk.kalman.tracks[ibesttrk].truth.pdg) == 13);
178  // });
179 
180  const Cut kRecoVtxDetectorCut([](const caf::SRProxy* sr)
181  {
182  if (!sr->vtx.elastic.IsValid)
183  return false;
184  return VtxInBounds(&sr->vtx.elastic.vtx, detector_vtx_min, detector_vtx_max);
185  });
186 
187  const Cut kRecoVtxLooseCut([](const caf::SRProxy* sr)
188  {
189  if (!sr->vtx.elastic.IsValid)
190  return false;
191  return VtxInBounds(&sr->vtx.elastic.vtx, loose_vtx_min, loose_vtx_max);
192  });
193 
194 
195  const Cut kRecoVtxNumuFiducialCut([](const caf::SRProxy* sr)
196  {
197  if (!sr->vtx.elastic.IsValid)
198  return false;
199  return VtxInBounds(&sr->vtx.elastic.vtx, numucc_fiducial_min, numucc_fiducial_max);
200  });
201 
203 
204  const float kMuonIDCutVal = 0.0;
205  const Cut kMuonIDCut([](const caf::SRProxy* sr)
206  {
207  return kMuonID(sr) >= kMuonIDCutVal;
208  });
209 
210  const Cut kActive([](const caf::SRProxy* sr)
211  {
212  int ibesttrk = kBestMuonIDIndex(sr);
213  if(sr->trk.kalman.ntracks < 1 || ibesttrk < 0 || ibesttrk >= int(sr->trk.kalman.ntracks))
214  return false;
215 
216  const caf::SRKalmanTrackProxy& bestmuon = sr->trk.kalman.tracks[ibesttrk];
217 
218  return (bestmuon.leninact > 0 &&
219  bestmuon.lenincat < 0);
220  });
221 
222  const Cut kActiveAndCatcher([](const caf::SRProxy* sr)
223  {
224  int ibesttrk = kBestMuonIDIndex(sr);
225  if(sr->trk.kalman.ntracks < 1 || ibesttrk < 0 || ibesttrk >= int(sr->trk.kalman.ntracks))
226  return false;
227 
228  const caf::SRKalmanTrackProxy& bestmuon = sr->trk.kalman.tracks[ibesttrk];
229 
230  return (bestmuon.leninact > 0 &&
231  bestmuon.lenincat > 0);
232  });
233 
234 } } } // End ana::xsec::numubar namespace
const NuTruthCut kIsNueCC_NT([](const caf::SRNeutrinoProxy *nu){return(nu->iscc &&(nu->pdg==12));})
Signal cuts.
const Cut kRecoVtxNumuFiducialCut([](const caf::SRProxy *sr){if(!sr->vtx.elastic.IsValid) return false;return VtxInBounds(&sr->vtx.elastic.vtx, numucc_fiducial_min, numucc_fiducial_max);})
caf::Proxy< double > leninact
Definition: SRProxy.h:1677
caf::Proxy< caf::SRFuzzyK > fuzzyk
Definition: SRProxy.h:2059
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
const NuTruthCut kTrueVtxDetectorCut_NT([](const caf::SRNeutrinoProxy *nu){return VtxInBounds(&nu->vtx, detector_vtx_min, detector_vtx_max);})
const NuTruthCut kTrueMuon_NT([](const caf::SRNeutrinoProxy *nu){for(int iprim=0;iprim< (int) nu->prim.size();++iprim) if(std::abs(nu->prim[iprim].pdg)==13) return true;return false;})
caf::Proxy< caf::SRVector3D > stop
Definition: SRProxy.h:1539
const TVector3 * detector_vtx_max
caf::Proxy< size_t > ntracks
Definition: SRProxy.h:1778
Proxy for caf::SRNeutrino.
Definition: SRProxy.h:510
const Cut kActive([](const caf::SRProxy *sr){int ibesttrk=kBestMuonIDIndex(sr);if(sr->trk.kalman.ntracks< 1||ibesttrk< 0||ibesttrk >=int(sr->trk.kalman.ntracks)) return false;const caf::SRKalmanTrackProxy &bestmuon=sr->trk.kalman.tracks[ibesttrk];return(bestmuon.leninact > 0 && bestmuon.lenincat< 0);})
Proxy for caf::StandardRecord.
Definition: SRProxy.h:2126
caf::Proxy< unsigned int > ncontplanes
Definition: SRProxy.h:1314
const TVector3 * containHigh
const TVector3 * numucc_fiducial_max
const NuTruthCut kIsNuebarCC_NT([](const caf::SRNeutrinoProxy *nu){return(nu->iscc &&(nu->pdg==-12));})
Proxy for caf::SRFuzzyKProng.
Definition: SRProxy.h:2003
const Cut kActiveAndCatcher([](const caf::SRProxy *sr){int ibesttrk=kBestMuonIDIndex(sr);if(sr->trk.kalman.ntracks< 1||ibesttrk< 0||ibesttrk >=int(sr->trk.kalman.ntracks)) return false;const caf::SRKalmanTrackProxy &bestmuon=sr->trk.kalman.tracks[ibesttrk];return(bestmuon.leninact > 0 && bestmuon.lenincat > 0);})
float abs(float number)
Definition: d0nt_math.hpp:39
_Cut< caf::SRNeutrinoProxy > NuTruthCut
Cut designed to be used over the nuTree, ie all neutrinos, not just those that got slices...
Definition: Cut.h:104
caf::Proxy< caf::SRElastic > elastic
Definition: SRProxy.h:2118
caf::Proxy< caf::SRVector3D > vtx
Definition: SRProxy.h:573
Track finder for cosmic rays.
caf::Proxy< std::vector< caf::SRFuzzyKProng > > png
Definition: SRProxy.h:2043
const TVector3 * loose_vtx_max
const Cut kQualityCut([](const caf::SRProxy *sr){return(sr->trk.kalman.ntracks > 0 &&sr->slc.nhit > 20 &&sr->slc.ncontplanes > 4);})
Preselection.
caf::Proxy< unsigned int > nhit
Definition: SRProxy.h:1315
const NuTruthCut kIsNueorbarCC_NT([](const caf::SRNeutrinoProxy *nu){return(nu->iscc &&(std::abs(nu->pdg)==12));})
Proxy for caf::SRVector3D.
Definition: SRProxy.h:78
caf::Proxy< double > lenincat
Definition: SRProxy.h:1678
caf::Proxy< caf::SRTrackBranch > trk
Definition: SRProxy.h:2145
const Cut kContainmentCut([](const caf::SRProxy *sr){if(!sr->vtx.elastic.IsValid) return false;int ibesttrk=kBestMuonIDIndex(sr);if(ibesttrk< 0||ibesttrk >=int(sr->trk.kalman.ntracks)) return false;for(const caf::SRFuzzyKProngProxy &prong:sr->vtx.elastic.fuzzyk.png) if(!VtxInBounds(&prong.shwlid.start, containLow, containHigh)||!VtxInBounds(&prong.shwlid.stop, containLow, containHigh)) return false;if(sr->trk.kalman.ntracks< 1) return false;const unsigned short muon_catcher_edge=1275;for(unsigned int i=0;i< sr->trk.kalman.ntracks;++i){if(int(i)==ibesttrk) continue;if(sr->trk.kalman.tracks[i].start.Z() > muon_catcher_edge|| sr->trk.kalman.tracks[i].stop.Z() > muon_catcher_edge) return false;}const caf::SRKalmanTrackProxy &besttrack=sr->trk.kalman.tracks[ibesttrk];return((besttrack.stop.Z()< muon_catcher_edge||besttrack.trkyposattrans< 55) &&besttrack.trkfwdcellnd > 5 &&besttrack.trkbakcellnd > 10);})
caf::Proxy< bool > iscc
Definition: SRProxy.h:538
caf::StandardRecord * sr
caf::Proxy< caf::SRVector3D > start
Definition: SRProxy.h:1538
const NuTruthCut kIsNumuNumubarCC_NT
const NuTruthCut kIsNC_NT
const NuTruthCut kIsNumuCC_NT([](const caf::SRNeutrinoProxy *nu){return(nu->iscc &&nu->pdg==+14);})
caf::Proxy< caf::SRShowerLID > shwlid
Definition: SRProxy.h:2023
const TVector3 * numucc_fiducial_min
Double_t xsec[nknots]
Definition: testXsec.C:47
const Cut kTrueKalmanMuon([](const caf::SRProxy *sr){ int ibesttrk=kBestMuonIDIndex(sr);if(ibesttrk< 0||ibesttrk >=int(sr->trk.kalman.ntracks)) return false;return std::abs(sr->trk.kalman.tracks[ibesttrk].truth.pdg)==13;})
const Cut kRecoVtxDetectorCut([](const caf::SRProxy *sr){if(!sr->vtx.elastic.IsValid) return false;return VtxInBounds(&sr->vtx.elastic.vtx, detector_vtx_min, detector_vtx_max);})
Main analysis reconstruction cuts.
caf::Proxy< caf::SRKalman > kalman
Definition: SRProxy.h:1797
const NuTruthCut kTrueVtxCut_NT
caf::Proxy< bool > IsValid
Definition: SRProxy.h:2058
const NuTruthCut kIsNumubarCC_NT([](const caf::SRNeutrinoProxy *nu){return(nu->iscc &&nu->pdg==-14);})
caf::Proxy< caf::SRSlice > slc
Definition: SRProxy.h:2142
const NuTruthCut kTrueVtxNumuFiducialCut_NT([](const caf::SRNeutrinoProxy *nu){return VtxInBounds(&nu->vtx, numucc_fiducial_min, numucc_fiducial_max);})
caf::Proxy< caf::SRVector3D > vtx
Definition: SRProxy.h:2073
enum BeamMode numubarcc
caf::Proxy< short int > pdg
Definition: SRProxy.h:552
const TVector3 * detector_vtx_min
Geometry.
caf::Proxy< std::vector< caf::SRKalmanTrack > > tracks
Definition: SRProxy.h:1780
const NuTruthCut kTrueVtxLooseCut_NT([](const caf::SRNeutrinoProxy *nu){return VtxInBounds(&nu->vtx, loose_vtx_min, loose_vtx_max);})
Template for Cut and SpillCut.
Definition: Cut.h:15
caf::Proxy< caf::SRVertexBranch > vtx
Definition: SRProxy.h:2146
const TVector3 * loose_vtx_min
const Cut kRecoVtxLooseCut([](const caf::SRProxy *sr){if(!sr->vtx.elastic.IsValid) return false;return VtxInBounds(&sr->vtx.elastic.vtx, loose_vtx_min, loose_vtx_max);})
Proxy for caf::SRKalmanTrack.
Definition: SRProxy.h:1735
const NuTruthCut kIsCC_NT([](const caf::SRNeutrinoProxy *nu){return(nu->iscc);})
const Cut kMuonIDCut([](const caf::SRProxy *sr){return kMuonID(sr) >=kMuonIDCutVal;})
Cut CutFromNuTruthCut(const NuTruthCut &stc)
Definition: Cut.cxx:7
caf::Proxy< std::vector< caf::SRTrueParticle > > prim
Definition: SRProxy.h:555
bool VtxInBounds(const caf::SRVector3DProxy *vec, const TVector3 *vmin, const TVector3 *vmax)
Main analysis truth cuts.
const TVector3 * containLow