Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:28:02 

0001 #include "SeedForPhotonConversionFromQuadruplets.h"
0002 
0003 #include "Conv4HitsReco2.h"
0004 
0005 #include "TrackingTools/KalmanUpdators/interface/KFUpdator.h"
0006 #include "RecoTracker/TkSeedGenerator/interface/FastHelix.h"
0007 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0008 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
0009 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0010 #include "MagneticField/Engine/interface/MagneticField.h"
0011 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
0012 #include "TrackingTools/Records/interface/TrackingComponentsRecord.h"
0013 #include "TrackingTools/Records/interface/TransientRecHitRecord.h"
0014 #include "TrackingTools/MaterialEffects/interface/PropagatorWithMaterial.h"
0015 #include "TrackingTools/GeomPropagators/interface/PropagationExceptions.h"
0016 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
0017 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
0018 #include "TrackingTools/Records/interface/TransientRecHitRecord.h"
0019 //ClusterShapeIncludes
0020 #include "RecoTracker/TkTrackingRegions/interface/OrderedHitsGenerator.h"
0021 #include "RecoTracker/TkTrackingRegions/interface/TrackingRegion.h"
0022 #include "RecoTracker/TkSeedingLayers/interface/SeedComparitor.h"
0023 #include "RecoTracker/TkSeedGenerator/interface/SeedCreator.h"
0024 #include "RecoTracker/TkSeedingLayers/interface/SeedComparitorFactory.h"
0025 #include "RecoTracker/TkSeedGenerator/interface/SeedCreatorFactory.h"
0026 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateClosestToBeamLine.h"
0027 #include "TrackingTools/PatternTools/interface/TSCBLBuilderNoMaterial.h"
0028 #include "DataFormats/BeamSpot/interface/BeamSpot.h"
0029 
0030 //#define mydebug_knuenz
0031 //#define mydebug_sguazz
0032 //#define quadDispLine
0033 template <class T>
0034 T sqr(T t) {
0035   return t * t;
0036 }
0037 
0038 SeedForPhotonConversionFromQuadruplets::SeedForPhotonConversionFromQuadruplets(
0039     const edm::ParameterSet& cfg, edm::ConsumesCollector& iC, const edm::ParameterSet& SeedComparitorPSet)
0040     : SeedForPhotonConversionFromQuadruplets(
0041           iC,
0042           SeedComparitorPSet,
0043           cfg.getParameter<std::string>("propagator"),
0044           cfg.existsAs<double>("SeedMomentumForBOFF") ? cfg.getParameter<double>("SeedMomentumForBOFF") : 5.0) {}
0045 
0046 SeedForPhotonConversionFromQuadruplets::SeedForPhotonConversionFromQuadruplets(
0047     edm::ConsumesCollector& iC,
0048     const edm::ParameterSet& SeedComparitorPSet,
0049     const std::string& propagator,
0050     double seedMomentumForBOFF)
0051     : theBfieldToken(iC.esConsumes()),
0052       theTTRHBuilderToken(iC.esConsumes(edm::ESInputTag("", "WithTrackAngle"))),  // FIXME: add to config
0053       theTrackerToken(iC.esConsumes()),
0054       thePropagatorToken(iC.esConsumes(edm::ESInputTag("", propagator))),
0055       theBOFFMomentum(seedMomentumForBOFF) {
0056   std::string comparitorName = SeedComparitorPSet.getParameter<std::string>("ComponentName");
0057   if (comparitorName != "none") {
0058     theComparitor = SeedComparitorFactory::get()->create(comparitorName, SeedComparitorPSet, iC);
0059   }
0060 }
0061 
0062 SeedForPhotonConversionFromQuadruplets::~SeedForPhotonConversionFromQuadruplets() {}
0063 
0064 const TrajectorySeed* SeedForPhotonConversionFromQuadruplets::trajectorySeed(TrajectorySeedCollection& seedCollection,
0065                                                                              const SeedingHitSet& phits,
0066                                                                              const SeedingHitSet& mhits,
0067                                                                              const TrackingRegion& region,
0068                                                                              const edm::Event& ev,
0069                                                                              const edm::EventSetup& es,
0070                                                                              std::stringstream& ss,
0071                                                                              std::vector<Quad>& quadV,
0072                                                                              edm::ParameterSet& QuadCutPSet) {
0073   //// CUT DEFINITIONS ////
0074 
0075   bool rejectAllQuads = QuadCutPSet.getParameter<bool>("rejectAllQuads");
0076   if (rejectAllQuads)
0077     return nullptr;
0078 
0079   bool applyDeltaPhiCuts = QuadCutPSet.getParameter<bool>("apply_DeltaPhiCuts");
0080   bool ClusterShapeFiltering = QuadCutPSet.getParameter<bool>("apply_ClusterShapeFilter");
0081   bool applyArbitration = QuadCutPSet.getParameter<bool>("apply_Arbitration");
0082   bool applydzCAcut = QuadCutPSet.getParameter<bool>("apply_zCACut");
0083   double CleaningmaxRadialDistance = QuadCutPSet.getParameter<double>("Cut_DeltaRho");  //cm
0084   double BeamPipeRadiusCut = QuadCutPSet.getParameter<double>("Cut_BeamPipeRadius");    //cm
0085   double CleaningMinLegPt = QuadCutPSet.getParameter<double>("Cut_minLegPt");           //GeV
0086   double maxLegPt = QuadCutPSet.getParameter<double>("Cut_maxLegPt");                   //GeV
0087   double dzcut = QuadCutPSet.getParameter<double>("Cut_zCA");                           //cm
0088 
0089   double toleranceFactorOnDeltaPhiCuts = 0.1;
0090 
0091   //  return 0; //FIXME, remove this line to make the code working.
0092 
0093   pss = &ss;
0094 
0095   if (phits.size() < 2)
0096     return nullptr;
0097   if (mhits.size() < 2)
0098     return nullptr;
0099 
0100     //PUT HERE THE QUADRUPLET ALGORITHM, AND IN CASE USE THE METHODS ALREADY DEVELOPED, ADAPTING THEM
0101 
0102     //
0103     // Rozzo ma efficace (per ora)
0104     //
0105 
0106 #ifdef mydebug_sguazz
0107   std::cout << " --------------------------------------------------------------------------"
0108             << "\n";
0109   std::cout << "  Starting a hit quad fast reco "
0110             << "\n";
0111   std::cout << " --------------------------------------------------------------------------"
0112             << "\n";
0113 #endif
0114 
0115   //
0116   // Let's build the 4 hits
0117   SeedingHitSet::ConstRecHitPointer ptth1 = phits[0];
0118   SeedingHitSet::ConstRecHitPointer ptth2 = phits[1];
0119   SeedingHitSet::ConstRecHitPointer mtth1 = mhits[0];
0120   SeedingHitSet::ConstRecHitPointer mtth2 = mhits[1];
0121 
0122   GlobalPoint vHit[4];
0123   vHit[0] = ptth2->globalPosition();
0124   vHit[1] = ptth1->globalPosition();
0125   vHit[2] = mtth1->globalPosition();
0126   vHit[3] = mtth2->globalPosition();
0127 
0128   //Photon source vertex primary vertex
0129   const GlobalPoint& vgPhotVertex = region.origin();
0130   math::XYZVector vPhotVertex(vgPhotVertex.x(), vgPhotVertex.y(), vgPhotVertex.z());
0131 
0132   math::XYZVector h1(vHit[0].x(), vHit[0].y(), vHit[0].z());
0133   math::XYZVector h2(vHit[1].x(), vHit[1].y(), vHit[1].z());
0134   math::XYZVector h3(vHit[2].x(), vHit[2].y(), vHit[2].z());
0135   math::XYZVector h4(vHit[3].x(), vHit[3].y(), vHit[3].z());
0136 
0137   // At this point implement cleaning cuts before building the seed:::
0138 
0139   /*
0140   Notes:
0141 
0142 h1, h2: positron
0143 h3, h4: electron
0144 
0145 P1, P2: positron, ordered with radius
0146 M1, M2: electron, ordered with radius
0147 
0148 Evan's notation:
0149 V1=P1
0150 V2=M1
0151 V3=P2
0152 V4=M2
0153 
0154 */
0155 
0156   math::XYZVector P1;
0157   math::XYZVector P2;
0158   math::XYZVector M1;
0159   math::XYZVector M2;
0160 
0161   if (h1.x() * h1.x() + h1.y() * h1.y() < h2.x() * h2.x() + h2.y() * h2.y()) {
0162     P1 = h1;
0163     P2 = h2;
0164   } else {
0165     P1 = h2;
0166     P2 = h1;
0167   }
0168 
0169   if (h3.x() * h3.x() + h3.y() * h3.y() < h4.x() * h4.x() + h4.y() * h4.y()) {
0170     M1 = h3;
0171     M2 = h4;
0172   } else {
0173     M1 = h4;
0174     M2 = h3;
0175   }
0176 
0177   ////////////////////////
0178   // Intersection-point:::
0179   ////////////////////////
0180   /*
0181 Calculate the intersection point of the lines P2-P1 and M2-M1.
0182 If this point is in the beam pipe, or if the distance of this
0183 point to the layer of the most inner hit of the seed is less
0184 than CleaningmaxRadialDistance cm, the combination is rejected.
0185 */
0186 
0187   math::XYZVector IP(0, 0, 0);
0188 
0189   //Line1:
0190   double kP = (P1.y() - P2.y()) / (P1.x() - P2.x());
0191   double dP = P1.y() - kP * P1.x();
0192   //Line2:
0193   double kM = (M1.y() - M2.y()) / (M1.x() - M2.x());
0194   double dM = M1.y() - kM * M1.x();
0195   //Intersection:
0196   double IPx = (dM - dP) / (kP - kM);
0197   double IPy = kP * IPx + dP;
0198 
0199   IP.SetXYZ(IPx, IPy, 0);
0200 
0201   double IPrho = std::sqrt(IP.x() * IP.x() + IP.y() * IP.y());
0202   double P1rho2 = P1.x() * P1.x() + P1.y() * P1.y();
0203   double M1rho2 = M1.x() * M1.x() + M1.y() * M1.y();
0204   double maxIPrho2 = IPrho + CleaningmaxRadialDistance;
0205   maxIPrho2 *= maxIPrho2;
0206 
0207   if (IPrho < BeamPipeRadiusCut || P1rho2 > maxIPrho2 || M1rho2 > maxIPrho2) {
0208     return nullptr;
0209   }
0210 
0211   if (applyDeltaPhiCuts) {
0212     kPI_ = std::atan(1.0) * 4;
0213 
0214     const auto& bfield = es.getData(theBfieldToken);
0215     math::XYZVector QuadMean(0, 0, 0);
0216     QuadMean.SetXYZ((M1.x() + M2.x() + P1.x() + P2.x()) / 4.,
0217                     (M1.y() + M2.y() + P1.y() + P2.y()) / 4.,
0218                     (M1.z() + M2.z() + P1.z() + P2.z()) / 4.);
0219 
0220     double fBField = bfield.inTesla(GlobalPoint(QuadMean.x(), QuadMean.y(), QuadMean.z())).z();
0221 
0222     double rMax = CleaningMinLegPt / (0.01 * 0.3 * fBField);
0223     double rMax_squared = rMax * rMax;
0224     double Mx = M1.x();
0225     double My = M1.y();
0226 
0227     math::XYZVector B(0, 0, 0);
0228     math::XYZVector C(0, 0, 0);
0229 
0230     if (rMax_squared * 4. > Mx * Mx + My * My) {
0231       ////////////////////////
0232       // Cleaning P1 points:::
0233       ////////////////////////
0234 
0235       //Cx, Cy = Coordinates of circle center
0236       //C_=line that contains the circle center
0237 
0238       //C_=k*x+d
0239       double k = -Mx / My;
0240       double d = My / 2. - k * Mx / 2.;
0241 
0242 #ifdef mydebug_knuenz
0243       std::cout << "k" << k << std::endl;
0244       std::cout << "d" << d << std::endl;
0245 #endif
0246 
0247       //Cx1,2 and Cy1,2 are the two points that have a distance of rMax to 0,0
0248       double CsolutionPart1 = -2 * k * d;
0249       double CsolutionPart2 = std::sqrt(4 * k * k * d * d - 4 * (1 + k * k) * (d * d - rMax_squared));
0250       double CsolutionPart3 = 2 * (1 + k * k);
0251       double Cx1 = (CsolutionPart1 + CsolutionPart2) / CsolutionPart3;
0252       double Cx2 = (CsolutionPart1 - CsolutionPart2) / CsolutionPart3;
0253       double Cy1 = k * Cx1 + d;
0254       double Cy2 = k * Cx2 + d;
0255 
0256       // Decide between solutions: phi(C) > phi(P)
0257       double Cx, Cy;
0258       math::XYZVector C1(Cx1, Cy1, 0);
0259       if (C1.x() * M1.y() - C1.y() * M1.x() < 0) {
0260         Cx = Cx1;
0261         Cy = Cy1;
0262       } else {
0263         Cx = Cx2;
0264         Cy = Cy2;
0265       }
0266       C.SetXYZ(Cx, Cy, 0);
0267 
0268 #ifdef mydebug_knuenz
0269       std::cout << "Cx1" << Cx1 << std::endl;
0270       std::cout << "Cx2" << Cx2 << std::endl;
0271       std::cout << "Cy1" << Cy1 << std::endl;
0272       std::cout << "Cy2" << Cy2 << std::endl;
0273       std::cout << "Cx" << Cx << std::endl;
0274       std::cout << "Cy" << Cy << std::endl;
0275 #endif
0276 
0277       // Find Tangent at 0,0 to minPtCircle and point (Bx,By) on the first layer which bisects the allowed angle
0278       k = -Cx / Cy;
0279       d = 0;
0280       double Bx1 = std::sqrt(Mx * Mx + My * My / (1 + k * k));
0281       double Bx2 = -Bx1;
0282       double By1 = k * Bx1 + d;
0283       double By2 = k * Bx2 + d;
0284 
0285 #ifdef mydebug_knuenz
0286       std::cout << "k" << k << std::endl;
0287       std::cout << "d" << d << std::endl;
0288 #endif
0289 
0290       // Decide between solutions: phi(B) < phi(P)
0291       double Bx, By;
0292       math::XYZVector B1(Bx1, By1, 0);
0293       if (M1.x() * B1.y() - M1.y() * B1.x() < 0) {
0294         Bx = Bx1;
0295         By = By1;
0296       } else {
0297         Bx = Bx2;
0298         By = By2;
0299       }
0300       B.SetXYZ(Bx, By, 0);
0301 
0302 #ifdef mydebug_knuenz
0303       std::cout << "Bx1" << Bx1 << std::endl;
0304       std::cout << "Bx2" << Bx2 << std::endl;
0305       std::cout << "By1" << By1 << std::endl;
0306       std::cout << "By2" << By2 << std::endl;
0307       std::cout << "Bx" << Bx << std::endl;
0308       std::cout << "By" << By << std::endl;
0309 #endif
0310 
0311       double DeltaPhiMaxM1P1 = DeltaPhiManual(M1, B) * 2;
0312 
0313 #ifdef mydebug_knuenz
0314       std::cout << "DeltaPhiMaxM1P1 " << DeltaPhiMaxM1P1 << std::endl;
0315       std::cout << "M1.DeltaPhi(P1) " << DeltaPhiManual(M1, P1) << std::endl;
0316       std::cout << "rho P1: " << std::sqrt(P1.x() * P1.x() + P1.y() * P1.y()) << "phi P1: " << P1.Phi() << std::endl;
0317       std::cout << "rho M1: " << std::sqrt(M1.x() * M1.x() + M1.y() * M1.y()) << "phi M1: " << M1.Phi() << std::endl;
0318 #endif
0319 
0320       //Finally Cut on DeltaPhi of P1 and M1
0321 
0322       double tol_DeltaPhiMaxM1P1 = DeltaPhiMaxM1P1 * toleranceFactorOnDeltaPhiCuts;
0323       double DeltaPhiManualM1P1 = DeltaPhiManual(M1, P1);
0324 
0325       if (DeltaPhiManualM1P1 > DeltaPhiMaxM1P1 + tol_DeltaPhiMaxM1P1 || DeltaPhiManualM1P1 < 0 - tol_DeltaPhiMaxM1P1) {
0326         return nullptr;
0327       }
0328 
0329     }  //if rMax > rLayerM1
0330 
0331     ////////////////////////
0332     // Cleaning M2 points:::
0333     ////////////////////////
0334 
0335     //  if(B.DeltaPhi(P1)>0){//normal algo (with minPt circle)
0336 
0337     double rM2_squared = M2.x() * M2.x() + M2.y() * M2.y();
0338     if (rMax_squared * 4. > rM2_squared) {  //if minPt circle is smaller than 2*M2-layer radius, algo makes no sense
0339 
0340       //Chordales equation (line containing the two intersection points of the two circles)
0341       double k = -C.x() / C.y();
0342       double d = (rM2_squared - rMax_squared + C.x() * C.x() + C.y() * C.y()) / (2 * C.y());
0343 
0344       double M2solutionPart1 = -2 * k * d;
0345       double M2solutionPart2 = std::sqrt(4 * k * k * d * d - 4 * (1 + k * k) * (d * d - rM2_squared));
0346       double M2solutionPart3 = 2 + 2 * k * k;
0347       double M2xMax1 = (M2solutionPart1 + M2solutionPart2) / M2solutionPart3;
0348       double M2xMax2 = (M2solutionPart1 - M2solutionPart2) / M2solutionPart3;
0349       double M2yMax1 = k * M2xMax1 + d;
0350       double M2yMax2 = k * M2xMax2 + d;
0351 
0352       //double M2xMax,M2yMax;
0353       math::XYZVector M2MaxVec1(M2xMax1, M2yMax1, 0);
0354       math::XYZVector M2MaxVec2(M2xMax2, M2yMax2, 0);
0355       math::XYZVector M2MaxVec(0, 0, 0);
0356       if (M2MaxVec1.x() * M2MaxVec2.y() - M2MaxVec1.y() * M2MaxVec2.x() < 0) {
0357         M2MaxVec.SetXYZ(M2xMax2, M2yMax2, 0);
0358       } else {
0359         M2MaxVec.SetXYZ(M2xMax1, M2yMax1, 0);
0360       }
0361 
0362       double DeltaPhiMaxM2 = DeltaPhiManual(M2MaxVec, M1);
0363 
0364 #ifdef mydebug_knuenz
0365       std::cout << "C.x() " << C.x() << std::endl;
0366       std::cout << "C.y() " << C.y() << std::endl;
0367       std::cout << "M1.x() " << M1.x() << std::endl;
0368       std::cout << "M1.y() " << M1.y() << std::endl;
0369       std::cout << "M2.x() " << M2.x() << std::endl;
0370       std::cout << "M2.y() " << M2.y() << std::endl;
0371       std::cout << "k " << k << std::endl;
0372       std::cout << "d " << d << std::endl;
0373       std::cout << "M2xMax1 " << M2xMax1 << std::endl;
0374       std::cout << "M2xMax2 " << M2xMax2 << std::endl;
0375       std::cout << "M2yMax1 " << M2yMax1 << std::endl;
0376       std::cout << "M2yMax2 " << M2yMax2 << std::endl;
0377       std::cout << "M2xMax " << M2MaxVec.x() << std::endl;
0378       std::cout << "M2yMax " << M2MaxVec.y() << std::endl;
0379       std::cout << "rM2_squared " << rM2_squared << std::endl;
0380       std::cout << "rMax " << rMax << std::endl;
0381       std::cout << "DeltaPhiMaxM2 " << DeltaPhiMaxM2 << std::endl;
0382 #endif
0383 
0384       double tol_DeltaPhiMaxM2 = DeltaPhiMaxM2 * toleranceFactorOnDeltaPhiCuts;
0385       double DeltaPhiManualM2M1 = DeltaPhiManual(M2, M1);
0386 
0387       if (DeltaPhiManualM2M1 > DeltaPhiMaxM2 + tol_DeltaPhiMaxM2 || DeltaPhiManualM2M1 < 0 - tol_DeltaPhiMaxM2) {
0388         return nullptr;
0389       }
0390 
0391       //Using the lazy solution for P2: DeltaPhiMaxP2=DeltaPhiMaxM2
0392       double DeltaPhiManualP1P2 = DeltaPhiManual(P1, P2);
0393       if (DeltaPhiManualP1P2 > DeltaPhiMaxM2 + tol_DeltaPhiMaxM2 || DeltaPhiManualP1P2 < 0 - tol_DeltaPhiMaxM2) {
0394         return nullptr;
0395       }
0396     }
0397   }
0398 
0399   //  }
0400 
0401   //  if(B.DeltaPhi(P1)<0){//different algo (without minPt circle)
0402 
0403   //  }
0404 
0405   // End of pre-seed cleaning
0406 
0407   Conv4HitsReco2 quad(vPhotVertex, h1, h2, h3, h4);
0408   quad.SetMaxNumberOfIterations(100);
0409 #ifdef mydebug_sguazz
0410   quad.Dump();
0411 #endif
0412   math::XYZVector candVtx;
0413   double candPtPlus, candPtMinus;
0414   //double rPlus, rMinus;
0415   int nite = quad.ConversionCandidate(candVtx, candPtPlus, candPtMinus);
0416 
0417   if (!(nite && abs(nite) < 25 && nite != -1000 && nite != -2000))
0418     return nullptr;
0419 
0420     //  math::XYZVector plusCenter = quad.GetPlusCenter(rPlus);
0421     //  math::XYZVector minusCenter = quad.GetMinusCenter(rMinus);
0422 
0423 #ifdef mydebug_sguazz
0424   std::cout << " >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>"
0425             << "\n";
0426   std::cout << " >>>>>>>>>>> Conv Cand: "
0427             << " Vertex X: " << candVtx.X() << " [cm] Y: " << candVtx.Y() << " [cm] pt+: " << candPtPlus
0428             << " [GeV] pt-: " << candPtMinus << " [GeV]; #its: " << nite << "\n";
0429 #endif
0430 
0431   //Add here cuts
0432   double minLegPt = CleaningMinLegPt;
0433   double maxRadialDistance = CleaningmaxRadialDistance;
0434 
0435   //
0436   // Cut on leg's transverse momenta
0437   if (candPtPlus < minLegPt)
0438     return nullptr;
0439   if (candPtMinus < minLegPt)
0440     return nullptr;
0441   //
0442   if (candPtPlus > maxLegPt)
0443     return nullptr;
0444   if (candPtMinus > maxLegPt)
0445     return nullptr;
0446   //
0447   // Cut on radial distance between estimated conversion vertex and inner hits
0448   double cr = std::sqrt(candVtx.Perp2());
0449   double maxr2 = (maxRadialDistance + cr);
0450   maxr2 *= maxr2;
0451   if (h2.Perp2() > maxr2)
0452     return nullptr;
0453   if (h3.Perp2() > maxr2)
0454     return nullptr;
0455 
0456   // At this point implement cleaning cuts after building the seed
0457 
0458   //ClusterShapeFilter_knuenz:::
0459   const auto& bfield = es.getData(theBfieldToken);
0460   float nomField = bfield.nominalValue();
0461 
0462   if (ClusterShapeFiltering) {
0463     if (theComparitor)
0464       theComparitor->init(ev, es);
0465 
0466     GlobalTrajectoryParameters pkine;
0467     GlobalTrajectoryParameters mkine;
0468 
0469     // SeedingHitSet::ConstRecHitPointer ptth1 = phits[0];  // (never used???)
0470     SeedingHitSet::ConstRecHitPointer ptth2 = phits[1];
0471     SeedingHitSet::ConstRecHitPointer mtth1 = mhits[0];
0472     SeedingHitSet::ConstRecHitPointer mtth2 = mhits[1];
0473 
0474     GlobalPoint vertexPos(candVtx.x(), candVtx.y(), candVtx.z());
0475 
0476     float ptMinReg = 0.1;
0477     GlobalTrackingRegion region(ptMinReg, vertexPos, 0, 0, true);
0478 
0479     FastHelix phelix(ptth2->globalPosition(), mtth1->globalPosition(), vertexPos, nomField, &bfield, vertexPos);
0480     pkine = phelix.stateAtVertex();
0481     FastHelix mhelix(mtth2->globalPosition(), mtth1->globalPosition(), vertexPos, nomField, &bfield, vertexPos);
0482     mkine = mhelix.stateAtVertex();
0483 
0484     if (theComparitor && !theComparitor->compatible(phits, pkine, phelix)) {
0485       return nullptr;
0486     }
0487     if (theComparitor && !theComparitor->compatible(mhits, mkine, mhelix)) {
0488       return nullptr;
0489     }
0490   }
0491 
0492   //Do a very simple fit to estimate the slope
0493   double quadPhotCotTheta = 0.;
0494   double quadZ0 = 0.;
0495   simpleGetSlope(ptth2, ptth1, mtth1, mtth2, region, quadPhotCotTheta, quadZ0);
0496 
0497   double quadPhotPhi = (candVtx - vPhotVertex).Phi();
0498 
0499   math::XYZVector fittedPrimaryVertex(vgPhotVertex.x(), vgPhotVertex.y(), quadZ0);
0500 
0501   candVtx.SetZ(std::sqrt(candVtx.Perp2()) * quadPhotCotTheta + quadZ0);
0502   GlobalPoint convVtxGlobalPoint(candVtx.X(), candVtx.Y(), candVtx.Z());
0503 
0504   //
0505   // Comparing new quad with old quad
0506   //
0507   //Arbitration
0508 
0509   Quad thisQuad;
0510   thisQuad.x = candVtx.X();
0511   thisQuad.y = candVtx.Y();
0512   thisQuad.z = candVtx.Z();
0513   thisQuad.ptPlus = candPtPlus;
0514   thisQuad.ptMinus = candPtMinus;
0515   thisQuad.cot = quadPhotCotTheta;
0516   if (similarQuadExist(thisQuad, quadV) && applyArbitration)
0517     return nullptr;
0518 
0519   // not able to get the mag field... doing the dirty way
0520   //
0521   // Plus
0522   FastHelix helixPlus(
0523       ptth2->globalPosition(), ptth1->globalPosition(), convVtxGlobalPoint, nomField, &bfield, convVtxGlobalPoint);
0524   GlobalTrajectoryParameters kinePlus = helixPlus.stateAtVertex();
0525   kinePlus = GlobalTrajectoryParameters(
0526       convVtxGlobalPoint,
0527       GlobalVector(candPtPlus * cos(quadPhotPhi), candPtPlus * sin(quadPhotPhi), candPtPlus * quadPhotCotTheta),
0528       1,  //1
0529       &kinePlus.magneticField());
0530 
0531   //
0532   // Minus
0533   FastHelix helixMinus(
0534       mtth2->globalPosition(), mtth1->globalPosition(), convVtxGlobalPoint, nomField, &bfield, convVtxGlobalPoint);
0535   GlobalTrajectoryParameters kineMinus = helixMinus.stateAtVertex();
0536   kineMinus = GlobalTrajectoryParameters(
0537       convVtxGlobalPoint,
0538       GlobalVector(candPtMinus * cos(quadPhotPhi), candPtMinus * sin(quadPhotPhi), candPtMinus * quadPhotCotTheta),
0539       -1,  //-1
0540       &kineMinus.magneticField());
0541 
0542   float sinThetaPlus = sin(kinePlus.momentum().theta());
0543   float sinThetaMinus = sin(kineMinus.momentum().theta());
0544   float ptmin = region.ptMin();
0545   //vertexBounds da region
0546   GlobalVector vertexBounds(region.originRBound(), region.originRBound(), region.originZBound());
0547 
0548   CurvilinearTrajectoryError errorPlus = initialError(vertexBounds, ptmin, sinThetaPlus);
0549   CurvilinearTrajectoryError errorMinus = initialError(vertexBounds, ptmin, sinThetaMinus);
0550   FreeTrajectoryState ftsPlus(kinePlus, errorPlus);
0551   FreeTrajectoryState ftsMinus(kineMinus, errorMinus);
0552 
0553   //FIXME: here probably you want to go in parallel with phits and mhits
0554   //NB: the seedCollection is filled (the important thing) the return of the last TrajectorySeed is not used, but is needed
0555   //to maintain the inheritance
0556 
0557 #ifdef quadDispLine
0558   double vError = region.originZBound();
0559   if (vError > 15.)
0560     vError = 1.;
0561   std::cout << "QuadDispLine " << vgPhotVertex.x() << " " << vgPhotVertex.y() << " " << vgPhotVertex.z() << " "
0562             << vError << " " << vHit[0].x() << " " << vHit[0].y() << " " << vHit[0].z() << " "
0563             << sqrt(getSqrEffectiveErrorOnZ(ptth2, region)) << " " << vHit[1].x() << " " << vHit[1].y() << " "
0564             << vHit[1].z() << " " << sqrt(getSqrEffectiveErrorOnZ(ptth1, region)) << " " << vHit[2].x() << " "
0565             << vHit[2].y() << " " << vHit[2].z() << " " << sqrt(getSqrEffectiveErrorOnZ(mtth1, region)) << " "
0566             << vHit[3].x() << " " << vHit[3].y() << " " << vHit[3].z() << " "
0567             << sqrt(getSqrEffectiveErrorOnZ(mtth2, region)) << " " << candVtx.X() << " " << candVtx.Y() << " "
0568             << candVtx.Z() << " " << fittedPrimaryVertex.X() << " " << fittedPrimaryVertex.Y() << " "
0569             << fittedPrimaryVertex.Z()
0570             << " "
0571             //        << candPtPlus  << " " << rPlus << " " << plusCenter.X() << " " << plusCenter.Y() << " "
0572             //        << candPtMinus << " " << rMinus << " " << minusCenter.X() << " " << minusCenter.Y() << " "
0573             << nite << " " << chi2 << "\n";
0574 #endif
0575 #ifdef mydebug_sguazz
0576   std::cout << " >>>>> Hit quad fast reco done >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>"
0577             << "\n";
0578   uint32_t detid;
0579   std::cout << "[SeedForPhotonConversionFromQuadruplets]\n ptth1 ";
0580   detid = ptth1->geographicalId().rawId();
0581   //    po.print(std::cout , detid );
0582   std::cout << " \t " << detid << " " << ptth1->localPosition() << " " << ptth1->globalPosition();
0583   detid = ptth2->geographicalId().rawId();
0584   std::cout << " \n\t ptth2 ";
0585   //    po.print(std::cout , detid );
0586   std::cout << " \t " << detid << " " << ptth2->localPosition() << " " << ptth2->globalPosition() << "\nhelix momentum "
0587             << kinePlus.momentum() << " pt " << kinePlus.momentum().perp() << " radius "
0588             << 1 / kinePlus.transverseCurvature() << " q " << kinePlus.charge();
0589   std::cout << " \n\t mtth1 ";
0590   detid = mtth1->geographicalId().rawId();
0591   std::cout << " \t " << detid << " " << mtth1->localPosition() << " " << mtth1->globalPosition();
0592   std::cout << " \n\t mtth2 ";
0593   detid = mtth2->geographicalId().rawId();
0594   //    po.print(std::cout , detid );
0595   std::cout << " \t " << detid << " " << mtth2->localPosition() << " " << mtth2->globalPosition() << "\nhelix momentum "
0596             << kineMinus.momentum() << " pt " << kineMinus.momentum().perp() << " radius "
0597             << 1 / kineMinus.transverseCurvature() << " q " << kineMinus.charge();
0598   std::cout << "\n <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<"
0599             << "\n";
0600 #endif
0601 
0602   // get cloner
0603   auto builder = static_cast<TkTransientTrackingRecHitBuilder const*>(&es.getData(theTTRHBuilderToken));
0604   cloner = (*builder).cloner();
0605 
0606   bool buildSeedBoolPos = buildSeedBool(seedCollection, phits, ftsPlus, es, applydzCAcut, region, dzcut);
0607   bool buildSeedBoolNeg = buildSeedBool(seedCollection, mhits, ftsMinus, es, applydzCAcut, region, dzcut);
0608 
0609   if (buildSeedBoolPos && buildSeedBoolNeg) {
0610     buildSeed(seedCollection, phits, ftsPlus, es, false, region);
0611     buildSeed(seedCollection, mhits, ftsMinus, es, false, region);
0612   }
0613 
0614   return nullptr;
0615 }
0616 
0617 GlobalTrajectoryParameters SeedForPhotonConversionFromQuadruplets::initialKinematic(const SeedingHitSet& hits,
0618                                                                                     const GlobalPoint& vertexPos,
0619                                                                                     const edm::EventSetup& es,
0620                                                                                     const float cotTheta) const {
0621   GlobalTrajectoryParameters kine;
0622 
0623   SeedingHitSet::ConstRecHitPointer tth1 = hits[0];
0624   SeedingHitSet::ConstRecHitPointer tth2 = hits[1];
0625 
0626   const auto& bfield = es.getData(theBfieldToken);
0627   float nomField = bfield.nominalValue();
0628   bool isBOFF = (0 == nomField);
0629 
0630   FastHelix helix(tth2->globalPosition(), tth1->globalPosition(), vertexPos, nomField, &bfield, vertexPos);
0631   kine = helix.stateAtVertex();
0632 
0633   //force the pz/pt equal to the measured one
0634   if (fabs(cotTheta) < cotTheta_Max)
0635     kine = GlobalTrajectoryParameters(
0636         kine.position(),
0637         GlobalVector(kine.momentum().x(), kine.momentum().y(), kine.momentum().perp() * cotTheta),
0638         kine.charge(),
0639         &kine.magneticField());
0640   else
0641     kine = GlobalTrajectoryParameters(
0642         kine.position(),
0643         GlobalVector(kine.momentum().x(), kine.momentum().y(), kine.momentum().perp() * cotTheta_Max),
0644         kine.charge(),
0645         &kine.magneticField());
0646 
0647 #ifdef mydebug_seed
0648   uint32_t detid;
0649   (*pss) << "[SeedForPhotonConversionFromQuadruplets] initialKinematic tth1 ";
0650   detid = tth1->geographicalId().rawId();
0651   po.print(*pss, detid);
0652   (*pss) << " \t " << detid << " " << tth1->localPosition() << " " << tth1->globalPosition();
0653   detid = tth2->geographicalId().rawId();
0654   (*pss) << " \n\t tth2 ";
0655   po.print(*pss, detid);
0656   (*pss) << " \t " << detid << " " << tth2->localPosition() << " " << tth2->globalPosition() << "\nhelix momentum "
0657          << kine.momentum() << " pt " << kine.momentum().perp() << " radius " << 1 / kine.transverseCurvature();
0658 #endif
0659 
0660   if (isBOFF && (theBOFFMomentum > 0)) {
0661     kine =
0662         GlobalTrajectoryParameters(kine.position(), kine.momentum().unit() * theBOFFMomentum, kine.charge(), &bfield);
0663   }
0664   return kine;
0665 }
0666 
0667 CurvilinearTrajectoryError SeedForPhotonConversionFromQuadruplets::initialError(const GlobalVector& vertexBounds,
0668                                                                                 float ptMin,
0669                                                                                 float sinTheta) const {
0670   // Set initial uncertainty on track parameters, using only P.V. constraint and no hit
0671   // information.
0672   GlobalError vertexErr(sqr(vertexBounds.x()), 0, sqr(vertexBounds.y()), 0, 0, sqr(vertexBounds.z()));
0673 
0674   AlgebraicSymMatrix55 C = ROOT::Math::SMatrixIdentity();
0675 
0676   // FIXME: minC00. Prevent apriori uncertainty in 1/P from being too small,
0677   // to avoid instabilities.
0678   // N.B. This parameter needs optimising ...
0679   float sin2th = sqr(sinTheta);
0680   float minC00 = 1.0;
0681   C[0][0] = std::max(sin2th / sqr(ptMin), minC00);
0682   float zErr = vertexErr.czz();
0683   float transverseErr = vertexErr.cxx();  // assume equal cxx cyy
0684   C[3][3] = transverseErr;
0685   C[4][4] = zErr * sin2th + transverseErr * (1 - sin2th);
0686 
0687   return CurvilinearTrajectoryError(C);
0688 }
0689 
0690 const TrajectorySeed* SeedForPhotonConversionFromQuadruplets::buildSeed(TrajectorySeedCollection& seedCollection,
0691                                                                         const SeedingHitSet& hits,
0692                                                                         const FreeTrajectoryState& fts,
0693                                                                         const edm::EventSetup& es,
0694                                                                         bool apply_dzCut,
0695                                                                         const TrackingRegion& region) const {
0696   // get tracker
0697   const auto& tracker = es.getData(theTrackerToken);
0698 
0699   // get propagator
0700   const Propagator* propagator = &es.getData(thePropagatorToken);
0701 
0702   // get updator
0703   KFUpdator updator;
0704 
0705   // Now update initial state track using information from seed hits.
0706 
0707   TrajectoryStateOnSurface updatedState;
0708   edm::OwnVector<TrackingRecHit> seedHits;
0709 
0710   const TrackingRecHit* hit = nullptr;
0711   for (unsigned int iHit = 0; iHit < hits.size() && iHit < 2; iHit++) {
0712     hit = hits[iHit];
0713     TrajectoryStateOnSurface state =
0714         (iHit == 0) ? propagator->propagate(fts, tracker.idToDet(hit->geographicalId())->surface())
0715                     : propagator->propagate(updatedState, tracker.idToDet(hit->geographicalId())->surface());
0716 
0717     SeedingHitSet::ConstRecHitPointer tth = hits[iHit];
0718 
0719     SeedingHitSet::RecHitPointer newtth = refitHit(tth, state);
0720 
0721     updatedState = updator.update(state, *newtth);
0722 
0723     seedHits.push_back(newtth);
0724 #ifdef mydebug_seed
0725     uint32_t detid = hit->geographicalId().rawId();
0726     (*pss) << "\n[SeedForPhotonConversionFromQuadruplets] hit " << iHit;
0727     po.print(*pss, detid);
0728     (*pss) << " " << detid << "\t lp " << hit->localPosition() << " tth " << tth->localPosition() << " newtth "
0729            << newtth->localPosition() << " state " << state.globalMomentum().perp();
0730 #endif
0731   }
0732 
0733   if (!hit)
0734     return nullptr;
0735 
0736   PTrajectoryStateOnDet const& PTraj =
0737       trajectoryStateTransform::persistentState(updatedState, hit->geographicalId().rawId());
0738 
0739   seedCollection.push_back(TrajectorySeed(PTraj, seedHits, alongMomentum));
0740   return &seedCollection.back();
0741 }
0742 
0743 bool SeedForPhotonConversionFromQuadruplets::buildSeedBool(TrajectorySeedCollection& seedCollection,
0744                                                            const SeedingHitSet& hits,
0745                                                            const FreeTrajectoryState& fts,
0746                                                            const edm::EventSetup& es,
0747                                                            bool apply_dzCut,
0748                                                            const TrackingRegion& region,
0749                                                            double dzcut) const {
0750   // get tracker
0751   const auto& tracker = es.getData(theTrackerToken);
0752 
0753   // get propagator
0754   const Propagator* propagator = &es.getData(thePropagatorToken);
0755 
0756   // get updator
0757   KFUpdator updator;
0758 
0759   // Now update initial state track using information from seed hits.
0760 
0761   TrajectoryStateOnSurface updatedState;
0762   edm::OwnVector<TrackingRecHit> seedHits;
0763 
0764   const TrackingRecHit* hit = nullptr;
0765   for (unsigned int iHit = 0; iHit < hits.size() && iHit < 2; iHit++) {
0766     hit = hits[iHit]->hit();
0767     TrajectoryStateOnSurface state =
0768         (iHit == 0) ? propagator->propagate(fts, tracker.idToDet(hit->geographicalId())->surface())
0769                     : propagator->propagate(updatedState, tracker.idToDet(hit->geographicalId())->surface());
0770     if (!state.isValid()) {
0771       return false;
0772     }
0773 
0774     SeedingHitSet::ConstRecHitPointer tth = hits[iHit];
0775 
0776     SeedingHitSet::RecHitPointer newtth = refitHit(tth, state);
0777 
0778     if (!checkHit(state, newtth, es)) {
0779       return false;
0780     }
0781 
0782     updatedState = updator.update(state, *newtth);
0783     if (!updatedState.isValid()) {
0784       return false;
0785     }
0786 
0787     seedHits.push_back(newtth->hit()->clone());
0788   }
0789 
0790   if (apply_dzCut) {
0791     /// Implement here the dz cut:::
0792 
0793     double zCA;
0794 
0795     math::XYZVector EstMomGam(
0796         updatedState.globalMomentum().x(), updatedState.globalMomentum().y(), updatedState.globalMomentum().z());
0797     math::XYZVector EstPosGam(
0798         updatedState.globalPosition().x(), updatedState.globalPosition().y(), updatedState.globalPosition().z());
0799 
0800     double EstMomGamLength =
0801         std::sqrt(EstMomGam.x() * EstMomGam.x() + EstMomGam.y() * EstMomGam.y() + EstMomGam.z() * EstMomGam.z());
0802     math::XYZVector EstMomGamNorm(
0803         EstMomGam.x() / EstMomGamLength, EstMomGam.y() / EstMomGamLength, EstMomGam.z() / EstMomGamLength);
0804 
0805     //Calculate dz of point of closest approach of the two lines (WA approach) -> cut on dz
0806 
0807     const GlobalPoint EstPosGamGlobalPoint(
0808         updatedState.globalPosition().x(), updatedState.globalPosition().y(), updatedState.globalPosition().z());
0809     const GlobalVector EstMomGamGlobalVector(
0810         updatedState.globalMomentum().x(), updatedState.globalMomentum().y(), updatedState.globalMomentum().z());
0811 
0812     const MagneticField* magField = &es.getData(theBfieldToken);
0813     TrackCharge qCharge = 0;
0814 
0815     const GlobalTrajectoryParameters myGlobalTrajectoryParameter(
0816         EstPosGamGlobalPoint, EstMomGamGlobalVector, qCharge, magField);
0817 
0818     AlgebraicSymMatrix66 aCovarianceMatrix;
0819 
0820     for (int i = 0; i < 6; ++i)
0821       for (int j = 0; j < 6; ++j)
0822         aCovarianceMatrix(i, j) = 1e-4;
0823 
0824     CartesianTrajectoryError myCartesianError(aCovarianceMatrix);
0825 
0826     const FreeTrajectoryState stateForProjectionToBeamLine(myGlobalTrajectoryParameter, myCartesianError);
0827 
0828     const GlobalPoint BeamSpotGlobalPoint(0, 0, 0);
0829 
0830     const reco::BeamSpot::Point BeamSpotPoint(region.origin().x(), region.origin().y(), region.origin().z());
0831 
0832     TSCBLBuilderNoMaterial tscblBuilder;
0833 
0834     double CovMatEntry = 0.;
0835     reco::BeamSpot::CovarianceMatrix cov;
0836     for (int i = 0; i < 3; ++i) {
0837       cov(i, i) = CovMatEntry;
0838     }
0839     reco::BeamSpot::BeamType BeamType_ = reco::BeamSpot::Unknown;
0840 
0841     reco::BeamSpot myBeamSpot(BeamSpotPoint, 0., 0., 0., 0., cov, BeamType_);
0842 
0843     TrajectoryStateClosestToBeamLine tscbl = tscblBuilder(stateForProjectionToBeamLine, myBeamSpot);
0844     if (tscbl.isValid() == false) {
0845       zCA = 0;
0846     } else {
0847       GlobalPoint v = tscbl.trackStateAtPCA().position();  // Position of closest approach to BS
0848       zCA = v.z();
0849     }
0850 
0851     /*  //Calculate dz of point of closest approach of the two lines -> cut on dz
0852 
0853   double newX,newY,newR;
0854   double Rbuff=std::sqrt(EstPosGam.x()*EstPosGam.x()+EstPosGam.y()*EstPosGam.y());
0855   double deltas,s,sbuff;
0856   double rMin=1e9;
0857 
0858   double InitXp=EstPosGam.x()+1*EstMomGamNorm.x();
0859   double InitXm=EstPosGam.x()-1*EstMomGamNorm.x();
0860   double InitYp=EstPosGam.y()+1*EstMomGamNorm.y();
0861   double InitYm=EstPosGam.y()-1*EstMomGamNorm.y();
0862 
0863   if(InitXp*InitXp+InitYp*InitYp < InitXm*InitXm+InitYm*InitYm)  {s=5; deltas=5;}
0864   else {s=-5; deltas=-5;}
0865 
0866   int nTurns=0;
0867   int nIterZ=0;
0868   for (int i_dz=0;i_dz<1000;i_dz++){
0869   newX=EstPosGam.x()+s*EstMomGamNorm.x();
0870   newY=EstPosGam.y()+s*EstMomGamNorm.y();
0871   newR=std::sqrt(newX*newX+newY*newY);
0872   if(newR>Rbuff) {deltas=-1*deltas/10;nTurns++;}
0873   else {Rbuff=newR;}
0874   if(newR<rMin) {rMin=newR; sbuff=s;}
0875   s=s+deltas;
0876   nIterZ++;
0877   if(nTurns>1) break;
0878   }
0879 
0880   zCA=EstPosGam.z()+sbuff*EstMomGamNorm.z();
0881 */
0882 
0883 #ifdef mydebug_knuenz
0884     std::cout << "zCA: " << zCA << std::endl;
0885 #endif
0886 
0887     if (std::fabs(zCA) > dzcut)
0888       return false;
0889   }
0890 
0891   return true;
0892 }
0893 
0894 SeedingHitSet::RecHitPointer SeedForPhotonConversionFromQuadruplets::refitHit(
0895     SeedingHitSet::ConstRecHitPointer hit, const TrajectoryStateOnSurface& state) const {
0896   //const TransientTrackingRecHit* a= hit.get();
0897   //return const_cast<TransientTrackingRecHit*> (a);
0898   //This was modified otherwise the rechit will have just the local x component and local y=0
0899   // To understand how to modify for pixels
0900 
0901   //const TSiStripRecHit2DLocalPos* b = dynamic_cast<const TSiStripRecHit2DLocalPos*>(a);
0902   //return const_cast<TSiStripRecHit2DLocalPos*>(b);
0903   return (SeedingHitSet::RecHitPointer)(cloner(*hit, state));
0904 }
0905 
0906 //
0907 // Below: stupid utils method by sguazz
0908 //
0909 //
0910 void SeedForPhotonConversionFromQuadruplets::stupidPrint(std::string s, float* d) {
0911   (*pss) << "\n" << s << "\t";
0912   for (size_t i = 0; i < 2; ++i)
0913     (*pss) << std::setw(60) << d[i] << std::setw(1) << " | ";
0914 }
0915 
0916 void SeedForPhotonConversionFromQuadruplets::stupidPrint(std::string s, double* d) {
0917   (*pss) << "\n" << s << "\t";
0918   for (size_t i = 0; i < 2; ++i)
0919     (*pss) << std::setw(60) << d[i] << std::setw(1) << " | ";
0920 }
0921 
0922 void SeedForPhotonConversionFromQuadruplets::stupidPrint(const char* s, GlobalPoint* d) {
0923   (*pss) << "\n" << s << "\t";
0924   for (size_t i = 0; i < 2; ++i)
0925     (*pss) << std::setw(20) << d[i] << " r " << d[i].perp() << " phi " << d[i].phi() << " | ";
0926 }
0927 
0928 void SeedForPhotonConversionFromQuadruplets::stupidPrint(const char* s, GlobalPoint* d, int n) {
0929   (*pss) << "\n" << s << "\n";
0930   for (int i = 0; i < n; ++i)
0931     (*pss) << std::setw(20) << d[i] << " r " << d[i].perp() << " phi " << d[i].phi() << "\n";
0932 }
0933 
0934 #include "DataFormats/Math/interface/deltaPhi.h"
0935 
0936 void SeedForPhotonConversionFromQuadruplets::bubbleSortVsPhi(GlobalPoint arr[], int n, GlobalPoint vtx) {
0937   bool swapped = true;
0938   int j = 0;
0939   GlobalPoint tmp;
0940   while (swapped) {
0941     swapped = false;
0942     j++;
0943     for (int i = 0; i < n - j; i++) {
0944       if (reco::deltaPhi((arr[i] - vtx).barePhi(), (arr[i + 1] - vtx).barePhi()) > 0.) {
0945         tmp = arr[i];
0946         arr[i] = arr[i + 1];
0947         arr[i + 1] = tmp;
0948         swapped = true;
0949       }
0950     }
0951   }
0952 }
0953 
0954 void SeedForPhotonConversionFromQuadruplets::bubbleReverseSortVsPhi(GlobalPoint arr[], int n, GlobalPoint vtx) {
0955   bool swapped = true;
0956   int j = 0;
0957   GlobalPoint tmp;
0958   while (swapped) {
0959     swapped = false;
0960     j++;
0961     for (int i = 0; i < n - j; i++) {
0962       if (reco::deltaPhi((arr[i] - vtx).barePhi(), (arr[i + 1] - vtx).barePhi()) < 0.) {
0963         tmp = arr[i];
0964         arr[i] = arr[i + 1];
0965         arr[i + 1] = tmp;
0966         swapped = true;
0967       }
0968     }
0969   }
0970 }
0971 
0972 double SeedForPhotonConversionFromQuadruplets::simpleGetSlope(const SeedingHitSet::ConstRecHitPointer& ohit,
0973                                                               const SeedingHitSet::ConstRecHitPointer& nohit,
0974                                                               const SeedingHitSet::ConstRecHitPointer& ihit,
0975                                                               const SeedingHitSet::ConstRecHitPointer& nihit,
0976                                                               const TrackingRegion& region,
0977                                                               double& cotTheta,
0978                                                               double& z0) {
0979   double x[5], y[5], e2y[5];
0980 
0981   //The fit is done filling x with r values, y with z values of the four hits and the vertex
0982   //
0983   //Hits
0984   x[0] = ohit->globalPosition().perp();
0985   y[0] = ohit->globalPosition().z();
0986   e2y[0] = getSqrEffectiveErrorOnZ(ohit, region);
0987   //
0988   x[1] = nohit->globalPosition().perp();
0989   y[1] = nohit->globalPosition().z();
0990   e2y[1] = getSqrEffectiveErrorOnZ(nohit, region);
0991   //
0992   x[2] = nihit->globalPosition().perp();
0993   y[2] = nihit->globalPosition().z();
0994   e2y[2] = getSqrEffectiveErrorOnZ(nihit, region);
0995   //
0996   x[3] = ihit->globalPosition().perp();
0997   y[3] = ihit->globalPosition().z();
0998   e2y[3] = getSqrEffectiveErrorOnZ(ihit, region);
0999   //
1000   //Vertex
1001   x[4] = region.origin().perp();
1002   y[4] = region.origin().z();
1003   double vError = region.originZBound();
1004   if (vError > 15.)
1005     vError = 1.;
1006   e2y[4] = sqr(vError);
1007 
1008   double e2z0;
1009   double chi2 = verySimpleFit(5, x, y, e2y, z0, e2z0, cotTheta);
1010 
1011   return chi2;
1012 }
1013 
1014 double SeedForPhotonConversionFromQuadruplets::verySimpleFit(
1015     int size, double* ax, double* ay, double* e2y, double& p0, double& e2p0, double& p1) {
1016   //#include "verySimpleFit.icc"
1017   return 0;
1018 }
1019 
1020 double SeedForPhotonConversionFromQuadruplets::getSqrEffectiveErrorOnZ(const SeedingHitSet::ConstRecHitPointer& hit,
1021                                                                        const TrackingRegion& region) {
1022   //
1023   //Fit-wise the effective error on Z is the sum in quadrature of the error on Z
1024   //and the error on R correctly projected by using hit-vertex direction
1025 
1026   double sqrProjFactor =
1027       sqr((hit->globalPosition().z() - region.origin().z()) / (hit->globalPosition().perp() - region.origin().perp()));
1028   return (hit->globalPositionError().czz() + sqrProjFactor * hit->globalPositionError().rerr(hit->globalPosition()));
1029 }
1030 
1031 bool SeedForPhotonConversionFromQuadruplets::similarQuadExist(Quad& thisQuad, std::vector<Quad>& quadV) {
1032   for (auto const& quad : quadV) {
1033     double dx = thisQuad.x - quad.x;
1034     double dy = thisQuad.y - quad.y;
1035     //double dz = abs(thisQuad.z-quad.z);
1036     //std::cout<<"thisQuad.x="<<thisQuad.x<<"  "<<"quad.x="<<quad.x<<"  "<<"thisQuad.y="<<thisQuad.y<<"  "<<"quad.y="<<quad.y<<"  "<<"thisQuad.z"<<thisQuad.z<<"  "<<"quad.z="<<quad.z<<"  "<<"thisQuad.ptPlus"<<thisQuad.ptPlus<<"  "<<"quad.ptPlus="<<quad.ptPlus<<"  "<<"thisQuad.ptMinus="<<thisQuad.ptMinus<<"  "<<"quad.ptMinus="<<quad.ptMinus<<"  "<<"thisQuad.cot="<<thisQuad.cot<<"  "<<"quad.cot="<<quad.cot<<std::endl; //ValDebug
1037     //std::cout<<"x1-x2="<<dx<<"y1-y2="<<dy<<"dx*dx+dy*dy="<<dx*dx+dy*dy<<std::endl; //ValDebug
1038     //std::cout<<"thisQuad.ptPlus-quad.ptPlus="<<thisQuad.ptPlus-quad.ptPlus<<"abs(thisQuad.ptPlus-quad.ptPlus)="<<abs(thisQuad.ptPlus-quad.ptPlus)<<std::endl; //ValDebug
1039     //std::cout <<sqrt(dx*dx+dy*dy)<<" <1? "<<dz<<" <3? "<<abs(thisQuad.ptPlus-quad.ptPlus)<<" <0.5? "<<abs(thisQuad.ptMinus-quad.ptMinus)<<" <0.5? "<<abs(thisQuad.cot-quad.cot)<<" <0.3? "<<std::endl; //ValDebug
1040     //if ( sqrt(dx*dx+dy*dy)<1.)
1041     //    std::cout <<sqrt(dx*dx+dy*dy)<<" <1? "<<dz<<" <3? "<<abs(thisQuad.ptPlus-quad.ptPlus)<<" <"<<0.5*quad.ptPlus<<"? "<<abs(thisQuad.ptMinus-quad.ptMinus)<<" <"<<0.5*quad.ptMinus<<"? "<<abs(thisQuad.cot-quad.cot)<<" <"<<0.3*quad.cot<<"? "<<std::endl; //ValDebug
1042     // ( sqrt(dx*dx+dy*dy)<1. &&
1043     //z<3. &&
1044     //bs(thisQuad.ptPlus-quad.ptPlus)<0.5*quad.ptPlus &&
1045     //bs(thisQuad.ptMinus-quad.ptMinus)<0.5*quad.ptMinus &&
1046     //bs(thisQuad.cot-quad.cot)<0.3*quad.cot
1047     //
1048     //  {
1049     //  //std::cout<<"Seed rejected due to arbitration"<<std::endl;
1050     //  return true;
1051     //  }
1052     if (sqrt(dx * dx + dy * dy) < 1. && fabs(thisQuad.ptPlus - quad.ptPlus) < 0.5 * quad.ptPlus &&
1053         fabs(thisQuad.ptMinus - quad.ptMinus) < 0.5 * quad.ptMinus) {
1054       //std::cout<<"Seed rejected due to arbitration"<<std::endl;
1055       return true;
1056     }
1057   }
1058   quadV.push_back(thisQuad);
1059   return false;
1060 }
1061 
1062 double SeedForPhotonConversionFromQuadruplets::DeltaPhiManual(const math::XYZVector& v1, const math::XYZVector& v2) {
1063   double kTWOPI = 2. * kPI_;
1064   double DeltaPhiMan = v1.Phi() - v2.Phi();
1065   while (DeltaPhiMan >= kPI_)
1066     DeltaPhiMan -= kTWOPI;
1067   while (DeltaPhiMan < -kPI_)
1068     DeltaPhiMan += kTWOPI;
1069 
1070   return DeltaPhiMan;
1071 }
[ 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!