CSCXonStrip_MatchGatti.cc

Go to the documentation of this file.

// This is CSCXonStrip_MatchGatti.cc

//---- Large part is copied from RecHitB
//---- author: Stoyan Stoynev - NU

#include <RecoLocalMuon/CSCRecHitD/src/CSCXonStrip_MatchGatti.h>
#include <RecoLocalMuon/CSCRecHitD/src/CSCStripHit.h>

#include <Geometry/CSCGeometry/interface/CSCLayer.h>
#include <Geometry/CSCGeometry/interface/CSCChamberSpecs.h>

#include <CondFormats/CSCObjects/interface/CSCDBGains.h>
#include <CondFormats/DataRecord/interface/CSCDBGainsRcd.h>
#include <CondFormats/CSCObjects/interface/CSCDBCrosstalk.h>
#include <CondFormats/DataRecord/interface/CSCDBCrosstalkRcd.h>
#include <CondFormats/CSCObjects/interface/CSCDBNoiseMatrix.h> 
#include <CondFormats/DataRecord/interface/CSCDBNoiseMatrixRcd.h>


#include <FWCore/MessageLogger/interface/MessageLogger.h>
#include <FWCore/Utilities/interface/Exception.h>


#include "DataFormats/CLHEP/interface/AlgebraicObjects.h"

#include <vector>
#include <cmath>
#include <iostream>
#include <fstream> 
//#include <iomanip.h> 
//#include <iomanip> 
                                                                                                 
#ifndef M_PI_2
#define M_PI_2 1.57079632679489661923
#endif
                                                                                                 

CSCXonStrip_MatchGatti::CSCXonStrip_MatchGatti(const edm::ParameterSet& ps) :
  recoConditions_( 0 )  {

  useCalib                   = ps.getParameter<bool>("CSCUseCalibrations");
  xtalksOffset               = ps.getParameter<double>("CSCStripxtalksOffset");
  noise_level_ME1a                 = ps.getParameter<double>("NoiseLevel_ME1a");
  xt_asymmetry_ME1a                = ps.getParameter<double>("XTasymmetry_ME1a");
  const_syst_ME1a                   = ps.getParameter<double>("ConstSyst_ME1a");
  noise_level_ME1b                 = ps.getParameter<double>("NoiseLevel_ME1b");
  xt_asymmetry_ME1b                = ps.getParameter<double>("XTasymmetry_ME1b");
  const_syst_ME1b                   = ps.getParameter<double>("ConstSyst_ME1b");
  noise_level_ME12                 = ps.getParameter<double>("NoiseLevel_ME12");
  xt_asymmetry_ME12                = ps.getParameter<double>("XTasymmetry_ME12");
  const_syst_ME12                   = ps.getParameter<double>("ConstSyst_ME12");
  noise_level_ME13                 = ps.getParameter<double>("NoiseLevel_ME13");
  xt_asymmetry_ME13                = ps.getParameter<double>("XTasymmetry_ME13");
  const_syst_ME13                   = ps.getParameter<double>("ConstSyst_ME13");
  noise_level_ME21                 = ps.getParameter<double>("NoiseLevel_ME21");
  xt_asymmetry_ME21                = ps.getParameter<double>("XTasymmetry_ME21");
  const_syst_ME21                   = ps.getParameter<double>("ConstSyst_ME21");
  noise_level_ME22                 = ps.getParameter<double>("NoiseLevel_ME22");
  xt_asymmetry_ME22                = ps.getParameter<double>("XTasymmetry_ME22");
  const_syst_ME22                   = ps.getParameter<double>("ConstSyst_ME22");
  noise_level_ME31                 = ps.getParameter<double>("NoiseLevel_ME31");
  xt_asymmetry_ME31                = ps.getParameter<double>("XTasymmetry_ME31");
  const_syst_ME31                   = ps.getParameter<double>("ConstSyst_ME31");
  noise_level_ME32                 = ps.getParameter<double>("NoiseLevel_ME32");
  xt_asymmetry_ME32                = ps.getParameter<double>("XTasymmetry_ME32");
  const_syst_ME32                   = ps.getParameter<double>("ConstSyst_ME32");
  noise_level_ME41                 = ps.getParameter<double>("NoiseLevel_ME41");
  xt_asymmetry_ME41                = ps.getParameter<double>("XTasymmetry_ME41");
  const_syst_ME41                   = ps.getParameter<double>("ConstSyst_ME41");
 
  getCorrectionValues("StringCurrentlyNotUsed");
}


