#include <DTMeantimerPatternReco.h>

Inheritance diagram for DTMeantimerPatternReco:

Public Member Functions
virtual std::string	algoName () const
	return the algo name More...

	DTMeantimerPatternReco (const edm::ParameterSet &pset)
	Constructor. More...

virtual edm::OwnVector < DTSLRecSegment2D >	reconstruct (const DTSuperLayer *sl, const std::vector< DTRecHit1DPair > &hits)
	this function is called in the producer More...

virtual void	setES (const edm::EventSetup &setup)

virtual	~DTMeantimerPatternReco ()
	Destructor. More...

Public Member Functions inherited from DTRecSegment2DBaseAlgo
	DTRecSegment2DBaseAlgo (const edm::ParameterSet &)
	Constructor. More...

virtual	~DTRecSegment2DBaseAlgo ()
	Destructor. More...

Private Member Functions
void	addHits (const DTSuperLayer sl, std::vector< DTSegmentCand::AssPoint > &assHits, const std::vector< std::shared_ptr< DTHitPairForFit >> &hits, std::vector< DTSegmentCand > &result)

std::vector< DTSegmentCand * >	buildSegments (const DTSuperLayer *sl, const std::vector< std::shared_ptr< DTHitPairForFit >> &hits)

bool	checkDoubleCandidates (std::vector< DTSegmentCand * > &segs, DTSegmentCand *seg)

std::unique_ptr< DTSegmentCand >	fitWithT0 (const DTSuperLayer *sl, const std::vector< DTSegmentCand::AssPoint > &assHits, double &chi2, double &t0_corr, const bool fitdebug)

bool	geometryFilter (const DTWireId first, const DTWireId second) const

std::vector< std::shared_ptr < DTHitPairForFit > >	initHits (const DTSuperLayer *sl, const std::vector< DTRecHit1DPair > &hits)

void	printPattern (std::vector< DTSegmentCand::AssPoint > &assHits, const DTHitPairForFit *hit)

Private Attributes
bool	debug

unsigned int	maxfound

std::string	theAlgoName

double	theAlphaMaxPhi

double	theAlphaMaxTheta

DTSegmentCleaner *	theCleaner

edm::ESHandle< DTGeometry >	theDTGeometry

DTLinearFit *	theFitter

unsigned int	theMaxAllowedHits

double	theMaxChi2

DTSegmentUpdator *	theUpdator

Friends
class	DTMeantimerPatternReco4D

Detailed Description

Algo for reconstructing 2d segment in DT using a combinatorial approach with a T0 estimation produced along the way

Author: Stefano Lacaprara - INFN Legnaro stefa.nosp@m.no.l.nosp@m.acapr.nosp@m.ara@.nosp@m.pd.in.nosp@m.fn.i.nosp@m.t; Riccardo Bellan - INFN TO ricca.nosp@m.rdo..nosp@m.bella.nosp@m.n@ce.nosp@m.rn.ch; Piotr Traczyk - SINS Warsaw ptrac.nosp@m.zyk@.nosp@m.fuw.e.nosp@m.du.p.nosp@m.l

Definition at line 44 of file DTMeantimerPatternReco.h.

Constructor & Destructor Documentation

DTMeantimerPatternReco::DTMeantimerPatternReco ( const edm::ParameterSet & pset )

Constructor.

Definition at line 39 of file DTMeantimerPatternReco.cc.

References debug, edm::ParameterSet::getParameter(), edm::ParameterSet::getUntrackedParameter(), theAlphaMaxPhi, theAlphaMaxTheta, theCleaner, theMaxAllowedHits, theMaxChi2, and theUpdator.

                                                                           : 
   DTRecSegment2DBaseAlgo(pset), 
   theFitter(new DTLinearFit()),
   theAlgoName("DTMeantimerPatternReco")
 {
   theMaxAllowedHits = pset.getParameter<unsigned int>("MaxAllowedHits"); // 100
   theAlphaMaxTheta = pset.getParameter<double>("AlphaMaxTheta");// 0.1 ;
   theAlphaMaxPhi = pset.getParameter<double>("AlphaMaxPhi");// 1.0 ;
   theMaxChi2 = pset.getParameter<double>("MaxChi2");// 8.0 ;
   debug = pset.getUntrackedParameter<bool>("debug");
   theUpdator = new DTSegmentUpdator(pset);
   theCleaner = new DTSegmentCleaner(pset);
 }

DTMeantimerPatternReco::~DTMeantimerPatternReco ( )

virtual

Destructor.

Definition at line 54 of file DTMeantimerPatternReco.cc.

