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
icarus::AnalysisTree Class Reference

Creates a simple ROOT tree with tracking and calorimetry information. More...

Inheritance diagram for icarus::AnalysisTree:

Public Member Functions

 AnalysisTree (fhicl::ParameterSet const &pset)
 
virtual ~AnalysisTree ()
 
void analyze (const art::Event &evt)
 read access to event More...
 
void beginSubRun (const art::SubRun &sr)
 

Private Member Functions

void HitsPurity (detinfo::DetectorClocksData const &clockData, std::vector< art::Ptr< recob::Hit > > const &hits, Int_t &trackid, Float_t &purity, double &maxe)
 
double length (const recob::Track &track)
 
double length (const simb::MCParticle &part, TVector3 &start, TVector3 &end)
 
double bdist (const recob::tracking::Point_t &pos)
 
size_t GetNTrackers () const
 Returns the number of trackers configured. More...
 
void CreateData (bool bClearData=false)
 Creates the structure for the tree data; optionally initializes it. More...
 
void SetAddresses ()
 Sets the addresses of all the tree branches, creating the missing ones. More...
 
void SetTrackerAddresses (size_t iTracker)
 
void CreateTree (bool bClearData=false)
 Create the output tree and the data structures, if needed. More...
 
void DestroyData ()
 Destroy the local buffers (existing branches will point to invalid address!) More...
 
void CheckData (std::string caller) const
 Helper function: throws if no data structure is available. More...
 
void CheckTree (std::string caller) const
 Helper function: throws if no tree is available. More...
 

Private Attributes

TTree * fTree
 
AnalysisTreeDataStruct * fData
 
AnalysisTreeDataStruct::SubRunData_t SubRunData
 
std::string fDigitModuleLabel
 
std::string fHitsModuleLabel
 
std::string fLArG4ModuleLabel
 
std::string fCalDataModuleLabel
 
std::string fGenieGenModuleLabel
 
std::string fCryGenModuleLabel
 
std::string fG4ModuleLabel
 
std::string fVertexModuleLabel
 
std::vector< std::string > fTrackModuleLabel
 
std::vector< std::string > fCalorimetryModuleLabel
 
std::vector< std::string > fParticleIDModuleLabel
 
std::string fPOTModuleLabel
 
bool fUseBuffer
 whether to use a permanent buffer (faster, huge memory) More...
 
bool fSaveAuxDetInfo
 whether to extract and save auxiliary detector data More...
 
bool fSaveCryInfo
 
bool fSaveGenieInfo
 whether to extract and save CRY particle data More...
 
bool fSaveGeantInfo
 whether to extract and save Genie information More...
 
bool fSaveHitInfo
 whether to extract and save Geant information More...
 
bool fSaveTrackInfo
 whether to extract and save Hit information More...
 
bool fSaveVertexInfo
 whether to extract and save Track information More...
 
std::vector< std::string > fCosmicTaggerAssocLabel
 whether to extract and save Vertex information More...
 
std::vector< std::string > fFlashMatchAssocLabel
 
bool isCosmics
 if it contains cosmics More...
 
bool fSaveCaloCosmics
 save calorimetry information for cosmics More...
 
float fG4minE
 

Detailed Description

Creates a simple ROOT tree with tracking and calorimetry information.

Configuration parameters

Definition at line 670 of file icaruscode/icaruscode/Analysis/AnalysisTree_module.cc.

Constructor & Destructor Documentation

icarus::AnalysisTree::AnalysisTree ( fhicl::ParameterSet const &  pset)
explicit

Definition at line 1661 of file icaruscode/icaruscode/Analysis/AnalysisTree_module.cc.

1661  :
1662  EDAnalyzer(pset),
1663  fTree(nullptr), fData(nullptr),
1664  fDigitModuleLabel (pset.get< std::string >("DigitModuleLabel") ),
1665  fHitsModuleLabel (pset.get< std::string >("HitsModuleLabel") ),
1666  fLArG4ModuleLabel (pset.get< std::string >("LArGeantModuleLabel") ),
1667  fCalDataModuleLabel (pset.get< std::string >("CalDataModuleLabel") ),
1668  fGenieGenModuleLabel (pset.get< std::string >("GenieGenModuleLabel") ),
1669  fCryGenModuleLabel (pset.get< std::string >("CryGenModuleLabel") ),
1670  fG4ModuleLabel (pset.get< std::string >("G4ModuleLabel") ),
1671  fVertexModuleLabel (pset.get< std::string> ("VertexModuleLabel") ),
1672  fTrackModuleLabel (pset.get< std::vector<std::string> >("TrackModuleLabel")),
1673  fCalorimetryModuleLabel (pset.get< std::vector<std::string> >("CalorimetryModuleLabel")),
1674  fParticleIDModuleLabel (pset.get< std::vector<std::string> >("ParticleIDModuleLabel") ),
1675  fPOTModuleLabel (pset.get< std::string >("POTModuleLabel") ),
1676  fUseBuffer (pset.get< bool >("UseBuffers", false)),
1677  fSaveAuxDetInfo (pset.get< bool >("SaveAuxDetInfo", false)),
1678  fSaveCryInfo (pset.get< bool >("SaveCryInfo", false)),
1679  fSaveGenieInfo (pset.get< bool >("SaveGenieInfo", false)),
1680  fSaveGeantInfo (pset.get< bool >("SaveGeantInfo", false)),
1681  fSaveHitInfo (pset.get< bool >("SaveHitInfo", false)),
1682  fSaveTrackInfo (pset.get< bool >("SaveTrackInfo", false)),
1683  fSaveVertexInfo (pset.get< bool >("SaveVertexInfo", false)),
1684  //fCosmicTaggerAssocLabel (pset.get<std::vector< std::string > >("CosmicTaggerAssocLabel") ),
1685  //fFlashMatchAssocLabel (pset.get<std::vector< std::string > >("FlashMatchAssocLabel") ),
1686  isCosmics(false),
1687  fSaveCaloCosmics (pset.get< bool >("SaveCaloCosmics",false)),
1688  fG4minE (pset.get< float>("G4minE",0.01))
1689 {
1690  if (fSaveAuxDetInfo == true) fSaveGeantInfo = true;
1691  mf::LogInfo("AnalysisTree") << "Configuration:"
1692  << "\n UseBuffers: " << std::boolalpha << fUseBuffer
1693  ;
1694  if (GetNTrackers() > kMaxTrackers) {
1695  throw art::Exception(art::errors::Configuration)
1696  << "AnalysisTree currently supports only up to " << kMaxTrackers
1697  << " tracking algorithms, but " << GetNTrackers() << " are specified."
1698  << "\nYou can increase kMaxTrackers and recompile.";
1699  } // if too many trackers
1700  if (fTrackModuleLabel.size() != fCalorimetryModuleLabel.size()){
1701  throw art::Exception(art::errors::Configuration)
1702  << "fTrackModuleLabel.size() = "<<fTrackModuleLabel.size()<<" does not match "
1703  << "fCalorimetryModuleLabel.size() = "<<fCalorimetryModuleLabel.size();
1704  }
1705  if (fTrackModuleLabel.size() != fParticleIDModuleLabel.size()){
1706  throw art::Exception(art::errors::Configuration)
1707  << "fTrackModuleLabel.size() = "<<fTrackModuleLabel.size()<<" does not match "
1708  << "fParticleIDModuleLabel.size() = "<<fParticleIDModuleLabel.size();
1709  }
1710 } // icarus::AnalysisTree::AnalysisTree()
icarus::AnalysisTree::fSaveTrackInfo
bool fSaveTrackInfo
whether to extract and save Hit information
Definition: icaruscode/icaruscode/Analysis/AnalysisTree_module.cc:716
kMaxTrackers
constexpr int kMaxTrackers
Definition: icaruscode/icaruscode/Analysis/AnalysisTree_module.cc:143
icarus::AnalysisTree::fSaveAuxDetInfo
bool fSaveAuxDetInfo
whether to extract and save auxiliary detector data
Definition: icaruscode/icaruscode/Analysis/AnalysisTree_module.cc:711
icarus::AnalysisTree::fUseBuffer
bool fUseBuffer
whether to use a permanent buffer (faster, huge memory)
Definition: icaruscode/icaruscode/Analysis/AnalysisTree_module.cc:710
icarus::AnalysisTree::fSaveHitInfo
bool fSaveHitInfo
whether to extract and save Geant information
Definition: icaruscode/icaruscode/Analysis/AnalysisTree_module.cc:715
icarus::AnalysisTree::GetNTrackers
size_t GetNTrackers() const
Returns the number of trackers configured.
Definition: icaruscode/icaruscode/Analysis/AnalysisTree_module.cc:726
icarus::AnalysisTree::fSaveCaloCosmics
bool fSaveCaloCosmics
save calorimetry information for cosmics
Definition: icaruscode/icaruscode/Analysis/AnalysisTree_module.cc:723
icarus::AnalysisTree::fSaveGenieInfo
bool fSaveGenieInfo
whether to extract and save CRY particle data
Definition: icaruscode/icaruscode/Analysis/AnalysisTree_module.cc:713
icarus::AnalysisTree::fSaveGeantInfo
bool fSaveGeantInfo
whether to extract and save Genie information
Definition: icaruscode/icaruscode/Analysis/AnalysisTree_module.cc:714
icarus::AnalysisTree::fSaveVertexInfo
bool fSaveVertexInfo
whether to extract and save Track information
Definition: icaruscode/icaruscode/Analysis/AnalysisTree_module.cc:717
icarus::AnalysisTree::~AnalysisTree ( )
virtual

Definition at line 1713 of file icaruscode/icaruscode/Analysis/AnalysisTree_module.cc.

1714 {
1715  DestroyData();
1716 }
icarus::AnalysisTree::DestroyData
void DestroyData()
Destroy the local buffers (existing branches will point to invalid address!)
Definition: icaruscode/icaruscode/Analysis/AnalysisTree_module.cc:772

Member Function Documentation

void icarus::AnalysisTree::analyze ( const art::Event &  evt)

read access to event

transfer the run and subrun data to the tree data object

Definition at line 1741 of file icaruscode/icaruscode/Analysis/AnalysisTree_module.cc.

