CSCRecHitD/src/CSCMake2DRecHit.cc

0001 // This is CSCMake2DRecHit
0002
0003 #include <RecoLocalMuon/CSCRecHitD/src/CSCMake2DRecHit.h>
0004 #include <RecoLocalMuon/CSCRecHitD/src/CSCXonStrip_MatchGatti.h>
0005 #include <RecoLocalMuon/CSCRecHitD/src/CSCStripHit.h>
0006 #include <RecoLocalMuon/CSCRecHitD/src/CSCWireHit.h>
0007 #include <RecoLocalMuon/CSCRecHitD/src/CSCRecoConditions.h>
0008
0009 #include <DataFormats/CSCRecHit/interface/CSCRecHit2D.h>
0010 #include <DataFormats/MuonDetId/interface/CSCDetId.h>
0011
0012 #include <Geometry/CSCGeometry/interface/CSCLayer.h>
0013 #include <Geometry/CSCGeometry/interface/CSCChamberSpecs.h>
0014 #include <Geometry/CSCGeometry/interface/CSCLayerGeometry.h>
0015
0016 #include <FWCore/MessageLogger/interface/MessageLogger.h>
0017 #include <FWCore/Utilities/interface/Exception.h>
0018
0019 #include <iostream>
0020 #include <algorithm>
0021 #include <cmath>
0022 #include <string>
0023
0024 CSCMake2DRecHit::CSCMake2DRecHit(const edm::ParameterSet& ps) : peakTimeFinder_(new CSCFindPeakTime(ps)) {
0025   useCalib = ps.getParameter<bool>("CSCUseCalibrations");
0026   stripWireDeltaTime = ps.getParameter<int>("CSCstripWireDeltaTime");  //@@ Non-standard  CSC*s*trip...
0027   useTimingCorrections = ps.getParameter<bool>("CSCUseTimingCorrections");
0028   useGasGainCorrections = ps.getParameter<bool>("CSCUseGasGainCorrections");
0029
0030   xMatchGatti_ = new CSCXonStrip_MatchGatti(ps);
0031 }
0032
0033 CSCMake2DRecHit::~CSCMake2DRecHit() { delete xMatchGatti_; }
0034
0035 CSCRecHit2D CSCMake2DRecHit::hitFromStripAndWire(const CSCDetId& id,
0036                                                  const CSCLayer* layer,
0037                                                  const CSCWireHit& wHit,
0038                                                  const CSCStripHit& sHit) {
0039   // Cache layer info for ease of access
0040   layer_ = layer;
0041   layergeom_ = layer_->geometry();
0042   specs_ = layer->chamber()->specs();
0043   id_ = id;
0044
0045   static constexpr float inv_sqrt_12 = 0.2886751;
0046
0047   float tpeak = -99.f;
0048
0049   CSCRecHit2D::ADCContainer adcMap;
0050   CSCRecHit2D::ChannelContainer wgroups;
0051
0052   // Find wire hit position and wire properties
0053   wgroups = wHit.wgroups();
0054
0055   int wg_left = wgroups[0];
0056   ;
0057   int wg_right = wgroups[wgroups.size() - 1];
0058
0059   int Nwires1 = layergeom_->numberOfWiresPerGroup(wg_left);
0060   int Nwires2 = layergeom_->numberOfWiresPerGroup(wg_right);
0061
0062   float Mwire1 = layergeom_->middleWireOfGroup(wg_left);
0063   float Mwire2 = layergeom_->middleWireOfGroup(wg_right);
0064
0065   int centerWire_left = (int)(Mwire1 - Nwires1 / 2. + 0.5);
0066   int centerWire_right = (int)(Mwire2 + Nwires2 / 2.);
0067
0068   float centerWire = (centerWire_left + centerWire_right) / 2.;
0069
0070   //---- WGs around dead HV segment regions may need special treatment...
0071   //---- This is not addressed here.
0072
0073   float sigmaWire = 0.;
0074   if (wHit.deadWG() > 0 || wgroups.size() > 2) {
0075     //---- worst possible case; take most conservative approach
0076     for (unsigned int iWG = 0; iWG < wgroups.size(); iWG++) {
0077       sigmaWire += layergeom_->yResolution(wgroups[iWG]);
0078     }
0079   } else if (2 == wgroups.size()) {
0080     //---- 2 WGs - get the larger error (overestimation if a single track is passing
0081     //---- between the WGs; underestimation if there are two separate signal sources)
0082     if (layergeom_->yResolution(wgroups[0]) > layergeom_->yResolution(wgroups[1])) {
0083       sigmaWire = layergeom_->yResolution(wgroups[0]);
0084     } else {
0085       sigmaWire = layergeom_->yResolution(wgroups[1]);
0086     }
0087   } else if (1 == wgroups.size()) {
0088     //---- simple - just 1 WG
0089     sigmaWire = layergeom_->yResolution(wgroups[0]);
0090   }
0091
0092   // Find strips position and properties
0093
0094   CSCRecHit2D::ChannelContainer const& strips = sHit.strips();
0095   int tmax = sHit.tmax();
0096   int nStrip = strips.size();
0097   int idCenterStrip = nStrip / 2;
0098   int centerStrip = strips[idCenterStrip];
0099
0100   // Retrieve strip pulseheights from the CSCStripHit
0101   const std::vector<float>& adc = sHit.s_adc();
0102   const std::vector<float>& adcRaw = sHit.s_adcRaw();
0103
0104   std::vector<float> adc2;
0105   std::vector<float> adc2Raw;
0106
0107   LogTrace("CSCMake2DRecHit") << "[CSCMake2DRecHit] dump of adc values to be added to rechit follows...";
0108
0109   for (int iStrip = 0; iStrip < nStrip; ++iStrip) {
0110     adc2.clear();
0111     adc2Raw.clear();
0112     adc2.reserve(4);
0113     adc2Raw.reserve(4);
0114     for (int t = 0; t < 4; ++t) {
0115       adc2.push_back(adc[t + iStrip * 4]);
0116       adc2Raw.push_back(adcRaw[t + iStrip * 4]);
0117     }
0118     //After CMSSW_5_0: ADC value is pedestal-subtracted and electronics-gain corrected
0119     adcMap.put(strips[iStrip], adc2.begin(), adc2.end());
0120     // Up to CMSSW_5_0, Rechit takes _raw_ adc values
0121     // adcMap.put( strips[iStrip], adc2Raw.begin(), adc2Raw.end() );
0122
0123     LogTrace("CSCMake2DRecHit") << "[CSCMake2DRecHit] strip = " << strips[iStrip] << " adcs= " << adc2Raw[0] << " "
0124                                 << adc2Raw[1] << " " << adc2Raw[2] << " " << adc2Raw[3];
0125   }
0126
0127   //The tpeak finding for both edge and non-edge strips has been moved to here
0128   //tpeak will be a const argument for xMatchGatti_->findXOnStrip
0129   float adcArray[4];
0130   for (int t = 0; t < 4; ++t) {
0131     int k = t + 4 * (idCenterStrip);
0132     adcArray[t] = adc[k];
0133   }
0134   tpeak = peakTimeFinder_->peakTime(tmax, adcArray, tpeak);
0135   // Just for completeness, the start time of the pulse is 133 ns earlier, according to Stan :)
0136   LogTrace("CSCRecHit") << "[CSCMake2DRecHit] " << id << " strip=" << centerStrip << ", t_zero=" << tpeak - 133.f
0137                         << ", tpeak=" << tpeak;
0138
0139   float positionWithinTheStrip = -99.;
0140   float sigmaWithinTheStrip = -99.;
0141   int quality = -1;
0142   LocalPoint lp0(0., 0.);
0143   int wglo = wg_left;
0144   int wghi = wg_right;
0145
0146   // First wire of first wiregroup in cluster
0147   wglo = layergeom_->middleWireOfGroup(wglo) + 0.5 - 0.5 * layergeom_->numberOfWiresPerGroup(wglo);
0148   // Last wire in last wiregroup of cluster (OK if only 1 wiregroup too)
0149   wghi = layergeom_->middleWireOfGroup(wghi) - 0.5 + 0.5 * layergeom_->numberOfWiresPerGroup(wghi);
0150
0151   float stripWidth = -99.f;
0152   // If at the edge, then used 1 strip cluster only
0153   if (centerStrip == 1 || centerStrip == specs_->nStrips() || nStrip < 2) {
0154     lp0 = (layergeom_->possibleRecHitPosition(float(centerStrip) - 0.5, wglo, wghi)).first;
0155     positionWithinTheStrip = 0.f;
0156     stripWidth = layergeom_->stripPitch(lp0);
0157     sigmaWithinTheStrip = stripWidth * inv_sqrt_12;
0158     quality = 2;
0159   } else {
0160     // If not at the edge, used cluster of size ClusterSize:
0161     LocalPoint lp11 = (layergeom_->possibleRecHitPosition(float(centerStrip) - 0.5, wglo, wghi)).first;
0162     float ymiddle = lp11.y();
0163     stripWidth = layergeom_->stripPitch(lp11);
0164
0165     //---- Calculate local position within the strip
0166     float xWithinChamber = lp11.x();
0167     quality = 0;
0168     //check if strip wire intersect is within the sensitive area of chamber
0169     if (layergeom_->inside(lp11)) {
0170       xMatchGatti_->findXOnStrip(id,
0171                                  layer_,
0172                                  sHit,
0173                                  centerStrip,
0174                                  xWithinChamber,
0175                                  stripWidth,
0176                                  tpeak,
0177                                  positionWithinTheStrip,
0178                                  sigmaWithinTheStrip,
0179                                  quality);
0180     }
0181     lp0 = LocalPoint(xWithinChamber, ymiddle);
0182   }
0183
0184   // compute the errors in local x and y
0185   LocalError localerr = layergeom_->localError(centerStrip, sigmaWithinTheStrip, sigmaWire);
0186
0187   // Before storing the recHit, take the opportunity to change its time
0188   if (useTimingCorrections) {
0189     float chipCorrection = recoConditions_->chipCorrection(id, centerStrip);
0190     float phaseCorrection = (sHit.stripsl1a()[0] >> (15 - 0) & 0x1) * 25.;
0191     float chamberCorrection = recoConditions_->chamberTimingCorrection(id);
0192
0193     GlobalPoint gp0 = layer_->toGlobal(lp0);
0194     float tofCorrection = gp0.mag() / 29.9792458;
0195     float signalPropagationSpeed[11] = {0.0, -78, -76, -188, -262, -97, -99, -90, -99, -99, -113};
0196     float position = lp0.y() / sin(layergeom_->stripAngle(centerStrip));
0197     float yCorrection = position / signalPropagationSpeed[id_.iChamberType()];
0198     //printf("RecHit in e:%d s:%d r:%d c:%d l:%d strip:%d \n",id.endcap(),id.station(), id.ring(),id.chamber(),id.layer(),centerStrip);
0199     //printf("\t tpeak before = %5.2f \t chipCorr %5.2f phaseCorr %5.2f chamberCorr %5.2f tofCorr %5.2f \n",
0200     //     tpeak,chipCorrection, phaseCorrection,chamberCorrection,tofCorrection);
0201     //printf("localy = %5.2f, yCorr = %5.2f \n",lp0.y(),yCorrection);
0202     tpeak = tpeak + chipCorrection + phaseCorrection + chamberCorrection - tofCorrection + yCorrection;
0203     //printf("\t tpeak after = %5.2f\n",tpeak);
0204   }
0205
0206   // Calculate wire time to the half bx level using time bins on
0207   // Store wire time with a precision of 0.01 as an int (multiply by 100)
0208   // Convert from bx to ns (multiply by 25)
0209   int scaledWireTime = 100 * findWireBx(wHit.timeBinsOn(), tpeak, id) * 25;
0210
0211   //Get the gas-gain correction for this rechit
0212   float gasGainCorrection = -999.f;
0213   if (useGasGainCorrections) {
0214     gasGainCorrection = recoConditions_->gasGainCorrection(id, centerStrip, centerWire);
0215   }
0216
0217   /// Correct the 3x3 ADC sum into the energy deposited in the layer.  Note:  this correction will
0218   /// give you dE.  In order to get the dE/dX, you will need to divide by the path length...
0219   /// If the algorithm to compute the corrected energy fails, flag it by a specific nonsense value:
0220   /// If the user has chosen not to use the gas gain correction --->  -998.
0221   /// If the gas gain correction from the database is a bad value ->  -997.
0222   /// If it is an edge strip -------------------------------------->  -996.
0223   /// If gas-gain is OK, but the ADC vector is the wrong size  ---->  -999.
0224   /// If the user has created the Rechit without the energy deposit>  -995.
0225   float energyDeposit = -999.f;
0226   if (gasGainCorrection < -998.f) {
0227     // if the user has chosen not to use the gas gain correction, set the energy to a different nonsense value
0228     energyDeposit = -998.f;
0229
0230   } else if (gasGainCorrection < 0.f) {
0231     // if the gas gain correction from the database is a bad value, set the energy to yet another nonsense value
0232     energyDeposit = -997.f;
0233
0234   } else {
0235     // gas-gain correction is OK, correct the 3x3 ADC sum
0236     if (adcMap.size() == 12) {
0237       energyDeposit = adcMap[0] * gasGainCorrection + adcMap[1] * gasGainCorrection + adcMap[2] * gasGainCorrection +
0238                       adcMap[4] * gasGainCorrection + adcMap[5] * gasGainCorrection + adcMap[6] * gasGainCorrection +
0239                       adcMap[8] * gasGainCorrection + adcMap[9] * gasGainCorrection + adcMap[10] * gasGainCorrection;
0240
0241     } else if (adcMap.size() == 4) {
0242       // if this is an edge strip, set the energy to yet another nonsense value
0243       energyDeposit = -996.f;
0244     }
0245   }
0246
0247   /// store rechit
0248
0249   /// Retrive the L1APhase+strips combination
0250   CSCRecHit2D::ChannelContainer const& L1A_and_strips = sHit.stripsTotal();  /// L1A
0251   /// Retrive the Bx + wgroups combination
0252   CSCRecHit2D::ChannelContainer const& BX_and_wgroups = wHit.wgroupsBXandWire();  /// BX
0253   // (sigmaWithinTheStrip/stripWidth) is in strip widths just like positionWithinTheStrip is!
0254
0255   CSCRecHit2D rechit(id,
0256                      lp0,
0257                      localerr,
0258                      L1A_and_strips,  /// L1A;
0259                      //adcMap, wgroups, tpeak, positionWithinTheStrip,
0260                      adcMap,
0261                      BX_and_wgroups,
0262                      tpeak,
0263                      positionWithinTheStrip,  /// BX
0264                      sigmaWithinTheStrip / stripWidth,
0265                      quality,
0266                      sHit.deadStrip(),
0267                      wHit.deadWG(),
0268                      scaledWireTime,
0269                      energyDeposit);
0270
0271   /// To see RecHit content (L1A feature included) (to be commented out)
0272   // rechit.print();
0273
0274   LogTrace("CSCRecHit") << "[CSCMake2DRecHit] rechit created in layer " << id << "... \n" << rechit;
0275
0276   return rechit;
0277 }
0278
0279 bool CSCMake2DRecHit::isHitInFiducial(const CSCLayer* layer, const CSCRecHit2D& rh) {
0280   bool isInFiducial = true;
0281   const CSCLayerGeometry* layergeom_ = layer->geometry();
0282   LocalPoint rhPosition = rh.localPosition();
0283   // is the rechit within the chamber?
0284   //(the problem occurs in ME11a/b otherwise it is OK)
0285   // we could use also
0286   //bool inside( const Local3DPoint&, const LocalError&, float scale=1.f ) const;
0287   if (!layergeom_->inside(rhPosition)) {
0288     isInFiducial = false;
0289   }
0290
0291   return isInFiducial;
0292 }
0293
0294 void CSCMake2DRecHit::setConditions(const CSCRecoConditions* reco) {
0295   xMatchGatti_->setConditions(reco);
0296   // And cache for use here
0297   recoConditions_ = reco;
0298 }
0299
0300 float CSCMake2DRecHit::findWireBx(const std::vector<int>& timeBinsOn, float tpeak, const CSCDetId& id) {
0301   // Determine the wire Bx from the vector of time bins on for the wire digi with peak time as an intial estimate.
0302   // Assumes that a single hit should create either one time bin on or two consecutive time bins on
0303   // so algorithm looks for bin on nearest to peak time and checks if it has a bin on consecutive with it
0304   float anode_bx_offset = recoConditions_->anodeBXoffset(id);
0305   float wireBx = -1;
0306   float timeGuess = tpeak / 25.f + anode_bx_offset;
0307   float diffMin = 9999.f;
0308   int bestMatch = -9;
0309   for (int j = 0; j < (int)timeBinsOn.size(); j++) {
0310     auto diff = timeGuess - timeBinsOn[j];
0311     // Find bin on closest to peak time
0312     if (std::abs(diff) < std::abs(diffMin)) {
0313       diffMin = diff;
0314       bestMatch = j;
0315       wireBx = timeBinsOn[j];
0316     }
0317   }
0318   int side = diffMin > 0 ? 1 : -1;  // diffMin/fabs(diffMin);
0319   bool unchanged = true;
0320   // First check if bin on the same side as peak time is on
0321   if ((bestMatch + side) > -1 &&
0322       (bestMatch + side) < (int)timeBinsOn.size()) {  // Make sure one next to it within vector limits
0323     if (timeBinsOn[bestMatch] ==
0324         (timeBinsOn[bestMatch + side] - side)) {  // See if next bin on in vector is consecutive in time
0325       // Set time to the average of the two bins
0326       wireBx = wireBx + 0.5f * side;
0327       unchanged = false;
0328     }
0329   }
0330   // If no match is found then check the other side
0331   if ((bestMatch - side) > -1 && (bestMatch - side) < (int)timeBinsOn.size() &&
0332       unchanged) {                                                         // Make sure one next to it exists
0333     if (timeBinsOn[bestMatch] == (timeBinsOn[bestMatch - side] + side)) {  // See if nextbin on is consecutive in time
0334       wireBx = wireBx - 0.5f * side;
0335     }
0336   }
0337   return wireBx - anode_bx_offset;  // expect collision muons to be centered near 0 bx
0338 }