All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Public Member Functions | Private Member Functions | Private Attributes | List of all members
cluster::ClusterAna Class Reference
Inheritance diagram for cluster::ClusterAna:

Public Member Functions

 ClusterAna (fhicl::ParameterSet const &pset)
 

Private Member Functions

void analyze (const art::Event &evt) override
 
void beginJob () override
 
void endJob () override
 

Private Attributes

TH1F * fNClusters
 
TH1F * fNHitInCluster
 
TH1F * fCREP2
 
TH1F * fCRE
 
TH1F * fCRP
 
TH1F * fNuKE_elec
 
TH1F * fNuKE_muon
 
TH1F * fNuKE_pion
 
TH1F * fNuKE_kaon
 
TH1F * fNuKE_prot
 
TH1F * fNuEP2_elec
 
TH1F * fNuEP2_muon
 
TH1F * fNuEP2_pion
 
TH1F * fNuEP2_kaon
 
TH1F * fNuEP2_prot
 
TH1F * fNuE_elec
 
TH1F * fNuE_muon
 
TH1F * fNuE_pion
 
TH1F * fNuE_kaon
 
TH1F * fNuE_prot
 
TH1F * fNuP_elec
 
TH1F * fNuP_muon
 
TH1F * fNuP_pion
 
TH1F * fNuP_kaon
 
TH1F * fNuP_prot
 
TH1F * fNuVtx_dx
 
TH1F * fNuVtx_dy
 
TH1F * fNuVtx_dz
 
TProfile * fNuEP2_KE_elec
 
TProfile * fNuEP2_KE_muon
 
TProfile * fNuEP2_KE_pion
 
TProfile * fNuEP2_KE_kaon
 
TProfile * fNuEP2_KE_prot
 
std::string fHitsModuleLabel
 
std::string fClusterModuleLabel
 
std::string fVertexModuleLabel
 
std::vector< float > fElecKERange
 
std::vector< float > fMuonKERange
 
std::vector< float > fPionKERange
 
std::vector< float > fKaonKERange
 
std::vector< float > fProtKERange
 
short fTrackWeightOption
 
bool fMergeDaughters
 
bool fSkipCosmics
 
short fPrintLevel
 
short moduleID
 
std::ofstream outFile
 

Detailed Description

Definition at line 44 of file ClusterAna_module.cc.

Constructor & Destructor Documentation

cluster::ClusterAna::ClusterAna ( fhicl::ParameterSet const &  pset)
explicit

Definition at line 116 of file ClusterAna_module.cc.

117  : EDAnalyzer(pset)
118  , fHitsModuleLabel(pset.get<std::string>("HitsModuleLabel"))
119  , fClusterModuleLabel(pset.get<std::string>("ClusterModuleLabel"))
120  , fVertexModuleLabel(pset.get<std::string>("VertexModuleLabel"))
121  , fElecKERange(pset.get<std::vector<float>>("ElecKERange"))
122  , fMuonKERange(pset.get<std::vector<float>>("MuonKERange"))
123  , fPionKERange(pset.get<std::vector<float>>("PionKERange"))
124  , fKaonKERange(pset.get<std::vector<float>>("KaonKERange"))
125  , fProtKERange(pset.get<std::vector<float>>("ProtKERange"))
126  , fTrackWeightOption(pset.get<short>("TrackWeightOption"))
127  , fMergeDaughters(pset.get<bool>("MergeDaughters"))
128  , fSkipCosmics(pset.get<bool>("SkipCosmics"))
129  , fPrintLevel(pset.get<short>("PrintLevel"))
130  {
131 
132  if (fPrintLevel == -1) {
133  // encode the clustermodule label into an integer
134  moduleID = 0;
135  size_t found = fClusterModuleLabel.find("traj");
136  if (found != std::string::npos) moduleID = 1;
137  found = fClusterModuleLabel.find("line");
138  if (found != std::string::npos) moduleID = 2;
139  found = fClusterModuleLabel.find("fuzz");
140  if (found != std::string::npos) moduleID = 3;
141  found = fClusterModuleLabel.find("pand");
142  if (found != std::string::npos) moduleID = 4;
143  std::cout << "fClusterModuleLabel " << fClusterModuleLabel << " ID " << moduleID << "\n";
144 
145  std::string fileName = fClusterModuleLabel + ".tru";
146  outFile.open(fileName);
147  }
148  }
cluster::ClusterAna::fHitsModuleLabel
std::string fHitsModuleLabel
Definition: ClusterAna_module.cc:92
cluster::ClusterAna::fElecKERange
std::vector< float > fElecKERange
Definition: ClusterAna_module.cc:95
cluster::ClusterAna::fProtKERange
std::vector< float > fProtKERange
Definition: ClusterAna_module.cc:99
cluster::ClusterAna::fKaonKERange
std::vector< float > fKaonKERange
Definition: ClusterAna_module.cc:98
cluster::ClusterAna::fVertexModuleLabel
std::string fVertexModuleLabel
Definition: ClusterAna_module.cc:94
fileName
services TFileService fileName
Definition: prodcorsika_genie_protononly_icarus_numi_nue_volDetEnclosure_tpc.fcl:18
cluster::ClusterAna::fClusterModuleLabel
std::string fClusterModuleLabel
Definition: ClusterAna_module.cc:93
cluster::ClusterAna::fPionKERange
std::vector< float > fPionKERange
Definition: ClusterAna_module.cc:97
cluster::ClusterAna::fMuonKERange
std::vector< float > fMuonKERange
Definition: ClusterAna_module.cc:96
cout
BEGIN_PROLOG could also be cout
Definition: messageservice.fcl:10

