Back to home page

Project CMSSW displayed by LXR
 
 

    

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

0001 /**
0002  *  See header file for a description of this class.
0003  *
0004  *  \author: Shih-Chuan Kao, Dominique Fortin - UCR
0005  *
0006  */
0007 
0008 #include "RecoMuon/MuonSeedGenerator/src/MuonSeedCreator.h"
0009 
0010 #include "RecoMuon/TransientTrackingRecHit/interface/MuonTransientTrackingRecHit.h"
0011 
0012 #include "RecoMuon/TrackingTools/interface/MuonPatternRecoDumper.h"
0013 
0014 #include "MagneticField/Engine/interface/MagneticField.h"
0015 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
0016 
0017 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
0018 #include "TrackingTools/DetLayers/interface/DetLayer.h"
0019 
0020 #include "DataFormats/TrajectoryState/interface/PTrajectoryStateOnDet.h"
0021 #include "DataFormats/Common/interface/OwnVector.h"
0022 #include "DataFormats/MuonDetId/interface/DTChamberId.h"
0023 #include "DataFormats/MuonDetId/interface/CSCDetId.h"
0024 #include "DataFormats/MuonDetId/interface/RPCDetId.h"
0025 
0026 #include "FWCore/Framework/interface/EventSetup.h"
0027 #include "FWCore/Framework/interface/ESHandle.h"
0028 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0029 
0030 #include "gsl/gsl_statistics.h"
0031 
0032 /*
0033  * Constructor
0034  */
0035 MuonSeedCreator::MuonSeedCreator(const edm::ParameterSet& pset) {
0036   theMinMomentum = pset.getParameter<double>("minimumSeedPt");
0037   theMaxMomentum = pset.getParameter<double>("maximumSeedPt");
0038   defaultMomentum = pset.getParameter<double>("defaultSeedPt");
0039   debug = pset.getParameter<bool>("DebugMuonSeed");
0040   sysError = pset.getParameter<double>("SeedPtSystematics");
0041   // load seed PT parameters
0042   DT12 = pset.getParameter<std::vector<double> >("DT_12");
0043   DT13 = pset.getParameter<std::vector<double> >("DT_13");
0044   DT14 = pset.getParameter<std::vector<double> >("DT_14");
0045   DT23 = pset.getParameter<std::vector<double> >("DT_23");
0046   DT24 = pset.getParameter<std::vector<double> >("DT_24");
0047   DT34 = pset.getParameter<std::vector<double> >("DT_34");
0048 
0049   CSC01 = pset.getParameter<std::vector<double> >("CSC_01");
0050   CSC12 = pset.getParameter<std::vector<double> >("CSC_12");
0051   CSC02 = pset.getParameter<std::vector<double> >("CSC_02");
0052   CSC13 = pset.getParameter<std::vector<double> >("CSC_13");
0053   CSC03 = pset.getParameter<std::vector<double> >("CSC_03");
0054   CSC14 = pset.getParameter<std::vector<double> >("CSC_14");
0055   CSC23 = pset.getParameter<std::vector<double> >("CSC_23");
0056   CSC24 = pset.getParameter<std::vector<double> >("CSC_24");
0057   CSC34 = pset.getParameter<std::vector<double> >("CSC_34");
0058 
0059   OL1213 = pset.getParameter<std::vector<double> >("OL_1213");
0060   OL1222 = pset.getParameter<std::vector<double> >("OL_1222");
0061   OL1232 = pset.getParameter<std::vector<double> >("OL_1232");
0062   OL1213 = pset.getParameter<std::vector<double> >("OL_1213");
0063   OL2222 = pset.getParameter<std::vector<double> >("OL_1222");
0064 
0065   SME11 = pset.getParameter<std::vector<double> >("SME_11");
0066   SME12 = pset.getParameter<std::vector<double> >("SME_12");
0067   SME13 = pset.getParameter<std::vector<double> >("SME_13");
0068   SME21 = pset.getParameter<std::vector<double> >("SME_21");
0069   SME22 = pset.getParameter<std::vector<double> >("SME_22");
0070   SME31 = pset.getParameter<std::vector<double> >("SME_31");
0071   SME32 = pset.getParameter<std::vector<double> >("SME_32");
0072   SME41 = pset.getParameter<std::vector<double> >("SME_41");
0073 
0074   SMB10 = pset.getParameter<std::vector<double> >("SMB_10");
0075   SMB11 = pset.getParameter<std::vector<double> >("SMB_11");
0076   SMB12 = pset.getParameter<std::vector<double> >("SMB_12");
0077   SMB20 = pset.getParameter<std::vector<double> >("SMB_20");
0078   SMB21 = pset.getParameter<std::vector<double> >("SMB_21");
0079   SMB22 = pset.getParameter<std::vector<double> >("SMB_22");
0080   SMB30 = pset.getParameter<std::vector<double> >("SMB_30");
0081   SMB31 = pset.getParameter<std::vector<double> >("SMB_31");
0082   SMB32 = pset.getParameter<std::vector<double> >("SMB_32");
0083 
0084   // Load dphi scale parameters
0085   CSC01_1 = pset.getParameter<std::vector<double> >("CSC_01_1_scale");
0086   CSC12_1 = pset.getParameter<std::vector<double> >("CSC_12_1_scale");
0087   CSC12_2 = pset.getParameter<std::vector<double> >("CSC_12_2_scale");
0088   CSC12_3 = pset.getParameter<std::vector<double> >("CSC_12_3_scale");
0089   CSC13_2 = pset.getParameter<std::vector<double> >("CSC_13_2_scale");
0090   CSC13_3 = pset.getParameter<std::vector<double> >("CSC_13_3_scale");
0091   CSC14_3 = pset.getParameter<std::vector<double> >("CSC_14_3_scale");
0092   CSC23_1 = pset.getParameter<std::vector<double> >("CSC_23_1_scale");
0093   CSC23_2 = pset.getParameter<std::vector<double> >("CSC_23_2_scale");
0094   CSC24_1 = pset.getParameter<std::vector<double> >("CSC_24_1_scale");
0095   CSC34_1 = pset.getParameter<std::vector<double> >("CSC_34_1_scale");
0096 
0097   DT12_1 = pset.getParameter<std::vector<double> >("DT_12_1_scale");
0098   DT12_2 = pset.getParameter<std::vector<double> >("DT_12_2_scale");
0099   DT13_1 = pset.getParameter<std::vector<double> >("DT_13_1_scale");
0100   DT13_2 = pset.getParameter<std::vector<double> >("DT_13_2_scale");
0101   DT14_1 = pset.getParameter<std::vector<double> >("DT_14_1_scale");
0102   DT14_2 = pset.getParameter<std::vector<double> >("DT_14_2_scale");
0103   DT23_1 = pset.getParameter<std::vector<double> >("DT_23_1_scale");
0104   DT23_2 = pset.getParameter<std::vector<double> >("DT_23_2_scale");
0105   DT24_1 = pset.getParameter<std::vector<double> >("DT_24_1_scale");
0106   DT24_2 = pset.getParameter<std::vector<double> >("DT_24_2_scale");
0107   DT34_1 = pset.getParameter<std::vector<double> >("DT_34_1_scale");
0108   DT34_2 = pset.getParameter<std::vector<double> >("DT_34_2_scale");
0109 
0110   OL_1213 = pset.getParameter<std::vector<double> >("OL_1213_0_scale");
0111   OL_1222 = pset.getParameter<std::vector<double> >("OL_1222_0_scale");
0112   OL_1232 = pset.getParameter<std::vector<double> >("OL_1232_0_scale");
0113   OL_2213 = pset.getParameter<std::vector<double> >("OL_2213_0_scale");
0114   OL_2222 = pset.getParameter<std::vector<double> >("OL_2222_0_scale");
0115 
0116   SMB_10S = pset.getParameter<std::vector<double> >("SMB_10_0_scale");
0117   SMB_11S = pset.getParameter<std::vector<double> >("SMB_11_0_scale");
0118   SMB_12S = pset.getParameter<std::vector<double> >("SMB_12_0_scale");
0119   SMB_20S = pset.getParameter<std::vector<double> >("SMB_20_0_scale");
0120   SMB_21S = pset.getParameter<std::vector<double> >("SMB_21_0_scale");
0121   SMB_22S = pset.getParameter<std::vector<double> >("SMB_22_0_scale");
0122   SMB_30S = pset.getParameter<std::vector<double> >("SMB_30_0_scale");
0123   SMB_31S = pset.getParameter<std::vector<double> >("SMB_31_0_scale");
0124   SMB_32S = pset.getParameter<std::vector<double> >("SMB_32_0_scale");
0125 
0126   SME_11S = pset.getParameter<std::vector<double> >("SME_11_0_scale");
0127   SME_12S = pset.getParameter<std::vector<double> >("SME_12_0_scale");
0128   SME_13S = pset.getParameter<std::vector<double> >("SME_13_0_scale");
0129   SME_21S = pset.getParameter<std::vector<double> >("SME_21_0_scale");
0130   SME_22S = pset.getParameter<std::vector<double> >("SME_22_0_scale");
0131 }
0132 
0133 /*
0134  * Destructor
0135  */
0136 MuonSeedCreator::~MuonSeedCreator() {}
0137 
0138 /*
0139  * createSeed
0140  *
0141  * Note type = 1 --> CSC
0142  *           = 2 --> Overlap
0143  *           = 3 --> DT
0144  */
0145 TrajectorySeed MuonSeedCreator::createSeed(
0146     int type, const SegmentContainer& seg, const std::vector<int>& layers, int NShowers, int NShowerSegments) {
0147   // The index of the station closest to the IP
0148   int last = 0;
0149 
0150   double ptmean = theMinMomentum;
0151   double sptmean = theMinMomentum;
0152 
0153   AlgebraicVector t(4);
0154   AlgebraicSymMatrix mat(5, 0);
0155   LocalPoint segPos;
0156 
0157   // Compute the pt according to station types used;
0158   if (type == 1)
0159     estimatePtCSC(seg, layers, ptmean, sptmean);
0160   if (type == 2)
0161     estimatePtOverlap(seg, layers, ptmean, sptmean);
0162   if (type == 3)
0163     estimatePtDT(seg, layers, ptmean, sptmean);
0164   if (type == 4)
0165     estimatePtSingle(seg, layers, ptmean, sptmean);
0166   // type 5 are the seeding for ME1/4
0167   if (type == 5)
0168     estimatePtCSC(seg, layers, ptmean, sptmean);
0169 
0170   // for certain clear showering case, set-up the minimum value
0171   if (NShowers > 0)
0172     estimatePtShowering(NShowers, NShowerSegments, ptmean, sptmean);
0173   //if ( NShowers > 0 ) std::cout<<" Showering happened "<<NShowers<<" times w/ "<< NShowerSegments<<std::endl; ;
0174 
0175   // Minimal pt
0176   double charge = 1.0;
0177   if (ptmean < 0.)
0178     charge = -1.0;
0179   if ((charge * ptmean) < theMinMomentum) {
0180     ptmean = theMinMomentum * charge;
0181     sptmean = theMinMomentum;
0182   } else if ((charge * ptmean) > theMaxMomentum) {
0183     ptmean = theMaxMomentum * charge;
0184     sptmean = theMaxMomentum * 0.25;
0185   }
0186 
0187   LocalTrajectoryParameters param;
0188 
0189   double p_err = 0.0;
0190   // determine the seed layer
0191   int best_seg = 0;
0192   double chi2_dof = 9999.0;
0193   unsigned int ini_seg = 0;
0194   // avoid generating seed from  1st layer(ME1/1)
0195   if (type == 5)
0196     ini_seg = 1;
0197   for (size_t i = ini_seg; i < seg.size(); i++) {
0198     double dof = static_cast<double>(seg[i]->degreesOfFreedom());
0199     if (chi2_dof < (seg[i]->chi2() / dof))
0200       continue;
0201     chi2_dof = seg[i]->chi2() / dof;
0202     best_seg = static_cast<int>(i);
0203   }
0204 
0205   if (type == 1 || type == 5 || type == 4) {
0206     // Fill the LocalTrajectoryParameters
0207     /// get the Global position
0208     last = best_seg;
0209     // last = 0;
0210     GlobalVector mom = seg[last]->globalPosition() - GlobalPoint();
0211     segPos = seg[last]->localPosition();
0212     /// get the Global direction
0213 
0214     GlobalVector polar(GlobalVector::Spherical(mom.theta(), seg[last]->globalDirection().phi(), 1.));
0215 
0216     /// scale the magnitude of total momentum
0217     polar *= fabs(ptmean) / polar.perp();
0218     /// Trasfer into local direction
0219     LocalVector segDirFromPos = seg[last]->det()->toLocal(polar);
0220     int chargeI = static_cast<int>(charge);
0221     LocalTrajectoryParameters param1(segPos, segDirFromPos, chargeI);
0222     param = param1;
0223     p_err = (sptmean * sptmean) / (polar.mag() * polar.mag() * ptmean * ptmean);
0224     mat = seg[last]->parametersError().similarityT(seg[last]->projectionMatrix());
0225     mat[0][0] = p_err;
0226     if (type == 5) {
0227       mat[0][0] = mat[0][0] / fabs(tan(mom.theta()));
0228       mat[1][1] = mat[1][1] / fabs(tan(mom.theta()));
0229       mat[3][3] = 2.25 * mat[3][3];
0230       mat[4][4] = 2.25 * mat[4][4];
0231     }
0232     if (type == 4) {
0233       mat[0][0] = mat[0][0] / fabs(tan(mom.theta()));
0234       mat[1][1] = mat[1][1] / fabs(tan(mom.theta()));
0235       mat[2][2] = 2.25 * mat[2][2];
0236       mat[3][3] = 2.25 * mat[3][3];
0237       mat[4][4] = 2.25 * mat[4][4];
0238     }
0239     double dh = fabs(seg[last]->globalPosition().eta()) - 1.6;
0240     if (fabs(dh) < 0.1 && type == 1) {
0241       mat[1][1] = 4. * mat[1][1];
0242       mat[2][2] = 4. * mat[2][2];
0243       mat[3][3] = 9. * mat[3][3];
0244       mat[4][4] = 9. * mat[4][4];
0245     }
0246 
0247     //if ( !highPt && type != 1 ) mat[1][1]= 2.25*mat[1][1];
0248     //if (  highPt && type != 1 ) mat[3][3]= 2.25*mat[1][1];
0249     //mat[2][2]= 3.*mat[2][2];
0250     //mat[3][3]= 2.*mat[3][3];
0251     //mat[4][4]= 2.*mat[4][4];
0252   } else {
0253     // Fill the LocalTrajectoryParameters
0254     /// get the Global position
0255     last = 0;
0256     segPos = seg[last]->localPosition();
0257     GlobalVector mom = seg[last]->globalPosition() - GlobalPoint();
0258     /// get the Global direction
0259     GlobalVector polar(GlobalVector::Spherical(mom.theta(), seg[last]->globalDirection().phi(), 1.));
0260     //GlobalVector polar(GlobalVector::Spherical(seg[last]->globalDirection().theta(),seg[last]->globalDirection().phi(),1.));
0261 
0262     /// count the energy loss - from parameterization
0263     //double ptRatio = 1.0 - (2.808/(fabs(ptmean) -1)) + (4.546/( (fabs(ptmean)-1)*(fabs(ptmean)-1)) );
0264     //ptmean = ptmean*ptRatio ;
0265 
0266     /// scale the magnitude of total momentum
0267     polar *= fabs(ptmean) / polar.perp();
0268     /// Trasfer into local direction
0269     LocalVector segDirFromPos = seg[last]->det()->toLocal(polar);
0270     int chargeI = static_cast<int>(charge);
0271     LocalTrajectoryParameters param1(segPos, segDirFromPos, chargeI);
0272     param = param1;
0273     p_err = (sptmean * sptmean) / (polar.mag() * polar.mag() * ptmean * ptmean);
0274     mat = seg[last]->parametersError().similarityT(seg[last]->projectionMatrix());
0275     //mat[0][0]= 1.44 * p_err;
0276     mat[0][0] = p_err;
0277   }
0278 
0279   if (debug) {
0280     GlobalPoint gp = seg[last]->globalPosition();
0281     float Geta = gp.eta();
0282 
0283     std::cout << "Type " << type << " Nsegments " << layers.size() << " ";
0284     std::cout << "pt " << ptmean << " errpt    " << sptmean << " eta " << Geta << " charge " << charge << std::endl;
0285   }
0286 
0287   // this perform H.T() * parErr * H, which is the projection of the
0288   // the measurement error (rechit rf) to the state error (TSOS rf)
0289   // Legend:
0290   // H => is the 4x5 projection matrix
0291   // parError the 4x4 parameter error matrix of the RecHit
0292 
0293   LocalTrajectoryError error(asSMatrix<5>(mat));
0294 
0295   // Create the TrajectoryStateOnSurface
0296   TrajectoryStateOnSurface tsos(param, error, seg[last]->det()->surface(), &*BField);
0297 
0298   // Take the DetLayer on which relies the segment
0299   DetId id = seg[last]->geographicalId();
0300 
0301   // Transform it in a TrajectoryStateOnSurface
0302 
0303   PTrajectoryStateOnDet seedTSOS = trajectoryStateTransform::persistentState(tsos, id.rawId());
0304 
0305   edm::OwnVector<TrackingRecHit> container;
0306   for (unsigned l = 0; l < seg.size(); l++) {
0307     container.push_back(seg[l]->hit()->clone());
0308     //container.push_back(seg[l]->hit());
0309   }
0310 
0311   TrajectorySeed theSeed(seedTSOS, container, alongMomentum);
0312   return theSeed;
0313 }
0314 
0315 /*
0316  * estimatePtCSC
0317  *
0318  * Look at delta phi to determine pt as:
0319  * pt = (c_1 * eta + c_2) / dphi
0320  *
0321  * Note that only segment pairs with one segment in the ME1 layers give sensitive results
0322  *
0323  */
0324 void MuonSeedCreator::estimatePtCSC(const SegmentContainer& seg,
0325                                     const std::vector<int>& layers,
0326                                     double& thePt,
0327                                     double& theSpt) {
0328   unsigned size = seg.size();
0329   if (size < 2)
0330     return;
0331 
0332   // reverse the segment and layer container first for pure CSC case
0333   //if ( layers[0] > layers[ layers.size()-1 ] ) {
0334   //   reverse( layers.begin(), layers.end() );
0335   //   reverse( seg.begin(), seg.end() );
0336   //}
0337 
0338   std::vector<double> ptEstimate;
0339   std::vector<double> sptEstimate;
0340 
0341   thePt = defaultMomentum;
0342   theSpt = defaultMomentum;
0343 
0344   double pt = 0.;
0345   double spt = 0.;
0346   GlobalPoint segPos[2];
0347 
0348   int layer0 = layers[0];
0349   segPos[0] = seg[0]->globalPosition();
0350   float eta = fabs(segPos[0].eta());
0351   //float corr = fabs( tan(segPos[0].theta()) );
0352   // use pt from vertex information
0353   /*
0354   if ( layer0 == 0 ) {
0355      SegmentContainer seg0;
0356      seg0.push_back(seg[0]);
0357      std::vector<int> lyr0(1,0);
0358      estimatePtSingle( seg0, lyr0, thePt, theSpt);
0359      ptEstimate.push_back( thePt );
0360      sptEstimate.push_back( theSpt );
0361   }
0362   */
0363 
0364   //std::cout<<" estimate CSC "<<std::endl;
0365 
0366   unsigned idx1 = size;
0367   if (size > 1) {
0368     while (idx1 > 1) {
0369       idx1--;
0370       int layer1 = layers[idx1];
0371       if (layer0 == layer1)
0372         continue;
0373       segPos[1] = seg[idx1]->globalPosition();
0374 
0375       double dphi = segPos[0].phi() - segPos[1].phi();
0376       //double temp_dphi = dphi/corr;
0377       double temp_dphi = dphi;
0378 
0379       double sign = 1.0;
0380       if (temp_dphi < 0.) {
0381         temp_dphi = -1.0 * temp_dphi;
0382         sign = -1.0;
0383       }
0384 
0385       // Ensure that delta phi is not too small to prevent pt from blowing up
0386       if (temp_dphi < 0.0001) {
0387         temp_dphi = 0.0001;
0388         pt = theMaxMomentum;
0389         spt = theMaxMomentum * 0.25;
0390         ptEstimate.push_back(pt * sign);
0391         sptEstimate.push_back(spt);
0392       }
0393       // ME1 is inner-most
0394       if (layer0 == 0 && temp_dphi >= 0.0001) {
0395         // ME1/2 is outer-most
0396         if (layer1 == 1) {
0397           //temp_dphi = scaledPhi(temp_dphi, CSC01_1[3] );
0398           pt = getPt(CSC01, eta, temp_dphi)[0];
0399           spt = getPt(CSC01, eta, temp_dphi)[1];
0400         }
0401         // ME2 is outer-most
0402         else if (layer1 == 2) {
0403           //temp_dphi = scaledPhi(temp_dphi, CSC12_3[3] );
0404           pt = getPt(CSC02, eta, temp_dphi)[0];
0405           spt = getPt(CSC02, eta, temp_dphi)[1];
0406         }
0407         // ME3 is outer-most
0408         else if (layer1 == 3) {
0409           //temp_dphi = scaledPhi(temp_dphi, CSC13_3[3] );
0410           pt = getPt(CSC03, eta, temp_dphi)[0];
0411           spt = getPt(CSC03, eta, temp_dphi)[1];
0412         }
0413         // ME4 is outer-most
0414         else {
0415           //temp_dphi = scaledPhi(temp_dphi, CSC14_3[3]);
0416           pt = getPt(CSC14, eta, temp_dphi)[0];
0417           spt = getPt(CSC14, eta, temp_dphi)[1];
0418         }
0419         ptEstimate.push_back(pt * sign);
0420         sptEstimate.push_back(spt);
0421       }
0422 
0423       // ME1/2,ME1/3 is inner-most
0424       if (layer0 == 1) {
0425         // ME2 is outer-most
0426         if (layer1 == 2) {
0427           //if ( eta <= 1.2 )  {   temp_dphi = scaledPhi(temp_dphi, CSC12_1[3]); }
0428           //if ( eta >  1.2 )  {   temp_dphi = scaledPhi(temp_dphi, CSC12_2[3]); }
0429           pt = getPt(CSC12, eta, temp_dphi)[0];
0430           spt = getPt(CSC12, eta, temp_dphi)[1];
0431         }
0432         // ME3 is outer-most
0433         else if (layer1 == 3) {
0434           temp_dphi = scaledPhi(temp_dphi, CSC13_2[3]);
0435           pt = getPt(CSC13, eta, temp_dphi)[0];
0436           spt = getPt(CSC13, eta, temp_dphi)[1];
0437         }
0438         // ME4 is outer-most
0439         else {
0440           temp_dphi = scaledPhi(temp_dphi, CSC14_3[3]);
0441           pt = getPt(CSC14, eta, temp_dphi)[0];
0442           spt = getPt(CSC14, eta, temp_dphi)[1];
0443         }
0444         ptEstimate.push_back(pt * sign);
0445         sptEstimate.push_back(spt);
0446       }
0447 
0448       // ME2 is inner-most
0449       if (layer0 == 2 && temp_dphi > 0.0001) {
0450         // ME4 is outer-most
0451         if (layer1 == 4) {
0452           temp_dphi = scaledPhi(temp_dphi, CSC24_1[3]);
0453           pt = getPt(CSC24, eta, temp_dphi)[0];
0454           spt = getPt(CSC24, eta, temp_dphi)[1];
0455         }
0456         // ME3 is outer-most
0457         else {
0458           // if ME2-4 is availabe , discard ME2-3
0459           if (eta <= 1.7) {
0460             temp_dphi = scaledPhi(temp_dphi, CSC23_1[3]);
0461           }
0462           if (eta > 1.7) {
0463             temp_dphi = scaledPhi(temp_dphi, CSC23_2[3]);
0464           }
0465           pt = getPt(CSC23, eta, temp_dphi)[0];
0466           spt = getPt(CSC23, eta, temp_dphi)[1];
0467         }
0468         ptEstimate.push_back(pt * sign);
0469         sptEstimate.push_back(spt);
0470       }
0471 
0472       // ME3 is inner-most
0473       if (layer0 == 3 && temp_dphi > 0.0001) {
0474         temp_dphi = scaledPhi(temp_dphi, CSC34_1[3]);
0475         pt = getPt(CSC34, eta, temp_dphi)[0];
0476         spt = getPt(CSC34, eta, temp_dphi)[1];
0477         ptEstimate.push_back(pt * sign);
0478         sptEstimate.push_back(spt);
0479       }
0480     }
0481   }
0482 
0483   // Compute weighted average if have more than one estimator
0484   if (!ptEstimate.empty())
0485     weightedPt(ptEstimate, sptEstimate, thePt, theSpt);
0486 }
0487 
0488 /*
0489  * estimatePtDT
0490  *
0491  * Look at delta phi between segments to determine pt as:
0492  * pt = (c_1 * eta + c_2) / dphi
0493  */
0494 void MuonSeedCreator::estimatePtDT(const SegmentContainer& seg,
0495                                    const std::vector<int>& layers,
0496                                    double& thePt,
0497                                    double& theSpt) {
0498   unsigned size = seg.size();
0499   if (size < 2)
0500     return;
0501 
0502   std::vector<double> ptEstimate;
0503   std::vector<double> sptEstimate;
0504 
0505   thePt = defaultMomentum;
0506   theSpt = defaultMomentum;
0507 
0508   double pt = 0.;
0509   double spt = 0.;
0510   GlobalPoint segPos[2];
0511 
0512   int layer0 = layers[0];
0513   segPos[0] = seg[0]->globalPosition();
0514   float eta = fabs(segPos[0].eta());
0515 
0516   //std::cout<<" estimate DT "<<std::endl;
0517   // inner-most layer
0518   //for ( unsigned idx0 = 0; idx0 < size-1; ++idx0 ) {
0519   //  layer0 = layers[idx0];
0520   //  segPos[0]  = seg[idx0]->globalPosition();
0521   // outer-most layer
0522   // for ( unsigned idx1 = idx0+1; idx1 < size; ++idx1 ) {
0523   for (unsigned idx1 = 1; idx1 < size; ++idx1) {
0524     int layer1 = layers[idx1];
0525     segPos[1] = seg[idx1]->globalPosition();
0526 
0527     //eta = fabs(segPos[1].eta());
0528     //if (layer1 == -4) eta = fabs(segPos[0].eta());
0529 
0530     double dphi = segPos[0].phi() - segPos[1].phi();
0531     double temp_dphi = dphi;
0532 
0533     // Ensure that delta phi is not too small to prevent pt from blowing up
0534 
0535     double sign = 1.0;
0536     if (temp_dphi < 0.) {
0537       temp_dphi = -temp_dphi;
0538       sign = -1.0;
0539     }
0540 
0541     if (temp_dphi < 0.0001) {
0542       temp_dphi = 0.0001;
0543       pt = theMaxMomentum;
0544       spt = theMaxMomentum * 0.25;
0545       ptEstimate.push_back(pt * sign);
0546       sptEstimate.push_back(spt);
0547     }
0548 
0549     // MB1 is inner-most
0550     if (layer0 == -1 && temp_dphi > 0.0001) {
0551       // MB2 is outer-most
0552       if (layer1 == -2) {
0553         if (eta <= 0.7) {
0554           temp_dphi = scaledPhi(temp_dphi, DT12_1[3]);
0555         }
0556         if (eta > 0.7) {
0557           temp_dphi = scaledPhi(temp_dphi, DT12_2[3]);
0558         }
0559         pt = getPt(DT12, eta, temp_dphi)[0];
0560         spt = getPt(DT12, eta, temp_dphi)[1];
0561       }
0562       // MB3 is outer-most
0563       else if (layer1 == -3) {
0564         if (eta <= 0.6) {
0565           temp_dphi = scaledPhi(temp_dphi, DT13_1[3]);
0566         }
0567         if (eta > 0.6) {
0568           temp_dphi = scaledPhi(temp_dphi, DT13_2[3]);
0569         }
0570         pt = getPt(DT13, eta, temp_dphi)[0];
0571         spt = getPt(DT13, eta, temp_dphi)[1];
0572       }
0573       // MB4 is outer-most
0574       else {
0575         if (eta <= 0.52) {
0576           temp_dphi = scaledPhi(temp_dphi, DT14_1[3]);
0577         }
0578         if (eta > 0.52) {
0579           temp_dphi = scaledPhi(temp_dphi, DT14_2[3]);
0580         }
0581         pt = getPt(DT14, eta, temp_dphi)[0];
0582         spt = getPt(DT14, eta, temp_dphi)[1];
0583       }
0584       ptEstimate.push_back(pt * sign);
0585       sptEstimate.push_back(spt);
0586     }
0587 
0588     // MB2 is inner-most
0589     if (layer0 == -2 && temp_dphi > 0.0001) {
0590       // MB3 is outer-most
0591       if (layer1 == -3) {
0592         if (eta <= 0.6) {
0593           temp_dphi = scaledPhi(temp_dphi, DT23_1[3]);
0594         }
0595         if (eta > 0.6) {
0596           temp_dphi = scaledPhi(temp_dphi, DT23_2[3]);
0597         }
0598         pt = getPt(DT23, eta, temp_dphi)[0];
0599         spt = getPt(DT23, eta, temp_dphi)[1];
0600       }
0601       // MB4 is outer-most
0602       else {
0603         if (eta <= 0.52) {
0604           temp_dphi = scaledPhi(temp_dphi, DT24_1[3]);
0605         }
0606         if (eta > 0.52) {
0607           temp_dphi = scaledPhi(temp_dphi, DT24_2[3]);
0608         }
0609         pt = getPt(DT24, eta, temp_dphi)[0];
0610         spt = getPt(DT24, eta, temp_dphi)[1];
0611       }
0612       ptEstimate.push_back(pt * sign);
0613       sptEstimate.push_back(spt);
0614     }
0615 
0616     // MB3 is inner-most    -> only marginally useful to pick up the charge
0617     if (layer0 == -3 && temp_dphi > 0.0001) {
0618       // MB4 is outer-most
0619 
0620       if (eta <= 0.51) {
0621         temp_dphi = scaledPhi(temp_dphi, DT34_1[3]);
0622       }
0623       if (eta > 0.51) {
0624         temp_dphi = scaledPhi(temp_dphi, DT34_2[3]);
0625       }
0626       pt = getPt(DT34, eta, temp_dphi)[0];
0627       spt = getPt(DT34, eta, temp_dphi)[1];
0628       ptEstimate.push_back(pt * sign);
0629       sptEstimate.push_back(spt);
0630     }
0631   }
0632   //}
0633 
0634   // Compute weighted average if have more than one estimator
0635   if (!ptEstimate.empty())
0636     weightedPt(ptEstimate, sptEstimate, thePt, theSpt);
0637 }
0638 
0639 /*
0640  * estimatePtOverlap
0641  *
0642  */
0643 void MuonSeedCreator::estimatePtOverlap(const SegmentContainer& seg,
0644                                         const std::vector<int>& layers,
0645                                         double& thePt,
0646                                         double& theSpt) {
0647   int size = layers.size();
0648 
0649   thePt = defaultMomentum;
0650   theSpt = defaultMomentum;
0651 
0652   SegmentContainer segCSC;
0653   std::vector<int> layersCSC;
0654   SegmentContainer segDT;
0655   std::vector<int> layersDT;
0656 
0657   // DT layers are numbered as -4 to -1, whereas CSC layers go from 0 to 4:
0658   for (unsigned j = 0; j < layers.size(); ++j) {
0659     if (layers[j] > -1) {
0660       segCSC.push_back(seg[j]);
0661       layersCSC.push_back(layers[j]);
0662     } else {
0663       segDT.push_back(seg[j]);
0664       layersDT.push_back(layers[j]);
0665     }
0666   }
0667 
0668   std::vector<double> ptEstimate;
0669   std::vector<double> sptEstimate;
0670 
0671   GlobalPoint segPos[2];
0672   int layer0 = layers[0];
0673   segPos[0] = seg[0]->globalPosition();
0674   float eta = fabs(segPos[0].eta());
0675   //std::cout<<" estimate OL "<<std::endl;
0676 
0677   if (!segDT.empty() && !segCSC.empty()) {
0678     int layer1 = layers[size - 1];
0679     segPos[1] = seg[size - 1]->globalPosition();
0680 
0681     double dphi = segPos[0].phi() - segPos[1].phi();
0682     double temp_dphi = dphi;
0683 
0684     // Ensure that delta phi is not too small to prevent pt from blowing up
0685 
0686     double sign = 1.0;
0687     if (temp_dphi < 0.) {
0688       temp_dphi = -temp_dphi;
0689       sign = -1.0;
0690     }
0691 
0692     if (temp_dphi < 0.0001) {
0693       temp_dphi = 0.0001;
0694       thePt = theMaxMomentum;
0695       theSpt = theMaxMomentum * 0.25;
0696       ptEstimate.push_back(thePt * sign);
0697       sptEstimate.push_back(theSpt);
0698     }
0699 
0700     // MB1 is inner-most
0701     if (layer0 == -1 && temp_dphi > 0.0001) {
0702       // ME1/3 is outer-most
0703       if (layer1 == 1) {
0704         temp_dphi = scaledPhi(temp_dphi, OL_1213[3]);
0705         thePt = getPt(OL1213, eta, temp_dphi)[0];
0706         theSpt = getPt(OL1213, eta, temp_dphi)[1];
0707       }
0708       // ME2 is outer-most
0709       else if (layer1 == 2) {
0710         temp_dphi = scaledPhi(temp_dphi, OL_1222[3]);
0711         thePt = getPt(OL1222, eta, temp_dphi)[0];
0712         theSpt = getPt(OL1222, eta, temp_dphi)[1];
0713       }
0714       // ME3 is outer-most
0715       else {
0716         temp_dphi = scaledPhi(temp_dphi, OL_1232[3]);
0717         thePt = getPt(OL1232, eta, temp_dphi)[0];
0718         theSpt = getPt(OL1232, eta, temp_dphi)[1];
0719       }
0720       ptEstimate.push_back(thePt * sign);
0721       sptEstimate.push_back(theSpt);
0722     }
0723     // MB2 is inner-most
0724     if (layer0 == -2 && temp_dphi > 0.0001) {
0725       // ME1/3 is outer-most
0726       if (layer1 == 1) {
0727         temp_dphi = scaledPhi(temp_dphi, OL_2213[3]);
0728         thePt = getPt(OL2213, eta, temp_dphi)[0];
0729         theSpt = getPt(OL2213, eta, temp_dphi)[1];
0730         ptEstimate.push_back(thePt * sign);
0731         sptEstimate.push_back(theSpt);
0732       }
0733       // ME2 is outer-most
0734       if (layer1 == 2) {
0735         temp_dphi = scaledPhi(temp_dphi, OL_2222[3]);
0736         thePt = getPt(OL2222, eta, temp_dphi)[0];
0737         theSpt = getPt(OL2222, eta, temp_dphi)[1];
0738       }
0739     }
0740   }
0741 
0742   if (segDT.size() > 1) {
0743     estimatePtDT(segDT, layersDT, thePt, theSpt);
0744     ptEstimate.push_back(thePt);
0745     sptEstimate.push_back(theSpt);
0746   }
0747 
0748   /*
0749   // not useful ....and pt estimation is bad
0750   if ( segCSC.size() > 1 ) {
0751     // don't estimate pt without ME1 information
0752     bool CSCLayer1=false;
0753     for (unsigned i=0; i< layersCSC.size(); i++) {
0754         if ( layersCSC[i]==0 || layersCSC[i]==1 ) CSCLayer1 = true;
0755     }
0756     if (CSCLayer1) {
0757        estimatePtCSC(segCSC, layersCSC, thePt, theSpt);
0758        ptEstimate.push_back(thePt);
0759        sptEstimate.push_back(theSpt);
0760     }
0761   }
0762   */
0763 
0764   // Compute weighted average if have more than one estimator
0765   if (!ptEstimate.empty())
0766     weightedPt(ptEstimate, sptEstimate, thePt, theSpt);
0767 }
0768 /*
0769  *
0770  *   estimate Pt for single segment events
0771  *
0772  */
0773 void MuonSeedCreator::estimatePtSingle(const SegmentContainer& seg,
0774                                        const std::vector<int>& layers,
0775                                        double& thePt,
0776                                        double& theSpt) {
0777   thePt = defaultMomentum;
0778   theSpt = defaultMomentum;
0779 
0780   GlobalPoint segPos = seg[0]->globalPosition();
0781   double eta = segPos.eta();
0782   GlobalVector gv = seg[0]->globalDirection();
0783 
0784   // Psi is angle between the segment origin and segment direction
0785   // Use dot product between two vectors to get Psi in global x-y plane
0786   double cosDpsi = (gv.x() * segPos.x() + gv.y() * segPos.y());
0787   cosDpsi /= sqrt(segPos.x() * segPos.x() + segPos.y() * segPos.y());
0788   cosDpsi /= sqrt(gv.x() * gv.x() + gv.y() * gv.y());
0789 
0790   double axb = (segPos.x() * gv.y()) - (segPos.y() * gv.x());
0791   double sign = (axb < 0.) ? 1.0 : -1.0;
0792 
0793   double dpsi = fabs(acos(cosDpsi));
0794   if (dpsi > 1.570796) {
0795     dpsi = 3.141592 - dpsi;
0796     sign = -1. * sign;
0797   }
0798   if (fabs(dpsi) < 0.00005) {
0799     dpsi = 0.00005;
0800   }
0801 
0802   // the 1st layer
0803   if (layers[0] == -1) {
0804     // MB10
0805     if (fabs(eta) < 0.3) {
0806       dpsi = scaledPhi(dpsi, SMB_10S[3]);
0807       thePt = getPt(SMB10, eta, dpsi)[0];
0808       theSpt = getPt(SMB10, eta, dpsi)[1];
0809     }
0810     // MB11
0811     if (fabs(eta) >= 0.3 && fabs(eta) < 0.82) {
0812       dpsi = scaledPhi(dpsi, SMB_11S[3]);
0813       thePt = getPt(SMB11, eta, dpsi)[0];
0814       theSpt = getPt(SMB11, eta, dpsi)[1];
0815     }
0816     // MB12
0817     if (fabs(eta) >= 0.82 && fabs(eta) < 1.2) {
0818       dpsi = scaledPhi(dpsi, SMB_12S[3]);
0819       thePt = getPt(SMB12, eta, dpsi)[0];
0820       theSpt = getPt(SMB12, eta, dpsi)[1];
0821     }
0822   }
0823   if (layers[0] == 1) {
0824     // ME13
0825     if (fabs(eta) > 0.92 && fabs(eta) < 1.16) {
0826       dpsi = scaledPhi(dpsi, SME_13S[3]);
0827       thePt = getPt(SME13, eta, dpsi)[0];
0828       theSpt = getPt(SME13, eta, dpsi)[1];
0829     }
0830     // ME12
0831     if (fabs(eta) >= 1.16 && fabs(eta) <= 1.6) {
0832       dpsi = scaledPhi(dpsi, SME_12S[3]);
0833       thePt = getPt(SME12, eta, dpsi)[0];
0834       theSpt = getPt(SME12, eta, dpsi)[1];
0835     }
0836   }
0837   if (layers[0] == 0) {
0838     // ME11
0839     if (fabs(eta) > 1.6) {
0840       dpsi = scaledPhi(dpsi, SMB_11S[3]);
0841       thePt = getPt(SME11, eta, dpsi)[0];
0842       theSpt = getPt(SME11, eta, dpsi)[1];
0843     }
0844   }
0845   // the 2nd layer
0846   if (layers[0] == -2) {
0847     // MB20
0848     if (fabs(eta) < 0.25) {
0849       dpsi = scaledPhi(dpsi, SMB_20S[3]);
0850       thePt = getPt(SMB20, eta, dpsi)[0];
0851       theSpt = getPt(SMB20, eta, dpsi)[1];
0852     }
0853     // MB21
0854     if (fabs(eta) >= 0.25 && fabs(eta) < 0.72) {
0855       dpsi = scaledPhi(dpsi, SMB_21S[3]);
0856       thePt = getPt(SMB21, eta, dpsi)[0];
0857       theSpt = getPt(SMB21, eta, dpsi)[1];
0858     }
0859     // MB22
0860     if (fabs(eta) >= 0.72 && fabs(eta) < 1.04) {
0861       dpsi = scaledPhi(dpsi, SMB_22S[3]);
0862       thePt = getPt(SMB22, eta, dpsi)[0];
0863       theSpt = getPt(SMB22, eta, dpsi)[1];
0864     }
0865   }
0866   if (layers[0] == 2) {
0867     // ME22
0868     if (fabs(eta) > 0.95 && fabs(eta) <= 1.6) {
0869       dpsi = scaledPhi(dpsi, SME_22S[3]);
0870       thePt = getPt(SME22, eta, dpsi)[0];
0871       theSpt = getPt(SME22, eta, dpsi)[1];
0872     }
0873     // ME21
0874     if (fabs(eta) > 1.6 && fabs(eta) < 2.45) {
0875       dpsi = scaledPhi(dpsi, SME_21S[3]);
0876       thePt = getPt(SME21, eta, dpsi)[0];
0877       theSpt = getPt(SME21, eta, dpsi)[1];
0878     }
0879   }
0880 
0881   // the 3rd layer
0882   if (layers[0] == -3) {
0883     // MB30
0884     if (fabs(eta) <= 0.22) {
0885       dpsi = scaledPhi(dpsi, SMB_30S[3]);
0886       thePt = getPt(SMB30, eta, dpsi)[0];
0887       theSpt = getPt(SMB30, eta, dpsi)[1];
0888     }
0889     // MB31
0890     if (fabs(eta) > 0.22 && fabs(eta) <= 0.6) {
0891       dpsi = scaledPhi(dpsi, SMB_31S[3]);
0892       thePt = getPt(SMB31, eta, dpsi)[0];
0893       theSpt = getPt(SMB31, eta, dpsi)[1];
0894     }
0895     // MB32
0896     if (fabs(eta) > 0.6 && fabs(eta) < 0.95) {
0897       dpsi = scaledPhi(dpsi, SMB_32S[3]);
0898       thePt = getPt(SMB32, eta, dpsi)[0];
0899       theSpt = getPt(SMB32, eta, dpsi)[1];
0900     }
0901   }
0902   thePt = fabs(thePt) * sign;
0903   theSpt = fabs(theSpt);
0904 
0905   return;
0906 }
0907 
0908 // setup the minimum value for obvious showering cases
0909 void MuonSeedCreator::estimatePtShowering(int& NShowers, int& NShowerSegments, double& thePt, double& theSpt) {
0910   if (NShowers > 2 && thePt < 300.) {
0911     thePt = 800.;
0912     theSpt = 200.;
0913   }
0914   if (NShowers == 2 && NShowerSegments > 11 && thePt < 150.) {
0915     thePt = 280.;
0916     theSpt = 70.;
0917   }
0918   if (NShowers == 2 && NShowerSegments <= 11 && thePt < 50.) {
0919     thePt = 80.;
0920     theSpt = 40.;
0921   }
0922   if (NShowers == 1 && NShowerSegments <= 5 && thePt < 10.) {
0923     thePt = 16.;
0924     theSpt = 8.;
0925   }
0926 }
0927 
0928 /*
0929  * weightedPt
0930  *
0931  * Look at delta phi between segments to determine pt as:
0932  * pt = (c_1 * eta + c_2) / dphi
0933  */
0934 void MuonSeedCreator::weightedPt(const std::vector<double>& ptEstimate,
0935                                  const std::vector<double>& sptEstimate,
0936                                  double& thePt,
0937                                  double& theSpt) {
0938   int size = ptEstimate.size();
0939 
0940   // If only one element, by-pass computation below
0941   if (size < 2) {
0942     thePt = ptEstimate[0];
0943     theSpt = sptEstimate[0];
0944     return;
0945   }
0946 
0947   double charge = 0.;
0948   // If have more than one pt estimator, first look if any estimator is biased
0949   // by comparing the charge of each estimator
0950 
0951   for (unsigned j = 0; j < ptEstimate.size(); j++) {
0952     //std::cout<<" weighting pt: "<< ptEstimate[j] <<std::endl;
0953     if (ptEstimate[j] < 0.) {
0954       // To prevent from blowing up, add 0.1
0955       // weight by relative error on pt
0956       charge -= 1. * (ptEstimate[j] * ptEstimate[j]) / (sptEstimate[j] * sptEstimate[j]);
0957     } else {
0958       // weight by relative error on pt
0959       charge += 1. * (ptEstimate[j] * ptEstimate[j]) / (sptEstimate[j] * sptEstimate[j]);
0960     }
0961   }
0962 
0963   // No need to normalize as we want to know only sign ( + or - )
0964   if (charge < 0.) {
0965     charge = -1.;
0966   } else {
0967     charge = 1.;
0968   }
0969 
0970   //int n = 0;
0971   double weightPtSum = 0.;
0972   double sigmaSqr_sum = 0.;
0973 
0974   // Now, we want to compute average Pt using estimators with "correct" charge
0975   // This is to remove biases
0976   for (unsigned j = 0; j < ptEstimate.size(); ++j) {
0977     //if ( (minpt_ratio < 0.5) && (fabs(ptEstimate[j]) < 5.0) ) continue;
0978     //if ( ptEstimate[j] * charge > 0. ) {
0979     //n++;
0980     sigmaSqr_sum += 1.0 / (sptEstimate[j] * sptEstimate[j]);
0981     weightPtSum += fabs(ptEstimate[j]) / (sptEstimate[j] * sptEstimate[j]);
0982     //}
0983   }
0984   /*  
0985   if (n < 1) {
0986     thePt  = defaultMomentum*charge;
0987     theSpt = defaultMomentum; 
0988     return;
0989   } 
0990   */
0991   // Compute weighted mean and error
0992 
0993   thePt = (charge * weightPtSum) / sigmaSqr_sum;
0994   theSpt = sqrt(1.0 / sigmaSqr_sum);
0995 
0996   //std::cout<<" final weighting : "<< thePt <<" ~ "<< fabs( theSpt/thePt ) <<std::endl;
0997 
0998   return;
0999 }
1000 
1001 std::vector<double> MuonSeedCreator::getPt(const std::vector<double>& vPara, double eta, double dPhi) {
1002   double h = fabs(eta);
1003   double estPt = (vPara[0] + vPara[1] * h + vPara[2] * h * h) / dPhi;
1004   double estSPt = (vPara[3] + vPara[4] * h + vPara[5] * h * h) * estPt;
1005   std::vector<double> paraPt;
1006   paraPt.push_back(estPt);
1007   paraPt.push_back(estSPt);
1008 
1009   //std::cout<<"      pt:"<<estPt<<" +/-"<< estSPt<<"   h:"<<eta<<"  df:"<<dPhi<<std::endl;
1010   return paraPt;
1011 }
1012 
1013 double MuonSeedCreator::scaledPhi(double dphi, double t1) {
1014   if (dphi != 0.) {
1015     double oPhi = 1. / dphi;
1016     dphi = dphi / (1. + t1 / (oPhi + 10.));
1017     return dphi;
1018 
1019   } else {
1020     return dphi;
1021   }
1022 }
[ 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!