1742 {
1743  //services
1744  art::ServiceHandle<geo::Geometry> geom;
1745  art::ServiceHandle<cheat::BackTrackerService> backTracker;
1746  art::ServiceHandle<cheat::ParticleInventoryService> partInventory;
1747  // auto const* detprop = lar::providerFrom<detinfo::DetectorPropertiesService>();
1748  // auto const* LArProp = lar::providerFrom<detinfo::LArPropertiesService>();
1749 
1750  // collect the sizes which might me needed to resize the tree data structure:
1751  bool isMC = !evt.isRealData();
1752 
1753  // * hits
1754  art::Handle< std::vector<recob::Hit> > hitListHandle;
1755  std::vector<art::Ptr<recob::Hit> > hitlist;
1756  if (evt.getByLabel(fHitsModuleLabel,hitListHandle))
1757  art::fill_ptr_vector(hitlist, hitListHandle);
1758 
1759  // * vertices
1760  art::Handle< std::vector<recob::Vertex> > vtxListHandle;
1761  std::vector<art::Ptr<recob::Vertex> > vtxlist;
1762  if (evt.getByLabel(fVertexModuleLabel,vtxListHandle))
1763  art::fill_ptr_vector(vtxlist, vtxListHandle);
1764 
1765 
1766  // * MC truth information
1767  art::Handle< std::vector<simb::MCTruth> > mctruthListHandle;
1768  std::vector<art::Ptr<simb::MCTruth> > mclist;
1769  if (evt.getByLabel(fGenieGenModuleLabel,mctruthListHandle))
1770  art::fill_ptr_vector(mclist, mctruthListHandle);
1771 
1772  // *MC truth cosmic generator information
1773  art::Handle< std::vector<simb::MCTruth> > mctruthcryListHandle;
1774  std::vector<art::Ptr<simb::MCTruth> > mclistcry;
1775  if (fSaveCryInfo){
1776  if (evt.getByLabel(fCryGenModuleLabel,mctruthcryListHandle))
1777  art::fill_ptr_vector(mclistcry, mctruthcryListHandle);
1778  }
1779 
1780  art::Ptr<simb::MCTruth> mctruthcry;
1781  int nCryPrimaries = 0;
1782 
1783  if (fSaveCryInfo){
1784  mctruthcry = mclistcry[0];
1785  nCryPrimaries = mctruthcry->NParticles();
1786  }
1787 
1788  int nGeniePrimaries = 0, nGEANTparticles = 0, nMCNeutrinos = 0;
1789 
1790  auto const clockData = art::ServiceHandle<detinfo::DetectorClocksService const>()->DataFor(evt);
1791  art::Ptr<simb::MCTruth> mctruth;
1792  //Brailsford 2017/10/16
1793  //Fix for issue 17917
1794  //With the code change to accept multiple neutrinos per TTree::Entry into the TTree, this int is no longer needed (it makes compilation fail due to a warning)
1795  //int imc = 0;
1796  if (isMC) { //is MC
1797  // GENIE
1798  if (!mclist.empty()){//at least one mc record
1799  //if (fSaveGenieInfo){
1800  //if (mclist[0]->NeutrinoSet()){//is neutrino
1801  //sometimes there can be multiple neutrino interactions in one spill
1802  //this is trying to find the primary interaction
1803  //by looking for the highest energy deposition
1804  //std::map<art::Ptr<simb::MCTruth>,double> mctruthemap;
1805  static bool isfirsttime = true;
1806  if (isfirsttime){
1807  for (size_t i = 0; i<hitlist.size(); i++){
1808  //if (hitlist[i]->View() == geo::kV){//collection view
1809  std::vector<sim::TrackIDE> eveIDs = backTracker->HitToEveTrackIDEs(clockData, hitlist[i]);
1810  for (size_t e = 0; e<eveIDs.size(); e++){
1811  art::Ptr<simb::MCTruth> ev_mctruth = partInventory->TrackIdToMCTruth_P(eveIDs[e].trackID);
1812  //mctruthemap[ev_mctruth]+=eveIDs[e].energy;
1813  if (ev_mctruth->Origin() == simb::kCosmicRay) isCosmics = true;
1814  }
1815  //}
1816  }
1817  isfirsttime = false;
1818  if (fSaveCaloCosmics) isCosmics = false; //override to save calo info
1819  }
1820 
1821 // double maxenergy = -1;
1822 // int imc0 = 0;
1823 // for (std::map<art::Ptr<simb::MCTruth>,double>::iterator ii=mctruthemap.begin(); ii!=mctruthemap.end(); ++ii){
1824 // if ((ii->second)>maxenergy){
1825 // maxenergy = ii->second;
1826 // mctruth = ii->first;
1827 // imc = imc0;
1828 // }
1829 // imc0++;
1830 // }
1831 
1832  //imc = 0; //set imc to 0 to solve a confusion for BNB+cosmic files where there are two MCTruth
1833  mctruth = mclist[0];
1834 
1835  if (mctruth->NeutrinoSet()) nGeniePrimaries = mctruth->NParticles();
1836  //} //end (fSaveGenieInfo)
1837 
1838  const sim::ParticleList& plist = partInventory->ParticleList();
1839  nGEANTparticles = plist.size();
1840 
1841  // to know the number of particles in AV would require
1842  // looking at all of them; so we waste some memory here
1843  } // if have MC truth
1844  MF_LOG_DEBUG("AnalysisTree") << "Expected "
1845  << nGEANTparticles << " GEANT particles, "
1846  << nGeniePrimaries << " GENIE particles";
1847  } // if MC
1848 
1849  //Brailsford 2017/10/16
1850  //Initially call the number of neutrinos to be stored the number of MCTruth objects. This is not strictly true i.e. BNB + cosmic overlay but we will count the number of neutrinos later
1851  nMCNeutrinos = mclist.size();
1852 
1853  CreateData(); // tracker data is created with default constructor
1854  if (fSaveGenieInfo){
1855  fData->ResizeGenie(nGeniePrimaries);
1856  fData->ResizeMCNeutrino(nMCNeutrinos);
1857  }
1858  if (fSaveCryInfo)
1859  fData->ResizeCry(nCryPrimaries);
1860  if (fSaveGeantInfo)
1861  fData->ResizeGEANT(nGEANTparticles);
1862  fData->ClearLocalData(); // don't bother clearing tracker data yet
1863 
1864 // const size_t Nplanes = 3; // number of wire planes; pretty much constant...
1865  const size_t NTrackers = GetNTrackers(); // number of trackers passed into fTrackModuleLabel
1866  const size_t NHits = hitlist.size(); // number of hits
1867  const size_t NVertices = vtxlist.size(); // number of vertices
1868  // make sure there is the data, the tree and everything;
1869  CreateTree();
1870 
1871  /// transfer the run and subrun data to the tree data object
1872 // fData->RunData = RunData;
1873  fData->SubRunData = SubRunData;
1874 
1875  fData->isdata = int(!isMC);
1876 
1877  std::vector< art::Handle< std::vector<recob::Track> > > trackListHandle(NTrackers);
1878  std::vector< std::vector<art::Ptr<recob::Track> > > tracklist(NTrackers);
1879  for (unsigned int it = 0; it < NTrackers; ++it){
1880  if (evt.getByLabel(fTrackModuleLabel[it],trackListHandle[it]))
1881  art::fill_ptr_vector(tracklist[it], trackListHandle[it]);
1882  }
1883 
1884  art::Handle< std::vector<simb::MCFlux> > mcfluxListHandle;
1885  std::vector<art::Ptr<simb::MCFlux> > fluxlist;
1886  if (evt.getByLabel(fGenieGenModuleLabel,mcfluxListHandle))
1887  art::fill_ptr_vector(fluxlist, mcfluxListHandle);
1888 
1889  std::vector<const sim::AuxDetSimChannel*> fAuxDetSimChannels;
1890  if (fSaveAuxDetInfo && fSaveGeantInfo){
1891  evt.getView(fLArG4ModuleLabel, fAuxDetSimChannels);
1892  }
1893 
1894  std::vector<const sim::SimChannel*> fSimChannels;
1895  if (isMC && fSaveGeantInfo){
1896  evt.getView(fLArG4ModuleLabel, fSimChannels);
1897  }
1898 
1899  fData->run = evt.run();
1900  fData->subrun = evt.subRun();
1901  fData->event = evt.id().event();
1902 
1903  art::Timestamp ts = evt.time();
1904  TTimeStamp tts(ts.timeHigh(), ts.timeLow());
1905  fData->evttime = tts.AsDouble();
1906 
1907  //copied from MergeDataPaddles.cxx
1908  art::Handle< raw::BeamInfo > beam;
1909  if (evt.getByLabel("beam",beam)){
1910  fData->beamtime = (double)beam->get_t_ms();
1911  fData->beamtime/=1000.; //in second
1912  }
1913 
1914 // std::cout<<detprop->NumberTimeSamples()<<" "<<detprop->ReadOutWindowSize()<<std::endl;
1915 // std::cout<<geom->DetHalfHeight()*2<<" "<<geom->DetHalfWidth()*2<<" "<<geom->DetLength()<<std::endl;
1916 // std::cout<<geom->Nwires(0)<<" "<<geom->Nwires(1)<<" "<<geom->Nwires(2)<<std::endl;
1917 
1918  //hit information
1919  if (fSaveHitInfo){
1920  fData->no_hits = (int) NHits;
1921  if (NHits > kMaxHits) {
1922  // got this error? consider increasing kMaxHits
1923  // (or ask for a redesign using vectors)
1924  mf::LogError("AnalysisTree:limits") << "event has " << NHits
1925  << " hits, only kMaxHits=" << kMaxHits << " stored in tree";
1926  }
1927  for (size_t i = 0; i < NHits && i < kMaxHits ; ++i){//loop over hits
1928  fData->hit_channel[i] = hitlist[i]->Channel();
1929  fData->hit_tpc[i] = hitlist[i]->WireID().TPC;
1930  fData->hit_plane[i] = hitlist[i]->WireID().Plane;
1931  fData->hit_wire[i] = hitlist[i]->WireID().Wire;
1932  fData->hit_peakT[i] = hitlist[i]->PeakTime();
1933  fData->hit_charge[i] = hitlist[i]->Integral();
1934  fData->hit_ph[i] = hitlist[i]->PeakAmplitude();
1935  fData->hit_startT[i] = hitlist[i]->PeakTimeMinusRMS();
1936  fData->hit_endT[i] = hitlist[i]->PeakTimePlusRMS();
1937  /*
1938  for (unsigned int it=0; it<fTrackModuleLabel.size();++it){
1939  art::FindManyP<recob::Track> fmtk(hitListHandle,evt,fTrackModuleLabel[it]);
1940  if (fmtk.at(i).size()!=0){
1941  hit_trkid[it][i] = fmtk.at(i)[0]->ID();
1942  }
1943  else
1944  hit_trkid[it][i] = 0;
1945  }
1946  */
1947 
1948  if (!evt.isRealData()){
1949  fData -> hit_nelec[i] = 0;
1950  fData -> hit_energy[i] = 0;
1951  const sim::SimChannel* chan = 0;
1952  for(size_t sc = 0; sc < fSimChannels.size(); ++sc){
1953  if(fSimChannels[sc]->Channel() == hitlist[i]->Channel()) chan = fSimChannels[sc];
1954  }
1955  if (chan){
1956  for (auto const& tdcide: chan->TDCIDEMap()) {
1957  // loop over the vector of IDE objects.
1958  std::vector<sim::IDE> const& idevec = tdcide.second;
1959  for(auto const& ide: idevec) {
1960  fData -> hit_nelec[i] += ide.numElectrons;
1961  fData -> hit_energy[i] += ide.energy;
1962  } // all IDEs in the TDC tick
1963  } // all TDC ticks
1964  }
1965  }
1966  }
1967 
1968  if (evt.getByLabel(fHitsModuleLabel,hitListHandle)){
1969  //Find tracks associated with hits
1970  art::FindManyP<recob::Track> fmtk(hitListHandle,evt,fTrackModuleLabel[0]);
1971  for (size_t i = 0; i < NHits && i < kMaxHits ; ++i){//loop over hits
1972  if (fmtk.isValid()){
1973  if (fmtk.at(i).size()!=0){
1974  fData->hit_trkid[i] = fmtk.at(i)[0]->ID();
1975  }
1976  else
1977  fData->hit_trkid[i] = -1;
1978  }
1979  }
1980  }
1981  }// end (fSaveHitInfo)
1982 
1983  //vertex information
1984  if (fSaveVertexInfo){
1985  fData->nvtx = NVertices;
1986  if (NVertices > kMaxVertices){
1987  // got this error? consider increasing kMaxVerticestra
1988  // (or ask for a redesign using vectors)
1989  mf::LogError("AnalysisTree:limits") << "event has " << NVertices
1990  << " vertices, only kMaxVertices=" << kMaxVertices << " stored in tree";
1991  }
1992  for (size_t i = 0; i < NVertices && i < kMaxVertices ; ++i){//loop over hits
1993  Double_t xyz[3] = {};
1994  vtxlist[i]->XYZ(xyz);
1995  for (size_t j = 0; j<3; ++j) fData->vtx[i][j] = xyz[j];
1996  }
1997  }// end (fSaveVertexInfo)
1998 
1999  //track information for multiple trackers
2000  if (fSaveTrackInfo){
2001  for (unsigned int iTracker=0; iTracker < NTrackers; ++iTracker){
2002  AnalysisTreeDataStruct::TrackDataStruct& TrackerData = fData->GetTrackerData(iTracker);
2003 
2004  size_t NTracks = tracklist[iTracker].size();
2005  // allocate enough space for this number of tracks (but at least for one of them!)
2006  TrackerData.SetMaxTracks(std::max(NTracks, (size_t) 1));
2007  TrackerData.Clear(); // clear all the data
2008 
2009  TrackerData.ntracks = (int) NTracks;
2010 
2011  // now set the tree addresses to the newly allocated memory;
2012  // this creates the tree branches in case they are not there yet
2013  SetTrackerAddresses(iTracker);
2014  if (NTracks > TrackerData.GetMaxTracks()) {
2015  // got this error? it might be a bug,
2016  // since we are supposed to have allocated enough space to fit all tracks
2017  mf::LogError("AnalysisTree:limits") << "event has " << NTracks
2018  << " " << fTrackModuleLabel[iTracker] << " tracks, only "
2019  << TrackerData.GetMaxTracks() << " stored in tree";
2020  }
2021 
2022  //call the track momentum algorithm that gives you momentum based on track range
2023  //trkf::TrackMomentumCalculator trkm;
2024 
2025  for(size_t iTrk=0; iTrk < NTracks; ++iTrk){//loop over tracks
2026 
2027  //Cosmic Tagger information
2028  if (fCosmicTaggerAssocLabel.size() > iTracker) {
2029  art::FindManyP<anab::CosmicTag> fmct(trackListHandle[iTracker],evt,fCosmicTaggerAssocLabel[iTracker]);
2030  if (fmct.isValid()){
2031  TrackerData.trkncosmictags_tagger[iTrk] = fmct.at(iTrk).size();
2032  if (fmct.at(iTrk).size()>0){
2033  if(fmct.at(iTrk).size()>1)
2034  std::cerr << "\n Warning : more than one cosmic tag per track in module! assigning the first tag to the track" << fCosmicTaggerAssocLabel[iTracker];
2035  TrackerData.trkcosmicscore_tagger[iTrk] = fmct.at(iTrk).at(0)->CosmicScore();
2036  TrackerData.trkcosmictype_tagger[iTrk] = fmct.at(iTrk).at(0)->CosmicType();
2037  }
2038  }
2039  } // if we have matching fCosmicTaggerAssocLabel
2040 
2041  //Flash match compatibility information
2042  if (fFlashMatchAssocLabel.size() > iTracker) {
2043  //Unlike CosmicTagger, Flash match doesn't assign a cosmic tag for every track. For those tracks, AnalysisTree initializes them with -9999 or -99999
2044  art::FindManyP<anab::CosmicTag> fmbfm(trackListHandle[iTracker],evt,fFlashMatchAssocLabel[iTracker]);
2045  if (fmbfm.isValid()){
2046  TrackerData.trkncosmictags_flashmatch[iTrk] = fmbfm.at(iTrk).size();
2047  if (fmbfm.at(iTrk).size()>0){
2048  if(fmbfm.at(iTrk).size()>1)
2049  std::cerr << "\n Warning : more than one cosmic tag per track in module! assigning the first tag to the track" << fFlashMatchAssocLabel[iTracker];
2050  TrackerData.trkcosmicscore_flashmatch[iTrk] = fmbfm.at(iTrk).at(0)->CosmicScore();
2051  TrackerData.trkcosmictype_flashmatch[iTrk] = fmbfm.at(iTrk).at(0)->CosmicType();
2052  //std::cout<<"\n"<<evt.event()<<"\t"<<iTrk<<"\t"<<fmbfm.at(iTrk).at(0)->CosmicScore()<<"\t"<<fmbfm.at(iTrk).at(0)->CosmicType();
2053  }
2054  }
2055  } // if we have matching fFlashMatchAssocLabel
2056 
2057  art::Ptr<recob::Track> ptrack(trackListHandle[iTracker], iTrk);
2058  const recob::Track& track = *ptrack;
2059 
2060  //TVector3 pos, dir_start, dir_end, end;
2061  recob::tracking::Point_t pos,end;
2062  recob::tracking::Vector_t dir_start,dir_end;
2063 
2064  double tlen = 0.; //mom = 0.;
2065  int TrackID = -1;
2066 
2067  int ntraj = track.NumberTrajectoryPoints();
2068  if (ntraj > 0) {
2069  pos = track.Vertex();
2070  dir_start = track.VertexDirection();
2071  dir_end = track.EndDirection();
2072  end = track.End();
2073 
2074  tlen = length(track);
2075  TrackID = track.ID();
2076 
2077  double theta_xz = std::atan2(dir_start.X(), dir_start.Z());
2078  double theta_yz = std::atan2(dir_start.Y(), dir_start.Z());
2079  double dpos = bdist(pos);
2080  double dend = bdist(end);
2081 
2082  TrackerData.trkId[iTrk] = TrackID;
2083  TrackerData.trkstartx[iTrk] = pos.X();
2084  TrackerData.trkstarty[iTrk] = pos.Y();
2085  TrackerData.trkstartz[iTrk] = pos.Z();
2086  TrackerData.trkstartd[iTrk] = dpos;
2087  TrackerData.trkendx[iTrk] = end.X();
2088  TrackerData.trkendy[iTrk] = end.Y();
2089  TrackerData.trkendz[iTrk] = end.Z();
2090  TrackerData.trkendd[iTrk] = dend;
2091  TrackerData.trktheta[iTrk] = dir_start.Theta();
2092  TrackerData.trkphi[iTrk] = dir_start.Phi();
2093  TrackerData.trkstartdcosx[iTrk] = dir_start.X();
2094  TrackerData.trkstartdcosy[iTrk] = dir_start.Y();
2095  TrackerData.trkstartdcosz[iTrk] = dir_start.Z();
2096  TrackerData.trkenddcosx[iTrk] = dir_end.X();
2097  TrackerData.trkenddcosy[iTrk] = dir_end.Y();
2098  TrackerData.trkenddcosz[iTrk] = dir_end.Z();
2099  TrackerData.trkthetaxz[iTrk] = theta_xz;
2100  TrackerData.trkthetayz[iTrk] = theta_yz;
2101  //TrackerData.trkmom[iTrk] = mom;
2102  TrackerData.trklen[iTrk] = tlen;
2103  //TrackerData.trkmomrange[iTrk] = trkm.GetTrackMomentum(tlen,13);
2104  //TrackerData.trkmommschi2[iTrk] = trkm.GetMomentumMultiScatterChi2(ptrack);
2105  //TrackerData.trkmommsllhd[iTrk] = trkm.GetMomentumMultiScatterLLHD(ptrack);
2106  } // if we have trajectory
2107 
2108  // find vertices associated with this track
2109  /*
2110  art::FindMany<recob::Vertex> fmvtx(trackListHandle[iTracker], evt, fVertexModuleLabel[iTracker]);
2111  if(fmvtx.isValid()) {
2112  std::vector<const recob::Vertex*> verts = fmvtx.at(iTrk);
2113  // should have two at most
2114  for(size_t ivx = 0; ivx < verts.size(); ++ivx) {
2115  verts[ivx]->XYZ(xyz);
2116  // find the vertex in TrackerData to get the index
2117  short theVtx = -1;
2118  for(short jvx = 0; jvx < TrackerData.nvtx; ++jvx) {
2119  if(TrackerData.vtx[jvx][2] == xyz[2]) {
2120  theVtx = jvx;
2121  break;
2122  }
2123  } // jvx
2124  // decide if it should be assigned to the track Start or End.
2125  // A simple dz test should suffice
2126  if(fabs(xyz[2] - TrackerData.trkstartz[iTrk]) <
2127  fabs(xyz[2] - TrackerData.trkendz[iTrk])) {
2128  TrackerData.trksvtxid[iTrk] = theVtx;
2129  } else {
2130  TrackerData.trkevtxid[iTrk] = theVtx;
2131  }
2132  } // vertices
2133  } // fmvtx.isValid()
2134  */
2135  Float_t minsdist = 10000;
2136  Float_t minedist = 10000;
2137  for (int ivx = 0; ivx < fData->nvtx && ivx < kMaxVertices; ++ivx){
2138  Float_t sdist = sqrt(pow(TrackerData.trkstartx[iTrk]-fData->vtx[ivx][0],2)+
2139  pow(TrackerData.trkstarty[iTrk]-fData->vtx[ivx][1],2)+
2140  pow(TrackerData.trkstartz[iTrk]-fData->vtx[ivx][2],2));
2141  Float_t edist = sqrt(pow(TrackerData.trkendx[iTrk]-fData->vtx[ivx][0],2)+
2142  pow(TrackerData.trkendy[iTrk]-fData->vtx[ivx][1],2)+
2143  pow(TrackerData.trkendz[iTrk]-fData->vtx[ivx][2],2));
2144  if (sdist<minsdist){
2145  minsdist = sdist;
2146  if (minsdist<10) TrackerData.trksvtxid[iTrk] = ivx;
2147  }
2148  if (edist<minedist){
2149  minedist = edist;
2150  if (minedist<10) TrackerData.trkevtxid[iTrk] = ivx;
2151  }
2152  }
2153  // find particle ID info
2154  art::FindMany<anab::ParticleID> fmpid(trackListHandle[iTracker], evt, fParticleIDModuleLabel[iTracker]);
2155  if(fmpid.isValid()) {
2156  std::vector<const anab::ParticleID*> pids = fmpid.at(iTrk);
2157  if (pids.size() == 0){
2158  mf::LogError("AnalysisTree:limits")
2159  << "No track-PID association found for " << fTrackModuleLabel[iTracker]
2160  << " track " << iTrk << ". Not saving particleID information.";
2161  }
2162  // Set dummy values
2163  double pidpdg[3] = {0,0,0};
2164  double pidchi[3] = {0.,0.,0.};
2165  for(size_t planenum=0; planenum<3; ++planenum) {
2166  TrackerData.trkpidchimu[iTrk][planenum] = 0.;
2167  TrackerData.trkpidchipi[iTrk][planenum] = 0.;
2168  TrackerData.trkpidchika[iTrk][planenum] = 0.;
2169  TrackerData.trkpidchipr[iTrk][planenum] = 0.;
2170  TrackerData.trkpidpida[iTrk][planenum] = 0.;
2171  }
2172  for (size_t ipid=0; ipid<pids.size(); ipid++){
2173  std::vector<anab::sParticleIDAlgScores> AlgScoresVec = pids[ipid]->ParticleIDAlgScores();
2174 
2175  // Loop though AlgScoresVec and find the variables we want
2176  for (size_t i_algscore=0; i_algscore<AlgScoresVec.size(); i_algscore++){
2177  anab::sParticleIDAlgScores AlgScore = AlgScoresVec.at(i_algscore);
2178 
2179  /*std::cout << "\n ParticleIDAlg " << AlgScore.fAlgName
2180  << "\n -- Variable type: " << AlgScore.fVariableType
2181  << "\n -- Track direction: " << AlgScore.fTrackDir
2182  << "\n -- Assuming PDG: " << AlgScore.fAssumedPdg
2183  << "\n -- Number of degrees of freedom: " << AlgScore.fNdf
2184  << "\n -- Value: " << AlgScore.fValue
2185  << "\n -- Using planeMask: " << AlgScore.fPlaneMask << std::endl;*/
2186 
2187  if (AlgScore.fPlaneMask.none() || AlgScore.fPlaneMask.count() > 1) continue;
2188  int planenum = -1;
2189  if (AlgScore.fPlaneMask.test(0)) planenum = 0;
2190  if (AlgScore.fPlaneMask.test(1)) planenum = 1;
2191  if (AlgScore.fPlaneMask.test(2)) planenum = 2;
2192  if (planenum<0 || planenum>2) continue;
2193 
2194  if (AlgScore.fAlgName == "Chi2"){
2195  if (TMath::Abs(AlgScore.fAssumedPdg) == 13){ // chi2mu
2196  TrackerData.trkpidchimu[iTrk][planenum] = AlgScore.fValue;
2197  if (AlgScore.fValue<pidchi[planenum] || (pidchi[planenum] == 0. && AlgScore.fValue>0.)){
2198  pidchi[planenum] = AlgScore.fValue;
2199  pidpdg[planenum] = TMath::Abs(AlgScore.fAssumedPdg);
2200  }
2201  }
2202  else if (TMath::Abs(AlgScore.fAssumedPdg) == 2212){ // chi2pr
2203  TrackerData.trkpidchipr[iTrk][planenum] = AlgScore.fValue;
2204  if (AlgScore.fValue<pidchi[planenum] || (pidchi[planenum] == 0. && AlgScore.fValue>0.)){
2205  pidchi[planenum] = AlgScore.fValue;
2206  pidpdg[planenum] = TMath::Abs(AlgScore.fAssumedPdg);
2207  }
2208  }
2209  else if (TMath::Abs(AlgScore.fAssumedPdg) == 211){ // chi2pi
2210  TrackerData.trkpidchipi[iTrk][planenum] = AlgScore.fValue;
2211  if (AlgScore.fValue<pidchi[planenum] || (pidchi[planenum] == 0. && AlgScore.fValue>0.)){
2212  pidchi[planenum] = AlgScore.fValue;
2213  pidpdg[planenum] = TMath::Abs(AlgScore.fAssumedPdg);
2214  }
2215  }
2216  else if (TMath::Abs(AlgScore.fAssumedPdg) == 321){ // chi2ka
2217  TrackerData.trkpidchika[iTrk][planenum] = AlgScore.fValue;
2218  if (AlgScore.fValue<pidchi[planenum] || (pidchi[planenum] == 0. && AlgScore.fValue>0.)){
2219  pidchi[planenum] = AlgScore.fValue;
2220  pidpdg[planenum] = TMath::Abs(AlgScore.fAssumedPdg);
2221  }
2222  }
2223 
2224  }
2225  else if (AlgScore.fVariableType==anab::kPIDA){
2226  TrackerData.trkpidpida[iTrk][planenum] = AlgScore.fValue;
2227  }
2228 
2229  } // end loop though AlgScoresVec
2230  } // end loop over pid[ipid]
2231 
2232  // Finally, set min chi2
2233  for (size_t planenum=0; planenum<3; planenum++){
2234  TrackerData.trkpidchi[iTrk][planenum] = pidchi[planenum];
2235  TrackerData.trkpidpdg[iTrk][planenum] = pidpdg[planenum];
2236  }
2237  } // fmpid.isValid()
2238 
2239 
2240 
2241  art::FindMany<anab::Calorimetry> fmcal(trackListHandle[iTracker], evt, fCalorimetryModuleLabel[iTracker]);
2242  if (fmcal.isValid()){
2243  std::vector<const anab::Calorimetry*> calos = fmcal.at(iTrk);
2244  if (calos.size() > TrackerData.GetMaxPlanesPerTrack(iTrk)) {
2245  // if you get this message, there is probably a bug somewhere since
2246  // the calorimetry planes should be 3.
2247  mf::LogError("AnalysisTree:limits")
2248  << "the " << fTrackModuleLabel[iTracker] << " track #" << iTrk
2249  << " has " << calos.size() << " planes for calorimetry , only "
2250  << TrackerData.GetMaxPlanesPerTrack(iTrk) << " stored in tree";
2251  }
2252  for (size_t ical = 0; ical<calos.size(); ++ical){
2253  if (!calos[ical]) continue;
2254  if (!calos[ical]->PlaneID().isValid) continue;
2255  int planenum = calos[ical]->PlaneID().Plane;
2256  if (planenum<0||planenum>2) continue;
2257  TrackerData.trkke[iTrk][planenum] = calos[ical]->KineticEnergy();
2258  TrackerData.trkrange[iTrk][planenum] = calos[ical]->Range();
2259  //For now make the second argument as 13 for muons.
2260  TrackerData.trkpitchc[iTrk][planenum]= calos[ical] -> TrkPitchC();
2261  const size_t NHits = calos[ical] -> dEdx().size();
2262  TrackerData.ntrkhits[iTrk][planenum] = (int) NHits;
2263  if (NHits > TrackerData.GetMaxHitsPerTrack(iTrk, planenum)) {
2264  // if you get this error, you'll have to increase kMaxTrackHits
2265  mf::LogError("AnalysisTree:limits")
2266  << "the " << fTrackModuleLabel[iTracker] << " track #" << iTrk
2267  << " has " << NHits << " hits on calorimetry plane #" << planenum
2268  <<", only "
2269  << TrackerData.GetMaxHitsPerTrack(iTrk, planenum) << " stored in tree";
2270  }
2271  if (!isCosmics){
2272  for(size_t iTrkHit = 0; iTrkHit < NHits && iTrkHit < TrackerData.GetMaxHitsPerTrack(iTrk, planenum); ++iTrkHit) {
2273  TrackerData.trkdedx[iTrk][planenum][iTrkHit] = (calos[ical] -> dEdx())[iTrkHit];
2274  TrackerData.trkdqdx[iTrk][planenum][iTrkHit] = (calos[ical] -> dQdx())[iTrkHit];
2275  TrackerData.trkresrg[iTrk][planenum][iTrkHit] = (calos[ical] -> ResidualRange())[iTrkHit];
2276  const auto& TrkPos = (calos[ical] -> XYZ())[iTrkHit];
2277  auto& TrkXYZ = TrackerData.trkxyz[iTrk][planenum][iTrkHit];
2278  TrkXYZ[0] = TrkPos.X();
2279  TrkXYZ[1] = TrkPos.Y();
2280  TrkXYZ[2] = TrkPos.Z();
2281  TrackerData.trkxp[iTrk][planenum][iTrkHit] = TrkPos.X();
2282  TrackerData.trkyp[iTrk][planenum][iTrkHit] = TrkPos.Y();
2283  TrackerData.trkzp[iTrk][planenum][iTrkHit] = TrkPos.Z();
2284  } // for track hits
2285  }
2286  } // for calorimetry info
2287  if(TrackerData.ntrkhits[iTrk][0] > TrackerData.ntrkhits[iTrk][1] && TrackerData.ntrkhits[iTrk][0] > TrackerData.ntrkhits[iTrk][2]) TrackerData.trkpidbestplane[iTrk] = 0;
2288  else if(TrackerData.ntrkhits[iTrk][1] > TrackerData.ntrkhits[iTrk][0] && TrackerData.ntrkhits[iTrk][1] > TrackerData.ntrkhits[iTrk][2]) TrackerData.trkpidbestplane[iTrk] = 1;
2289  else if(TrackerData.ntrkhits[iTrk][2] > TrackerData.ntrkhits[iTrk][0] && TrackerData.ntrkhits[iTrk][2] > TrackerData.ntrkhits[iTrk][1]) TrackerData.trkpidbestplane[iTrk] = 2;
2290  else if(TrackerData.ntrkhits[iTrk][2] == TrackerData.ntrkhits[iTrk][0] && TrackerData.ntrkhits[iTrk][2] > TrackerData.ntrkhits[iTrk][1]) TrackerData.trkpidbestplane[iTrk] = 2;
2291  else if(TrackerData.ntrkhits[iTrk][2] == TrackerData.ntrkhits[iTrk][1] && TrackerData.ntrkhits[iTrk][2] > TrackerData.ntrkhits[iTrk][0]) TrackerData.trkpidbestplane[iTrk] = 2;
2292  else if(TrackerData.ntrkhits[iTrk][1] == TrackerData.ntrkhits[iTrk][0] && TrackerData.ntrkhits[iTrk][1] > TrackerData.ntrkhits[iTrk][2]) TrackerData.trkpidbestplane[iTrk] = 0;
2293  else if(TrackerData.ntrkhits[iTrk][1] == TrackerData.ntrkhits[iTrk][0] && TrackerData.ntrkhits[iTrk][1] == TrackerData.ntrkhits[iTrk][2]) TrackerData.trkpidbestplane[iTrk] = 2;
2294  } // if has calorimetry info
2295 
2296  //track truth information
2297  if (isMC){
2298  //get the hits on each plane
2299  art::FindManyP<recob::Hit> fmht(trackListHandle[iTracker], evt, fTrackModuleLabel[iTracker]);
2300  std::vector< art::Ptr<recob::Hit> > allHits = fmht.at(iTrk);
2301  std::vector< art::Ptr<recob::Hit> > hits[kNplanes];
2302 
2303  for(size_t ah = 0; ah < allHits.size(); ++ah){
2304  if (/* allHits[ah]->WireID().Plane >= 0 && */ // always true
2305  allHits[ah]->WireID().Plane < 3){
2306  hits[allHits[ah]->WireID().Plane].push_back(allHits[ah]);
2307  }
2308  }
2309 
2310  for (size_t ipl = 0; ipl < 3; ++ipl){
2311  TrackerData.trkidtruth_recoutils_totaltrueenergy[iTrk][ipl] = RecoUtils::TrueParticleIDFromTotalTrueEnergy(clockData, hits[ipl]);
2312  TrackerData.trkidtruth_recoutils_totalrecocharge[iTrk][ipl] = RecoUtils::TrueParticleIDFromTotalRecoCharge(clockData, hits[ipl]);
2313  TrackerData.trkidtruth_recoutils_totalrecohits[iTrk][ipl] = RecoUtils::TrueParticleIDFromTotalRecoHits(clockData, hits[ipl]);
2314  double maxe = 0;
2315  HitsPurity(clockData, hits[ipl],TrackerData.trkidtruth[iTrk][ipl],TrackerData.trkpurtruth[iTrk][ipl],maxe);
2316  //std::cout<<"\n"<<iTracker<<"\t"<<iTrk<<"\t"<<ipl<<"\t"<<trkidtruth[iTracker][iTrk][ipl]<<"\t"<<trkpurtruth[iTracker][iTrk][ipl]<<"\t"<<maxe;
2317  if (TrackerData.trkidtruth[iTrk][ipl]>0){
2318  const art::Ptr<simb::MCTruth> mc = partInventory->TrackIdToMCTruth_P(TrackerData.trkidtruth[iTrk][ipl]);
2319  TrackerData.trkorigin[iTrk][ipl] = mc->Origin();
2320  const simb::MCParticle *particle = partInventory->TrackIdToParticle_P(TrackerData.trkidtruth[iTrk][ipl]);
2321  double tote = 0;
2322  const std::vector<const sim::IDE*> vide(backTracker->TrackIdToSimIDEs_Ps(TrackerData.trkidtruth[iTrk][ipl]));
2323  for (auto ide: vide) {
2324  tote += ide->energy;
2325  TrackerData.trksimIDEenergytruth[iTrk][ipl] = ide->energy;
2326  TrackerData.trksimIDExtruth[iTrk][ipl] = ide->x;
2327  TrackerData.trksimIDEytruth[iTrk][ipl] = ide->y;
2328  TrackerData.trksimIDEztruth[iTrk][ipl] = ide->z;
2329  }
2330  TrackerData.trkpdgtruth[iTrk][ipl] = particle->PdgCode();
2331  TrackerData.trkefftruth[iTrk][ipl] = maxe/(tote/kNplanes); //tote include both induction and collection energies
2332  //std::cout<<"\n"<<trkpdgtruth[iTracker][iTrk][ipl]<<"\t"<<trkefftruth[iTracker][iTrk][ipl];
2333  }
2334  }
2335  }//end if (isMC)
2336  }//end loop over track
2337  }//end loop over track module labels
2338  }// end (fSaveTrackInfo)
2339 
2340  /*trkf::TrackMomentumCalculator trkm;
2341  std::cout<<"\t"<<trkm.GetTrackMomentum(200,2212)<<"\t"<<trkm.GetTrackMomentum(-10, 13)<<"\t"<<trkm.GetTrackMomentum(300,-19)<<"\n";
2342 */
2343  //mc truth information
2344  if (isMC){
2345  if (fSaveCryInfo){
2346  //store cry (cosmic generator information)
2347  fData->mcevts_truthcry = mclistcry.size();
2348  fData->cry_no_primaries = nCryPrimaries;
2349  //fData->cry_no_primaries;
2350  for(Int_t iPartc = 0; iPartc < mctruthcry->NParticles(); ++iPartc){
2351  const simb::MCParticle& partc(mctruthcry->GetParticle(iPartc));
2352  fData->cry_primaries_pdg[iPartc]=partc.PdgCode();
2353  fData->cry_Eng[iPartc]=partc.E();
2354  fData->cry_Px[iPartc]=partc.Px();
2355  fData->cry_Py[iPartc]=partc.Py();
2356  fData->cry_Pz[iPartc]=partc.Pz();
2357  fData->cry_P[iPartc]=partc.P();
2358  fData->cry_StartPointx[iPartc] = partc.Vx();
2359  fData->cry_StartPointy[iPartc] = partc.Vy();
2360  fData->cry_StartPointz[iPartc] = partc.Vz();
2361  fData->cry_status_code[iPartc]=partc.StatusCode();
2362  fData->cry_mass[iPartc]=partc.Mass();
2363  fData->cry_trackID[iPartc]=partc.TrackId();
2364  fData->cry_ND[iPartc]=partc.NumberDaughters();
2365  fData->cry_mother[iPartc]=partc.Mother();
2366  } // for cry particles
2367  }// end fSaveCryInfo
2368 
2369  fData->mcevts_truth = mclist.size();
2370  //Brailsford 2017/10/16
2371  //Issue 17917
2372  //To keep a 1:1 between neutrinos and 'flux' we need the assns
2373  art::FindOneP<simb::MCFlux> fmFluxNeutrino(mctruthListHandle, evt, fGenieGenModuleLabel);
2374 
2375  if (fData->mcevts_truth > 0){//at least one mc record
2376  if (fSaveGenieInfo){
2377 
2378  //Brailsford 2017/10/16
2379  //Issue 17917
2380  //Loop over every truth in the spill rather than just looking at the first one.
2381  //Because MCTruth could be a neutrino OR something else (e.g. cosmics) we are going to have to count up how many neutrinos there are
2382  fData->mcevts_truth = 0;
2383  for (unsigned int i_mctruth = 0; i_mctruth < mclist.size(); i_mctruth++){
2384  //fetch an mctruth
2385  art::Ptr<simb::MCTruth> curr_mctruth = mclist[i_mctruth];
2386  //Check if it's a neutrino
2387  if (!curr_mctruth->NeutrinoSet()) continue;
2388  fData->nuPDG_truth[i_mctruth] = curr_mctruth->GetNeutrino().Nu().PdgCode();
2389  fData->ccnc_truth[i_mctruth] = curr_mctruth->GetNeutrino().CCNC();
2390  fData->mode_truth[i_mctruth] = curr_mctruth->GetNeutrino().Mode();
2391  fData->Q2_truth[i_mctruth] = curr_mctruth->GetNeutrino().QSqr();
2392  fData->W_truth[i_mctruth] = curr_mctruth->GetNeutrino().W();
2393  fData->hitnuc_truth[i_mctruth] = curr_mctruth->GetNeutrino().HitNuc();
2394  fData->enu_truth[i_mctruth] = curr_mctruth->GetNeutrino().Nu().E();
2395  fData->nuvtxx_truth[i_mctruth] = curr_mctruth->GetNeutrino().Nu().Vx();
2396  fData->nuvtxy_truth[i_mctruth] = curr_mctruth->GetNeutrino().Nu().Vy();
2397  fData->nuvtxz_truth[i_mctruth] = curr_mctruth->GetNeutrino().Nu().Vz();
2398  if (curr_mctruth->GetNeutrino().Nu().P()){
2399  fData->nu_dcosx_truth[i_mctruth] = curr_mctruth->GetNeutrino().Nu().Px()/curr_mctruth->GetNeutrino().Nu().P();
2400  fData->nu_dcosy_truth[i_mctruth] = curr_mctruth->GetNeutrino().Nu().Py()/curr_mctruth->GetNeutrino().Nu().P();
2401  fData->nu_dcosz_truth[i_mctruth] = curr_mctruth->GetNeutrino().Nu().Pz()/curr_mctruth->GetNeutrino().Nu().P();
2402  }
2403  fData->lep_mom_truth[i_mctruth] = curr_mctruth->GetNeutrino().Lepton().P();
2404  if (curr_mctruth->GetNeutrino().Lepton().P()){
2405  fData->lep_dcosx_truth[i_mctruth] = curr_mctruth->GetNeutrino().Lepton().Px()/curr_mctruth->GetNeutrino().Lepton().P();
2406  fData->lep_dcosy_truth[i_mctruth] = curr_mctruth->GetNeutrino().Lepton().Py()/curr_mctruth->GetNeutrino().Lepton().P();
2407  fData->lep_dcosz_truth[i_mctruth] = curr_mctruth->GetNeutrino().Lepton().Pz()/curr_mctruth->GetNeutrino().Lepton().P();
2408  }
2409  //Brailsford
2410  //2017/10/17
2411  //Issue 12918
2412  //Use the art::Ptr key as the neutrino's unique ID
2413  fData->nuID_truth[i_mctruth] = curr_mctruth.key();
2414  //We need to also store N 'flux' neutrinos per event so now check that the FindOneP is valid and, if so, use it!
2415  if (fmFluxNeutrino.isValid()){
2416  art::Ptr<simb::MCFlux> curr_mcflux = fmFluxNeutrino.at(i_mctruth);
2417  fData->tpx_flux[i_mctruth] = curr_mcflux->ftpx;
2418  fData->tpy_flux[i_mctruth] = curr_mcflux->ftpy;
2419  fData->tpz_flux[i_mctruth] = curr_mcflux->ftpz;
2420  fData->tptype_flux[i_mctruth] = curr_mcflux->ftptype;
2421  }
2422 
2423  //Let's increase the neutrino count
2424  fData->mcevts_truth++;
2425  }
2426 
2427  if (mctruth->NeutrinoSet()){
2428  //Brailsford 2017/10/16
2429  //Issue 17917
2430  //Bypass this section of filling code as there can now be multiple neutrinos per TTree::entry stored in the TTree
2431  /*
2432  fData->nuPDG_truth = mctruth->GetNeutrino().Nu().PdgCode();
2433  fData->ccnc_truth = mctruth->GetNeutrino().CCNC();
2434  fData->mode_truth = mctruth->GetNeutrino().Mode();
2435  fData->Q2_truth = mctruth->GetNeutrino().QSqr();
2436  fData->W_truth = mctruth->GetNeutrino().W();
2437  fData->hitnuc_truth = mctruth->GetNeutrino().HitNuc();
2438  fData->enu_truth = mctruth->GetNeutrino().Nu().E();
2439  fData->nuvtxx_truth = mctruth->GetNeutrino().Nu().Vx();
2440  fData->nuvtxy_truth = mctruth->GetNeutrino().Nu().Vy();
2441  fData->nuvtxz_truth = mctruth->GetNeutrino().Nu().Vz();
2442  if (mctruth->GetNeutrino().Nu().P()){
2443  fData->nu_dcosx_truth = mctruth->GetNeutrino().Nu().Px()/mctruth->GetNeutrino().Nu().P();
2444  fData->nu_dcosy_truth = mctruth->GetNeutrino().Nu().Py()/mctruth->GetNeutrino().Nu().P();
2445  fData->nu_dcosz_truth = mctruth->GetNeutrino().Nu().Pz()/mctruth->GetNeutrino().Nu().P();
2446  }
2447  fData->lep_mom_truth = mctruth->GetNeutrino().Lepton().P();
2448  if (mctruth->GetNeutrino().Lepton().P()){
2449  fData->lep_dcosx_truth = mctruth->GetNeutrino().Lepton().Px()/mctruth->GetNeutrino().Lepton().P();
2450  fData->lep_dcosy_truth = mctruth->GetNeutrino().Lepton().Py()/mctruth->GetNeutrino().Lepton().P();
2451  fData->lep_dcosz_truth = mctruth->GetNeutrino().Lepton().Pz()/mctruth->GetNeutrino().Lepton().P();
2452  }
2453  //flux information
2454  art::Ptr<simb::MCFlux> mcflux = fluxlist[imc];
2455  fData->tpx_flux = mcflux->ftpx;
2456  fData->tpy_flux = mcflux->ftpy;
2457  fData->tpz_flux = mcflux->ftpz;
2458  fData->tptype_flux = mcflux->ftptype;
2459  */
2460 
2461  //genie particles information
2462  fData->genie_no_primaries = mctruth->NParticles();
2463 
2464  size_t StoreParticles = std::min((size_t) fData->genie_no_primaries, fData->GetMaxGeniePrimaries());
2465  if (fData->genie_no_primaries > (int) StoreParticles) {
2466  // got this error? it might be a bug,
2467  // since the structure should have enough room for everything
2468  mf::LogError("AnalysisTree:limits") << "event has "
2469  << fData->genie_no_primaries << " MC particles, only "
2470  << StoreParticles << " stored in tree";
2471  }
2472  for(size_t iPart = 0; iPart < StoreParticles; ++iPart){
2473  const simb::MCParticle& part(mctruth->GetParticle(iPart));
2474  fData->genie_primaries_pdg[iPart]=part.PdgCode();
2475  fData->genie_Eng[iPart]=part.E();
2476  fData->genie_Px[iPart]=part.Px();
2477  fData->genie_Py[iPart]=part.Py();
2478  fData->genie_Pz[iPart]=part.Pz();
2479  fData->genie_P[iPart]=part.P();
2480  fData->genie_status_code[iPart]=part.StatusCode();
2481  fData->genie_mass[iPart]=part.Mass();
2482  fData->genie_trackID[iPart]=part.TrackId();
2483  fData->genie_ND[iPart]=part.NumberDaughters();
2484  fData->genie_mother[iPart]=part.Mother();
2485  } // for particle
2486  } //if neutrino set
2487  }// end (fSaveGenieInfo)
2488 
2489  //GEANT particles information
2490  if (fSaveGeantInfo){
2491  const sim::ParticleList& plist = partInventory->ParticleList();
2492 
2493  std::string pri("primary");
2494  int primary=0;
2495  int active = 0;
2496  int geant_particle=0;
2497  sim::ParticleList::const_iterator itPart = plist.begin(),
2498  pend = plist.end(); // iterator to pairs (track id, particle)
2499 
2500  for(size_t iPart = 0; (iPart < plist.size()) && (itPart != pend); ++iPart){
2501  const simb::MCParticle* pPart = (itPart++)->second;
2502  if (!pPart) {
2503  throw art::Exception(art::errors::LogicError)
2504  << "GEANT particle #" << iPart << " returned a null pointer";
2505  }
2506 
2507  ++geant_particle;
2508  bool isPrimary = pPart->Process() == pri;
2509  if (isPrimary) ++primary;
2510 
2511  int TrackID = pPart->TrackId();
2512 
2513  TVector3 mcstart, mcend;
2514  double plen = length(*pPart, mcstart, mcend);
2515 
2516  bool isActive = plen != 0;
2517  if (plen) active++;
2518 
2519  if (iPart < fData->GetMaxGEANTparticles()) {
2520  if (pPart->E()>fG4minE||isPrimary){
2521  fData->process_primary[iPart] = int(isPrimary);
2522  fData->processname[iPart]= pPart->Process();
2523  fData->Mother[iPart]=pPart->Mother();
2524  fData->TrackId[iPart]=TrackID;
2525  fData->pdg[iPart]=pPart->PdgCode();
2526  fData->status[iPart] = pPart->StatusCode();
2527  fData->Eng[iPart]=pPart->E();
2528  fData->EndE[iPart]=pPart->EndE();
2529  fData->Mass[iPart]=pPart->Mass();
2530  fData->Px[iPart]=pPart->Px();
2531  fData->Py[iPart]=pPart->Py();
2532  fData->Pz[iPart]=pPart->Pz();
2533  fData->P[iPart]=pPart->Momentum().Vect().Mag();
2534  fData->StartPointx[iPart]=pPart->Vx();
2535  fData->StartPointy[iPart]=pPart->Vy();
2536  fData->StartPointz[iPart]=pPart->Vz();
2537  fData->StartT[iPart] = pPart->T();
2538  fData->EndPointx[iPart]=pPart->EndPosition()[0];
2539  fData->EndPointy[iPart]=pPart->EndPosition()[1];
2540  fData->EndPointz[iPart]=pPart->EndPosition()[2];
2541  fData->EndT[iPart] = pPart->EndT();
2542  fData->theta[iPart] = pPart->Momentum().Theta();
2543  fData->phi[iPart] = pPart->Momentum().Phi();
2544  fData->theta_xz[iPart] = std::atan2(pPart->Px(), pPart->Pz());
2545  fData->theta_yz[iPart] = std::atan2(pPart->Py(), pPart->Pz());
2546  fData->pathlen[iPart] = plen;
2547  fData->NumberDaughters[iPart]=pPart->NumberDaughters();
2548  fData->inTPCActive[iPart] = int(isActive);
2549  if (isActive){
2550  fData->StartPointx_tpcAV[iPart] = mcstart.X();
2551  fData->StartPointy_tpcAV[iPart] = mcstart.Y();
2552  fData->StartPointz_tpcAV[iPart] = mcstart.Z();
2553  fData->EndPointx_tpcAV[iPart] = mcend.X();
2554  fData->EndPointy_tpcAV[iPart] = mcend.Y();
2555  fData->EndPointz_tpcAV[iPart] = mcend.Z();
2556  }
2557  //Brailsford
2558  //2017/10/17
2559  //Issue 17918
2560  //Get the mother neutrino of this particle and use the hosting art::Ptr key as the matched ID
2561  const art::Ptr<simb::MCTruth>& matched_mctruth = partInventory->ParticleToMCTruth_P(pPart);
2562  fData->MotherNuId[iPart] = matched_mctruth.key();
2563  }
2564  //access auxiliary detector parameters
2565  if (fSaveAuxDetInfo) {
2566  unsigned short nAD = 0; // number of cells that particle hit
2567 
2568  // find deposit of this particle in each of the detector cells
2569  for (const sim::AuxDetSimChannel* c: fAuxDetSimChannels) {
2570 
2571  // find if this cell has a contribution (IDE) from this particle,
2572  // and which one
2573  const std::vector<sim::AuxDetIDE>& setOfIDEs = c->AuxDetIDEs();
2574  // using a C++ "lambda" function here; this one:
2575  // - sees only TrackID from the current scope
2576  // - takes one parameter: the AuxDetIDE to be tested
2577  // - returns if that IDE belongs to the track we are looking for
2578  std::vector<sim::AuxDetIDE>::const_iterator iIDE
2579  = std::find_if(
2580  setOfIDEs.begin(), setOfIDEs.end(),
2581  [TrackID](const sim::AuxDetIDE& IDE){ return IDE.trackID == TrackID; }
2582  );
2583  if (iIDE == setOfIDEs.end()) continue;
2584 
2585  // now iIDE points to the energy released by the track #i (TrackID)
2586 
2587  // look for IDE with matching trackID
2588  // find trackIDs stored in setOfIDEs with the same trackID, but negative,
2589  // this is an untracked particle who's energy should be added as deposited by this original trackID
2590  float totalE = 0.; // total energy deposited around by the GEANT particle in this cell
2591  for(const auto& adtracks: setOfIDEs) {
2592  if( fabs(adtracks.trackID) == TrackID )
2593  totalE += adtracks.energyDeposited;
2594  } // for
2595 
2596  // fill the structure
2597  if (nAD < kMaxAuxDets) {
2598  fData->AuxDetID[iPart][nAD] = c->AuxDetID();
2599  fData->entryX[iPart][nAD] = iIDE->entryX;
2600  fData->entryY[iPart][nAD] = iIDE->entryY;
2601  fData->entryZ[iPart][nAD] = iIDE->entryZ;
2602  fData->entryT[iPart][nAD] = iIDE->entryT;
2603  fData->exitX[iPart][nAD] = iIDE->exitX;
2604  fData->exitY[iPart][nAD] = iIDE->exitY;
2605  fData->exitZ[iPart][nAD] = iIDE->exitZ;
2606  fData->exitT[iPart][nAD] = iIDE->exitT;
2607  fData->exitPx[iPart][nAD] = iIDE->exitMomentumX;
2608  fData->exitPy[iPart][nAD] = iIDE->exitMomentumY;
2609  fData->exitPz[iPart][nAD] = iIDE->exitMomentumZ;
2610  fData->CombinedEnergyDep[iPart][nAD] = totalE;
2611  }
2612  ++nAD;
2613  } // for aux det sim channels
2614  fData->NAuxDets[iPart] = nAD;
2615 
2616  if (nAD > kMaxAuxDets) {
2617  // got this error? consider increasing kMaxAuxDets
2618  mf::LogError("AnalysisTree:limits") << "particle #" << iPart
2619  << " touches " << nAD << " auxiliary detector cells, only "
2620  << kMaxAuxDets << " of them are saved in the tree";
2621  } // if too many detector cells
2622  } // if (fSaveAuxDetInfo)
2623  }
2624  else if (iPart == fData->GetMaxGEANTparticles()) {
2625  // got this error? it might be a bug,
2626  // since the structure should have enough room for everything
2627  mf::LogError("AnalysisTree:limits") << "event has "
2628  << plist.size() << " MC particles, only "
2629  << fData->GetMaxGEANTparticles() << " will be stored in tree";
2630  }
2631  } // for particles
2632 
2633  fData->geant_list_size_in_tpcAV = active;
2634  fData->no_primaries = primary;
2635  fData->geant_list_size = geant_particle;
2636 
2637  MF_LOG_DEBUG("AnalysisTree") << "Counted "
2638  << fData->geant_list_size << " GEANT particles ("
2639  << fData->geant_list_size_in_tpcAV << " in AV), "
2640  << fData->no_primaries << " primaries, "
2641  << fData->genie_no_primaries << " GENIE particles";
2642 
2643  FillWith(fData->MergedId, 0);
2644 
2645  // helper map track ID => index
2646  std::map<int, size_t> TrackIDtoIndex;
2647  const size_t nTrackIDs = fData->TrackId.size();
2648  for (size_t index = 0; index < nTrackIDs; ++index)
2649  TrackIDtoIndex.emplace(fData->TrackId[index], index);
2650 
2651  // for each particle, consider all the direct ancestors with the same
2652  // PDG ID, and mark them as belonging to the same "group"
2653  // (having the same MergedId)
2654  int currentMergedId = 1;
2655  for(size_t iPart = fData->geant_list_size; iPart-- > 0; ) {
2656  // if the particle already belongs to a group, don't bother
2657  if (fData->MergedId[iPart]) continue;
2658  // the particle starts its own group
2659  fData->MergedId[iPart] = currentMergedId;
2660 
2661  // look in the ancestry, one by one
2662  int currentMotherTrackId = fData->Mother[iPart];
2663  while(currentMotherTrackId > 0) {
2664  // find the mother (we have its track ID in currentMotherTrackId)
2665  std::map<int, size_t>::const_iterator iMother
2666  = TrackIDtoIndex.find(currentMotherTrackId);
2667  if (iMother == TrackIDtoIndex.end()) break; // no mother found
2668  size_t currentMotherIndex = iMother->second;
2669  // if the mother particle is of a different type,
2670  // don't bother with iPart ancestry any further
2671  if (fData->pdg[iPart] != fData->pdg[currentMotherIndex]) break;
2672 
2673  // group the "current mother" (actually, ancestor) with iPart
2674  fData->MergedId[currentMotherIndex] = currentMergedId;
2675  // then consider the grandmother (one generation earlier)
2676  currentMotherTrackId = fData->Mother[currentMotherIndex];
2677  } // while ancestry exists
2678  ++currentMergedId;
2679  } // for merging check
2680  } // if (fSaveGeantInfo)
2681 
2682  }//if (mcevts_truth)
2683  }//if (isMC){
2684  //fData->taulife = detprop->ElectronLifetime();
2685  fTree->Fill();
2686 
2687  if (mf::isDebugEnabled()) {
2688  // use mf::LogDebug instead of MF_LOG_DEBUG because we reuse it in many lines;
2689  // thus, we protect this part of the code with the line above
2690  mf::LogDebug logStream("AnalysisTreeStructure");
2691  logStream
2692  << "Tree data structure contains:"
2693  << "\n - " << fData->no_hits << " hits (" << fData->GetMaxHits() << ")"
2694  << "\n - " << fData->genie_no_primaries << " genie primaries (" << fData->GetMaxGeniePrimaries() << ")"
2695  << "\n - " << fData->geant_list_size << " GEANT particles (" << fData->GetMaxGEANTparticles() << "), "
2696  << fData->no_primaries << " primaries"
2697  << "\n - " << fData->geant_list_size_in_tpcAV << " GEANT particles in AV "
2698  << "\n - " << ((int) fData->kNTracker) << " trackers:"
2699  ;
2700 
2701  size_t iTracker = 0;
2702  for (auto tracker = fData->TrackData.cbegin();
2703  tracker != fData->TrackData.cend(); ++tracker, ++iTracker
2704  ) {
2705  logStream
2706  << "\n -> " << tracker->ntracks << " " << fTrackModuleLabel[iTracker]
2707  << " tracks (" << tracker->GetMaxTracks() << ")"
2708  ;
2709  for (int iTrk = 0; iTrk < tracker->ntracks; ++iTrk) {
2710  logStream << "\n [" << iTrk << "] "<< tracker->ntrkhits[iTrk][0];
2711  for (size_t ipl = 1; ipl < tracker->GetMaxPlanesPerTrack(iTrk); ++ipl)
2712  logStream << " + " << tracker->ntrkhits[iTrk][ipl];
2713  logStream << " hits (" << tracker->GetMaxHitsPerTrack(iTrk, 0);
2714  for (size_t ipl = 1; ipl < tracker->GetMaxPlanesPerTrack(iTrk); ++ipl)
2715  logStream << " + " << tracker->GetMaxHitsPerTrack(iTrk, ipl);
2716  logStream << ")";
2717  } // for tracks
2718  } // for trackers
2719  } // if logging enabled
2720 
2721  // if we don't use a permanent buffer (which can be huge),
2722  // delete the current buffer, and we'll create a new one on the next event
2723  if (!fUseBuffer) {
2724  MF_LOG_DEBUG("AnalysisTreeStructure") << "Freeing the tree data structure";
2725  DestroyData();
2726  }
2727 } // icarus::AnalysisTree::analyze()
icarus::AnalysisTreeDataStruct::ClearLocalData
void ClearLocalData()
Clear all fields if this object (not the tracker algorithm data)
Definition: icaruscode/icaruscode/Analysis/AnalysisTree_module.cc:1189
kNplanes
constexpr int kNplanes
Definition: icaruscode/icaruscode/Analysis/AnalysisTree_module.cc:140
second
second
Definition: detectorproperties_icarus.fcl:8
tca::TCEvent::event
unsigned int event
Definition: DataStructs.h:634
icarus::AnalysisTreeDataStruct::NAuxDets
std::vector< UShort_t > NAuxDets
Number of AuxDets crossed by this particle.
Definition: icaruscode/icaruscode/Analysis/AnalysisTree_module.cc:488
sim::SimChannel
Energy deposited on a readout channel by simulated tracks.
Definition: SimChannel.h:145
RecoUtils::TrueParticleIDFromTotalTrueEnergy
int TrueParticleIDFromTotalTrueEnergy(detinfo::DetectorClocksData const &clockData, const std::vector< art::Ptr< recob::Hit > > &hits, bool rollup_unsaved_ids=1)
Definition: icaruscode/icaruscode/CRT/CRTUtils/RecoUtils.cc:30
icarus::AnalysisTreeDataStruct::SubRunData
SubRunData_t SubRunData
subrun data collected at begin of subrun
Definition: icaruscode/icaruscode/Analysis/AnalysisTree_module.cc:343
tca::TCEvent::run
unsigned int run
Definition: DataStructs.h:635
icarus::AnalysisTreeDataStruct::ResizeGenie
void ResizeGenie(int nPrimaries)
Resize the data strutcure for Genie primaries.
Definition: icaruscode/icaruscode/Analysis/AnalysisTree_module.cc:1434
cerr
BEGIN_PROLOG could also be cerr
Definition: messageservice.fcl:10
icarus::AnalysisTreeDataStruct::ResizeMCNeutrino
void ResizeMCNeutrino(int nNeutrinos)
Resize the data structure for MCNeutrino particles.
Definition: icaruscode/icaruscode/Analysis/AnalysisTree_module.cc:1342
gen_crt_frags.beam
def beam
Definition: gen_crt_frags.py:212
icarus::AnalysisTree::fSaveTrackInfo
bool fSaveTrackInfo
whether to extract and save Hit information
Definition: icaruscode/icaruscode/Analysis/AnalysisTree_module.cc:716
icarus::AnalysisTreeDataStruct::GetMaxGEANTparticles
size_t GetMaxGEANTparticles() const
Returns the number of GEANT particles for which memory is allocated.
Definition: icaruscode/icaruscode/Analysis/AnalysisTree_module.cc:576
icarus::AnalysisTreeDataStruct::entryY
AuxDetMCData_t< Float_t > entryY
Entry Y position of particle into AuxDet.
Definition: icaruscode/icaruscode/Analysis/AnalysisTree_module.cc:491
kMaxVertices
constexpr unsigned short kMaxVertices
Definition: icaruscode/icaruscode/Analysis/AnalysisTree_module.cc:144
icarus::AnalysisTreeDataStruct::GetMaxGeniePrimaries
size_t GetMaxGeniePrimaries() const
Returns the number of GENIE primaries for which memory is allocated.
Definition: icaruscode/icaruscode/Analysis/AnalysisTree_module.cc:579
icarus::AnalysisTreeDataStruct::exitPz
AuxDetMCData_t< Float_t > exitPz
Exit z momentum of particle out of AuxDet.
Definition: icaruscode/icaruscode/Analysis/AnalysisTree_module.cc:500
recob::Track::NumberTrajectoryPoints
size_t NumberTrajectoryPoints() const
Various functions related to the presence and the number of (valid) points.
Definition: bj/lardataobj/RecoBase/Track.h:102
recob::Track::VertexDirection
Vector_t VertexDirection() const
Definition: bj/lardataobj/RecoBase/Track.h:132
icarus::AnalysisTree::bdist
double bdist(const recob::tracking::Point_t &pos)
Definition: icaruscode/icaruscode/Analysis/AnalysisTree_module.cc:2770
icarus::AnalysisTree::fSaveAuxDetInfo
bool fSaveAuxDetInfo
whether to extract and save auxiliary detector data
Definition: icaruscode/icaruscode/Analysis/AnalysisTree_module.cc:711
track
process_name use argoneut_mc_hitfinder track
Definition: standard_reco.fcl:36
icarus::AnalysisTreeDataStruct::entryX
AuxDetMCData_t< Float_t > entryX
Entry X position of particle into AuxDet.
Definition: icaruscode/icaruscode/Analysis/AnalysisTree_module.cc:490
icarus::AnalysisTree::fUseBuffer
bool fUseBuffer
whether to use a permanent buffer (faster, huge memory)
Definition: icaruscode/icaruscode/Analysis/AnalysisTree_module.cc:710
anab::sParticleIDAlgScores::fValue
float fValue
Result of Particle ID algorithm/test.
Definition: ParticleID.h:28
icarus::AnalysisTree::DestroyData
void DestroyData()
Destroy the local buffers (existing branches will point to invalid address!)
Definition: icaruscode/icaruscode/Analysis/AnalysisTree_module.cc:772
icarus::AnalysisTree::fCosmicTaggerAssocLabel
std::vector< std::string > fCosmicTaggerAssocLabel
whether to extract and save Vertex information
Definition: icaruscode/icaruscode/Analysis/AnalysisTree_module.cc:719
icarus::AnalysisTreeDataStruct::entryZ
AuxDetMCData_t< Float_t > entryZ
Entry Z position of particle into AuxDet.
Definition: icaruscode/icaruscode/Analysis/AnalysisTree_module.cc:492
icarus::AnalysisTreeDataStruct::exitPx
AuxDetMCData_t< Float_t > exitPx
Exit x momentum of particle out of AuxDet.
Definition: icaruscode/icaruscode/Analysis/AnalysisTree_module.cc:498
anab::sParticleIDAlgScores::fPlaneMask
std::bitset< 8 > fPlaneMask
Bitset for PlaneID used by algorithm, allowing for multiple planes and up to 8 total planes...
Definition: ParticleID.h:29
icarus::AnalysisTreeDataStruct::entryT
AuxDetMCData_t< Float_t > entryT
Entry T position of particle into AuxDet.
Definition: icaruscode/icaruscode/Analysis/AnalysisTree_module.cc:493
geo::fhicl::PlaneID
IDparameter< geo::PlaneID > PlaneID
Member type of validated geo::PlaneID parameter.
Definition: geo_types_fhicl.h:375
anab::sParticleIDAlgScores::fAlgName
std::string fAlgName
&lt; determined particle ID
Definition: ParticleID.h:23
kMaxHits
constexpr int kMaxHits
Definition: icaruscode/icaruscode/Analysis/AnalysisTree_module.cc:141
sim::AuxDetSimChannel
Collection of particles crossing one auxiliary detector cell.
Definition: AuxDetSimChannel.h:60
icarus::AnalysisTree::HitsPurity
void HitsPurity(detinfo::DetectorClocksData const &clockData, std::vector< art::Ptr< recob::Hit > > const &hits, Int_t &trackid, Float_t &purity, double &maxe)
Definition: icaruscode/icaruscode/Analysis/AnalysisTree_module.cc:2729
icarus::AnalysisTreeDataStruct::ResizeGEANT
void ResizeGEANT(int nParticles)
Resize the data strutcure for GEANT particles.
Definition: icaruscode/icaruscode/Analysis/AnalysisTree_module.cc:1373
anab::sParticleIDAlgScores::fVariableType
kVariableType fVariableType
Variable type enum: defined in ParticleID_VariableTypeEnums.h. Set to kNotSet by default.
Definition: ParticleID.h:24
recob::tracking::Vector_t
ROOT::Math::DisplacementVector3D< ROOT::Math::Cartesian3D< Coord_t >, ROOT::Math::GlobalCoordinateSystemTag > Vector_t
Type for representation of momenta in 3D space. See recob::tracking::Coord_t for more details on the ...
Definition: TrackingTypes.h:29
icarus::AnalysisTreeDataStruct::no_primaries
Int_t no_primaries
! how many particles there is currently room for
Definition: icaruscode/icaruscode/Analysis/AnalysisTree_module.cc:441
icarus::AnalysisTreeDataStruct::mcevts_truth
Int_t mcevts_truth
! The number of MCNeutrinos there is currently room for
Definition: icaruscode/icaruscode/Analysis/AnalysisTree_module.cc:380
icarus::AnalysisTree::fSaveHitInfo
bool fSaveHitInfo
whether to extract and save Geant information
Definition: icaruscode/icaruscode/Analysis/AnalysisTree_module.cc:715
geo::fhicl::WireID
IDparameter< geo::WireID > WireID
Member type of validated geo::WireID parameter.
Definition: geo_types_fhicl.h:407
icarus::AnalysisTreeDataStruct::exitT
AuxDetMCData_t< Float_t > exitT
Exit T position of particle out of AuxDet.
Definition: icaruscode/icaruscode/Analysis/AnalysisTree_module.cc:497
icarus::AnalysisTreeDataStruct::CombinedEnergyDep
AuxDetMCData_t< Float_t > CombinedEnergyDep
Sum energy of all particles with this trackID (+ID or -ID) in AuxDet.
Definition: icaruscode/icaruscode/Analysis/AnalysisTree_module.cc:501
recob::Track::Vertex
Point_t const & Vertex() const
Definition: bj/lardataobj/RecoBase/Track.h:124
icarus::AnalysisTreeDataStruct::exitX
AuxDetMCData_t< Float_t > exitX
Exit X position of particle out of AuxDet.
Definition: icaruscode/icaruscode/Analysis/AnalysisTree_module.cc:494
icarus::ns::util::end
auto end(FixedBins< T, C > const &) noexcept
Definition: FixedBins.h:585
kCosmicRay
const Cut kCosmicRay
Definition: PAC2020/SBND/helper.h:54
RecoUtils::TrueParticleIDFromTotalRecoHits
int TrueParticleIDFromTotalRecoHits(detinfo::DetectorClocksData const &clockData, const std::vector< art::Ptr< recob::Hit > > &hits, bool rollup_unsaved_ids=1)
Definition: icaruscode/icaruscode/CRT/CRTUtils/RecoUtils.cc:83
TrackID
TrackID
Definition: ClusteringValidation_module.cc:57
icarus::AnalysisTreeDataStruct::GetMaxHits
size_t GetMaxHits() const
Returns the number of hits for which memory is allocated.
Definition: icaruscode/icaruscode/Analysis/AnalysisTree_module.cc:570
tca::dEdx
float dEdx(detinfo::DetectorClocksData const &clockData, detinfo::DetectorPropertiesData const &detProp, const TCSlice &slc, TP3D &tp3d)
Definition: PFPUtils.cxx:2687
icarus::AnalysisTreeDataStruct::exitZ
AuxDetMCData_t< Float_t > exitZ
Exit Z position of particle out of AuxDet.
Definition: icaruscode/icaruscode/Analysis/AnalysisTree_module.cc:496
icarus::AnalysisTreeDataStruct::ResizeCry
void ResizeCry(int nPrimaries)
Resize the data strutcure for Cry primaries.
Definition: icaruscode/icaruscode/Analysis/AnalysisTree_module.cc:1452
icarus::AnalysisTree::GetNTrackers
size_t GetNTrackers() const
Returns the number of trackers configured.
Definition: icaruscode/icaruscode/Analysis/AnalysisTree_module.cc:726
icarus::AnalysisTree::CreateTree
void CreateTree(bool bClearData=false)
Create the output tree and the data structures, if needed.
Definition: icaruscode/icaruscode/Analysis/AnalysisTree_module.cc:1718
icarus::AnalysisTreeDataStruct::exitY
AuxDetMCData_t< Float_t > exitY
Exit Y position of particle out of AuxDet.
Definition: icaruscode/icaruscode/Analysis/AnalysisTree_module.cc:495
icarus::AnalysisTreeDataStruct::AuxDetID
AuxDetMCData_t< Short_t > AuxDetID
Which AuxDet this particle went through.
Definition: icaruscode/icaruscode/Analysis/AnalysisTree_module.cc:489
icarus::AnalysisTree::fSaveCaloCosmics
bool fSaveCaloCosmics
save calorimetry information for cosmics
Definition: icaruscode/icaruscode/Analysis/AnalysisTree_module.cc:723
recob::Track::EndDirection
Vector_t EndDirection() const
Definition: bj/lardataobj/RecoBase/Track.h:133
sim::AuxDetIDE
MC truth information to make RawDigits and do back tracking.
Definition: AuxDetSimChannel.h:27
e
do i e
Definition: updateVersionList.sh:68
icarus::AnalysisTree::fSaveGenieInfo
bool fSaveGenieInfo
whether to extract and save CRY particle data
Definition: icaruscode/icaruscode/Analysis/AnalysisTree_module.cc:713
recob::Track::End
Point_t const & End() const
Definition: bj/lardataobj/RecoBase/Track.h:125
tca::TCEvent::subRun
unsigned int subRun
Definition: DataStructs.h:636
sim::SimChannel::TDCIDEMap
TDCIDEs_t const & TDCIDEMap() const
Returns all the deposited energy information as stored.
Definition: SimChannel.h:334
RecoUtils::TrueParticleIDFromTotalRecoCharge
int TrueParticleIDFromTotalRecoCharge(detinfo::DetectorClocksData const &clockData, const std::vector< art::Ptr< recob::Hit > > &hits, bool rollup_unsaved_ids=1)
Definition: icaruscode/icaruscode/CRT/CRTUtils/RecoUtils.cc:60
icarus::AnalysisTree::fSaveGeantInfo
bool fSaveGeantInfo
whether to extract and save Genie information
Definition: icaruscode/icaruscode/Analysis/AnalysisTree_module.cc:714
tca::evt
TCEvent evt
Definition: DataStructs.cxx:8
icarus::AnalysisTreeDataStruct::exitPy
AuxDetMCData_t< Float_t > exitPy
Exit y momentum of particle out of AuxDet.
Definition: icaruscode/icaruscode/Analysis/AnalysisTree_module.cc:499
recob::tracking::Point_t
ROOT::Math::PositionVector3D< ROOT::Math::Cartesian3D< Coord_t >, ROOT::Math::GlobalCoordinateSystemTag > Point_t
Type for representation of position in physical 3D space. See recob::tracking::Coord_t for more detai...
Definition: TrackingTypes.h:26
icarus::AnalysisTree::CreateData
void CreateData(bool bClearData=false)
Creates the structure for the tree data; optionally initializes it.
Definition: icaruscode/icaruscode/Analysis/AnalysisTree_module.cc:729
icarus::AnalysisTree::fSaveVertexInfo
bool fSaveVertexInfo
whether to extract and save Track information
Definition: icaruscode/icaruscode/Analysis/AnalysisTree_module.cc:717
recob::Track
Track from a non-cascading particle.A recob::Track consists of a recob::TrackTrajectory, plus additional members relevant for a &quot;fitted&quot; track:
Definition: bj/lardataobj/RecoBase/Track.h:49
anab::sParticleIDAlgScores::fAssumedPdg
int fAssumedPdg
PDG of particle hypothesis assumed by algorithm, if applicable. Set to 0 by default.
Definition: ParticleID.h:27
kMaxAuxDets
constexpr unsigned short kMaxAuxDets
max number of auxiliary detector cells per MC particle
Definition: icaruscode/icaruscode/Analysis/AnalysisTree_module.cc:145
double icarus::AnalysisTree::bdist ( const recob::tracking::Point_t &  pos)
private

