CSCSegAlgoST Class Reference

#include <CSCSegAlgoST.h>

Inheritance diagram for CSCSegAlgoST:

Public Types
typedef std::deque< bool >	BoolContainer
typedef std::vector< const CSCRecHit2D * >	ChamberHitContainer
	Typedefs. More...
typedef std::vector< std::vector< const CSCRecHit2D * > >	Segments

Public Member Functions
std::vector< CSCSegment >	buildSegments (const ChamberHitContainer &rechits)
std::vector< CSCSegment >	buildSegments2 (const ChamberHitContainer &rechits)
std::vector< std::vector< const CSCRecHit2D * > >	chainHits (const CSCChamber *aChamber, const ChamberHitContainer &rechits)
std::vector< std::vector< const CSCRecHit2D * > >	clusterHits (const CSCChamber *aChamber, const ChamberHitContainer &rechits)
	CSCSegAlgoST (const edm::ParameterSet &ps)
	Constructor. More...
std::vector< CSCSegment >	prune_bad_hits (const CSCChamber *aChamber, std::vector< CSCSegment > &segments)
std::vector< CSCSegment >	run (const CSCChamber *aChamber, const ChamberHitContainer &rechits) override
	~CSCSegAlgoST () override
	Destructor. More...
Public Member Functions inherited from CSCSegmentAlgorithm
	CSCSegmentAlgorithm (const edm::ParameterSet &)
	Constructor. More...
virtual std::vector< CSCSegment >	run (const CSCChamber *chamber, const std::vector< const CSCRecHit2D * > &rechits)=0
virtual	~CSCSegmentAlgorithm ()
	Destructor. More...

Private Member Functions
bool	adjustCovariance (void)
const CSCChamber *	chamber () const
void	ChooseSegments (void)
void	ChooseSegments2 (int best_seg)
void	ChooseSegments2a (std::vector< ChamberHitContainer > &best_segments, int best_seg)
void	ChooseSegments3 (int best_seg)
void	ChooseSegments3 (std::vector< ChamberHitContainer > &best_segments, std::vector< float > &best_weight, int best_seg)
void	dumpSegment (const CSCSegment &seg) const
void	findDuplicates (std::vector< CSCSegment > &segments)
bool	isGoodToMerge (bool isME11a, ChamberHitContainer &newChain, ChamberHitContainer &oldChain)
const edm::ParameterSet &	pset (void) const
double	theWeight (double coordinate_1, double coordinate_2, double coordinate_3, float layer_1, float layer_2, float layer_3)
	Utility functions. More...

Private Attributes
float	a_yweightPenaltyThreshold [5][5]
bool	adjustCovariance_
double	BPMinImprovement
bool	BrutePruning
double	chi2Norm_3D_
std::vector< ChamberHitContainer >	chosen_Psegments
std::vector< float >	chosen_weight_A
bool	condpass1
	Flag whether to 'improve' covariance matrix. More...
bool	condpass2
std::vector< float >	curv_A
std::vector< float >	curv_noL1_A
std::vector< float >	curv_noL2_A
std::vector< float >	curv_noL3_A
std::vector< float >	curv_noL4_A
std::vector< float >	curv_noL5_A
std::vector< float >	curv_noL6_A
double	curvePenalty
double	curvePenaltyThreshold
bool	debug
double	dXclusBoxMax
double	dYclusBoxMax
Segments	GoodSegments
double	hitDropLimit4Hits
double	hitDropLimit5Hits
double	hitDropLimit6Hits
int	maxRecHitsInCluster
int	minHitsPerSegment
const std::string	myName_
bool	onlyBestSegment
ChamberHitContainer	PAhits_onLayer [6]
bool	preClustering
bool	preClustering_useChaining
bool	prePrun_
	Chi^2 normalization for the initial fit. More...
double	prePrunLimit_
ChamberHitContainer	protoSegment
bool	Pruning
const edm::ParameterSet	ps_
std::vector< ChamberHitContainer >	Psegments
ChamberHitContainer	Psegments_hits
std::vector< ChamberHitContainer >	Psegments_noL1
std::vector< ChamberHitContainer >	Psegments_noL2
std::vector< ChamberHitContainer >	Psegments_noL3
std::vector< ChamberHitContainer >	Psegments_noL4
std::vector< ChamberHitContainer >	Psegments_noL5
std::vector< ChamberHitContainer >	Psegments_noL6
std::vector< ChamberHitContainer >	Psegments_noLx
CSCSegAlgoShowering *	showering_
const CSCChamber *	theChamber
bool	useShowering
std::vector< float >	weight_A
std::vector< float >	weight_B
std::vector< float >	weight_noL1_A
std::vector< float >	weight_noL1_B
std::vector< float >	weight_noL2_A
std::vector< float >	weight_noL2_B
std::vector< float >	weight_noL3_A
std::vector< float >	weight_noL3_B
std::vector< float >	weight_noL4_A
std::vector< float >	weight_noL4_B
std::vector< float >	weight_noL5_A
std::vector< float >	weight_noL5_B
std::vector< float >	weight_noL6_A
std::vector< float >	weight_noL6_B
std::vector< float >	weight_noLx_A
double	yweightPenalty
double	yweightPenaltyThreshold

Detailed Description

This algorithm is based on the Minimum Spanning Tree (ST) approach for building endcap muon track segments out of the rechit's in a CSCChamber.

A CSCSegment is a RecSegment4D, and is built from CSCRecHit2D objects, each of which is a RecHit2DLocalPos.

This builds segments consisting of at least 3 hits. Segments can share a common rechit, but only one.

Authors

S. Stoynev - NWU I. Bloch - FNAL E. James - FNAL A. Sakharov - WSU (extensive revision to handle weird segments) ... ... ... T. Cox - UC Davis (struggling to handle this monster)

Definition at line 32 of file CSCSegAlgoST.h.

Member Typedef Documentation

typedef std::deque<bool> CSCSegAlgoST::BoolContainer

Definition at line 41 of file CSCSegAlgoST.h.

typedef std::vector<const CSCRecHit2D*> CSCSegAlgoST::ChamberHitContainer

Typedefs.

Definition at line 39 of file CSCSegAlgoST.h.

typedef std::vector< std::vector<const CSCRecHit2D* > > CSCSegAlgoST::Segments

Definition at line 40 of file CSCSegAlgoST.h.

Constructor & Destructor Documentation

CSCSegAlgoST::CSCSegAlgoST ( const edm::ParameterSet &  ps )

explicit

Constructor.

CSCSegAlgoST.cc

Authors

: S. Stoynev - NU I. Bloch - FNAL E. James - FNAL A. Sakharov - WSU T. Cox - UC Davis - segment fit factored out of entangled code - Jan 2015

Improve the covariance matrix?

Definition at line 33 of file CSCSegAlgoST.cc.

References adjustCovariance_, BPMinImprovement, BrutePruning, chi2Norm_3D_, curvePenalty, curvePenaltyThreshold, debug, dXclusBoxMax, dYclusBoxMax, edm::ParameterSet::getParameter(), edm::ParameterSet::getUntrackedParameter(), hitDropLimit4Hits, hitDropLimit5Hits, hitDropLimit6Hits, maxRecHitsInCluster, minHitsPerSegment, onlyBestSegment, preClustering, preClustering_useChaining, prePrun_, prePrunLimit_, Pruning, showering_, useShowering, yweightPenalty, and yweightPenaltyThreshold.

                                                    : 
  CSCSegmentAlgorithm(ps), myName_("CSCSegAlgoST"), ps_(ps),  showering_(nullptr) {
    
  debug                  = ps.getUntrackedParameter<bool>("CSCDebug");
  //  minLayersApart         = ps.getParameter<int>("minLayersApart");
  //  nSigmaFromSegment      = ps.getParameter<double>("nSigmaFromSegment");
  minHitsPerSegment      = ps.getParameter<int>("minHitsPerSegment");
  //  muonsPerChamberMax     = ps.getParameter<int>("CSCSegmentPerChamberMax");      
  //  chi2Max                = ps.getParameter<double>("chi2Max");
  dXclusBoxMax           = ps.getParameter<double>("dXclusBoxMax");
  dYclusBoxMax           = ps.getParameter<double>("dYclusBoxMax");
  preClustering          = ps.getParameter<bool>("preClustering");
  preClustering_useChaining    = ps.getParameter<bool>("preClusteringUseChaining");
  Pruning                = ps.getParameter<bool>("Pruning");
  BrutePruning           = ps.getParameter<bool>("BrutePruning");
  BPMinImprovement        = ps.getParameter<double>("BPMinImprovement");
  // maxRecHitsInCluster is the maximal number of hits in a precluster that is being processed
  // This cut is intended to remove messy events. Currently nothing is returned if there are
  // more that maxRecHitsInCluster hits. It could be useful to return an estimate of the 
  // cluster position, which is available.
  maxRecHitsInCluster    = ps.getParameter<int>("maxRecHitsInCluster");
  onlyBestSegment        = ps.getParameter<bool>("onlyBestSegment");

  hitDropLimit4Hits      = ps.getParameter<double>("hitDropLimit4Hits");
  hitDropLimit5Hits      = ps.getParameter<double>("hitDropLimit5Hits");
  hitDropLimit6Hits      = ps.getParameter<double>("hitDropLimit6Hits");
  
  yweightPenaltyThreshold      = ps.getParameter<double>("yweightPenaltyThreshold");
  yweightPenalty               = ps.getParameter<double>("yweightPenalty");
                                             
  curvePenaltyThreshold        = ps.getParameter<double>("curvePenaltyThreshold");
  curvePenalty                 = ps.getParameter<double>("curvePenalty");

  useShowering = ps.getParameter<bool>("useShowering");
  if (useShowering) showering_   = new CSCSegAlgoShowering( ps );

  adjustCovariance_ = ps.getParameter<bool>("CorrectTheErrors");

  chi2Norm_3D_      = ps.getParameter<double>("NormChi2Cut3D");
  prePrun_          = ps.getParameter<bool>("prePrun");
  prePrunLimit_     = ps.getParameter<double>("prePrunLimit");

  if (debug) edm::LogVerbatim("CSCSegment") << "CSCSegAlgoST: with factored conditioned segment fit";
}

CSCSegAlgoST::~CSCSegAlgoST ( )

override

Destructor.

Definition at line 82 of file CSCSegAlgoST.cc.

References showering_.

                            {
  delete showering_;
}

Member Function Documentation

bool CSCSegAlgoST::adjustCovariance ( void  )

inlineprivate

Definition at line 89 of file CSCSegAlgoST.h.

References adjustCovariance_, ChooseSegments(), ChooseSegments2(), ChooseSegments2a(), ChooseSegments3(), dumpSegment(), findDuplicates(), isGoodToMerge(), and theWeight().

Referenced by buildSegments(), and prune_bad_hits().

{ return adjustCovariance_;}

std::vector< CSCSegment > CSCSegAlgoST::buildSegments

(

const ChamberHitContainer &

rechits

)

Build track segments in this chamber (this is where the actual segment-building algorithm hides.)

Definition at line 666 of file CSCSegAlgoST.cc.

References a_yweightPenaltyThreshold, adjustCovariance(), chamber(), CSCSegment::chi2(), CSCSegFit::chi2(), chi2Norm_3D_, ChooseSegments2a(), ChooseSegments3(), chosen_Psegments, chosen_weight_A, condpass1, condpass2, CSCSegFit::covarianceMatrix(), curv_A, curv_noL1_A, curv_noL2_A, curv_noL3_A, curv_noL4_A, curv_noL5_A, curv_noL6_A, curvePenalty, curvePenaltyThreshold, debug, dumpSegment(), relativeConstraints::empty, end, CSCCondSegFit::fit(), GoodSegments, hitDropLimit4Hits, hitDropLimit5Hits, hitDropLimit6Hits, mps_fire::i, triggerObjects_cff::id, createfilelist::int, CSCSegFit::intercept(), CSCSegFit::localdir(), LogTrace, maxRecHitsInCluster, minHitsPerSegment, CSCSegFit::ndof(), CSCSegFit::nhits(), CSCSegment::nRecHits(), onlyBestSegment, PAhits_onLayer, prePrun_, prePrunLimit_, protoSegment, Psegments, Psegments_hits, Psegments_noL1, Psegments_noL2, Psegments_noL3, Psegments_noL4, Psegments_noL5, Psegments_noL6, Psegments_noLx, pset(), CSCSegFit::scaleXError(), CSCSegFit::setScaleXError(), showering_, CSCSegAlgoShowering::showerSeg(), findQualityFiles::size, mathSSE::sqrt(), groupFilesInBlocks::temp, theChamber, theWeight(), useShowering, weight_A, weight_B, weight_noL1_A, weight_noL1_B, weight_noL2_A, weight_noL2_B, weight_noL3_A, weight_noL3_B, weight_noL4_A, weight_noL4_B, weight_noL5_A, weight_noL5_B, weight_noL6_A, weight_noL6_B, weight_noLx_A, CSCCondSegFit::worstHit(), x, y, and yweightPenalty.

