32 {
35
37 mcCollptra->
push_back(
new HepMC::GenEvent(*evt));
38 }
39
41#ifdef HEPMC3
42 std::vector<HepMC::GenVertexPtr> lambda_vertices_to_remove;
43 for (auto part: **evt){
44 const int abspid = std::abs(
part->pdg_id());
45 if (abspid != 310 && abspid != 3122 && abspid != 3222 && abspid != 3112 && abspid != 3322 && abspid != 3312 && abspid != 3334 ) continue;
46 if (!
part->end_vertex())
continue;
47 lambda_vertices_to_remove.push_back(
part->end_vertex());
49 }
50 for (auto v: lambda_vertices_to_remove) (*evt)->remove_vertex(std::move(v));
51#else
52
53 std::set<int> Lambdas;
54 std::set<int> vtx_to_delete;
55 for (HepMC::GenEvent::vertex_const_iterator vtx = (*evt)->vertices_begin(); vtx != (*evt)->vertices_end(); ++vtx) {
56
57 if ( vtx_to_delete.find((*vtx)->barcode()) == vtx_to_delete.end() &&
58 (*vtx)->particles_in_const_begin() != (*vtx)->particles_in_const_end() ) {
59 HepMC::GenVertex::particles_in_const_iterator
part = (*vtx)->particles_in_const_begin();
60 bool lambda_not_found = true;
61 const int abspid = std::abs((*part)->pdg_id());
62 do {
63 if (abspid == 310 || abspid == 3122 || abspid == 3222 || abspid == 3112 ||
64 abspid == 3322 || abspid == 3312 || abspid == 3334 ) {
65 lambda_not_found = false;
66 Lambdas.insert((*part)->barcode());
67 vtx_to_delete.insert((*vtx)->barcode());
68
69
70 for (HepMC::GenVertex::vertex_iterator desc = (*vtx)->vertices_begin(HepMC::descendants);
71 desc != (*vtx)->vertices_end(HepMC::descendants); ++desc) {
72 vtx_to_delete.insert((*desc)->barcode());
73 }
74 }
76 } while (part != (*vtx)->particles_in_const_end() && lambda_not_found);
77 }
78 }
79
80
81 for (std::set<int>::iterator l = Lambdas.begin(); l != Lambdas.end(); ++l) {
82 HepMC::GenParticle* lam = (*evt)->barcode_to_particle(*l);
83 lam->set_status(1);
84 }
85
86 for (std::set<int>::iterator v = vtx_to_delete.begin(); v != vtx_to_delete.end(); ++v) {
87 HepMC::GenVertex* vdel = (*evt)->barcode_to_vertex(*v);
88 (*evt)->remove_vertex(vdel);
89 delete vdel;
90 }
91#endif
92 }
93
96 return StatusCode::FAILURE;
97 }
98 }
99
100
102 auto pdfinfo =
evt->pdf_info();
103 auto ion =
evt->heavy_ion();
104#ifdef HEPMC3
105 if (pdfinfo) {
106 std::cout << "PdfInfo: "
107 << pdfinfo->parton_id[0] << ", "
108 << pdfinfo->parton_id[1] << ", "
109 << pdfinfo->x[0] << ", "
110 << pdfinfo->x[1] << ", "
111 << pdfinfo->scale << ", "
112 << pdfinfo->xf[0] << ", "
113 << pdfinfo->xf[1] << ", "
114 << pdfinfo->pdf_id[0] << ", "
115 << pdfinfo->pdf_id[1]
116 << std::endl;
117 }
118
119 if (ion) {
120 std::cout << std::endl;
121 std::cout << "Heavy Ion: "
122 << ion->Ncoll_hard <<", "
123 << ion->Npart_proj <<" , "
124 << ion->Npart_targ<< ", "
125 << ion->Ncoll<< ", "
126 << ion->spectator_neutrons << ", "
127 << ion->spectator_protons << ", "
128 << ion->N_Nwounded_collisions << ", "
129 << ion->Nwounded_N_collisions << ", "
130 << ion->Nwounded_Nwounded_collisions << ", "
131 << ion->impact_parameter << ", "
132 << ion->event_plane_angle << ", "
133 << ion->eccentricity << ", "
134 << ion->sigma_inel_NN
135 << std::endl;
136 }
137#else
138 if (pdfinfo) {
139 std::cout << "PdfInfo: "
140 << pdfinfo->id1() << ", "
141 << pdfinfo->id2() << ", "
142 << pdfinfo->x1() << ", "
143 << pdfinfo->x2() << ", "
144 << pdfinfo->scalePDF() << ", "
145 << pdfinfo->pdf1() << ", "
146 << pdfinfo->pdf2() << ", "
147 << pdfinfo->pdf_id1() << ", "
148 << pdfinfo->pdf_id2()
149 << std::endl;
150 }
151
152 if (ion) {
153 std::cout << std::endl;
154 std::cout << "Heavy Ion: "
155 << ion->Ncoll_hard() <<", "
156 << ion->Npart_proj() <<" , "
157 << ion->Npart_targ()<< ", "
158 << ion->Ncoll()<< ", "
159 << ion->spectator_neutrons() << ", "
160 << ion->spectator_protons() << ", "
161 << ion->N_Nwounded_collisions() << ", "
162 << ion->Nwounded_N_collisions() << ", "
163 << ion->Nwounded_Nwounded_collisions() << ", "
164 << ion->impact_parameter() << ", "
165 << ion->event_plane_angle() << ", "
166 << ion->eccentricity() << ", "
167 << ion->sigma_inel_NN()
168 << std::endl;
169 }
170#endif
174 } else {
175
176 for (auto part: *evt) {
177 double rapid =
part->momentum().pseudoRapidity();
178 double phi =
part->momentum().phi();
179 double et=
part->momentum().perp();
180 int p_stat=
part->status();
181 int p_id =
part->pdg_id();
182 std::cout <<
" eta = " << rapid<<
" Phi = " <<
phi <<
" Et = " <<
et/
GeV <<
" Status= " << p_stat <<
" PDG ID= "<< p_id << std::endl;
183 }
184 }
186 for (const auto& pitr: *evt) {
187 int p_stat=pitr->status();
188 if(p_stat==2 && pitr->production_vertex() && pitr->end_vertex()) {
189 const auto& prodVtx = pitr->production_vertex()->position();
190 const auto& endVtx = pitr->end_vertex()->position();
191 const CLHEP::HepLorentzVector lv0( prodVtx.x(), prodVtx.y(), prodVtx.z(), prodVtx.t() );
192 const CLHEP::HepLorentzVector lv1( endVtx.x(), endVtx.y(), endVtx.z(), endVtx.t() );
193
194 CLHEP::HepLorentzVector dist4D(lv1);
195 dist4D-=lv0;
196 CLHEP::Hep3Vector dist3D=dist4D.vect();
197 if(dist3D.mag()>1*Gaudi::Units::mm) {
198 const auto& GenMom = pitr->momentum();
199 const CLHEP::HepLorentzVector
mom( GenMom.x(), GenMom.y(), GenMom.z(), GenMom.t() );
203 ATH_MSG_INFO(
" gamma(Momentum)="<<
mom.gamma()<<
" gamma(Vertices)="<<dist4D.gamma());
204 }
205 }
206 }
207 }
208 }
209 }
210 return StatusCode::SUCCESS;
211}
Scalar phi() const
phi method
float et(const xAOD::jFexSRJetRoI *j)
ServiceHandle< StoreGateSvc > & evtStore()
value_type push_back(value_type pElem)
Add an element to the end of the collection.
DataModel_detail::iterator< DataVector > iterator
const_iterator end() const noexcept
Return a const_iterator pointing past the end of the collection.
void line(std::ostream &os, const GenEvent &e)
1.14.0