#include <CalibrationTreeBuilder.h>

Inheritance diagram for CalibrationTreeBuilder::CalibrationTreeBuilder:

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

Observer &	operator= (Observer const &)=delete

Observer &	operator= (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, bool >	EmplaceEve (const simb::MCParticle *new_eve)

std::pair< std::vector< CalibTreeRecord::ParticleRecord >::iterator, bool >	EmplaceParticle (const simb::MCParticle *new_part)

const art::Ptr< sim::OpDetDivRec >	FindDivRec (int const &opDetNum) const

const std::vector< sim::SDP >	OpHitToChannelWeightedSimSDPs (art::Ptr< recob::OpHit > const &opHit_P) const

Private Attributes
art::ServiceHandle< cheat::ParticleInventoryService >	PIS

art::ServiceHandle< cheat::BackTrackerService >	BTS

art::ServiceHandle< cheat::PhotonBackTrackerService >	PBS

art::ServiceHandle< geo::Geometry >	GS

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< TriggerResults >	getTriggerResults (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
ConsumesCollector &	consumesCollector ()

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.

     :EDAnalyzer(pSet),
     private_HitLabel(pSet.get<art::InputTag>("HitLabel","gaushit")),
     private_OpHitLabel(pSet.get<art::InputTag>("OpHitLabel","ophit")),
     fWavLabel(pSet.get<art::InputTag>("WavLabel", "opdigi"))
   {  }

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.

                                                                                         {
     art::ServiceHandle<geo::Geometry> geom;
     std::vector<const sim::IDE*> ides_ps = BTS->HitToSimIDEs_Ps(clockData, hit);
     std::vector<std::pair<int, CalibTreeRecord::PartialHit>> track_partials;
     Double_t total_charge = 0.0;
     if(ides_ps.empty()){return false;}
     for(auto ide_p : ides_ps){
       int track_id = ide_p->trackID;
       CalibTreeRecord::PartialHit tmp;
       total_charge += ide_p->numElectrons;
       //        tmp.charge = ide_p->numElectrons;
       tmp.num_electrons = ide_p->numElectrons;
       tmp.time   = hit->PeakTime();
       tmp.width  = hit->SigmaPeakTime();
       tmp.split  = 0.0;
       tmp.wire   = hit->Channel();
       tmp.is_collection_wire = (geom->SignalType(tmp.wire)==geo::kCollection)?(true):(false);
       tmp.energy = (tmp.is_collection_wire)?(ide_p->energy):(0.0);
       tmp.index  = counter;
       track_partials.push_back(std::make_pair(track_id, tmp));
     }
     for(auto track_partial : track_partials){
       track_partial.second.split  = track_partial.second.num_electrons / total_charge; //What percentage of the hit is in this partial.
       track_partial.second.charge = hit->Integral()*track_partial.second.split; //Charge now normalized by deposition;
       //Add partial to CalibTreeRecord;
       const simb::MCParticle* part = PIS->TrackIdToParticle_P(track_partial.first);
       //Emplace Particle. (Emplace Eve will be handeled in emplace particle
       std::pair<std::vector<CalibTreeRecord::ParticleRecord>::iterator, bool> part_marker = EmplaceParticle(part);
       part_marker.first->partial_hits.push_back(track_partial.second); //Add post normalized partial hit.
       //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.
       mf::LogDebug(__FUNCTION__)<<"This partialhit is being added to the record with charge "<<track_partial.second.charge<<"\n";
     }//end for track partial.
 
     return true;
   }

bool CalibrationTreeBuilder::CalibrationTreeBuilder::AddHit	(	const art::Ptr< recob::OpHit >	hit,
		unsigned int &	counter
	)

private

Definition at line 440 of file CalibrationTreeBuilder_module.cc.

                                                                                         {
     {//make each loop local so the vectors don't persist
 //      try{
         const std::vector< sim::SDP> sdps = this->OpHitToChannelWeightedSimSDPs(hit);
 //      }catch(const std::exception &e){
 //        return false;
 //      }
         std::vector<std::pair<int, CalibTreeRecord::PartialOpHit>> track_partials;
         Double_t total_charge = 0.0;
         if(sdps.empty()){return false;}
         for(auto sdp : sdps){
           int track_id = sdp.trackID;
           CalibTreeRecord::PartialOpHit tmp;
           total_charge += sdp.numPhotons;
           //        tmp.pes = //Not used yet.
           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.
           tmp.energy = sdp.energy;
           //tmp.energy = 0.0;// sdp.energy;
           tmp.time   = hit->PeakTime();
           tmp.width  = hit->Width();
           tmp.split  = 0.0;
           tmp.opchan = hit->OpChannel();
           tmp.opdet  = GS->OpDetFromOpChannel(hit->OpChannel());
           tmp.index  = counter;
           track_partials.push_back(std::make_pair(track_id, tmp));
         }
         for(auto track_partial : track_partials){
           track_partial.second.split  = track_partial.second.num_photons / total_charge; //What percentage of the hit is in this partial.
           track_partial.second.pes = hit->PE()*track_partial.second.split; //PEs now normalized by deposition;
           //Add partial to CalibTreeRecord;
           const simb::MCParticle* part = PIS->TrackIdToParticle_P(track_partial.first);
           //Emplace Particle. (Emplace Eve will be handeled in emplace particle
           std::pair<std::vector<CalibTreeRecord::ParticleRecord>::iterator, bool> part_marker = EmplaceParticle(part);
           part_marker.first->partial_ophits.push_back(track_partial.second); //Add post normalized partial hit.
           //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.
         }//end for track partial.
 
       }
       return true;
 //    }
 //    catch(...){
 //      return false;
 //    }
   }

void CalibrationTreeBuilder::CalibrationTreeBuilder::analyze ( const art::Event & evt )

overridevirtual

Definition at line 83 of file CalibrationTreeBuilder_module.cc.

                                                           {
     //eprV.clear();
     private_eventPrep.Clear();
     private_eventBuffer.Clear();
     priv_DivRecs.clear();
     this->PrepDivRec(evt);
 
 
     //Set up BackTracking //Now in header
     /*    art::ServiceHandle<cheat::ParticleInventoryService> PIS;
           art::ServiceHandle<cheat::BackTrackerService> BTS;
           art::ServiceHandle<cheat::PhotonBackTrackerService> PBS;*/
     //Get a list of generator names.
     auto const mcHandles = evt.getMany<std::vector<simb::MCTruth>>();
     std::vector< std::pair<int, std::string>> track_id_to_label;
 
     for( auto const& mcHandle : mcHandles ){
       const std::string& sModuleLabel = mcHandle.provenance()->moduleLabel();
       art::FindManyP<simb::MCParticle> findMCParts(mcHandle, evt, "largeant");
       std::vector<art::Ptr<simb::MCParticle> > mcParts = findMCParts.at(0);
       for( const art::Ptr<simb::MCParticle> ptr : mcParts){
         int track_id = ptr->TrackId();
         gf->add(track_id, sModuleLabel);
       }
     }
 
     std::vector<art::Ptr<recob::Hit>> hitList;
     auto hitHandle = evt.getHandle<std::vector<recob::Hit>>(private_HitLabel);
     if ( hitHandle )
       art::fill_ptr_vector(hitList, hitHandle);
     std::vector<art::Ptr<recob::OpHit>> opHitList;
     auto opHitHandle = evt.getHandle<std::vector<recob::OpHit>>(private_OpHitLabel);
     if ( opHitHandle )
       art::fill_ptr_vector(opHitList, opHitHandle);
 
     //private_eventBuffer is the event record
     private_eventBuffer.bunch_hits=true;
     private_eventBuffer.run=evt.id().run();
     private_eventBuffer.subrun=evt.id().subRun();
     private_eventBuffer.event=evt.id().event();
     //private_eventPrep is the event record
     private_eventPrep.bunch_hits=true;
     private_eventPrep.run=evt.id().run();
     private_eventPrep.subrun=evt.id().subRun();
     private_eventPrep.event=evt.id().event();
 
     auto const clockData = art::ServiceHandle<detinfo::DetectorClocksService const>()->DataFor(evt);
     unsigned int hitCounter=0;
     for ( auto const& hit : hitList){
       hitCounter++;
       if(!AddHit(clockData, hit, hitCounter)){
         mf::LogDebug(__FUNCTION__)<<"Hit "<<hitCounter<<" failed to add to the record. Is this hit caused by noise?";
       }
     }//end for hit in hitlist
 
     hitCounter=0;
     for ( auto const& hit : opHitList ){
       ++hitCounter;
       if(!AddHit(hit, hitCounter)){
         mf::LogDebug(__FUNCTION__)<<"OpHit "<<hitCounter<<" failed to add to the record. Is this hit caused by noise?";
       }
     } //Add ophits to hitmap.
 
     //Buffer normalized per hit as it is built.
     //Make the event record
     /*
        if(!private_eventPrep.eves.empty()){
        mf::LogDebug(__FUNCTION__)<<"EventBuffer for CalibrationTree contains "<<private_eventPrep.eves.size()<<" eves.\n";
        if(!private_eventPrep.eves.at(0).particles.empty()){
        mf::LogDebug(__FUNCTION__)<<"The first particle has "<<private_eventPrep.eves.at(0).particles.size()<<" particles.\n";
        if(!private_eventPrep.eves.at(0).particles.at(0).partial_hits.empty()){
        mf::LogDebug(__FUNCTION__)<<"The first of these particles has "<<private_eventPrep.eves.at(0).particles.at(0).partial_hits.size()<<" partial hits.\n";
        }else{
        mf::LogDebug(__FUNCTION__)<<"The particle contains no hits.\n";
        }
        }else{
        mf::LogDebug(__FUNCTION__)<<"Eve contains no particles.\n";
        }
        }else{
        mf::LogDebug(__FUNCTION__)<<"EventBuffer for CalibrationTree contains no eves.\n";
        }
        */
 
     {//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.
       int eve_pos=0;
       for(auto eve : private_eventPrep.eves){
         if(!private_eventPrep.eves.empty()){ //This is not needed and is redundant.
           int part_pos=0;
           for(auto part : eve.particles){
             { //keep partial_hit_pos local scope for the loop
               int partial_hit_pos=0;
               for(auto partial_hit : part.partial_hits){
                 int ind = partial_hit.index;
                 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;});
                 if(hc_pos==private_eventPrep.hits.end()){
                   private_eventPrep.hits.emplace_back(ind);
                   if(hc_pos!=(private_eventPrep.hits.end()-1)){
                     hc_pos=private_eventPrep.hits.end()-1;
                   }
                 }
                 hc_pos->locations.emplace_back(eve_pos,part_pos,partial_hit_pos);
                 //                private_eventPrep.hits.emplace_back(eve_pos,part_pos,partial_hit_pos);
                 ++partial_hit_pos;
               }
             }
             {
               int partial_ophit_pos=0;
               for(auto partial_ophit : part.partial_ophits){
                 int ind = partial_ophit.index;
                 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;});
                 if(hc_pos==private_eventPrep.ophits.end()){
                   private_eventPrep.ophits.emplace_back(ind);
                   if(hc_pos!=(private_eventPrep.ophits.end()-1)){
                     hc_pos=private_eventPrep.ophits.end()-1;
                   }
                 }
                 hc_pos->locations.emplace_back(eve_pos,part_pos,partial_ophit_pos);
                 ++partial_ophit_pos;
               }
             }
             ++part_pos;
           }
         }
         ++eve_pos;
       }
     }
 
     //This is where the partial hit bunching should happen.
     //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.
     //Now we loop the private_eventPrep and build private_eventBuffer.
 
     if( private_eventPrep.bunch_hits==true){
       for(size_t eve_num=0; eve_num<private_eventPrep.eves.size(); ++eve_num ){
         auto eve = private_eventPrep.eves.at(eve_num);
         private_eventBuffer.eves.push_back(eve); //We are gettign them in order, we are filling them in order. This should not break anything.
         //private_eventBuffer.eves.at(eve_num).particles.clear();//clear the particles list from the copeid eve. We will now refill it.
         for(size_t part_num=0; part_num<eve.particles.size(); ++part_num ){
           auto part = eve.particles.at(part_num);
           auto& npart = private_eventBuffer.eves.at(eve_num).particles.at(part_num);
           npart.partial_hits.clear();
           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.
           //add part to particle buffer.
           std::map<int64_t, CalibTreeRecord::PartialHit> bunched_hits;
           std::map<int64_t, std::list<Double_t>> remembered_energies;
           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);
             auto emp_h = bunched_hits.emplace(std::make_pair(partialhit.index,partialhit));
             std::list<Double_t> dummy;
             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.
             if(emp_h.second==true){ //emp_e.second and emp_h.second are always the same.
               emp_e.first->second.emplace_back(partialhit.energy);
               emp_h.first->second=partialhit;
             }else{
               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";}
               emp_e.first->second.emplace_back(partialhit.energy);
               emp_h.first->second.charge+=partialhit.charge;
               emp_h.first->second.split+=partialhit.split;
               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).
             }
             emp_e.first->second.unique();//Throw out duplicates.
             Double_t energy=0.0;
             for( auto en : emp_e.first->second){
               energy += en;
             }
             emp_h.first->second.energy=energy;
           }
           remembered_energies.clear();
           std::map<int64_t, CalibTreeRecord::PartialOpHit> bunched_ophits;
           for(auto partialop : part.partial_ophits){
             auto emp_o = bunched_ophits.emplace(std::make_pair(partialop.index,partialop));
             std::list<Double_t> dummy;
             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.
             if(emp_o.second==true){
               emp_e.first->second.emplace_back(partialop.energy);
               emp_o.first->second=partialop;
             }else{
               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";}
               emp_e.first->second.emplace_back(partialop.energy);
               emp_o.first->second.pes+=partialop.pes;
               emp_o.first->second.split+=partialop.split;
               emp_o.first->second.num_photons+=partialop.num_photons;
               //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.
             }
             emp_e.first->second.unique();//Throw out duplicates.
             Double_t energy=0.0;
             for( auto en : emp_e.first->second ){
               energy += en;
             }
             emp_o.first->second.energy=energy;
           }
           for(auto entry : bunched_hits){
             npart.partial_hits.push_back(entry.second);
           }
           for(auto entry : bunched_ophits){
             //This is where I should make the combined per detector records.
             //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.
             npart.partial_ophits.push_back(entry.second);
           }
         }
       }
 
       //Reding the bunched record is a good time to make the flattened record
       /*
          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";
          */
       {//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.
         int eve_pos=0;
         fl.run           = private_eventBuffer.run;
         fl.subrun        = private_eventBuffer.subrun;
         fl.event_n       = private_eventBuffer.event;
         for(auto eve : private_eventBuffer.eves){
           if(!private_eventBuffer.eves.empty()){
             fl.eve_x         = eve.x_pos;
             fl.eve_y         = eve.y_pos;
             fl.eve_z         = eve.z_pos;
             fl.eve_t         = eve.t_pos;
             fl.eve_trackid   = eve.trackId;
             fl.eve_pdgid     = eve.pdgid;
             fl.eve_index     = eve_pos;
             int part_pos=0;
             for(auto part : eve.particles){
               fl.part_x        = part.x_pos;
               fl.part_y        = part.y_pos;
               fl.part_z        = part.z_pos;
               fl.part_t        = part.t_pos;
               fl.part_trackid  = part.trackId;
               fl.part_pdgid    = part.pdgid;
               fl.part_iseve    = part.isEve;
               fl.part_index    = part_pos;
               { //keep partial_hit_pos local scope for the loop
                 fl.ophit_pes     = 0.0;
                 fl.ophit_energy  = 0.0;
                 fl.ophit_time    = 0.0;
                 fl.ophit_width   = 0.0;
                 fl.ophit_split   = 0.0;
                 fl.ophit_index   = 0.0;
                 fl.ophit_opchan  = 0.0;
                 int partial_hit_pos=0;
                 for(auto partial_hit : part.partial_hits){
                   //This section is to fill the flat tree
                   fl.hit_charge    = partial_hit.charge;
                   fl.hit_energy    = partial_hit.energy;
                   fl.hit_time      = partial_hit.time;
                   fl.hit_width     = partial_hit.width;
                   fl.hit_split     = partial_hit.split;
                   fl.hit_index     = partial_hit.index;
                   fl.hit_wire      = partial_hit.wire;
                   private_FlatCalibrationTree->Fill();
                   //Flat tree finished.
                   int ind = partial_hit.index;
                   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;});
                   if(hc_pos==private_eventBuffer.hits.end()){
                     private_eventBuffer.hits.emplace_back(ind);
                     if(hc_pos!=(private_eventBuffer.hits.end()-1)){
                       hc_pos=private_eventBuffer.hits.end()-1;
                     }
                   }
                   hc_pos->locations.emplace_back(eve_pos,part_pos,partial_hit_pos);
                   //                private_eventBuffer.hits.emplace_back(eve_pos,part_pos,partial_hit_pos);
                   ++partial_hit_pos;
                 }
               }
               fl.hit_charge    = 0.0;
               fl.hit_energy    = 0.0;
               fl.hit_time      = 0.0;
               fl.hit_width     = 0.0;
               fl.hit_split     = 0.0;
               fl.hit_index     = 0.0;
               fl.hit_wire      = 0.0;
               {
                 int partial_ophit_pos=0;
                 for(auto partial_ophit : part.partial_ophits){
                   int ind = partial_ophit.index;
                   fl.ophit_pes     = partial_ophit.pes;
                   fl.ophit_energy  = partial_ophit.energy;
                   fl.ophit_time    = partial_ophit.time;
                   fl.ophit_width   = partial_ophit.width;
                   fl.ophit_split   = partial_ophit.split;
                   fl.ophit_index   = partial_ophit.index;
                   fl.ophit_opchan  = partial_ophit.opchan;
                   private_FlatCalibrationTree->Fill();
 
                   //Seg fault here. I suspect it is because of an empty private_eventBuffer.ophits
                   std::vector<CalibTreeRecord::HitContributor>::iterator hc_pos = std::find_if(private_eventBuffer.ophits.begin(), private_eventBuffer.ophits.end(), [&ind](CalibTreeRecord::HitContributor& hc){ return ind==hc.index;});
                   if(hc_pos==private_eventBuffer.ophits.end()){
                     private_eventBuffer.ophits.emplace_back(ind);
                     if(hc_pos!=(private_eventBuffer.ophits.end()-1)){
                       hc_pos=private_eventBuffer.ophits.end()-1;
                     }
                   }
                   hc_pos->locations.emplace_back(eve_pos,part_pos,partial_ophit_pos);
                   ++partial_ophit_pos;
                 }
               }
               ++part_pos;
             }
           }
           ++eve_pos;
         }
       }
       //Private event buffer eves are now done for this particle.
       //Loop private event buffer to make the hit and ophit partial positions list.
     }else{ private_eventBuffer=private_eventPrep;}
 
 
     private_CalibrationTree->Fill();
 
 
 
   } //end analyze

