Public Types | Public Member Functions | Static Public Member Functions | Protected Member Functions | Private Attributes | List of all members
cheat::MCCheater Class Reference
Inheritance diagram for cheat::MCCheater:
art::EDProducer art::ProducerBase art::Consumer art::EngineCreator art::ProductRegistryHelper

Public Types

using ModuleType = EDProducer
 
using WorkerType = WorkerT< EDProducer >
 
template<typename UserConfig , typename KeysToIgnore = void>
using Table = ProducerBase::Table< UserConfig, KeysToIgnore >
 

Public Member Functions

 MCCheater (fhicl::ParameterSet const &pset)
 
virtual ~MCCheater ()
 
void produce (art::Event &evt)
 
void beginJob ()
 
template<typename PROD , BranchType B = InEvent>
ProductID getProductID (std::string const &instanceName={}) const
 
template<typename PROD , BranchType B>
ProductID getProductID (ModuleDescription const &moduleDescription, std::string const &instanceName) const
 
bool modifiesEvent () const
 
template<typename T , BranchType = InEvent>
ProductToken< Tconsumes (InputTag const &)
 
template<typename T , art::BranchType BT>
art::ProductToken< Tconsumes (InputTag const &it)
 
template<typename T , BranchType = InEvent>
void consumesMany ()
 
template<typename Element , BranchType = InEvent>
ViewToken< Element > consumesView (InputTag const &)
 
template<typename T , art::BranchType BT>
art::ViewToken< TconsumesView (InputTag const &it)
 
template<typename T , BranchType = InEvent>
ProductToken< TmayConsume (InputTag const &)
 
template<typename T , art::BranchType BT>
art::ProductToken< TmayConsume (InputTag const &it)
 
template<typename T , BranchType = InEvent>
void mayConsumeMany ()
 
template<typename Element , BranchType = InEvent>
ViewToken< Element > mayConsumeView (InputTag const &)
 
template<typename T , art::BranchType BT>
art::ViewToken< TmayConsumeView (InputTag const &it)
 
base_engine_tcreateEngine (seed_t seed)
 
base_engine_tcreateEngine (seed_t seed, std::string const &kind_of_engine_to_make)
 
base_engine_tcreateEngine (seed_t seed, std::string const &kind_of_engine_to_make, label_t const &engine_label)
 
seed_t get_seed_value (fhicl::ParameterSet const &pset, char const key[]="seed", seed_t const implicit_seed=-1)
 

Static Public Member Functions

static cet::exempt_ptr< Consumernon_module_context ()
 

Protected Member Functions

CurrentProcessingContext const * currentContext () const
 
void validateConsumedProduct (BranchType const bt, ProductInfo const &pi)
 
void prepareForJob (fhicl::ParameterSet const &pset)
 
void showMissingConsumes () const
 

Private Attributes

int ntest
 
std::string fMCTruthListLabel
 
std::string fFLSHitListLabel
 
TH1F * fMuE
 
TH1F * fElE
 
TH1F * felL
 
TH1F * fmuL
 
TH1F * ftotalE
 
TH1F * ftotEXV
 
TH1F * ftotEYV
 
TH1F * fdcosX
 
TH1F * fdcosY
 
TH1F * fdcosZ
 
TH1F * fsimhitE
 
TH1F * fFLSE
 
TH1F * fnumsh
 
TH1F * fNumFLS
 
TH1F * ftlen
 
TH1F * fthreshE
 
TH1F * fnplanes
 
TH1F * fparte
 
TH1F * fnpart
 
TH1F * fntrajpts
 
TH1F * fcorx
 
TH1F * fcory
 
TH1F * fclhit
 
TH1F * fzstart
 

Detailed Description

Definition at line 53 of file MCCheater_module.cc.

Member Typedef Documentation

using art::EDProducer::ModuleType = EDProducer
inherited

Definition at line 34 of file EDProducer.h.

template<typename UserConfig , typename KeysToIgnore = void>
using art::EDProducer::Table = ProducerBase::Table<UserConfig, KeysToIgnore>
inherited

Definition at line 43 of file EDProducer.h.

using art::EDProducer::WorkerType = WorkerT<EDProducer>
inherited

Definition at line 35 of file EDProducer.h.

Constructor & Destructor Documentation

cheat::MCCheater::MCCheater ( fhicl::ParameterSet const &  pset)
explicit

Definition at line 122 of file MCCheater_module.cc.

122  :
123  fMCTruthListLabel (pset.get< std::string >("MCTruthListLabel" )),
124  fFLSHitListLabel (pset.get< std::string >("FLSHitListLabel" ))
125  {
126  // tell module what its making
127  produces< std::vector<cheat::SimHit> >();
128  produces< std::vector<cheat::SimTrack> >();
129  // produces< std::vector<cheat::SimVertex> >();
130  // produces< std::vector<cheat::SimCluster> >();
131 
132  // Be noisy to demonstrate what's happening
133  mf::LogInfo("MCCheater") << " MCCheater::MCCheater()\n";
134 
135  // print this to keep track of version number
136  mf::LogInfo("MCCheater") << " ****************************************\n"
137  << " * MCCheater Version Number 4 *\n"
138  << " ****************************************\n";
139 
140  }
MaybeLogger_< ELseverityLevel::ELsev_info, false > LogInfo
std::string fMCTruthListLabel
std::string fFLSHitListLabel
enum BeamMode string
cheat::MCCheater::~MCCheater ( )
virtual

Definition at line 143 of file MCCheater_module.cc.

144  {
145  //======================================================================
146  // Clean up any memory allocated by your module
147  //======================================================================
148  }

Member Function Documentation

void cheat::MCCheater::beginJob ( )
virtual

Reimplemented from art::EDProducer.

Definition at line 151 of file MCCheater_module.cc.

References fclhit, fcorx, fcory, fdcosX, fdcosY, fdcosZ, fElE, felL, fFLSE, fMuE, fmuL, fnpart, fnplanes, fntrajpts, fNumFLS, fnumsh, fparte, fsimhitE, fthreshE, ftlen, ftotalE, ftotEXV, ftotEYV, fzstart, and art::TFileDirectory::make().