Definition at line 2770 of file icaruscode/icaruscode/Analysis/AnalysisTree_module.cc.

2771 {
2772  // Get geometry.
2773  art::ServiceHandle<geo::Geometry> geom;
2774 
2775  geo::CryostatGeo const& cryo0 = geom->Cryostat(0);
2776  geo::CryostatGeo const& cryo1 = geom->Cryostat(1);
2777 
2778  geo::TPCGeo const& tpc00 = cryo0.TPC(0);
2779  TVector3 xyzcenter00 = tpc00.GetActiveVolumeCenter();
2780  std::cout << xyzcenter00[0] << " " << xyzcenter00[1] << " " << xyzcenter00[2] << std::endl;
2781 
2782  geo::TPCGeo const& tpc01 = cryo0.TPC(1);
2783  TVector3 xyzcenter01 = tpc01.GetActiveVolumeCenter();
2784  std::cout << xyzcenter01[0] << " " << xyzcenter01[1] << " " << xyzcenter01[2] << std::endl;
2785 
2786  geo::TPCGeo const& tpc10 = cryo1.TPC(0);
2787  TVector3 xyzcenter10 = tpc10.GetActiveVolumeCenter();
2788  std::cout << xyzcenter10[0] << " " << xyzcenter10[1] << " " << xyzcenter10[2] << std::endl;
2789 
2790  geo::TPCGeo const& tpc11 = cryo1.TPC(1);
2791  TVector3 xyzcenter11 = tpc11.GetActiveVolumeCenter();
2792  std::cout << xyzcenter11[0] << " " << xyzcenter11[1] << " " << xyzcenter11[2] << std::endl;
2793 
2794  double h00=tpc00.ActiveHalfHeight();
2795  double w00=tpc00.ActiveHalfWidth();
2796  double l00=tpc00.ActiveLength();
2797  std::cout << h00 << " " << w00 << " " << l00 << std::endl;
2798 
2799 
2800  double xmin = xyzcenter10[0]-w00;
2801  double xmax = xyzcenter11[0]+w00;
2802  double ymin = xyzcenter00[1]-h00;
2803  double ymax = xyzcenter00[1]+h00;
2804  double zmin = xyzcenter00[2]-l00/2;
2805  double zmax = xyzcenter00[2]+l00/2;
2806 
2807  double d1; // Distance to right side (wires).
2808  double d2; // Distance to left side (cathode).
2809 
2810  if (pos.X()>0) {
2811  d1=pos.X()-xmin;
2812  d2=xmax-pos.X();
2813  }
2814  if (pos.X()<0) {
2815  d1=fabs(pos.X())-xmin;
2816  d2=xmax-fabs(pos.X());
2817  }
2818 
2819  double d3 = pos.Y() -ymin; // Distance to bottom.
2820  double d4 = ymax - pos.Y(); // Distance to top.
2821  double d5 = pos.Z()-zmin; // Distance to front.
2822  double d6 = zmax - pos.Z(); // Distance to back.
2823 
2824  double result = std::min(std::min(std::min(std::min(std::min(d1, d2), d3), d4), d5), d6);
2825  return result;
2826 }
geo::TPCGeo::GetActiveVolumeCenter
Point GetActiveVolumeCenter() const
Returns the center of the TPC active volume in world coordinates [cm].
Definition: TPCGeo.h:288
geo::TPCGeo::ActiveHalfHeight
double ActiveHalfHeight() const
Half height (associated with y coordinate) of active TPC volume [cm].
Definition: TPCGeo.h:99
geo::TPCGeo
Geometry information for a single TPC.
Definition: TPCGeo.h:38
geo::CryostatGeo
Geometry information for a single cryostat.
Definition: CryostatGeo.h:43
xmin
process_name pandoraGausCryo1 vertexChargeCryo1 vertexStubCryo1 xmin
Definition: run_vertexcharge_icarus.fcl:109
geo::TPCGeo::ActiveHalfWidth
double ActiveHalfWidth() const
Half width (associated with x coordinate) of active TPC volume [cm].
Definition: TPCGeo.h:95
geo::TPCGeo::ActiveLength
double ActiveLength() const
Length (associated with z coordinate) of active TPC volume [cm].
Definition: TPCGeo.h:103
geo::CryostatGeo::TPC
const TPCGeo & TPC(unsigned int itpc) const
Return the itpc&#39;th TPC in the cryostat.
Definition: CryostatGeo.cxx:93
cout
BEGIN_PROLOG could also be cout
Definition: messageservice.fcl:10
void icarus::AnalysisTree::beginSubRun ( const art::SubRun &  sr)