References theCleaner, theFitter, and theUpdator.

                                                 {
   delete theUpdator;
   delete theCleaner;
   delete theFitter;
 }

Member Function Documentation

void DTMeantimerPatternReco::addHits	(	const DTSuperLayer *	sl,
		std::vector< DTSegmentCand::AssPoint > &	assHits,
		const std::vector< std::shared_ptr< DTHitPairForFit >> &	hits,
		std::vector< DTSegmentCand * > &	result
	)

private

Definition at line 200 of file DTMeantimerPatternReco.cc.

References checkDoubleCandidates(), gather_cfg::cout, debug, fitWithT0(), geometryFilter(), DTEnums::Left, maxfound, and DTEnums::Right.

Referenced by buildSegments().

                                                                 {
 
   double chi2l,chi2r,t0_corrl,t0_corrr;
   bool foundSomething = false;
 
 //  if (debug) 
 //    cout << "DTMeantimerPatternReco::addHit " << endl << "   Picked " << assHits.size() << " hits, " << hits.size() << " left." << endl;
   
   if (assHits.size()+hits.size()<maxfound) return;
 
   // loop over the remaining hits
   for (hitCont::const_iterator hit=hits.begin(); hit!=hits.end(); ++hit) {
 
 //    if (debug) {
 //      int nHits=assHits.size()+1;
 //      cout << "     Trying B: " << **hit<< " wire: " << (*hit)->id() << endl;
 //      printPattern(assHits,*hit);
 //    }
 
     assHits.push_back(DTSegmentCand::AssPoint(*hit, DTEnums::Left));
     std::unique_ptr<DTSegmentCand> left_seg = fitWithT0(sl,assHits, chi2l, t0_corrl,0);
     assHits.pop_back();
 //    if (debug) 
 //      cout << "    Left:  t0= " << t0_corrl << "  chi2/nHits= " << chi2l << "/" << nHits << "  ok: " << left_ok << endl;
 
     assHits.push_back(DTSegmentCand::AssPoint(*hit, DTEnums::Right));
     std::unique_ptr<DTSegmentCand> right_seg = fitWithT0(sl,assHits, chi2r, t0_corrr,0);
     assHits.pop_back();
 //    if (debug) 
 //      cout << "   Right:  t0= " << t0_corrr << "  chi2/nHits= " << chi2r << "/" << nHits << "  ok: " << right_ok << endl;
 
     bool left_ok=(left_seg)?true:false;
     bool right_ok=(right_seg)?true:false;
 
     if (!left_ok && !right_ok) continue;
 
     foundSomething = true;    
 
     // prepare the hit set for the next search, start from the other side                        
     hitCont hitsForFit;
     for (hitCont::const_iterator tmpHit=hit+1; tmpHit!=hits.end(); tmpHit++) 
       if (geometryFilter((*tmpHit)->id(),(*hit)->id())) hitsForFit.push_back(*tmpHit); 
 
     reverse(hitsForFit.begin(),hitsForFit.end());
 
     // choose only one - left or right
     if (assHits.size()>3 && left_ok && right_ok) {
       if (chi2l<chi2r-0.1) right_ok=false; else
         if (chi2r<chi2l-0.1) left_ok=false;
     }
     if (left_ok) { 
       assHits.push_back(DTSegmentCand::AssPoint(*hit, DTEnums::Left));
       addHits(sl,assHits,hitsForFit,result);
       assHits.pop_back();
     }
     if (right_ok) { 
       assHits.push_back(DTSegmentCand::AssPoint(*hit, DTEnums::Right));
       addHits(sl,assHits,hitsForFit,result);
       assHits.pop_back();
     }
   }
 
   if (foundSomething) return;
   // if we didn't find any new hits compatible with the current candidate, we proceed to check and store the candidate
 
   // If we already have a segment with more hits from this hit pair - don't save this one.  
   if (assHits.size()<maxfound) return;
 
   // Check if semgent Ok, calculate chi2
   std::unique_ptr<DTSegmentCand> seg = fitWithT0(sl,assHits, chi2l, t0_corrl,debug);
   if (!seg) return;
   
   if (!seg->good()) {
 //    if (debug) cout << "   Segment not good() - skipping" << endl;
     return;
   }
 
   if (assHits.size()>maxfound) maxfound = assHits.size();
   if (debug) cout << endl << "   Seg t0= " << t0_corrl << endl << *seg<< endl;
   
   if (checkDoubleCandidates(result,seg.get())) {
     result.push_back(seg.release());
     if (debug) cout << "   Result is now " << result.size() << endl;
   } else {
     if (debug) cout << "   Exists - skipping" << endl;
   }
 }