152  {
153 
154  //
155  // Book histograms, ntuples, initialize counts etc., etc., ...
156  //
158 
159  fMuE = tfs->make<TH1F>("fMuE", ";True Muon Energy (GeV)", 50, 0., 10.);
160  fElE = tfs->make<TH1F>("fElE", ";True Electron Energy (GeV)", 50, 0., 10.);
161  felL = tfs->make<TH1F>("felL", ";length(cm)", 60, 0., 1200.);
162  fmuL = tfs->make<TH1F>("fmuL", ";length(cm)", 60, 0., 1200.);
163  ftotalE = tfs->make<TH1F>("ftotalE", ";E of SimHit (MeV)", 100, 0., 50.);
164  ftotEXV = tfs->make<TH1F>("ftotEXV", ";E of Xview SimHit (MeV)", 100, 0., 50.);
165  ftotEYV = tfs->make<TH1F>("ftotEYV", ";E of Yview SimHit (MeV)", 100, 0., 50.);
166  fdcosX = tfs->make<TH1F>("fdcosX", ";Dir cos X", 100, -1., 1.);
167  fdcosY = tfs->make<TH1F>("fdcosY", ";Dir cos Y", 100, -1., 1.);
168  fdcosZ = tfs->make<TH1F>("fdcosZ", ";Dir cos Z", 100, -1., 1.);
169  fthreshE = tfs->make<TH1F>("fthreshE", ";E of threshold SimHit (MeV)", 100, 0., 50.);
170  fsimhitE = tfs->make<TH1F>("fsimhitE", ";Total SimHit Energy per Event (GeV)", 200, 0., 20.);
171  // fshlen = tfs->make<TH1F>("fshlen", "Shower Length (cm)", 250, 0., 1000.);
172  ftlen = tfs->make<TH1F>("ftlen", ";Track Length (cm)", 250, 0., 500.);
173  // fEpar = tfs->make<TH1F>("fEpar", ";SimHit E per Event Tr 0 (GeV)", 250, 0., 5.);
174  fFLSE = tfs->make<TH1F>("fFLSE", ";Total FLSHit Energy per Event (GeV)", 200, 0., 20.);
175  fnumsh = tfs->make<TH1F>("fnumsh", ";Number of SimHits per Event", 50, -0.5, 999.5);
176  fnplanes = tfs->make<TH1F>("fnplanes", ";Number of Planes in Track", 100, 0., 200.);
177  // fVertx = tfs->make<TH2F>("fVertx", ";x (cm);y (cm)", 400, -800., 800., 400, -800., 800.);
178  // fstsxy = tfs->make<TH2F>("fstsxy", ";SimTrack start y vs x", 150, -150., 150., 240, -220., 220.);
179  // fstexy = tfs->make<TH2F>("fstexy", ";SimTrack end y vs x", 150, -150., 150., 240, -220., 220.);
180  // fFidVert = tfs->make<TH2F>("fFidfVert", ";x (cm);y (cm)", 400, -800., 800., 400, -800., 800.);
181  // fzvtx = tfs->make<TH1F>("fzvtx", ";Zvtx (cm)", 350, 0., 7000.);
182  fNumFLS = tfs->make<TH1F>("fNumFLS", ";Number of FLSHits in Event", 100, -0.5, 1999.5);
183  fnpart = tfs->make<TH1F>("fnpart", ";Number of Particles in Event", 50, -0.5, 49.5);
184  fparte = tfs->make<TH1F>("fparte", ";Total Energy of Particles in Event (GeV)", 200, 0., 20.);
185  fntrajpts = tfs->make<TH1F>("fntrajpts", ";Number of Trajectory Points per Particle", 50, 0., 50.);
186  // fXRange = tfs->make<TH1F>("fXRange", ";Range of SimHit X Coords", 120, -210., 210.);
187  // fYRange = tfs->make<TH1F>("fYRange", ";Range of SimHit Y Coords", 120, -210., 210.);
188  // fZRange = tfs->make<TH1F>("fZRange", ";Range of SimHit Z Coords", 100, 0., 1000.);
189  // fnshptr = tfs->make<TH1F>("fnshptr", ";Number of SimHits per Track", 50, 0., 100.);
190  fcorx = tfs->make<TH1F>("fcorx", ";X Coord of xview plane", 120, -810., 810.);
191  fcory = tfs->make<TH1F>("fcory", ";Y Coord of yview plane", 120, -810., 810.);
192  fclhit = tfs->make<TH1F>("fclhit", ";Penetration Depth with SimHits", 125, 0., 250.);
193  fzstart = tfs->make<TH1F>("fzstart", ";Plane of particle start", 100, -1.0, 199.0);
194 
195  }
T * make(ARGS...args) const
template<typename T , BranchType = InEvent>
ProductToken<T> art::Consumer::consumes ( InputTag const &  )
inherited
template<typename T , art::BranchType BT>
art::ProductToken<T> art::Consumer::consumes ( InputTag const &  it)
inherited

Definition at line 146 of file Consumer.h.

References art::InputTag::instance(), PandAna.reco_validation.prod5_pid_validation::invalid, art::InputTag::label(), art::InputTag::process(), and T.

147 {
148  if (!moduleContext_)
149  return ProductToken<T>::invalid();
150 
151  consumables_[BT].emplace_back(ConsumableType::Product,
152  TypeID{typeid(T)},
153  it.label(),
154  it.instance(),
155  it.process());
156  return ProductToken<T>{it};
157 }
set< int >::iterator it
static ProductToken< T > invalid()
Definition: ProductToken.h:47
ConsumableProducts consumables_
Definition: Consumer.h:137
double T
Definition: Xdiff_gwt.C:5
bool moduleContext_
Definition: Consumer.h:135
template<typename T , art::BranchType BT>
void art::Consumer::consumesMany ( )
inherited

Definition at line 161 of file Consumer.h.

References T.

162 {
163  if (!moduleContext_)
164  return;
165 
166  consumables_[BT].emplace_back(ConsumableType::Many, TypeID{typeid(T)});
167 }
ConsumableProducts consumables_
Definition: Consumer.h:137
double T
Definition: Xdiff_gwt.C:5
bool moduleContext_
Definition: Consumer.h:135
template<typename Element , BranchType = InEvent>
ViewToken<Element> art::Consumer::consumesView ( InputTag const &  )
inherited
template<typename T , art::BranchType BT>
art::ViewToken<T> art::Consumer::consumesView ( InputTag const &  it)
inherited

Definition at line 171 of file Consumer.h.

References art::InputTag::instance(), PandAna.reco_validation.prod5_pid_validation::invalid, art::InputTag::label(), art::InputTag::process(), and T.

172 {
173  if (!moduleContext_)
174  return ViewToken<T>::invalid();
175 
176  consumables_[BT].emplace_back(ConsumableType::ViewElement,
177  TypeID{typeid(T)},
178  it.label(),
179  it.instance(),
180  it.process());
181  return ViewToken<T>{it};
182 }
set< int >::iterator it
static ViewToken< Element > invalid()
Definition: ProductToken.h:75
ConsumableProducts consumables_
Definition: Consumer.h:137
double T
Definition: Xdiff_gwt.C:5
bool moduleContext_
Definition: Consumer.h:135
base_engine_t& art::EngineCreator::createEngine ( seed_t  seed)
inherited
base_engine_t& art::EngineCreator::createEngine ( seed_t  seed,
std::string const &  kind_of_engine_to_make 
)
inherited
base_engine_t& art::EngineCreator::createEngine ( seed_t  seed,
std::string const &  kind_of_engine_to_make,
label_t const &  engine_label 
)
inherited
CurrentProcessingContext const* art::EDProducer::currentContext ( ) const
protectedinherited
seed_t art::EngineCreator::get_seed_value ( fhicl::ParameterSet const &  pset,
char const  key[] = "seed",
seed_t const  implicit_seed = -1 
)
inherited
template<typename PROD , BranchType B>
ProductID art::EDProducer::getProductID ( std::string const &  instanceName = {}) const
inlineinherited

Definition at line 123 of file EDProducer.h.

References art::EDProducer::moduleDescription_.

Referenced by skim::NueSkimmer::CopyMichelSlice(), and skim::NueSkimmer::CopyMichelTrack().

124  {
125  return ProducerBase::getProductID<PROD, B>(moduleDescription_,
126  instanceName);
127  }
ModuleDescription moduleDescription_
Definition: EDProducer.h:115
template<typename PROD , BranchType B>
ProductID art::ProducerBase::getProductID ( ModuleDescription const &  moduleDescription,
std::string const &  instanceName 
) const
inherited

Definition at line 56 of file ProducerBase.h.

References art::ModuleDescription::moduleLabel().

Referenced by art::ProducerBase::modifiesEvent().

58  {
59  auto const& pd =
60  get_ProductDescription<PROD>(B, md.moduleLabel(), instanceName);
61  return pd.productID();
62  }
template<typename T , BranchType = InEvent>
ProductToken<T> art::Consumer::mayConsume ( InputTag const &  )
inherited
template<typename T , art::BranchType BT>
art::ProductToken<T> art::Consumer::mayConsume ( InputTag const &  it)
inherited

Definition at line 189 of file Consumer.h.

References art::InputTag::instance(), PandAna.reco_validation.prod5_pid_validation::invalid, art::InputTag::label(), art::InputTag::process(), and T.

190 {
191  if (!moduleContext_)
192  return ProductToken<T>::invalid();
193 
194  consumables_[BT].emplace_back(ConsumableType::Product,
195  TypeID{typeid(T)},
196  it.label(),
197  it.instance(),
198  it.process());
199  return ProductToken<T>{it};
200 }
set< int >::iterator it
static ProductToken< T > invalid()
Definition: ProductToken.h:47
ConsumableProducts consumables_
Definition: Consumer.h:137
double T
Definition: Xdiff_gwt.C:5
bool moduleContext_
Definition: Consumer.h:135
template<typename T , art::BranchType BT>
void art::Consumer::mayConsumeMany ( )
inherited

Definition at line 204 of file Consumer.h.

References T.