Referenced by run().

                                                                                    {

  // Clear buffer for segment vector
  std::vector<CSCSegment> segmentInChamber;
  segmentInChamber.clear(); // list of final segments

  // CSC Ring;
  unsigned int thering    = 999;
  unsigned int thestation = 999;
  //unsigned int thecham    = 999;

  std::vector<int> hits_onLayerNumber(6);

  unsigned int UpperLimit = maxRecHitsInCluster;
  if (int(rechits.size()) < minHitsPerSegment) return segmentInChamber;
 
  for(int iarray = 0; iarray <6; ++iarray) { // magic number 6: number of layers in CSC chamber - not gonna change :)
    PAhits_onLayer[iarray].clear();
    hits_onLayerNumber[iarray] = 0;    
  }

  chosen_Psegments.clear();
  chosen_weight_A.clear();

  Psegments.clear();
  Psegments_noLx.clear();
  Psegments_noL1.clear();
  Psegments_noL2.clear();
  Psegments_noL3.clear();
  Psegments_noL4.clear();
  Psegments_noL5.clear();
  Psegments_noL6.clear();

  Psegments_hits.clear();
  
  weight_A.clear();
  weight_noLx_A.clear();
  weight_noL1_A.clear();
  weight_noL2_A.clear();
  weight_noL3_A.clear();
  weight_noL4_A.clear();
  weight_noL5_A.clear();
  weight_noL6_A.clear();

  weight_B.clear();
  weight_noL1_B.clear();
  weight_noL2_B.clear();
  weight_noL3_B.clear();
  weight_noL4_B.clear();
  weight_noL5_B.clear();
  weight_noL6_B.clear();

  curv_A.clear();
  curv_noL1_A.clear();
  curv_noL2_A.clear();
  curv_noL3_A.clear();
  curv_noL4_A.clear();
  curv_noL5_A.clear();
  curv_noL6_A.clear();

  // definition of middle layer for n-hit segment
  int midlayer_pointer[6] = {0,0,2,3,3,4};
  
  // int n_layers_missed_tot = 0;
  int n_layers_occupied_tot = 0;
  int n_layers_processed = 0;

  float min_weight_A = 99999.9;
  float min_weight_noLx_A = 99999.9;

  //float best_weight_B = -1.;
  //float best_weight_noLx_B = -1.;

  //float best_curv_A = -1.;
  //float best_curv_noLx_A = -1.;

  int best_pseg = -1;
  int best_noLx_pseg = -1;
  int best_Layer_noLx = -1;

  //************************************************************************;    
  //***   Start segment building   *****************************************;    
  //************************************************************************;    
  
  // Determine how many layers with hits we have
  // Fill all hits into the layer hit container:
  
  // Have 2 standard arrays: one giving the number of hits per layer. 
  // The other the corresponding hits. 
  
  // Loop all available hits, count hits per layer and fill the hits into array by layer
  for(size_t M = 0; M < rechits.size(); ++M) {
    // add hits to array per layer and count hits per layer:
    hits_onLayerNumber[ rechits[M]->cscDetId().layer()-1 ] += 1;
    if(hits_onLayerNumber[ rechits[M]->cscDetId().layer()-1 ] == 1 ) n_layers_occupied_tot += 1;
    // add hits to vector in array
    PAhits_onLayer[rechits[M]->cscDetId().layer()-1]    .push_back(rechits[M]);    
  }
 
  // We have now counted the hits per layer and filled pointers to the hits into an array
  
  int tothits = 0;
  int maxhits = 0;
  int nexthits = 0;
  int maxlayer = -1;
  int nextlayer = -1;

  for(size_t i = 0; i< hits_onLayerNumber.size(); ++i){
    //std::cout<<"We have "<<hits_onLayerNumber[i]<<" hits on layer "<<i+1<<std::endl;
    tothits += hits_onLayerNumber[i];
    if (hits_onLayerNumber[i] > maxhits) {
      nextlayer = maxlayer;
      nexthits = maxhits;
      maxlayer = i;
      maxhits = hits_onLayerNumber[i];
    }
    else if (hits_onLayerNumber[i] > nexthits) {
      nextlayer = i;
      nexthits = hits_onLayerNumber[i];
    }
  }


  if (tothits > (int)UpperLimit) {
    if (n_layers_occupied_tot > 4) {
      tothits = tothits - hits_onLayerNumber[maxlayer];
      n_layers_occupied_tot = n_layers_occupied_tot - 1;
      PAhits_onLayer[maxlayer].clear();
      hits_onLayerNumber[maxlayer] = 0;
    }
  }

  if (tothits > (int)UpperLimit) {
    if (n_layers_occupied_tot > 4) {
      tothits = tothits - hits_onLayerNumber[nextlayer];
      n_layers_occupied_tot = n_layers_occupied_tot - 1;
      PAhits_onLayer[nextlayer].clear();
      hits_onLayerNumber[nextlayer] = 0;
    }
  }

  if (tothits > (int)UpperLimit){ 

  //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  // Showering muon - returns nothing if chi2 == -1 (see comment in SegAlgoShowering)
  //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++  
    if (useShowering) {
      CSCSegment segShower = showering_->showerSeg(theChamber, rechits);
      if( debug ) dumpSegment( segShower );
      // Make sure have at least 3 hits...
      if ( segShower.nRecHits() < 3 ) return segmentInChamber;
      if ( segShower.chi2() == -1 ) return segmentInChamber;
      segmentInChamber.push_back(segShower);
      return segmentInChamber;  
    } 
    else {
    //  LogTrace("CSCSegment|CSC") << "[CSCSegAlgoST::buildSegments] No. of rechits in the cluster/chamber > " 
    //  << UpperLimit << " ... Segment finding in the cluster/chamber canceled!";
      CSCDetId id = rechits[0]->cscDetId();
      edm::LogVerbatim("CSCSegment|CSC") << "[CSCSegAlgoST::buildSegments] No. of rechits in the cluster/chamber > " 
                     << UpperLimit << " ... Segment finding canceled in CSC" << id;
      return segmentInChamber;  
    }
  }

  // Find out which station, ring and chamber we are in 
  // Used to choose station/ring dependant y-weight cuts

  if( !rechits.empty() ) {
    thering = rechits[0]->cscDetId().ring();
    thestation = rechits[0]->cscDetId().station();
    //thecham = rechits[0]->cscDetId().chamber();
  }

  // std::cout<<"We are in Station/ring/chamber: "<<thestation <<" "<< thering<<" "<< thecham<<std::endl;

  // Cut-off parameter - don't reconstruct segments with less than X hits
  if( n_layers_occupied_tot < minHitsPerSegment ) { 
    return segmentInChamber;
  }
  
  // Start building all possible hit combinations:

  // loop over the six chamber layers and form segment candidates from the available hits:

  for(int layer = 0; layer < 6; ++layer) {

    // *****************************************************************
    // *** Set missed layer counter here (not currently implemented) ***
    // *****************************************************************
    // if( PAhits_onLayer[layer].size() == 0 ) {
    //   n_layers_missed_tot += 1;
    // }

    if( !PAhits_onLayer[layer].empty() ) {
      n_layers_processed += 1;
    }

    // Save the size of the protosegment before hits were added on the current layer
    int orig_number_of_psegs = Psegments.size();
    int orig_number_of_noL1_psegs = Psegments_noL1.size();
    int orig_number_of_noL2_psegs = Psegments_noL2.size();
    int orig_number_of_noL3_psegs = Psegments_noL3.size();
    int orig_number_of_noL4_psegs = Psegments_noL4.size();
    int orig_number_of_noL5_psegs = Psegments_noL5.size();
    int orig_number_of_noL6_psegs = Psegments_noL6.size();

    // loop over the hits on the layer and initiate protosegments or add hits to protosegments
    for(int hit = 0; hit < int(PAhits_onLayer[layer].size()); ++hit) { // loop all hits on the Layer number "layer"

      // create protosegments from all hits on the first layer with hits
      if( orig_number_of_psegs == 0 ) { // would be faster to turn this around - ask for "orig_number_of_psegs != 0"

    Psegments_hits.push_back(PAhits_onLayer[layer][hit]);

    Psegments.push_back(Psegments_hits); 
    Psegments_noL6.push_back(Psegments_hits); 
    Psegments_noL5.push_back(Psegments_hits); 
    Psegments_noL4.push_back(Psegments_hits); 
    Psegments_noL3.push_back(Psegments_hits); 
    Psegments_noL2.push_back(Psegments_hits); 

    // Initialize weights corresponding to this segment for first hit (with 0)

    curv_A.push_back(0.0);
    curv_noL6_A.push_back(0.0); 
    curv_noL5_A.push_back(0.0); 
    curv_noL4_A.push_back(0.0); 
    curv_noL3_A.push_back(0.0); 
    curv_noL2_A.push_back(0.0); 

    weight_A.push_back(0.0);
    weight_noL6_A.push_back(0.0); 
    weight_noL5_A.push_back(0.0); 
    weight_noL4_A.push_back(0.0); 
    weight_noL3_A.push_back(0.0); 
    weight_noL2_A.push_back(0.0); 

    weight_B.push_back(0.0);
    weight_noL6_B.push_back(0.0); 
    weight_noL5_B.push_back(0.0); 
    weight_noL4_B.push_back(0.0); 
    weight_noL3_B.push_back(0.0); 
    weight_noL2_B.push_back(0.0); 
    
    // reset array for next hit on next layer
    Psegments_hits    .clear();
      }
      else {
    if( orig_number_of_noL1_psegs == 0 ) {

      Psegments_hits.push_back(PAhits_onLayer[layer][hit]);

      Psegments_noL1.push_back(Psegments_hits); 

      // Initialize weight corresponding to this segment for first hit (with 0)

      curv_noL1_A.push_back(0.0);

      weight_noL1_A.push_back(0.0);

      weight_noL1_B.push_back(0.0);
    
      // reset array for next hit on next layer
      Psegments_hits    .clear();

    }

    // loop over the protosegments and create a new protosegments for each hit-1 on this layer
    
        for( int pseg = 0; pseg < orig_number_of_psegs; ++pseg ) { 

      int pseg_pos = (pseg)+((hit)*orig_number_of_psegs);
      int pseg_noL1_pos = (pseg)+((hit)*orig_number_of_noL1_psegs);
      int pseg_noL2_pos = (pseg)+((hit)*orig_number_of_noL2_psegs);
      int pseg_noL3_pos = (pseg)+((hit)*orig_number_of_noL3_psegs);
      int pseg_noL4_pos = (pseg)+((hit)*orig_number_of_noL4_psegs);
      int pseg_noL5_pos = (pseg)+((hit)*orig_number_of_noL5_psegs);
      int pseg_noL6_pos = (pseg)+((hit)*orig_number_of_noL6_psegs);

      // - Loop all psegs. 
      // - If not last hit, clone  existing protosegments  (PAhits_onLayer[layer].size()-1) times
      // - then add the new hits

      if( ! (hit == int(PAhits_onLayer[layer].size()-1)) ) { // not the last hit - prepare (copy) new protosegments for the following hits
        // clone psegs (to add next hits or last hit on layer):

        Psegments.push_back(Psegments[ pseg_pos ]); 
        if (n_layers_processed != 2 && pseg < orig_number_of_noL1_psegs) Psegments_noL1.push_back(Psegments_noL1[ pseg_noL1_pos ]); 
        if (n_layers_processed != 2 && pseg < orig_number_of_noL2_psegs) Psegments_noL2.push_back(Psegments_noL2[ pseg_noL2_pos ]); 
        if (n_layers_processed != 3 && pseg < orig_number_of_noL3_psegs) Psegments_noL3.push_back(Psegments_noL3[ pseg_noL3_pos ]); 
        if (n_layers_processed != 4 && pseg < orig_number_of_noL4_psegs) Psegments_noL4.push_back(Psegments_noL4[ pseg_noL4_pos ]); 
        if (n_layers_processed != 5 && pseg < orig_number_of_noL5_psegs) Psegments_noL5.push_back(Psegments_noL5[ pseg_noL5_pos ]); 
        if (n_layers_processed != 6 && pseg < orig_number_of_noL6_psegs) Psegments_noL6.push_back(Psegments_noL6[ pseg_noL6_pos ]); 
        // clone weight corresponding to this segment too
        weight_A.push_back(weight_A[ pseg_pos ]);
        if (n_layers_processed != 2 && pseg < orig_number_of_noL1_psegs) weight_noL1_A.push_back(weight_noL1_A[ pseg_noL1_pos ]);
        if (n_layers_processed != 2 && pseg < orig_number_of_noL2_psegs) weight_noL2_A.push_back(weight_noL2_A[ pseg_noL2_pos ]);
        if (n_layers_processed != 3 && pseg < orig_number_of_noL3_psegs) weight_noL3_A.push_back(weight_noL3_A[ pseg_noL3_pos ]);
        if (n_layers_processed != 4 && pseg < orig_number_of_noL4_psegs) weight_noL4_A.push_back(weight_noL4_A[ pseg_noL4_pos ]);
        if (n_layers_processed != 5 && pseg < orig_number_of_noL5_psegs) weight_noL5_A.push_back(weight_noL5_A[ pseg_noL5_pos ]);
        if (n_layers_processed != 6 && pseg < orig_number_of_noL6_psegs) weight_noL6_A.push_back(weight_noL6_A[ pseg_noL6_pos ]);
        // clone curvature variable corresponding to this segment too
        curv_A.push_back(curv_A[ pseg_pos ]);
        if (n_layers_processed != 2 && pseg < orig_number_of_noL1_psegs) curv_noL1_A.push_back(curv_noL1_A[ pseg_noL1_pos ]);
        if (n_layers_processed != 2 && pseg < orig_number_of_noL2_psegs) curv_noL2_A.push_back(curv_noL2_A[ pseg_noL2_pos ]);
        if (n_layers_processed != 3 && pseg < orig_number_of_noL3_psegs) curv_noL3_A.push_back(curv_noL3_A[ pseg_noL3_pos ]);
        if (n_layers_processed != 4 && pseg < orig_number_of_noL4_psegs) curv_noL4_A.push_back(curv_noL4_A[ pseg_noL4_pos ]);
        if (n_layers_processed != 5 && pseg < orig_number_of_noL5_psegs) curv_noL5_A.push_back(curv_noL5_A[ pseg_noL5_pos ]);
        if (n_layers_processed != 6 && pseg < orig_number_of_noL6_psegs) curv_noL6_A.push_back(curv_noL6_A[ pseg_noL6_pos ]);
        // clone "y"-weight corresponding to this segment too
        weight_B.push_back(weight_B[ pseg_pos ]);
        if (n_layers_processed != 2 && pseg < orig_number_of_noL1_psegs) weight_noL1_B.push_back(weight_noL1_B[ pseg_noL1_pos ]);
        if (n_layers_processed != 2 && pseg < orig_number_of_noL2_psegs) weight_noL2_B.push_back(weight_noL2_B[ pseg_noL2_pos ]);
        if (n_layers_processed != 3 && pseg < orig_number_of_noL3_psegs) weight_noL3_B.push_back(weight_noL3_B[ pseg_noL3_pos ]);
        if (n_layers_processed != 4 && pseg < orig_number_of_noL4_psegs) weight_noL4_B.push_back(weight_noL4_B[ pseg_noL4_pos ]);
        if (n_layers_processed != 5 && pseg < orig_number_of_noL5_psegs) weight_noL5_B.push_back(weight_noL5_B[ pseg_noL5_pos ]);
        if (n_layers_processed != 6 && pseg < orig_number_of_noL6_psegs) weight_noL6_B.push_back(weight_noL6_B[ pseg_noL6_pos ]);
      }
      // add hits to original pseg:
      Psegments[ pseg_pos ].push_back(PAhits_onLayer[ layer ][ hit ]);
      if (n_layers_processed != 2 && pseg < orig_number_of_noL1_psegs) Psegments_noL1[ pseg_noL1_pos ].push_back(PAhits_onLayer[ layer ][ hit ]);
      if (n_layers_processed != 2 && pseg < orig_number_of_noL2_psegs) Psegments_noL2[ pseg_noL2_pos ].push_back(PAhits_onLayer[ layer ][ hit ]);
      if (n_layers_processed != 3 && pseg < orig_number_of_noL3_psegs) Psegments_noL3[ pseg_noL3_pos ].push_back(PAhits_onLayer[ layer ][ hit ]);
      if (n_layers_processed != 4 && pseg < orig_number_of_noL4_psegs) Psegments_noL4[ pseg_noL4_pos ].push_back(PAhits_onLayer[ layer ][ hit ]);
      if (n_layers_processed != 5 && pseg < orig_number_of_noL5_psegs) Psegments_noL5[ pseg_noL5_pos ].push_back(PAhits_onLayer[ layer ][ hit ]);
      if (n_layers_processed != 6 && pseg < orig_number_of_noL6_psegs) Psegments_noL6[ pseg_noL6_pos ].push_back(PAhits_onLayer[ layer ][ hit ]);
            
      // calculate/update the weight (only for >2 hits on psegment):

      if( Psegments[ pseg_pos ].size() > 2 ) {
              
        // looks more exciting than it is. Here the weight is calculated. It is the difference in x of the last two and one but the last two hits, 
            // divided by the distance of the corresponding hits. Please refer to twiki page XXXX or CMS Note YYY (and use layer_distance)

        weight_A[ pseg_pos ] += theWeight(
                          (*(Psegments[ pseg_pos ].end()-1 ))->localPosition().x(), 
                          (*(Psegments[ pseg_pos ].end()-2 ))->localPosition().x(),
                          (*(Psegments[ pseg_pos ].end()-3 ))->localPosition().x(),
                          float((*(Psegments[ pseg_pos ].end()-1))->cscDetId().layer()),
                          float((*(Psegments[ pseg_pos ].end()-2))->cscDetId().layer()),
                          float((*(Psegments[ pseg_pos ].end()-3))->cscDetId().layer())
                          );

        weight_B[ pseg_pos ] += theWeight(
                          (*(Psegments[ pseg_pos ].end()-1 ))->localPosition().y(), 
                          (*(Psegments[ pseg_pos ].end()-2 ))->localPosition().y(),
                          (*(Psegments[ pseg_pos ].end()-3 ))->localPosition().y(),
                          float((*(Psegments[ pseg_pos ].end()-1))->cscDetId().layer()),
                          float((*(Psegments[ pseg_pos ].end()-2))->cscDetId().layer()),
                          float((*(Psegments[ pseg_pos ].end()-3))->cscDetId().layer())
                          );

        // if we have picked up the last hit go looking for pseg with the lowest (and second lowest?) weight

        if(int(Psegments[ pseg_pos ].size()) == n_layers_occupied_tot) {

          curv_A[ pseg_pos ] += theWeight(
                          (*(Psegments[ pseg_pos ].end()-1 ))->localPosition().x(), 
                          (*(Psegments[ pseg_pos ].end()-midlayer_pointer[n_layers_occupied_tot-1] ))->localPosition().x(),
                          (*(Psegments[ pseg_pos ].end()-n_layers_occupied_tot ))->localPosition().x(),
                          float((*(Psegments[ pseg_pos ].end()-1))->cscDetId().layer()),
                          float((*(Psegments[ pseg_pos ].end()-midlayer_pointer[n_layers_occupied_tot-1] ))->cscDetId().layer()),
                          float((*(Psegments[ pseg_pos ].end()-n_layers_occupied_tot ))->cscDetId().layer())
                          );

              if (curv_A[ pseg_pos ] > curvePenaltyThreshold) weight_A[ pseg_pos ] = weight_A[ pseg_pos ] * curvePenalty;

          if (weight_B[ pseg_pos ] > a_yweightPenaltyThreshold[thestation][thering]) weight_A[ pseg_pos ] = weight_A[ pseg_pos ] * yweightPenalty;
               
              if (weight_A[ pseg_pos ] < min_weight_A ) {
            min_weight_A = weight_A[ pseg_pos ];
                //best_weight_B = weight_B[ pseg_pos ];
                //best_curv_A = curv_A[ pseg_pos ];
            best_pseg = pseg_pos ;
              }

        }

        // alternative: fill map with weight and pseg (which is already ordered)? Seems a very good tool to go looking for segments from. 
            // As I understand, the segments would be inserted according to their weight, so the list would "automatically" be sorted.

      }

          if ( n_layers_occupied_tot > 3 ) {
        if (pseg < orig_number_of_noL1_psegs && (n_layers_processed != 2)) {
          if(( Psegments_noL1[ pseg_noL1_pos ].size() > 2 ) ) {
        
        // looks more exciting than it is. Here the weight is calculated. It is the difference in x of the last two and one but the last two hits, 
        // divided by the distance of the corresponding hits. Please refer to twiki page XXXX or CMS Note YYY (and use layer_distance)
        
        weight_noL1_A[ pseg_noL1_pos ] += theWeight(
                                (*(Psegments_noL1[ pseg_noL1_pos ].end()-1 ))->localPosition().x(), 
                                (*(Psegments_noL1[ pseg_noL1_pos ].end()-2 ))->localPosition().x(),
                                (*(Psegments_noL1[ pseg_noL1_pos ].end()-3 ))->localPosition().x(),
                                float((*(Psegments_noL1[ pseg_noL1_pos ].end()-1))->cscDetId().layer()),
                                float((*(Psegments_noL1[ pseg_noL1_pos ].end()-2))->cscDetId().layer()),
                                float((*(Psegments_noL1[ pseg_noL1_pos ].end()-3))->cscDetId().layer())
                                );

        weight_noL1_B[ pseg_noL1_pos ] += theWeight(
                                (*(Psegments_noL1[ pseg_noL1_pos ].end()-1 ))->localPosition().y(), 
                                (*(Psegments_noL1[ pseg_noL1_pos ].end()-2 ))->localPosition().y(),
                                (*(Psegments_noL1[ pseg_noL1_pos ].end()-3 ))->localPosition().y(),
                                float((*(Psegments_noL1[ pseg_noL1_pos ].end()-1))->cscDetId().layer()),
                                float((*(Psegments_noL1[ pseg_noL1_pos ].end()-2))->cscDetId().layer()),
                                float((*(Psegments_noL1[ pseg_noL1_pos ].end()-3))->cscDetId().layer())
                                );

        //if we have picked up the last hit go looking for pseg with the lowest (and second lowest?) weight

        if(int(Psegments_noL1[ pseg_noL1_pos ].size()) == n_layers_occupied_tot -1 ) {

          curv_noL1_A[ pseg_noL1_pos ] += theWeight(
                                (*(Psegments_noL1[ pseg_noL1_pos ].end()-1 ))->localPosition().x(), 
                                (*(Psegments_noL1[ pseg_noL1_pos ].end()-midlayer_pointer[n_layers_occupied_tot-2] ))->localPosition().x(),
                                (*(Psegments_noL1[ pseg_noL1_pos ].end()-(n_layers_occupied_tot-1) ))->localPosition().x(),
                                float((*(Psegments_noL1[ pseg_noL1_pos ].end()-1 ))->cscDetId().layer()),
                                float((*(Psegments_noL1[ pseg_noL1_pos ].end()-midlayer_pointer[n_layers_occupied_tot-2] ))->cscDetId().layer()),
                                float((*(Psegments_noL1[ pseg_noL1_pos ].end()-(n_layers_occupied_tot-1) ))->cscDetId().layer())
                                );

          if (curv_noL1_A[ pseg_noL1_pos ] > curvePenaltyThreshold) weight_noL1_A[ pseg_noL1_pos ] = weight_noL1_A[ pseg_noL1_pos ] * curvePenalty;

          if (weight_noL1_B[ pseg_noL1_pos ] > a_yweightPenaltyThreshold[thestation][thering]) 
            weight_noL1_A[ pseg_noL1_pos ] = weight_noL1_A[ pseg_noL1_pos ] * yweightPenalty;

          if (weight_noL1_A[ pseg_noL1_pos ] < min_weight_noLx_A ) {
            min_weight_noLx_A = weight_noL1_A[ pseg_noL1_pos ];
            //best_weight_noLx_B = weight_noL1_B[ pseg_noL1_pos ];
            //best_curv_noLx_A = curv_noL1_A[ pseg_noL1_pos ];
            best_noLx_pseg = pseg_noL1_pos;
                    best_Layer_noLx = 1;
          }

        }

        // alternative: fill map with weight and pseg (which is already ordered)? Seems a very good tool to go looking for segments from. 
        // As I understand, the segments would be inserted according to their weight, so the list would "automatically" be sorted.
        
          }
        }
      }

          if ( n_layers_occupied_tot > 3 ) {
        if (pseg < orig_number_of_noL2_psegs && ( n_layers_processed != 2 )) {
          if(( Psegments_noL2[ pseg_noL2_pos ].size() > 2 )) {
              
        // looks more exciting than it is. Here the weight is calculated. It is the difference in x of the last two and one but the last two hits, 
        // divided by the distance of the corresponding hits. Please refer to twiki page XXXX or CMS Note YYY (and use layer_distance)

        weight_noL2_A[ pseg_noL2_pos ] += theWeight(
                                (*(Psegments_noL2[ pseg_noL2_pos ].end()-1 ))->localPosition().x(), 
                                (*(Psegments_noL2[ pseg_noL2_pos ].end()-2 ))->localPosition().x(),
                                (*(Psegments_noL2[ pseg_noL2_pos ].end()-3 ))->localPosition().x(),
                                float((*(Psegments_noL2[ pseg_noL2_pos ].end()-1))->cscDetId().layer()),
                                float((*(Psegments_noL2[ pseg_noL2_pos ].end()-2))->cscDetId().layer()),
                                float((*(Psegments_noL2[ pseg_noL2_pos ].end()-3))->cscDetId().layer())
                                );

        weight_noL2_B[ pseg_noL2_pos ] += theWeight(
                                (*(Psegments_noL2[ pseg_noL2_pos ].end()-1 ))->localPosition().y(), 
                                (*(Psegments_noL2[ pseg_noL2_pos ].end()-2 ))->localPosition().y(),
                                (*(Psegments_noL2[ pseg_noL2_pos ].end()-3 ))->localPosition().y(),
                                float((*(Psegments_noL2[ pseg_noL2_pos ].end()-1))->cscDetId().layer()),
                                float((*(Psegments_noL2[ pseg_noL2_pos ].end()-2))->cscDetId().layer()),
                                float((*(Psegments_noL2[ pseg_noL2_pos ].end()-3))->cscDetId().layer())
                                );

        //if we have picked up the last hit go looking for pseg with the lowest (and second lowest?) weight

        if(int(Psegments_noL2[ pseg_noL2_pos ].size()) == n_layers_occupied_tot -1 ) {

          curv_noL2_A[ pseg_noL2_pos ] += theWeight(
                                (*(Psegments_noL2[ pseg_noL2_pos ].end()-1 ))->localPosition().x(), 
                                (*(Psegments_noL2[ pseg_noL2_pos ].end()-midlayer_pointer[n_layers_occupied_tot-2] ))->localPosition().x(),
                                (*(Psegments_noL2[ pseg_noL2_pos ].end()-(n_layers_occupied_tot-1) ))->localPosition().x(),
                                float((*(Psegments_noL2[ pseg_noL2_pos ].end()-1 ))->cscDetId().layer()),
                                float((*(Psegments_noL2[ pseg_noL2_pos ].end()-midlayer_pointer[n_layers_occupied_tot-2] ))->cscDetId().layer()),
                                float((*(Psegments_noL2[ pseg_noL2_pos ].end()-(n_layers_occupied_tot-1) ))->cscDetId().layer())
                                );

          if (curv_noL2_A[ pseg_noL2_pos ] > curvePenaltyThreshold) weight_noL2_A[ pseg_noL2_pos ] = weight_noL2_A[ pseg_noL2_pos ] * curvePenalty;

          if (weight_noL2_B[ pseg_noL2_pos ] > a_yweightPenaltyThreshold[thestation][thering]) 
            weight_noL2_A[ pseg_noL2_pos ] = weight_noL2_A[ pseg_noL2_pos ] * yweightPenalty;

          if (weight_noL2_A[ pseg_noL2_pos ] < min_weight_noLx_A ) {
            min_weight_noLx_A = weight_noL2_A[ pseg_noL2_pos ];
            //best_weight_noLx_B = weight_noL2_B[ pseg_noL2_pos ];
            //best_curv_noLx_A = curv_noL2_A[ pseg_noL2_pos ];
            best_noLx_pseg = pseg_noL2_pos;
                    best_Layer_noLx = 2;
          }

        }

        // alternative: fill map with weight and pseg (which is already ordered)? Seems a very good tool to go looking for segments from. 
        // As I understand, the segments would be inserted according to their weight, so the list would "automatically" be sorted.

          }
        }
      }

          if ( n_layers_occupied_tot > 3 ) {
        if (pseg < orig_number_of_noL3_psegs && ( n_layers_processed != 3 )) {
          if(( Psegments_noL3[ pseg_noL3_pos ].size() > 2 )) {
              
        // looks more exciting than it is. Here the weight is calculated. It is the difference in x of the last two and one but the last two hits, 
        // divided by the distance of the corresponding hits. Please refer to twiki page XXXX or CMS Note YYY (and use layer_distance)

        weight_noL3_A[ pseg_noL3_pos ] += theWeight(
                                (*(Psegments_noL3[ pseg_noL3_pos ].end()-1 ))->localPosition().x(), 
                                (*(Psegments_noL3[ pseg_noL3_pos ].end()-2 ))->localPosition().x(),
                                (*(Psegments_noL3[ pseg_noL3_pos ].end()-3 ))->localPosition().x(),
                                float((*(Psegments_noL3[ pseg_noL3_pos ].end()-1))->cscDetId().layer()),
                                float((*(Psegments_noL3[ pseg_noL3_pos ].end()-2))->cscDetId().layer()),
                                float((*(Psegments_noL3[ pseg_noL3_pos ].end()-3))->cscDetId().layer())
                                );

        weight_noL3_B[ pseg_noL3_pos ] += theWeight(
                                (*(Psegments_noL3[ pseg_noL3_pos ].end()-1 ))->localPosition().y(), 
                                (*(Psegments_noL3[ pseg_noL3_pos ].end()-2 ))->localPosition().y(),
                                (*(Psegments_noL3[ pseg_noL3_pos ].end()-3 ))->localPosition().y(),
                                float((*(Psegments_noL3[ pseg_noL3_pos ].end()-1))->cscDetId().layer()),
                                float((*(Psegments_noL3[ pseg_noL3_pos ].end()-2))->cscDetId().layer()),
                                float((*(Psegments_noL3[ pseg_noL3_pos ].end()-3))->cscDetId().layer())
                                );

        //if we have picked up the last hit go looking for pseg with the lowest (and second lowest?) weight

        if(int(Psegments_noL3[ pseg_noL3_pos ].size()) == n_layers_occupied_tot -1 ) {

          curv_noL3_A[ pseg_noL3_pos ] += theWeight(
                                (*(Psegments_noL3[ pseg_noL3_pos ].end()-1 ))->localPosition().x(), 
                                (*(Psegments_noL3[ pseg_noL3_pos ].end()-midlayer_pointer[n_layers_occupied_tot-2] ))->localPosition().x(),
                                (*(Psegments_noL3[ pseg_noL3_pos ].end()-(n_layers_occupied_tot-1) ))->localPosition().x(),
                                float((*(Psegments_noL3[ pseg_noL3_pos ].end()-1 ))->cscDetId().layer()),
                                float((*(Psegments_noL3[ pseg_noL3_pos ].end()-midlayer_pointer[n_layers_occupied_tot-2] ))->cscDetId().layer()),
                                float((*(Psegments_noL3[ pseg_noL3_pos ].end()-(n_layers_occupied_tot-1) ))->cscDetId().layer())
                                );

          if (curv_noL3_A[ pseg_noL3_pos ] > curvePenaltyThreshold) weight_noL3_A[ pseg_noL3_pos ] = weight_noL3_A[ pseg_noL3_pos ] * curvePenalty;

          if (weight_noL3_B[ pseg_noL3_pos ] > a_yweightPenaltyThreshold[thestation][thering]) 
            weight_noL3_A[ pseg_noL3_pos ] = weight_noL3_A[ pseg_noL3_pos ] * yweightPenalty;

          if (weight_noL3_A[ pseg_noL3_pos ] < min_weight_noLx_A ) {
            min_weight_noLx_A = weight_noL3_A[ pseg_noL3_pos ];
            //best_weight_noLx_B = weight_noL3_B[ pseg_noL3_pos ];
            //best_curv_noLx_A = curv_noL3_A[ pseg_noL3_pos ];
            best_noLx_pseg = pseg_noL3_pos;
                    best_Layer_noLx = 3;
          }

        }

        // alternative: fill map with weight and pseg (which is already ordered)? Seems a very good tool to go looking for segments from. 
        // As I understand, the segments would be inserted according to their weight, so the list would "automatically" be sorted.

          }
        }
      }

          if ( n_layers_occupied_tot > 3 ) {
        if (pseg < orig_number_of_noL4_psegs && ( n_layers_processed != 4 )) {
          if(( Psegments_noL4[ pseg_noL4_pos ].size() > 2 )) {
              
        // looks more exciting than it is. Here the weight is calculated. It is the difference in x of the last two and one but the last two hits, 
        // divided by the distance of the corresponding hits. Please refer to twiki page XXXX or CMS Note YYY (and use layer_distance)
          
        weight_noL4_A[ pseg_noL4_pos ] += theWeight(
                                (*(Psegments_noL4[ pseg_noL4_pos ].end()-1 ))->localPosition().x(), 
                                (*(Psegments_noL4[ pseg_noL4_pos ].end()-2 ))->localPosition().x(),
                                (*(Psegments_noL4[ pseg_noL4_pos ].end()-3 ))->localPosition().x(),
                                float((*(Psegments_noL4[ pseg_noL4_pos ].end()-1))->cscDetId().layer()),
                                float((*(Psegments_noL4[ pseg_noL4_pos ].end()-2))->cscDetId().layer()),
                                float((*(Psegments_noL4[ pseg_noL4_pos ].end()-3))->cscDetId().layer())
                                );

        weight_noL4_B[ pseg_noL4_pos ] += theWeight(
                                (*(Psegments_noL4[ pseg_noL4_pos ].end()-1 ))->localPosition().y(), 
                                (*(Psegments_noL4[ pseg_noL4_pos ].end()-2 ))->localPosition().y(),
                                (*(Psegments_noL4[ pseg_noL4_pos ].end()-3 ))->localPosition().y(),
                                float((*(Psegments_noL4[ pseg_noL4_pos ].end()-1))->cscDetId().layer()),
                                float((*(Psegments_noL4[ pseg_noL4_pos ].end()-2))->cscDetId().layer()),
                                float((*(Psegments_noL4[ pseg_noL4_pos ].end()-3))->cscDetId().layer())
                                );

        //if we have picked up the last hit go looking for pseg with the lowest (and second lowest?) weight

        if(int(Psegments_noL4[ pseg_noL4_pos ].size()) == n_layers_occupied_tot -1 ) {

          curv_noL4_A[ pseg_noL4_pos ] += theWeight(
                                (*(Psegments_noL4[ pseg_noL4_pos ].end()-1 ))->localPosition().x(), 
                                (*(Psegments_noL4[ pseg_noL4_pos ].end()-midlayer_pointer[n_layers_occupied_tot-2] ))->localPosition().x(),
                                (*(Psegments_noL4[ pseg_noL4_pos ].end()-(n_layers_occupied_tot-1) ))->localPosition().x(),
                                float((*(Psegments_noL4[ pseg_noL4_pos ].end()-1 ))->cscDetId().layer()),
                                float((*(Psegments_noL4[ pseg_noL4_pos ].end()-midlayer_pointer[n_layers_occupied_tot-2] ))->cscDetId().layer()),
                                float((*(Psegments_noL4[ pseg_noL4_pos ].end()-(n_layers_occupied_tot-1) ))->cscDetId().layer())
                                );

          if (curv_noL4_A[ pseg_noL4_pos ] > curvePenaltyThreshold) weight_noL4_A[ pseg_noL4_pos ] = weight_noL4_A[ pseg_noL4_pos ] * curvePenalty;

          if (weight_noL4_B[ pseg_noL4_pos ] > a_yweightPenaltyThreshold[thestation][thering]) 
            weight_noL4_A[ pseg_noL4_pos ] = weight_noL4_A[ pseg_noL4_pos ] * yweightPenalty;

          if (weight_noL4_A[ pseg_noL4_pos ] < min_weight_noLx_A ) {
            min_weight_noLx_A = weight_noL4_A[ pseg_noL4_pos ];
            //best_weight_noLx_B = weight_noL4_B[ pseg_noL4_pos ];
            //best_curv_noLx_A = curv_noL4_A[ pseg_noL4_pos ];
            best_noLx_pseg = pseg_noL4_pos;
                    best_Layer_noLx = 4;
          }

        }

        // alternative: fill map with weight and pseg (which is already ordered)? Seems a very good tool to go looking for segments from. 
        // As I understand, the segments would be inserted according to their weight, so the list would "automatically" be sorted.

          }
        }
      }

          if ( n_layers_occupied_tot > 4 ) {
        if (pseg < orig_number_of_noL5_psegs && ( n_layers_processed != 5 )) {
          if(( Psegments_noL5[ pseg_noL5_pos ].size() > 2 )){
              
        // looks more exciting than it is. Here the weight is calculated. It is the difference in x of the last two and one but the last two hits, 
        // divided by the distance of the corresponding hits. Please refer to twiki page XXXX or CMS Note YYY (and use layer_distance)

        weight_noL5_A[ pseg_noL5_pos ] += theWeight(
                                (*(Psegments_noL5[ pseg_noL5_pos ].end()-1 ))->localPosition().x(), 
                                (*(Psegments_noL5[ pseg_noL5_pos ].end()-2 ))->localPosition().x(),
                                (*(Psegments_noL5[ pseg_noL5_pos ].end()-3 ))->localPosition().x(),
                                float((*(Psegments_noL5[ pseg_noL5_pos ].end()-1))->cscDetId().layer()),
                                float((*(Psegments_noL5[ pseg_noL5_pos ].end()-2))->cscDetId().layer()),
                                float((*(Psegments_noL5[ pseg_noL5_pos ].end()-3))->cscDetId().layer())
                                );

        weight_noL5_B[ pseg_noL5_pos ] += theWeight(
                                (*(Psegments_noL5[ pseg_noL5_pos ].end()-1 ))->localPosition().y(), 
                                (*(Psegments_noL5[ pseg_noL5_pos ].end()-2 ))->localPosition().y(),
                                (*(Psegments_noL5[ pseg_noL5_pos ].end()-3 ))->localPosition().y(),
                                float((*(Psegments_noL5[ pseg_noL5_pos ].end()-1))->cscDetId().layer()),
                                float((*(Psegments_noL5[ pseg_noL5_pos ].end()-2))->cscDetId().layer()),
                                float((*(Psegments_noL5[ pseg_noL5_pos ].end()-3))->cscDetId().layer())
                                );

        //if we have picked up the last hit go looking for pseg with the lowest (and second lowest?) weight

        if(int(Psegments_noL5[ pseg_noL5_pos ].size()) == n_layers_occupied_tot -1 ) {

          curv_noL5_A[ pseg_noL5_pos ] += theWeight(
                                (*(Psegments_noL5[ pseg_noL5_pos ].end()-1 ))->localPosition().x(), 
                                (*(Psegments_noL5[ pseg_noL5_pos ].end()-midlayer_pointer[n_layers_occupied_tot-2] ))->localPosition().x(),
                                (*(Psegments_noL5[ pseg_noL5_pos ].end()-(n_layers_occupied_tot-1) ))->localPosition().x(),
                                float((*(Psegments_noL5[ pseg_noL5_pos ].end()-1 ))->cscDetId().layer()),
                                float((*(Psegments_noL5[ pseg_noL5_pos ].end()-midlayer_pointer[n_layers_occupied_tot-2] ))->cscDetId().layer()),
                                float((*(Psegments_noL5[ pseg_noL5_pos ].end()-(n_layers_occupied_tot-1) ))->cscDetId().layer())
                                );

          if (curv_noL5_A[ pseg_noL5_pos ] > curvePenaltyThreshold) weight_noL5_A[ pseg_noL5_pos ] = weight_noL5_A[ pseg_noL5_pos ] * curvePenalty;

          if (weight_noL5_B[ pseg_noL5_pos ] > a_yweightPenaltyThreshold[thestation][thering]) 
            weight_noL5_A[ pseg_noL5_pos ] = weight_noL5_A[ pseg_noL5_pos ] * yweightPenalty;

          if (weight_noL5_A[ pseg_noL5_pos ] < min_weight_noLx_A ) {
            min_weight_noLx_A = weight_noL5_A[ pseg_noL5_pos ];
            //best_weight_noLx_B = weight_noL5_B[ pseg_noL5_pos ];
            //best_curv_noLx_A = curv_noL5_A[ pseg_noL5_pos ];
            best_noLx_pseg = pseg_noL5_pos;
                    best_Layer_noLx = 5;
          }

        }

        // alternative: fill map with weight and pseg (which is already ordered)? Seems a very good tool to go looking for segments from. 
        // As I understand, the segments would be inserted according to their weight, so the list would "automatically" be sorted.

          }
        }
      }

          if ( n_layers_occupied_tot > 5 ) {
        if (pseg < orig_number_of_noL6_psegs && ( n_layers_processed != 6 )) {
          if(( Psegments_noL6[ pseg_noL6_pos ].size() > 2 )){
              
        // looks more exciting than it is. Here the weight is calculated. It is the difference in x of the last two and one but the last two hits, 
        // divided by the distance of the corresponding hits. Please refer to twiki page XXXX or CMS Note YYY (and use layer_distance)

        weight_noL6_A[ pseg_noL6_pos ] += theWeight(
                                (*(Psegments_noL6[ pseg_noL6_pos ].end()-1 ))->localPosition().x(), 
                                (*(Psegments_noL6[ pseg_noL6_pos ].end()-2 ))->localPosition().x(),
                                (*(Psegments_noL6[ pseg_noL6_pos ].end()-3 ))->localPosition().x(),
                                float((*(Psegments_noL6[ pseg_noL6_pos ].end()-1))->cscDetId().layer()),
                                float((*(Psegments_noL6[ pseg_noL6_pos ].end()-2))->cscDetId().layer()),
                                float((*(Psegments_noL6[ pseg_noL6_pos ].end()-3))->cscDetId().layer())
                                );

        weight_noL6_B[ pseg_noL6_pos ] += theWeight(
                                (*(Psegments_noL6[ pseg_noL6_pos ].end()-1 ))->localPosition().y(), 
                                (*(Psegments_noL6[ pseg_noL6_pos ].end()-2 ))->localPosition().y(),
                                (*(Psegments_noL6[ pseg_noL6_pos ].end()-3 ))->localPosition().y(),
                                float((*(Psegments_noL6[ pseg_noL6_pos ].end()-1))->cscDetId().layer()),
                                float((*(Psegments_noL6[ pseg_noL6_pos ].end()-2))->cscDetId().layer()),
                                float((*(Psegments_noL6[ pseg_noL6_pos ].end()-3))->cscDetId().layer())
                                );

        //if we have picked up the last hit go looking for pseg with the lowest (and second lowest?) weight

        if(int(Psegments_noL6[ pseg_noL6_pos ].size()) == n_layers_occupied_tot -1 ) {

          curv_noL6_A[ pseg_noL6_pos ] += theWeight(
                                (*(Psegments_noL6[ pseg_noL6_pos ].end()-1 ))->localPosition().x(), 
                                (*(Psegments_noL6[ pseg_noL6_pos ].end()-midlayer_pointer[n_layers_occupied_tot-2] ))->localPosition().x(),
                                (*(Psegments_noL6[ pseg_noL6_pos ].end()-(n_layers_occupied_tot-1) ))->localPosition().x(),
                                float((*(Psegments_noL6[ pseg_noL6_pos ].end()-1 ))->cscDetId().layer()),
                                float((*(Psegments_noL6[ pseg_noL6_pos ].end()-midlayer_pointer[n_layers_occupied_tot-2] ))->cscDetId().layer()),
                                float((*(Psegments_noL6[ pseg_noL6_pos ].end()-(n_layers_occupied_tot-1) ))->cscDetId().layer())
                                );

          if (curv_noL6_A[ pseg_noL6_pos ] > curvePenaltyThreshold) weight_noL6_A[ pseg_noL6_pos ] = weight_noL6_A[ pseg_noL6_pos ] * curvePenalty;

          if (weight_noL6_B[ pseg_noL6_pos ] > a_yweightPenaltyThreshold[thestation][thering]) 
            weight_noL6_A[ pseg_noL6_pos ] = weight_noL6_A[ pseg_noL6_pos ] * yweightPenalty;

          if (weight_noL6_A[ pseg_noL6_pos ] < min_weight_noLx_A ) {
            min_weight_noLx_A = weight_noL6_A[ pseg_noL6_pos ];
            //best_weight_noLx_B = weight_noL6_B[ pseg_noL6_pos ];
            //best_curv_noLx_A = curv_noL6_A[ pseg_noL6_pos ];
            best_noLx_pseg = pseg_noL6_pos;
                    best_Layer_noLx = 6;
          }

        }

        // alternative: fill map with weight and pseg (which is already ordered)? Seems a very good tool to go looking for segments from. 
        // As I understand, the segments would be inserted according to their weight, so the list would "automatically" be sorted.

          }
        }
      }

    }
      }
    }
  }

  //************************************************************************;    
  //***   End segment building   *******************************************;    
  //************************************************************************;    

  // Important part! Here segment(s) are actually chosen. All the good segments
  // could be chosen or some (best) ones only (in order to save time).

  // Check if there is a segment with n-1 hits that has a signifcantly better 
  // weight than the best n hit segment

  // IBL 070828: implicit assumption here is that there will always only be one segment per 
  // cluster - if there are >1 we will need to find out which segment the alternative n-1 hit 
  // protosegment belongs to!


  //float chosen_weight = min_weight_A;
  //float chosen_ywgt = best_weight_B;
  //float chosen_curv = best_curv_A;
  //int chosen_nlayers = n_layers_occupied_tot;
  int chosen_pseg = best_pseg;
  if (best_pseg<0) { 
    return segmentInChamber; 
  }
  chosen_Psegments = (Psegments);
  chosen_weight_A = (weight_A);

  float hit_drop_limit = -999999.999;

  // define different weight improvement requirements depending on how many layers are in the segment candidate
  switch ( n_layers_processed ) {
  case 1 : 
    // do nothing;
    break;
  case 2 :
    // do nothing;
    break;
  case 3 : 
    // do nothing;
    break;
  case 4 : 
    hit_drop_limit =  hitDropLimit6Hits * (1./2.) * hitDropLimit4Hits;
    if ((best_Layer_noLx < 1) || (best_Layer_noLx > 4)) {
      //      std::cout<<"CSCSegAlgoST: For four layers, best_Layer_noLx = "<< best_Layer_noLx << std::endl;
    }
    if ((best_Layer_noLx == 2) || (best_Layer_noLx == 3)) hit_drop_limit = hit_drop_limit * (1./2.); 
    break;
  case 5 : 
    hit_drop_limit =  hitDropLimit6Hits * (2./3.) * hitDropLimit5Hits;
    if ((best_Layer_noLx < 1) || (best_Layer_noLx > 5)) {
      //      std::cout<<"CSCSegAlgoST: For five layers, best_Layer_noLx = "<< best_Layer_noLx << std::endl;
    }
    if ((best_Layer_noLx == 2) || (best_Layer_noLx == 4)) hit_drop_limit = hit_drop_limit * (1./2.); 
    if (best_Layer_noLx == 3) hit_drop_limit = hit_drop_limit * (1./3.); 
    break;
  case 6 : 
    hit_drop_limit =  hitDropLimit6Hits * (3./4.);
    if ((best_Layer_noLx < 1) || (best_Layer_noLx > 6)) {
      //      std::cout<<"CSCSegAlgoST: For six layers, best_Layer_noLx = "<< best_Layer_noLx << std::endl;
    }
    if ((best_Layer_noLx == 2) || (best_Layer_noLx == 5)) hit_drop_limit = hit_drop_limit * (1./2.); 
    if ((best_Layer_noLx == 3) || (best_Layer_noLx == 4)) hit_drop_limit = hit_drop_limit * (1./3.); 
    break;
    
  default : 
    // Fallback - should never occur.
    LogTrace("CSCSegment|CSC") << "[CSCSegAlgoST::buildSegments] Unexpected number of layers with hits - please inform CSC DPG.";
    hit_drop_limit = 0.1;
  }

  // choose the NoLx collection (the one that contains the best N-1 candidate)
  switch ( best_Layer_noLx ) {
  case 1 : 
    Psegments_noLx.clear();
    Psegments_noLx = Psegments_noL1;
    weight_noLx_A.clear();
    weight_noLx_A = weight_noL1_A;
    break;
  case 2 :
    Psegments_noLx.clear();
    Psegments_noLx = Psegments_noL2;
    weight_noLx_A.clear();
    weight_noLx_A = weight_noL2_A;
    break;
  case 3 : 
    Psegments_noLx.clear();
    Psegments_noLx = Psegments_noL3;
    weight_noLx_A.clear();
    weight_noLx_A = weight_noL3_A;
    break;
  case 4 : 
    Psegments_noLx.clear();
    Psegments_noLx = Psegments_noL4;
    weight_noLx_A.clear();
    weight_noLx_A = weight_noL4_A;
    break;
  case 5 : 
    Psegments_noLx.clear();
    Psegments_noLx = Psegments_noL5;
    weight_noLx_A.clear();
    weight_noLx_A = weight_noL5_A;
    break;
  case 6 : 
    Psegments_noLx.clear();
    Psegments_noLx = Psegments_noL6;
    weight_noLx_A.clear();
    weight_noLx_A = weight_noL6_A;
    break;
    
  default : 
    // Fallback - should occur only for preclusters with only 3 layers with hits.
    Psegments_noLx.clear();
    weight_noLx_A.clear();
  }
  
  if( min_weight_A > 0. ) {
    if ( min_weight_noLx_A/min_weight_A < hit_drop_limit ) {
      //chosen_weight = min_weight_noLx_A;
      //chosen_ywgt = best_weight_noLx_B;
      //chosen_curv = best_curv_noLx_A;
      //chosen_nlayers = n_layers_occupied_tot-1;
      chosen_pseg = best_noLx_pseg;
      chosen_Psegments.clear();
      chosen_weight_A.clear();
      chosen_Psegments = (Psegments_noLx);
      chosen_weight_A = (weight_noLx_A);
    }
  }

  if(onlyBestSegment) {
    ChooseSegments2a( chosen_Psegments, chosen_pseg );
  }
  else {
    ChooseSegments3( chosen_Psegments, chosen_weight_A, chosen_pseg ); 
  }

  for(unsigned int iSegment=0; iSegment<GoodSegments.size();++iSegment){
    protoSegment = GoodSegments[iSegment];

    // Create new fitter object
    CSCCondSegFit* segfit = new CSCCondSegFit( pset(), chamber(), protoSegment );
    condpass1 = false;
    condpass2 = false;
    segfit->setScaleXError( 1.0 );
    segfit->fit(condpass1, condpass2);

    // Attempt to handle numerical instability of the fit.
    // (Any segment with chi2/dof > chi2Norm_3D_ is considered
    // as potentially suffering from numerical instability in fit.)
    if( adjustCovariance() ){
    // Call the fit with prefitting option:
    // First fit a straight line to X-Z coordinates and calculate chi2
    // This is done in CSCCondSegFit::correctTheCovX()
    // Scale up errors in X if this chi2 is too big (default 'too big' is >20);
    // Then refit XY-Z with the scaled-up X errors 
      if(segfit->chi2()/segfit->ndof()>chi2Norm_3D_){
    condpass1 = true;
    segfit->fit(condpass1, condpass2);
      }
      if(segfit->scaleXError() < 1.00005){
        LogTrace("CSCWeirdSegment") << "[CSCSegAlgoST::buildSegments] Segment ErrXX scaled and refit " << std::endl;
    if(segfit->chi2()/segfit->ndof()>chi2Norm_3D_){
     // Call the fit with direct adjustment of condition number;
     // If the attempt to improve fit by scaling up X error fails
     // call the procedure to make the condition number of M compatible with
     // the precision of X and Y measurements;
     // Achieved by decreasing abs value of the Covariance
          LogTrace("CSCWeirdSegment") << "[CSCSegAlgoST::buildSegments] Segment ErrXY changed to match cond. number and refit " << std::endl;
      condpass2 = true;
      segfit->fit(condpass1, condpass2);
        }
      }
      // Call the pre-pruning procedure;
      // If the attempt to improve fit by scaling up X error is successfull,
      // while scale factor for X errors is too big.
      // Prune the recHit inducing the biggest contribution into X-Z chi^2
      // and refit;
      if(prePrun_ && (sqrt(segfit->scaleXError())>prePrunLimit_) &&
     (segfit->nhits()>3)){   
        LogTrace("CSCWeirdSegment") << "[CSCSegAlgoST::buildSegments] Scale factor chi2uCorrection too big, pre-Prune and refit " << std::endl;
    protoSegment.erase(protoSegment.begin() + segfit->worstHit(),
               protoSegment.begin() + segfit->worstHit()+1 );                 

    // Need to create new fitter object to repeat fit with fewer hits
    // Original code maintained current values of condpass1, condpass2, scaleXError - calc in CorrectTheCovX()
    //@@ DO THE SAME THING HERE, BUT IS THAT CORRECT?! It does make a difference. 
    double tempcorr = segfit->scaleXError(); // save current value
    delete segfit;   
    segfit = new CSCCondSegFit( pset(), chamber(), protoSegment );           
        segfit->setScaleXError( tempcorr ); // reset to previous value (rather than init to 1)
    segfit->fit(condpass1, condpass2);
      }
    }
    
    // calculate covariance matrix 
    //    AlgebraicSymMatrix temp2 = segfit->covarianceMatrix();

    // build an actual CSC segment
    CSCSegment temp(protoSegment, segfit->intercept(), segfit->localdir(), 
                       segfit->covarianceMatrix(), segfit->chi2() );
    delete segfit;

    if( debug ) dumpSegment( temp );
 
    segmentInChamber.push_back(temp); 
  }
  return segmentInChamber;
}