Member Function Documentation

void cluster::ClusterAna::analyze ( const art::Event &  evt)
overrideprivate

Definition at line 211 of file ClusterAna_module.cc.

212  {
213  // code stolen from TrackAna_module.cc
214  art::ServiceHandle<geo::Geometry const> geom;
215  if (geom->Nplanes() > 3) {
216  mf::LogError("ClusterAna") << "Too many planes for this code...";
217  return;
218  }
219 
220  // get all hits in the event
221  art::Handle<std::vector<recob::Hit>> hitListHandle;
222  evt.getByLabel(fHitsModuleLabel, hitListHandle);
223  std::vector<art::Ptr<recob::Hit>> allhits;
224  art::fill_ptr_vector(allhits, hitListHandle);
225  if (allhits.size() == 0) return;
226 
227  // get clusters and cluster-hit associations
228  art::Handle<std::vector<recob::Cluster>> clusterListHandle;
229  evt.getByLabel(fClusterModuleLabel, clusterListHandle);
230  art::FindManyP<recob::Hit> fmh(clusterListHandle, evt, fClusterModuleLabel);
231  if (clusterListHandle->size() == 0) {
232  std::cout << "No clusters in event. Hits size " << allhits.size() << "\n";
233  // don't return or the statistics will be off
234  }
235 
236  // get 3D vertices
237  art::Handle<std::vector<recob::Vertex>> vertexListHandle;
238  double xyz[3] = {0, 0, 0};
239  art::PtrVector<recob::Vertex> recoVtxList;
240  if (vertexListHandle.isValid()) {
241  evt.getByLabel(fVertexModuleLabel, vertexListHandle);
242  for (unsigned int ii = 0; ii < vertexListHandle->size(); ++ii) {
243  art::Ptr<recob::Vertex> vertex(vertexListHandle, ii);
244  recoVtxList.push_back(vertex);
245  vertex->XYZ(xyz);
246  } // ii
247  } // vertexListHandle
248 
249  // list of all true particles
250  art::ServiceHandle<cheat::BackTrackerService const> bt_serv;
251  art::ServiceHandle<cheat::ParticleInventoryService const> pi_serv;
252  sim::ParticleList const& plist = pi_serv->ParticleList();
253  // list of all true particles that will be considered
254  std::vector<const simb::MCParticle*> plist2;
255  // true (reconstructed) hits for each particle in plist2
256  std::vector<std::vector<art::Ptr<recob::Hit>>> hlist2;
257  // index of cluster matched to each particle in plist2 in each plane
258  std::vector<std::vector<short>> truToCl;
259  // number of true hits in each plane and cluster
260  std::vector<std::vector<unsigned short>> nTruHitInCl;
261  //number of reconstructed hits in all clusters
262  std::vector<unsigned short> nRecHitInCl;
263 
264  // calculate average EP2 for every event to facilitate code development
265  // Beam Neutrinos - muons and not-muons
266 
267  float aveNuEP2mu = 0.;
268  float numNuEP2mu = 0.;
269  float aveNuEP2nm = 0.;
270  float numNuEP2nm = 0.;
271  // Cosmic Rays
272  float aveCREP2 = 0.;
273  float numCREP2 = 0.;
274 
275  // track ID of the neutrino (or single particle)
276  int neutTrackID = -1;
277  std::vector<int> tidlist;
278  float neutEnergy = -1.;
279  int neutIntType = -1;
280  int neutCCNC = -1;
281  bool isNeutrino = false;
282  bool isSingleParticle = false;
283 
284  auto const clockData = art::ServiceHandle<detinfo::DetectorClocksService const>()->DataFor(evt);
285 
286  for (sim::ParticleList::const_iterator ipart = plist.begin(); ipart != plist.end(); ++ipart) {
287  const simb::MCParticle* part = (*ipart).second;
288  assert(part != 0);
289  int pdg = abs(part->PdgCode());
290  int trackID = part->TrackId();
291  art::Ptr<simb::MCTruth> theTruth = pi_serv->TrackIdToMCTruth_P(trackID);
292  if (fSkipCosmics && theTruth->Origin() == simb::kCosmicRay) continue;
293  if (theTruth->Origin() == simb::kBeamNeutrino) isNeutrino = true;
294  if (theTruth->Origin() == simb::kSingleParticle) isSingleParticle = true;
295  float KE = 1000 * (part->E() - part->Mass());
296 
297  // Get the neutrino track ID. Assume that there is only one neutrino
298  // interaction and it is first in the list of BeamNeutrino particles
299  if ((isNeutrino || isSingleParticle) && neutTrackID < 0) {
300  neutTrackID = trackID;
301  if (isNeutrino) {
302  simb::MCNeutrino theNeutrino = theTruth->GetNeutrino();
303  neutEnergy = 1000. * theNeutrino.Nu().E();
304  neutIntType = theNeutrino.InteractionType();
305  neutCCNC = theNeutrino.CCNC();
306  }
307 
308  for (unsigned short iv = 0; iv < recoVtxList.size(); ++iv) {
309  recoVtxList[iv]->XYZ(xyz);
310  fNuVtx_dx->Fill(part->Vx() - xyz[0]);
311  fNuVtx_dy->Fill(part->Vy() - xyz[1]);
312  fNuVtx_dz->Fill(part->Vz() - xyz[2]);
313  } // iv
314  } // theTruth->Origin() == simb::kBeamNeutrino && neutTrackID <
315 
316  if (fPrintLevel > 3)
317  mf::LogVerbatim("ClusterAna") << "Origin " << theTruth->Origin() << " trackID " << trackID
318  << " PDG " << part->PdgCode() << " KE " << (int)KE << " MeV "
319  << " neutTrackID " << neutTrackID << " Mother "
320  << part->Mother() << " Proc " << part->Process();
321 
322  bool isCharged = (pdg == 11) || (pdg == 13) || (pdg == 211) || (pdg == 321) || (pdg == 2212);
323 
324  if (!isCharged) continue;
325 
326  // KE (MeV) cuts
327  if (pdg == 11) {
328  if (fElecKERange[0] < 0) continue;
329  // only allow primary electrons
330  if (part->Process() != "primary") continue;
331  if (KE < fElecKERange[0] || KE > fElecKERange[1]) continue;
332  }
333  if (pdg == 13) {
334  if (fMuonKERange[0] < 0) continue;
335  if (KE < fMuonKERange[0] || KE > fMuonKERange[1]) continue;
336  }
337  if (pdg == 211) {
338  if (fPionKERange[0] < 0) continue;
339  if (KE < fPionKERange[0] || KE > fPionKERange[1]) continue;
340  }
341  if (pdg == 321) {
342  if (fKaonKERange[0] < 0) continue;
343  if (KE < fKaonKERange[0] || KE > fKaonKERange[1]) continue;
344  }
345  if (pdg == 2212) {
346  if (fProtKERange[0] < 0) continue;
347  if (KE < fProtKERange[0] || KE > fProtKERange[1]) continue;
348  }
349  // ignore secondaries from neutron interactions
350  if (part->Process() == "NeutronInelastic") continue;
351  plist2.push_back(part);
352  tidlist.push_back(trackID);
353  // initialize the true->(cluster,plane) association
354  std::vector<short> temp{-1, -1, -1};
355  truToCl.push_back(temp);
356  // initialize the true hit count
357  std::vector<unsigned short> temp2(3);
358  nTruHitInCl.push_back(temp2);
359 
360  if (fPrintLevel > 2)
361  mf::LogVerbatim("ClusterAna")
362  << "plist2[" << plist2.size() - 1 << "]"
363  << " Origin " << theTruth->Origin() << " trackID " << trackID << " PDG "
364  << part->PdgCode() << " KE " << (int)KE << " Mother " << part->Mother() + neutTrackID - 1
365  << " Proc " << part->Process();
366  }
367 
368  if (plist2.size() == 0) return;
369 
370  // get the hits (in all planes) that are matched to the true tracks
371  hlist2 = bt_serv->TrackIdsToHits_Ps(clockData, tidlist, allhits);
372  if (hlist2.size() != plist2.size()) {
373  mf::LogError("ClusterAna") << "MC particle list size " << plist2.size()
374  << " != size of MC particle true hits lists " << hlist2.size();
375  return;
376  }
377  tidlist.clear();
378 
379  // vector of (mother, daughter) pairs
380  std::vector<std::pair<unsigned short, unsigned short>> moda;
381  // Deal with mother-daughter tracks
382  if (fMergeDaughters && neutTrackID >= 0) {
383  // Assume that daughters appear later in the list. Step backwards
384  // to accumulate all generations of daughters
385  unsigned short ii, dpl, jj, jpl, kpl;
386  for (ii = 0; ii < plist2.size(); ++ii) {
387  dpl = plist2.size() - 1 - ii;
388  // no mother
389  if (plist2[dpl]->Mother() == 0) continue;
390  // electron
391  if (abs(plist2[dpl]->PdgCode()) == 11) continue;
392  // the actual mother trackID is offset from the neutrino trackID
393  int motherID = neutTrackID + plist2[dpl]->Mother() - 1;
394  // ensure that we are only looking at BeamNeutrino or single particle daughters
395  if (motherID != neutTrackID) continue;
396  // count the number of daughters
397  int ndtr = 0;
398  for (kpl = 0; kpl < plist2.size(); ++kpl) {
399  if (plist2[kpl]->Mother() == motherID) ++ndtr;
400  }
401  // require only one daughter
402  if (ndtr > 1) continue;
403  // find the mother in the list
404  int mpl = -1;
405  for (jj = 0; jj < plist2.size(); ++jj) {
406  jpl = plist2.size() - 1 - jj;
407  if (plist2[jpl]->TrackId() == motherID) {
408  mpl = jpl;
409  break;
410  }
411  } // jpl
412  // mother not found for some reason
413  if (mpl < 0) {
414  mf::LogVerbatim("ClusterAna")
415  << " mother of daughter " << dpl << " not found. mpl = " << mpl;
416  continue;
417  }
418  // ensure that PDG code for mother and daughter are the same
419  if (plist2[dpl]->PdgCode() != plist2[mpl]->PdgCode()) continue;
420  moda.push_back(std::make_pair(mpl, dpl));
421  } // dpl
422  } // MergeDaughters
423 
424  // Now match reconstructed clusters to true particles.
425  art::PtrVector<recob::Cluster> clusters;
426  for (unsigned int ii = 0; ii < clusterListHandle->size(); ++ii) {
427  art::Ptr<recob::Cluster> clusterHolder(clusterListHandle, ii);
428  clusters.push_back(clusterHolder);
429  }
430 
431  fNClusters->Fill(clusterListHandle->size());
432  nRecHitInCl.resize(clusters.size());
433 
434  // get the plane from the view. Perhaps there is a method that does
435  // this somewhere...
436  std::map<geo::View_t, unsigned int> ViewToPlane;
437  for (unsigned int plane = 0; plane < geom->Nplanes(); ++plane) {
438  geo::View_t view = geom->Plane(plane).View();
439  ViewToPlane[view] = plane;
440  }
441  for (size_t icl = 0; icl < clusters.size(); ++icl) {
442  unsigned int plane = ViewToPlane[clusters[icl]->View()];
443  std::vector<art::Ptr<recob::Hit>> cluhits = fmh.at(icl);
444  fNHitInCluster->Fill(cluhits.size());
445  nRecHitInCl[icl] = cluhits.size();
446  // count the number of hits matched to each true particle in plist2
447  std::vector<unsigned short> nHitInPl2(plist2.size());
448  for (size_t iht = 0; iht < cluhits.size(); ++iht) {
449 
450  // look for this hit in all of the truth hit lists
451  short hitInPl2 = -1;
452  for (unsigned short ipl = 0; ipl < plist2.size(); ++ipl) {
453  unsigned short imat = 0;
454  for (imat = 0; imat < hlist2[ipl].size(); ++imat) {
455  if (cluhits[iht] == hlist2[ipl][imat]) break;
456  } // imat
457  if (imat < hlist2[ipl].size()) {
458  hitInPl2 = ipl;
459  break;
460  }
461  } // ipl
462  if (hitInPl2 < 0) continue;
463  // Assign the hit count to the mother if this is a daughter.
464  // Mother-daughter pairs are entered in the moda vector in reverse
465  // order, so assign daughter hits to the highest generation mother.
466  for (unsigned short imd = 0; imd < moda.size(); ++imd) {
467  if (moda[imd].second == hitInPl2) hitInPl2 = moda[imd].first;
468  }
469  // count
470  ++nHitInPl2[hitInPl2];
471  } // iht
472  // Associate the cluster with the truth particle that has the highest
473  // number of cluster hits
474  unsigned short nhit = 0;
475  short imtru = -1;
476  for (unsigned int ipl = 0; ipl < nHitInPl2.size(); ++ipl) {
477  if (nHitInPl2[ipl] > nhit) {
478  nhit = nHitInPl2[ipl];
479  imtru = ipl;
480  }
481  } // ipl
482  // make the cluster->(true,plane) association and save the
483  // number of true hits in the cluster
484  if (imtru != -1) {
485  // clobber a previously made association?
486  if (nhit > nTruHitInCl[imtru][plane]) {
487  truToCl[imtru][plane] = icl;
488  nTruHitInCl[imtru][plane] = nhit;
489  }
490  } // imtru != 1
491  } // icl
492 
493  // ready to calculate Efficiency, Purity in each plane and EP2
494  for (unsigned short ipl = 0; ipl < plist2.size(); ++ipl) {
495  // ignore daughters
496  bool skipit = false;
497  for (unsigned short ii = 0; ii < moda.size(); ++ii) {
498  if (moda[ii].second == ipl) {
499  skipit = true;
500  break;
501  }
502  } // ii
503  if (skipit) continue;
504  // ignore true particles with few true hits. Outside the detector
505  // or not reconstructable?
506  if (hlist2[ipl].size() < 3) continue;
507 
508  int trackID = plist2[ipl]->TrackId();
509  art::Ptr<simb::MCTruth> theTruth = pi_serv->TrackIdToMCTruth_P(trackID);
510  bool isCosmic = (theTruth->Origin() == simb::kCosmicRay);
511  float KE = 1000 * (plist2[ipl]->E() - plist2[ipl]->Mass());
512  int PDG = abs(plist2[ipl]->PdgCode());
513 
514  std::vector<short> nTru(geom->Nplanes());
515  std::vector<short> nRec(geom->Nplanes());
516  std::vector<short> nTruRec(geom->Nplanes());
517  std::vector<float> eff(geom->Nplanes());
518  std::vector<float> pur(geom->Nplanes());
519  std::vector<float> ep(geom->Nplanes());
520  for (unsigned int plane = 0; plane < geom->Nplanes(); ++plane) {
521  // count the number of true hits in this plane for the true particle.
522  // First count the mother hits
523  for (unsigned short ii = 0; ii < hlist2[ipl].size(); ++ii) {
524  if (ViewToPlane[hlist2[ipl][ii]->View()] == plane) ++nTru[plane];
525  } // ii
526 
527  // next look for daughters and count those hits in all generations
528  unsigned short mom = ipl;
529  std::vector<std::pair<unsigned short, unsigned short>>::reverse_iterator rit =
530  moda.rbegin();
531  while (rit != moda.rend()) {
532  if ((*rit).first == mom) {
533  unsigned short dau = (*rit).second;
534  for (unsigned short jj = 0; jj < hlist2[dau].size(); ++jj) {
535  if (ViewToPlane[hlist2[dau][jj]->View()] == plane) ++nTru[plane];
536  } // jj
537  // It is likely that one hit appears in the mother list
538  // as well as the daughter list, so subtract one from the count
539  mom = (*rit).second;
540  } // (*rit).first == mom
541  ++rit;
542  } // rit
543 
544  if (nTru[plane] == 0) continue;
545  // Ignore short truth clusters
546  if (nTru[plane] < 3) continue;
547  short icl = truToCl[ipl][plane];
548  if (icl < 0) { nRec[plane] = 0; }
549  else {
550  nRec[plane] = nRecHitInCl[icl];
551  }
552  nTruRec[plane] = nTruHitInCl[ipl][plane];
553  eff[plane] = 0;
554  pur[plane] = 0;
555  if (nTru[plane] > 0) eff[plane] = (float)nTruRec[plane] / (float)nTru[plane];
556  if (nRec[plane] > 0) pur[plane] = (float)nTruRec[plane] / (float)nRec[plane];
557  // do this to prevent histogram under/over flow
558  if (eff[plane] == 0) { eff[plane] = 0.001; }
559  if (pur[plane] == 0) { pur[plane] = 0.001; }
560  if (eff[plane] >= 1) { eff[plane] = 0.999; }
561  if (pur[plane] >= 1) { pur[plane] = 0.999; }
562  ep[plane] = eff[plane] * pur[plane];
563  if (fPrintLevel == -1)
564  outFile << moduleID << " " << evt.id().event() << " " << trackID << " " << PDG << " "
565  << (int)KE << " " << plane << " " << nTru[plane] << " " << std::setprecision(3)
566  << eff[plane] << " " << pur[plane] << "\n";
567  } // plane
568  // sort the ep values in ascending order
569  std::vector<float> temp;
570  temp = ep;
571  std::sort(temp.begin(), temp.end());
572  // EP2 is the second highest value
573  unsigned short ii = temp.size() - 2;
574  float ep2 = temp[ii];
575  // find the plane that defined EP2
576  short ep2Plane = 0;
577  short ep2Cluster = 0;
578  for (unsigned short jj = 0; jj < temp.size(); ++jj) {
579  if (ep[jj] == ep2) {
580  ep2Plane = jj;
581  ep2Cluster = truToCl[ipl][ep2Plane];
582  break;
583  }
584  } // jj
585  // find the US and DS ends of the cluster for printing
586  std::array<double, 2> clBeg, clEnd;
587  if (ep2Cluster >= 0) {
588  clBeg[0] = clusters[ep2Cluster]->StartWire();
589  clBeg[1] = clusters[ep2Cluster]->StartTick();
590  clEnd[0] = clusters[ep2Cluster]->EndWire();
591  clEnd[1] = clusters[ep2Cluster]->EndTick();
592  }
593  else {
594  clBeg.fill(0.);
595  clEnd.fill(0.);
596  }
597  // fill histograms
598  if (isCosmic) {
599  fCREP2->Fill(ep2);
600  fCRE->Fill(eff[ep2Plane]);
601  fCRP->Fill(pur[ep2Plane]);
602  aveCREP2 += ep2;
603  numCREP2 += 1.;
604  if (fPrintLevel > 1)
605  mf::LogVerbatim("ClusterAna")
606  << ">>>CREP2 " << std::fixed << std::setprecision(2) << ep2 << " E " << eff[ep2Plane]
607  << std::setprecision(2) << " P " << pur[ep2Plane] << " P:W:T " << ep2Plane << ":"
608  << (int)clBeg[0] << ":" << (int)clBeg[1] << "-" << ep2Plane << ":" << (int)clEnd[0]
609  << ":" << (int)clEnd[1] << " PDG " << PDG << " KE " << (int)KE << " MeV";
610  } // isCosmic
611  else {
612  float wght = 1.;
613  if (fTrackWeightOption == 1) wght = KE;
614  // accumulate statistics for muons and not-muons
615  if (PDG == 13) {
616  aveNuEP2mu += ep2 * wght;
617  numNuEP2mu += wght;
618  }
619  else {
620  aveNuEP2nm += ep2 * wght;
621  numNuEP2nm += wght;
622  }
623  if (PDG == 11) {
624  fNuKE_elec->Fill(KE, wght);
625  fNuE_elec->Fill(eff[ep2Plane], wght);
626  fNuP_elec->Fill(pur[ep2Plane], wght);
627  fNuEP2_elec->Fill(ep2, wght);
628  fNuEP2_KE_elec->Fill(KE, ep2, wght);
629  }
630  else if (PDG == 13) {
631  fNuKE_muon->Fill(KE, wght);
632  fNuE_muon->Fill(eff[ep2Plane], wght);
633  fNuP_muon->Fill(pur[ep2Plane], wght);
634  fNuEP2_muon->Fill(ep2, wght);
635  fNuEP2_KE_muon->Fill(KE, ep2, wght);
636  }
637  else if (PDG == 211) {
638  fNuKE_pion->Fill(KE, wght);
639  fNuE_pion->Fill(eff[ep2Plane], wght);
640  fNuP_pion->Fill(pur[ep2Plane], wght);
641  fNuEP2_pion->Fill(ep2, wght);
642  fNuEP2_KE_pion->Fill(KE, ep2, wght);
643  }
644  else if (PDG == 321) {
645  fNuKE_kaon->Fill(KE, wght);
646  fNuE_kaon->Fill(eff[ep2Plane], wght);
647  fNuP_kaon->Fill(pur[ep2Plane], wght);
648  fNuEP2_kaon->Fill(ep2, wght);
649  fNuEP2_KE_kaon->Fill(KE, ep2, wght);
650  }
651  else if (PDG == 2212) {
652  fNuKE_prot->Fill(KE, wght);
653  fNuE_prot->Fill(eff[ep2Plane], wght);
654  fNuP_prot->Fill(pur[ep2Plane], wght);
655  fNuEP2_prot->Fill(ep2, wght);
656  fNuEP2_KE_prot->Fill(KE, ep2, wght);
657  }
658  if (fPrintLevel > 1)
659  mf::LogVerbatim("ClusterAna")
660  << ">>>NuEP2 " << std::fixed << std::setprecision(2) << ep2 << " E " << eff[ep2Plane]
661  << std::setprecision(2) << " P " << pur[ep2Plane] << " P:W:T " << ep2Plane << ":"
662  << (int)clBeg[0] << ":" << (int)clBeg[1] << "-" << ep2Plane << ":" << (int)clEnd[0]
663  << ":" << (int)clEnd[1] << " PDG " << PDG << " KE " << (int)KE << " MeV ";
664  if (fPrintLevel > 2) {
665  // print out the begin/end true hits
666  mf::LogVerbatim mfp("ClusterAna");
667  mfp << " Truth P:W:T ";
668  for (unsigned int plane = 0; plane < geom->Nplanes(); ++plane) {
669  unsigned short loW = 9999;
670  unsigned short loT = 0;
671  unsigned short hiW = 0;
672  unsigned short hiT = 0;
673  for (unsigned short ii = 0; ii < hlist2[ipl].size(); ++ii) {
674  if (ViewToPlane[hlist2[ipl][ii]->View()] == plane) {
675  art::Ptr<recob::Hit> theHit = hlist2[ipl][ii];
676  if (theHit->WireID().Wire < loW) {
677  loW = theHit->WireID().Wire;
678  loT = theHit->PeakTime();
679  }
680  if (theHit->WireID().Wire > hiW) {
681  hiW = theHit->WireID().Wire;
682  hiT = theHit->PeakTime();
683  }
684  } // correct view
685  } // ii
686  mfp << plane << ":" << loW << ":" << loT << "-" << plane << ":" << hiW << ":" << hiT
687  << " ";
688  } // plane
689  } // fPrintLevel > 2
690  } // !isCosmic
691  } // ipl
692 
693  float ave1 = -1.;
694  if (numNuEP2mu > 0.) ave1 = aveNuEP2mu / numNuEP2mu;
695 
696  float ave2 = -1.;
697  if (numNuEP2nm > 0.) ave2 = aveNuEP2nm / numNuEP2nm;
698 
699  float ave3 = -1.;
700  if (numCREP2 > 0.) ave3 = aveCREP2 / numCREP2;
701 
702  if (fPrintLevel > 0) {
703  std::string nuType = "Other";
704  if (neutCCNC == simb::kCC) {
705  if (neutIntType == 1001) nuType = "CCQE";
706  if (neutIntType == 1091) nuType = "DIS";
707  if (neutIntType == 1097) nuType = "COH";
708  if (neutIntType > 1002 && neutIntType < 1091) nuType = "RES";
709  }
710  else if (neutCCNC == simb::kNC) {
711  nuType = "NC";
712  }
713  else {
714  nuType = "Unknown";
715  }
716  mf::LogVerbatim("ClusterAna")
717  << "EvtEP2 " << evt.id().event() << " NuType " << nuType << " Enu " << std::fixed
718  << std::setprecision(0) << neutEnergy << std::right << std::fixed << std::setprecision(2)
719  << " NuMuons " << ave1 << " NuPiKp " << ave2 << " CosmicRays " << ave3 << " CCNC "
720  << neutCCNC << " IntType " << neutIntType;
721  }
722  } // analyze
vertex
process_name vertex
Definition: cheaterreco.fcl:51
PDG
process_name stream1 can override from command line with o or output services DetectorPropertiesService services DetectorPropertiesService services DetectorPropertiesService services DetectorPropertiesService physics analyzers pmtresponse NeutronTrackingCut services LArG4Parameters gaussian physics producers generator PDG
Definition: prodsingle_libraryana.fcl:81
if
then if[["$THISISATEST"==1]]
Definition: neoSmazza.sh:95
cluster::ClusterAna::fHitsModuleLabel
std::string fHitsModuleLabel
Definition: ClusterAna_module.cc:92
second
second
Definition: detectorproperties_icarus.fcl:8
tca::TCEvent::event
unsigned int event
Definition: DataStructs.h:634
pdg
var pdg
Definition: selectors.fcl:14
geo::View_t
enum geo::_plane_proj View_t
Enumerate the possible plane projections.
right
walls no right
Definition: selectors.fcl:105
icarus::ns::util::size
std::size_t size(FixedBins< T, C > const &) noexcept
Definition: FixedBins.h:561
cluster::ClusterAna::fElecKERange
std::vector< float > fElecKERange
Definition: ClusterAna_module.cc:95
abs
T abs(T value)
Definition: sparse_vector_test.cc:30
ana::Current::kCC
Charged-current interactions.
Definition: IPrediction.h:38
util::quantities::concepts::Quantity
A value measured in the specified unit.
Definition: quantities.h:566
cluster::ClusterAna::fProtKERange
std::vector< float > fProtKERange
Definition: ClusterAna_module.cc:99
cluster::ClusterAna::fKaonKERange
std::vector< float > fKaonKERange
Definition: ClusterAna_module.cc:98
kCosmicRay
const Cut kCosmicRay
Definition: PAC2020/SBND/helper.h:54
cluster::ClusterAna::fVertexModuleLabel
std::string fVertexModuleLabel
Definition: ClusterAna_module.cc:94
ana::Current::kNC
Neutral-current interactions.
Definition: IPrediction.h:39
cluster::ClusterAna::fClusterModuleLabel
std::string fClusterModuleLabel
Definition: ClusterAna_module.cc:93
tca::evt
TCEvent evt
Definition: DataStructs.cxx:8
else
else
Definition: generate_icarus_geometry_splitwire.pl:39
cluster::ClusterAna::fPionKERange
std::vector< float > fPionKERange
Definition: ClusterAna_module.cc:97
cluster::ClusterAna::fMuonKERange
std::vector< float > fMuonKERange
Definition: ClusterAna_module.cc:96
cout
BEGIN_PROLOG could also be cout
Definition: messageservice.fcl:10
void cluster::ClusterAna::beginJob ( )
overrideprivate

