Public Member Functions | Private Member Functions | Private Attributes | List of all members
CalibrationTreeBuilder::CalibrationTreeBuilder Class Reference

#include <CalibrationTreeBuilder.h>

Inheritance diagram for CalibrationTreeBuilder::CalibrationTreeBuilder:
art::EDAnalyzer art::detail::Analyzer art::detail::LegacyModule art::Observer art::ModuleBase

Public Member Functions

 CalibrationTreeBuilder (fhicl::ParameterSet const &pSet)
 
virtual void beginJob () override
 
virtual void analyze (const art::Event &evt) override
 
virtual void endJob () override
 
- Public Member Functions inherited from art::EDAnalyzer
 EDAnalyzer (fhicl::ParameterSet const &pset)
 
template<typename Config >
 EDAnalyzer (Table< Config > const &config)
 
std::string workerType () const
 
- Public Member Functions inherited from art::detail::Analyzer
virtual ~Analyzer () noexcept
 
 Analyzer (fhicl::ParameterSet const &pset)
 
template<typename Config >
 Analyzer (Table< Config > const &config)
 
void doBeginJob (SharedResources const &resources)
 
void doEndJob ()
 
void doRespondToOpenInputFile (FileBlock const &fb)
 
void doRespondToCloseInputFile (FileBlock const &fb)
 
void doRespondToOpenOutputFiles (FileBlock const &fb)
 
void doRespondToCloseOutputFiles (FileBlock const &fb)
 
bool doBeginRun (RunPrincipal &rp, ModuleContext const &mc)
 
bool doEndRun (RunPrincipal &rp, ModuleContext const &mc)
 
bool doBeginSubRun (SubRunPrincipal &srp, ModuleContext const &mc)
 
bool doEndSubRun (SubRunPrincipal &srp, ModuleContext const &mc)
 
bool doEvent (EventPrincipal &ep, ModuleContext const &mc, std::atomic< std::size_t > &counts_run, std::atomic< std::size_t > &counts_passed, std::atomic< std::size_t > &counts_failed)
 
- Public Member Functions inherited from art::Observer
 ~Observer () noexcept
 
 Observer (Observer const &)=delete
 
 Observer (Observer &&)=delete
 
Observeroperator= (Observer const &)=delete
 
Observeroperator= (Observer &&)=delete
 
void registerProducts (ProductDescriptions &, ModuleDescription const &)
 
void fillDescriptions (ModuleDescription const &)
 
fhicl::ParameterSetID selectorConfig () const
 
- Public Member Functions inherited from art::ModuleBase
virtual ~ModuleBase () noexcept
 
 ModuleBase ()
 
ModuleDescription const & moduleDescription () const
 
void setModuleDescription (ModuleDescription const &)
 
std::array< std::vector< ProductInfo >, NumBranchTypes > const & getConsumables () const
 
void sortConsumables (std::string const &current_process_name)
 
template<typename T , BranchType BT>
ViewToken< T > consumesView (InputTag const &tag)
 
template<typename T , BranchType BT>
ViewToken< T > mayConsumeView (InputTag const &tag)
 

Private Member Functions

bool AddHit (detinfo::DetectorClocksData const &clockData, const art::Ptr< recob::Hit > hit, unsigned int &counter)
 
bool AddHit (const art::Ptr< recob::OpHit > hit, unsigned int &counter)
 
void PrepDivRec (const art::Event &evt)
 
std::pair< std::vector< CalibTreeRecord::EveRecord >::iterator, boolEmplaceEve (const simb::MCParticle *new_eve)
 
std::pair< std::vector< CalibTreeRecord::ParticleRecord >::iterator, boolEmplaceParticle (const simb::MCParticle *new_part)
 
const art::Ptr< sim::OpDetDivRecFindDivRec (int const &opDetNum) const
 
const std::vector< sim::SDPOpHitToChannelWeightedSimSDPs (art::Ptr< recob::OpHit > const &opHit_P) const
 

Private Attributes

art::ServiceHandle< cheat::ParticleInventoryServicePIS
 
art::ServiceHandle< cheat::BackTrackerServiceBTS
 
art::ServiceHandle< cheat::PhotonBackTrackerServicePBS
 
art::ServiceHandle< geo::GeometryGS
 
art::InputTag private_HitLabel
 
art::InputTag private_OpHitLabel
 
art::InputTag fWavLabel
 
art::ServiceHandle< art::TFileService > private_service_tfs
 
TTree * private_CalibrationTree
 
TTree * private_FlatCalibrationTree
 
TBranch * private_CalibrationRecord
 
TBranch * private_FlatCalibrationRecord
 
CalibTreeRecord::CalibTreeRecord private_eventBuffer
 
CalibTreeRecord::CalibTreeRecord private_eventPrep
 
flatfiller fl
 
std::vector< art::Ptr< sim::OpDetDivRec > > priv_DivRecs
 

Additional Inherited Members

- Public Types inherited from art::EDAnalyzer
using WorkerType = WorkerT< EDAnalyzer >
 
using ModuleType = EDAnalyzer
 
- Protected Member Functions inherited from art::Observer
std::string const & processName () const
 
bool wantAllEvents () const noexcept
 
bool wantEvent (ScheduleID id, Event const &e) const
 
Handle< TriggerResultsgetTriggerResults (Event const &e) const
 
 Observer (fhicl::ParameterSet const &config)
 
 Observer (std::vector< std::string > const &select_paths, std::vector< std::string > const &reject_paths, fhicl::ParameterSet const &config)
 
- Protected Member Functions inherited from art::ModuleBase
ConsumesCollectorconsumesCollector ()
 
template<typename T , BranchType = InEvent>
ProductToken< T > consumes (InputTag const &)
 
template<typename Element , BranchType = InEvent>
ViewToken< Element > consumesView (InputTag const &)
 
template<typename T , BranchType = InEvent>
void consumesMany ()
 
template<typename T , BranchType = InEvent>
ProductToken< T > mayConsume (InputTag const &)
 
template<typename Element , BranchType = InEvent>
ViewToken< Element > mayConsumeView (InputTag const &)
 
template<typename T , BranchType = InEvent>
void mayConsumeMany ()
 