std::vector<CSCSegment> CSCSegAlgoST::buildSegments2

(

const ChamberHitContainer &

rechits

)

Build track segments in this chamber (this is where the actual segment-building algorithm hides.)

std::vector< std::vector< const CSCRecHit2D * > > CSCSegAlgoST::chainHits	(	const CSCChamber *	aChamber,
		const ChamberHitContainer &	rechits
	)

Definition at line 512 of file CSCSegAlgoST.cc.

References begin, CSCChamberSpecs::chamberTypeName(), end, CSCChamberSpecs::gangedStrips(), mps_fire::i, isGoodToMerge(), CSCChamber::specs(), and groupFilesInBlocks::temp.

Referenced by run().

                                                                                                                                  {

  std::vector<ChamberHitContainer> rechits_chains; // this is a collection of groups of rechits


  std::vector<const CSCRecHit2D*> temp;

  std::vector< ChamberHitContainer > seeds;

  std::vector <bool> usedCluster;

  // split rechits into subvectors and return vector of vectors:
  // Loop over rechits
  // Create one seed per hit
  //std::cout<<" rechits.size() = "<<rechits.size()<<std::endl;
  for(unsigned int i = 0; i < rechits.size(); ++i) {
    temp.clear();
    temp.push_back(rechits[i]);
    seeds.push_back(temp);
    usedCluster.push_back(false);
  }
  // Only ME1/1A can have ganged strips so no need to test name
  bool gangedME11a = false;
  if ( ("ME1/a" == aChamber->specs()->chamberTypeName()) && aChamber->specs()->gangedStrips() ){
  //  if ( aChamber->specs()->gangedStrips() ){
    gangedME11a = true;
  }
  // merge chains that are too close ("touch" each other)
  for(size_t NNN = 0; NNN < seeds.size(); ++NNN) {
    for(size_t MMM = NNN+1; MMM < seeds.size(); ++MMM) {
      if(usedCluster[MMM] || usedCluster[NNN]){
        continue;
      }
      // all is in the way we define "good";
      // try not to "cluster" the hits but to "chain" them;
      // it does the clustering but also does a better job
      // for inclined tracks (not clustering them together;
      // crossed tracks would be still clustered together) 
      // 22.12.09: In fact it is not much more different 
      // than the "clustering", we just introduce another
      // variable in the game - Z. And it makes sense 
      // to re-introduce Y (or actually wire group mumber)
      // in a similar way as for the strip number - see
      // the code below.
      bool goodToMerge  = isGoodToMerge(gangedME11a, seeds[NNN], seeds[MMM]);
      if(goodToMerge){
        // merge chains!
        // merge by adding seed NNN to seed MMM and erasing seed NNN

        // add seed NNN to MMM (lower to larger number)
        seeds[MMM].insert(seeds[MMM].end(),seeds[NNN].begin(),seeds[NNN].end());

        // mark seed NNN as used
        usedCluster[NNN] = true;
        // we have merged a seed (NNN) to the highter seed (MMM) - need to contimue to
        // next seed (NNN+1)
        break;
      }

    }
  }

  // hand over the final seeds to the output
  // would be more elegant if we could do the above step with
  // erasing the merged ones, rather than the

  for(size_t NNN = 0; NNN < seeds.size(); ++NNN) {
    if(usedCluster[NNN]) continue; //skip seeds that have been marked as used up in merging
    rechits_chains.push_back(seeds[NNN]);
  }

  //***************************************************************

      return rechits_chains;
}

