Back to home page

Project CMSSW displayed by LXR
 
 

    

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

0001 /*
0002  *  See header file for a description of this class.
0003  *
0004  *
0005  *  \author Haiyun.Teng - Peking University
0006  *
0007  */
0008 
0009 #include "RecoMuon/MuonSeedGenerator/src/RPCSeedPattern.h"
0010 #include <RecoMuon/TrackingTools/interface/MuonPatternRecoDumper.h>
0011 #include <MagneticField/Engine/interface/MagneticField.h>
0012 #include <TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h>
0013 #include <TrackingTools/DetLayers/interface/DetLayer.h>
0014 #include <DataFormats/TrajectoryState/interface/PTrajectoryStateOnDet.h>
0015 #include <DataFormats/Common/interface/OwnVector.h>
0016 #include <DataFormats/MuonDetId/interface/RPCDetId.h>
0017 #include <FWCore/Framework/interface/ESHandle.h>
0018 #include <FWCore/MessageLogger/interface/MessageLogger.h>
0019 #include <Geometry/CommonDetUnit/interface/GeomDet.h>
0020 #include <Geometry/RPCGeometry/interface/RPCChamber.h>
0021 #include <Geometry/RPCGeometry/interface/RPCGeometry.h>
0022 
0023 #include "gsl/gsl_statistics.h"
0024 #include "TH1F.h"
0025 #include <cmath>
0026 
0027 using namespace std;
0028 using namespace edm;
0029 
0030 RPCSeedPattern::RPCSeedPattern() {
0031   isPatternChecked = false;
0032   isConfigured = false;
0033   MagnecticFieldThreshold = 0.5;
0034 }
0035 
0036 RPCSeedPattern::~RPCSeedPattern() {}
0037 
0038 void RPCSeedPattern::configure(const edm::ParameterSet& iConfig) {
0039   MaxRSD = iConfig.getParameter<double>("MaxRSD");
0040   deltaRThreshold = iConfig.getParameter<double>("deltaRThreshold");
0041   AlgorithmType = iConfig.getParameter<unsigned int>("AlgorithmType");
0042   autoAlgorithmChoose = iConfig.getParameter<bool>("autoAlgorithmChoose");
0043   ZError = iConfig.getParameter<double>("ZError");
0044   MinDeltaPhi = iConfig.getParameter<double>("MinDeltaPhi");
0045   MagnecticFieldThreshold = iConfig.getParameter<double>("MagnecticFieldThreshold");
0046   stepLength = iConfig.getParameter<double>("stepLength");
0047   sampleCount = iConfig.getParameter<unsigned int>("sampleCount");
0048   isConfigured = true;
0049 }
0050 
0051 RPCSeedPattern::weightedTrajectorySeed RPCSeedPattern::seed(const MagneticField& Field,
0052                                                             const RPCGeometry& rpcGeom,
0053                                                             int& isGoodSeed) {
0054   if (isConfigured == false) {
0055     cout << "Configuration not set yet" << endl;
0056     return createFakeSeed(isGoodSeed, Field);
0057   }
0058 
0059   // Check recHit number, if fail we return a fake seed and set pattern to "wrong"
0060   unsigned int NumberofHitsinSeed = nrhit();
0061   if (NumberofHitsinSeed < 3) {
0062     return createFakeSeed(isGoodSeed, Field);
0063   }
0064   // If only three recHits, we don't have other choice
0065   if (NumberofHitsinSeed == 3)
0066     ThreePointsAlgorithm();
0067 
0068   if (NumberofHitsinSeed > 3) {
0069     if (autoAlgorithmChoose == false) {
0070       cout << "computePtWithmorerecHits" << endl;
0071       if (AlgorithmType == 0)
0072         ThreePointsAlgorithm();
0073       if (AlgorithmType == 1)
0074         MiddlePointsAlgorithm();
0075       if (AlgorithmType == 2)
0076         SegmentAlgorithm();
0077       if (AlgorithmType == 3) {
0078         if (checkSegment()) {
0079           SegmentAlgorithmSpecial(Field);
0080         } else {
0081           cout << "Not enough recHits for Special Segment Algorithm" << endl;
0082           return createFakeSeed(isGoodSeed, Field);
0083         }
0084       }
0085     } else {
0086       if (checkSegment()) {
0087         AlgorithmType = 3;
0088         SegmentAlgorithmSpecial(Field);
0089       } else {
0090         AlgorithmType = 2;
0091         SegmentAlgorithm();
0092       }
0093     }
0094   }
0095 
0096   // Check the pattern
0097   if (isPatternChecked == false) {
0098     if (AlgorithmType != 3) {
0099       checkSimplePattern(Field);
0100     } else {
0101       checkSegmentAlgorithmSpecial(Field, rpcGeom);
0102     }
0103   }
0104 
0105   return createSeed(isGoodSeed, Field);
0106 }
0107 
0108 void RPCSeedPattern::ThreePointsAlgorithm() {
0109   cout << "computePtWith3recHits" << endl;
0110   unsigned int NumberofHitsinSeed = nrhit();
0111   // Check recHit number, if fail we set the pattern to "wrong"
0112   if (NumberofHitsinSeed < 3) {
0113     isPatternChecked = true;
0114     isGoodPattern = -1;
0115     return;
0116   }
0117   // Choose every 3 recHits to form a part
0118   unsigned int NumberofPart = NumberofHitsinSeed * (NumberofHitsinSeed - 1) * (NumberofHitsinSeed - 2) / (3 * 2);
0119   ;
0120   double* pt = new double[NumberofPart];
0121   double* pt_err = new double[NumberofPart];
0122   // Loop for each three-recHits part
0123   ConstMuonRecHitPointer precHit[3];
0124   unsigned int n = 0;
0125   unsigned int NumberofStraight = 0;
0126   for (unsigned int i = 0; i < (NumberofHitsinSeed - 2); i++)
0127     for (unsigned int j = (i + 1); j < (NumberofHitsinSeed - 1); j++)
0128       for (unsigned int k = (j + 1); k < NumberofHitsinSeed; k++) {
0129         precHit[0] = theRecHits[i];
0130         precHit[1] = theRecHits[j];
0131         precHit[2] = theRecHits[k];
0132         bool checkStraight = checkStraightwithThreerecHits(precHit, MinDeltaPhi);
0133         if (!checkStraight) {
0134           GlobalVector Center_temp = computePtwithThreerecHits(pt[n], pt_err[n], precHit);
0135           // For simple pattern
0136           Center += Center_temp;
0137         } else {
0138           // For simple pattern
0139           NumberofStraight++;
0140           pt[n] = upper_limit_pt;
0141           pt_err[n] = 0;
0142         }
0143         n++;
0144       }
0145   // For simple pattern, only one general parameter for pattern
0146   if (NumberofStraight == NumberofPart) {
0147     isStraight = true;
0148     meanRadius = -1;
0149   } else {
0150     isStraight = false;
0151     Center /= (NumberofPart - NumberofStraight);
0152     double meanR = 0.;
0153     for (ConstMuonRecHitContainer::const_iterator iter = theRecHits.begin(); iter != theRecHits.end(); iter++)
0154       meanR += getDistance(*iter, Center);
0155     meanR /= NumberofHitsinSeed;
0156     meanRadius = meanR;
0157   }
0158 
0159   // Unset the pattern estimation signa
0160   isPatternChecked = false;
0161 
0162   //double ptmean0 = 0;
0163   //double sptmean0 = 0;
0164   //computeBestPt(pt, pt_err, ptmean0, sptmean0, (NumberofPart - NumberofStraight));
0165 
0166   delete[] pt;
0167   delete[] pt_err;
0168 }
0169 
0170 void RPCSeedPattern::MiddlePointsAlgorithm() {
0171   cout << "Using middle points algorithm" << endl;
0172   unsigned int NumberofHitsinSeed = nrhit();
0173   // Check recHit number, if fail we set the pattern to "wrong"
0174   if (NumberofHitsinSeed < 4) {
0175     isPatternChecked = true;
0176     isGoodPattern = -1;
0177     return;
0178   }
0179   double* X = new double[NumberofHitsinSeed];
0180   double* Y = new double[NumberofHitsinSeed];
0181   unsigned int n = 0;
0182   for (ConstMuonRecHitContainer::const_iterator iter = theRecHits.begin(); iter != theRecHits.end(); iter++) {
0183     X[n] = (*iter)->globalPosition().x();
0184     Y[n] = (*iter)->globalPosition().y();
0185     cout << "X[" << n << "] = " << X[n] << ", Y[" << n << "]= " << Y[n] << endl;
0186     n++;
0187   }
0188   unsigned int NumberofPoints = NumberofHitsinSeed;
0189   while (NumberofPoints > 3) {
0190     for (unsigned int i = 0; i <= (NumberofPoints - 2); i++) {
0191       X[i] = (X[i] + X[i + 1]) / 2;
0192       Y[i] = (Y[i] + Y[i + 1]) / 2;
0193     }
0194     NumberofPoints--;
0195   }
0196   double x[3], y[3];
0197   for (unsigned int i = 0; i < 3; i++) {
0198     x[i] = X[i];
0199     y[i] = Y[i];
0200   }
0201   double pt = 0;
0202   double pt_err = 0;
0203   bool checkStraight = checkStraightwithThreerecHits(x, y, MinDeltaPhi);
0204   if (!checkStraight) {
0205     GlobalVector Center_temp = computePtWithThreerecHits(pt, pt_err, x, y);
0206     double meanR = 0.;
0207     for (ConstMuonRecHitContainer::const_iterator iter = theRecHits.begin(); iter != theRecHits.end(); iter++)
0208       meanR += getDistance(*iter, Center_temp);
0209     meanR /= NumberofHitsinSeed;
0210     // For simple pattern
0211     isStraight = false;
0212     Center = Center_temp;
0213     meanRadius = meanR;
0214   } else {
0215     // For simple pattern
0216     isStraight = true;
0217     meanRadius = -1;
0218   }
0219 
0220   // Unset the pattern estimation signa
0221   isPatternChecked = false;
0222 
0223   delete[] X;
0224   delete[] Y;
0225 }
0226 
0227 void RPCSeedPattern::SegmentAlgorithm() {
0228   cout << "Using segments algorithm" << endl;
0229   unsigned int NumberofHitsinSeed = nrhit();
0230   // Check recHit number, if fail we set the pattern to "wrong"
0231   if (NumberofHitsinSeed < 4) {
0232     isPatternChecked = true;
0233     isGoodPattern = -1;
0234     return;
0235   }
0236 
0237   RPCSegment* Segment;
0238   unsigned int NumberofSegment = NumberofHitsinSeed - 2;
0239   Segment = new RPCSegment[NumberofSegment];
0240   unsigned int n = 0;
0241   for (ConstMuonRecHitContainer::const_iterator iter = theRecHits.begin(); iter != (theRecHits.end() - 2); iter++) {
0242     Segment[n].first = (*iter);
0243     Segment[n].second = (*(iter + 2));
0244     n++;
0245   }
0246   unsigned int NumberofStraight = 0;
0247   for (unsigned int i = 0; i < NumberofSegment - 1; i++) {
0248     bool checkStraight = checkStraightwithSegment(Segment[i], Segment[i + 1], MinDeltaPhi);
0249     if (checkStraight == true) {
0250       // For simple patterm
0251       NumberofStraight++;
0252     } else {
0253       GlobalVector Center_temp = computePtwithSegment(Segment[i], Segment[i + 1]);
0254       // For simple patterm
0255       Center += Center_temp;
0256     }
0257   }
0258   // For simple pattern, only one general parameter for pattern
0259   if ((NumberofSegment - 1 - NumberofStraight) > 0) {
0260     isStraight = false;
0261     Center /= (NumberofSegment - 1 - NumberofStraight);
0262     double meanR = 0.;
0263     for (ConstMuonRecHitContainer::const_iterator iter = theRecHits.begin(); iter != theRecHits.end(); iter++)
0264       meanR += getDistance(*iter, Center);
0265     meanR /= NumberofHitsinSeed;
0266     meanRadius = meanR;
0267   } else {
0268     isStraight = true;
0269     meanRadius = -1;
0270   }
0271 
0272   // Unset the pattern estimation signal
0273   isPatternChecked = false;
0274 
0275   delete[] Segment;
0276 }
0277 
0278 void RPCSeedPattern::SegmentAlgorithmSpecial(const MagneticField& Field) {
0279   //unsigned int NumberofHitsinSeed = nrhit();
0280   if (!checkSegment()) {
0281     isPatternChecked = true;
0282     isGoodPattern = -1;
0283     return;
0284   }
0285 
0286   // Get magnetice field sampling information, recHit's position is not the border of Chamber and Iron
0287   for (ConstMuonRecHitContainer::const_iterator iter = theRecHits.begin(); iter != (theRecHits.end() - 1); iter++) {
0288     GlobalPoint gpFirst = (*iter)->globalPosition();
0289     GlobalPoint gpLast = (*(iter + 1))->globalPosition();
0290     GlobalPoint* gp = new GlobalPoint[sampleCount];
0291     double dx = (gpLast.x() - gpFirst.x()) / (sampleCount + 1);
0292     double dy = (gpLast.y() - gpFirst.y()) / (sampleCount + 1);
0293     double dz = (gpLast.z() - gpFirst.z()) / (sampleCount + 1);
0294     for (unsigned int index = 0; index < sampleCount; index++) {
0295       gp[index] = GlobalPoint(
0296           (gpFirst.x() + dx * (index + 1)), (gpFirst.y() + dy * (index + 1)), (gpFirst.z() + dz * (index + 1)));
0297       GlobalVector MagneticVec_temp = Field.inTesla(gp[index]);
0298       cout << "Sampling magnetic field : " << MagneticVec_temp << endl;
0299       //BValue.push_back(MagneticVec_temp);
0300     }
0301     delete[] gp;
0302   }
0303 
0304   // form two segments
0305   ConstMuonRecHitContainer::const_iterator iter = theRecHits.begin();
0306   for (unsigned int n = 0; n <= 1; n++) {
0307     SegmentRB[n].first = (*iter);
0308     cout << "SegmentRB " << n << " recHit: " << (*iter)->globalPosition() << endl;
0309     iter++;
0310     SegmentRB[n].second = (*iter);
0311     cout << "SegmentRB " << n << " recHit: " << (*iter)->globalPosition() << endl;
0312     iter++;
0313   }
0314   GlobalVector segvec1 = (SegmentRB[0].second)->globalPosition() - (SegmentRB[0].first)->globalPosition();
0315   GlobalVector segvec2 = (SegmentRB[1].second)->globalPosition() - (SegmentRB[1].first)->globalPosition();
0316 
0317   // extrapolate the segment to find the Iron border which magnetic field is at large value
0318   entryPosition = (SegmentRB[0].second)->globalPosition();
0319   leavePosition = (SegmentRB[1].first)->globalPosition();
0320   while (fabs(Field.inTesla(entryPosition).z()) < MagnecticFieldThreshold) {
0321     cout << "Entry position is : " << entryPosition << ", and stepping into next point" << endl;
0322     entryPosition += segvec1.unit() * stepLength;
0323   }
0324   // Loop back for more accurate by stepLength/10
0325   while (fabs(Field.inTesla(entryPosition).z()) >= MagnecticFieldThreshold) {
0326     cout << "Entry position is : " << entryPosition << ", and stepping back into next point" << endl;
0327     entryPosition -= segvec1.unit() * stepLength / 10;
0328   }
0329   entryPosition += 0.5 * segvec1.unit() * stepLength / 10;
0330   cout << "Final entry position is : " << entryPosition << endl;
0331 
0332   while (fabs(Field.inTesla(leavePosition).z()) < MagnecticFieldThreshold) {
0333     cout << "Leave position is : " << leavePosition << ", and stepping into next point" << endl;
0334     leavePosition -= segvec2.unit() * stepLength;
0335   }
0336   // Loop back for more accurate by stepLength/10
0337   while (fabs(Field.inTesla(leavePosition).z()) >= MagnecticFieldThreshold) {
0338     cout << "Leave position is : " << leavePosition << ", and stepping back into next point" << endl;
0339     leavePosition += segvec2.unit() * stepLength / 10;
0340   }
0341   leavePosition -= 0.5 * segvec2.unit() * stepLength / 10;
0342   cout << "Final leave position is : " << leavePosition << endl;
0343 
0344   // Sampling magnetic field in Iron region
0345   GlobalPoint* gp = new GlobalPoint[sampleCount];
0346   double dx = (leavePosition.x() - entryPosition.x()) / (sampleCount + 1);
0347   double dy = (leavePosition.y() - entryPosition.y()) / (sampleCount + 1);
0348   double dz = (leavePosition.z() - entryPosition.z()) / (sampleCount + 1);
0349   std::vector<GlobalVector> BValue;
0350   BValue.clear();
0351   for (unsigned int index = 0; index < sampleCount; index++) {
0352     gp[index] = GlobalPoint((entryPosition.x() + dx * (index + 1)),
0353                             (entryPosition.y() + dy * (index + 1)),
0354                             (entryPosition.z() + dz * (index + 1)));
0355     GlobalVector MagneticVec_temp = Field.inTesla(gp[index]);
0356     cout << "Sampling magnetic field : " << MagneticVec_temp << endl;
0357     BValue.push_back(MagneticVec_temp);
0358   }
0359   delete[] gp;
0360   GlobalVector meanB2(0, 0, 0);
0361   for (std::vector<GlobalVector>::const_iterator BIter = BValue.begin(); BIter != BValue.end(); BIter++)
0362     meanB2 += (*BIter);
0363   meanB2 /= BValue.size();
0364   cout << "Mean B field is " << meanB2 << endl;
0365   meanMagneticField2 = meanB2;
0366 
0367   double meanBz2 = meanB2.z();
0368   double deltaBz2 = 0.;
0369   for (std::vector<GlobalVector>::const_iterator BIter = BValue.begin(); BIter != BValue.end(); BIter++)
0370     deltaBz2 += (BIter->z() - meanBz2) * (BIter->z() - meanBz2);
0371   deltaBz2 /= BValue.size();
0372   deltaBz2 = sqrt(deltaBz2);
0373   cout << "delta Bz is " << deltaBz2 << endl;
0374 
0375   // Distance of the initial 3 segment
0376   S = 0;
0377   bool checkStraight = checkStraightwithSegment(SegmentRB[0], SegmentRB[1], MinDeltaPhi);
0378   if (checkStraight == true) {
0379     // Just for complex pattern
0380     isStraight2 = checkStraight;
0381     Center2 = GlobalVector(0, 0, 0);
0382     meanRadius2 = -1;
0383     GlobalVector MomentumVec = (SegmentRB[1].second)->globalPosition() - (SegmentRB[0].first)->globalPosition();
0384     S += MomentumVec.perp();
0385     lastPhi = MomentumVec.phi().value();
0386   } else {
0387     GlobalVector seg1 = entryPosition - (SegmentRB[0].first)->globalPosition();
0388     S += seg1.perp();
0389     GlobalVector seg2 = (SegmentRB[1].second)->globalPosition() - leavePosition;
0390     S += seg2.perp();
0391     GlobalVector vecZ(0, 0, 1);
0392     GlobalVector gvec1 = seg1.cross(vecZ);
0393     GlobalVector gvec2 = seg2.cross(vecZ);
0394     double A1 = gvec1.x();
0395     double B1 = gvec1.y();
0396     double A2 = gvec2.x();
0397     double B2 = gvec2.y();
0398     double X1 = entryPosition.x();
0399     double Y1 = entryPosition.y();
0400     double X2 = leavePosition.x();
0401     double Y2 = leavePosition.y();
0402     double XO = (A1 * A2 * (Y2 - Y1) + A2 * B1 * X1 - A1 * B2 * X2) / (A2 * B1 - A1 * B2);
0403     double YO = (B1 * B2 * (X2 - X1) + B2 * A1 * Y1 - B1 * A2 * Y2) / (B2 * A1 - B1 * A2);
0404     GlobalVector Center_temp(XO, YO, 0);
0405     // Just for complex pattern
0406     isStraight2 = checkStraight;
0407     Center2 = Center_temp;
0408 
0409     cout << "entryPosition: " << entryPosition << endl;
0410     cout << "leavePosition: " << leavePosition << endl;
0411     cout << "Center2 is : " << Center_temp << endl;
0412 
0413     double R1 = GlobalVector((entryPosition.x() - Center_temp.x()),
0414                              (entryPosition.y() - Center_temp.y()),
0415                              (entryPosition.z() - Center_temp.z()))
0416                     .perp();
0417     double R2 = GlobalVector((leavePosition.x() - Center_temp.x()),
0418                              (leavePosition.y() - Center_temp.y()),
0419                              (leavePosition.z() - Center_temp.z()))
0420                     .perp();
0421     double meanR = (R1 + R2) / 2;
0422     double deltaR = sqrt(((R1 - meanR) * (R1 - meanR) + (R2 - meanR) * (R2 - meanR)) / 2);
0423     meanRadius2 = meanR;
0424     cout << "R1 is " << R1 << ", R2 is " << R2 << endl;
0425     cout << "Mean radius is " << meanR << endl;
0426     cout << "Delta R is " << deltaR << endl;
0427     double deltaPhi =
0428         fabs(((leavePosition - GlobalPoint(XO, YO, 0)).phi() - (entryPosition - GlobalPoint(XO, YO, 0)).phi()).value());
0429     S += meanR * deltaPhi;
0430     lastPhi = seg2.phi().value();
0431   }
0432 
0433   // Unset the pattern estimation signa
0434   isPatternChecked = false;
0435 }
0436 
0437 bool RPCSeedPattern::checkSegment() const {
0438   bool isFit = true;
0439   unsigned int count = 0;
0440   // first 4 recHits should be located in RB1 and RB2
0441   for (ConstMuonRecHitContainer::const_iterator iter = theRecHits.begin(); iter != theRecHits.end(); iter++) {
0442     count++;
0443     const GeomDet* Detector = (*iter)->det();
0444     if (dynamic_cast<const RPCChamber*>(Detector) != nullptr) {
0445       const RPCChamber* RPCCh = dynamic_cast<const RPCChamber*>(Detector);
0446       RPCDetId RPCId = RPCCh->id();
0447       int Region = RPCId.region();
0448       unsigned int Station = RPCId.station();
0449       //int Layer = RPCId.layer();
0450       if (count <= 4) {
0451         if (Region != 0)
0452           isFit = false;
0453         if (Station > 2)
0454           isFit = false;
0455       }
0456     }
0457   }
0458   // more than 4 recHits for pattern building
0459   if (count <= 4)
0460     isFit = false;
0461   cout << "Check for segment fit: " << isFit << endl;
0462   return isFit;
0463 }
0464 
0465 MuonTransientTrackingRecHit::ConstMuonRecHitPointer RPCSeedPattern::FirstRecHit() const { return theRecHits.front(); }
0466 
0467 MuonTransientTrackingRecHit::ConstMuonRecHitPointer RPCSeedPattern::BestRefRecHit() const {
0468   ConstMuonRecHitPointer best;
0469   int index = 0;
0470   // Use the last one for recHit on last layer has minmum delta Z for barrel or delta R for endcap while calculating the momentum
0471   // But for Algorithm 3 we use the 4th recHit on the 2nd segment for more accurate
0472   for (ConstMuonRecHitContainer::const_iterator iter = theRecHits.begin(); iter != theRecHits.end(); iter++) {
0473     if (AlgorithmType != 3)
0474       best = (*iter);
0475     else if (index < 4)
0476       best = (*iter);
0477     index++;
0478   }
0479   return best;
0480 }
0481 
0482 double RPCSeedPattern::getDistance(const ConstMuonRecHitPointer& precHit, const GlobalVector& Center) const {
0483   return sqrt((precHit->globalPosition().x() - Center.x()) * (precHit->globalPosition().x() - Center.x()) +
0484               (precHit->globalPosition().y() - Center.y()) * (precHit->globalPosition().y() - Center.y()));
0485 }
0486 
0487 bool RPCSeedPattern::checkStraightwithThreerecHits(ConstMuonRecHitPointer (&precHit)[3], double MinDeltaPhi) const {
0488   GlobalVector segvec1 = precHit[1]->globalPosition() - precHit[0]->globalPosition();
0489   GlobalVector segvec2 = precHit[2]->globalPosition() - precHit[1]->globalPosition();
0490   double dPhi = (segvec2.phi() - segvec1.phi()).value();
0491   if (fabs(dPhi) > MinDeltaPhi) {
0492     cout << "Part is estimate to be not straight" << endl;
0493     return false;
0494   } else {
0495     cout << "Part is estimate to be straight" << endl;
0496     return true;
0497   }
0498 }
0499 
0500 GlobalVector RPCSeedPattern::computePtwithThreerecHits(double& pt,
0501                                                        double& pt_err,
0502                                                        ConstMuonRecHitPointer (&precHit)[3]) const {
0503   double x[3], y[3];
0504   x[0] = precHit[0]->globalPosition().x();
0505   y[0] = precHit[0]->globalPosition().y();
0506   x[1] = precHit[1]->globalPosition().x();
0507   y[1] = precHit[1]->globalPosition().y();
0508   x[2] = precHit[2]->globalPosition().x();
0509   y[2] = precHit[2]->globalPosition().y();
0510   double A = (y[2] - y[1]) / (x[2] - x[1]) - (y[1] - y[0]) / (x[1] - x[0]);
0511   double TYO = (x[2] - x[0]) / A + (y[2] * y[2] - y[1] * y[1]) / ((x[2] - x[1]) * A) -
0512                (y[1] * y[1] - y[0] * y[0]) / ((x[1] - x[0]) * A);
0513   double TXO = (x[2] + x[1]) + (y[2] * y[2] - y[1] * y[1]) / (x[2] - x[1]) - TYO * (y[2] - y[1]) / (x[2] - x[1]);
0514   double XO = 0.5 * TXO;
0515   double YO = 0.5 * TYO;
0516   double R2 = (x[0] - XO) * (x[0] - XO) + (y[0] - YO) * (y[0] - YO);
0517   cout << "R2 is " << R2 << endl;
0518   // How this algorithm get the pt without magnetic field??
0519   pt = 0.01 * sqrt(R2) * 2 * 0.3;
0520   cout << "pt is " << pt << endl;
0521   GlobalVector Center(XO, YO, 0);
0522   return Center;
0523 }
0524 
0525 bool RPCSeedPattern::checkStraightwithSegment(const RPCSegment& Segment1,
0526                                               const RPCSegment& Segment2,
0527                                               double MinDeltaPhi) const {
0528   GlobalVector segvec1 = (Segment1.second)->globalPosition() - (Segment1.first)->globalPosition();
0529   GlobalVector segvec2 = (Segment2.second)->globalPosition() - (Segment2.first)->globalPosition();
0530   GlobalVector segvec3 = (Segment2.first)->globalPosition() - (Segment1.first)->globalPosition();
0531   // compare segvec 1&2 for paralle, 1&3 for straight
0532   double dPhi1 = (segvec2.phi() - segvec1.phi()).value();
0533   double dPhi2 = (segvec3.phi() - segvec1.phi()).value();
0534   cout << "Checking straight with 2 segments. dPhi1: " << dPhi1 << ", dPhi2: " << dPhi2 << endl;
0535   cout << "Checking straight with 2 segments. dPhi1 in degree: " << dPhi1 * 180 / 3.1415926
0536        << ", dPhi2 in degree: " << dPhi2 * 180 / 3.1415926 << endl;
0537   if (fabs(dPhi1) > MinDeltaPhi || fabs(dPhi2) > MinDeltaPhi) {
0538     cout << "Segment is estimate to be not straight" << endl;
0539     return false;
0540   } else {
0541     cout << "Segment is estimate to be straight" << endl;
0542     return true;
0543   }
0544 }
0545 
0546 GlobalVector RPCSeedPattern::computePtwithSegment(const RPCSegment& Segment1, const RPCSegment& Segment2) const {
0547   GlobalVector segvec1 = (Segment1.second)->globalPosition() - (Segment1.first)->globalPosition();
0548   GlobalVector segvec2 = (Segment2.second)->globalPosition() - (Segment2.first)->globalPosition();
0549   GlobalPoint Point1(((Segment1.second)->globalPosition().x() + (Segment1.first)->globalPosition().x()) / 2,
0550                      ((Segment1.second)->globalPosition().y() + (Segment1.first)->globalPosition().y()) / 2,
0551                      ((Segment1.second)->globalPosition().z() + (Segment1.first)->globalPosition().z()) / 2);
0552   GlobalPoint Point2(((Segment2.second)->globalPosition().x() + (Segment2.first)->globalPosition().x()) / 2,
0553                      ((Segment2.second)->globalPosition().y() + (Segment2.first)->globalPosition().y()) / 2,
0554                      ((Segment2.second)->globalPosition().z() + (Segment2.first)->globalPosition().z()) / 2);
0555   GlobalVector vecZ(0, 0, 1);
0556   GlobalVector gvec1 = segvec1.cross(vecZ);
0557   GlobalVector gvec2 = segvec2.cross(vecZ);
0558   double A1 = gvec1.x();
0559   double B1 = gvec1.y();
0560   double A2 = gvec2.x();
0561   double B2 = gvec2.y();
0562   double X1 = Point1.x();
0563   double Y1 = Point1.y();
0564   double X2 = Point2.x();
0565   double Y2 = Point2.y();
0566   double XO = (A1 * A2 * (Y2 - Y1) + A2 * B1 * X1 - A1 * B2 * X2) / (A2 * B1 - A1 * B2);
0567   double YO = (B1 * B2 * (X2 - X1) + B2 * A1 * Y1 - B1 * A2 * Y2) / (B2 * A1 - B1 * A2);
0568   GlobalVector Center(XO, YO, 0);
0569   return Center;
0570 }
0571 
0572 bool RPCSeedPattern::checkStraightwithThreerecHits(double (&x)[3], double (&y)[3], double MinDeltaPhi) const {
0573   GlobalVector segvec1((x[1] - x[0]), (y[1] - y[0]), 0);
0574   GlobalVector segvec2((x[2] - x[1]), (y[2] - y[1]), 0);
0575   double dPhi = (segvec2.phi() - segvec1.phi()).value();
0576   if (fabs(dPhi) > MinDeltaPhi) {
0577     cout << "Part is estimate to be not straight" << endl;
0578     return false;
0579   } else {
0580     cout << "Part is estimate to be straight" << endl;
0581     return true;
0582   }
0583 }
0584 
0585 GlobalVector RPCSeedPattern::computePtWithThreerecHits(double& pt,
0586                                                        double& pt_err,
0587                                                        double (&x)[3],
0588                                                        double (&y)[3]) const {
0589   double A = (y[2] - y[1]) / (x[2] - x[1]) - (y[1] - y[0]) / (x[1] - x[0]);
0590   double TYO = (x[2] - x[0]) / A + (y[2] * y[2] - y[1] * y[1]) / ((x[2] - x[1]) * A) -
0591                (y[1] * y[1] - y[0] * y[0]) / ((x[1] - x[0]) * A);
0592   double TXO = (x[2] + x[1]) + (y[2] * y[2] - y[1] * y[1]) / (x[2] - x[1]) - TYO * (y[2] - y[1]) / (x[2] - x[1]);
0593   double XO = 0.5 * TXO;
0594   double YO = 0.5 * TYO;
0595   double R2 = (x[0] - XO) * (x[0] - XO) + (y[0] - YO) * (y[0] - YO);
0596   cout << "R2 is " << R2 << endl;
0597   // How this algorithm get the pt without magnetic field??
0598   pt = 0.01 * sqrt(R2) * 2 * 0.3;
0599   cout << "pt is " << pt << endl;
0600   GlobalVector Center(XO, YO, 0);
0601   return Center;
0602 }
0603 
0604 void RPCSeedPattern::checkSimplePattern(const MagneticField& Field) {
0605   if (isPatternChecked == true)
0606     return;
0607 
0608   unsigned int NumberofHitsinSeed = nrhit();
0609 
0610   // Print the recHit's position
0611   for (ConstMuonRecHitContainer::const_iterator iter = theRecHits.begin(); iter != theRecHits.end(); iter++)
0612     cout << "Position of recHit is: " << (*iter)->globalPosition() << endl;
0613 
0614   // Get magnetice field information
0615   std::vector<double> BzValue;
0616   for (ConstMuonRecHitContainer::const_iterator iter = theRecHits.begin(); iter != (theRecHits.end() - 1); iter++) {
0617     GlobalPoint gpFirst = (*iter)->globalPosition();
0618     GlobalPoint gpLast = (*(iter + 1))->globalPosition();
0619     GlobalPoint* gp = new GlobalPoint[sampleCount];
0620     double dx = (gpLast.x() - gpFirst.x()) / (sampleCount + 1);
0621     double dy = (gpLast.y() - gpFirst.y()) / (sampleCount + 1);
0622     double dz = (gpLast.z() - gpFirst.z()) / (sampleCount + 1);
0623     for (unsigned int index = 0; index < sampleCount; index++) {
0624       gp[index] = GlobalPoint(
0625           (gpFirst.x() + dx * (index + 1)), (gpFirst.y() + dy * (index + 1)), (gpFirst.z() + dz * (index + 1)));
0626       GlobalVector MagneticVec_temp = Field.inTesla(gp[index]);
0627       cout << "Sampling magnetic field : " << MagneticVec_temp << endl;
0628       BzValue.push_back(MagneticVec_temp.z());
0629     }
0630     delete[] gp;
0631   }
0632   meanBz = 0.;
0633   for (unsigned int index = 0; index < BzValue.size(); index++)
0634     meanBz += BzValue[index];
0635   meanBz /= BzValue.size();
0636   cout << "Mean Bz is " << meanBz << endl;
0637   deltaBz = 0.;
0638   for (unsigned int index = 0; index < BzValue.size(); index++)
0639     deltaBz += (BzValue[index] - meanBz) * (BzValue[index] - meanBz);
0640   deltaBz /= BzValue.size();
0641   deltaBz = sqrt(deltaBz);
0642   cout << "delata Bz is " << deltaBz << endl;
0643 
0644   // Set isGoodPattern to default true and check the failure
0645   isGoodPattern = 1;
0646 
0647   // Check the Z direction
0648   if (fabs((*(theRecHits.end() - 1))->globalPosition().z() - (*(theRecHits.begin()))->globalPosition().z()) > ZError) {
0649     if (((*(theRecHits.end() - 1))->globalPosition().z() - (*(theRecHits.begin()))->globalPosition().z()) > ZError)
0650       isParralZ = 1;
0651     else
0652       isParralZ = -1;
0653   } else
0654     isParralZ = 0;
0655 
0656   cout << " Check isParralZ is :" << isParralZ << endl;
0657   for (ConstMuonRecHitContainer::const_iterator iter = theRecHits.begin(); iter != (theRecHits.end() - 1); iter++) {
0658     if (isParralZ == 0) {
0659       if (fabs((*(iter + 1))->globalPosition().z() - (*iter)->globalPosition().z()) > ZError) {
0660         cout << "Pattern find error in Z direction: wrong perpendicular direction" << endl;
0661         isGoodPattern = 0;
0662       }
0663     } else {
0664       if ((int)(((*(iter + 1))->globalPosition().z() - (*iter)->globalPosition().z()) / ZError) * isParralZ < 0) {
0665         cout << "Pattern find error in Z direction: wrong Z direction" << endl;
0666         isGoodPattern = 0;
0667       }
0668     }
0669   }
0670 
0671   // Check pattern
0672   if (isStraight == false) {
0673     // Check clockwise direction
0674     GlobalVector* vec = new GlobalVector[NumberofHitsinSeed];
0675     unsigned int index = 0;
0676     for (ConstMuonRecHitContainer::const_iterator iter = theRecHits.begin(); iter != theRecHits.end(); iter++) {
0677       GlobalVector vec_temp(((*iter)->globalPosition().x() - Center.x()),
0678                             ((*iter)->globalPosition().y() - Center.y()),
0679                             ((*iter)->globalPosition().z() - Center.z()));
0680       vec[index] = vec_temp;
0681       index++;
0682     }
0683     isClockwise = 0;
0684     for (unsigned int index = 0; index < (NumberofHitsinSeed - 1); index++) {
0685       // Check phi direction, all sub-dphi direction should be the same
0686       if ((vec[index + 1].phi() - vec[index].phi()) > 0)
0687         isClockwise--;
0688       else
0689         isClockwise++;
0690       cout << "Current isClockwise is : " << isClockwise << endl;
0691     }
0692     cout << "Check isClockwise is : " << isClockwise << endl;
0693     if ((unsigned int)abs(isClockwise) != (NumberofHitsinSeed - 1)) {
0694       cout << "Pattern find error in Phi direction" << endl;
0695       isGoodPattern = 0;
0696       isClockwise = 0;
0697     } else
0698       isClockwise /= abs(isClockwise);
0699     delete[] vec;
0700 
0701     // Get meanPt and meanSpt
0702     double deltaRwithBz = fabs(deltaBz * meanRadius / meanBz);
0703     cout << "deltaR with Bz is " << deltaRwithBz << endl;
0704 
0705     if (isClockwise == 0)
0706       meanPt = upper_limit_pt;
0707     else
0708       meanPt = 0.01 * meanRadius * meanBz * 0.3 * isClockwise;
0709     if (fabs(meanPt) > upper_limit_pt)
0710       meanPt = upper_limit_pt * meanPt / fabs(meanPt);
0711     cout << " meanRadius is " << meanRadius << endl;
0712     cout << " meanPt is " << meanPt << endl;
0713 
0714     double deltaR = 0.;
0715     for (ConstMuonRecHitContainer::const_iterator iter = theRecHits.begin(); iter != theRecHits.end(); iter++) {
0716       deltaR += (getDistance(*iter, Center) - meanRadius) * (getDistance(*iter, Center) - meanRadius);
0717     }
0718     deltaR = deltaR / NumberofHitsinSeed;
0719     deltaR = sqrt(deltaR);
0720     //meanSpt = 0.01 * deltaR * meanBz * 0.3;
0721     meanSpt = deltaR;
0722     cout << "DeltaR is " << deltaR << endl;
0723     if (deltaR > deltaRThreshold) {
0724       cout << "Pattern find error: deltaR over threshold" << endl;
0725       isGoodPattern = 0;
0726     }
0727   } else {
0728     // Just set pattern to be straight
0729     isClockwise = 0;
0730     meanPt = upper_limit_pt;
0731     // Set the straight pattern with lowest priority among good pattern
0732     meanSpt = deltaRThreshold;
0733   }
0734   cout << "III--> Seed Pt : " << meanPt << endl;
0735   cout << "III--> Pattern is: " << isGoodPattern << endl;
0736 
0737   // Set the pattern estimation signal
0738   isPatternChecked = true;
0739 }
0740 
0741 void RPCSeedPattern::checkSegmentAlgorithmSpecial(MagneticField const& Field, RPCGeometry const& rpcGeom) {
0742   if (isPatternChecked == true)
0743     return;
0744 
0745   if (!checkSegment()) {
0746     isPatternChecked = true;
0747     isGoodPattern = -1;
0748     return;
0749   }
0750 
0751   // Set isGoodPattern to default true and check the failure
0752   isGoodPattern = 1;
0753 
0754   // Print the recHit's position
0755   for (ConstMuonRecHitContainer::const_iterator iter = theRecHits.begin(); iter != theRecHits.end(); iter++)
0756     cout << "Position of recHit is: " << (*iter)->globalPosition() << endl;
0757 
0758   // Check the Z direction
0759   if (fabs((*(theRecHits.end() - 1))->globalPosition().z() - (*(theRecHits.begin()))->globalPosition().z()) > ZError) {
0760     if (((*(theRecHits.end() - 1))->globalPosition().z() - (*(theRecHits.begin()))->globalPosition().z()) > ZError)
0761       isParralZ = 1;
0762     else
0763       isParralZ = -1;
0764   } else
0765     isParralZ = 0;
0766 
0767   cout << " Check isParralZ is :" << isParralZ << endl;
0768   for (ConstMuonRecHitContainer::const_iterator iter = theRecHits.begin(); iter != (theRecHits.end() - 1); iter++) {
0769     if (isParralZ == 0) {
0770       if (fabs((*(iter + 1))->globalPosition().z() - (*iter)->globalPosition().z()) > ZError) {
0771         cout << "Pattern find error in Z direction: wrong perpendicular direction" << endl;
0772         isGoodPattern = 0;
0773       }
0774     } else {
0775       if ((int)(((*(iter + 1))->globalPosition().z() - (*iter)->globalPosition().z()) / ZError) * isParralZ < 0) {
0776         cout << "Pattern find error in Z direction: wrong Z direction" << endl;
0777         isGoodPattern = 0;
0778       }
0779     }
0780   }
0781 
0782   // Check the pattern
0783   if (isStraight2 == true) {
0784     // Set pattern to be straight
0785     isClockwise = 0;
0786     meanPt = upper_limit_pt;
0787     // Set the straight pattern with lowest priority among good pattern
0788     meanSpt = deltaRThreshold;
0789 
0790     // Extrapolate to other recHits and check deltaR
0791     GlobalVector startSegment = (SegmentRB[1].second)->globalPosition() - (SegmentRB[1].first)->globalPosition();
0792     GlobalPoint startPosition = (SegmentRB[1].first)->globalPosition();
0793     GlobalVector startMomentum = startSegment * (meanPt / startSegment.perp());
0794     unsigned int index = 0;
0795 
0796     for (ConstMuonRecHitContainer::const_iterator iter = theRecHits.begin(); iter != theRecHits.end(); iter++) {
0797       if (index < 4) {
0798         index++;
0799         continue;
0800       }
0801       double tracklength = 0;
0802       cout << "Now checking recHit " << index << endl;
0803       double Distance = extropolateStep(startPosition, startMomentum, iter, isClockwise, tracklength, Field, rpcGeom);
0804       cout << "Final distance is " << Distance << endl;
0805       if (Distance > MaxRSD) {
0806         cout << "Pattern find error in distance for other recHits: " << Distance << endl;
0807         isGoodPattern = 0;
0808       }
0809       index++;
0810     }
0811   } else {
0812     // Get clockwise direction
0813     GlobalVector vec[2];
0814     vec[0] = GlobalVector(
0815         (entryPosition.x() - Center2.x()), (entryPosition.y() - Center2.y()), (entryPosition.z() - Center2.z()));
0816     vec[1] = GlobalVector(
0817         (leavePosition.x() - Center2.x()), (leavePosition.y() - Center2.y()), (leavePosition.z() - Center2.z()));
0818     isClockwise = 0;
0819     if ((vec[1].phi() - vec[0].phi()).value() > 0)
0820       isClockwise = -1;
0821     else
0822       isClockwise = 1;
0823 
0824     cout << "Check isClockwise is : " << isClockwise << endl;
0825 
0826     // Get meanPt
0827     meanPt = 0.01 * meanRadius2 * meanMagneticField2.z() * 0.3 * isClockwise;
0828     //meanPt = 0.01 * meanRadius2[0] * (-3.8) * 0.3 * isClockwise;
0829     cout << " meanRadius is " << meanRadius2 << ", with meanBz " << meanMagneticField2.z() << endl;
0830     cout << " meanPt is " << meanPt << endl;
0831     if (fabs(meanPt) > upper_limit_pt)
0832       meanPt = upper_limit_pt * meanPt / fabs(meanPt);
0833 
0834     // Check the initial 3 segments
0835     cout << "entryPosition: " << entryPosition << endl;
0836     cout << "leavePosition: " << leavePosition << endl;
0837     cout << "Center2 is : " << Center2 << endl;
0838     double R1 = vec[0].perp();
0839     double R2 = vec[1].perp();
0840     double deltaR = sqrt(((R1 - meanRadius2) * (R1 - meanRadius2) + (R2 - meanRadius2) * (R2 - meanRadius2)) / 2);
0841     meanSpt = deltaR;
0842     cout << "R1 is " << R1 << ", R2 is " << R2 << endl;
0843     cout << "Delta R for the initial 3 segments is " << deltaR << endl;
0844     if (deltaR > deltaRThreshold) {
0845       cout << "Pattern find error in delta R for the initial 3 segments" << endl;
0846       isGoodPattern = 0;
0847     }
0848 
0849     // Extrapolate to other recHits and check deltaR
0850     GlobalVector startSegment = (SegmentRB[1].second)->globalPosition() - (SegmentRB[1].first)->globalPosition();
0851     GlobalPoint startPosition = (SegmentRB[1].first)->globalPosition();
0852     GlobalVector startMomentum = startSegment * (fabs(meanPt) / startSegment.perp());
0853     unsigned int index = 0;
0854     for (ConstMuonRecHitContainer::const_iterator iter = theRecHits.begin(); iter != theRecHits.end(); iter++) {
0855       if (index < 4) {
0856         index++;
0857         continue;
0858       }
0859       double tracklength = 0;
0860       cout << "Now checking recHit " << index << endl;
0861       double Distance = extropolateStep(startPosition, startMomentum, iter, isClockwise, tracklength, Field, rpcGeom);
0862       cout << "Final distance is " << Distance << endl;
0863       if (Distance > MaxRSD) {
0864         cout << "Pattern find error in distance for other recHits: " << Distance << endl;
0865         isGoodPattern = 0;
0866       }
0867       index++;
0868     }
0869   }
0870 
0871   cout << "Checking finish, isGoodPattern now is " << isGoodPattern << endl;
0872   // Set the pattern estimation signal
0873   isPatternChecked = true;
0874 }
0875 
0876 double RPCSeedPattern::extropolateStep(const GlobalPoint& startPosition,
0877                                        const GlobalVector& startMomentum,
0878                                        ConstMuonRecHitContainer::const_iterator iter,
0879                                        const int ClockwiseDirection,
0880                                        double& tracklength,
0881                                        const MagneticField& Field,
0882                                        const RPCGeometry& rpcGeometry) {
0883   cout << "Extrapolating the track to check the pattern" << endl;
0884   tracklength = 0;
0885   // Get the iter recHit's detector geometry
0886   DetId hitDet = (*iter)->hit()->geographicalId();
0887   RPCDetId RPCId = RPCDetId(hitDet.rawId());
0888   //const RPCChamber* hitRPC = dynamic_cast<const RPCChamber*>(hitDet);
0889 
0890   const BoundPlane RPCSurface = rpcGeometry.chamber(RPCId)->surface();
0891   double startSide = RPCSurface.localZ(startPosition);
0892   cout << "Start side : " << startSide;
0893 
0894   GlobalPoint currentPosition = startPosition;
0895   GlobalVector currentMomentum = startMomentum;
0896   GlobalVector ZDirection(0, 0, 1);
0897 
0898   // Use the perp other than mag, since initial segment might have small value while final recHit have large difference value at Z direction
0899   double currentDistance = ((GlobalVector)(currentPosition - (*iter)->globalPosition())).perp();
0900   cout << "Start current position is : " << currentPosition << endl;
0901   cout << "Start current Momentum is: " << currentMomentum.mag() << ", in vector: " << currentMomentum << endl;
0902   cout << "Start current distance is " << currentDistance << endl;
0903   cout << "Start current radius is " << currentPosition.perp() << endl;
0904   cout << "Destination radius is " << (*iter)->globalPosition().perp() << endl;
0905 
0906   // Judge roughly if the stepping cross the Det surface of the recHit
0907   //while((currentPosition.perp() < ((*iter)->globalPosition().perp())))
0908   double currentDistance_next = currentDistance;
0909   do {
0910     currentDistance = currentDistance_next;
0911     if (ClockwiseDirection == 0) {
0912       currentPosition += currentMomentum.unit() * stepLength;
0913     } else {
0914       double Bz = Field.inTesla(currentPosition).z();
0915       double Radius = currentMomentum.perp() / fabs(Bz * 0.01 * 0.3);
0916       double deltaPhi = (stepLength * currentMomentum.perp() / currentMomentum.mag()) / Radius;
0917 
0918       // Get the center for current step
0919       GlobalVector currentPositiontoCenter = currentMomentum.unit().cross(ZDirection);
0920       currentPositiontoCenter *= Radius;
0921       // correction of ClockwiseDirection correction
0922       currentPositiontoCenter *= ClockwiseDirection;
0923       // continue to get the center for current step
0924       GlobalPoint currentCenter = currentPosition;
0925       currentCenter += currentPositiontoCenter;
0926 
0927       // Get the next step position
0928       GlobalVector CentertocurrentPosition = (GlobalVector)(currentPosition - currentCenter);
0929       double Phi = CentertocurrentPosition.phi().value();
0930       Phi += deltaPhi * (-1) * ClockwiseDirection;
0931       double deltaZ = stepLength * currentMomentum.z() / currentMomentum.mag();
0932       GlobalVector CentertonewPosition(GlobalVector::Cylindrical(CentertocurrentPosition.perp(), Phi, deltaZ));
0933       double PtPhi = currentMomentum.phi().value();
0934       PtPhi += deltaPhi * (-1) * ClockwiseDirection;
0935       currentMomentum = GlobalVector(GlobalVector::Cylindrical(currentMomentum.perp(), PtPhi, currentMomentum.z()));
0936       currentPosition = currentCenter + CentertonewPosition;
0937     }
0938 
0939     // count the total step length
0940     tracklength += stepLength * currentMomentum.perp() / currentMomentum.mag();
0941 
0942     // Get the next step distance
0943     double currentSide = RPCSurface.localZ(currentPosition);
0944     cout << "Stepping current side : " << currentSide << endl;
0945     cout << "Stepping current position is: " << currentPosition << endl;
0946     cout << "Stepping current Momentum is: " << currentMomentum.mag() << ", in vector: " << currentMomentum << endl;
0947     currentDistance_next = ((GlobalVector)(currentPosition - (*iter)->globalPosition())).perp();
0948     cout << "Stepping current distance is " << currentDistance << endl;
0949     cout << "Stepping current radius is " << currentPosition.perp() << endl;
0950   } while (currentDistance_next < currentDistance);
0951 
0952   return currentDistance;
0953 }
0954 
0955 RPCSeedPattern::weightedTrajectorySeed RPCSeedPattern::createFakeSeed(int& isGoodSeed, const MagneticField& Field) {
0956   // Create a fake seed and return
0957   cout << "Now create a fake seed" << endl;
0958   isPatternChecked = true;
0959   isGoodPattern = -1;
0960   isStraight = true;
0961   meanPt = upper_limit_pt;
0962   meanSpt = 0;
0963   Charge = 0;
0964   isClockwise = 0;
0965   isParralZ = 0;
0966   meanRadius = -1;
0967   //return createSeed(isGoodSeed, eSetup);
0968 
0969   // Get the reference recHit, DON'T use the recHit on 1st layer(inner most layer)
0970   const ConstMuonRecHitPointer best = BestRefRecHit();
0971 
0972   Momentum = GlobalVector(0, 0, 0);
0973   LocalPoint segPos = best->localPosition();
0974   LocalVector segDirFromPos = best->det()->toLocal(Momentum);
0975   LocalTrajectoryParameters param(segPos, segDirFromPos, Charge);
0976 
0977   //AlgebraicVector t(4);
0978   AlgebraicSymMatrix mat(5, 0);
0979   mat = best->parametersError().similarityT(best->projectionMatrix());
0980   mat[0][0] = meanSpt;
0981   LocalTrajectoryError error(asSMatrix<5>(mat));
0982 
0983   TrajectoryStateOnSurface tsos(param, error, best->det()->surface(), &Field);
0984 
0985   DetId id = best->geographicalId();
0986 
0987   PTrajectoryStateOnDet seedTSOS = trajectoryStateTransform::persistentState(tsos, id.rawId());
0988 
0989   edm::OwnVector<TrackingRecHit> container;
0990   for (ConstMuonRecHitContainer::const_iterator iter = theRecHits.begin(); iter != theRecHits.end(); iter++)
0991     container.push_back((*iter)->hit()->clone());
0992 
0993   TrajectorySeed theSeed(seedTSOS, container, alongMomentum);
0994   weightedTrajectorySeed theweightedSeed;
0995   theweightedSeed.first = theSeed;
0996   theweightedSeed.second = meanSpt;
0997   isGoodSeed = isGoodPattern;
0998 
0999   return theweightedSeed;
1000 }
1001 
1002 RPCSeedPattern::weightedTrajectorySeed RPCSeedPattern::createSeed(int& isGoodSeed, const MagneticField& Field) {
1003   if (isPatternChecked == false || isGoodPattern == -1) {
1004     cout << "Pattern is not yet checked! Create a fake seed instead!" << endl;
1005     return createFakeSeed(isGoodSeed, Field);
1006   }
1007 
1008   MuonPatternRecoDumper debug;
1009 
1010   //double theMinMomentum = 3.0;
1011   //if(fabs(meanPt) < lower_limit_pt)
1012   //meanPt = lower_limit_pt * meanPt / fabs(meanPt);
1013 
1014   // For pattern we use is Clockwise other than isStraight to estimate charge
1015   if (isClockwise == 0)
1016     Charge = 0;
1017   else
1018     Charge = (int)(meanPt / fabs(meanPt));
1019 
1020   // Get the reference recHit, DON'T use the recHit on 1st layer(inner most layer)
1021   const ConstMuonRecHitPointer best = BestRefRecHit();
1022   const ConstMuonRecHitPointer first = FirstRecHit();
1023 
1024   if (isClockwise != 0) {
1025     if (AlgorithmType != 3) {
1026       // Get the momentum on reference recHit
1027       GlobalVector vecRef1((first->globalPosition().x() - Center.x()),
1028                            (first->globalPosition().y() - Center.y()),
1029                            (first->globalPosition().z() - Center.z()));
1030       GlobalVector vecRef2((best->globalPosition().x() - Center.x()),
1031                            (best->globalPosition().y() - Center.y()),
1032                            (best->globalPosition().z() - Center.z()));
1033 
1034       double deltaPhi = (vecRef2.phi() - vecRef1.phi()).value();
1035       double deltaS = meanRadius * fabs(deltaPhi);
1036       double deltaZ = best->globalPosition().z() - first->globalPosition().z();
1037 
1038       GlobalVector vecZ(0, 0, 1);
1039       GlobalVector vecPt = (vecRef2.unit()).cross(vecZ);
1040       if (isClockwise == -1)
1041         vecPt *= -1;
1042       vecPt *= deltaS;
1043       Momentum = GlobalVector(0, 0, deltaZ);
1044       Momentum += vecPt;
1045       Momentum *= fabs(meanPt / deltaS);
1046     } else {
1047       double deltaZ = best->globalPosition().z() - first->globalPosition().z();
1048       Momentum = GlobalVector(GlobalVector::Cylindrical(S, lastPhi, deltaZ));
1049       Momentum *= fabs(meanPt / S);
1050     }
1051   } else {
1052     Momentum = best->globalPosition() - first->globalPosition();
1053     double deltaS = Momentum.perp();
1054     Momentum *= fabs(meanPt / deltaS);
1055   }
1056   LocalPoint segPos = best->localPosition();
1057   LocalVector segDirFromPos = best->det()->toLocal(Momentum);
1058   LocalTrajectoryParameters param(segPos, segDirFromPos, Charge);
1059 
1060   LocalTrajectoryError error = getSpecialAlgorithmErrorMatrix(first, best);
1061 
1062   TrajectoryStateOnSurface tsos(param, error, best->det()->surface(), &Field);
1063   cout << "Trajectory State on Surface before the extrapolation" << endl;
1064   cout << debug.dumpTSOS(tsos);
1065   DetId id = best->geographicalId();
1066   cout << "The RecSegment relies on: " << endl;
1067   cout << debug.dumpMuonId(id);
1068   cout << debug.dumpTSOS(tsos);
1069 
1070   PTrajectoryStateOnDet const& seedTSOS = trajectoryStateTransform::persistentState(tsos, id.rawId());
1071 
1072   edm::OwnVector<TrackingRecHit> container;
1073   for (ConstMuonRecHitContainer::const_iterator iter = theRecHits.begin(); iter != theRecHits.end(); iter++) {
1074     // This casting withou clone will cause memory overflow when used in push_back
1075     // Since container's deconstructor functiion free the pointer menber!
1076     //TrackingRecHit* pt = dynamic_cast<TrackingRecHit*>(&*(*iter));
1077     //cout << "Push recHit type " << pt->getType() << endl;
1078     container.push_back((*iter)->hit()->clone());
1079   }
1080 
1081   TrajectorySeed theSeed(seedTSOS, container, alongMomentum);
1082   weightedTrajectorySeed theweightedSeed;
1083   theweightedSeed.first = theSeed;
1084   theweightedSeed.second = meanSpt;
1085   isGoodSeed = isGoodPattern;
1086 
1087   return theweightedSeed;
1088 }
1089 
1090 LocalTrajectoryError RPCSeedPattern::getSpecialAlgorithmErrorMatrix(const ConstMuonRecHitPointer& first,
1091                                                                     const ConstMuonRecHitPointer& best) {
1092   LocalTrajectoryError Error;
1093   double dXdZ = 0;
1094   double dYdZ = 0;
1095   double dP = 0;
1096   AlgebraicSymMatrix mat(5, 0);
1097   mat = best->parametersError().similarityT(best->projectionMatrix());
1098   if (AlgorithmType != 3) {
1099     GlobalVector vecRef1((first->globalPosition().x() - Center.x()),
1100                          (first->globalPosition().y() - Center.y()),
1101                          (first->globalPosition().z() - Center.z()));
1102     GlobalVector vecRef2((best->globalPosition().x() - Center.x()),
1103                          (best->globalPosition().y() - Center.y()),
1104                          (best->globalPosition().z() - Center.z()));
1105     double deltaPhi = (vecRef2.phi() - vecRef1.phi()).value();
1106     double L = meanRadius * fabs(deltaPhi);
1107     double N = nrhit();
1108     double A_N = 180 * N * N * N / ((N - 1) * (N + 1) * (N + 2) * (N + 3));
1109     double sigma_x = sqrt(mat[3][3]);
1110     double betaovergame = Momentum.mag() / 0.1066;
1111     double beta = sqrt((betaovergame * betaovergame) / (1 + betaovergame * betaovergame));
1112     double dPt = meanPt * (0.0136 * sqrt(1 / 100) * sqrt(4 * A_N / N) / (beta * 0.3 * meanBz * L) +
1113                            sigma_x * meanPt * sqrt(4 * A_N) / (0.3 * meanBz * L * L));
1114     double dP = dPt * Momentum.mag() / meanPt;
1115     mat[0][0] = (dP * dP) / (Momentum.mag() * Momentum.mag() * Momentum.mag() * Momentum.mag());
1116     mat[1][1] = dXdZ * dXdZ;
1117     mat[2][2] = dYdZ * dYdZ;
1118     Error = LocalTrajectoryError(asSMatrix<5>(mat));
1119   } else {
1120     AlgebraicSymMatrix mat0(5, 0);
1121     mat0 = (SegmentRB[0].first)->parametersError().similarityT((SegmentRB[0].first)->projectionMatrix());
1122     double dX0 = sqrt(mat0[3][3]);
1123     double dY0 = sqrt(mat0[4][4]);
1124     AlgebraicSymMatrix mat1(5, 0);
1125     mat1 = (SegmentRB[0].second)->parametersError().similarityT((SegmentRB[0].second)->projectionMatrix());
1126     double dX1 = sqrt(mat1[3][3]);
1127     double dY1 = sqrt(mat1[4][4]);
1128     AlgebraicSymMatrix mat2(5, 0);
1129     mat2 = (SegmentRB[1].first)->parametersError().similarityT((SegmentRB[1].first)->projectionMatrix());
1130     double dX2 = sqrt(mat2[3][3]);
1131     double dY2 = sqrt(mat2[4][4]);
1132     AlgebraicSymMatrix mat3(5, 0);
1133     mat3 = (SegmentRB[1].second)->parametersError().similarityT((SegmentRB[1].second)->projectionMatrix());
1134     double dX3 = sqrt(mat3[3][3]);
1135     double dY3 = sqrt(mat3[4][4]);
1136     cout << "Local error for 4 recHits are: " << dX0 << ", " << dY0 << ", " << dX1 << ", " << dY1 << ", " << dX2 << ", "
1137          << dY2 << ", " << dX3 << ", " << dY3 << endl;
1138     const GeomDetUnit* refRPC1 = (SegmentRB[0].second)->detUnit();
1139     LocalPoint recHit0 = refRPC1->toLocal((SegmentRB[0].first)->globalPosition());
1140     LocalPoint recHit1 = refRPC1->toLocal((SegmentRB[0].second)->globalPosition());
1141     LocalVector localSegment00 = (LocalVector)(recHit1 - recHit0);
1142     LocalVector localSegment01 = LocalVector(localSegment00.x() + dX0 + dX1, localSegment00.y(), localSegment00.z());
1143     LocalVector localSegment02 = LocalVector(localSegment00.x() - dX0 - dX1, localSegment00.y(), localSegment00.z());
1144     GlobalVector globalSegment00 = refRPC1->toGlobal(localSegment00);
1145     GlobalVector globalSegment01 = refRPC1->toGlobal(localSegment01);
1146     GlobalVector globalSegment02 = refRPC1->toGlobal(localSegment02);
1147 
1148     const GeomDetUnit* refRPC2 = (SegmentRB[1].first)->detUnit();
1149     LocalPoint recHit2 = refRPC2->toLocal((SegmentRB[1].first)->globalPosition());
1150     LocalPoint recHit3 = refRPC2->toLocal((SegmentRB[1].second)->globalPosition());
1151     LocalVector localSegment10 = (LocalVector)(recHit3 - recHit2);
1152     LocalVector localSegment11 = LocalVector(localSegment10.x() + dX2 + dX3, localSegment10.y(), localSegment10.z());
1153     LocalVector localSegment12 = LocalVector(localSegment10.x() - dX2 - dX3, localSegment10.y(), localSegment10.z());
1154     GlobalVector globalSegment10 = refRPC2->toGlobal(localSegment10);
1155     GlobalVector globalSegment11 = refRPC2->toGlobal(localSegment11);
1156     GlobalVector globalSegment12 = refRPC2->toGlobal(localSegment12);
1157 
1158     if (isClockwise != 0) {
1159       GlobalVector vec[2];
1160       vec[0] = GlobalVector(
1161           (entryPosition.x() - Center2.x()), (entryPosition.y() - Center2.y()), (entryPosition.z() - Center2.z()));
1162       vec[1] = GlobalVector(
1163           (leavePosition.x() - Center2.x()), (leavePosition.y() - Center2.y()), (leavePosition.z() - Center2.z()));
1164       double halfPhiCenter = fabs((vec[1].phi() - vec[0].phi()).value()) / 2;
1165       // dPhi0 shoule be the same clockwise direction, while dPhi1 should be opposite clockwise direction, w.r.t to the track clockwise
1166       double dPhi0 = (((globalSegment00.phi() - globalSegment01.phi()).value() * isClockwise) > 0)
1167                          ? fabs((globalSegment00.phi() - globalSegment01.phi()).value())
1168                          : fabs((globalSegment00.phi() - globalSegment02.phi()).value());
1169       double dPhi1 = (((globalSegment10.phi() - globalSegment11.phi()).value() * isClockwise) < 0)
1170                          ? fabs((globalSegment10.phi() - globalSegment11.phi()).value())
1171                          : fabs((globalSegment10.phi() - globalSegment12.phi()).value());
1172       // For the deltaR should be kept small, we assume the delta Phi0/Phi1 should be in a same limit value
1173       double dPhi = (dPhi0 <= dPhi1) ? dPhi0 : dPhi1;
1174       cout << "DPhi for new Segment is about " << dPhi << endl;
1175       // Check the variance of halfPhiCenter
1176       double newhalfPhiCenter = ((halfPhiCenter - dPhi) > 0 ? (halfPhiCenter - dPhi) : 0);
1177       if (newhalfPhiCenter != 0) {
1178         double newmeanPt = meanPt * halfPhiCenter / newhalfPhiCenter;
1179         if (fabs(newmeanPt) > upper_limit_pt)
1180           newmeanPt = upper_limit_pt * meanPt / fabs(meanPt);
1181         cout << "The error is inside range. Max meanPt could be " << newmeanPt << endl;
1182         dP = fabs(Momentum.mag() * (newmeanPt - meanPt) / meanPt);
1183       } else {
1184         double newmeanPt = upper_limit_pt * meanPt / fabs(meanPt);
1185         cout << "The error is outside range. Max meanPt could be " << newmeanPt << endl;
1186         dP = fabs(Momentum.mag() * (newmeanPt - meanPt) / meanPt);
1187       }
1188     } else {
1189       double dPhi0 = (fabs((globalSegment00.phi() - globalSegment01.phi()).value()) <=
1190                       fabs((globalSegment00.phi() - globalSegment02.phi()).value()))
1191                          ? fabs((globalSegment00.phi() - globalSegment01.phi()).value())
1192                          : fabs((globalSegment00.phi() - globalSegment02.phi()).value());
1193       double dPhi1 = (fabs((globalSegment10.phi() - globalSegment11.phi()).value()) <=
1194                       fabs((globalSegment10.phi() - globalSegment12.phi()).value()))
1195                          ? fabs((globalSegment10.phi() - globalSegment11.phi()).value())
1196                          : fabs((globalSegment10.phi() - globalSegment12.phi()).value());
1197       double dPhi = (dPhi0 <= dPhi1) ? dPhi0 : dPhi1;
1198       GlobalVector middleSegment = leavePosition - entryPosition;
1199       double halfDistance = middleSegment.perp() / 2;
1200       double newmeanPt = halfDistance / dPhi;
1201       cout << "The error is for straight. Max meanPt could be " << newmeanPt << endl;
1202       dP = fabs(Momentum.mag() * (newmeanPt - meanPt) / meanPt);
1203     }
1204 
1205     double dXdZ1 = globalSegment11.x() / globalSegment11.z() - globalSegment10.x() / globalSegment10.z();
1206     double dXdZ2 = globalSegment12.x() / globalSegment12.z() - globalSegment10.x() / globalSegment10.z();
1207     dXdZ = (fabs(dXdZ1) >= fabs(dXdZ2)) ? dXdZ1 : dXdZ2;
1208 
1209     LocalVector localSegment13 = LocalVector(localSegment10.x(), localSegment10.y() + dY2 + dY3, localSegment10.z());
1210     LocalVector localSegment14 = LocalVector(localSegment10.x(), localSegment10.y() - dY2 - dY3, localSegment10.z());
1211     GlobalVector globalSegment13 = refRPC2->toGlobal(localSegment13);
1212     GlobalVector globalSegment14 = refRPC2->toGlobal(localSegment14);
1213     double dYdZ1 = globalSegment13.y() / globalSegment13.z() - globalSegment10.y() / globalSegment10.z();
1214     double dYdZ2 = globalSegment14.y() / globalSegment14.z() - globalSegment10.y() / globalSegment10.z();
1215     dYdZ = (fabs(dYdZ1) >= fabs(dYdZ2)) ? dYdZ1 : dYdZ2;
1216 
1217     mat[0][0] = (dP * dP) / (Momentum.mag() * Momentum.mag() * Momentum.mag() * Momentum.mag());
1218     mat[1][1] = dXdZ * dXdZ;
1219     mat[2][2] = dYdZ * dYdZ;
1220     Error = LocalTrajectoryError(asSMatrix<5>(mat));
1221   }
1222   return Error;
1223 }
[ 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!