CSCXonStrip_MatchGatti::~CSCXonStrip_MatchGatti(){
}


/* findPosition
 *
 */
void CSCXonStrip_MatchGatti::findXOnStrip( const CSCDetId& id, const CSCLayer* layer, const CSCStripHit& stripHit, 
                                            int centralStrip, float& xWithinChamber, float& sWidth, 
                                            const float& tpeak, float& xWithinStrip, float& sigma, int & quality_flag) {
  quality_flag = 0; 
  // Initialize Gatti parameters using chamberSpecs
  // Cache specs_ info for ease of access
  specs_ = layer->chamber()->specs();
  stripWidth = sWidth;
  //initChamberSpecs();
  // Initialize output parameters  
  xWithinStrip = xWithinChamber;  

  CSCStripHit::ChannelContainer const & strips = stripHit.strips();
  int nStrips = strips.size();
  int centStrip = nStrips/2 + 1;   
  std::vector<float> const & adcs  = stripHit.s_adc();
  int tmax = stripHit.tmax();

  //float t_peak = tpeak;
  //float t_zero = 0.;
  //float adc[4];
  //
  //if ( useCalib ) {
  //
  //  for ( int t = 0; t < 4; ++t ) {
  //    int k  = t + 4 * (centStrip-1);
  //    adc[t] = adcs[k];
  //  }
  //
  //  // t_peak from peak finder is now 'absolute' i.e. in ns from start of sca time bin 0
  //  t_peak = peakTimeFinder_->peakTime( tmax, adc, t_peak );
  //  // Just for completeness, the start time of the pulse is 133 ns earlier, according to Stan :)
  //  t_zero = t_peak - 133.;
  //  // and reset tpeak since that's the way it gets passed out of this function (Argh!)
  //  tpeak = t_peak;
  //  LogTrace("CSCRecHit|CSCXonStrip_MatchGatti") << "CSCXonStrip_MatchGatti: " << 
  //     id << " strip=" << centralStrip << ", t_zero=" << t_zero << ", t_peak=" << t_peak;
  //}
      
  //---- fill the charge matrix (3x3)
  float adc[4];
  int j = 0;
  for ( int i = 1; i <= nStrips; ++i ) {
    if ( i > (centStrip-2) && i < (centStrip+2) ) {
      std::vector<float> adcsFit;
      for ( int t = 0; t < 4; ++t ) {
        int k  = t + 4*(i-1);
        adc[t] = adcs[k];
        if ( t < 3) chargeSignal[j][t] = adc[t];
      }
      j++;
    }
  }


  // Load in x-talks:

  if ( useCalib ) {
    std::vector<float> xtalks;
    recoConditions_->crossTalk( id, centralStrip, xtalks );
    float dt = 50.f * tmax - tpeak;
    // XTalks; l,r are for left, right XTalk; lr0,1,2 are for what charge "remains" in the strip 
    for ( int t = 0; t < 3; ++t ) {
      xt_l[0][t] = xtalks[0] * (50.f* (t-1) + dt) + xtalks[1] + xtalksOffset;
      xt_r[0][t] = xtalks[2] * (50.f* (t-1) + dt) + xtalks[3] + xtalksOffset;
      xt_l[1][t] = xtalks[4] * (50.f* (t-1) + dt) + xtalks[5] + xtalksOffset;
      xt_r[1][t] = xtalks[6] * (50.f* (t-1) + dt) + xtalks[7] + xtalksOffset;
      xt_l[2][t] = xtalks[8] * (50.f* (t-1) + dt) + xtalks[9] + xtalksOffset;
      xt_r[2][t] = xtalks[10]* (50.f* (t-1) + dt) + xtalks[11]+ xtalksOffset;

      xt_lr0[t] = (1. - xt_l[0][t] - xt_r[0][t]);
      xt_lr1[t] = (1. - xt_l[1][t] - xt_r[1][t]);
      xt_lr2[t] = (1. - xt_l[2][t] - xt_r[2][t]);
    }
  } else { 
    for ( int t = 0; t < 3; ++t ) {
      xt_l[0][t] = xtalksOffset;
      xt_r[0][t] = xtalksOffset;
      xt_l[1][t] = xtalksOffset; 
      xt_r[1][t] = xtalksOffset; 
      xt_l[2][t] = xtalksOffset; 
      xt_r[2][t] = xtalksOffset; 

      xt_lr0[t] = (1. - xt_l[0][t] - xt_r[0][t]);
      xt_lr1[t] = (1. - xt_l[1][t] - xt_r[1][t]);
      xt_lr2[t] = (1. - xt_l[2][t] - xt_r[2][t]);
    } 
  }
  

  // vector containing noise starts at tmax - 1, and tmax > 3, but....
  int tbin = tmax - 4;

  // .... originally, suppose to have tmax in tbin 4 or 5, but noticed in MTCC lots of 
  // hits with tmax == 3, so let's allow these, and shift down noise matrix by one element...
  // This is a patch because the calibration database doesn't have elements for tbin = 2, 
  // e.g. there is no element e[tmax-1,tmax+1] = e[2,4].

  if (tmax < 4) tbin = 0;    // patch

  // Load in auto-correlation noise matrices
  if ( useCalib ) {
    std::vector<float> nmatrix;
    recoConditions_->noiseMatrix( id, centralStrip, nmatrix );
    for ( int istrip =0; istrip < 3; ++istrip ) {
      a11[istrip] = nmatrix[0+tbin*3+istrip*15];
      a12[istrip] = nmatrix[1+tbin*3+istrip*15];
      a13[istrip] = nmatrix[2+tbin*3+istrip*15];
      a22[istrip] = nmatrix[3+tbin*3+istrip*15];
      a23[istrip] = nmatrix[4+tbin*3+istrip*15];
      a33[istrip] = nmatrix[6+tbin*3+istrip*15];
    }
  } else {
    // FIXME:  NO HARD WIRED VALUES !!!
    for ( int istrip =0; istrip < 3; ++istrip ) {
      a11[istrip] = 10.0;
      a12[istrip] = 0.0;
      a13[istrip] = 0.0;
      a22[istrip] = 10.0;
      a23[istrip] = 0.0;
      a33[istrip] = 10.0;
    }
  }
  
  //---- Set up noise, XTalk matrices 
  setupMatrix();
  //---- Calculate the coordinate within the strip and associate uncertainty

  static const std::string ME1a("ME1/a");
  static const std::string ME1b("ME1/b");

  bool ME1_1 = (ME1a == specs_->chamberTypeName() || ME1b == specs_->chamberTypeName());


  // due to cross talks and 3 time bin sum it is in principe possible that the center strip is not the maximum strip
  // in some cases the consequences could be quite extreme
  // take some measures against the extreme cases
  bool peakMismatch = false;
  std::vector <float> charges(3);
  charges[0] = q_sumL;
  charges[1] = q_sumC;
  charges[2] = q_sumR;
  int min_index = min_element(charges.begin(),charges.end()) - charges.begin();
  int max_index = max_element(charges.begin(),charges.end()) - charges.begin();
  if(1!=max_index && (1==min_index ||
              // the condition below means that if the initial position estimate within strip (|xF|)
              // is above  1.1/2 = 0.55 "strip widths" peakMismatch is set to true (so special case);
              // in normal cases |xF|<=0.5 (0.5 is at the edge of the "central" strip)
              (charges[max_index] - charges[min_index]) > 1.1*(q_sumC - charges[min_index]) )){
    peakMismatch = true;
    switch (max_index){
    case 0:
      xWithinStrip = -1;
      break;
    case 2:
      xWithinStrip = 1;
      break;
    default:
      // should be an error message here
      xWithinStrip = 0;// in case?
      break;
    }
  }
  // we don't have the needed information (it's similar to the "edge" strip case)  
  //else if(stripHit.isNearDeadStrip()){
  else if(stripHit.deadStrip()>0){
    xWithinStrip = 0;
  }
  else{
    //
    xWithinStrip = float(calculateXonStripPosition(stripWidth, ME1_1));
  }
  xWithinChamber = xWithinChamber + (xWithinStrip * stripWidth);
  if(peakMismatch){
    sigma = stripWidth/std::sqrt(12.f);
  }
  else{
    
    //---- error estimation
    int factorStripWidth = int( std::sqrt(stripWidth/0.38f) );
    int maxConsecutiveStrips = 8;
    if(factorStripWidth){
      maxConsecutiveStrips /=  factorStripWidth ;
    }
    maxConsecutiveStrips++;

    struct ChamberTypes {
    std::map <std::string, int> chamberTypes;
      int  operator()(std::string const & s) const { auto p = chamberTypes.find(s);  return p!=chamberTypes.end() ? (*p).second : 0; }
      ChamberTypes() {
    chamberTypes["ME1/a"] = 1;
    chamberTypes["ME1/b"] = 2;
    chamberTypes["ME1/2"] = 3;
    chamberTypes["ME1/3"] = 4;
    chamberTypes["ME2/1"] = 5;
    chamberTypes["ME2/2"] = 6;
    chamberTypes["ME3/1"] = 7;
    chamberTypes["ME3/2"] = 8;
    chamberTypes["ME4/1"] = 9;
    chamberTypes["ME4/2"] = 8;
      }
    };
    static const ChamberTypes chamberTypes;
    
    switch(chamberTypes(specs_->chamberTypeName())){
    case 1:
      noise_level  = noise_level_ME1a;
      xt_asymmetry = xt_asymmetry_ME1a;
      const_syst = const_syst_ME1a;
      break;
      
    case 2:
      noise_level  = noise_level_ME1b;
      xt_asymmetry = xt_asymmetry_ME1b;
      const_syst = const_syst_ME1b;
      
    case 3:
      noise_level  = noise_level_ME12;
      xt_asymmetry = xt_asymmetry_ME12;
      const_syst = const_syst_ME12;
      break;
      
    case 4:
      noise_level  = noise_level_ME13;
      xt_asymmetry = xt_asymmetry_ME13;
      const_syst = const_syst_ME13;
      break;
      
    case 5:
      noise_level  = noise_level_ME21;
      xt_asymmetry = xt_asymmetry_ME21;
      const_syst = const_syst_ME21;
      break;
      
    case 6:
      noise_level  = noise_level_ME22;
      xt_asymmetry = xt_asymmetry_ME22;
      const_syst = const_syst_ME22;
      break;

    case 7:
      noise_level  = noise_level_ME31;
      xt_asymmetry = xt_asymmetry_ME31;
      const_syst = const_syst_ME31;
      break;

    case 8:
      noise_level  = noise_level_ME32;
      xt_asymmetry = xt_asymmetry_ME32;
      const_syst = const_syst_ME32;
      break;

    case 9:
      noise_level  = noise_level_ME41;
      xt_asymmetry = xt_asymmetry_ME41;
      const_syst = const_syst_ME41;
      break;

    default:
      noise_level  = noise_level_ME22;
      xt_asymmetry = xt_asymmetry_ME22;
      const_syst = const_syst_ME22;

    }
    if(0==stripHit.deadStrip() &&
       stripHit.numberOfConsecutiveStrips()<maxConsecutiveStrips &&
       std::abs(stripHit.closestMaximum())>maxConsecutiveStrips/2 ){
      sigma =  float(calculateXonStripError(stripWidth, ME1_1));
    }
    else{ //---- near dead strip or too close maxima or too wide strip cluster
      sigma = stripWidth/std::sqrt(12.f);
    }
  }
  quality_flag = 1;
}

