#include <DTMeantimerPatternReco.h>

Inheritance diagram for DTMeantimerPatternReco:

Public Member Functions
virtual std::string	algoName () const
	return the algo name
	DTMeantimerPatternReco (const edm::ParameterSet &pset)
	Constructor.
virtual edm::OwnVector < DTSLRecSegment2D >	reconstruct (const DTSuperLayer *sl, const std::vector< DTRecHit1DPair > &hits)
	this function is called in the producer
virtual void	setES (const edm::EventSetup &setup)
virtual	~DTMeantimerPatternReco ()
	Destructor.
Private Types
typedef std::pair < DTHitPairForFit *, DTEnums::DTCellSide >	AssPoint
Private Member Functions
void	addHits (const DTSuperLayer sl, std::vector< AssPoint > &assHits, const std::vector< DTHitPairForFit > &hits, std::vector< DTSegmentCand * > &result, std::vector< AssPoint > &usedHits)
std::vector< DTSegmentCand * >	buildSegments (const DTSuperLayer sl, const std::vector< DTHitPairForFit > &hits)
bool	checkDoubleCandidates (std::vector< DTSegmentCand * > &segs, DTSegmentCand *seg)
bool	fitWithT0 (const std::vector< AssPoint > &assHits, double &chi2, double &t0_corr, const bool fitdebug)
bool	geometryFilter (const DTWireId first, const DTWireId second) const
std::vector< DTHitPairForFit * >	initHits (const DTSuperLayer *sl, const std::vector< DTRecHit1DPair > &hits)
void	rawFit (double &a, double &b, const std::vector< std::pair< double, double > > &hits)
Private Attributes
bool	debug
unsigned int	maxfound
std::string	theAlgoName
double	theAlphaMaxPhi
double	theAlphaMaxTheta
DTSegmentCleaner *	theCleaner
edm::ESHandle< DTGeometry >	theDTGeometry
unsigned int	theMaxAllowedHits
double	theMaxChi2
double	theMaxT0
double	theMinT0
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

Date:: 2008/03/10 11:18:20

Revision:: 1.2

Author:: Stefano Lacaprara - INFN Legnaro <stefano.lacaprara@pd.infn.it>; Riccardo Bellan - INFN TO <riccardo.bellan@cern.ch>; Piotr Traczyk - SINS Warsaw <ptraczyk@fuw.edu.pl>

Definition at line 43 of file DTMeantimerPatternReco.h.

Member Typedef Documentation

typedef std::pair<DTHitPairForFit*, DTEnums::DTCellSide> DTMeantimerPatternReco::AssPoint [private]

Definition at line 72 of file DTMeantimerPatternReco.h.

Constructor & Destructor Documentation

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

Constructor.

Definition at line 37 of file DTMeantimerPatternReco.cc.

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

                                                                          : 
DTRecSegment2DBaseAlgo(pset), theAlgoName("DTMeantimerPatternReco")
{

  theMaxAllowedHits = pset.getParameter<unsigned int>("MaxAllowedHits"); // 100
  theAlphaMaxTheta = pset.getParameter<double>("AlphaMaxTheta");// 0.1 ;
  theAlphaMaxPhi = pset.getParameter<double>("AlphaMaxPhi");// 1.0 ;
  theMaxT0 = pset.getParameter<double>("MaxT0");
  theMinT0 = pset.getParameter<double>("MinT0");
  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 53 of file DTMeantimerPatternReco.cc.

References theCleaner, and theUpdator.

                                                {
  delete theUpdator;
  delete theCleaner;
}

Member Function Documentation

void DTMeantimerPatternReco::addHits	(	const DTSuperLayer *	sl,
		std::vector< AssPoint > &	assHits,
		const std::vector< DTHitPairForFit * > &	hits,
		std::vector< DTSegmentCand * > &	result,
		std::vector< AssPoint > &	usedHits
	)		`[private]`

Referenced by CSCNeutronReader::addHits(), and buildSegments().

virtual std::string DTMeantimerPatternReco::algoName ( void ) const [inline, virtual]

return the algo name

Implements DTRecSegment2DBaseAlgo.

Definition at line 61 of file DTMeantimerPatternReco.h.

References theAlgoName.

{ return theAlgoName; }

vector< DTSegmentCand * > DTMeantimerPatternReco::buildSegments	(	const DTSuperLayer *	sl,
		const std::vector< DTHitPairForFit * > &	hits
	)		`[private]`

Definition at line 104 of file DTMeantimerPatternReco.cc.

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

Referenced by reconstruct().

                                                                                  {

  typedef vector<DTHitPairForFit*> hitCont;
  typedef hitCont::const_iterator  hitIter;
  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;
  }

  // 10-Mar-2004 SL
  // put a protection against heavily populated chambers, for which the segment
  // building could lead to infinite memory usage...
  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;
  
  vector<AssPoint> usedHits;

  // 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<AssPoint> assHits;
          // create a candidate hit list
          assHits.push_back(AssPoint(*firstHit,codes[firstLR]));
          assHits.push_back(AssPoint(*lastHit,codes[lastLR]));

          // run hit adding/segment building 
          maxfound = 3;
          addHits(sl,assHits,hitsForFit,result, usedHits);
        }
      }
    }
  }

  // now I have a couple of segment hypotheses, should check for ghost
  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 433 of file DTCombinatorialExtendedPatternReco.cc.

