Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:27:08 

0001 #include "RecoMuon/MuonSeedGenerator/src/MuonSeedPtExtractor.h"
0002 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0003 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0004 #include "DataFormats/MuonDetId/interface/CSCDetId.h"
0005 #include "DataFormats/MuonDetId/interface/ME0DetId.h"
0006 
0007 #include "TMath.h"
0008 #include <sstream>
0009 
0010 MuonSeedPtExtractor::MuonSeedPtExtractor(const edm::ParameterSet& par)
0011     : theBeamSpot(0., 0., 0.), scaleDT_(par.getParameter<bool>("scaleDT")) {
0012   init(par);
0013 }
0014 
0015 void MuonSeedPtExtractor::init(const edm::ParameterSet& par) {
0016   // load pT seed parameters
0017   // DT combinations
0018   fillParametersForCombo("DT_12", par);
0019   fillParametersForCombo("DT_13", par);
0020   fillParametersForCombo("DT_14", par);
0021   fillParametersForCombo("DT_23", par);
0022   fillParametersForCombo("DT_24", par);
0023   fillParametersForCombo("DT_34", par);
0024   // CSC combinations
0025   fillParametersForCombo("CSC_01", par);
0026   fillParametersForCombo("CSC_12", par);
0027   fillParametersForCombo("CSC_02", par);
0028   fillParametersForCombo("CSC_13", par);
0029   fillParametersForCombo("CSC_03", par);
0030   fillParametersForCombo("CSC_14", par);
0031   fillParametersForCombo("CSC_23", par);
0032   fillParametersForCombo("CSC_24", par);
0033   fillParametersForCombo("CSC_34", par);
0034 
0035   // Overlap combinations
0036   fillParametersForCombo("OL_1213", par);
0037   fillParametersForCombo("OL_1222", par);
0038   fillParametersForCombo("OL_1232", par);
0039   fillParametersForCombo("OL_2213", par);
0040   fillParametersForCombo("OL_2222", par);
0041 
0042   // Single segments (CSC)
0043   fillParametersForCombo("SME_11", par);
0044   fillParametersForCombo("SME_12", par);
0045   fillParametersForCombo("SME_13", par);
0046   fillParametersForCombo("SME_21", par);
0047   fillParametersForCombo("SME_22", par);
0048   fillParametersForCombo("SME_31", par);
0049   fillParametersForCombo("SME_32", par);
0050   fillParametersForCombo("SME_41", par);
0051   fillParametersForCombo("SME_42", par);
0052 
0053   // Single segments (DT)
0054   fillParametersForCombo("SMB_10", par);
0055   fillParametersForCombo("SMB_11", par);
0056   fillParametersForCombo("SMB_12", par);
0057   fillParametersForCombo("SMB_20", par);
0058   fillParametersForCombo("SMB_21", par);
0059   fillParametersForCombo("SMB_22", par);
0060   fillParametersForCombo("SMB_30", par);
0061   fillParametersForCombo("SMB_31", par);
0062   fillParametersForCombo("SMB_32", par);
0063 
0064   fillScalesForCombo("CSC_01_1_scale", par);
0065   fillScalesForCombo("CSC_12_1_scale", par);
0066   fillScalesForCombo("CSC_12_2_scale", par);
0067   fillScalesForCombo("CSC_12_3_scale", par);
0068   fillScalesForCombo("CSC_13_2_scale", par);
0069   fillScalesForCombo("CSC_13_3_scale", par);
0070   fillScalesForCombo("CSC_14_3_scale", par);
0071   fillScalesForCombo("CSC_23_1_scale", par);
0072   fillScalesForCombo("CSC_23_2_scale", par);
0073   fillScalesForCombo("CSC_24_1_scale", par);
0074   fillScalesForCombo("CSC_34_1_scale", par);
0075   fillScalesForCombo("DT_12_1_scale", par);
0076   fillScalesForCombo("DT_12_2_scale", par);
0077   fillScalesForCombo("DT_13_1_scale", par);
0078   fillScalesForCombo("DT_13_2_scale", par);
0079   fillScalesForCombo("DT_14_1_scale", par);
0080   fillScalesForCombo("DT_14_2_scale", par);
0081   fillScalesForCombo("DT_23_1_scale", par);
0082   fillScalesForCombo("DT_23_2_scale", par);
0083   fillScalesForCombo("DT_24_1_scale", par);
0084   fillScalesForCombo("DT_24_2_scale", par);
0085   fillScalesForCombo("DT_34_1_scale", par);
0086   fillScalesForCombo("DT_34_2_scale", par);
0087   fillScalesForCombo("OL_1213_0_scale", par);
0088   fillScalesForCombo("OL_1222_0_scale", par);
0089   fillScalesForCombo("OL_1232_0_scale", par);
0090   fillScalesForCombo("OL_2213_0_scale", par);
0091   fillScalesForCombo("OL_2222_0_scale", par);
0092   fillScalesForCombo("SMB_10_0_scale", par);
0093   fillScalesForCombo("SMB_11_0_scale", par);
0094   fillScalesForCombo("SMB_12_0_scale", par);
0095   fillScalesForCombo("SMB_20_0_scale", par);
0096   fillScalesForCombo("SMB_21_0_scale", par);
0097   fillScalesForCombo("SMB_22_0_scale", par);
0098   fillScalesForCombo("SMB_30_0_scale", par);
0099   fillScalesForCombo("SMB_31_0_scale", par);
0100   fillScalesForCombo("SMB_32_0_scale", par);
0101   fillScalesForCombo("SME_11_0_scale", par);
0102   fillScalesForCombo("SME_12_0_scale", par);
0103   fillScalesForCombo("SME_13_0_scale", par);
0104   fillScalesForCombo("SME_21_0_scale", par);
0105   fillScalesForCombo("SME_22_0_scale", par);
0106 }
0107 
0108 MuonSeedPtExtractor::~MuonSeedPtExtractor() {}
0109 
0110 void MuonSeedPtExtractor::fillParametersForCombo(const std::string& name, const edm::ParameterSet& pset) {
0111   theParametersForCombo[name] = pset.getParameter<std::vector<double> >(name);
0112 }
0113 
0114 void MuonSeedPtExtractor::fillScalesForCombo(const std::string& name, const edm::ParameterSet& pset) {
0115   theScalesForCombo[name] = pset.getParameter<std::vector<double> >(name);
0116 }
0117 
0118 std::vector<double> MuonSeedPtExtractor::pT_extract(
0119     MuonTransientTrackingRecHit::ConstMuonRecHitPointer firstHit,
0120     MuonTransientTrackingRecHit::ConstMuonRecHitPointer secondHit) const {
0121   GlobalPoint innerPoint = firstHit->globalPosition() - theBeamSpot;
0122   GlobalPoint outerPoint = secondHit->globalPosition() - theBeamSpot;
0123   MuonTransientTrackingRecHit::ConstMuonRecHitPointer innerHit = firstHit;
0124   MuonTransientTrackingRecHit::ConstMuonRecHitPointer outerHit = secondHit;
0125 
0126   // ways in which the hits could be in the wrong order
0127   if ((outerHit->isDT() && innerHit->isCSC()) ||
0128       (outerHit->isDT() && innerHit->isDT() && (innerPoint.perp() > outerPoint.perp())) ||
0129       (outerHit->isCSC() && innerHit->isCSC() && (fabs(innerPoint.z()) > fabs(outerPoint.z())))) {
0130     innerHit = secondHit;
0131     outerHit = firstHit;
0132     innerPoint = innerHit->globalPosition() - theBeamSpot;
0133     outerPoint = outerHit->globalPosition() - theBeamSpot;
0134   }
0135 
0136   double phiInner = innerPoint.phi();
0137   double phiOuter = outerPoint.phi();
0138 
0139   double etaInner = innerPoint.eta();
0140   double etaOuter = outerPoint.eta();
0141   //std::cout<<" inner pos = "<< innerPoint << " phi eta " << phiInner << " " << etaInner << std::endl;
0142   //std::cout<<" outer pos = "<< outerPoint << " phi eta " << phiOuter << " " << etaOuter << std::endl;
0143   //double thetaInner = firstHit->globalPosition().theta();
0144   // if some of the segments is missing r-phi measurement then we should
0145   // use only the 4D phi estimate (also use 4D eta estimate only)
0146   // the direction is not so important (it will be corrected)
0147   /*
0148     bool firstOK = (4==allValidSets[iSet][firstMeasurement]->hit->dimension());
0149     bool lastOK = (4==allValidSets[iSet][lastMeasurement]->hit->dimension());
0150     if(!(firstOK * lastOK)){
0151     if(!firstOK){
0152     }
0153     if(!firstOK){
0154     }
0155     }
0156   */
0157   double dPhi = phiInner - phiOuter;
0158   if (dPhi < -TMath::Pi()) {
0159     dPhi += 2 * TMath::Pi();
0160   } else if (dPhi > TMath::Pi()) {
0161     dPhi -= 2 * TMath::Pi();
0162   }
0163   int sign = 1;
0164   if (dPhi < 0.) {
0165     dPhi = -dPhi;
0166     sign = -1;
0167   }
0168 
0169   if (dPhi < 1.0e-6) {
0170     dPhi = 1.0e-6;
0171   }
0172   double eta = fabs(etaOuter);  // what if it is 2D segment? use inner?
0173 
0174   std::vector<int> stationCoded(2);
0175   stationCoded[0] = stationCoded[1] = 999;
0176 
0177   DetId detId_inner = innerHit->hit()->geographicalId();
0178   DetId detId_outer = outerHit->hit()->geographicalId();
0179 
0180   stationCoded[0] = stationCode(innerHit);
0181   stationCoded[1] = stationCode(outerHit);
0182 
0183   std::ostringstream os0;
0184   std::ostringstream os1;
0185   os0 << abs(stationCoded[0]);
0186   os1 << abs(stationCoded[1]);
0187 
0188   //std::cout<<" st1 = "<<stationCoded[0]<<" st2 = "<<stationCoded[1]<<std::endl;
0189   //std::cout<<" detId_inner = "<<detId_inner.rawId()<<"  detId_outer = "<< detId_outer.rawId()<<std::endl;
0190   std::string combination = "0";
0191   std::string init_combination = combination;
0192   bool singleSegment = false;
0193   //if(detId_first == detId_second){
0194   if (stationCoded[0] == stationCoded[1]) {
0195     // single segment - DT or CSC
0196     singleSegment = true;
0197     //eta = innerPoint.eta();
0198     GlobalVector gv = innerHit->globalDirection();
0199     double gvPerpPos = gv.x() * gv.x() + gv.y() * gv.y();
0200     if (gvPerpPos < 1e-32)
0201       gvPerpPos = 1e-32;
0202     gvPerpPos = sqrt(gvPerpPos);
0203     // Psi is angle between the segment origin and segment direction
0204     // Use dot product between two vectors to get Psi in global x-y plane
0205     double cosDpsi = (gv.x() * innerPoint.x() + gv.y() * innerPoint.y());
0206     if (cosDpsi != 0) {
0207       cosDpsi /= sqrt(innerPoint.x() * innerPoint.x() + innerPoint.y() * innerPoint.y());
0208       cosDpsi /= gvPerpPos;
0209       cosDpsi = cosDpsi > 1 ? 1 : cosDpsi;
0210       cosDpsi = cosDpsi < -1 ? -1 : cosDpsi;
0211     }
0212 
0213     double axb = (innerPoint.x() * gv.y()) - (innerPoint.y() * gv.x());
0214     sign = (axb < 0.) ? 1 : -1;
0215 
0216     double dpsi = fabs(acos(cosDpsi));
0217     if (dpsi > TMath::Pi() / 2.) {
0218       dpsi = TMath::Pi() - dpsi;
0219     }
0220     if (fabs(dpsi) < 0.00005) {
0221       dpsi = 0.00005;
0222     }
0223     dPhi = dpsi;
0224 
0225     if (innerHit->isDT()) {
0226       DTChamberId dtCh(detId_inner);
0227       std::ostringstream os;
0228       os << dtCh.station() << abs(dtCh.wheel());
0229       combination = "SMB_" + os.str();
0230     } else if (innerHit->isCSC()) {
0231       CSCDetId cscId(detId_inner);
0232       std::ostringstream os;
0233       int ring = cscId.ring();
0234       if (ring == 4)
0235         ring = 1;
0236       os << cscId.station() << ring;
0237       combination = "SME_" + os.str();
0238     } else if (innerHit->isME0()) {
0239       ME0DetId me0Id(detId_inner);
0240       std::ostringstream os;
0241       int ring = 1;  //me0Id.ring(); me0 only in ring 1
0242       os << me0Id.station() << ring;
0243       combination = "SME_" + os.str();
0244     } else {
0245       throw cms::Exception("MuonSeedPtExtractor") << "Bad hit DetId";
0246     }
0247   } else {
0248     if (stationCoded[0] < 0) {
0249       if (stationCoded[1] < 0) {
0250         // DT-DT
0251         combination = "DT_" + os0.str() + os1.str();
0252       } else {
0253         // DT-CSC
0254         eta = fabs(etaInner);
0255         if (-1 == stationCoded[0]) {
0256           switch (stationCoded[1]) {
0257             case 1:
0258               combination = "OL_1213";
0259               break;
0260             case 2:
0261               combination = "OL_1222";
0262               break;
0263             case 3:
0264               combination = "OL_1232";
0265               break;
0266             default:
0267               // can not be
0268               break;
0269           }
0270         } else if (-2 == stationCoded[0]) {
0271           if (1 == stationCoded[1]) {
0272             combination = "OL_2213";
0273           } else {
0274             // can not be (not coded?)
0275             combination = "OL_2222";  // in case
0276           }
0277         } else {
0278           // can not be
0279         }
0280       }
0281     } else {
0282       if (stationCoded[1] < 0) {
0283         // CSC-DT
0284         // can not be
0285       } else {
0286         // CSC-CSC
0287         combination = "CSC_" + os0.str() + os1.str();
0288         if ("CSC_04" == combination) {
0289           combination = "CSC_14";
0290         }
0291       }
0292     }
0293   }
0294 
0295   std::vector<double> pTestimate(2);  //
0296   //std::cout<<" combination = "<<combination<<std::endl;
0297   if (init_combination != combination) {
0298     //std::cout<<" combination = "<<combination<<" eta = "<<eta<<" dPhi = "<<dPhi<<std::endl;
0299     ParametersMap::const_iterator parametersItr = theParametersForCombo.find(combination);
0300     if (parametersItr == theParametersForCombo.end()) {
0301       //      edm::LogWarning("RecoMuon|MuonSeedGenerator|MuonSeedPtExtractor") << "Cannot find parameters for combo " << combination;
0302       edm::LogWarning("BadSegmentCombination") << "Cannot find parameters for combo " << combination;
0303       pTestimate[0] = pTestimate[1] = 100;
0304       //       throw cms::Exception("MuonSeedPtEstimator") << "Cannot find parameters for combo " << combination;
0305     } else {
0306       if (scaleDT_ && outerHit->isDT()) {
0307         pTestimate = getPt(parametersItr->second, eta, dPhi, combination, detId_outer);
0308       } else {
0309         pTestimate = getPt(parametersItr->second, eta, dPhi);
0310       }
0311 
0312       if (singleSegment) {
0313         pTestimate[0] = fabs(pTestimate[0]);
0314         pTestimate[1] = fabs(pTestimate[1]);
0315       }
0316       pTestimate[0] *= double(sign);
0317     }
0318   } else {
0319     // often a MB3 - ME1/3 seed
0320     pTestimate[0] = pTestimate[1] = 100;
0321     // hmm
0322   }
0323   // std::cout<<" pTestimate[0] = "<<pTestimate[0]<<" pTestimate[1] = "<<pTestimate[1]<<std::endl;
0324   /*
0325                               MuonRecHitContainer_clusters[cluster][iHit+1]->isDT());
0326           if(specialCase){
0327             DTChamberId dtCh(detId);
0328             DTChamberId dtCh_2(detId_2);
0329             specialCase =  (dtCh.station() == dtCh_2.station());
0330           }
0331   */
0332   //return vPara;
0333   return pTestimate;
0334 }
0335 
0336 int MuonSeedPtExtractor::stationCode(MuonTransientTrackingRecHit::ConstMuonRecHitPointer hit) const {
0337   DetId detId(hit->hit()->geographicalId());
0338   int result = -999;
0339   if (hit->isDT()) {
0340     DTChamberId dtCh(detId);
0341     //std::cout<<"first (DT) St/W/S = "<<dtCh.station()<<"/"<<dtCh.wheel()<<"/"<<dtCh.sector()<<"/"<<std::endl;
0342     result = -1 * dtCh.station();
0343   } else if (hit->isCSC()) {
0344     CSCDetId cscID(detId);
0345     //std::cout<<"first (CSC) E/S/R/C = "<<cscID.endcap()<<"/"<<cscID.station()<<"/"<<cscID.ring()<<"/"<<cscID.chamber()<<std::endl;
0346     result = cscID.station();
0347     if (result == 1 && (1 == cscID.ring() || 4 == cscID.ring()))
0348       result = 0;
0349   } else if (hit->isME0()) {
0350     result = 0;
0351   } else if (hit->isRPC()) {
0352   }
0353   return result;
0354 }
0355 
0356 std::vector<double> MuonSeedPtExtractor::getPt(const std::vector<double>& vPara, double eta, double dPhi) const {
0357   //std::cout<<" eta = "<<eta<<" dPhi = "<<dPhi<<" vPara[0] = "<<vPara[0]<<" vPara[1] = "<<vPara[1]<<" vPara[2] = "<<vPara[2]<<std::endl;
0358   double h = fabs(eta);
0359   double estPt = (vPara[0] + vPara[1] * h + vPara[2] * h * h) / dPhi;
0360   double estSPt = (vPara[3] + vPara[4] * h + vPara[5] * h * h) / dPhi;
0361   // std::cout<<"estPt = "<<estPt<<std::endl;
0362   std::vector<double> paraPt;
0363   paraPt.push_back(estPt);
0364   paraPt.push_back(estSPt);
0365   return paraPt;
0366 }
0367 
0368 // combined pt parameterization and pt scale
0369 std::vector<double> MuonSeedPtExtractor::getPt(const std::vector<double>& vPara,
0370                                                double eta,
0371                                                double dPhi,
0372                                                const std::string& combination,
0373                                                const DTChamberId& outerDetId) const {
0374   double h = fabs(eta);
0375   double estPt = (vPara[0] + vPara[1] * h + vPara[2] * h * h) / dPhi;
0376   // changed by S.C.
0377   double estSPt = (vPara[3] + vPara[4] * h + vPara[5] * h * h) / dPhi;
0378 
0379   // scale the pt and spt , changed by S.C.
0380   int wheel = 0;
0381   if (combination[0] == 'D') {
0382     wheel = abs(outerDetId.wheel());
0383   }
0384 
0385   std::ostringstream os;
0386   os << combination << "_" << wheel << "_scale";
0387 
0388   ScalesMap::const_iterator scalesItr = theScalesForCombo.find(os.str());
0389   if (scalesItr != theScalesForCombo.end()) {
0390     double t1 = scalesItr->second[0];
0391     double scaleFactor = 1. / (1. + t1 / (estPt + 10.));
0392 
0393     estSPt = estSPt * scaleFactor;
0394     estPt = estPt * scaleFactor;
0395   }
0396 
0397   // std::cout<<"estPt = "<<estPt<<std::endl;
0398   std::vector<double> paraPt;
0399   paraPt.push_back(estPt);
0400   paraPt.push_back(estSPt);
0401   return paraPt;
0402 }
[ 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!