205 {
206  if (!moduleContext_)
207  return;
208 
209  consumables_[BT].emplace_back(ConsumableType::Many, TypeID{typeid(T)});
210 }
ConsumableProducts consumables_
Definition: Consumer.h:137
double T
Definition: Xdiff_gwt.C:5
bool moduleContext_
Definition: Consumer.h:135
template<typename Element , BranchType = InEvent>
ViewToken<Element> art::Consumer::mayConsumeView ( InputTag const &  )
inherited
template<typename T , art::BranchType BT>
art::ViewToken<T> art::Consumer::mayConsumeView ( InputTag const &  it)
inherited

Definition at line 214 of file Consumer.h.

References art::InputTag::instance(), PandAna.reco_validation.prod5_pid_validation::invalid, art::InputTag::label(), art::InputTag::process(), and T.

215 {
216  if (!moduleContext_)
217  return ViewToken<T>::invalid();
218 
219  consumables_[BT].emplace_back(ConsumableType::ViewElement,
220  TypeID{typeid(T)},
221  it.label(),
222  it.instance(),
223  it.process());
224  return ViewToken<T>{it};
225 }
set< int >::iterator it
static ViewToken< Element > invalid()
Definition: ProductToken.h:75
ConsumableProducts consumables_
Definition: Consumer.h:137
double T
Definition: Xdiff_gwt.C:5
bool moduleContext_
Definition: Consumer.h:135
bool art::ProducerBase::modifiesEvent ( ) const
inlineinherited

Definition at line 40 of file ProducerBase.h.

References art::ProducerBase::getProductID(), and string.

41  {
42  return true;
43  }
static cet::exempt_ptr<Consumer> art::Consumer::non_module_context ( )
staticinherited
void art::Consumer::prepareForJob ( fhicl::ParameterSet const &  pset)
protectedinherited
void cheat::MCCheater::produce ( art::Event evt)
virtual

Implements art::EDProducer.

Definition at line 197 of file MCCheater_module.cc.

References abs(), rb::Track::AppendTrajectoryPoint(), sim::ParticleNavigator::begin(), geo::PlaneGeo::Cell(), om::cout, DEFINE_ART_MODULE(), geo::GeometryBase::DetHalfHeight(), geo::GeometryBase::DetHalfWidth(), geo::GeometryBase::DetLength(), dir, dist, simb::MCParticle::E(), febshutoff_auto::end, sim::ParticleNavigator::end(), sim::ParticleNavigator::EveId(), fcorx, fcory, fElE, fFLSE, fFLSHitListLabel, sim::FLSHitList::fHits, fMuE, cheat::SimHit::fNP, fnpart, fntrajpts, fNumFLS, fnumsh, fparte, fsimhitE, fthreshE, ftotalE, cheat::SimHit::fTotE, ftotEXV, ftotEYV, cheat::SimHit::fTrID, art::DataViewImpl::getByLabel(), geo::CellGeo::GetCenter(), make_syst_table_plots::h, hits(), MECModelEnuComparisons::i, it, calib::j, geo::kX, geo::kY, geo::CellGeo::LocalToWorld(), simb::MCParticle::Mother(), ntest, simb::MCParticle::NumberTrajectoryPoints(), part, cheat::BackTracker::ParticleNavigator(), simb::MCParticle::PdgCode(), BlessedPlots::pid, geo::GeometryBase::Plane(), elec2geo::pos, simb::MCParticle::Position(), art::Event::put(), moon_position_table_new3::second, rb::CellHit::SetCell(), rawdata::RawDigit::SetChannel(), cheat::SimHit::SetCoorX(), cheat::SimHit::SetCoorY(), cheat::SimHit::SetCoorZ(), rb::Prong::SetDir(), cheat::SimHit::SetEnD(), cheat::SimHit::SetNP(), rb::CellHit::SetPE(), cheat::SimTrack::SetPID(), cheat::SimHit::SetPID(), rb::CellHit::SetPlane(), rb::Track::SetStart(), cheat::SimHit::SetTim(), cheat::SimHit::SetTotE(), cheat::SimTrack::SetTrackID(), cheat::SimHit::SetTrID(), rb::CellHit::SetView(), std::sqrt(), febshutoff_auto::start, strack, nd_projection_maker::tpp, simb::MCParticle::TrackId(), POTSpillRate::view, rb::CellHit::View(), geo::PlaneGeo::View(), X, Y, and Z.