const CSCChamber* CSCSegAlgoST::chamber ( ) const

inlineprivate

Definition at line 112 of file CSCSegAlgoST.h.

References theChamber.

Referenced by buildSegments(), dumpSegment(), geometryXMLparser.CSCAlignable::index(), and prune_bad_hits().

{return theChamber;}

void CSCSegAlgoST::ChooseSegments ( void  )

private

Referenced by adjustCovariance().

void CSCSegAlgoST::ChooseSegments2 ( int  best_seg )

private

Definition at line 1681 of file CSCSegAlgoST.cc.

References GoodSegments, createfilelist::int, LogTrace, Psegments, findQualityFiles::size, and weight_A.

Referenced by adjustCovariance().

                                               {
  //  std::vector <int> CommonHits(6); // nice  concept :)
  std::vector <unsigned int> BadCandidate;
  int SumCommonHits =0;
  GoodSegments.clear();
  BadCandidate.clear();
  for(unsigned int iCand=0;iCand<Psegments.size();++iCand) {
    // skip here if segment was marked bad
    for(unsigned int iiCand=iCand+1;iiCand<Psegments.size();++iiCand){
      // skip here too if segment was marked bad
      SumCommonHits =0;
      if( Psegments[iCand].size() != Psegments[iiCand].size() ) {
    LogTrace("CSCSegment|CSC") << "[CSCSegAlgoST::ChooseSegments2] ALARM!! Should not happen - please inform CSC DPG";
      }
      else {
    for( int ihits = 0; ihits < int(Psegments[iCand].size()); ++ihits ) { // iCand and iiCand NEED to have same nr of hits! (alsways have by construction)
      if( Psegments[iCand][ihits] == Psegments[iiCand][ihits]) {
        ++SumCommonHits;
      }
    }
      }
      if(SumCommonHits>1) {
    if( weight_A[iCand]>weight_A[iiCand] ) { // use weight_A here
      BadCandidate.push_back(iCand);
      // rather mark segment bad by an array which is in sync with protosegments!! e.g. set weight = weight*1000 or have an addidional array or set it to weight *= -1
    }
    else{
      BadCandidate.push_back(iiCand);
      // rather mark segment bad by an array which is in sync with protosegments!! e.g. set weight = weight*1000 or have an addidional array or set it to weight *= -1
    }
      }
    }
  }
  bool discard;
  for(unsigned int isegm=0;isegm<Psegments.size();++isegm) {
    // For best results another iteration/comparison over Psegments 
    //should be applied here... It would make the program much slower.
    discard = false;
    for(unsigned int ibad=0;ibad<BadCandidate.size();++ibad) {
      // can save this loop if we used an array in sync with Psegments!!!!
      if(isegm == BadCandidate[ibad]) {
    discard = true;
      }
    }
    if(!discard) {
      GoodSegments.push_back( Psegments[isegm] );
    }
  }
}

