Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
DISKinematicsReco.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file DISKinematicsReco.C
1 #include "DISKinematicsReco.h"
2 #include "PidCandidatev1.h"
3 #include "TruthTrackerHepMC.h"
4 #include "TrackProjectionTools.h"
5 #include "TrackProjectorPlaneECAL.h"
6 /* STL includes */
7 #include <cassert>
8 
9 /* Fun4All includes */
10 #include <trackbase_historic/SvtxTrack.h>
11 #include <trackbase_historic/SvtxTrackMap.h>
12 #include <trackbase_historic/SvtxTrack_FastSim.h>
13 
14 #include <trackbase_historic/SvtxTrackMap_v1.h>
15 
16 #include <phool/getClass.h>
17 
18 #include <fun4all/Fun4AllReturnCodes.h>
19 
20 #include <g4jets/JetMap.h>
21 
22 #include <calobase/RawTowerGeomContainer.h>
23 #include <calobase/RawTowerContainer.h>
24 #include <calobase/RawTowerGeom.h>
25 #include <calobase/RawTowerv1.h>
26 #include <calobase/RawTowerDefs.h>
27 
28 #include <calobase/RawClusterContainer.h>
29 #include <calobase/RawCluster.h>
30 
31 #include <g4vertex/GlobalVertexMap.h>
32 #include <g4vertex/GlobalVertex.h>
33 
34 #include <g4main/PHG4Particle.h>
35 
36 #include <g4eval/CaloEvalStack.h>
37 #include <g4eval/CaloRawTowerEval.h>
38 
39 /* ROOT includes */
40 #include <TLorentzVector.h>
41 #include <TString.h>
42 #include <TNtuple.h>
43 #include <TTree.h>
44 #include <TFile.h>
45 #include <TDatabasePDG.h>
46 
47 using namespace std;
48 
49 DISKinematicsReco::DISKinematicsReco(std::string filename) :
50  SubsysReco("DISKinematicsReco" ),
51  _mproton( 0.938272 ),
52  _save_towers(false),
53  _save_tracks(false),
54  _do_process_geant4_cluster(false),
55  _do_process_truth(false),
56  _ievent(0),
57  _filename(filename),
58  _tfile(nullptr),
59  _tree_event_cluster(nullptr),
60  _tree_event_truth(nullptr),
61  _beam_electron_ptotal(10),
62  _beam_hadron_ptotal(250),
63  _trackproj(nullptr)
64 {
65 
66 }
67 
68 int
69 DISKinematicsReco::Init(PHCompositeNode *topNode)
70 {
71  _topNode = topNode;
72 
73  _ievent = 0;
74 
75  _tfile = new TFile(_filename.c_str(), "RECREATE");
76 
77  /* Add PidCandidate properties to map that defines output tree */
78  float dummy = 0;
79  vector< float > vdummy;
80 
81  _map_em_candidate_branches.insert( make_pair( PidCandidate::em_cluster_id , vdummy ) );
82  _map_em_candidate_branches.insert( make_pair( PidCandidate::em_cluster_prob , vdummy ) );
83  _map_em_candidate_branches.insert( make_pair( PidCandidate::em_cluster_posx , vdummy ) );
84  _map_em_candidate_branches.insert( make_pair( PidCandidate::em_cluster_posy , vdummy ) );
85  _map_em_candidate_branches.insert( make_pair( PidCandidate::em_cluster_posz , vdummy ) );
86  _map_em_candidate_branches.insert( make_pair( PidCandidate::em_cluster_e , vdummy ) );
87  _map_em_candidate_branches.insert( make_pair( PidCandidate::em_cluster_ecore , vdummy ) );
88  _map_em_candidate_branches.insert( make_pair( PidCandidate::em_cluster_et_iso , vdummy ) );
89  _map_em_candidate_branches.insert( make_pair( PidCandidate::em_cluster_theta , vdummy ) );
90  _map_em_candidate_branches.insert( make_pair( PidCandidate::em_cluster_phi , vdummy ) );
91  _map_em_candidate_branches.insert( make_pair( PidCandidate::em_cluster_eta , vdummy ) );
92  _map_em_candidate_branches.insert( make_pair( PidCandidate::em_cluster_pt , vdummy ) );
93  _map_em_candidate_branches.insert( make_pair( PidCandidate::em_cluster_ntower , vdummy ) );
94  _map_em_candidate_branches.insert( make_pair( PidCandidate::em_cluster_caloid , vdummy ) );
95 
96  _map_em_candidate_branches.insert( make_pair( PidCandidate::em_track_id , vdummy ) );
97  _map_em_candidate_branches.insert( make_pair( PidCandidate::em_track_quality , vdummy ) );
98  _map_em_candidate_branches.insert( make_pair( PidCandidate::em_track_theta , vdummy ) );
99  _map_em_candidate_branches.insert( make_pair( PidCandidate::em_track_phi , vdummy ) );
100  _map_em_candidate_branches.insert( make_pair( PidCandidate::em_track_eta , vdummy ) );
101  _map_em_candidate_branches.insert( make_pair( PidCandidate::em_track_ptotal , vdummy ) );
102  _map_em_candidate_branches.insert( make_pair( PidCandidate::em_track_ptrans , vdummy ) );
103  _map_em_candidate_branches.insert( make_pair( PidCandidate::em_track_charge , vdummy ) );
104  _map_em_candidate_branches.insert( make_pair( PidCandidate::em_track_dca , vdummy ) );
105  _map_em_candidate_branches.insert( make_pair( PidCandidate::em_track_section , vdummy ) );
106  _map_em_candidate_branches.insert( make_pair( PidCandidate::em_track_e3x3_cemc , vdummy ) );
107  _map_em_candidate_branches.insert( make_pair( PidCandidate::em_track_e3x3_femc , vdummy ) );
108  _map_em_candidate_branches.insert( make_pair( PidCandidate::em_track_e3x3_eemc , vdummy ) );
109  _map_em_candidate_branches.insert( make_pair( PidCandidate::em_track_e3x3_ihcal , vdummy ) );
110  _map_em_candidate_branches.insert( make_pair( PidCandidate::em_track_e3x3_ohcal , vdummy ) );
111  _map_em_candidate_branches.insert( make_pair( PidCandidate::em_track_e3x3_fhcal , vdummy ) );
112  _map_em_candidate_branches.insert( make_pair( PidCandidate::em_track_e3x3_ehcal , vdummy ) );
113  _map_em_candidate_branches.insert( make_pair( PidCandidate::em_track_cluster_dr , vdummy ) );
114 
115  _map_em_candidate_branches.insert( make_pair( PidCandidate::em_track_theta2cluster , vdummy ) );
116  _map_em_candidate_branches.insert( make_pair( PidCandidate::em_track_eta2cluster , vdummy ) );
117  _map_em_candidate_branches.insert( make_pair( PidCandidate::em_track_phi2cluster , vdummy ) );
118  _map_em_candidate_branches.insert( make_pair( PidCandidate::em_track_p2cluster , vdummy ) );
119  _map_em_candidate_branches.insert( make_pair( PidCandidate::em_track_x2cluster , vdummy ) );
120  _map_em_candidate_branches.insert( make_pair( PidCandidate::em_track_y2cluster , vdummy ) );
121  _map_em_candidate_branches.insert( make_pair( PidCandidate::em_track_z2cluster , vdummy ) );
122 
123  _map_em_candidate_branches.insert( make_pair( PidCandidate::em_pid_prob_electron , vdummy ) );
124  _map_em_candidate_branches.insert( make_pair( PidCandidate::em_pid_prob_pion , vdummy ) );
125  _map_em_candidate_branches.insert( make_pair( PidCandidate::em_pid_prob_kaon , vdummy ) );
126  _map_em_candidate_branches.insert( make_pair( PidCandidate::em_pid_prob_proton , vdummy ) );
127 
128  _map_em_candidate_branches.insert( make_pair( PidCandidate::em_evtgen_pid , vdummy ) );
129  _map_em_candidate_branches.insert( make_pair( PidCandidate::em_evtgen_ptotal , vdummy ) );
130  _map_em_candidate_branches.insert( make_pair( PidCandidate::em_evtgen_theta , vdummy ) );
131  _map_em_candidate_branches.insert( make_pair( PidCandidate::em_evtgen_phi , vdummy ) );
132  _map_em_candidate_branches.insert( make_pair( PidCandidate::em_evtgen_eta , vdummy ) );
133  _map_em_candidate_branches.insert( make_pair( PidCandidate::em_evtgen_charge , vdummy ) );
134  _map_em_candidate_branches.insert( make_pair( PidCandidate::em_evtgen_is_scattered_lepton , vdummy ) );
135 
136  _map_em_candidate_branches.insert( make_pair( PidCandidate::em_reco_x_e , vdummy ) );
137  _map_em_candidate_branches.insert( make_pair( PidCandidate::em_reco_y_e , vdummy ) );
138  _map_em_candidate_branches.insert( make_pair( PidCandidate::em_reco_q2_e , vdummy ) );
139  _map_em_candidate_branches.insert( make_pair( PidCandidate::em_reco_w_e , vdummy ) );
140 
141  /* Add branches to map that defines output tree for event-wise properties */
142  /*
143  runnr
144  eventnr
145  eventweight
146  cal_px
147  cal_py
148  cal_pz
149  cal_etotal
150  cal_et
151  cal_empz
152  cal_pt --> "pt_miss" ?
153  cal_pt_phi = atan2(sum(py)/sum(px)) --> "pt_miss_angle" ?
154  */
155  _map_event_branches.insert( make_pair( "em_ncandidates" , dummy ) );
156 
157  _map_event_branches.insert( make_pair( "Et_miss" , dummy ) );
158  _map_event_branches.insert( make_pair( "Et_miss_phi" , dummy ) );
159 
160  _map_event_branches.insert( make_pair( "evtgen_process_id" , dummy ) );
161  _map_event_branches.insert( make_pair( "evtgen_s" , dummy ) );
162  _map_event_branches.insert( make_pair( "evtgen_x" , dummy ) );
163  _map_event_branches.insert( make_pair( "evtgen_Q2" , dummy ) );
164  _map_event_branches.insert( make_pair( "evtgen_y" , dummy ) );
165  _map_event_branches.insert( make_pair( "evtgen_W" , dummy ) );
166 
167  /* Create tree for information about full event */
168  _tree_event_cluster = new TTree("event_cluster", "a Tree with global event information and EM candidates");
169  _tree_event_cluster->Branch("event", &_ievent, "event/I");
170 
171  /* Add event branches */
172  for ( map< string , float >::iterator iter = _map_event_branches.begin();
173  iter != _map_event_branches.end();
174  ++iter)
175  {
176  _tree_event_cluster->Branch( (iter->first).c_str(),
177  &(iter->second) );
178  }
179 
180  /* Add EM candidate branches */
181  for ( map< PidCandidate::PROPERTY , vector<float> >::iterator iter = _map_em_candidate_branches.begin();
182  iter != _map_em_candidate_branches.end();
183  ++iter)
184  {
185  _tree_event_cluster->Branch( PidCandidate::get_property_info( (iter->first) ).first.c_str(),
186  &(iter->second) );
187  }
188 
189  /* create clone of tree for truth particle candidates */
190  _tree_event_truth = _tree_event_cluster->CloneTree();
191  _tree_event_truth->SetName("event_truth");
192  _tree_event_truth->SetTitle("a Tree with global event information and truth particle candidates");
193 
194 
195  return 0;
196 }
197 
198 int
199 DISKinematicsReco::InitRun(PHCompositeNode *topNode)
200 {
201  if(_do_process_geant4_cluster)
202  _trackproj=new TrackProjectorPlaneECAL( topNode );
203  return 0;
204 }
205 
206 int
207 DISKinematicsReco::process_event(PHCompositeNode *topNode)
208 {
209  /* Find electron candidates based on calorimeter cluster from full Geant4 detector simulation */
210  if ( _do_process_geant4_cluster )
211  {
212  /* Reset branch map */
213  ResetBranchMap();
214 
215  /* Create map to collect electron candidates.
216  * Use energy as 'key' to the map because energy is unique for each jet, while there are sometimes multiple jets (same energy,
217  * different jet ID) in the input jet collection. Also, map automatically sorts entries by key, i.e. this gives list of tau candidates
218  * sorted by energy. */
219  type_map_tcan electronCandidateMap;
220 
221  /* Collect EM candidates from calorimeter cluster */
222  CollectEmCandidatesFromCluster( electronCandidateMap );
223 
224  /* Calculate kinematics for each em candidate */
225  AddReconstructedKinematics( electronCandidateMap, "cluster" );
226 
227  /* Add information about em candidats to output tree */
228  WriteCandidatesToTree( electronCandidateMap );
229 
230  /* Add global event information */
231  AddGlobalCalorimeterInformation();
232  AddTruthEventInformation();
233 
234  /* fill event information tree */
235  _tree_event_cluster->Fill();
236  }
237 
238  /* Find electron candidates based on truth particle inforamtion */
239  if ( _do_process_truth )
240  {
241  /* reset branches, fill with truth particle candidates, and fill different tree */
242  ResetBranchMap();
243 
244  /* Collect electron candidates based on TRUTH information */
245  type_map_tcan electronTruthCandidateMap;
246 
247  CollectEmCandidatesFromTruth( electronTruthCandidateMap );
248  AddReconstructedKinematics( electronTruthCandidateMap, "truth" );
249  WriteCandidatesToTree( electronTruthCandidateMap );
250  AddTruthEventInformation();
251 
252  /* fill event information tree */
253  _tree_event_truth->Fill();
254  }
255 
256  /* count up event number */
257  _ievent ++;
258 
259  return 0;
260 }
261 
262 
263 
264 int
265 DISKinematicsReco::CollectEmCandidatesFromCluster( type_map_tcan& electronCandidateMap )
266 {
267 
268  /* Loop over all clusters in EEMC, CEMC, and FEMC to find electron candidates */
269  vector< string > v_ecals;
270  v_ecals.push_back("EEMC");
271  v_ecals.push_back("CEMC");
272  v_ecals.push_back("FEMC");
273  for ( unsigned idx = 0; idx < v_ecals.size(); idx++ )
274  {
275  CaloEvalStack * caloevalstack = new CaloEvalStack(_topNode, v_ecals.at( idx ) );
276 
277  string clusternodename = "CLUSTER_" + v_ecals.at( idx );
278  RawClusterContainer *clusterList = findNode::getClass<RawClusterContainer>(_topNode,clusternodename.c_str());
279  SvtxTrackMap *trackmap = findNode::getClass<SvtxTrackMap>(_topNode,"TrackMap");
280 
281  if (!clusterList) {
282  cerr << PHWHERE << " ERROR: Can't find node " << clusternodename << endl;
283  return false;
284  }
285  if(!trackmap) {
286  cerr << PHWHERE << " ERROR: Can't find node SvtxTrackMap" << endl;
287  return false;
288  }
289 
290  for (unsigned int k = 0; k < clusterList->size(); ++k)
291  {
292  RawCluster *cluster = clusterList->getCluster(k);
293  /* Check if cluster energy is below threshold */
294  float e_cluster_threshold = 0.3;
295  if ( cluster->get_energy() < e_cluster_threshold )
296  continue;
297 
298  InsertCandidateFromCluster( electronCandidateMap , cluster , caloevalstack , trackmap);
299  }
300  }
301 
302  return 0;
303 }
304 
305 
306 int
307 DISKinematicsReco::CollectEmCandidatesFromTruth( type_map_tcan& candidateMap )
308 {
309  /* Get collection of truth particles from event generator */
310  PHHepMCGenEventMap *geneventmap = findNode::getClass<PHHepMCGenEventMap>(_topNode,"PHHepMCGenEventMap");
311  if (!geneventmap) {
312  //std::cout << PHWHERE << " WARNING: Can't find requested PHHepMCGenEventMap" << endl;
313  return -1;
314  }
315 
316  /* Add truth kinematics */
317  int embedding_id = 1;
318  PHHepMCGenEvent *genevt = geneventmap->get(embedding_id);
319  if (!genevt)
320  {
321  std::cout << PHWHERE << "WARNING: Node PHHepMCGenEventMap missing subevent with embedding ID "<< embedding_id;
322  std::cout <<". Print PHHepMCGenEventMap:";
323  geneventmap->identify();
324  return -1;
325  }
326 
327  HepMC::GenEvent* theEvent = genevt->getEvent();
328 
329  if ( !theEvent )
330  {
331  std::cout << PHWHERE << "WARNING: Missing requested GenEvent!" << endl;
332  return -1;
333  }
334 
335  /* Look for beam particles and scattered lepton */
336  TruthTrackerHepMC truth;
337  truth.set_hepmc_geneventmap( geneventmap );
338  //HepMC::GenParticle* particle_beam_l = truth.FindBeamLepton();
339  //HepMC::GenParticle* particle_beam_h = truth.FindBeamHadron();
340  HepMC::GenParticle* particle_scattered_l = truth.FindScatteredLepton();
341  /* loop over all particles */
342  for (HepMC::GenEvent::particle_const_iterator p = theEvent->particles_begin();
343  p != theEvent->particles_end(); ++p)
344  {
345  TParticlePDG * pdg_p = TDatabasePDG::Instance()->GetParticle( (*p)->pdg_id() );
346  int charge = -999;
347  if ( pdg_p )
348  {
349  /* NOTE: TParticlePDG::Charge() returns charge in units of |e|/3 (see ROOT documentation) */
350  charge = pdg_p->Charge() / 3;
351  }
352 
353  /* skip particles that are not stable final state particles (status 1) */
354  if ( (*p)->status() != 1 )
355  continue;
356  /* create new pid candidate */
357  PidCandidatev1 *tc = new PidCandidatev1();
358  tc->set_candidate_id( candidateMap.size()+1 );
359 
360  float mom_eta = -1 * log ( tan( (*p)->momentum().theta() / 2.0 ) );
361 
362  tc->set_property( PidCandidate::em_evtgen_pid, (*p)->pdg_id() );
363  tc->set_property( PidCandidate::em_evtgen_ptotal, (float) (*p)->momentum().e() );
364  tc->set_property( PidCandidate::em_evtgen_theta, (float) (*p)->momentum().theta() );
365  tc->set_property( PidCandidate::em_evtgen_phi, (float) (*p)->momentum().phi() );
366  tc->set_property( PidCandidate::em_evtgen_eta, mom_eta );
367  tc->set_property( PidCandidate::em_evtgen_charge, charge );
368  tc->set_property( PidCandidate::em_evtgen_is_scattered_lepton, (uint)0 );
369 
370  /* is scattered lepton? */
371  if ( particle_scattered_l &&
372  (*p) == particle_scattered_l )
373  tc->set_property( PidCandidate::em_evtgen_is_scattered_lepton, (uint)1 );
374 
375  /* add pid candidate to collection */
376  candidateMap.insert( make_pair( tc->get_candidate_id(), tc ) );
377  }
378 
379  return 0;
380 }
381 
382 
383 int
384 DISKinematicsReco::InsertCandidateFromCluster( type_map_tcan& candidateMap , RawCluster *cluster , CaloEvalStack *caloevalstack, SvtxTrackMap *trackmap )
385 {
386 
387  /* create new electron candidate */
388  PidCandidatev1 *tc = new PidCandidatev1();
389  tc->set_candidate_id( candidateMap.size()+1 );
390 
391  /* set some initial cluster properties */
392  float theta = atan2( cluster->get_r() , cluster->get_z() );
393  float eta = -1 * log( tan( theta / 2.0 ) );
394  float pt = cluster->get_energy() * sin( theta );
395 
396  /* get calorimeter ID where towers of cluster are found */
397  unsigned caloid = 0;
398  RawCluster::TowerConstIterator rtiter = cluster->get_towers().first;
399  caloid = RawTowerDefs::decode_caloid( rtiter->first );
400  _trackproj->set_detector(_trackproj->get_detector_from_cluster(cluster));
401  //cout << "Calo ID: " << caloid << " -> " << RawTowerDefs::convert_caloid_to_name( RawTowerDefs::decode_caloid( rtiter->first ) ) << endl;
402 
403  tc->set_property( PidCandidate::em_cluster_id, cluster->get_id() );
404  tc->set_property( PidCandidate::em_cluster_prob, cluster->get_prob() );
405  tc->set_property( PidCandidate::em_cluster_posx, cluster->get_x() );
406  tc->set_property( PidCandidate::em_cluster_posy, cluster->get_y() );
407  tc->set_property( PidCandidate::em_cluster_posz, cluster->get_z() );
408  tc->set_property( PidCandidate::em_cluster_e, cluster->get_energy() );
409  tc->set_property( PidCandidate::em_cluster_ecore, cluster->get_ecore() );
410  tc->set_property( PidCandidate::em_cluster_et_iso, cluster->get_et_iso() );
411  tc->set_property( PidCandidate::em_cluster_theta, theta );
412  tc->set_property( PidCandidate::em_cluster_phi, cluster->get_phi() );
413  tc->set_property( PidCandidate::em_cluster_eta, eta );
414  tc->set_property( PidCandidate::em_cluster_pt, pt );
415  tc->set_property( PidCandidate::em_cluster_ntower, (unsigned)cluster->getNTowers() );
416  tc->set_property( PidCandidate::em_cluster_caloid, caloid );
417 
418  /* get track projection helper class */
419  TrackProjectionTools tpt( _topNode );
420 
421  /* find matching reco track */
422  float best_track_dr = NAN;
423  SvtxTrack* best_track = NULL; //tpt.FindClosestTrack( cluster, best_track_dr ); /* @TODO switch track finding back on as soon as we understand it better */
424 
425  /* IF matching track found: set track properties */
426  if ( best_track )
427  {
428  /* set some initial track properties */
429  float theta = 2.0 * atan( exp( -1 * best_track->get_eta() ) );
430 
431  tc->set_property( PidCandidate::em_track_id, (uint)best_track->get_id() );
432  tc->set_property( PidCandidate::em_track_quality, best_track->get_quality() );
433  tc->set_property( PidCandidate::em_track_theta, theta );
434  tc->set_property( PidCandidate::em_track_phi, best_track->get_phi() );
435  tc->set_property( PidCandidate::em_track_eta, best_track->get_eta() );
436  tc->set_property( PidCandidate::em_track_ptotal, best_track->get_p() );
437  tc->set_property( PidCandidate::em_track_ptrans, best_track->get_pt() );
438  tc->set_property( PidCandidate::em_track_charge, best_track->get_charge() );
439  tc->set_property( PidCandidate::em_track_dca, best_track->get_dca() );
440  tc->set_property( PidCandidate::em_track_section, (uint)0 );
441  tc->set_property( PidCandidate::em_track_e3x3_cemc, best_track->get_cal_energy_3x3(SvtxTrack::CEMC) );
442  tc->set_property( PidCandidate::em_track_e3x3_ihcal, best_track->get_cal_energy_3x3(SvtxTrack::HCALIN) );
443  tc->set_property( PidCandidate::em_track_e3x3_ohcal, best_track->get_cal_energy_3x3(SvtxTrack::HCALOUT) );
444  tc->set_property( PidCandidate::em_track_cluster_dr, best_track_dr );
445 
446  // Use the track states to project to the FEMC / FHCAL / EEMC and generate
447  // energy sums.
448  float e3x3_femc = NAN;
449  float e3x3_fhcal = NAN;
450  float e3x3_eemc = NAN;
451 
452  for (SvtxTrack::ConstStateIter state_itr = best_track->begin_states();
453  state_itr != best_track->end_states(); state_itr++) {
454 
455  SvtxTrackState *temp = dynamic_cast<SvtxTrackState*>(state_itr->second);
456 
457  if( (temp->get_name()=="FHCAL") )
458  e3x3_fhcal = tpt.getE33Forward( "FHCAL" , temp->get_x() , temp->get_y() );
459 
460  if( (temp->get_name()=="FEMC") )
461  e3x3_femc = tpt.getE33Forward( "FEMC" , temp->get_x() , temp->get_y() );
462 
463  if( (temp->get_name()=="EEMC") )
464  e3x3_eemc = tpt.getE33Forward( "EEMC" , temp->get_x() , temp->get_y() );
465  }
466 
467  /* set candidate properties */
468  tc->set_property( PidCandidate::em_track_e3x3_fhcal, e3x3_fhcal );
469  tc->set_property( PidCandidate::em_track_e3x3_femc, e3x3_femc );
470  tc->set_property( PidCandidate::em_track_e3x3_eemc, e3x3_eemc );
471 
472  /* Get information on track from PID detectors */
473  tc->set_property( PidCandidate::em_pid_prob_electron, (float)0.0 );
474  tc->set_property( PidCandidate::em_pid_prob_pion, (float)0.0 );
475  tc->set_property( PidCandidate::em_pid_prob_kaon, (float)0.0 );
476  tc->set_property( PidCandidate::em_pid_prob_proton, (float)0.0 );
477  }
478 
479  /* set em candidate MC truth properties */
480  tc->set_property( PidCandidate::em_evtgen_pid, (int)NAN );
481  tc->set_property( PidCandidate::em_evtgen_ptotal, (float)NAN );
482  tc->set_property( PidCandidate::em_evtgen_theta, (float)NAN );
483  tc->set_property( PidCandidate::em_evtgen_phi, (float)NAN );
484  tc->set_property( PidCandidate::em_evtgen_eta, (float)NAN );
485  tc->set_property( PidCandidate::em_evtgen_charge, (int)NAN );
486 
487  /* If matching truth primary particle found: update truth information */
488  CaloRawClusterEval* clustereval = caloevalstack->get_rawcluster_eval();
489  PHG4Particle* primary = clustereval->max_truth_primary_particle_by_energy(cluster);
490 
491  if ( primary )
492  {
493  /* get particle momenta and theta, phi angles */
494  float gpx = primary->get_px();
495  float gpy = primary->get_py();
496  float gpz = primary->get_pz();
497  float gpt = sqrt(gpx * gpx + gpy * gpy);
498  float gptotal = sqrt(gpx * gpx + gpy * gpy + gpz * gpz);
499  //float ge = (float)primary->get_e();
500 
501  float gphi = NAN;
502  float gtheta = NAN;
503  float geta = NAN;
504 
505  if (gpt != 0.0)
506  {
507  gtheta = atan2( gpt, gpz );
508  geta = -1 * log( tan( gtheta / 2.0 ) );
509  }
510 
511  gphi = atan2(gpy, gpx);
512 
513  /* get charge based on PDG code of particle */
514  TParticlePDG * pdg_p = TDatabasePDG::Instance()->GetParticle( primary->get_pid() );
515  int gcharge = -999;
516  if ( pdg_p )
517  {
518  /* NOTE: TParticlePDG::Charge() returns charge in units of |e|/3 (see ROOT documentation) */
519  gcharge = pdg_p->Charge() / 3;
520  }
521 
522  tc->set_property( PidCandidate::em_evtgen_pid, primary->get_pid() );
523  tc->set_property( PidCandidate::em_evtgen_ptotal, gptotal );
524  tc->set_property( PidCandidate::em_evtgen_theta, gtheta );
525  tc->set_property( PidCandidate::em_evtgen_phi, gphi );
526  tc->set_property( PidCandidate::em_evtgen_eta, geta );
527  tc->set_property( PidCandidate::em_evtgen_charge, gcharge );
528 
529 
530  /* Track pointing to cluster */
531 
532  bool fill_in_truth = false;
533 
534  /* Try to find a track which best points to the current cluster */
535  /* TODO: Depending on cluster theta, deltaR changes (currently -1) */
536 
537  SvtxTrack_FastSim* the_track = _trackproj->get_best_track(trackmap, cluster);
538 
539  /* Test if a best track was found */
540  if(the_track!=NULL) //Fill in reconstructed info
541  {
542  /* Get track position and momentum */
543 
544  SvtxTrackState *the_state = _trackproj->get_best_state(the_track, cluster);
545  double posv[3] = {0.,0.,0.};
546  posv[0] = the_state->get_x();
547  posv[1] = the_state->get_y();
548  posv[2] = the_state->get_z();
549  /* Extrapolate track2cluster values */
550  float track2cluster_theta = atan(sqrt(posv[0]*posv[0]+
551  posv[1]*posv[1]) / posv[2]);
552  float track2cluster_eta = -log(abs(tan(track2cluster_theta/2)));
553  if(tan(track2cluster_theta/2)<0)
554  track2cluster_eta*=-1;
555  float track2cluster_phi = atan(posv[1]/posv[0]);
556  float track2cluster_x = posv[0];
557  float track2cluster_y = posv[1];
558  float track2cluster_z = posv[2];
559 
560  /* Fills in weird values right now */
561  float track2cluster_p = 0;
562 
563 
564  tc->set_property( PidCandidate::em_track_theta2cluster, track2cluster_theta );
565  tc->set_property( PidCandidate::em_track_eta2cluster, track2cluster_eta );
566  tc->set_property( PidCandidate::em_track_phi2cluster, track2cluster_phi );
567  tc->set_property( PidCandidate::em_track_p2cluster, track2cluster_p );
568  tc->set_property( PidCandidate::em_track_x2cluster, track2cluster_x );
569  tc->set_property( PidCandidate::em_track_y2cluster, track2cluster_y);
570  tc->set_property( PidCandidate::em_track_z2cluster, track2cluster_z);
571 
572 
573  tc->set_property( PidCandidate::em_track_id, (uint)primary->get_track_id() );
574  tc->set_property( PidCandidate::em_track_quality, (float)100.0 );
575  tc->set_property( PidCandidate::em_track_theta,2*atan( exp(-the_track->get_eta()) ));
576  tc->set_property( PidCandidate::em_track_phi, the_track->get_phi() );
577  tc->set_property( PidCandidate::em_track_eta, the_track->get_eta() );
578  tc->set_property( PidCandidate::em_track_ptotal, the_track->get_p() );
579  tc->set_property( PidCandidate::em_track_ptrans, the_track->get_pt() );
580  tc->set_property( PidCandidate::em_track_charge, the_track->get_charge() );
581  tc->set_property( PidCandidate::em_track_dca, NAN );
582  tc->set_property( PidCandidate::em_track_section, (uint)0 );
583  tc->set_property( PidCandidate::em_track_e3x3_cemc, NAN );
584  tc->set_property( PidCandidate::em_track_e3x3_ihcal, NAN );
585  tc->set_property( PidCandidate::em_track_e3x3_ohcal, NAN );
586  tc->set_property( PidCandidate::em_track_cluster_dr, (float)0.0 );
587 
588  tc->set_property( PidCandidate::em_pid_prob_electron, (float)0.0 );
589  tc->set_property( PidCandidate::em_pid_prob_pion, (float)0.0 );
590  tc->set_property( PidCandidate::em_pid_prob_kaon, (float)0.0 );
591  tc->set_property( PidCandidate::em_pid_prob_proton, (float)0.0 );
592 
593 
594  /*
595  cout << "------------------------------------------" << endl;
596  cout << "Track Projection Successful!" << endl;
597  cout << "Cluster Detector = " << _trackproj->get_detector() << "EMC" << endl;
598  cout << "Cluster X: " << cluster->get_x() << " | Track X: " << posv[0] << endl;
599  cout << "Cluster Y: " << cluster->get_y() << " | Track Y: " << posv[1] << endl;
600  cout << "Cluster Z: " << cluster->get_z() << " | Track Z: " << posv[2] << endl;
601  */
602 
603 
604  }
605  else if(fill_in_truth) //Fill in truth
606  {
607  tc->set_property( PidCandidate::em_track_theta2cluster, NAN);
608  tc->set_property( PidCandidate::em_track_eta2cluster, NAN);
609  tc->set_property( PidCandidate::em_track_phi2cluster, NAN );
610  tc->set_property( PidCandidate::em_track_p2cluster, NAN);
611  tc->set_property( PidCandidate::em_track_x2cluster, NAN );
612  tc->set_property( PidCandidate::em_track_y2cluster, NAN);
613  tc->set_property( PidCandidate::em_track_z2cluster, NAN);
614  tc->set_property( PidCandidate::em_track_id, (uint)primary->get_track_id() );
615  tc->set_property( PidCandidate::em_track_quality, (float)100.0 );
616  tc->set_property( PidCandidate::em_track_theta, gtheta );
617  tc->set_property( PidCandidate::em_track_phi, gphi );
618  tc->set_property( PidCandidate::em_track_eta, geta );
619  tc->set_property( PidCandidate::em_track_ptotal, gptotal );
620  tc->set_property( PidCandidate::em_track_ptrans, gpt );
621  tc->set_property( PidCandidate::em_track_charge, gcharge );
622  tc->set_property( PidCandidate::em_track_dca, NAN );
623  tc->set_property( PidCandidate::em_track_section, (uint)0 );
624  tc->set_property( PidCandidate::em_track_e3x3_cemc, NAN );
625  tc->set_property( PidCandidate::em_track_e3x3_ihcal, NAN );
626  tc->set_property( PidCandidate::em_track_e3x3_ohcal, NAN );
627  tc->set_property( PidCandidate::em_track_cluster_dr, (float)0.0 );
628 
629  /* Get information on track from PID detectors */
630  tc->set_property( PidCandidate::em_pid_prob_electron, (float)0.0 );
631  tc->set_property( PidCandidate::em_pid_prob_pion, (float)0.0 );
632  tc->set_property( PidCandidate::em_pid_prob_kaon, (float)0.0 );
633  tc->set_property( PidCandidate::em_pid_prob_proton, (float)0.0 );
634 
635  if ( abs( primary->get_pid() ) == 11 )
636  tc->set_property( PidCandidate::em_pid_prob_electron, (float)1.0 );
637 
638  if ( abs( primary->get_pid() ) == 211 )
639  tc->set_property( PidCandidate::em_pid_prob_pion, (float)1.0 );
640 
641  if ( abs( primary->get_pid() ) == 321 )
642  tc->set_property( PidCandidate::em_pid_prob_kaon, (float)1.0 );
643 
644  if ( abs( primary->get_pid() ) == 2212 )
645  tc->set_property( PidCandidate::em_pid_prob_proton, (float)1.0 );
646  }
647  else
648  {
649  tc->set_property( PidCandidate::em_track_theta2cluster, NAN);
650  tc->set_property( PidCandidate::em_track_eta2cluster, NAN);
651  tc->set_property( PidCandidate::em_track_phi2cluster, NAN );
652  tc->set_property( PidCandidate::em_track_p2cluster, NAN);
653  tc->set_property( PidCandidate::em_track_x2cluster, NAN );
654  tc->set_property( PidCandidate::em_track_y2cluster, NAN);
655  tc->set_property( PidCandidate::em_track_z2cluster, NAN);
656  tc->set_property( PidCandidate::em_track_id,(uint)0);
657  tc->set_property( PidCandidate::em_track_quality, (float)0.0 );
658  tc->set_property( PidCandidate::em_track_theta, NAN );
659  tc->set_property( PidCandidate::em_track_phi, NAN );
660  tc->set_property( PidCandidate::em_track_eta, NAN );
661  tc->set_property( PidCandidate::em_track_ptotal, NAN );
662  tc->set_property( PidCandidate::em_track_ptrans, NAN );
663  tc->set_property( PidCandidate::em_track_charge, 0 );
664  tc->set_property( PidCandidate::em_track_dca, NAN );
665  tc->set_property( PidCandidate::em_track_section, (uint)0 );
666  tc->set_property( PidCandidate::em_track_e3x3_cemc, NAN );
667  tc->set_property( PidCandidate::em_track_e3x3_ihcal, NAN );
668  tc->set_property( PidCandidate::em_track_e3x3_ohcal, NAN );
669  tc->set_property( PidCandidate::em_track_cluster_dr, (float)0.0 );
670 
671  /* Get information on track from PID detectors */
672  tc->set_property( PidCandidate::em_pid_prob_electron, (float)0.0 );
673  tc->set_property( PidCandidate::em_pid_prob_pion, (float)0.0 );
674  tc->set_property( PidCandidate::em_pid_prob_kaon, (float)0.0 );
675  tc->set_property( PidCandidate::em_pid_prob_proton, (float)0.0 );
676 
677  }
678  }
679 
680  /* add tau candidate to collection */
681  candidateMap.insert( make_pair( tc->get_candidate_id(), tc ) );
682  return 0;
683 }
684 
685 
686 int
687 DISKinematicsReco::AddGlobalCalorimeterInformation()
688 {
689  /* Missing transverse energy, transverse energy direction, and Energy sums, x and y components separately */
690  float Et_miss = 0;
691  float Et_miss_phi = 0;
692  float Ex_sum = 0;
693  float Ey_sum = 0;
694 
695  /* energy threshold for considering tower */
696  float tower_emin = 0.0;
697 
698  /* Loop over all tower (and geometry) collections */
699  vector < string > calo_names;
700  calo_names.push_back( "CEMC" );
701  calo_names.push_back( "HCALIN" );
702  calo_names.push_back( "HCALOUT" );
703  calo_names.push_back( "FEMC" );
704  calo_names.push_back( "FHCAL" );
705  calo_names.push_back( "EEMC" );
706 
707  for ( unsigned i = 0; i < calo_names.size(); i++ )
708  {
709 
710  string towernodename = "TOWER_CALIB_" + calo_names.at( i );
711  RawTowerContainer* towers = findNode::getClass<RawTowerContainer>(_topNode, towernodename.c_str());
712  if (!towers)
713  {
714  std::cout << PHWHERE << ": Could not find node " << towernodename.c_str() << std::endl;
715  return -1;
716  }
717 
718  string towergeomnodename = "TOWERGEOM_" + calo_names.at( i );
719  RawTowerGeomContainer *towergeoms = findNode::getClass<RawTowerGeomContainer>(_topNode, towergeomnodename.c_str());
720  if (! towergeoms)
721  {
722  cout << PHWHERE << ": Could not find node " << towergeomnodename.c_str() << endl;
723  return -1;
724  }
725 
726  /* define tower iterator */
727  RawTowerContainer::ConstRange begin_end = towers->getTowers();
728  RawTowerContainer::ConstIterator rtiter;
729 
730  /* loop over all tower in CEMC calorimeter */
731  for (rtiter = begin_end.first; rtiter != begin_end.second; ++rtiter)
732  {
733  /* get tower energy */
734  RawTower *tower = rtiter->second;
735  float tower_energy = tower->get_energy();
736 
737  /* check if tower above energy treshold */
738  if ( tower_energy < tower_emin )
739  continue;
740 
741  /* get eta and phi of tower and check angle delta_R w.r.t. jet axis */
742  RawTowerGeom * tower_geom = towergeoms->get_tower_geometry(tower -> get_key());
743  float tower_eta = tower_geom->get_eta();
744  float tower_phi = tower_geom->get_phi();
745 
746  /* from https://en.wikipedia.org/wiki/Pseudorapidity:
747  p_x = p_T * cos( phi )
748  p_y = p_T * sin( phi )
749  p_z = p_T * sinh( eta )
750  |p| = p_T * cosh( eta )
751  */
752 
753  /* calculate 'transverse' tower energy */
754  float tower_energy_t = tower_energy / cosh( tower_eta );
755 
756  /* add energy components of this tower to total energy components */
757  Ex_sum += tower_energy_t * cos( tower_phi );
758  Ey_sum += tower_energy_t * sin( tower_phi );
759  }
760  }
761 
762  /* calculate Et_miss */
763  Et_miss = sqrt( Ex_sum * Ex_sum + Ey_sum * Ey_sum );
764  Et_miss_phi = atan( Ey_sum / Ex_sum );
765 
766  /* set event branch paraemters */
767  ( _map_event_branches.find( "Et_miss" ) )->second = Et_miss;
768  ( _map_event_branches.find( "Et_miss_phi" ) )->second = Et_miss_phi;
769 
770  return 0;
771 }
772 
773 int
774 DISKinematicsReco::WriteCandidatesToTree( type_map_tcan& electronCandidateMap )
775 {
776  /* Get number of electron candidates */
777  ( _map_event_branches.find( "em_ncandidates" ) )->second = electronCandidateMap.size();
778 
779  /* Loop over all EM candidates and add them to tree */
780  for (type_map_tcan::iterator iter = electronCandidateMap.begin();
781  iter != electronCandidateMap.end();
782  ++iter)
783  {
784  /* update information in map and fill tree */
785  for ( map< PidCandidate::PROPERTY , vector<float> >::iterator iter_prop = _map_em_candidate_branches.begin();
786  iter_prop != _map_em_candidate_branches.end();
787  ++iter_prop)
788  {
789  switch ( PidCandidate::get_property_info( (iter_prop->first) ).second ) {
790 
791  case PidCandidate::type_float :
792  (iter_prop->second).push_back( (iter->second)->get_property_float( (iter_prop->first) ) );
793  break;
794 
795  case PidCandidate::type_int :
796  (iter_prop->second).push_back( (iter->second)->get_property_int( (iter_prop->first) ) );
797  break;
798 
799  case PidCandidate::type_uint :
800  (iter_prop->second).push_back( (iter->second)->get_property_uint( (iter_prop->first) ) );
801  break;
802 
803  case PidCandidate::type_unknown :
804  break;
805  }
806  }
807  }
808 
809  return 0;
810 }
811 
812 
813 PidCandidate*
814 DISKinematicsReco::FindMinDeltaRCandidate( type_map_tcan *candidates, const float eta_ref, const float phi_ref )
815 {
816  /* PidCandidate with eta, phi closest to reference */
817  PidCandidate* best_candidate = NULL;
818 
819  return best_candidate;
820 }
821 
822 
823 float
824 DISKinematicsReco::CalculateDeltaR( float eta1, float phi1, float eta2, float phi2 )
825 {
826  /* Particles at phi = PI+x and phi = PI-x are actually close to each other in phi, but simply calculating
827  * the difference in phi would give a large distance (because phi ranges from -PI to +PI in the convention
828  * used. Account for this by subtracting 2PI is particles fall within this border area. */
829  if( ( phi1 < -0.9*TMath::Pi()) && ( phi2 > 0.9*TMath::Pi() ) ) phi2 = phi2 - 2*TMath::Pi();
830  if( ( phi1 > 0.9*TMath::Pi()) && ( phi2 < -0.9*TMath::Pi() ) ) phi1 = phi1 - 2*TMath::Pi();
831 
832  float delta_R = sqrt( pow( eta2 - eta1, 2 ) + pow( phi2 - phi1, 2 ) );
833 
834  return delta_R;
835 }
836 
837 int
838 DISKinematicsReco::AddReconstructedKinematics( type_map_tcan& em_candidates , string mode )
839 {
840  /* Loop over all EM candidates */
841  for (type_map_tcan::iterator iter = em_candidates.begin();
842  iter != em_candidates.end();
843  ++iter)
844  {
845  /* Add reco kinematics based on scattered electron data */
846  PidCandidate* the_electron = iter->second;
847 
848  float e0_E = _beam_electron_ptotal;
849  float p0_E = _beam_hadron_ptotal;
850 
851  /* get scattered particle kinematicsbased on chosen mode */
852  float e1_E = NAN;
853  float e1_theta = NAN;
854 
855  if ( mode == "cluster" )
856  {
857  e1_E = the_electron->get_property_float( PidCandidate::em_cluster_e );
858  e1_theta = the_electron->get_property_float( PidCandidate::em_cluster_theta );
859  }
860  else if ( mode == "truth" )
861  {
862  e1_E = the_electron->get_property_float( PidCandidate::em_evtgen_ptotal );
863  e1_theta = the_electron->get_property_float( PidCandidate::em_evtgen_theta );
864  }
865  else
866  {
867  cout << "WARNING: Unknown mode " << mode << " selected." << endl;
868  return -1;
869  }
870 
871  /* for purpose of calculations, 'theta' angle of the scattered electron is defined as angle
872  between 'scattered electron' and 'direction of incoming electron'. Since initial electron
873  has 'theta = Pi' in coordinate system, we need to use 'theta_rel = Pi - theta' for calculations
874  of kinematics. */
875  float e1_theta_rel = M_PI - e1_theta;
876 
877  /* event kinematics */
878  float dis_s = 4.0 * e0_E * p0_E;
879 
880  float dis_Q2 = 2.0 * e0_E * e1_E * ( 1 - cos( e1_theta_rel ) );
881 
882  /* ePHENIX LOI definition of y: */
883  float dis_y = 1.0 - ( e1_E / e0_E ) + ( dis_Q2 / ( 4.0 * e0_E * e0_E ) );
884 
885  /* G. Wolf Hera Physics definitions of y: */
886  //float dis_y = 1 - ( e1_E / (2*e0_E) ) * ( 1 - cos( e1_theta_rel ) );
887 
888  float dis_x = dis_Q2 / ( dis_s * dis_y );
889 
890  float dis_W2 = _mproton*_mproton + dis_Q2 * ( ( 1 / dis_x ) - 1 );
891  float dis_W = sqrt( dis_W2 );
892 
893  the_electron->set_property( PidCandidate::em_reco_x_e, dis_x );
894  the_electron->set_property( PidCandidate::em_reco_y_e, dis_y );
895  the_electron->set_property( PidCandidate::em_reco_q2_e, dis_Q2 );
896  the_electron->set_property( PidCandidate::em_reco_w_e, dis_W );
897  }
898 
899  return 0;
900 }
901 
902 int
903 DISKinematicsReco::AddTruthEventInformation()
904 {
905  /* Get collection of truth particles from event generator */
906  PHHepMCGenEventMap *geneventmap = findNode::getClass<PHHepMCGenEventMap>(_topNode,"PHHepMCGenEventMap");
907  if (!geneventmap) {
908  std::cout << PHWHERE << " WARNING: Can't find requested PHHepMCGenEventMap" << endl;
909  return -1;
910  }
911 
912  /* Add truth kinematics */
913  int embedding_id = 1;
914  PHHepMCGenEvent *genevt = geneventmap->get(embedding_id);
915  if (!genevt)
916  {
917  std::cout << PHWHERE << "WARNING: Node PHHepMCGenEventMap missing subevent with embedding ID "<< embedding_id;
918  std::cout <<". Print PHHepMCGenEventMap:";
919  geneventmap->identify();
920  return -1;
921  }
922 
923  HepMC::GenEvent* theEvent = genevt->getEvent();
924 
925  if ( !theEvent )
926  {
927  std::cout << PHWHERE << "WARNING: Missing requested GenEvent!" << endl;
928  return -1;
929  }
930 
931  int true_process_id = theEvent->signal_process_id();
932  float ev_x = -NAN;
933  float ev_Q2 = -NAN;
934  float ev_W = -NAN;
935  float ev_y = -NAN;
936 
937  float ev_s = 4 * _beam_electron_ptotal * _beam_hadron_ptotal;
938 
939 
940  if ( theEvent->pdf_info() )
941  {
942  //ev_x = theEvent->pdf_info()->x1();
943  ev_x = theEvent->pdf_info()->x2();
944  ev_Q2 = theEvent->pdf_info()->scalePDF();
945  ev_y = ev_Q2 / ( ev_x * ev_s );
946  ev_W = sqrt( _mproton*_mproton + ev_Q2 * ( 1 / ev_x - 1 ) );
947  }
948 
949  ( _map_event_branches.find( "evtgen_process_id" ) )->second = true_process_id;
950  ( _map_event_branches.find( "evtgen_s" ) )->second = ev_s;
951  ( _map_event_branches.find( "evtgen_x" ) )->second = ev_x;
952  ( _map_event_branches.find( "evtgen_Q2" ) )->second = ev_Q2;
953  ( _map_event_branches.find( "evtgen_y" ) )->second = ev_y;
954  ( _map_event_branches.find( "evtgen_W" ) )->second = ev_W;
955 
956  return 0;
957 }
958 
959 void
960 DISKinematicsReco::ResetBranchMap()
961 {
962  /* Event branches */
963  for ( map< string , float >::iterator iter = _map_event_branches.begin();
964  iter != _map_event_branches.end();
965  ++iter)
966  {
967  (iter->second) = NAN;
968  }
969 
970  /* EM candidate branches */
971  for ( map< PidCandidate::PROPERTY , vector<float> >::iterator iter = _map_em_candidate_branches.begin();
972  iter != _map_em_candidate_branches.end();
973  ++iter)
974  {
975  (iter->second).clear();
976  }
977 
978  return;
979 }
980 
981 int
982 DISKinematicsReco::End(PHCompositeNode *topNode)
983 {
984  _tfile->cd();
985 
986  if ( _tree_event_cluster )
987  _tree_event_cluster->Write();
988 
989  if ( _tree_event_truth )
990  _tree_event_truth->Write();
991 
992  _tfile->Close();
993 
994  return 0;
995 }