Inheritance diagram for gar::rec::dayoneconverter:

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

	dayoneconverter (dayoneconverter const &)=delete

	dayoneconverter (dayoneconverter &&)=delete

dayoneconverter &	operator= (dayoneconverter const &)=delete

dayoneconverter &	operator= (dayoneconverter &&)=delete

void	produce (art::Event &e) override

Public Member Functions inherited from art::EDProducer
	EDProducer (fhicl::ParameterSet const &pset)

template<typename Config >
	EDProducer (Table< Config > const &config)

std::string	workerType () const

Public Member Functions inherited from art::detail::Producer
virtual	~Producer () noexcept

	Producer (fhicl::ParameterSet const &)

	Producer (Producer const &)=delete

	Producer (Producer &&)=delete

Producer &	operator= (Producer const &)=delete

Producer &	operator= (Producer &&)=delete

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::Modifier
	~Modifier () noexcept

	Modifier ()

	Modifier (Modifier const &)=delete

	Modifier (Modifier &&)=delete

Modifier &	operator= (Modifier const &)=delete

Modifier &	operator= (Modifier &&)=delete

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
int	makepatrectrack (std::vector< gar::rec::TPCCluster > &trackTPCClusters, gar::rec::TrackPar &trackpar)

void	digitizeCaloHitsSimple (gar::sdp::CaloDeposit cd, float *fcpos, float &energy, float &time)

void	digitizeCaloHitsMu2e (gar::sdp::CaloDeposit cd, float *fcpos, float &energy, float &time)

Private Attributes
std::string	fInputEdepLabel
	Input label for edeps. More...

std::string	fInputEdepInstanceTPC
	Input instance for TPC edeps. More...

std::string	fInputEdepInstanceMuID
	Input instance for MuID edeps. More...

bool	fIncludeMuIDhits
	Include MuID hits as TPCClusters. More...

float	fSmearX
	amount by which to smear X, in cm More...

float	fSmearY
	amount by which to smear Y, in cm More...

float	fSmearT
	amount by which to smear T, in ns More...

float	fPECm
	conversion factor from cm step length to pe More...

float	fSmearLY
	amount by which to smear the LY More...

float	fThrPE
	threshold cut in pe More...

float	fZCut1
	Cut to ensure TPC Clusters are on different planes. More...

float	fZCut2
	Cut to ensure TPC clusters are on the same plane. More...

float	fRCut
	Road in the YZ plane to add hits on a circle. More...

float	fTimeBunch
	Time bunch spread, in ns. More...

bool	fUseOnlyTrueMuonHits
	whether to cheat and use true muon hits More...

std::string	fG4ModuleLabel
	Use this to get the MCParticles. More...

CLHEP::HepRandomEngine &	fEngine

const gar::geo::GeometryCore *	fGeo
	geometry information More...

gar::geo::BitFieldCoder *	fFieldDecoderTrk

gar::geo::BitFieldCoder *	fFieldDecoderMuID

TF1 *	fMu2e

Additional Inherited Members
Public Types inherited from art::EDProducer
using	ModuleType = EDProducer

using	WorkerType = WorkerT< EDProducer >

Public Types inherited from art::detail::Producer
template<typename UserConfig , typename KeysToIgnore = void>
using	Table = Modifier::Table< UserConfig, KeysToIgnore >

Public Types inherited from art::Modifier
template<typename UserConfig , typename UserKeysToIgnore = void>
using	Table = ProducerTable< UserConfig, detail::ModuleConfig, UserKeysToIgnore >

Static Public Member Functions inherited from art::EDProducer
static void	commitEvent (EventPrincipal &ep, Event &e)

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 54 of file dayoneconverter_module.cc.

Constructor & Destructor Documentation

gar::rec::dayoneconverter::dayoneconverter ( fhicl::ParameterSet const & p )

explicit

