55 for (
unsigned i=0;
i<seedingAlgo.size(); ++
i )
56 produces<TrajectorySeedCollection>(seedingAlgo[
i]);
60 if (
pTMin.size() != seedingAlgo.size() )
62 <<
" WARNING : pTMin does not have the proper size "
65 for (
unsigned i=0; i<
pTMin.size(); ++
i )
72 <<
" WARNING : minRecHits does not have the proper size "
76 for (
unsigned ialgo=0; ialgo<
minRecHits.size(); ++ialgo )
81 if (
maxD0.size() != seedingAlgo.size() )
83 <<
" WARNING : maxD0 does not have the proper size "
87 if (
maxZ0.size() != seedingAlgo.size() )
89 <<
" WARNING : maxZ0 does not have the proper size "
100 if ( numberOfHits.size() != seedingAlgo.size() )
102 <<
" WARNING : numberOfHits does not have the proper size "
114 for(std::vector<std::string>::const_iterator it=layerList.begin(); it < layerList.end(); ++it) {
115 std::vector<LayerSpec> tempResult;
118 while (pos != std::string::npos) {
123 layerSpec.
name = line.substr(0, pos);;
126 if (layerSpec.
name.substr(0,4)==
"BPix" ) {
129 layerSpec.
idLayer = std::atoi(layerSpec.
name.substr(4,1).c_str());
132 <<
"Bad data naming in IterativeInitialStep_cff.py iterativeInitialSeeds.layerList" << std::endl
133 <<
"Layers: " << layerSpec.
name <<
", number needs to be in range 1-3" << std::endl;
137 else if (layerSpec.
name.substr(0,4)==
"FPix" ) {
139 if(layerSpec.
name.substr(layerSpec.
name.size()-3)==
"pos")
143 layerSpec.
idLayer = std::atoi(layerSpec.
name.substr(4,1).c_str());
146 <<
"Bad data naming in IterativeInitialStep_cff.py iterativeInitialSeeds.layerList" << std::endl
147 <<
"Layers: " << layerSpec.
name <<
", number needs to be in range 1-2" << std::endl;
151 else if (layerSpec.
name.substr(0,3)==
"TIB" ) {
154 layerSpec.
idLayer = std::atoi(layerSpec.
name.substr(3,1).c_str());
157 <<
"Bad data naming in IterativeInitialStep_cff.py iterativeInitialSeeds.layerList" << std::endl
158 <<
"Layers: " << layerSpec.
name <<
", number needs to be in range 1-4" << std::endl;
162 else if (layerSpec.
name.substr(0,3)==
"TID" ) {
164 if(layerSpec.
name.substr(layerSpec.
name.size()-3)==
"pos")
168 layerSpec.
idLayer = std::atoi(layerSpec.
name.substr(3,1).c_str());
171 <<
"Bad data naming in IterativeInitialStep_cff.py iterativeInitialSeeds.layerList" << std::endl
172 <<
"Layers: " << layerSpec.
name <<
", number needs to be in range 1-3" << std::endl;
176 else if (layerSpec.
name.substr(0,3)==
"TOB" ) {
179 layerSpec.
idLayer = std::atoi(layerSpec.
name.substr(3,1).c_str());
182 <<
"Bad data naming in IterativeInitialStep_cff.py iterativeInitialSeeds.layerList" << std::endl
183 <<
"Layers: " << layerSpec.
name <<
", number needs to be in range 1-6" << std::endl;
187 else if (layerSpec.
name.substr(0,3)==
"TEC" ) {
189 if(layerSpec.
name.substr(layerSpec.
name.size()-3)==
"pos")
193 layerSpec.
idLayer = std::atoi(layerSpec.
name.substr(3,1).c_str());
196 <<
"Bad data naming in IterativeInitialStep_cff.py iterativeInitialSeeds.layerList" << std::endl
197 <<
"Layers: " << layerSpec.
name <<
", number needs to be in range 1-9" << std::endl;
202 <<
"Bad data naming in IterativeInitialStep_cff.py iterativeInitialSeeds.layerList" << std::endl
203 <<
"Layer: " << layerSpec.
name <<
", shouldn't exist" << std::endl
204 <<
"Case sensitive names: BPix FPix TIB TID TOB TEC" << std::endl;
207 tempResult.push_back(layerSpec);
208 line=line.substr(pos+1,std::string::npos);
223 conf.
getParameter<std::vector<unsigned int> >(
"firstHitSubDetectorNumber");
226 <<
" WARNING : firstHitSubDetectorNumber does not have the proper size "
229 std::vector<unsigned int> firstSubDets =
230 conf.
getParameter<std::vector<unsigned int> >(
"firstHitSubDetectors");
240 if ( firstSubDets.size() != check1 )
242 <<
" WARNING : firstHitSubDetectors does not have the proper size (should be " << check1 <<
")"
247 conf.
getParameter<std::vector<unsigned int> >(
"secondHitSubDetectorNumber");
250 <<
" WARNING : secondHitSubDetectorNumber does not have the proper size "
253 std::vector<unsigned int> secondSubDets =
254 conf.
getParameter<std::vector<unsigned int> >(
"secondHitSubDetectors");
264 if ( secondSubDets.size() != check2 )
266 <<
" WARNING : secondHitSubDetectors does not have the proper size (should be " << check2 <<
")"
270 conf.
getParameter<std::vector<unsigned int> >(
"thirdHitSubDetectorNumber");
273 <<
" WARNING : thirdHitSubDetectorNumber does not have the proper size "
276 std::vector<unsigned int> thirdSubDets =
277 conf.
getParameter<std::vector<unsigned int> >(
"thirdHitSubDetectors");
287 if ( thirdSubDets.size() != check3 )
289 <<
" WARNING : thirdHitSubDetectors does not have the proper size (should be " << check3 <<
")"
296 <<
" WARNING : originRadius does not have the proper size "
302 <<
" WARNING : originHalfLength does not have the proper size "
308 <<
" WARNING : originpTMin does not have the proper size "
314 <<
" WARNING : primaryVertices does not have the proper size "
320 <<
" WARNING : zVertexConstraint does not have the proper size "
332 std::cout <<
"TrajectorySeedProducer destructed" << std::endl;
370 std::cout <<
"################################################################" << std::endl;
371 std::cout <<
" TrajectorySeedProducer produce init " << std::endl;
380 unsigned nSimTracks = 0;
381 unsigned nTracksWithHits = 0;
382 unsigned nTracksWithPT = 0;
383 unsigned nTracksWithD0Z0 = 0;
388 std::vector<TrajectorySeedCollection*>
390 for (
unsigned ialgo=0; ialgo<
seedingAlgo.size(); ++ialgo ) {
405 x0 = BSPosition_.X();
406 y0 = BSPosition_.Y();
407 z0 = BSPosition_.Z();
417 std::cout <<
" Step A: SimTracks found " << theSimTracks->size() << std::endl;
426 std::cout <<
" Step B: Full GS RecHits found " << theGSRecHits->size() << std::endl;
428 if(theGSRecHits->size() == 0) {
429 for (
unsigned ialgo=0; ialgo<
seedingAlgo.size(); ++ialgo ) {
430 std::auto_ptr<TrajectorySeedCollection>
p(output[ialgo]);
437 vertices = std::vector<const reco::VertexCollection*>
439 for (
unsigned ialgo=0; ialgo<
seedingAlgo.size(); ++ialgo ) {
444 if (!isVertexCollection )
continue;
449 std::cout <<
" Step C: Loop over the RecHits, track by track " << std::endl;
453 const std::vector<unsigned> theSimTrackIds = theGSRecHits->ids();
456 for (
unsigned tkId=0; tkId != theSimTrackIds.size(); ++tkId ) {
459 std::cout <<
"Track number " << tkId << std::endl;
463 unsigned simTrackId = theSimTrackIds[tkId];
464 const SimTrack& theSimTrack = (*theSimTracks)[simTrackId];
467 <<
" eta " << theSimTrack.
momentum().Eta()
468 <<
" pdg ID " << theSimTrack.
type()
476 int vertexIndex = theSimTrack.
vertIndex();
477 const SimVertex& theSimVertex = (*theSimVtx)[vertexIndex];
479 std::cout <<
" o SimTrack " << theSimTrack << std::endl;
480 std::cout <<
" o SimVertex " << theSimVertex << std::endl;
505 unsigned numberOfRecHits = 0;
507 for ( iterRecHit = theRecHitRangeIteratorBegin;
508 iterRecHit != theRecHitRangeIteratorEnd;
510 previousHit = currentHit;
518 for (
unsigned int ialgo = 0; ialgo <
seedingAlgo.size(); ++ialgo ) {
525 std::cout <<
"The number of RecHits = " << numberOfRecHits << std::endl;
527 if ( numberOfRecHits <
minRecHits[ialgo] )
continue;
531 if ( theSimTrack.
momentum().Perp2() <
pTMin[ialgo] )
continue;
539 std::vector<TrackerRecHit >
546 bool compatible =
false;
548 for ( iterRecHit1 = theRecHitRangeIteratorBegin; iterRecHit1 != theRecHitRangeIteratorEnd; ++iterRecHit1) {
553 std::cout <<
"The first hit subDetId = " << theSeedHits0.
subDetId() << std::endl;
557 bool isInside =
true;
565 if ( isInside )
continue;
579 if ( !isOndet )
continue;
583 std::cout <<
"Apparently the first hit is on the requested detector! " << std::endl;
585 for ( iterRecHit2 = iterRecHit1+1; iterRecHit2 != theRecHitRangeIteratorEnd; ++iterRecHit2) {
589 std::cout <<
"The second hit subDetId = " << theSeedHits1.
subDetId() << std::endl;
599 if ( isInside )
continue;
603 isOndet = theSeedHits[0].isOnRequestedDet(
theLayersInSets, theSeedHits[1]);
606 if ( !isOndet )
break;
612 std::cout <<
"Apparently the second hit is on the requested detector! " << std::endl;
628 std::cout <<
"Algo" <<
seedingAlgo[0] <<
"\t Are the two hits compatible with the PV? " << compatible << std::endl;
636 compatible = compatible && theSeedHits[0].makesAPairWith3rd(theSeedHits[1]);
638 compatible = compatible && theSeedHits[0].makesAPairWith(theSeedHits[1]);
651 if ( !compatible )
continue;
662 for ( iterRecHit3 = iterRecHit2+1; iterRecHit3 != theRecHitRangeIteratorEnd; ++iterRecHit3) {
667 std::cout <<
"The third hit subDetId = " << theSeedHits2.
subDetId() << std::endl;
677 if ( isInside )
continue;
681 isOndet = theSeedHits[0].isOnRequestedDet(
theLayersInSets, theSeedHits[1], theSeedHits[2]);
685 if ( !isOndet )
continue;
692 if (!
selectMuons) compatible = compatible && theSeedHits[0].makesATripletWith(theSeedHits[1],theSeedHits[2]);
696 std::cout <<
"Apparently the third hit is on the requested detector! " << std::endl;
699 if ( compatible )
break;
703 if ( compatible )
break;
707 if ( compatible )
break;
713 if ( !compatible )
continue;
716 std::cout <<
"@@@ There is at least a compatible seed" << std::endl;
718 std::cout <<
"@@@ There is no compatible seed" << std::endl;
723 std::cout <<
"Preparing to create the TrajectorySeed" << std::endl;
728 for (
unsigned ih=0; ih<theSeedHits.size(); ++ih ) {
733 std::cout <<
"with " << recHits.
size() <<
" hits." << std::endl;
739 (*theSimVtx)[vertexIndex].
position().
y(),
740 (*theSimVtx)[vertexIndex].
position().
z());
743 GlobalVector momentum( (*theSimTracks)[simTrackId].momentum().
x() ,
744 (*theSimTracks)[simTrackId].momentum().
y() ,
745 (*theSimTracks)[simTrackId].momentum().
z() );
747 float charge = (*theSimTracks)[simTrackId].charge();
756 if(theSeedHits0.
subDetId() !=1 || theSeedHits0.
subDetId() !=2) errorMatrix = errorMatrix * 0.0000001;
761 std::cout <<
"TrajectorySeedProducer: SimTrack parameters " << std::endl;
762 std::cout <<
"\t\t pT = " << (*theSimTracks)[simTrackId].momentum().Pt() << std::endl;
763 std::cout <<
"\t\t eta = " << (*theSimTracks)[simTrackId].momentum().Eta() << std::endl;
764 std::cout <<
"\t\t phi = " << (*theSimTracks)[simTrackId].momentum().Phi() << std::endl;
765 std::cout <<
"TrajectorySeedProducer: AlgebraicSymMatrix " << errorMatrix << std::endl;
771 std::cout <<
"TrajectorySeedProducer: FTS momentum " << initialFTS.
momentum() << std::endl;
780 if (!initialTSOS.
isValid())
continue;
790 std::cout <<
"TrajectorySeedProducer: TSOS surface side " << initialTSOS.
surfaceSide() << std::endl;
798 std::cout <<
"Trajectory seed created ! " << std::endl;
806 for (
unsigned ialgo=0; ialgo<
seedingAlgo.size(); ++ialgo ) {
807 std::auto_ptr<TrajectorySeedCollection>
p(output[ialgo]);
829 float localErrors[15];
831 for (
int i=0;
i<dim; ++
i) {
832 for (
int j=0;
j<=
i; ++
j) {
833 localErrors[k++] =
m(
i,
j);
836 int surfaceSide =
static_cast<int>(ts.
surfaceSide());
848 unsigned algo)
const {
856 std::cout <<
"ThePos1 = " << thePos1 << std::endl;
857 std::cout <<
"ThePos2 = " << thePos2 << std::endl;
875 if ( !intersect )
return false;
878 std::cout <<
"MyPart R = " << myPart.
R() <<
"\t Z = " << myPart.
Z()
879 <<
"\t pT = " << myPart.Pt() << std::endl;
884 if ( myPart.Pt() <
originpTMin[algo] )
return false;
891 if (!theVertices)
return true;
892 unsigned nVertices = theVertices->size();
896 + (thePos1.Y()-
y0)*(thePos1.Y()-
y0) );
898 + (thePos2.Y()-
y0)*(thePos2.Y()-
y0) );
900 for (
unsigned iv=0; iv<nVertices; ++iv ) {
902 double zV = (*theVertices)[iv].z();
904 double checkRZ1 = forward ?
907 double checkRZ2 = forward ?
913 bool compat = forward ?
914 checkRZmin < R1 && R1 < checkRZmax :
915 checkRZmin < thePos1.Z()-zV && thePos1.Z()-zV < checkRZmax;
917 if ( compat )
return compat;
bool compatibleWithBeamAxis(GlobalPoint &gpos1, GlobalPoint &gpos2, double error, bool forward, unsigned algo) const
void setCharge(float q)
set the MEASURED charge
T getParameter(std::string const &) const
std::vector< unsigned > minRecHits
std::pair< const_iterator, const_iterator > range
iterator range
double zImpactParameter(double x0=0, double y0=0.) const
Longitudinal impact parameter.
std::vector< unsigned int > firstHitSubDetectorNumber
virtual void beginRun(edm::Run const &run, const edm::EventSetup &es) override
std::vector< std::vector< LayerSpec > > theLayersInSets
TrajectorySeedProducer(const edm::ParameterSet &conf)
const LocalTrajectoryParameters & localParameters() const
ROOT::Math::SMatrixIdentity AlgebraicMatrixID
std::vector< double > pTMin
Geom::Phi< T > phi() const
bool propagateToBeamCylinder(const XYZTLorentzVector &v, double radius=0.)
std::vector< Vertex > VertexCollection
collection of Vertex objects
std::vector< double > originRadius
ROOT::Math::SMatrix< double, 5, 5, ROOT::Math::MatRepSym< double, 5 > > AlgebraicSymMatrix55
const MagneticFieldMap * theFieldMap
void stateOnDet(const TrajectoryStateOnSurface &ts, unsigned int detid, PTrajectoryStateOnDet &pts) const
A mere copy (without memory leak) of an existing tracking method.
std::vector< double > originpTMin
const Plane & surface() const
The nominal surface of the GeomDet.
GlobalPoint globalPosition() const
The global position.
LocalError positionError() const
The Signals That Services Can Subscribe To This is based on ActivityRegistry and is current per Services can connect to the signals distributed by the ActivityRegistry in order to monitor the activity of the application Each possible callback has some defined which we here list in angle e g
static int position[TOTALCHAMBERS][3]
edm::InputTag hitProducer
bool isOnTheSameLayer(const TrackerRecHit &other) const
Check if two hits are on the same layer of the same subdetector.
std::vector< unsigned int > numberOfHits
uint32_t rawId() const
get the raw id
LocalVector localMomentum() const
PropagatorWithMaterial * thePropagator
double R() const
vertex radius
static const double pts[33]
unsigned int absMinRecHits
C::const_iterator const_iterator
constant access iterator type
std::vector< std::vector< unsigned int > > secondHitSubDetectors
virtual void produce(edm::Event &e, const edm::EventSetup &es) override
std::vector< TrajectorySeed > TrajectorySeedCollection
const T & max(const T &a, const T &b)
std::vector< std::vector< unsigned int > > thirdHitSubDetectors
OrphanHandle< PROD > put(std::auto_ptr< PROD > product)
Put a new product.
std::vector< double > zVertexConstraint
SurfaceSide surfaceSide() const
Position relative to material, defined relative to momentum vector.
virtual ~TrajectorySeedProducer()
double Z() const
z of vertex
Abs< T >::type abs(const T &t)
const AlgebraicSymMatrix55 & matrix() const
const math::XYZTLorentzVectorD & position() const
const LocalTrajectoryError & localError() const
const TrackerGeometry * theGeometry
virtual const GeomDet * idToDet(DetId) const
bool getByLabel(InputTag const &tag, Handle< PROD > &result) const
GlobalVector momentum() const
std::vector< unsigned int > thirdHitSubDetectorNumber
virtual TrackingRecHit * clone() const =0
int vertIndex() const
index of the vertex in the Event container (-1 if no vertex)
bool isForward() const
Is it a forward hit ?
XYZPointD XYZPoint
point in space with cartesian internal representation
std::vector< unsigned int > secondHitSubDetectorNumber
T const * product() const
unsigned int layerNumber() const
The Layer Number.
edm::InputTag theBeamSpot
ESHandle< TrackerGeometry > geometry
int type() const
particle type (HEP PDT convension)
bool isOnRequestedDet(const std::vector< unsigned int > &whichDet, const std::string &seedingAlgo) const
Check if the hit is on one of the requested detector.
GlobalVector globalMomentum() const
std::vector< const reco::VertexCollection * > vertices
std::vector< double > originHalfLength
const math::XYZTLorentzVectorD & momentum() const
particle info...
std::vector< double > maxZ0
unsigned int subDetId() const
The subdet Id.
std::vector< std::string > seedingAlgo
std::vector< edm::InputTag > primaryVertices
std::vector< std::vector< unsigned int > > firstHitSubDetectors
std::vector< double > maxD0
DetId geographicalId() const
double xyImpactParameter(double x0=0., double y0=0.) const
Transverse impact parameter.
math::XYZTLorentzVector XYZTLorentzVector
const MagneticField * theMagField