void CalibrationTreeBuilder::CalibrationTreeBuilder::beginJob ( )

overridevirtual

Reimplemented from art::EDAnalyzer.

Definition at line 74 of file CalibrationTreeBuilder_module.cc.

                                        {
     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");
     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");
     private_CalibrationRecord = private_CalibrationTree->Branch("event_records", &private_eventBuffer );
     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");
 
   }

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.

                                                                                                                         {
     std::vector<CalibTreeRecord::EveRecord>& eves = private_eventPrep.eves;
     int eve_id = eve->TrackId();
     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;});
     if( eve_ptr != eves.end() ){ return std::make_pair(eve_ptr, false); }else{
       CalibTreeRecord::EveRecord tmpRec;
       tmpRec.trackId = eve_id;
       tmpRec.generator = gf->get_gen(eve_id);
       tmpRec.pdgid = eve->PdgCode();
       mf::LogDebug(__FUNCTION__)<<"The eve now enetered into the particle record has PDGID "<<eve->PdgCode()<<"\n";
       tmpRec.x_pos   = eve->Position(0).X();
       tmpRec.y_pos   = eve->Position(0).Y();
       tmpRec.z_pos   = eve->Position(0).Z();      tmpRec.t_pos   = eve->Position(0).T();
       tmpRec.t_pos   = eve->Position(0).T();
       eves.push_back(tmpRec);
       if(eve_ptr == (eves.end()-1) ){
         return std::make_pair(eve_ptr, true);
       }else{
         return std::make_pair(eves.end()-1, true);
       }
     }
   }

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.

                                                                                                                                    {
     int track_id = part->TrackId();
     int eve_id   = PIS->TrackIdToEveTrackId(part->TrackId());
     const simb::MCParticle* eve = PIS->TrackIdToParticle_P(eve_id);
     //This search can be put into EmplaceEve
     //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;});
     std::pair<std::vector<CalibTreeRecord::EveRecord>::iterator, bool> eve_marker = EmplaceEve(eve);
     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;});
     if(part_ptr!=eve_marker.first->particles.end()){
       return std::make_pair(part_ptr, false);
     }else{
       CalibTreeRecord::ParticleRecord tmpRec;
       tmpRec.isEve   = (track_id==eve_id)?true:false; //The conditional operator here is unneccessary, but I think it is more clear to include it.
       tmpRec.trackId = track_id;
       tmpRec.pdgid   = part->PdgCode();
       tmpRec.dP      = ( part->EndMomentum().P() - part->Momentum(0).P() );
       tmpRec.dE      = ( part->EndMomentum().E() - part->Momentum(0).E() );
       mf::LogDebug(__FUNCTION__)<<"The particle now enetered into the particle record has PDGID "<<part->PdgCode()<<"\n";
       tmpRec.x_pos   = part->Position(0).X();
       tmpRec.y_pos   = part->Position(0).Y();
       tmpRec.z_pos   = part->Position(0).Z();
       tmpRec.t_pos   = part->Position(0).T();
       eve_marker.first->particles.push_back(tmpRec);
       if(part_ptr == (eve_marker.first->particles.end()-1) ){
         return std::make_pair(part_ptr, true); //If the vector hasn't been moved (can happen with expanding), just return what we have made).
       }else{//The vector has moved. Return the right answer (this could honestly be the only return, the previous being unnecessary).
         return std::make_pair(eve_marker.first->particles.end()-1, true);
       }
     }//end else (part_ptr!=eve_marker.first->particles.end());
   }//end EmplaceParticle

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.

   {
     art::Ptr< sim::OpDetDivRec > opDet;
     for(size_t detnum = 0; detnum < priv_DivRecs.size(); ++detnum){
       if(priv_DivRecs.at(detnum)->OpDetNum() == opDetNum)
         opDet = priv_DivRecs.at(detnum);
     }
     if(!opDet)
     {
       throw cet::exception("CalibTreeBuilder") << "No sim:: OpDetDivRec corresponding "
         << "to opDetNum: " << opDetNum << "\n";
     }
     return opDet;
   }

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.

   {
     const double fDelay = PBS->GetDelay();
     std::vector< sim::SDP > retVec;
     double fPeakTime = opHit_P->PeakTime();
     double fWidth = opHit_P->Width();
     UInt_t fChan = opHit_P->OpChannel();
     art::ServiceHandle<geo::Geometry> geom;
     int fDet  = geom->OpDetFromOpChannel(fChan);
     //I should use the timing service for these time conversions.
     sim::OpDetBacktrackerRecord::timePDclock_t start_time = ((fPeakTime- fWidth)*1000.0)-fDelay;
     sim::OpDetBacktrackerRecord::timePDclock_t end_time = ((fPeakTime+ fWidth)*1000.0)-fDelay;
     if(start_time > end_time){throw;}//This is bad. Give a reasonable error message here, and use cet::except
 
     //BUG!!!fGeom->OpDetFromOpChannel(channel)
     art::Ptr<sim::OpDetBacktrackerRecord> fBTR;
     try{
       fBTR = PBS->FindOpDetBTR(fDet);
     }catch(cet::exception &e){
       if(true){
         std::vector<sim::SDP> ret;
         return ret;
       }else{
         throw e;
       }
     }
     const std::vector<std::pair<double, std::vector<sim::SDP>> >& timeSDPMap
       = fBTR->timePDclockSDPsMap(); //Not guranteed to be sorted.
     art::Ptr<sim::OpDetDivRec> div_rec = this->FindDivRec(fDet);//This is an OpDetDivRec collected from this BTR.
 
 
     std::vector<const std::pair<double, std::vector<sim::SDP>>*> timePDclockSDPMap_SortedPointers;
     std::vector<const sim::OpDet_Time_Chans*> div_rec_SortedPointers;
     for ( auto& pair : timeSDPMap ){ timePDclockSDPMap_SortedPointers.push_back(&pair); }
     auto pairSort = [](auto& a, auto& b) { return a->first < b->first ; } ;
     if( !std::is_sorted( timePDclockSDPMap_SortedPointers.begin(), timePDclockSDPMap_SortedPointers.end(), pairSort)){
       std::sort(timePDclockSDPMap_SortedPointers.begin(), timePDclockSDPMap_SortedPointers.end(), pairSort);
     }
 
     //This section is a hack to make comparisons work right.
     std::vector<sim::SDP> dummyVec;
     std::pair<double, std::vector<sim::SDP>> start_timePair = std::make_pair(start_time, dummyVec);
     std::pair<double, std::vector<sim::SDP>> end_timePair = std::make_pair(end_time, dummyVec);
     auto start_timePair_P = &start_timePair;
     auto end_timePair_P = &end_timePair;
 
     //First interesting iterator.
     auto map_pdsdp_itr = std::lower_bound(timePDclockSDPMap_SortedPointers.begin(), timePDclockSDPMap_SortedPointers.end(), start_timePair_P, pairSort);
     //Last interesting iterator.
     auto map_pdsdp_last = std::upper_bound(map_pdsdp_itr, timePDclockSDPMap_SortedPointers.end(), end_timePair_P, pairSort);
 
     //retvec.push_back(map_pdsdp_first.second[0]);
     //This code screams fragile. Really. If you read this, feel free to clean up the methods.
     for(auto& sdp_time_pair =  map_pdsdp_itr; sdp_time_pair != map_pdsdp_last; ++sdp_time_pair){
       //cut div_rec by time
       auto time=(*sdp_time_pair)->first;
       for(auto& sdp : (*sdp_time_pair)->second){
         //      auto sdp = (*sdp_time_pair)->second; //This is a vector
         //auto slice = div_rec->GetSlice(start_time, end_time);
         auto time_divr_pair = div_rec->FindClosestTimeChan(time);
         auto time_divr=time_divr_pair.first;
         auto time_divr_found=time_divr_pair.second;
         sim::SDP tmp = sdp;
         if( time_divr_found && time_divr->time == time){ //This will break if time_divr is the end value of the vector
           tmp.energy = time_divr->GetFrac(fChan) * tmp.energy;
           // ((*map_divrec_itr)->DivChans.eScaleFrac(fChan))*tmp.energy;
         }else{
           //          std::cout<<"I am having trouble reconciling "<<time<<" and "<<time_divr->time<<"\n";
           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.
         }
         retVec.emplace_back(tmp);
 
       }
     }
 
     //while( map_pdsdp_itr != map_pdsdp_last /*&& map_divrec_itr != map_divrec_last*/){ //Stepping through.
     /*for(auto const& sdp :  (*map_pdsdp_itr)->second){
       sim::SDP tmp = sdp;
       tmp.energy = ((*map_divrec_itr)->DivChans.eScaleFrac(fChan))*tmp.energy;
     //              second[fChan].tick_photons_frac)*tmp.energy;
     retVec.emplace_back(tmp);
     }
     ++map_pdsdp_itr;
     ++map_divrec_itr;
     }*/
 
     return retVec;
   }