void CSCSegAlgoST::ChooseSegments2a ( std::vector< ChamberHitContainer > &  best_segments, int  best_seg  )

private

Definition at line 1623 of file CSCSegAlgoST.cc.

References GoodSegments.

Referenced by adjustCovariance(), and buildSegments().

                                                                                                      {
  // just return best segment
  GoodSegments.clear();
  GoodSegments.push_back( chosen_segments[chosen_seg] );
}

void CSCSegAlgoST::ChooseSegments3 ( int  best_seg )

private

Referenced by adjustCovariance(), and buildSegments().

void CSCSegAlgoST::ChooseSegments3 ( std::vector< ChamberHitContainer > &  best_segments, std::vector< float > &  best_weight, int  best_seg  )

private

Definition at line 1629 of file CSCSegAlgoST.cc.

References GoodSegments, createfilelist::int, and findQualityFiles::size.

                                                                                                                                         {

  int SumCommonHits = 0;
  GoodSegments.clear();
  int nr_remaining_candidates;
  unsigned int nr_of_segment_candidates;
  
  nr_remaining_candidates = nr_of_segment_candidates = chosen_segments.size();

  // always select and return best protosegment:  
  GoodSegments.push_back( chosen_segments[ chosen_seg ] );

  float chosen_weight_temp = 999999.;
  int chosen_seg_temp = -1;

  // try to find further segment candidates:
  while( nr_remaining_candidates > 0 ) {

    for(unsigned int iCand=0; iCand < nr_of_segment_candidates; ++iCand) {
      //only compare current best to psegs that have not been marked bad:
      if( chosen_weight[iCand] < 0. ) continue;
      SumCommonHits = 0;

      for( int ihits = 0; ihits < int(chosen_segments[iCand].size()); ++ihits ) { // iCand and iiCand NEED to have same nr of hits! (always have by construction)
    if( chosen_segments[iCand][ihits] == chosen_segments[chosen_seg][ihits]) {
      ++SumCommonHits;
    }
      }

      //mark a pseg bad:
      if(SumCommonHits>1) { // needs to be a card; should be investigated first
    chosen_weight[iCand] = -1.;
    nr_remaining_candidates -= 1;
      }
      else {
    // save the protosegment with the smallest weight
    if( chosen_weight[ iCand ] < chosen_weight_temp ) {
      chosen_weight_temp = chosen_weight[ iCand ];
      chosen_seg_temp = iCand ;
    }
      }
    }

    if( chosen_seg_temp > -1 ) GoodSegments.push_back( chosen_segments[ chosen_seg_temp ] );

    chosen_seg = chosen_seg_temp;
    // re-initialze temporary best parameters
    chosen_weight_temp = 999999;
    chosen_seg_temp = -1;
  }
}