Referenced by CSCNeutronReader::addHits().

                                                                      {
  for (vector<DTSegmentCand*>::iterator cand=cands.begin();
       cand!=cands.end(); ++cand) 
    if (*(*cand)==*seg) return false;
  return true;
}

bool DTMeantimerPatternReco::fitWithT0	(	const std::vector< AssPoint > &	assHits,
		double &	chi2,
		double &	t0_corr,
		const bool	fitdebug
	)		`[private]`

Definition at line 336 of file DTMeantimerPatternReco.cc.

References a, b, gather_cfg::cout, debug, delta, DTEnums::Left, asciidump::s, theMaxChi2, theMinT0, ExpressReco_HICollisions_FallBack::x, and ExpressReco_HICollisions_FallBack::y.

Referenced by CSCNeutronReader::addHits().

{
  typedef vector < pair<double,double> > hitCoord;
  double a,b,coordError,x,y;
  double sx=0,sy=0,sxy=0,sxx=0,ssx=0,ssy=0,s=0,ss=0;
  int leftHits=0,rightHits=0;

  if (assHits.size()<3) return false;

  // I'm assuming the single hit error is the same for all hits...
  coordError=((*(assHits.begin())).first)->localPositionError().xx();

  for (vector<AssPoint>::const_iterator assHit=assHits.begin(); assHit!=assHits.end(); ++assHit) {
    if (coordError!=((*(assHits.begin())).first)->localPositionError().xx()) 
      cout << "   DTMeantimerPatternReco: Warning! Hit errors different within one segment!" << endl;

    x=((*assHit).first)->localPosition((*assHit).second).z();
    y=((*assHit).first)->localPosition((*assHit).second).x();

    sx+=x;
    sy+=y;
    sxy+=x*y;
    sxx+=x*x;
    s++;
    
    if ((*assHit).second==DTEnums::Left) {
      leftHits++;
      ssx+=x;
      ssy+=y;
      ss++;
    } else {
      rightHits++;
      ssx-=x;
      ssy-=y;
      ss--;
    }  
  }      
                    
  if (fitdebug && debug)
    cout  << "   DTMeantimerPatternReco::fitWithT0   Left hits: " << leftHits << "  Right hits: " << rightHits << endl;

  if (leftHits && rightHits) {

    double delta = ss*ss*sxx+s*sx*sx+s*ssx*ssx-s*s*sxx-2*ss*sx*ssx;
    t0_corr=0.;

    if (delta) {
      a=(ssy*s*ssx+sxy*ss*ss+sy*sx*s-sy*ss*ssx-ssy*sx*ss-sxy*s*s)/delta;
      b=(ssx*sy*ssx+sxx*ssy*ss+sx*sxy*s-sxx*sy*s-ssx*sxy*ss-sx*ssy*ssx)/delta;
      t0_corr=(ssx*s*sxy+sxx*ss*sy+sx*sx*ssy-sxx*s*ssy-sx*ss*sxy-ssx*sx*sy)/delta;
    } else return false;
  } else {
    double d = s*sxx - sx*sx;
    b = (sxx*sy- sx*sxy)/ d;
    a = (s*sxy - sx*sy) / d;
    t0_corr=0;
  }

// Calculate the chi^2 of the hits AFTER t0 correction

  double chi,chi2_not0;  
  chi2=0;
  chi2_not0=0;
  
  for (vector<AssPoint>::const_iterator assHit=assHits.begin(); assHit!=assHits.end(); ++assHit) {
    x=((*assHit).first)->localPosition((*assHit).second).z();
    y=((*assHit).first)->localPosition((*assHit).second).x();

    chi=y-a*x-b;
    chi2_not0+=chi*chi/coordError;

    if ((*assHit).second==DTEnums::Left) chi-=t0_corr;
      else chi+=t0_corr;
      
    chi2+=chi*chi/coordError;
  }

  // For 3-hit segments ignore timing information
  if (assHits.size()<4) {
    chi2=chi2_not0;
//    if (chi2<theMaxChi2) return true;
    if (chi2<200.) return true;
      else return false;
  }

  t0_corr/=-0.00543; // convert drift distance to time

  if (debug && fitdebug) 
  {
    cout << "      t0_corr = " << t0_corr << "ns  chi2/nHits = " << chi2 << "/" << assHits.size() << endl;
    
    if (((chi2/(assHits.size()-2)<theMaxChi2) && (t0_corr<theMaxT0) && (t0_corr>theMinT0)) || (assHits.size()==4)) {
      cout << "      a  = " << a  << "       b  = " << b  << endl;
      for (vector<AssPoint>::const_iterator assHit=assHits.begin(); assHit!=assHits.end(); ++assHit) {
        x=((*assHit).first)->localPosition((*assHit).second).z();
        y=((*assHit).first)->localPosition((*assHit).second).x();
        cout << "   z= " << x << "   x= " << y;
        if ((*assHit).second==DTEnums::Left) cout << "   x_corr= " << y+t0_corr*0.00543;
                                       else  cout << "   x_corr= " << y-t0_corr*0.00543;
        cout << "   seg= " << a*x+b << endl;
      }
    }
  }

  if ((chi2/(assHits.size()-2)<theMaxChi2) && (t0_corr<theMaxT0) && (t0_corr>theMinT0)) return true;
    else return false;
}