Definition at line 1728 of file icaruscode/icaruscode/Analysis/AnalysisTree_module.cc.

1729 {
1730 
1731  art::Handle< sumdata::POTSummary > potListHandle;
1732  //sr.getByLabel(fPOTModuleLabel,potListHandle);
1733 
1734  if(sr.getByLabel(fPOTModuleLabel,potListHandle))
1735  SubRunData.pot=potListHandle->totpot;
1736  else
1737  SubRunData.pot=0.;
1738 
1739 }
icarus::AnalysisTreeDataStruct::SubRunData_t::pot
Double_t pot
Definition: icaruscode/icaruscode/Analysis/AnalysisTree_module.cc:339
void icarus::AnalysisTree::CheckData ( std::string  caller) const
inlineprivate

Helper function: throws if no data structure is available.

Definition at line 775 of file icaruscode/icaruscode/Analysis/AnalysisTree_module.cc.

776  {
777  if (fData) return;
778  throw art::Exception(art::errors::LogicError)
779  << "AnalysisTree::" << caller << ": no data";
780  } // CheckData()
void icarus::AnalysisTree::CheckTree ( std::string  caller) const
inlineprivate

Helper function: throws if no tree is available.

Definition at line 782 of file icaruscode/icaruscode/Analysis/AnalysisTree_module.cc.