Detailed Description

Definition at line 89 of file CalibrationTreeBuilder.h.

Constructor & Destructor Documentation

CalibrationTreeBuilder::CalibrationTreeBuilder::CalibrationTreeBuilder ( fhicl::ParameterSet const &  pSet)

Definition at line 57 of file CalibrationTreeBuilder_module.cc.

58  :EDAnalyzer(pSet),
59  private_HitLabel(pSet.get<art::InputTag>("HitLabel","gaushit")),
60  private_OpHitLabel(pSet.get<art::InputTag>("OpHitLabel","ophit")),
61  fWavLabel(pSet.get<art::InputTag>("WavLabel", "opdigi"))
62  { }
EDAnalyzer(fhicl::ParameterSet const &pset)
Definition: EDAnalyzer.h:25

Member Function Documentation

bool CalibrationTreeBuilder::CalibrationTreeBuilder::AddHit ( detinfo::DetectorClocksData const &  clockData,
const art::Ptr< recob::Hit hit,
unsigned int &  counter 
)
private

Definition at line 403 of file CalibrationTreeBuilder_module.cc.

404  {
406  std::vector<const sim::IDE*> ides_ps = BTS->HitToSimIDEs_Ps(clockData, hit);
407  std::vector<std::pair<int, CalibTreeRecord::PartialHit>> track_partials;
408  Double_t total_charge = 0.0;
409  if(ides_ps.empty()){return false;}
410  for(auto ide_p : ides_ps){
411  int track_id = ide_p->trackID;
413  total_charge += ide_p->numElectrons;
414  // tmp.charge = ide_p->numElectrons;
415  tmp.num_electrons = ide_p->numElectrons;
416  tmp.time = hit->PeakTime();
417  tmp.width = hit->SigmaPeakTime();
418  tmp.split = 0.0;
419  tmp.wire = hit->Channel();
420  tmp.is_collection_wire = (geom->SignalType(tmp.wire)==geo::kCollection)?(true):(false);
421  tmp.energy = (tmp.is_collection_wire)?(ide_p->energy):(0.0);
422  tmp.index = counter;
423  track_partials.push_back(std::make_pair(track_id, tmp));
424  }
425  for(auto track_partial : track_partials){
426  track_partial.second.split = track_partial.second.num_electrons / total_charge; //What percentage of the hit is in this partial.
427  track_partial.second.charge = hit->Integral()*track_partial.second.split; //Charge now normalized by deposition;
428  //Add partial to CalibTreeRecord;
429  const simb::MCParticle* part = PIS->TrackIdToParticle_P(track_partial.first);
430  //Emplace Particle. (Emplace Eve will be handeled in emplace particle
432  part_marker.first->partial_hits.push_back(track_partial.second); //Add post normalized partial hit.
433  //Don't forget to add the indicies to the hit vector. Actually. To keep it easy, we can just loop the eve records and add the indicies then. This is inefficient, but easier to code since I don't immediately have access to eve_pos here.
434  mf::LogDebug(__FUNCTION__)<<"This partialhit is being added to the record with charge "<<track_partial.second.charge<<"\n";
435  }//end for track partial.
436 
437  return true;
438  }
intermediate_table::iterator iterator
const simb::MCParticle * TrackIdToParticle_P(int id) const
std::pair< std::vector< CalibTreeRecord::ParticleRecord >::iterator, bool > EmplaceParticle(const simb::MCParticle *new_part)
float Integral() const
Integral under the calibrated signal waveform of the hit, in tick x ADC units.
Definition: Hit.h:224
SigType_t SignalType(geo::PlaneID const &pid) const
Returns the type of signal on the channels of specified TPC plane.
string tmp
Definition: languages.py:63
art::ServiceHandle< cheat::BackTrackerService > BTS
art::ServiceHandle< cheat::ParticleInventoryService > PIS
float PeakTime() const
Time of the signal peak, in tick units.
Definition: Hit.h:218
MaybeLogger_< ELseverityLevel::ELsev_success, false > LogDebug
float SigmaPeakTime() const
Uncertainty for the signal peak, in tick units.
Definition: Hit.h:219
std::vector< const sim::IDE * > HitToSimIDEs_Ps(detinfo::DetectorClocksData const &clockData, recob::Hit const &hit) const
raw::ChannelID_t Channel() const
ID of the readout channel the hit was extracted from.
Definition: Hit.h:230
Signal from collection planes.
Definition: geo_types.h:146
bool CalibrationTreeBuilder::CalibrationTreeBuilder::AddHit ( const art::Ptr< recob::OpHit hit,
unsigned int &  counter 
)
private

Definition at line 440 of file CalibrationTreeBuilder_module.cc.

