Inheritance diagram for sbnd::CRTTrackProducer:

Public Member Functions
	CRTTrackProducer (fhicl::ParameterSet const &p)

	CRTTrackProducer (CRTTrackProducer const &)=delete

	CRTTrackProducer (CRTTrackProducer &&)=delete

CRTTrackProducer &	operator= (CRTTrackProducer const &)=delete

CRTTrackProducer &	operator= (CRTTrackProducer &&)=delete

void	produce (art::Event &e) override

void	beginJob () override

void	endJob () override

Private Attributes
std::string	fDataLabelHits

std::string	fDataLabelTZeros

int	fTrackMethodType

int	fStoreTrack

CRTTrackRecoAlg	trackAlg

Detailed Description

Definition at line 47 of file sbndcode/sbndcode/CRT/CRTTrackProducer_module.cc.

Constructor & Destructor Documentation

sbnd::CRTTrackProducer::CRTTrackProducer ( fhicl::ParameterSet const & p )

explicit

Definition at line 119 of file sbndcode/sbndcode/CRT/CRTTrackProducer_module.cc.

   : EDProducer(p), trackAlg(p.get<fhicl::ParameterSet>("TrackAlg"))
 {  
 
   // Initialize member data here.
   fDataLabelHits      = p.get<std::string>("DataLabelHits");      // CRTHit producer module name
   fDataLabelTZeros    = p.get<std::string>("DataLabelTZeros");    // CRTTzero producer module name
   fStoreTrack         = p.get<int>        ("StoreTrack");         // method 1 = all, method 2 = ave, method 3 = pure, method 4 = top plane
   fTrackMethodType    = p.get<int>        ("TrackMethodType");    // Print stuff
   
   // Call appropriate produces<>() functions here.
   if(fStoreTrack == 1){ 
     produces< std::vector<sbn::crt::CRTTrack>   >();
     produces< art::Assns<sbn::crt::CRTTrack , sbn::crt::CRTHit> >();
   } 
 
 } // CRTTrackProducer()

sbnd::CRTTrackProducer::CRTTrackProducer ( CRTTrackProducer const & )

delete

sbnd::CRTTrackProducer::CRTTrackProducer ( CRTTrackProducer && )

delete

Member Function Documentation

void sbnd::CRTTrackProducer::beginJob ( )

override

Definition at line 361 of file sbndcode/sbndcode/CRT/CRTTrackProducer_module.cc.

 {
 
 } // CRTTrackProducer::beginJob()

void sbnd::CRTTrackProducer::endJob ( )

override

Definition at line 367 of file sbndcode/sbndcode/CRT/CRTTrackProducer_module.cc.

 {
 
 } // CRTTrackProducer::endJob()

CRTTrackProducer& sbnd::CRTTrackProducer::operator= ( CRTTrackProducer const & )

delete

CRTTrackProducer& sbnd::CRTTrackProducer::operator= ( CRTTrackProducer && )

delete

void sbnd::CRTTrackProducer::produce ( art::Event & e )

override