Definition at line 152 of file ClusterAna_module.cc.

153  {
154 
155  // get access to the TFile service
156  art::ServiceHandle<art::TFileService const> tfs;
157 
158  fNClusters = tfs->make<TH1F>("fNoClustersInEvent", "Number of Clusters", 40, 0, 400);
159  fNHitInCluster = tfs->make<TH1F>("fNHitInCluster", "NHitInCluster", 100, 0, 100);
160 
161  if (!fSkipCosmics) {
162  // Cosmic ray histos
163  fCREP2 = tfs->make<TH1F>("CREP2", "CREP2", 50, 0, 1);
164  fCRE = tfs->make<TH1F>("CRE", "CR Efficiency", 50, 0, 1);
165  fCRP = tfs->make<TH1F>("CRP", "CR Purity", 50, 0, 1);
166  }
167  // Neutrino Int histos
168  fNuKE_elec = tfs->make<TH1F>("NuKE_elec", "NuKE electron", 100, 0, 4000);
169  fNuKE_muon = tfs->make<TH1F>("NuKE_muon", "NuKE muon", 100, 0, 4000);
170  fNuKE_pion = tfs->make<TH1F>("NuKE_pion", "NuKE pion", 100, 0, 4000);
171  fNuKE_kaon = tfs->make<TH1F>("NuKE_kaon", "NuKE kaon", 100, 0, 4000);
172  fNuKE_prot = tfs->make<TH1F>("NuKE_prot", "NuKE proton", 100, 0, 4000);
173 
174  fNuEP2_elec = tfs->make<TH1F>("NuEP2_elec", "NuEP2 electron", 50, 0, 1);
175  fNuEP2_muon = tfs->make<TH1F>("NuEP2_muon", "NuEP2 muon", 50, 0, 1);
176  fNuEP2_pion = tfs->make<TH1F>("NuEP2_pion", "NuEP2 pion", 50, 0, 1);
177  fNuEP2_kaon = tfs->make<TH1F>("NuEP2_kaon", "NuEP2 kaon", 50, 0, 1);
178  fNuEP2_prot = tfs->make<TH1F>("NuEP2_prot", "NuEP2 proton", 50, 0, 1);
179 
180  fNuE_elec = tfs->make<TH1F>("NuE_elec", "Nu Efficiency electron", 50, 0, 1);
181  fNuE_muon = tfs->make<TH1F>("NuE_muon", "Nu Efficiency muon", 50, 0, 1);
182  fNuE_pion = tfs->make<TH1F>("NuE_pion", "Nu Efficiency pion", 50, 0, 1);
183  fNuE_kaon = tfs->make<TH1F>("NuE_kaon", "Nu Efficiency kaon", 50, 0, 1);
184  fNuE_prot = tfs->make<TH1F>("NuE_prot", "Nu Efficiency proton", 50, 0, 1);
185 
186  fNuP_elec = tfs->make<TH1F>("NuP_elec", "Nu Purity electron", 50, 0, 1);
187  fNuP_muon = tfs->make<TH1F>("NuP_muon", "Nu Purity muon", 50, 0, 1);
188  fNuP_pion = tfs->make<TH1F>("NuP_pion", "Nu Purity pion", 50, 0, 1);
189  fNuP_kaon = tfs->make<TH1F>("NuP_kaon", "Nu Purity kaon", 50, 0, 1);
190  fNuP_prot = tfs->make<TH1F>("NuP_prot", "Nu Purity proton", 50, 0, 1);
191 
192  // True - Reco vertex difference
193  fNuVtx_dx = tfs->make<TH1F>("Vtx dx", "Vtx dx", 80, -10, 10);
194  fNuVtx_dy = tfs->make<TH1F>("Vtx dy", "Vtx dy", 80, -10, 10);
195  fNuVtx_dz = tfs->make<TH1F>("Vtx dz", "Vtx dz", 80, -10, 10);
196 
197  fNuEP2_KE_elec = tfs->make<TProfile>("NuEP2_KE_elec", "NuEP2 electron vs KE", 200, 0, 2000);
198  fNuEP2_KE_muon = tfs->make<TProfile>("NuEP2_KE_muon", "NuEP2 muon vs KE", 200, 0, 2000);
199  fNuEP2_KE_pion = tfs->make<TProfile>("NuEP2_KE_pion", "NuEP2 pion vs KE", 200, 0, 2000);
200  fNuEP2_KE_kaon = tfs->make<TProfile>("NuEP2_KE_kaon", "NuEP2 kaon vs KE", 200, 0, 2000);
201  fNuEP2_KE_prot = tfs->make<TProfile>("NuEP2_KE_prot", "NuEP2 proton vs KE", 200, 0, 2000);
202  }
tfs
art::ServiceHandle< art::TFileService > tfs
Definition: NuVertexChargeTree_module.cc:254
void cluster::ClusterAna::endJob ( )
overrideprivate