440  {
441  {//make each loop local so the vectors don't persist
442 // try{
443  const std::vector< sim::SDP> sdps = this->OpHitToChannelWeightedSimSDPs(hit);
444 // }catch(const std::exception &e){
445 // return false;
446 // }
447  std::vector<std::pair<int, CalibTreeRecord::PartialOpHit>> track_partials;
448  Double_t total_charge = 0.0;
449  if(sdps.empty()){return false;}
450  for(auto sdp : sdps){
451  int track_id = sdp.trackID;
453  total_charge += sdp.numPhotons;
454  // tmp.pes = //Not used yet.
455  tmp.num_photons = sdp.numPhotons; //This is not PEs. This is incident photons. I am just storing it here for later logic (see final assignment in the next for loop.
456  tmp.energy = sdp.energy;
457  //tmp.energy = 0.0;// sdp.energy;
458  tmp.time = hit->PeakTime();
459  tmp.width = hit->Width();
460  tmp.split = 0.0;
461  tmp.opchan = hit->OpChannel();
462  tmp.opdet = GS->OpDetFromOpChannel(hit->OpChannel());
463  tmp.index = counter;
464  track_partials.push_back(std::make_pair(track_id, tmp));
465  }
466  for(auto track_partial : track_partials){
467  track_partial.second.split = track_partial.second.num_photons / total_charge; //What percentage of the hit is in this partial.
468  track_partial.second.pes = hit->PE()*track_partial.second.split; //PEs now normalized by deposition;
469  //Add partial to CalibTreeRecord;
470  const simb::MCParticle* part = PIS->TrackIdToParticle_P(track_partial.first);
471  //Emplace Particle. (Emplace Eve will be handeled in emplace particle
473  part_marker.first->partial_ophits.push_back(track_partial.second); //Add post normalized partial hit.
474  //Don't forget to add the indicies to the hit vector. Actually. To keep it easy, we can just loop the eve records and add the indicies then. This is inefficient, but easier to code since I don't immediately have access to eve_pos here.
475  }//end for track partial.
476 
477  }
478  return true;
479 // }
480 // catch(...){
481 // return false;
482 // }
483  }
intermediate_table::iterator iterator
const simb::MCParticle * TrackIdToParticle_P(int id) const
double PeakTime() const
Definition: OpHit.h:64
std::pair< std::vector< CalibTreeRecord::ParticleRecord >::iterator, bool > EmplaceParticle(const simb::MCParticle *new_part)
unsigned int OpDetFromOpChannel(int opChannel) const
Convert unique channel to detector number.
const std::vector< sim::SDP > OpHitToChannelWeightedSimSDPs(art::Ptr< recob::OpHit > const &opHit_P) const
double Width() const
Definition: OpHit.h:66
string tmp
Definition: languages.py:63
art::ServiceHandle< cheat::ParticleInventoryService > PIS
double PE() const
Definition: OpHit.h:69
int OpChannel() const
Definition: OpHit.h:62
void CalibrationTreeBuilder::CalibrationTreeBuilder::analyze ( const art::Event evt)
overridevirtual

Definition at line 83 of file CalibrationTreeBuilder_module.cc.