void CalibrationTreeBuilder::CalibrationTreeBuilder::PrepDivRec ( const art::Event & evt )

private

Definition at line 653 of file CalibrationTreeBuilder_module.cc.

 {
   if( 0 ){ return;} //Insert check for DivRecs here, or don't use validHandle below.
   auto const& divrecHandle = evt.getValidHandle <std::vector<sim::OpDetDivRec>>(fWavLabel);
   if(divrecHandle.failedToGet()){
     return;
   }
   art::fill_ptr_vector(priv_DivRecs, divrecHandle);
   auto compareDivReclambda = [](art::Ptr<sim::OpDetDivRec> a, art::Ptr<sim::OpDetDivRec> b) {return(a->OpDetNum() < b->OpDetNum());};
   if (!std::is_sorted(priv_DivRecs.begin(), priv_DivRecs.end(), compareDivReclambda))
     std::sort(priv_DivRecs.begin(), priv_DivRecs.end(), compareDivReclambda);
 
 }

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:

duneana/duneana/DAQSimAna/CalibrationTreeBuilder.h
duneana/duneana/DAQSimAna/CalibrationTreeBuilder_module.cc

Public Member Functions

Private Member Functions

Private Attributes

Additional Inherited Members

Detailed Description

Constructor & Destructor Documentation

Member Function Documentation

Member Data Documentation