Definition at line 205 of file ClusterAna_module.cc.

206  {
207  if (fPrintLevel == -1) outFile.close();
208  }

Member Data Documentation

std::string cluster::ClusterAna::fClusterModuleLabel
private

Definition at line 93 of file ClusterAna_module.cc.

TH1F* cluster::ClusterAna::fCRE
private

Definition at line 57 of file ClusterAna_module.cc.

TH1F* cluster::ClusterAna::fCREP2
private

Definition at line 56 of file ClusterAna_module.cc.

TH1F* cluster::ClusterAna::fCRP
private

Definition at line 58 of file ClusterAna_module.cc.

std::vector<float> cluster::ClusterAna::fElecKERange
private

Definition at line 95 of file ClusterAna_module.cc.

std::string cluster::ClusterAna::fHitsModuleLabel
private

Definition at line 92 of file ClusterAna_module.cc.

std::vector<float> cluster::ClusterAna::fKaonKERange
private

Definition at line 98 of file ClusterAna_module.cc.

bool cluster::ClusterAna::fMergeDaughters
private

Definition at line 101 of file ClusterAna_module.cc.

std::vector<float> cluster::ClusterAna::fMuonKERange
private

Definition at line 96 of file ClusterAna_module.cc.

TH1F* cluster::ClusterAna::fNClusters
private