Definition at line 103 of file dayoneconverter_module.cc.

       : EDProducer{p}
       , fEngine(art::ServiceHandle<rndm::NuRandomService>()->createEngine(*this,p,"Seed"))
         // More initializers here.
       {
         // Call appropriate produces<>() functions here.
         // Call appropriate consumes<>() for any products to be retrieved by this module.
 
         fInputEdepLabel = p.get<std::string>("InputLabel","edepconvert");
         fInputEdepInstanceTPC = p.get<std::string>("InputInstanceTPC","TrackerSc");
         fInputEdepInstanceMuID = p.get<std::string>("InputInstanceMuID","MuID");
         fIncludeMuIDhits = p.get<bool>("IncludeMuIDhits", false);
         fSmearX = p.get<float>("SmearX",0.3);  // in cm
         fSmearY = p.get<float>("SmearY",0.3);  // in cm
         fSmearT = p.get<float>("SmearT",1.0);  // in ns
         fPECm  = p.get<float>("PeCm", 50.);  // in pe/cm
         fSmearLY  = p.get<float>("SmearLY", 10.0);  // in pe
         fThrPE    = p.get<float>("ThrPE", 5.0);  // in pe
         fZCut1    = p.get<float>("ZCut1",10.0); // in cm
         fZCut2    = p.get<float>("ZCut2",0.5); // in cm
         fRCut     = p.get<float>("RCut" ,5.0);  // in cm
         fTimeBunch = p.get<float>("TimeBunch",150000); // in ns            
         fUseOnlyTrueMuonHits = p.get<bool>("UseOnlyTrueMuonHits",false);  // whether to cheat or not
         fG4ModuleLabel = p.get<std::string>("G4ModuleLabel","geant");  // for retrieving MCParticles
 
         art::InputTag tpcedeptag(fInputEdepLabel,fInputEdepInstanceTPC);
         art::InputTag muidedeptag(fInputEdepLabel,fInputEdepInstanceMuID);
         consumes<std::vector<gar::sdp::CaloDeposit> >(tpcedeptag);
         consumes<std::vector<gar::sdp::CaloDeposit> >(muidedeptag);
         if (fUseOnlyTrueMuonHits)
           {
             consumes<std::vector<simb::MCParticle> >(fG4ModuleLabel);
           }
         produces<std::vector<gar::rec::Track> >();
         produces<std::vector<gar::rec::TPCCluster> >();
         produces< art::Assns<gar::rec::TPCCluster, gar::rec::Track> >();
 
         fGeo = gar::providerFrom<gar::geo::GeometryGAr>();
         std::string fEncoding = fGeo->GetMinervaCellIDEncoding();
         fFieldDecoderTrk = new gar::geo::BitFieldCoder( fEncoding );
         std::string fEncodingMuID = fGeo->GetMuIDCellIDEncoding();
         fFieldDecoderMuID = new gar::geo::BitFieldCoder( fEncodingMuID );
 
         fMu2e = new TF1("mu2e_pe", "[0] + [1]*x", 0, 600);
         fMu2e->SetParameter(0, 50.);//pe
         fMu2e->SetParameter(1, -0.06);//pe/cm
       }

gar::rec::dayoneconverter::dayoneconverter ( dayoneconverter const & )

delete

gar::rec::dayoneconverter::dayoneconverter ( dayoneconverter && )

delete

Member Function Documentation

void gar::rec::dayoneconverter::digitizeCaloHitsMu2e	(	gar::sdp::CaloDeposit	cd,
		float *	fcpos,
		float &	energy,
		float &	time
	)

private