783  {
784  if (fTree) return;
785  throw art::Exception(art::errors::LogicError)
786  << "AnalysisTree::" << caller << ": no tree";
787  } // CheckData()
void icarus::AnalysisTree::CreateData ( bool  bClearData = false)
inlineprivate

Creates the structure for the tree data; optionally initializes it.

Definition at line 729 of file icaruscode/icaruscode/Analysis/AnalysisTree_module.cc.

730  {
731  if (!fData) {
732  fData = new AnalysisTreeDataStruct(GetNTrackers());
733  fData->SetBits(AnalysisTreeDataStruct::tdAuxDet, !fSaveAuxDetInfo);
734  fData->SetBits(AnalysisTreeDataStruct::tdCry, !fSaveCryInfo);
735  fData->SetBits(AnalysisTreeDataStruct::tdGenie, !fSaveGenieInfo);
736  fData->SetBits(AnalysisTreeDataStruct::tdGeant, !fSaveGeantInfo);
737  }
738  else {
739  fData->SetBits(AnalysisTreeDataStruct::tdHit, !fSaveHitInfo);
740  fData->SetBits(AnalysisTreeDataStruct::tdTrack, !fSaveTrackInfo);
741  fData->SetBits(AnalysisTreeDataStruct::tdVtx, !fSaveVertexInfo);
742  fData->SetTrackers(GetNTrackers());
743  if (bClearData) fData->Clear();
744  }
745  } // CreateData()
icarus::AnalysisTree::fSaveTrackInfo
bool fSaveTrackInfo
whether to extract and save Hit information
Definition: icaruscode/icaruscode/Analysis/AnalysisTree_module.cc:716
icarus::AnalysisTree::fSaveAuxDetInfo
bool fSaveAuxDetInfo
whether to extract and save auxiliary detector data
Definition: icaruscode/icaruscode/Analysis/AnalysisTree_module.cc:711
icarus::AnalysisTreeDataStruct::SetBits
void SetBits(unsigned int setbits, bool unset=false)
Sets the specified bits.
Definition: icaruscode/icaruscode/Analysis/AnalysisTree_module.cc:527
icarus::AnalysisTreeDataStruct::SetTrackers
void SetTrackers(size_t nTrackers)
Allocates data structures for the given number of trackers (no Clear())
Definition: icaruscode/icaruscode/Analysis/AnalysisTree_module.cc:548
icarus::AnalysisTree::fSaveHitInfo
bool fSaveHitInfo
whether to extract and save Geant information
Definition: icaruscode/icaruscode/Analysis/AnalysisTree_module.cc:715
icarus::AnalysisTree::GetNTrackers
size_t GetNTrackers() const
Returns the number of trackers configured.
Definition: icaruscode/icaruscode/Analysis/AnalysisTree_module.cc:726
icarus::AnalysisTree::fSaveGenieInfo
bool fSaveGenieInfo
whether to extract and save CRY particle data
Definition: icaruscode/icaruscode/Analysis/AnalysisTree_module.cc:713
icarus::AnalysisTree::fSaveGeantInfo
bool fSaveGeantInfo
whether to extract and save Genie information
Definition: icaruscode/icaruscode/Analysis/AnalysisTree_module.cc:714
icarus::AnalysisTree::fSaveVertexInfo
bool fSaveVertexInfo
whether to extract and save Track information
Definition: icaruscode/icaruscode/Analysis/AnalysisTree_module.cc:717
void icarus::AnalysisTree::CreateTree ( bool  bClearData = false)
private