Definition at line 53 of file ClusterAna_module.cc.

TH1F* cluster::ClusterAna::fNHitInCluster
private

Definition at line 54 of file ClusterAna_module.cc.

TH1F* cluster::ClusterAna::fNuE_elec
private

Definition at line 71 of file ClusterAna_module.cc.

TH1F* cluster::ClusterAna::fNuE_kaon
private

Definition at line 74 of file ClusterAna_module.cc.

TH1F* cluster::ClusterAna::fNuE_muon
private

Definition at line 72 of file ClusterAna_module.cc.

TH1F* cluster::ClusterAna::fNuE_pion
private

Definition at line 73 of file ClusterAna_module.cc.

TH1F* cluster::ClusterAna::fNuE_prot
private

Definition at line 75 of file ClusterAna_module.cc.

TH1F* cluster::ClusterAna::fNuEP2_elec
private

Definition at line 66 of file ClusterAna_module.cc.

TH1F* cluster::ClusterAna::fNuEP2_kaon
private

Definition at line 69 of file ClusterAna_module.cc.

TProfile* cluster::ClusterAna::fNuEP2_KE_elec
private

Definition at line 86 of file ClusterAna_module.cc.

TProfile* cluster::ClusterAna::fNuEP2_KE_kaon
private

Definition at line 89 of file ClusterAna_module.cc.