virtual std::string DTMeantimerPatternReco::algoName ( void ) const

inlinevirtual

return the algo name

Implements DTRecSegment2DBaseAlgo.

Definition at line 62 of file DTMeantimerPatternReco.h.

References theAlgoName.

62 { return theAlgoName; }

DTMeantimerPatternReco::theAlgoName

std::string theAlgoName

Definition: DTMeantimerPatternReco.h:109

vector< DTSegmentCand * > DTMeantimerPatternReco::buildSegments	(	const DTSuperLayer *	sl,
		const std::vector< std::shared_ptr< DTHitPairForFit >> &	hits
	)

private

Definition at line 105 of file DTMeantimerPatternReco.cc.

References addHits(), DTSegmentCleaner::clean(), gather_cfg::cout, debug, geometryFilter(), DTSuperLayer::id(), hit::id, listHistos::IP, DTEnums::Left, maxfound, query::result, DTEnums::Right, theAlphaMaxPhi, theAlphaMaxTheta, theCleaner, theMaxAllowedHits, PV3DBase< T, PVType, FrameType >::theta(), and GeomDet::toGlobal().

Referenced by DTMeantimerPatternReco4D::buildPhiSuperSegmentsCandidates(), and reconstruct().

                                                                                             {
 
   vector<DTSegmentCand*> result;
   DTEnums::DTCellSide codes[2]={DTEnums::Left, DTEnums::Right};
 
   if(debug) {
     cout << "buildSegments: " << sl->id() << " nHits " << hits.size() << endl;
     for (hitIter hit=hits.begin(); hit!=hits.end(); ++hit) cout << **hit<< " wire: " << (*hit)->id() << " DigiTime: " << (*hit)->digiTime() << endl;
   }
 
   if (hits.size() > theMaxAllowedHits ) {
     if(debug) {
       cout << "Warning: this SuperLayer " << sl->id() << " has too many hits : "
         << hits.size() << " max allowed is " << theMaxAllowedHits << endl;
       cout << "Skipping segment reconstruction... " << endl;
     }
     return result;
   }
 
   GlobalPoint IP;
   float DAlphaMax;
   if ((sl->id()).superlayer()==2)  // Theta SL
     DAlphaMax=theAlphaMaxTheta;
   else // Phi SL
     DAlphaMax=theAlphaMaxPhi;
   
   // get two hits in different layers and see if there are other hits
   //  compatible with them
   for (hitCont::const_iterator firstHit=hits.begin(); firstHit!=hits.end();
        ++firstHit) {
     for (hitCont::const_reverse_iterator lastHit=hits.rbegin(); 
          (*lastHit)!=(*firstHit); ++lastHit) {
       
       // a geometrical sensibility cut for the two hits
       if (!geometryFilter((*firstHit)->id(),(*lastHit)->id())) continue;
      
       // create a set of hits for the fit (only the hits between the two selected ones)
       hitCont hitsForFit;
       for (hitCont::const_iterator tmpHit=firstHit+1; (*tmpHit)!=(*lastHit); tmpHit++) 
         if ((geometryFilter((*tmpHit)->id(),(*lastHit)->id())) 
             && (geometryFilter((*tmpHit)->id(),(*firstHit)->id()))) hitsForFit.push_back(*tmpHit);
 
       for (int firstLR=0; firstLR<2; ++firstLR) {
         for (int lastLR=0; lastLR<2; ++lastLR) {
 
       // TODO move the global transformation in the DTHitPairForFit class
       // when it will be moved I will able to remove the sl from the input parameter
       GlobalPoint gposFirst=sl->toGlobal( (*firstHit)->localPosition(codes[firstLR]) );
       GlobalPoint gposLast= sl->toGlobal( (*lastHit)->localPosition(codes[lastLR]) );
           GlobalVector gvec=gposLast-gposFirst;
           GlobalVector gvecIP=gposLast-IP;
 
           // difference in angle measured
           float DAlpha=fabs(gvec.theta()-gvecIP.theta());
           if (DAlpha>DAlphaMax) continue;
       
 //    if(debug) {
 //              cout << "Selected hit pair:" << endl;
 //              cout << "  First " << *(*firstHit) << " Layer Id: " << (*firstHit)->id().layerId() << " Side: " << firstLR << " DigiTime: " << (*firstHit)->digiTime() << endl;
 //              cout << "  Last "  << *(*lastHit)  << " Layer Id: " << (*lastHit)->id().layerId()  << " Side: " << lastLR << " DigiTime: " << (*lastHit)->digiTime() <<  endl;
 //          }
         
           vector<DTSegmentCand::AssPoint> assHits;
           // create a candidate hit list
           assHits.push_back(DTSegmentCand::AssPoint(*firstHit,codes[firstLR]));
           assHits.push_back(DTSegmentCand::AssPoint(*lastHit,codes[lastLR]));
 
           // run hit adding/segment building 
           maxfound = 3;
           addHits(sl,assHits,hitsForFit,result);
         }
       }
     }
    }
 
   // now I have a couple of segment hypotheses, should check for ghosts
   if (debug) {
     cout << "Result (before cleaning): " << result.size() << endl;
     for (vector<DTSegmentCand*>::const_iterator seg=result.begin(); seg!=result.end(); ++seg) cout << *(*seg) << endl;
   }        
 
   result = theCleaner->clean(result);
 
   if (debug) {
     cout << "Result (after cleaning): " << result.size() << endl;
     for (vector<DTSegmentCand*>::const_iterator seg=result.begin(); seg!=result.end(); ++seg) cout << *(*seg) << endl;
   }
 
   return result;
 }