Definition at line 580 of file dayoneconverter_module.cc.

     {
       //Need to check in which direction to smear (depends on the segmenetation)
       //Can use the cellID decoder to know looking at cellX and cellY values
       CLHEP::RandGauss GaussRand(fEngine);
 
       gar::raw::CellID_t cID = cd.CellID();
       int cellX = fFieldDecoderTrk->get(cID, "cellX");
       int cellY = fFieldDecoderTrk->get(cID, "cellY");
 
       fcpos[0] = cd.X(); // default values
       fcpos[1] = cd.Y();
       fcpos[2] = cd.Z();
       energy = cd.Energy();
       time = cd.Time();
 
       if(cellX == 0) {
         //Segmented in Y
         fcpos[0] += GaussRand.fire(0., fSmearX);
       }
 
       if(cellY == 0) {
         //Segmented in X
         fcpos[1] += GaussRand.fire(0., fSmearY);
       }
 
       time += GaussRand.fire(0., fSmearT);
 
       //convert energy to pe based on the distance in the strip
       //Curve as "linear" with parameters a = 0.06 pe/cm, b = 50 pe (80*0.6 from Mu2e)
       double local_distance = 0.;
       std::array<double, 3> point = {cd.X(), cd.Y(), cd.Z()};
       std::array<double, 3> pointLocal;
       gar::geo::LocalTransformation<TGeoHMatrix> trans;
       fGeo->WorldToLocal(point, pointLocal, trans);
       TVector3 tpoint(point[0], point[1], point[2]);
       std::string name = fGeo->VolumeName(tpoint);
 
       if(cellX == 0) {
         //Segmented in Y -> get the x pos
         local_distance = 300. - pointLocal[0];
       }
 
       if(cellY == 0) {
         //Segmented in X -> get the y pos
         local_distance = 300. - pointLocal[1];
       }
 
       energy = fMu2e->Eval(local_distance);
       // std::cout << "Volume " << name << std::endl;
       // std::cout << "CellX " << cellX << " CellY " << cellY << std::endl;
       // std::cout << "Global Position " << cd.X() << ", " << cd.Y() << ", " << cd.Z() << std::endl;
       // std::cout << "Local Position " << pointLocal[0] << ", " << pointLocal[1] << ", " << pointLocal[2] << std::endl;
       // std::cout << "local distance in strip " << local_distance << " cm has " << energy << " pe detected" << std::endl;
       energy += GaussRand.fire(0., fSmearLY);
 
       if(energy < fThrPE) energy = 0.;
 
       return;
     }

void gar::rec::dayoneconverter::digitizeCaloHitsSimple	(	gar::sdp::CaloDeposit	cd,
		float *	fcpos,
		float &	energy,
		float &	time
	)

private

Definition at line 541 of file dayoneconverter_module.cc.

     {
       //Need to check in which direction to smear (depends on the segmenetation)
       //Can use the cellID decoder to know looking at cellX and cellY values
       CLHEP::RandGauss GaussRand(fEngine);
 
       gar::raw::CellID_t cID = cd.CellID();
       int cellX = fFieldDecoderTrk->get(cID, "cellX");
       int cellY = fFieldDecoderTrk->get(cID, "cellY");
 
       fcpos[0] = cd.X(); // default values
       fcpos[1] = cd.Y();
       fcpos[2] = cd.Z();
       energy = cd.Energy();
       time = cd.Time();
 
       if(cellX == 0) {
         //Segmented in Y
         fcpos[0] += GaussRand.fire(0., fSmearX);
       }
 
       if(cellY == 0) {
         //Segmented in X
         fcpos[1] += GaussRand.fire(0., fSmearY);
       }
 
       time += GaussRand.fire(0., fSmearT);
 
       //convert energy to pe
       energy = fPECm * cd.StepLength();
       energy += GaussRand.fire(0., fSmearLY);
 
       if(energy < fThrPE) energy = 0.;
 
       return;
     }

int gar::rec::dayoneconverter::makepatrectrack	(	std::vector< gar::rec::TPCCluster > &	trackTPCClusters,
		gar::rec::TrackPar &	trackpar
	)

private

