151 using namespace reco;
161 int ZprodVtx, ZdecayVtx = 0;
162 bool vtxFound =
false;
163 HepMC::GenEvent::particle_const_iterator parIt, parItE;
164 parIt = prodIn->GetEvent()->particles_begin();
165 parItE = prodIn->GetEvent()->particles_end();
166 for ( ;parIt != parItE; ++parIt){
167 int pdg =
abs( (*parIt)->pdg_id() ) ;
169 ZdecayVtx = (*parIt)->end_vertex()->barcode();
170 ZprodVtx = (*parIt)->production_vertex()->barcode();
184 std::set<int> barcodesZ;
185 barcodesZ.insert(ZdecayVtx);
187 HepMC::GenEvent::vertex_const_iterator itVtx, itVtxE;
188 itVtx = prodIn->GetEvent()->vertices_begin();
189 itVtxE = prodIn->GetEvent()->vertices_end();
191 for(;itVtx!=itVtxE;++itVtx){
192 if ( barcodesZ.find( (*itVtx)->barcode() )!=barcodesZ.end()){
195 HepMC::GenVertex::particles_out_const_iterator parIt, parItE;
196 parIt=(*itVtx)->particles_out_const_begin();
197 parItE=(*itVtx)->particles_out_const_end();
198 for ( ; parIt!=parItE; ++parIt){
199 if ((*parIt)->end_vertex()){
200 barcodesZ.insert((*parIt)->end_vertex()->barcode());
208 barcodesZ.insert(ZprodVtx);
209 barcodesZ.insert(ZdecayVtx);
216 std::auto_ptr<HepMC::GenEvent> evtUE(
new HepMC::GenEvent());
217 std::auto_ptr<HepMC::GenEvent> evtTauTau(
new HepMC::GenEvent());
220 itVtx = prodIn->GetEvent()->vertices_begin();
221 itVtxE = prodIn->GetEvent()->vertices_end();
224 std::map<int, HepMC::GenVertex*> barcodeMapUE;
225 std::map<int, HepMC::GenVertex*> barcodeMapTauTau;
226 for(;itVtx!=itVtxE;++itVtx){
227 bool isZ = barcodesZ.find( (*itVtx)->barcode() )!=barcodesZ.end();
229 bool leadToZ =
false;
230 if ( (*parIt)->production_vertex() ) {
236 if (isZ || (*itVtx)->barcode()==ZprodVtx || leadToZ ){
237 HepMC::GenVertex* newvertex =
new HepMC::GenVertex((*itVtx)->position(),
239 (*itVtx)->weights() );
241 newvertex->suggest_barcode( -(++newBarcode) );
243 barcodeMapTauTau[(*itVtx)->barcode()] = newvertex;
244 evtTauTau->add_vertex(newvertex);
247 if (!isZ || (*itVtx)->barcode()==ZprodVtx){
248 HepMC::GenVertex* newvertex =
new HepMC::GenVertex((*itVtx)->position(),
250 (*itVtx)->weights() );
252 newvertex->suggest_barcode( -(++newBarcode) );
253 barcodeMapUE[ (*itVtx)->barcode()] = newvertex;
254 evtUE->add_vertex(newvertex);
265 parIt = prodIn->GetEvent()->particles_begin();
266 parItE = prodIn->GetEvent()->particles_end();
267 for ( ;parIt != parItE; ++parIt){
270 bool fromZprodVtx =
false;
271 if ( (*parIt)->production_vertex() ){
272 fromZprodVtx = ( (*parIt)->production_vertex()->barcode() == ZprodVtx);
277 if ( (*parIt)->production_vertex() ){
278 fromZ = barcodesZ.find( (*parIt)->production_vertex()->barcode() )!=barcodesZ.end();
283 bool leadToZ =
false;
285 if ( (*parIt)->end_vertex()) {
287 leadToZ =
std::abs((*parIt)->end_vertex()->barcode()) <
std::abs(ZdecayVtx);
294 newparticle->suggest_barcode(newBarcode++);
295 if ( (*parIt)->end_vertex()){
297 barcodeMapUE[ (*parIt)->end_vertex()->barcode()]->add_particle_in(newparticle);
303 newparticle->suggest_barcode(newBarcode++);
304 if ( (*parIt)->end_vertex()){
307 barcodeMapTauTau[ (*parIt)->end_vertex()->barcode()]->add_particle_in(newparticle);
312 }
else if (fromZprodVtx) {
313 bool isZ = (*parIt)->pdg_id();
315 newparticle->suggest_barcode(newBarcode++);
318 if ( (*parIt)->end_vertex())
319 barcodeMapTauTau[(*parIt)->end_vertex()->barcode()]->add_particle_in(newparticle);
320 if ((*parIt)->production_vertex())
321 barcodeMapTauTau[(*parIt)->production_vertex()->barcode()]->add_particle_out(newparticle);
323 if ((*parIt)->end_vertex())
324 barcodeMapUE[(*parIt)->end_vertex()->barcode()]->add_particle_in(newparticle);
325 if ((*parIt)->production_vertex())
326 barcodeMapUE[(*parIt)->production_vertex()->barcode()]->add_particle_out(newparticle);
331 newparticle->suggest_barcode(newBarcode++);
334 if ((*parIt)->end_vertex())
335 barcodeMapTauTau[(*parIt)->end_vertex()->barcode()]->add_particle_in(newparticle);
336 if ((*parIt)->production_vertex())
337 barcodeMapTauTau[(*parIt)->production_vertex()->barcode()]->add_particle_out(newparticle);
339 if ((*parIt)->end_vertex())
341 barcodeMapUE[(*parIt)->end_vertex()->barcode()]->add_particle_in(newparticle);
342 if ((*parIt)->production_vertex())
343 barcodeMapUE[(*parIt)->production_vertex()->barcode()]->add_particle_out(newparticle);
353 if (orgCnt != tautauCnt + UEcnt){
354 std::cout <<
" XXX Warning - wrong num particles - ORG: "
356 <<
" tautau " << tautauCnt
374 prodUE->addHepMCData( evtUE.release() );
376 else iEvent.
put( prodUE);
383 prodTauTau->addHepMCData( evtTauTau.release() );
384 if (
_doUe ) iEvent.
put( prodTauTau,
"Ztautau");
385 else iEvent.
put( prodTauTau);
int cntStableParticles(const HepMC::GenEvent *evtUE)
OrphanHandle< PROD > put(std::auto_ptr< PROD > product)
Put a new product.
Abs< T >::type abs(const T &t)
bool getByLabel(InputTag const &tag, Handle< PROD > &result) const