Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
LeptoquarksReco.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file LeptoquarksReco.C
1 #include "LeptoquarksReco.h"
2 #include "PidCandidatev1.h"
3 #include "TruthTrackerHepMC.h"
4 
5 /* STL includes */
6 #include <cassert>
7 
8 /* Fun4All includes */
9 #include <trackbase_historic/SvtxTrackMap_v1.h>
10 #include <trackbase_historic/SvtxVertexMap_v1.h>
11 //#include <trackbase_historic/SvtxClusterMap_v1.h>
12 
13 #include <phool/getClass.h>
14 
15 #include <fun4all/Fun4AllReturnCodes.h>
16 
17 #include <g4jets/JetMap.h>
18 
19 #include <calobase/RawTowerGeomContainer.h>
20 #include <calobase/RawTowerContainer.h>
21 #include <calobase/RawTowerGeom.h>
22 #include <calobase/RawTowerv1.h>
23 
24 #include <g4vertex/GlobalVertexMap.h>
25 #include <g4vertex/GlobalVertex.h>
26 
27 #include <g4main/PHG4VtxPoint.h>
28 #include <g4main/PHG4Particle.h>
29 
30 #include <g4eval/CaloRawTowerEval.h>
31 #include <g4eval/SvtxEvalStack.h>
32 
33 
34 /* ROOT includes */
35 #include <TLorentzVector.h>
36 #include <TString.h>
37 #include <TNtuple.h>
38 #include <TTree.h>
39 #include <TFile.h>
40 
41 using namespace std;
42 
43 LeptoquarksReco::LeptoquarksReco(std::string filename) :
44  SubsysReco("LeptoquarksReco" ),
45  _save_towers(false),
46  _save_tracks(false),
47  _ievent(0),
48  _filename(filename),
49  _tfile(nullptr),
50  _t_event(nullptr),
51  _ntp_tower(nullptr),
52  _ntp_track(nullptr),
53  _ebeam_E(0),
54  _pbeam_E(0),
55  _tau_jet_emin(5.0),
56  _jetcolname("AntiKt_Tower_r05")
57 {
58 
59 }
60 
61 int
62 LeptoquarksReco::Init(PHCompositeNode *topNode)
63 {
64  _ievent = 0;
65 
66  _tfile = new TFile(_filename.c_str(), "RECREATE");
67 
68  /* Add PidCandidate properties to map that defines output tree */
69  float dummy = 0;
70  vector< float > vdummy;
71 
72  _map_tau_candidate_branches.insert( make_pair( PidCandidate::jet_id , vdummy ) );
73  _map_tau_candidate_branches.insert( make_pair( PidCandidate::evtgen_pid , vdummy ) );
74  _map_tau_candidate_branches.insert( make_pair( PidCandidate::evtgen_etotal , vdummy ) );
75  _map_tau_candidate_branches.insert( make_pair( PidCandidate::evtgen_ptotal , vdummy ) );
76  _map_tau_candidate_branches.insert( make_pair( PidCandidate::evtgen_theta , vdummy ) );
77  _map_tau_candidate_branches.insert( make_pair( PidCandidate::evtgen_eta , vdummy ) );
78  _map_tau_candidate_branches.insert( make_pair( PidCandidate::evtgen_phi , vdummy ) );
79  _map_tau_candidate_branches.insert( make_pair( PidCandidate::evtgen_decay_prong , vdummy ) );
80  _map_tau_candidate_branches.insert( make_pair( PidCandidate::evtgen_decay_hcharged , vdummy ) );
81  _map_tau_candidate_branches.insert( make_pair( PidCandidate::evtgen_decay_lcharged , vdummy ) );
82  _map_tau_candidate_branches.insert( make_pair( PidCandidate::jet_eta , vdummy ) );
83  _map_tau_candidate_branches.insert( make_pair( PidCandidate::jet_phi , vdummy ) );
84  _map_tau_candidate_branches.insert( make_pair( PidCandidate::jet_etotal , vdummy ) );
85  _map_tau_candidate_branches.insert( make_pair( PidCandidate::jet_etrans , vdummy ) );
86  _map_tau_candidate_branches.insert( make_pair( PidCandidate::jet_ptotal , vdummy ) );
87  _map_tau_candidate_branches.insert( make_pair( PidCandidate::jet_ptrans , vdummy ) );
88  _map_tau_candidate_branches.insert( make_pair( PidCandidate::jet_minv , vdummy ) );
89  _map_tau_candidate_branches.insert( make_pair( PidCandidate::jet_mtrans , vdummy ) );
90  _map_tau_candidate_branches.insert( make_pair( PidCandidate::jet_ncomp , vdummy ) );
91  _map_tau_candidate_branches.insert( make_pair( PidCandidate::jet_ncomp_above_0p1 , vdummy ) );
92  _map_tau_candidate_branches.insert( make_pair( PidCandidate::jet_ncomp_above_1 , vdummy ) );
93  _map_tau_candidate_branches.insert( make_pair( PidCandidate::jet_ncomp_above_10 , vdummy ) );
94  _map_tau_candidate_branches.insert( make_pair( PidCandidate::jet_ncomp_emcal , vdummy ) );
95  _map_tau_candidate_branches.insert( make_pair( PidCandidate::jetshape_econe_r01 , vdummy ) );
96  _map_tau_candidate_branches.insert( make_pair( PidCandidate::jetshape_econe_r02 , vdummy ) );
97  _map_tau_candidate_branches.insert( make_pair( PidCandidate::jetshape_econe_r03 , vdummy ) );
98  _map_tau_candidate_branches.insert( make_pair( PidCandidate::jetshape_econe_r04 , vdummy ) );
99  _map_tau_candidate_branches.insert( make_pair( PidCandidate::jetshape_econe_r05 , vdummy ) );
100  _map_tau_candidate_branches.insert( make_pair( PidCandidate::jetshape_r90 , vdummy ) );
101  _map_tau_candidate_branches.insert( make_pair( PidCandidate::jetshape_rms , vdummy ) );
102  _map_tau_candidate_branches.insert( make_pair( PidCandidate::jetshape_radius , vdummy ) );
103  _map_tau_candidate_branches.insert( make_pair( PidCandidate::jetshape_emcal_econe_r01 , vdummy ) );
104  _map_tau_candidate_branches.insert( make_pair( PidCandidate::jetshape_emcal_econe_r02 , vdummy ) );
105  _map_tau_candidate_branches.insert( make_pair( PidCandidate::jetshape_emcal_econe_r03 , vdummy ) );
106  _map_tau_candidate_branches.insert( make_pair( PidCandidate::jetshape_emcal_econe_r04 , vdummy ) );
107  _map_tau_candidate_branches.insert( make_pair( PidCandidate::jetshape_emcal_econe_r05 , vdummy ) );
108  _map_tau_candidate_branches.insert( make_pair( PidCandidate::jetshape_emcal_r90 , vdummy ) );
109  _map_tau_candidate_branches.insert( make_pair( PidCandidate::jetshape_emcal_rms , vdummy ) );
110  _map_tau_candidate_branches.insert( make_pair( PidCandidate::jetshape_emcal_radius , vdummy ) );
111  _map_tau_candidate_branches.insert( make_pair( PidCandidate::tracks_count_r02 , vdummy ) );
112  _map_tau_candidate_branches.insert( make_pair( PidCandidate::tracks_chargesum_r02 , vdummy ) );
113  _map_tau_candidate_branches.insert( make_pair( PidCandidate::tracks_rmax_r02 , vdummy ) );
114  _map_tau_candidate_branches.insert( make_pair( PidCandidate::tracks_count_r04 , vdummy ) );
115  _map_tau_candidate_branches.insert( make_pair( PidCandidate::tracks_chargesum_r04 , vdummy ) );
116  _map_tau_candidate_branches.insert( make_pair( PidCandidate::tracks_rmax_r04 , vdummy ) );
117  _map_tau_candidate_branches.insert( make_pair( PidCandidate::tracks_count_R , vdummy ) );
118  _map_tau_candidate_branches.insert( make_pair( PidCandidate::tracks_chargesum_R , vdummy ) );
119  _map_tau_candidate_branches.insert( make_pair( PidCandidate::tracks_rmax_R , vdummy ) );
120  _map_tau_candidate_branches.insert( make_pair( PidCandidate::tracks_vertex , vdummy ) );
121 
122  /* Add branches to map that defines output tree for event-wise properties */
123  _map_event_branches.insert( make_pair( "Et_miss" , dummy ) );
124  _map_event_branches.insert( make_pair( "Et_miss_phi" , dummy ) );
125  _map_event_branches.insert( make_pair( "reco_tau_found" , dummy ) );
126  _map_event_branches.insert( make_pair( "reco_tau_is_tau" , dummy ) );
127  _map_event_branches.insert( make_pair( "reco_tau_eta" , dummy ) );
128  _map_event_branches.insert( make_pair( "reco_tau_phi" , dummy ) );
129  _map_event_branches.insert( make_pair( "reco_tau_ptotal" , dummy ) );
130  _map_event_branches.insert( make_pair( "is_lq_event" , dummy ) );
131 
132 
133  /* Create tree for information about full event */
134  _t_event = new TTree("event", "a Tree with global event information and tau candidates");
135  _t_event->Branch("event", &_ievent, "event/I");
136 
137  /* Add event branches */
138  for ( map< string , float >::iterator iter = _map_event_branches.begin();
139  iter != _map_event_branches.end();
140  ++iter)
141  {
142  _t_event->Branch( (iter->first).c_str(),
143  &(iter->second) );
144  }
145 
146  /* Add tau candidate branches */
147  for ( map< PidCandidate::PROPERTY , vector< float > >::iterator iter = _map_tau_candidate_branches.begin();
148  iter != _map_tau_candidate_branches.end();
149  ++iter)
150  {
151  _t_event->Branch(PidCandidate::get_property_info( (iter->first) ).first.c_str(),
152  &(iter->second) );
153  }
154 
155 
156  /* Create additional tree for towers */
157  if ( _save_towers )
158  {
159  _ntp_tower = new TNtuple("ntp_tower","towers from all all tau candidates",
160  "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");
161  }
162 
163  /* Create additional tree for tracks */
164  if ( _save_tracks )
165  {
166  _ntp_track = new TNtuple("ntp_track","tracks from all all tau candidates",
167  "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");
168  }
169 
170  return 0;
171 }
172 
173 int
174 LeptoquarksReco::process_event(PHCompositeNode *topNode)
175 {
176  /* Reset branch map */
177  ResetBranchMap();
178 
179  /* Create map to collect tau candidates.
180  * Use energy as 'key' to the map because energy is unique for each jet, while there are sometimes multiple jets (same energy,
181  * different jet ID) in the input jet collection. Also, map automatically sorts entries by key, i.e. this gives list of tau candidates
182  * sorted by energy. */
183  type_map_tcan tauCandidateMap;
184 
185  /* Get reco jets collection */
186  JetMap* recojets = findNode::getClass<JetMap>(topNode,_jetcolname.c_str());
187  if (!recojets)
188  {
189  cerr << PHWHERE << " ERROR: Can't find " << _jetcolname << endl;
190  return Fun4AllReturnCodes::ABORTEVENT;
191  }
192 
193  /* Get track collection with all tracks in this event */
194  SvtxTrackMap* trackmap =
195  findNode::getClass<SvtxTrackMap>(topNode,"SvtxTrackMap");
196  if (!trackmap)
197  {
198  cout<< PHWHERE << "SvtxTrackMap node not found on node tree"
199  << endl;
200  return Fun4AllReturnCodes::ABORTEVENT;
201  }
202 
203  /* Get vertex collection with all vertices in this event */
204  SvtxVertexMap* vertexmap =
205  findNode::getClass<SvtxVertexMap>(topNode,"SvtxVertexMap");
206  if (!vertexmap)
207  {
208  cout<< PHWHERE << "SvtxVertexMap node not found on node tree"
209  << endl;
210  return Fun4AllReturnCodes::ABORTEVENT;
211  }
212 
213 
214  /* Get collection of truth particles from event generator */
215  PHHepMCGenEventMap *genevtmap = findNode::getClass<PHHepMCGenEventMap>(topNode,"PHHepMCGenEventMap");
216  if (!genevtmap) {
217  cerr << PHWHERE << " ERROR: Can't find PHHepMCGenEventMap" << endl;
218  return Fun4AllReturnCodes::ABORTEVENT;
219  }
220 
221  /* create map of all tower collections and tower geometry collections used */
222  type_map_cdata map_calotower;
223 
224  /* Select calorimeter to include */
225  vector< RawTowerDefs::CalorimeterId > v_caloids;
226  v_caloids.push_back( RawTowerDefs::CEMC );
227  v_caloids.push_back( RawTowerDefs::HCALIN );
228  v_caloids.push_back( RawTowerDefs::HCALOUT );
229  v_caloids.push_back( RawTowerDefs::FEMC );
230  v_caloids.push_back( RawTowerDefs::FHCAL );
231  v_caloids.push_back( RawTowerDefs::EEMC );
232 
233  /* Fill map with calorimeter data */
234  for ( unsigned i = 0; i < v_caloids.size(); i++ )
235  {
236  /* Get collection of towers */
237  string towers_name = "TOWER_CALIB_";
238  towers_name.append( RawTowerDefs::convert_caloid_to_name( v_caloids.at( i ) ) );
239 
240  RawTowerContainer *towers = findNode::getClass<RawTowerContainer>(topNode,towers_name);
241  if (!towers) {
242  cerr << PHWHERE << " ERROR: Can't find " << towers_name << endl;
243  return Fun4AllReturnCodes::ABORTEVENT;
244  }
245 
246  /* get geometry collection */
247  string towergeom_name = "TOWERGEOM_";
248  towergeom_name.append( RawTowerDefs::convert_caloid_to_name( v_caloids.at( i ) ) );
249 
250  RawTowerGeomContainer *geom = findNode::getClass<RawTowerGeomContainer>(topNode,towergeom_name);
251  if (!geom) {
252  cerr << PHWHERE << " ERROR: Can't find " << towergeom_name << endl;
253  return Fun4AllReturnCodes::ABORTEVENT;
254  }
255 
256  /* Insert tower and geometry collections in map */
257  map_calotower.insert( make_pair( v_caloids.at( i ), make_pair( towers, geom ) ) );
258  }
259 
260 
261  /* Loop over all jets to collect tau candidate objects */
262  for (JetMap::Iter iter = recojets->begin();
263  iter != recojets->end();
264  ++iter)
265  {
266  /* skip jets below energy threshold */
267  if ( (iter->second)->get_e() < _tau_jet_emin )
268  continue;
269 
270  /* create new tau candidate */
271  PidCandidatev1 *tc = new PidCandidatev1();
272  tc->set_candidate_id( (iter->second)->get_id() );
273  tc->set_property( PidCandidate::jet_id , (iter->second)->get_id() );
274 
275  tc->set_property( PidCandidate::evtgen_pid, (int)0 );
276 
277  /* add tau candidate to collection */
278  tauCandidateMap.insert( make_pair( (iter->second)->get_e(), tc ) );
279  }
280 
281  /* vertex information */
282  SvtxEvalStack *svtxevalstack = new SvtxEvalStack(topNode);
283 
284  /* Add jet information to tau candidates */
285  AddJetInformation( tauCandidateMap, recojets, &map_calotower );
286 
287  /* Add jet structure information to tau candidates */
288  AddJetStructureInformation( tauCandidateMap, &map_calotower );
289 
290  /* Add track information to tau candidates */
291  AddTrackInformation( tauCandidateMap, trackmap, vertexmap, svtxevalstack, recojets->get_par());
292 
293  /* Add tag for true Tau particle jet to tau candidates */
294  AddTrueTauTag( tauCandidateMap, genevtmap );
295 
296  /* Add information about tau candidats to output tree */
297  WritePidCandidatesToTree( tauCandidateMap );
298 
299  /* Add global event information to separate tree */
300  AddGlobalEventInformation( tauCandidateMap, &map_calotower );
301 
302  /* fill event information tree */
303  _t_event->Fill();
304 
305  /* count up event number */
306  _ievent ++;
307 
308  return 0;
309 }
310 
311 // Returns average of values while excluding outliers //
312 float Average(vector<float> v)
313 {
314  double sum = 0;
315  int med_i = v.size()/2;
316  vector<float> temp_v;
317  for(int i=0;(unsigned)i<v.size();i++)
318  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]);
319 
320  for(int j=0;(unsigned)j<temp_v.size();j++) sum = sum + temp_v[j];
321 
322  return sum/temp_v.size();
323 }
324 
325 int
326 LeptoquarksReco::AddTrueTauTag( type_map_tcan& tauCandidateMap, PHHepMCGenEventMap *genevtmap )
327 {
328  /* Look for leptoquark and tau particle in the event */
329  TruthTrackerHepMC truth;
330  truth.set_hepmc_geneventmap( genevtmap );
331 
332  int pdg_lq = 39; // leptoquark
333  int pdg_tau = 15; // tau lepton
334  int pdg_parton = 0;
335  int pdg_electron = 11; // electron
336 
337  // Check if electron exists //
338  HepMC::GenParticle* particle_electron = NULL;
339 
340  for (PHHepMCGenEventMap::ReverseIter iter = genevtmap->rbegin(); iter != genevtmap->rend(); ++iter)
341  {
342  PHHepMCGenEvent *genevt = iter->second;
343  HepMC::GenEvent *theEvent = genevt->getEvent();
344 
345  // Loop over all truth particles in event generator collection
346  for ( HepMC::GenEvent::particle_iterator p = theEvent->particles_begin();
347  p != theEvent->particles_end(); ++p ) {
348  // If final state electron then tag it //
349  if((*p)->pdg_id() == pdg_electron && (*p)->status() == 1) particle_electron = (*p);
350  }
351  }
352 
353  /* Search for leptoquark in event */
354  HepMC::GenParticle* particle_lq = truth.FindParticle( pdg_lq );
355 
356  // Tag event as LQ event or not //
357  if (particle_lq) ( _map_event_branches.find( "is_lq_event" ) )->second = 1;
358  else ( _map_event_branches.find( "is_lq_event" ) )->second = 0;
359 
360  /* Search for lq->tau decay in event */
361  HepMC::GenParticle* particle_tau = NULL;
362 
363  // Find tau that comes from LQ or just any tau in event //
364  if(particle_lq) particle_tau = truth.FindDaughterParticle( pdg_tau, particle_lq );
365  else particle_tau = truth.FindParticle(pdg_tau);
366 
367  /* Search for lq->quark decay in event.
368  * Loop over all quark PDG codes until finding a matching quark. */
369  HepMC::GenParticle* particle_quark = NULL;
370  for ( int pdg_quark = 1; pdg_quark < 7; pdg_quark++ )
371  {
372  /* try quark */
373  particle_quark = truth.FindDaughterParticle( pdg_quark, particle_lq );
374  if (particle_quark)
375  {
376  pdg_parton = pdg_quark;
377  break;
378  }
379 
380  /* try anti-quark */
381  particle_quark = truth.FindDaughterParticle( -pdg_quark, particle_lq );
382  if (particle_quark)
383  {
384  pdg_parton = -pdg_quark;
385  break;
386  }
387  }
388 
389  /* If TAU in event: Tag the tau candidate (i.e. jet) with smalles delta_R from this tau */
390  if( particle_tau )
391  {
392  PidCandidate* best_match = FindMinDeltaRCandidate( &tauCandidateMap,
393  particle_tau->momentum().eta(),
394  particle_tau->momentum().phi() );
395 
396 
397  /* set is_tau = TRUE for PidCandiate with smallest delta_R within reasonable range*/
398  if ( best_match )
399  {
400  /* Check: If PidCandidate::Evtgen_pid has already been set to a value != 0, exit function. */
401  if( best_match->get_property_int( PidCandidate::evtgen_pid ) != 0 )
402  {
403  cout << "ERROR: Try to set PidCandidate::evtgen_pid for PidCandidate with evtgen_pid != 0" << endl;
404  return -1;
405  }
406 
407  /* Update PidCandidate entry */
408  best_match->set_property( PidCandidate::evtgen_pid, pdg_tau );
409  best_match->set_property( PidCandidate::evtgen_etotal, (float)particle_tau->momentum().e() );
410  best_match->set_property( PidCandidate::evtgen_eta, (float)particle_tau->momentum().eta() );
411  best_match->set_property( PidCandidate::evtgen_phi, (float)particle_tau->momentum().phi() );
412 
413  /* Check particle decay if end-vertex found */
414  if ( particle_tau->end_vertex() )
415  {
416  /* Add information about tau decay */
417  uint tau_decay_prong = 0;
418  uint tau_decay_hcharged = 0;
419  uint tau_decay_lcharged = 0;
420 
421  /* Count how many charged particles (hadrons and leptons) a given particle decays into. */
422  truth.FindDecayParticles( particle_tau, tau_decay_prong, tau_decay_hcharged, tau_decay_lcharged );
423 
424 
425  /* Update tau candidate entry */
426  best_match->set_property( PidCandidate::evtgen_decay_prong, tau_decay_prong );
427  best_match->set_property( PidCandidate::evtgen_decay_hcharged, tau_decay_hcharged );
428  best_match->set_property( PidCandidate::evtgen_decay_lcharged, tau_decay_lcharged );
429  }
430  }
431  }
432 
433  /* If QUARK (->jet) in event: Tag the tau candidate (i.e. jet) with smalles delta_R from this quark */
434  if( particle_quark )
435  {
436  PidCandidate* best_match = FindMinDeltaRCandidate( &tauCandidateMap,
437  particle_quark->momentum().eta(),
438  particle_quark->momentum().phi() );
439 
440  /* set is_uds = TRUE for PidCandiate with smallest delta_R if found */
441  if ( best_match )
442  {
443  /* Check: If PidCandidate::Evtgen_pid has already been set to a value != 0, exit function. */
444  if( best_match->get_property_int( PidCandidate::evtgen_pid ) != 0 )
445  {
446  cout << "ERROR: Try to set PidCandidate::evtgen_pid for PidCandidate with evtgen_pid != 0" << endl;
447  return -1;
448  }
449  /* Set properties of PidCandidate */
450  best_match->set_property( PidCandidate::evtgen_pid, pdg_parton );
451  best_match->set_property( PidCandidate::evtgen_etotal, (float)particle_quark->momentum().e() );
452  best_match->set_property( PidCandidate::evtgen_eta, (float)particle_quark->momentum().eta() );
453  best_match->set_property( PidCandidate::evtgen_phi, (float)particle_quark->momentum().phi() );
454  }
455  }
456 
457  if( particle_electron )
458  {
459  //cout<<"ELECTRON"<<endl;
460  PidCandidate* best_match = FindMinDeltaRCandidate( &tauCandidateMap,
461  particle_electron->momentum().eta(),
462  particle_electron->momentum().phi() );
463 
464  /* set electron = TRUE for PidCandiate with smallest delta_R if found */
465  if ( best_match )
466  {
467  /* Check: If PidCandidate::Evtgen_pid has already been set to a value != 0, exit function. */
468  if( best_match->get_property_int( PidCandidate::evtgen_pid ) != 0 )
469  {
470  cout << "ERROR: Try to set PidCandidate::evtgen_pid for PidCandidate with evtgen_pid != 0" << endl;
471  return -1;
472  }
473 
474  /* Set properties of PidCandidate */
475  best_match->set_property( PidCandidate::evtgen_pid, pdg_electron );
476  best_match->set_property( PidCandidate::evtgen_etotal, (float)particle_electron->momentum().e() );
477  best_match->set_property( PidCandidate::evtgen_eta, (float)particle_electron->momentum().eta() );
478  best_match->set_property( PidCandidate::evtgen_phi, (float)particle_electron->momentum().phi() );
479  }
480  }
481 
482 
483  return 0;
484 }
485 
486 
487 int
488 LeptoquarksReco::AddJetInformation( type_map_tcan& tauCandidateMap, JetMap* recojets, type_map_cdata* map_calotower )
489 {
490  /* Loop over tau candidates */
491  for (type_map_tcan::iterator iter = tauCandidateMap.begin();
492  iter != tauCandidateMap.end();
493  ++iter)
494  {
495  Jet* jetx = recojets->get( (iter->second)->get_property_uint( PidCandidate::jet_id ) );
496 
497  /* calculate transverse mass of jet */
498  float jet_mtrans = sqrt( pow( jetx->get_mass(), 2 ) +
499  pow( jetx->get_pt(), 2 ) );
500 
501  /* count jet ncomp above thresholds */
502  unsigned int jet_ncomp_above_0p1 = 0;
503  unsigned int jet_ncomp_above_1 = 0;
504  unsigned int jet_ncomp_above_10 = 0;
505 
506  for (Jet::ConstIter jcompiter = jetx->begin_comp(); jcompiter != jetx->end_comp(); ++jcompiter)
507  {
508  RawTowerDefs::CalorimeterId calo_id = RawTowerDefs::NONE;
509 
510  switch ( jcompiter->first )
511  {
512  case Jet::CEMC_TOWER:
513  calo_id = RawTowerDefs::CEMC;
514  break;
515  case Jet::HCALIN_TOWER:
516  calo_id = RawTowerDefs::HCALIN;
517  break;
518  case Jet::HCALOUT_TOWER:
519  calo_id = RawTowerDefs::HCALOUT;
520  break;
521  default:
522  break;
523  }
524 
525  /* continue if no calorimeter id found */
526  if ( calo_id == RawTowerDefs::NONE )
527  continue;
528 
529  /* get tower container from map, find tower in tower container, get tower energy */
530  float e_component = 0;
531  if ( map_calotower->find( calo_id ) != map_calotower->end() )
532  e_component = ( ( ( map_calotower->find( calo_id ) )->second ).first )->getTower( jcompiter->second )->get_energy();
533 
534  /* check if energy is above threshold and count up matching counters accordingly */
535  if ( e_component > 0.1 )
536  jet_ncomp_above_0p1++;
537  if ( e_component > 1 )
538  jet_ncomp_above_1++;
539  if ( e_component > 10 )
540  jet_ncomp_above_10++;
541  }
542 
543  /* set tau candidate jet properties */
544  (iter->second)->set_property( PidCandidate::jet_eta , jetx->get_eta() );
545  (iter->second)->set_property( PidCandidate::jet_phi , jetx->get_phi() );
546  (iter->second)->set_property( PidCandidate::jet_etotal , jetx->get_e() );
547  (iter->second)->set_property( PidCandidate::jet_etrans , jetx->get_et() );
548  (iter->second)->set_property( PidCandidate::jet_ptotal , jetx->get_p() );
549  (iter->second)->set_property( PidCandidate::jet_ptrans , jetx->get_pt() );
550  (iter->second)->set_property( PidCandidate::jet_minv , jetx->get_mass() );
551  (iter->second)->set_property( PidCandidate::jet_mtrans , jet_mtrans );
552  (iter->second)->set_property( PidCandidate::jet_ncomp , (uint)jetx->size_comp() );
553  (iter->second)->set_property( PidCandidate::jet_ncomp_above_0p1 , jet_ncomp_above_0p1 );
554  (iter->second)->set_property( PidCandidate::jet_ncomp_above_1 , jet_ncomp_above_1 );
555  (iter->second)->set_property( PidCandidate::jet_ncomp_above_10 , jet_ncomp_above_10 );
556  (iter->second)->set_property( PidCandidate::jet_ncomp_emcal , (uint)jetx->count_comp( Jet::CEMC_TOWER ) );
557  }
558 
559  return 0;
560 }
561 
562 
563 
564 int
565 LeptoquarksReco::AddJetStructureInformation( type_map_tcan& tauCandidateMap, type_map_cdata* map_towers )
566 {
567  /* Cone size around jet axis within which to look for tracks */
568  float delta_R_cutoff_r1 = 0.1;
569  float delta_R_cutoff_r2 = 0.2;
570  float delta_R_cutoff_r3 = 0.3;
571  float delta_R_cutoff_r4 = 0.4;
572  float delta_R_cutoff_r5 = 0.5;
573 
574  /* energy threshold for considering tower */
575  float tower_emin = 0.0;
576 
577  /* number of steps for finding r90 */
578  int n_steps = 50;
579 
580  /* Loop over tau candidates */
581  for (type_map_tcan::iterator iter = tauCandidateMap.begin();
582  iter != tauCandidateMap.end();
583  ++iter)
584  {
585  /* get jet axis */
586  float jet_eta = (iter->second)->get_property_float( PidCandidate::jet_eta );
587  float jet_phi = (iter->second)->get_property_float( PidCandidate::jet_phi );
588  float jet_e = (iter->second)->get_property_float( PidCandidate::jet_etotal );
589 
590  /* collect jet structure properties */
591  float er1 = 0;
592  float er2 = 0;
593  float er3 = 0;
594  float er4 = 0;
595  float er5 = 0;
596  float r90 = 0;
597  float radius = 0;
598  float rms = 0;
599  float rms_esum = 0;
600 
601  /* collect jet structure properties- EMCal ONLY */
602  float emcal_er1 = 0;
603  float emcal_er2 = 0;
604  float emcal_er3 = 0;
605  float emcal_er4 = 0;
606  float emcal_er5 = 0;
607  float emcal_r90 = 0;
608  float emcal_radius = 0;
609  float emcal_rms = 0;
610  float emcal_rms_esum = 0;
611 
612  /* Loop over all tower (and geometry) collections */
613  for (type_map_cdata::iterator iter_calo = map_towers->begin();
614  iter_calo != map_towers->end();
615  ++iter_calo)
616  {
617  /* define tower iterator */
618  RawTowerContainer::ConstRange begin_end = ((iter_calo->second).first)->getTowers();
619  RawTowerContainer::ConstIterator rtiter;
620 
621  /* loop over all tower in CEMC calorimeter */
622  for (rtiter = begin_end.first; rtiter != begin_end.second; ++rtiter)
623  {
624  /* get tower energy */
625  RawTower *tower = rtiter->second;
626  float tower_energy = tower->get_energy();
627 
628  /* check if tower above energy treshold */
629  if ( tower_energy < tower_emin )
630  continue;
631 
632  /* get eta and phi of tower and check angle delta_R w.r.t. jet axis */
633  RawTowerGeom * tower_geom = ((iter_calo->second).second)->get_tower_geometry(tower -> get_key());
634  float tower_eta = tower_geom->get_eta();
635  float tower_phi = tower_geom->get_phi();
636 
637  /* If accounting for displaced vertex, need to calculate eta and phi:
638  double r = tower_geom->get_center_radius();
639  double phi = atan2(tower_geom->get_center_y(), tower_geom->get_center_x());
640  double z0 = tower_geom->get_center_z();
641  double z = z0 - vtxz;
642  double eta = asinh(z/r); // eta after shift from vertex
643  */
644 
645  float delta_R = CalculateDeltaR( tower_eta , tower_phi, jet_eta, jet_phi );
646 
647  /* if save_towers set true: add tower to tree */
648  if ( _save_towers )
649  {
650  float tower_data[17] = {(float) _ievent,
651  (float) (iter->second)->get_property_uint( PidCandidate::jet_id ),
652  (float) (iter->second)->get_property_int( PidCandidate::evtgen_pid ),
653  (float) (iter->second)->get_property_float( PidCandidate::evtgen_etotal ),
654  (float) (iter->second)->get_property_float( PidCandidate::evtgen_eta ),
655  (float) (iter->second)->get_property_float( PidCandidate::evtgen_phi ),
656  (float) (iter->second)->get_property_uint( PidCandidate::evtgen_decay_prong ),
657  (float) (iter->second)->get_property_uint( PidCandidate::evtgen_decay_hcharged ),
658  (float) (iter->second)->get_property_uint( PidCandidate::evtgen_decay_lcharged ),
659  (float) (iter->second)->get_property_float( PidCandidate::jet_eta ),
660  (float) (iter->second)->get_property_float( PidCandidate::jet_phi ),
661  (float) (iter->second)->get_property_float( PidCandidate::jet_etotal ),
662  (float) (iter_calo->first),
663  (float) tower_eta,
664  (float) tower_phi,
665  (float) delta_R,
666  (float) tower_energy
667  };
668 
669  _ntp_tower->Fill(tower_data);
670  }
671 
672  /* check delta R distance tower from jet axis */
673  if ( delta_R <= delta_R_cutoff_r1 )
674  {
675  er1 += tower_energy;
676  if ( (iter_calo->first) == RawTowerDefs::CEMC )
677  emcal_er1 += tower_energy;
678  }
679  if ( delta_R <= delta_R_cutoff_r2 )
680  {
681  er2 += tower_energy;
682  if ( (iter_calo->first) == RawTowerDefs::CEMC )
683  emcal_er2 += tower_energy;
684  }
685  if ( delta_R <= delta_R_cutoff_r3 )
686  {
687  er3 += tower_energy;
688  if ( (iter_calo->first) == RawTowerDefs::CEMC )
689  emcal_er3 += tower_energy;
690  }
691  if ( delta_R <= delta_R_cutoff_r4 )
692  {
693  er4 += tower_energy;
694  if ( (iter_calo->first) == RawTowerDefs::CEMC )
695  emcal_er4 += tower_energy;
696  }
697  if ( delta_R <= delta_R_cutoff_r5 )
698  {
699  er5 += tower_energy;
700  if ( (iter_calo->first) == RawTowerDefs::CEMC )
701  emcal_er5 += tower_energy;
702  }
703 
704  if ( delta_R <= delta_R_cutoff_r5 )
705  {
706  rms += tower_energy*delta_R*delta_R;
707  rms_esum += tower_energy;
708 
709  radius += tower_energy*delta_R;
710 
711  if ( (iter_calo->first) == RawTowerDefs::CEMC )
712  {
713  emcal_rms += tower_energy*delta_R*delta_R;
714  emcal_rms_esum += tower_energy;
715  emcal_radius += tower_energy*delta_R;
716  }
717  }
718  }
719  }
720 
721  /* finalize calculation of rms and radius */
722  if ( rms_esum > 0 )
723  {
724  radius /= rms_esum;
725  rms /= rms_esum;
726  rms = sqrt( rms );
727  }
728  else
729  {
730  radius = -1;
731  rms = -1;
732  }
733  if ( emcal_rms_esum > 0 )
734  {
735  emcal_radius /= emcal_rms_esum;
736  emcal_rms /= emcal_rms_esum;
737  emcal_rms = sqrt( emcal_rms );
738  }
739  else
740  {
741  emcal_radius = -1;
742  emcal_rms = -1;
743  }
744 
745  /* Search for cone angle that contains 90% of jet energy */
746  for(int r_i = 1; r_i<n_steps+1; r_i++){
747  float e_tower_sum = 0;
748 
749  /* Loop over all tower (and geometry) collections */
750  for (type_map_cdata::iterator iter_calo = map_towers->begin();
751  iter_calo != map_towers->end();
752  ++iter_calo)
753  {
754  /* define tower iterator */
755  RawTowerContainer::ConstRange begin_end = ((iter_calo->second).first)->getTowers();
756  RawTowerContainer::ConstIterator rtiter;
757 
758  if (e_tower_sum >= 0.9*jet_e) break;
759 
760  for (rtiter = begin_end.first; rtiter != begin_end.second; ++rtiter)
761  {
762  RawTower *tower = rtiter->second;
763  float tower_energy = tower->get_energy();
764 
765  /* check if tower is above minimum tower energy required */
766  if ( tower_energy < tower_emin )
767  continue;
768 
769  RawTowerGeom * tower_geom = ((iter_calo->second).second)->get_tower_geometry(tower -> get_key());
770  float tower_eta = tower_geom->get_eta();
771  float tower_phi = tower_geom->get_phi();
772 
773  float delta_R = CalculateDeltaR( tower_eta , tower_phi, jet_eta, jet_phi );
774 
775  if(delta_R < r_i*delta_R_cutoff_r5/n_steps) {
776  e_tower_sum = e_tower_sum + tower_energy;
777  r90 = r_i*delta_R_cutoff_r5/n_steps;
778  }
779  }
780 
781  if (e_tower_sum >= 0.9*jet_e) break;
782  }
783  }
784 
785  /* set tau candidate properties */
786  (iter->second)->set_property( PidCandidate::jetshape_econe_r01, er1 );
787  (iter->second)->set_property( PidCandidate::jetshape_econe_r02, er2 );
788  (iter->second)->set_property( PidCandidate::jetshape_econe_r03, er3 );
789  (iter->second)->set_property( PidCandidate::jetshape_econe_r04, er4 );
790  (iter->second)->set_property( PidCandidate::jetshape_econe_r05, er5 );
791  (iter->second)->set_property( PidCandidate::jetshape_r90, r90 );
792  (iter->second)->set_property( PidCandidate::jetshape_rms, rms );
793  (iter->second)->set_property( PidCandidate::jetshape_radius, radius );
794  (iter->second)->set_property( PidCandidate::jetshape_emcal_econe_r01, emcal_er1 );
795  (iter->second)->set_property( PidCandidate::jetshape_emcal_econe_r02, emcal_er2 );
796  (iter->second)->set_property( PidCandidate::jetshape_emcal_econe_r03, emcal_er3 );
797  (iter->second)->set_property( PidCandidate::jetshape_emcal_econe_r04, emcal_er4 );
798  (iter->second)->set_property( PidCandidate::jetshape_emcal_econe_r05, emcal_er5 );
799  (iter->second)->set_property( PidCandidate::jetshape_emcal_r90, emcal_r90 );
800  (iter->second)->set_property( PidCandidate::jetshape_emcal_rms, emcal_rms );
801  (iter->second)->set_property( PidCandidate::jetshape_emcal_radius, emcal_radius );
802  }
803 
804  return 0;
805 }
806 
807 int
808 LeptoquarksReco::AddTrackInformation( type_map_tcan& tauCandidateMap, SvtxTrackMap* trackmap, SvtxVertexMap* vertexmap, SvtxEvalStack *svtxevalstack, double R_max )
809 {
810 
811  // Pointers for tracks //
812  SvtxTrackEval* trackeval = svtxevalstack->get_track_eval();
813  SvtxTruthEval* trutheval = svtxevalstack->get_truth_eval();
814 
815  /* Loop over tau candidates */
816  for (type_map_tcan::iterator iter = tauCandidateMap.begin();
817  iter != tauCandidateMap.end();
818  ++iter)
819  {
820  uint tracks_count_r02 = 0;
821  int tracks_chargesum_r02 = 0;
822  float tracks_rmax_r02 = 0;
823 
824  uint tracks_count_r04 = 0;
825  int tracks_chargesum_r04 = 0;
826  float tracks_rmax_r04 = 0;
827 
828  uint tracks_count_R = 0;
829  int tracks_chargesum_R = 0;
830  float tracks_rmax_R = 0;
831 
832  vector<float> tracks_vertex;
833  vector<float> temp_vertex;
834 
835  float jet_eta = (iter->second)->get_property_float( PidCandidate::jet_eta );
836  float jet_phi = (iter->second)->get_property_float( PidCandidate::jet_phi );
837 
838  /* Loop over tracks
839  * (float) track->get_eta(), //eta of the track
840  * (float) track->get_phi(), //phi of the track
841  * (float) track->get_pt(), //transverse momentum of track
842  * (float) track->get_p(), //total momentum of track
843  * (float) track->get_charge(), //electric charge of track
844  * (float) track->get_quality() //track quality */
845 
846  //Loop over tracks in event //
847  for (SvtxTrackMap::ConstIter track_itr = trackmap->begin();
848  track_itr != trackmap->end(); track_itr++) {
849 
850  SvtxTrack* track = dynamic_cast<SvtxTrack*>(track_itr->second);
851 
852  // Get track variables //
853  float track_eta = track->get_eta();
854  float track_phi = track->get_phi();
855  int track_charge = track->get_charge();
856  double gvx,gvy,gvz;
857 
858  float delta_R = CalculateDeltaR( track_eta, track_phi, jet_eta, jet_phi );
859 
860  PHG4Particle* g4particle = trackeval->max_truth_particle_by_nclusters(track);
861 
862  // Get true vertex distances //
863  PHG4VtxPoint* vtx = trutheval->get_vertex(g4particle);
864  gvx = vtx->get_x();
865  gvy = vtx->get_y();
866  gvz = vtx->get_z();
867 
868 
869 
870  // If charged track is within jet then use its vertex //
871  if(delta_R < 0.5 && trutheval->is_primary(g4particle)) tracks_vertex.push_back(sqrt(pow(gvx,2)+pow(gvy,2)+pow(gvz,2)));
872 
873  /* if save_tracks set true: add track to tree */
874  if ( _save_tracks )
875  {
876  float track_data[17] = {(float) _ievent,
877  (float) (iter->second)->get_property_uint( PidCandidate::jet_id ),
878  (float) (iter->second)->get_property_int( PidCandidate::evtgen_pid ),
879  (float) (iter->second)->get_property_float( PidCandidate::evtgen_etotal ),
880  (float) (iter->second)->get_property_float( PidCandidate::evtgen_eta ),
881  (float) (iter->second)->get_property_float( PidCandidate::evtgen_phi ),
882  (float) (iter->second)->get_property_uint( PidCandidate::evtgen_decay_prong ),
883  (float) (iter->second)->get_property_uint( PidCandidate::evtgen_decay_hcharged ),
884  (float) (iter->second)->get_property_uint( PidCandidate::evtgen_decay_lcharged ),
885  (float) (iter->second)->get_property_float( PidCandidate::jet_eta ),
886  (float) (iter->second)->get_property_float( PidCandidate::jet_phi ),
887  (float) (iter->second)->get_property_float( PidCandidate::jet_etotal ),
888  (float) track->get_quality(),
889  (float) track_eta,
890  (float) track_phi,
891  (float) delta_R,
892  (float) track->get_p()
893  };
894 
895  _ntp_track->Fill(track_data);
896  }
897 
898  /* If track within search cone, update track information for tau candidate */
899  if ( delta_R < 0.2 )
900  {
901  tracks_count_r02++;
902  tracks_chargesum_r02 += track_charge;
903 
904  if ( delta_R > tracks_rmax_r02 )
905  tracks_rmax_r02 = delta_R;
906  }
907 
908  if ( delta_R < 0.4 )
909  {
910  tracks_count_r04++;
911  tracks_chargesum_r04 += track_charge;
912 
913  if ( delta_R > tracks_rmax_r04 )
914  tracks_rmax_r04 = delta_R;
915  }
916 
917  if ( delta_R < R_max )
918  {
919  tracks_count_R++;
920  tracks_chargesum_R += track_charge;
921 
922  if ( delta_R > tracks_rmax_R )
923  tracks_rmax_R = delta_R;
924  }
925 
926  } // end loop over reco tracks //
927  // sort vertex array in increasing order //
928  std::sort(tracks_vertex.begin(),tracks_vertex.end());
929 
930  // Compute average vertex distance of tracks in jet //
931  float avg = Average(tracks_vertex);
932 
933  /* Set track-based properties for tau candidate */
934  (iter->second)->set_property( PidCandidate::tracks_count_r02, tracks_count_r02 );
935  (iter->second)->set_property( PidCandidate::tracks_chargesum_r02, tracks_chargesum_r02 );
936  (iter->second)->set_property( PidCandidate::tracks_rmax_r02, tracks_rmax_r02 );
937  (iter->second)->set_property( PidCandidate::tracks_count_r04, tracks_count_r04 );
938  (iter->second)->set_property( PidCandidate::tracks_chargesum_r04, tracks_chargesum_r04 );
939  (iter->second)->set_property( PidCandidate::tracks_rmax_r04, tracks_rmax_r04 );
940  (iter->second)->set_property( PidCandidate::tracks_count_R, tracks_count_R );
941  (iter->second)->set_property( PidCandidate::tracks_chargesum_R, tracks_chargesum_R );
942  (iter->second)->set_property( PidCandidate::tracks_rmax_R, tracks_rmax_R );
943  if(avg == avg) (iter->second)->set_property( PidCandidate::tracks_vertex, avg);
944 
945  } // end loop over tau candidates
946 
947  return 0;
948 }
949 
950 
951 int
952 LeptoquarksReco::WritePidCandidatesToTree( type_map_tcan& tauCandidateMap )
953 {
954  /* Loop over all tau candidates and add them to tree*/
955  for (type_map_tcan::iterator iter = tauCandidateMap.begin();
956  iter != tauCandidateMap.end();
957  ++iter)
958  {
959  /* update information in map and fill tree */
960  for ( map< PidCandidate::PROPERTY , vector<float> >::iterator iter_prop = _map_tau_candidate_branches.begin();
961  iter_prop != _map_tau_candidate_branches.end();
962  ++iter_prop)
963  {
964  switch ( PidCandidate::get_property_info( (iter_prop->first) ).second ) {
965 
966  case PidCandidate::type_float :
967  (iter_prop->second).push_back( (iter->second)->get_property_float( (iter_prop->first) ) );
968  break;
969 
970  case PidCandidate::type_int :
971  (iter_prop->second).push_back( (iter->second)->get_property_int( (iter_prop->first) ) );
972  break;
973 
974  case PidCandidate::type_uint :
975  (iter_prop->second).push_back( (iter->second)->get_property_uint( (iter_prop->first) ) );
976  break;
977 
978  case PidCandidate::type_unknown :
979  break;
980  }
981  }
982  }
983 
984  return 0;
985 }
986 
987 
988 PidCandidate*
989 LeptoquarksReco::FindMinDeltaRCandidate( type_map_tcan *candidates, const float eta_ref, const float phi_ref )
990 {
991  /* PidCandidate with eta, phi closest to reference */
992  PidCandidate* best_candidate = NULL;
993 
994  float eta_ref_local = eta_ref;
995  float phi_ref_local = phi_ref;
996 
997  /* Event record uses 0 < phi < 2Pi, while Fun4All uses -Pi < phi < Pi.
998  * Therefore, correct angle for 'wraparound' at phi == Pi */
999  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 }