analysis/EICAnalysis/TruthTrackerHepMC.C

0001 #include "TruthTrackerHepMC.h"
0002 #include "DISKinematicsReco.h"
0003
0004 /* ROOT includes */
0005 #include <TDatabasePDG.h>
0006
0007 /* STL includes */
0008 #include <iostream>
0009
0010 using namespace std;
0011
0012
0013 TruthTrackerHepMC::TruthTrackerHepMC()
0014 {
0015
0016 }
0017
0018
0019 HepMC::GenParticle* TruthTrackerHepMC::FindParticle( int pid )
0020 {
0021
0022   HepMC::GenParticle *particle = NULL;
0023
0024   /* --> Loop over all truth events in event generator collection */
0025   for (PHHepMCGenEventMap::ReverseIter iter = _genevtmap->rbegin(); iter != _genevtmap->rend(); ++iter)
0026     {
0027       PHHepMCGenEvent *genevt = iter->second;
0028       HepMC::GenEvent *theEvent = genevt->getEvent();
0029
0030       /* check if GenEvent object found */
0031       if ( !theEvent )
0032         {
0033           cout << "ERROR: Missing GenEvent!" << endl;
0034           return NULL;
0035         }
0036
0037       /* Loop over all truth particles in event generator collection */
0038       for ( HepMC::GenEvent::particle_iterator p = theEvent->particles_begin();
0039             p != theEvent->particles_end(); ++p ) {
0040
0041         /* check particle ID */
0042         if ( (*p)->pdg_id() == pid )
0043           {
0044             particle = (*p);
0045
0046             /* end loop if matching particle found */
0047             break;
0048           } // end if matching PID //
0049       } // end loop over all particles in event //
0050     }// end loop over genevtmap //
0051
0052   return particle;
0053 }
0054
0055
0056 HepMC::GenParticle* TruthTrackerHepMC::FindDaughterParticle( int pid, HepMC::GenParticle* particle_mother )
0057 {
0058   HepMC::GenParticle* particle_daughter = NULL;
0059
0060   /*Check if particle mother is found*/
0061   if(!particle_mother);
0062   /* Where does mother particle end (=decay)? */
0063   else if ( particle_mother->end_vertex() )
0064     {
0065       /* Loop over child particles at mother's end vertex */
0066       for ( HepMC::GenVertex::particle_iterator child
0067               = particle_mother->end_vertex()->
0068               particles_begin(HepMC::children);
0069             child != particle_mother->end_vertex()->
0070               particles_end(HepMC::children);
0071             ++child )
0072         {
0073
0074           /* Has child correct PDG code? */
0075           if ( (*child)->pdg_id() == pid )
0076             {
0077               particle_daughter = (*child);
0078               UpdateFinalStateParticle( particle_daughter );
0079         }
0080         }
0081     }
0082
0083   return particle_daughter;
0084 }
0085
0086
0087 HepMC::GenParticle* TruthTrackerHepMC::FindBeamLepton( )
0088 {
0089   HepMC::GenParticle *particle = NULL;
0090
0091   int embedding_id = 1;
0092   PHHepMCGenEvent *theEvent = _genevtmap->get(embedding_id);
0093
0094   /* check if GenEvent object found */
0095   if ( !theEvent )
0096     {
0097       cout << "ERROR: Missing GenEvent!" << endl;
0098       return NULL;
0099     }
0100
0101   HepMC::GenEvent* genEvent = theEvent->getEvent();
0102   /* Find beam lepton */
0103   particle = genEvent->beam_particles().first;
0104   return particle;
0105 }
0106
0107
0108 HepMC::GenParticle* TruthTrackerHepMC::FindBeamHadron( )
0109 {
0110   HepMC::GenParticle *particle = NULL;
0111
0112   int embedding_id = 1;
0113   PHHepMCGenEvent *theEvent = _genevtmap->get(embedding_id);
0114
0115   /* check if GenEvent object found */
0116   if ( !theEvent )
0117     {
0118       cout << "ERROR: Missing GenEvent!" << endl;
0119       return NULL;
0120     }
0121
0122   HepMC::GenEvent* genEvent = theEvent->getEvent();
0123
0124   particle = genEvent->beam_particles().second;
0125
0126   return particle;
0127 }
0128
0129
0130 HepMC::GenParticle* TruthTrackerHepMC::FindScatteredLepton( )
0131 {
0132   HepMC::GenParticle *particle = NULL;
0133
0134   /* @TODO How to select scattered lepton in an unambiguous way for
0135      DIS and Exclusive Processes? (Pythia, Sartre, ...)
0136      Return NULL pointer for now. */
0137   /* In some event generators (SARTRE), no beam particles are created,
0138      if so, assume beam lepton is an electron for now */
0139   int pid_beam_lepton=0;
0140   if(FindBeamLepton()==NULL)
0141     {
0142       /* Assume electron beam */
0143       pid_beam_lepton=11;
0144     }
0145   else
0146     {
0147       pid_beam_lepton = FindBeamLepton()->pdg_id();
0148     }
0149   int embedding_id = 1;
0150
0151   PHHepMCGenEvent *genevt = _genevtmap->get(embedding_id);
0152   HepMC::GenEvent *theEvent = genevt->getEvent();
0153
0154   /* check if GenEvent object found */
0155   if ( !theEvent )
0156     {
0157       cout << "ERROR: Missing GenEvent!" << endl;
0158       return NULL;
0159     }
0160   /* Loop over all truth particles in event generator collection */
0161   for ( HepMC::GenEvent::particle_iterator p = theEvent->particles_begin();
0162     p != theEvent->particles_end(); ++p ) {
0163
0164     /* check particle status and ID */
0165     if ( (*p)->status() == 1 &&
0166      (*p)->pdg_id()  == pid_beam_lepton )
0167       {
0168     particle = (*p);
0169
0170     /* end loop if matching particle found */
0171     break;
0172       } // end if matching status and PID //
0173   } // end loop over all particles in event //
0174   return particle;
0175 }
0176
0177
0178 void
0179 TruthTrackerHepMC::UpdateFinalStateParticle( HepMC::GenParticle *&particle )
0180 {
0181   bool final_state = false;
0182   while ( !final_state )
0183     {
0184       /* if particle has no vertex it is final state */
0185       if(!particle->end_vertex()) break;
0186
0187       /* Loop over all children at the end vertex of this particle */
0188       for ( HepMC::GenVertex::particle_iterator child
0189               = particle->end_vertex()->particles_begin(HepMC::children);
0190             child != particle->end_vertex()->particles_end(HepMC::children);
0191             ++child )
0192         {
0193
0194       /* If there is a child of same particle ID, this was not the final state particle- update pointer to particle and repeat */
0195           if ( (*child)->pdg_id() == particle->pdg_id() )
0196             {
0197               particle = (*child);
0198               break;
0199             }
0200           final_state = true;
0201         }
0202     }
0203   return;
0204 }
0205
0206
0207 void
0208 TruthTrackerHepMC::FindDecayParticles( HepMC::GenParticle *particle_mother, uint &decay_prong, uint &decay_hcharged, uint &decay_lcharged)
0209 {
0210   /* Reset counter variables */
0211   decay_prong = 0;
0212   decay_hcharged = 0;
0213   decay_lcharged = 0;
0214
0215   /* Loop over particles at end vertex of particle */
0216   for ( HepMC::GenVertex::particle_iterator decay
0217           = particle_mother->end_vertex()->
0218           particles_begin(HepMC::children);
0219         decay != particle_mother->end_vertex()->
0220           particles_end(HepMC::children);
0221         ++decay )
0222     {
0223       /* check if particle decays further */
0224       if(!(*decay)->end_vertex()){
0225
0226         /* Get entry from TParticlePDG because HepMC::GenPArticle does not provide charge or class of particle */
0227         TParticlePDG * pdg_p = TDatabasePDG::Instance()->GetParticle( (*decay)->pdg_id() );
0228
0229         /* Check if particle is charged */
0230         if ( pdg_p->Charge() != 0 )
0231           {
0232             decay_prong += 1;
0233
0234             /* Check if particle is lepton */
0235             if ( string( pdg_p->ParticleClass() ) == "Lepton" )
0236               decay_lcharged += 1;
0237
0238             /* Check if particle is hadron, i.e. Meson or Baryon */
0239             else if ( ( string( pdg_p->ParticleClass() ) == "Meson"  ) ||
0240                       ( string( pdg_p->ParticleClass() ) == "Baryon" ) )
0241               decay_hcharged += 1;
0242           }
0243       }
0244
0245       /* loop over decay if particle decays further */
0246       else if((*decay)->end_vertex()){
0247
0248         /*further decay loop */
0249         for ( HepMC::GenVertex::particle_iterator second_decay
0250                 = (*decay)->end_vertex()->
0251                 particles_begin(HepMC::children);
0252               second_decay != (*decay)->end_vertex()->
0253                 particles_end(HepMC::children);
0254               ++second_decay )
0255           {
0256
0257             // Get entry from TParticlePDG because HepMC::GenPArticle does not provide charge or class of particle
0258             TParticlePDG * pdg_p2 = TDatabasePDG::Instance()->GetParticle( (*second_decay)->pdg_id() );
0259
0260             // Check if particle is charged
0261             if ( pdg_p2->Charge() != 0 )
0262               {
0263                 decay_prong += 1;
0264
0265                 // Check if particle is lepton
0266                 if ( string( pdg_p2->ParticleClass() ) == "Lepton" )
0267                   decay_lcharged += 1;
0268
0269                 // Check if particle is hadron, i.e. Meson or Baryon
0270                 else if ( ( string( pdg_p2->ParticleClass() ) == "Meson"  ) ||
0271                           ( string( pdg_p2->ParticleClass() ) == "Baryon" ) )
0272                   decay_hcharged += 1;
0273               }
0274           }
0275       }
0276     }
0277
0278   return;
0279
0280   /* @TODO: Below = Function that calls itself instead of doing the same loop. Still need to debug */
0281
0282   //  for ( HepMC::GenVertex::particle_iterator decay
0283   //          = particle_mother->end_vertex()->
0284   //          particles_begin(HepMC::children);
0285   //        decay != particle_mother->end_vertex()->
0286   //          particles_end(HepMC::children);
0287   //        ++decay )
0288   //    {
0289   //
0290   //      if(!(*decay)->end_vertex()){
0291   //
0292   //        // Get entry from TParticlePDG because HepMC::GenPArticle does not provide charge or class of particle
0293   //        TParticlePDG * pdg_p = TDatabasePDG::Instance()->GetParticle( (*decay)->pdg_id() );
0294   //
0295   //        // Check if particle is charged
0296   //        if ( pdg_p->Charge() != 0 )
0297   //          {
0298   //            decay_prong += 1;
0299   //
0300   //            // Check if particle is lepton
0301   //            if ( string( pdg_p->ParticleClass() ) == "Lepton" )
0302   //              decay_lcharged += 1;
0303   //
0304   //            // Check if particle is hadron, i.e. Meson or Baryon
0305   //            else if ( ( string( pdg_p->ParticleClass() ) == "Meson"  ) ||
0306   //                      ( string( pdg_p->ParticleClass() ) == "Baryon" ) )
0307   //              decay_hcharged += 1;
0308   //          }
0309   //      }
0310   //
0311   //      //HepMC::GenVertex::particle_iterator secondary_decay = (*decay)->end_vertex;
0312   //      else if((*decay)->end_vertex()) FindDecayProducts((*decay), decay_prong, decay_hcharged, decay_lcharged);
0313   //    }
0314
0315 }
0316
0317
0318 void
0319 TruthTrackerHepMC::FindMissingPt( float &pt_miss, float &pt_miss_phi )
0320 {
0321   /* set output values to '0' */
0322   pt_miss = 0;
0323   pt_miss_phi = 0;
0324
0325   /* variables to keep track of components */
0326   float px_sum = 0;
0327   float py_sum = 0;
0328
0329   /* --> Loop over all truth events in event generator collection */
0330   for (PHHepMCGenEventMap::ReverseIter iter = _genevtmap->rbegin(); iter != _genevtmap->rend(); ++iter)
0331     {
0332       PHHepMCGenEvent *genevt = iter->second;
0333       HepMC::GenEvent *theEvent = genevt->getEvent();
0334
0335       /* check if GenEvent object found */
0336       if ( !theEvent )
0337         {
0338           cout << "ERROR: Missing GenEvent!" << endl;
0339           return;
0340         }
0341
0342       /* Loop over all truth particles in event generator collection */
0343       for ( HepMC::GenEvent::particle_iterator p = theEvent->particles_begin();
0344             p != theEvent->particles_end(); ++p ) {
0345
0346         /* skip particles that are not stable final state particles */
0347         if ( (*p)->status() != 1 )
0348           continue;
0349
0350         /* Get entry from TParticlePDG because HepMC::GenPArticle does not provide charge or class of particle */
0351         TParticlePDG * pdg_p = TDatabasePDG::Instance()->GetParticle( (*p)->pdg_id() );
0352
0353         /* skip neutral leptons = neutrinos (invisible to detector) */
0354         if ( ( string( pdg_p->ParticleClass() ) == "Lepton" ) &&
0355              ( pdg_p->Charge() == 0 ) )
0356           continue;
0357
0358         /* update momentum component sum */
0359         px_sum += (*p)->momentum().px();
0360         py_sum += (*p)->momentum().py();
0361
0362       } // end loop over all particles in event //
0363     }// end loop over genevtmap //
0364
0365   /* calculate pt_miss and phi */
0366   pt_miss = sqrt( px_sum * px_sum + py_sum * py_sum );
0367   pt_miss_phi = atan( py_sum / px_sum );
0368
0369   return;
0370 }