198  {
199  //======================================================================
200  // Called for every event. "Reco" implies that you are adding
201  // information to the event
202  //======================================================================
203 
204  mf::LogInfo("MCCheater") << "MCCheater::produce() " << " New Reco Event\n";
205 
206  // increment counter for M matrix
207  ntest++;
208 
209  /*
210  // get MCTruth info - should loop over this list since there may be more
211  // than 1 interaction, a neutrino event and a cosmic event for example
212  // other than that, don't need this list for the products of the cheater
213  art::Handle< std::vector<simb::MCTruth> > mclist;
214  evt.getByLabel(fMCTruthListLabel, mclist);
215  if (mclist->empty())
216  {
217  mf::LogError("MCCheater") << "No MCTruth List - bummer\n";
218  return;
219  }
220 
221  // for now, only have 1 MCTruth per event, so only need mclist[0]
222  // art::PtrVector<simb::MCTruth> mctruth; don't need this either?
223  // for (unsigned int i=0; i<mclist->size(); ++i)
224  // {
225  art::Ptr<simb::MCTruth> mct(mclist, 0);
226  std::cout << " Number of Particles from MCTruth = " << mct->NParticles() << "\n";
227  // }
228 
229  mf::LogInfo("MCCheater") << "Number of MCTruth lists = "<< mclist->size() << "\n";
230 
231 
232  // fill a histogram with the size of each MCTruth
233  // for now, there is only 1 MCTruth list per MC event, so only need to use i=0
234 
235  for(unsigned int i=0; i < mclist->size(); ++i)
236  {
237  std::cout << " Number of particles from MCTruth = " << mctruth[i]->NParticles() << "\n";
238  }
239  */
240 
241  // PARTICLE LIST NEEDED FOR CHEATER
242  // get Particle info - needed for Sim Objects
244  const sim::ParticleNavigator &pnav = bt->ParticleNavigator();
245 
246  // First list needed for Cheater - the particle list
247  // using the particle list, make some plots,
248  // make a map linking trackID with particles
249  // and a map linking trackID with its particleID - this is all that's needed for cheater
250  std::map<int, const sim::Particle*> trpartmap; // dont really need this map
251  std::map<int, int> trpidmap; // TrackID ParticleID map
252  double PartE = 0.;
253  unsigned int NPart = 0;
254  int npar = 0;
255  int pcount = 0;
256  for (sim::ParticleNavigator::const_iterator j = pnav.begin(); j != pnav.end(); ++j){
257 
258  const sim::Particle p = *((*j).second);
259  NPart += 1;
260  npar += 1;
261 
262  int pid = p.PdgCode(); // what be this particle
263  // energy of particle
264  PartE += p.E();
265  unsigned int ntrajpts = p.NumberTrajectoryPoints();
266  fntrajpts->Fill(ntrajpts);
267  // for now, just check length of particle
268  // double plength = sqrt((enx-stx)*(enx-stx)+(eny-sty)*(eny-sty)+(enz-stz)*(enz-stz));
269  if (abs(pid)==13) fMuE->Fill(p.E()); // true muon energy for some reason?
270  if (abs(pid)==11) fElE->Fill(p.E()); // true electron energy
271  int tid = p.TrackId();
272  // reduce number of particles with some obvious? cuts, no neutrinos, pi0s, photons, nuclei
273  if (abs(pid)==12 || abs(pid)==14) continue; // no neutrinos is only cut so far
274  pcount++;
275  trpartmap[tid] = &p;
276  trpidmap[tid]=pid;
277  } // end loop over particles in this list
278  mf::LogInfo("MCCheater") << "Final particle count after cuts in particle list = " << pcount << "\n";
279  fparte->Fill(PartE);
280  fnpart->Fill(NPart);
281 
282  // mf::LogInfo("MCCheater") << "TrackID Particle Map size " << trpartmap.size() << "\n";
283  // mf::LogInfo("MCCheater") << "TrackID ParticleID Map size " << trpidmap.size() << "\n";
284 
285 
286  // FLSHITS ALSO NEEDED FOR CHEATER
287  // need FLSHits also - get FLSHit info
289  evt.getByLabel(fFLSHitListLabel, hitlist);
290  if (hitlist->empty()){
291  mf::LogError("MCCheater") << "No FLSHit Lists - bummer\n";
292  return;
293  }
294  // mf::LogInfo("MCCheater") << "Number of FLSHit hitlists = " << hitlist->size() << "\n";
295 
296  // Make maps if needed - these maps not necessary for the Cheater
297  // std::map<unsigned int, std::pair<int, int> > tidmgmmap; // map linking track ID to its mother and grandmother
298  std::map<unsigned int, std::vector<std::pair<unsigned int, unsigned int> > > tidpcpmap;
299  // std::map<std::pair<unsigned int, unsigned int>, std::vector<int> > pcptidmap; // map linking pair plane cell to vector of TrackIDs
300  // std::map<std::pair<unsigned int, unsigned int>, std::vector<int> > pcpmotmap; // map linking pair plane cell to vector of TrackID mothers
301  // std::map<std::pair<unsigned int, unsigned int>, std::vector<int> > pcppidmap; // map linking pair plane cell to vector of particleIDs
302  // std::map<std::pair<unsigned int, unsigned int>, std::vector<double> > pcpendmap; // map linking pair plane cell to vector of Edeps
303  // std::map<std::pair<unsigned int, unsigned int>, std::vector<double> > pcptimmap; // map linking pair plane cell to vector of times
304 
305  // FLSHits are formed for each particle that deposits energy in a cell, so there
306  // can be more than 1 FLSHit in a particular cell corresponding to each particle
307  // that contributed to the energy. And the cell IDs in FLSHits are not unique, when
308  // combined with the PlaneId, then a unique channel ID for the whole detector can be determined.
309  int nFLShits = 0;
310  double flshitE = 0.;
311  // geo::Geometry& geo = util::EDMUtils::GetGeometry(evt);
313  // Make a map linking FLSHits to TrackID
314  std::map<int, std::vector<const sim::FLSHit*> > trhitmap; // map linking TrackID to FLSHits
315  // can define a map with a pair of keys and a value, so can make a map with plane and cell as key and
316  // the value a vector of eg, FLSHits
317  std::map<std::pair<unsigned int, unsigned int>, std::vector<const sim::FLSHit*> > pcpflsmap; // map linking plane cell pair to FLShits
318 
319  for (unsigned int i=0; i<hitlist->size(); ++i) // loop over FLSHit Lists - 1 per process interaction
320  {
321 
322  art::Ptr<sim::FLSHitList> hlist(hitlist, i);
323  // mf::LogInfo("MCCheater") << "Num of FLSHits in list = " << hlist->fHits.size() << "\n";
324 
325  const std::vector<sim::FLSHit>& hits = hlist->fHits;
326 
327  nFLShits += hits.size();
328  // loop over FLSHits in this list
329  for (unsigned int j=0; j<hits.size(); ++j)
330  {
331  int tid = hits[j].GetTrackID();
332  if (tid<0) continue; // just protection - shouldn't be here this is incoming particle (neutrino)
333  // here, loop over particle list to get mother of this trackID
334  int tmoth = -2; // set mother to -2 for now
335  int tgmoth = -2; // set grandmother to -2 for now
336  int tpid = 0;
337  // need to use only the particle id from particle list for this track - FLSHit pid can be different one
338  for (sim::ParticleNavigator::const_iterator k = pnav.begin(); k != pnav.end(); ++k)
339  {
340  const sim::Particle part = *((*k).second);
341  if (part.TrackId()==tid)
342  {
343  tpid = part.PdgCode();
344  tmoth = part.Mother();
345  }
346  // have to find grandma myself -
347  tgmoth = pnav.EveId(part.TrackId());
348  }
349  int smoth = tmoth;
350  if (smoth > -2) smoth = 1;
351  if (tmoth==-2) mf::LogInfo("MCCheater") << " THIS TrackID NOT IN PARTICLE LIST = " << tid << " " << hits[i].GetPDG() <<"\n";
352  // make a pair out of tmoth and tgmoth
353  std::pair<int, int> mgmpair(tmoth, tgmoth);
354  // if (tpid>6000) continue;
355  unsigned int planeid = hits[j].GetPlaneID();
356  unsigned int cellid = hits[j].GetCellID();
357  // define the key as a plane cell pair here
358  std::pair<unsigned int, unsigned int> pcpair(planeid, cellid);
359  // this is pair of trackID and particleID
360  std::pair<unsigned int, unsigned int> tppair(tid, tpid);
361  double xyzloc[3];
362  double xyzwor[3];
363  xyzloc[0] = 0.5*(hits[j].GetEntryX()+hits[j].GetExitX());
364  xyzloc[1] = 0.5*(hits[j].GetEntryY()+hits[j].GetExitY());
365  xyzloc[2] = 0.5*(hits[j].GetEntryZ()+hits[j].GetExitZ());
366  fgeo->Plane(planeid)->Cell(cellid)->LocalToWorld(xyzloc,xyzwor);
367  // if (h.GetEdep()>0.001) trhitmap[tid].push_back(& h);
368  // no selection - fill maps with everything
369  trhitmap[tid].push_back(& hits[j]);
370  pcpflsmap[pcpair].push_back(& hits[j]);
371  double end = hits[j].GetEdep();
372  flshitE += end;
373  // this is time, but don't use it yet
374  // double tim = hits[j].fT;
375  // geo::CellUniqueId cid = hits[j].fId;
376  // fill maps which will be used to make SimHits, so don't make any cuts yet
377  // if tpid is 22, a photon, don't use
378 
379  if (end>0.000000000001)
380  {
381  tidpcpmap[tid].push_back(pcpair);
382  // tidmgmmap[tid].push_back(mgmpair);
383  // tidmgmmap[tid]=mgmpair;
384  // pcptidmap[pcpair].push_back(tid);
385  // pcpmotmap[pcpair].push_back(tmoth);
386  // pcppidmap[pcpair].push_back(tpid);
387  // pcpendmap[pcpair].push_back(end);
388  // pcptimmap[pcpair].push_back(tim);
389  }
390 
391  } // end of FLSHit hits loop
392  } // end of FLSHit lists loop
393  fFLSE->Fill(flshitE); // histogram of total FLSHit E
394  fNumFLS->Fill(nFLShits); // histogram of number of FLSHits
395  /*
396  std::cout << " TrackID FLSHit Map size " << trhitmap.size() << "\n";
397  std::cout << " TrackID PlaneCell pair Map size " << tidpcpmap.size() << "\n";
398  std::cout << " TrackID MotherGrandmother pair Map size = " << tidmgmmap.size() << "\n";
399  std::cout << " PlaneCell pair FLSHit Map size " << pcpflsmap.size() << "\n";
400  std::cout << " PlaneCell pair TrackID Map size " << pcptidmap.size() << "\n";
401  std::cout << " PlaneCell pair TID Mother Map size " << pcpmotmap.size() << "\n";
402  std::cout << " PlaneCell pair ParticleID Map size " << pcppidmap.size() << "\n";
403  std::cout << " PlaneCell pair Edep Map size " << pcpendmap.size() << "\n";
404  std::cout << " PlaneCell pair Tdep Map size " << pcptimmap.size() << "\n";
405  */
406  // test some maps
407 
408  // mf::LogInfo("MCCheater") << "Testing PlaneCell pair FLSHit map\n";
409  /*
410  for (std::map<std::pair<unsigned int, unsigned int>, std::vector<const sim::FLSHit*> >::iterator itm = pcpflsmap.begin(); itm != pcpflsmap.end(); ++itm)
411  {
412  std::pair<unsigned int, unsigned int> planecell = itm->first;
413  std::vector<const sim::FLSHit*> harr = itm->second;
414  std::cout <<" Number of FLSHits = " << harr.size() << " For Plane and Cell = " << planecell.first << " " << planecell.second << "\n";
415  }
416  */
417 
418  // NEED THIS MAP FOR SIMTRACKS
419  // find unique planecells for trackIDs, so can tell whether a track is longer than n cells
420  // make a map linking each trackID with its unique planecells and the size of the value
421  // planecell array is the number of unique planecells on this track
422  std::map<unsigned int, std::vector<std::pair<unsigned int, unsigned int> > > tidupcpmap;
423  for (std::map<unsigned int, std::vector<std::pair<unsigned int, unsigned int> > >::iterator itm = tidpcpmap.begin(); itm != tidpcpmap.end(); ++itm)
424  {
425  unsigned int trackid = itm->first;
426  std::vector<std::pair<unsigned int, unsigned int> > pcarray = itm->second;
427 
428  // eliminate duplicates, resize array with erase
429  // stl unique moves duplicates to the end of the array, doesn't eliminate
430  // erase deletes the items after the new location from unique
431  // this is part of that magic of C++ (R Rechenmacher)
432  std::sort(pcarray.begin(), pcarray.end());
433  std::vector<std::pair<unsigned int, unsigned int> >::iterator endloc;
434  endloc = std::unique(pcarray.begin(), pcarray.end());
435  pcarray.erase(endloc, pcarray.end());
436  for (unsigned int i=0; i<pcarray.size(); ++i)
437  {
438  tidupcpmap[trackid].push_back(pcarray[i]);
439  }
440  }
441 
442 
443  // MAKE SIMHITS HERE
444  // make SimHits using the pc pair, FLSHit map
445  double Etotev = 0.;
446  unsigned int NumSH = 0;
447  std::vector<cheat::SimHit> simhits;
448  for (std::map<std::pair<unsigned int, unsigned int>, std::vector<const sim::FLSHit*> >::iterator itm = pcpflsmap.begin(); itm != pcpflsmap.end(); ++itm)
449  {
450  // cheat::SimHit* hit = new cheat::SimHit();
452  std::pair<unsigned int, unsigned int> pcp = itm->first; // plane-cell pair of this hit
453  std::vector<const sim::FLSHit*> fh = itm->second; // get the FLSHits in this pc pair
454  unsigned int hsize = fh.size(); // number of FLSHits in this cell
455 
456  // can set plane cell and view for this SimHit now
457  hit.SetPlane(pcp.first); // set plane of this hit
458  unsigned int planen = pcp.first; // can use this later to set fiducial cut on fd
459  hit.SetCell(pcp.second); // set cell of this hit
460  unsigned int celln = pcp.second; // can use this later for fiducial cuts on fd
461  // next 2 lines set view for this plane
462  geo::View_t view = fgeo->Plane(pcp.first)->View();
463  hit.SetView(view); // set view of this hit
464  uint32_t ichan = 1; // have to find this somewhere?
465  hit.SetChannel(ichan);
466 
467  // make the arrays of E, T, PID, and TrID
468  std::vector<double> earr;
469  std::vector<double> carr;
470  std::vector<int> tarr;
471  std::vector<int> parr;
472 
473  std::vector<std::pair<int, int> > tparr; // trackID, PID pair
474  std::vector<std::pair<int, double> > tearr;
475  std::vector<std::pair<int, double> > tcarr;
476  double ene = 0.;
477  double coorx = 0.;
478  double coory = 0.;
479  double coorz = 0.;
480 
481  double_t xxyyzz[3];
482  fgeo->Plane(planen)->Cell(celln)->GetCenter(xxyyzz);
483  coorz = xxyyzz[2];
484  if (view==geo::kX)
485  {
486  coorx = xxyyzz[0];
487  fcorx->Fill(coorx);
488  // need to get y from flshit
489  } else if (view==geo::kY)
490  {
491  coory = xxyyzz[1];
492  fcory->Fill(coory);
493  // need to get x from flshit
494  }
495 
496  for (unsigned int i = 0; i<hsize; ++i)
497  {
498  earr.push_back(fh[i]->GetEdep()); //energy of this hit
499  std::pair<int, double> tep(fh[i]->GetTrackID(), fh[i]->GetEdep()); // pair of trackID, hitE
500  tearr.push_back(tep);
501  carr.push_back(fh[i]->GetEntryT()); // time of this hit
502  std::pair<int, double> ttp(fh[i]->GetTrackID(), fh[i]->GetEntryT()); // pair of trackID, time
503  tcarr.push_back(ttp);
504 
505  tarr.push_back(fh[i]->GetTrackID()); // trackID of this hit
506  // get particle ID from TrackID ParticleID map, not from FLSHit since its different!
507  for (std::map<int, int>::iterator itm = trpidmap.begin(); itm != trpidmap.end(); ++itm)
508  {
509  int trid = itm->first;
510  if (trid == fh[i]->GetTrackID())
511  {
512  parr.push_back(itm->second); // particle ID of this trackID
513  std::pair<int, int> tpp(fh[i]->GetTrackID(), itm->second); // pair of trackID, particleID
514  tparr.push_back(tpp);
515  }
516  }
517  ene += fh[i]->GetEdep(); // energy sum of all contributions to this hit
518 
519  // this uses first FLS contribution to get the non-view coordinate for a SimHit
520  // replace this with an array of the actual values for each contribution
521  if (i==0)
522  {
523  double loc0[3];
524  double wor0[3];
525  double zexo = fh[0]->GetExitZ();
526  double zeno = fh[0]->GetEntryZ();
527  loc0[2] = 0.5*(zexo+zeno);
528  double xexo = fh[0]->GetExitX();
529  double xeno = fh[0]->GetEntryX();
530  loc0[0] = 0.5*(xexo+xeno);
531  double yexo = fh[0]->GetExitY();
532  double yeno = fh[0]->GetEntryY();
533  loc0[1] = 0.5*(yexo+yeno);
534  // change from local to global
535  fgeo->Plane(planen)->Cell(celln)->LocalToWorld(loc0,wor0);
536  if (view==geo::kX) coory = wor0[1];
537  if (view==geo::kY) coorx = wor0[0];
538  }
539  }
540  hit.SetTotE(ene); // set E of this hit
541  ftotalE->Fill(ene * 1000.); // fill this hist also in units of MeV
542  hit.SetPE(ene*5000.); // took a guess to set the number of photoelectrons
543  // hit.SetPECorr(ene*5000.); // took another guess to seet the number of corrected photoelectrons
544  hit.SetCoorX(coorx);
545  hit.SetCoorY(coory);
546  hit.SetCoorZ(coorz);
547 
548  // now need to finish definition based on number of FLSHits in this cell
549  if (hsize == 1) // there is only 1 hit in this cell, can make SimHit directly
550  {
551  hit.SetNP(1); // just 1 here
552  hit.SetEnD(earr);
553  hit.SetTim(carr);
554  hit.SetTrID(tarr);
555  hit.SetPID(parr);
556  }
557  else // size greater than 1 combine if duplicates or mother-daughters
558  {
559  // without sorting, can do with remove - then can adjust other arrays as well
560  std::vector<int>::iterator startloc(tarr.begin());
561  std::vector<int>::iterator endloc(tarr.end());
562  int value(*startloc);
563 
564  while (++startloc < endloc)
565  {
566  endloc = std::remove(startloc, endloc, value);
567  value = *startloc;
568  }
569  tarr.erase(endloc,tarr.end());
570  if (tarr.size()==1)
571  {
572  hit.SetNP(1);
573  hit.SetTrID(tarr);
574  std::vector<int> mparr (1,parr[0]);
575  hit.SetPID(mparr);
576  std::vector<double> mcarr (1,carr[0]);
577  hit.SetTim(mcarr);
578  double esum = 0.;
579  for (unsigned int i=0; i<hsize; ++i)
580  {
581  esum += earr[i];
582  }
583  std::vector<double> mearr (1,esum);
584  hit.SetEnD(mearr);
585  }
586  else // tarr size greater than 1
587  {
588  // if the size of the tarr array is the same as the size of tpp, for now, all hits are
589  // valid and unique and the arrays should be used to make a SimHit
590  if (hsize==tarr.size())
591  {
592  hit.SetNP(hsize);
593  hit.SetPID(parr);
594  hit.SetEnD(earr);
595  hit.SetTrID(tarr);
596  hit.SetTim(carr);
597  }
598  else
599  {
600  // tarr size less than hsize so must analyze
601  // can set the size to tarr size
602  hit.SetNP(tarr.size());
603  hit.SetTrID(tarr);
604 
605  // use pairs to re-generate PID, Energy, time arrays
606  std::vector<int> pvec(tarr.size());
607  std::vector<double> evec(tarr.size());
608  std::vector<double> tvec(tarr.size());
609  for (unsigned int i=0; i<tarr.size(); ++i)
610  {
611  for (unsigned int j=0; j<tparr.size(); ++j)
612  {
613  double lowt = 9999999999.;
614  if (tparr[j].first == tarr[i])
615  {
616  pvec[i] = tparr[j].second;
617  evec[i] += tearr[j].second;
618  if (tcarr[j].second > 0. && tcarr[j].second < lowt)
619  {
620  tvec[i] = tcarr[j].second;
621  lowt = tcarr[j].second;
622  }
623  }
624  }
625  } // end of tarr loop
626  hit.SetPID(pvec);
627  hit.SetEnD(evec);
628  hit.SetTim(tvec);
629  }
630  }
631  }
632 
633  // can put threshold and timing cuts here if wanted before making hit list
634  // 0.0005 is half MeV threshold cut
635  if( hit.fTotE > 0.00000005) // for now, just use a threshold cut
636  {
637  // finally, put this hit into list
638  simhits.push_back(hit);
639  Etotev += hit.fTotE;
640  NumSH++;
641  }
642  }
643 
644  fsimhitE->Fill(Etotev); //fill histogram of total SimHit E in event
645  fnumsh->Fill(NumSH); // fill histogram of number of SimHits in Event
646  std::cout << " Total E (GeV) in SimHits = " << Etotev << "\n";
647 
648  // check some SimHit stuff
649  std::cout << "Size of SimHit List = " << simhits.size() << "\n";
650  // loop over all SimHits to get Info
651  for (unsigned int i=0; i<simhits.size(); ++i)
652  {
653  cheat::SimHit h=simhits[i];
654  // std::cout <<" Num of particles for "<< i << " th hit = " << h.fNP << "\n"
655  // <<" Plane and Cell of this hit = " << h.fPlane << " " << h.fCell << "\n"
656  // <<" View of this hit = " << h.fView << "\n"
657  // <<" Total E of this hit = "<< h.fTotE << "\n";
658  // fill some histograms
659  fthreshE->Fill(h.fTotE * 1000.); // this is energy per hit in MeV
660  if (h.View()==geo::kX) ftotEXV->Fill(h.fTotE * 1000.); // this is energy per hit in MeV per X view
661  if (h.View()==geo::kY) ftotEYV->Fill(h.fTotE * 1000.); // this is energy per hti in MeV per Y view
662  // for (unsigned int j=0; j<h.fNP; ++j)
663  // {
664  // std::cout <<" Particle Edep, ID, TrID, Time = "<< h.fEnD[j] << " " << h.fPID[j] << " " << h.fTrID[j] << " " << h.fTim[j] << "\n";
665  // }
666  }
667 
668  // MAKE SIMTRACKS HERE
669  // to make SimTracks, use TrackID unique PlaneCell pair map - if a TrackID has more than
670  // 2 unique planecell pairs, both views represented, its a SimTrack
671  // will need a map linking trackID to SimHits for SimTracks
672  // std::map<unsigned int, art::Ptr<cheat::SimHit> > simtrhitmap;
673  std::map<unsigned int, std::vector<cheat::SimHit*> > simtrhitmap;
674  std::vector<cheat::SimTrack> simtracks; // initialize SimTracks
675  for (std::map<unsigned int, std::vector<std::pair<unsigned int, unsigned int> > >::iterator itm = tidupcpmap.begin(); itm != tidupcpmap.end(); ++itm)
676  {
677  unsigned int tid = itm->first;
678  unsigned int upcsize = itm->second.size();
679  // if only 1 or 2 unique pc pairs, this is not reconstructable as a track
680  // also eventually probably want to add view requirement for both views
681  if (upcsize>2)
682  {
683  // this is a SimTrack
685  // first, get start and end coords from true particle length using trackID particle map
686  int trpid = 0;
687  // make some variables to hold the trajectory information
688  TVector3 scoord; // Track start point
689  bool foundStart = false;
690  unsigned int startpoint = 0;
691  TVector3 ecoord; // Track end point
692  bool foundEnd = false;
693  unsigned int endpoint = 0;
694  std::vector<TLorentzVector> trpoints; // vector holding the particle trajectory points between the start and end pont.
695 
696 
697  std::map<int, const sim::Particle*>::const_iterator it;
698  for (it = trpartmap.begin(); it != trpartmap.end(); ++it)
699  {
700  unsigned int trid = it->first;
701  const sim::Particle* p=it->second;
702  if (trid==tid)
703  {
704 
705  trpid = p->PdgCode();
706 
707  // Get the starting position of the track as the first trajectory point inside the detector volume
708  unsigned int trajPoint = 0;
709  while(!foundStart && trajPoint<p->NumberTrajectoryPoints()){
710  TLorentzVector start = p->Position(trajPoint);
711  if(TMath::Abs(start.X())<fgeo->DetHalfWidth() && TMath::Abs(start.Y())<fgeo->DetHalfHeight() && start.Z()>0.0 && start.Z()<fgeo->DetLength()){
712  scoord.SetXYZ(start.X(),start.Y(),start.Z());
713  foundStart = true;
714  startpoint = trajPoint;
715  // std::cout<<"Found the start point: "<<startpoint<<std::endl;
716  }
717  else{
718  ++trajPoint;
719  }
720  }
721  //initiate endpoint to same as startpoint
722  endpoint = startpoint;
723  trajPoint = p->NumberTrajectoryPoints()-1;
724  // Get the endpoint as the last trajectory point that is inside the detector volume
725  while(!foundEnd && trajPoint>=startpoint){
726  TLorentzVector end = p->Position(trajPoint);
727  if(TMath::Abs(end.X())<fgeo->DetHalfWidth() && TMath::Abs(end.Y())<fgeo->DetHalfHeight() && end.Z()>0.0 && end.Z()<fgeo->DetLength()){
728  ecoord.SetXYZ(end.X(),end.Y(),end.Z());
729  foundEnd = true;
730  endpoint = trajPoint;
731  // std::cout<<"Found the end point: "<<endpoint<<std::endl;
732  }
733  else{
734  --trajPoint;
735  }
736  }
737  // Save all the trajectory points after start point up to the end point.
738  if(foundStart && foundEnd && endpoint>startpoint){
739  for(unsigned int trpoint = startpoint+1; trpoint <= endpoint;++trpoint){
740  // std::cout<<"Start point: "<<startpoint<<" End point: "<<endpoint<<std::endl;
741  // std::cout<<"Putting trajectory point into track list: "<<trpoint<<std::endl;
742  trpoints.push_back(p->Position(trpoint));
743  }
744  }
745 
746  }
747  }
748  // now only make a simtrack object if the track start and stop where found
749  if(foundStart && foundEnd && endpoint>=startpoint){
751 
752  // calculate direction cosines from start and end of this MC track
753  strack.SetPID(trpid);
754  strack.SetTrackID(tid);
755  strack.SetStart(scoord);
756  for(unsigned int trpoint = 0; trpoint < trpoints.size();++trpoint){
757  TVector3 pos(trpoints[trpoint].X(),trpoints[trpoint].Y(),trpoints[trpoint].Z());
758  strack.AppendTrajectoryPoint(pos);
759  // Get initial track direction
760  if(trpoint == 0){
761  // calculate direction cosines from the first and second trajectory points of this MC track
762 
763  double dist = sqrt((ecoord.X()-trpoints[trpoint].X())*(ecoord.X()-trpoints[trpoint].X())+
764  (ecoord.Y()-trpoints[trpoint].Y())*(ecoord.Y()-trpoints[trpoint].Y())+
765  (ecoord.Z()-trpoints[trpoint].Z())*(ecoord.Z()-trpoints[trpoint].Z()));
766  TVector3 dir((ecoord.X()-trpoints[trpoint].X())/dist,(ecoord.Y()-trpoints[trpoint].Y())/dist,(ecoord.Z()-trpoints[trpoint].Z())/dist);
767  strack.SetDir(dir);
768 
769  }
770  }
771  // add to SimTracks
772  simtracks.push_back(strack);
773  }
774 
775  // loop over SimHits, adding to the TrackID, SimHit map
776  // simtrhitmap[tid].resize(0); // Jon suggested this to fix map?
777  for (unsigned int i=0; i<simhits.size(); ++i)
778  {
779  bool ismatch = false;
780  // const cheat::SimHit h=simhits[i];
781  cheat::SimHit h=simhits[i];
782  for (unsigned int j=0; j<h.fNP; ++j)
783  {
784  if (h.fTrID[j]==(int) tid) ismatch = true;
785  }
786  if (ismatch) simtrhitmap[tid].push_back(& h);
787  }
788  } // end of upc greater than 2 loop
789  } // end of loop over trackID planecell pair map
790 
791  // write SimTracks to event
792  std::unique_ptr<std::vector<cheat::SimTrack> > simtrackcol(new std::vector<cheat::SimTrack>(simtracks) );
793  evt.put(std::move(simtrackcol));
794 
795  /*
796  // ********************** Diagnostics for MCCheater ***********************
797  // setup a test for single particles - if trajectory along z axis, test variable is true
798  bool zdir = false;
799  // report some SimTrack stuff
800  std::cout << "Size of SimTrack List = " << simtracks.size() << "\n";
801  double tlensum = 0.;
802  // define min and max coords for NDOS
803  double xmin = -130.1;
804  double xmax = 130.1;
805  double ymin = -195.1;
806  double ymax = 195.1;
807  double zmin = 0.;
808  double zmax = 1200.;
809  for (unsigned int i=0; i<simtracks.size(); ++i)
810  {
811  fzstart->Fill(simtracks[i].StartPos()[2]/6.0);
812  double xlen = simtracks[i].EndPos()[0]-simtracks[i].StartPos()[0];
813  double ylen = simtracks[i].EndPos()[1]-simtracks[i].StartPos()[1];
814  double zlen = simtracks[i].EndPos()[2]-simtracks[i].StartPos()[2];
815  fnplanes->Fill(zlen/6.);
816  double tlength = sqrt(xlen*xlen+ylen*ylen+zlen*zlen);
817  tlensum += tlength;
818  fdcosX->Fill(simtracks[i].Dir()[0]);
819  fdcosY->Fill(simtracks[i].Dir()[1]);
820  fdcosZ->Fill(simtracks[i].Dir()[2]);
821  if (simtracks[i].Dir()[2]>0.5) zdir = true;
822  if (zdir && simtracks[i].TrackID()==0) ftlen->Fill(tlength);
823  if (simtracks[i].Dir()[2]>0.5 && simtracks[i].PID()==13) fmuL->Fill(tlength);
824  if (simtracks[i].Dir()[2]>0.5 && simtracks[i].PID()==11) felL->Fill(tlength);
825  std::cout << " Track Number " << i << "\n";
826  std::cout << " Track ID = " << simtracks[i].TrackID() << " Pdg ID = " << simtracks[i].PID() << "\n";
827  std::cout << " Start Coord = " << simtracks[i].StartPos()[0] << " " << simtracks[i].StartPos()[1] << " " << simtracks[i].StartPos()[2] << "\n";
828  std::cout << " End Coord = " << simtracks[i].EndPos()[0] << " " << simtracks[i].EndPos()[1] << " " << simtracks[i].EndPos()[2] << "\n";
829  std::cout << " Direction = " << simtracks[i].Dir()[0] << " " << simtracks[i].Dir()[1] << " " << simtracks[i].Dir()[2] << "\n";
830  // find point of intersection with plane of NDOS
831  // first, use top plane
832  double ttop = (195.-simtracks[i].StartPos()[1])/ylen;
833  double xint = simtracks[i].StartPos()[0]+ttop*(xlen);
834  double yint = simtracks[i].StartPos()[1]+ttop*(ylen);
835  double zint = simtracks[i].StartPos()[2]+ttop*(zlen);
836  std::cout << " TP ttrack = " << ttop << " Coords = " << xint << " " << yint << " " << zint << "\n";
837  if (xint>xmin && xint<xmax && yint>ymin && yint<ymax && zint>zmin && zint<zmax)
838  {
839  // loop over all simhits to find closest to this point
840  double sdist = 999.;
841  for (unsigned int j=0; j<simhits.size(); ++j)
842  {
843  cheat::SimHit h=simhits[j];
844  // does this hit contain contribution from this track?
845  bool htid = false;
846  for (unsigned int k=0; k<h.fNP; ++k)
847  {
848  if (h.fTrID[k]==simtracks[i].TrackID()) htid = true;
849  }
850  if (htid)
851  {
852  double hdist = sqrt((h.fCoorX-xint)*(h.fCoorX-xint)+(h.fCoorY-yint)*(h.fCoorY-yint)+(h.fCoorZ-zint)*(h.fCoorZ-zint));
853  // std::cout << " Hit distance = " << hdist << "\n";
854  if (hdist<sdist) sdist = hdist;
855  } // satisfies trackID
856  } // end of SimHit loop
857  std::cout << " Closest distance to edge = " << sdist << " for TrackID = " << simtracks[i].TrackID() << "\n";
858  fclhit->Fill(sdist);
859  }
860  }
861  */
862 
863  /*
864  double E0SHTot = 0.;
865  if (zdir)
866  {
867  for (unsigned int i=0; i<simhits.size(); ++i)
868  {
869  // cheat::SimHit hit=simhits[i];
870  E0SHTot += simhits[i].fTotE;
871  }
872  fEpar->Fill(E0SHTot); // histogram of forward-going energy sum
873  }
874  */
875 
876  // check map of SimTrack TrackID SimHit
877  // this map doesn't work?
878  /*
879  for (std::map<unsigned int, std::vector<cheat::SimHit*> >::iterator itm = simtrhitmap.begin(); itm != simtrhitmap.end(); ++itm)
880  {
881  std::vector<cheat::SimHit*> sharr = itm->second;
882  mf::LogInfo("MCCheater") <<" Number of SimHits = " << sharr.size() << " For SimTrackID = " << itm->first << "\n";
883  fnshptr->Fill(sharr.size());
884  double stc[3];
885  double enc[3];
886  stc[1] = -210.;
887  enc[1] = 210.;
888 
889  for (unsigned int i=0; i<sharr.size(); ++i)
890  {
891  cheat::SimHit* h=sharr[i];
892  double xcor = h->fCoorX;
893  double ycor = h->fCoorY;
894  double zcor = h->fCoorZ;
895  // histogram these coordinates
896  fXRange->Fill(xcor);
897  fYRange->Fill(ycor);
898  fZRange->Fill(zcor);
899 
900  if (ycor > stc[1])
901  {
902  stc[1] = ycor;
903  stc[2] = zcor;
904  stc[0] = xcor;
905  continue;
906  } else if (ycor < stc[1])
907  {
908  if (ycor < enc[1])
909  {
910  enc[1] = ycor;
911  enc[2] = zcor;
912  enc[0] = xcor;
913  continue;
914  }
915  }
916  }
917  fstsxy->Fill(stc[0],stc[1]);
918  fstexy->Fill(enc[0],enc[1]);
919  }
920  */
921 
922  // write SimHits to event
923  std::unique_ptr<std::vector<cheat::SimHit> > simhitcol(new std::vector<cheat::SimHit>(simhits) );
924  std::cout << " Wrote SimHits to event of size = " << simhits.size() << "\n";
925  evt.put(std::move(simhitcol));
926 
927  return;
928  }
double E(const int i=0) const
Definition: MCParticle.h:232
unsigned int NumberTrajectoryPoints() const
Definition: MCParticle.h:217
const TLorentzVector & Position(const int i=0) const
Definition: MCParticle.h:218
int PdgCode() const
Definition: MCParticle.h:211
virtual void SetStart(TVector3 start)
Definition: Track.cxx:168
void GetCenter(double *xyz, double localz=0.0) const
Definition: CellGeo.cxx:159
set< int >::iterator it
const sim::ParticleNavigator & ParticleNavigator() const
Get a reference to the ParticleNavigator.
Definition: BackTracker.h:744
enum geo::_plane_proj View_t
Enumerate the possible plane projections.
MaybeLogger_< ELseverityLevel::ELsev_info, false > LogInfo
void LocalToWorld(const double *local, double *world) const
Definition: CellGeo.cxx:80
int Mother() const
Definition: MCParticle.h:212
void SetPID(int pid)
Definition: SimTrack.h:28
geo::View_t View() const
Definition: CellHit.h:41
const CellGeo * Cell(int icell) const
Definition: PlaneGeo.h:48
const char * p
Definition: xmltok.h:285
T sqrt(T number)
Definition: d0nt_math.hpp:156
list_type::const_iterator const_iterator
Vertical planes which measure X.
Definition: PlaneGeo.h:28
double DetLength() const
void SetChannel(uint32_t iChan)
Definition: RawDigit.h:99
void abs(TH1 *hist)
const PlaneGeo * Plane(unsigned int i) const
MaybeLogger_< ELseverityLevel::ELsev_error, false > LogError
void SetTotE(double tote)
Definition: SimHit.h:50
void SetTrackID(int trackid)
Definition: SimTrack.h:29
Horizontal planes which measure Y.
Definition: PlaneGeo.h:29
int TrackId() const
Definition: MCParticle.h:209
Float_t Y
Definition: plot.C:38
void SetView(geo::View_t view)
Definition: CellHit.h:54
void SetPlane(unsigned short plane)
Definition: CellHit.h:53
double dist
Definition: runWimpSim.h:113
Float_t Z
Definition: plot.C:38
ProductID put(std::unique_ptr< PROD > &&product)
Definition: Event.h:102
void SetCell(unsigned short cell)
Definition: CellHit.h:52
TString part[npart]
Definition: Style.C:32
View_t View() const
Which coordinate does this plane measure.
Definition: PlaneGeo.h:53
void hits()
Definition: readHits.C:15
void SetPE(float pe)
Definition: CellHit.h:55
const XML_Char int const XML_Char * value
Definition: expat.h:331
void SetNP(unsigned short np)
Definition: SimHit.h:51
void SetCoorZ(double coorz)
Definition: SimHit.h:49
virtual void SetDir(TVector3 dir)
Definition: Prong.cxx:104
Float_t strack
Definition: plot.C:35
void AppendTrajectoryPoint(TVector3 pt)
Definition: Track.cxx:112
void SetPID(std::vector< int > PID)
Definition: SimHit.h:53
void SetTim(std::vector< double > Tim)
Definition: SimHit.h:55
const double j
Definition: BetheBloch.cxx:29
double DetHalfHeight() const
void SetTrID(std::vector< int > TrID)
Definition: SimHit.h:52
OStream cout
Definition: OStream.cxx:6
void SetCoorX(double coorx)
Definition: SimHit.h:47
void SetEnD(std::vector< double > EnD)
Definition: SimHit.h:54
double DetHalfWidth() const
TDirectory * dir
Definition: macro.C:5
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
std::string fFLSHitListLabel
Definition: structs.h:12
double fTotE
Total Energy dep from all particles.
Definition: SimHit.h:65
Float_t X
Definition: plot.C:38
Definition: fwd.h:28
void SetCoorY(double coory)
Definition: SimHit.h:48
unsigned short fNP
size of vectors
Definition: SimHit.h:66
std::vector< int > fTrID
vector of sim track IDs
Definition: SimHit.h:67
int EveId(const int trackID) const
void art::Consumer::showMissingConsumes ( ) const
protectedinherited