Definition at line 645 of file dayoneconverter_module.cc.

     {
       // track parameters:  x is the independent variable
       // 0: y
       // 1: z
       // 2: curvature
       // 3: phi
       // 4: lambda = angle from the cathode plane
       // 5: x   /// added on to the end
 
       float fSortTransWeight = 1.0;
       int fInitialTPNTPCClusters = 100;
       float fSortDistBack = 2.0;
       int fPrintLevel = 0;
 
       float lengthforwards = 0;
       std::vector<int> hlf;
       float lengthbackwards = 0;
       std::vector<int> hlb;
 
       gar::rec::sort_TPCClusters_along_track(trackTPCClusters,hlf,hlb,fPrintLevel,lengthforwards,lengthbackwards,fSortTransWeight,fSortDistBack);
 
       std::vector<float> tparbeg(6,0);
       float xother = 0;
       if ( gar::rec::initial_trackpar_estimate(trackTPCClusters, hlf, tparbeg[2], tparbeg[4],
                                                tparbeg[3], tparbeg[5], tparbeg[0], tparbeg[1], xother, fInitialTPNTPCClusters, fPrintLevel) != 0)
         {
           return 1;
         }
       float chisqforwards = 0;
       float tvpos[3] = {tparbeg[5],tparbeg[0],tparbeg[1]};
       for (size_t i=0; i<trackTPCClusters.size(); ++i)
         {
           float dist = 0;
           const float *pos = trackTPCClusters.at(i).Position();
           int retcode = util::TrackPropagator::DistXYZ(tparbeg.data(),tvpos,pos,dist);
           if (retcode == 0)
             {
               chisqforwards += dist*dist/(fSmearX*fSmearX);   // crude parameterization of errors
             }
         }
 
       std::vector<float> tparend(6,0);
       if ( gar::rec::initial_trackpar_estimate(trackTPCClusters, hlb, tparend[2], tparend[4],
                                                tparend[3], tparend[5], tparend[0], tparend[1], xother, fInitialTPNTPCClusters, fPrintLevel) != 0)
         {
           return 1;
         }
       float chisqbackwards = 0;
       float tepos[3] = {tparend[5],tparend[0],tparend[1]};
       for (size_t i=0; i<trackTPCClusters.size(); ++i)
         {
           float dist = 0;
           const float *pos = trackTPCClusters.at(i).Position();
           int retcode = util::TrackPropagator::DistXYZ(tparend.data(),tepos,pos,dist);
           if (retcode == 0)
             {
               chisqbackwards += dist*dist/(fSmearX*fSmearX);   // crude parameterization of errors
             }
         }
 
 
       // no covariance matrix in patrec tracks
       float covmatbeg[25];
       float covmatend[25];
       for (size_t i=0; i<25; ++i) // no covmat in patrec tracks
         {
           covmatend[i] = 0;
           covmatbeg[i] = 0;
         }
 
       trackpar.setNTPCClusters(trackTPCClusters.size());
       trackpar.setTime(0);
       trackpar.setChisqForwards(chisqforwards);
       trackpar.setChisqBackwards(chisqbackwards);
       trackpar.setLengthForwards(lengthforwards);
       trackpar.setLengthBackwards(lengthbackwards);
       trackpar.setCovMatBeg(covmatbeg);
       trackpar.setCovMatEnd(covmatend);
       trackpar.setTrackParametersBegin(tparbeg.data());
       trackpar.setXBeg(tparbeg[5]);
       trackpar.setTrackParametersEnd(tparend.data());
       trackpar.setXEnd(tparend[5]);
 
       return 0;
     }

dayoneconverter& gar::rec::dayoneconverter::operator= ( dayoneconverter const & )

delete

dayoneconverter& gar::rec::dayoneconverter::operator= ( dayoneconverter && )

delete

void gar::rec::dayoneconverter::produce ( art::Event & e )

overridevirtual

Implements art::EDProducer.