TProfile* cluster::ClusterAna::fNuEP2_KE_muon
private

Definition at line 87 of file ClusterAna_module.cc.

TProfile* cluster::ClusterAna::fNuEP2_KE_pion
private

Definition at line 88 of file ClusterAna_module.cc.

TProfile* cluster::ClusterAna::fNuEP2_KE_prot
private

Definition at line 90 of file ClusterAna_module.cc.

TH1F* cluster::ClusterAna::fNuEP2_muon
private

Definition at line 67 of file ClusterAna_module.cc.

TH1F* cluster::ClusterAna::fNuEP2_pion
private

Definition at line 68 of file ClusterAna_module.cc.

TH1F* cluster::ClusterAna::fNuEP2_prot
private

Definition at line 70 of file ClusterAna_module.cc.

TH1F* cluster::ClusterAna::fNuKE_elec
private

Definition at line 60 of file ClusterAna_module.cc.

TH1F* cluster::ClusterAna::fNuKE_kaon
private

Definition at line 63 of file ClusterAna_module.cc.

TH1F* cluster::ClusterAna::fNuKE_muon
private

Definition at line 61 of file ClusterAna_module.cc.

TH1F* cluster::ClusterAna::fNuKE_pion
private

Definition at line 62 of file ClusterAna_module.cc.

