CMSSW Reference Manual

00001 // This is CSCXonStrip_MatchGatti.cc
00002 
00003 //---- Large part is copied from RecHitB
00004 //---- author: Stoyan Stoynev - NU
00005 
00006 #include <RecoLocalMuon/CSCRecHitD/src/CSCXonStrip_MatchGatti.h>
00007 //#include <RecoLocalMuon/CSCRecHitD/src/CSCStripCrosstalk.h>
00008 //#include <RecoLocalMuon/CSCRecHitD/src/CSCStripNoiseMatrix.h>
00009 #include <RecoLocalMuon/CSCRecHitD/src/CSCFindPeakTime.h>
00010 #include <RecoLocalMuon/CSCRecHitD/src/CSCStripHit.h>
00011 
00012 #include <Geometry/CSCGeometry/interface/CSCLayer.h>
00013 #include <Geometry/CSCGeometry/interface/CSCChamberSpecs.h>
00014 
00015 #include <CondFormats/CSCObjects/interface/CSCDBGains.h>
00016 #include <CondFormats/DataRecord/interface/CSCDBGainsRcd.h>
00017 #include <CondFormats/CSCObjects/interface/CSCDBCrosstalk.h>
00018 #include <CondFormats/DataRecord/interface/CSCDBCrosstalkRcd.h>
00019 #include <CondFormats/CSCObjects/interface/CSCDBNoiseMatrix.h> 
00020 #include <CondFormats/DataRecord/interface/CSCDBNoiseMatrixRcd.h>
00021 
00022 
00023 #include <FWCore/MessageLogger/interface/MessageLogger.h>
00024 #include <FWCore/Utilities/interface/Exception.h>
00025 
00026 
00027 #include "DataFormats/CLHEP/interface/AlgebraicObjects.h"
00028 
00029 #include <vector>
00030 #include <cmath>
00031 #include <iostream>
00032 #include <fstream> 
00033 //#include <iomanip.h> 
00034 //#include <iomanip> 
00035                                                                                                  
00036 #ifndef M_PI_2
00037 #define M_PI_2 1.57079632679489661923
00038 #endif
00039                                                                                                  
00040 
00041 CSCXonStrip_MatchGatti::CSCXonStrip_MatchGatti(const edm::ParameterSet& ps) :
00042    recoConditions_( 0 ){
00043 
00044   useCalib                   = ps.getUntrackedParameter<bool>("CSCUseCalibrations");
00045   xtalksOffset               = ps.getUntrackedParameter<double>("CSCStripxtalksOffset");
00046   noise_level_ME1a                 = ps.getUntrackedParameter<double>("NoiseLevel_ME1a");
00047   xt_asymmetry_ME1a                = ps.getUntrackedParameter<double>("XTasymmetry_ME1a");
00048   const_syst_ME1a                   = ps.getUntrackedParameter<double>("ConstSyst_ME1a");
00049   noise_level_ME1b                 = ps.getUntrackedParameter<double>("NoiseLevel_ME1b");
00050   xt_asymmetry_ME1b                = ps.getUntrackedParameter<double>("XTasymmetry_ME1b");
00051   const_syst_ME1b                   = ps.getUntrackedParameter<double>("ConstSyst_ME1b");
00052   noise_level_ME12                 = ps.getUntrackedParameter<double>("NoiseLevel_ME12");
00053   xt_asymmetry_ME12                = ps.getUntrackedParameter<double>("XTasymmetry_ME12");
00054   const_syst_ME12                   = ps.getUntrackedParameter<double>("ConstSyst_ME12");
00055   noise_level_ME13                 = ps.getUntrackedParameter<double>("NoiseLevel_ME13");
00056   xt_asymmetry_ME13                = ps.getUntrackedParameter<double>("XTasymmetry_ME13");
00057   const_syst_ME13                   = ps.getUntrackedParameter<double>("ConstSyst_ME13");
00058   noise_level_ME21                 = ps.getUntrackedParameter<double>("NoiseLevel_ME21");
00059   xt_asymmetry_ME21                = ps.getUntrackedParameter<double>("XTasymmetry_ME21");
00060   const_syst_ME21                   = ps.getUntrackedParameter<double>("ConstSyst_ME21");
00061   noise_level_ME22                 = ps.getUntrackedParameter<double>("NoiseLevel_ME22");
00062   xt_asymmetry_ME22                = ps.getUntrackedParameter<double>("XTasymmetry_ME22");
00063   const_syst_ME22                   = ps.getUntrackedParameter<double>("ConstSyst_ME22");
00064   noise_level_ME31                 = ps.getUntrackedParameter<double>("NoiseLevel_ME31");
00065   xt_asymmetry_ME31                = ps.getUntrackedParameter<double>("XTasymmetry_ME31");
00066   const_syst_ME31                   = ps.getUntrackedParameter<double>("ConstSyst_ME31");
00067   noise_level_ME32                 = ps.getUntrackedParameter<double>("NoiseLevel_ME32");
00068   xt_asymmetry_ME32                = ps.getUntrackedParameter<double>("XTasymmetry_ME32");
00069   const_syst_ME32                   = ps.getUntrackedParameter<double>("ConstSyst_ME32");
00070   noise_level_ME41                 = ps.getUntrackedParameter<double>("NoiseLevel_ME41");
00071   xt_asymmetry_ME41                = ps.getUntrackedParameter<double>("XTasymmetry_ME41");
00072   const_syst_ME41                   = ps.getUntrackedParameter<double>("ConstSyst_ME41");
00073   peakTimeFinder_            = new CSCFindPeakTime();
00074   getCorrectionValues("StringCurrentlyNotUsed");
00075 }
00076 
00077 
00078 CSCXonStrip_MatchGatti::~CSCXonStrip_MatchGatti(){
00079   delete peakTimeFinder_;
00080 }
00081 
00082 
00083 /* findPosition
00084  *
00085  */
00086 void CSCXonStrip_MatchGatti::findXOnStrip( const CSCDetId& id, const CSCLayer* layer, const CSCStripHit& stripHit, 
00087                                             int centralStrip, float& xWithinChamber, float& sWidth, 
00088                                             float& tpeak,  float& xWithinStrip, float& sigma, int & quality_flag) {
00089   quality_flag = 0; 
00090   // Initialize Gatti parameters using chamberSpecs
00091   // Cache specs_ info for ease of access
00092   specs_ = layer->chamber()->specs();
00093   stripWidth = sWidth;
00094   //initChamberSpecs();
00095   // Initialize output parameters  
00096   xWithinStrip = xWithinChamber;  
00097 
00098   CSCStripHit::ChannelContainer strips = stripHit.strips();
00099   int nStrips = strips.size();
00100   int centStrip = nStrips/2 + 1;   
00101   std::vector<float> adcs = stripHit.s_adc();
00102   int tmax = stripHit.tmax();
00103 
00104   // Fit peaking time only if using calibrations
00105   float t_peak = tpeak;
00106   float t_zero = 0.;
00107   float adc[4];
00108 
00109   if ( useCalib ) {
00110     bool useFittedCharge = false;
00111     for ( int t = 0; t < 4; ++t ) {
00112       int k  = t + 4 * (centStrip-1);
00113       adc[t] = adcs[k];
00114     }
00115     useFittedCharge = peakTimeFinder_->FindPeakTime( tmax, adc, t_zero, t_peak );  // Clumsy...  can remove t_peak ?
00116     tpeak = t_peak+t_zero;
00117   }
00118       
00119   //---- fill the charge matrix (3x3)
00120   int j = 0;
00121   for ( int i = 1; i <= nStrips; ++i ) {
00122     if ( i > (centStrip-2) && i < (centStrip+2) ) {
00123       std::vector<float> adcsFit;
00124       for ( int t = 0; t < 4; ++t ) {
00125         int k  = t + 4*(i-1);
00126         adc[t] = adcs[k];
00127         if ( t < 3) chargeSignal[j][t] = adc[t];
00128       }
00129       j++;
00130     }
00131   }
00132 
00133 
00134   // Load in x-talks:
00135 
00136   if ( useCalib ) {
00137     std::vector<float> xtalks;
00138     recoConditions_->crossTalk( id, centralStrip, xtalks );
00139     float dt = 50. * tmax - (t_peak + t_zero);  // QUESTION:  should it be only - (t_peak) ???
00140     // XTalks; l,r are for left, right XTalk; lr0,1,2 are for what charge "remains" in the strip 
00141     for ( int t = 0; t < 3; ++t ) {
00142       xt_l[0][t] = xtalks[0] * (50.* (t-1) + dt) + xtalks[1] + xtalksOffset;
00143       xt_r[0][t] = xtalks[2] * (50.* (t-1) + dt) + xtalks[3] + xtalksOffset;
00144       xt_l[1][t] = xtalks[4] * (50.* (t-1) + dt) + xtalks[5] + xtalksOffset;
00145       xt_r[1][t] = xtalks[6] * (50.* (t-1) + dt) + xtalks[7] + xtalksOffset;
00146       xt_l[2][t] = xtalks[8] * (50.* (t-1) + dt) + xtalks[9] + xtalksOffset;
00147       xt_r[2][t] = xtalks[10]* (50.* (t-1) + dt) + xtalks[11]+ xtalksOffset;
00148 
00149       xt_lr0[t] = (1. - xt_l[0][t] - xt_r[0][t]);
00150       xt_lr1[t] = (1. - xt_l[1][t] - xt_r[1][t]);
00151       xt_lr2[t] = (1. - xt_l[2][t] - xt_r[2][t]);
00152     }
00153   } else { 
00154     for ( int t = 0; t < 3; ++t ) {
00155       xt_l[0][t] = xtalksOffset;
00156       xt_r[0][t] = xtalksOffset;
00157       xt_l[1][t] = xtalksOffset; 
00158       xt_r[1][t] = xtalksOffset; 
00159       xt_l[2][t] = xtalksOffset; 
00160       xt_r[2][t] = xtalksOffset; 
00161 
00162       xt_lr0[t] = (1. - xt_l[0][t] - xt_r[0][t]);
00163       xt_lr1[t] = (1. - xt_l[1][t] - xt_r[1][t]);
00164       xt_lr2[t] = (1. - xt_l[2][t] - xt_r[2][t]);
00165     } 
00166   }
00167   
00168 
00169   // vector containing noise starts at tmax - 1, and tmax > 3, but....
00170   int tbin = tmax - 4;
00171 
00172   // .... originally, suppose to have tmax in tbin 4 or 5, but noticed in MTCC lots of 
00173   // hits with tmax == 3, so let's allow these, and shift down noise matrix by one element...
00174   // This is a patch because the calibration database doesn't have elements for tbin = 2, 
00175   // e.g. there is no element e[tmax-1,tmax+1] = e[2,4].
00176 
00177   if (tmax < 4) tbin = 0;    // patch
00178 
00179   // Load in auto-correlation noise matrices
00180   if ( useCalib ) {
00181     std::vector<float> nmatrix;
00182     recoConditions_->noiseMatrix( id, centralStrip, nmatrix );
00183     for ( int istrip =0; istrip < 3; ++istrip ) {
00184       a11[istrip] = nmatrix[0+tbin*3+istrip*15];
00185       a12[istrip] = nmatrix[1+tbin*3+istrip*15];
00186       a13[istrip] = nmatrix[2+tbin*3+istrip*15];
00187       a22[istrip] = nmatrix[3+tbin*3+istrip*15];
00188       a23[istrip] = nmatrix[4+tbin*3+istrip*15];
00189       a33[istrip] = nmatrix[6+tbin*3+istrip*15];
00190     }
00191   } else {
00192     // FIXME:  NO HARD WIRED VALUES !!!
00193     for ( int istrip =0; istrip < 3; ++istrip ) {
00194       a11[istrip] = 10.0;
00195       a12[istrip] = 0.0;
00196       a13[istrip] = 0.0;
00197       a22[istrip] = 10.0;
00198       a23[istrip] = 0.0;
00199       a33[istrip] = 10.0;
00200     }
00201   }
00202   
00203   //---- Set up noise, XTalk matrices 
00204   setupMatrix();
00205   //---- Calculate the coordinate within the strip and associate uncertainty
00206   bool ME1_1 = false;
00207   if("ME1/a" == specs_->chamberTypeName() || "ME1/b" == specs_->chamberTypeName()){
00208     ME1_1 = true;
00209   } 
00210 
00211   // due to cross talks and 3 time bin sum it is in principe possible that the center strip is not the maximum strip
00212   // in some cases the consequences could be quite extreme
00213   // take some measures against the extreme cases
00214   bool peakMismatch = false;
00215   std::vector <float> charges(3);
00216   charges[0] = q_sumL;
00217   charges[1] = q_sumC;
00218   charges[2] = q_sumR;
00219   int min_index = min_element(charges.begin(),charges.end()) - charges.begin();
00220   int max_index = max_element(charges.begin(),charges.end()) - charges.begin();
00221   if(1!=max_index && (1==min_index || charges[max_index]/q_sumC > 1.1)){// 10% below the maximum charge
00222     peakMismatch = true;
00223     switch (max_index){
00224     case 0:
00225       xWithinStrip = -1;
00226       break;
00227     case 2:
00228       xWithinStrip = 1;
00229       break;
00230     default:
00231       // should be an error message here
00232       xWithinStrip = 0;// in case?
00233       break;
00234     }
00235   }
00236   else{
00237     //
00238     xWithinStrip = float(calculateXonStripPosition(stripWidth, ME1_1));
00239   }
00240   xWithinChamber = xWithinChamber + (xWithinStrip * stripWidth);
00241   if(peakMismatch){
00242     sigma = stripWidth/sqrt(12);
00243   }
00244   else{
00245     
00246     //---- error estimation
00247     int factorStripWidth = int( sqrt(stripWidth/0.38) );
00248     int maxConsecutiveStrips = 8;
00249     if(factorStripWidth){
00250       maxConsecutiveStrips /=  factorStripWidth ;
00251     }
00252     maxConsecutiveStrips++;
00253     
00254     std::map <std::string, int> chamberTypes;
00255     chamberTypes["ME1/a"] = 1;
00256     chamberTypes["ME1/b"] = 2;
00257     chamberTypes["ME1/2"] = 3;
00258     chamberTypes["ME1/3"] = 4;
00259     chamberTypes["ME2/1"] = 5;
00260     chamberTypes["ME2/2"] = 6;
00261     chamberTypes["ME3/1"] = 7;
00262     chamberTypes["ME3/2"] = 8;
00263     chamberTypes["ME4/1"] = 9;
00264     chamberTypes["ME4/2"] = 8;
00265     
00266     switch(chamberTypes[specs_->chamberTypeName()]){
00267     case 1:
00268       noise_level  = noise_level_ME1a;
00269       xt_asymmetry = xt_asymmetry_ME1a;
00270       const_syst = const_syst_ME1a;
00271       break;
00272       
00273     case 2:
00274       noise_level  = noise_level_ME1b;
00275       xt_asymmetry = xt_asymmetry_ME1b;
00276       const_syst = const_syst_ME1b;
00277       
00278     case 3:
00279       noise_level  = noise_level_ME12;
00280       xt_asymmetry = xt_asymmetry_ME12;
00281       const_syst = const_syst_ME12;
00282       break;
00283       
00284     case 4:
00285       noise_level  = noise_level_ME13;
00286       xt_asymmetry = xt_asymmetry_ME13;
00287       const_syst = const_syst_ME13;
00288       break;
00289       
00290     case 5:
00291       noise_level  = noise_level_ME21;
00292       xt_asymmetry = xt_asymmetry_ME21;
00293       const_syst = const_syst_ME21;
00294       break;
00295       
00296     case 6:
00297       noise_level  = noise_level_ME22;
00298       xt_asymmetry = xt_asymmetry_ME22;
00299       const_syst = const_syst_ME22;
00300       break;
00301 
00302     case 7:
00303       noise_level  = noise_level_ME31;
00304       xt_asymmetry = xt_asymmetry_ME31;
00305       const_syst = const_syst_ME31;
00306       break;
00307 
00308     case 8:
00309       noise_level  = noise_level_ME32;
00310       xt_asymmetry = xt_asymmetry_ME32;
00311       const_syst = const_syst_ME32;
00312       break;
00313 
00314     case 9:
00315       noise_level  = noise_level_ME41;
00316       xt_asymmetry = xt_asymmetry_ME41;
00317       const_syst = const_syst_ME41;
00318       break;
00319 
00320     default:
00321       noise_level  = noise_level_ME22;
00322       xt_asymmetry = xt_asymmetry_ME22;
00323       const_syst = const_syst_ME22;
00324 
00325     }
00326     if(false==stripHit.isNearDeadStrip() &&
00327        stripHit.numberOfConsecutiveStrips()<maxConsecutiveStrips &&
00328        fabs(stripHit.closestMaximum())>maxConsecutiveStrips/2 ){
00329       sigma =  float(calculateXonStripError(stripWidth, ME1_1));
00330     }
00331     else{ //---- near dead strip or too close maxima or too wide strip cluster
00332       sigma = stripWidth/sqrt(12);
00333     }
00334   }
00335   quality_flag = 1;
00336 }
00337 
00338 /* setupMatrix
00339  *
00340  */
00341 void CSCXonStrip_MatchGatti::setupMatrix() {
00342   //---- a??? and v??[] could be skipped for now...; not used yet
00343 
00344   /*
00345   double dd, a11t, a12t, a13t, a22t, a23t, a33t;
00346   double syserr = adcSystematics;
00347   double xtlk_err = xtalksSystematics;
00348   // Left strip
00349   a11t = a11[0] + syserr*syserr * ChargeSignal[0][0]*ChargeSignal[0][0] + xtlk_err*xtlk_err*ChargeSignal[1][0]*ChargeSignal[1][0];
00350   a12t = a12[0] + syserr*syserr * ChargeSignal[0][0]*ChargeSignal[0][1];
00351   a13t = a13[0] + syserr*syserr * ChargeSignal[0][0]*ChargeSignal[0][2];
00352   a22t = a22[0] + syserr*syserr * ChargeSignal[0][1]*ChargeSignal[0][1] + xtlk_err*xtlk_err*ChargeSignal[1][1]*ChargeSignal[1][1];
00353   a23t = a23[0] + syserr*syserr * ChargeSignal[0][1]*ChargeSignal[0][2];
00354   a33t = a33[0] + syserr*syserr * ChargeSignal[0][2]*ChargeSignal[0][2] + xtlk_err*xtlk_err*ChargeSignal[1][2]*ChargeSignal[1][2];
00355 
00356   dd     = (a11t*a33t*a22t - a11t*a23t*a23t - a33t*a12t*a12t 
00357                        + 2.* a12t*a13t*a23t - a13t*a13t*a22t );
00358 
00359   v11[0] = (a33t*a22t - a23t*a23t)/dd;
00360   v12[0] =-(a33t*a12t - a13t*a23t)/dd;
00361   v13[0] = (a12t*a23t - a13t*a22t)/dd;
00362   v22[0] = (a33t*a11t - a13t*a13t)/dd;
00363   v23[0] =-(a23t*a11t - a12t*a13t)/dd;
00364   v33[0] = (a22t*a11t - a12t*a12t)/dd;
00365      
00366   // Center strip
00367   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]);
00368   a12t = a12[1] + syserr*syserr * ChargeSignal[1][0]*ChargeSignal[1][1];
00369   a13t = a13[1] + syserr*syserr * ChargeSignal[1][0]*ChargeSignal[1][2];
00370   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]);
00371   a23t = a23[1] + syserr*syserr * ChargeSignal[1][1]*ChargeSignal[1][2];
00372   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]);
00373 
00374   dd     = (a11t*a33t*a22t - a11t*a23t*a23t - a33t*a12t*a12t
00375                        + 2.* a12t*a13t*a23t - a13t*a13t*a22t );
00376 
00377   v11[1] = (a33t*a22t - a23t*a23t)/dd;
00378   v12[1] =-(a33t*a12t - a13t*a23t)/dd;
00379   v13[1] = (a12t*a23t - a13t*a22t)/dd;
00380   v22[1] = (a33t*a11t - a13t*a13t)/dd;
00381   v23[1] =-(a23t*a11t - a12t*a13t)/dd;
00382   v33[1] = (a22t*a11t - a12t*a12t)/dd;
00383 
00384   // Right strip
00385   a11t = a11[2] + syserr*syserr * ChargeSignal[2][0]*ChargeSignal[2][0] + xtlk_err*xtlk_err*ChargeSignal[1][0]*ChargeSignal[1][0];
00386   a12t = a12[2] + syserr*syserr * ChargeSignal[2][0]*ChargeSignal[2][1];
00387   a13t = a13[2] + syserr*syserr * ChargeSignal[2][0]*ChargeSignal[2][2];
00388   a22t = a22[2] + syserr*syserr * ChargeSignal[2][1]*ChargeSignal[2][1] + xtlk_err*xtlk_err*ChargeSignal[1][1]*ChargeSignal[1][1];
00389   a23t = a23[2] + syserr*syserr * ChargeSignal[2][1]*ChargeSignal[2][2];
00390   a33t = a33[2] + syserr*syserr * ChargeSignal[2][2]*ChargeSignal[2][2] + xtlk_err*xtlk_err*ChargeSignal[1][2]*ChargeSignal[1][2];
00391 
00392   dd     = (a11t*a33t*a22t - a11t*a23t*a23t - a33t*a12t*a12t
00393                         +2.* a12t*a13t*a23t - a13t*a13t*a22t );
00394 
00395   v11[2] = (a33t*a22t - a23t*a23t)/dd;
00396   v12[2] =-(a33t*a12t - a13t*a23t)/dd;
00397   v13[2] = (a12t*a23t - a13t*a22t)/dd;
00398   v22[2] = (a33t*a11t - a13t*a13t)/dd;
00399   v23[2] =-(a23t*a11t - a12t*a13t)/dd;
00400   v33[2] = (a22t*a11t - a12t*a12t)/dd;
00401 */
00402   //---- Find the inverted XTalk matrix and apply it to the charge (3x3)
00403   //---- Thus the charge before the XTalk is obtained
00404   HepMatrix cross_talks(3,3);
00405   HepMatrix cross_talks_inv(3,3);
00406   int err = 0;
00407   //---- q_sum is 3 time bins summed; L, C, R - left, central, right strips
00408   q_sum = q_sumL = q_sumC = q_sumR = 0.;
00409   double charge = 0.;
00410   for(int iTime=0;iTime<3;iTime++){
00411     cross_talks_inv(1,1) = cross_talks(1,1) = xt_lr0[iTime];
00412     cross_talks_inv(1,2) = cross_talks(1,2) = xt_l[1][iTime];
00413     cross_talks_inv(1,3) = cross_talks(1,3) = 0.;
00414     cross_talks_inv(2,1) = cross_talks(2,1) =  xt_r[0][iTime];
00415     cross_talks_inv(2,2) = cross_talks(2,2) = xt_lr1[iTime];
00416     cross_talks_inv(2,3) = cross_talks(2,3) = xt_l[2][iTime];
00417     cross_talks_inv(3,1) = cross_talks(3,1) = 0.;
00418     cross_talks_inv(3,2) = cross_talks(3,2) = xt_r[1][iTime];
00419     cross_talks_inv(3,3) = cross_talks(3,3) = xt_lr2[iTime];
00420     cross_talks_inv.invert(err);
00421     if (err != 0) {
00422       LogTrace("CSCRecHit")<<" Failed to invert XTalks matrix. Inaccurate cross-talk correction..."<<"\n";
00423     }
00424     //---- "charge" is XT-corrected charge
00425     charge = chargeSignal[0][iTime]*cross_talks_inv(1,1) + chargeSignal[1][iTime]*cross_talks_inv(1,2) + 
00426       chargeSignal[2][iTime]*cross_talks_inv(1,3);
00427     //---- Negative charge? According to studies (and logic) - better use 0 charge
00428     if(charge<0.){
00429       charge = 0.;
00430     }
00431     q_sum+=charge;
00432     q_sumL+=charge;
00433     charge = chargeSignal[0][iTime]*cross_talks_inv(2,1) + chargeSignal[1][iTime]*cross_talks_inv(2,2) + 
00434       chargeSignal[2][iTime]* cross_talks_inv(2,3);
00435     if(charge<0.){
00436       charge = 0.;
00437     }
00438     q_sum+=charge;
00439     q_sumC+=charge;
00440     charge = chargeSignal[0][iTime]*cross_talks_inv(3,1) + chargeSignal[1][iTime]*cross_talks_inv(3,2) + 
00441       chargeSignal[2][iTime]*cross_talks_inv(3,3);
00442     if(charge<0.){
00443       charge = 0.;
00444     }
00445     q_sum+=charge;
00446     q_sumR+=charge;
00447   }
00448 }
00449 
00450 
00451 /* initChamberSpecs
00452  *
00453  */
00454 void CSCXonStrip_MatchGatti::initChamberSpecs() {
00455   // Not used directly but these are parameters used for extracting the correction values
00456   // in coordinate and error estimators
00457 
00458   // Distance between anode and cathode
00459   h = specs_->anodeCathodeSpacing();
00460   r = h / stripWidth;
00461 
00462   // Wire spacing
00463   double wspace = specs_->wireSpacing();
00464 
00465   // Wire radius
00466   double wradius = specs_->wireRadius();
00467 
00468   // Accepted parameters in Gatti function
00469   const double parm[5] = {.1989337e-02, -.6901542e-04, .8665786, 154.6177, -.680163e-03 };
00470 
00471   k_3 = ( parm[0]*wspace/h + parm[1] )
00472       * ( parm[2]*wspace/wradius + parm[3] + parm[4]*(wspace/wradius)*(wspace/wradius) );
00473 
00474   sqrt_k_3 = sqrt( k_3 );
00475   norm     = r * (0.5 / atan( sqrt_k_3 )); // changed from norm to r * norm
00476   k_2      = M_PI_2 * ( 1. - sqrt_k_3 /2. );
00477   k_1      = 0.25 * k_2 * sqrt_k_3 / atan( sqrt_k_3 );
00478 }
00479 
00480 
00481 void CSCXonStrip_MatchGatti::getCorrectionValues(std::string estimator){
00482   hardcodedCorrectionInitialization();
00483 }
00484 
00485 double CSCXonStrip_MatchGatti::estimated2GattiCorrection(double x_estimated, float stripWidth, bool ME1_1) {
00486   //---- 11 "nominal" strip widths : 0.6 - 1.6 cm; for ME1_1 just 6 "nominal" strip widths : 0.3 - 0.8 cm; see HardCodedCorrectionInitialization() 
00487   //---- Calculate corrections at specific  Xestimated (linear interpolation between points)
00488   int n_SW;
00489   int min_SW;
00490   if(ME1_1){
00491     n_SW = n_SW_ME1_1;
00492     min_SW = 3; // min SW calculated is 0.3 cm
00493   }
00494   else{
00495     n_SW = n_SW_noME1_1;
00496     min_SW = 6;// min SW calculated is 0.6 cm
00497   }
00498   int stripDown = int(10.*stripWidth) - min_SW; // 0 is at min strip width calculated
00499   int stripUp = stripDown + 1;
00500   if(stripUp>n_SW-1){
00501     //---- to be checked...
00502     //if(stripUp>n_SW){
00503       //std::cout<<" Is strip width = "<<stripWidth<<" OK?" <<std::endl;
00504     //}
00505     stripUp = n_SW-1;
00506   }
00507 
00508   double half_strip_width = 0.5;
00509   //const int Nbins = 501;
00510   const int n_bins = n_val;
00511   double corr_2_xestim = 999.;
00512   if(stripDown<0){
00513     corr_2_xestim = 1;
00514   }
00515   else{
00516     //---- Parametrized x_gatti minus x_estimated differences
00517 
00518     int xestim_bin = -999;
00519     double delta_strip_width = 999.;
00520     double delta_strip_widthUpDown = 999.;
00521     double diff_2_strip_width = 999.;
00522     delta_strip_width = stripWidth - int(stripWidth*10)/10.;
00523     delta_strip_widthUpDown = 0.1;
00524 
00525     if(fabs(x_estimated)>0.5){
00526       if(fabs(x_estimated)>1.){
00527         corr_2_xestim = 1.;// for now; to be investigated
00528       }
00529       else{      
00530         //if(fabs(Xestimated)>0.55){
00531           //std::cout<<"X position from the estimated position above 0.55 (safty margin)?! "<<std::endl;
00532           //CorrToXc = 999.;
00533         //}
00534         xestim_bin = int((1.- fabs(x_estimated))/half_strip_width * n_bins);
00535         if(ME1_1){
00536           diff_2_strip_width = x_correction_ME1_1[stripUp][xestim_bin]-x_correction_ME1_1[stripDown][xestim_bin];
00537           corr_2_xestim =  x_correction_ME1_1[stripDown][xestim_bin] +
00538             (delta_strip_width/delta_strip_widthUpDown)*diff_2_strip_width ;
00539         }
00540         else{
00541           diff_2_strip_width = x_correction_noME1_1[stripUp][xestim_bin]-x_correction_noME1_1[stripDown][xestim_bin];
00542           corr_2_xestim =  x_correction_noME1_1[stripDown][xestim_bin] +
00543             (delta_strip_width/delta_strip_widthUpDown)*diff_2_strip_width ;
00544         }
00545         corr_2_xestim = -corr_2_xestim;
00546       }
00547     }
00548     else{
00549       xestim_bin = int((fabs(x_estimated)/half_strip_width) * n_bins);
00550       if(ME1_1){
00551         diff_2_strip_width = x_correction_ME1_1[stripUp][xestim_bin] - x_correction_ME1_1[stripDown][xestim_bin];
00552         corr_2_xestim =  x_correction_ME1_1[stripDown][xestim_bin] +
00553           (delta_strip_width/delta_strip_widthUpDown)*diff_2_strip_width ;
00554       }
00555       else{
00556         diff_2_strip_width = x_correction_noME1_1[stripUp][xestim_bin] - x_correction_noME1_1[stripDown][xestim_bin];
00557         corr_2_xestim =  x_correction_noME1_1[stripDown][xestim_bin] +
00558           (delta_strip_width/delta_strip_widthUpDown)*diff_2_strip_width ;
00559       }
00560     }
00561     if(x_estimated<0.){
00562       corr_2_xestim = -corr_2_xestim;
00563     }
00564   }
00565   
00566   return corr_2_xestim;
00567 }
00568 
00569 
00570 double CSCXonStrip_MatchGatti::estimated2Gatti(double x_estimated, float stripWidth, bool ME1_1) {
00571 
00572  int sign;
00573  if(x_estimated>0.){
00574    sign = 1;
00575  }
00576  else{
00577    sign = - 1;
00578  }
00579  double x_corr = estimated2GattiCorrection(x_estimated, stripWidth, ME1_1);
00580  double x_gatti = x_estimated + x_corr;
00581 
00582  return x_gatti;
00583 }
00584 
00585 double CSCXonStrip_MatchGatti::xfError_Noise(double noise){
00586 
00587   double min, max;
00588   if(q_sumR>q_sumL){
00589     min = q_sumL;
00590     max = q_sumR;
00591   }
00592   else{
00593     min = q_sumR;
00594     max = q_sumL;
00595   }
00596   //---- Error propagation...
00597   //---- Names here are fake! Due to technical features
00598   double dr_L2 = pow(q_sumR-q_sumL,2);
00599   double dr_C2 = pow(q_sumC-min,2);
00600   double dr_R2 = pow(q_sumC-max,2);
00601   double error = sqrt(dr_L2 + dr_C2 + dr_R2)*noise/pow(q_sumC-min,2)/2;
00602 
00603   return error;
00604 }
00605 
00606 double CSCXonStrip_MatchGatti::xfError_XTasym(double xt_asym){
00607 
00608   double min;
00609   if(q_sumR>q_sumL){
00610     min = q_sumL;
00611   }
00612   else{
00613     min = q_sumR;
00614   }
00615   //---- Error propagation
00616   double dXTL = (pow(q_sumC,2)+pow(q_sumR,2)-q_sumL*q_sumR-q_sumR*q_sumC);
00617   double dXTR = (pow(q_sumC,2)+pow(q_sumL,2)-q_sumL*q_sumR-q_sumL*q_sumC);
00618   double dXT = sqrt(pow(dXTL,2) + pow(dXTR,2))/pow((q_sumC-min),2)/2;
00619   double error = dXT * xt_asym;
00620 
00621   return error;
00622 }
00623 
00624 
00625 double CSCXonStrip_MatchGatti::calculateXonStripError(float stripWidth, bool ME1_1){
00626   double min;
00627   if(q_sumR>q_sumL){
00628     min = q_sumL;
00629   }
00630   else{
00631     min = q_sumR;
00632   }
00633   
00634   double xf = (q_sumR - q_sumL)/(q_sumC - min)/2;
00635   double xf_ErrorNoise = xfError_Noise(noise_level);
00636   double xf_ErrorXTasym = xfError_XTasym(xt_asymmetry);
00637   // x_G = x_F + correction_functon(x_F)
00638   // 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)]
00639   double d_xf = sqrt( pow( xf_ErrorNoise, 2) + pow( xf_ErrorXTasym, 2));
00640         double d_corr = estimated2GattiCorrection(xf+d_xf, stripWidth, ME1_1) - estimated2GattiCorrection(xf, stripWidth, ME1_1);
00641   double x_shift = d_xf + d_corr;
00642   //  double x_shift = sqrt( pow( xf_ErrorNoise, 2) + pow( xf_ErrorXTasym, 2)) * 
00643   //(1 + (estimated2GattiCorrection(xf+0.001, stripWidth, ME1_1) -
00644   //  estimated2GattiCorrection(xf, stripWidth, ME1_1))*1000.);
00645   double x_error =   sqrt( pow( fabs(x_shift)*stripWidth, 2) + pow(const_syst, 2) );
00646   return  x_error; 
00647 }
00648 
00649 double CSCXonStrip_MatchGatti::calculateXonStripPosition(float stripWidth, bool ME1_1){ 
00650 
00651   double  x_estimated = -99.;
00652   double min;
00653   if(q_sumR>q_sumL){
00654     min = q_sumL;
00655   }
00656   else{
00657     min = q_sumR;
00658   }
00659   //---- This is XF ( X Florida - after the first group that used it)  
00660   x_estimated = (q_sumR - q_sumL)/(q_sumC - min)/2;
00661   double x_gatti = estimated2Gatti(x_estimated, stripWidth, ME1_1);
00662   return x_gatti;
00663 }
00664 
00665