bool DTMeantimerPatternReco::geometryFilter	(	const DTWireId	first,
		const DTWireId	second
	)		const `[private]`

Definition at line 300 of file DTMeantimerPatternReco.cc.

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

Referenced by CSCNeutronReader::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< DTHitPairForFit * > DTMeantimerPatternReco::initHits	(	const DTSuperLayer *	sl,
		const std::vector< DTRecHit1DPair > &	hits
	)		`[private]`

Definition at line 163 of file DTClusterer.cc.

References DTSuperLayer::layer(), DTLayerId::layer(), query::result, python::multivaluedict::sort(), DTLayer::specificTopology(), GeomDet::toGlobal(), GeomDet::toLocal(), DTWireId::wire(), DTTopology::wirePosition(), and PV3DBase< T, PVType, FrameType >::x().

Referenced by reconstruct().

                                                                                          {
  vector<pair<float, DTRecHit1DPair> > result;
  for (vector<DTRecHit1DPair>::const_iterator pair=pairs.begin();
       pair!=pairs.end(); ++pair ) {

    // get wire
    DTWireId wid = (*pair).wireId();
    // get Layer
    const DTLayer* lay= sl->layer(wid.layer());
    // get wire position in SL (only x)
    LocalPoint
      posInLayer(lay->specificTopology().wirePosition(wid.wire()),0.,0.);
    LocalPoint posInSL = sl->toLocal(lay->toGlobal(posInLayer));
    // put the pair into result
    result.push_back(make_pair(posInSL.x(), *pair));
  }
  // sorted by x
  sort(result.begin(), result.end(), sortClusterByX());

  return result;
}

void DTMeantimerPatternReco::rawFit	(	double &	a,
		double &	b,
		const std::vector< std::pair< double, double > > &	hits
	)		`[private]`

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 60 of file DTMeantimerPatternReco.cc.

References buildSegments(), 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;
  vector<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);

    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 85 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);
}