/* setupMatrix
 *
 */
void CSCXonStrip_MatchGatti::setupMatrix() {
  //---- a??? and v??[] could be skipped for now...; not used yet

  /*
  double dd, a11t, a12t, a13t, a22t, a23t, a33t;
  double syserr = adcSystematics;
  double xtlk_err = xtalksSystematics;
  // Left strip
  a11t = a11[0] + syserr*syserr * ChargeSignal[0][0]*ChargeSignal[0][0] + xtlk_err*xtlk_err*ChargeSignal[1][0]*ChargeSignal[1][0];
  a12t = a12[0] + syserr*syserr * ChargeSignal[0][0]*ChargeSignal[0][1];
  a13t = a13[0] + syserr*syserr * ChargeSignal[0][0]*ChargeSignal[0][2];
  a22t = a22[0] + syserr*syserr * ChargeSignal[0][1]*ChargeSignal[0][1] + xtlk_err*xtlk_err*ChargeSignal[1][1]*ChargeSignal[1][1];
  a23t = a23[0] + syserr*syserr * ChargeSignal[0][1]*ChargeSignal[0][2];
  a33t = a33[0] + syserr*syserr * ChargeSignal[0][2]*ChargeSignal[0][2] + xtlk_err*xtlk_err*ChargeSignal[1][2]*ChargeSignal[1][2];

  dd     = (a11t*a33t*a22t - a11t*a23t*a23t - a33t*a12t*a12t 
                       + 2.* a12t*a13t*a23t - a13t*a13t*a22t );

  v11[0] = (a33t*a22t - a23t*a23t)/dd;
  v12[0] =-(a33t*a12t - a13t*a23t)/dd;
  v13[0] = (a12t*a23t - a13t*a22t)/dd;
  v22[0] = (a33t*a11t - a13t*a13t)/dd;
  v23[0] =-(a23t*a11t - a12t*a13t)/dd;
  v33[0] = (a22t*a11t - a12t*a12t)/dd;
     
  // Center strip
  a11t = a11[1] + syserr*syserr * ChargeSignal[1][0]*ChargeSignal[1][0] + xtlk_err*xtlk_err*(ChargeSignal[0][0]*ChargeSignal[0][0]+ChargeSignal[2][0]*ChargeSignal[2][0]);
  a12t = a12[1] + syserr*syserr * ChargeSignal[1][0]*ChargeSignal[1][1];
  a13t = a13[1] + syserr*syserr * ChargeSignal[1][0]*ChargeSignal[1][2];
  a22t = a22[1] + syserr*syserr * ChargeSignal[1][1]*ChargeSignal[1][1] + xtlk_err*xtlk_err*(ChargeSignal[0][1]*ChargeSignal[0][1]+ChargeSignal[2][1]*ChargeSignal[2][1]);
  a23t = a23[1] + syserr*syserr * ChargeSignal[1][1]*ChargeSignal[1][2];
  a33t = a33[1] + syserr*syserr * ChargeSignal[1][2]*ChargeSignal[1][2] + xtlk_err*xtlk_err*(ChargeSignal[0][2]*ChargeSignal[0][2]+ChargeSignal[2][2]*ChargeSignal[2][2]);

  dd     = (a11t*a33t*a22t - a11t*a23t*a23t - a33t*a12t*a12t
                       + 2.* a12t*a13t*a23t - a13t*a13t*a22t );

  v11[1] = (a33t*a22t - a23t*a23t)/dd;
  v12[1] =-(a33t*a12t - a13t*a23t)/dd;
  v13[1] = (a12t*a23t - a13t*a22t)/dd;
  v22[1] = (a33t*a11t - a13t*a13t)/dd;
  v23[1] =-(a23t*a11t - a12t*a13t)/dd;
  v33[1] = (a22t*a11t - a12t*a12t)/dd;

  // Right strip
  a11t = a11[2] + syserr*syserr * ChargeSignal[2][0]*ChargeSignal[2][0] + xtlk_err*xtlk_err*ChargeSignal[1][0]*ChargeSignal[1][0];
  a12t = a12[2] + syserr*syserr * ChargeSignal[2][0]*ChargeSignal[2][1];
  a13t = a13[2] + syserr*syserr * ChargeSignal[2][0]*ChargeSignal[2][2];
  a22t = a22[2] + syserr*syserr * ChargeSignal[2][1]*ChargeSignal[2][1] + xtlk_err*xtlk_err*ChargeSignal[1][1]*ChargeSignal[1][1];
  a23t = a23[2] + syserr*syserr * ChargeSignal[2][1]*ChargeSignal[2][2];
  a33t = a33[2] + syserr*syserr * ChargeSignal[2][2]*ChargeSignal[2][2] + xtlk_err*xtlk_err*ChargeSignal[1][2]*ChargeSignal[1][2];

  dd     = (a11t*a33t*a22t - a11t*a23t*a23t - a33t*a12t*a12t
                        +2.* a12t*a13t*a23t - a13t*a13t*a22t );

  v11[2] = (a33t*a22t - a23t*a23t)/dd;
  v12[2] =-(a33t*a12t - a13t*a23t)/dd;
  v13[2] = (a12t*a23t - a13t*a22t)/dd;
  v22[2] = (a33t*a11t - a13t*a13t)/dd;
  v23[2] =-(a23t*a11t - a12t*a13t)/dd;
  v33[2] = (a22t*a11t - a12t*a12t)/dd;
*/
  //---- Find the inverted XTalk matrix and apply it to the charge (3x3)
  //---- Thus the charge before the XTalk is obtained
  CLHEP::HepMatrix cross_talks_inv(3,3);
  int err = 0;
  //---- q_sum is 3 time bins summed; L, C, R - left, central, right strips
  q_sum = q_sumL = q_sumC = q_sumR = 0.;
  double charge = 0.;
  for(int iTime=0;iTime<3;iTime++){
    cross_talks_inv(1,1) = xt_lr0[iTime];
    cross_talks_inv(1,2) = xt_l[1][iTime];
    cross_talks_inv(1,3) = 0.;
    cross_talks_inv(2,1) =  xt_r[0][iTime];
    cross_talks_inv(2,2) = xt_lr1[iTime];
    cross_talks_inv(2,3)  = xt_l[2][iTime];
    cross_talks_inv(3,1) = 0.;
    cross_talks_inv(3,2) = xt_r[1][iTime];
    cross_talks_inv(3,3) = xt_lr2[iTime];
    cross_talks_inv.invert(err);
    if (err != 0) {
      edm::LogWarning("FailedXTalkiInversionNoCrosstalkCorrection") <<"Failed to invert XTalks matrix. No cross-talk correction for this rechit."; 
      //edm::LogError("CSCRecHit") << "Failed to invert XTalks matrix. No cross-talk correction for this rechit.";
      return;
    }
    //---- "charge" is XT-corrected charge
    charge = chargeSignal[0][iTime]*cross_talks_inv(1,1) + chargeSignal[1][iTime]*cross_talks_inv(1,2) + 
      chargeSignal[2][iTime]*cross_talks_inv(1,3);
    //---- Negative charge? According to studies (and logic) - better use 0 charge
    //----- Later studies suggest that this only do harm. I am still worried about 
    // charges of -50 ADC and below (0.5% of the cases) but let see 
    //if(charge<0.){
      //charge = 0.;
    //}
    q_sum+=charge;
    q_sumL+=charge;
    charge = chargeSignal[0][iTime]*cross_talks_inv(2,1) + chargeSignal[1][iTime]*cross_talks_inv(2,2) + 
      chargeSignal[2][iTime]* cross_talks_inv(2,3);
    //if(charge<0.){
      //charge = 0.;
    //}
    q_sum+=charge;
    q_sumC+=charge;
    charge = chargeSignal[0][iTime]*cross_talks_inv(3,1) + chargeSignal[1][iTime]*cross_talks_inv(3,2) + 
      chargeSignal[2][iTime]*cross_talks_inv(3,3);
    //if(charge<0.){
      //charge = 0.;
    //}
    q_sum+=charge;
    q_sumR+=charge;
  }
}