Definition at line 151 of file dayoneconverter_module.cc.

     {
       // Implementation of required member function here.
       std::unique_ptr< std::vector<gar::rec::TPCCluster> > TPCClusterCol(new std::vector<gar::rec::TPCCluster>);
       std::unique_ptr< std::vector<gar::rec::Track> > trkCol(new std::vector<gar::rec::Track>);
       std::unique_ptr< art::Assns<gar::rec::TPCCluster,gar::rec::Track> > TPCClusterTrkAssns(new ::art::Assns<gar::rec::TPCCluster,gar::rec::Track>);
 
       art::InputTag tpcedeptag(fInputEdepLabel,fInputEdepInstanceTPC);
       auto tccdHandle = e.getValidHandle< std::vector<gar::sdp::CaloDeposit> >(tpcedeptag);
       auto const& tccds = *tccdHandle;
 
       // put these back when we have tracks and can associate them
       auto const trackPtrMaker = art::PtrMaker<gar::rec::Track>(e);
       auto const tpcclusPtrMaker = art::PtrMaker<gar::rec::TPCCluster>(e);
 
       auto const *magFieldService = gar::providerFrom<mag::MagneticFieldService>();
       G4ThreeVector zerovec(0,0,0);
       G4ThreeVector magfield = magFieldService->FieldAtPoint(zerovec);
 
       art::Handle<std::vector<simb::MCParticle > > mcphandle;
       std::unordered_map<int, const simb::MCParticle * > TrkIDMap;
       if (fUseOnlyTrueMuonHits)
         {
           mcphandle = e.getHandle<std::vector<simb::MCParticle> >(fG4ModuleLabel);
           if (!mcphandle)
             {
                throw cet::exception("DayOneConverter::Produce")
                     << "Attempting to identify muon hits without MC particles";
             }
           for (auto const &mcp : *mcphandle)
             {
               TrkIDMap[mcp.TrackId()] = &mcp;
             }
         }
 
       // first stab, just make all the TPC clusters in the event, and then worry later
       // about pattern recognition
 
       for (auto const& cd : tccds)
         {
           const int trackid = cd.TrackID();
           if (fUseOnlyTrueMuonHits)
             {
               auto mcpt = TrkIDMap.find(trackid);
               if (mcpt == TrkIDMap.end()) continue;
               if (std::abs(mcpt->second->PdgCode()) != 13) continue;
             }
 
           float energy = 0.;
           float time = 0.;
           float fcpos[3] = {0., 0., 0.};
           
           // digitizeCaloHitsSimple(cd, fcpos, energy, time);
           digitizeCaloHitsMu2e(cd, fcpos, energy, time);
 
           float covmat[6] = {0,0,0,0,0,0};  // TODO -- fill this in with something reasonble
 
           if(energy <= 0.) continue;
 
           TPCClusterCol->emplace_back(energy,
                                       fcpos,
                                       time,     // time is in ns
                                       time,
                                       time,
                                       fSmearX,
                                       covmat);
         }
 
       if(fIncludeMuIDhits)
         {
           art::InputTag muidedeptag(fInputEdepLabel,fInputEdepInstanceMuID);
           auto muIDHandle = e.getValidHandle< std::vector<gar::sdp::CaloDeposit> >(muidedeptag);
           auto const& muids = *muIDHandle;
 
           for (auto const& cd : muids)
             {
 
               const int trackid = cd.TrackID();
               if (fUseOnlyTrueMuonHits)
                 {
                   auto mcpt = TrkIDMap.find(trackid);
                   if (mcpt == TrkIDMap.end()) continue;
                   if (std::abs(mcpt->second->PdgCode()) != 13) continue;
                 }
 
               float energy = 0.;
               float time = 0.;
               float fcpos[3] = {0., 0., 0.};
 
               // digitizeCaloHitsSimple(cd, fcpos, energy, time);
               digitizeCaloHitsMu2e(cd, fcpos, energy, time);
 
               if(energy <= 0.) continue;
 
               float covmat[6] = {0,0,0,0,0,0};  // TODO -- fill this in with something reasonble
 
               TPCClusterCol->emplace_back(energy,
                                           fcpos,
                                           time,     // time is in ns
                                           time,
                                           time,
                                           fSmearX,
                                           covmat);
             }
         }
 
       if (true)
         {
 
           // sort the TPC Clusters based on time
           bool debug=false;
 
           std::vector<std::vector<gar::rec::TPCCluster>> TPCClusterColTime;
           std::vector<gar::rec::TPCCluster> TPCClusterColTimeBlock;
           std::sort(TPCClusterCol->begin(),TPCClusterCol->end(),
                     [](const gar::rec::TPCCluster &a, const gar::rec::TPCCluster &b)->bool
                     { return a.Time() < b.Time(); } );
           size_t ntpcclus = TPCClusterCol->size();
 
           if (ntpcclus!=0&&debug) std::cout<<"Time start and end: "<<TPCClusterCol->at(0).Time()<<" "<<TPCClusterCol->at(TPCClusterCol->size()-1).Time()<<std::endl;
 
           float startt= 0;
           if (ntpcclus!=0)
             {
               startt=TPCClusterCol->at(0).Time();
               for(size_t i=0; i<ntpcclus; i++)
                 {
                         
                   if(TPCClusterCol->at(i).Time()-startt<=fTimeBunch)
                     {
                       TPCClusterColTimeBlock.push_back(TPCClusterCol->at(i));
                     }
                   else
                     {
                       TPCClusterColTime.push_back(TPCClusterColTimeBlock);
                       startt=TPCClusterCol->at(i).Time();
                       TPCClusterColTimeBlock.clear();     
                       TPCClusterColTimeBlock.push_back(TPCClusterCol->at(i));             
                     }
                         
                 }
 
               if(TPCClusterColTimeBlock.size()!=0) TPCClusterColTime.push_back(TPCClusterColTimeBlock);
             }
 
 
           if(debug)
             {
               for(size_t c=0; c<TPCClusterColTime.size(); c++)
                 {
                     
                   std::cout << "l"<<c<<" = { ";
                   for (size_t d=0; d<TPCClusterColTime.at(c).size(); d++) {
                     std::cout << TPCClusterColTime.at(c).at(d).Time() << ", ";
                   }
                   std::cout << "};\n";
                 
                 }
             }
                 
           for(size_t t=0; t<TPCClusterColTime.size(); t++)
             {
 
               ////Sort the clusters along Z
               std::sort(TPCClusterColTime.at(t).begin(),TPCClusterColTime.at(t).end(),
                         [](const gar::rec::TPCCluster &a, const gar::rec::TPCCluster &b)->bool
                         { return a.Position()[2] < b.Position()[2]; } );
 
 
               size_t ntpcclust = TPCClusterColTime.at(t).size();
                     
                     
               // look for best-fit tracks in the list of TPC clusters
               // find the best triplet of TPC Clusters with spacing at least fZCut that fit on a circle
               // to think about -- time and energy cuts
 
                                   //create list of unused TPC clusters: a total one, and one for each plane
                     
                     std::list<size_t> unusedTPC;
                     std::list<size_t> TPCplane;
                     std::vector<std::list<size_t>> unusedTPCplanes;
                     float startz= 0;
                     float fZCut3= 4; //in cm thickness of the tracking plane
 
                     if (ntpcclust!=0)
                     {
                         startz=TPCClusterColTime.at(t).at(0).Position()[2];
                         for(size_t i=0; i<ntpcclust; i++)
                         {
                             unusedTPC.push_back(i);
                             
                             if(TPCClusterColTime.at(t).at(i).Position()[2]-startz<=fZCut3)
                             {
                                 TPCplane.push_back(i);
                             }
                             else
                             {
                                 unusedTPCplanes.push_back(TPCplane);
                                 startz=TPCClusterColTime.at(t).at(i).Position()[2];
                                 TPCplane.clear();     
                                 TPCplane.push_back(i);             
                             }
                             
                         }
 
                         if(TPCplane.size()!=0) unusedTPCplanes.push_back(TPCplane);
                     }
                     else
                     {
                     for(size_t i=0; i<ntpcclust; i++)
                         {
                             unusedTPC.push_back(i);
                         }  
                     }
 
                     if(debug){
                     std::cout << "Plane division: p = { ";
                     for (size_t n : unusedTPC) {
                     std::cout << TPCClusterColTime.at(t).at(n).Time() << ", ";
                     }
                     std::cout << "};\n";
 
                     for(size_t c=0; c<unusedTPCplanes.size(); c++)
                     {
                     std::cout << "p"<<c<<" = { ";
                         for (size_t n : unusedTPCplanes.at(c)) {
                         std::cout << TPCClusterColTime.at(t).at(n).Time() << ", ";
                         }
                     std::cout << "};\n";
                     }
                     }
                     
                     int done= 0;
                         
                     while (done==0)
                     {
                         size_t bestnpts = 0;
                         float bestsumr2 = 0;
                         std::vector<size_t> besttpcclusindex;
                     
                         for (size_t i=0; i<unusedTPCplanes.size(); ++i)
                         {
                         for (size_t it : unusedTPCplanes.at(i))
                         {
                             //const float *ipos = TPCClusterCol->at(it).Position();
                             for (size_t j=i+1; j<unusedTPCplanes.size(); ++j)
                             {
                             for (size_t jt : unusedTPCplanes.at(j))
                             {   
                                 //const float *jpos = TPCClusterCol->at(jt).Position();
                                 //if (TMath::Abs( ipos[2] - jpos[2] ) < fZCut1) continue;         Now superfluous using the lists
                                 for (size_t k=j+1; k<unusedTPCplanes.size(); ++k)
                                 {
                                 for (size_t kt : unusedTPCplanes.at(k))
                                 {
                                     //const float *kpos = TPCClusterCol->at(kt).Position();
                                     //if (TMath::Abs( ipos[2] - kpos[2] ) < fZCut1) continue;       Now superfluous using the lists
                                     std::vector<gar::rec::TPCCluster> triplet;
                                     triplet.push_back(TPCClusterCol->at(it));
                                     triplet.push_back(TPCClusterCol->at(jt));
                                     triplet.push_back(TPCClusterCol->at(kt));
                                     gar::rec::TrackPar triplettrack;
                                     makepatrectrack(triplet,triplettrack);
 
                                     // pick the best TPC clusters at each Z position.  Save their
                                     // indices in tpcclusindex and sum the squares of distances in sum
                                     std::vector<size_t> tpcclusindex;
                                     float sumr2 = 0;
 
                                     float zcur = -2E9;
                                     int tpcclusindexb = -1;
                                     float dbest=1E9;
                                     gar::rec::Track tpt;
                                     gar::rec::TPCCluster tpcclusb;
 
                                     for (size_t kt2 : unusedTPC)
                                     {
                                         const float *k2pos = TPCClusterCol->at(kt2).Position();
 
                                         // clusters are sorted along Z.  If we found a new Z, put the best point on the list
 
                                         if ((TMath::Abs(zcur - k2pos[2]) > fZCut2) &&
                                         (tpcclusindexb > -1) &&
                                         (dbest < fRCut) )
                                         {
                                             tpcclusindex.push_back(tpcclusindexb);
                                             //TPCtrial.push_back(tpcclusb);
                                             sumr2 += dbest*dbest;
                                             dbest = 1E9;
                                             tpcclusindexb = -1;
                                             zcur = k2pos[2];
                                         }
 
                                         float dist=0;
                                         tpt = triplettrack.CreateTrack();
                                         int retcode = util::TrackPropagator::DistXYZ(tpt.TrackParBeg(),tpt.Vertex(),k2pos,dist);
                                         if (retcode != 0) continue;
                                         if (dist > fRCut) continue;
                                         if (dist<dbest)
                                         {
                                             dbest = dist;
                                             tpcclusindexb = kt2;
                                         }
                                         // last point -- check to see if it gets added.
                                         if (kt2 == *std::prev(unusedTPC.end(),1) && tpcclusindexb > -1)
                                         {
                                             tpcclusindex.push_back(tpcclusindexb);
                                             sumr2 += dbest*dbest;
                                         }
                                     }  // end loop over kt2 -- assigning clusters to this track
 
                                     if (tpcclusindex.size() > bestnpts ||
                                     ((tpcclusindex.size() == bestnpts) &&
                                     (sumr2 < bestsumr2)))
                                     {
                                         bestnpts = tpcclusindex.size();
                                         bestsumr2 = sumr2;
                                         besttpcclusindex = tpcclusindex;
                                     }
                                 } // end loop over k in triplet
                                 }
                             } // end loop over j in triplet
                             }
                         } // end loop over i in triplet
                         }
 
                         for(size_t i=0;i<bestnpts;i++)
                         {
                         unusedTPC.remove(besttpcclusindex.at(i));
                         for(size_t j=0;j<unusedTPCplanes.size();j++)
                         {
                             unusedTPCplanes.at(j).remove(besttpcclusindex.at(i));
                         }
                         }
                         
                         if (debug)
                         {  
                             std::cout<<"bestnpts: "<<bestnpts<<std::endl;
                             for(size_t j=0;j<unusedTPCplanes.size();j++)
                                 {
                                     std::cout<<"unusedTPCplane"<<j<<" : "<<unusedTPCplanes.at(j).size()<<std::endl;
                                 }
                             std::cout<<"unusedTPC: "<<unusedTPC.size()<<std::endl<<std::endl;
                         }
 
                         if (bestnpts > 0)
                         {
                             // make track from all the TPC clusters
                             //trkCol->push_back(besttrack);
                             std::vector<gar::rec::TPCCluster> tcv;
                             for (size_t i=0;i<besttpcclusindex.size(); ++i)
                             {
                                 tcv.push_back(TPCClusterCol->at(besttpcclusindex.at(i)));
                             }
 
                             gar::rec::TrackPar btp;
                             makepatrectrack(tcv,btp);
                             gar::rec::Track btt = btp.CreateTrack();
                             trkCol->push_back(btt);
 
                             auto const trackpointer = trackPtrMaker(trkCol->size()-1);
                             for (size_t i=0; i<besttpcclusindex.size(); ++i)
                             {
                                 auto const tpccluspointer = tpcclusPtrMaker(besttpcclusindex.at(i));
                                 TPCClusterTrkAssns->addSingle(tpccluspointer,trackpointer);
                             }
                         }
                         else
                         {
                         done= 1;
                         if (debug) std::cout<<"Not able to find a new track, stop cycle, done="<<done<<std::endl;
                         }
 
                     
                     }
 
                 }
 
             if(debug) std::cout<<"Found this many tracks: "<<trkCol->size()<<std::endl<<std::endl<<std::endl;
         }  
         
             
             e.put(std::move(trkCol));
             e.put(std::move(TPCClusterCol));
             e.put(std::move(TPCClusterTrkAssns));
             
     }