std::vector< std::vector< const CSCRecHit2D * > > CSCSegAlgoST::clusterHits	(	const CSCChamber *	aChamber,
		const ChamberHitContainer &	rechits
	)

Build groups of rechits that are separated in x and y to save time on the segment finding

Definition at line 396 of file CSCSegAlgoST.cc.

References begin, dXclusBoxMax, dYclusBoxMax, end, mps_fire::i, LogTrace, findQualityFiles::size, groupFilesInBlocks::temp, theChamber, x, and y.

Referenced by run().

                                                                                                                                    {
  theChamber = aChamber; 

  std::vector<ChamberHitContainer> rechits_clusters; // this is a collection of groups of rechits
  //   const float dXclus_box_cut       = 4.; // seems to work reasonably 070116
  //   const float dYclus_box_cut       = 8.; // seems to work reasonably 070116

  //float dXclus = 0.0;
  //float dYclus = 0.0;
  float dXclus_box = 0.0;
  float dYclus_box = 0.0;

  std::vector<const CSCRecHit2D*> temp;

  std::vector< ChamberHitContainer > seeds;

  std::vector<float> running_meanX;
  std::vector<float> running_meanY;

  std::vector<float> seed_minX;
  std::vector<float> seed_maxX;
  std::vector<float> seed_minY;
  std::vector<float> seed_maxY;

  //std::cout<<"*************************************************************"<<std::endl;
  //std::cout<<"Called clusterHits in Chamber "<< theChamber->specs()->chamberTypeName()<<std::endl;
  //std::cout<<"*************************************************************"<<std::endl;

  // split rechits into subvectors and return vector of vectors:
  // Loop over rechits 
  // Create one seed per hit
  for(unsigned int i = 0; i < rechits.size(); ++i) {

    temp.clear();

    temp.push_back(rechits[i]);

    seeds.push_back(temp);

    // First added hit in seed defines the mean to which the next hit is compared
    // for this seed.

    running_meanX.push_back( rechits[i]->localPosition().x() );
    running_meanY.push_back( rechits[i]->localPosition().y() );
    
    // set min/max X and Y for box containing the hits in the precluster:
    seed_minX.push_back( rechits[i]->localPosition().x() );
    seed_maxX.push_back( rechits[i]->localPosition().x() );
    seed_minY.push_back( rechits[i]->localPosition().y() );
    seed_maxY.push_back( rechits[i]->localPosition().y() );
  }
    
  // merge clusters that are too close
  // measure distance between final "running mean"
  for(size_t NNN = 0; NNN < seeds.size(); ++NNN) {
    
    for(size_t MMM = NNN+1; MMM < seeds.size(); ++MMM) {
      if(running_meanX[MMM] == 999999. || running_meanX[NNN] == 999999. ) {
    LogTrace("CSCSegment|CSC") << "[CSCSegAlgoST::clusterHits] ALARM! Skipping used seeds, this should not happen - inform CSC DPG";
    //  std::cout<<"We should never see this line now!!!"<<std::endl;
    continue; //skip seeds that have been used 
      }
      
      // calculate cut criteria for simple running mean distance cut:
      //dXclus = fabs(running_meanX[NNN] - running_meanX[MMM]);
      //dYclus = fabs(running_meanY[NNN] - running_meanY[MMM]);

      // calculate minmal distance between precluster boxes containing the hits:
      if ( running_meanX[NNN] > running_meanX[MMM] ) dXclus_box = seed_minX[NNN] - seed_maxX[MMM];
      else                                           dXclus_box = seed_minX[MMM] - seed_maxX[NNN];
      if ( running_meanY[NNN] > running_meanY[MMM] ) dYclus_box = seed_minY[NNN] - seed_maxY[MMM];
      else                                           dYclus_box = seed_minY[MMM] - seed_maxY[NNN];
      
      
      if( dXclus_box < dXclusBoxMax && dYclus_box < dYclusBoxMax ) {
    // merge clusters!
    // merge by adding seed NNN to seed MMM and erasing seed NNN
        
    // calculate running mean for the merged seed:
    running_meanX[MMM] = (running_meanX[NNN]*seeds[NNN].size() + running_meanX[MMM]*seeds[MMM].size()) / (seeds[NNN].size()+seeds[MMM].size());
    running_meanY[MMM] = (running_meanY[NNN]*seeds[NNN].size() + running_meanY[MMM]*seeds[MMM].size()) / (seeds[NNN].size()+seeds[MMM].size());
        
    // update min/max X and Y for box containing the hits in the merged cluster:
    if ( seed_minX[NNN] <= seed_minX[MMM] ) seed_minX[MMM] = seed_minX[NNN];
    if ( seed_maxX[NNN] >  seed_maxX[MMM] ) seed_maxX[MMM] = seed_maxX[NNN];
    if ( seed_minY[NNN] <= seed_minY[MMM] ) seed_minY[MMM] = seed_minY[NNN];
    if ( seed_maxY[NNN] >  seed_maxY[MMM] ) seed_maxY[MMM] = seed_maxY[NNN];
        
    // add seed NNN to MMM (lower to larger number)
    seeds[MMM].insert(seeds[MMM].end(),seeds[NNN].begin(),seeds[NNN].end());
        
    // mark seed NNN as used (at the moment just set running mean to 999999.)
    running_meanX[NNN] = 999999.;
    running_meanY[NNN] = 999999.;
    // we have merged a seed (NNN) to the highter seed (MMM) - need to contimue to 
    // next seed (NNN+1)
    break;
      }

    }
  }

  // hand over the final seeds to the output
  // would be more elegant if we could do the above step with 
  // erasing the merged ones, rather than the 
  for(size_t NNN = 0; NNN < seeds.size(); ++NNN) {
    if(running_meanX[NNN] == 999999.) continue; //skip seeds that have been marked as used up in merging
    rechits_clusters.push_back(seeds[NNN]);
  }

  //***************************************************************

      return rechits_clusters; 
}

void CSCSegAlgoST::dumpSegment ( const CSCSegment &  seg ) const

private

Definition at line 1759 of file CSCSegAlgoST.cc.

References chamber(), CSCSegment::chi2(), CSCSegment::degreesOfFreedom(), CSCChamber::id(), CSCSegment::localDirection(), CSCSegment::localDirectionError(), CSCSegment::localPosition(), CSCSegment::localPositionError(), CSCSegment::parametersError(), CSCSegment::specificRecHits(), and CSCSegment::time().

Referenced by adjustCovariance(), and buildSegments().

                                                            {

  // Only called if pset value 'CSCDebug' is set in config

  edm::LogVerbatim("CSCSegment") << "CSCSegment in " << chamber()->id()
                         << "\nlocal position = " << seg.localPosition()
             << "\nerror = " << seg.localPositionError() 
             << "\nlocal direction = " << seg.localDirection()
             << "\nerror =" << seg.localDirectionError()
             << "\ncovariance matrix" 
             << seg.parametersError()
             << "chi2/ndf = " << seg.chi2() << "/" << seg.degreesOfFreedom()
                         << "\n#rechits = " << seg.specificRecHits().size()
                         << "\ntime = " << seg.time();
}

void CSCSegAlgoST::findDuplicates ( std::vector< CSCSegment > &  segments )

private

Definition at line 1732 of file CSCSegAlgoST.cc.

Referenced by adjustCovariance(), and run().

                                                                   {
  // this is intended for ME1/1a only - we have ghost segments because of the strips ganging 
  // this function finds them (first the rechits by sharesInput() )
  // if a segment shares all the rechits with another segment it is a duplicate (even if
  // it has less rechits) 
  
  for(std::vector<CSCSegment>::iterator it=segments.begin(); it != segments.end(); ++it) {
    std::vector<CSCSegment*> duplicateSegments;
    for(std::vector<CSCSegment>::iterator it2=segments.begin(); it2 != segments.end(); ++it2) {
      //
      bool allShared = true;
      if(it!=it2){
    allShared = it->sharesRecHits(*it2);
      }
      else{
        allShared = false;
      }
      //
      if(allShared){
        duplicateSegments.push_back(&(*it2));
      }
    }
    it->setDuplicateSegments(duplicateSegments);
  }

}

bool CSCSegAlgoST::isGoodToMerge ( bool  isME11a, ChamberHitContainer &  newChain, ChamberHitContainer &  oldChain  )

private

Definition at line 588 of file CSCSegAlgoST.cc.

Referenced by adjustCovariance(), and chainHits().

                                                                                                                 {
  for(size_t iRH_new = 0;iRH_new<newChain.size();++iRH_new){
    int layer_new = newChain[iRH_new]->cscDetId().layer()-1;     
    int middleStrip_new = newChain[iRH_new]->nStrips()/2;
    int centralStrip_new = newChain[iRH_new]->channels(middleStrip_new);
    int centralWire_new = newChain[iRH_new]->hitWire();
    bool layerRequirementOK = false;
    bool stripRequirementOK = false;
    bool wireRequirementOK = false;
    bool goodToMerge = false;
    for(size_t iRH_old = 0;iRH_old<oldChain.size();++iRH_old){      
      int layer_old = oldChain[iRH_old]->cscDetId().layer()-1;
      int middleStrip_old = oldChain[iRH_old]->nStrips()/2;
      int centralStrip_old = oldChain[iRH_old]->channels(middleStrip_old);
      int centralWire_old = oldChain[iRH_old]->hitWire();

      // to be chained, two hits need to be in neighbouring layers...
      // or better allow few missing layers (upto 3 to avoid inefficiencies);
      // however we'll not make an angle correction because it
      // worsen the situation in some of the "regular" cases 
      // (not making the correction means that the conditions for
      // forming a cluster are different if we have missing layers -
      // this could affect events at the boundaries ) 
      if(layer_new==layer_old+1 || layer_new==layer_old-1 ||
     layer_new==layer_old+2 || layer_new==layer_old-2 ||
     layer_new==layer_old+3 || layer_new==layer_old-3 ||
     layer_new==layer_old+4 || layer_new==layer_old-4 ){
        layerRequirementOK = true;
      }
      int allStrips = 48;
      //to be chained, two hits need to be "close" in strip number (can do it in phi
      // but it doesn't really matter); let "close" means upto 2 strips (3?) - 
      // this is more compared to what CLCT readout patterns allow 
      if(centralStrip_new==centralStrip_old ||
         centralStrip_new==centralStrip_old+1 || centralStrip_new==centralStrip_old-1 ||
         centralStrip_new==centralStrip_old+2 || centralStrip_new==centralStrip_old-2){
        stripRequirementOK = true;
      }
      // same for wires (and ALCT patterns)
      if(centralWire_new==centralWire_old ||
         centralWire_new==centralWire_old+1 || centralWire_new==centralWire_old-1 ||
         centralWire_new==centralWire_old+2 || centralWire_new==centralWire_old-2){
        wireRequirementOK = true;
      }

      if(gangedME11a){
    if(centralStrip_new==centralStrip_old+1-allStrips || centralStrip_new==centralStrip_old-1-allStrips ||
       centralStrip_new==centralStrip_old+2-allStrips || centralStrip_new==centralStrip_old-2-allStrips ||
       centralStrip_new==centralStrip_old+1+allStrips || centralStrip_new==centralStrip_old-1+allStrips ||
       centralStrip_new==centralStrip_old+2+allStrips || centralStrip_new==centralStrip_old-2+allStrips){
      stripRequirementOK = true;
    }
      }
      if(layerRequirementOK && stripRequirementOK && wireRequirementOK){
        goodToMerge = true;
        return goodToMerge;
      }
    }
  }
  return false;
}