/* initChamberSpecs
 *
 */
void CSCXonStrip_MatchGatti::initChamberSpecs() {
  // Not used directly but these are parameters used for extracting the correction values
  // in coordinate and error estimators

  // Distance between anode and cathode
  h = specs_->anodeCathodeSpacing();
  r = h / stripWidth;

  // Wire spacing
  double wspace = specs_->wireSpacing();

  // Wire radius
  double wradius = specs_->wireRadius();

  // Accepted parameters in Gatti function
  const double parm[5] = {.1989337e-02, -.6901542e-04, .8665786, 154.6177, -.680163e-03 };

  k_3 = ( parm[0]*wspace/h + parm[1] )
      * ( parm[2]*wspace/wradius + parm[3] + parm[4]*(wspace/wradius)*(wspace/wradius) );

  sqrt_k_3 = std::sqrt( k_3 );
  norm     = r * (0.5 / std::atan( sqrt_k_3 )); // changed from norm to r * norm
  k_2      = M_PI_2 * ( 1. - sqrt_k_3 /2. );
  k_1      = 0.25 * k_2 * sqrt_k_3 / std::atan( sqrt_k_3 );
}


void CSCXonStrip_MatchGatti::getCorrectionValues(std::string estimator){
  hardcodedCorrectionInitialization();
}