Definition at line 138 of file sbndcode/sbndcode/CRT/CRTTrackProducer_module.cc.

 {
 
   // For validation
   int nTrack = 0;
   int nCompTrack = 0;
   int nIncTrack = 0;
 
   // CRTTrack collection on this event                                                                         
   std::unique_ptr<std::vector<sbn::crt::CRTTrack> > CRTTrackCol(new std::vector<sbn::crt::CRTTrack>);
   std::unique_ptr< art::Assns<sbn::crt::CRTTrack, sbn::crt::CRTHit> > Trackassn( new art::Assns<sbn::crt::CRTTrack, sbn::crt::CRTHit>);
   art::PtrMaker<sbn::crt::CRTTrack> makeTrackPtr(evt);
 
   // Implementation of required member function here.
   art::Handle< std::vector<sbn::crt::CRTHit> > rawHandle;
   evt.getByLabel(fDataLabelHits, rawHandle); //what is the product instance name? no BernZMQ
 
   // Check to make sure the data we asked for is valid                                                                         
   if(!rawHandle.isValid()){
     mf::LogWarning("CRTTrackProducer")
       <<"No CRTHits from producer module "<<fDataLabelHits;
     return;
   }
 
   // Track method 4 = SBND method with top plane (doesn't use CRTTzero)
   if(fTrackMethodType == 4){
 
     //Get the CRT hits from the event
     std::vector<art::Ptr<sbn::crt::CRTHit> > hitlist;
     if (evt.getByLabel(fDataLabelHits, rawHandle))
       art::fill_ptr_vector(hitlist, rawHandle);
 
     std::map<art::Ptr<sbn::crt::CRTHit>, int> hitIds;
     for(size_t i = 0; i<hitlist.size(); i++){
       hitIds[hitlist[i]] = i;
     }
 
     std::vector<std::vector<art::Ptr<sbn::crt::CRTHit>>> CRTTzeroVect = trackAlg.CreateCRTTzeros(hitlist);
 
     // Loop over tzeros
     for(size_t i = 0; i<CRTTzeroVect.size(); i++){
 
       //loop over hits for this tzero, sort by tagger
       std::map<std::string, std::vector<art::Ptr<sbn::crt::CRTHit>>> hits;
       for (size_t ah = 0; ah< CRTTzeroVect[i].size(); ++ah){        
         std::string ip = CRTTzeroVect[i][ah]->tagger;       
         hits[ip].push_back(CRTTzeroVect[i][ah]);
       } // loop over hits
       
       //loop over planes and calculate average hits
       std::vector<std::pair<sbn::crt::CRTHit, std::vector<int>>> allHits;
       for (auto &keyVal : hits){
         std::string ip = keyVal.first;
         std::vector<std::pair<sbn::crt::CRTHit, std::vector<int>>> ahits = trackAlg.AverageHits(hits[ip], hitIds);
         allHits.insert(allHits.end(), ahits.begin(), ahits.end());
       }
 
       //Create tracks with hits at the same tzero
       std::vector<std::pair<sbn::crt::CRTTrack, std::vector<int>>> trackCandidates = trackAlg.CreateTracks(allHits);
       nTrack += trackCandidates.size();
       for(size_t j = 0; j < trackCandidates.size(); j++){
         CRTTrackCol->emplace_back(trackCandidates[j].first);
 
         art::Ptr<sbn::crt::CRTTrack> trackPtr = makeTrackPtr(CRTTrackCol->size()-1);
         for (size_t ah = 0; ah< CRTTzeroVect[i].size(); ++ah){        
           Trackassn->addSingle(trackPtr, CRTTzeroVect[i][ah]);
         }
         if(trackCandidates[j].first.complete) nCompTrack++;
         else nIncTrack++;
       }
     }
 
   //Older track reconstruction methods from MicroBooNE
   }else{
     //Get list of tzeros             
     //  std::vector<sbn::crt::CRTHit> const& CRTHitCollection(*rawHandle);
     art::Handle< std::vector<sbn::crt::CRTTzero> > rawHandletzero;
     evt.getByLabel(fDataLabelTZeros, rawHandletzero); //what is the product instance name? no BernZMQ
     
     //check to make sure the data we asked for is valid                                                                            
     if(!rawHandletzero.isValid()){
       mf::LogWarning("CRTTrackProducer")
         <<"No CRTTzeros from producer module "<<fDataLabelTZeros;
       return;
     }
     
     std::vector<art::Ptr<sbn::crt::CRTTzero> > tzerolist;
     if (evt.getByLabel(fDataLabelTZeros,rawHandletzero))
       art::fill_ptr_vector(tzerolist, rawHandletzero);
  
     art::FindManyP<sbn::crt::CRTHit> fmht(rawHandletzero, evt, fDataLabelTZeros);
  
     //loop over tzeros
     for(size_t tzIter = 0; tzIter < tzerolist.size(); ++tzIter){   
       
       //count planes with hits for this tzero
       int np =0 ; //int ipflag[7] = {}; // CHANGED FROM 4 TO 7
       int tothits =0;
       for (int ip=0;ip<7;++ip) { // CHANGED FROM 4 TO 7
         if (tzerolist[tzIter]->nhits[ip]>0)  { np++; /*ipflag[ip]=1;*/ tothits+=tzerolist[tzIter]->nhits[ip];}  
       }
  
       if (np<2) continue;
       std::vector<art::Ptr<sbn::crt::CRTHit> > hitlist=fmht.at(tzIter);
       //for(size_t hit_i = 0; hit_i < hitlist.size(); hit
       if (fTrackMethodType==1) {
         double time_s_A = hitlist[0]->ts0_s;
  
         // find pairs of hits in different planes
         for (size_t ah = 0; ah< hitlist.size()-1; ++ah){        
           sbn::crt::CRTHit temphit=*hitlist[ah];
           CRTavehit Ahit = copyme(temphit);
           int planeA = hitlist[ah]->plane;
  
           for (size_t bh = ah+1; bh< hitlist.size(); ++bh){        
             int planeB = hitlist[bh]->plane;
  
             if (planeB!=planeA && !((planeA==3&&planeB==4)||(planeA==4&&planeB==3))) {  // make a track               
               temphit=*hitlist[bh];
               CRTavehit Bhit = copyme(temphit);
               sbn::crt::CRTTrack CRTcanTrack=shcut(Ahit,Bhit,time_s_A,0);
               CRTTrackCol->emplace_back(CRTcanTrack);
             }
  
           }
  
         }
  
       }
  
       else if ((fTrackMethodType==2) || (fTrackMethodType==3 && np==2 && tothits==2)) {        
  
         //loop over hits and get average x,y,z,pe for each plane CHANGED FROM 4 TO 7
         std::vector<float> thittime0[7];
         std::vector<float> thittime1[7];
         std::vector<float> tx[7];
         std::vector<float> ty[7];
         std::vector<float> tz[7];
         std::vector<float> pe[7];
         
         double time_s_A = hitlist[0]->ts0_s;
         //      double time_s_err = hitlist[0]->ts0_s_err;
         double time_s_err = 0.;
         double time1_ns_A = hitlist[0]->ts1_ns;
         double time0_ns_A = hitlist[0]->ts0_ns;
        
         //loop over hits for this tzero, sort by plane
         for (size_t ah = 0; ah< hitlist.size(); ++ah){        
           int ip = hitlist[ah]->plane;       
           thittime0[ip].push_back(hitlist[ah]->ts0_ns-time0_ns_A);
           thittime1[ip].push_back(hitlist[ah]->ts1_ns-time1_ns_A);
           tx[ip].push_back(hitlist[ah]->x_pos);
           ty[ip].push_back(hitlist[ah]->y_pos);
           tz[ip].push_back(hitlist[ah]->z_pos);
           pe[ip].push_back(hitlist[ah]->peshit);        
         } // loop over hits
  
         CRTavehit aveHits[7];
         //loop over planes and calculate average hits
         for (int ip = 0; ip < 7; ip++){
           if (tx[ip].size()>0){
             uint32_t at0; int32_t at1; uint16_t rt0,rt1;
             float totpe=0.0;
             float avet1=0.0; float rmst1 =0.0; 
             float avet0=0.0; float rmst0 =0.0; 
             float avex=0.0; float rmsx =0.0; 
             float avey=0.0; float rmsy =0.0; 
             float avez=0.0; float rmsz =0.0;
             vmanip(thittime0[ip],&avet0,&rmst0);
             vmanip(thittime1[ip],&avet1,&rmst1);
             at0 = (uint32_t)(avet0+time0_ns_A); rt0 = (uint16_t)rmst0;   
             at1 = (int32_t)(avet1+time1_ns_A); rt1 = (uint16_t)rmst1;
             vmanip(tx[ip],&avex,&rmsx);
             vmanip(ty[ip],&avey,&rmsy);
             vmanip(tz[ip],&avez,&rmsz);
             totpe=std::accumulate(pe[ip].begin(), pe[ip].end(), 0.0);
             CRTavehit aveHit = fillme(at0,rt0,at1,rt1,avex,rmsx,avey,rmsy,avez,rmsz,totpe,ip,"");
             aveHits[ip] = aveHit;
           }
           else {
             CRTavehit aveHit = fillme(0,0,0,0,-99999,-99999,-99999,-99999,-99999,-99999,-99999,ip,"");
             aveHits[ip] = aveHit;
           }
         }  
  
         // find pairs of hits in different planes
         for (size_t ah = 0; ah< 6; ++ah){        
           CRTavehit Ahit = aveHits[ah];
           if( Ahit.x_pos==-99999 ) continue;
  
           for (size_t bh = ah+1; bh< 7; ++bh){        
             if ( aveHits[bh].x_pos==-99999 ) continue;
  
             if (ah!=bh && !(ah==3&&bh==4)) {  // make a track               
               CRTavehit Bhit = aveHits[bh];
               sbn::crt::CRTTrack CRTcanTrack=shcut(Ahit,Bhit,time_s_A,time_s_err);
               CRTTrackCol->emplace_back(CRTcanTrack);
               nTrack++;
             }
  
           }
  
         }
  
       }
       
     }// loop over tzeros
   }
   
   mf::LogInfo("CRTTrackProducer")
     <<"Number of tracks            = "<<CRTTrackCol->size()<<"\n"
     <<"Number of complete tracks   = "<<CRTTrackCol->size()-nIncTrack<<"\n"
     <<"Number of incomplete tracks = "<<nIncTrack;
 
   //store track collection into event
   if(fStoreTrack == 1){
     evt.put(std::move(CRTTrackCol));
     evt.put(std::move(Trackassn));
   }
   
 } // CRTTrackProducer::produce()