43 ROOT::Math::SMatrixIdentity
id;
45 const float uncert = 0.001;
49 void setPropagationTargetZ(
const float &
z);
51 bool propagate(
const double px,
const double py,
const double pz,
const double x,
const double y,
69 targetPlaneBackward_ =
74 const double x,
const double y,
const double z,
89 trackStateP = RKProp().propagate(startingStateP, *targetPlaneForward_);
91 trackStateP = RKProp().propagate(startingStateP, *targetPlaneBackward_);
107 return propagate(momentum.px(), momentum.py(), momentum.pz(), position.x(), position.y(),
128 outsideEESurface = 1,
140 float gen_TrueNumInt_;
150 std::vector<
float> genpart_eta_;
151 std::vector<
float> genpart_phi_;
152 std::vector<
float> genpart_pt_;
153 std::vector<
float> genpart_energy_;
154 std::vector<
float> genpart_dvx_;
155 std::vector<
float> genpart_dvy_;
156 std::vector<
float> genpart_dvz_;
157 std::vector<
float> genpart_ovx_;
158 std::vector<
float> genpart_ovy_;
159 std::vector<
float> genpart_ovz_;
160 std::vector<
float> genpart_exx_;
161 std::vector<
float> genpart_exy_;
162 std::vector<
int> genpart_mother_;
163 std::vector<
float> genpart_exphi_;
164 std::vector<
float> genpart_exeta_;
165 std::vector<
float> genpart_fbrem_;
166 std::vector<
int> genpart_pid_;
167 std::vector<
int> genpart_gen_;
168 std::vector<
int> genpart_reachedEE_;
169 std::vector<
bool> genpart_fromBeamPipe_;
170 std::vector<
std::vector<
float>> genpart_posx_;
171 std::vector<
std::vector<
float>> genpart_posy_;
172 std::vector<
std::vector<
float>> genpart_posz_;
177 std::vector<
float> gen_eta_;
178 std::vector<
float> gen_phi_;
179 std::vector<
float> gen_pt_;
180 std::vector<
float> gen_energy_;
181 std::vector<
int> gen_charge_;
182 std::vector<
int> gen_pdgid_;
183 std::vector<
int> gen_status_;
184 std::vector<
std::vector<
int>> gen_daughters_;
208 "HGCalTriggerNtupleGen" );
211 HGCalTriggerNtupleGen::
223 mySimEvent_ = std::make_unique<FSimEvent>(particleFilter_);
227 tree.Branch(
"gen_n", &gen_n_,
"gen_n/I");
228 tree.Branch(
"gen_PUNumInt", &gen_PUNumInt_ ,
"gen_PUNumInt/I");
229 tree.Branch(
"gen_TrueNumInt", &gen_TrueNumInt_ ,
"gen_TrueNumInt/F");
236 tree.Branch(
"vtx_x", &vtx_x_);
237 tree.Branch(
"vtx_y", &vtx_y_);
238 tree.Branch(
"vtx_z", &vtx_z_);
241 tree.Branch(
"gen_eta", &gen_eta_);
242 tree.Branch(
"gen_phi", &gen_phi_);
243 tree.Branch(
"gen_pt", &gen_pt_);
244 tree.Branch(
"gen_energy", &gen_energy_);
245 tree.Branch(
"gen_charge", &gen_charge_);
246 tree.Branch(
"gen_pdgid", &gen_pdgid_);
247 tree.Branch(
"gen_status", &gen_status_);
248 tree.Branch(
"gen_daughters", &gen_daughters_);
250 tree.Branch(
"genpart_eta", &genpart_eta_);
251 tree.Branch(
"genpart_phi", &genpart_phi_);
252 tree.Branch(
"genpart_pt", &genpart_pt_);
253 tree.Branch(
"genpart_energy", &genpart_energy_);
254 tree.Branch(
"genpart_dvx", &genpart_dvx_);
255 tree.Branch(
"genpart_dvy", &genpart_dvy_);
256 tree.Branch(
"genpart_dvz", &genpart_dvz_);
257 tree.Branch(
"genpart_ovx", &genpart_ovx_);
258 tree.Branch(
"genpart_ovy", &genpart_ovy_);
259 tree.Branch(
"genpart_ovz", &genpart_ovz_);
260 tree.Branch(
"genpart_mother", &genpart_mother_);
261 tree.Branch(
"genpart_exphi", &genpart_exphi_);
262 tree.Branch(
"genpart_exeta", &genpart_exeta_);
263 tree.Branch(
"genpart_exx", &genpart_exx_);
264 tree.Branch(
"genpart_exy", &genpart_exy_);
265 tree.Branch(
"genpart_fbrem", &genpart_fbrem_);
266 tree.Branch(
"genpart_pid", &genpart_pid_);
267 tree.Branch(
"genpart_gen", &genpart_gen_);
268 tree.Branch(
"genpart_reachedEE", &genpart_reachedEE_);
269 tree.Branch(
"genpart_fromBeamPipe", &genpart_fromBeamPipe_);
270 tree.Branch(
"genpart_posx", &genpart_posx_);
271 tree.Branch(
"genpart_posy", &genpart_posy_);
272 tree.Branch(
"genpart_posz", &genpart_posz_);
285 const std::vector< PileupSummaryInfo >& PupInfo = *PupInfo_h;
287 if(pdt_watcher_.check(es))
291 mySimEvent_->initializePdt(&(*pdt));
294 if(magfield_watcher_.check(es))
298 aField_ = &(*magfield);
301 triggerTools_.eventSetup(es);
310 iEvent.
getByToken(simTracks_token_, simTracksHandle);
311 iEvent.
getByToken(simVertices_token_, simVerticesHandle);
312 mySimEvent_->fill(*simTracksHandle, *simVerticesHandle);
314 HepMC::GenVertex *
primaryVertex = *(hevH)->GetEvent()->vertices_begin();
315 const float mm2cm = 0.1;
316 vtx_x_ = primaryVertex->position().x() * mm2cm;
317 vtx_y_ = primaryVertex->position().y() * mm2cm;
318 vtx_z_ = primaryVertex->position().z() * mm2cm;
319 Point sim_pv(vtx_x_, vtx_y_, vtx_z_);
324 std::vector<FSimTrack *> allselectedgentracks;
325 const float eeInnerRadius = 25.;
326 const float eeOuterRadius = 160.;
327 unsigned int npart = mySimEvent_->nTracks();
328 for (
unsigned int i = 0;
i <
npart; ++
i) {
329 std::vector<float> xp, yp, zp;
333 int reachedEE = ReachHGCal::notReach;
337 if (
std::abs(myTrack.vertex().position().z()) >= triggerTools_.getLayerZ(1))
continue;
339 const unsigned nlayers = triggerTools_.lastLayerBH();
340 if (myTrack.noEndVertex())
343 bool reachesHGCal = toHGCalPropagator.
propagate(
344 myTrack.momentum(), myTrack.vertex().position(), myTrack.charge(), propcoords);
347 if (reachesHGCal && vtx.Rho() < eeOuterRadius && vtx.Rho() > eeInnerRadius) {
348 reachedEE = ReachHGCal::onEESurface;
351 for (
int i = 0;
i < myTrack.nDaughters(); ++
i) dpt += myTrack.daughter(
i).momentum().pt();
352 if (
abs(myTrack.type()) == 11) fbrem = dpt / myTrack.momentum().pt();
353 }
else if (reachesHGCal && vtx.Rho() > eeOuterRadius)
354 reachedEE = ReachHGCal::outsideEESurface;
357 for (
unsigned il = 1; il <=
nlayers; ++il) {
358 const float charge = myTrack.charge();
361 indiv_particleProp.
propagate(myTrack.momentum(), myTrack.vertex().position(),
charge,
364 xp.push_back(propCoords.
x);
365 yp.push_back(propCoords.
y);
366 zp.push_back(propCoords.
z);
369 vtx = myTrack.endVertex().position();
371 auto orig_vtx = myTrack.vertex().position();
373 allselectedgentracks.push_back(&mySimEvent_->track(
i));
375 genpart_eta_.push_back(myTrack.momentum().eta());
376 genpart_phi_.push_back(myTrack.momentum().phi());
377 genpart_pt_.push_back(myTrack.momentum().pt());
378 genpart_energy_.push_back(myTrack.momentum().energy());
379 genpart_dvx_.push_back(vtx.x());
380 genpart_dvy_.push_back(vtx.y());
381 genpart_dvz_.push_back(vtx.z());
383 genpart_ovx_.push_back(orig_vtx.x());
384 genpart_ovy_.push_back(orig_vtx.y());
385 genpart_ovz_.push_back(orig_vtx.z());
388 toHGCalPropagator.
propagate(myTrack.momentum(), orig_vtx, myTrack.charge(), hitsHGCal);
390 genpart_exphi_.push_back(hitsHGCal.
phi);
391 genpart_exeta_.push_back(hitsHGCal.
eta);
392 genpart_exx_.push_back(hitsHGCal.
x);
393 genpart_exy_.push_back(hitsHGCal.
y);
395 genpart_fbrem_.push_back(fbrem);
396 genpart_pid_.push_back(myTrack.type());
397 genpart_gen_.push_back(myTrack.genpartIndex());
398 genpart_reachedEE_.push_back(reachedEE);
399 genpart_fromBeamPipe_.push_back(
true);
401 genpart_posx_.push_back(xp);
402 genpart_posy_.push_back(yp);
403 genpart_posz_.push_back(zp);
408 iEvent.
getByToken(gen_token_, genParticlesHandle);
409 gen_n_ = genParticlesHandle->size();
411 for (
const auto& particle : *genParticlesHandle) {
412 gen_eta_.push_back(particle.eta());
413 gen_phi_.push_back(particle.phi());
414 gen_pt_.push_back(particle.pt());
415 gen_energy_.push_back(particle.energy());
416 gen_charge_.push_back(particle.charge());
417 gen_pdgid_.push_back(particle.pdgId());
418 gen_status_.push_back(particle.status());
419 std::vector<int>
daughters(particle.daughterRefVector().size(), 0);
420 for (
unsigned j = 0; j < particle.daughterRefVector().size(); ++j) {
421 daughters[j] =
static_cast<int>(particle.daughterRefVector().at(j).key());
428 genpart_mother_.resize(genpart_posz_.size(), -1);
429 for (
size_t i = 0;
i < allselectedgentracks.size();
i++) {
430 const auto tracki = allselectedgentracks.at(
i);
432 for (
size_t j =
i + 1; j < allselectedgentracks.size(); j++) {
433 const auto trackj = allselectedgentracks.at(j);
435 if (!tracki->noMother()) {
436 if (&tracki->mother() == trackj) genpart_mother_.at(
i) = j;
438 if (!trackj->noMother()) {
439 if (&trackj->mother() == tracki) genpart_mother_.at(j) =
i;
444 for(
const auto& PVI : PupInfo)
446 if(PVI.getBunchCrossing() == 0)
448 gen_PUNumInt_ = PVI.getPU_NumInteractions();
449 gen_TrueNumInt_ = PVI.getTrueNumInteractions();
463 gen_TrueNumInt_ = 0.;
471 genpart_eta_.clear();
472 genpart_phi_.clear();
474 genpart_energy_.clear();
475 genpart_dvx_.clear();
476 genpart_dvy_.clear();
477 genpart_dvz_.clear();
478 genpart_ovx_.clear();
479 genpart_ovy_.clear();
480 genpart_ovz_.clear();
481 genpart_exx_.clear();
482 genpart_exy_.clear();
483 genpart_mother_.clear();
484 genpart_exphi_.clear();
485 genpart_exeta_.clear();
486 genpart_fbrem_.clear();
487 genpart_pid_.clear();
488 genpart_gen_.clear();
489 genpart_reachedEE_.clear();
490 genpart_fromBeamPipe_.clear();
491 genpart_posx_.clear();
492 genpart_posy_.clear();
493 genpart_posz_.clear();
505 gen_daughters_.clear();
std::vector< GenParticle > GenParticleCollection
collection of GenParticles
T getParameter(std::string const &) const
static AlgebraicMatrix initialize()
const MagneticField * field_
ROOT::Math::LorentzVector< ROOT::Math::PxPyPzE4D< double > > XYZTLorentzVectorD
Lorentz vector with cylindrical internal representation using pseudorapidity.
bool getByToken(EDGetToken token, Handle< PROD > &result) const
Geom::Phi< T > phi() const
CurvilinearTrajectoryError err_
GlobalPoint globalPosition() const
ROOT::Math::SMatrix< double, 5, 5, ROOT::Math::MatRepSym< double, 5 > > AlgebraicSymMatrix55
void initialize(TTree &, const edm::ParameterSet &, edm::ConsumesCollector &&) final
void clear(CLHEP::HepGenMatrix &m)
Helper function: Reset all elements of a matrix to 0.
static PlanePointer build(Args &&...args)
void fill(const edm::Event &, const edm::EventSetup &) final
Abs< T >::type abs(const T &t)
const RKPropagatorInS & RKProp() const
ROOT::Math::Transform3DPJ Transform3D
Plane::PlanePointer targetPlaneForward_
TrajectoryStateOnSurface TSOS
SimpleTrackPropagator(const MagneticField *f)
bool propagate(const double px, const double py, const double pz, const double x, const double y, const double z, const float charge, Coordinates &coords) const
void setPropagationTargetZ(const float &z)
math::XYZTLorentzVectorD toVector()
static int position[264][3]
#define DEFINE_EDM_PLUGIN(factory, type, name)
ROOT::Math::Transform3DPJ::Point Point
defaultRKPropagator::Product prod_