analysis/EICAnalysis/DISKinematicsReco.C

0001 #include "DISKinematicsReco.h"
0002 #include "PidCandidatev1.h"
0003 #include "TruthTrackerHepMC.h"
0004 #include "TrackProjectionTools.h"
0005 #include "TrackProjectorPlaneECAL.h"
0006 /* STL includes */
0007 #include <cassert>
0008
0009 /* Fun4All includes */
0010 #include <trackbase_historic/SvtxTrack.h>
0011 #include <trackbase_historic/SvtxTrackMap.h>
0012 #include <trackbase_historic/SvtxTrack_FastSim.h>
0013
0014 #include <trackbase_historic/SvtxTrackMap_v1.h>
0015
0016 #include <phool/getClass.h>
0017
0018 #include <fun4all/Fun4AllReturnCodes.h>
0019
0020 #include <g4jets/JetMap.h>
0021
0022 #include <calobase/RawTowerGeomContainer.h>
0023 #include <calobase/RawTowerContainer.h>
0024 #include <calobase/RawTowerGeom.h>
0025 #include <calobase/RawTowerv1.h>
0026 #include <calobase/RawTowerDefs.h>
0027
0028 #include <calobase/RawClusterContainer.h>
0029 #include <calobase/RawCluster.h>
0030
0031 #include <g4vertex/GlobalVertexMap.h>
0032 #include <g4vertex/GlobalVertex.h>
0033
0034 #include <g4main/PHG4Particle.h>
0035
0036 #include <g4eval/CaloEvalStack.h>
0037 #include <g4eval/CaloRawTowerEval.h>
0038
0039 /* ROOT includes */
0040 #include <TLorentzVector.h>
0041 #include <TString.h>
0042 #include <TNtuple.h>
0043 #include <TTree.h>
0044 #include <TFile.h>
0045 #include <TDatabasePDG.h>
0046
0047 using namespace std;
0048
0049 DISKinematicsReco::DISKinematicsReco(std::string filename) :
0050   SubsysReco("DISKinematicsReco" ),
0051   _mproton( 0.938272 ),
0052   _save_towers(false),
0053   _save_tracks(false),
0054   _do_process_geant4_cluster(false),
0055   _do_process_truth(false),
0056   _ievent(0),
0057   _filename(filename),
0058   _tfile(nullptr),
0059   _tree_event_cluster(nullptr),
0060   _tree_event_truth(nullptr),
0061   _beam_electron_ptotal(10),
0062   _beam_hadron_ptotal(250),
0063   _trackproj(nullptr)
0064 {
0065
0066 }
0067
0068 int
0069 DISKinematicsReco::Init(PHCompositeNode *topNode)
0070 {
0071   _topNode = topNode;
0072
0073   _ievent = 0;
0074
0075   _tfile = new TFile(_filename.c_str(), "RECREATE");
0076
0077   /* Add PidCandidate properties to map that defines output tree */
0078   float dummy = 0;
0079   vector< float > vdummy;
0080
0081   _map_em_candidate_branches.insert( make_pair( PidCandidate::em_cluster_id , vdummy ) );
0082   _map_em_candidate_branches.insert( make_pair( PidCandidate::em_cluster_prob , vdummy ) );
0083   _map_em_candidate_branches.insert( make_pair( PidCandidate::em_cluster_posx , vdummy ) );
0084   _map_em_candidate_branches.insert( make_pair( PidCandidate::em_cluster_posy , vdummy ) );
0085   _map_em_candidate_branches.insert( make_pair( PidCandidate::em_cluster_posz , vdummy ) );
0086   _map_em_candidate_branches.insert( make_pair( PidCandidate::em_cluster_e , vdummy ) );
0087   _map_em_candidate_branches.insert( make_pair( PidCandidate::em_cluster_ecore , vdummy ) );
0088   _map_em_candidate_branches.insert( make_pair( PidCandidate::em_cluster_et_iso , vdummy ) );
0089   _map_em_candidate_branches.insert( make_pair( PidCandidate::em_cluster_theta , vdummy ) );
0090   _map_em_candidate_branches.insert( make_pair( PidCandidate::em_cluster_phi , vdummy ) );
0091   _map_em_candidate_branches.insert( make_pair( PidCandidate::em_cluster_eta , vdummy ) );
0092   _map_em_candidate_branches.insert( make_pair( PidCandidate::em_cluster_pt , vdummy ) );
0093   _map_em_candidate_branches.insert( make_pair( PidCandidate::em_cluster_ntower , vdummy ) );
0094   _map_em_candidate_branches.insert( make_pair( PidCandidate::em_cluster_caloid , vdummy ) );
0095
0096   _map_em_candidate_branches.insert( make_pair( PidCandidate::em_track_id , vdummy ) );
0097   _map_em_candidate_branches.insert( make_pair( PidCandidate::em_track_quality , vdummy ) );
0098   _map_em_candidate_branches.insert( make_pair( PidCandidate::em_track_theta , vdummy ) );
0099   _map_em_candidate_branches.insert( make_pair( PidCandidate::em_track_phi , vdummy ) );
0100   _map_em_candidate_branches.insert( make_pair( PidCandidate::em_track_eta , vdummy ) );
0101   _map_em_candidate_branches.insert( make_pair( PidCandidate::em_track_ptotal , vdummy ) );
0102   _map_em_candidate_branches.insert( make_pair( PidCandidate::em_track_ptrans , vdummy ) );
0103   _map_em_candidate_branches.insert( make_pair( PidCandidate::em_track_charge , vdummy ) );
0104   _map_em_candidate_branches.insert( make_pair( PidCandidate::em_track_dca , vdummy ) );
0105   _map_em_candidate_branches.insert( make_pair( PidCandidate::em_track_section , vdummy ) );
0106   _map_em_candidate_branches.insert( make_pair( PidCandidate::em_track_e3x3_cemc , vdummy ) );
0107   _map_em_candidate_branches.insert( make_pair( PidCandidate::em_track_e3x3_femc , vdummy ) );
0108   _map_em_candidate_branches.insert( make_pair( PidCandidate::em_track_e3x3_eemc , vdummy ) );
0109   _map_em_candidate_branches.insert( make_pair( PidCandidate::em_track_e3x3_ihcal , vdummy ) );
0110   _map_em_candidate_branches.insert( make_pair( PidCandidate::em_track_e3x3_ohcal , vdummy ) );
0111   _map_em_candidate_branches.insert( make_pair( PidCandidate::em_track_e3x3_fhcal , vdummy ) );
0112   _map_em_candidate_branches.insert( make_pair( PidCandidate::em_track_e3x3_ehcal , vdummy ) );
0113   _map_em_candidate_branches.insert( make_pair( PidCandidate::em_track_cluster_dr , vdummy ) );
0114
0115    _map_em_candidate_branches.insert( make_pair( PidCandidate::em_track_theta2cluster , vdummy ) );
0116    _map_em_candidate_branches.insert( make_pair( PidCandidate::em_track_eta2cluster , vdummy ) );
0117   _map_em_candidate_branches.insert( make_pair( PidCandidate::em_track_phi2cluster , vdummy ) );
0118   _map_em_candidate_branches.insert( make_pair( PidCandidate::em_track_p2cluster , vdummy ) );
0119   _map_em_candidate_branches.insert( make_pair( PidCandidate::em_track_x2cluster , vdummy ) );
0120   _map_em_candidate_branches.insert( make_pair( PidCandidate::em_track_y2cluster , vdummy ) );
0121   _map_em_candidate_branches.insert( make_pair( PidCandidate::em_track_z2cluster , vdummy ) );
0122
0123   _map_em_candidate_branches.insert( make_pair( PidCandidate::em_pid_prob_electron , vdummy ) );
0124   _map_em_candidate_branches.insert( make_pair( PidCandidate::em_pid_prob_pion , vdummy ) );
0125   _map_em_candidate_branches.insert( make_pair( PidCandidate::em_pid_prob_kaon , vdummy ) );
0126   _map_em_candidate_branches.insert( make_pair( PidCandidate::em_pid_prob_proton , vdummy ) );
0127
0128   _map_em_candidate_branches.insert( make_pair( PidCandidate::em_evtgen_pid , vdummy ) );
0129   _map_em_candidate_branches.insert( make_pair( PidCandidate::em_evtgen_ptotal , vdummy ) );
0130   _map_em_candidate_branches.insert( make_pair( PidCandidate::em_evtgen_theta , vdummy ) );
0131   _map_em_candidate_branches.insert( make_pair( PidCandidate::em_evtgen_phi , vdummy ) );
0132   _map_em_candidate_branches.insert( make_pair( PidCandidate::em_evtgen_eta , vdummy ) );
0133   _map_em_candidate_branches.insert( make_pair( PidCandidate::em_evtgen_charge , vdummy ) );
0134   _map_em_candidate_branches.insert( make_pair( PidCandidate::em_evtgen_is_scattered_lepton , vdummy ) );
0135
0136   _map_em_candidate_branches.insert( make_pair( PidCandidate::em_reco_x_e , vdummy ) );
0137   _map_em_candidate_branches.insert( make_pair( PidCandidate::em_reco_y_e , vdummy ) );
0138   _map_em_candidate_branches.insert( make_pair( PidCandidate::em_reco_q2_e , vdummy ) );
0139   _map_em_candidate_branches.insert( make_pair( PidCandidate::em_reco_w_e , vdummy ) );
0140
0141   /* Add branches to map that defines output tree for event-wise properties */
0142   /*
0143     runnr
0144     eventnr
0145     eventweight
0146     cal_px
0147     cal_py
0148     cal_pz
0149     cal_etotal
0150     cal_et
0151     cal_empz
0152     cal_pt --> "pt_miss" ?
0153     cal_pt_phi = atan2(sum(py)/sum(px)) --> "pt_miss_angle" ?
0154   */
0155   _map_event_branches.insert( make_pair( "em_ncandidates" , dummy ) );
0156
0157   _map_event_branches.insert( make_pair( "Et_miss" , dummy ) );
0158   _map_event_branches.insert( make_pair( "Et_miss_phi" , dummy ) );
0159
0160   _map_event_branches.insert( make_pair( "evtgen_process_id" , dummy ) );
0161   _map_event_branches.insert( make_pair( "evtgen_s" , dummy ) );
0162   _map_event_branches.insert( make_pair( "evtgen_x" , dummy ) );
0163   _map_event_branches.insert( make_pair( "evtgen_Q2" , dummy ) );
0164   _map_event_branches.insert( make_pair( "evtgen_y" , dummy ) );
0165   _map_event_branches.insert( make_pair( "evtgen_W" , dummy ) );
0166
0167   /* Create tree for information about full event */
0168   _tree_event_cluster = new TTree("event_cluster", "a Tree with global event information and EM candidates");
0169   _tree_event_cluster->Branch("event", &_ievent, "event/I");
0170
0171   /* Add event branches */
0172   for ( map< string , float >::iterator iter = _map_event_branches.begin();
0173         iter != _map_event_branches.end();
0174         ++iter)
0175     {
0176       _tree_event_cluster->Branch( (iter->first).c_str(),
0177                                    &(iter->second) );
0178     }
0179
0180   /* Add EM candidate branches */
0181   for ( map< PidCandidate::PROPERTY , vector<float> >::iterator iter = _map_em_candidate_branches.begin();
0182         iter != _map_em_candidate_branches.end();
0183         ++iter)
0184     {
0185       _tree_event_cluster->Branch( PidCandidate::get_property_info( (iter->first) ).first.c_str(),
0186                                    &(iter->second) );
0187     }
0188
0189   /* create clone of tree for truth particle candidates */
0190   _tree_event_truth = _tree_event_cluster->CloneTree();
0191   _tree_event_truth->SetName("event_truth");
0192   _tree_event_truth->SetTitle("a Tree with global event information and truth particle candidates");
0193
0194
0195   return 0;
0196 }
0197
0198 int
0199 DISKinematicsReco::InitRun(PHCompositeNode *topNode)
0200 {
0201   if(_do_process_geant4_cluster)
0202     _trackproj=new TrackProjectorPlaneECAL( topNode );
0203    return 0;
0204 }
0205
0206 int
0207 DISKinematicsReco::process_event(PHCompositeNode *topNode)
0208 {
0209   /* Find electron candidates based on calorimeter cluster from full Geant4 detector simulation */
0210   if ( _do_process_geant4_cluster )
0211     {
0212       /* Reset branch map */
0213       ResetBranchMap();
0214
0215       /* Create map to collect electron candidates.
0216        * Use energy as 'key' to the map because energy is unique for each jet, while there are sometimes multiple jets (same energy,
0217        * different jet ID) in the input jet collection. Also, map automatically sorts entries by key, i.e. this gives list of tau candidates
0218        * sorted by energy. */
0219       type_map_tcan electronCandidateMap;
0220
0221       /* Collect EM candidates from calorimeter cluster */
0222       CollectEmCandidatesFromCluster( electronCandidateMap );
0223
0224       /* Calculate kinematics for each em candidate */
0225       AddReconstructedKinematics( electronCandidateMap, "cluster" );
0226
0227       /* Add information about em candidats to output tree */
0228       WriteCandidatesToTree( electronCandidateMap );
0229
0230       /* Add global event information */
0231       AddGlobalCalorimeterInformation();
0232       AddTruthEventInformation();
0233
0234       /* fill event information tree */
0235       _tree_event_cluster->Fill();
0236     }
0237
0238   /* Find electron candidates based on truth particle inforamtion */
0239   if ( _do_process_truth )
0240     {
0241       /* reset branches, fill with truth particle candidates, and fill different tree */
0242       ResetBranchMap();
0243
0244       /* Collect electron candidates based on TRUTH information */
0245       type_map_tcan electronTruthCandidateMap;
0246
0247       CollectEmCandidatesFromTruth( electronTruthCandidateMap );
0248       AddReconstructedKinematics( electronTruthCandidateMap, "truth" );
0249       WriteCandidatesToTree( electronTruthCandidateMap );
0250       AddTruthEventInformation();
0251
0252       /* fill event information tree */
0253       _tree_event_truth->Fill();
0254     }
0255
0256   /* count up event number */
0257   _ievent ++;
0258
0259   return 0;
0260 }
0261
0262
0263
0264 int
0265 DISKinematicsReco::CollectEmCandidatesFromCluster( type_map_tcan& electronCandidateMap )
0266 {
0267
0268   /* Loop over all clusters in EEMC, CEMC, and FEMC to find electron candidates */
0269   vector< string > v_ecals;
0270   v_ecals.push_back("EEMC");
0271   v_ecals.push_back("CEMC");
0272   v_ecals.push_back("FEMC");
0273   for ( unsigned idx = 0; idx < v_ecals.size(); idx++ )
0274     {
0275       CaloEvalStack * caloevalstack = new CaloEvalStack(_topNode, v_ecals.at( idx ) );
0276
0277       string clusternodename = "CLUSTER_" + v_ecals.at( idx );
0278       RawClusterContainer *clusterList = findNode::getClass<RawClusterContainer>(_topNode,clusternodename.c_str());
0279       SvtxTrackMap *trackmap = findNode::getClass<SvtxTrackMap>(_topNode,"TrackMap");
0280
0281       if (!clusterList) {
0282         cerr << PHWHERE << " ERROR: Can't find node " << clusternodename << endl;
0283         return false;
0284       }
0285       if(!trackmap) {
0286     cerr << PHWHERE << " ERROR: Can't find node SvtxTrackMap" << endl;
0287     return false;
0288       }
0289
0290       for (unsigned int k = 0; k < clusterList->size(); ++k)
0291         {
0292           RawCluster *cluster = clusterList->getCluster(k);
0293           /* Check if cluster energy is below threshold */
0294           float e_cluster_threshold = 0.3;
0295           if ( cluster->get_energy() < e_cluster_threshold )
0296             continue;
0297
0298           InsertCandidateFromCluster( electronCandidateMap , cluster , caloevalstack , trackmap);
0299         }
0300     }
0301
0302   return 0;
0303 }
0304
0305
0306 int
0307 DISKinematicsReco::CollectEmCandidatesFromTruth( type_map_tcan& candidateMap )
0308 {
0309   /* Get collection of truth particles from event generator */
0310   PHHepMCGenEventMap *geneventmap = findNode::getClass<PHHepMCGenEventMap>(_topNode,"PHHepMCGenEventMap");
0311   if (!geneventmap) {
0312     //std::cout << PHWHERE << " WARNING: Can't find requested PHHepMCGenEventMap" << endl;
0313     return -1;
0314   }
0315
0316   /* Add truth kinematics */
0317   int embedding_id = 1;
0318   PHHepMCGenEvent *genevt = geneventmap->get(embedding_id);
0319   if (!genevt)
0320     {
0321       std::cout << PHWHERE << "WARNING: Node PHHepMCGenEventMap missing subevent with embedding ID "<< embedding_id;
0322       std::cout <<". Print PHHepMCGenEventMap:";
0323       geneventmap->identify();
0324       return -1;
0325     }
0326
0327   HepMC::GenEvent* theEvent = genevt->getEvent();
0328
0329   if ( !theEvent )
0330     {
0331       std::cout << PHWHERE << "WARNING: Missing requested GenEvent!" << endl;
0332       return -1;
0333     }
0334
0335   /* Look for beam particles and scattered lepton */
0336   TruthTrackerHepMC truth;
0337   truth.set_hepmc_geneventmap( geneventmap );
0338   //HepMC::GenParticle* particle_beam_l = truth.FindBeamLepton();
0339   //HepMC::GenParticle* particle_beam_h = truth.FindBeamHadron();
0340   HepMC::GenParticle* particle_scattered_l = truth.FindScatteredLepton();
0341   /* loop over all particles */
0342   for (HepMC::GenEvent::particle_const_iterator p = theEvent->particles_begin();
0343        p != theEvent->particles_end(); ++p)
0344     {
0345       TParticlePDG * pdg_p = TDatabasePDG::Instance()->GetParticle( (*p)->pdg_id() );
0346       int charge = -999;
0347       if ( pdg_p )
0348     {
0349       /* NOTE: TParticlePDG::Charge() returns charge in units of |e|/3 (see ROOT documentation) */
0350       charge = pdg_p->Charge() / 3;
0351     }
0352
0353       /* skip particles that are not stable final state particles (status 1) */
0354       if ( (*p)->status() != 1 )
0355         continue;
0356       /* create new pid candidate */
0357       PidCandidatev1 *tc = new PidCandidatev1();
0358       tc->set_candidate_id( candidateMap.size()+1 );
0359
0360       float mom_eta = -1 * log ( tan( (*p)->momentum().theta() / 2.0 ) );
0361
0362       tc->set_property( PidCandidate::em_evtgen_pid, (*p)->pdg_id() );
0363       tc->set_property( PidCandidate::em_evtgen_ptotal, (float) (*p)->momentum().e() );
0364       tc->set_property( PidCandidate::em_evtgen_theta, (float) (*p)->momentum().theta() );
0365       tc->set_property( PidCandidate::em_evtgen_phi, (float) (*p)->momentum().phi() );
0366       tc->set_property( PidCandidate::em_evtgen_eta, mom_eta );
0367       tc->set_property( PidCandidate::em_evtgen_charge, charge );
0368       tc->set_property( PidCandidate::em_evtgen_is_scattered_lepton, (uint)0 );
0369
0370       /* is scattered lepton? */
0371       if ( particle_scattered_l &&
0372        (*p) == particle_scattered_l )
0373     tc->set_property( PidCandidate::em_evtgen_is_scattered_lepton, (uint)1 );
0374
0375       /* add pid candidate to collection */
0376       candidateMap.insert( make_pair( tc->get_candidate_id(), tc ) );
0377     }
0378
0379   return 0;
0380 }
0381
0382
0383 int
0384 DISKinematicsReco::InsertCandidateFromCluster( type_map_tcan& candidateMap , RawCluster *cluster , CaloEvalStack *caloevalstack, SvtxTrackMap *trackmap )
0385 {
0386
0387   /* create new electron candidate */
0388   PidCandidatev1 *tc = new PidCandidatev1();
0389   tc->set_candidate_id( candidateMap.size()+1 );
0390
0391   /* set some initial cluster properties */
0392   float theta = atan2( cluster->get_r() , cluster->get_z() );
0393   float eta =  -1 * log( tan( theta / 2.0 ) );
0394   float pt = cluster->get_energy() * sin( theta );
0395
0396   /* get calorimeter ID where towers of cluster are found */
0397   unsigned caloid = 0;
0398   RawCluster::TowerConstIterator rtiter = cluster->get_towers().first;
0399   caloid = RawTowerDefs::decode_caloid( rtiter->first );
0400   _trackproj->set_detector(_trackproj->get_detector_from_cluster(cluster));
0401   //cout << "Calo ID: " << caloid << " -> " << RawTowerDefs::convert_caloid_to_name( RawTowerDefs::decode_caloid( rtiter->first ) ) << endl;
0402
0403   tc->set_property( PidCandidate::em_cluster_id, cluster->get_id() );
0404   tc->set_property( PidCandidate::em_cluster_prob, cluster->get_prob() );
0405   tc->set_property( PidCandidate::em_cluster_posx, cluster->get_x() );
0406   tc->set_property( PidCandidate::em_cluster_posy, cluster->get_y() );
0407   tc->set_property( PidCandidate::em_cluster_posz, cluster->get_z() );
0408   tc->set_property( PidCandidate::em_cluster_e, cluster->get_energy() );
0409   tc->set_property( PidCandidate::em_cluster_ecore, cluster->get_ecore() );
0410   tc->set_property( PidCandidate::em_cluster_et_iso, cluster->get_et_iso() );
0411   tc->set_property( PidCandidate::em_cluster_theta, theta );
0412   tc->set_property( PidCandidate::em_cluster_phi, cluster->get_phi() );
0413   tc->set_property( PidCandidate::em_cluster_eta, eta );
0414   tc->set_property( PidCandidate::em_cluster_pt, pt );
0415   tc->set_property( PidCandidate::em_cluster_ntower, (unsigned)cluster->getNTowers() );
0416   tc->set_property( PidCandidate::em_cluster_caloid, caloid );
0417
0418   /* get track projection helper class */
0419   TrackProjectionTools tpt( _topNode );
0420
0421   /* find matching reco track */
0422   float best_track_dr = NAN;
0423   SvtxTrack* best_track = NULL; //tpt.FindClosestTrack( cluster, best_track_dr ); /* @TODO switch track finding back on as soon as we understand it better */
0424
0425   /* IF matching track found: set track properties */
0426   if ( best_track )
0427     {
0428       /* set some initial track properties */
0429       float theta = 2.0 * atan( exp( -1 * best_track->get_eta() ) );
0430
0431       tc->set_property( PidCandidate::em_track_id, (uint)best_track->get_id() );
0432       tc->set_property( PidCandidate::em_track_quality, best_track->get_quality() );
0433       tc->set_property( PidCandidate::em_track_theta, theta );
0434       tc->set_property( PidCandidate::em_track_phi, best_track->get_phi() );
0435       tc->set_property( PidCandidate::em_track_eta, best_track->get_eta() );
0436       tc->set_property( PidCandidate::em_track_ptotal, best_track->get_p() );
0437       tc->set_property( PidCandidate::em_track_ptrans, best_track->get_pt() );
0438       tc->set_property( PidCandidate::em_track_charge, best_track->get_charge() );
0439       tc->set_property( PidCandidate::em_track_dca, best_track->get_dca() );
0440       tc->set_property( PidCandidate::em_track_section, (uint)0 );
0441       tc->set_property( PidCandidate::em_track_e3x3_cemc, best_track->get_cal_energy_3x3(SvtxTrack::CEMC) );
0442       tc->set_property( PidCandidate::em_track_e3x3_ihcal, best_track->get_cal_energy_3x3(SvtxTrack::HCALIN) );
0443       tc->set_property( PidCandidate::em_track_e3x3_ohcal, best_track->get_cal_energy_3x3(SvtxTrack::HCALOUT) );
0444       tc->set_property( PidCandidate::em_track_cluster_dr, best_track_dr );
0445
0446       // Use the track states to project to the FEMC / FHCAL / EEMC and generate
0447       // energy sums.
0448       float e3x3_femc = NAN;
0449       float e3x3_fhcal = NAN;
0450       float e3x3_eemc = NAN;
0451
0452       for (SvtxTrack::ConstStateIter state_itr = best_track->begin_states();
0453            state_itr != best_track->end_states(); state_itr++) {
0454
0455         SvtxTrackState *temp = dynamic_cast<SvtxTrackState*>(state_itr->second);
0456
0457         if( (temp->get_name()=="FHCAL") )
0458           e3x3_fhcal = tpt.getE33Forward( "FHCAL" , temp->get_x() , temp->get_y() );
0459
0460         if( (temp->get_name()=="FEMC") )
0461           e3x3_femc = tpt.getE33Forward( "FEMC" , temp->get_x() , temp->get_y() );
0462
0463         if( (temp->get_name()=="EEMC") )
0464           e3x3_eemc = tpt.getE33Forward( "EEMC" , temp->get_x() , temp->get_y() );
0465       }
0466
0467       /* set candidate properties */
0468       tc->set_property( PidCandidate::em_track_e3x3_fhcal, e3x3_fhcal );
0469       tc->set_property( PidCandidate::em_track_e3x3_femc, e3x3_femc );
0470       tc->set_property( PidCandidate::em_track_e3x3_eemc, e3x3_eemc );
0471
0472       /* Get information on track from PID detectors */
0473       tc->set_property( PidCandidate::em_pid_prob_electron, (float)0.0 );
0474       tc->set_property( PidCandidate::em_pid_prob_pion, (float)0.0 );
0475       tc->set_property( PidCandidate::em_pid_prob_kaon, (float)0.0 );
0476       tc->set_property( PidCandidate::em_pid_prob_proton, (float)0.0 );
0477     }
0478
0479   /* set em candidate MC truth properties */
0480   tc->set_property( PidCandidate::em_evtgen_pid, (int)NAN );
0481   tc->set_property( PidCandidate::em_evtgen_ptotal, (float)NAN );
0482   tc->set_property( PidCandidate::em_evtgen_theta, (float)NAN );
0483   tc->set_property( PidCandidate::em_evtgen_phi, (float)NAN );
0484   tc->set_property( PidCandidate::em_evtgen_eta, (float)NAN );
0485   tc->set_property( PidCandidate::em_evtgen_charge, (int)NAN );
0486
0487   /* If matching truth primary particle found: update truth information */
0488   CaloRawClusterEval* clustereval = caloevalstack->get_rawcluster_eval();
0489   PHG4Particle* primary = clustereval->max_truth_primary_particle_by_energy(cluster);
0490
0491   if ( primary )
0492     {
0493       /* get particle momenta and theta, phi angles */
0494       float gpx = primary->get_px();
0495       float gpy = primary->get_py();
0496       float gpz = primary->get_pz();
0497       float gpt = sqrt(gpx * gpx + gpy * gpy);
0498       float gptotal = sqrt(gpx * gpx + gpy * gpy + gpz * gpz);
0499       //float ge = (float)primary->get_e();
0500
0501       float gphi = NAN;
0502       float gtheta = NAN;
0503       float geta = NAN;
0504
0505       if (gpt != 0.0)
0506     {
0507       gtheta = atan2( gpt, gpz );
0508       geta = -1 * log( tan( gtheta / 2.0 ) );
0509     }
0510
0511       gphi = atan2(gpy, gpx);
0512
0513       /* get charge based on PDG code of particle */
0514       TParticlePDG * pdg_p = TDatabasePDG::Instance()->GetParticle( primary->get_pid() );
0515       int gcharge = -999;
0516       if ( pdg_p )
0517     {
0518       /* NOTE: TParticlePDG::Charge() returns charge in units of |e|/3 (see ROOT documentation) */
0519       gcharge = pdg_p->Charge() / 3;
0520     }
0521
0522       tc->set_property( PidCandidate::em_evtgen_pid, primary->get_pid() );
0523       tc->set_property( PidCandidate::em_evtgen_ptotal, gptotal );
0524       tc->set_property( PidCandidate::em_evtgen_theta, gtheta );
0525       tc->set_property( PidCandidate::em_evtgen_phi, gphi );
0526       tc->set_property( PidCandidate::em_evtgen_eta, geta );
0527       tc->set_property( PidCandidate::em_evtgen_charge, gcharge );
0528
0529
0530       /* Track pointing to cluster */
0531
0532       bool fill_in_truth = false;
0533
0534       /* Try to find a track which best points to the current cluster */
0535       /* TODO: Depending on cluster theta, deltaR changes (currently -1) */
0536
0537       SvtxTrack_FastSim* the_track = _trackproj->get_best_track(trackmap, cluster);
0538
0539       /* Test if a best track was found */
0540       if(the_track!=NULL) //Fill in reconstructed info
0541     {
0542       /* Get track position and momentum */
0543
0544       SvtxTrackState *the_state = _trackproj->get_best_state(the_track, cluster);
0545       double posv[3] = {0.,0.,0.};
0546       posv[0] = the_state->get_x();
0547       posv[1] = the_state->get_y();
0548       posv[2] = the_state->get_z();
0549       /* Extrapolate track2cluster values */
0550       float track2cluster_theta = atan(sqrt(posv[0]*posv[0]+
0551                         posv[1]*posv[1]) / posv[2]);
0552       float track2cluster_eta = -log(abs(tan(track2cluster_theta/2)));
0553       if(tan(track2cluster_theta/2)<0)
0554         track2cluster_eta*=-1;
0555       float track2cluster_phi = atan(posv[1]/posv[0]);
0556       float track2cluster_x = posv[0];
0557       float track2cluster_y = posv[1];
0558       float track2cluster_z = posv[2];
0559
0560       /* Fills in weird values right now */
0561       float track2cluster_p = 0;
0562
0563
0564       tc->set_property( PidCandidate::em_track_theta2cluster, track2cluster_theta );
0565       tc->set_property( PidCandidate::em_track_eta2cluster, track2cluster_eta );
0566       tc->set_property( PidCandidate::em_track_phi2cluster, track2cluster_phi );
0567       tc->set_property( PidCandidate::em_track_p2cluster, track2cluster_p );
0568       tc->set_property( PidCandidate::em_track_x2cluster, track2cluster_x );
0569       tc->set_property( PidCandidate::em_track_y2cluster, track2cluster_y);
0570       tc->set_property( PidCandidate::em_track_z2cluster, track2cluster_z);
0571
0572
0573       tc->set_property( PidCandidate::em_track_id, (uint)primary->get_track_id() );
0574       tc->set_property( PidCandidate::em_track_quality, (float)100.0 );
0575       tc->set_property( PidCandidate::em_track_theta,2*atan( exp(-the_track->get_eta()) ));
0576       tc->set_property( PidCandidate::em_track_phi, the_track->get_phi() );
0577       tc->set_property( PidCandidate::em_track_eta, the_track->get_eta() );
0578       tc->set_property( PidCandidate::em_track_ptotal, the_track->get_p() );
0579       tc->set_property( PidCandidate::em_track_ptrans, the_track->get_pt() );
0580       tc->set_property( PidCandidate::em_track_charge, the_track->get_charge() );
0581       tc->set_property( PidCandidate::em_track_dca, NAN );
0582       tc->set_property( PidCandidate::em_track_section, (uint)0 );
0583       tc->set_property( PidCandidate::em_track_e3x3_cemc, NAN );
0584       tc->set_property( PidCandidate::em_track_e3x3_ihcal, NAN );
0585       tc->set_property( PidCandidate::em_track_e3x3_ohcal, NAN );
0586       tc->set_property( PidCandidate::em_track_cluster_dr, (float)0.0 );
0587
0588       tc->set_property( PidCandidate::em_pid_prob_electron, (float)0.0 );
0589       tc->set_property( PidCandidate::em_pid_prob_pion, (float)0.0 );
0590       tc->set_property( PidCandidate::em_pid_prob_kaon, (float)0.0 );
0591       tc->set_property( PidCandidate::em_pid_prob_proton, (float)0.0 );
0592
0593
0594       /*
0595          cout << "------------------------------------------" << endl;
0596          cout << "Track Projection Successful!" << endl;
0597          cout << "Cluster Detector = " << _trackproj->get_detector() << "EMC" << endl;
0598          cout << "Cluster X: " << cluster->get_x() << " | Track X: " << posv[0] << endl;
0599          cout << "Cluster Y: " << cluster->get_y() << " | Track Y: " << posv[1] << endl;
0600          cout << "Cluster Z: " << cluster->get_z() << " | Track Z: " << posv[2] << endl;
0601       */
0602
0603
0604     }
0605       else if(fill_in_truth) //Fill in truth
0606     {
0607       tc->set_property( PidCandidate::em_track_theta2cluster, NAN);
0608       tc->set_property( PidCandidate::em_track_eta2cluster, NAN);
0609       tc->set_property( PidCandidate::em_track_phi2cluster, NAN );
0610       tc->set_property( PidCandidate::em_track_p2cluster, NAN);
0611       tc->set_property( PidCandidate::em_track_x2cluster, NAN );
0612       tc->set_property( PidCandidate::em_track_y2cluster, NAN);
0613       tc->set_property( PidCandidate::em_track_z2cluster, NAN);
0614       tc->set_property( PidCandidate::em_track_id, (uint)primary->get_track_id() );
0615       tc->set_property( PidCandidate::em_track_quality, (float)100.0 );
0616       tc->set_property( PidCandidate::em_track_theta, gtheta );
0617       tc->set_property( PidCandidate::em_track_phi, gphi );
0618       tc->set_property( PidCandidate::em_track_eta, geta );
0619       tc->set_property( PidCandidate::em_track_ptotal, gptotal );
0620       tc->set_property( PidCandidate::em_track_ptrans, gpt );
0621       tc->set_property( PidCandidate::em_track_charge, gcharge );
0622       tc->set_property( PidCandidate::em_track_dca, NAN );
0623       tc->set_property( PidCandidate::em_track_section, (uint)0 );
0624       tc->set_property( PidCandidate::em_track_e3x3_cemc, NAN );
0625       tc->set_property( PidCandidate::em_track_e3x3_ihcal, NAN );
0626       tc->set_property( PidCandidate::em_track_e3x3_ohcal, NAN );
0627       tc->set_property( PidCandidate::em_track_cluster_dr, (float)0.0 );
0628
0629       /* Get information on track from PID detectors */
0630       tc->set_property( PidCandidate::em_pid_prob_electron, (float)0.0 );
0631       tc->set_property( PidCandidate::em_pid_prob_pion, (float)0.0 );
0632       tc->set_property( PidCandidate::em_pid_prob_kaon, (float)0.0 );
0633       tc->set_property( PidCandidate::em_pid_prob_proton, (float)0.0 );
0634
0635       if ( abs( primary->get_pid() ) == 11 )
0636         tc->set_property( PidCandidate::em_pid_prob_electron, (float)1.0 );
0637
0638       if ( abs( primary->get_pid() ) == 211 )
0639         tc->set_property( PidCandidate::em_pid_prob_pion, (float)1.0 );
0640
0641       if ( abs( primary->get_pid() ) == 321 )
0642         tc->set_property( PidCandidate::em_pid_prob_kaon, (float)1.0 );
0643
0644       if ( abs( primary->get_pid() ) == 2212 )
0645         tc->set_property( PidCandidate::em_pid_prob_proton, (float)1.0 );
0646     }
0647       else
0648     {
0649       tc->set_property( PidCandidate::em_track_theta2cluster, NAN);
0650       tc->set_property( PidCandidate::em_track_eta2cluster, NAN);
0651       tc->set_property( PidCandidate::em_track_phi2cluster, NAN );
0652       tc->set_property( PidCandidate::em_track_p2cluster, NAN);
0653       tc->set_property( PidCandidate::em_track_x2cluster, NAN );
0654       tc->set_property( PidCandidate::em_track_y2cluster, NAN);
0655       tc->set_property( PidCandidate::em_track_z2cluster, NAN);
0656       tc->set_property( PidCandidate::em_track_id,(uint)0);
0657       tc->set_property( PidCandidate::em_track_quality, (float)0.0 );
0658       tc->set_property( PidCandidate::em_track_theta, NAN );
0659       tc->set_property( PidCandidate::em_track_phi, NAN );
0660       tc->set_property( PidCandidate::em_track_eta, NAN );
0661       tc->set_property( PidCandidate::em_track_ptotal, NAN );
0662       tc->set_property( PidCandidate::em_track_ptrans, NAN );
0663       tc->set_property( PidCandidate::em_track_charge, 0 );
0664       tc->set_property( PidCandidate::em_track_dca, NAN );
0665       tc->set_property( PidCandidate::em_track_section, (uint)0 );
0666       tc->set_property( PidCandidate::em_track_e3x3_cemc, NAN );
0667       tc->set_property( PidCandidate::em_track_e3x3_ihcal, NAN );
0668       tc->set_property( PidCandidate::em_track_e3x3_ohcal, NAN );
0669       tc->set_property( PidCandidate::em_track_cluster_dr, (float)0.0 );
0670
0671       /* Get information on track from PID detectors */
0672       tc->set_property( PidCandidate::em_pid_prob_electron, (float)0.0 );
0673       tc->set_property( PidCandidate::em_pid_prob_pion, (float)0.0 );
0674       tc->set_property( PidCandidate::em_pid_prob_kaon, (float)0.0 );
0675       tc->set_property( PidCandidate::em_pid_prob_proton, (float)0.0 );
0676
0677     }
0678     }
0679
0680   /* add tau candidate to collection */
0681   candidateMap.insert( make_pair( tc->get_candidate_id(), tc ) );
0682   return 0;
0683 }
0684
0685
0686 int
0687 DISKinematicsReco::AddGlobalCalorimeterInformation()
0688 {
0689   /* Missing transverse energy, transverse energy direction, and Energy sums, x and y components separately */
0690   float Et_miss = 0;
0691   float Et_miss_phi = 0;
0692   float Ex_sum = 0;
0693   float Ey_sum = 0;
0694
0695   /* energy threshold for considering tower */
0696   float tower_emin = 0.0;
0697
0698   /* Loop over all tower (and geometry) collections */
0699   vector < string > calo_names;
0700   calo_names.push_back( "CEMC" );
0701   calo_names.push_back( "HCALIN" );
0702   calo_names.push_back( "HCALOUT" );
0703   calo_names.push_back( "FEMC" );
0704   calo_names.push_back( "FHCAL" );
0705   calo_names.push_back( "EEMC" );
0706
0707   for ( unsigned i = 0; i < calo_names.size(); i++ )
0708     {
0709
0710       string towernodename = "TOWER_CALIB_" + calo_names.at( i );
0711       RawTowerContainer* towers = findNode::getClass<RawTowerContainer>(_topNode, towernodename.c_str());
0712       if (!towers)
0713         {
0714           std::cout << PHWHERE << ": Could not find node " << towernodename.c_str() << std::endl;
0715           return -1;
0716         }
0717
0718       string towergeomnodename = "TOWERGEOM_" + calo_names.at( i );
0719       RawTowerGeomContainer *towergeoms = findNode::getClass<RawTowerGeomContainer>(_topNode, towergeomnodename.c_str());
0720       if (! towergeoms)
0721         {
0722           cout << PHWHERE << ": Could not find node " << towergeomnodename.c_str() << endl;
0723           return -1;
0724         }
0725
0726       /* define tower iterator */
0727       RawTowerContainer::ConstRange begin_end = towers->getTowers();
0728       RawTowerContainer::ConstIterator rtiter;
0729
0730       /* loop over all tower in CEMC calorimeter */
0731       for (rtiter = begin_end.first; rtiter !=  begin_end.second; ++rtiter)
0732         {
0733           /* get tower energy */
0734           RawTower *tower = rtiter->second;
0735           float tower_energy = tower->get_energy();
0736
0737           /* check if tower above energy treshold */
0738           if ( tower_energy < tower_emin )
0739             continue;
0740
0741           /* get eta and phi of tower and check angle delta_R w.r.t. jet axis */
0742           RawTowerGeom * tower_geom = towergeoms->get_tower_geometry(tower -> get_key());
0743           float tower_eta = tower_geom->get_eta();
0744           float tower_phi = tower_geom->get_phi();
0745
0746           /* from https://en.wikipedia.org/wiki/Pseudorapidity:
0747              p_x = p_T * cos( phi )
0748              p_y = p_T * sin( phi )
0749              p_z = p_T * sinh( eta )
0750              |p| = p_T * cosh( eta )
0751           */
0752
0753           /* calculate 'transverse' tower energy */
0754           float tower_energy_t = tower_energy / cosh( tower_eta );
0755
0756           /* add energy components of this tower to total energy components */
0757           Ex_sum += tower_energy_t * cos( tower_phi );
0758           Ey_sum += tower_energy_t * sin( tower_phi );
0759         }
0760     }
0761
0762   /* calculate Et_miss */
0763   Et_miss = sqrt( Ex_sum * Ex_sum + Ey_sum * Ey_sum );
0764   Et_miss_phi = atan( Ey_sum / Ex_sum );
0765
0766   /* set event branch paraemters */
0767   ( _map_event_branches.find( "Et_miss" ) )->second = Et_miss;
0768   ( _map_event_branches.find( "Et_miss_phi" ) )->second = Et_miss_phi;
0769
0770   return 0;
0771 }
0772
0773 int
0774 DISKinematicsReco::WriteCandidatesToTree( type_map_tcan& electronCandidateMap )
0775 {
0776   /* Get number of electron candidates */
0777   ( _map_event_branches.find( "em_ncandidates" ) )->second = electronCandidateMap.size();
0778
0779   /* Loop over all EM candidates and add them to tree */
0780   for (type_map_tcan::iterator iter = electronCandidateMap.begin();
0781        iter != electronCandidateMap.end();
0782        ++iter)
0783     {
0784       /* update information in map and fill tree */
0785       for ( map< PidCandidate::PROPERTY , vector<float> >::iterator iter_prop = _map_em_candidate_branches.begin();
0786             iter_prop != _map_em_candidate_branches.end();
0787             ++iter_prop)
0788         {
0789           switch ( PidCandidate::get_property_info( (iter_prop->first) ).second ) {
0790
0791           case PidCandidate::type_float :
0792             (iter_prop->second).push_back( (iter->second)->get_property_float( (iter_prop->first) ) );
0793             break;
0794
0795           case PidCandidate::type_int :
0796             (iter_prop->second).push_back( (iter->second)->get_property_int( (iter_prop->first) ) );
0797             break;
0798
0799           case PidCandidate::type_uint :
0800             (iter_prop->second).push_back( (iter->second)->get_property_uint( (iter_prop->first) ) );
0801             break;
0802
0803           case PidCandidate::type_unknown :
0804             break;
0805           }
0806         }
0807     }
0808
0809   return 0;
0810 }
0811
0812
0813 PidCandidate*
0814 DISKinematicsReco::FindMinDeltaRCandidate( type_map_tcan *candidates, const float eta_ref, const float phi_ref )
0815 {
0816   /* PidCandidate with eta, phi closest to reference */
0817   PidCandidate* best_candidate = NULL;
0818
0819   return best_candidate;
0820 }
0821
0822
0823 float
0824 DISKinematicsReco::CalculateDeltaR( float eta1, float phi1, float eta2, float phi2 )
0825 {
0826   /* Particles at phi = PI+x and phi = PI-x are actually close to each other in phi, but simply calculating
0827    * the difference in phi would give a large distance (because phi ranges from -PI to +PI in the convention
0828    * used. Account for this by subtracting 2PI is particles fall within this border area. */
0829   if( ( phi1 < -0.9*TMath::Pi()) && ( phi2 >  0.9*TMath::Pi() ) ) phi2 = phi2 - 2*TMath::Pi();
0830   if( ( phi1 >  0.9*TMath::Pi()) && ( phi2 < -0.9*TMath::Pi() ) ) phi1 = phi1 - 2*TMath::Pi();
0831
0832   float delta_R = sqrt( pow( eta2 - eta1, 2 ) + pow( phi2 - phi1, 2 ) );
0833
0834   return delta_R;
0835 }
0836
0837 int
0838 DISKinematicsReco::AddReconstructedKinematics( type_map_tcan& em_candidates , string mode )
0839 {
0840   /* Loop over all EM candidates */
0841   for (type_map_tcan::iterator iter = em_candidates.begin();
0842        iter != em_candidates.end();
0843        ++iter)
0844     {
0845       /* Add reco kinematics based on scattered electron data */
0846       PidCandidate* the_electron = iter->second;
0847
0848       float e0_E = _beam_electron_ptotal;
0849       float p0_E = _beam_hadron_ptotal;
0850
0851       /* get scattered particle kinematicsbased on chosen mode */
0852       float e1_E = NAN;
0853       float e1_theta = NAN;
0854
0855       if ( mode == "cluster" )
0856         {
0857           e1_E = the_electron->get_property_float( PidCandidate::em_cluster_e );
0858           e1_theta = the_electron->get_property_float( PidCandidate::em_cluster_theta );
0859         }
0860       else if ( mode == "truth" )
0861         {
0862           e1_E = the_electron->get_property_float( PidCandidate::em_evtgen_ptotal );
0863           e1_theta = the_electron->get_property_float( PidCandidate::em_evtgen_theta );
0864         }
0865       else
0866         {
0867           cout << "WARNING: Unknown mode " << mode << " selected." << endl;
0868           return -1;
0869         }
0870
0871       /* for purpose of calculations, 'theta' angle of the scattered electron is defined as angle
0872          between 'scattered electron' and 'direction of incoming electron'. Since initial electron
0873          has 'theta = Pi' in coordinate system, we need to use 'theta_rel = Pi - theta' for calculations
0874          of kinematics. */
0875       float e1_theta_rel = M_PI - e1_theta;
0876
0877       /* event kinematics */
0878       float dis_s = 4.0 * e0_E * p0_E;
0879
0880       float dis_Q2 = 2.0 * e0_E * e1_E * ( 1 - cos( e1_theta_rel ) );
0881
0882       /* ePHENIX LOI definition of y: */
0883       float dis_y = 1.0 - ( e1_E / e0_E ) + ( dis_Q2 / ( 4.0 * e0_E * e0_E ) );
0884
0885       /* G. Wolf Hera Physics definitions of y: */
0886       //float dis_y = 1 - ( e1_E / (2*e0_E) ) * ( 1 - cos( e1_theta_rel ) );
0887
0888       float dis_x = dis_Q2 / ( dis_s * dis_y );
0889
0890       float dis_W2 = _mproton*_mproton + dis_Q2 * ( ( 1 / dis_x ) - 1 );
0891       float dis_W = sqrt( dis_W2 );
0892
0893       the_electron->set_property( PidCandidate::em_reco_x_e, dis_x );
0894       the_electron->set_property( PidCandidate::em_reco_y_e, dis_y );
0895       the_electron->set_property( PidCandidate::em_reco_q2_e, dis_Q2 );
0896       the_electron->set_property( PidCandidate::em_reco_w_e, dis_W );
0897     }
0898
0899   return 0;
0900 }
0901
0902 int
0903 DISKinematicsReco::AddTruthEventInformation()
0904 {
0905   /* Get collection of truth particles from event generator */
0906   PHHepMCGenEventMap *geneventmap = findNode::getClass<PHHepMCGenEventMap>(_topNode,"PHHepMCGenEventMap");
0907   if (!geneventmap) {
0908     std::cout << PHWHERE << " WARNING: Can't find requested PHHepMCGenEventMap" << endl;
0909     return -1;
0910   }
0911
0912   /* Add truth kinematics */
0913   int embedding_id = 1;
0914   PHHepMCGenEvent *genevt = geneventmap->get(embedding_id);
0915   if (!genevt)
0916     {
0917       std::cout << PHWHERE << "WARNING: Node PHHepMCGenEventMap missing subevent with embedding ID "<< embedding_id;
0918       std::cout <<". Print PHHepMCGenEventMap:";
0919       geneventmap->identify();
0920       return -1;
0921     }
0922
0923   HepMC::GenEvent* theEvent = genevt->getEvent();
0924
0925   if ( !theEvent )
0926     {
0927       std::cout << PHWHERE << "WARNING: Missing requested GenEvent!" << endl;
0928       return -1;
0929     }
0930
0931   int true_process_id = theEvent->signal_process_id();
0932   float ev_x = -NAN;
0933   float ev_Q2 = -NAN;
0934   float ev_W = -NAN;
0935   float ev_y = -NAN;
0936
0937   float ev_s = 4 * _beam_electron_ptotal * _beam_hadron_ptotal;
0938
0939
0940   if ( theEvent->pdf_info() )
0941     {
0942       //ev_x = theEvent->pdf_info()->x1();
0943       ev_x = theEvent->pdf_info()->x2();
0944       ev_Q2 = theEvent->pdf_info()->scalePDF();
0945       ev_y = ev_Q2 / ( ev_x * ev_s );
0946       ev_W = sqrt( _mproton*_mproton + ev_Q2 * ( 1 / ev_x - 1 ) );
0947     }
0948
0949   ( _map_event_branches.find( "evtgen_process_id" ) )->second = true_process_id;
0950   ( _map_event_branches.find( "evtgen_s" ) )->second = ev_s;
0951   ( _map_event_branches.find( "evtgen_x" ) )->second = ev_x;
0952   ( _map_event_branches.find( "evtgen_Q2" ) )->second = ev_Q2;
0953   ( _map_event_branches.find( "evtgen_y" ) )->second = ev_y;
0954   ( _map_event_branches.find( "evtgen_W" ) )->second = ev_W;
0955
0956   return 0;
0957 }
0958
0959 void
0960 DISKinematicsReco::ResetBranchMap()
0961 {
0962   /* Event branches */
0963   for ( map< string , float >::iterator iter = _map_event_branches.begin();
0964         iter != _map_event_branches.end();
0965         ++iter)
0966     {
0967       (iter->second) = NAN;
0968     }
0969
0970   /* EM candidate branches */
0971   for ( map< PidCandidate::PROPERTY , vector<float> >::iterator iter = _map_em_candidate_branches.begin();
0972         iter != _map_em_candidate_branches.end();
0973         ++iter)
0974     {
0975       (iter->second).clear();
0976     }
0977
0978   return;
0979 }
0980
0981 int
0982 DISKinematicsReco::End(PHCompositeNode *topNode)
0983 {
0984   _tfile->cd();
0985
0986   if ( _tree_event_cluster )
0987     _tree_event_cluster->Write();
0988
0989   if ( _tree_event_truth )
0990     _tree_event_truth->Write();
0991
0992   _tfile->Close();
0993
0994   return 0;
0995 }