analysis/EICAnalysis/Leptoquarks.C

0001 #include "Leptoquarks.h"
0002 #include "TruthTrackerHepMC.h"
0003
0004 #include <phool/getClass.h>
0005 #include <fun4all/Fun4AllServer.h>
0006 #include <fun4all/Fun4AllReturnCodes.h>
0007
0008 #include <phool/PHCompositeNode.h>
0009
0010 #include <TLorentzVector.h>
0011 #include <TNtuple.h>
0012 #include <TFile.h>
0013 #include <TString.h>
0014 #include <TH2D.h>
0015 #include <TDatabasePDG.h>
0016
0017 #include <phhepmc/PHHepMCGenEvent.h>
0018 #include <phhepmc/PHHepMCGenEventMap.h>
0019 #include <HepMC/GenEvent.h>
0020 #include <HepMC/GenVertex.h>
0021
0022 using namespace std;
0023
0024 Leptoquarks::Leptoquarks(std::string filename) :
0025   SubsysReco("Leptoquarks" ),
0026   _foutname(filename),
0027   _fout_root(NULL),
0028   _tree_event(NULL),
0029   _ebeam_E(10),
0030   _pbeam_E(250),
0031   _embedding_id(1)
0032 {
0033
0034 }
0035
0036 int
0037 Leptoquarks::Init(PHCompositeNode *topNode)
0038 {
0039   _ievent = 0;
0040
0041   _fout_root = new TFile(_foutname.c_str(), "RECREATE");
0042
0043   _tree_event = new TNtuple("ntp_event","event information from LQ events",
0044                             "ievent:pt_miss:pt_miss_phi:has_tau:has_uds:tau_eta:tau_phi:tau_e:tau_p:tau_pt:tau_decay_prong:tau_decay_hcharged:tau_decay_lcharged:uds_eta:uds_phi:uds_e:uds_p:uds_pt");
0045
0046
0047   return 0;
0048 }
0049
0050 int
0051 Leptoquarks::process_event(PHCompositeNode *topNode)
0052 {
0053   /* Access GenEventMap */
0054   PHHepMCGenEventMap * geneventmap = findNode::getClass<PHHepMCGenEventMap>(topNode, "PHHepMCGenEventMap");
0055   if (!geneventmap)
0056     {
0057       std::cout <<PHWHERE<<" - Fatal error - missing node PHHepMCGenEventMap"<<std::endl;
0058       return Fun4AllReturnCodes::ABORTRUN;
0059     }
0060
0061   /* Look for leptoquark and tau particle in the event */
0062   TruthTrackerHepMC truth;
0063   truth.set_hepmc_geneventmap( geneventmap );
0064
0065   int pdg_lq = 39; // leptoquark
0066   int pdg_tau = 15; // tau lepton
0067
0068   /* Search for leptoquark in event */
0069   HepMC::GenParticle* particle_lq = truth.FindParticle( pdg_lq );
0070
0071   /* Search for lq->tau decay in event */
0072   HepMC::GenParticle* particle_tau = truth.FindDaughterParticle( pdg_tau, particle_lq );
0073
0074   /* Search for lq->quark decay in event.
0075    * Loop over all quark PDG codes until finding a matching quark. */
0076   HepMC::GenParticle* particle_quark = NULL;
0077   for ( int pdg_quark = 1; pdg_quark < 7; pdg_quark++ )
0078     {
0079       /* try quark */
0080       particle_quark = truth.FindDaughterParticle( pdg_quark, particle_lq );
0081       if (particle_quark)
0082         break;
0083
0084       /* try anti-quark */
0085       particle_quark = truth.FindDaughterParticle( -pdg_quark, particle_lq );
0086       if (particle_quark)
0087         break;
0088     }
0089
0090   /* extract numbers from truth particles */
0091   bool tau_found = false;
0092   bool quark_found = false;
0093
0094   double tau_eta = 0;
0095   double tau_phi = 0;
0096   double tau_e = 0;
0097   double tau_p = 0;
0098   double tau_pt = 0;
0099
0100   uint tau_decay_prong = 0;
0101   uint tau_decay_hcharged = 0;
0102   uint tau_decay_lcharged = 0;
0103
0104   double quark_eta = 0;
0105   double quark_phi = 0;
0106   double quark_e = 0;
0107   double quark_p = 0;
0108   double quark_pt = 0;
0109
0110   float pt_miss = 0;
0111   float pt_miss_phi = 0;
0112
0113   /* Calculate missing transverse momentum in event */
0114   truth.FindMissingPt( pt_miss, pt_miss_phi );
0115
0116   /* If TAU in event: update tau information */
0117   if( particle_tau )
0118     {
0119       tau_found = true;
0120       tau_eta = particle_tau->momentum().eta();
0121       tau_phi = particle_tau->momentum().phi();
0122
0123       /* Event record uses 0 < phi < 2Pi, while Fun4All uses -Pi < phi < Pi.
0124        * Therefore, correct angle for 'wraparound' at phi == Pi */
0125       if(tau_phi>TMath::Pi()) tau_phi = tau_phi-2*TMath::Pi();
0126
0127       tau_e = particle_tau->momentum().e();
0128       tau_p = sqrt( pow( particle_tau->momentum().px(), 2 ) +
0129                     pow( particle_tau->momentum().py(), 2 ) +
0130                     pow( particle_tau->momentum().pz(), 2 ) );
0131       tau_pt = particle_tau->momentum().perp();
0132
0133       /* Count how many charged particles (hadrons and leptons) the tau particle decays into. */
0134       truth.FindDecayParticles( particle_tau, tau_decay_prong, tau_decay_hcharged, tau_decay_lcharged );
0135
0136     }
0137
0138   /* If QUARK (->jet) in event: Tag the tau candidate (i.e. jet) with smalles delta_R from this quark */
0139   /* @TODO: This is a copy from the loop to tag the tau candidate with smallest delta_R w.r.t. final state tau-
0140    * instead of copying the code, make a function that can be called for both (GetMinDeltaRElement or similar) */
0141   if( particle_quark )
0142     {
0143       quark_found = true;
0144       quark_eta = particle_quark->momentum().eta();
0145       quark_phi = particle_quark->momentum().phi();
0146
0147       /* Event record uses 0 < phi < 2Pi, while Fun4All uses -Pi < phi < Pi.
0148        * Therefore, correct angle for 'wraparound' at phi == Pi */
0149       if(quark_phi>TMath::Pi()) quark_phi = quark_phi-2*TMath::Pi();
0150
0151       quark_e  = particle_quark->momentum().e();
0152       quark_p  = sqrt( pow( particle_quark->momentum().px(), 2 ) +
0153                        pow( particle_quark->momentum().py(), 2 ) +
0154                        pow( particle_quark->momentum().pz(), 2 ) );
0155       quark_pt = particle_quark->momentum().perp();
0156
0157     }
0158
0159   /* Fill event information to ntupple */
0160   float event_data[18] = { (float) _ievent,
0161                            (float) pt_miss,
0162                            (float) pt_miss_phi,
0163                            (float) tau_found,
0164                            (float) quark_found,
0165                            (float) tau_eta,
0166                            (float) tau_phi,
0167                            (float) tau_e,
0168                            (float) tau_p,
0169                            (float) tau_pt,
0170                            (float) tau_decay_prong,
0171                            (float) tau_decay_hcharged,
0172                            (float) tau_decay_lcharged,
0173                            (float) quark_eta,
0174                            (float) quark_phi,
0175                            (float) quark_e,
0176                            (float) quark_p,
0177                            (float) quark_pt
0178   };
0179
0180   _tree_event->Fill( event_data );
0181
0182   /* count up event number */
0183   _ievent ++;
0184
0185   return 0;
0186 }
0187
0188 int
0189 Leptoquarks::End(PHCompositeNode *topNode)
0190 {
0191   _fout_root->cd();
0192   _tree_event->Write();
0193   _fout_root->Close();
0194
0195   return 0;
0196 }
0197
0198
0199 void
0200 Leptoquarks::UpdateFinalStateParticle( HepMC::GenParticle *&particle )
0201 {
0202   bool final_state = false;
0203   while ( !final_state )
0204     {
0205       /* Loop over all children at the end vertex of this particle */
0206       for ( HepMC::GenVertex::particle_iterator child
0207               = particle->end_vertex()->particles_begin(HepMC::children);
0208             child != particle->end_vertex()->particles_end(HepMC::children);
0209             ++child )
0210         {
0211           /* If there is a child of same particle ID, this was not the final state particle- update pointer to particle and repeat */
0212           if ( (*child)->pdg_id() == particle->pdg_id() )
0213             {
0214               particle = (*child);
0215               break;
0216             }
0217           final_state = true;
0218         }
0219     }
0220   return;
0221 }