Referenced by art::RootOutput::endJob().

void art::Consumer::validateConsumedProduct ( BranchType const  bt,
ProductInfo const &  pi 
)
protectedinherited

Member Data Documentation

TH1F* cheat::MCCheater::fclhit
private

Definition at line 117 of file MCCheater_module.cc.

Referenced by beginJob().

TH1F* cheat::MCCheater::fcorx
private

Definition at line 113 of file MCCheater_module.cc.

Referenced by beginJob(), and produce().

TH1F* cheat::MCCheater::fcory
private

Definition at line 114 of file MCCheater_module.cc.

Referenced by beginJob(), and produce().

TH1F* cheat::MCCheater::fdcosX
private

Definition at line 90 of file MCCheater_module.cc.

Referenced by beginJob().

TH1F* cheat::MCCheater::fdcosY
private

Definition at line 91 of file MCCheater_module.cc.

Referenced by beginJob().

TH1F* cheat::MCCheater::fdcosZ
private

Definition at line 92 of file MCCheater_module.cc.

Referenced by beginJob().

TH1F* cheat::MCCheater::fElE
private

Definition at line 84 of file MCCheater_module.cc.

Referenced by beginJob(), and produce().

TH1F* cheat::MCCheater::felL
private

Definition at line 85 of file MCCheater_module.cc.

