CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
List of all members | Public Member Functions | Static Public Attributes | Protected Member Functions | Protected Attributes
SeedForPhotonConversionFromQuadruplets Class Reference

#include <SeedForPhotonConversionFromQuadruplets.h>

Public Member Functions

void bubbleReverseSortVsPhi (GlobalPoint arr[], int n, GlobalPoint vtx)
 
void bubbleSortVsPhi (GlobalPoint arr[], int n, GlobalPoint vtx)
 
double getSqrEffectiveErrorOnZ (const TransientTrackingRecHit::ConstRecHitPointer &hit, const TrackingRegion &region)
 
 SeedForPhotonConversionFromQuadruplets (const edm::ParameterSet &cfg)
 
 SeedForPhotonConversionFromQuadruplets (const std::string &propagator="PropagatorWithMaterial", double seedMomentumForBOFF=-5.0)
 
double simpleGetSlope (const TransientTrackingRecHit::ConstRecHitPointer &ohit, const TransientTrackingRecHit::ConstRecHitPointer &nohit, const TransientTrackingRecHit::ConstRecHitPointer &ihit, const TransientTrackingRecHit::ConstRecHitPointer &nihit, const TrackingRegion &region, double &cotTheta, double &z0)
 
void stupidPrint (std::string s, float *d)
 
void stupidPrint (std::string s, double *d)
 
void stupidPrint (const char *s, GlobalPoint *d)
 
void stupidPrint (const char *s, GlobalPoint *d, int n)
 
const TrajectorySeedtrajectorySeed (TrajectorySeedCollection &seedCollection, const SeedingHitSet &phits, const SeedingHitSet &mhits, const TrackingRegion &region, const edm::EventSetup &es, std::stringstream &ss, std::vector< Quad > &quadV, edm::ParameterSet &SeedComparitorPSet, edm::ParameterSet &QuadCutPSet)
 
double verySimpleFit (int size, double *ax, double *ay, double *e2y, double &p0, double &e2p0, double &p1)
 
 ~SeedForPhotonConversionFromQuadruplets ()
 

Static Public Attributes

static const int cotTheta_Max =99999
 

Protected Member Functions

const TrajectorySeedbuildSeed (TrajectorySeedCollection &seedCollection, const SeedingHitSet &hits, const FreeTrajectoryState &fts, const edm::EventSetup &es, bool apply_dzCut, const TrackingRegion &region) const
 
bool buildSeedBool (TrajectorySeedCollection &seedCollection, const SeedingHitSet &hits, const FreeTrajectoryState &fts, const edm::EventSetup &es, bool apply_dzCut, const TrackingRegion &region, double dzcut) const
 
bool checkHit (const TrajectoryStateOnSurface &, const TransientTrackingRecHit::ConstRecHitPointer &hit, const edm::EventSetup &es) const
 
double DeltaPhiManual (const math::XYZVector &v1, const math::XYZVector &v2)
 
CurvilinearTrajectoryError initialError (const GlobalVector &vertexBounds, float ptMin, float sinTheta) const
 
GlobalTrajectoryParameters initialKinematic (const SeedingHitSet &hits, const GlobalPoint &vertexPos, const edm::EventSetup &es, const float cotTheta) const
 
TransientTrackingRecHit::RecHitPointer refitHit (const TransientTrackingRecHit::ConstRecHitPointer &hit, const TrajectoryStateOnSurface &state) const
 
bool similarQuadExist (Quad &thisQuad, std::vector< Quad > &quadV)
 

Protected Attributes

double kPI_
 
PrintRecoObjects po
 
std::stringstream * pss
 
double theBOFFMomentum
 
std::string thePropagatorLabel
 

Detailed Description

Definition at line 14 of file SeedForPhotonConversionFromQuadruplets.h.

Constructor & Destructor Documentation

SeedForPhotonConversionFromQuadruplets::SeedForPhotonConversionFromQuadruplets ( const edm::ParameterSet cfg)
inline

Definition at line 18 of file SeedForPhotonConversionFromQuadruplets.h.

18  :
19  thePropagatorLabel(cfg.getParameter<std::string>("propagator")),
20  theBOFFMomentum(cfg.existsAs<double>("SeedMomentumForBOFF") ? cfg.getParameter<double>("SeedMomentumForBOFF") : 5.0)
21  {}
T getParameter(std::string const &) const
bool existsAs(std::string const &parameterName, bool trackiness=true) const
checks if a parameter exists as a given type
Definition: ParameterSet.h:184
SeedForPhotonConversionFromQuadruplets::SeedForPhotonConversionFromQuadruplets ( const std::string &  propagator = "PropagatorWithMaterial",
double  seedMomentumForBOFF = -5.0 
)
inline
SeedForPhotonConversionFromQuadruplets::~SeedForPhotonConversionFromQuadruplets ( )
inline

Definition at line 28 of file SeedForPhotonConversionFromQuadruplets.h.

28 {}

Member Function Documentation

void SeedForPhotonConversionFromQuadruplets::bubbleReverseSortVsPhi ( GlobalPoint  arr[],
int  n,
GlobalPoint  vtx 
)

Definition at line 958 of file SeedForPhotonConversionFromQuadruplets.cc.

References reco::deltaPhi(), i, j, phi, and tmp.

958  {
959  bool swapped = true;
960  int j = 0;
962  while (swapped) {
963  swapped = false;
964  j++;
965  for (int i = 0; i < n - j; i++) {
966  if ( reco::deltaPhi( (arr[i]-vtx).phi(), (arr[i + 1]-vtx).phi() ) < 0. ) {
967  tmp = arr[i];
968  arr[i] = arr[i + 1];
969  arr[i + 1] = tmp;
970  swapped = true;
971  }
972  }
973  }
974 }
int i
Definition: DBlmapReader.cc:9
int j
Definition: DBlmapReader.cc:9
double deltaPhi(double phi1, double phi2)
Definition: deltaPhi.h:12
std::vector< std::vector< double > > tmp
Definition: MVATrainer.cc:100
Definition: DDAxes.h:10
void SeedForPhotonConversionFromQuadruplets::bubbleSortVsPhi ( GlobalPoint  arr[],
int  n,
GlobalPoint  vtx 
)

Definition at line 939 of file SeedForPhotonConversionFromQuadruplets.cc.

References reco::deltaPhi(), i, j, phi, and tmp.

939  {
940  bool swapped = true;
941  int j = 0;
943  while (swapped) {
944  swapped = false;
945  j++;
946  for (int i = 0; i < n - j; i++) {
947  if ( reco::deltaPhi( (arr[i]-vtx).phi(), (arr[i + 1]-vtx).phi() ) > 0. ) {
948  tmp = arr[i];
949  arr[i] = arr[i + 1];
950  arr[i + 1] = tmp;
951  swapped = true;
952  }
953  }
954  }
955 }
int i
Definition: DBlmapReader.cc:9
int j
Definition: DBlmapReader.cc:9
double deltaPhi(double phi1, double phi2)
Definition: deltaPhi.h:12
std::vector< std::vector< double > > tmp
Definition: MVATrainer.cc:100
Definition: DDAxes.h:10
const TrajectorySeed * SeedForPhotonConversionFromQuadruplets::buildSeed ( TrajectorySeedCollection seedCollection,
const SeedingHitSet hits,
const FreeTrajectoryState fts,
const edm::EventSetup es,
bool  apply_dzCut,
const TrackingRegion region 
) const
protected

Definition at line 672 of file SeedForPhotonConversionFromQuadruplets.cc.

References alongMomentum, cond::rpcobgas::detid, TrackingRecHit::geographicalId(), edm::EventSetup::get(), TrajectoryStateOnSurface::globalMomentum(), TrackingRecHit::localPosition(), PV3DBase< T, PVType, FrameType >::perp(), trajectoryStateTransform::persistentState(), po, PrintRecoObjects::print(), Propagator::propagate(), LargeD0_PixelPairStep_cff::propagator, edm::OwnVector< T, P >::push_back(), DetId::rawId(), refitHit(), SeedingHitSet::size(), thePropagatorLabel, and patCandidatesForDimuonsSequences_cff::tracker.

Referenced by trajectorySeed().

679 {
680  // get tracker
682  es.get<TrackerDigiGeometryRecord>().get(tracker);
683 
684  // get propagator
685  edm::ESHandle<Propagator> propagatorHandle;
686  es.get<TrackingComponentsRecord>().get(thePropagatorLabel, propagatorHandle);
687  const Propagator* propagator = &(*propagatorHandle);
688 
689  // get updator
690  KFUpdator updator;
691 
692  // Now update initial state track using information from seed hits.
693 
694  TrajectoryStateOnSurface updatedState;
696 
697  const TrackingRecHit* hit = 0;
698  for ( unsigned int iHit = 0; iHit < hits.size() && iHit<2; iHit++) {
699  hit = hits[iHit]->hit();
700  TrajectoryStateOnSurface state = (iHit==0) ?
701  propagator->propagate(fts,tracker->idToDet(hit->geographicalId())->surface())
702  : propagator->propagate(updatedState, tracker->idToDet(hit->geographicalId())->surface());
703 
705 
706  TransientTrackingRecHit::RecHitPointer newtth = refitHit( tth, state);
707 
708  updatedState = updator.update(state, *newtth);
709 
710  seedHits.push_back(newtth->hit()->clone());
711 #ifdef mydebug_seed
712  uint32_t detid = hit->geographicalId().rawId();
713  (*pss) << "\n[SeedForPhotonConversionFromQuadruplets] hit " << iHit;
714  po.print(*pss, detid);
715  (*pss) << " " << detid << "\t lp " << hit->localPosition()
716  << " tth " << tth->localPosition() << " newtth " << newtth->localPosition() << " state " << state.globalMomentum().perp();
717 #endif
718  }
719 
720  PTrajectoryStateOnDet const & PTraj =
722 
723 
724  seedCollection.push_back( TrajectorySeed(PTraj,seedHits,alongMomentum));
725  return &seedCollection.back();
726 }
std::pair< ALIstring, ALIstring > pss
Definition: Fit.h:27
T perp() const
Definition: PV3DBase.h:72
PTrajectoryStateOnDet persistentState(const TrajectoryStateOnSurface &ts, unsigned int detid)
uint32_t rawId() const
get the raw id
Definition: DetId.h:43
void push_back(D *&d)
Definition: OwnVector.h:273
void print(std::stringstream &ss, const SiStripCluster &clus)
TransientTrackingRecHit::RecHitPointer refitHit(const TransientTrackingRecHit::ConstRecHitPointer &hit, const TrajectoryStateOnSurface &state) const
virtual TrajectoryStateOnSurface propagate(const FreeTrajectoryState &, const Surface &) const
Definition: Propagator.cc:12
const T & get() const
Definition: EventSetup.h:55
GlobalVector globalMomentum() const
unsigned int size() const
Definition: SeedingHitSet.h:41
DetId geographicalId() const
virtual LocalPoint localPosition() const =0
bool SeedForPhotonConversionFromQuadruplets::buildSeedBool ( TrajectorySeedCollection seedCollection,
const SeedingHitSet hits,
const FreeTrajectoryState fts,
const edm::EventSetup es,
bool  apply_dzCut,
const TrackingRegion region,
double  dzcut 
) const
protected

Implement here the dz cut:::

Definition at line 732 of file SeedForPhotonConversionFromQuadruplets.cc.

References checkHit(), gather_cfg::cout, alignCSCRings::e, TrackingRecHit::geographicalId(), edm::EventSetup::get(), TrajectoryStateOnSurface::globalMomentum(), TrajectoryStateOnSurface::globalPosition(), i, TrajectoryStateClosestToBeamLine::isValid(), TrajectoryStateOnSurface::isValid(), j, TrackingRegion::origin(), FreeTrajectoryState::position(), Propagator::propagate(), LargeD0_PixelPairStep_cff::propagator, edm::OwnVector< T, P >::push_back(), refitHit(), SeedingHitSet::size(), mathSSE::sqrt(), thePropagatorLabel, patCandidatesForDimuonsSequences_cff::tracker, TrajectoryStateClosestToBeamLine::trackStateAtPCA(), reco::BeamSpot::Unknown, findQualityFiles::v, PV3DBase< T, PVType, FrameType >::x(), PV3DBase< T, PVType, FrameType >::y(), and PV3DBase< T, PVType, FrameType >::z().

Referenced by trajectorySeed().

740 {
741  // get tracker
743  es.get<TrackerDigiGeometryRecord>().get(tracker);
744 
745  // get propagator
746  edm::ESHandle<Propagator> propagatorHandle;
747  es.get<TrackingComponentsRecord>().get(thePropagatorLabel, propagatorHandle);
748  const Propagator* propagator = &(*propagatorHandle);
749 
750  // get updator
751  KFUpdator updator;
752 
753  // Now update initial state track using information from seed hits.
754 
755  TrajectoryStateOnSurface updatedState;
757 
758  const TrackingRecHit* hit = 0;
759  for ( unsigned int iHit = 0; iHit < hits.size() && iHit<2; iHit++) {
760  hit = hits[iHit]->hit();
761  TrajectoryStateOnSurface state = (iHit==0) ?
762  propagator->propagate(fts,tracker->idToDet(hit->geographicalId())->surface())
763  : propagator->propagate(updatedState, tracker->idToDet(hit->geographicalId())->surface());
764  if (!state.isValid()) {
765  return false;}
766 
768 
769  TransientTrackingRecHit::RecHitPointer newtth = refitHit( tth, state);
770 
771 
772  if (!checkHit(state,newtth,es)){
773  return false;
774  }
775 
776  updatedState = updator.update(state, *newtth);
777  if (!updatedState.isValid()){
778  return false;
779  }
780 
781  seedHits.push_back(newtth->hit()->clone());
782  }
783 
784  if(apply_dzCut){
786 
787  double zCA;
788 
789  math::XYZVector EstMomGam(updatedState.globalMomentum().x(),updatedState.globalMomentum().y(),updatedState.globalMomentum().z());
790  math::XYZVector EstPosGam(updatedState.globalPosition().x(),updatedState.globalPosition().y(),updatedState.globalPosition().z());
791 
792  double EstMomGamLength=std::sqrt(EstMomGam.x()*EstMomGam.x()+EstMomGam.y()*EstMomGam.y()+EstMomGam.z()*EstMomGam.z());
793  math::XYZVector EstMomGamNorm(EstMomGam.x()/EstMomGamLength,EstMomGam.y()/EstMomGamLength,EstMomGam.z()/EstMomGamLength);
794 
795  //Calculate dz of point of closest approach of the two lines (WA approach) -> cut on dz
796 
797 
798  const GlobalPoint EstPosGamGlobalPoint(updatedState.globalPosition().x(),updatedState.globalPosition().y(),updatedState.globalPosition().z());
799  const GlobalVector EstMomGamGlobalVector(updatedState.globalMomentum().x(),updatedState.globalMomentum().y(),updatedState.globalMomentum().z());
800 
801 
803  es.get<IdealMagneticFieldRecord>().get(bfield);
804  const MagneticField* magField = bfield.product();
805  TrackCharge qCharge = 0;
806 
807  const GlobalTrajectoryParameters myGlobalTrajectoryParameter(EstPosGamGlobalPoint, EstMomGamGlobalVector, qCharge, magField);
808 
809  AlgebraicSymMatrix66 aCovarianceMatrix;
810 
811  for (int i =0;i<6;++i)
812  for (int j =0;j<6;++j)
813  aCovarianceMatrix(i, j) = 1e-4;
814 
815  CartesianTrajectoryError myCartesianError (aCovarianceMatrix);
816 
817  const FreeTrajectoryState stateForProjectionToBeamLine(myGlobalTrajectoryParameter,myCartesianError);
818 
819  const GlobalPoint BeamSpotGlobalPoint(0,0,0);
820 
821  const reco::BeamSpot::Point BeamSpotPoint(region.origin().x(),region.origin().y(),region.origin().z());
822 
823  TSCBLBuilderNoMaterial tscblBuilder;
824 
825  double CovMatEntry=0.;
827  for (int i=0;i<3;++i) {
828  cov(i,i) = CovMatEntry;
829  }
831 
832  reco::BeamSpot myBeamSpot(BeamSpotPoint, 0.,0.,0.,0., cov,BeamType_);
833 
834  TrajectoryStateClosestToBeamLine tscbl = tscblBuilder(stateForProjectionToBeamLine,myBeamSpot);
835  if (tscbl.isValid()==false) {
836  zCA=0;
837  }
838  else{
839  GlobalPoint v = tscbl.trackStateAtPCA().position(); // Position of closest approach to BS
840  zCA=v.z();
841  }
842 
843 /* //Calculate dz of point of closest approach of the two lines -> cut on dz
844 
845  double newX,newY,newR;
846  double Rbuff=std::sqrt(EstPosGam.x()*EstPosGam.x()+EstPosGam.y()*EstPosGam.y());
847  double deltas,s,sbuff;
848  double rMin=1e9;
849 
850  double InitXp=EstPosGam.x()+1*EstMomGamNorm.x();
851  double InitXm=EstPosGam.x()-1*EstMomGamNorm.x();
852  double InitYp=EstPosGam.y()+1*EstMomGamNorm.y();
853  double InitYm=EstPosGam.y()-1*EstMomGamNorm.y();
854 
855  if(InitXp*InitXp+InitYp*InitYp < InitXm*InitXm+InitYm*InitYm) {s=5; deltas=5;}
856  else {s=-5; deltas=-5;}
857 
858  int nTurns=0;
859  int nIterZ=0;
860  for (int i_dz=0;i_dz<1000;i_dz++){
861  newX=EstPosGam.x()+s*EstMomGamNorm.x();
862  newY=EstPosGam.y()+s*EstMomGamNorm.y();
863  newR=std::sqrt(newX*newX+newY*newY);
864  if(newR>Rbuff) {deltas=-1*deltas/10;nTurns++;}
865  else {Rbuff=newR;}
866  if(newR<rMin) {rMin=newR; sbuff=s;}
867  s=s+deltas;
868  nIterZ++;
869  if(nTurns>1) break;
870  }
871 
872  zCA=EstPosGam.z()+sbuff*EstMomGamNorm.z();
873 */
874 
875 #ifdef mydebug_knuenz
876  std::cout<< "zCA: " << zCA <<std::endl;
877 #endif
878 
879  if(std::fabs(zCA)>dzcut) return false;
880 
881 
882  }
883 
884 
885  return true;
886 }
math::Error< dimension >::type CovarianceMatrix
Definition: BeamSpot.h:31
BeamType
beam spot flags
Definition: BeamSpot.h:26
int i
Definition: DBlmapReader.cc:9
GlobalPoint const & origin() const
ROOT::Math::SMatrix< double, 6, 6, ROOT::Math::MatRepSym< double, 6 > > AlgebraicSymMatrix66
T y() const
Definition: PV3DBase.h:63
GlobalPoint globalPosition() const
math::XYZPoint Point
point in the space
Definition: BeamSpot.h:29
bool checkHit(const TrajectoryStateOnSurface &, const TransientTrackingRecHit::ConstRecHitPointer &hit, const edm::EventSetup &es) const
int TrackCharge
Definition: TrackCharge.h:4
void push_back(D *&d)
Definition: OwnVector.h:273
T sqrt(T t)
Definition: SSEVec.h:48
T z() const
Definition: PV3DBase.h:64
int j
Definition: DBlmapReader.cc:9
TransientTrackingRecHit::RecHitPointer refitHit(const TransientTrackingRecHit::ConstRecHitPointer &hit, const TrajectoryStateOnSurface &state) const
virtual TrajectoryStateOnSurface propagate(const FreeTrajectoryState &, const Surface &) const
Definition: Propagator.cc:12
GlobalPoint position() const
XYZVectorD XYZVector
spatial vector with cartesian internal representation
Definition: Vector3D.h:30
const T & get() const
Definition: EventSetup.h:55
GlobalVector globalMomentum() const
unsigned int size() const
Definition: SeedingHitSet.h:41
tuple cout
Definition: gather_cfg.py:121
DetId geographicalId() const
T x() const
Definition: PV3DBase.h:62
bool SeedForPhotonConversionFromQuadruplets::checkHit ( const TrajectoryStateOnSurface ,
const TransientTrackingRecHit::ConstRecHitPointer hit,
const edm::EventSetup es 
) const
inlineprotected

Definition at line 60 of file SeedForPhotonConversionFromQuadruplets.h.

Referenced by buildSeedBool().

63  { return true; }
double SeedForPhotonConversionFromQuadruplets::DeltaPhiManual ( const math::XYZVector v1,
const math::XYZVector v2 
)
protected

Definition at line 1069 of file SeedForPhotonConversionFromQuadruplets.cc.

References kPI_.

Referenced by trajectorySeed().

1069  {
1070 
1071  double kTWOPI = 2.*kPI_;
1072  double DeltaPhiMan=v1.Phi()-v2.Phi();
1073  while (DeltaPhiMan >= kPI_) DeltaPhiMan -= kTWOPI;
1074  while (DeltaPhiMan < -kPI_) DeltaPhiMan += kTWOPI;
1075 
1076  return DeltaPhiMan;
1077 
1078 }
double SeedForPhotonConversionFromQuadruplets::getSqrEffectiveErrorOnZ ( const TransientTrackingRecHit::ConstRecHitPointer hit,
const TrackingRegion region 
)

Definition at line 1021 of file SeedForPhotonConversionFromQuadruplets.cc.

References TrackingRegion::origin(), PV3DBase< T, PVType, FrameType >::perp(), funct::sqr(), and PV3DBase< T, PVType, FrameType >::z().

Referenced by simpleGetSlope(), and trajectorySeed().

1021  {
1022 
1023  //
1024  //Fit-wise the effective error on Z is the sum in quadrature of the error on Z
1025  //and the error on R correctly projected by using hit-vertex direction
1026 
1027  double sqrProjFactor = 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 }
T perp() const
Definition: PV3DBase.h:72
GlobalPoint const & origin() const
T z() const
Definition: PV3DBase.h:64
Square< F >::type sqr(const F &f)
Definition: Square.h:13
CurvilinearTrajectoryError SeedForPhotonConversionFromQuadruplets::initialError ( const GlobalVector vertexBounds,
float  ptMin,
float  sinTheta 
) const
protected

Definition at line 643 of file SeedForPhotonConversionFromQuadruplets.cc.

References funct::C, GlobalErrorBase< T, ErrorWeightType >::cxx(), GlobalErrorBase< T, ErrorWeightType >::czz(), max(), funct::sqr(), PV3DBase< T, PVType, FrameType >::x(), PV3DBase< T, PVType, FrameType >::y(), and PV3DBase< T, PVType, FrameType >::z().

Referenced by trajectorySeed().

647 {
648  // Set initial uncertainty on track parameters, using only P.V. constraint and no hit
649  // information.
650  GlobalError vertexErr( sqr(vertexBounds.x()), 0,
651  sqr(vertexBounds.y()), 0, 0,
652  sqr(vertexBounds.z())
653  );
654 
655 
656  AlgebraicSymMatrix55 C = ROOT::Math::SMatrixIdentity();
657 
658 // FIXME: minC00. Prevent apriori uncertainty in 1/P from being too small,
659 // to avoid instabilities.
660 // N.B. This parameter needs optimising ...
661  float sin2th = sqr(sinTheta);
662  float minC00 = 1.0;
663  C[0][0] = std::max(sin2th/sqr(ptMin), minC00);
664  float zErr = vertexErr.czz();
665  float transverseErr = vertexErr.cxx(); // assume equal cxx cyy
666  C[3][3] = transverseErr;
667  C[4][4] = zErr*sin2th + transverseErr*(1-sin2th);
668 
669  return CurvilinearTrajectoryError(C);
670 }
T y() const
Definition: PV3DBase.h:63
ROOT::Math::SMatrix< double, 5, 5, ROOT::Math::MatRepSym< double, 5 > > AlgebraicSymMatrix55
const T & max(const T &a, const T &b)
T z() const
Definition: PV3DBase.h:64
Square< F >::type sqr(const F &f)
Definition: Square.h:13
T x() const
Definition: PV3DBase.h:62
GlobalTrajectoryParameters SeedForPhotonConversionFromQuadruplets::initialKinematic ( const SeedingHitSet hits,
const GlobalPoint vertexPos,
const edm::EventSetup es,
const float  cotTheta 
) const
protected

Definition at line 584 of file SeedForPhotonConversionFromQuadruplets.cc.

References GlobalTrajectoryParameters::charge(), cotTheta_Max, cond::rpcobgas::detid, edm::EventSetup::get(), GlobalTrajectoryParameters::magneticField(), GlobalTrajectoryParameters::momentum(), PV3DBase< T, PVType, FrameType >::perp(), po, GlobalTrajectoryParameters::position(), PrintRecoObjects::print(), theBOFFMomentum, GlobalTrajectoryParameters::transverseCurvature(), Vector3DBase< T, FrameTag >::unit(), PV3DBase< T, PVType, FrameType >::x(), and PV3DBase< T, PVType, FrameType >::y().

589 {
591 
594 
596  es.get<IdealMagneticFieldRecord>().get(bfield);
597  float nomField = bfield->nominalValue();
598  bool isBOFF = (0==nomField);
599 
600 
601  FastHelix helix(tth2->globalPosition(), tth1->globalPosition(), vertexPos, nomField,&*bfield, vertexPos);
602  kine = helix.stateAtVertex();
603 
604  //force the pz/pt equal to the measured one
605  if(fabs(cotTheta)<cotTheta_Max)
606  kine = GlobalTrajectoryParameters(kine.position(),
607  GlobalVector(kine.momentum().x(),kine.momentum().y(),kine.momentum().perp()*cotTheta),
608  kine.charge(),
609  & kine.magneticField()
610  );
611  else
612  kine = GlobalTrajectoryParameters(kine.position(),
613  GlobalVector(kine.momentum().x(),kine.momentum().y(),kine.momentum().perp()*cotTheta_Max),
614  kine.charge(),
615  & kine.magneticField()
616  );
617 
618 #ifdef mydebug_seed
619  uint32_t detid;
620  (*pss) << "[SeedForPhotonConversionFromQuadruplets] initialKinematic tth1 " ;
621  detid=tth1->geographicalId().rawId();
622  po.print(*pss, detid );
623  (*pss) << " \t " << detid << " " << tth1->localPosition() << " " << tth1->globalPosition() ;
624  detid= tth2->geographicalId().rawId();
625  (*pss) << " \n\t tth2 ";
626  po.print(*pss, detid );
627  (*pss) << " \t " << detid << " " << tth2->localPosition() << " " << tth2->globalPosition()
628  << "\nhelix momentum " << kine.momentum() << " pt " << kine.momentum().perp() << " radius " << 1/kine.transverseCurvature();
629 #endif
630 
631  if (isBOFF && (theBOFFMomentum > 0)) {
632  kine = GlobalTrajectoryParameters(kine.position(),
633  kine.momentum().unit() * theBOFFMomentum,
634  kine.charge(),
635  &*bfield);
636  }
637  return kine;
638 }
std::pair< ALIstring, ALIstring > pss
Definition: Fit.h:27
T perp() const
Definition: PV3DBase.h:72
T y() const
Definition: PV3DBase.h:63
void print(std::stringstream &ss, const SiStripCluster &clus)
Vector3DBase unit() const
Definition: Vector3DBase.h:57
const T & get() const
Definition: EventSetup.h:55
const MagneticField & magneticField() const
T x() const
Definition: PV3DBase.h:62
Global3DVector GlobalVector
Definition: GlobalVector.h:10
TransientTrackingRecHit::RecHitPointer SeedForPhotonConversionFromQuadruplets::refitHit ( const TransientTrackingRecHit::ConstRecHitPointer hit,
const TrajectoryStateOnSurface state 
) const
protected

Definition at line 890 of file SeedForPhotonConversionFromQuadruplets.cc.

Referenced by buildSeed(), and buildSeedBool().

893 {
894  //const TransientTrackingRecHit* a= hit.get();
895  //return const_cast<TransientTrackingRecHit*> (a);
896  //This was modified otherwise the rechit will have just the local x component and local y=0
897  // To understand how to modify for pixels
898 
899  //const TSiStripRecHit2DLocalPos* b = dynamic_cast<const TSiStripRecHit2DLocalPos*>(a);
900  //return const_cast<TSiStripRecHit2DLocalPos*>(b);
901  return hit->clone(state);
902 }
bool SeedForPhotonConversionFromQuadruplets::similarQuadExist ( Quad thisQuad,
std::vector< Quad > &  quadV 
)
protected

Definition at line 1032 of file SeedForPhotonConversionFromQuadruplets.cc.

References Quad::ptMinus, Quad::ptPlus, mathSSE::sqrt(), Quad::x, and Quad::y.

Referenced by trajectorySeed().

1032  {
1033 
1034  BOOST_FOREACH( Quad quad, quadV )
1035  {
1036  double dx = thisQuad.x-quad.x;
1037  double dy = thisQuad.y-quad.y;
1038  //double dz = abs(thisQuad.z-quad.z);
1039  //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
1040  //std::cout<<"x1-x2="<<dx<<"y1-y2="<<dy<<"dx*dx+dy*dy="<<dx*dx+dy*dy<<std::endl; //ValDebug
1041  //std::cout<<"thisQuad.ptPlus-quad.ptPlus="<<thisQuad.ptPlus-quad.ptPlus<<"abs(thisQuad.ptPlus-quad.ptPlus)="<<abs(thisQuad.ptPlus-quad.ptPlus)<<std::endl; //ValDebug
1042  //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
1043  //if ( sqrt(dx*dx+dy*dy)<1.)
1044  // 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
1045  // ( sqrt(dx*dx+dy*dy)<1. &&
1046  //z<3. &&
1047  //bs(thisQuad.ptPlus-quad.ptPlus)<0.5*quad.ptPlus &&
1048  //bs(thisQuad.ptMinus-quad.ptMinus)<0.5*quad.ptMinus &&
1049  //bs(thisQuad.cot-quad.cot)<0.3*quad.cot
1050  //
1051  // {
1052  // //std::cout<<"Seed rejected due to arbitration"<<std::endl;
1053  // return true;
1054  // }
1055  if ( sqrt(dx*dx+dy*dy)<1. &&
1056  fabs(thisQuad.ptPlus-quad.ptPlus)<0.5*quad.ptPlus &&
1057  fabs(thisQuad.ptMinus-quad.ptMinus)<0.5*quad.ptMinus
1058  )
1059  {
1060  //std::cout<<"Seed rejected due to arbitration"<<std::endl;
1061  return true;
1062  }
1063  }
1064  quadV.push_back(thisQuad);
1065  return false;
1066 
1067 }
double ptMinus
Definition: Quad.h:9
double ptPlus
Definition: Quad.h:8
double x
Definition: Quad.h:5
T sqrt(T t)
Definition: SSEVec.h:48
double y
Definition: Quad.h:6
Definition: Quad.h:4
double SeedForPhotonConversionFromQuadruplets::simpleGetSlope ( const TransientTrackingRecHit::ConstRecHitPointer ohit,
const TransientTrackingRecHit::ConstRecHitPointer nohit,
const TransientTrackingRecHit::ConstRecHitPointer ihit,
const TransientTrackingRecHit::ConstRecHitPointer nihit,
const TrackingRegion region,
double &  cotTheta,
double &  z0 
)

Definition at line 978 of file SeedForPhotonConversionFromQuadruplets.cc.

References getSqrEffectiveErrorOnZ(), TrackingRegion::origin(), TrackingRegion::originZBound(), PV3DBase< T, PVType, FrameType >::perp(), funct::sqr(), verySimpleFit(), x, detailsBasic3DVector::y, and PV3DBase< T, PVType, FrameType >::z().

Referenced by trajectorySeed().

978  {
979 
980  double x[5], y[5], e2y[5];
981 
982  //The fit is done filling x with r values, y with z values of the four hits and the vertex
983  //
984  //Hits
985  x[0] = ohit->globalPosition().perp();
986  y[0] = ohit->globalPosition().z();
987  e2y[0] = getSqrEffectiveErrorOnZ(ohit, region);
988  //
989  x[1] = nohit->globalPosition().perp();
990  y[1] = nohit->globalPosition().z();
991  e2y[1] = getSqrEffectiveErrorOnZ(nohit, region);
992  //
993  x[2] = nihit->globalPosition().perp();
994  y[2] = nihit->globalPosition().z();
995  e2y[2] = getSqrEffectiveErrorOnZ(nihit, region);
996  //
997  x[3] = ihit->globalPosition().perp();
998  y[3] = ihit->globalPosition().z();
999  e2y[3] = getSqrEffectiveErrorOnZ(ihit, region);
1000  //
1001  //Vertex
1002  x[4] = region.origin().perp();
1003  y[4] = region.origin().z();
1004  double vError = region.originZBound();
1005  if ( vError > 15. ) 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 }
T perp() const
Definition: PV3DBase.h:72
GlobalPoint const & origin() const
double getSqrEffectiveErrorOnZ(const TransientTrackingRecHit::ConstRecHitPointer &hit, const TrackingRegion &region)
double verySimpleFit(int size, double *ax, double *ay, double *e2y, double &p0, double &e2p0, double &p1)
T z() const
Definition: PV3DBase.h:64
float originZBound() const
bounds the particle vertex in the longitudinal plane
Square< F >::type sqr(const F &f)
Definition: Square.h:13
Definition: DDAxes.h:10
void SeedForPhotonConversionFromQuadruplets::stupidPrint ( std::string  s,
float *  d 
)

Definition at line 909 of file SeedForPhotonConversionFromQuadruplets.cc.

References i.

909  {
910  (*pss) << "\n" << s << "\t";
911  for(size_t i=0;i<2;++i)
912  (*pss) << std::setw (60) << d[i] << std::setw(1) << " | ";
913 }
int i
Definition: DBlmapReader.cc:9
std::pair< ALIstring, ALIstring > pss
Definition: Fit.h:27
void SeedForPhotonConversionFromQuadruplets::stupidPrint ( std::string  s,
double *  d 
)

Definition at line 916 of file SeedForPhotonConversionFromQuadruplets.cc.

References i.

916  {
917  (*pss) << "\n" << s << "\t";
918  for(size_t i=0;i<2;++i)
919  (*pss) << std::setw (60) << d[i] << std::setw(1) << " | ";
920 }
int i
Definition: DBlmapReader.cc:9
std::pair< ALIstring, ALIstring > pss
Definition: Fit.h:27
void SeedForPhotonConversionFromQuadruplets::stupidPrint ( const char *  s,
GlobalPoint d 
)

Definition at line 923 of file SeedForPhotonConversionFromQuadruplets.cc.

References i, PV3DBase< T, PVType, FrameType >::perp(), and PV3DBase< T, PVType, FrameType >::phi().

923  {
924  (*pss) << "\n" << s << "\t";
925  for(size_t i=0;i<2;++i)
926  (*pss) << std::setw(20) << d[i] << " r " << d[i].perp() << " phi " << d[i].phi() << " | ";
927 }
int i
Definition: DBlmapReader.cc:9
std::pair< ALIstring, ALIstring > pss
Definition: Fit.h:27
T perp() const
Definition: PV3DBase.h:72
Geom::Phi< T > phi() const
Definition: PV3DBase.h:69
void SeedForPhotonConversionFromQuadruplets::stupidPrint ( const char *  s,
GlobalPoint d,
int  n 
)

Definition at line 930 of file SeedForPhotonConversionFromQuadruplets.cc.

References i, n, PV3DBase< T, PVType, FrameType >::perp(), and PV3DBase< T, PVType, FrameType >::phi().

930  {
931  (*pss) << "\n" << s << "\n";
932  for(int i=0;i<n;++i)
933  (*pss) << std::setw(20) << d[i] << " r " << d[i].perp() << " phi " << d[i].phi() << "\n";
934 }
int i
Definition: DBlmapReader.cc:9
std::pair< ALIstring, ALIstring > pss
Definition: Fit.h:27
T perp() const
Definition: PV3DBase.h:72
Geom::Phi< T > phi() const
Definition: PV3DBase.h:69
const TrajectorySeed * SeedForPhotonConversionFromQuadruplets::trajectorySeed ( TrajectorySeedCollection seedCollection,
const SeedingHitSet phits,
const SeedingHitSet mhits,
const TrackingRegion region,
const edm::EventSetup es,
std::stringstream &  ss,
std::vector< Quad > &  quadV,
edm::ParameterSet SeedComparitorPSet,
edm::ParameterSet QuadCutPSet 
)

Definition at line 37 of file SeedForPhotonConversionFromQuadruplets.cc.

References funct::abs(), buildSeed(), buildSeedBool(), funct::C, Conv4HitsReco2::ConversionCandidate(), funct::cos(), gather_cfg::cout, DeltaPhiManual(), cond::rpcobgas::detid, Conv4HitsReco2::Dump(), reco::get(), edm::EventSetup::get(), edm::ParameterSet::getParameter(), getSqrEffectiveErrorOnZ(), initialError(), listHistos::IP, gen::k, kPI_, TrackingRegion::origin(), TrackingRegion::originRBound(), TrackingRegion::originZBound(), colinearityKinematic::Phi, TrackingRegion::ptMin(), ptmin, Conv4HitsReco2::SetMaxNumberOfIterations(), similarQuadExist(), simpleGetSlope(), funct::sin(), SeedingHitSet::size(), mathSSE::sqrt(), AlCaHLTBitMon_QueryRunRegistry::string, Quad::x, x, PV3DBase< T, PVType, FrameType >::x(), detailsBasic3DVector::y, PV3DBase< T, PVType, FrameType >::y(), detailsBasic3DVector::z, and PV3DBase< T, PVType, FrameType >::z().

Referenced by PhotonConversionTrajectorySeedProducerFromQuadrupletsAlgo::inspect().

47 {
48 
49 
51 
52 
53  bool rejectAllQuads=QuadCutPSet.getParameter<bool>("rejectAllQuads");
54  if(rejectAllQuads) return 0;
55 
56  bool applyDeltaPhiCuts=QuadCutPSet.getParameter<bool>("apply_DeltaPhiCuts");
57  bool ClusterShapeFiltering=QuadCutPSet.getParameter<bool>("apply_ClusterShapeFilter");
58  bool applyArbitration=QuadCutPSet.getParameter<bool>("apply_Arbitration");
59  bool applydzCAcut=QuadCutPSet.getParameter<bool>("apply_zCACut");
60  double CleaningmaxRadialDistance=QuadCutPSet.getParameter<double>("Cut_DeltaRho");//cm
61  double BeamPipeRadiusCut=QuadCutPSet.getParameter<double>("Cut_BeamPipeRadius");//cm
62  double CleaningMinLegPt = QuadCutPSet.getParameter<double>("Cut_minLegPt"); //GeV
63  double maxLegPt = QuadCutPSet.getParameter<double>("Cut_maxLegPt"); //GeV
64  double dzcut=QuadCutPSet.getParameter<double>("Cut_zCA");//cm
65 
66  double toleranceFactorOnDeltaPhiCuts=0.1;
67 
68  // return 0; //FIXME, remove this line to make the code working.
69 
70  pss = &ss;
71 
72  if ( phits.size() < 2) return 0;
73  if ( mhits.size() < 2) return 0;
74 
75  //PUT HERE THE QUADRUPLET ALGORITHM, AND IN CASE USE THE METHODS ALREADY DEVELOPED, ADAPTING THEM
76 
77  //
78  // Rozzo ma efficace (per ora)
79  //
80 
81 #ifdef mydebug_sguazz
82  std::cout << " --------------------------------------------------------------------------" << "\n";
83  std::cout << " Starting a hit quad fast reco " << "\n";
84  std::cout << " --------------------------------------------------------------------------" << "\n";
85 #endif
86 
87  //
88  // Let's build the 4 hits
93 
94  GlobalPoint vHit[4];
95  vHit[0]=ptth2->globalPosition();
96  vHit[1]=ptth1->globalPosition();
97  vHit[2]=mtth1->globalPosition();
98  vHit[3]=mtth2->globalPosition();
99 
100  //Photon source vertex primary vertex
101  GlobalPoint vgPhotVertex=region.origin();
102  math::XYZVector vPhotVertex(vgPhotVertex.x(), vgPhotVertex.y(), vgPhotVertex.z());
103 
104  math::XYZVector h1(vHit[0].x(),vHit[0].y(),vHit[0].z());
105  math::XYZVector h2(vHit[1].x(),vHit[1].y(),vHit[1].z());
106  math::XYZVector h3(vHit[2].x(),vHit[2].y(),vHit[2].z());
107  math::XYZVector h4(vHit[3].x(),vHit[3].y(),vHit[3].z());
108 
109 // At this point implement cleaning cuts before building the seed:::
110 
111 /*
112  Notes:
113 
114 h1, h2: positron
115 h3, h4: electron
116 
117 P1, P2: positron, ordered with radius
118 M1, M2: electron, ordered with radius
119 
120 Evan's notation:
121 V1=P1
122 V2=M1
123 V3=P2
124 V4=M2
125 
126 */
127 
128  math::XYZVector P1;
129  math::XYZVector P2;
130  math::XYZVector M1;
131  math::XYZVector M2;
132 
133  if(h1.x()*h1.x()+h1.y()*h1.y() < h2.x()*h2.x()+h2.y()*h2.y()){
134  P1=h1;
135  P2=h2;
136  }
137  else{
138  P1=h2;
139  P2=h1;
140  }
141 
142  if(h3.x()*h3.x()+h3.y()*h3.y() < h4.x()*h4.x()+h4.y()*h4.y()){
143  M1=h3;
144  M2=h4;
145  }
146  else{
147  M1=h4;
148  M2=h3;
149  }
150 
152 // Intersection-point:::
154 /*
155 Calculate the intersection point of the lines P2-P1 and M2-M1.
156 If this point is in the beam pipe, or if the distance of this
157 point to the layer of the most inner hit of the seed is less
158 than CleaningmaxRadialDistance cm, the combination is rejected.
159 */
160 
161 
162  math::XYZVector IP(0,0,0);
163 
164  //Line1:
165  double kP=(P1.y()-P2.y())/(P1.x()-P2.x());
166  double dP=P1.y()-kP*P1.x();
167  //Line2:
168  double kM=(M1.y()-M2.y())/(M1.x()-M2.x());
169  double dM=M1.y()-kM*M1.x();
170  //Intersection:
171  double IPx=(dM-dP)/(kP-kM);
172  double IPy=kP*IPx+dP;
173 
174  IP.SetXYZ(IPx,IPy,0);
175 
176  double IPrho=std::sqrt(IP.x()*IP.x()+IP.y()*IP.y());
177  double P1rho2=P1.x()*P1.x()+P1.y()*P1.y();
178  double M1rho2=M1.x()*M1.x()+M1.y()*M1.y();
179  double maxIPrho2=IPrho+CleaningmaxRadialDistance; maxIPrho2*=maxIPrho2;
180 
181  if( IPrho<BeamPipeRadiusCut || P1rho2>maxIPrho2 || M1rho2>maxIPrho2){
182  return 0;
183  }
184 
185  if(applyDeltaPhiCuts) {
186 
187  kPI_ = std::atan(1.0)*4;
188 
190  es.get<IdealMagneticFieldRecord>().get(bfield);
191  math::XYZVector QuadMean(0,0,0);
192  QuadMean.SetXYZ((M1.x()+M2.x()+P1.x()+P2.x())/4.,(M1.y()+M2.y()+P1.y()+P2.y())/4.,(M1.z()+M2.z()+P1.z()+P2.z())/4.);
193 
194  double fBField = bfield->inTesla(GlobalPoint(QuadMean.x(),QuadMean.y(),QuadMean.z())).z();
195 
196  double rMax=CleaningMinLegPt/(0.01*0.3*fBField);
197  double rMax_squared=rMax*rMax;
198  double Mx=M1.x();
199  double My=M1.y();
200 
201  math::XYZVector B(0,0,0);
202  math::XYZVector C(0,0,0);
203 
204  if(rMax_squared*4. > Mx*Mx+My*My){
205 
207 // Cleaning P1 points:::
209 
210  //Cx, Cy = Coordinates of circle center
211  //C_=line that contains the circle center
212 
213  //C_=k*x+d
214  double k=-Mx/My;
215  double d=My/2.-k*Mx/2.;
216 
217 #ifdef mydebug_knuenz
218 std::cout << "k" << k << std::endl;
219 std::cout << "d" << d << std::endl;
220 #endif
221 
222  //Cx1,2 and Cy1,2 are the two points that have a distance of rMax to 0,0
223  double CsolutionPart1=-2*k*d;
224  double CsolutionPart2=std::sqrt(4*k*k*d*d-4*(1+k*k)*(d*d-rMax_squared));
225  double CsolutionPart3=2*(1+k*k);
226  double Cx1=(CsolutionPart1+CsolutionPart2)/CsolutionPart3;
227  double Cx2=(CsolutionPart1-CsolutionPart2)/CsolutionPart3;
228  double Cy1=k*Cx1+d;
229  double Cy2=k*Cx2+d;
230 
231 
232  // Decide between solutions: phi(C) > phi(P)
233  double Cx,Cy;
234  math::XYZVector C1(Cx1,Cy1,0);
235  if(C1.x()*M1.y()-C1.y()*M1.x()<0){
236  Cx=Cx1;
237  Cy=Cy1;
238  }
239  else{
240  Cx=Cx2;
241  Cy=Cy2;
242  }
243  C.SetXYZ(Cx,Cy,0);
244 
245 #ifdef mydebug_knuenz
246  std::cout << "Cx1" << Cx1 << std::endl;
247  std::cout << "Cx2" << Cx2 << std::endl;
248  std::cout << "Cy1" << Cy1 << std::endl;
249  std::cout << "Cy2" << Cy2 << std::endl;
250  std::cout << "Cx" << Cx << std::endl;
251  std::cout << "Cy" << Cy << std::endl;
252 #endif
253 
254 // Find Tangent at 0,0 to minPtCircle and point (Bx,By) on the first layer which bisects the allowed angle
255  k=-Cx/Cy;
256  d=0;
257  double Bx1=std::sqrt(Mx*Mx+My*My/(1+k*k));
258  double Bx2=-Bx1;
259  double By1=k*Bx1+d;
260  double By2=k*Bx2+d;
261 
262 #ifdef mydebug_knuenz
263  std::cout << "k" << k << std::endl;
264  std::cout << "d" << d << std::endl;
265 #endif
266 
267 // Decide between solutions: phi(B) < phi(P)
268  double Bx,By;
269  math::XYZVector B1(Bx1,By1,0);
270  if(M1.x()*B1.y()-M1.y()*B1.x()<0){
271  Bx=Bx1;
272  By=By1;
273  }
274  else{
275  Bx=Bx2;
276  By=By2;
277  }
278  B.SetXYZ(Bx,By,0);
279 
280 #ifdef mydebug_knuenz
281  std::cout << "Bx1" << Bx1 << std::endl;
282  std::cout << "Bx2" << Bx2 << std::endl;
283  std::cout << "By1" << By1 << std::endl;
284  std::cout << "By2" << By2 << std::endl;
285  std::cout << "Bx" << Bx << std::endl;
286  std::cout << "By" << By << std::endl;
287 #endif
288 
289  double DeltaPhiMaxM1P1=DeltaPhiManual(M1,B)*2;
290 
291 #ifdef mydebug_knuenz
292  std::cout << "DeltaPhiMaxM1P1 " << DeltaPhiMaxM1P1 << std::endl;
293  std::cout << "M1.DeltaPhi(P1) " << DeltaPhiManual(M1,P1) << std::endl;
294  std::cout << "rho P1: " << std::sqrt(P1.x()*P1.x()+P1.y()*P1.y()) << "phi P1: " << P1.Phi() << std::endl;
295  std::cout << "rho M1: " << std::sqrt(M1.x()*M1.x()+M1.y()*M1.y()) << "phi M1: " << M1.Phi() << std::endl;
296 #endif
297 
298 //Finally Cut on DeltaPhi of P1 and M1
299 
300  double tol_DeltaPhiMaxM1P1=DeltaPhiMaxM1P1*toleranceFactorOnDeltaPhiCuts;
301  double DeltaPhiManualM1P1=DeltaPhiManual(M1,P1);
302 
303 if(DeltaPhiManualM1P1>DeltaPhiMaxM1P1+tol_DeltaPhiMaxM1P1 || DeltaPhiManualM1P1<0-tol_DeltaPhiMaxM1P1){
304  return 0;
305 }
306 
307  }//if rMax > rLayerM1
308 
309 
311 // Cleaning M2 points:::
313 
314 // if(B.DeltaPhi(P1)>0){//normal algo (with minPt circle)
315 
316  double rM2_squared=M2.x()*M2.x()+M2.y()*M2.y();
317  if(rMax_squared*4. > rM2_squared){//if minPt circle is smaller than 2*M2-layer radius, algo makes no sense
318 
319  //Chordales equation (line containing the two intersection points of the two circles)
320  double k=-C.x()/C.y();
321  double d=(rM2_squared-rMax_squared+C.x()*C.x()+C.y()*C.y())/(2*C.y());
322 
323  double M2solutionPart1=-2*k*d;
324  double M2solutionPart2=std::sqrt(4*k*k*d*d-4*(1+k*k)*(d*d-rM2_squared));
325  double M2solutionPart3=2+2*k*k;
326  double M2xMax1=(M2solutionPart1+M2solutionPart2)/M2solutionPart3;
327  double M2xMax2=(M2solutionPart1-M2solutionPart2)/M2solutionPart3;
328  double M2yMax1=k*M2xMax1+d;
329  double M2yMax2=k*M2xMax2+d;
330 
331  //double M2xMax,M2yMax;
332  math::XYZVector M2MaxVec1(M2xMax1,M2yMax1,0);
333  math::XYZVector M2MaxVec2(M2xMax2,M2yMax2,0);
334  math::XYZVector M2MaxVec(0,0,0);
335  if(M2MaxVec1.x()*M2MaxVec2.y()-M2MaxVec1.y()*M2MaxVec2.x()<0){
336  M2MaxVec.SetXYZ(M2xMax2,M2yMax2,0);
337  }
338  else{
339  M2MaxVec.SetXYZ(M2xMax1,M2yMax1,0);
340  }
341 
342  double DeltaPhiMaxM2=DeltaPhiManual(M2MaxVec,M1);
343 
344 #ifdef mydebug_knuenz
345  std::cout << "C.x() " << C.x() << std::endl;
346  std::cout << "C.y() " << C.y() << std::endl;
347  std::cout << "M1.x() " << M1.x() << std::endl;
348  std::cout << "M1.y() " << M1.y() << std::endl;
349  std::cout << "M2.x() " << M2.x() << std::endl;
350  std::cout << "M2.y() " << M2.y() << std::endl;
351  std::cout << "k " << k << std::endl;
352  std::cout << "d " << d << std::endl;
353  std::cout << "M2xMax1 " << M2xMax1 << std::endl;
354  std::cout << "M2xMax2 " << M2xMax2 << std::endl;
355  std::cout << "M2yMax1 " << M2yMax1 << std::endl;
356  std::cout << "M2yMax2 " << M2yMax2 << std::endl;
357  std::cout << "M2xMax " << M2MaxVec.x() << std::endl;
358  std::cout << "M2yMax " << M2MaxVec.y() << std::endl;
359  std::cout << "rM2_squared " << rM2_squared << std::endl;
360  std::cout << "rMax " << rMax << std::endl;
361  std::cout << "DeltaPhiMaxM2 " << DeltaPhiMaxM2 << std::endl;
362 #endif
363 
364 
365  double tol_DeltaPhiMaxM2=DeltaPhiMaxM2*toleranceFactorOnDeltaPhiCuts;
366  double DeltaPhiManualM2M1=DeltaPhiManual(M2,M1);
367 
368  if(DeltaPhiManualM2M1>DeltaPhiMaxM2+tol_DeltaPhiMaxM2 || DeltaPhiManualM2M1<0-tol_DeltaPhiMaxM2){
369  return 0;
370  }
371 
372  //Using the lazy solution for P2: DeltaPhiMaxP2=DeltaPhiMaxM2
373  double DeltaPhiManualP1P2=DeltaPhiManual(P1,P2);
374  if(DeltaPhiManualP1P2>DeltaPhiMaxM2+tol_DeltaPhiMaxM2 || DeltaPhiManualP1P2<0-tol_DeltaPhiMaxM2){
375  return 0;
376  }
377 
378  }
379 
380 }
381 
382 // }
383 
384 // if(B.DeltaPhi(P1)<0){//different algo (without minPt circle)
385 
386 // }
387 
388 // End of pre-seed cleaning
389 
390 
391 
392  Conv4HitsReco2 quad(vPhotVertex, h1, h2, h3, h4);
393  quad.SetMaxNumberOfIterations(100);
394 #ifdef mydebug_sguazz
395  quad.Dump();
396 #endif
397  math::XYZVector candVtx;
398  double candPtPlus, candPtMinus;
399  //double rPlus, rMinus;
400  int nite = quad.ConversionCandidate(candVtx, candPtPlus, candPtMinus);
401 
402  if ( ! (nite && abs(nite) < 25 && nite != -1000 && nite != -2000) ) return 0;
403 
404 // math::XYZVector plusCenter = quad.GetPlusCenter(rPlus);
405 // math::XYZVector minusCenter = quad.GetMinusCenter(rMinus);
406 
407 #ifdef mydebug_sguazz
408  std::cout << " >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>" << "\n";
409  std::cout << " >>>>>>>>>>> Conv Cand: " << " Vertex X: " << candVtx.X() << " [cm] Y: " << candVtx.Y() << " [cm] pt+: " << candPtPlus<< " [GeV] pt-: " << candPtMinus << " [GeV]; #its: " << nite << "\n";
410 #endif
411 
412  //Add here cuts
413  double minLegPt = CleaningMinLegPt;
414  double maxRadialDistance = CleaningmaxRadialDistance;
415 
416  //
417  // Cut on leg's transverse momenta
418  if ( candPtPlus < minLegPt ) return 0;
419  if ( candPtMinus < minLegPt ) return 0;
420  //
421  if ( candPtPlus > maxLegPt ) return 0;
422  if ( candPtMinus > maxLegPt ) return 0;
423  //
424  // Cut on radial distance between estimated conversion vertex and inner hits
425  double cr = std::sqrt(candVtx.Perp2());
426  double maxr2 = (maxRadialDistance + cr); maxr2*=maxr2;
427  if (h2.Perp2() > maxr2) return 0;
428  if (h3.Perp2() > maxr2) return 0;
429 
430 
431 // At this point implement cleaning cuts after building the seed
432 
433  //ClusterShapeFilter_knuenz:::
434  std::string comparitorName = SeedComparitorPSet.getParameter<std::string>("ComponentName");
435  SeedComparitor * theComparitor = (comparitorName == "none") ? 0 : SeedComparitorFactory::get()->create( comparitorName, SeedComparitorPSet);
437  es.get<IdealMagneticFieldRecord>().get(bfield);
438  float nomField = bfield->nominalValue();
439 
440 
441  if(ClusterShapeFiltering){
442  if (theComparitor) theComparitor->init(es);
443 
446 
451 
452  GlobalPoint vertexPos(candVtx.x(),candVtx.y(),candVtx.z());
453 
454  float ptMinReg=0.1;
455  GlobalTrackingRegion region(ptMinReg,vertexPos,0,0,true);
456 
457  FastHelix phelix(ptth2->globalPosition(), mtth1->globalPosition(), vertexPos, nomField,&*bfield, vertexPos);
458  pkine = phelix.stateAtVertex();
459  FastHelix mhelix(mtth2->globalPosition(), mtth1->globalPosition(), vertexPos, nomField,&*bfield, vertexPos);
460  mkine = mhelix.stateAtVertex();
461 
462  if(theComparitor&&!theComparitor->compatible(phits, pkine, phelix, region)) { return 0; }
463  if(theComparitor&&!theComparitor->compatible(mhits, mkine, mhelix, region)) { return 0; }
464  }
465 
466 
467  //Do a very simple fit to estimate the slope
468  double quadPhotCotTheta = 0.;
469  double quadZ0 = 0.;
470  simpleGetSlope(ptth2, ptth1, mtth1, mtth2, region, quadPhotCotTheta, quadZ0);
471 
472  double quadPhotPhi = (candVtx-vPhotVertex).Phi();
473 
474  math::XYZVector fittedPrimaryVertex(vgPhotVertex.x(), vgPhotVertex.y(),quadZ0);
475 
476  candVtx.SetZ(std::sqrt(candVtx.Perp2())*quadPhotCotTheta+quadZ0);
477  GlobalPoint convVtxGlobalPoint(candVtx.X(),candVtx.Y(),candVtx.Z());
478 
479  //
480  // Comparing new quad with old quad
481  //
482  //Arbitration
483 
484  Quad thisQuad;
485  thisQuad.x = candVtx.X();
486  thisQuad.y = candVtx.Y();
487  thisQuad.z = candVtx.Z();
488  thisQuad.ptPlus = candPtPlus;
489  thisQuad.ptMinus = candPtMinus;
490  thisQuad.cot = quadPhotCotTheta;
491  if ( similarQuadExist(thisQuad, quadV) && applyArbitration ) return 0;
492 
493  // not able to get the mag field... doing the dirty way
494  //
495  // Plus
496  FastHelix helixPlus(ptth2->globalPosition(), ptth1->globalPosition(), convVtxGlobalPoint, nomField,&*bfield, convVtxGlobalPoint);
497  GlobalTrajectoryParameters kinePlus = helixPlus.stateAtVertex();
498  kinePlus = GlobalTrajectoryParameters(convVtxGlobalPoint,
499  GlobalVector(candPtPlus*cos(quadPhotPhi),candPtPlus*sin(quadPhotPhi),candPtPlus*quadPhotCotTheta),
500  1,//1
501  & kinePlus.magneticField()
502  );
503 
504  //
505  // Minus
506  FastHelix helixMinus(mtth2->globalPosition(), mtth1->globalPosition(), convVtxGlobalPoint, nomField,&*bfield, convVtxGlobalPoint);
507  GlobalTrajectoryParameters kineMinus = helixMinus.stateAtVertex();
508  kineMinus = GlobalTrajectoryParameters(convVtxGlobalPoint,
509  GlobalVector(candPtMinus*cos(quadPhotPhi),candPtMinus*sin(quadPhotPhi),candPtMinus*quadPhotCotTheta),
510  -1,//-1
511  & kineMinus.magneticField()
512  );
513 
514  float sinThetaPlus = sin(kinePlus.momentum().theta());
515  float sinThetaMinus = sin(kineMinus.momentum().theta());
516  float ptmin = region.ptMin();
517  //vertexBounds da region
518  GlobalVector vertexBounds(region.originRBound(),region.originRBound(),region.originZBound());
519 
520  CurvilinearTrajectoryError errorPlus = initialError(vertexBounds, ptmin, sinThetaPlus);
521  CurvilinearTrajectoryError errorMinus = initialError(vertexBounds, ptmin, sinThetaMinus);
522  FreeTrajectoryState ftsPlus(kinePlus, errorPlus);
523  FreeTrajectoryState ftsMinus(kineMinus, errorMinus);
524 
525  //FIXME: here probably you want to go in parallel with phits and mhits
526  //NB: the seedCollection is filled (the important thing) the return of the last TrajectorySeed is not used, but is needed
527  //to maintain the inheritance
528 
529 #ifdef quadDispLine
530  double vError = region.originZBound();
531  if ( vError > 15. ) vError = 1.;
532  std::cout << "QuadDispLine "
533  << vgPhotVertex.x() << " " << vgPhotVertex.y() << " " << vgPhotVertex.z() << " " << vError << " "
534  << vHit[0].x() << " " << vHit[0].y() << " " << vHit[0].z() << " " << sqrt(getSqrEffectiveErrorOnZ(ptth2, region)) << " "
535  << vHit[1].x() << " " << vHit[1].y() << " " << vHit[1].z() << " " << sqrt(getSqrEffectiveErrorOnZ(ptth1, region)) << " "
536  << vHit[2].x() << " " << vHit[2].y() << " " << vHit[2].z() << " " << sqrt(getSqrEffectiveErrorOnZ(mtth1, region)) << " "
537  << vHit[3].x() << " " << vHit[3].y() << " " << vHit[3].z() << " " << sqrt(getSqrEffectiveErrorOnZ(mtth2, region)) << " "
538  << candVtx.X() << " " << candVtx.Y() << " " << candVtx.Z() << " "
539  << fittedPrimaryVertex.X() << " " << fittedPrimaryVertex.Y() << " " << fittedPrimaryVertex.Z() << " "
540 // << candPtPlus << " " << rPlus << " " << plusCenter.X() << " " << plusCenter.Y() << " "
541 // << candPtMinus << " " << rMinus << " " << minusCenter.X() << " " << minusCenter.Y() << " "
542  << nite << " " << chi2 << "\n";
543 #endif
544 #ifdef mydebug_sguazz
545  std::cout << " >>>>> Hit quad fast reco done >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>" << "\n";
546  uint32_t detid;
547  std::cout << "[SeedForPhotonConversionFromQuadruplets]\n ptth1 " ;
548  detid=ptth1->geographicalId().rawId();
549  // po.print(std::cout , detid );
550  std::cout << " \t " << detid << " " << ptth1->localPosition() << " " << ptth1->globalPosition() ;
551  detid=ptth2->geographicalId().rawId();
552  std::cout << " \n\t ptth2 ";
553  // po.print(std::cout , detid );
554  std::cout << " \t " << detid << " " << ptth2->localPosition() << " " << ptth2->globalPosition()
555  << "\nhelix momentum " << kinePlus.momentum() << " pt " << kinePlus.momentum().perp() << " radius " << 1/kinePlus.transverseCurvature() << " q " << kinePlus.charge();
556  std::cout << " \n\t mtth1 ";
557  detid=mtth1->geographicalId().rawId();
558  std::cout << " \t " << detid << " " << mtth1->localPosition() << " " << mtth1->globalPosition() ;
559  std::cout << " \n\t mtth2 ";
560  detid=mtth2->geographicalId().rawId();
561  // po.print(std::cout , detid );
562  std::cout << " \t " << detid << " " << mtth2->localPosition() << " " << mtth2->globalPosition()
563  << "\nhelix momentum " << kineMinus.momentum() << " pt " << kineMinus.momentum().perp() << " radius " << 1/kineMinus.transverseCurvature() << " q " << kineMinus.charge();
564  std::cout << "\n <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<" << "\n";
565 #endif
566 
567 
568 
569  bool buildSeedBoolPos = buildSeedBool(seedCollection,phits,ftsPlus,es,applydzCAcut,region, dzcut);
570  bool buildSeedBoolNeg = buildSeedBool(seedCollection,mhits,ftsMinus,es,applydzCAcut,region, dzcut);
571 
572 
573 
574  if( buildSeedBoolPos && buildSeedBoolNeg ){
575  buildSeed(seedCollection,phits,ftsPlus,es,false,region);
576  buildSeed(seedCollection,mhits,ftsMinus,es,false,region);
577  }
578 
579  return 0;
580 
581 }
float originRBound() const
bounds the particle vertex in the transverse plane
T getParameter(std::string const &) const
CurvilinearTrajectoryError initialError(const GlobalVector &vertexBounds, float ptMin, float sinTheta) const
bool buildSeedBool(TrajectorySeedCollection &seedCollection, const SeedingHitSet &hits, const FreeTrajectoryState &fts, const edm::EventSetup &es, bool apply_dzCut, const TrackingRegion &region, double dzcut) const
std::pair< ALIstring, ALIstring > pss
Definition: Fit.h:27
GlobalPoint const & origin() const
double DeltaPhiManual(const math::XYZVector &v1, const math::XYZVector &v2)
double_binary B
Definition: DDStreamer.cc:234
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
Global3DPoint GlobalPoint
Definition: GlobalPoint.h:10
T y() const
Definition: PV3DBase.h:63
float float float z
double getSqrEffectiveErrorOnZ(const TransientTrackingRecHit::ConstRecHitPointer &hit, const TrackingRegion &region)
double x
Definition: Quad.h:5
T sqrt(T t)
Definition: SSEVec.h:48
T z() const
Definition: PV3DBase.h:64
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
tuple IP
Definition: listHistos.py:76
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
const TrajectorySeed * buildSeed(TrajectorySeedCollection &seedCollection, const SeedingHitSet &hits, const FreeTrajectoryState &fts, const edm::EventSetup &es, bool apply_dzCut, const TrackingRegion &region) const
float originZBound() const
bounds the particle vertex in the longitudinal plane
int k[5][pyjets_maxn]
double simpleGetSlope(const TransientTrackingRecHit::ConstRecHitPointer &ohit, const TransientTrackingRecHit::ConstRecHitPointer &nohit, const TransientTrackingRecHit::ConstRecHitPointer &ihit, const TransientTrackingRecHit::ConstRecHitPointer &nihit, const TrackingRegion &region, double &cotTheta, double &z0)
XYZVectorD XYZVector
spatial vector with cartesian internal representation
Definition: Vector3D.h:30
const T & get() const
Definition: EventSetup.h:55
float ptMin() const
minimal pt of interest
double ptmin
Definition: HydjetWrapper.h:85
bool similarQuadExist(Quad &thisQuad, std::vector< Quad > &quadV)
unsigned int size() const
Definition: SeedingHitSet.h:41
tuple cout
Definition: gather_cfg.py:121
Definition: DDAxes.h:10
Definition: Quad.h:4
T x() const
Definition: PV3DBase.h:62
T get(const Candidate &c)
Definition: component.h:55
Global3DVector GlobalVector
Definition: GlobalVector.h:10
double SeedForPhotonConversionFromQuadruplets::verySimpleFit ( int  size,
double *  ax,
double *  ay,
double *  e2y,
double &  p0,
double &  e2p0,
double &  p1 
)

Definition at line 1015 of file SeedForPhotonConversionFromQuadruplets.cc.

Referenced by simpleGetSlope().

1015  {
1016 
1017  //#include "RecoTracker/ConversionSeedGenerators/interface/verySimpleFit.icc"
1018  return 0;
1019 }

Member Data Documentation

const int SeedForPhotonConversionFromQuadruplets::cotTheta_Max =99999
static

Definition at line 16 of file SeedForPhotonConversionFromQuadruplets.h.

Referenced by initialKinematic().

double SeedForPhotonConversionFromQuadruplets::kPI_
protected

Definition at line 105 of file SeedForPhotonConversionFromQuadruplets.h.

Referenced by DeltaPhiManual(), and trajectorySeed().

PrintRecoObjects SeedForPhotonConversionFromQuadruplets::po
protected

Definition at line 108 of file SeedForPhotonConversionFromQuadruplets.h.

Referenced by buildSeed(), and initialKinematic().

std::stringstream* SeedForPhotonConversionFromQuadruplets::pss
protected

Definition at line 107 of file SeedForPhotonConversionFromQuadruplets.h.

double SeedForPhotonConversionFromQuadruplets::theBOFFMomentum
protected

Definition at line 104 of file SeedForPhotonConversionFromQuadruplets.h.

Referenced by initialKinematic().

std::string SeedForPhotonConversionFromQuadruplets::thePropagatorLabel
protected

Definition at line 103 of file SeedForPhotonConversionFromQuadruplets.h.

Referenced by buildSeed(), and buildSeedBool().