Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 11:56:44 

0001 #include <iomanip>
0002 #include "TrackingTools/TrackFitters/interface/TrajectoryStateCombiner.h"
0003 #include "FWCore/ServiceRegistry/interface/Service.h"
0004 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0005 #include "TrackingTools/PatternTools/interface/Trajectory.h"
0006 #include "Alignment/MuonAlignmentAlgorithms/plugins/CSCOverlapsAlignmentAlgorithm.h"
0007 #include "Alignment/MuonAlignmentAlgorithms/plugins/CSCPairResidualsConstraint.h"
0008 
0009 double CSCPairResidualsConstraint::value() const {
0010   double delta = (m_sum1 * m_sumxx) - (m_sumx * m_sumx);
0011   assert(delta > 0.);
0012   if (m_parent->m_mode == kModePhiy || m_parent->m_mode == kModePhiPos || m_parent->m_mode == kModeRadius) {
0013     return ((m_sumxx * m_sumy) - (m_sumx * m_sumxy)) / delta;
0014   } else if (m_parent->m_mode == kModePhiz) {
0015     return ((m_sum1 * m_sumxy) - (m_sumx * m_sumy)) / delta;
0016   } else
0017     assert(false);
0018 }
0019 
0020 double CSCPairResidualsConstraint::error() const {
0021   if (m_parent->m_errorFromRMS) {
0022     assert(m_sum1 > 0.);
0023     return sqrt((m_sumyy / m_sum1) - pow(m_sumy / m_sum1, 2)) / sqrt(m_sumN);
0024   } else {
0025     double delta = (m_sum1 * m_sumxx) - (m_sumx * m_sumx);
0026     assert(delta > 0.);
0027     if (m_parent->m_mode == kModePhiy || m_parent->m_mode == kModePhiPos || m_parent->m_mode == kModeRadius) {
0028       return sqrt(m_sumxx / delta);
0029     } else if (m_parent->m_mode == kModePhiz) {
0030       return sqrt(m_sum1 / delta);
0031     } else
0032       assert(false);
0033   }
0034 }
0035 
0036 bool CSCPairResidualsConstraint::valid() const { return (m_sumN >= m_parent->m_minTracksPerOverlap); }
0037 
0038 void CSCPairResidualsConstraint::configure(CSCOverlapsAlignmentAlgorithm *parent) {
0039   m_parent = parent;
0040 
0041   if (m_parent->m_makeHistograms) {
0042     edm::Service<TFileService> tFileService;
0043 
0044     std::stringstream name, name2, name3, title;
0045     title << "i =" << m_id_i << " j =" << m_id_j;
0046 
0047     name << "slopeResiduals_" << m_identifier;
0048     m_slopeResiduals = tFileService->make<TH1F>(name.str().c_str(), title.str().c_str(), 300, -30., 30.);
0049 
0050     name2 << "offsetResiduals_" << m_identifier;
0051     m_offsetResiduals = tFileService->make<TH1F>(name2.str().c_str(), title.str().c_str(), 300, -30., 30.);
0052 
0053     name3 << "radial_" << m_identifier;
0054     m_radial = tFileService->make<TH1F>(name3.str().c_str(), title.str().c_str(), 700, 0., 700.);
0055   } else {
0056     m_slopeResiduals = nullptr;
0057     m_offsetResiduals = nullptr;
0058     m_radial = nullptr;
0059   }
0060 }
0061 
0062 void CSCPairResidualsConstraint::setZplane(const CSCGeometry *cscGeometry) {
0063   m_cscGeometry = cscGeometry;
0064 
0065   m_Zplane = (m_cscGeometry->idToDet(m_id_i)->surface().position().z() +
0066               m_cscGeometry->idToDet(m_id_j)->surface().position().z()) /
0067              2.;
0068   m_averageRadius = (m_cscGeometry->idToDet(m_id_i)->surface().position().perp() +
0069                      m_cscGeometry->idToDet(m_id_j)->surface().position().perp()) /
0070                     2.;
0071 
0072   m_iZ = m_cscGeometry->idToDet(m_id_i)->surface().position().z();
0073   m_jZ = m_cscGeometry->idToDet(m_id_j)->surface().position().z();
0074 
0075   CSCDetId i1(m_id_i.endcap(), m_id_i.station(), m_id_i.ring(), m_id_i.chamber(), 1);
0076   CSCDetId i6(m_id_i.endcap(), m_id_i.station(), m_id_i.ring(), m_id_i.chamber(), 6);
0077   CSCDetId j1(m_id_j.endcap(), m_id_j.station(), m_id_j.ring(), m_id_j.chamber(), 1);
0078   CSCDetId j6(m_id_j.endcap(), m_id_j.station(), m_id_j.ring(), m_id_j.chamber(), 6);
0079   m_iZ1 = m_cscGeometry->idToDet(m_id_i)->surface().toLocal(m_cscGeometry->idToDet(i1)->surface().position()).z();
0080   m_iZ6 = m_cscGeometry->idToDet(m_id_i)->surface().toLocal(m_cscGeometry->idToDet(i6)->surface().position()).z();
0081   m_jZ1 = m_cscGeometry->idToDet(m_id_j)->surface().toLocal(m_cscGeometry->idToDet(j1)->surface().position()).z();
0082   m_jZ6 = m_cscGeometry->idToDet(m_id_j)->surface().toLocal(m_cscGeometry->idToDet(j6)->surface().position()).z();
0083 
0084   m_Zsurface = Plane::build(Plane::PositionType(0., 0., m_Zplane), Plane::RotationType());
0085 }
0086 
0087 void CSCPairResidualsConstraint::setPropagator(const Propagator *propagator) { m_propagator = propagator; }
0088 
0089 bool CSCPairResidualsConstraint::addTrack(const std::vector<TrajectoryMeasurement> &measurements,
0090                                           const reco::TransientTrack &track,
0091                                           const TrackTransformer *trackTransformer) {
0092   std::vector<const TransientTrackingRecHit *> hits_i;
0093   std::vector<const TransientTrackingRecHit *> hits_j;
0094 
0095   for (std::vector<TrajectoryMeasurement>::const_iterator measurement = measurements.begin();
0096        measurement != measurements.end();
0097        ++measurement) {
0098     const TransientTrackingRecHit *hit = &*(measurement->recHit());
0099 
0100     DetId id = hit->geographicalId();
0101     if (id.det() == DetId::Muon && id.subdetId() == MuonSubdetId::CSC) {
0102       CSCDetId cscid(id.rawId());
0103       CSCDetId chamberId(cscid.endcap(), cscid.station(), cscid.ring(), cscid.chamber(), 0);
0104       if (m_parent->m_combineME11 && cscid.station() == 1 && cscid.ring() == 4)
0105         chamberId = CSCDetId(cscid.endcap(), 1, 1, cscid.chamber(), 0);
0106 
0107       if (chamberId == m_id_i)
0108         hits_i.push_back(hit);
0109       if (chamberId == m_id_j)
0110         hits_j.push_back(hit);
0111     }
0112   }
0113 
0114   if (m_parent->m_makeHistograms) {
0115     m_parent->m_hitsPerChamber->Fill(hits_i.size());
0116     m_parent->m_hitsPerChamber->Fill(hits_j.size());
0117   }
0118 
0119   // require minimum number of hits (if the requirement is too low (~2), some NANs might result...)
0120   if (int(hits_i.size()) < m_parent->m_minHitsPerChamber || int(hits_j.size()) < m_parent->m_minHitsPerChamber)
0121     return false;
0122 
0123   // maybe require segments to be fiducial
0124   if (m_parent->m_fiducial && !(isFiducial(hits_i, true) && isFiducial(hits_j, false)))
0125     return false;
0126 
0127   double intercept_i = 0.;
0128   double interceptError2_i = 0.;
0129   double slope_i = 0.;
0130   double slopeError2_i = 0.;
0131   double intercept_j = 0.;
0132   double interceptError2_j = 0.;
0133   double slope_j = 0.;
0134   double slopeError2_j = 0.;
0135 
0136   // if slopeFromTrackRefit, then you'll need to refit the whole track without this station's hits;
0137   // need at least two other stations for that to be reliable
0138   if (m_parent->m_slopeFromTrackRefit) {
0139     double dphidz;
0140     if (dphidzFromTrack(measurements, track, trackTransformer, dphidz)) {
0141       double sum1_i = 0.;
0142       double sumy_i = 0.;
0143       for (std::vector<const TransientTrackingRecHit *>::const_iterator hit = hits_i.begin(); hit != hits_i.end();
0144            ++hit) {
0145         double phi, phierr2;
0146         calculatePhi(*hit, phi, phierr2, false, true);
0147         double z = (*hit)->globalPosition().z() - m_Zplane;
0148 
0149         double weight = 1.;
0150         if (m_parent->m_useHitWeights)
0151           weight = 1. / phierr2;
0152         sum1_i += weight;
0153         sumy_i += weight * (phi - z * dphidz);
0154       }
0155 
0156       double sum1_j = 0.;
0157       double sumy_j = 0.;
0158       for (std::vector<const TransientTrackingRecHit *>::const_iterator hit = hits_j.begin(); hit != hits_j.end();
0159            ++hit) {
0160         double phi, phierr2;
0161         calculatePhi(*hit, phi, phierr2, false, true);
0162         double z = (*hit)->globalPosition().z() - m_Zplane;
0163 
0164         double weight = 1.;
0165         if (m_parent->m_useHitWeights)
0166           weight = 1. / phierr2;
0167         sum1_j += weight;
0168         sumy_j += weight * (phi - z * dphidz);
0169       }
0170 
0171       if (sum1_i != 0. && sum1_j != 0.) {
0172         slope_i = slope_j = dphidz;
0173 
0174         intercept_i = sumy_i / sum1_i;
0175         interceptError2_i = 1. / sum1_i;
0176 
0177         intercept_j = sumy_j / sum1_j;
0178         interceptError2_j = 1. / sum1_j;
0179       } else
0180         return false;
0181     }
0182   }
0183 
0184   else {  // not slopeFromTrackRefit
0185     double sum1_i = 0.;
0186     double sumx_i = 0.;
0187     double sumy_i = 0.;
0188     double sumxx_i = 0.;
0189     double sumxy_i = 0.;
0190     for (std::vector<const TransientTrackingRecHit *>::const_iterator hit = hits_i.begin(); hit != hits_i.end();
0191          ++hit) {
0192       double phi, phierr2;
0193       calculatePhi(*hit, phi, phierr2, false, true);
0194       double z = (*hit)->globalPosition().z() - m_Zplane;
0195 
0196       double weight = 1.;
0197       if (m_parent->m_useHitWeights)
0198         weight = 1. / phierr2;
0199       sum1_i += weight;
0200       sumx_i += weight * z;
0201       sumy_i += weight * phi;
0202       sumxx_i += weight * z * z;
0203       sumxy_i += weight * z * phi;
0204     }
0205 
0206     double sum1_j = 0.;
0207     double sumx_j = 0.;
0208     double sumy_j = 0.;
0209     double sumxx_j = 0.;
0210     double sumxy_j = 0.;
0211     for (std::vector<const TransientTrackingRecHit *>::const_iterator hit = hits_j.begin(); hit != hits_j.end();
0212          ++hit) {
0213       double phi, phierr2;
0214       calculatePhi(*hit, phi, phierr2, false, true);
0215       double z = (*hit)->globalPosition().z() - m_Zplane;
0216 
0217       double weight = 1.;
0218       if (m_parent->m_useHitWeights)
0219         weight = 1. / phierr2;
0220       sum1_j += weight;
0221       sumx_j += weight * z;
0222       sumy_j += weight * phi;
0223       sumxx_j += weight * z * z;
0224       sumxy_j += weight * z * phi;
0225     }
0226 
0227     double delta_i = (sum1_i * sumxx_i) - (sumx_i * sumx_i);
0228     double delta_j = (sum1_j * sumxx_j) - (sumx_j * sumx_j);
0229     if (delta_i != 0. && delta_j != 0.) {
0230       intercept_i = ((sumxx_i * sumy_i) - (sumx_i * sumxy_i)) / delta_i;
0231       interceptError2_i = sumxx_i / delta_i;
0232       slope_i = ((sum1_i * sumxy_i) - (sumx_i * sumy_i)) / delta_i;
0233       slopeError2_i = sum1_i / delta_i;
0234 
0235       intercept_j = ((sumxx_j * sumy_j) - (sumx_j * sumxy_j)) / delta_j;
0236       interceptError2_j = sumxx_j / delta_j;
0237       slope_j = ((sum1_j * sumxy_j) - (sumx_j * sumy_j)) / delta_j;
0238       slopeError2_j = sum1_j / delta_j;
0239     } else
0240       return false;
0241   }
0242 
0243   // from hits on the two chambers, determine radial_intercepts separately and radial_slope together
0244   double sum1_ri = 0.;
0245   double sumx_ri = 0.;
0246   double sumy_ri = 0.;
0247   double sumxx_ri = 0.;
0248   double sumxy_ri = 0.;
0249   for (std::vector<const TransientTrackingRecHit *>::const_iterator hit = hits_i.begin(); hit != hits_i.end(); ++hit) {
0250     double r = (*hit)->globalPosition().perp();
0251     double z = (*hit)->globalPosition().z() - m_Zplane;
0252     sum1_ri += 1.;
0253     sumx_ri += z;
0254     sumy_ri += r;
0255     sumxx_ri += z * z;
0256     sumxy_ri += z * r;
0257   }
0258   double radial_delta_i = (sum1_ri * sumxx_ri) - (sumx_ri * sumx_ri);
0259   if (radial_delta_i == 0.)
0260     return false;
0261   double radial_slope_i = ((sum1_ri * sumxy_ri) - (sumx_ri * sumy_ri)) / radial_delta_i;
0262   double radial_intercept_i =
0263       ((sumxx_ri * sumy_ri) - (sumx_ri * sumxy_ri)) / radial_delta_i + radial_slope_i * (m_iZ - m_Zplane);
0264 
0265   double sum1_rj = 0.;
0266   double sumx_rj = 0.;
0267   double sumy_rj = 0.;
0268   double sumxx_rj = 0.;
0269   double sumxy_rj = 0.;
0270   for (std::vector<const TransientTrackingRecHit *>::const_iterator hit = hits_j.begin(); hit != hits_j.end(); ++hit) {
0271     double r = (*hit)->globalPosition().perp();
0272     double z = (*hit)->globalPosition().z() - m_Zplane;
0273     sum1_rj += 1.;
0274     sumx_rj += z;
0275     sumy_rj += r;
0276     sumxx_rj += z * z;
0277     sumxy_rj += z * r;
0278   }
0279   double radial_delta_j = (sum1_rj * sumxx_rj) - (sumx_rj * sumx_rj);
0280   if (radial_delta_j == 0.)
0281     return false;
0282   double radial_slope_j = ((sum1_rj * sumxy_rj) - (sumx_rj * sumy_rj)) / radial_delta_j;
0283   double radial_intercept_j =
0284       ((sumxx_rj * sumy_rj) - (sumx_rj * sumxy_rj)) / radial_delta_j + radial_slope_j * (m_jZ - m_Zplane);
0285 
0286   double radial_delta = ((sum1_ri + sum1_rj) * (sumxx_ri + sumxx_rj)) - ((sumx_ri + sumx_rj) * (sumx_ri + sumx_rj));
0287   if (radial_delta == 0.)
0288     return false;
0289   double radial_intercept =
0290       (((sumxx_ri + sumxx_rj) * (sumy_ri + sumy_rj)) - ((sumx_ri + sumx_rj) * (sumxy_ri + sumxy_rj))) / radial_delta;
0291   double radial_slope =
0292       (((sum1_ri + sum1_rj) * (sumxy_ri + sumxy_rj)) - ((sumx_ri + sumx_rj) * (sumy_ri + sumy_rj))) / radial_delta;
0293 
0294   if (m_parent->m_makeHistograms) {
0295     m_parent->m_drdz->Fill(radial_slope);
0296   }
0297   if (m_parent->m_maxdrdz > 0. && fabs(radial_slope) > m_parent->m_maxdrdz)
0298     return false;
0299 
0300   double quantity = 0.;
0301   double quantityError2 = 0.;
0302   if (m_parent->m_mode == kModePhiy) {  // phiy comes from track d(rphi)/dz
0303     quantity = (slope_i * radial_intercept_i) - (slope_j * radial_intercept_j);
0304     quantityError2 = (slopeError2_i)*pow(radial_intercept_i, 2) + (slopeError2_j)*pow(radial_intercept_j, 2);
0305   } else if (m_parent->m_mode == kModePhiPos || m_parent->m_mode == kModeRadius) {  // phipos comes from phi intercepts
0306     quantity = intercept_i - intercept_j;
0307     quantityError2 = interceptError2_i + interceptError2_j;
0308   } else if (m_parent->m_mode == kModePhiz) {  // phiz comes from the slope of rphi intercepts
0309     quantity = (intercept_i - intercept_j) * radial_intercept;
0310     quantityError2 = (interceptError2_i + interceptError2_j) * pow(radial_intercept, 2);
0311   } else
0312     assert(false);
0313 
0314   if (quantityError2 == 0.)
0315     return false;
0316 
0317   double slopeResid = ((slope_i * radial_intercept_i) - (slope_j * radial_intercept_j)) * 1000.;
0318   double slopeResidError2 =
0319       ((slopeError2_i)*pow(radial_intercept_i, 2) + (slopeError2_j)*pow(radial_intercept_j, 2)) * 1000. * 1000.;
0320   double offsetResid = (intercept_i - intercept_j) * radial_intercept * 10.;
0321   double offsetResidError2 = (interceptError2_i + interceptError2_j) * pow(radial_intercept, 2) * 10. * 10.;
0322 
0323   if (m_parent->m_truncateSlopeResid > 0. && fabs(slopeResid) > m_parent->m_truncateSlopeResid)
0324     return false;
0325   if (m_parent->m_truncateOffsetResid > 0. && fabs(offsetResid) > m_parent->m_truncateOffsetResid)
0326     return false;
0327 
0328   double weight = 1.;
0329   if (m_parent->m_useTrackWeights)
0330     weight = 1. / quantityError2;
0331 
0332   // fill the running sums for this CSCPairResidualsConstraint
0333   m_sumN += 1;
0334   m_sum1 += weight;
0335   m_sumx += weight * (radial_intercept - m_averageRadius);
0336   m_sumy += weight * quantity;
0337   m_sumxx += weight * pow(radial_intercept - m_averageRadius, 2);
0338   m_sumyy += weight * quantity * quantity;
0339   m_sumxy += weight * (radial_intercept - m_averageRadius) * quantity;
0340 
0341   if (m_parent->m_makeHistograms) {
0342     double rphi_slope_i = slope_i * radial_intercept_i;
0343     double rphi_slope_j = slope_j * radial_intercept_j;
0344 
0345     if (m_parent->m_slopeFromTrackRefit) {
0346       m_parent->m_slope->Fill(rphi_slope_i);  // == rphi_slope_j
0347 
0348       if (m_id_i.endcap() == 1 && m_id_i.station() == 4)
0349         m_parent->m_slope_MEp4->Fill(rphi_slope_i);
0350       if (m_id_i.endcap() == 1 && m_id_i.station() == 3)
0351         m_parent->m_slope_MEp3->Fill(rphi_slope_i);
0352       if (m_id_i.endcap() == 1 && m_id_i.station() == 2)
0353         m_parent->m_slope_MEp2->Fill(rphi_slope_i);
0354       if (m_id_i.endcap() == 1 && m_id_i.station() == 1)
0355         m_parent->m_slope_MEp1->Fill(rphi_slope_i);
0356       if (m_id_i.endcap() == 2 && m_id_i.station() == 1)
0357         m_parent->m_slope_MEm1->Fill(rphi_slope_i);
0358       if (m_id_i.endcap() == 2 && m_id_i.station() == 2)
0359         m_parent->m_slope_MEm2->Fill(rphi_slope_i);
0360       if (m_id_i.endcap() == 2 && m_id_i.station() == 3)
0361         m_parent->m_slope_MEm3->Fill(rphi_slope_i);
0362       if (m_id_i.endcap() == 2 && m_id_i.station() == 4)
0363         m_parent->m_slope_MEm4->Fill(rphi_slope_i);
0364     } else {
0365       m_parent->m_slope->Fill(rphi_slope_i);
0366       m_parent->m_slope->Fill(rphi_slope_j);
0367 
0368       if (m_id_i.endcap() == 1 && m_id_i.station() == 4) {
0369         m_parent->m_slope_MEp4->Fill(rphi_slope_i);
0370         m_parent->m_slope_MEp4->Fill(rphi_slope_j);
0371       }
0372       if (m_id_i.endcap() == 1 && m_id_i.station() == 3) {
0373         m_parent->m_slope_MEp3->Fill(rphi_slope_i);
0374         m_parent->m_slope_MEp3->Fill(rphi_slope_j);
0375       }
0376       if (m_id_i.endcap() == 1 && m_id_i.station() == 2) {
0377         m_parent->m_slope_MEp2->Fill(rphi_slope_i);
0378         m_parent->m_slope_MEp2->Fill(rphi_slope_j);
0379       }
0380       if (m_id_i.endcap() == 1 && m_id_i.station() == 1) {
0381         m_parent->m_slope_MEp1->Fill(rphi_slope_i);
0382         m_parent->m_slope_MEp1->Fill(rphi_slope_j);
0383       }
0384       if (m_id_i.endcap() == 2 && m_id_i.station() == 1) {
0385         m_parent->m_slope_MEm1->Fill(rphi_slope_i);
0386         m_parent->m_slope_MEm1->Fill(rphi_slope_j);
0387       }
0388       if (m_id_i.endcap() == 2 && m_id_i.station() == 2) {
0389         m_parent->m_slope_MEm2->Fill(rphi_slope_i);
0390         m_parent->m_slope_MEm2->Fill(rphi_slope_j);
0391       }
0392       if (m_id_i.endcap() == 2 && m_id_i.station() == 3) {
0393         m_parent->m_slope_MEm3->Fill(rphi_slope_i);
0394         m_parent->m_slope_MEm3->Fill(rphi_slope_j);
0395       }
0396       if (m_id_i.endcap() == 2 && m_id_i.station() == 4) {
0397         m_parent->m_slope_MEm4->Fill(rphi_slope_i);
0398         m_parent->m_slope_MEm4->Fill(rphi_slope_j);
0399       }
0400     }
0401 
0402     m_slopeResiduals->Fill(slopeResid);
0403     m_offsetResiduals->Fill(offsetResid);
0404     m_radial->Fill(radial_intercept);
0405 
0406     m_parent->m_slopeResiduals->Fill(slopeResid);
0407     m_parent->m_slopeResiduals_weighted->Fill(slopeResid, 1. / slopeResidError2);
0408     m_parent->m_slopeResiduals_normalized->Fill(slopeResid / sqrt(slopeResidError2));
0409 
0410     m_parent->m_offsetResiduals->Fill(offsetResid);
0411     m_parent->m_offsetResiduals_weighted->Fill(offsetResid, 1. / offsetResidError2);
0412     m_parent->m_offsetResiduals_normalized->Fill(offsetResid / sqrt(offsetResidError2));
0413 
0414     double ringbin = 0;
0415     if (m_id_i.endcap() == 2 && m_id_i.station() == 4 && m_id_i.ring() == 2)
0416       ringbin = 1.5;
0417     else if (m_id_i.endcap() == 2 && m_id_i.station() == 4 && m_id_i.ring() == 1)
0418       ringbin = 2.5;
0419     else if (m_id_i.endcap() == 2 && m_id_i.station() == 3 && m_id_i.ring() == 2)
0420       ringbin = 3.5;
0421     else if (m_id_i.endcap() == 2 && m_id_i.station() == 3 && m_id_i.ring() == 1)
0422       ringbin = 4.5;
0423     else if (m_id_i.endcap() == 2 && m_id_i.station() == 2 && m_id_i.ring() == 2)
0424       ringbin = 5.5;
0425     else if (m_id_i.endcap() == 2 && m_id_i.station() == 2 && m_id_i.ring() == 1)
0426       ringbin = 6.5;
0427     else if (m_id_i.endcap() == 2 && m_id_i.station() == 1 && m_id_i.ring() == 3)
0428       ringbin = 7.5;
0429     else if (m_id_i.endcap() == 2 && m_id_i.station() == 1 && m_id_i.ring() == 2)
0430       ringbin = 8.5;
0431     else if (m_id_i.endcap() == 2 && m_id_i.station() == 1 && m_id_i.ring() == 1)
0432       ringbin = 9.5;
0433     else if (m_id_i.endcap() == 2 && m_id_i.station() == 1 && m_id_i.ring() == 4)
0434       ringbin = 10.5;
0435     else if (m_id_i.endcap() == 1 && m_id_i.station() == 1 && m_id_i.ring() == 4)
0436       ringbin = 11.5;
0437     else if (m_id_i.endcap() == 1 && m_id_i.station() == 1 && m_id_i.ring() == 1)
0438       ringbin = 12.5;
0439     else if (m_id_i.endcap() == 1 && m_id_i.station() == 1 && m_id_i.ring() == 2)
0440       ringbin = 13.5;
0441     else if (m_id_i.endcap() == 1 && m_id_i.station() == 1 && m_id_i.ring() == 3)
0442       ringbin = 14.5;
0443     else if (m_id_i.endcap() == 1 && m_id_i.station() == 2 && m_id_i.ring() == 1)
0444       ringbin = 15.5;
0445     else if (m_id_i.endcap() == 1 && m_id_i.station() == 2 && m_id_i.ring() == 2)
0446       ringbin = 16.5;
0447     else if (m_id_i.endcap() == 1 && m_id_i.station() == 3 && m_id_i.ring() == 1)
0448       ringbin = 17.5;
0449     else if (m_id_i.endcap() == 1 && m_id_i.station() == 3 && m_id_i.ring() == 2)
0450       ringbin = 18.5;
0451     else if (m_id_i.endcap() == 1 && m_id_i.station() == 4 && m_id_i.ring() == 1)
0452       ringbin = 19.5;
0453     else if (m_id_i.endcap() == 1 && m_id_i.station() == 4 && m_id_i.ring() == 2)
0454       ringbin = 20.5;
0455     m_parent->m_occupancy->Fill(m_id_i.chamber() + 0.5, ringbin);
0456   }
0457 
0458   return true;
0459 }
0460 
0461 bool CSCPairResidualsConstraint::dphidzFromTrack(const std::vector<TrajectoryMeasurement> &measurements,
0462                                                  const reco::TransientTrack &track,
0463                                                  const TrackTransformer *trackTransformer,
0464                                                  double &dphidz) {
0465   // make a list of hits on all chambers *other* than the ones associated with this constraint
0466   std::map<int, int> stations;
0467   TransientTrackingRecHit::ConstRecHitContainer cscHits;
0468   for (std::vector<TrajectoryMeasurement>::const_iterator measurement = measurements.begin();
0469        measurement != measurements.end();
0470        ++measurement) {
0471     DetId id = measurement->recHit()->geographicalId();
0472     if (id.det() == DetId::Muon && id.subdetId() == MuonSubdetId::CSC) {
0473       CSCDetId cscid(id.rawId());
0474       CSCDetId chamberId(cscid.endcap(), cscid.station(), cscid.ring(), cscid.chamber(), 0);
0475       if (m_parent->m_combineME11 && cscid.station() == 1 && cscid.ring() == 4)
0476         chamberId = CSCDetId(cscid.endcap(), 1, 1, cscid.chamber(), 0);
0477 
0478       if (chamberId != m_id_i && chamberId != m_id_j) {
0479         int station = (cscid.endcap() == 1 ? 1 : -1) * cscid.station();
0480         if (stations.find(station) == stations.end()) {
0481           stations[station] = 0;
0482         }
0483         stations[station]++;
0484 
0485         cscHits.push_back(measurement->recHit());
0486       }
0487     }
0488   }
0489 
0490   // for the fit to be reliable, it needs to cross multiple stations
0491   int numStations = 0;
0492   for (std::map<int, int>::const_iterator station = stations.begin(); station != stations.end(); ++station) {
0493     if (station->second >= m_parent->m_minHitsPerChamber) {
0494       numStations++;
0495     }
0496   }
0497 
0498   if (numStations >= m_parent->m_minStationsInTrackRefits) {
0499     // refit the track with these hits
0500     std::vector<Trajectory> trajectories = trackTransformer->transform(track, cscHits);
0501 
0502     if (!trajectories.empty()) {
0503       const std::vector<TrajectoryMeasurement> &measurements2 = trajectories.begin()->measurements();
0504 
0505       // find the closest TSOS to the Z plane (on both sides)
0506       bool found_plus = false;
0507       bool found_minus = false;
0508       TrajectoryStateOnSurface tsos_plus, tsos_minus;
0509       for (std::vector<TrajectoryMeasurement>::const_iterator measurement = measurements2.begin();
0510            measurement != measurements2.end();
0511            ++measurement) {
0512         double z = measurement->recHit()->globalPosition().z();
0513         if (z > m_Zplane) {
0514           if (!found_plus || fabs(z - m_Zplane) < fabs(tsos_plus.globalPosition().z() - m_Zplane)) {
0515             tsos_plus = TrajectoryStateCombiner().combine(measurement->forwardPredictedState(),
0516                                                           measurement->backwardPredictedState());
0517           }
0518           if (tsos_plus.isValid())
0519             found_plus = true;
0520         } else {
0521           if (!found_minus || fabs(z - m_Zplane) < fabs(tsos_minus.globalPosition().z() - m_Zplane)) {
0522             tsos_minus = TrajectoryStateCombiner().combine(measurement->forwardPredictedState(),
0523                                                            measurement->backwardPredictedState());
0524           }
0525           if (tsos_minus.isValid())
0526             found_minus = true;
0527         }
0528       }
0529 
0530       // propagate from the closest TSOS to the Z plane (from both sides, if possible)
0531       TrajectoryStateOnSurface from_plus, from_minus;
0532       if (found_plus) {
0533         from_plus = m_propagator->propagate(tsos_plus, *m_Zsurface);
0534       }
0535       if (found_minus) {
0536         from_minus = m_propagator->propagate(tsos_minus, *m_Zsurface);
0537       }
0538 
0539       // if you have two sides, merge them
0540       TrajectoryStateOnSurface merged;
0541       if (found_plus && from_plus.isValid() && found_minus && from_minus.isValid()) {
0542         merged = TrajectoryStateCombiner().combine(from_plus, from_minus);
0543       } else if (found_plus && from_plus.isValid()) {
0544         merged = from_plus;
0545       } else if (found_minus && from_minus.isValid()) {
0546         merged = from_minus;
0547       } else
0548         return false;
0549 
0550       // if, after all that, we have a good fit-and-propagation, report the direction
0551       if (merged.isValid()) {
0552         double angle = merged.globalPosition().phi() + M_PI / 2.;
0553         GlobalVector direction = merged.globalDirection();
0554         double dxdz = direction.x() / direction.z();
0555         double dydz = direction.y() / direction.z();
0556         dphidz = (dxdz * cos(angle) + dydz * sin(angle)) / merged.globalPosition().perp();
0557         return true;
0558       }
0559 
0560     }  // end if refit successful
0561   }    // end if enough hits
0562   return false;
0563 }
0564 
0565 void CSCPairResidualsConstraint::write(std::ofstream &output) {
0566   output << std::setprecision(14) << std::fixed;
0567   output << "CSCPairResidualsConstraint " << m_identifier << " " << i() << " " << j() << " " << m_sumN << " " << m_sum1
0568          << " " << m_sumx << " " << m_sumy << " " << m_sumxx << " " << m_sumyy << " " << m_sumxy << " EOLN"
0569          << std::endl;
0570 }
0571 
0572 void CSCPairResidualsConstraint::read(std::vector<std::ifstream *> &input, std::vector<std::string> &filenames) {
0573   m_sumN = 0;
0574   m_sum1 = 0.;
0575   m_sumx = 0.;
0576   m_sumy = 0.;
0577   m_sumxx = 0.;
0578   m_sumyy = 0.;
0579   m_sumxy = 0.;
0580 
0581   std::vector<std::ifstream *>::const_iterator inputiter = input.begin();
0582   std::vector<std::string>::const_iterator filename = filenames.begin();
0583   for (; inputiter != input.end(); ++inputiter, ++filename) {
0584     int linenumber = 0;
0585     bool touched = false;
0586     while (!(*inputiter)->eof()) {
0587       linenumber++;
0588       std::string name, eoln;
0589       unsigned int identifier;
0590       int i, j;
0591       int sumN;
0592       double sum1, sumx, sumy, sumxx, sumyy, sumxy;
0593 
0594       (**inputiter) >> name >> identifier >> i >> j >> sumN >> sum1 >> sumx >> sumy >> sumxx >> sumyy >> sumxy >> eoln;
0595 
0596       if (!(*inputiter)->eof() && (name != "CSCPairResidualsConstraint" || eoln != "EOLN"))
0597         throw cms::Exception("CorruptTempFile")
0598             << "Temporary file " << *filename << " is incorrectly formatted on line " << linenumber << std::endl;
0599 
0600       if (identifier == m_identifier) {
0601         if (i != m_i || j != m_j)
0602           throw cms::Exception("CorruptTempFile")
0603               << "Wrong (i,j) for CSCPairResidualsConstraint " << m_identifier << " (" << m_i << "," << m_j
0604               << ") in file " << *filename << " on line " << linenumber << std::endl;
0605         touched = true;
0606 
0607         m_sumN += sumN;
0608         m_sum1 += sum1;
0609         m_sumx += sumx;
0610         m_sumy += sumy;
0611         m_sumxx += sumxx;
0612         m_sumyy += sumyy;
0613         m_sumxy += sumxy;
0614       }
0615     }
0616 
0617     (*inputiter)->clear();
0618     (*inputiter)->seekg(0, std::ios::beg);
0619 
0620     if (!touched)
0621       throw cms::Exception("CorruptTempFile")
0622           << "CSCPairResidualsConstraint " << m_identifier << " is missing from file " << *filename << std::endl;
0623   }
0624 }
0625 
0626 void CSCPairResidualsConstraint::calculatePhi(
0627     const TransientTrackingRecHit *hit, double &phi, double &phierr2, bool doRphi, bool globalPhi) {
0628   align::LocalPoint pos = hit->localPosition();
0629   DetId id = hit->geographicalId();
0630   CSCDetId cscid = CSCDetId(id.rawId());
0631 
0632   double r = 0.;
0633   if (globalPhi) {
0634     phi = hit->globalPosition().phi();
0635     r = hit->globalPosition().perp();
0636 
0637     //     double sinAngle = sin(phi);
0638     //     double cosAngle = cos(phi);
0639     //     double xx = hit->globalPositionError().cxx();
0640     //     double xy = hit->globalPositionError().cyx();
0641     //     double yy = hit->globalPositionError().cyy();
0642     //     phierr2 = (xx*cosAngle*cosAngle + 2.*xy*sinAngle*cosAngle + yy*sinAngle*sinAngle) / (r*r);
0643   } else {
0644     // these constants are related to the way CSC chambers are built--- really constant!
0645     const double R_ME11 = 181.5;
0646     const double R_ME12 = 369.7;
0647     const double R_ME21 = 242.7;
0648     const double R_ME31 = 252.7;
0649     const double R_ME41 = 262.65;
0650     const double R_MEx2 = 526.5;
0651 
0652     double R = 0.;
0653     if (cscid.station() == 1 && (cscid.ring() == 1 || cscid.ring() == 4))
0654       R = R_ME11;
0655     else if (cscid.station() == 1 && cscid.ring() == 2)
0656       R = R_ME12;
0657     else if (cscid.station() == 2 && cscid.ring() == 1)
0658       R = R_ME21;
0659     else if (cscid.station() == 3 && cscid.ring() == 1)
0660       R = R_ME31;
0661     else if (cscid.station() == 4 && cscid.ring() == 1)
0662       R = R_ME41;
0663     else if (cscid.station() > 1 && cscid.ring() == 2)
0664       R = R_MEx2;
0665     else
0666       assert(false);
0667     r = (pos.y() + R);
0668 
0669     phi = atan2(pos.x(), r);
0670 
0671     if (cscid.endcap() == 1 && cscid.station() >= 3)
0672       phi *= -1;
0673     else if (cscid.endcap() == 2 && cscid.station() <= 2)
0674       phi *= -1;
0675   }
0676 
0677   int strip = m_cscGeometry->layer(id)->geometry()->nearestStrip(pos);
0678   double angle = m_cscGeometry->layer(id)->geometry()->stripAngle(strip) - M_PI / 2.;
0679   double sinAngle = sin(angle);
0680   double cosAngle = cos(angle);
0681   double xx = hit->localPositionError().xx();
0682   double xy = hit->localPositionError().xy();
0683   double yy = hit->localPositionError().yy();
0684   phierr2 = (xx * cosAngle * cosAngle + 2. * xy * sinAngle * cosAngle + yy * sinAngle * sinAngle) / (r * r);
0685 
0686   if (doRphi) {
0687     phi *= r;
0688     phierr2 *= r * r;
0689   }
0690 }
0691 
0692 bool CSCPairResidualsConstraint::isFiducial(std::vector<const TransientTrackingRecHit *> &hits, bool is_i) {
0693   // these constants are related to the way CSC chambers are built--- really constant!
0694   const double cut_ME11 = 0.086;
0695   const double cut_ME12 = 0.090;
0696   const double cut_MEx1 = 0.180;
0697   const double cut_MEx2 = 0.090;
0698 
0699   double sum1 = 0.;
0700   double sumx = 0.;
0701   double sumy = 0.;
0702   double sumxx = 0.;
0703   double sumxy = 0.;
0704   for (std::vector<const TransientTrackingRecHit *>::const_iterator hit = hits.begin(); hit != hits.end(); ++hit) {
0705     double phi, phierr2;
0706     calculatePhi(*hit, phi, phierr2);
0707     double z = (is_i ? m_cscGeometry->idToDet(m_id_i)->surface() : m_cscGeometry->idToDet(m_id_j)->surface())
0708                    .toLocal((*hit)->globalPosition())
0709                    .z();
0710 
0711     if (m_parent->m_makeHistograms) {
0712       if (m_id_i.station() == 1 && (m_id_i.ring() == 1 || m_id_i.ring() == 4)) {
0713         m_parent->m_fiducial_ME11->Fill(fabs(phi), sqrt(phierr2));
0714       } else if (m_id_i.station() == 1 && m_id_i.ring() == 2) {
0715         m_parent->m_fiducial_ME12->Fill(fabs(phi), sqrt(phierr2));
0716       } else if (m_id_i.station() > 1 && m_id_i.ring() == 1) {
0717         m_parent->m_fiducial_MEx1->Fill(fabs(phi), sqrt(phierr2));
0718       } else if (m_id_i.station() > 1 && m_id_i.ring() == 2) {
0719         m_parent->m_fiducial_MEx2->Fill(fabs(phi), sqrt(phierr2));
0720       }
0721     }
0722 
0723     double weight = 1.;
0724     if (m_parent->m_useHitWeights)
0725       weight = 1. / phierr2;
0726     sum1 += weight;
0727     sumx += weight * z;
0728     sumy += weight * phi;
0729     sumxx += weight * z * z;
0730     sumxy += weight * z * phi;
0731   }
0732   double delta = (sum1 * sumxx) - (sumx * sumx);
0733   if (delta == 0.)
0734     return false;
0735   double intercept = ((sumxx * sumy) - (sumx * sumxy)) / delta;
0736   double slope = ((sum1 * sumxy) - (sumx * sumy)) / delta;
0737 
0738   double phi1 = intercept + slope * (is_i ? m_iZ1 : m_jZ1);
0739   double phi6 = intercept + slope * (is_i ? m_iZ6 : m_jZ6);
0740 
0741   if (m_id_i.station() == 1 && (m_id_i.ring() == 1 || m_id_i.ring() == 4)) {
0742     return (fabs(phi1) < cut_ME11 && fabs(phi6) < cut_ME11);
0743   } else if (m_id_i.station() == 1 && m_id_i.ring() == 2) {
0744     return (fabs(phi1) < cut_ME12 && fabs(phi6) < cut_ME12);
0745   } else if (m_id_i.station() > 1 && m_id_i.ring() == 1) {
0746     return (fabs(phi1) < cut_MEx1 && fabs(phi6) < cut_MEx1);
0747   } else if (m_id_i.station() > 1 && m_id_i.ring() == 2) {
0748     return (fabs(phi1) < cut_MEx2 && fabs(phi6) < cut_MEx2);
0749   } else
0750     assert(false);
0751 }
[ 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!