Referenced by beginJob().

TH1F* cheat::MCCheater::fFLSE
private

Definition at line 96 of file MCCheater_module.cc.

Referenced by beginJob(), and produce().

std::string cheat::MCCheater::fFLSHitListLabel
private

Definition at line 78 of file MCCheater_module.cc.

Referenced by produce().

std::string cheat::MCCheater::fMCTruthListLabel
private

Definition at line 77 of file MCCheater_module.cc.

TH1F* cheat::MCCheater::fMuE
private

Definition at line 83 of file MCCheater_module.cc.

Referenced by beginJob(), and produce().

TH1F* cheat::MCCheater::fmuL
private

Definition at line 86 of file MCCheater_module.cc.

Referenced by beginJob().

TH1F* cheat::MCCheater::fnpart
private

Definition at line 107 of file MCCheater_module.cc.

Referenced by beginJob(), and produce().

TH1F* cheat::MCCheater::fnplanes
private

Definition at line 102 of file MCCheater_module.cc.

Referenced by beginJob().

TH1F* cheat::MCCheater::fntrajpts
private

Definition at line 108 of file MCCheater_module.cc.

Referenced by beginJob(), and produce().

TH1F* cheat::MCCheater::fNumFLS
private

Definition at line 98 of file MCCheater_module.cc.