Member Data Documentation

CLHEP::HepRandomEngine& gar::rec::dayoneconverter::fEngine

private

Definition at line 91 of file dayoneconverter_module.cc.

gar::geo::BitFieldCoder* gar::rec::dayoneconverter::fFieldDecoderMuID

private

Definition at line 95 of file dayoneconverter_module.cc.

gar::geo::BitFieldCoder* gar::rec::dayoneconverter::fFieldDecoderTrk

private

Definition at line 94 of file dayoneconverter_module.cc.

std::string gar::rec::dayoneconverter::fG4ModuleLabel

private

Use this to get the MCParticles.

Definition at line 89 of file dayoneconverter_module.cc.

const gar::geo::GeometryCore* gar::rec::dayoneconverter::fGeo

private

geometry information

Definition at line 93 of file dayoneconverter_module.cc.

bool gar::rec::dayoneconverter::fIncludeMuIDhits

private

Include MuID hits as TPCClusters.

Definition at line 76 of file dayoneconverter_module.cc.

std::string gar::rec::dayoneconverter::fInputEdepInstanceMuID

private

Input instance for MuID edeps.

Definition at line 75 of file dayoneconverter_module.cc.

std::string gar::rec::dayoneconverter::fInputEdepInstanceTPC

private

Input instance for TPC edeps.