double CSCXonStrip_MatchGatti::estimated2GattiCorrection(double x_estimated, float stripWidth, bool ME1_1) {
  //---- 11 "nominal" strip widths : 0.6 - 1.6 cm; for ME1_1 just 6 "nominal" strip widths : 0.3 - 0.8 cm; see HardCodedCorrectionInitialization() 
  //---- Calculate corrections at specific  Xestimated (linear interpolation between points)
  int n_SW;
  int min_SW;
  if(ME1_1){
    n_SW = n_SW_ME1_1;
    min_SW = 3; // min SW calculated is 0.3 cm
  }
  else{
    n_SW = n_SW_noME1_1;
    min_SW = 6;// min SW calculated is 0.6 cm
  }
  int stripDown = int(10.*stripWidth) - min_SW; // 0 is at min strip width calculated
  int stripUp = stripDown + 1;
  if(stripUp>n_SW-1){
    //---- to be checked...
    //if(stripUp>n_SW){
      //std::cout<<" Is strip width = "<<stripWidth<<" OK?" <<std::endl;
    //}
    stripUp = n_SW-1;
  }

  double half_strip_width = 0.5;
  //const int Nbins = 501;
  const int n_bins = n_val;
  double corr_2_xestim = 999.;
  if(stripDown<0){
    corr_2_xestim = 1;
  }
  else{
    //---- Parametrized x_gatti minus x_estimated differences

    int xestim_bin = -999;
    double delta_strip_width = 999.;
    double delta_strip_widthUpDown = 999.;
    double diff_2_strip_width = 999.;
    delta_strip_width = stripWidth - int(stripWidth*10)/10.;
    delta_strip_widthUpDown = 0.1;

    if(std::abs(x_estimated)>0.5){
      if(std::abs(x_estimated)>1.){
        corr_2_xestim = 1.;// for now; to be investigated
      }
      else{  
    //if(std::abs(Xestimated)>0.55){
      //std::cout<<"X position from the estimated position above 0.55 (safty margin)?! "<<std::endl;
      //CorrToXc = 999.;
    //}
    xestim_bin = int((1.- std::abs(x_estimated))/half_strip_width * n_bins);
    if(ME1_1){
      diff_2_strip_width = x_correction_ME1_1[stripUp][xestim_bin]-x_correction_ME1_1[stripDown][xestim_bin];
      corr_2_xestim =  x_correction_ME1_1[stripDown][xestim_bin] +
        (delta_strip_width/delta_strip_widthUpDown)*diff_2_strip_width ;
    }
    else{
      diff_2_strip_width = x_correction_noME1_1[stripUp][xestim_bin]-x_correction_noME1_1[stripDown][xestim_bin];
      corr_2_xestim =  x_correction_noME1_1[stripDown][xestim_bin] +
        (delta_strip_width/delta_strip_widthUpDown)*diff_2_strip_width ;
    }
    corr_2_xestim = -corr_2_xestim;
      }
    }
    else{
      xestim_bin = int((std::abs(x_estimated)/half_strip_width) * n_bins);
      if(ME1_1){
    diff_2_strip_width = x_correction_ME1_1[stripUp][xestim_bin] - x_correction_ME1_1[stripDown][xestim_bin];
    corr_2_xestim =  x_correction_ME1_1[stripDown][xestim_bin] +
      (delta_strip_width/delta_strip_widthUpDown)*diff_2_strip_width ;
      }
      else{
    diff_2_strip_width = x_correction_noME1_1[stripUp][xestim_bin] - x_correction_noME1_1[stripDown][xestim_bin];
    corr_2_xestim =  x_correction_noME1_1[stripDown][xestim_bin] +
      (delta_strip_width/delta_strip_widthUpDown)*diff_2_strip_width ;
      }
    }
    if(x_estimated<0.){
      corr_2_xestim = -corr_2_xestim;
    }
  }
  
  return corr_2_xestim;
}


