CSCSegAlgoST.cc

Go to the documentation of this file.

#include "CSCSegAlgoST.h"
#include "CSCCondSegFit.h"
#include "CSCSegAlgoShowering.h"

#include "DataFormats/GeometryVector/interface/GlobalPoint.h"

#include "Geometry/CSCGeometry/interface/CSCLayer.h"

#include "FWCore/MessageLogger/interface/MessageLogger.h"

#include "CommonTools/Statistics/interface/ChiSquaredProbability.h"

#include <algorithm>
#include <cmath>
#include <iostream>
#include <string>

/* constructor
 *
 */
CSCSegAlgoST::CSCSegAlgoST(const edm::ParameterSet& ps)
    : 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";
}

/* Destructor
 *
 */
CSCSegAlgoST::~CSCSegAlgoST() { delete showering_; }

std::vector<CSCSegment> CSCSegAlgoST::run(const CSCChamber* aChamber, const ChamberHitContainer& rechits) {
  // 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);
  }
}

// ********************************************************************;
// *** This method is meant to remove clearly bad hits, using as    ***;
// *** much information from the chamber as possible (e.g. charge,  ***;
// *** hit position, timing, etc.)                                  ***;
// ********************************************************************;
std::vector<CSCSegment> CSCSegAlgoST::prune_bad_hits(const CSCChamber* aChamber, std::vector<CSCSegment>& segments) {
  //   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;
}

// ********************************************************************;
std::vector<std::vector<const CSCRecHit2D*> > CSCSegAlgoST::clusterHits(const CSCChamber* aChamber,
                                                                        const ChamberHitContainer& rechits) {
  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;
}

std::vector<std::vector<const CSCRecHit2D*> > CSCSegAlgoST::chainHits(const CSCChamber* aChamber,
                                                                      const ChamberHitContainer& rechits) {
  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;
}

bool CSCSegAlgoST::isGoodToMerge(bool gangedME11a, ChamberHitContainer& newChain, ChamberHitContainer& oldChain) {
  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;
}

double CSCSegAlgoST::theWeight(
    double coordinate_1, double coordinate_2, double coordinate_3, float layer_1, float layer_2, float layer_3) {
  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;
}

/* 
 * This algorithm uses a Minimum Spanning Tree (ST) approach to build
 * endcap muon track segments from the rechits in the 6 layers of a CSC.
 */
std::vector<CSCSegment> CSCSegAlgoST::buildSegments(const ChamberHitContainer& rechits) {
  // 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;
}

void CSCSegAlgoST::ChooseSegments2a(std::vector<ChamberHitContainer>& chosen_segments, int chosen_seg) {
  // just return best segment
  GoodSegments.clear();
  GoodSegments.push_back(chosen_segments[chosen_seg]);
}

void CSCSegAlgoST::ChooseSegments3(std::vector<ChamberHitContainer>& chosen_segments,
                                   std::vector<float>& chosen_weight,
                                   int chosen_seg) {
  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;
  }
}

void CSCSegAlgoST::ChooseSegments2(int best_seg) {
  //  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::findDuplicates(std::vector<CSCSegment>& segments) {
  // 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);
  }
}

void CSCSegAlgoST::dumpSegment(const CSCSegment& seg) const {
  // 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();
}