83  {
84  //eprV.clear();
87  priv_DivRecs.clear();
88  this->PrepDivRec(evt);
89 
90 
91  //Set up BackTracking //Now in header
92  /* art::ServiceHandle<cheat::ParticleInventoryService> PIS;
93  art::ServiceHandle<cheat::BackTrackerService> BTS;
94  art::ServiceHandle<cheat::PhotonBackTrackerService> PBS;*/
95  //Get a list of generator names.
96  auto const mcHandles = evt.getMany<std::vector<simb::MCTruth>>();
97  std::vector< std::pair<int, std::string>> track_id_to_label;
98 
99  for( auto const& mcHandle : mcHandles ){
100  const std::string& sModuleLabel = mcHandle.provenance()->moduleLabel();
101  art::FindManyP<simb::MCParticle> findMCParts(mcHandle, evt, "largeant");
102  std::vector<art::Ptr<simb::MCParticle> > mcParts = findMCParts.at(0);
103  for( const art::Ptr<simb::MCParticle> ptr : mcParts){
104  int track_id = ptr->TrackId();
105  gf->add(track_id, sModuleLabel);
106  }
107  }
108 
109  std::vector<art::Ptr<recob::Hit>> hitList;
110  auto hitHandle = evt.getHandle<std::vector<recob::Hit>>(private_HitLabel);
111  if ( hitHandle )
112  art::fill_ptr_vector(hitList, hitHandle);
113  std::vector<art::Ptr<recob::OpHit>> opHitList;
114  auto opHitHandle = evt.getHandle<std::vector<recob::OpHit>>(private_OpHitLabel);
115  if ( opHitHandle )
116  art::fill_ptr_vector(opHitList, opHitHandle);
117 
118  //private_eventBuffer is the event record
120  private_eventBuffer.run=evt.id().run();
123  //private_eventPrep is the event record
125  private_eventPrep.run=evt.id().run();
128 
129  auto const clockData = art::ServiceHandle<detinfo::DetectorClocksService const>()->DataFor(evt);
130  unsigned int hitCounter=0;
131  for ( auto const& hit : hitList){
132  hitCounter++;
133  if(!AddHit(clockData, hit, hitCounter)){
134  mf::LogDebug(__FUNCTION__)<<"Hit "<<hitCounter<<" failed to add to the record. Is this hit caused by noise?";
135  }
136  }//end for hit in hitlist
137 
138  hitCounter=0;
139  for ( auto const& hit : opHitList ){
140  ++hitCounter;
141  if(!AddHit(hit, hitCounter)){
142  mf::LogDebug(__FUNCTION__)<<"OpHit "<<hitCounter<<" failed to add to the record. Is this hit caused by noise?";
143  }
144  } //Add ophits to hitmap.
145 
146  //Buffer normalized per hit as it is built.
147  //Make the event record
148  /*
149  if(!private_eventPrep.eves.empty()){
150  mf::LogDebug(__FUNCTION__)<<"EventBuffer for CalibrationTree contains "<<private_eventPrep.eves.size()<<" eves.\n";
151  if(!private_eventPrep.eves.at(0).particles.empty()){
152  mf::LogDebug(__FUNCTION__)<<"The first particle has "<<private_eventPrep.eves.at(0).particles.size()<<" particles.\n";
153  if(!private_eventPrep.eves.at(0).particles.at(0).partial_hits.empty()){
154  mf::LogDebug(__FUNCTION__)<<"The first of these particles has "<<private_eventPrep.eves.at(0).particles.at(0).partial_hits.size()<<" partial hits.\n";
155  }else{
156  mf::LogDebug(__FUNCTION__)<<"The particle contains no hits.\n";
157  }
158  }else{
159  mf::LogDebug(__FUNCTION__)<<"Eve contains no particles.\n";
160  }
161  }else{
162  mf::LogDebug(__FUNCTION__)<<"EventBuffer for CalibrationTree contains no eves.\n";
163  }
164  */
165 
166  {//Forced local scope for a few variables. This is a good indication that the code I am writing doesn't belong here, but in another function.
167  int eve_pos=0;
168  for(auto eve : private_eventPrep.eves){
169  if(!private_eventPrep.eves.empty()){ //This is not needed and is redundant.
170  int part_pos=0;
171  for(auto part : eve.particles){
172  { //keep partial_hit_pos local scope for the loop
173  int partial_hit_pos=0;
174  for(auto partial_hit : part.partial_hits){
175  int ind = partial_hit.index;
176  std::vector<CalibTreeRecord::HitContributor>::iterator hc_pos = std::find_if(private_eventPrep.hits.begin(), private_eventPrep.hits.end(), [&ind](CalibTreeRecord::HitContributor& hc){ return ind==hc.index;});
177  if(hc_pos==private_eventPrep.hits.end()){
178  private_eventPrep.hits.emplace_back(ind);
179  if(hc_pos!=(private_eventPrep.hits.end()-1)){
180  hc_pos=private_eventPrep.hits.end()-1;
181  }
182  }
183  hc_pos->locations.emplace_back(eve_pos,part_pos,partial_hit_pos);
184  // private_eventPrep.hits.emplace_back(eve_pos,part_pos,partial_hit_pos);
185  ++partial_hit_pos;
186  }
187  }
188  {
189  int partial_ophit_pos=0;
190  for(auto partial_ophit : part.partial_ophits){
191  int ind = partial_ophit.index;
192  std::vector<CalibTreeRecord::HitContributor>::iterator hc_pos = std::find_if(private_eventPrep.ophits.begin(), private_eventPrep.ophits.end(), [&ind](CalibTreeRecord::HitContributor& hc){ return ind==hc.index;});
193  if(hc_pos==private_eventPrep.ophits.end()){
194  private_eventPrep.ophits.emplace_back(ind);
195  if(hc_pos!=(private_eventPrep.ophits.end()-1)){
196  hc_pos=private_eventPrep.ophits.end()-1;
197  }
198  }
199  hc_pos->locations.emplace_back(eve_pos,part_pos,partial_ophit_pos);
200  ++partial_ophit_pos;
201  }
202  }
203  ++part_pos;
204  }
205  }
206  ++eve_pos;
207  }
208  }
209 
210  //This is where the partial hit bunching should happen.
211  //I can use the hits records to bunch partials on the same particle. This however would mean I would then need to recalculate the positions for all hits post bunching.
212  //Now we loop the private_eventPrep and build private_eventBuffer.
213 
214  if( private_eventPrep.bunch_hits==true){
215  for(size_t eve_num=0; eve_num<private_eventPrep.eves.size(); ++eve_num ){
216  auto eve = private_eventPrep.eves.at(eve_num);
217  private_eventBuffer.eves.push_back(eve); //We are gettign them in order, we are filling them in order. This should not break anything.
218  //private_eventBuffer.eves.at(eve_num).particles.clear();//clear the particles list from the copeid eve. We will now refill it.
219  for(size_t part_num=0; part_num<eve.particles.size(); ++part_num ){
220  auto part = eve.particles.at(part_num);
221  auto& npart = private_eventBuffer.eves.at(eve_num).particles.at(part_num);
222  npart.partial_hits.clear();
223  npart.partial_ophits.clear(); //I did not clear the particles from the new list because I don't need to. They will be the same particles. The only thing that chages are the partial records, as we are bunching them.
224  //add part to particle buffer.
225  std::map<int64_t, CalibTreeRecord::PartialHit> bunched_hits;
226  std::map<int64_t, std::list<Double_t>> remembered_energies;
227  for(auto partialhit : part.partial_hits){ //No value in using an index loop here because I don't need the hit number (it will be internal, and won't match the counter);
228  auto emp_h = bunched_hits.emplace(std::make_pair(partialhit.index,partialhit));
229  std::list<Double_t> dummy;
230  auto emp_e = remembered_energies.emplace(std::make_pair(partialhit.index, dummy));//very memory inefficient, but I am not in a position to care about that right now. One could easily improve on my code here.
231  if(emp_h.second==true){ //emp_e.second and emp_h.second are always the same.
232  emp_e.first->second.emplace_back(partialhit.energy);
233  emp_h.first->second=partialhit;
234  }else{
235  if(emp_h.first->second.time != partialhit.time){throw cet::exception("CalibrationTreeBuilder")<<"Buncher failed. Trying to merge two of the same hit at different times?\n";}
236  emp_e.first->second.emplace_back(partialhit.energy);
237  emp_h.first->second.charge+=partialhit.charge;
238  emp_h.first->second.split+=partialhit.split;
239  emp_h.first->second.num_electrons+=partialhit.num_electrons;//numelectrons came from the ides that contributed to the given hit. They should stack because the IDEs stack(which is why we are adding them in the first place).
240  }
241  emp_e.first->second.unique();//Throw out duplicates.
242  Double_t energy=0.0;
243  for( auto en : emp_e.first->second){
244  energy += en;
245  }
246  emp_h.first->second.energy=energy;
247  }
248  remembered_energies.clear();
249  std::map<int64_t, CalibTreeRecord::PartialOpHit> bunched_ophits;
250  for(auto partialop : part.partial_ophits){
251  auto emp_o = bunched_ophits.emplace(std::make_pair(partialop.index,partialop));
252  std::list<Double_t> dummy;
253  auto emp_e = remembered_energies.emplace(std::make_pair(partialop.index, dummy));//very memory inefficient, but I am not in a position to care about that right now. One could easily improve on my code here.
254  if(emp_o.second==true){
255  emp_e.first->second.emplace_back(partialop.energy);
256  emp_o.first->second=partialop;
257  }else{
258  if(emp_o.first->second.time != partialop.time){throw cet::exception("CalibrationTreeBuilder")<<"Buncher failed. Trying to merge two of the same ophit at different times?\n";}
259  emp_e.first->second.emplace_back(partialop.energy);
260  emp_o.first->second.pes+=partialop.pes;
261  emp_o.first->second.split+=partialop.split;
262  emp_o.first->second.num_photons+=partialop.num_photons;
263  //numphotons came from the sdps that contributed to the given hit. They should stack because the SDPs stack(which is why we are adding them in the first place). NOTE!! We will not do this when we combine records together for the detectors.
264  }
265  emp_e.first->second.unique();//Throw out duplicates.
266  Double_t energy=0.0;
267  for( auto en : emp_e.first->second ){
268  energy += en;
269  }
270  emp_o.first->second.energy=energy;
271  }
272  for(auto entry : bunched_hits){
273  npart.partial_hits.push_back(entry.second);
274  }
275  for(auto entry : bunched_ophits){
276  //This is where I should make the combined per detector records.
277  //This new record should find any existing record of the same OpDetNum and Time, and add the PEs, and Photons to that record, if it doesn't exist, make a new one.
278  npart.partial_ophits.push_back(entry.second);
279  }
280  }
281  }
282 
283  //Reding the bunched record is a good time to make the flattened record
284  /*
285  mf::LogDebug(__FUNCTION__)<<"EventBuffer for CalibrationTree contains "<<private_eventBuffer.eves.size()<<" eves, the first of which has "<<private_eventBuffer.eves.at(0).particles.size()<<" particles. The first of these particles has "<<private_eventBuffer.eves.at(0).particles.at(0).partial_hits.size()<<" partial hits.\n";
286  */
287  {//Forced local scope for a few variables. This is a good indication that the code I am writing doesn't belong here, but in another function.
288  int eve_pos=0;
292  for(auto eve : private_eventBuffer.eves){
293  if(!private_eventBuffer.eves.empty()){
294  fl.eve_x = eve.x_pos;
295  fl.eve_y = eve.y_pos;
296  fl.eve_z = eve.z_pos;
297  fl.eve_t = eve.t_pos;
298  fl.eve_trackid = eve.trackId;
299  fl.eve_pdgid = eve.pdgid;
300  fl.eve_index = eve_pos;
301  int part_pos=0;
302  for(auto part : eve.particles){
303  fl.part_x = part.x_pos;
304  fl.part_y = part.y_pos;
305  fl.part_z = part.z_pos;
306  fl.part_t = part.t_pos;
307  fl.part_trackid = part.trackId;
308  fl.part_pdgid = part.pdgid;
309  fl.part_iseve = part.isEve;
310  fl.part_index = part_pos;
311  { //keep partial_hit_pos local scope for the loop
312  fl.ophit_pes = 0.0;
313  fl.ophit_energy = 0.0;
314  fl.ophit_time = 0.0;
315  fl.ophit_width = 0.0;
316  fl.ophit_split = 0.0;
317  fl.ophit_index = 0.0;
318  fl.ophit_opchan = 0.0;
319  int partial_hit_pos=0;
320  for(auto partial_hit : part.partial_hits){
321  //This section is to fill the flat tree
322  fl.hit_charge = partial_hit.charge;
323  fl.hit_energy = partial_hit.energy;
324  fl.hit_time = partial_hit.time;
325  fl.hit_width = partial_hit.width;
326  fl.hit_split = partial_hit.split;
327  fl.hit_index = partial_hit.index;
328  fl.hit_wire = partial_hit.wire;
330  //Flat tree finished.
331  int ind = partial_hit.index;
332  std::vector<CalibTreeRecord::HitContributor>::iterator hc_pos = std::find_if(private_eventBuffer.hits.begin(), private_eventBuffer.hits.end(), [&ind](CalibTreeRecord::HitContributor& hc){ return ind==hc.index;});
333  if(hc_pos==private_eventBuffer.hits.end()){
334  private_eventBuffer.hits.emplace_back(ind);
335  if(hc_pos!=(private_eventBuffer.hits.end()-1)){
336  hc_pos=private_eventBuffer.hits.end()-1;
337  }
338  }
339  hc_pos->locations.emplace_back(eve_pos,part_pos,partial_hit_pos);
340  // private_eventBuffer.hits.emplace_back(eve_pos,part_pos,partial_hit_pos);
341  ++partial_hit_pos;
342  }
343  }
344  fl.hit_charge = 0.0;
345  fl.hit_energy = 0.0;
346  fl.hit_time = 0.0;
347  fl.hit_width = 0.0;
348  fl.hit_split = 0.0;
349  fl.hit_index = 0.0;
350  fl.hit_wire = 0.0;
351  {
352  int partial_ophit_pos=0;
353  for(auto partial_ophit : part.partial_ophits){
354  int ind = partial_ophit.index;
355  fl.ophit_pes = partial_ophit.pes;
356  fl.ophit_energy = partial_ophit.energy;
357  fl.ophit_time = partial_ophit.time;
358  fl.ophit_width = partial_ophit.width;
359  fl.ophit_split = partial_ophit.split;
360  fl.ophit_index = partial_ophit.index;
361  fl.ophit_opchan = partial_ophit.opchan;
363 
364  //Seg fault here. I suspect it is because of an empty private_eventBuffer.ophits
366  if(hc_pos==private_eventBuffer.ophits.end()){
367  private_eventBuffer.ophits.emplace_back(ind);
368  if(hc_pos!=(private_eventBuffer.ophits.end()-1)){
369  hc_pos=private_eventBuffer.ophits.end()-1;
370  }
371  }
372  hc_pos->locations.emplace_back(eve_pos,part_pos,partial_ophit_pos);
373  ++partial_ophit_pos;
374  }
375  }
376  ++part_pos;
377  }
378  }
379  ++eve_pos;
380  }
381  }
382  //Private event buffer eves are now done for this particle.
383  //Loop private event buffer to make the hit and ophit partial positions list.
385 
386 
387  private_CalibrationTree->Fill();
388 
389 
390 
391  } //end analyze
intermediate_table::iterator iterator
std::vector< HitContributor > hits
QList< Entry > entry
Handle< PROD > getHandle(SelectorBase const &) const
Definition: DataViewImpl.h:382
std::string string
Definition: nybbler.cc:12
RunNumber_t run() const
Definition: EventID.h:98
CalibTreeRecord::CalibTreeRecord private_eventBuffer
void add(const int &track_id, const std::string &gname)
std::vector< Handle< PROD > > getMany(SelectorBase const &selector=MatchAllSelector{}) const
Definition: DataViewImpl.h:479
std::vector< HitContributor > ophits
Detector simulation of raw signals on wires.
std::vector< art::Ptr< sim::OpDetDivRec > > priv_DivRecs
bool AddHit(detinfo::DetectorClocksData const &clockData, const art::Ptr< recob::Hit > hit, unsigned int &counter)
CalibTreeRecord::CalibTreeRecord private_eventPrep
MaybeLogger_< ELseverityLevel::ELsev_success, false > LogDebug
std::vector< EveRecord > eves
cet::LibraryManager dummy("noplugin")
EventNumber_t event() const
Definition: EventID.h:116
void fill_ptr_vector(std::vector< Ptr< T >> &ptrs, H const &h)
Definition: Ptr.h:297
SubRunNumber_t subRun() const
Definition: EventID.h:110
EventID id() const
Definition: Event.cc:34
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
void CalibrationTreeBuilder::CalibrationTreeBuilder::beginJob ( )
overridevirtual

Reimplemented from art::EDAnalyzer.

Definition at line 74 of file CalibrationTreeBuilder_module.cc.

74  {
75  private_CalibrationTree = private_service_tfs->make<TTree>("CalibrationRecordTree"," A TTree for storing hit and ophit values under their particles for each event in the simulation");
76  private_FlatCalibrationTree = private_service_tfs->make<TTree>("FlatCalibrationRecordTree"," A TTree for storing hit and ophit values with their particles for each event in the simulation");
78  private_FlatCalibrationRecord = private_FlatCalibrationTree->Branch("flat_event_records", &fl, "eve_x/D:eve_y:eve_z:eve_t:part_x:part_y:part_z:part_t:hit_charge:hit_energy:hit_time:hit_width:hit_split:ophit_pes:ophit_energy:ophit_time:ophit_width:ophit_split:hit_index/L:ophit_index:run/i:subrun:event_n:hit_wire:ophit_opchan:eve_index:part_index:eve_trackid/I:eve_pdgid:part_trackid:part_pdgid:part_iseve/O");
79 
80  }
CalibTreeRecord::CalibTreeRecord private_eventBuffer
art::ServiceHandle< art::TFileService > private_service_tfs
std::pair< std::vector< CalibTreeRecord::EveRecord >::iterator, bool > CalibrationTreeBuilder::CalibrationTreeBuilder::EmplaceEve ( const simb::MCParticle new_eve)
private

Definition at line 518 of file CalibrationTreeBuilder_module.cc.

518  {
519  std::vector<CalibTreeRecord::EveRecord>& eves = private_eventPrep.eves;
520  int eve_id = eve->TrackId();
521  std::vector<CalibTreeRecord::EveRecord>::iterator eve_ptr = std::find_if(eves.begin(), eves.end(), [&eve_id](CalibTreeRecord::EveRecord& eve_rec){return eve_id == eve_rec.trackId;});
522  if( eve_ptr != eves.end() ){ return std::make_pair(eve_ptr, false); }else{
524  tmpRec.trackId = eve_id;
525  tmpRec.generator = gf->get_gen(eve_id);
526  tmpRec.pdgid = eve->PdgCode();
527  mf::LogDebug(__FUNCTION__)<<"The eve now enetered into the particle record has PDGID "<<eve->PdgCode()<<"\n";
528  tmpRec.x_pos = eve->Position(0).X();
529  tmpRec.y_pos = eve->Position(0).Y();
530  tmpRec.z_pos = eve->Position(0).Z(); tmpRec.t_pos = eve->Position(0).T();
531  tmpRec.t_pos = eve->Position(0).T();
532  eves.push_back(tmpRec);
533  if(eve_ptr == (eves.end()-1) ){
534  return std::make_pair(eve_ptr, true);
535  }else{
536  return std::make_pair(eves.end()-1, true);
537  }
538  }
539  }
intermediate_table::iterator iterator
CalibTreeRecord::CalibTreeRecord private_eventPrep
MaybeLogger_< ELseverityLevel::ELsev_success, false > LogDebug
std::vector< EveRecord > eves
std::pair< std::vector< CalibTreeRecord::ParticleRecord >::iterator, bool > CalibrationTreeBuilder::CalibrationTreeBuilder::EmplaceParticle ( const simb::MCParticle new_part)
private

Definition at line 487 of file CalibrationTreeBuilder_module.cc.

487  {
488  int track_id = part->TrackId();
489  int eve_id = PIS->TrackIdToEveTrackId(part->TrackId());
490  const simb::MCParticle* eve = PIS->TrackIdToParticle_P(eve_id);
491  //This search can be put into EmplaceEve
492  //std::vector<EvenRecord::ParticleRecord>::iterator = std::find_if(private_eventPrep.eves.begin(), private_eventPrep.eves.end(), [&eve_id](CalibTreeRecord::EveRecord& eve_rec){ return eve_id==eve_rec.trackId;});
494  std::vector<CalibTreeRecord::ParticleRecord>::iterator part_ptr = std::find_if(eve_marker.first->particles.begin(), eve_marker.first->particles.end(), [&track_id](CalibTreeRecord::ParticleRecord& part_rec){return track_id == part_rec.trackId;});
495  if(part_ptr!=eve_marker.first->particles.end()){
496  return std::make_pair(part_ptr, false);
497  }else{
499  tmpRec.isEve = (track_id==eve_id)?true:false; //The conditional operator here is unneccessary, but I think it is more clear to include it.
500  tmpRec.trackId = track_id;
501  tmpRec.pdgid = part->PdgCode();
502  tmpRec.dP = ( part->EndMomentum().P() - part->Momentum(0).P() );
503  tmpRec.dE = ( part->EndMomentum().E() - part->Momentum(0).E() );
504  mf::LogDebug(__FUNCTION__)<<"The particle now enetered into the particle record has PDGID "<<part->PdgCode()<<"\n";
505  tmpRec.x_pos = part->Position(0).X();
506  tmpRec.y_pos = part->Position(0).Y();
507  tmpRec.z_pos = part->Position(0).Z();
508  tmpRec.t_pos = part->Position(0).T();
509  eve_marker.first->particles.push_back(tmpRec);
510  if(part_ptr == (eve_marker.first->particles.end()-1) ){
511  return std::make_pair(part_ptr, true); //If the vector hasn't been moved (can happen with expanding), just return what we have made).
512  }else{//The vector has moved. Return the right answer (this could honestly be the only return, the previous being unnecessary).
513  return std::make_pair(eve_marker.first->particles.end()-1, true);
514  }
515  }//end else (part_ptr!=eve_marker.first->particles.end());
516  }//end EmplaceParticle
intermediate_table::iterator iterator
const simb::MCParticle * TrackIdToParticle_P(int id) const
std::pair< std::vector< CalibTreeRecord::EveRecord >::iterator, bool > EmplaceEve(const simb::MCParticle *new_eve)
art::ServiceHandle< cheat::ParticleInventoryService > PIS
MaybeLogger_< ELseverityLevel::ELsev_success, false > LogDebug
void CalibrationTreeBuilder::CalibrationTreeBuilder::endJob ( )
overridevirtual

Reimplemented from art::EDAnalyzer.

Definition at line 394 of file CalibrationTreeBuilder_module.cc.

394  {
395  }
const art::Ptr< sim::OpDetDivRec > CalibrationTreeBuilder::CalibrationTreeBuilder::FindDivRec ( int const &  opDetNum) const
private

Definition at line 545 of file CalibrationTreeBuilder_module.cc.

546  {
548  for(size_t detnum = 0; detnum < priv_DivRecs.size(); ++detnum){
549  if(priv_DivRecs.at(detnum)->OpDetNum() == opDetNum)
550  opDet = priv_DivRecs.at(detnum);
551  }
552  if(!opDet)
553  {
554  throw cet::exception("CalibTreeBuilder") << "No sim:: OpDetDivRec corresponding "
555  << "to opDetNum: " << opDetNum << "\n";
556  }
557  return opDet;
558  }
std::vector< art::Ptr< sim::OpDetDivRec > > priv_DivRecs
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
const std::vector< sim::SDP > CalibrationTreeBuilder::CalibrationTreeBuilder::OpHitToChannelWeightedSimSDPs ( art::Ptr< recob::OpHit > const &  opHit_P) const
private

Definition at line 564 of file CalibrationTreeBuilder_module.cc.

565  {
566  const double fDelay = PBS->GetDelay();
567  std::vector< sim::SDP > retVec;
568  double fPeakTime = opHit_P->PeakTime();
569  double fWidth = opHit_P->Width();
570  UInt_t fChan = opHit_P->OpChannel();
572  int fDet = geom->OpDetFromOpChannel(fChan);
573  //I should use the timing service for these time conversions.
574  sim::OpDetBacktrackerRecord::timePDclock_t start_time = ((fPeakTime- fWidth)*1000.0)-fDelay;
575  sim::OpDetBacktrackerRecord::timePDclock_t end_time = ((fPeakTime+ fWidth)*1000.0)-fDelay;
576  if(start_time > end_time){throw;}//This is bad. Give a reasonable error message here, and use cet::except
577 
578  //BUG!!!fGeom->OpDetFromOpChannel(channel)
580  try{
581  fBTR = PBS->FindOpDetBTR(fDet);
582  }catch(cet::exception &e){
583  if(true){
584  std::vector<sim::SDP> ret;
585  return ret;
586  }else{
587  throw e;
588  }
589  }
590  const std::vector<std::pair<double, std::vector<sim::SDP>> >& timeSDPMap
591  = fBTR->timePDclockSDPsMap(); //Not guranteed to be sorted.
592  art::Ptr<sim::OpDetDivRec> div_rec = this->FindDivRec(fDet);//This is an OpDetDivRec collected from this BTR.
593 
594 
595  std::vector<const std::pair<double, std::vector<sim::SDP>>*> timePDclockSDPMap_SortedPointers;
596  std::vector<const sim::OpDet_Time_Chans*> div_rec_SortedPointers;
597  for ( auto& pair : timeSDPMap ){ timePDclockSDPMap_SortedPointers.push_back(&pair); }
598  auto pairSort = [](auto& a, auto& b) { return a->first < b->first ; } ;
599  if( !std::is_sorted( timePDclockSDPMap_SortedPointers.begin(), timePDclockSDPMap_SortedPointers.end(), pairSort)){
600  std::sort(timePDclockSDPMap_SortedPointers.begin(), timePDclockSDPMap_SortedPointers.end(), pairSort);
601  }
602 
603  //This section is a hack to make comparisons work right.
604  std::vector<sim::SDP> dummyVec;
605  std::pair<double, std::vector<sim::SDP>> start_timePair = std::make_pair(start_time, dummyVec);
606  std::pair<double, std::vector<sim::SDP>> end_timePair = std::make_pair(end_time, dummyVec);
607  auto start_timePair_P = &start_timePair;
608  auto end_timePair_P = &end_timePair;
609 
610  //First interesting iterator.
611  auto map_pdsdp_itr = std::lower_bound(timePDclockSDPMap_SortedPointers.begin(), timePDclockSDPMap_SortedPointers.end(), start_timePair_P, pairSort);
612  //Last interesting iterator.
613  auto map_pdsdp_last = std::upper_bound(map_pdsdp_itr, timePDclockSDPMap_SortedPointers.end(), end_timePair_P, pairSort);
614 
615  //retvec.push_back(map_pdsdp_first.second[0]);
616  //This code screams fragile. Really. If you read this, feel free to clean up the methods.
617  for(auto& sdp_time_pair = map_pdsdp_itr; sdp_time_pair != map_pdsdp_last; ++sdp_time_pair){
618  //cut div_rec by time
619  auto time=(*sdp_time_pair)->first;
620  for(auto& sdp : (*sdp_time_pair)->second){
621  // auto sdp = (*sdp_time_pair)->second; //This is a vector
622  //auto slice = div_rec->GetSlice(start_time, end_time);
623  auto time_divr_pair = div_rec->FindClosestTimeChan(time);
624  auto time_divr=time_divr_pair.first;
625  auto time_divr_found=time_divr_pair.second;
626  sim::SDP tmp = sdp;
627  if( time_divr_found && time_divr->time == time){ //This will break if time_divr is the end value of the vector
628  tmp.energy = time_divr->GetFrac(fChan) * tmp.energy;
629  // ((*map_divrec_itr)->DivChans.eScaleFrac(fChan))*tmp.energy;
630  }else{
631  // std::cout<<"I am having trouble reconciling "<<time<<" and "<<time_divr->time<<"\n";
632  tmp.energy = 0; //This does not corespond to an energy detected. A photon may or may not have been incident on the detector, but it didn't get picked up.
633  }
634  retVec.emplace_back(tmp);
635 
636  }
637  }
638 
639  //while( map_pdsdp_itr != map_pdsdp_last /*&& map_divrec_itr != map_divrec_last*/){ //Stepping through.
640  /*for(auto const& sdp : (*map_pdsdp_itr)->second){
641  sim::SDP tmp = sdp;
642  tmp.energy = ((*map_divrec_itr)->DivChans.eScaleFrac(fChan))*tmp.energy;
643  // second[fChan].tick_photons_frac)*tmp.energy;
644  retVec.emplace_back(tmp);
645  }
646  ++map_pdsdp_itr;
647  ++map_divrec_itr;
648  }*/
649 
650  return retVec;
651  }
art::Ptr< sim::OpDetBacktrackerRecord > FindOpDetBTR(int const &opDetNum)
double PeakTime() const
Definition: OpHit.h:64
unsigned int OpDetFromOpChannel(int opChannel) const
Convert unique channel to detector number.
const double e
double Width() const
Definition: OpHit.h:66
const double a
double timePDclock_t
Type for iTimePDclock tick used in the interface.
string tmp
Definition: languages.py:63
float energy
energy deposited by ionization
const art::Ptr< sim::OpDetDivRec > FindDivRec(int const &opDetNum) const
static bool * b
Definition: config.cpp:1043
art::ServiceHandle< cheat::PhotonBackTrackerService > PBS
int OpChannel() const
Definition: OpHit.h:62
timePDclockSDPs_t const & timePDclockSDPsMap() const
Returns all the deposited energy information as stored.
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
std::pair< OpDetDivRec::Time_Chans_t::const_iterator, bool > FindClosestTimeChan(OpDet_Time_Chans::stored_time_t pdTime) const
Definition: OpDetDivRec.cxx:73
void CalibrationTreeBuilder::CalibrationTreeBuilder::PrepDivRec ( const art::Event evt)
private

Definition at line 653 of file CalibrationTreeBuilder_module.cc.

654 {
655  if( 0 ){ return;} //Insert check for DivRecs here, or don't use validHandle below.
656  auto const& divrecHandle = evt.getValidHandle <std::vector<sim::OpDetDivRec>>(fWavLabel);
657  if(divrecHandle.failedToGet()){
658  return;
659  }
660  art::fill_ptr_vector(priv_DivRecs, divrecHandle);
661  auto compareDivReclambda = [](art::Ptr<sim::OpDetDivRec> a, art::Ptr<sim::OpDetDivRec> b) {return(a->OpDetNum() < b->OpDetNum());};
662  if (!std::is_sorted(priv_DivRecs.begin(), priv_DivRecs.end(), compareDivReclambda))
663  std::sort(priv_DivRecs.begin(), priv_DivRecs.end(), compareDivReclambda);
664 
665 }
const double a
ValidHandle< PROD > getValidHandle(InputTag const &tag) const
Definition: DataViewImpl.h:441
int OpDetNum() const
Definition: OpDetDivRec.h:102
std::vector< art::Ptr< sim::OpDetDivRec > > priv_DivRecs
static bool * b
Definition: config.cpp:1043
void fill_ptr_vector(std::vector< Ptr< T >> &ptrs, H const &h)
Definition: Ptr.h:297

Member Data Documentation

art::ServiceHandle<cheat::BackTrackerService> CalibrationTreeBuilder::CalibrationTreeBuilder::BTS
private

Definition at line 111 of file CalibrationTreeBuilder.h.

flatfiller CalibrationTreeBuilder::CalibrationTreeBuilder::fl
private

Definition at line 131 of file CalibrationTreeBuilder.h.

art::InputTag CalibrationTreeBuilder::CalibrationTreeBuilder::fWavLabel
private

Definition at line 117 of file CalibrationTreeBuilder.h.

art::ServiceHandle<geo::Geometry> CalibrationTreeBuilder::CalibrationTreeBuilder::GS
private

Definition at line 113 of file CalibrationTreeBuilder.h.

art::ServiceHandle<cheat::PhotonBackTrackerService> CalibrationTreeBuilder::CalibrationTreeBuilder::PBS
private

Definition at line 112 of file CalibrationTreeBuilder.h.

art::ServiceHandle<cheat::ParticleInventoryService> CalibrationTreeBuilder::CalibrationTreeBuilder::PIS
private

Definition at line 110 of file CalibrationTreeBuilder.h.

std::vector<art::Ptr<sim::OpDetDivRec> > CalibrationTreeBuilder::CalibrationTreeBuilder::priv_DivRecs
mutableprivate

Definition at line 133 of file CalibrationTreeBuilder.h.

TBranch* CalibrationTreeBuilder::CalibrationTreeBuilder::private_CalibrationRecord
private

Definition at line 123 of file CalibrationTreeBuilder.h.

TTree* CalibrationTreeBuilder::CalibrationTreeBuilder::private_CalibrationTree
private

Definition at line 121 of file CalibrationTreeBuilder.h.

CalibTreeRecord::CalibTreeRecord CalibrationTreeBuilder::CalibrationTreeBuilder::private_eventBuffer
private

Definition at line 129 of file CalibrationTreeBuilder.h.

CalibTreeRecord::CalibTreeRecord CalibrationTreeBuilder::CalibrationTreeBuilder::private_eventPrep
private

Definition at line 130 of file CalibrationTreeBuilder.h.

TBranch* CalibrationTreeBuilder::CalibrationTreeBuilder::private_FlatCalibrationRecord
private

Definition at line 124 of file CalibrationTreeBuilder.h.

TTree* CalibrationTreeBuilder::CalibrationTreeBuilder::private_FlatCalibrationTree
private

Definition at line 122 of file CalibrationTreeBuilder.h.

art::InputTag CalibrationTreeBuilder::CalibrationTreeBuilder::private_HitLabel
private

Definition at line 115 of file CalibrationTreeBuilder.h.

art::InputTag CalibrationTreeBuilder::CalibrationTreeBuilder::private_OpHitLabel
private

Definition at line 116 of file CalibrationTreeBuilder.h.

art::ServiceHandle<art::TFileService> CalibrationTreeBuilder::CalibrationTreeBuilder::private_service_tfs
private

Definition at line 119 of file CalibrationTreeBuilder.h.


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