24 #include "TDatabasePDG.h" 48 #include "Utilities/func/MathUtil.h" 193 for(
unsigned int hitIdx = 0; hitIdx < hitHandle->size(); hitIdx++){
201 const unsigned int sliceMax = sliceHandle->size();
204 for(
unsigned int sliceIdx = 0; sliceIdx < sliceMax; ++sliceIdx){
208 if ( slice->
IsNoise() )
continue;
212 if(SlicePurr.size()==0)
continue;
217 double sliceEfficiency = SliceNuPurr[0].efficiency;
222 double E = nu.
Nu().
E();
226 int intCurrent = nu.
CCNC();
229 if (
fSkipNuMu && nuPDG != 12 && nuPDG != -12 )
return;
232 if ( lepPDG == 12 || lepPDG == 14 || lepPDG == 16 || lepPDG == 13 || lepPDG == -13 )
235 std::vector<int> hitsNfls;
236 std::vector<double> hitsIsEve;
237 std::vector<double> hitsIsEve_e;
238 std::vector<double> hitsIsEve_mu;
239 std::vector<double> hitsIsEve_tau;
240 std::vector<double> hitsIsEM;
241 std::vector<double> hitsEveIsEM;
242 std::vector<double> hitsIsMu;
243 std::vector<double> hitsEveIsMu;
244 std::vector<double> hitsADC;
245 std::vector<double> hitsPECorr;
246 std::vector<double> hitsTrueE;
247 std::vector<double> hitsRecoE;
250 bool selected =
true;
251 float xmin[2] = { 9999, 9999};
252 float ymin[2] = { 9999, 9999};
253 float zmin[2] = { 9999, 9999};
254 float xmax[2] = {-9999,-9999};
255 float ymax[2] = {-9999,-9999};
256 float zmax[2] = {-9999,-9999};
258 const int interactionNhits = slice->
NCell();
262 for(
unsigned int j = 0;
j < slice->
NCell(geoview); ++
j) {
269 float hitRecoE = rh.
GeV();
270 float hitPECorr = rh.
PECorr();
271 double hitADC = ch->
ADC();
273 const std::vector< sim::FLSHit > CellFLSs = bt->
HitToFLSHit(chptr);
274 int nFLSs = CellFLSs.size();
278 if ( nFLSs == 0 )
continue;
279 if ( nFLSs == 1 ) hitTrueE = CellFLSs[0].GetEdepBirks();
281 for (
int FlsIdx = 0; FlsIdx < nFLSs; ++FlsIdx)
282 hitTrueE += CellFLSs[FlsIdx].GetEdepBirks();
287 double hitIsEve_e = 0;
288 double hitIsEve_mu = 0;
289 double hitIsEve_tau = 0;
291 double hitEveIsEM = 0;
293 double hitEveIsMu = 0;
295 for (
int flsIdx = 0; flsIdx < nFLSs; ++flsIdx){
299 int EveID = pnav.
EveId(trkID);
305 double thisIsEve_tau;
306 double thisIsEM, thisEveIsEM;
307 double thisIsMu, thisEveIsMu;
310 if ( trkID == EveID ) thisIsEve = 1;
312 if ( trkID == EveID && (Eve_p->
PdgCode() == 11 || Eve_p->
PdgCode() == -11) ) thisIsEve_e = 1;
313 else thisIsEve_e = 0;
314 if ( trkID == EveID && (Eve_p->
PdgCode() == 13 || Eve_p->
PdgCode() == -13) ) thisIsEve_mu = 1;
315 else thisIsEve_mu = 0;
316 if ( trkID == EveID && (Eve_p->
PdgCode() == 15 || Eve_p->
PdgCode() == -15) ) thisIsEve_tau = 1;
317 else thisIsEve_tau = 0;
320 Eve_p->
PdgCode() == 22 || Eve_p->
PdgCode() == 111 ) thisEveIsEM = 1;
321 else thisEveIsEM = 0;
326 if ( Eve_p->
PdgCode() == 13 || Eve_p->
PdgCode() == -13 ) thisEveIsMu = 1;
327 else thisEveIsMu = 0;
328 if ( fls_p->
PdgCode() == 13 || fls_p->
PdgCode() == -13 ) thisIsMu = 1;
334 hitIsEve += (thisflsE/hitTrueE)*thisIsEve;
335 hitIsEve_e += (thisflsE/hitTrueE)*thisIsEve_e;
336 hitIsEve_mu += (thisflsE/hitTrueE)*thisIsEve_mu;
337 hitIsEve_tau += (thisflsE/hitTrueE)*thisIsEve_tau;
338 hitIsEM += (thisflsE/hitTrueE)*thisIsEM;
339 hitEveIsEM += (thisflsE/hitTrueE)*thisEveIsEM;
340 hitIsMu += (thisflsE/hitTrueE)*thisIsMu;
341 hitEveIsMu += (thisflsE/hitTrueE)*thisEveIsMu;
345 hitsADC.push_back(hitADC);
346 hitsPECorr.push_back((
double)hitPECorr);
347 hitsNfls.push_back(nFLSs);
348 hitsTrueE.push_back(hitTrueE);
349 hitsRecoE.push_back((
double)hitRecoE);
350 hitsIsEve.push_back(hitIsEve);
351 hitsIsEve_e.push_back(hitIsEve_e);
352 hitsIsEve_mu.push_back(hitIsEve_mu);
353 hitsIsEve_tau.push_back(hitIsEve_tau);
354 hitsIsEM.push_back(hitIsEM);
355 hitsEveIsEM.push_back(hitEveIsEM);
356 hitsIsMu.push_back(hitIsMu);
357 hitsEveIsMu.push_back(hitEveIsMu);
362 const double x = xyz[0];
363 const double y = xyz[1];
364 const double z = xyz[2];
367 if ( x < xmin[0] ) { xmin[1] = xmin[0]; xmin[0] =
x;}
368 else if ( x < xmin[1] ) xmin[1] =
x;
369 if ( z < zmin[0] ) { zmin[1] = zmin[0]; zmin[0] =
z;}
370 else if ( z < zmin[1] ) zmin[1] =
z;
371 if ( x > xmax[0] ) { xmax[1] = xmax[0]; xmax[0] =
x;}
372 else if ( x > xmax[1] ) xmax[1] =
x;
373 if ( z > zmax[0] ) { zmax[1] = zmax[0]; zmax[0] =
z;}
374 else if ( z > zmax[1] ) zmax[1] = zmax[0];
377 if ( y < ymin[0] ) { ymin[1] = ymin[0]; ymin[0] =
y;}
378 else if ( y < ymin[1] ) ymin[1] =
y;
379 if ( z < zmin[0] ) { zmin[1] = zmin[0]; zmin[0] =
z;}
380 else if ( z < zmin[1] ) zmin[1] =
z;
381 if ( y > ymax[0] ) { ymax[1] = ymax[0]; ymax[0] =
y;}
382 else if ( y > ymax[1] ) ymax[1] =
y;
383 if ( z > zmax[0] ) { zmax[1] = zmax[0]; zmax[0] =
z;}
384 else if ( z > zmax[1] ) zmax[1] = zmax[0];
400 double TotalRecoE =0;
401 double TotalRecoE_EM =0;
421 TotalRecoE += (1.0 -EveMuness)*(
HitsRecoE[
h]);
425 double InteractionHitRecoE = TotalRecoE;
427 double R_emRE_HitsRE = TotalRecoE_EM/InteractionHitRecoE;
439 if( !selected )
continue;
double E(const int i=0) const
back track the reconstruction to the simulation
double IntSliceEfficiency
unsigned int NCell(geo::View_t view) const
Number of cells in view view.
const simb::MCNeutrino & GetNeutrino() const
std::map< std::string, double > xmax
ShowerEnergyAna(const fhicl::ParameterSet &pset)
std::vector< NeutrinoEffPur > SliceToNeutrinoInteractions(const std::vector< const rb::CellHit * > &sliceHits, const std::vector< const rb::CellHit * > &allHits, bool sortPur=false) const
Given a collection of hits (often a slice), returns vector of structures of neutrino interactions...
const sim::ParticleNavigator & ParticleNavigator() const
Get a reference to the ParticleNavigator.
enum geo::_plane_proj View_t
Enumerate the possible plane projections.
unsigned short Plane() const
const simb::MCParticle & Nu() const
Vertical planes which measure X.
A single unit of energy deposition in the liquid scintillator.
const PlaneGeo * Plane(unsigned int i) const
DEFINE_ART_MODULE(TestTMapFile)
void analyze(const art::Event &evt)
Horizontal planes which measure Y.
art::PtrVector< rb::CellHit > AllCells() const
Get all cells from both views.
Calibrated quantities relying on position in the orthogonal view. To generate a rb::CellHit from a rb...
unsigned short Cell() const
int InteractionType() const
void reconfigure(const fhicl::ParameterSet &pset)
const simb::MCParticle & Lepton() const
void push_back(Ptr< U > const &p)
std::vector< sim::FLSHit > HitToFLSHit(const rb::CellHit &hit) const
All the FLSHits that contributed to this hit, sorted from most to least light.
T get(std::string const &key) const
rb::RecoHit RecoHit(const art::Ptr< rb::CellHit > &chit) const
Return calibrated hit based on assumed W coordinate.
Collect Geo headers and supply basic geometry functions.
bool IsCalibrated() const
You MUST check here before accessing PECorr, MIP or GeV.
bool IsFiltered(const art::Event &evt, art::Ptr< T > x, const std::vector< std::string > &labels)
Is this Ptr marked "filtered out"?
EDAnalyzer(Table< Config > const &config)
float GetEdepBirks() const
Get total Energy with Birks suppression deposited into the cell for the whole FLSHit.
const sim::Particle * TrackIDToParticle(int const &id) const
Returns a pointer to the sim::Particle object corresponding to the given TrackID. ...
double HitsIsEve_mu[1000]
art::Ptr< rb::CellHit > Cell(geo::View_t view, unsigned int viewIdx) const
Get the ith cell from view view.
A rawdata::RawDigit with channel information decoded.
T * make(ARGS...args) const
double IntSliceEMfraction
int GetTrackID() const
Track ID.
int16_t ADC(uint32_t i) const
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
bool IsNoise() const
Is the noise flag set?
Event generator information.
std::vector< NeutrinoEffPur > SliceToMCTruth(const std::vector< const rb::CellHit * > &sliceHits, const std::vector< const rb::CellHit * > &allHits, bool sortPur=false) const
Given a collection of hits (often a slice), returns vector of structures of MCTruth, efficiency, and purity of that neutrino interaction relative to the slice. Efficiency is defined as FLS energy from neutrino interaction in slice / total FLS energy from neutrino interaction in event. This vector is sorted from the highest efficiency interaction to lowest. This function returns all MCTruth, including those without neutrino interactions, i.e. cosmics.
double HitsIsEve_tau[1000]
Encapsulate the geometry of one entire detector (near, far, ndos)
int EveId(const int trackID) const