analysis/EICAnalysis/LeptoquarksReco.C

0001 #include "LeptoquarksReco.h"
0002 #include "PidCandidatev1.h"
0003 #include "TruthTrackerHepMC.h"
0004
0005 /* STL includes */
0006 #include <cassert>
0007
0008 /* Fun4All includes */
0009 #include <trackbase_historic/SvtxTrackMap_v1.h>
0010 #include <trackbase_historic/SvtxVertexMap_v1.h>
0011 //#include <trackbase_historic/SvtxClusterMap_v1.h>
0012
0013 #include <phool/getClass.h>
0014
0015 #include <fun4all/Fun4AllReturnCodes.h>
0016
0017 #include <g4jets/JetMap.h>
0018
0019 #include <calobase/RawTowerGeomContainer.h>
0020 #include <calobase/RawTowerContainer.h>
0021 #include <calobase/RawTowerGeom.h>
0022 #include <calobase/RawTowerv1.h>
0023
0024 #include <g4vertex/GlobalVertexMap.h>
0025 #include <g4vertex/GlobalVertex.h>
0026
0027 #include <g4main/PHG4VtxPoint.h>
0028 #include <g4main/PHG4Particle.h>
0029
0030 #include <g4eval/CaloRawTowerEval.h>
0031 #include <g4eval/SvtxEvalStack.h>
0032
0033
0034 /* ROOT includes */
0035 #include <TLorentzVector.h>
0036 #include <TString.h>
0037 #include <TNtuple.h>
0038 #include <TTree.h>
0039 #include <TFile.h>
0040
0041 using namespace std;
0042
0043 LeptoquarksReco::LeptoquarksReco(std::string filename) :
0044   SubsysReco("LeptoquarksReco" ),
0045   _save_towers(false),
0046   _save_tracks(false),
0047   _ievent(0),
0048   _filename(filename),
0049   _tfile(nullptr),
0050   _t_event(nullptr),
0051   _ntp_tower(nullptr),
0052   _ntp_track(nullptr),
0053   _ebeam_E(0),
0054   _pbeam_E(0),
0055   _tau_jet_emin(5.0),
0056   _jetcolname("AntiKt_Tower_r05")
0057 {
0058
0059 }
0060
0061 int
0062 LeptoquarksReco::Init(PHCompositeNode *topNode)
0063 {
0064   _ievent = 0;
0065
0066   _tfile = new TFile(_filename.c_str(), "RECREATE");
0067
0068   /* Add PidCandidate properties to map that defines output tree */
0069   float dummy = 0;
0070   vector< float > vdummy;
0071
0072   _map_tau_candidate_branches.insert( make_pair( PidCandidate::jet_id , vdummy ) );
0073   _map_tau_candidate_branches.insert( make_pair( PidCandidate::evtgen_pid , vdummy ) );
0074   _map_tau_candidate_branches.insert( make_pair( PidCandidate::evtgen_etotal , vdummy ) );
0075   _map_tau_candidate_branches.insert( make_pair( PidCandidate::evtgen_ptotal , vdummy ) );
0076   _map_tau_candidate_branches.insert( make_pair( PidCandidate::evtgen_theta , vdummy ) );
0077   _map_tau_candidate_branches.insert( make_pair( PidCandidate::evtgen_eta , vdummy ) );
0078   _map_tau_candidate_branches.insert( make_pair( PidCandidate::evtgen_phi , vdummy ) );
0079   _map_tau_candidate_branches.insert( make_pair( PidCandidate::evtgen_decay_prong , vdummy ) );
0080   _map_tau_candidate_branches.insert( make_pair( PidCandidate::evtgen_decay_hcharged , vdummy ) );
0081   _map_tau_candidate_branches.insert( make_pair( PidCandidate::evtgen_decay_lcharged , vdummy ) );
0082   _map_tau_candidate_branches.insert( make_pair( PidCandidate::jet_eta , vdummy ) );
0083   _map_tau_candidate_branches.insert( make_pair( PidCandidate::jet_phi , vdummy ) );
0084   _map_tau_candidate_branches.insert( make_pair( PidCandidate::jet_etotal , vdummy ) );
0085   _map_tau_candidate_branches.insert( make_pair( PidCandidate::jet_etrans , vdummy ) );
0086   _map_tau_candidate_branches.insert( make_pair( PidCandidate::jet_ptotal , vdummy ) );
0087   _map_tau_candidate_branches.insert( make_pair( PidCandidate::jet_ptrans , vdummy ) );
0088   _map_tau_candidate_branches.insert( make_pair( PidCandidate::jet_minv , vdummy ) );
0089   _map_tau_candidate_branches.insert( make_pair( PidCandidate::jet_mtrans , vdummy ) );
0090   _map_tau_candidate_branches.insert( make_pair( PidCandidate::jet_ncomp , vdummy ) );
0091   _map_tau_candidate_branches.insert( make_pair( PidCandidate::jet_ncomp_above_0p1 , vdummy ) );
0092   _map_tau_candidate_branches.insert( make_pair( PidCandidate::jet_ncomp_above_1 , vdummy ) );
0093   _map_tau_candidate_branches.insert( make_pair( PidCandidate::jet_ncomp_above_10 , vdummy ) );
0094   _map_tau_candidate_branches.insert( make_pair( PidCandidate::jet_ncomp_emcal , vdummy ) );
0095   _map_tau_candidate_branches.insert( make_pair( PidCandidate::jetshape_econe_r01 , vdummy ) );
0096   _map_tau_candidate_branches.insert( make_pair( PidCandidate::jetshape_econe_r02 , vdummy ) );
0097   _map_tau_candidate_branches.insert( make_pair( PidCandidate::jetshape_econe_r03 , vdummy ) );
0098   _map_tau_candidate_branches.insert( make_pair( PidCandidate::jetshape_econe_r04 , vdummy ) );
0099   _map_tau_candidate_branches.insert( make_pair( PidCandidate::jetshape_econe_r05 , vdummy ) );
0100   _map_tau_candidate_branches.insert( make_pair( PidCandidate::jetshape_r90 , vdummy ) );
0101   _map_tau_candidate_branches.insert( make_pair( PidCandidate::jetshape_rms , vdummy ) );
0102   _map_tau_candidate_branches.insert( make_pair( PidCandidate::jetshape_radius , vdummy ) );
0103   _map_tau_candidate_branches.insert( make_pair( PidCandidate::jetshape_emcal_econe_r01 , vdummy ) );
0104   _map_tau_candidate_branches.insert( make_pair( PidCandidate::jetshape_emcal_econe_r02 , vdummy ) );
0105   _map_tau_candidate_branches.insert( make_pair( PidCandidate::jetshape_emcal_econe_r03 , vdummy ) );
0106   _map_tau_candidate_branches.insert( make_pair( PidCandidate::jetshape_emcal_econe_r04 , vdummy ) );
0107   _map_tau_candidate_branches.insert( make_pair( PidCandidate::jetshape_emcal_econe_r05 , vdummy ) );
0108   _map_tau_candidate_branches.insert( make_pair( PidCandidate::jetshape_emcal_r90 , vdummy ) );
0109   _map_tau_candidate_branches.insert( make_pair( PidCandidate::jetshape_emcal_rms , vdummy ) );
0110   _map_tau_candidate_branches.insert( make_pair( PidCandidate::jetshape_emcal_radius , vdummy ) );
0111   _map_tau_candidate_branches.insert( make_pair( PidCandidate::tracks_count_r02 , vdummy ) );
0112   _map_tau_candidate_branches.insert( make_pair( PidCandidate::tracks_chargesum_r02 , vdummy ) );
0113   _map_tau_candidate_branches.insert( make_pair( PidCandidate::tracks_rmax_r02 , vdummy ) );
0114   _map_tau_candidate_branches.insert( make_pair( PidCandidate::tracks_count_r04 , vdummy ) );
0115   _map_tau_candidate_branches.insert( make_pair( PidCandidate::tracks_chargesum_r04 , vdummy ) );
0116   _map_tau_candidate_branches.insert( make_pair( PidCandidate::tracks_rmax_r04 , vdummy ) );
0117   _map_tau_candidate_branches.insert( make_pair( PidCandidate::tracks_count_R , vdummy ) );
0118   _map_tau_candidate_branches.insert( make_pair( PidCandidate::tracks_chargesum_R , vdummy ) );
0119   _map_tau_candidate_branches.insert( make_pair( PidCandidate::tracks_rmax_R , vdummy ) );
0120   _map_tau_candidate_branches.insert( make_pair( PidCandidate::tracks_vertex , vdummy ) );
0121
0122   /* Add branches to map that defines output tree for event-wise properties */
0123   _map_event_branches.insert( make_pair( "Et_miss" , dummy ) );
0124   _map_event_branches.insert( make_pair( "Et_miss_phi" , dummy ) );
0125   _map_event_branches.insert( make_pair( "reco_tau_found" , dummy ) );
0126   _map_event_branches.insert( make_pair( "reco_tau_is_tau" , dummy ) );
0127   _map_event_branches.insert( make_pair( "reco_tau_eta" , dummy ) );
0128   _map_event_branches.insert( make_pair( "reco_tau_phi" , dummy ) );
0129   _map_event_branches.insert( make_pair( "reco_tau_ptotal" , dummy ) );
0130   _map_event_branches.insert( make_pair( "is_lq_event" , dummy ) );
0131
0132
0133   /* Create tree for information about full event */
0134   _t_event = new TTree("event", "a Tree with global event information and tau candidates");
0135   _t_event->Branch("event", &_ievent, "event/I");
0136
0137   /* Add event branches */
0138   for ( map< string , float >::iterator iter = _map_event_branches.begin();
0139         iter != _map_event_branches.end();
0140         ++iter)
0141     {
0142       _t_event->Branch( (iter->first).c_str(),
0143                         &(iter->second) );
0144     }
0145
0146   /* Add tau candidate branches */
0147   for ( map< PidCandidate::PROPERTY , vector< float > >::iterator iter = _map_tau_candidate_branches.begin();
0148         iter != _map_tau_candidate_branches.end();
0149         ++iter)
0150     {
0151       _t_event->Branch(PidCandidate::get_property_info( (iter->first) ).first.c_str(),
0152                &(iter->second) );
0153     }
0154
0155
0156   /* Create additional tree for towers */
0157   if ( _save_towers )
0158     {
0159       _ntp_tower = new TNtuple("ntp_tower","towers from all all tau candidates",
0160                                "ievent:jet_id:evtgen_pid:evtgen_etotal:evtgen_eta:evtgen_phi:evtgen_decay_prong:evtgen_decay_hcharged:evtgen_decay_lcharged:jet_eta:jet_phi:jet_etotal:tower_calo_id:tower_eta:tower_phi:tower_delta_r:tower_e");
0161     }
0162
0163   /* Create additional tree for tracks */
0164   if ( _save_tracks )
0165     {
0166       _ntp_track = new TNtuple("ntp_track","tracks from all all tau candidates",
0167                                "ievent:jet_id:evtgen_pid:evtgen_etotal:evtgen_eta:evtgen_phi:evtgen_decay_prong:evtgen_decay_hcharged:evtgen_decay_lcharged:jet_eta:jet_phi:jet_etotal:track_quality:track_eta:track_phi:track_delta_r:track_p");
0168     }
0169
0170   return 0;
0171 }
0172
0173 int
0174 LeptoquarksReco::process_event(PHCompositeNode *topNode)
0175 {
0176   /* Reset branch map */
0177   ResetBranchMap();
0178
0179   /* Create map to collect tau candidates.
0180    * Use energy as 'key' to the map because energy is unique for each jet, while there are sometimes multiple jets (same energy,
0181    * different jet ID) in the input jet collection. Also, map automatically sorts entries by key, i.e. this gives list of tau candidates
0182    * sorted by energy. */
0183   type_map_tcan tauCandidateMap;
0184
0185   /* Get reco jets collection */
0186   JetMap* recojets = findNode::getClass<JetMap>(topNode,_jetcolname.c_str());
0187   if (!recojets)
0188     {
0189       cerr << PHWHERE << " ERROR: Can't find " << _jetcolname << endl;
0190       return Fun4AllReturnCodes::ABORTEVENT;
0191     }
0192
0193   /* Get track collection with all tracks in this event */
0194   SvtxTrackMap* trackmap =
0195     findNode::getClass<SvtxTrackMap>(topNode,"SvtxTrackMap");
0196   if (!trackmap)
0197     {
0198       cout<< PHWHERE << "SvtxTrackMap node not found on node tree"
0199            << endl;
0200       return Fun4AllReturnCodes::ABORTEVENT;
0201     }
0202
0203   /* Get vertex collection with all vertices in this event */
0204   SvtxVertexMap* vertexmap =
0205     findNode::getClass<SvtxVertexMap>(topNode,"SvtxVertexMap");
0206   if (!vertexmap)
0207     {
0208       cout<< PHWHERE << "SvtxVertexMap node not found on node tree"
0209           << endl;
0210       return Fun4AllReturnCodes::ABORTEVENT;
0211     }
0212
0213
0214     /* Get collection of truth particles from event generator */
0215   PHHepMCGenEventMap *genevtmap = findNode::getClass<PHHepMCGenEventMap>(topNode,"PHHepMCGenEventMap");
0216   if (!genevtmap) {
0217     cerr << PHWHERE << " ERROR: Can't find PHHepMCGenEventMap" << endl;
0218     return Fun4AllReturnCodes::ABORTEVENT;
0219   }
0220
0221   /* create map of all tower collections and tower geometry collections used */
0222   type_map_cdata map_calotower;
0223
0224   /* Select calorimeter to include */
0225   vector< RawTowerDefs::CalorimeterId > v_caloids;
0226   v_caloids.push_back( RawTowerDefs::CEMC );
0227   v_caloids.push_back( RawTowerDefs::HCALIN );
0228   v_caloids.push_back( RawTowerDefs::HCALOUT );
0229   v_caloids.push_back( RawTowerDefs::FEMC );
0230   v_caloids.push_back( RawTowerDefs::FHCAL );
0231   v_caloids.push_back( RawTowerDefs::EEMC );
0232
0233   /* Fill map with calorimeter data */
0234   for ( unsigned i = 0; i < v_caloids.size(); i++ )
0235     {
0236       /* Get collection of towers */
0237       string towers_name = "TOWER_CALIB_";
0238       towers_name.append( RawTowerDefs::convert_caloid_to_name( v_caloids.at( i ) ) );
0239
0240       RawTowerContainer *towers = findNode::getClass<RawTowerContainer>(topNode,towers_name);
0241       if (!towers) {
0242         cerr << PHWHERE << " ERROR: Can't find " << towers_name << endl;
0243         return Fun4AllReturnCodes::ABORTEVENT;
0244       }
0245
0246       /* get geometry collection */
0247       string towergeom_name = "TOWERGEOM_";
0248       towergeom_name.append( RawTowerDefs::convert_caloid_to_name( v_caloids.at( i ) ) );
0249
0250       RawTowerGeomContainer *geom = findNode::getClass<RawTowerGeomContainer>(topNode,towergeom_name);
0251       if (!geom) {
0252         cerr << PHWHERE << " ERROR: Can't find " << towergeom_name << endl;
0253         return Fun4AllReturnCodes::ABORTEVENT;
0254       }
0255
0256       /* Insert tower and geometry collections in map */
0257       map_calotower.insert( make_pair( v_caloids.at( i ), make_pair( towers, geom ) ) );
0258     }
0259
0260
0261   /* Loop over all jets to collect tau candidate objects */
0262   for (JetMap::Iter iter = recojets->begin();
0263        iter != recojets->end();
0264        ++iter)
0265     {
0266       /* skip jets below energy threshold */
0267       if ( (iter->second)->get_e() < _tau_jet_emin )
0268         continue;
0269
0270       /* create new tau candidate */
0271       PidCandidatev1 *tc = new PidCandidatev1();
0272       tc->set_candidate_id( (iter->second)->get_id() );
0273       tc->set_property( PidCandidate::jet_id , (iter->second)->get_id() );
0274
0275       tc->set_property( PidCandidate::evtgen_pid, (int)0 );
0276
0277       /* add tau candidate to collection */
0278       tauCandidateMap.insert( make_pair( (iter->second)->get_e(), tc ) );
0279     }
0280
0281   /* vertex information */
0282   SvtxEvalStack *svtxevalstack = new SvtxEvalStack(topNode);
0283
0284   /* Add jet information to tau candidates */
0285   AddJetInformation( tauCandidateMap, recojets, &map_calotower );
0286
0287   /* Add jet structure information to tau candidates */
0288   AddJetStructureInformation( tauCandidateMap, &map_calotower );
0289
0290   /* Add track information to tau candidates */
0291   AddTrackInformation( tauCandidateMap, trackmap, vertexmap, svtxevalstack, recojets->get_par());
0292
0293   /* Add tag for true Tau particle jet to tau candidates */
0294   AddTrueTauTag( tauCandidateMap, genevtmap );
0295
0296   /* Add information about tau candidats to output tree */
0297   WritePidCandidatesToTree( tauCandidateMap );
0298
0299   /* Add global event information to separate tree */
0300   AddGlobalEventInformation( tauCandidateMap, &map_calotower );
0301
0302   /* fill event information tree */
0303   _t_event->Fill();
0304
0305   /* count up event number */
0306   _ievent ++;
0307
0308   return 0;
0309 }
0310
0311 // Returns average of values while excluding outliers //
0312 float Average(vector<float> v)
0313 {
0314   double sum = 0;
0315   int med_i = v.size()/2;
0316   vector<float> temp_v;
0317   for(int i=0;(unsigned)i<v.size();i++)
0318     if (v[i] > (v[med_i] - v[med_i]*0.2) && v[i] < (v[med_i] + v[med_i]*0.2)) temp_v.push_back(v[i]);
0319
0320   for(int j=0;(unsigned)j<temp_v.size();j++) sum = sum + temp_v[j];
0321
0322   return sum/temp_v.size();
0323 }
0324
0325 int
0326 LeptoquarksReco::AddTrueTauTag( type_map_tcan& tauCandidateMap, PHHepMCGenEventMap *genevtmap )
0327 {
0328   /* Look for leptoquark and tau particle in the event */
0329   TruthTrackerHepMC truth;
0330   truth.set_hepmc_geneventmap( genevtmap );
0331
0332   int pdg_lq = 39; // leptoquark
0333   int pdg_tau = 15; // tau lepton
0334   int pdg_parton = 0;
0335   int pdg_electron = 11; // electron
0336
0337   // Check if electron exists  //
0338   HepMC::GenParticle* particle_electron = NULL;
0339
0340   for (PHHepMCGenEventMap::ReverseIter iter = genevtmap->rbegin(); iter != genevtmap->rend(); ++iter)
0341     {
0342       PHHepMCGenEvent *genevt = iter->second;
0343       HepMC::GenEvent *theEvent = genevt->getEvent();
0344
0345       // Loop over all truth particles in event generator collection
0346       for ( HepMC::GenEvent::particle_iterator p = theEvent->particles_begin();
0347         p != theEvent->particles_end(); ++p ) {
0348     // If final state electron then tag it //
0349     if((*p)->pdg_id() == pdg_electron && (*p)->status() == 1) particle_electron = (*p);
0350       }
0351     }
0352
0353   /* Search for leptoquark in event */
0354   HepMC::GenParticle* particle_lq = truth.FindParticle( pdg_lq );
0355
0356   // Tag event as LQ event or not //
0357   if (particle_lq) ( _map_event_branches.find( "is_lq_event" ) )->second = 1;
0358   else ( _map_event_branches.find( "is_lq_event" ) )->second = 0;
0359
0360   /* Search for lq->tau decay in event */
0361   HepMC::GenParticle* particle_tau = NULL;
0362
0363   // Find tau that comes from LQ or just any tau in event //
0364   if(particle_lq) particle_tau = truth.FindDaughterParticle( pdg_tau, particle_lq );
0365   else particle_tau = truth.FindParticle(pdg_tau);
0366
0367   /* Search for lq->quark decay in event.
0368    * Loop over all quark PDG codes until finding a matching quark. */
0369   HepMC::GenParticle* particle_quark = NULL;
0370   for ( int pdg_quark = 1; pdg_quark < 7; pdg_quark++ )
0371     {
0372       /* try quark */
0373       particle_quark = truth.FindDaughterParticle( pdg_quark, particle_lq );
0374       if (particle_quark)
0375         {
0376           pdg_parton = pdg_quark;
0377           break;
0378         }
0379
0380       /* try anti-quark */
0381       particle_quark = truth.FindDaughterParticle( -pdg_quark, particle_lq );
0382       if (particle_quark)
0383         {
0384           pdg_parton = -pdg_quark;
0385           break;
0386         }
0387     }
0388
0389   /* If TAU in event: Tag the tau candidate (i.e. jet) with smalles delta_R from this tau */
0390   if( particle_tau )
0391     {
0392       PidCandidate* best_match = FindMinDeltaRCandidate( &tauCandidateMap,
0393                                                          particle_tau->momentum().eta(),
0394                                                          particle_tau->momentum().phi() );
0395
0396
0397       /* set is_tau = TRUE for PidCandiate with smallest delta_R within reasonable range*/
0398       if ( best_match )
0399         {
0400       /* Check: If PidCandidate::Evtgen_pid has already been set to a value != 0, exit function. */
0401       if( best_match->get_property_int( PidCandidate::evtgen_pid ) != 0 )
0402         {
0403           cout << "ERROR: Try to set PidCandidate::evtgen_pid for PidCandidate with evtgen_pid != 0" << endl;
0404           return -1;
0405         }
0406
0407           /* Update PidCandidate entry */
0408           best_match->set_property( PidCandidate::evtgen_pid, pdg_tau );
0409           best_match->set_property( PidCandidate::evtgen_etotal, (float)particle_tau->momentum().e() );
0410           best_match->set_property( PidCandidate::evtgen_eta, (float)particle_tau->momentum().eta() );
0411           best_match->set_property( PidCandidate::evtgen_phi, (float)particle_tau->momentum().phi() );
0412
0413           /* Check particle decay if end-vertex found */
0414           if ( particle_tau->end_vertex() )
0415             {
0416               /* Add information about tau decay */
0417               uint tau_decay_prong = 0;
0418               uint tau_decay_hcharged = 0;
0419               uint tau_decay_lcharged = 0;
0420
0421           /* Count how many charged particles (hadrons and leptons) a given particle decays into. */
0422               truth.FindDecayParticles( particle_tau, tau_decay_prong, tau_decay_hcharged, tau_decay_lcharged );
0423
0424
0425               /* Update tau candidate entry */
0426               best_match->set_property( PidCandidate::evtgen_decay_prong, tau_decay_prong );
0427               best_match->set_property( PidCandidate::evtgen_decay_hcharged, tau_decay_hcharged );
0428               best_match->set_property( PidCandidate::evtgen_decay_lcharged, tau_decay_lcharged );
0429             }
0430         }
0431     }
0432
0433   /* If QUARK (->jet) in event: Tag the tau candidate (i.e. jet) with smalles delta_R from this quark */
0434   if( particle_quark )
0435     {
0436       PidCandidate* best_match = FindMinDeltaRCandidate( &tauCandidateMap,
0437                                                          particle_quark->momentum().eta(),
0438                                                          particle_quark->momentum().phi() );
0439
0440       /* set is_uds = TRUE for PidCandiate with smallest delta_R if found */
0441       if ( best_match )
0442         {
0443       /* Check: If PidCandidate::Evtgen_pid has already been set to a value != 0, exit function. */
0444       if( best_match->get_property_int( PidCandidate::evtgen_pid ) != 0 )
0445         {
0446           cout << "ERROR: Try to set PidCandidate::evtgen_pid for PidCandidate with evtgen_pid != 0" << endl;
0447           return -1;
0448         }
0449       /* Set properties of PidCandidate */
0450           best_match->set_property( PidCandidate::evtgen_pid, pdg_parton );
0451           best_match->set_property( PidCandidate::evtgen_etotal, (float)particle_quark->momentum().e() );
0452           best_match->set_property( PidCandidate::evtgen_eta, (float)particle_quark->momentum().eta() );
0453           best_match->set_property( PidCandidate::evtgen_phi, (float)particle_quark->momentum().phi() );
0454         }
0455     }
0456
0457   if( particle_electron )
0458     {
0459       //cout<<"ELECTRON"<<endl;
0460       PidCandidate* best_match = FindMinDeltaRCandidate( &tauCandidateMap,
0461                                                          particle_electron->momentum().eta(),
0462                                                          particle_electron->momentum().phi() );
0463
0464       /* set electron = TRUE for PidCandiate with smallest delta_R if found */
0465       if ( best_match )
0466         {
0467           /* Check: If PidCandidate::Evtgen_pid has already been set to a value != 0, exit function. */
0468           if( best_match->get_property_int( PidCandidate::evtgen_pid ) != 0 )
0469             {
0470               cout << "ERROR: Try to set PidCandidate::evtgen_pid for PidCandidate with evtgen_pid != 0" << endl;
0471               return -1;
0472             }
0473
0474           /* Set properties of PidCandidate */
0475           best_match->set_property( PidCandidate::evtgen_pid, pdg_electron );
0476           best_match->set_property( PidCandidate::evtgen_etotal, (float)particle_electron->momentum().e() );
0477           best_match->set_property( PidCandidate::evtgen_eta, (float)particle_electron->momentum().eta() );
0478           best_match->set_property( PidCandidate::evtgen_phi, (float)particle_electron->momentum().phi() );
0479         }
0480     }
0481
0482
0483   return 0;
0484 }
0485
0486
0487 int
0488 LeptoquarksReco::AddJetInformation( type_map_tcan& tauCandidateMap, JetMap* recojets, type_map_cdata* map_calotower )
0489 {
0490   /* Loop over tau candidates */
0491   for (type_map_tcan::iterator iter = tauCandidateMap.begin();
0492        iter != tauCandidateMap.end();
0493        ++iter)
0494     {
0495       Jet* jetx = recojets->get( (iter->second)->get_property_uint( PidCandidate::jet_id ) );
0496
0497       /* calculate transverse mass of jet */
0498       float jet_mtrans = sqrt( pow( jetx->get_mass(), 2 ) +
0499                                pow( jetx->get_pt(), 2 ) );
0500
0501       /* count jet ncomp above thresholds */
0502       unsigned int jet_ncomp_above_0p1 = 0;
0503       unsigned int jet_ncomp_above_1 = 0;
0504       unsigned int jet_ncomp_above_10 = 0;
0505
0506       for (Jet::ConstIter jcompiter = jetx->begin_comp(); jcompiter != jetx->end_comp(); ++jcompiter)
0507         {
0508           RawTowerDefs::CalorimeterId calo_id = RawTowerDefs::NONE;
0509
0510           switch ( jcompiter->first )
0511             {
0512             case Jet::CEMC_TOWER:
0513               calo_id = RawTowerDefs::CEMC;
0514               break;
0515             case Jet::HCALIN_TOWER:
0516               calo_id = RawTowerDefs::HCALIN;
0517               break;
0518             case Jet::HCALOUT_TOWER:
0519               calo_id = RawTowerDefs::HCALOUT;
0520               break;
0521             default:
0522               break;
0523             }
0524
0525           /* continue if no calorimeter id found */
0526           if ( calo_id == RawTowerDefs::NONE )
0527             continue;
0528
0529           /* get tower container from map, find tower in tower container, get tower energy */
0530           float e_component = 0;
0531           if ( map_calotower->find( calo_id ) != map_calotower->end() )
0532             e_component = ( ( ( map_calotower->find( calo_id ) )->second ).first )->getTower( jcompiter->second )->get_energy();
0533
0534           /* check if energy is above threshold and count up matching counters accordingly */
0535           if ( e_component > 0.1 )
0536             jet_ncomp_above_0p1++;
0537           if ( e_component > 1 )
0538             jet_ncomp_above_1++;
0539           if ( e_component > 10 )
0540             jet_ncomp_above_10++;
0541         }
0542
0543       /* set tau candidate jet properties */
0544       (iter->second)->set_property( PidCandidate::jet_eta , jetx->get_eta() );
0545       (iter->second)->set_property( PidCandidate::jet_phi , jetx->get_phi() );
0546       (iter->second)->set_property( PidCandidate::jet_etotal , jetx->get_e() );
0547       (iter->second)->set_property( PidCandidate::jet_etrans , jetx->get_et() );
0548       (iter->second)->set_property( PidCandidate::jet_ptotal , jetx->get_p() );
0549       (iter->second)->set_property( PidCandidate::jet_ptrans , jetx->get_pt() );
0550       (iter->second)->set_property( PidCandidate::jet_minv , jetx->get_mass() );
0551       (iter->second)->set_property( PidCandidate::jet_mtrans , jet_mtrans );
0552       (iter->second)->set_property( PidCandidate::jet_ncomp , (uint)jetx->size_comp() );
0553       (iter->second)->set_property( PidCandidate::jet_ncomp_above_0p1 , jet_ncomp_above_0p1 );
0554       (iter->second)->set_property( PidCandidate::jet_ncomp_above_1 , jet_ncomp_above_1 );
0555       (iter->second)->set_property( PidCandidate::jet_ncomp_above_10 , jet_ncomp_above_10 );
0556       (iter->second)->set_property( PidCandidate::jet_ncomp_emcal , (uint)jetx->count_comp( Jet::CEMC_TOWER ) );
0557     }
0558
0559   return 0;
0560 }
0561
0562
0563
0564 int
0565 LeptoquarksReco::AddJetStructureInformation( type_map_tcan& tauCandidateMap, type_map_cdata* map_towers )
0566 {
0567   /* Cone size around jet axis within which to look for tracks */
0568   float delta_R_cutoff_r1 = 0.1;
0569   float delta_R_cutoff_r2 = 0.2;
0570   float delta_R_cutoff_r3 = 0.3;
0571   float delta_R_cutoff_r4 = 0.4;
0572   float delta_R_cutoff_r5 = 0.5;
0573
0574   /* energy threshold for considering tower */
0575   float tower_emin = 0.0;
0576
0577   /* number of steps for finding r90 */
0578   int n_steps = 50;
0579
0580   /* Loop over tau candidates */
0581   for (type_map_tcan::iterator iter = tauCandidateMap.begin();
0582        iter != tauCandidateMap.end();
0583        ++iter)
0584     {
0585       /* get jet axis */
0586       float jet_eta = (iter->second)->get_property_float( PidCandidate::jet_eta );
0587       float jet_phi = (iter->second)->get_property_float( PidCandidate::jet_phi );
0588       float jet_e =   (iter->second)->get_property_float( PidCandidate::jet_etotal );
0589
0590       /* collect jet structure properties */
0591       float er1 = 0;
0592       float er2 = 0;
0593       float er3 = 0;
0594       float er4 = 0;
0595       float er5 = 0;
0596       float r90 = 0;
0597       float radius = 0;
0598       float rms = 0;
0599       float rms_esum = 0;
0600
0601       /* collect jet structure properties- EMCal ONLY */
0602       float emcal_er1 = 0;
0603       float emcal_er2 = 0;
0604       float emcal_er3 = 0;
0605       float emcal_er4 = 0;
0606       float emcal_er5 = 0;
0607       float emcal_r90 = 0;
0608       float emcal_radius = 0;
0609       float emcal_rms = 0;
0610       float emcal_rms_esum = 0;
0611
0612       /* Loop over all tower (and geometry) collections */
0613       for (type_map_cdata::iterator iter_calo = map_towers->begin();
0614            iter_calo != map_towers->end();
0615            ++iter_calo)
0616         {
0617           /* define tower iterator */
0618           RawTowerContainer::ConstRange begin_end = ((iter_calo->second).first)->getTowers();
0619           RawTowerContainer::ConstIterator rtiter;
0620
0621           /* loop over all tower in CEMC calorimeter */
0622           for (rtiter = begin_end.first; rtiter !=  begin_end.second; ++rtiter)
0623             {
0624               /* get tower energy */
0625               RawTower *tower = rtiter->second;
0626               float tower_energy = tower->get_energy();
0627
0628               /* check if tower above energy treshold */
0629               if ( tower_energy < tower_emin )
0630                 continue;
0631
0632               /* get eta and phi of tower and check angle delta_R w.r.t. jet axis */
0633               RawTowerGeom * tower_geom = ((iter_calo->second).second)->get_tower_geometry(tower -> get_key());
0634               float tower_eta = tower_geom->get_eta();
0635               float tower_phi = tower_geom->get_phi();
0636
0637               /* If accounting for displaced vertex, need to calculate eta and phi:
0638                  double r = tower_geom->get_center_radius();
0639                  double phi = atan2(tower_geom->get_center_y(), tower_geom->get_center_x());
0640                  double z0 = tower_geom->get_center_z();
0641                  double z = z0 - vtxz;
0642                  double eta = asinh(z/r); // eta after shift from vertex
0643               */
0644
0645               float delta_R = CalculateDeltaR( tower_eta , tower_phi, jet_eta, jet_phi );
0646
0647               /* if save_towers set true: add tower to tree */
0648               if ( _save_towers )
0649                 {
0650                   float tower_data[17] = {(float) _ievent,
0651                                           (float) (iter->second)->get_property_uint( PidCandidate::jet_id ),
0652                                           (float) (iter->second)->get_property_int( PidCandidate::evtgen_pid ),
0653                                           (float) (iter->second)->get_property_float( PidCandidate::evtgen_etotal ),
0654                                           (float) (iter->second)->get_property_float( PidCandidate::evtgen_eta ),
0655                                           (float) (iter->second)->get_property_float( PidCandidate::evtgen_phi ),
0656                                           (float) (iter->second)->get_property_uint( PidCandidate::evtgen_decay_prong ),
0657                                           (float) (iter->second)->get_property_uint( PidCandidate::evtgen_decay_hcharged ),
0658                                           (float) (iter->second)->get_property_uint( PidCandidate::evtgen_decay_lcharged ),
0659                                           (float) (iter->second)->get_property_float( PidCandidate::jet_eta ),
0660                                           (float) (iter->second)->get_property_float( PidCandidate::jet_phi ),
0661                                           (float) (iter->second)->get_property_float( PidCandidate::jet_etotal ),
0662                                           (float) (iter_calo->first),
0663                                           (float) tower_eta,
0664                                           (float) tower_phi,
0665                                           (float) delta_R,
0666                                           (float) tower_energy
0667                   };
0668
0669                   _ntp_tower->Fill(tower_data);
0670                 }
0671
0672               /* check delta R distance tower from jet axis */
0673               if ( delta_R <= delta_R_cutoff_r1 )
0674                 {
0675                   er1 += tower_energy;
0676                   if ( (iter_calo->first) == RawTowerDefs::CEMC )
0677                     emcal_er1 += tower_energy;
0678                 }
0679               if ( delta_R <= delta_R_cutoff_r2 )
0680                 {
0681                   er2 += tower_energy;
0682                   if ( (iter_calo->first) == RawTowerDefs::CEMC )
0683                     emcal_er2 += tower_energy;
0684                 }
0685               if ( delta_R <= delta_R_cutoff_r3 )
0686                 {
0687                   er3 += tower_energy;
0688                   if ( (iter_calo->first) == RawTowerDefs::CEMC )
0689                     emcal_er3 += tower_energy;
0690                 }
0691               if ( delta_R <= delta_R_cutoff_r4 )
0692                 {
0693                   er4 += tower_energy;
0694                   if ( (iter_calo->first) == RawTowerDefs::CEMC )
0695                     emcal_er4 += tower_energy;
0696                 }
0697               if ( delta_R <= delta_R_cutoff_r5 )
0698                 {
0699                   er5 += tower_energy;
0700                   if ( (iter_calo->first) == RawTowerDefs::CEMC )
0701                     emcal_er5 += tower_energy;
0702                 }
0703
0704               if ( delta_R <= delta_R_cutoff_r5 )
0705                 {
0706                   rms += tower_energy*delta_R*delta_R;
0707                   rms_esum += tower_energy;
0708
0709                   radius += tower_energy*delta_R;
0710
0711                   if ( (iter_calo->first) == RawTowerDefs::CEMC )
0712                     {
0713                       emcal_rms += tower_energy*delta_R*delta_R;
0714                       emcal_rms_esum += tower_energy;
0715                       emcal_radius += tower_energy*delta_R;
0716                     }
0717                 }
0718             }
0719         }
0720
0721       /* finalize calculation of rms and radius */
0722       if ( rms_esum > 0 )
0723         {
0724           radius /= rms_esum;
0725           rms    /= rms_esum;
0726           rms     = sqrt( rms );
0727         }
0728       else
0729         {
0730           radius = -1;
0731           rms = -1;
0732         }
0733       if ( emcal_rms_esum > 0 )
0734         {
0735           emcal_radius /= emcal_rms_esum;
0736           emcal_rms    /= emcal_rms_esum;
0737           emcal_rms     = sqrt( emcal_rms );
0738         }
0739       else
0740         {
0741           emcal_radius = -1;
0742           emcal_rms = -1;
0743         }
0744
0745       /* Search for cone angle that contains 90% of jet energy */
0746       for(int r_i = 1; r_i<n_steps+1; r_i++){
0747         float e_tower_sum = 0;
0748
0749         /* Loop over all tower (and geometry) collections */
0750         for (type_map_cdata::iterator iter_calo = map_towers->begin();
0751              iter_calo != map_towers->end();
0752              ++iter_calo)
0753           {
0754             /* define tower iterator */
0755             RawTowerContainer::ConstRange begin_end = ((iter_calo->second).first)->getTowers();
0756             RawTowerContainer::ConstIterator rtiter;
0757
0758             if (e_tower_sum >= 0.9*jet_e) break;
0759
0760             for (rtiter = begin_end.first; rtiter !=  begin_end.second; ++rtiter)
0761               {
0762                 RawTower *tower = rtiter->second;
0763                 float tower_energy = tower->get_energy();
0764
0765                 /* check if tower is above minimum tower energy required */
0766                 if ( tower_energy < tower_emin )
0767                   continue;
0768
0769                 RawTowerGeom * tower_geom = ((iter_calo->second).second)->get_tower_geometry(tower -> get_key());
0770                 float tower_eta = tower_geom->get_eta();
0771                 float tower_phi = tower_geom->get_phi();
0772
0773                 float delta_R = CalculateDeltaR( tower_eta , tower_phi, jet_eta, jet_phi );
0774
0775                 if(delta_R < r_i*delta_R_cutoff_r5/n_steps) {
0776                   e_tower_sum = e_tower_sum + tower_energy;
0777                   r90 = r_i*delta_R_cutoff_r5/n_steps;
0778                 }
0779               }
0780
0781             if (e_tower_sum >= 0.9*jet_e) break;
0782           }
0783       }
0784
0785       /* set tau candidate properties */
0786       (iter->second)->set_property( PidCandidate::jetshape_econe_r01, er1 );
0787       (iter->second)->set_property( PidCandidate::jetshape_econe_r02, er2 );
0788       (iter->second)->set_property( PidCandidate::jetshape_econe_r03, er3 );
0789       (iter->second)->set_property( PidCandidate::jetshape_econe_r04, er4 );
0790       (iter->second)->set_property( PidCandidate::jetshape_econe_r05, er5 );
0791       (iter->second)->set_property( PidCandidate::jetshape_r90, r90 );
0792       (iter->second)->set_property( PidCandidate::jetshape_rms, rms );
0793       (iter->second)->set_property( PidCandidate::jetshape_radius, radius );
0794       (iter->second)->set_property( PidCandidate::jetshape_emcal_econe_r01, emcal_er1 );
0795       (iter->second)->set_property( PidCandidate::jetshape_emcal_econe_r02, emcal_er2 );
0796       (iter->second)->set_property( PidCandidate::jetshape_emcal_econe_r03, emcal_er3 );
0797       (iter->second)->set_property( PidCandidate::jetshape_emcal_econe_r04, emcal_er4 );
0798       (iter->second)->set_property( PidCandidate::jetshape_emcal_econe_r05, emcal_er5 );
0799       (iter->second)->set_property( PidCandidate::jetshape_emcal_r90, emcal_r90 );
0800       (iter->second)->set_property( PidCandidate::jetshape_emcal_rms, emcal_rms );
0801       (iter->second)->set_property( PidCandidate::jetshape_emcal_radius, emcal_radius );
0802     }
0803
0804   return 0;
0805 }
0806
0807 int
0808 LeptoquarksReco::AddTrackInformation( type_map_tcan& tauCandidateMap, SvtxTrackMap* trackmap, SvtxVertexMap* vertexmap, SvtxEvalStack *svtxevalstack, double R_max )
0809 {
0810
0811   // Pointers for tracks //
0812   SvtxTrackEval* trackeval = svtxevalstack->get_track_eval();
0813   SvtxTruthEval* trutheval = svtxevalstack->get_truth_eval();
0814
0815   /* Loop over tau candidates */
0816   for (type_map_tcan::iterator iter = tauCandidateMap.begin();
0817        iter != tauCandidateMap.end();
0818        ++iter)
0819     {
0820       uint tracks_count_r02 = 0;
0821       int tracks_chargesum_r02 = 0;
0822       float tracks_rmax_r02 = 0;
0823
0824       uint tracks_count_r04 = 0;
0825       int tracks_chargesum_r04 = 0;
0826       float tracks_rmax_r04 = 0;
0827
0828       uint tracks_count_R = 0;
0829       int tracks_chargesum_R = 0;
0830       float tracks_rmax_R = 0;
0831
0832       vector<float> tracks_vertex;
0833       vector<float> temp_vertex;
0834
0835       float jet_eta = (iter->second)->get_property_float( PidCandidate::jet_eta );
0836       float jet_phi = (iter->second)->get_property_float( PidCandidate::jet_phi );
0837
0838       /* Loop over tracks
0839        * (float) track->get_eta(),     //eta of the track
0840        * (float) track->get_phi(),     //phi of the track
0841        * (float) track->get_pt(),      //transverse momentum of track
0842        * (float) track->get_p(),       //total momentum of track
0843        * (float) track->get_charge(),  //electric charge of track
0844        * (float) track->get_quality()  //track quality */
0845
0846       //Loop over tracks in event //
0847     for (SvtxTrackMap::ConstIter track_itr = trackmap->begin();
0848            track_itr != trackmap->end(); track_itr++) {
0849
0850         SvtxTrack* track = dynamic_cast<SvtxTrack*>(track_itr->second);
0851
0852     // Get track variables //
0853         float track_eta = track->get_eta();
0854         float track_phi = track->get_phi();
0855         int track_charge = track->get_charge();
0856     double gvx,gvy,gvz;
0857
0858         float delta_R = CalculateDeltaR( track_eta, track_phi, jet_eta, jet_phi );
0859
0860     PHG4Particle* g4particle = trackeval->max_truth_particle_by_nclusters(track);
0861
0862     // Get true vertex distances //
0863     PHG4VtxPoint* vtx = trutheval->get_vertex(g4particle);
0864     gvx = vtx->get_x();
0865     gvy = vtx->get_y();
0866     gvz = vtx->get_z();
0867
0868
0869
0870     // If charged track is within jet then use its vertex //
0871     if(delta_R < 0.5 && trutheval->is_primary(g4particle)) tracks_vertex.push_back(sqrt(pow(gvx,2)+pow(gvy,2)+pow(gvz,2)));
0872
0873     /* if save_tracks set true: add track to tree */
0874         if ( _save_tracks )
0875           {
0876             float track_data[17] = {(float) _ievent,
0877                                     (float) (iter->second)->get_property_uint( PidCandidate::jet_id ),
0878                                     (float) (iter->second)->get_property_int( PidCandidate::evtgen_pid ),
0879                                     (float) (iter->second)->get_property_float( PidCandidate::evtgen_etotal ),
0880                                     (float) (iter->second)->get_property_float( PidCandidate::evtgen_eta ),
0881                                     (float) (iter->second)->get_property_float( PidCandidate::evtgen_phi ),
0882                                     (float) (iter->second)->get_property_uint( PidCandidate::evtgen_decay_prong ),
0883                                     (float) (iter->second)->get_property_uint( PidCandidate::evtgen_decay_hcharged ),
0884                                     (float) (iter->second)->get_property_uint( PidCandidate::evtgen_decay_lcharged ),
0885                                     (float) (iter->second)->get_property_float( PidCandidate::jet_eta ),
0886                                     (float) (iter->second)->get_property_float( PidCandidate::jet_phi ),
0887                                     (float) (iter->second)->get_property_float( PidCandidate::jet_etotal ),
0888                                     (float) track->get_quality(),
0889                                     (float) track_eta,
0890                                     (float) track_phi,
0891                                     (float) delta_R,
0892                                     (float) track->get_p()
0893             };
0894
0895             _ntp_track->Fill(track_data);
0896           }
0897
0898         /* If track within search cone, update track information for tau candidate */
0899         if ( delta_R < 0.2 )
0900           {
0901             tracks_count_r02++;
0902             tracks_chargesum_r02 += track_charge;
0903
0904             if ( delta_R > tracks_rmax_r02 )
0905               tracks_rmax_r02 = delta_R;
0906           }
0907
0908         if ( delta_R < 0.4 )
0909           {
0910             tracks_count_r04++;
0911             tracks_chargesum_r04 += track_charge;
0912
0913             if ( delta_R > tracks_rmax_r04 )
0914               tracks_rmax_r04 = delta_R;
0915           }
0916
0917     if ( delta_R < R_max )
0918           {
0919             tracks_count_R++;
0920             tracks_chargesum_R += track_charge;
0921
0922             if ( delta_R > tracks_rmax_R )
0923               tracks_rmax_R = delta_R;
0924           }
0925
0926     } // end loop over reco tracks //
0927     // sort vertex array in increasing order //
0928     std::sort(tracks_vertex.begin(),tracks_vertex.end());
0929
0930     // Compute average vertex distance of tracks in jet //
0931     float avg = Average(tracks_vertex);
0932
0933     /* Set track-based properties for tau candidate */
0934       (iter->second)->set_property( PidCandidate::tracks_count_r02, tracks_count_r02 );
0935       (iter->second)->set_property( PidCandidate::tracks_chargesum_r02, tracks_chargesum_r02 );
0936       (iter->second)->set_property( PidCandidate::tracks_rmax_r02, tracks_rmax_r02 );
0937       (iter->second)->set_property( PidCandidate::tracks_count_r04, tracks_count_r04 );
0938       (iter->second)->set_property( PidCandidate::tracks_chargesum_r04, tracks_chargesum_r04 );
0939       (iter->second)->set_property( PidCandidate::tracks_rmax_r04, tracks_rmax_r04 );
0940       (iter->second)->set_property( PidCandidate::tracks_count_R, tracks_count_R );
0941       (iter->second)->set_property( PidCandidate::tracks_chargesum_R, tracks_chargesum_R );
0942       (iter->second)->set_property( PidCandidate::tracks_rmax_R, tracks_rmax_R );
0943       if(avg == avg) (iter->second)->set_property( PidCandidate::tracks_vertex, avg);
0944
0945     } // end loop over tau  candidates
0946
0947   return 0;
0948 }
0949
0950
0951 int
0952 LeptoquarksReco::WritePidCandidatesToTree( type_map_tcan& tauCandidateMap )
0953 {
0954   /* Loop over all tau candidates and add them to tree*/
0955   for (type_map_tcan::iterator iter = tauCandidateMap.begin();
0956        iter != tauCandidateMap.end();
0957        ++iter)
0958     {
0959       /* update information in map and fill tree */
0960       for ( map< PidCandidate::PROPERTY , vector<float> >::iterator iter_prop = _map_tau_candidate_branches.begin();
0961             iter_prop != _map_tau_candidate_branches.end();
0962             ++iter_prop)
0963         {
0964           switch ( PidCandidate::get_property_info( (iter_prop->first) ).second ) {
0965
0966           case PidCandidate::type_float :
0967             (iter_prop->second).push_back( (iter->second)->get_property_float( (iter_prop->first) ) );
0968             break;
0969
0970           case PidCandidate::type_int :
0971             (iter_prop->second).push_back( (iter->second)->get_property_int( (iter_prop->first) ) );
0972             break;
0973
0974           case PidCandidate::type_uint :
0975             (iter_prop->second).push_back( (iter->second)->get_property_uint( (iter_prop->first) ) );
0976             break;
0977
0978           case PidCandidate::type_unknown :
0979             break;
0980           }
0981         }
0982     }
0983
0984   return 0;
0985 }
0986
0987
0988 PidCandidate*
0989 LeptoquarksReco::FindMinDeltaRCandidate( type_map_tcan *candidates, const float eta_ref, const float phi_ref )
0990 {
0991   /* PidCandidate with eta, phi closest to reference */
0992   PidCandidate* best_candidate = NULL;
0993
0994   float eta_ref_local = eta_ref;
0995   float phi_ref_local = phi_ref;
0996
0997   /* Event record uses 0 < phi < 2Pi, while Fun4All uses -Pi < phi < Pi.
0998    * Therefore, correct angle for 'wraparound' at phi == Pi */
0999   if( phi_ref_local > TMath::Pi() ) phi_ref_local = phi_ref_local - 2*TMath::Pi();
1000
1001   /* find jet with smallest delta_R from quark */
1002   float min_delta_R = 100;
1003   type_map_tcan::iterator min_delta_R_iter = candidates->end();
1004
1005   for (type_map_tcan::iterator iter = candidates->begin();
1006        iter != candidates->end();
1007        ++iter)
1008     {
1009       float eta = (iter->second)->get_property_float( PidCandidate::jet_eta );
1010       float phi = (iter->second)->get_property_float( PidCandidate::jet_phi );
1011
1012       float delta_R = CalculateDeltaR( eta, phi, eta_ref_local, phi_ref_local );
1013       //cout<<delta_R<<endl;
1014       if ( delta_R < min_delta_R )
1015         {
1016           min_delta_R_iter = iter;            ;
1017           min_delta_R = delta_R;
1018         }
1019     }
1020
1021   /* set best_candidate to PidCandiate with smallest delta_R within reasonable range*/
1022   if ( min_delta_R_iter != candidates->end() && min_delta_R < 0.5 )
1023     best_candidate = min_delta_R_iter->second;
1024   //cout<<"Min R: "<<min_delta_R<<endl;
1025   return best_candidate;
1026 }
1027
1028
1029 float
1030 LeptoquarksReco::CalculateDeltaR( float eta1, float phi1, float eta2, float phi2 )
1031 {
1032   /* Particles at phi = PI+x and phi = PI-x are actually close to each other in phi, but simply calculating
1033    * the difference in phi would give a large distance (because phi ranges from -PI to +PI in the convention
1034    * used. Account for this by subtracting 2PI is particles fall within this border area. */
1035   if( ( phi1 < -0.9*TMath::Pi()) && ( phi2 >  0.9*TMath::Pi() ) ) phi2 = phi2 - 2*TMath::Pi();
1036   if( ( phi1 >  0.9*TMath::Pi()) && ( phi2 < -0.9*TMath::Pi() ) ) phi1 = phi1 - 2*TMath::Pi();
1037
1038   float delta_R = sqrt( pow( eta2 - eta1, 2 ) + pow( phi2 - phi1, 2 ) );
1039
1040   return delta_R;
1041 }
1042
1043
1044 int
1045 LeptoquarksReco::AddGlobalEventInformation( type_map_tcan& tauCandidateMap, type_map_cdata* map_towers )
1046 {
1047   /* missing transverse momentum */
1048   //float pt_miss = 0;
1049
1050   /* direction of missing transverse momentum */
1051   //float pt_miss_phi = 0;
1052
1053   /* Missing transverse energy, transverse energy direction, and Energy sums, x and y components separately */
1054   float Et_miss = 0;
1055   float Et_miss_phi = 0;
1056   float Ex_sum = 0;
1057   float Ey_sum = 0;
1058
1059   /* energy threshold for considering tower */
1060   float tower_emin = 0.0;
1061
1062   /* Loop over all tower (and geometry) collections */
1063   for (type_map_cdata::iterator iter_calo = map_towers->begin();
1064        iter_calo != map_towers->end();
1065        ++iter_calo)
1066     {
1067
1068       /* define tower iterator */
1069       RawTowerContainer::ConstRange begin_end = ((iter_calo->second).first)->getTowers();
1070       RawTowerContainer::ConstIterator rtiter;
1071
1072       /* loop over all tower in CEMC calorimeter */
1073       for (rtiter = begin_end.first; rtiter !=  begin_end.second; ++rtiter)
1074         {
1075           /* get tower energy */
1076           RawTower *tower = rtiter->second;
1077           float tower_energy = tower->get_energy();
1078
1079           /* check if tower above energy treshold */
1080           if ( tower_energy < tower_emin )
1081       continue;
1082
1083           /* get eta and phi of tower and check angle delta_R w.r.t. jet axis */
1084           RawTowerGeom * tower_geom = ((iter_calo->second).second)->get_tower_geometry(tower -> get_key());
1085           float tower_eta = tower_geom->get_eta();
1086           float tower_phi = tower_geom->get_phi();
1087
1088           /* from https://en.wikipedia.org/wiki/Pseudorapidity:
1089              p_x = p_T * cos( phi )
1090              p_y = p_T * sin( phi )
1091              p_z = p_T * sinh( eta )
1092              |p| = p_T * cosh( eta )
1093           */
1094
1095           /* calculate 'transverse' tower energy */
1096           float tower_energy_t = tower_energy / cosh( tower_eta );
1097
1098           /* add energy components of this tower to total energy components */
1099           Ex_sum += tower_energy_t * cos( tower_phi );
1100           Ey_sum += tower_energy_t * sin( tower_phi );
1101         }
1102       }
1103
1104   /* calculate Et_miss and phi angle*/
1105   Et_miss = sqrt( Ex_sum * Ex_sum + Ey_sum * Ey_sum );
1106   Et_miss_phi = atan2( Ey_sum , Ex_sum );
1107
1108   /* Loop over tau candidates and find tau jet*/
1109   PidCandidate* the_tau = NULL;
1110   for (type_map_tcan::iterator iter = tauCandidateMap.begin();
1111        iter != tauCandidateMap.end();
1112        ++iter)
1113     {
1114       if ( ( iter->second)->get_property_int( PidCandidate::evtgen_pid ) == 15 )
1115         the_tau = iter->second;
1116     }
1117
1118   /* update event information tree variables */
1119   /* @TODO make this better protected against errors- if 'find' returns NULL pointer,
1120      this will lead to a SEGMENTATION FAULT */
1121   ( _map_event_branches.find( "Et_miss" ) )->second = Et_miss;
1122   ( _map_event_branches.find( "Et_miss_phi" ) )->second = Et_miss_phi;
1123
1124   // If tau is found, write variables //
1125   if ( the_tau )
1126     {
1127       ( _map_event_branches.find( "reco_tau_found" ) )->second = 1;
1128       ( _map_event_branches.find( "reco_tau_is_tau" ) )->second =
1129         the_tau->get_property_int( PidCandidate::evtgen_pid );
1130       ( _map_event_branches.find( "reco_tau_eta" ) )->second =
1131         the_tau->get_property_float( PidCandidate::jet_eta );
1132       ( _map_event_branches.find( "reco_tau_phi" ) )->second =
1133         the_tau->get_property_float( PidCandidate::jet_phi );
1134       ( _map_event_branches.find( "reco_tau_ptotal" ) )->second =
1135         the_tau->get_property_float( PidCandidate::jet_ptotal );
1136     }
1137   else
1138     {
1139       ( _map_event_branches.find( "reco_tau_found" ) )->second = 0;
1140       ( _map_event_branches.find( "reco_tau_is_tau" ) )->second = NAN;
1141       ( _map_event_branches.find( "reco_tau_eta" ) )->second = NAN;
1142       ( _map_event_branches.find( "reco_tau_phi" ) )->second = NAN;
1143       ( _map_event_branches.find( "reco_tau_ptotal" ) )->second = NAN;
1144     }
1145
1146   return 0;
1147 }
1148
1149
1150 void
1151 LeptoquarksReco::ResetBranchMap()
1152 {
1153   /* Event branches */
1154   for ( map< string , float >::iterator iter = _map_event_branches.begin();
1155         iter != _map_event_branches.end();
1156         ++iter)
1157     {
1158       (iter->second) = NAN;
1159     }
1160
1161   /* Tau candidate branches */
1162   for ( map< PidCandidate::PROPERTY , vector<float> >::iterator iter = _map_tau_candidate_branches.begin();
1163         iter != _map_tau_candidate_branches.end();
1164         ++iter)
1165     {
1166       (iter->second).clear();
1167     }
1168
1169   return;
1170 }
1171
1172
1173 int
1174 LeptoquarksReco::End(PHCompositeNode *topNode)
1175 {
1176   _tfile->cd();
1177
1178   if ( _t_event )
1179     _t_event->Write();
1180
1181   if ( _ntp_tower )
1182     _ntp_tower->Write();
1183
1184   if ( _ntp_track )
1185     _ntp_track->Write();
1186
1187   _tfile->Close();
1188
1189   return 0;
1190 }