std::vector< CSCSegment > CSCSegAlgoST::prune_bad_hits	(	const CSCChamber *	aChamber,
		std::vector< CSCSegment > &	segments
	)

Remove bad hits from found segments based not only on chi2, but also on charge and further "low level" chamber information.

Definition at line 181 of file CSCSegAlgoST.cc.

References adjustCovariance(), begin, BPMinImprovement, BrutePruning, edmScanValgrind::buffer, chamber(), CSCSegFit::chi2(), CSCSegmentAlgorithmTC_cfi::chi2ndfProbMin, chi2Norm_3D_, ChiSquaredProbability(), condpass1, condpass2, CSCSegFit::covarianceMatrix(), MillePedeFileConverter_cfg::e, CSCCondSegFit::fit(), CSCSegFit::intercept(), CSCChamber::layer(), CSCSegFit::localdir(), funct::m, minHitsPerSegment, CSCSegFit::ndof(), nhits, protoSegment, pset(), CSCSegFit::scaleXError(), CSCSegFit::setScaleXError(), groupFilesInBlocks::temp, theChamber, GeomDet::toGlobal(), PV3DBase< T, PVType, FrameType >::x(), LocalError::xx(), and PV3DBase< T, PVType, FrameType >::z().

Referenced by run().

                                                                                                               {
  
  //   std::cout<<"*************************************************************"<<std::endl;
  //   std::cout<<"Called prune_bad_hits in Chamber "<< theChamber->specs()->chamberTypeName()<<std::endl;
  //   std::cout<<"*************************************************************"<<std::endl;
  
  std::vector<CSCSegment>          segments_temp;
  std::vector<ChamberHitContainer> rechits_clusters; // this is a collection of groups of rechits
  
  const float chi2ndfProbMin = 1.0e-4;
  bool   use_brute_force = BrutePruning;

  int hit_nr = 0;
  int hit_nr_worst = -1;
  //int hit_nr_2ndworst = -1;
  
  for(std::vector<CSCSegment>::iterator it=segments.begin(); it != segments.end(); ++it) {
    
    // do nothing for nhit <= minHitPerSegment
    if( (*it).nRecHits() <= minHitsPerSegment ) continue;
    
    if( !use_brute_force ) {// find worst hit
      
      float chisq    = (*it).chi2();
      int nhits      = (*it).nRecHits();
      LocalPoint localPos = (*it).localPosition();
      LocalVector segDir = (*it).localDirection();
      const CSCChamber* cscchamber = theChamber;
      float globZ       ;
      
      GlobalPoint globalPosition = cscchamber->toGlobal(localPos);
      globZ = globalPosition.z();
      
      
      if( ChiSquaredProbability((double)chisq,(double)(2*nhits-4)) < chi2ndfProbMin  ) {

    // find (rough) "residuals" (NOT excluding the hit from the fit - speed!) of hits on segment
    std::vector<CSCRecHit2D> theseRecHits = (*it).specificRecHits();
    std::vector<CSCRecHit2D>::const_iterator iRH_worst;
    //float xdist_local       = -99999.;

    float xdist_local_worst_sig = -99999.;
    float xdist_local_2ndworst_sig = -99999.;
    float xdist_local_sig       = -99999.;

    hit_nr = 0;
    hit_nr_worst = -1;
    //hit_nr_2ndworst = -1;

    for ( std::vector<CSCRecHit2D>::const_iterator iRH = theseRecHits.begin(); iRH != theseRecHits.end(); ++iRH) {
      //mark "worst" hit:
      
      //float z_at_target ;
      //float radius      ;
      float loc_x_at_target;
      //float loc_y_at_target;
      //float loc_z_at_target;

      //z_at_target  = 0.;
      //radius       = 0.;
      
      // set the z target in CMS global coordinates:
      const CSCLayer* csclayerRH = theChamber->layer((*iRH).cscDetId().layer());
      LocalPoint localPositionRH = (*iRH).localPosition();
      GlobalPoint globalPositionRH = csclayerRH->toGlobal(localPositionRH); 
      
      LocalError rerrlocal = (*iRH).localPositionError();  
      float xxerr = rerrlocal.xx();
      
      float target_z     = globalPositionRH.z();  // target z position in cm (z pos of the hit)
      
      if(target_z > 0.) {
        loc_x_at_target = localPos.x() + (segDir.x()/fabs(segDir.z())*( target_z - globZ ));
        //loc_y_at_target = localPos.y() + (segDir.y()/fabs(segDir.z())*( target_z - globZ ));
        //loc_z_at_target = target_z;
      }
      else {
        loc_x_at_target = localPos.x() + ((-1)*segDir.x()/fabs(segDir.z())*( target_z - globZ ));
        //loc_y_at_target = localPos.y() + ((-1)*segDir.y()/fabs(segDir.z())*( target_z - globZ ));
        //loc_z_at_target = target_z;
      }
      // have to transform the segments coordinates back to the local frame... how?!!!!!!!!!!!!
      
      //xdist_local  = fabs(localPositionRH.x() - loc_x_at_target);
      xdist_local_sig  = fabs((localPositionRH.x() -loc_x_at_target)/(xxerr));
      
      if( xdist_local_sig > xdist_local_worst_sig ) {
        xdist_local_2ndworst_sig = xdist_local_worst_sig;
        xdist_local_worst_sig    = xdist_local_sig;
        iRH_worst            = iRH;
        //hit_nr_2ndworst = hit_nr_worst;
        hit_nr_worst = hit_nr;
      }
      else if(xdist_local_sig > xdist_local_2ndworst_sig) {
        xdist_local_2ndworst_sig = xdist_local_sig;
        //hit_nr_2ndworst = hit_nr;
      }
      ++hit_nr;
    }

    // reset worst hit number if certain criteria apply.
    // Criteria: 2nd worst hit must be at least a factor of
    // 1.5 better than the worst in terms of sigma:
    if ( xdist_local_worst_sig / xdist_local_2ndworst_sig < 1.5 ) {
      hit_nr_worst    = -1;
      //hit_nr_2ndworst = -1;
    }
      }
    }

    // if worst hit was found, refit without worst hit and select if considerably better than original fit.
    // Can also use brute force: refit all n-1 hit segments and choose one over the n hit if considerably "better"
   
    std::vector< CSCRecHit2D > buffer;
    std::vector< std::vector< CSCRecHit2D > > reduced_segments;
    std::vector< CSCRecHit2D > theseRecHits = (*it).specificRecHits();
    float best_red_seg_prob = 0.0;
    // usefor chi2 1 diff   float best_red_seg_prob = 99999.;
    buffer.clear();

    if( ChiSquaredProbability((double)(*it).chi2(),(double)((2*(*it).nRecHits())-4)) < chi2ndfProbMin  ) {
    
      buffer = theseRecHits;

      // Dirty switch: here one can select to refit all possible subsets or just the one without the 
      // tagged worst hit:
      if( use_brute_force ) { // Brute force method: loop over all possible segments:
    for(size_t bi = 0; bi < buffer.size(); ++bi) {
      reduced_segments.push_back(buffer);
      reduced_segments[bi].erase(reduced_segments[bi].begin()+(bi),reduced_segments[bi].begin()+(bi+1));
    }
      }
      else { // More elegant but still biased: erase only worst hit
    // Comment: There is not a very strong correlation of the worst hit with the one that one should remove... 
    if( hit_nr_worst >= 0 && hit_nr_worst <= int(buffer.size())  ) {
      // fill segment in buffer, delete worst hit
      buffer.erase(buffer.begin()+(hit_nr_worst),buffer.begin()+(hit_nr_worst+1));
      reduced_segments.push_back(buffer);
    }
    else {
      // only fill segment in array, do not delete anything
      reduced_segments.push_back(buffer);
    }
      }
    }
      
    // Loop over the subsegments and fit (only one segment if "use_brute_force" is false):
    for(size_t iSegment=0; iSegment<reduced_segments.size(); ++iSegment) {
      // loop over hits on given segment and push pointers to hits into protosegment
      protoSegment.clear();
      for(size_t m = 0; m<reduced_segments[iSegment].size(); ++m ) {
    protoSegment.push_back(&reduced_segments[iSegment][m]);
      }

      // Create fitter object
      CSCCondSegFit* segfit = new CSCCondSegFit( pset(), chamber(), protoSegment );
      condpass1 = false;
      condpass2 = false;
      segfit->setScaleXError( 1.0 );
      segfit->fit(condpass1, condpass2);

      // Attempt to handle numerical instability of the fit;
      // The same as in the build method;
      // Preprune is not applied;
      if( adjustCovariance() ){
    if(segfit->chi2()/segfit->ndof()>chi2Norm_3D_){
      condpass1 = true;
      segfit->fit(condpass1, condpass2);
    }
    if( (segfit->scaleXError() < 1.00005)&&(segfit->chi2()/segfit->ndof()>chi2Norm_3D_) ){
      condpass2 = true;
      segfit->fit(condpass1, condpass2);
    }
      }
    
      // calculate error matrix 
      //      AlgebraicSymMatrix temp2 = segfit->covarianceMatrix();

      // build an actual segment
      CSCSegment temp(protoSegment, segfit->intercept(), segfit->localdir(), 
                         segfit->covarianceMatrix(), segfit->chi2() );

      // and finished with this fit
      delete segfit;

      // n-hit segment is (*it)
      // (n-1)-hit segment is temp
      // replace n-hit segment with (n-1)-hit segment if segment probability is BPMinImprovement better
      double oldchi2 = (*it).chi2();
      double olddof  =  2 * (*it).nRecHits() - 4;
      double newchi2 = temp.chi2();
      double newdof  = 2 * temp.nRecHits() - 4;
      if( ( ChiSquaredProbability(oldchi2,olddof) < (1./BPMinImprovement)*
        ChiSquaredProbability(newchi2,newdof) ) 
      && ( ChiSquaredProbability(newchi2,newdof) > best_red_seg_prob )
      && ( ChiSquaredProbability(newchi2,newdof) > 1e-10 )
    ) {
    best_red_seg_prob = ChiSquaredProbability(newchi2,newdof);
        // The (n-1)- segment is better than the n-hit segment.
    // If it has at least  minHitsPerSegment  hits replace the n-hit segment
        // with this  better (n-1)-hit segment:
        if( temp.nRecHits() >= minHitsPerSegment ) {
          (*it) = temp;
        }
      }
    } // end of loop over subsegments (iSegment)
    
  } // end loop over segments (it)
  
  return segments;
  
}

const edm::ParameterSet& CSCSegAlgoST::pset ( void  ) const

inlineprivate

Definition at line 86 of file CSCSegAlgoST.h.

References ps_.

Referenced by buildSegments(), and prune_bad_hits().

{ return ps_;}

std::vector< CSCSegment > CSCSegAlgoST::run ( const CSCChamber *  aChamber, const ChamberHitContainer &  rechits  )

override

Build segments for all desired groups of hits

Definition at line 87 of file CSCSegAlgoST.cc.