bool DTMeantimerPatternReco::checkDoubleCandidates	(	std::vector< DTSegmentCand * > &	segs,
		DTSegmentCand *	seg
	)

private

Definition at line 362 of file DTMeantimerPatternReco.cc.

References DTSegmentCand::chi2ndof(), and DTSegmentCand::nHits().

Referenced by addHits().

                                                                   {
   for (vector<DTSegmentCand*>::iterator cand=cands.begin();
        cand!=cands.end(); ++cand) {
     if (*(*cand)==*seg) return false;
     if (((*cand)->nHits()>=seg->nHits()) && ((*cand)->chi2ndof()<seg->chi2ndof()))
       if ((*cand)->nSharedHitPairs(*seg)>int(seg->nHits()-2)) return false;
   }
   return true;
 }

std::unique_ptr< DTSegmentCand > DTMeantimerPatternReco::fitWithT0	(	const DTSuperLayer *	sl,
		const std::vector< DTSegmentCand::AssPoint > &	assHits,
		double &	chi2,
		double &	t0_corr,
		const bool	fitdebug
	)

private

Definition at line 327 of file DTMeantimerPatternReco.cc.

References DTSegmentUpdator::fit(), theMaxChi2, and theUpdator.

Referenced by addHits().

 {
   // create a DTSegmentCand 
   DTSegmentCand::AssPointCont pointsSet;
   pointsSet.insert(assHits.begin(),assHits.end());
   std::unique_ptr<DTSegmentCand> seg(new DTSegmentCand(pointsSet,sl));
   
   // perform the 3 parameter fit on the segment candidate
   theUpdator->fit(seg.get(),1,fitdebug);
 
   chi2 = seg->chi2();
   t0_corr = seg->t0();
   
   // sanity check - drop segment candidates with a failed fit
   // for a 3-par fit this includes segments with hits after the calculated t0 correction ending up
   // beyond the chamber walls or on the other side of the wire
   if (chi2==-1.) return nullptr;
 
   // at this point we keep all 3-hit segments that passed the above check
   if (assHits.size()==3) return seg;
 
   // for segments with no t0 information we impose a looser chi2 cut
   if (t0_corr==0) {
     if (chi2<100.) return seg;
       else return nullptr; 
   }
 
   // cut on chi2/ndof of the segment candidate
   if ((chi2/(assHits.size()-3)<theMaxChi2)) return seg;
     else return nullptr;
 }

bool DTMeantimerPatternReco::geometryFilter	(	const DTWireId	first,
		const DTWireId	second
	)		const

private

Definition at line 291 of file DTMeantimerPatternReco.cc.

References funct::abs(), gather_cfg::cout, DTLayerId::layer(), DTWireId::layerId(), DTSuperLayerId::superLayer(), DTLayerId::superlayerId(), and DTWireId::wire().

Referenced by addHits(), and buildSegments().

 {
 //  return true;
 
   const int layerLowerCut[4]={0,-1,-2,-2};
   const int layerUpperCut[4]={0, 2, 2, 3};
 //  const int layerLowerCut[4]={0,-2,-4,-5};
 //  const int layerUpperCut[4]={0, 3, 4, 6};
 
   // deal only with hits that are in the same SL
   if (first.layerId().superlayerId().superLayer()!=second.layerId().superlayerId().superLayer()) 
     return true;
     
   int deltaLayer=abs(first.layerId().layer()-second.layerId().layer());
 
   // drop hits in the same layer
   if (!deltaLayer) return false;
 
   // protection against unexpected layer numbering
   if (deltaLayer>3) { 
     cout << "*** WARNING! DT Layer numbers differ by more than 3! for hits: " << endl;
     cout << "             " << first << endl;
     cout << "             " << second << endl;
     return false;
   }
 
   // accept only hits in cells "not too far away"
   int deltaWire=first.wire()-second.wire();
   if (second.layerId().layer()%2==0) deltaWire=-deltaWire; // yet another trick to get it right...
   if ((deltaWire<=layerLowerCut[deltaLayer]) || (deltaWire>=layerUpperCut[deltaLayer])) return false;
 
   return true;
 }