Definition at line 74 of file dayoneconverter_module.cc.

std::string gar::rec::dayoneconverter::fInputEdepLabel

private

Input label for edeps.

Definition at line 73 of file dayoneconverter_module.cc.

TF1* gar::rec::dayoneconverter::fMu2e

private

Definition at line 96 of file dayoneconverter_module.cc.

float gar::rec::dayoneconverter::fPECm

private

conversion factor from cm step length to pe

Definition at line 80 of file dayoneconverter_module.cc.

float gar::rec::dayoneconverter::fRCut

private

Road in the YZ plane to add hits on a circle.

Definition at line 85 of file dayoneconverter_module.cc.

float gar::rec::dayoneconverter::fSmearLY

private

amount by which to smear the LY

Definition at line 81 of file dayoneconverter_module.cc.

float gar::rec::dayoneconverter::fSmearT

private

amount by which to smear T, in ns

Definition at line 79 of file dayoneconverter_module.cc.

float gar::rec::dayoneconverter::fSmearX

private

amount by which to smear X, in cm

Definition at line 77 of file dayoneconverter_module.cc.

float gar::rec::dayoneconverter::fSmearY

private

amount by which to smear Y, in cm

Definition at line 78 of file dayoneconverter_module.cc.

float gar::rec::dayoneconverter::fThrPE

private

threshold cut in pe

Definition at line 82 of file dayoneconverter_module.cc.

float gar::rec::dayoneconverter::fTimeBunch

private

Time bunch spread, in ns.

Definition at line 86 of file dayoneconverter_module.cc.

bool gar::rec::dayoneconverter::fUseOnlyTrueMuonHits

private

whether to cheat and use true muon hits

Definition at line 88 of file dayoneconverter_module.cc.

float gar::rec::dayoneconverter::fZCut1

private

Cut to ensure TPC Clusters are on different planes.

Definition at line 83 of file dayoneconverter_module.cc.

float gar::rec::dayoneconverter::fZCut2

private

Cut to ensure TPC clusters are on the same plane.

Definition at line 84 of file dayoneconverter_module.cc.

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

garsoft/Reco/dayoneconverter_module.cc

Public Member Functions

Private Member Functions

Private Attributes

Additional Inherited Members

Detailed Description

Constructor & Destructor Documentation

Member Function Documentation

Member Data Documentation