double CSCXonStrip_MatchGatti::estimated2Gatti(double x_estimated, float stripWidth, bool ME1_1) {

 double x_corr = estimated2GattiCorrection(x_estimated, stripWidth, ME1_1);
 double x_gatti = x_estimated + x_corr;

 return x_gatti;
}

double CSCXonStrip_MatchGatti::xfError_Noise(double noise){

  double min, max;
  if(q_sumR>q_sumL){
    min = q_sumL;
    max = q_sumR;
  }
  else{
    min = q_sumR;
    max = q_sumL;
  }
  //---- Error propagation...
  //---- Names here are fake! Due to technical features
  double dr_L2 = pow(q_sumR-q_sumL,2);
  double dr_C2 = pow(q_sumC-min,2);
  double dr_R2 = pow(q_sumC-max,2);
  double error = std::sqrt(dr_L2 + dr_C2 + dr_R2)*noise/std::pow(q_sumC-min,2)/2;

  return error;
}

double CSCXonStrip_MatchGatti::xfError_XTasym(double xt_asym){

  double min;
  if(q_sumR>q_sumL){
    min = q_sumL;
  }
  else{
    min = q_sumR;
  }
  //---- Error propagation
  double dXTL = (std::pow(q_sumC,2)+std::pow(q_sumR,2)-q_sumL*q_sumR-q_sumR*q_sumC);
  double dXTR = (std::pow(q_sumC,2)+std::pow(q_sumL,2)-q_sumL*q_sumR-q_sumL*q_sumC);
  double dXT = std::sqrt(std::pow(dXTL,2) + std::pow(dXTR,2))/std::pow((q_sumC-min),2)/2;
  double error = dXT * xt_asym;

  return error;
}


