15 #include "NovaDAQConventions/DAQConventions.h" 27 #include "Utilities/AssociationUtil.h" 28 #include "Utilities/func/MathUtil.h" 33 #include "TMVA/Reader.h" 69 produces< std::vector<numusand::NumuSandObj> >();
70 produces< art::Assns<numusand::NumuSandObj, rb::Cluster> >();
92 if(slicevec->empty()) {
106 std::unique_ptr< std::vector<numusand::NumuSandObj> > outputObjects(
new std::vector<numusand::NumuSandObj>);
109 for(
size_t i = 0;
i < slicevec->size(); ++
i){
130 float xmin2 = 9999, ymin2 = 9999, zmin2 = 9999, xmint = 9999, ymint = 9999, zmint = 9999;
131 float xmax2 = -9999, ymax2 = -9999, zmax2 = -9999, xmaxt = -9999, ymaxt = -9999, zmaxt = -9999;
132 int ntx = 0, nty = 0;
138 if(fmElastic.isValid()){
139 std::vector< art::Ptr<rb::Vertex> > elastics = fmElastic.at(
i);
140 if (elastics.size() == 1){
142 TVector3 vtx = elastics[0]->GetXYZ();
147 for (
unsigned int hIdx = 0; hIdx < slice->
NCell(); hIdx++){
158 vtxE20 += rhit.
GeV();
159 vtxE40 += rhit.
GeV();
160 vtxE60 += rhit.
GeV();
163 vtxE40 += rhit.
GeV();
164 vtxE60 += rhit.
GeV();
167 vtxE60 += rhit.
GeV();
177 for(
unsigned int j = 0;
j < slice->
NCell(); ++
j) {
181 const double x = xyz[0];
182 const double y = xyz[1];
183 const double z = xyz[2];
186 if(x < xmin) { xmin2 =
xmin; xmin =
x;}
187 else if(x < xmin2) xmin2 =
x;
188 if(z < zmin) { zmin2 = zmin; zmin =
z;}
189 else if(z < zmin2) zmin2 =
z;
191 else if(x > xmax2) xmax2 =
x;
193 else if (z > zmax2) zmax2 =
zmax;
197 if(x < xmint) xmint =
x;
198 if(z < zmint) zmint =
z;
200 if(x > xmaxt) xmaxt =
x;
201 if(z > zmaxt) zmaxt =
z;
207 else if (y < ymin2) ymin2 =
y;
208 if(z < zmin) { zmin2 = zmin; zmin =
z;}
209 else if (z < zmin2) zmin2 =
z;
212 else if(y > ymax2) ymax2 =
y;
214 else if(z > zmax2) zmax2 =
z;
218 if(y < ymint) ymint =
y;
219 if(z < zmint) zmint =
z;
221 if(y > ymaxt) ymaxt =
y;
222 if(z > zmaxt) zmaxt =
z;
234 if(ntx > 0 && nty > 0) {
253 if(!trackAssnList.isValid()) {
255 outputObjects->push_back(sandbox);
260 const std::vector< art::Ptr<rb::Track> > sliceTracks = trackAssnList.at(
i);
263 if(!remidVec.isValid()) {
265 outputObjects->push_back(sandbox);
275 double bestpimuE = -999.;
279 double bestval = -999., val2 = -999;
285 float avededxtrk1 = 0.0;
286 float avededxtrk2 = 0.0;
287 float avededxtrk1Last4Cells = 0.0;
288 float avededxtrk2Last4Cells = 0.0;
289 float avededxtrk1Last6Cells = 0.0;
290 float avededxtrk2Last6Cells = 0.0;
291 float avededxtrk1Last8Cells = 0.0;
292 float avededxtrk2Last8Cells = 0.0;
293 float scattangletrk1 = 0.0;
294 float scattangletrk2 = 0.0;
298 if((sliceTracks.size() == 2)){
300 for(
size_t j = 0;
j < sliceTracks.size(); ++
j){
304 if(!track->
Is3D())
continue;
306 if(trackRemid->
Value() > bestval) {
309 bestval = trackRemid->
Value();
312 else if(trackRemid->
Value() > val2) {
313 val2 = trackRemid->
Value();
320 if((!remidVec.isValid()) ||
321 (bestval < 0.0) || (val2 < 0.0)){
399 float RecoMuonE = -5.0;
400 float Reco2ndTrkE = -5.0;
403 std::vector<art::Ptr<numue::NumuE> > numuEnergy = eAss.at(
i);
404 if(numuEnergy.size() > 0){
406 trkCCE = numuEnergy[0]->TrkCCE();
410 if(trackEnergy.isValid()){
411 RecoMuonE = trackEnergy.at(bestIdx)->E();
416 float tmpHadE = trkCCE - RecoMuonE;
418 float offTrkEnFrac = (tmpHadE - Reco2ndTrkE)/trkCCE;
434 for(
size_t j = 0;
j < sliceTracks.size(); ++
j){
444 if(trackRemid->
Value() > bestval) {
447 bestval = trackRemid->
Value();
450 else if(trackRemid->
Value() > val2) {
451 val2 = trackRemid->
Value();
455 if(trackRemid->
Value() > 0.6) nmutrks++;
460 const std::vector< cheat::TrackIDE > trids = bt->
HitsToTrackIDE(trkhits);
463 if( !parts.empty() && !trids.empty() ) {
464 if(
abs(parts[0]->
PdgCode())==13 && parts[0]->Process()==
"Decay") {
465 double tempE = trids[0].energyFrac * track->
TotalGeV();
466 if(tempE > bestpimuE) {
474 if(parts.size() > 1 && trids.size() > 1) {
475 if(
abs(parts[1]->
PdgCode())==13 && parts[1]->Process()==
"Decay") {
476 double tempE = trids[1].energyFrac * track->
TotalGeV();
477 if(tempE > bestpimuE) {
498 for (
unsigned int sliceHit = 0; sliceHit < allslchits.
size(); ++sliceHit){
502 for (
unsigned int trackHit = 0; trackHit < mutrkhits.
size(); ++trackHit){
503 match = (allslchits[sliceHit] == mutrkhits[trackHit]);
507 if (!match) vertexCluster.
Add(allslchits[sliceHit]);
511 nhadhits = vertexCluster.
NCell();
517 unsigned int minCellsFromEdge = 99999999;
522 for(
unsigned int hitIdx = 0; hitIdx < vertexCluster.
NCell(); ++hitIdx){
524 const int planeNum = chit->
Plane();
527 if(cellsFromEdge < minCellsFromEdge){
528 minCellsFromEdge = cellsFromEdge;
566 if(pimuIdx == -5) pimuIdx = 0;
569 const std::vector< cheat::TrackIDE > trids = bt->
HitsToTrackIDE(slchits);
570 if(!parts.empty() && !trids.empty()) {
571 for(
unsigned int i = 0;
i < parts.size();
i++) {
572 if(trids.size() >
i) {
573 if(parts[
i]->
PdgCode()==2212 && trids[
i].energyFrac * slice->
TotalGeV() > 0.05) nprotons++;
582 outputObjects->push_back(sandbox);
586 evt.
put(std::move(outputObjects));
587 evt.
put(std::move(assoc));
void SetNHadHits(int nhadhits)
void SetXMaxt(float dxmaxt)
void SetPDGSecondBest(int pdgsecondbest)
back track the reconstruction to the simulation
static bool CreateAssn(art::EDProducer const &prod, art::Event &evt, std::vector< T > &a, art::Ptr< U > b, art::Assns< T, U > &assn, size_t indx=UINT_MAX, std::string const &instance=std::string())
Create a 1 to 1 association between a new product and one already in the event.
unsigned int NCell(geo::View_t view) const
Number of cells in view view.
std::map< std::string, double > xmax
void beginRun(art::Run &run)
void SetZMaxt(float dzmaxt)
void SetXMint(float dxmint)
void SetAvedEdxTrack1Last4Cells(float avededxtrk1Last4Cells)
unsigned short Plane() const
void SetYMax2(float dymax2)
void SetAvedEdxTrack2Last4Cells(float avededxtrk2Last4Cells)
void SetOffTrkFra(float offtrk)
proxy for off track energy
void getAveTrackdEdx(const art::Ptr< rb::Track > track, float &avededx)
void getAveTrackdEdxLast4Cells(const art::Ptr< rb::Track > track, float &avededxlast4)
Vertical planes which measure X.
FillSandbox(fhicl::ParameterSet const &pset)
unsigned int Ncells() const
Number of cells in this plane.
A collection of associated CellHits.
float getScatt(const art::Ptr< rb::Track > track)
void SetAvedEdxTrack2Last6Cells(float avededxtrk2Last6Cells)
void SetScattAngleTrack1(float scattangletrk1)
void SetAvedEdxTrack2Last8Cells(float avededxtrk2Last8Cells)
const PlaneGeo * Plane(unsigned int i) const
DEFINE_ART_MODULE(TestTMapFile)
void SetPDGBest(int pdgbest)
T sqr(T x)
More efficient square function than pow(x,2)
Horizontal planes which measure Y.
void SetXMin2(float dxmin2)
minimum position of hit second farthest from edge in that dimension
art::PtrVector< rb::CellHit > AllCells() const
Get all cells from both views.
void SetNMuTrks(int nmutrks)
Calibrated quantities relying on position in the orthogonal view. To generate a rb::CellHit from a rb...
void SetPiMuDecay(int pitomu)
virtual void Add(const art::Ptr< rb::CellHit > &cell, double weight=1)
void SetVtxE20(float vtxE20)
ProductID put(std::unique_ptr< PROD > &&product)
virtual double TotalLength() const
Length (cm) of all the track segments.
unsigned short Cell() const
void produce(art::Event &evt)
int PdgCode(Resonance_t res, int Q)
(resonance id, charge) -> PDG code
void getActivityVtx(const art::PtrVector< rb::CellHit > track1Hits, const art::PtrVector< rb::CellHit > track2Hits, const art::PtrVector< rb::CellHit > sliceHits, float &actvtx)
void SetYMaxt(float dymaxt)
minimum position of hit closest to edge in that dimension that has more than 100 photoelectons (not p...
double CalorimetricEnergy(EEnergyCalcScheme escheme=kRecomputeEnergy) const
Simple estimate of neutrino energy.
void SetHadCellsFromEdge(int hadcellsfromedge)
void SetYMint(float dymint)
minimum position of hit closest to edge in that dimension that has more than 100 photoelectons (not p...
void SetAvedEdxTrack1Last6Cells(float avededxtrk1Last6Cells)
void SetAvedEdxTrack1(float avededxtrk1)
rb::RecoHit RecoHit(const art::Ptr< rb::CellHit > &chit) const
Return calibrated hit based on assumed W coordinate.
void getAveTrackdEdxLast8Cells(const art::Ptr< rb::Track > track, float &avededxlast8)
void SetZMint(float dzmint)
bool IsCalibrated() const
You MUST check here before accessing PECorr, MIP or GeV.
void SetAvedEdxTrack1Last8Cells(float avededxtrk1Last8Cells)
void SetAvedEdxTrack2(float avededxtrk2)
Vertex location in position and time.
Sandbox to make objects for numu analysis in CAF.
void SetZMax2(float dzmax2)
void SetMissingE(float misse)
total E of true particles leaving detector (>1 MeV, not neutron)
void SetYMin2(float dymin2)
numusand::NumuSandFxs * fNumuSandFxs
art::Ptr< rb::CellHit > Cell(geo::View_t view, unsigned int viewIdx) const
Get the ith cell from view view.
std::vector< TrackIDE > HitsToTrackIDE(const std::vector< const rb::CellHit * > &hits, bool useBirksE=false) const
Returns vector of TrackIDE structs contributing to the given collection of hits.
void SetNDropped(int ndropped)
A rawdata::RawDigit with channel information decoded.
void getAveTrackdEdxLast6Cells(const art::Ptr< rb::Track > track, float &avededxlast6)
void SetXMax2(float dxmax2)
maximum position of hit second farthest from edge in that dimension
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
bool HaveTruthInfo() const
Is this a file with truth info in? (Is BackTracker going to be any use to you?)
double TotalGeV(EEnergyCalcScheme escheme=kRecomputeEnergy) const
Simple sum of the estimated GeV of all the hits.
void SetActVtx(float startact)
proxy for hadE[all hits not in main track] - calE[2nd reco trk]
void SetScattAngleTrack2(float scattangletrk2)
bool getMissingE(const art::Ptr< rb::Cluster > slice, float &missE1)
void SetZMin2(float dzmin2)
T min(const caf::Proxy< T > &a, T b)
bool IsNoise() const
Is the noise flag set?
Encapsulate the geometry of one entire detector (near, far, ndos)
virtual bool Is3D() const
void SetVtxE60(float vtxE60)
void SetNProtons(int nprotons)
void SetVtxE40(float vtxE40)
std::vector< const sim::Particle * > HitsToParticle(const std::vector< const rb::CellHit * > &hits) const
Returns vector of sim::Particle objects contributing to the given collection of hits.