vector< std::shared_ptr< DTHitPairForFit > > DTMeantimerPatternReco::initHits	(	const DTSuperLayer *	sl,
		const std::vector< DTRecHit1DPair > &	hits
	)

private

Definition at line 93 of file DTMeantimerPatternReco.cc.

References query::result, and theDTGeometry.

Referenced by DTMeantimerPatternReco4D::buildPhiSuperSegmentsCandidates(), and reconstruct().

                                                                            {  
   
   hitCont result;
   for (vector<DTRecHit1DPair>::const_iterator hit=hits.begin();
        hit!=hits.end(); ++hit) {
     result.push_back(std::make_shared<DTHitPairForFit>(*hit, *sl, theDTGeometry));
   }
   return result;
 }

void DTMeantimerPatternReco::printPattern	(	std::vector< DTSegmentCand::AssPoint > &	assHits,
		const DTHitPairForFit *	hit
	)

private

Definition at line 376 of file DTMeantimerPatternReco.cc.

References gather_cfg::cout, i, DTEnums::Left, and cuy::mark.

                                                                                                         {
 
   char mark[26]={". . . . . . . . . . . . "};
   int wire[12]={0,0,0,0,0,0,0,0,0,0,0,0};
 
   for (vector<DTSegmentCand::AssPoint>::const_iterator assHit=assHits.begin(); assHit!=assHits.end(); ++assHit) {
     int lay  = (((*assHit).first)->id().superlayerId().superLayer()-1)*4 + ((*assHit).first)->id().layerId().layer()-1;
     wire[lay]= ((*assHit).first)->id().wire();
     if ((*assHit).second==DTEnums::Left) mark[lay*2]='L'; else mark[lay*2]='R';
   }
 
   int lay  = ((*hit).id().superlayerId().superLayer()-1)*4 + (*hit).id().layerId().layer() - 1;
   wire[lay]= (*hit).id().wire();
   mark[lay*2]='*';
 
   cout << "   " << mark << endl << "  ";
   for (int i=0; i<12; i++) if (wire[i]) cout << setw(2) << wire[i]; else cout << "  ";
   cout << endl;
 }

edm::OwnVector< DTSLRecSegment2D > DTMeantimerPatternReco::reconstruct	(	const DTSuperLayer *	sl,
		const std::vector< DTRecHit1DPair > &	hits
	)

virtual

this function is called in the producer

Implements DTRecSegment2DBaseAlgo.

Definition at line 62 of file DTMeantimerPatternReco.cc.

References buildSegments(), objects.IsoTrackAnalyzer::candidates, gather_cfg::cout, debug, initHits(), edm::OwnVector< T, P >::push_back(), query::result, theUpdator, and DTSegmentUpdator::update().

Referenced by DTMeantimerPatternReco4D::reconstruct().

                                                                                {
 
   edm::OwnVector<DTSLRecSegment2D> result;
   std::vector<std::shared_ptr<DTHitPairForFit>> hitsForFit = initHits(sl, pairs);
 
   vector<DTSegmentCand*> candidates = buildSegments(sl, hitsForFit);
 
   vector<DTSegmentCand*>::const_iterator cand=candidates.begin();
   while (cand<candidates.end()) {
     
     DTSLRecSegment2D *segment = (**cand);
     theUpdator->update(segment,1);
 
     if (debug) 
       cout<<"Reconstructed 2D segments "<<*segment<<endl;
     result.push_back(segment);
     
     delete *(cand++); // delete the candidate!
   }
 
   return result;
 }

void DTMeantimerPatternReco::setES ( const edm::EventSetup & setup )

virtual

Through this function the EventSetup is percolated to the objs which request it

Implements DTRecSegment2DBaseAlgo.

Definition at line 86 of file DTMeantimerPatternReco.cc.

References edm::EventSetup::get(), DTSegmentUpdator::setES(), theDTGeometry, and theUpdator.

Referenced by DTMeantimerPatternReco4D::setES().

                                                             {
   // Get the DT Geometry
   setup.get<MuonGeometryRecord>().get(theDTGeometry);
   theUpdator->setES(setup);
 }