Create the output tree and the data structures, if needed.

Definition at line 1718 of file icaruscode/icaruscode/Analysis/AnalysisTree_module.cc.

1718  {
1719  if (!fTree) {
1720  art::ServiceHandle<art::TFileService> tfs;
1721  fTree = tfs->make<TTree>("anatree","analysis tree");
1722  }
1723  CreateData(bClearData);
1724  SetAddresses();
1725 } // icarus::AnalysisTree::CreateTree()
icarus::AnalysisTree::SetAddresses
void SetAddresses()
Sets the addresses of all the tree branches, creating the missing ones.
Definition: icaruscode/icaruscode/Analysis/AnalysisTree_module.cc:748
tfs
art::ServiceHandle< art::TFileService > tfs
Definition: NuVertexChargeTree_module.cc:254
icarus::AnalysisTree::CreateData
void CreateData(bool bClearData=false)
Creates the structure for the tree data; optionally initializes it.
Definition: icaruscode/icaruscode/Analysis/AnalysisTree_module.cc:729
void icarus::AnalysisTree::DestroyData ( )
inlineprivate

Destroy the local buffers (existing branches will point to invalid address!)

Definition at line 772 of file icaruscode/icaruscode/Analysis/AnalysisTree_module.cc.

772 { if (fData) { delete fData; fData = nullptr; } }
size_t icarus::AnalysisTree::GetNTrackers ( ) const
inlineprivate

