Reference Manual

Go to the documentation of this file.
00001 #include <iomanip>
00002 #include "TrackingTools/TrackFitters/interface/TrajectoryStateCombiner.h"
00003 #include "FWCore/ServiceRegistry/interface/Service.h"
00004 #include "CommonTools/UtilAlgos/interface/TFileService.h"
00005 #include "TrackingTools/PatternTools/interface/Trajectory.h"
00006 #include "Alignment/MuonAlignmentAlgorithms/plugins/CSCOverlapsAlignmentAlgorithm.h"
00007 #include "Alignment/MuonAlignmentAlgorithms/plugins/CSCPairResidualsConstraint.h"
00008 
00009 double CSCPairResidualsConstraint::value() const {
00010   double delta = (m_sum1*m_sumxx) - (m_sumx*m_sumx);
00011   assert(delta > 0.);
00012   if (m_parent->m_mode == kModePhiy  ||  m_parent->m_mode == kModePhiPos  ||  m_parent->m_mode == kModeRadius) {
00013     return ((m_sumxx*m_sumy) - (m_sumx*m_sumxy))/delta;
00014   }
00015   else if (m_parent->m_mode == kModePhiz) {
00016     return ((m_sum1*m_sumxy) - (m_sumx*m_sumy))/delta;
00017   }
00018   else assert(false);
00019 }
00020 
00021 double CSCPairResidualsConstraint::error() const {
00022   if (m_parent->m_errorFromRMS) {
00023     assert(m_sum1 > 0.);
00024     return sqrt((m_sumyy/m_sum1) - pow(m_sumy/m_sum1, 2))/sqrt(m_sumN);
00025   }
00026   else {
00027     double delta = (m_sum1*m_sumxx) - (m_sumx*m_sumx);
00028     assert(delta > 0.);
00029     if (m_parent->m_mode == kModePhiy  ||  m_parent->m_mode == kModePhiPos  ||  m_parent->m_mode == kModeRadius) {
00030       return sqrt(m_sumxx/delta);
00031     }
00032     else if (m_parent->m_mode == kModePhiz) {
00033       return sqrt(m_sum1/delta);
00034     }
00035     else assert(false);
00036   }
00037 }
00038 
00039 bool CSCPairResidualsConstraint::valid() const {
00040    return (m_sumN >= m_parent->m_minTracksPerOverlap);
00041 }
00042 
00043 void CSCPairResidualsConstraint::configure(CSCOverlapsAlignmentAlgorithm *parent) {
00044   m_parent = parent;
00045 
00046   if (m_parent->m_makeHistograms) {
00047     edm::Service<TFileService> tFileService;
00048 
00049     std::stringstream name, name2, name3, title;
00050     title << "i =" << m_id_i << " j =" << m_id_j;
00051 
00052     name << "slopeResiduals_" << m_identifier;
00053     m_slopeResiduals = tFileService->make<TH1F>(name.str().c_str(), title.str().c_str(), 300, -30., 30.);
00054 
00055     name2 << "offsetResiduals_" << m_identifier;
00056     m_offsetResiduals = tFileService->make<TH1F>(name2.str().c_str(), title.str().c_str(), 300, -30., 30.);
00057 
00058     name3 << "radial_" << m_identifier;
00059     m_radial = tFileService->make<TH1F>(name3.str().c_str(), title.str().c_str(), 700, 0., 700.);
00060   }
00061   else {
00062     m_slopeResiduals = NULL;
00063     m_offsetResiduals = NULL;
00064     m_radial = NULL;
00065   }
00066 }
00067 
00068 void CSCPairResidualsConstraint::setZplane(const CSCGeometry *cscGeometry) {
00069   m_cscGeometry = cscGeometry;
00070 
00071   m_Zplane = (m_cscGeometry->idToDet(m_id_i)->surface().position().z() + m_cscGeometry->idToDet(m_id_j)->surface().position().z()) / 2.;
00072   m_averageRadius = (m_cscGeometry->idToDet(m_id_i)->surface().position().perp() + m_cscGeometry->idToDet(m_id_j)->surface().position().perp()) / 2.;
00073 
00074   m_iZ = m_cscGeometry->idToDet(m_id_i)->surface().position().z();
00075   m_jZ = m_cscGeometry->idToDet(m_id_j)->surface().position().z();
00076 
00077   CSCDetId i1(m_id_i.endcap(), m_id_i.station(), m_id_i.ring(), m_id_i.chamber(), 1);
00078   CSCDetId i6(m_id_i.endcap(), m_id_i.station(), m_id_i.ring(), m_id_i.chamber(), 6);
00079   CSCDetId j1(m_id_j.endcap(), m_id_j.station(), m_id_j.ring(), m_id_j.chamber(), 1);
00080   CSCDetId j6(m_id_j.endcap(), m_id_j.station(), m_id_j.ring(), m_id_j.chamber(), 6);
00081   m_iZ1 = m_cscGeometry->idToDet(m_id_i)->surface().toLocal(m_cscGeometry->idToDet(i1)->surface().position()).z();
00082   m_iZ6 = m_cscGeometry->idToDet(m_id_i)->surface().toLocal(m_cscGeometry->idToDet(i6)->surface().position()).z();
00083   m_jZ1 = m_cscGeometry->idToDet(m_id_j)->surface().toLocal(m_cscGeometry->idToDet(j1)->surface().position()).z();
00084   m_jZ6 = m_cscGeometry->idToDet(m_id_j)->surface().toLocal(m_cscGeometry->idToDet(j6)->surface().position()).z();
00085 
00086   m_Zsurface = Plane::build(Plane::PositionType(0., 0., m_Zplane), Plane::RotationType());
00087 }
00088 
00089 void CSCPairResidualsConstraint::setPropagator(const Propagator *propagator) {
00090   m_propagator = propagator;
00091 }
00092 
00093 bool CSCPairResidualsConstraint::addTrack(const std::vector<TrajectoryMeasurement> &measurements, const reco::TransientTrack &track, const TrackTransformer *trackTransformer) {
00094   std::vector<const TransientTrackingRecHit*> hits_i;
00095   std::vector<const TransientTrackingRecHit*> hits_j;
00096 
00097   for (std::vector<TrajectoryMeasurement>::const_iterator measurement = measurements.begin();  measurement != measurements.end();  ++measurement) {
00098     const TransientTrackingRecHit *hit = &*(measurement->recHit());
00099 
00100     DetId id = hit->geographicalId();
00101     if (id.det() == DetId::Muon  &&  id.subdetId() == MuonSubdetId::CSC) {
00102       CSCDetId cscid(id.rawId());
00103       CSCDetId chamberId(cscid.endcap(), cscid.station(), cscid.ring(), cscid.chamber(), 0);
00104       if (m_parent->m_combineME11  &&  cscid.station() == 1  &&  cscid.ring() == 4) chamberId = CSCDetId(cscid.endcap(), 1, 1, cscid.chamber(), 0);
00105 
00106       if (chamberId == m_id_i) hits_i.push_back(hit);
00107       if (chamberId == m_id_j) hits_j.push_back(hit);
00108     }
00109   }
00110 
00111   if (m_parent->m_makeHistograms) {
00112     m_parent->m_hitsPerChamber->Fill(hits_i.size());
00113     m_parent->m_hitsPerChamber->Fill(hits_j.size());
00114   }
00115 
00116   // require minimum number of hits (if the requirement is too low (~2), some NANs might result...)
00117   if (int(hits_i.size()) < m_parent->m_minHitsPerChamber  ||  int(hits_j.size()) < m_parent->m_minHitsPerChamber) return false;
00118 
00119   // maybe require segments to be fiducial
00120   if (m_parent->m_fiducial  &&  !(isFiducial(hits_i, true)  &&  isFiducial(hits_j, false))) return false;
00121 
00122   double intercept_i = 0.;
00123   double interceptError2_i = 0.;
00124   double slope_i = 0.;
00125   double slopeError2_i = 0.;
00126   double intercept_j = 0.;
00127   double interceptError2_j = 0.;
00128   double slope_j = 0.;
00129   double slopeError2_j = 0.;
00130 
00131   // if slopeFromTrackRefit, then you'll need to refit the whole track without this station's hits;
00132   // need at least two other stations for that to be reliable
00133   if (m_parent->m_slopeFromTrackRefit) {
00134     double dphidz;
00135     if (dphidzFromTrack(measurements, track, trackTransformer, dphidz)) {
00136       double sum1_i = 0.;
00137       double sumy_i = 0.;
00138       for (std::vector<const TransientTrackingRecHit*>::const_iterator hit = hits_i.begin();  hit != hits_i.end();  ++hit) {
00139         double phi, phierr2;
00140         calculatePhi(*hit, phi, phierr2, false, true);
00141         double z = (*hit)->globalPosition().z() - m_Zplane;
00142 
00143         double weight = 1.;
00144         if (m_parent->m_useHitWeights) weight = 1./phierr2;
00145         sum1_i += weight;
00146         sumy_i += weight * (phi - z*dphidz);
00147       }
00148 
00149       double sum1_j = 0.;
00150       double sumy_j = 0.;
00151       for (std::vector<const TransientTrackingRecHit*>::const_iterator hit = hits_j.begin();  hit != hits_j.end();  ++hit) {
00152         double phi, phierr2;
00153         calculatePhi(*hit, phi, phierr2, false, true);
00154         double z = (*hit)->globalPosition().z() - m_Zplane;
00155 
00156         double weight = 1.;
00157         if (m_parent->m_useHitWeights) weight = 1./phierr2;
00158         sum1_j += weight;
00159         sumy_j += weight * (phi - z*dphidz);
00160       }
00161 
00162       if (sum1_i != 0.  &&  sum1_j != 0.) {
00163         slope_i = slope_j = dphidz;
00164 
00165         intercept_i = sumy_i / sum1_i;
00166         interceptError2_i = 1. / sum1_i;
00167 
00168         intercept_j = sumy_j / sum1_j;
00169         interceptError2_j = 1. / sum1_j;
00170       }
00171       else return false;
00172     }
00173   }
00174 
00175   else { // not slopeFromTrackRefit
00176     double sum1_i = 0.;
00177     double sumx_i = 0.;
00178     double sumy_i = 0.;
00179     double sumxx_i = 0.;
00180     double sumxy_i = 0.;
00181     for (std::vector<const TransientTrackingRecHit*>::const_iterator hit = hits_i.begin();  hit != hits_i.end();  ++hit) {
00182       double phi, phierr2;
00183       calculatePhi(*hit, phi, phierr2, false, true);
00184       double z = (*hit)->globalPosition().z() - m_Zplane;
00185 
00186       double weight = 1.;
00187       if (m_parent->m_useHitWeights) weight = 1./phierr2;
00188       sum1_i += weight;
00189       sumx_i += weight * z;
00190       sumy_i += weight * phi;
00191       sumxx_i += weight * z*z;
00192       sumxy_i += weight * z*phi;
00193     }
00194 
00195     double sum1_j = 0.;
00196     double sumx_j = 0.;
00197     double sumy_j = 0.;
00198     double sumxx_j = 0.;
00199     double sumxy_j = 0.;
00200     for (std::vector<const TransientTrackingRecHit*>::const_iterator hit = hits_j.begin();  hit != hits_j.end();  ++hit) {
00201       double phi, phierr2;
00202       calculatePhi(*hit, phi, phierr2, false, true);
00203       double z = (*hit)->globalPosition().z() - m_Zplane;
00204         
00205       double weight = 1.;
00206       if (m_parent->m_useHitWeights) weight = 1./phierr2;
00207       sum1_j += weight;
00208       sumx_j += weight * z;
00209       sumy_j += weight * phi;
00210       sumxx_j += weight * z*z;
00211       sumxy_j += weight * z*phi;
00212     }
00213 
00214     double delta_i = (sum1_i*sumxx_i) - (sumx_i*sumx_i);    
00215     double delta_j = (sum1_j*sumxx_j) - (sumx_j*sumx_j);    
00216     if (delta_i != 0.  &&  delta_j != 0.) {
00217       intercept_i = ((sumxx_i*sumy_i) - (sumx_i*sumxy_i))/delta_i;
00218       interceptError2_i = sumxx_i/delta_i;
00219       slope_i = ((sum1_i*sumxy_i) - (sumx_i*sumy_i))/delta_i;
00220       slopeError2_i = sum1_i/delta_i;
00221 
00222       intercept_j = ((sumxx_j*sumy_j) - (sumx_j*sumxy_j))/delta_j;
00223       interceptError2_j = sumxx_j/delta_j;
00224       slope_j = ((sum1_j*sumxy_j) - (sumx_j*sumy_j))/delta_j;
00225       slopeError2_j = sum1_j/delta_j;
00226     }
00227     else return false;
00228   }
00229 
00230   // from hits on the two chambers, determine radial_intercepts separately and radial_slope together
00231   double sum1_ri = 0.;
00232   double sumx_ri = 0.;
00233   double sumy_ri = 0.;
00234   double sumxx_ri = 0.;
00235   double sumxy_ri = 0.;
00236   for (std::vector<const TransientTrackingRecHit*>::const_iterator hit = hits_i.begin();  hit != hits_i.end();  ++hit) {
00237       double r = (*hit)->globalPosition().perp();
00238       double z = (*hit)->globalPosition().z() - m_Zplane;
00239       sum1_ri += 1.;
00240       sumx_ri += z;
00241       sumy_ri += r;
00242       sumxx_ri += z*z;
00243       sumxy_ri += z*r;
00244   }
00245   double radial_delta_i = (sum1_ri*sumxx_ri) - (sumx_ri*sumx_ri);
00246   if (radial_delta_i == 0.) return false;
00247   double radial_slope_i = ((sum1_ri*sumxy_ri) - (sumx_ri*sumy_ri))/radial_delta_i;
00248   double radial_intercept_i = ((sumxx_ri*sumy_ri) - (sumx_ri*sumxy_ri))/radial_delta_i + radial_slope_i*(m_iZ - m_Zplane);
00249 
00250   double sum1_rj = 0.;
00251   double sumx_rj = 0.;
00252   double sumy_rj = 0.;
00253   double sumxx_rj = 0.;
00254   double sumxy_rj = 0.;
00255   for (std::vector<const TransientTrackingRecHit*>::const_iterator hit = hits_j.begin();  hit != hits_j.end();  ++hit) {
00256       double r = (*hit)->globalPosition().perp();
00257       double z = (*hit)->globalPosition().z() - m_Zplane;
00258       sum1_rj += 1.;
00259       sumx_rj += z;
00260       sumy_rj += r;
00261       sumxx_rj += z*z;
00262       sumxy_rj += z*r;
00263   }
00264   double radial_delta_j = (sum1_rj*sumxx_rj) - (sumx_rj*sumx_rj);
00265   if (radial_delta_j == 0.) return false;
00266   double radial_slope_j = ((sum1_rj*sumxy_rj) - (sumx_rj*sumy_rj))/radial_delta_j;
00267   double radial_intercept_j = ((sumxx_rj*sumy_rj) - (sumx_rj*sumxy_rj))/radial_delta_j + radial_slope_j*(m_jZ - m_Zplane);
00268 
00269   double radial_delta = ((sum1_ri + sum1_rj)*(sumxx_ri + sumxx_rj)) - ((sumx_ri + sumx_rj)*(sumx_ri + sumx_rj));
00270   if (radial_delta == 0.) return false;
00271   double radial_intercept = (((sumxx_ri + sumxx_rj)*(sumy_ri + sumy_rj)) - ((sumx_ri + sumx_rj)*(sumxy_ri + sumxy_rj)))/radial_delta;
00272   double radial_slope = (((sum1_ri + sum1_rj)*(sumxy_ri + sumxy_rj)) - ((sumx_ri + sumx_rj)*(sumy_ri + sumy_rj)))/radial_delta;
00273 
00274   if (m_parent->m_makeHistograms) {
00275     m_parent->m_drdz->Fill(radial_slope);
00276   }
00277   if (m_parent->m_maxdrdz > 0.  &&  fabs(radial_slope) > m_parent->m_maxdrdz) return false;
00278 
00279   double quantity = 0.;
00280   double quantityError2 = 0.;
00281   if (m_parent->m_mode == kModePhiy) {  // phiy comes from track d(rphi)/dz
00282     quantity = (slope_i*radial_intercept_i) - (slope_j*radial_intercept_j);
00283     quantityError2 = (slopeError2_i)*pow(radial_intercept_i, 2) + (slopeError2_j)*pow(radial_intercept_j, 2);
00284   }
00285   else if (m_parent->m_mode == kModePhiPos  ||  m_parent->m_mode == kModeRadius) {  // phipos comes from phi intercepts
00286     quantity = intercept_i - intercept_j;
00287     quantityError2 = interceptError2_i + interceptError2_j;
00288   }
00289   else if (m_parent->m_mode == kModePhiz) {  // phiz comes from the slope of rphi intercepts
00290     quantity = (intercept_i - intercept_j) * radial_intercept;
00291     quantityError2 = (interceptError2_i + interceptError2_j) * pow(radial_intercept, 2);
00292   }
00293   else assert(false);
00294 
00295   if (quantityError2 == 0.) return false;
00296 
00297   double slopeResid = ((slope_i*radial_intercept_i) - (slope_j*radial_intercept_j)) * 1000.;
00298   double slopeResidError2 = ((slopeError2_i)*pow(radial_intercept_i, 2) + (slopeError2_j)*pow(radial_intercept_j, 2)) * 1000. * 1000.;
00299   double offsetResid = (intercept_i - intercept_j) * radial_intercept * 10.;
00300   double offsetResidError2 = (interceptError2_i + interceptError2_j) * pow(radial_intercept, 2) * 10. * 10.;
00301 
00302   if (m_parent->m_truncateSlopeResid > 0.  &&  fabs(slopeResid) > m_parent->m_truncateSlopeResid) return false;
00303   if (m_parent->m_truncateOffsetResid > 0.  &&  fabs(offsetResid) > m_parent->m_truncateOffsetResid) return false;
00304 
00305   double weight = 1.;
00306   if (m_parent->m_useTrackWeights) weight = 1./quantityError2;
00307 
00308   // fill the running sums for this CSCPairResidualsConstraint
00309   m_sumN += 1;
00310   m_sum1 += weight;
00311   m_sumx += weight * (radial_intercept - m_averageRadius);
00312   m_sumy += weight * quantity;
00313   m_sumxx += weight * pow(radial_intercept - m_averageRadius, 2);
00314   m_sumyy += weight * quantity*quantity;
00315   m_sumxy += weight * (radial_intercept - m_averageRadius)*quantity;
00316 
00317   if (m_parent->m_makeHistograms) {
00318     double rphi_slope_i = slope_i * radial_intercept_i;
00319     double rphi_slope_j = slope_j * radial_intercept_j;
00320 
00321     if (m_parent->m_slopeFromTrackRefit) {
00322       m_parent->m_slope->Fill(rphi_slope_i);  // == rphi_slope_j
00323 
00324       if (m_id_i.endcap() == 1  &&  m_id_i.station() == 4) m_parent->m_slope_MEp4->Fill(rphi_slope_i);
00325       if (m_id_i.endcap() == 1  &&  m_id_i.station() == 3) m_parent->m_slope_MEp3->Fill(rphi_slope_i);
00326       if (m_id_i.endcap() == 1  &&  m_id_i.station() == 2) m_parent->m_slope_MEp2->Fill(rphi_slope_i);
00327       if (m_id_i.endcap() == 1  &&  m_id_i.station() == 1) m_parent->m_slope_MEp1->Fill(rphi_slope_i);
00328       if (m_id_i.endcap() == 2  &&  m_id_i.station() == 1) m_parent->m_slope_MEm1->Fill(rphi_slope_i);
00329       if (m_id_i.endcap() == 2  &&  m_id_i.station() == 2) m_parent->m_slope_MEm2->Fill(rphi_slope_i);
00330       if (m_id_i.endcap() == 2  &&  m_id_i.station() == 3) m_parent->m_slope_MEm3->Fill(rphi_slope_i);
00331       if (m_id_i.endcap() == 2  &&  m_id_i.station() == 4) m_parent->m_slope_MEm4->Fill(rphi_slope_i);
00332     }
00333     else {
00334       m_parent->m_slope->Fill(rphi_slope_i); m_parent->m_slope->Fill(rphi_slope_j);
00335 
00336       if (m_id_i.endcap() == 1  &&  m_id_i.station() == 4) { m_parent->m_slope_MEp4->Fill(rphi_slope_i); m_parent->m_slope_MEp4->Fill(rphi_slope_j); }
00337       if (m_id_i.endcap() == 1  &&  m_id_i.station() == 3) { m_parent->m_slope_MEp3->Fill(rphi_slope_i); m_parent->m_slope_MEp3->Fill(rphi_slope_j); }
00338       if (m_id_i.endcap() == 1  &&  m_id_i.station() == 2) { m_parent->m_slope_MEp2->Fill(rphi_slope_i); m_parent->m_slope_MEp2->Fill(rphi_slope_j); }
00339       if (m_id_i.endcap() == 1  &&  m_id_i.station() == 1) { m_parent->m_slope_MEp1->Fill(rphi_slope_i); m_parent->m_slope_MEp1->Fill(rphi_slope_j); }
00340       if (m_id_i.endcap() == 2  &&  m_id_i.station() == 1) { m_parent->m_slope_MEm1->Fill(rphi_slope_i); m_parent->m_slope_MEm1->Fill(rphi_slope_j); }
00341       if (m_id_i.endcap() == 2  &&  m_id_i.station() == 2) { m_parent->m_slope_MEm2->Fill(rphi_slope_i); m_parent->m_slope_MEm2->Fill(rphi_slope_j); }
00342       if (m_id_i.endcap() == 2  &&  m_id_i.station() == 3) { m_parent->m_slope_MEm3->Fill(rphi_slope_i); m_parent->m_slope_MEm3->Fill(rphi_slope_j); }
00343       if (m_id_i.endcap() == 2  &&  m_id_i.station() == 4) { m_parent->m_slope_MEm4->Fill(rphi_slope_i); m_parent->m_slope_MEm4->Fill(rphi_slope_j); }
00344     }
00345 
00346     m_slopeResiduals->Fill(slopeResid);
00347     m_offsetResiduals->Fill(offsetResid);
00348     m_radial->Fill(radial_intercept);
00349 
00350     m_parent->m_slopeResiduals->Fill(slopeResid);
00351     m_parent->m_slopeResiduals_weighted->Fill(slopeResid, 1./slopeResidError2);
00352     m_parent->m_slopeResiduals_normalized->Fill(slopeResid/sqrt(slopeResidError2));
00353 
00354     m_parent->m_offsetResiduals->Fill(offsetResid);
00355     m_parent->m_offsetResiduals_weighted->Fill(offsetResid, 1./offsetResidError2);
00356     m_parent->m_offsetResiduals_normalized->Fill(offsetResid/sqrt(offsetResidError2));
00357 
00358     double ringbin = 0;
00359     if (m_id_i.endcap() == 2  &&  m_id_i.station() == 4  &&  m_id_i.ring() == 2) ringbin = 1.5;
00360     else if (m_id_i.endcap() == 2  &&  m_id_i.station() == 4  &&  m_id_i.ring() == 1) ringbin = 2.5;
00361     else if (m_id_i.endcap() == 2  &&  m_id_i.station() == 3  &&  m_id_i.ring() == 2) ringbin = 3.5;
00362     else if (m_id_i.endcap() == 2  &&  m_id_i.station() == 3  &&  m_id_i.ring() == 1) ringbin = 4.5;
00363     else if (m_id_i.endcap() == 2  &&  m_id_i.station() == 2  &&  m_id_i.ring() == 2) ringbin = 5.5;
00364     else if (m_id_i.endcap() == 2  &&  m_id_i.station() == 2  &&  m_id_i.ring() == 1) ringbin = 6.5;
00365     else if (m_id_i.endcap() == 2  &&  m_id_i.station() == 1  &&  m_id_i.ring() == 3) ringbin = 7.5;
00366     else if (m_id_i.endcap() == 2  &&  m_id_i.station() == 1  &&  m_id_i.ring() == 2) ringbin = 8.5;
00367     else if (m_id_i.endcap() == 2  &&  m_id_i.station() == 1  &&  m_id_i.ring() == 1) ringbin = 9.5;
00368     else if (m_id_i.endcap() == 2  &&  m_id_i.station() == 1  &&  m_id_i.ring() == 4) ringbin = 10.5;
00369     else if (m_id_i.endcap() == 1  &&  m_id_i.station() == 1  &&  m_id_i.ring() == 4) ringbin = 11.5;
00370     else if (m_id_i.endcap() == 1  &&  m_id_i.station() == 1  &&  m_id_i.ring() == 1) ringbin = 12.5;
00371     else if (m_id_i.endcap() == 1  &&  m_id_i.station() == 1  &&  m_id_i.ring() == 2) ringbin = 13.5;
00372     else if (m_id_i.endcap() == 1  &&  m_id_i.station() == 1  &&  m_id_i.ring() == 3) ringbin = 14.5;
00373     else if (m_id_i.endcap() == 1  &&  m_id_i.station() == 2  &&  m_id_i.ring() == 1) ringbin = 15.5;
00374     else if (m_id_i.endcap() == 1  &&  m_id_i.station() == 2  &&  m_id_i.ring() == 2) ringbin = 16.5;
00375     else if (m_id_i.endcap() == 1  &&  m_id_i.station() == 3  &&  m_id_i.ring() == 1) ringbin = 17.5;
00376     else if (m_id_i.endcap() == 1  &&  m_id_i.station() == 3  &&  m_id_i.ring() == 2) ringbin = 18.5;
00377     else if (m_id_i.endcap() == 1  &&  m_id_i.station() == 4  &&  m_id_i.ring() == 1) ringbin = 19.5;
00378     else if (m_id_i.endcap() == 1  &&  m_id_i.station() == 4  &&  m_id_i.ring() == 2) ringbin = 20.5;
00379     m_parent->m_occupancy->Fill(m_id_i.chamber() + 0.5, ringbin);
00380   }
00381 
00382   return true;
00383 }
00384 
00385 bool CSCPairResidualsConstraint::dphidzFromTrack(const std::vector<TrajectoryMeasurement> &measurements, const reco::TransientTrack &track, const TrackTransformer *trackTransformer, double &dphidz) {
00386   // make a list of hits on all chambers *other* than the ones associated with this constraint
00387   std::map<int,int> stations;
00388   int total = 0;
00389   TransientTrackingRecHit::ConstRecHitContainer cscHits;
00390   for (std::vector<TrajectoryMeasurement>::const_iterator measurement = measurements.begin();  measurement != measurements.end();  ++measurement) {
00391     DetId id = measurement->recHit()->geographicalId();
00392     if (id.det() == DetId::Muon  &&  id.subdetId() == MuonSubdetId::CSC) {
00393       CSCDetId cscid(id.rawId());
00394       CSCDetId chamberId(cscid.endcap(), cscid.station(), cscid.ring(), cscid.chamber(), 0);
00395       if (m_parent->m_combineME11  &&  cscid.station() == 1  &&  cscid.ring() == 4) chamberId = CSCDetId(cscid.endcap(), 1, 1, cscid.chamber(), 0);
00396 
00397       if (chamberId != m_id_i  &&  chamberId != m_id_j) {
00398         int station = (cscid.endcap() == 1 ? 1 : -1) * cscid.station();
00399         if (stations.find(station) == stations.end()) {
00400           stations[station] = 0;
00401         }
00402         stations[station]++;
00403         total++;
00404 
00405         cscHits.push_back(measurement->recHit());
00406       }
00407     }
00408   }
00409 
00410   // for the fit to be reliable, it needs to cross multiple stations
00411   int numStations = 0;
00412   for (std::map<int,int>::const_iterator station = stations.begin();  station != stations.end();  ++station) {
00413     if (station->second >= m_parent->m_minHitsPerChamber) {
00414       numStations++;
00415     }
00416   }
00417 
00418   if (numStations >= m_parent->m_minStationsInTrackRefits) {
00419     // refit the track with these hits
00420     std::vector<Trajectory> trajectories = trackTransformer->transform(track, cscHits);
00421 
00422     if (!trajectories.empty()) {
00423       const std::vector<TrajectoryMeasurement> &measurements2 = trajectories.begin()->measurements();
00424 
00425       // find the closest TSOS to the Z plane (on both sides)
00426       bool found_plus = false;
00427       bool found_minus = false;
00428       TrajectoryStateOnSurface tsos_plus, tsos_minus;
00429       for (std::vector<TrajectoryMeasurement>::const_iterator measurement = measurements2.begin();  measurement != measurements2.end();  ++measurement) {
00430         double z = measurement->recHit()->globalPosition().z();
00431         if (z > m_Zplane) {
00432           if (!found_plus  ||  fabs(z - m_Zplane) < fabs(tsos_plus.globalPosition().z() - m_Zplane)) {
00433             tsos_plus = TrajectoryStateCombiner().combine(measurement->forwardPredictedState(), measurement->backwardPredictedState());
00434           }
00435           if (tsos_plus.isValid()) found_plus = true;
00436         }
00437         else {
00438           if (!found_minus  ||  fabs(z - m_Zplane) < fabs(tsos_minus.globalPosition().z() - m_Zplane)) {
00439             tsos_minus = TrajectoryStateCombiner().combine(measurement->forwardPredictedState(), measurement->backwardPredictedState());
00440           }
00441           if (tsos_minus.isValid()) found_minus = true;
00442         }
00443       }
00444 
00445       // propagate from the closest TSOS to the Z plane (from both sides, if possible)
00446       TrajectoryStateOnSurface from_plus, from_minus;
00447       if (found_plus) {
00448         from_plus = m_propagator->propagate(tsos_plus, *m_Zsurface);
00449       }
00450       if (found_minus) {
00451         from_minus = m_propagator->propagate(tsos_minus, *m_Zsurface);
00452       }
00453 
00454       // if you have two sides, merge them
00455       TrajectoryStateOnSurface merged;
00456       if (found_plus  &&  from_plus.isValid()  &&  found_minus  &&  from_minus.isValid()) {
00457         merged = TrajectoryStateCombiner().combine(from_plus, from_minus);
00458       }
00459       else if (found_plus  &&  from_plus.isValid()) {
00460         merged = from_plus;
00461       }
00462       else if (found_minus  &&  from_minus.isValid()) {
00463         merged = from_minus;
00464       }
00465       else return false;
00466 
00467       // if, after all that, we have a good fit-and-propagation, report the direction
00468       if (merged.isValid()) {
00469         double angle = merged.globalPosition().phi() + M_PI/2.;
00470         GlobalVector direction = merged.globalDirection();
00471         double dxdz = direction.x() / direction.z();
00472         double dydz = direction.y() / direction.z();
00473         dphidz = (dxdz*cos(angle) + dydz*sin(angle)) / merged.globalPosition().perp();
00474         return true;
00475       }
00476 
00477     } // end if refit successful
00478   } // end if enough hits
00479   return false;
00480 }
00481 
00482 void CSCPairResidualsConstraint::write(std::ofstream &output) {
00483   output << std::setprecision(14) << std::fixed;
00484   output << "CSCPairResidualsConstraint " << m_identifier << " " << i() << " " << j() << " "
00485          << m_sumN << " " << m_sum1 << " " << m_sumx << " " << m_sumy << " " << m_sumxx << " " << m_sumyy << " " << m_sumxy << " EOLN" << std::endl;
00486 }
00487 
00488 void CSCPairResidualsConstraint::read(std::vector<std::ifstream*> &input, std::vector<std::string> &filenames) {
00489   m_sumN = 0;
00490   m_sum1 = 0.;
00491   m_sumx = 0.;
00492   m_sumy = 0.;
00493   m_sumxx = 0.;
00494   m_sumyy = 0.;
00495   m_sumxy = 0.;
00496 
00497   std::vector<std::ifstream*>::const_iterator inputiter = input.begin();
00498   std::vector<std::string>::const_iterator filename = filenames.begin();
00499   for (;  inputiter != input.end();  ++inputiter, ++filename) {
00500     int linenumber = 0;
00501     bool touched = false;
00502     while (!(*inputiter)->eof()) {
00503       linenumber++;
00504       std::string name, eoln;
00505       unsigned int identifier;
00506       int i, j;
00507       int sumN;
00508       double sum1, sumx, sumy, sumxx, sumyy, sumxy;
00509     
00510       (**inputiter) >> name >> identifier >> i >> j >> sumN >> sum1 >> sumx >> sumy >> sumxx >> sumyy >> sumxy >> eoln;
00511 
00512       if (!(*inputiter)->eof()  &&  (name != "CSCPairResidualsConstraint"  ||  eoln != "EOLN")) throw cms::Exception("CorruptTempFile") << "Temporary file " << *filename << " is incorrectly formatted on line " << linenumber << std::endl;
00513 
00514       if (identifier == m_identifier) {
00515         if (i != m_i  ||  j != m_j) throw cms::Exception("CorruptTempFile") << "Wrong (i,j) for CSCPairResidualsConstraint " << m_identifier << " (" << m_i << "," << m_j << ") in file " << *filename << " on line " << linenumber << std::endl;
00516         touched = true;
00517 
00518         m_sumN += sumN;
00519         m_sum1 += sum1;
00520         m_sumx += sumx;
00521         m_sumy += sumy;
00522         m_sumxx += sumxx;
00523         m_sumyy += sumyy;
00524         m_sumxy += sumxy;
00525       }
00526     }
00527 
00528     (*inputiter)->clear();
00529     (*inputiter)->seekg(0, std::ios::beg);
00530 
00531     if (!touched) throw cms::Exception("CorruptTempFile") << "CSCPairResidualsConstraint " << m_identifier << " is missing from file " << *filename << std::endl;
00532   }
00533 }
00534 
00535 void CSCPairResidualsConstraint::calculatePhi(const TransientTrackingRecHit *hit, double &phi, double &phierr2, bool doRphi, bool globalPhi) {
00536   align::LocalPoint pos = hit->localPosition();
00537   DetId id = hit->geographicalId();
00538   CSCDetId cscid = CSCDetId(id.rawId());
00539 
00540   double r = 0.;
00541   if (globalPhi) {
00542     phi = hit->globalPosition().phi();
00543     r = hit->globalPosition().perp();
00544 
00545 //     double sinAngle = sin(phi);
00546 //     double cosAngle = cos(phi);
00547 //     double xx = hit->globalPositionError().cxx();
00548 //     double xy = hit->globalPositionError().cyx();
00549 //     double yy = hit->globalPositionError().cyy();
00550 //     phierr2 = (xx*cosAngle*cosAngle + 2.*xy*sinAngle*cosAngle + yy*sinAngle*sinAngle) / (r*r);
00551   }
00552   else {
00553     // these constants are related to the way CSC chambers are built--- really constant!
00554     const double R_ME11 = 181.5;
00555     const double R_ME12 = 369.7;
00556     const double R_ME21 = 242.7;
00557     const double R_ME31 = 252.7;
00558     const double R_ME41 = 262.65;
00559     const double R_MEx2 = 526.5;
00560 
00561     double R = 0.;
00562     if (cscid.station() == 1  &&  (cscid.ring() == 1  ||  cscid.ring() == 4)) R = R_ME11;
00563     else if (cscid.station() == 1  &&  cscid.ring() == 2) R = R_ME12;
00564     else if (cscid.station() == 2  &&  cscid.ring() == 1) R = R_ME21;
00565     else if (cscid.station() == 3  &&  cscid.ring() == 1) R = R_ME31;
00566     else if (cscid.station() == 4  &&  cscid.ring() == 1) R = R_ME41;
00567     else if (cscid.station() > 1  &&  cscid.ring() == 2) R = R_MEx2;
00568     else assert(false);
00569     r = (pos.y() + R);
00570   
00571     phi = atan2(pos.x(), r);
00572 
00573     if (cscid.endcap() == 1  &&  cscid.station() >= 3) phi *= -1;
00574     else if (cscid.endcap() == 2  &&  cscid.station() <= 2) phi *= -1;
00575   }
00576 
00577   int strip = m_cscGeometry->layer(id)->geometry()->nearestStrip(pos);
00578   double angle = m_cscGeometry->layer(id)->geometry()->stripAngle(strip) - M_PI/2.;
00579   double sinAngle = sin(angle);
00580   double cosAngle = cos(angle);
00581   double xx = hit->localPositionError().xx();
00582   double xy = hit->localPositionError().xy();
00583   double yy = hit->localPositionError().yy();
00584   phierr2 = (xx*cosAngle*cosAngle + 2.*xy*sinAngle*cosAngle + yy*sinAngle*sinAngle) / (r*r);
00585 
00586   if (doRphi) {
00587     phi *= r;
00588     phierr2 *= r*r;
00589   }
00590 }
00591 
00592 bool CSCPairResidualsConstraint::isFiducial(std::vector<const TransientTrackingRecHit*> &hits, bool is_i) {
00593   // these constants are related to the way CSC chambers are built--- really constant!
00594   const double cut_ME11 = 0.086;
00595   const double cut_ME12 = 0.090;
00596   const double cut_MEx1 = 0.180;
00597   const double cut_MEx2 = 0.090;
00598 
00599   double sum1 = 0.;
00600   double sumx = 0.;
00601   double sumy = 0.;
00602   double sumxx = 0.;
00603   double sumxy = 0.;
00604   for (std::vector<const TransientTrackingRecHit*>::const_iterator hit = hits.begin();  hit != hits.end();  ++hit) {
00605     double phi, phierr2;
00606     calculatePhi(*hit, phi, phierr2);
00607     double z = (is_i ? m_cscGeometry->idToDet(m_id_i)->surface() : m_cscGeometry->idToDet(m_id_j)->surface()).toLocal((*hit)->globalPosition()).z();
00608 
00609     if (m_parent->m_makeHistograms) {
00610       if (m_id_i.station() == 1  &&  (m_id_i.ring() == 1  ||  m_id_i.ring() == 4)) {
00611         m_parent->m_fiducial_ME11->Fill(fabs(phi), sqrt(phierr2));
00612       }
00613       else if (m_id_i.station() == 1  &&  m_id_i.ring() == 2) {
00614         m_parent->m_fiducial_ME12->Fill(fabs(phi), sqrt(phierr2));
00615       }
00616       else if (m_id_i.station() > 1  &&  m_id_i.ring() == 1) {
00617         m_parent->m_fiducial_MEx1->Fill(fabs(phi), sqrt(phierr2));
00618       }
00619       else if (m_id_i.station() > 1  &&  m_id_i.ring() == 2) {
00620         m_parent->m_fiducial_MEx2->Fill(fabs(phi), sqrt(phierr2));
00621       }
00622     }
00623     
00624     double weight = 1.;
00625     if (m_parent->m_useHitWeights) weight = 1./phierr2;
00626     sum1 += weight;
00627     sumx += weight * z;
00628     sumy += weight * phi;
00629     sumxx += weight * z*z;
00630     sumxy += weight * z*phi;
00631   }
00632   double delta = (sum1*sumxx) - (sumx*sumx);
00633   if (delta == 0.) return false;
00634   double intercept = ((sumxx*sumy) - (sumx*sumxy))/delta;
00635   double slope = ((sum1*sumxy) - (sumx*sumy))/delta;
00636 
00637   double phi1 = intercept + slope*(is_i ? m_iZ1 : m_jZ1);
00638   double phi6 = intercept + slope*(is_i ? m_iZ6 : m_jZ6);
00639 
00640   if (m_id_i.station() == 1  &&  (m_id_i.ring() == 1  ||  m_id_i.ring() == 4)) {
00641     return (fabs(phi1) < cut_ME11  &&  fabs(phi6) < cut_ME11);
00642   }
00643   else if (m_id_i.station() == 1  &&  m_id_i.ring() == 2) {
00644     return (fabs(phi1) < cut_ME12  &&  fabs(phi6) < cut_ME12);
00645   }
00646   else if (m_id_i.station() > 1  &&  m_id_i.ring() == 1) {
00647     return (fabs(phi1) < cut_MEx1  &&  fabs(phi6) < cut_MEx1);
00648   }
00649   else if (m_id_i.station() > 1  &&  m_id_i.ring() == 2) {
00650     return (fabs(phi1) < cut_MEx2  &&  fabs(phi6) < cut_MEx2);
00651   }
00652   else assert(false);
00653 }
/data/refman/pasoursint/CMSSW_5_3_9_patch3/src/Alignment/MuonAlignmentAlgorithms/plugins/CSCPairResidualsConstraint.cc