TH1F* cluster::ClusterAna::fNuKE_prot
private

Definition at line 64 of file ClusterAna_module.cc.

TH1F* cluster::ClusterAna::fNuP_elec
private

Definition at line 76 of file ClusterAna_module.cc.

TH1F* cluster::ClusterAna::fNuP_kaon
private

Definition at line 79 of file ClusterAna_module.cc.

TH1F* cluster::ClusterAna::fNuP_muon
private

Definition at line 77 of file ClusterAna_module.cc.

TH1F* cluster::ClusterAna::fNuP_pion
private

Definition at line 78 of file ClusterAna_module.cc.

TH1F* cluster::ClusterAna::fNuP_prot
private

Definition at line 80 of file ClusterAna_module.cc.

TH1F* cluster::ClusterAna::fNuVtx_dx
private

Definition at line 82 of file ClusterAna_module.cc.

TH1F* cluster::ClusterAna::fNuVtx_dy
private

Definition at line 83 of file ClusterAna_module.cc.

TH1F* cluster::ClusterAna::fNuVtx_dz
private

Definition at line 84 of file ClusterAna_module.cc.

std::vector<float> cluster::ClusterAna::fPionKERange
private

Definition at line 97 of file ClusterAna_module.cc.

short cluster::ClusterAna::fPrintLevel
private

Definition at line 104 of file ClusterAna_module.cc.

std::vector<float> cluster::ClusterAna::fProtKERange
private

Definition at line 99 of file ClusterAna_module.cc.

bool cluster::ClusterAna::fSkipCosmics
private

Definition at line 103 of file ClusterAna_module.cc.

short cluster::ClusterAna::fTrackWeightOption
private

Definition at line 100 of file ClusterAna_module.cc.

std::string cluster::ClusterAna::fVertexModuleLabel
private

Definition at line 94 of file ClusterAna_module.cc.

short cluster::ClusterAna::moduleID
private

Definition at line 105 of file ClusterAna_module.cc.

std::ofstream cluster::ClusterAna::outFile
private

Definition at line 107 of file ClusterAna_module.cc.

The documentation for this class was generated from the following file:
Generated by   doxygen 1.8.5