Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:25:59 

0001 // This is CSCXonStrip_MatchGatti.cc
0002 
0003 //---- Large part is copied from RecHitB
0004 //---- author: Stoyan Stoynev - NU
0005 
0006 #include <RecoLocalMuon/CSCRecHitD/src/CSCXonStrip_MatchGatti.h>
0007 #include <RecoLocalMuon/CSCRecHitD/src/CSCStripHit.h>
0008 
0009 #include <Geometry/CSCGeometry/interface/CSCLayer.h>
0010 #include <Geometry/CSCGeometry/interface/CSCChamberSpecs.h>
0011 
0012 #include <CondFormats/CSCObjects/interface/CSCDBGains.h>
0013 #include <CondFormats/DataRecord/interface/CSCDBGainsRcd.h>
0014 #include <CondFormats/CSCObjects/interface/CSCDBCrosstalk.h>
0015 #include <CondFormats/DataRecord/interface/CSCDBCrosstalkRcd.h>
0016 #include <CondFormats/CSCObjects/interface/CSCDBNoiseMatrix.h>
0017 #include <CondFormats/DataRecord/interface/CSCDBNoiseMatrixRcd.h>
0018 
0019 #include <FWCore/MessageLogger/interface/MessageLogger.h>
0020 #include <FWCore/Utilities/interface/Exception.h>
0021 
0022 #include "DataFormats/CLHEP/interface/AlgebraicObjects.h"
0023 
0024 #include <vector>
0025 #include <cmath>
0026 #include <iostream>
0027 #include <fstream>
0028 //#include <iomanip.h>
0029 //#include <iomanip>
0030 
0031 #ifndef M_PI_2
0032 #define M_PI_2 1.57079632679489661923
0033 #endif
0034 
0035 CSCXonStrip_MatchGatti::CSCXonStrip_MatchGatti(const edm::ParameterSet& ps) : recoConditions_(nullptr) {
0036   useCalib = ps.getParameter<bool>("CSCUseCalibrations");
0037   xtalksOffset = ps.getParameter<double>("CSCStripxtalksOffset");
0038   noise_level_ME1a = ps.getParameter<double>("NoiseLevel_ME1a");
0039   xt_asymmetry_ME1a = ps.getParameter<double>("XTasymmetry_ME1a");
0040   const_syst_ME1a = ps.getParameter<double>("ConstSyst_ME1a");
0041   noise_level_ME1b = ps.getParameter<double>("NoiseLevel_ME1b");
0042   xt_asymmetry_ME1b = ps.getParameter<double>("XTasymmetry_ME1b");
0043   const_syst_ME1b = ps.getParameter<double>("ConstSyst_ME1b");
0044   noise_level_ME12 = ps.getParameter<double>("NoiseLevel_ME12");
0045   xt_asymmetry_ME12 = ps.getParameter<double>("XTasymmetry_ME12");
0046   const_syst_ME12 = ps.getParameter<double>("ConstSyst_ME12");
0047   noise_level_ME13 = ps.getParameter<double>("NoiseLevel_ME13");
0048   xt_asymmetry_ME13 = ps.getParameter<double>("XTasymmetry_ME13");
0049   const_syst_ME13 = ps.getParameter<double>("ConstSyst_ME13");
0050   noise_level_ME21 = ps.getParameter<double>("NoiseLevel_ME21");
0051   xt_asymmetry_ME21 = ps.getParameter<double>("XTasymmetry_ME21");
0052   const_syst_ME21 = ps.getParameter<double>("ConstSyst_ME21");
0053   noise_level_ME22 = ps.getParameter<double>("NoiseLevel_ME22");
0054   xt_asymmetry_ME22 = ps.getParameter<double>("XTasymmetry_ME22");
0055   const_syst_ME22 = ps.getParameter<double>("ConstSyst_ME22");
0056   noise_level_ME31 = ps.getParameter<double>("NoiseLevel_ME31");
0057   xt_asymmetry_ME31 = ps.getParameter<double>("XTasymmetry_ME31");
0058   const_syst_ME31 = ps.getParameter<double>("ConstSyst_ME31");
0059   noise_level_ME32 = ps.getParameter<double>("NoiseLevel_ME32");
0060   xt_asymmetry_ME32 = ps.getParameter<double>("XTasymmetry_ME32");
0061   const_syst_ME32 = ps.getParameter<double>("ConstSyst_ME32");
0062   noise_level_ME41 = ps.getParameter<double>("NoiseLevel_ME41");
0063   xt_asymmetry_ME41 = ps.getParameter<double>("XTasymmetry_ME41");
0064   const_syst_ME41 = ps.getParameter<double>("ConstSyst_ME41");
0065 
0066   getCorrectionValues("StringCurrentlyNotUsed");
0067 }
0068 
0069 CSCXonStrip_MatchGatti::~CSCXonStrip_MatchGatti() {}
0070 
0071 /* findPosition
0072  *
0073  */
0074 void CSCXonStrip_MatchGatti::findXOnStrip(const CSCDetId& id,
0075                                           const CSCLayer* layer,
0076                                           const CSCStripHit& stripHit,
0077                                           int centralStrip,
0078                                           float& xWithinChamber,
0079                                           float& sWidth,
0080                                           const float& tpeak,
0081                                           float& xWithinStrip,
0082                                           float& sigma,
0083                                           int& quality_flag) {
0084   quality_flag = 0;
0085   // Initialize Gatti parameters using chamberSpecs
0086   // Cache specs_ info for ease of access
0087   specs_ = layer->chamber()->specs();
0088   stripWidth = sWidth;
0089   //initChamberSpecs();
0090   // Initialize output parameters
0091   xWithinStrip = xWithinChamber;
0092 
0093   CSCStripHit::ChannelContainer const& strips = stripHit.strips();
0094   int nStrips = strips.size();
0095   int centStrip = nStrips / 2 + 1;
0096   std::vector<float> const& adcs = stripHit.s_adc();
0097   int tmax = stripHit.tmax();
0098 
0099   //// Fit peaking time only if using calibrations
0100   //float t_peak = tpeak;
0101   //float t_zero = 0.;
0102   //float adc[4];
0103   //
0104   //if ( useCalib ) {
0105   //
0106   //  for ( int t = 0; t < 4; ++t ) {
0107   //    int k  = t + 4 * (centStrip-1);
0108   //    adc[t] = adcs[k];
0109   //  }
0110   //
0111   //  // t_peak from peak finder is now 'absolute' i.e. in ns from start of sca time bin 0
0112   //  t_peak = peakTimeFinder_->peakTime( tmax, adc, t_peak );
0113   //  // Just for completeness, the start time of the pulse is 133 ns earlier, according to Stan :)
0114   //  t_zero = t_peak - 133.;
0115   //  // and reset tpeak since that's the way it gets passed out of this function (Argh!)
0116   //  tpeak = t_peak;
0117   //  LogTrace("CSCRecHit|CSCXonStrip_MatchGatti") << "CSCXonStrip_MatchGatti: " <<
0118   //     id << " strip=" << centralStrip << ", t_zero=" << t_zero << ", t_peak=" << t_peak;
0119   //}
0120 
0121   //---- fill the charge matrix (3x3)
0122   float adc[4];
0123   int j = 0;
0124   for (int i = 1; i <= nStrips; ++i) {
0125     if (i > (centStrip - 2) && i < (centStrip + 2)) {
0126       std::vector<float> adcsFit;
0127       for (int t = 0; t < 4; ++t) {
0128         int k = t + 4 * (i - 1);
0129         adc[t] = adcs[k];
0130         if (t < 3)
0131           chargeSignal[j][t] = adc[t];
0132       }
0133       j++;
0134     }
0135   }
0136 
0137   // Load in x-talks:
0138 
0139   if (useCalib) {
0140     std::vector<float> xtalks;
0141     recoConditions_->crossTalk(id, centralStrip, xtalks);
0142     float dt = 50.f * tmax - tpeak;
0143     // XTalks; l,r are for left, right XTalk; lr0,1,2 are for what charge "remains" in the strip
0144     for (int t = 0; t < 3; ++t) {
0145       xt_l[0][t] = xtalks[0] * (50.f * (t - 1) + dt) + xtalks[1] + xtalksOffset;
0146       xt_r[0][t] = xtalks[2] * (50.f * (t - 1) + dt) + xtalks[3] + xtalksOffset;
0147       xt_l[1][t] = xtalks[4] * (50.f * (t - 1) + dt) + xtalks[5] + xtalksOffset;
0148       xt_r[1][t] = xtalks[6] * (50.f * (t - 1) + dt) + xtalks[7] + xtalksOffset;
0149       xt_l[2][t] = xtalks[8] * (50.f * (t - 1) + dt) + xtalks[9] + xtalksOffset;
0150       xt_r[2][t] = xtalks[10] * (50.f * (t - 1) + dt) + xtalks[11] + xtalksOffset;
0151 
0152       xt_lr0[t] = (1. - xt_l[0][t] - xt_r[0][t]);
0153       xt_lr1[t] = (1. - xt_l[1][t] - xt_r[1][t]);
0154       xt_lr2[t] = (1. - xt_l[2][t] - xt_r[2][t]);
0155     }
0156   } else {
0157     for (int t = 0; t < 3; ++t) {
0158       xt_l[0][t] = xtalksOffset;
0159       xt_r[0][t] = xtalksOffset;
0160       xt_l[1][t] = xtalksOffset;
0161       xt_r[1][t] = xtalksOffset;
0162       xt_l[2][t] = xtalksOffset;
0163       xt_r[2][t] = xtalksOffset;
0164 
0165       xt_lr0[t] = (1. - xt_l[0][t] - xt_r[0][t]);
0166       xt_lr1[t] = (1. - xt_l[1][t] - xt_r[1][t]);
0167       xt_lr2[t] = (1. - xt_l[2][t] - xt_r[2][t]);
0168     }
0169   }
0170 
0171   // vector containing noise starts at tmax - 1, and tmax > 3, but....
0172   int tbin = tmax - 4;
0173 
0174   // .... originally, suppose to have tmax in tbin 4 or 5, but noticed in MTCC lots of
0175   // hits with tmax == 3, so let's allow these, and shift down noise matrix by one element...
0176   // This is a patch because the calibration database doesn't have elements for tbin = 2,
0177   // e.g. there is no element e[tmax-1,tmax+1] = e[2,4].
0178 
0179   if (tmax < 4)
0180     tbin = 0;  // patch
0181 
0182   // Load in auto-correlation noise matrices
0183   if (useCalib) {
0184     std::vector<float> nmatrix;
0185     recoConditions_->noiseMatrix(id, centralStrip, nmatrix);
0186     for (int istrip = 0; istrip < 3; ++istrip) {
0187       a11[istrip] = nmatrix[0 + tbin * 3 + istrip * 15];
0188       a12[istrip] = nmatrix[1 + tbin * 3 + istrip * 15];
0189       a13[istrip] = nmatrix[2 + tbin * 3 + istrip * 15];
0190       a22[istrip] = nmatrix[3 + tbin * 3 + istrip * 15];
0191       a23[istrip] = nmatrix[4 + tbin * 3 + istrip * 15];
0192       a33[istrip] = nmatrix[6 + tbin * 3 + istrip * 15];
0193     }
0194   } else {
0195     // FIXME:  NO HARD WIRED VALUES !!!
0196     for (int istrip = 0; istrip < 3; ++istrip) {
0197       a11[istrip] = 10.0;
0198       a12[istrip] = 0.0;
0199       a13[istrip] = 0.0;
0200       a22[istrip] = 10.0;
0201       a23[istrip] = 0.0;
0202       a33[istrip] = 10.0;
0203     }
0204   }
0205 
0206   //---- Set up noise, XTalk matrices
0207   setupMatrix();
0208   //---- Calculate the coordinate within the strip and associate uncertainty
0209 
0210   static const std::string ME1a("ME1/a");
0211   static const std::string ME1b("ME1/b");
0212 
0213   bool ME1_1 = (ME1a == specs_->chamberTypeName() || ME1b == specs_->chamberTypeName());
0214 
0215   // due to cross talks and 3 time bin sum it is in principe possible that the center strip is not the maximum strip
0216   // in some cases the consequences could be quite extreme
0217   // take some measures against the extreme cases
0218   bool peakMismatch = false;
0219   std::vector<float> charges(3);
0220   charges[0] = q_sumL;
0221   charges[1] = q_sumC;
0222   charges[2] = q_sumR;
0223   int min_index = min_element(charges.begin(), charges.end()) - charges.begin();
0224   int max_index = max_element(charges.begin(), charges.end()) - charges.begin();
0225   if (1 != max_index && (1 == min_index ||
0226                          // the condition below means that if the initial position estimate within strip (|xF|)
0227                          // is above  1.1/2 = 0.55 "strip widths" peakMismatch is set to true (so special case);
0228                          // in normal cases |xF|<=0.5 (0.5 is at the edge of the "central" strip)
0229                          (charges[max_index] - charges[min_index]) > 1.1 * (q_sumC - charges[min_index]))) {
0230     peakMismatch = true;
0231     switch (max_index) {
0232       case 0:
0233         xWithinStrip = -1;
0234         break;
0235       case 2:
0236         xWithinStrip = 1;
0237         break;
0238       default:
0239         // should be an error message here
0240         xWithinStrip = 0;  // in case?
0241         break;
0242     }
0243   }
0244   // we don't have the needed information (it's similar to the "edge" strip case)
0245   //else if(stripHit.isNearDeadStrip()){
0246   else if (stripHit.deadStrip() > 0) {
0247     xWithinStrip = 0;
0248   } else {
0249     //
0250     xWithinStrip = float(calculateXonStripPosition(stripWidth, ME1_1));
0251   }
0252   xWithinChamber = xWithinChamber + (xWithinStrip * stripWidth);
0253   if (peakMismatch) {
0254     sigma = stripWidth / std::sqrt(12.f);
0255   } else {
0256     //---- error estimation
0257     int factorStripWidth = int(std::sqrt(stripWidth / 0.38f));
0258     int maxConsecutiveStrips = 8;
0259     if (factorStripWidth) {
0260       maxConsecutiveStrips /= factorStripWidth;
0261     }
0262     maxConsecutiveStrips++;
0263 
0264     struct ChamberTypes {
0265       std::map<std::string, int> chamberTypes;
0266       int operator()(std::string const& s) const {
0267         auto p = chamberTypes.find(s);
0268         return p != chamberTypes.end() ? (*p).second : 0;
0269       }
0270       ChamberTypes() {
0271         chamberTypes["ME1/a"] = 1;
0272         chamberTypes["ME1/b"] = 2;
0273         chamberTypes["ME1/2"] = 3;
0274         chamberTypes["ME1/3"] = 4;
0275         chamberTypes["ME2/1"] = 5;
0276         chamberTypes["ME2/2"] = 6;
0277         chamberTypes["ME3/1"] = 7;
0278         chamberTypes["ME3/2"] = 8;
0279         chamberTypes["ME4/1"] = 9;
0280         chamberTypes["ME4/2"] = 8;
0281       }
0282     };
0283     static const ChamberTypes chamberTypes;
0284 
0285     switch (chamberTypes(specs_->chamberTypeName())) {
0286       case 1:
0287         noise_level = noise_level_ME1a;
0288         xt_asymmetry = xt_asymmetry_ME1a;
0289         const_syst = const_syst_ME1a;
0290         break;
0291 
0292       case 2:
0293         noise_level = noise_level_ME1b;
0294         xt_asymmetry = xt_asymmetry_ME1b;
0295         const_syst = const_syst_ME1b;
0296         break;
0297 
0298       case 3:
0299         noise_level = noise_level_ME12;
0300         xt_asymmetry = xt_asymmetry_ME12;
0301         const_syst = const_syst_ME12;
0302         break;
0303 
0304       case 4:
0305         noise_level = noise_level_ME13;
0306         xt_asymmetry = xt_asymmetry_ME13;
0307         const_syst = const_syst_ME13;
0308         break;
0309 
0310       case 5:
0311         noise_level = noise_level_ME21;
0312         xt_asymmetry = xt_asymmetry_ME21;
0313         const_syst = const_syst_ME21;
0314         break;
0315 
0316       case 6:
0317         noise_level = noise_level_ME22;
0318         xt_asymmetry = xt_asymmetry_ME22;
0319         const_syst = const_syst_ME22;
0320         break;
0321 
0322       case 7:
0323         noise_level = noise_level_ME31;
0324         xt_asymmetry = xt_asymmetry_ME31;
0325         const_syst = const_syst_ME31;
0326         break;
0327 
0328       case 8:
0329         noise_level = noise_level_ME32;
0330         xt_asymmetry = xt_asymmetry_ME32;
0331         const_syst = const_syst_ME32;
0332         break;
0333 
0334       case 9:
0335         noise_level = noise_level_ME41;
0336         xt_asymmetry = xt_asymmetry_ME41;
0337         const_syst = const_syst_ME41;
0338         break;
0339 
0340       default:
0341         noise_level = noise_level_ME22;
0342         xt_asymmetry = xt_asymmetry_ME22;
0343         const_syst = const_syst_ME22;
0344     }
0345     if (0 == stripHit.deadStrip() && stripHit.numberOfConsecutiveStrips() < maxConsecutiveStrips &&
0346         std::abs(stripHit.closestMaximum()) > maxConsecutiveStrips / 2) {
0347       sigma = float(calculateXonStripError(stripWidth, ME1_1));
0348     } else {  //---- near dead strip or too close maxima or too wide strip cluster
0349       sigma = stripWidth / std::sqrt(12.f);
0350     }
0351   }
0352   quality_flag = 1;
0353 }
0354 
0355 /* setupMatrix
0356  *
0357  */
0358 void CSCXonStrip_MatchGatti::setupMatrix() {
0359   //---- a??? and v??[] could be skipped for now...; not used yet
0360 
0361   /*
0362   double dd, a11t, a12t, a13t, a22t, a23t, a33t;
0363   double syserr = adcSystematics;
0364   double xtlk_err = xtalksSystematics;
0365   // Left strip
0366   a11t = a11[0] + syserr*syserr * ChargeSignal[0][0]*ChargeSignal[0][0] + xtlk_err*xtlk_err*ChargeSignal[1][0]*ChargeSignal[1][0];
0367   a12t = a12[0] + syserr*syserr * ChargeSignal[0][0]*ChargeSignal[0][1];
0368   a13t = a13[0] + syserr*syserr * ChargeSignal[0][0]*ChargeSignal[0][2];
0369   a22t = a22[0] + syserr*syserr * ChargeSignal[0][1]*ChargeSignal[0][1] + xtlk_err*xtlk_err*ChargeSignal[1][1]*ChargeSignal[1][1];
0370   a23t = a23[0] + syserr*syserr * ChargeSignal[0][1]*ChargeSignal[0][2];
0371   a33t = a33[0] + syserr*syserr * ChargeSignal[0][2]*ChargeSignal[0][2] + xtlk_err*xtlk_err*ChargeSignal[1][2]*ChargeSignal[1][2];
0372 
0373   dd     = (a11t*a33t*a22t - a11t*a23t*a23t - a33t*a12t*a12t 
0374                        + 2.* a12t*a13t*a23t - a13t*a13t*a22t );
0375 
0376   v11[0] = (a33t*a22t - a23t*a23t)/dd;
0377   v12[0] =-(a33t*a12t - a13t*a23t)/dd;
0378   v13[0] = (a12t*a23t - a13t*a22t)/dd;
0379   v22[0] = (a33t*a11t - a13t*a13t)/dd;
0380   v23[0] =-(a23t*a11t - a12t*a13t)/dd;
0381   v33[0] = (a22t*a11t - a12t*a12t)/dd;
0382      
0383   // Center strip
0384   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]);
0385   a12t = a12[1] + syserr*syserr * ChargeSignal[1][0]*ChargeSignal[1][1];
0386   a13t = a13[1] + syserr*syserr * ChargeSignal[1][0]*ChargeSignal[1][2];
0387   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]);
0388   a23t = a23[1] + syserr*syserr * ChargeSignal[1][1]*ChargeSignal[1][2];
0389   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]);
0390 
0391   dd     = (a11t*a33t*a22t - a11t*a23t*a23t - a33t*a12t*a12t
0392                        + 2.* a12t*a13t*a23t - a13t*a13t*a22t );
0393 
0394   v11[1] = (a33t*a22t - a23t*a23t)/dd;
0395   v12[1] =-(a33t*a12t - a13t*a23t)/dd;
0396   v13[1] = (a12t*a23t - a13t*a22t)/dd;
0397   v22[1] = (a33t*a11t - a13t*a13t)/dd;
0398   v23[1] =-(a23t*a11t - a12t*a13t)/dd;
0399   v33[1] = (a22t*a11t - a12t*a12t)/dd;
0400 
0401   // Right strip
0402   a11t = a11[2] + syserr*syserr * ChargeSignal[2][0]*ChargeSignal[2][0] + xtlk_err*xtlk_err*ChargeSignal[1][0]*ChargeSignal[1][0];
0403   a12t = a12[2] + syserr*syserr * ChargeSignal[2][0]*ChargeSignal[2][1];
0404   a13t = a13[2] + syserr*syserr * ChargeSignal[2][0]*ChargeSignal[2][2];
0405   a22t = a22[2] + syserr*syserr * ChargeSignal[2][1]*ChargeSignal[2][1] + xtlk_err*xtlk_err*ChargeSignal[1][1]*ChargeSignal[1][1];
0406   a23t = a23[2] + syserr*syserr * ChargeSignal[2][1]*ChargeSignal[2][2];
0407   a33t = a33[2] + syserr*syserr * ChargeSignal[2][2]*ChargeSignal[2][2] + xtlk_err*xtlk_err*ChargeSignal[1][2]*ChargeSignal[1][2];
0408 
0409   dd     = (a11t*a33t*a22t - a11t*a23t*a23t - a33t*a12t*a12t
0410                         +2.* a12t*a13t*a23t - a13t*a13t*a22t );
0411 
0412   v11[2] = (a33t*a22t - a23t*a23t)/dd;
0413   v12[2] =-(a33t*a12t - a13t*a23t)/dd;
0414   v13[2] = (a12t*a23t - a13t*a22t)/dd;
0415   v22[2] = (a33t*a11t - a13t*a13t)/dd;
0416   v23[2] =-(a23t*a11t - a12t*a13t)/dd;
0417   v33[2] = (a22t*a11t - a12t*a12t)/dd;
0418 */
0419   //---- Find the inverted XTalk matrix and apply it to the charge (3x3)
0420   //---- Thus the charge before the XTalk is obtained
0421   CLHEP::HepMatrix cross_talks_inv(3, 3);
0422   int err = 0;
0423   //---- q_sum is 3 time bins summed; L, C, R - left, central, right strips
0424   q_sum = q_sumL = q_sumC = q_sumR = 0.;
0425   double charge = 0.;
0426   for (int iTime = 0; iTime < 3; iTime++) {
0427     cross_talks_inv(1, 1) = xt_lr0[iTime];
0428     cross_talks_inv(1, 2) = xt_l[1][iTime];
0429     cross_talks_inv(1, 3) = 0.;
0430     cross_talks_inv(2, 1) = xt_r[0][iTime];
0431     cross_talks_inv(2, 2) = xt_lr1[iTime];
0432     cross_talks_inv(2, 3) = xt_l[2][iTime];
0433     cross_talks_inv(3, 1) = 0.;
0434     cross_talks_inv(3, 2) = xt_r[1][iTime];
0435     cross_talks_inv(3, 3) = xt_lr2[iTime];
0436     cross_talks_inv.invert(err);
0437     if (err != 0) {
0438       edm::LogWarning("FailedXTalkiInversionNoCrosstalkCorrection")
0439           << "Failed to invert XTalks matrix. No cross-talk correction for this rechit.";
0440       //edm::LogError("CSCRecHit") << "Failed to invert XTalks matrix. No cross-talk correction for this rechit.";
0441       return;
0442     }
0443     //---- "charge" is XT-corrected charge
0444     charge = chargeSignal[0][iTime] * cross_talks_inv(1, 1) + chargeSignal[1][iTime] * cross_talks_inv(1, 2) +
0445              chargeSignal[2][iTime] * cross_talks_inv(1, 3);
0446     //---- Negative charge? According to studies (and logic) - better use 0 charge
0447     //----- Later studies suggest that this only do harm. I am still worried about
0448     // charges of -50 ADC and below (0.5% of the cases) but let see
0449     //if(charge<0.){
0450     //charge = 0.;
0451     //}
0452     q_sum += charge;
0453     q_sumL += charge;
0454     charge = chargeSignal[0][iTime] * cross_talks_inv(2, 1) + chargeSignal[1][iTime] * cross_talks_inv(2, 2) +
0455              chargeSignal[2][iTime] * cross_talks_inv(2, 3);
0456     //if(charge<0.){
0457     //charge = 0.;
0458     //}
0459     q_sum += charge;
0460     q_sumC += charge;
0461     charge = chargeSignal[0][iTime] * cross_talks_inv(3, 1) + chargeSignal[1][iTime] * cross_talks_inv(3, 2) +
0462              chargeSignal[2][iTime] * cross_talks_inv(3, 3);
0463     //if(charge<0.){
0464     //charge = 0.;
0465     //}
0466     q_sum += charge;
0467     q_sumR += charge;
0468   }
0469 }
0470 
0471 /* initChamberSpecs
0472  *
0473  */
0474 void CSCXonStrip_MatchGatti::initChamberSpecs() {
0475   // Not used directly but these are parameters used for extracting the correction values
0476   // in coordinate and error estimators
0477 
0478   // Distance between anode and cathode
0479   h = specs_->anodeCathodeSpacing();
0480   r = h / stripWidth;
0481 
0482   // Wire spacing
0483   double wspace = specs_->wireSpacing();
0484 
0485   // Wire radius
0486   double wradius = specs_->wireRadius();
0487 
0488   // Accepted parameters in Gatti function
0489   const double parm[5] = {.1989337e-02, -.6901542e-04, .8665786, 154.6177, -.680163e-03};
0490 
0491   k_3 = (parm[0] * wspace / h + parm[1]) *
0492         (parm[2] * wspace / wradius + parm[3] + parm[4] * (wspace / wradius) * (wspace / wradius));
0493 
0494   sqrt_k_3 = std::sqrt(k_3);
0495   norm = r * (0.5 / std::atan(sqrt_k_3));  // changed from norm to r * norm
0496   k_2 = M_PI_2 * (1. - sqrt_k_3 / 2.);
0497   k_1 = 0.25 * k_2 * sqrt_k_3 / std::atan(sqrt_k_3);
0498 }
0499 
0500 void CSCXonStrip_MatchGatti::getCorrectionValues(std::string estimator) { hardcodedCorrectionInitialization(); }
0501 
0502 double CSCXonStrip_MatchGatti::estimated2GattiCorrection(double x_estimated, float stripWidth, bool ME1_1) {
0503   //---- 11 "nominal" strip widths : 0.6 - 1.6 cm; for ME1_1 just 6 "nominal" strip widths : 0.3 - 0.8 cm; see HardCodedCorrectionInitialization()
0504   //---- Calculate corrections at specific  Xestimated (linear interpolation between points)
0505   int n_SW;
0506   int min_SW;
0507   if (ME1_1) {
0508     n_SW = n_SW_ME1_1;
0509     min_SW = 3;  // min SW calculated is 0.3 cm
0510   } else {
0511     n_SW = n_SW_noME1_1;
0512     min_SW = 6;  // min SW calculated is 0.6 cm
0513   }
0514   int stripDown = int(10. * stripWidth) - min_SW;  // 0 is at min strip width calculated
0515   int stripUp = stripDown + 1;
0516   if (stripUp > n_SW - 1) {
0517     //---- to be checked...
0518     //if(stripUp>n_SW){
0519     //std::cout<<" Is strip width = "<<stripWidth<<" OK?" <<std::endl;
0520     //}
0521     stripUp = n_SW - 1;
0522   }
0523 
0524   double half_strip_width = 0.5;
0525   //const int Nbins = 501;
0526   const int n_bins = n_val;
0527   double corr_2_xestim = 999.;
0528   if (stripDown < 0) {
0529     corr_2_xestim = 1;
0530   } else {
0531     //---- Parametrized x_gatti minus x_estimated differences
0532 
0533     int xestim_bin = -999;
0534     double delta_strip_width = 999.;
0535     double delta_strip_widthUpDown = 999.;
0536     double diff_2_strip_width = 999.;
0537     delta_strip_width = stripWidth - int(stripWidth * 10) / 10.;
0538     delta_strip_widthUpDown = 0.1;
0539 
0540     if (std::abs(x_estimated) > 0.5) {
0541       if (std::abs(x_estimated) > 1.) {
0542         corr_2_xestim = 1.;  // for now; to be investigated
0543       } else {
0544         //if(std::abs(Xestimated)>0.55){
0545         //std::cout<<"X position from the estimated position above 0.55 (safty margin)?! "<<std::endl;
0546         //CorrToXc = 999.;
0547         //}
0548         xestim_bin = int((1. - std::abs(x_estimated)) / half_strip_width * n_bins);
0549         if (ME1_1) {
0550           diff_2_strip_width = x_correction_ME1_1[stripUp][xestim_bin] - x_correction_ME1_1[stripDown][xestim_bin];
0551           corr_2_xestim = x_correction_ME1_1[stripDown][xestim_bin] +
0552                           (delta_strip_width / delta_strip_widthUpDown) * diff_2_strip_width;
0553         } else {
0554           diff_2_strip_width = x_correction_noME1_1[stripUp][xestim_bin] - x_correction_noME1_1[stripDown][xestim_bin];
0555           corr_2_xestim = x_correction_noME1_1[stripDown][xestim_bin] +
0556                           (delta_strip_width / delta_strip_widthUpDown) * diff_2_strip_width;
0557         }
0558         corr_2_xestim = -corr_2_xestim;
0559       }
0560     } else {
0561       xestim_bin = int((std::abs(x_estimated) / half_strip_width) * n_bins);
0562       if (ME1_1) {
0563         diff_2_strip_width = x_correction_ME1_1[stripUp][xestim_bin] - x_correction_ME1_1[stripDown][xestim_bin];
0564         corr_2_xestim = x_correction_ME1_1[stripDown][xestim_bin] +
0565                         (delta_strip_width / delta_strip_widthUpDown) * diff_2_strip_width;
0566       } else {
0567         diff_2_strip_width = x_correction_noME1_1[stripUp][xestim_bin] - x_correction_noME1_1[stripDown][xestim_bin];
0568         corr_2_xestim = x_correction_noME1_1[stripDown][xestim_bin] +
0569                         (delta_strip_width / delta_strip_widthUpDown) * diff_2_strip_width;
0570       }
0571     }
0572     if (x_estimated < 0.) {
0573       corr_2_xestim = -corr_2_xestim;
0574     }
0575   }
0576 
0577   return corr_2_xestim;
0578 }
0579 
0580 double CSCXonStrip_MatchGatti::estimated2Gatti(double x_estimated, float stripWidth, bool ME1_1) {
0581   double x_corr = estimated2GattiCorrection(x_estimated, stripWidth, ME1_1);
0582   double x_gatti = x_estimated + x_corr;
0583 
0584   return x_gatti;
0585 }
0586 
0587 double CSCXonStrip_MatchGatti::xfError_Noise(double noise) {
0588   double min, max;
0589   if (q_sumR > q_sumL) {
0590     min = q_sumL;
0591     max = q_sumR;
0592   } else {
0593     min = q_sumR;
0594     max = q_sumL;
0595   }
0596   //---- Error propagation...
0597   //---- Names here are fake! Due to technical features
0598   double dr_L2 = pow(q_sumR - q_sumL, 2);
0599   double dr_C2 = pow(q_sumC - min, 2);
0600   double dr_R2 = pow(q_sumC - max, 2);
0601   double error = std::sqrt(dr_L2 + dr_C2 + dr_R2) * noise / std::pow(q_sumC - min, 2) / 2;
0602 
0603   return error;
0604 }
0605 
0606 double CSCXonStrip_MatchGatti::xfError_XTasym(double xt_asym) {
0607   double min;
0608   if (q_sumR > q_sumL) {
0609     min = q_sumL;
0610   } else {
0611     min = q_sumR;
0612   }
0613   //---- Error propagation
0614   double dXTL = (std::pow(q_sumC, 2) + std::pow(q_sumR, 2) - q_sumL * q_sumR - q_sumR * q_sumC);
0615   double dXTR = (std::pow(q_sumC, 2) + std::pow(q_sumL, 2) - q_sumL * q_sumR - q_sumL * q_sumC);
0616   double dXT = std::sqrt(std::pow(dXTL, 2) + std::pow(dXTR, 2)) / std::pow((q_sumC - min), 2) / 2;
0617   double error = dXT * xt_asym;
0618 
0619   return error;
0620 }
0621 
0622 double CSCXonStrip_MatchGatti::calculateXonStripError(float stripWidth, bool ME1_1) {
0623   double min;
0624   if (q_sumR > q_sumL) {
0625     min = q_sumL;
0626   } else {
0627     min = q_sumR;
0628   }
0629 
0630   double xf = (q_sumR - q_sumL) / (q_sumC - min) / 2;
0631   double xf_ErrorNoise = xfError_Noise(noise_level);
0632   double xf_ErrorXTasym = xfError_XTasym(xt_asymmetry);
0633   // x_G = x_F + correction_functon(x_F)
0634   // 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)]
0635   double d_xf = std::sqrt(std::pow(xf_ErrorNoise, 2) + std::pow(xf_ErrorXTasym, 2));
0636   double d_corr =
0637       estimated2GattiCorrection(xf + d_xf, stripWidth, ME1_1) - estimated2GattiCorrection(xf, stripWidth, ME1_1);
0638   double x_shift = d_xf + d_corr;
0639   //  double x_shift = sqrt( pow( xf_ErrorNoise, 2) + pow( xf_ErrorXTasym, 2)) *
0640   //(1 + (estimated2GattiCorrection(xf+0.001, stripWidth, ME1_1) -
0641   //  estimated2GattiCorrection(xf, stripWidth, ME1_1))*1000.);
0642   double x_error = std::sqrt(std::pow(std::abs(x_shift) * stripWidth, 2) + std::pow(const_syst, 2));
0643   return x_error;
0644 }
0645 
0646 double CSCXonStrip_MatchGatti::calculateXonStripPosition(float stripWidth, bool ME1_1) {
0647   double x_estimated = -99.;
0648   double min;
0649   if (q_sumR > q_sumL) {
0650     min = q_sumL;
0651   } else {
0652     min = q_sumR;
0653   }
0654   //---- This is XF ( X Florida - after the first group that used it)
0655   x_estimated = (q_sumR - q_sumL) / (q_sumC - min) / 2;
0656   double x_gatti = estimated2Gatti(x_estimated, stripWidth, ME1_1);
0657   return x_gatti;
0658 }
0659 
0660 // Define space for statics
0661 const int CSCXonStrip_MatchGatti::n_val;
0662 const int CSCXonStrip_MatchGatti::n_SW_noME1_1;
0663 const int CSCXonStrip_MatchGatti::n_SW_ME1_1;
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!