Referenced by beginJob(), and produce().

TH1F* cheat::MCCheater::fnumsh
private

Definition at line 97 of file MCCheater_module.cc.

Referenced by beginJob(), and produce().

TH1F* cheat::MCCheater::fparte
private

Definition at line 106 of file MCCheater_module.cc.

Referenced by beginJob(), and produce().

TH1F* cheat::MCCheater::fsimhitE
private

Definition at line 94 of file MCCheater_module.cc.

Referenced by beginJob(), and produce().

TH1F* cheat::MCCheater::fthreshE
private

Definition at line 101 of file MCCheater_module.cc.

Referenced by beginJob(), and produce().

TH1F* cheat::MCCheater::ftlen
private

Definition at line 99 of file MCCheater_module.cc.

Referenced by beginJob().

TH1F* cheat::MCCheater::ftotalE
private

Definition at line 87 of file MCCheater_module.cc.

Referenced by beginJob(), and produce().

TH1F* cheat::MCCheater::ftotEXV
private

Definition at line 88 of file MCCheater_module.cc.

Referenced by beginJob(), and produce().

TH1F* cheat::MCCheater::ftotEYV
private

Definition at line 89 of file MCCheater_module.cc.

Referenced by beginJob(), and produce().

TH1F* cheat::MCCheater::fzstart
private

Definition at line 118 of file MCCheater_module.cc.

Referenced by beginJob().

int cheat::MCCheater::ntest
private

Definition at line 75 of file MCCheater_module.cc.

Referenced by produce().


The documentation for this class was generated from the following file: