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 SeedingHitSet::ConstRecHitPointer &hit, const TrackingRegion &region)
 
 SeedForPhotonConversionFromQuadruplets (const edm::ParameterSet &cfg, edm::ConsumesCollector &iC, const edm::ParameterSet &SeedComparitorPSet)
 
 SeedForPhotonConversionFromQuadruplets (edm::ConsumesCollector &iC, const edm::ParameterSet &SeedComparitorPSet, const std::string &propagator="PropagatorWithMaterial", double seedMomentumForBOFF=-5.0)
 
double simpleGetSlope (const SeedingHitSet::ConstRecHitPointer &ohit, const SeedingHitSet::ConstRecHitPointer &nohit, const SeedingHitSet::ConstRecHitPointer &ihit, const SeedingHitSet::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::Event &ev, const edm::EventSetup &es, std::stringstream &ss, std::vector< Quad > &quadV, 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 SeedingHitSet::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
 
SeedingHitSet::RecHitPointer refitHit (SeedingHitSet::ConstRecHitPointer hit, const TrajectoryStateOnSurface &state) const
 
bool similarQuadExist (Quad &thisQuad, std::vector< Quad > &quadV)
 

Protected Attributes

TkClonerImpl cloner
 
double kPI_
 
PrintRecoObjects po
 
std::stringstream * pss
 
double theBOFFMomentum
 
std::unique_ptr< SeedComparitortheComparitor
 
std::string thePropagatorLabel
 

Detailed Description

Definition at line 20 of file SeedForPhotonConversionFromQuadruplets.h.

Constructor & Destructor Documentation

SeedForPhotonConversionFromQuadruplets::SeedForPhotonConversionFromQuadruplets ( const edm::ParameterSet cfg,
edm::ConsumesCollector iC,
const edm::ParameterSet SeedComparitorPSet 
)

Definition at line 37 of file SeedForPhotonConversionFromQuadruplets.cc.

37  :
39  SeedComparitorPSet,
40  cfg.getParameter<std::string>("propagator"),
41  cfg.existsAs<double>("SeedMomentumForBOFF") ? cfg.getParameter<double>("SeedMomentumForBOFF") : 5.0)
42 {}
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(const edm::ParameterSet &cfg, edm::ConsumesCollector &iC, const edm::ParameterSet &SeedComparitorPSet)
SeedForPhotonConversionFromQuadruplets::SeedForPhotonConversionFromQuadruplets ( edm::ConsumesCollector iC,
const edm::ParameterSet SeedComparitorPSet,
const std::string &  propagator = "PropagatorWithMaterial",
double  seedMomentumForBOFF = -5.0 
)

Definition at line 44 of file SeedForPhotonConversionFromQuadruplets.cc.

References SurfaceDeformationFactory::create(), reco::get(), edm::ParameterSet::getParameter(), AlCaHLTBitMon_QueryRunRegistry::string, and theComparitor.

45  : thePropagatorLabel(propagator), theBOFFMomentum(seedMomentumForBOFF)
46 {
47  std::string comparitorName = SeedComparitorPSet.getParameter<std::string>("ComponentName");
48  if(comparitorName != "none") {
49  theComparitor.reset(SeedComparitorFactory::get()->create( comparitorName, SeedComparitorPSet, iC));
50  }
51 }
T getParameter(std::string const &) const
SurfaceDeformation * create(int type, const std::vector< double > &params)
T get(const Candidate &c)
Definition: component.h:55
SeedForPhotonConversionFromQuadruplets::~SeedForPhotonConversionFromQuadruplets ( )

Definition at line 53 of file SeedForPhotonConversionFromQuadruplets.cc.

53 {}

Member Function Documentation

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

Definition at line 982 of file SeedForPhotonConversionFromQuadruplets.cc.

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

982  {
983  bool swapped = true;
984  int j = 0;
986  while (swapped) {
987  swapped = false;
988  j++;
989  for (int i = 0; i < n - j; i++) {
990  if ( reco::deltaPhi( (arr[i]-vtx).phi(), (arr[i + 1]-vtx).phi() ) < 0. ) {
991  tmp = arr[i];
992  arr[i] = arr[i + 1];
993  arr[i + 1] = tmp;
994  swapped = true;
995  }
996  }
997  }
998 }
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 963 of file SeedForPhotonConversionFromQuadruplets.cc.

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

963  {
964  bool swapped = true;
965  int j = 0;
967  while (swapped) {
968  swapped = false;
969  j++;
970  for (int i = 0; i < n - j; i++) {
971  if ( reco::deltaPhi( (arr[i]-vtx).phi(), (arr[i + 1]-vtx).phi() ) > 0. ) {
972  tmp = arr[i];
973  arr[i] = arr[i + 1];
974  arr[i + 1] = tmp;
975  swapped = true;
976  }
977  }
978  }
979 }
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 689 of file SeedForPhotonConversionFromQuadruplets.cc.

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

Referenced by trajectorySeed().

696 {
697  // get tracker
699  es.get<TrackerDigiGeometryRecord>().get(tracker);
700 
701  // get propagator
702  edm::ESHandle<Propagator> propagatorHandle;
703  es.get<TrackingComponentsRecord>().get(thePropagatorLabel, propagatorHandle);
704  const Propagator* propagator = &(*propagatorHandle);
705 
706  // get updator
707  KFUpdator updator;
708 
709  // Now update initial state track using information from seed hits.
710 
711  TrajectoryStateOnSurface updatedState;
713 
714  const TrackingRecHit* hit = 0;
715  for ( unsigned int iHit = 0; iHit < hits.size() && iHit<2; iHit++) {
716  hit = hits[iHit];
717  TrajectoryStateOnSurface state = (iHit==0) ?
718  propagator->propagate(fts,tracker->idToDet(hit->geographicalId())->surface())
719  : propagator->propagate(updatedState, tracker->idToDet(hit->geographicalId())->surface());
720 
721  SeedingHitSet::ConstRecHitPointer tth = hits[iHit];
722 
723  SeedingHitSet::RecHitPointer newtth = refitHit( tth, state);
724 
725  updatedState = updator.update(state, *newtth);
726 
727  seedHits.push_back(newtth);
728 #ifdef mydebug_seed
729  uint32_t detid = hit->geographicalId().rawId();
730  (*pss) << "\n[SeedForPhotonConversionFromQuadruplets] hit " << iHit;
731  po.print(*pss, detid);
732  (*pss) << " " << detid << "\t lp " << hit->localPosition()
733  << " tth " << tth->localPosition() << " newtth " << newtth->localPosition() << " state " << state.globalMomentum().perp();
734 #endif
735  }
736 
737  PTrajectoryStateOnDet const & PTraj =
739 
740 
741  seedCollection.push_back( TrajectorySeed(PTraj,seedHits,alongMomentum));
742  return &seedCollection.back();
743 }
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:274
BaseTrackerRecHit const * ConstRecHitPointer
Definition: SeedingHitSet.h:11
void print(std::stringstream &ss, const SiStripCluster &clus)
const T & get() const
Definition: EventSetup.h:55
TrajectoryStateOnSurface propagate(STA const &state, SUR const &surface) const
Definition: Propagator.h:56
GlobalVector globalMomentum() const
unsigned int size() const
Definition: SeedingHitSet.h:44
DetId geographicalId() const
SeedingHitSet::RecHitPointer refitHit(SeedingHitSet::ConstRecHitPointer hit, const TrajectoryStateOnSurface &state) const
virtual LocalPoint localPosition() const GCC11_FINAL
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 749 of file SeedForPhotonConversionFromQuadruplets.cc.

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

Referenced by trajectorySeed().

757 {
758  // get tracker
760  es.get<TrackerDigiGeometryRecord>().get(tracker);
761 
762  // get propagator
763  edm::ESHandle<Propagator> propagatorHandle;
764  es.get<TrackingComponentsRecord>().get(thePropagatorLabel, propagatorHandle);
765  const Propagator* propagator = &(*propagatorHandle);
766 
767  // get cloner (FIXME: add to config)
768  auto TTRHBuilder = "WithTrackAngle";
770  es.get<TransientRecHitRecord>().get(TTRHBuilder, builderH);
771  auto builder = (TkTransientTrackingRecHitBuilder const *)(builderH.product());
772  cloner = (*builder).cloner();
773 
774  // get updator
775  KFUpdator updator;
776 
777  // Now update initial state track using information from seed hits.
778 
779  TrajectoryStateOnSurface updatedState;
781 
782  const TrackingRecHit* hit = 0;
783  for ( unsigned int iHit = 0; iHit < hits.size() && iHit<2; iHit++) {
784  hit = hits[iHit]->hit();
785  TrajectoryStateOnSurface state = (iHit==0) ?
786  propagator->propagate(fts,tracker->idToDet(hit->geographicalId())->surface())
787  : propagator->propagate(updatedState, tracker->idToDet(hit->geographicalId())->surface());
788  if (!state.isValid()) {
789  return false;}
790 
791  SeedingHitSet::ConstRecHitPointer tth = hits[iHit];
792 
793  SeedingHitSet::RecHitPointer newtth = refitHit( tth, state);
794 
795 
796  if (!checkHit(state,newtth,es)){
797  return false;
798  }
799 
800  updatedState = updator.update(state, *newtth);
801  if (!updatedState.isValid()){
802  return false;
803  }
804 
805  seedHits.push_back(newtth->hit()->clone());
806  }
807 
808  if(apply_dzCut){
810 
811  double zCA;
812 
813  math::XYZVector EstMomGam(updatedState.globalMomentum().x(),updatedState.globalMomentum().y(),updatedState.globalMomentum().z());
814  math::XYZVector EstPosGam(updatedState.globalPosition().x(),updatedState.globalPosition().y(),updatedState.globalPosition().z());
815 
816  double EstMomGamLength=std::sqrt(EstMomGam.x()*EstMomGam.x()+EstMomGam.y()*EstMomGam.y()+EstMomGam.z()*EstMomGam.z());
817  math::XYZVector EstMomGamNorm(EstMomGam.x()/EstMomGamLength,EstMomGam.y()/EstMomGamLength,EstMomGam.z()/EstMomGamLength);
818 
819  //Calculate dz of point of closest approach of the two lines (WA approach) -> cut on dz
820 
821 
822  const GlobalPoint EstPosGamGlobalPoint(updatedState.globalPosition().x(),updatedState.globalPosition().y(),updatedState.globalPosition().z());
823  const GlobalVector EstMomGamGlobalVector(updatedState.globalMomentum().x(),updatedState.globalMomentum().y(),updatedState.globalMomentum().z());
824 
825 
827  es.get<IdealMagneticFieldRecord>().get(bfield);
828  const MagneticField* magField = bfield.product();
829  TrackCharge qCharge = 0;
830 
831  const GlobalTrajectoryParameters myGlobalTrajectoryParameter(EstPosGamGlobalPoint, EstMomGamGlobalVector, qCharge, magField);
832 
833  AlgebraicSymMatrix66 aCovarianceMatrix;
834 
835  for (int i =0;i<6;++i)
836  for (int j =0;j<6;++j)
837  aCovarianceMatrix(i, j) = 1e-4;
838 
839  CartesianTrajectoryError myCartesianError (aCovarianceMatrix);
840 
841  const FreeTrajectoryState stateForProjectionToBeamLine(myGlobalTrajectoryParameter,myCartesianError);
842 
843  const GlobalPoint BeamSpotGlobalPoint(0,0,0);
844 
845  const reco::BeamSpot::Point BeamSpotPoint(region.origin().x(),region.origin().y(),region.origin().z());
846 
847  TSCBLBuilderNoMaterial tscblBuilder;
848 
849  double CovMatEntry=0.;
851  for (int i=0;i<3;++i) {
852  cov(i,i) = CovMatEntry;
853  }
855 
856  reco::BeamSpot myBeamSpot(BeamSpotPoint, 0.,0.,0.,0., cov,BeamType_);
857 
858  TrajectoryStateClosestToBeamLine tscbl = tscblBuilder(stateForProjectionToBeamLine,myBeamSpot);
859  if (tscbl.isValid()==false) {
860  zCA=0;
861  }
862  else{
863  GlobalPoint v = tscbl.trackStateAtPCA().position(); // Position of closest approach to BS
864  zCA=v.z();
865  }
866 
867 /* //Calculate dz of point of closest approach of the two lines -> cut on dz
868 
869  double newX,newY,newR;
870  double Rbuff=std::sqrt(EstPosGam.x()*EstPosGam.x()+EstPosGam.y()*EstPosGam.y());
871  double deltas,s,sbuff;
872  double rMin=1e9;
873 
874  double InitXp=EstPosGam.x()+1*EstMomGamNorm.x();
875  double InitXm=EstPosGam.x()-1*EstMomGamNorm.x();
876  double InitYp=EstPosGam.y()+1*EstMomGamNorm.y();
877  double InitYm=EstPosGam.y()-1*EstMomGamNorm.y();
878 
879  if(InitXp*InitXp+InitYp*InitYp < InitXm*InitXm+InitYm*InitYm) {s=5; deltas=5;}
880  else {s=-5; deltas=-5;}
881 
882  int nTurns=0;
883  int nIterZ=0;
884  for (int i_dz=0;i_dz<1000;i_dz++){
885  newX=EstPosGam.x()+s*EstMomGamNorm.x();
886  newY=EstPosGam.y()+s*EstMomGamNorm.y();
887  newR=std::sqrt(newX*newX+newY*newY);
888  if(newR>Rbuff) {deltas=-1*deltas/10;nTurns++;}
889  else {Rbuff=newR;}
890  if(newR<rMin) {rMin=newR; sbuff=s;}
891  s=s+deltas;
892  nIterZ++;
893  if(nTurns>1) break;
894  }
895 
896  zCA=EstPosGam.z()+sbuff*EstMomGamNorm.z();
897 */
898 
899 #ifdef mydebug_knuenz
900  std::cout<< "zCA: " << zCA <<std::endl;
901 #endif
902 
903  if(std::fabs(zCA)>dzcut) return false;
904 
905 
906  }
907 
908 
909  return true;
910 }
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
BaseTrackerRecHit const * hit() const GCC11_FINAL
math::XYZPoint Point
point in the space
Definition: BeamSpot.h:29
int TrackCharge
Definition: TrackCharge.h:4
void push_back(D *&d)
Definition: OwnVector.h:274
BaseTrackerRecHit const * ConstRecHitPointer
Definition: SeedingHitSet.h:11
T sqrt(T t)
Definition: SSEVec.h:48
T z() const
Definition: PV3DBase.h:64
int j
Definition: DBlmapReader.cc:9
bool checkHit(const TrajectoryStateOnSurface &, const SeedingHitSet::ConstRecHitPointer &hit, const edm::EventSetup &es) const
virtual TrackingRecHit * clone() const =0
GlobalPoint position() const
XYZVectorD XYZVector
spatial vector with cartesian internal representation
Definition: Vector3D.h:30
const T & get() const
Definition: EventSetup.h:55
T const * product() const
Definition: ESHandle.h:62
TrajectoryStateOnSurface propagate(STA const &state, SUR const &surface) const
Definition: Propagator.h:56
GlobalVector globalMomentum() const
unsigned int size() const
Definition: SeedingHitSet.h:44
tuple cout
Definition: gather_cfg.py:121
DetId geographicalId() const
SeedingHitSet::RecHitPointer refitHit(SeedingHitSet::ConstRecHitPointer hit, const TrajectoryStateOnSurface &state) const
T x() const
Definition: PV3DBase.h:62
bool SeedForPhotonConversionFromQuadruplets::checkHit ( const TrajectoryStateOnSurface ,
const SeedingHitSet::ConstRecHitPointer hit,
const edm::EventSetup es 
) const
inlineprotected

Definition at line 62 of file SeedForPhotonConversionFromQuadruplets.h.

Referenced by buildSeedBool().

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

Definition at line 1093 of file SeedForPhotonConversionFromQuadruplets.cc.

References kPI_.

Referenced by trajectorySeed().

1093  {
1094 
1095  double kTWOPI = 2.*kPI_;
1096  double DeltaPhiMan=v1.Phi()-v2.Phi();
1097  while (DeltaPhiMan >= kPI_) DeltaPhiMan -= kTWOPI;
1098  while (DeltaPhiMan < -kPI_) DeltaPhiMan += kTWOPI;
1099 
1100  return DeltaPhiMan;
1101 
1102 }
double SeedForPhotonConversionFromQuadruplets::getSqrEffectiveErrorOnZ ( const SeedingHitSet::ConstRecHitPointer hit,
const TrackingRegion region 
)

Definition at line 1045 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().

1045  {
1046 
1047  //
1048  //Fit-wise the effective error on Z is the sum in quadrature of the error on Z
1049  //and the error on R correctly projected by using hit-vertex direction
1050 
1051  double sqrProjFactor = sqr((hit->globalPosition().z()-region.origin().z())/(hit->globalPosition().perp()-region.origin().perp()));
1052  return (hit->globalPositionError().czz()+sqrProjFactor*hit->globalPositionError().rerr(hit->globalPosition()));
1053 
1054 }
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 660 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().

664 {
665  // Set initial uncertainty on track parameters, using only P.V. constraint and no hit
666  // information.
667  GlobalError vertexErr( sqr(vertexBounds.x()), 0,
668  sqr(vertexBounds.y()), 0, 0,
669  sqr(vertexBounds.z())
670  );
671 
672 
673  AlgebraicSymMatrix55 C = ROOT::Math::SMatrixIdentity();
674 
675 // FIXME: minC00. Prevent apriori uncertainty in 1/P from being too small,
676 // to avoid instabilities.
677 // N.B. This parameter needs optimising ...
678  float sin2th = sqr(sinTheta);
679  float minC00 = 1.0;
680  C[0][0] = std::max(sin2th/sqr(ptMin), minC00);
681  float zErr = vertexErr.czz();
682  float transverseErr = vertexErr.cxx(); // assume equal cxx cyy
683  C[3][3] = transverseErr;
684  C[4][4] = zErr*sin2th + transverseErr*(1-sin2th);
685 
686  return CurvilinearTrajectoryError(C);
687 }
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 601 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().

606 {
608 
609  SeedingHitSet::ConstRecHitPointer tth1 = hits[0];
610  SeedingHitSet::ConstRecHitPointer tth2 = hits[1];
611 
613  es.get<IdealMagneticFieldRecord>().get(bfield);
614  float nomField = bfield->nominalValue();
615  bool isBOFF = (0==nomField);
616 
617 
618  FastHelix helix(tth2->globalPosition(), tth1->globalPosition(), vertexPos, nomField,&*bfield, vertexPos);
619  kine = helix.stateAtVertex();
620 
621  //force the pz/pt equal to the measured one
622  if(fabs(cotTheta)<cotTheta_Max)
623  kine = GlobalTrajectoryParameters(kine.position(),
624  GlobalVector(kine.momentum().x(),kine.momentum().y(),kine.momentum().perp()*cotTheta),
625  kine.charge(),
626  & kine.magneticField()
627  );
628  else
629  kine = GlobalTrajectoryParameters(kine.position(),
630  GlobalVector(kine.momentum().x(),kine.momentum().y(),kine.momentum().perp()*cotTheta_Max),
631  kine.charge(),
632  & kine.magneticField()
633  );
634 
635 #ifdef mydebug_seed
636  uint32_t detid;
637  (*pss) << "[SeedForPhotonConversionFromQuadruplets] initialKinematic tth1 " ;
638  detid=tth1->geographicalId().rawId();
639  po.print(*pss, detid );
640  (*pss) << " \t " << detid << " " << tth1->localPosition() << " " << tth1->globalPosition() ;
641  detid= tth2->geographicalId().rawId();
642  (*pss) << " \n\t tth2 ";
643  po.print(*pss, detid );
644  (*pss) << " \t " << detid << " " << tth2->localPosition() << " " << tth2->globalPosition()
645  << "\nhelix momentum " << kine.momentum() << " pt " << kine.momentum().perp() << " radius " << 1/kine.transverseCurvature();
646 #endif
647 
648  if (isBOFF && (theBOFFMomentum > 0)) {
649  kine = GlobalTrajectoryParameters(kine.position(),
650  kine.momentum().unit() * theBOFFMomentum,
651  kine.charge(),
652  &*bfield);
653  }
654  return kine;
655 }
std::pair< ALIstring, ALIstring > pss
Definition: Fit.h:27
T perp() const
Definition: PV3DBase.h:72
T y() const
Definition: PV3DBase.h:63
BaseTrackerRecHit const * ConstRecHitPointer
Definition: SeedingHitSet.h:11
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
SeedingHitSet::RecHitPointer SeedForPhotonConversionFromQuadruplets::refitHit ( SeedingHitSet::ConstRecHitPointer  hit,
const TrajectoryStateOnSurface state 
) const
protected

Definition at line 914 of file SeedForPhotonConversionFromQuadruplets.cc.

References cloner.

Referenced by buildSeed(), and buildSeedBool().

917 {
918  //const TransientTrackingRecHit* a= hit.get();
919  //return const_cast<TransientTrackingRecHit*> (a);
920  //This was modified otherwise the rechit will have just the local x component and local y=0
921  // To understand how to modify for pixels
922 
923  //const TSiStripRecHit2DLocalPos* b = dynamic_cast<const TSiStripRecHit2DLocalPos*>(a);
924  //return const_cast<TSiStripRecHit2DLocalPos*>(b);
925  return (SeedingHitSet::RecHitPointer)(cloner(*hit,state));
926 }
bool SeedForPhotonConversionFromQuadruplets::similarQuadExist ( Quad thisQuad,
std::vector< Quad > &  quadV 
)
protected

Definition at line 1056 of file SeedForPhotonConversionFromQuadruplets.cc.

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

Referenced by trajectorySeed().

1056  {
1057 
1058  BOOST_FOREACH( Quad quad, quadV )
1059  {
1060  double dx = thisQuad.x-quad.x;
1061  double dy = thisQuad.y-quad.y;
1062  //double dz = abs(thisQuad.z-quad.z);
1063  //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
1064  //std::cout<<"x1-x2="<<dx<<"y1-y2="<<dy<<"dx*dx+dy*dy="<<dx*dx+dy*dy<<std::endl; //ValDebug
1065  //std::cout<<"thisQuad.ptPlus-quad.ptPlus="<<thisQuad.ptPlus-quad.ptPlus<<"abs(thisQuad.ptPlus-quad.ptPlus)="<<abs(thisQuad.ptPlus-quad.ptPlus)<<std::endl; //ValDebug
1066  //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
1067  //if ( sqrt(dx*dx+dy*dy)<1.)
1068  // 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
1069  // ( sqrt(dx*dx+dy*dy)<1. &&
1070  //z<3. &&
1071  //bs(thisQuad.ptPlus-quad.ptPlus)<0.5*quad.ptPlus &&
1072  //bs(thisQuad.ptMinus-quad.ptMinus)<0.5*quad.ptMinus &&
1073  //bs(thisQuad.cot-quad.cot)<0.3*quad.cot
1074  //
1075  // {
1076  // //std::cout<<"Seed rejected due to arbitration"<<std::endl;
1077  // return true;
1078  // }
1079  if ( sqrt(dx*dx+dy*dy)<1. &&
1080  fabs(thisQuad.ptPlus-quad.ptPlus)<0.5*quad.ptPlus &&
1081  fabs(thisQuad.ptMinus-quad.ptMinus)<0.5*quad.ptMinus
1082  )
1083  {
1084  //std::cout<<"Seed rejected due to arbitration"<<std::endl;
1085  return true;
1086  }
1087  }
1088  quadV.push_back(thisQuad);
1089  return false;
1090 
1091 }
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 SeedingHitSet::ConstRecHitPointer ohit,
const SeedingHitSet::ConstRecHitPointer nohit,
const SeedingHitSet::ConstRecHitPointer ihit,
const SeedingHitSet::ConstRecHitPointer nihit,
const TrackingRegion region,
double &  cotTheta,
double &  z0 
)

Definition at line 1002 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().

1002  {
1003 
1004  double x[5], y[5], e2y[5];
1005 
1006  //The fit is done filling x with r values, y with z values of the four hits and the vertex
1007  //
1008  //Hits
1009  x[0] = ohit->globalPosition().perp();
1010  y[0] = ohit->globalPosition().z();
1011  e2y[0] = getSqrEffectiveErrorOnZ(ohit, region);
1012  //
1013  x[1] = nohit->globalPosition().perp();
1014  y[1] = nohit->globalPosition().z();
1015  e2y[1] = getSqrEffectiveErrorOnZ(nohit, region);
1016  //
1017  x[2] = nihit->globalPosition().perp();
1018  y[2] = nihit->globalPosition().z();
1019  e2y[2] = getSqrEffectiveErrorOnZ(nihit, region);
1020  //
1021  x[3] = ihit->globalPosition().perp();
1022  y[3] = ihit->globalPosition().z();
1023  e2y[3] = getSqrEffectiveErrorOnZ(ihit, region);
1024  //
1025  //Vertex
1026  x[4] = region.origin().perp();
1027  y[4] = region.origin().z();
1028  double vError = region.originZBound();
1029  if ( vError > 15. ) vError = 1.;
1030  e2y[4] = sqr(vError);
1031 
1032  double e2z0;
1033  double chi2 = verySimpleFit(5, x, y, e2y, z0, e2z0, cotTheta);
1034 
1035  return chi2;
1036 
1037 }
T perp() const
Definition: PV3DBase.h:72
GlobalPoint const & origin() const
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
double getSqrEffectiveErrorOnZ(const SeedingHitSet::ConstRecHitPointer &hit, const TrackingRegion &region)
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 933 of file SeedForPhotonConversionFromQuadruplets.cc.

References i.

933  {
934  (*pss) << "\n" << s << "\t";
935  for(size_t i=0;i<2;++i)
936  (*pss) << std::setw (60) << d[i] << std::setw(1) << " | ";
937 }
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 940 of file SeedForPhotonConversionFromQuadruplets.cc.

References i.

940  {
941  (*pss) << "\n" << s << "\t";
942  for(size_t i=0;i<2;++i)
943  (*pss) << std::setw (60) << d[i] << std::setw(1) << " | ";
944 }
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 947 of file SeedForPhotonConversionFromQuadruplets.cc.

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

947  {
948  (*pss) << "\n" << s << "\t";
949  for(size_t i=0;i<2;++i)
950  (*pss) << std::setw(20) << d[i] << " r " << d[i].perp() << " phi " << d[i].phi() << " | ";
951 }
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 954 of file SeedForPhotonConversionFromQuadruplets.cc.

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

954  {
955  (*pss) << "\n" << s << "\n";
956  for(int i=0;i<n;++i)
957  (*pss) << std::setw(20) << d[i] << " r " << d[i].perp() << " phi " << d[i].phi() << "\n";
958 }
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::Event ev,
const edm::EventSetup es,
std::stringstream &  ss,
std::vector< Quad > &  quadV,
edm::ParameterSet QuadCutPSet 
)

Definition at line 56 of file SeedForPhotonConversionFromQuadruplets.cc.

References funct::abs(), buildSeed(), buildSeedBool(), funct::C, Conv4HitsReco2::ConversionCandidate(), funct::cos(), gather_cfg::cout, DeltaPhiManual(), cond::rpcobgas::detid, Conv4HitsReco2::Dump(), 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(), contentValuesCheck::ss, theComparitor, Quad::x, x, PV3DBase< T, PVType, FrameType >::x(), detailsBasic3DVector::y, PV3DBase< T, PVType, FrameType >::y(), detailsBasic3DVector::z, and PV3DBase< T, PVType, FrameType >::z().

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

Referenced by simpleGetSlope().

1039  {
1040 
1041  //#include "RecoTracker/ConversionSeedGenerators/interface/verySimpleFit.icc"
1042  return 0;
1043 }

Member Data Documentation

TkClonerImpl SeedForPhotonConversionFromQuadruplets::cloner
mutableprotected

Definition at line 110 of file SeedForPhotonConversionFromQuadruplets.h.

Referenced by buildSeedBool(), and refitHit().

const int SeedForPhotonConversionFromQuadruplets::cotTheta_Max =99999
static

Definition at line 22 of file SeedForPhotonConversionFromQuadruplets.h.

Referenced by initialKinematic().

double SeedForPhotonConversionFromQuadruplets::kPI_
protected

Definition at line 107 of file SeedForPhotonConversionFromQuadruplets.h.

Referenced by DeltaPhiManual(), and trajectorySeed().

PrintRecoObjects SeedForPhotonConversionFromQuadruplets::po
protected

Definition at line 113 of file SeedForPhotonConversionFromQuadruplets.h.

Referenced by buildSeed(), and initialKinematic().

std::stringstream* SeedForPhotonConversionFromQuadruplets::pss
protected

Definition at line 112 of file SeedForPhotonConversionFromQuadruplets.h.

double SeedForPhotonConversionFromQuadruplets::theBOFFMomentum
protected

Definition at line 106 of file SeedForPhotonConversionFromQuadruplets.h.

Referenced by initialKinematic().

std::unique_ptr<SeedComparitor> SeedForPhotonConversionFromQuadruplets::theComparitor
protected
std::string SeedForPhotonConversionFromQuadruplets::thePropagatorLabel
protected

Definition at line 105 of file SeedForPhotonConversionFromQuadruplets.h.

Referenced by buildSeed(), and buildSeedBool().