Returns the number of trackers configured.

Definition at line 726 of file icaruscode/icaruscode/Analysis/AnalysisTree_module.cc.

726 { return fTrackModuleLabel.size(); }
void icarus::AnalysisTree::HitsPurity ( detinfo::DetectorClocksData const &  clockData,
std::vector< art::Ptr< recob::Hit > > const &  hits,
Int_t &  trackid,
Float_t &  purity,
double &  maxe 
)
private

Definition at line 2729 of file icaruscode/icaruscode/Analysis/AnalysisTree_module.cc.

2729  {
2730 
2731  trackid = -1;
2732  purity = -1;
2733 
2734  art::ServiceHandle<cheat::BackTrackerService> backTracker;
2735 
2736  std::map<int,double> trkide;
2737 
2738  for(size_t h = 0; h < hits.size(); ++h){
2739 
2740  art::Ptr<recob::Hit> hit = hits[h];
2741  std::vector<sim::IDE> ides;
2742 
2743  std::vector<sim::TrackIDE> eveIDs = backTracker->HitToEveTrackIDEs(clockData, hit);
2744 
2745  for(size_t e = 0; e < eveIDs.size(); ++e){
2746  //std::cout<<h<<" "<<e<<" "<<eveIDs[e].trackID<<" "<<eveIDs[e].energy<<" "<<eveIDs[e].energyFrac<<std::endl;
2747  trkide[eveIDs[e].trackID] += eveIDs[e].energy;
2748  }
2749  }
2750 
2751  maxe = -1;
2752  double tote = 0;
2753  for (std::map<int,double>::iterator ii = trkide.begin(); ii!=trkide.end(); ++ii){
2754  tote += ii->second;
2755  if ((ii->second)>maxe){
2756  maxe = ii->second;
2757  trackid = ii->first;
2758  }
2759  }
2760 
2761  //std::cout << "the total energy of this reco track is: " << tote << std::endl;
2762 
2763 
2764  if (tote>0){
2765  purity = maxe/tote;
2766  }
2767 }
hit
process_name hit
Definition: cheaterreco.fcl:51
h
while getopts h
Definition: getStatistics.sh:13
e
do i e
Definition: updateVersionList.sh:68
double icarus::AnalysisTree::length ( const recob::Track &  track)
private

Definition at line 2829 of file icaruscode/icaruscode/Analysis/AnalysisTree_module.cc.

