20 #include "TDatabasePDG.h" 38 #include "NovaDAQConventions/DAQConventions.h" 39 #include "DAQDataFormats/NanoSliceVersionConvention.h" 57 #include "Utilities/AssociationUtil.h" 105 std::unique_ptr< std::vector<sim::FLSHitList> >& mrccFLSHits,
106 std::unique_ptr< std::vector<sim::PhotonSignal> >& mrccPhotSig,
107 std::unique_ptr< std::vector<sim::Particle> >& mrccPartList,
108 std::unique_ptr< std::vector<simb::MCParticle> >& muonInfo) ;
209 produces< std::vector <rawdata::RawDigit> >();
210 produces< std::vector< simb::MCParticle> >();
211 produces< art::Assns<simb::MCParticle, rb::Cluster> >();
215 if(
fMuonID ==
"removeByTruth"){
216 produces< std::vector <sim::FLSHitList> >();
217 produces< std::vector <sim::PhotonSignal> >();
218 produces< std::vector <sim::Particle> >();
221 mf::LogInfo(
"MuonRemove") <<
" MuonRemove::MuonRemove()\n";
289 std::unique_ptr< std::vector<rawdata::RawDigit> >
290 muonRemovedDigits(
new std::vector<rawdata::RawDigit>);
291 std::unique_ptr< std::vector<simb::MCParticle > >
292 muonInfo(
new std::vector< simb::MCParticle>);
293 std::unique_ptr< art::Assns<simb::MCParticle, rb::Cluster> >
296 std::unique_ptr< std::vector<sim::Particle > >
297 mrccPartList(
new std::vector< sim::Particle>);
298 std::unique_ptr< std::vector<sim::FLSHitList> >
299 mrccFLSHits(
new std::vector<sim::FLSHitList>);
300 std::unique_ptr< std::vector<sim::PhotonSignal> >
301 mrccPhotSig(
new std::vector<sim::PhotonSignal>);
303 std::vector< rawdata::RawDigit> newDigits;
304 std::vector< rb::WeightedHit> allTrackHits;
309 int nRaw = digitcol->size();
310 std::list< rawdata::RawDigit >digs;
311 for(
int iRaw = 0; iRaw < nRaw; ++iRaw ){
313 digs.push_back( *(iDig));
317 std::list< rawdata::RawDigit >::iterator digIter, digIterBegin, digIterEnd;
323 for(
unsigned int i = 0;
i < slicecol->size(); ++
i){
341 int nSlices = slices.
size();
345 for(
int iSli = 0; iSli < nSlices; iSli++ ){
348 std::vector< art::Ptr<rb::Track> > alltracks = fmt.at(iSli);
349 std::vector< art::Ptr<rb::Track> > tracks;
354 for (
unsigned int trackId = 0; trackId < alltracks.size(); ++trackId){
355 if (fmBPFFitSums.isValid()){
356 std::vector<art::Ptr<rb::FitSum>> bpfFitSums = fmBPFFitSums.at(trackId);
357 if (
abs(bpfFitSums[0]->PDG()) == 13 ) tracks.push_back(alltracks[trackId]);
364 std::vector< rb::WeightedHit> trackHits;
368 double MuEnergy = -5;
370 if(
fMuonID ==
"removeByLength")
373 if(
fMuonID ==
"removeByEfficiency")
377 if(
fMuonID ==
"removeByTruth"){
381 std::cerr<<
" No MC truth info in this file.\n" 382 <<
" Choose a different method for muon removal.\n";
394 if(
fMuonID !=
"removeByTruth" && muonTrack.
NCell() > 0){
395 muonInfo->push_back(
MuonInfo(muonTrack, MuEnergy) );
397 slices[iSli], *assns);
400 if( muonTrack.
NCell() != 0 &&
fMuonID !=
"removeByTruth" )
403 allTrackHits.insert( allTrackHits.end(), trackHits.begin(), trackHits.end() );
409 std::vector< rb::WeightedHit >::iterator trkIterStart = allTrackHits.begin(), trkIter;
410 for(
auto digIter = digs.begin(); digIter != digs.end(); ++digIter){
414 int digPlane = cHit.
Plane();
415 int digCell = cHit.Cell();
416 int digTDC = cHit.TDC();
418 bool matched =
false;
420 for( trkIter = trkIterStart; trkIter != allTrackHits.end(); ++trkIter ){
424 if( digPlane < hit->Plane() ){
425 trkIterStart = trkIter;
428 else if( digPlane == hit->
Plane() &&
429 digCell < hit->
Cell() ){
430 trkIterStart = trkIter;
434 if(digPlane == hit->
Plane() &&
435 digCell == hit->
Cell() &&
436 digTDC == hit->
TDC() ){
438 float muonWeight = trkIter->weight;
439 float newADC = (1-muonWeight)*cHit.ADC();
442 trkIterStart = trkIter;
446 else if(muonWeight > 0 && muonWeight < 1){
452 dig.
SetADC(1, (1-muonWeight)*cHit.ADC(1));
453 dig.
SetADC(2, (1-muonWeight)*cHit.ADC(2));
462 int nadc = dig.
NADC();
463 for(
int iadc = 0; iadc < nadc; ++iadc){
464 int16_t newadc = ((dig.
ADC(iadc) -
baseline )*( 1-muonWeight) ) + baseline;
465 dig.
SetADC( iadc, newadc);
469 muonRemovedDigits->push_back( dig );
470 trkIterStart = trkIter;
476 muonRemovedDigits->push_back( *digIter );
480 std::cout<<
"raw digits before and after removal: " 481 <<nRaw<<
" "<<muonRemovedDigits->size()<<
std::endl;
484 if(
fMuonID ==
"removeByTruth"){
485 FillTruthInfo(evt, mrccFLSHits, mrccPhotSig, mrccPartList, muonInfo);
488 evt.
put(std::move(mrccPartList));
489 evt.
put(std::move(mrccPhotSig));
490 evt.
put(std::move(mrccFLSHits));
494 evt.
put(std::move(muonRemovedDigits));
495 evt.
put(std::move(muonInfo));
496 evt.
put(std::move(assns));
532 bool vtxValid =
false;
534 bool NProngValid =
false;
536 double muonEmax = -1;
540 if (fmElastic.isValid()){
541 std::vector<art::Ptr<rb::Vertex>> elastics = fmElastic.at(0);
542 if(elastics.size() > 0) {
546 std::vector< art::Ptr<rb::Prong> > prongs = fmFuzzyProng3D.at(0);
547 if (prongs.size() > 0) NProngValid =
true;
555 std::vector< art::Ptr<rb::Track> > tracks =
fmt.at(0);
556 for (
unsigned int i = 0;
i<tracks.size();
i++){
560 double muonE = muE.at(0)->E();
561 if (muonE > muonEmax) muonEmax = muonE;
566 if (fmBPFFitSums.isValid()) {
567 std::vector<art::Ptr<rb::FitSum>> bpfFitSums = fmBPFFitSums.at(0);
568 if (
abs(bpfFitSums[0]->PDG()) == 13) {
569 double muonE = bpfFitSums[0]->FourMom().E();
570 if (muonE > muonEmax) muonEmax = muonE;
578 bool ret = vtxValid && NProngValid && (muonEmax > 0);
592 double maxLength = 0;
594 for(
unsigned int i = 0;
i < tracks.size();
i++){
595 length = tracks[
i]->TotalLength();
596 if( length > maxLength){
598 longestTrack = *tracks[
i];
604 muonEnergy = muE.at(0)->E();
608 if (fmBPFFitSums.isValid()) {
609 std::vector<art::Ptr<rb::FitSum>> bpfFitSums = fmBPFFitSums.at(0);
610 muonEnergy = bpfFitSums[0]->FourMom().E();
632 std::vector<rb::WeightedHit> muonHits;
637 std::vector< int> primMuTrackIDs;
639 for(
int p = 0;
p < nPrimaries; ++
p){
643 int muMotherID = prim->
TrackId();
644 primMuTrackIDs.push_back( muMotherID);
649 if(part->
Mother() == muMotherID){
650 primMuTrackIDs.push_back( part->
TrackId() );
665 for(
int p = 0;
p < (
int)pSignal.size();
p++){
666 totPhotons += pSignal[
p].NPhoton();
667 int partID = pSignal[
p].TrackId();
670 for(
int i = 0;
i < (
int)primMuTrackIDs.size();
i++ ){
671 if( partID == primMuTrackIDs[
i] ){
672 muonPhotons += pSignal[
p].NPhoton();
678 if( muonPhotons == totPhotons ){
680 muonHits.push_back(newHit);
684 muonHits.push_back(newHit);
718 std::vector< std::set <int> > primMuTrackIDs;
720 for(
int n = 0;
n < nPrimaries; ++
n){
721 std::set <int> thisMuIDs;
726 int muMotherID = prim->
TrackId();
727 thisMuIDs.insert(muMotherID);
733 if (part->
Mother() == muMotherID )
734 thisMuIDs.insert( part->
TrackId() );
738 primMuTrackIDs.push_back(thisMuIDs);
744 for(
int p = 0;
p< (
int)primMuTrackIDs.size();
p++){
748 for(
int i = 0;
i< (
int)tracks.size();
i++){
757 muonTrack = *tracks[
i];
759 muonEnergy = muE.at(0)->E();
780 static TDatabasePDG* pdgt = TDatabasePDG::Instance();
781 TParticlePDG* pdgpm = pdgt->GetParticle(13);
783 const double muonMass = pdgpm->Mass();
788 double muonMom =
TMath::Sqrt((muonEnergy*muonEnergy) - (muonMass*muonMass));
790 TLorentzVector startPosition( muonTrack.
Start().X(), muonTrack.
Start().Y(),
793 TLorentzVector startMomentum( muonTrack.
Dir().X() * muonMom,
794 muonTrack.
Dir().Y() * muonMom ,
795 muonTrack.
Dir().Z() * muonMom, muonEnergy );
797 TLorentzVector stopPosition( muonTrack.
Stop().X(),
798 muonTrack.
Stop().Y(),
799 muonTrack.
Stop().Z(),
800 muonTrack.
Cell( muonTrack.
NCell() -1 )->TNS() );
801 TLorentzVector stopMomentum( 0,0,0,0 );
820 std::unique_ptr< std::vector<sim::FLSHitList> >& mrccFLSHits,
821 std::unique_ptr< std::vector<sim::PhotonSignal> >& mrccPhotSig,
822 std::unique_ptr< std::vector<sim::Particle> >& mrccPartList,
823 std::unique_ptr< std::vector<simb::MCParticle> >& muonInfo)
830 if (MCTruthList->empty())
834 std::vector< int> primMuTrackIDs;
842 muonInfo->push_back( *prim);
843 int muMotherID = prim->
TrackId();
844 primMuTrackIDs.push_back(muMotherID);
848 if( part->
Mother() == muMotherID ){
849 primMuTrackIDs.push_back( part->
TrackId() );
861 for(
size_t l = 0;
l < flslistHandle->size(); ++
l){
864 const std::vector<sim::FLSHit> &
hits = flslistHandle->at(
l).fHits;
866 for(
size_t f = 0;
f < hits.size(); ++
f){
867 bool hitIsFromMuon =
false;
869 for(
int id = 0;
id< (
int)primMuTrackIDs.size();
id++){
870 if(hits[
f].GetTrackID() == primMuTrackIDs[
id])
871 hitIsFromMuon =
true;
874 newFLSList.
fHits.push_back(hits[
f]);
878 mrccFLSHits->push_back(newFLSList);
887 for(
size_t i = 0;
i < phots->size(); ++
i){
889 bool sigIsFromMuon =
false;
891 for(
int id = 0;
id<(
int)primMuTrackIDs.size();
id++){
892 if( phot.
TrackId() == primMuTrackIDs[
id])
893 sigIsFromMuon =
true;
897 mrccPhotSig->push_back(phot);
903 for (
unsigned int s = 0;
s < MCTruthList->size();
s++){
910 for (
unsigned int i = 0;
i < particleList.size(); ++
i){
911 bool partIsFromMuon =
false;
912 for(
int id = 0;
id<(
int)primMuTrackIDs.size();
id++){
913 if(particleList[
i]->TrackId() == primMuTrackIDs[
id])
914 partIsFromMuon =
true;
917 mrccPartList->push_back(*particleList[
i]);
929 float maxPidVal = -1;
933 for(
int i = 0;
i < (
int)tracks.size();
i++){
936 float pidVal = pid->
Value();
937 if( pidVal > maxPidVal ){
939 mostMuLikeTrack = *tracks[
i];
941 muonEnergy = muE.at(0)->E();
946 std::cerr<<
" There are no remid objects in this file.\n" 947 <<
" Choose another method of muon removal.\n";
951 return mostMuLikeTrack;
966 float maxPidVal = -1;
971 std::vector<art::Ptr<rb::Vertex>> elastics = fmElastic.at(0);
972 if( elastics.size() > 0 ){
975 std::vector< art::Ptr<rb::Prong> > prongs = fmFuzzyProng3D.at(0);
978 for (
unsigned int iPng = 0; iPng < prongs.size(); ++iPng ) {
980 std::vector<art::Ptr<rb::PID>> cvnparts = fmCVNParticleResult.at(iPng);
981 for (
unsigned int j=0;
j<cvnparts.size(); ++
j ){
982 if (
abs(cvnparts[
j]->Pdg()) == 13 ){
983 double thisPid = cvnparts[
j]->Value();
984 if ( thisPid > maxPidVal) {
986 mostMuLikeProng = prongs[iPng];
991 if ( cvnparts.size() == 0 && fmCVNParticleResult.isValid() && prong->
TotalLength() > 500){
993 mostMuLikeProng = prongs[iPng];
1000 std::vector< art::Ptr<rb::Track> > MuLikeTracks = fmMuLikeTracks.at(0);
1003 for (
unsigned int trackId = 0; trackId < MuLikeTracks.size(); ++trackId){
1004 if (fmBPFFitSums.isValid()){
1005 std::vector<art::Ptr<rb::FitSum>> bpfFitSums = fmBPFFitSums.at(trackId);
1006 if (
abs(bpfFitSums[0]->PDG()) == 13 ) {
1007 muonEnergy = bpfFitSums[0]->FourMom().E();
1008 mostMuLikeTrack=*MuLikeTracks[trackId];
1013 return mostMuLikeTrack;
::xsd::cxx::tree::id< char, ncname > id
back track the reconstruction to the simulation
static bool CreateAssn(art::EDProducer const &prod, art::Event &evt, std::vector< T > &a, art::Ptr< U > b, art::Assns< T, U > &assn, size_t indx=UINT_MAX, std::string const &instance=std::string())
Create a 1 to 1 association between a new product and one already in the event.
bool fTrackCleanUp
A FHiCL flag to turn of track clean up. Default if true.
std::vector< sim::FLSHit > fHits
std::vector< rb::WeightedHit > RemoveByTruth(const art::Ptr< rb::Cluster > &slice)
Returns a vector of Cellhits attributed to the muon by checking with truth. For implementing muon rem...
unsigned int NCell(geo::View_t view) const
Number of cells in view view.
void produce(art::Event &evt)
void AddTrajectoryPoint(TLorentzVector const &position, TLorentzVector const &momentum)
int32_t TDC() const
The time of the last baseline sample.
const sim::ParticleNavigator & ParticleNavigator() const
Get a reference to the ParticleNavigator.
std::string fFuzzy3DInput
Label of module that produced prongs.
MaybeLogger_< ELseverityLevel::ELsev_info, false > LogInfo
unsigned short Plane() const
std::string fCVNLabel
Label of module that produced slice CVN PID.
A collection of associated CellHits.
rb::Track RemoveByProngCVN(std::vector< art::Ptr< rb::Track > > &tracks, const art::Ptr< rb::Cluster > &slice, art::Event &evt, double &muonEnergy)
Returns the BPF muon track with the highest Prong CVN value.
std::string fKalmanTrackModuleLabel
Label of module that produced tracks.
A rb::Prong with full reconstructed trajectory.
std::string fEnergyLabelBPF
Label of module that produced Energy object (in case of BPF tracks)
DEFINE_ART_MODULE(TestTMapFile)
virtual TVector3 Start() const
std::string fCVNParticleLabel
Label of module that produced Prong CVN PIDs.
double fNumuPIDCut
Numu PID cut.
rb::Track LongestTrack(std::vector< art::Ptr< rb::Track > > &tracks, art::Event &evt, double &muonEnergy)
A function that returns the longest track amongst a vector of tracks.
rb::CellHit MakeCellHit(const rawdata::RawDigit *rawdigit)
void FillTruthInfo(const art::Event &evt, std::unique_ptr< std::vector< sim::FLSHitList > > &mrccFLSHits, std::unique_ptr< std::vector< sim::PhotonSignal > > &mrccPhotSig, std::unique_ptr< std::vector< sim::Particle > > &mrccPartList, std::unique_ptr< std::vector< simb::MCParticle > > &muonInfo)
Saves muon removed FLSHitLists, PhotonSignals, Particle list and information about the removed muon i...
art::PtrVector< rb::CellHit > AllCells() const
Get all cells from both views.
bool IsNoise(const art::Ptr< rb::CellHit > &hit) const
Is this hit not associated with any particles?
double HitCollectionEfficiency(const std::set< int > &trackIDs, const std::vector< rb::WeightedHit > &whits, const std::vector< const rb::CellHit * > &allhits, const geo::View_t &view, std::map< int, double > *effMap=0, bool energyEff=false, double *desiredEnergy=0, double *totalEnergy=0, int *desiredHits=0, int *totalHits=0) const
Returns the fraction of all energy in an event from a specific set of Geant4 track IDs that are repre...
bool SliceSelection(const art::Event &evt, const art::Ptr< rb::Cluster > &slice)
Function for preselecting the slices for the further mrcc process.
fhicl::ParameterSet fTrkCleanUpPSet
ParameterSet of TrackCleanUpAlg.
ProductID put(std::unique_ptr< PROD > &&product)
unsigned short Cell() const
void SetADC(uint32_t i, int16_t iADC)
virtual double TotalLength() const
Distance along prong to reach last cell hit.
Far Detector at Ash River, MN.
int fDCSDistance
DCS Distance: multipoint readout look back.
std::string fClusterModuleLabel
Label of module that produced slice clusters.
const sim::Particle * Primary(const int) const
void beginRun(art::Run &run)
void push_back(Ptr< U > const &p)
bool CompareDigitsByPlaneAndCell(rawdata::RawDigit ra, rawdata::RawDigit rb)
Compare two rawdigits by plane. Return true if plane of ra is greater then the plane of rb...
std::string fBPFTrackModuleLabel
Label of module that produced BPF tracks.
Prototype Near Detector on the surface at FNAL.
T get(std::string const &key) const
novadaq::cnv::DetId DetId() const
Prefer ds::DetectorService::DetId() instead.
virtual TVector3 Dir() const
Unit vector describing prong direction.
An algorithm to determine the energy of a slice.
Near Detector in the NuMI cavern.
unsigned short GetPlane(const rawdata::RawDigit *dig)
bool IsFiltered(const art::Event &evt, art::Ptr< T > x, const std::vector< std::string > &labels)
Is this Ptr marked "filtered out"?
float fADCThreshold_Far
ADC thresholds.
void SortDigitsByPlaneAndCell(std::list< rawdata::RawDigit > &digitcol)
Takes in a vector of raw digits and sorts it by plane.
Vertex location in position and time.
Perform a "2 point" Hough transform on a collection of hits.
rb::Track RemoveByEfficiency(std::vector< art::Ptr< rb::Track > > &tracks, const rb::Cluster &slice, art::Event &evt, double &muonEnergy)
Returns the highest efficiency (by truth) muon track for removal;.
std::vector< rb::WeightedHit > CleanUpTrack(const rb::Track &muonTrack, const rb::Cluster &slice, double *mipRangeHigh=NULL, double *mipRangeLow=NULL, double *mipValue=NULL, double *vertexRegionDeDxCutOff=NULL)
A function to disentangle the muon track from the hadronic energy close to the vertex. Returns a vector of rb::WeightedHits where the hit is the cellhits that were on the given track and weight is the fraction of the energy in the hit that came from the track. mipRangeLow, mipRangeHigh, mipValue and vertexRegionDeDxCutOff are pointers that point to the value of mip range limits, peak values andthe dedx cut off used to find the vertex region. If these pointers are not provided, the default values are used. The default values are good for muons, not for other particles.
bool CompareByPlaneAndCell(rb::WeightedHit whita, rb::WeightedHit whitb)
Compare two rawdigits by plane. Return true if plane of ra is greater then the plane of rb...
art::Ptr< rb::CellHit > Cell(geo::View_t view, unsigned int viewIdx) const
Get the ith cell from view view.
const std::vector< sim::PhotonSignal > HitToPhotonSignal(const art::Ptr< rb::CellHit > &hit) const
Returns the PhotonSignals contributing the signal in the specified hit.
A rawdata::RawDigit with channel information decoded.
std::string fEnergyLabel
Label of module that produced Energy object (in case of Kalman tracks)
numue::NumuEAlg * fNumuEAlg
An object of NumuEAlg.
rb::Track RemoveByReMId(std::vector< art::Ptr< rb::Track > > &tracks, art::Event &evt, double &muonEnergy)
Returns the track with the highest PID value returned by the ReMId package.
int16_t ADC(uint32_t i) const
A vector of FLSHit from single neutrino interaction.
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Var Sqrt(const Var &v)
Use to take sqrt of a var.
MuonRemove(fhicl::ParameterSet const &pset)
fhicl::ParameterSet fNumuEAlgPSet
ParameterSet of NumuEAlg.
art::Ptr< rb::CellHit > hit
bool HaveTruthInfo() const
Is this a file with truth info in? (Is BackTracker going to be any use to you?)
std::string fMCTruthModuleLabel
Label of module that produced MC truth information.
Simple hit+weight pair, returned from rb::Cluster::WeightedHits.
std::string fMuonID
Label of module that produced muon pid object.
assert(nhit_max >=nhit_nbins)
simb::MCParticle MuonInfo(rb::Track muonTrack, double muonEnergy)
Saves the information about the removed muon in a simb::MCParticle object The energy of the muon is o...
std::string fElasticInput
Label of module that produced verteces.
unsigned short GetCell(const rawdata::RawDigit *dig)
std::string fNumuReMIdLabel
Label of module that produced Numu PID.
void reconfigure(const fhicl::ParameterSet &pset)
void SortByPlaneAndCell(std::vector< rb::WeightedHit > &whits)
TVector3 Stop() const
Position of the final trajectory point.
TrackCleanUpAlg * fTrkCleanUpAlg
An object of TrackCleanUpAlg.
Encapsulate the geometry of one entire detector (near, far, ndos)
void SortByPlane(std::vector< art::Ptr< rb::CellHit > > &c)
Sort c in plane order (low to high).
std::string fWhichTrackModule
Label of module which will be used for track loading.
std::string fRawDataLabel
Label of module that produced raw digits.
art::ServiceHandle< geo::Geometry > fgeom
Handle to Geometry.
int NumberOfPrimaries() const
mapped_type Particle(const size_type) const
std::vector< const sim::Particle * > MCTruthToParticles(art::Ptr< simb::MCTruth > const &mct) const