28 const std::vector<GlobalPoint>& surfaces)
31 std::vector<GlobalPoint> trajectory;
35 for(std::vector<GlobalPoint>::const_iterator surface_iter = surfaces.begin();
36 surface_iter != surfaces.end(); surface_iter++) {
38 std::unique_ptr<Cylinder> cylinder = std::make_unique<Cylinder>(surface_iter->perp(),
41 std::unique_ptr<Plane> forwardEndcap = std::make_unique<Plane>(
Surface::PositionType(0,0,surface_iter->z()),
43 std::unique_ptr<Plane> backwardEndcap = std::make_unique<Plane>(
Surface::PositionType(0,0,-surface_iter->z()),
47 LogTrace(
"StartingPoint")<<
"Propagate from "<<
"\n" 48 <<
"\tx: " << ftsStart.
position().
x()<<
"\n" 49 <<
"\ty: " << ftsStart.
position().
y()<<
"\n" 50 <<
"\tz: " << ftsStart.
position().
z()<<
"\n" 51 <<
"\tmomentum eta: " << ftsStart.
momentum().
eta()<<
"\n" 52 <<
"\tmomentum phi: " << ftsStart.
momentum().
phi()<<
"\n" 53 <<
"\tmomentum: " << ftsStart.
momentum().
mag()<<
"\n";
56 float corner = surface_iter->perp()/surface_iter->z();
71 if (fabs(tanTheta) > corner)
76 else if(tanTheta > 0.)
87 if (tanTheta < 0 ) tSOSDest =
ivProp_->
propagate( ftsCurrent,*forwardEndcap);
88 if (tanTheta >= 0 ) tSOSDest =
ivProp_->
propagate( ftsCurrent,*backwardEndcap);
106 if (tanTheta < 0 ) tSOSDest =
ivProp_->
propagate( ftsCurrent,*forwardEndcap);
107 if (tanTheta >= 0 ) tSOSDest =
ivProp_->
propagate( ftsCurrent,*backwardEndcap);
124 if (! tSOSDest.
isValid())
return trajectory;
127 LogTrace(
"SuccessfullPropagation") <<
"Great, I reached something." <<
"\n" 139 trajectory.push_back(point);
152 LogTrace(
"MatchPoint") <<
"point (eta,phi): " << direction.
eta() <<
"," << direction.
phi() <<
"\n";
153 int ieta =
iEta(direction);
154 int iphi =
iPhi(direction);
156 LogTrace(
"MatchPoint") <<
"(ieta,iphi): " << ieta <<
"," << iphi <<
"\n";
158 if (ieta>=0 && ieta<nEta_ && iphi>=0 && iphi<
nPhi_){
159 set = (*theMap_)[ieta][iphi];
162 LogTrace(
"MatchPoint") <<
"Add neighbors (ieta,iphi): " << ieta <<
"," << iphi <<
"\n";
164 int maxIEta = ieta+idR;
165 int minIEta = ieta-idR;
167 if(minIEta<0) minIEta = 0;
168 int maxIPhi = iphi+idR;
169 int minIPhi = iphi-idR;
174 LogTrace(
"MatchPoint") <<
"\tieta (min,max): " << minIEta <<
"," << maxIEta<<
"\n";
175 LogTrace(
"MatchPoint") <<
"\tiphi (min,max): " << minIPhi <<
"," << maxIPhi<<
"\n";
176 for (
int i=minIEta;
i<=maxIEta;
i++)
177 for (
int j=minIPhi;j<=maxIPhi;j++) {
178 if(
i==ieta && j==iphi)
continue;
203 if (point.
eta()>1.827 && point.
eta()<=1.830)
return iEta1-1;
204 else if (point.
eta()>1.914 && point.
eta()<=1.930)
return iEta1-1;
205 else if (point.
eta()>2.001 && point.
eta()<=2.043)
return iEta1-1;
206 else if (point.
eta()>2.088 && point.
eta()<=2.172)
return iEta1-1;
207 else if (point.
eta()>2.175 && point.
eta()<=2.262)
return iEta1-1;
208 else if (point.
eta()>2.262 && point.
eta()<=2.332)
return iEta1-2;
209 else if (point.
eta()>2.332 && point.
eta()<=2.349)
return iEta1-1;
210 else if (point.
eta()>2.349 && point.
eta()<=2.436)
return iEta1-2;
211 else if (point.
eta()>2.436 && point.
eta()<=2.500)
return iEta1-3;
212 else if (point.
eta()>2.500 && point.
eta()<=2.523)
return iEta1-2;
213 else if (point.
eta()>2.523 && point.
eta()<=2.610)
return iEta1-3;
214 else if (point.
eta()>2.610 && point.
eta()<=2.650)
return iEta1-4;
215 else if (point.
eta()>2.650 && point.
eta()<=2.697)
return iEta1-3;
216 else if (point.
eta()>2.697 && point.
eta()<=2.784)
return iEta1-4;
217 else if (point.
eta()>2.784 && point.
eta()<=2.868)
return iEta1-5;
218 else if (point.
eta()>2.868 && point.
eta()<=2.871)
return iEta1-4;
219 else if (point.
eta()>2.871 && point.
eta()<=2.958)
return iEta1-5;
220 else if (point.
eta()>2.958)
return iEta1-6;
221 else if (point.
eta()<-1.827 && point.
eta()>=-1.830)
return iEta1+1;
222 else if (point.
eta()<-1.914 && point.
eta()>=-1.930)
return iEta1+1;
223 else if (point.
eta()<-2.001 && point.
eta()>=-2.043)
return iEta1+1;
224 else if (point.
eta()<-2.088 && point.
eta()>=-2.172)
return iEta1+1;
225 else if (point.
eta()<-2.175 && point.
eta()>=-2.262)
return iEta1+1;
226 else if (point.
eta()<-2.262 && point.
eta()>=-2.332)
return iEta1+2;
227 else if (point.
eta()<-2.332 && point.
eta()>=-2.349)
return iEta1+1;
228 else if (point.
eta()<-2.349 && point.
eta()>=-2.436)
return iEta1+2;
229 else if (point.
eta()<-2.436 && point.
eta()>=-2.500)
return iEta1+3;
230 else if (point.
eta()<-2.500 && point.
eta()>=-2.523)
return iEta1+2;
231 else if (point.
eta()<-2.523 && point.
eta()>=-2.610)
return iEta1+3;
232 else if (point.
eta()<-2.610 && point.
eta()>=-2.650)
return iEta1+4;
233 else if (point.
eta()<-2.650 && point.
eta()>=-2.697)
return iEta1+3;
234 else if (point.
eta()<-2.697 && point.
eta()>=-2.784)
return iEta1+4;
235 else if (point.
eta()<-2.784 && point.
eta()>=-2.868)
return iEta1+5;
236 else if (point.
eta()<-2.868 && point.
eta()>=-2.871)
return iEta1+4;
237 else if (point.
eta()<-2.871 && point.
eta()>=-2.958)
return iEta1+5;
238 else if (point.
eta()<-2.349)
return iEta1+6;
245 double pi=4*atan(1.);
255 LogTrace(
"HDetIdAssociator")<<
"building map" <<
"\n";
257 theMap_ =
new std::vector<std::vector<std::set<DetId> > >(
nEta_,std::vector<std::set<DetId> >(
nPhi_));
258 int numberOfDetIdsOutsideEtaRange = 0;
259 int numberOfDetIdsActive = 0;
261 for (std::set<DetId>::const_iterator id_itr = validIds.begin(); id_itr!=validIds.end(); id_itr++) {
265 if(point.
eta()==0)
continue;
267 int ieta =
iEta(point);
268 int iphi =
iPhi(point);
275 if ( etaMin > ieta) etaMin = ieta;
276 if ( etaMax < ieta) etaMax = ieta;
277 if ( phiMin > iphi) phiMin = iphi;
278 if ( phiMax < iphi) phiMax = iphi;
280 if ((ieta>54||ieta<15) && iphi%2==0) phiMax++;
281 if ((ieta>54||ieta<15) && iphi%2==1) phiMin--;
283 if (etaMax<0||phiMax<0||etaMin>=nEta_||phiMin>=
nPhi_) {
284 LogTrace(
"HDetIdAssociator")<<
"Out of range: DetId:" << id_itr->rawId() <<
285 "\n\teta (min,max): " << etaMin <<
"," << etaMax <<
286 "\n\tphi (min,max): " << phiMin <<
"," << phiMax <<
287 "\nTower id: " << id_itr->rawId() <<
"\n";
288 numberOfDetIdsOutsideEtaRange++;
292 if (phiMax-phiMin > phiMin+
nPhi_-phiMax){
296 for(
int ieta = etaMin; ieta <=
etaMax; ieta++)
297 for(
int iphi = phiMin; iphi <=
phiMax; iphi++)
299 numberOfDetIdsActive++;
301 LogTrace(
"HDetIdAssociator") <<
"Number of elements outside the allowed range ( |eta|>"<<
302 nEta_/2*
etaBinSize_ <<
"): " << numberOfDetIdsOutsideEtaRange <<
"\n";
303 LogTrace(
"HDetIdAssociator") <<
"Number of active DetId's mapped: " <<
304 numberOfDetIdsActive <<
"\n";
309 const std::vector<GlobalPoint>& trajectory,
314 std::set<DetId> outset;
317 for(std::set<DetId>::const_iterator id_iter = inset.begin(); id_iter != inset.end(); id_iter++)
318 for(std::vector<GlobalPoint>::const_iterator point_iter = trajectory.begin(); point_iter != trajectory.end(); point_iter++)
319 if (
nearElement(*point_iter,*id_iter,dR)) outset.insert(*id_iter);
322 if (inset.size()!=1)
return outset;
323 std::set<DetId>::const_iterator id_inp = inset.begin();
329 for (
int i=ieta-1;
i<=ieta+1;
i++) {
330 for (
int j=iphi-1;j<=iphi+1;j++) {
332 if( i<0 || i>=
nEta_)
continue;
333 int j2fill = j%
nPhi_;
334 if(j2fill<0) j2fill+=
nPhi_;
354 const std::vector<GlobalPoint>& trajectory)
357 std::set<DetId> outset;
358 for(std::set<DetId>::const_iterator id_iter = inset.begin(); id_iter != inset.end(); id_iter++)
359 for(std::vector<GlobalPoint>::const_iterator point_iter = trajectory.begin(); point_iter != trajectory.end(); point_iter++)
360 if (
insideElement(*point_iter, *id_iter)) outset.insert(*id_iter);
370 std::set<DetId> outset;
371 std::set<DetId>::const_iterator id_max = inset.begin();
374 for(std::set<DetId>::const_iterator id_iter = inset.begin(); id_iter != inset.end(); id_iter++) {
380 if(tower != (*caloTowers).end() && tower->hadEnergy()>Ehadmax) {
382 Ehadmax = tower->hadEnergy();
386 if (Ehadmax > 0) outset.insert(*id_max);
405 std::set<DetId> outset;
406 std::set<DetId>::const_iterator id_max = inset.begin();
409 for(std::set<DetId>::const_iterator id_iter = inset.begin(); id_iter != inset.end(); id_iter++) {
415 if(hit != (*recHits).end() && hit->energy()>Ehadmax) {
417 Ehadmax = hit->energy();
421 if (Ehadmax > 0) outset.insert(*id_max);
virtual void check_setup()
virtual std::vector< GlobalPoint > getTrajectory(const FreeTrajectoryState &, const std::vector< GlobalPoint > &)
virtual std::set< DetId > getDetIdsCloseToAPoint(const GlobalPoint &, const int idR=0)
virtual GlobalPoint getPosition(const DetId &)=0
std::vector< std::vector< std::set< DetId > > > * theMap_
virtual bool insideElement(const GlobalPoint &, const DetId &)=0
Geom::Phi< T > phi() const
std::vector< CaloTower >::const_iterator const_iterator
GlobalPoint globalPosition() const
virtual bool nearElement(const GlobalPoint &point, const DetId &id, const double distance)
void swap(edm::DataFrameContainer &lhs, edm::DataFrameContainer &rhs)
FreeTrajectoryState const * freeState(bool withErrors=true) const
virtual int iPhi(const GlobalPoint &)
Abs< T >::type abs(const T &t)
virtual std::set< DetId > getDetIdsInACone(const std::set< DetId > &, const std::vector< GlobalPoint > &trajectory, const double)
virtual std::set< DetId > getCrossedDetIds(const std::set< DetId > &, const std::vector< GlobalPoint > &trajectory)
virtual std::set< DetId > getASetOfValidDetIds()=0
virtual std::set< DetId > getMaxEDetId(const std::set< DetId > &, edm::Handle< CaloTowerCollection > caloTowers)
GlobalVector momentum() const
virtual int iEta(const GlobalPoint &)
GlobalPoint position() const
Point3DBase< float, GlobalTag > PositionType
TrajectoryStateOnSurface propagate(STA const &state, SUR const &surface) const
TkRotation< float > RotationType
*vegas h *****************************************************used in the default bin number in original ***version of VEGAS is ***a higher bin number might help to derive a more precise ***grade subtle point