2830 {
2831  return track.Length();
2832 }
recob::Track::Length
double Length(size_t p=0) const
Access to various track properties.
Definition: bj/lardataobj/RecoBase/Track.h:167
double icarus::AnalysisTree::length ( const simb::MCParticle &  part,
TVector3 &  start,
TVector3 &  end 
)
private

Definition at line 2835 of file icaruscode/icaruscode/Analysis/AnalysisTree_module.cc.

2836 {
2837  // Get geometry.
2838  art::ServiceHandle<geo::Geometry> geom;
2839  // auto const* detprop = lar::providerFrom<detinfo::DetectorPropertiesService>();
2840 
2841  // Get active volume boundary.
2842  //double xmin = 0.;
2843 
2844 
2845  geo::CryostatGeo const& cryo0 = geom->Cryostat(0);
2846  geo::CryostatGeo const& cryo1 = geom->Cryostat(1);
2847 
2848  geo::TPCGeo const& tpc00 = cryo0.TPC(0);
2849  TVector3 xyzcenter00 = tpc00.GetActiveVolumeCenter();
2850  std::cout << xyzcenter00[0] << " " << xyzcenter00[1] << " " << xyzcenter00[2] << std::endl;
2851 
2852  geo::TPCGeo const& tpc01 = cryo0.TPC(1);
2853  TVector3 xyzcenter01 = tpc01.GetActiveVolumeCenter();
2854  std::cout << xyzcenter01[0] << " " << xyzcenter01[1] << " " << xyzcenter01[2] << std::endl;
2855 
2856  geo::TPCGeo const& tpc10 = cryo1.TPC(0);
2857  TVector3 xyzcenter10 = tpc10.GetActiveVolumeCenter();
2858  std::cout << xyzcenter10[0] << " " << xyzcenter10[1] << " " << xyzcenter10[2] << std::endl;
2859 
2860  geo::TPCGeo const& tpc11 = cryo1.TPC(1);
2861  TVector3 xyzcenter11 = tpc11.GetActiveVolumeCenter();
2862  std::cout << xyzcenter11[0] << " " << xyzcenter11[1] << " " << xyzcenter11[2] << std::endl;
2863 
2864  double h00=tpc00.ActiveHalfHeight();
2865  double w00=tpc00.ActiveHalfWidth();
2866  double l00=tpc00.ActiveLength();
2867  std::cout << h00 << " " << w00 << " " << l00 << std::endl;
2868 
2869 
2870  double xmin = xyzcenter10[0]-w00;
2871  double xmax = xyzcenter11[0]+w00;
2872  double ymin = xyzcenter00[1]-h00;
2873  double ymax = xyzcenter00[1]+h00;
2874  double zmin = xyzcenter00[2]-l00/2;
2875  double zmax = xyzcenter00[2]+l00/2;
2876  //double vDrift = 160*pow(10,-6);
2877 
2878  std::cout << "DET DIMENSIONS: xmin = " << xmin << " xmax = " << xmax << " ymin = " << ymin << " ymax = " << ymax << " zmin = " << zmin << " zmax = " << zmax << std::endl;
2879 
2880  double result = 0.;
2881  TVector3 disp;
2882  int n = part.NumberTrajectoryPoints();
2883  bool first = true;
2884 
2885  for(int i = 0; i < n; ++i) {
2886  // check if the particle is inside a TPC
2887  double mypos[3] = {part.Vx(i), part.Vy(i), part.Vz(i)};
2888  if (fabs(mypos[0]) >= xmin && fabs(mypos[0]) <= xmax && mypos[1] >= ymin && mypos[1] <= ymax && mypos[2] >= zmin && mypos[2] <= zmax){
2889  //if (mypos[0] >= -300.0 && mypos[0] <= 200.0 && mypos[1] >= ymin && mypos[1] <= ymax && mypos[2] >= zmin && mypos[2] <= zmax){
2890  double xGen = part.Vx(i);
2891  //double tGen = part.T(i);
2892  // Doing some manual shifting to account for
2893  // an interaction not occuring with the beam dump
2894  // we will reconstruct an x distance different from
2895  // where the particle actually passed to to the time
2896  // being different from in-spill interactions
2897  //double newX = xGen+(tGen*vDrift);
2898  double newX = xGen;
2899  //if (newX < -xmax || newX > (2*xmax)) continue;
2900  //if (newX < (2.*xmin) || newX > (2*xmax)) continue;
2901 
2902  TVector3 pos(newX,part.Vy(i),part.Vz(i));
2903  if(first){
2904  start = pos;
2905  }
2906  else {
2907  disp -= pos;
2908  result += disp.Mag();
2909  }
2910  first = false;
2911  disp = pos;
2912  end = pos;
2913  }
2914  }
2915  return result;
2916 }
geo::TPCGeo::GetActiveVolumeCenter
Point GetActiveVolumeCenter() const
Returns the center of the TPC active volume in world coordinates [cm].
Definition: TPCGeo.h:288
geo::TPCGeo::ActiveHalfHeight
double ActiveHalfHeight() const
Half height (associated with y coordinate) of active TPC volume [cm].
Definition: TPCGeo.h:99
geo::TPCGeo
Geometry information for a single TPC.
Definition: TPCGeo.h:38
geo::CryostatGeo
Geometry information for a single cryostat.
Definition: CryostatGeo.h:43
xmin
process_name pandoraGausCryo1 vertexChargeCryo1 vertexStubCryo1 xmin
Definition: run_vertexcharge_icarus.fcl:109
geo::TPCGeo::ActiveHalfWidth
double ActiveHalfWidth() const
Half width (associated with x coordinate) of active TPC volume [cm].
Definition: TPCGeo.h:95
icarus::ns::util::end
auto end(FixedBins< T, C > const &) noexcept
Definition: FixedBins.h:585
geo::TPCGeo::ActiveLength
double ActiveLength() const
Length (associated with z coordinate) of active TPC volume [cm].
Definition: TPCGeo.h:103
geo::CryostatGeo::TPC
const TPCGeo & TPC(unsigned int itpc) const
Return the itpc&#39;th TPC in the cryostat.
Definition: CryostatGeo.cxx:93
cout
BEGIN_PROLOG could also be cout
Definition: messageservice.fcl:10
void icarus::AnalysisTree::SetAddresses ( )
inlineprivate

Sets the addresses of all the tree branches, creating the missing ones.

Definition at line 748 of file icaruscode/icaruscode/Analysis/AnalysisTree_module.cc.

749  {
750  CheckData("SetAddress()"); CheckTree("SetAddress()");
751  fData->SetAddresses(fTree, fTrackModuleLabel, isCosmics);
752  } // SetAddresses()
icarus::AnalysisTreeDataStruct::SetAddresses
void SetAddresses(TTree *pTree, const std::vector< std::string > &trackers, bool isCosmics)
Connect this object with a tree.
Definition: icaruscode/icaruscode/Analysis/AnalysisTree_module.cc:1471
icarus::AnalysisTree::CheckTree
void CheckTree(std::string caller) const
Helper function: throws if no tree is available.
Definition: icaruscode/icaruscode/Analysis/AnalysisTree_module.cc:782
icarus::AnalysisTree::CheckData
void CheckData(std::string caller) const
Helper function: throws if no data structure is available.
Definition: icaruscode/icaruscode/Analysis/AnalysisTree_module.cc:775
void icarus::AnalysisTree::SetTrackerAddresses ( size_t  iTracker)
inlineprivate

Sets the addresses of all the tree branches of the specified tracking algo, creating the missing ones

Definition at line 756 of file icaruscode/icaruscode/Analysis/AnalysisTree_module.cc.

757  {
758  CheckData("SetTrackerAddresses()"); CheckTree("SetTrackerAddresses()");
759  if (iTracker >= fData->GetNTrackers()) {
760  throw art::Exception(art::errors::LogicError)
761  << "AnalysisTree::SetTrackerAddresses(): no tracker #" << iTracker
762  << " (" << fData->GetNTrackers() << " available)";
763  }
764  fData->GetTrackerData(iTracker) \
765  .SetAddresses(fTree, fTrackModuleLabel[iTracker], isCosmics);
766  } // SetTrackerAddresses()
icarus::AnalysisTreeDataStruct::TrackDataStruct::SetAddresses
void SetAddresses(TTree *pTree, std::string tracker, bool isCosmics)
Definition: icaruscode/icaruscode/Analysis/AnalysisTree_module.cc:974
icarus::AnalysisTree::CheckTree
void CheckTree(std::string caller) const
Helper function: throws if no tree is available.
Definition: icaruscode/icaruscode/Analysis/AnalysisTree_module.cc:782
icarus::AnalysisTree::CheckData
void CheckData(std::string caller) const
Helper function: throws if no data structure is available.
Definition: icaruscode/icaruscode/Analysis/AnalysisTree_module.cc:775
icarus::AnalysisTreeDataStruct::GetNTrackers
size_t GetNTrackers() const
Returns the number of trackers for which data structures are allocated.
Definition: icaruscode/icaruscode/Analysis/AnalysisTree_module.cc:567

Member Data Documentation

std::string icarus::AnalysisTree::fCalDataModuleLabel
private

Definition at line 701 of file icaruscode/icaruscode/Analysis/AnalysisTree_module.cc.

std::vector<std::string> icarus::AnalysisTree::fCalorimetryModuleLabel
private

Definition at line 707 of file icaruscode/icaruscode/Analysis/AnalysisTree_module.cc.

std::vector<std::string> icarus::AnalysisTree::fCosmicTaggerAssocLabel
private

whether to extract and save Vertex information

Definition at line 719 of file icaruscode/icaruscode/Analysis/AnalysisTree_module.cc.

std::string icarus::AnalysisTree::fCryGenModuleLabel
private

Definition at line 703 of file icaruscode/icaruscode/Analysis/AnalysisTree_module.cc.

AnalysisTreeDataStruct* icarus::AnalysisTree::fData
private

Definition at line 694 of file icaruscode/icaruscode/Analysis/AnalysisTree_module.cc.

std::string icarus::AnalysisTree::fDigitModuleLabel
private

Definition at line 698 of file icaruscode/icaruscode/Analysis/AnalysisTree_module.cc.

std::vector<std::string> icarus::AnalysisTree::fFlashMatchAssocLabel
private

Definition at line 720 of file icaruscode/icaruscode/Analysis/AnalysisTree_module.cc.

float icarus::AnalysisTree::fG4minE
private

Energy threshold to save g4 particle info

Definition at line 724 of file icaruscode/icaruscode/Analysis/AnalysisTree_module.cc.

std::string icarus::AnalysisTree::fG4ModuleLabel
private

Definition at line 704 of file icaruscode/icaruscode/Analysis/AnalysisTree_module.cc.

std::string icarus::AnalysisTree::fGenieGenModuleLabel
private

Definition at line 702 of file icaruscode/icaruscode/Analysis/AnalysisTree_module.cc.

std::string icarus::AnalysisTree::fHitsModuleLabel
private

Definition at line 699 of file icaruscode/icaruscode/Analysis/AnalysisTree_module.cc.

std::string icarus::AnalysisTree::fLArG4ModuleLabel
private

Definition at line 700 of file icaruscode/icaruscode/Analysis/AnalysisTree_module.cc.

std::vector<std::string> icarus::AnalysisTree::fParticleIDModuleLabel
private

Definition at line 708 of file icaruscode/icaruscode/Analysis/AnalysisTree_module.cc.

std::string icarus::AnalysisTree::fPOTModuleLabel
private

Definition at line 709 of file icaruscode/icaruscode/Analysis/AnalysisTree_module.cc.

bool icarus::AnalysisTree::fSaveAuxDetInfo
private

whether to extract and save auxiliary detector data

Definition at line 711 of file icaruscode/icaruscode/Analysis/AnalysisTree_module.cc.

bool icarus::AnalysisTree::fSaveCaloCosmics
private

save calorimetry information for cosmics

Definition at line 723 of file icaruscode/icaruscode/Analysis/AnalysisTree_module.cc.

bool icarus::AnalysisTree::fSaveCryInfo
private

Definition at line 712 of file icaruscode/icaruscode/Analysis/AnalysisTree_module.cc.

bool icarus::AnalysisTree::fSaveGeantInfo
private

whether to extract and save Genie information

Definition at line 714 of file icaruscode/icaruscode/Analysis/AnalysisTree_module.cc.

bool icarus::AnalysisTree::fSaveGenieInfo
private

whether to extract and save CRY particle data

Definition at line 713 of file icaruscode/icaruscode/Analysis/AnalysisTree_module.cc.

bool icarus::AnalysisTree::fSaveHitInfo
private

whether to extract and save Geant information

Definition at line 715 of file icaruscode/icaruscode/Analysis/AnalysisTree_module.cc.

bool icarus::AnalysisTree::fSaveTrackInfo
private

whether to extract and save Hit information

Definition at line 716 of file icaruscode/icaruscode/Analysis/AnalysisTree_module.cc.

bool icarus::AnalysisTree::fSaveVertexInfo
private

whether to extract and save Track information

Definition at line 717 of file icaruscode/icaruscode/Analysis/AnalysisTree_module.cc.

std::vector<std::string> icarus::AnalysisTree::fTrackModuleLabel
private

Definition at line 706 of file icaruscode/icaruscode/Analysis/AnalysisTree_module.cc.

TTree* icarus::AnalysisTree::fTree
private

Definition at line 689 of file icaruscode/icaruscode/Analysis/AnalysisTree_module.cc.

bool icarus::AnalysisTree::fUseBuffer
private

whether to use a permanent buffer (faster, huge memory)

Definition at line 710 of file icaruscode/icaruscode/Analysis/AnalysisTree_module.cc.

std::string icarus::AnalysisTree::fVertexModuleLabel
private

Definition at line 705 of file icaruscode/icaruscode/Analysis/AnalysisTree_module.cc.

bool icarus::AnalysisTree::isCosmics
private

if it contains cosmics

Definition at line 722 of file icaruscode/icaruscode/Analysis/AnalysisTree_module.cc.

AnalysisTreeDataStruct::SubRunData_t icarus::AnalysisTree::SubRunData
private

Definition at line 696 of file icaruscode/icaruscode/Analysis/AnalysisTree_module.cc.

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