Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 14:26:07

0001 /** \class SETSeedFinder
0002     I. Bloch, E. James, S. Stoynev
0003  */
0004 
0005 #include "RecoMuon/MuonSeedGenerator/src/SETSeedFinder.h"
0006 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0007 #include "FWCore/Framework/interface/Event.h"
0008 #include "RecoMuon/TrackingTools/interface/MuonPatternRecoDumper.h"
0009 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0010 #include "TrackingTools/DetLayers/interface/DetLayer.h"
0011 #include "Geometry/CommonDetUnit/interface/GeomDet.h"
0012 #include "DataFormats/MuonDetId/interface/MuonSubdetId.h"
0013 #include "DataFormats/MuonDetId/interface/CSCDetId.h"
0014 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
0015 
0016 #include "TMath.h"
0017 
0018 using namespace edm;
0019 using namespace std;
0020 
0021 const string metname = "Muon|RecoMuon|SETMuonSeedFinder";
0022 
0023 SETSeedFinder::SETSeedFinder(const ParameterSet& parameterSet) : MuonSeedVFinder() {
0024   // Parameter set for the Builder
0025   ParameterSet trajectoryBuilderParameters = parameterSet.getParameter<ParameterSet>("SETTrajBuilderParameters");
0026   apply_prePruning = trajectoryBuilderParameters.getParameter<bool>("Apply_prePruning");
0027   // load pT seed parameters
0028   thePtExtractor = new MuonSeedPtExtractor(trajectoryBuilderParameters);
0029 }
0030 
0031 void SETSeedFinder::seeds(const MuonRecHitContainer& cluster, std::vector<TrajectorySeed>& result) {}
0032 
0033 // there is an existing sorter somewhere in the CMSSW code (I think) - delete that
0034 namespace {
0035   struct sorter {
0036     sorter() {}
0037     bool operator()(MuonTransientTrackingRecHit::MuonRecHitPointer const& hit_1,
0038                     MuonTransientTrackingRecHit::MuonRecHitPointer const& hit_2) const {
0039       return (hit_1->globalPosition().mag2() < hit_2->globalPosition().mag2());
0040     }
0041   };  // smaller first
0042   const sorter sortSegRadius;
0043 }  // namespace
0044 
0045 std::vector<SETSeedFinder::MuonRecHitContainer> SETSeedFinder::sortByLayer(MuonRecHitContainer& cluster) const {
0046   stable_sort(cluster.begin(), cluster.end(), sortSegRadius);
0047   //---- group hits in detector layers (if in same layer); the idea is that
0048   //---- some hits could not belong to a track simultaneously - these will be in a
0049   //---- group; two hits from one and the same group will not go to the same track
0050   std::vector<MuonRecHitContainer> MuonRecHitContainer_perLayer;
0051   if (!cluster.empty()) {
0052     int iHit = 0;
0053     MuonRecHitContainer hitsInThisLayer;
0054     hitsInThisLayer.push_back(cluster[iHit]);
0055     DetId detId = cluster[iHit]->hit()->geographicalId();
0056     const GeomDet* geomDet = theService->trackingGeometry()->idToDet(detId);
0057     while (iHit < int(cluster.size()) - 1) {
0058       DetId detId_2 = cluster[iHit + 1]->hit()->geographicalId();
0059       const GlobalPoint gp_nextHit = cluster[iHit + 1]->globalPosition();
0060 
0061       // this is the distance of the "second" hit to the "first" detector (containing the "first hit")
0062       float distanceToDetector = fabs(geomDet->surface().localZ(gp_nextHit));
0063 
0064       //---- hits from DT and CSC  could be very close in angle but incosistent with
0065       //---- belonging to a common track (and these are different surfaces);
0066       //---- also DT (and CSC now - 090822) hits from a station (in a pre-cluster) should be always in a group together;
0067       //---- take this into account and put such hits in a group together
0068 
0069       bool specialCase = false;
0070       if (detId.subdetId() == MuonSubdetId::DT && detId_2.subdetId() == MuonSubdetId::DT) {
0071         DTChamberId dtCh(detId);
0072         DTChamberId dtCh_2(detId_2);
0073         specialCase = (dtCh.station() == dtCh_2.station());
0074       } else if (detId.subdetId() == MuonSubdetId::CSC && detId_2.subdetId() == MuonSubdetId::CSC) {
0075         CSCDetId cscCh(detId);
0076         CSCDetId cscCh_2(detId_2);
0077         specialCase = (cscCh.station() == cscCh_2.station() && cscCh.ring() == cscCh_2.ring());
0078       }
0079 
0080       if (distanceToDetector < 0.001 || true == specialCase) {  // hardcoded value - remove!
0081         hitsInThisLayer.push_back(cluster[iHit + 1]);
0082       } else {
0083         specialCase = false;
0084         if (((cluster[iHit]->isDT() && cluster[iHit + 1]->isCSC()) ||
0085              (cluster[iHit]->isCSC() && cluster[iHit + 1]->isDT())) &&
0086             //---- what is the minimal distance between a DT and a CSC hit belonging
0087             //---- to a common track? (well, with "reasonable" errors; put 10 cm for now)
0088             fabs(cluster[iHit + 1]->globalPosition().mag() - cluster[iHit]->globalPosition().mag()) < 10.) {
0089           hitsInThisLayer.push_back(cluster[iHit + 1]);
0090           // change to Stoyan - now we also update the detID here... give it a try. IBL 080905
0091           detId = cluster[iHit + 1]->hit()->geographicalId();
0092           geomDet = theService->trackingGeometry()->idToDet(detId);
0093         } else if (!specialCase) {
0094           //---- put the group of hits in the vector (containing the groups of hits)
0095           //---- and continue with next layer (group)
0096           MuonRecHitContainer_perLayer.push_back(hitsInThisLayer);
0097           hitsInThisLayer.clear();
0098           hitsInThisLayer.push_back(cluster[iHit + 1]);
0099           detId = cluster[iHit + 1]->hit()->geographicalId();
0100           geomDet = theService->trackingGeometry()->idToDet(detId);
0101         }
0102       }
0103       ++iHit;
0104     }
0105     MuonRecHitContainer_perLayer.push_back(hitsInThisLayer);
0106   }
0107   return MuonRecHitContainer_perLayer;
0108 }
0109 //
0110 void SETSeedFinder::limitCombinatorics(std::vector<MuonRecHitContainer>& MuonRecHitContainer_perLayer) {
0111   const int maximumNumberOfCombinations = 1000000;
0112   unsigned nLayers = MuonRecHitContainer_perLayer.size();
0113   if (1 == nLayers) {
0114     return;
0115   }
0116   // maximal number of (segment) layers would be upto ~12; see next function
0117   // below is just a quick fix for a rare "overflow"
0118   if (MuonRecHitContainer_perLayer.size() > 15) {
0119     MuonRecHitContainer_perLayer.resize(1);
0120     return;
0121   }
0122 
0123   std::vector<double> sizeOfLayer(nLayers);
0124   //std::cout<<" nLayers = "<<nLayers<<std::endl;
0125   double nAllCombinations = 1.;
0126   for (unsigned int i = 0; i < nLayers; ++i) {
0127     //std::cout<<" i = "<<i<<" size = "<<MuonRecHitContainer_perLayer.at(i).size()<<std::endl;
0128     sizeOfLayer.at(i) = MuonRecHitContainer_perLayer.at(i).size();
0129     nAllCombinations *= MuonRecHitContainer_perLayer.at(i).size();
0130   }
0131   //std::cout<<"nAllCombinations = "<<nAllCombinations<<std::endl;
0132   //---- Erase most busy detector layers until we get less than maximumNumberOfCombinations combinations
0133   int iCycle = 0;
0134   while (nAllCombinations > float(maximumNumberOfCombinations)) {
0135     ++iCycle;
0136     std::vector<double>::iterator maxEl_it = max_element(sizeOfLayer.begin(), sizeOfLayer.end());
0137     int maxEl = maxEl_it - sizeOfLayer.begin();
0138     nAllCombinations /= MuonRecHitContainer_perLayer.at(maxEl).size();
0139     //std::cout<<" iCycle = "<<iCycle<<" nAllCombinations = "<<nAllCombinations<<std::endl;
0140     MuonRecHitContainer_perLayer.erase(MuonRecHitContainer_perLayer.begin() + maxEl);
0141     sizeOfLayer.erase(sizeOfLayer.begin() + maxEl);
0142   }
0143   return;
0144 }
0145 //
0146 std::vector<SETSeedFinder::MuonRecHitContainer> SETSeedFinder::findAllValidSets(
0147     const std::vector<SETSeedFinder::MuonRecHitContainer>& MuonRecHitContainer_perLayer) {
0148   std::vector<MuonRecHitContainer> allValidSets;
0149   // build all possible combinations (i.e valid sets; the algorithm name is after this feature -
0150   // SET algorithm)
0151   //
0152   // ugly... use recursive function?!
0153   // or implement Ingo's suggestion (a la ST)
0154   unsigned nLayers = MuonRecHitContainer_perLayer.size();
0155   if (1 == nLayers) {
0156     return allValidSets;
0157   }
0158   MuonRecHitContainer validSet;
0159   unsigned int iPos0 = 0;
0160   std::vector<unsigned int> iLayer(12);  // could there be more than 11 layers?
0161   std::vector<unsigned int> size(12);
0162   if (iPos0 < nLayers) {
0163     size.at(iPos0) = MuonRecHitContainer_perLayer.at(iPos0).size();
0164     for (iLayer[iPos0] = 0; iLayer[iPos0] < size[iPos0]; ++iLayer[iPos0]) {
0165       validSet.clear();
0166       validSet.push_back(MuonRecHitContainer_perLayer[iPos0][iLayer[iPos0]]);
0167       unsigned int iPos1 = 1;
0168       if (iPos1 < nLayers) {
0169         size.at(iPos1) = MuonRecHitContainer_perLayer.at(iPos1).size();
0170         for (iLayer[iPos1] = 0; iLayer[iPos1] < size[iPos1]; ++iLayer[iPos1]) {
0171           validSet.resize(iPos1);
0172           validSet.push_back(MuonRecHitContainer_perLayer[iPos1][iLayer[iPos1]]);
0173           unsigned int iPos2 = 2;
0174           if (iPos2 < nLayers) {
0175             size.at(iPos2) = MuonRecHitContainer_perLayer.at(iPos2).size();
0176             for (iLayer[iPos2] = 0; iLayer[iPos2] < size[iPos2]; ++iLayer[iPos2]) {
0177               validSet.resize(iPos2);
0178               validSet.push_back(MuonRecHitContainer_perLayer[iPos2][iLayer[iPos2]]);
0179               unsigned int iPos3 = 3;
0180               if (iPos3 < nLayers) {
0181                 size.at(iPos3) = MuonRecHitContainer_perLayer.at(iPos3).size();
0182                 for (iLayer[iPos3] = 0; iLayer[iPos3] < size[iPos3]; ++iLayer[iPos3]) {
0183                   validSet.resize(iPos3);
0184                   validSet.push_back(MuonRecHitContainer_perLayer[iPos3][iLayer[iPos3]]);
0185                   unsigned int iPos4 = 4;
0186                   if (iPos4 < nLayers) {
0187                     size.at(iPos4) = MuonRecHitContainer_perLayer.at(iPos4).size();
0188                     for (iLayer[iPos4] = 0; iLayer[iPos4] < size[iPos4]; ++iLayer[iPos4]) {
0189                       validSet.resize(iPos4);
0190                       validSet.push_back(MuonRecHitContainer_perLayer[iPos4][iLayer[iPos4]]);
0191                       unsigned int iPos5 = 5;
0192                       if (iPos5 < nLayers) {
0193                         size.at(iPos5) = MuonRecHitContainer_perLayer.at(iPos5).size();
0194                         for (iLayer[iPos5] = 0; iLayer[iPos5] < size[iPos5]; ++iLayer[iPos5]) {
0195                           validSet.resize(iPos5);
0196                           validSet.push_back(MuonRecHitContainer_perLayer[iPos5][iLayer[iPos5]]);
0197                           unsigned int iPos6 = 6;
0198                           if (iPos6 < nLayers) {
0199                             size.at(iPos6) = MuonRecHitContainer_perLayer.at(iPos6).size();
0200                             for (iLayer[iPos6] = 0; iLayer[iPos6] < size[iPos6]; ++iLayer[iPos6]) {
0201                               validSet.resize(iPos6);
0202                               validSet.push_back(MuonRecHitContainer_perLayer[iPos6][iLayer[iPos6]]);
0203                               unsigned int iPos7 = 7;
0204                               if (iPos7 < nLayers) {
0205                                 size.at(iPos7) = MuonRecHitContainer_perLayer.at(iPos7).size();
0206                                 for (iLayer[iPos7] = 0; iLayer[iPos7] < size[iPos7]; ++iLayer[iPos7]) {
0207                                   validSet.resize(iPos7);
0208                                   validSet.push_back(MuonRecHitContainer_perLayer[iPos7][iLayer[iPos7]]);
0209                                   unsigned int iPos8 = 8;
0210                                   if (iPos8 < nLayers) {
0211                                     size.at(iPos8) = MuonRecHitContainer_perLayer.at(iPos8).size();
0212                                     for (iLayer[iPos8] = 0; iLayer[iPos8] < size[iPos8]; ++iLayer[iPos8]) {
0213                                       validSet.resize(iPos8);
0214                                       validSet.push_back(MuonRecHitContainer_perLayer[iPos8][iLayer[iPos8]]);
0215                                       unsigned int iPos9 = 9;
0216                                       if (iPos9 < nLayers) {
0217                                         size.at(iPos9) = MuonRecHitContainer_perLayer.at(iPos9).size();
0218                                         for (iLayer[iPos9] = 0; iLayer[iPos9] < size[iPos9]; ++iLayer[iPos9]) {
0219                                           validSet.resize(iPos9);
0220                                           validSet.push_back(MuonRecHitContainer_perLayer[iPos9][iLayer[iPos9]]);
0221                                           unsigned int iPos10 = 10;
0222                                           if (iPos10 < nLayers) {
0223                                             size.at(iPos10) = MuonRecHitContainer_perLayer.at(iPos10).size();
0224                                             for (iLayer[iPos10] = 0; iLayer[iPos10] < size[iPos10]; ++iLayer[iPos10]) {
0225                                               validSet.resize(iPos10);
0226                                               validSet.push_back(MuonRecHitContainer_perLayer[iPos10][iLayer[iPos10]]);
0227                                               unsigned int iPos11 = 11;  // more?
0228                                               if (iPos11 < nLayers) {
0229                                                 size.at(iPos11) = MuonRecHitContainer_perLayer.at(iPos11).size();
0230                                                 for (iLayer[iPos11] = 0; iLayer[iPos11] < size[iPos11];
0231                                                      ++iLayer[iPos11]) {
0232                                                 }
0233                                               } else {
0234                                                 allValidSets.push_back(validSet);
0235                                               }
0236                                             }
0237                                           } else {
0238                                             allValidSets.push_back(validSet);
0239                                           }
0240                                         }
0241                                       } else {
0242                                         allValidSets.push_back(validSet);
0243                                       }
0244                                     }
0245                                   } else {
0246                                     allValidSets.push_back(validSet);
0247                                   }
0248                                 }
0249                               } else {
0250                                 allValidSets.push_back(validSet);
0251                               }
0252                             }
0253                           } else {
0254                             allValidSets.push_back(validSet);
0255                           }
0256                         }
0257                       } else {
0258                         allValidSets.push_back(validSet);
0259                       }
0260                     }
0261                   } else {
0262                     allValidSets.push_back(validSet);
0263                   }
0264                 }
0265               } else {
0266                 allValidSets.push_back(validSet);
0267               }
0268             }
0269           } else {
0270             allValidSets.push_back(validSet);
0271           }
0272         }
0273       } else {
0274         allValidSets.push_back(validSet);
0275       }
0276     }
0277   } else {
0278     allValidSets.push_back(validSet);
0279   }
0280   return allValidSets;
0281 }
0282 
0283 std::pair<int, int>  // or <bool, bool>
0284 SETSeedFinder::checkAngleDeviation(double dPhi_1, double dPhi_2) const {
0285   // Two conditions:
0286   // a) deviations should be only to one side (above some absolute value cut to avoid
0287   //    material effects; this should be refined)
0288   // b) deviatiation in preceding steps should be bigger due to higher magnetic field
0289   //    (again - a minimal value cut should be in place; this also should account for
0290   //     the small (Z) distances in overlaping CSC chambers)
0291 
0292   double mult = dPhi_1 * dPhi_2;
0293   int signVal = 1;
0294   if (fabs(dPhi_1) < fabs(dPhi_2)) {
0295     signVal = -1;
0296   }
0297   int signMult = -1;
0298   if (mult > 0)
0299     signMult = 1;
0300   std::pair<int, int> sign;
0301   sign = make_pair(signVal, signMult);
0302 
0303   return sign;
0304 }
0305 
0306 void SETSeedFinder::validSetsPrePruning(std::vector<SETSeedFinder::MuonRecHitContainer>& allValidSets) {
0307   //---- this actually is a pre-pruning; it does not include any fit information;
0308   //---- it is intended to remove only very "wild" segments from a set;
0309   //---- no "good" segment is to be lost (otherwise - widen the parameters)
0310 
0311   for (unsigned int iSet = 0; iSet < allValidSets.size(); ++iSet) {
0312     pre_prune(allValidSets[iSet]);
0313   }
0314 }
0315 
0316 void SETSeedFinder::pre_prune(SETSeedFinder::MuonRecHitContainer& validSet) const {
0317   unsigned nHits = validSet.size();
0318   if (nHits > 3) {  // to decide we need at least 4 measurements
0319     // any information could be used to make a decision for pruning
0320     // maybe dPhi (delta Phi) is enough
0321     std::vector<double> dPhi;
0322     double dPhi_tmp;
0323     bool wildCandidate;
0324     int pruneHit_tmp;
0325 
0326     for (unsigned int iHit = 1; iHit < nHits; ++iHit) {
0327       dPhi_tmp = validSet[iHit]->globalPosition().phi() - validSet[iHit - 1]->globalPosition().phi();
0328       dPhi.push_back(dPhi_tmp);
0329     }
0330     std::vector<int> pruneHit;
0331     //---- loop over all the hits in a set
0332 
0333     for (unsigned int iHit = 0; iHit < nHits; ++iHit) {
0334       double dPHI_MIN = 0.02;  //?? hardcoded - remove it
0335       if (iHit) {
0336         // if we have to remove the very first hit (iHit == 0) then
0337         // we'll probably be in trouble already
0338         wildCandidate = false;
0339         // actually 2D is bad only if not r-phi... Should I refine it?
0340         // a hit is a candidate for pruning only if dPhi > dPHI_MIN;
0341         // pruning decision is based on combination of hits characteristics
0342         if (4 == validSet[iHit - 1]->dimension() && 4 == validSet[iHit]->dimension() &&
0343             fabs(validSet[iHit]->globalPosition().phi() - validSet[iHit - 1]->globalPosition().phi()) > dPHI_MIN) {
0344           wildCandidate = true;
0345         }
0346         pruneHit_tmp = -1;
0347         if (wildCandidate) {
0348           // OK - this couple doesn't look good (and is from 4D segments); proceed...
0349           if (1 == iHit) {  // the first  and the last hits are special case
0350             if (4 == validSet[iHit + 1]->dimension() && 4 == validSet[iHit + 2]->dimension()) {  //4D?
0351               // is the picture better if we remove the second hit?
0352               dPhi_tmp = validSet[iHit + 1]->globalPosition().phi() - validSet[iHit - 1]->globalPosition().phi();
0353               // is the deviation what we expect (sign, not magnitude)?
0354               std::pair<int, int> sign = checkAngleDeviation(dPhi_tmp, dPhi[2]);
0355               if (1 == sign.first && 1 == sign.second) {
0356                 pruneHit_tmp = iHit;  // mark the hit 1 for removing
0357               }
0358             }
0359           } else if (iHit > 1 && iHit < validSet.size() - 1) {
0360             if (4 == validSet[0]->dimension() &&  // we rely on the first (most important) couple
0361                 4 == validSet[1]->dimension() && pruneHit.back() != int(iHit - 1) &&
0362                 pruneHit.back() != 1) {  // check if hits are already marked
0363               // decide which of the two hits should be removed (if any; preferably the outer one i.e.
0364               // iHit rather than iHit-1); here - check what we get by removing iHit
0365               dPhi_tmp = validSet[iHit + 1]->globalPosition().phi() - validSet[iHit - 1]->globalPosition().phi();
0366               // first couple is most important anyway so again compare to it
0367               std::pair<int, int> sign = checkAngleDeviation(dPhi[0], dPhi_tmp);
0368               if (1 == sign.first && 1 == sign.second) {
0369                 pruneHit_tmp = iHit;  // mark the hit iHit for removing
0370               } else {                // iHit is not to be removed; proceed...
0371                 // what if we remove (iHit - 1) instead of iHit?
0372                 dPhi_tmp = validSet[iHit + 1]->globalPosition().phi() - validSet[iHit]->globalPosition().phi();
0373                 std::pair<int, int> sign = checkAngleDeviation(dPhi[0], dPhi_tmp);
0374                 if (1 == sign.first && 1 == sign.second) {
0375                   pruneHit_tmp = iHit - 1;  // mark the hit (iHit -1) for removing
0376                 }
0377               }
0378             }
0379           } else {
0380             // the last hit: if picture is not good - remove it
0381             if (pruneHit.size() > 1 && pruneHit[pruneHit.size() - 1] < 0 && pruneHit[pruneHit.size() - 2] < 0) {
0382               std::pair<int, int> sign = checkAngleDeviation(dPhi[dPhi.size() - 2], dPhi[dPhi.size() - 1]);
0383               if (-1 == sign.first && -1 == sign.second) {  // here logic is a bit twisted
0384                 pruneHit_tmp = iHit;                        // mark the last hit for removing
0385               }
0386             }
0387           }
0388         }
0389         pruneHit.push_back(pruneHit_tmp);
0390       }
0391     }
0392     // }
0393     // actual pruning
0394     for (unsigned int iHit = 1; iHit < nHits; ++iHit) {
0395       int count = 0;
0396       if (pruneHit[iHit - 1] > 0) {
0397         validSet.erase(validSet.begin() + pruneHit[iHit - 1] - count);
0398         ++count;
0399       }
0400     }
0401   }
0402 }
0403 
0404 std::vector<SeedCandidate> SETSeedFinder::fillSeedCandidates(std::vector<MuonRecHitContainer>& allValidSets) {
0405   //---- we have the valid sets constructed; transform the information in an
0406   //---- apropriate form; meanwhile - estimate the momentum for a given set
0407 
0408   // RPCs should not be used (no parametrization)
0409   std::vector<SeedCandidate> seedCandidates_inCluster;
0410   // calculate and fill the inputs needed
0411   // loop over all valid sets
0412   for (unsigned int iSet = 0; iSet < allValidSets.size(); ++iSet) {
0413     //
0414     //std::cout<<"  This is SET number : "<<iSet<<std::endl;
0415     //for(unsigned int iHit = 0;iHit<allValidSets[iSet].size();++iHit){
0416     //std::cout<<"   measurements in the SET:  iHit = "<<iHit<<" pos = "<<allValidSets[iSet][iHit]->globalPosition()<<
0417     //" dim = "<<allValidSets[iSet][iHit]->dimension()<<std::endl;
0418     //}
0419 
0420     CLHEP::Hep3Vector momEstimate;
0421     int chargeEstimate;
0422     estimateMomentum(allValidSets[iSet], momEstimate, chargeEstimate);
0423     MuonRecHitContainer MuonRecHitContainer_theSet_prep;
0424     // currently hardcoded - will be in proper loop of course:
0425 
0426     SeedCandidate seedCandidates_inCluster_prep;
0427     seedCandidates_inCluster_prep.theSet = allValidSets[iSet];
0428     seedCandidates_inCluster_prep.momentum = momEstimate;
0429     seedCandidates_inCluster_prep.charge = chargeEstimate;
0430     seedCandidates_inCluster.push_back(seedCandidates_inCluster_prep);
0431     // END estimateMomentum
0432   }
0433   return seedCandidates_inCluster;
0434 }
0435 
0436 void SETSeedFinder::estimateMomentum(const MuonRecHitContainer& validSet,
0437                                      CLHEP::Hep3Vector& momEstimate,
0438                                      int& charge) const {
0439   int firstMeasurement = -1;
0440   int lastMeasurement = -1;
0441 
0442   // don't use 2D measurements for momentum estimation
0443 
0444   //if( 4==allValidSets[iSet].front()->dimension() &&
0445   //(allValidSets[iSet].front()->isCSC() || allValidSets[iSet].front()->isDT())){
0446   //firstMeasurement = 0;
0447   //}
0448   //else{
0449   // which is the "first" hit (4D)?
0450   for (unsigned int iMeas = 0; iMeas < validSet.size(); ++iMeas) {
0451     if (4 == validSet[iMeas]->dimension() && (validSet[iMeas]->isCSC() || validSet[iMeas]->isDT())) {
0452       firstMeasurement = iMeas;
0453       break;
0454     }
0455   }
0456   //}
0457 
0458   std::vector<double> momentum_estimate;
0459   double pT = 0.;
0460   MuonTransientTrackingRecHit::ConstMuonRecHitPointer firstHit;
0461   MuonTransientTrackingRecHit::ConstMuonRecHitPointer secondHit;
0462   // which is the second hit?
0463   for (int loop = 0; loop < 2; ++loop) {  // it is actually not used; to be removed
0464     // this is the last measurement
0465     if (!loop) {  // this is what is used currently
0466       // 23.04.09 : it becomes a problem with introduction of ME42 chambers -
0467       // the initial pT parametrization is incorrect for them
0468       for (int iMeas = validSet.size() - 1; iMeas > -1; --iMeas) {
0469         if (4 == validSet[iMeas]->dimension() && (validSet[iMeas]->isCSC() || validSet[iMeas]->isDT()) &&
0470             // below is a fix saying "don't use ME4 chambers for initial pT estimation";
0471             // not using ME41 should not be a big loss too (and is more "symmetric" solution)
0472             fabs(validSet[iMeas]->globalPosition().z()) < 1000.) {
0473           lastMeasurement = iMeas;
0474           break;
0475         }
0476       }
0477     } else {
0478       // this is the second measurement
0479       for (unsigned int iMeas = 1; iMeas < validSet.size(); ++iMeas) {
0480         if (4 == validSet[iMeas]->dimension() && (validSet[iMeas]->isCSC() || validSet[iMeas]->isDT())) {
0481           lastMeasurement = iMeas;
0482           break;
0483         }
0484       }
0485     }
0486     // only 2D measurements (it should have been already abandoned)
0487     if (-1 == lastMeasurement && -1 == firstMeasurement) {
0488       firstMeasurement = 0;
0489       lastMeasurement = validSet.size() - 1;
0490     }
0491     // because of the ME42 above lastMeasurement could be -1
0492     else if (-1 == lastMeasurement) {
0493       lastMeasurement = firstMeasurement;
0494     } else if (-1 == firstMeasurement) {
0495       firstMeasurement = lastMeasurement;
0496     }
0497 
0498     firstHit = validSet[firstMeasurement];
0499     secondHit = validSet[lastMeasurement];
0500     if (firstHit->isRPC() && secondHit->isRPC()) {  // remove all RPCs from here?
0501       momentum_estimate.push_back(300.);
0502       momentum_estimate.push_back(300.);
0503     } else {
0504       if (firstHit->isRPC()) {
0505         firstHit = secondHit;
0506       } else if (secondHit->isRPC()) {
0507         secondHit = firstHit;
0508       }
0509       //---- estimate pT given two hits
0510       //std::cout<<"   hits for initial pT estimate: first -> dim = "<<firstHit->dimension()<<" pos = "<<firstHit->globalPosition()<<
0511       //" , second -> "<<" dim = "<<secondHit->dimension()<<" pos = "<<secondHit->globalPosition()<<std::endl;
0512       //---- pT throws exception if hits are MB4
0513       // (no coding for them - 2D hits in the outer station)
0514       if (2 == firstHit->dimension() && 2 == secondHit->dimension()) {
0515         momentum_estimate.push_back(999999999.);
0516         momentum_estimate.push_back(999999999.);
0517       } else {
0518         momentum_estimate = thePtExtractor->pT_extract(firstHit, secondHit);
0519       }
0520     }
0521     pT = fabs(momentum_estimate[0]);
0522     if (true || pT > 40.) {  //it is skipped; we have to look at least into number of hits in the chamber actually...
0523       // and then decide which segment to use
0524       // use the last measurement, otherwise use the second; this is to be investigated
0525       break;
0526     }
0527   }
0528 
0529   const float pT_min = 1.99;  // many hardcoded - remove them!
0530   if (pT > 3000.) {
0531     pT = 3000.;
0532   } else if (pT < pT_min) {
0533     if (pT > 0) {
0534       pT = pT_min;
0535     } else if (pT > (-1) * pT_min) {
0536       pT = (-1) * pT_min;
0537     } else if (pT < -3000.) {
0538       pT = -3000;
0539     }
0540   }
0541   //std::cout<<"  THE pT from the parametrization: "<<momentum_estimate[0]<<std::endl;
0542   // estimate the charge of the track candidate from the delta phi of two segments:
0543   //int charge      = dPhi > 0 ? 1 : -1; // what we want is: dphi < 0 => charge = -1
0544   charge = momentum_estimate[0] > 0 ? 1 : -1;
0545 
0546   // we have the pT - get the 3D momentum estimate as well
0547 
0548   // this is already final info:
0549   double xHit = validSet[firstMeasurement]->globalPosition().x();
0550   double yHit = validSet[firstMeasurement]->globalPosition().y();
0551   double rHit = TMath::Sqrt(pow(xHit, 2) + pow(yHit, 2));
0552 
0553   double thetaInner = validSet[firstMeasurement]->globalPosition().theta();
0554   // if some of the segments is missing r-phi measurement then we should
0555   // use only the 4D phi estimate (also use 4D eta estimate only)
0556   // the direction is not so important (it will be corrected)
0557 
0558   double rTrack = (pT / (0.3 * 3.8)) * 100.;  //times 100 for conversion to cm!
0559 
0560   double par = -1. * (2. / charge) * (TMath::ASin(rHit / (2 * rTrack)));
0561   double sinPar = TMath::Sin(par);
0562   double cosPar = TMath::Cos(par);
0563 
0564   // calculate phi at coordinate origin (0,0,0).
0565   double sinPhiH = 1. / (2. * charge * rTrack) * (xHit + ((sinPar) / (cosPar - 1.)) * yHit);
0566   double cosPhiH = -1. / (2. * charge * rTrack) * (((sinPar) / (1. - cosPar)) * xHit + yHit);
0567 
0568   // finally set the return vector
0569 
0570   // try out the reco info:
0571   // should used into to theta directly here (rather than tan(atan2(...)))
0572   momEstimate = CLHEP::Hep3Vector(pT * cosPhiH, pT * sinPhiH, pT / TMath::Tan(thetaInner));
0573   //Hep3Vector momEstimate(6.97961,      5.89732,     -50.0855);
0574   const float minMomenum = 5.;  //hardcoded - remove it! same in SETFilter
0575   if (momEstimate.mag() < minMomenum) {
0576     int sign = (pT < 0.) ? -1 : 1;
0577     pT = sign * (fabs(pT) + 1);
0578     CLHEP::Hep3Vector momEstimate2(pT * cosPhiH, pT * sinPhiH, pT / TMath::Tan(thetaInner));
0579     momEstimate = momEstimate2;
0580     if (momEstimate.mag() < minMomenum) {
0581       pT = sign * (fabs(pT) + 1);
0582       CLHEP::Hep3Vector momEstimate3(pT * cosPhiH, pT * sinPhiH, pT / TMath::Tan(thetaInner));
0583       momEstimate = momEstimate3;
0584       if (momEstimate.mag() < minMomenum) {
0585         pT = sign * (fabs(pT) + 1);
0586         CLHEP::Hep3Vector momEstimate4(pT * cosPhiH, pT * sinPhiH, pT / TMath::Tan(thetaInner));
0587         momEstimate = momEstimate4;
0588       }
0589     }
0590   }
0591 }
0592 
0593 TrajectorySeed SETSeedFinder::makeSeed(const TrajectoryStateOnSurface& firstTSOS,
0594                                        const TransientTrackingRecHit::ConstRecHitContainer& hits) const {
0595   edm::OwnVector<TrackingRecHit> recHitsContainer;
0596   for (unsigned int iHit = 0; iHit < hits.size(); ++iHit) {
0597     recHitsContainer.push_back(hits.at(iHit)->hit()->clone());
0598   }
0599   PropagationDirection dir = oppositeToMomentum;
0600   if (useSegmentsInTrajectory) {
0601     dir = alongMomentum;  // why forward (for rechits) later?
0602   }
0603 
0604   PTrajectoryStateOnDet const& seedTSOS =
0605       trajectoryStateTransform::persistentState(firstTSOS, hits.at(0)->geographicalId().rawId());
0606   TrajectorySeed seed(seedTSOS, recHitsContainer, dir);
0607 
0608   //MuonPatternRecoDumper debug;
0609   //std::cout<<" firstTSOS = "<<debug.dumpTSOS(firstTSOS)<<std::endl;
0610   //std::cout<<" iTraj = ???"<<" hits = "<<range.second-range.first<<std::endl;
0611   //std::cout<<" nhits = "<<hits.size()<<std::endl;
0612   //for(unsigned int iRH=0;iRH<hits.size();++iRH){
0613   //std::cout<<" RH = "<<iRH+1<<" globPos = "<<hits.at(iRH)->globalPosition()<<std::endl;
0614   //}
0615   return seed;
0616 }