References a, a_yweightPenaltyThreshold, b, buildSegments(), chainHits(), CSCChamberSpecs::chamberTypeName(), clusterHits(), findDuplicates(), CSCChamberSpecs::gangedStrips(), CSCChamber::id(), LogTrace, preClustering, preClustering_useChaining, prune_bad_hits(), Pruning, CSCChamber::specs(), theChamber, and yweightPenaltyThreshold.

                                                                                                      {

  // Set member variable
  theChamber = aChamber; 

  LogTrace("CSCSegAlgoST") << "[CSCSegAlgoST::run] Start building segments in chamber " << theChamber->id();

  // pre-cluster rechits and loop over all sub clusters seperately
  std::vector<CSCSegment>          segments_temp;
  std::vector<CSCSegment>          segments;
  std::vector<ChamberHitContainer> rechits_clusters; // a collection of clusters of rechits

  // Define yweight penalty depending on chamber. 
  // We fixed the relative ratios, but they can be scaled by parameters:
  
  for(int a = 0; a<5; ++a) {
    for(int b = 0; b<5; ++b) {
      a_yweightPenaltyThreshold[a][b] = 0.0;
    }
  }
  
  a_yweightPenaltyThreshold[1][1] = yweightPenaltyThreshold * 10.20;
  a_yweightPenaltyThreshold[1][2] = yweightPenaltyThreshold * 14.00;
  a_yweightPenaltyThreshold[1][3] = yweightPenaltyThreshold * 20.40;
  a_yweightPenaltyThreshold[1][4] = yweightPenaltyThreshold * 10.20;
  a_yweightPenaltyThreshold[2][1] = yweightPenaltyThreshold *  7.60;
  a_yweightPenaltyThreshold[2][2] = yweightPenaltyThreshold * 20.40;
  a_yweightPenaltyThreshold[3][1] = yweightPenaltyThreshold *  7.60;
  a_yweightPenaltyThreshold[3][2] = yweightPenaltyThreshold * 20.40;
  a_yweightPenaltyThreshold[4][1] = yweightPenaltyThreshold *  6.75;
  a_yweightPenaltyThreshold[4][2] = yweightPenaltyThreshold * 20.40;
  
  if(preClustering) {
    // run a pre-clusterer on the given rechits to split clearly-separated segment seeds:
    if(preClustering_useChaining){
      // it uses X,Y,Z information; there are no configurable parameters used;
      // the X, Y, Z "cuts" are just (much) wider than the LCT readout ones
      // (which are actually not step functions); this new code could accomodate
      // the clusterHits one below but we leave it for security and backward 
      // comparison reasons 
      rechits_clusters = chainHits( theChamber, rechits );
    }
    else{
      // it uses X,Y information + configurable parameters
      rechits_clusters = clusterHits( theChamber, rechits );
    }
    // loop over the found clusters:
    for(std::vector<ChamberHitContainer>::iterator sub_rechits = rechits_clusters.begin(); sub_rechits !=  rechits_clusters.end(); ++sub_rechits ) {
      // clear the buffer for the subset of segments:
      segments_temp.clear();
      // build the subset of segments:
      segments_temp = buildSegments( (*sub_rechits) );
      // add the found subset of segments to the collection of all segments in this chamber:
      segments.insert( segments.end(), segments_temp.begin(), segments_temp.end() );
    }
    // Any pruning of hits?
    if( Pruning ) {
      segments_temp.clear(); // segments_temp needed?!?!
      segments_temp = prune_bad_hits( theChamber, segments );
      segments.clear(); // segments_temp needed?!?!
      segments.swap(segments_temp); // segments_temp needed?!?!
    }
  
    // Ganged strips in ME1/1A?
    if ( ("ME1/a" == aChamber->specs()->chamberTypeName()) && aChamber->specs()->gangedStrips() ){
    //  if ( aChamber->specs()->gangedStrips() ){
      findDuplicates(segments);
    }
    return segments;
  }
  else {
    segments = buildSegments(rechits);
    if( Pruning ) {
      segments_temp.clear(); // segments_temp needed?!?!
      segments_temp = prune_bad_hits( theChamber, segments );
      segments.clear(); // segments_temp needed?!?!
      segments.swap(segments_temp); // segments_temp needed?!?!
    }

    // Ganged strips in ME1/1A?
    if ( ("ME1/a" == aChamber->specs()->chamberTypeName()) && aChamber->specs()->gangedStrips() ){
    //  if ( aChamber->specs()->gangedStrips() ){
      findDuplicates(segments);
    }
    return segments;
    //return buildSegments(rechits); 
  }
}

double CSCSegAlgoST::theWeight ( double  coordinate_1, double  coordinate_2, double  coordinate_3, float  layer_1, float  layer_2, float  layer_3  )

private

Utility functions.

Definition at line 653 of file CSCSegAlgoST.cc.

Referenced by adjustCovariance(), and buildSegments().

                                                                                                                                         {
  double sub_weight = 0;
  sub_weight = fabs( 
            ( (coordinate_2 - coordinate_3) / (layer_2  - layer_3) ) - 
            ( (coordinate_1 - coordinate_2) / (layer_1  - layer_2) ) 
            );
  return sub_weight;
}

Member Data Documentation

float CSCSegAlgoST::a_yweightPenaltyThreshold[5][5]

private

Definition at line 182 of file CSCSegAlgoST.h.

Referenced by buildSegments(), and run().

bool CSCSegAlgoST::adjustCovariance_

private

Definition at line 191 of file CSCSegAlgoST.h.

Referenced by adjustCovariance(), and CSCSegAlgoST().

double CSCSegAlgoST::BPMinImprovement

private

Definition at line 174 of file CSCSegAlgoST.h.

Referenced by CSCSegAlgoST(), and prune_bad_hits().

bool CSCSegAlgoST::BrutePruning

private

Definition at line 173 of file CSCSegAlgoST.h.

Referenced by CSCSegAlgoST(), and prune_bad_hits().

double CSCSegAlgoST::chi2Norm_3D_

private

Definition at line 195 of file CSCSegAlgoST.h.

Referenced by buildSegments(), CSCSegAlgoST(), and prune_bad_hits().

std::vector< ChamberHitContainer > CSCSegAlgoST::chosen_Psegments

private

Definition at line 133 of file CSCSegAlgoST.h.

Referenced by buildSegments().

std::vector< float > CSCSegAlgoST::chosen_weight_A

private

Definition at line 142 of file CSCSegAlgoST.h.

Referenced by buildSegments().

bool CSCSegAlgoST::condpass1

private

Flag whether to 'improve' covariance matrix.

Definition at line 193 of file CSCSegAlgoST.h.

Referenced by buildSegments(), and prune_bad_hits().

bool CSCSegAlgoST::condpass2

private

Definition at line 193 of file CSCSegAlgoST.h.

Referenced by buildSegments(), and prune_bad_hits().

std::vector< float > CSCSegAlgoST::curv_A

private

Definition at line 143 of file CSCSegAlgoST.h.

Referenced by buildSegments().

std::vector< float > CSCSegAlgoST::curv_noL1_A

private

Definition at line 144 of file CSCSegAlgoST.h.

Referenced by buildSegments().

std::vector< float > CSCSegAlgoST::curv_noL2_A

private

Definition at line 145 of file CSCSegAlgoST.h.

Referenced by buildSegments().

std::vector< float > CSCSegAlgoST::curv_noL3_A

private

Definition at line 146 of file CSCSegAlgoST.h.

Referenced by buildSegments().

std::vector< float > CSCSegAlgoST::curv_noL4_A

private

Definition at line 147 of file CSCSegAlgoST.h.

Referenced by buildSegments().

std::vector< float > CSCSegAlgoST::curv_noL5_A

private

Definition at line 148 of file CSCSegAlgoST.h.

Referenced by buildSegments().

std::vector< float > CSCSegAlgoST::curv_noL6_A

private

Definition at line 149 of file CSCSegAlgoST.h.

Referenced by buildSegments().

double CSCSegAlgoST::curvePenalty

private

Definition at line 188 of file CSCSegAlgoST.h.

Referenced by buildSegments(), and CSCSegAlgoST().

double CSCSegAlgoST::curvePenaltyThreshold

private

Definition at line 187 of file CSCSegAlgoST.h.

Referenced by buildSegments(), and CSCSegAlgoST().

bool CSCSegAlgoST::debug

private

Definition at line 161 of file CSCSegAlgoST.h.

Referenced by buildSegments(), and CSCSegAlgoST().

double CSCSegAlgoST::dXclusBoxMax

private

Definition at line 167 of file CSCSegAlgoST.h.

Referenced by clusterHits(), and CSCSegAlgoST().

double CSCSegAlgoST::dYclusBoxMax

private

Definition at line 168 of file CSCSegAlgoST.h.

Referenced by clusterHits(), and CSCSegAlgoST().

Segments CSCSegAlgoST::GoodSegments

private

Definition at line 120 of file CSCSegAlgoST.h.

Referenced by buildSegments(), ChooseSegments2(), ChooseSegments2a(), and ChooseSegments3().

double CSCSegAlgoST::hitDropLimit4Hits

private

Definition at line 178 of file CSCSegAlgoST.h.

Referenced by buildSegments(), and CSCSegAlgoST().

double CSCSegAlgoST::hitDropLimit5Hits

private

Definition at line 179 of file CSCSegAlgoST.h.

Referenced by buildSegments(), and CSCSegAlgoST().

double CSCSegAlgoST::hitDropLimit6Hits

private

Definition at line 180 of file CSCSegAlgoST.h.

Referenced by buildSegments(), and CSCSegAlgoST().

int CSCSegAlgoST::maxRecHitsInCluster

private

Definition at line 169 of file CSCSegAlgoST.h.

Referenced by buildSegments(), and CSCSegAlgoST().

int CSCSegAlgoST::minHitsPerSegment

private

Definition at line 164 of file CSCSegAlgoST.h.

Referenced by buildSegments(), CSCSegAlgoST(), and prune_bad_hits().

const std::string CSCSegAlgoST::myName_

private

Definition at line 115 of file CSCSegAlgoST.h.

bool CSCSegAlgoST::onlyBestSegment

private

Definition at line 175 of file CSCSegAlgoST.h.

Referenced by buildSegments(), and CSCSegAlgoST().

ChamberHitContainer CSCSegAlgoST::PAhits_onLayer[6]

private

Definition at line 122 of file CSCSegAlgoST.h.

Referenced by buildSegments().

bool CSCSegAlgoST::preClustering

private

Definition at line 170 of file CSCSegAlgoST.h.

Referenced by CSCSegAlgoST(), and run().

bool CSCSegAlgoST::preClustering_useChaining

private

Definition at line 171 of file CSCSegAlgoST.h.

Referenced by CSCSegAlgoST(), and run().

bool CSCSegAlgoST::prePrun_

private

Chi^2 normalization for the initial fit.

Definition at line 197 of file CSCSegAlgoST.h.

Referenced by buildSegments(), and CSCSegAlgoST().

double CSCSegAlgoST::prePrunLimit_

private

Allow to prune a (rechit in a) segment in segment buld method once it passed through Chi^2-X and chi2uCorrection is big.

Definition at line 199 of file CSCSegAlgoST.h.

Referenced by buildSegments(), and CSCSegAlgoST().

ChamberHitContainer CSCSegAlgoST::protoSegment

private

Definition at line 158 of file CSCSegAlgoST.h.

Referenced by buildSegments(), and prune_bad_hits().

bool CSCSegAlgoST::Pruning

private

Definition at line 172 of file CSCSegAlgoST.h.

Referenced by CSCSegAlgoST(), and run().

const edm::ParameterSet CSCSegAlgoST::ps_

private

Definition at line 116 of file CSCSegAlgoST.h.

Referenced by pset().

std::vector< ChamberHitContainer > CSCSegAlgoST::Psegments

private

Definition at line 125 of file CSCSegAlgoST.h.

Referenced by buildSegments(), and ChooseSegments2().

ChamberHitContainer CSCSegAlgoST::Psegments_hits

private

Definition at line 123 of file CSCSegAlgoST.h.

Referenced by buildSegments().

std::vector< ChamberHitContainer > CSCSegAlgoST::Psegments_noL1

private

Definition at line 127 of file CSCSegAlgoST.h.

Referenced by buildSegments().

std::vector< ChamberHitContainer > CSCSegAlgoST::Psegments_noL2

private

Definition at line 128 of file CSCSegAlgoST.h.

Referenced by buildSegments().

std::vector< ChamberHitContainer > CSCSegAlgoST::Psegments_noL3

private

Definition at line 129 of file CSCSegAlgoST.h.

Referenced by buildSegments().

std::vector< ChamberHitContainer > CSCSegAlgoST::Psegments_noL4

private

Definition at line 130 of file CSCSegAlgoST.h.

Referenced by buildSegments().

std::vector< ChamberHitContainer > CSCSegAlgoST::Psegments_noL5

private

Definition at line 131 of file CSCSegAlgoST.h.

Referenced by buildSegments().

std::vector< ChamberHitContainer > CSCSegAlgoST::Psegments_noL6

private

Definition at line 132 of file CSCSegAlgoST.h.

Referenced by buildSegments().

std::vector< ChamberHitContainer > CSCSegAlgoST::Psegments_noLx

private

Definition at line 126 of file CSCSegAlgoST.h.

Referenced by buildSegments().

CSCSegAlgoShowering* CSCSegAlgoST::showering_

private

Definition at line 117 of file CSCSegAlgoST.h.

Referenced by buildSegments(), CSCSegAlgoST(), and ~CSCSegAlgoST().

const CSCChamber* CSCSegAlgoST::theChamber

private

Definition at line 119 of file CSCSegAlgoST.h.

Referenced by buildSegments(), chamber(), clusterHits(), prune_bad_hits(), and run().

bool CSCSegAlgoST::useShowering

private

Definition at line 176 of file CSCSegAlgoST.h.

Referenced by buildSegments(), and CSCSegAlgoST().

std::vector< float > CSCSegAlgoST::weight_A

private

Definition at line 134 of file CSCSegAlgoST.h.

Referenced by buildSegments(), and ChooseSegments2().

std::vector< float > CSCSegAlgoST::weight_B

private

Definition at line 150 of file CSCSegAlgoST.h.

Referenced by buildSegments().

std::vector< float > CSCSegAlgoST::weight_noL1_A

private

Definition at line 136 of file CSCSegAlgoST.h.

Referenced by buildSegments().

std::vector< float > CSCSegAlgoST::weight_noL1_B

private

Definition at line 151 of file CSCSegAlgoST.h.

Referenced by buildSegments().

std::vector< float > CSCSegAlgoST::weight_noL2_A

private

Definition at line 137 of file CSCSegAlgoST.h.

Referenced by buildSegments().

std::vector< float > CSCSegAlgoST::weight_noL2_B

private

Definition at line 152 of file CSCSegAlgoST.h.

Referenced by buildSegments().

std::vector< float > CSCSegAlgoST::weight_noL3_A

private

Definition at line 138 of file CSCSegAlgoST.h.

Referenced by buildSegments().

std::vector< float > CSCSegAlgoST::weight_noL3_B

private

Definition at line 153 of file CSCSegAlgoST.h.

Referenced by buildSegments().

std::vector< float > CSCSegAlgoST::weight_noL4_A

private

Definition at line 139 of file CSCSegAlgoST.h.

Referenced by buildSegments().

std::vector< float > CSCSegAlgoST::weight_noL4_B

private

Definition at line 154 of file CSCSegAlgoST.h.

Referenced by buildSegments().

std::vector< float > CSCSegAlgoST::weight_noL5_A

private

Definition at line 140 of file CSCSegAlgoST.h.

Referenced by buildSegments().

std::vector< float > CSCSegAlgoST::weight_noL5_B

private

Definition at line 155 of file CSCSegAlgoST.h.

Referenced by buildSegments().

std::vector< float > CSCSegAlgoST::weight_noL6_A

private

Definition at line 141 of file CSCSegAlgoST.h.

Referenced by buildSegments().

std::vector< float > CSCSegAlgoST::weight_noL6_B

private

Definition at line 156 of file CSCSegAlgoST.h.

Referenced by buildSegments().

std::vector< float > CSCSegAlgoST::weight_noLx_A

private

Definition at line 135 of file CSCSegAlgoST.h.

Referenced by buildSegments().

double CSCSegAlgoST::yweightPenalty

private

Definition at line 185 of file CSCSegAlgoST.h.

Referenced by buildSegments(), and CSCSegAlgoST().

double CSCSegAlgoST::yweightPenaltyThreshold

private

Definition at line 184 of file CSCSegAlgoST.h.

Referenced by CSCSegAlgoST(), and run().