double CSCXonStrip_MatchGatti::calculateXonStripError(float stripWidth, bool ME1_1){
  double min;
  if(q_sumR>q_sumL){
    min = q_sumL;
  }
  else{
    min = q_sumR;
  }
  
  double xf = (q_sumR - q_sumL)/(q_sumC - min)/2;
  double xf_ErrorNoise = xfError_Noise(noise_level);
  double xf_ErrorXTasym = xfError_XTasym(xt_asymmetry);
  // x_G = x_F + correction_functon(x_F)
  // as these are correlated the error should be simply d(x_G) = |d(x_F)| + [correction_functon(x_F+|d(x_F)|) - correction_functon(x_F)]
  double d_xf = std::sqrt( std::pow( xf_ErrorNoise, 2) + std::pow( xf_ErrorXTasym, 2));
  double d_corr = estimated2GattiCorrection(xf+d_xf, stripWidth, ME1_1) - estimated2GattiCorrection(xf, stripWidth, ME1_1);
  double x_shift = d_xf + d_corr;
  //  double x_shift = sqrt( pow( xf_ErrorNoise, 2) + pow( xf_ErrorXTasym, 2)) * 
  //(1 + (estimated2GattiCorrection(xf+0.001, stripWidth, ME1_1) -
  //  estimated2GattiCorrection(xf, stripWidth, ME1_1))*1000.);
  double x_error =   std::sqrt( std::pow( std::abs(x_shift)*stripWidth, 2) + std::pow(const_syst, 2) );
  return  x_error; 
}

double CSCXonStrip_MatchGatti::calculateXonStripPosition(float stripWidth, bool ME1_1){ 

  double  x_estimated = -99.;
  double min;
  if(q_sumR>q_sumL){
    min = q_sumL;
  }
  else{
    min = q_sumR;
  }
  //---- This is XF ( X Florida - after the first group that used it)  
  x_estimated = (q_sumR - q_sumL)/(q_sumC - min)/2;
  double x_gatti = estimated2Gatti(x_estimated, stripWidth, ME1_1);
  return x_gatti;
}

// Define space for statics
const int CSCXonStrip_MatchGatti::n_val;
const int CSCXonStrip_MatchGatti::n_SW_noME1_1;
const int CSCXonStrip_MatchGatti::n_SW_ME1_1;