analysis/EICAnalysis/ExclusiveReco.C

0001 #include "ExclusiveReco.h"
0002 #include "TruthTrackerHepMC.h"
0003 #include "TrackProjectorPlaneECAL.h"
0004 #include "DVMPHelper.h"
0005 /* STL includes */
0006 #include <cassert>
0007
0008 /* Fun4All includes */
0009 #include <trackbase_historic/SvtxTrackMap_v1.h>
0010
0011 #include <phool/getClass.h>
0012
0013 #include <fun4all/Fun4AllReturnCodes.h>
0014
0015 #include <calobase/RawTowerGeomContainer.h>
0016 #include <calobase/RawTowerContainer.h>
0017 #include <calobase/RawTowerGeom.h>
0018 #include <calobase/RawTowerv1.h>
0019 #include <calobase/RawTowerDefs.h>
0020
0021 #include <calobase/RawClusterContainer.h>
0022 #include <calobase/RawCluster.h>
0023
0024 #include <g4eval/CaloEvalStack.h>
0025
0026 #include <g4main/PHG4Particle.h>
0027
0028
0029 /* ROOT includes */
0030 #include <TString.h>
0031 #include <TTree.h>
0032 #include <TFile.h>
0033
0034 using namespace std;
0035
0036 ExclusiveReco::ExclusiveReco(std::string filename) :
0037   SubsysReco("ExclusiveReco" ),
0038   _mproton( 0.938272 ),
0039   _save_towers(false),
0040   _save_tracks(false),
0041   _ievent(0),
0042   _filename(filename),
0043   _tfile(nullptr),
0044   _tree_invariant_mass(nullptr),
0045   _tree_event_reco(nullptr),
0046   _tree_event_truth(nullptr),
0047   _beam_electron_ptotal(10),
0048   _beam_hadron_ptotal(250),
0049   _trackproj(nullptr)
0050 {
0051
0052 }
0053
0054 int
0055 ExclusiveReco::Init(PHCompositeNode *topNode)
0056 {
0057   _topNode = topNode;
0058
0059   _ievent = 0;
0060
0061   _tfile = new TFile(_filename.c_str(), "RECREATE");
0062
0063   /* Create tree for information about full event */
0064   _tree_invariant_mass = new TTree("event_exclusive", "A tree with exclusive event information");
0065   _tree_invariant_mass->Branch("event", &_ievent, "event/I");
0066
0067   _tree_invariant_mass->Branch("reco_inv", &_vect1);
0068   _tree_invariant_mass->Branch("true_inv", &_vect2);
0069   _tree_invariant_mass->Branch("reco_inv_decay", &_vect3);
0070   _tree_invariant_mass->Branch("reco_inv_scatter", &_vect4);
0071   _tree_invariant_mass->Branch("true_inv_decay", &_vect5);
0072   _tree_invariant_mass->Branch("true_inv_scatter", &_vect6);
0073
0074   /* Create tree for information about reco event */
0075   _tree_event_reco = new TTree("event_reco", "A tree with reco particle information");
0076   _tree_event_reco->Branch("event", &_ievent, "event/I");
0077   _tree_event_reco->Branch("eta", &reco_eta);
0078   _tree_event_reco->Branch("phi", &reco_phi);
0079   _tree_event_reco->Branch("ptotal", &reco_ptotal);
0080   _tree_event_reco->Branch("charge", &reco_charge);
0081   _tree_event_reco->Branch("cluster_e", &reco_cluster_e);
0082   _tree_event_reco->Branch("is_scattered_lepton", &reco_is_scattered_lepton);
0083
0084
0085   /* Create tree for information about reco event */
0086   _tree_event_truth = new TTree("event_truth", "A tree with truth particle information");
0087   _tree_event_truth->Branch("event", &_ievent, "event/I");
0088   _tree_event_truth->Branch("geta", &true_eta);
0089   _tree_event_truth->Branch("gphi", &true_phi);
0090   _tree_event_truth->Branch("gptotal", &true_ptotal);
0091   _tree_event_truth->Branch("pid",&true_pid);
0092   _tree_event_truth->Branch("is_scattered_lepton",&is_scattered_lepton);
0093
0094   return 0;
0095 }
0096
0097 int
0098 ExclusiveReco::InitRun(PHCompositeNode *topNode)
0099 {
0100    _trackproj=new TrackProjectorPlaneECAL( topNode );
0101    return 0;
0102 }
0103
0104 int
0105 ExclusiveReco::process_event(PHCompositeNode *topNode)
0106 {
0107   /* Calculate invariant mass of a DVMP event */
0108   if ( _do_process_dvmp )
0109     {
0110       AddInvariantMassInformation();
0111     }
0112   /* count up event number */
0113   _ievent ++;
0114
0115   return 0;
0116 }
0117
0118 int
0119 ExclusiveReco::AddInvariantMassInformation()
0120 {
0121   /* First, add truth particle information */
0122   // ------------------------------------------------------------------------//
0123   //true vectors are defined in header file
0124   true_eta.clear(); true_phi.clear(); true_ptotal.clear(); true_pid.clear(); is_scattered_lepton.clear();
0125   // -----------------------------------------------------------------------//
0126   /* Get collection of truth particles from event generator */
0127   PHHepMCGenEventMap *geneventmap = findNode::getClass<PHHepMCGenEventMap>(_topNode,"PHHepMCGenEventMap");
0128   if (!geneventmap) {
0129     std::cout << PHWHERE << " WARNING: Can't find requested PHHepMCGenEventMap" << endl;
0130     return -1;
0131   }
0132
0133   /* Add truth kinematics */
0134   int embedding_id = 1;
0135   PHHepMCGenEvent *genevt = geneventmap->get(embedding_id);
0136   if (!genevt)
0137     {
0138       std::cout << PHWHERE << "WARNING: Node PHHepMCGenEventMap missing subevent with embedding ID "<< embedding_id;
0139       std::cout <<". Print PHHepMCGenEventMap:";
0140       geneventmap->identify();
0141       return -1;
0142     }
0143
0144   HepMC::GenEvent* theEvent = genevt->getEvent();
0145
0146   if ( !theEvent )
0147     {
0148       std::cout << PHWHERE << "WARNING: Missing requested GenEvent!" << endl;
0149       return -1;
0150     }
0151
0152   /* Look for scattered lepton */
0153   TruthTrackerHepMC truth;
0154   truth.set_hepmc_geneventmap( geneventmap );
0155   HepMC::GenParticle* particle_scattered_l = truth.FindScatteredLepton();
0156   /* loop over all particles */
0157   for (HepMC::GenEvent::particle_const_iterator p = theEvent->particles_begin();
0158        p != theEvent->particles_end(); ++p)
0159     {
0160       /* skip particles that are not stable final state particles (status 1) */
0161       if ( (*p)->status() != 1 )
0162         continue;
0163
0164       float mom_eta = -1 * log ( tan( (*p)->momentum().theta() / 2.0 ) );
0165
0166       true_eta.push_back(mom_eta);
0167       true_phi.push_back( (*p)->momentum().phi() );
0168       true_ptotal.push_back( (*p)->momentum().e() );
0169       true_pid.push_back( (*p)->pdg_id() );
0170
0171       if ( particle_scattered_l &&
0172        (*p) == particle_scattered_l )
0173     is_scattered_lepton.push_back(true);
0174       else
0175     is_scattered_lepton.push_back(false);
0176     }
0177   /* Second, add reconstructed particle information */
0178   // --------------------------------------------------------------------------
0179   //reco vectors are defined in header file
0180   reco_eta.clear(); reco_phi.clear(); reco_ptotal.clear(); reco_cluster_e.clear(); reco_charge.clear(); reco_is_scattered_lepton.clear();
0181   // --------------------------------------------------------------------------
0182   vector< string > v_ecals;
0183   v_ecals.push_back("EEMC");
0184   v_ecals.push_back("CEMC");
0185   v_ecals.push_back("FEMC");
0186   for ( unsigned idx = 0; idx < v_ecals.size(); idx++ )
0187     {
0188       string clusternodename = "CLUSTER_" + v_ecals.at( idx );
0189       RawClusterContainer *clusterList = findNode::getClass<RawClusterContainer>(_topNode,clusternodename.c_str());
0190       SvtxTrackMap *trackmap = findNode::getClass<SvtxTrackMap>(_topNode,"SvtxTrackMap");
0191       if (!clusterList) {
0192         cerr << PHWHERE << " ERROR: Can't find node " << clusternodename << endl;
0193         return false;
0194       }
0195       if(!trackmap) {
0196     cerr << PHWHERE << " ERROR: Can't find node SvtxTrackMap" << endl;
0197     return false;
0198       }
0199
0200       for (unsigned int k = 0; k < clusterList->size(); ++k)
0201         {
0202           RawCluster *cluster = clusterList->getCluster(k);
0203           /* Check if cluster energy is below threshold */
0204           float e_cluster_threshold = 0.3;
0205           if ( cluster->get_energy() < e_cluster_threshold )
0206             continue;
0207
0208       _trackproj->set_detector(_trackproj->get_detector_from_cluster(cluster));
0209
0210       SvtxTrack *best_track = _trackproj->get_best_track(trackmap,cluster);
0211       if(best_track!=NULL)
0212         {
0213           reco_eta.push_back(best_track->get_eta());
0214           reco_phi.push_back(best_track->get_phi());
0215           // See if particle is in EEMC, if so, swap ptotal with cluster e
0216           RawCluster::TowerConstIterator rtiter = cluster->get_towers().first;
0217           const char * cc = RawTowerDefs::convert_caloid_to_name( RawTowerDefs::decode_caloid( rtiter->first )).c_str();
0218           if(strcmp(cc,"EEMC")==0)
0219         {
0220           reco_ptotal.push_back(cluster->get_energy());
0221         }
0222           else
0223         {
0224           reco_ptotal.push_back(best_track->get_p());
0225         }
0226           reco_charge.push_back(best_track->get_charge());
0227           reco_cluster_e.push_back(cluster->get_energy());
0228
0229           // Look at every truth particle
0230           for(unsigned i = 0 ; i < true_eta.size() ; i++)
0231         {
0232           // Compute difference in reco particle eta and phi w/ true
0233           float eta_diff = abs(best_track->get_eta()-true_eta.at(i));
0234           float phi_diff = abs(best_track->get_phi()-true_phi.at(i));
0235
0236           // If the differences in these values is very small
0237           if(eta_diff<0.1&&phi_diff<0.1)
0238             {
0239               bool b = is_scattered_lepton.at(i)&&best_track->get_charge()==-1;
0240               if(b)
0241             reco_is_scattered_lepton.push_back(true);
0242               else
0243             reco_is_scattered_lepton.push_back(false);
0244             }
0245           else
0246             {
0247               reco_is_scattered_lepton.push_back(false);
0248             }
0249         }
0250         }
0251       else
0252         {
0253           reco_eta.push_back(NAN);
0254           reco_phi.push_back(NAN);
0255           reco_ptotal.push_back(NAN);
0256           reco_charge.push_back(0);
0257           reco_cluster_e.push_back(NAN);
0258           reco_is_scattered_lepton.push_back(false);
0259         }
0260     }
0261     }
0262
0263   // At this point, we have all the truth and reco event information we need to fiddle around with measuring the invariant mass //
0264   DVMPHelper * dvmp = new DVMPHelper(reco_eta,reco_phi,reco_ptotal,reco_charge,reco_cluster_e,reco_is_scattered_lepton,true_eta,true_phi,true_ptotal,true_pid,is_scattered_lepton);
0265
0266   // 1) Invariant Mass of all reconstructed e- e+ pairs
0267   // 2) Invariant Mass of all truth e- e+ pairs
0268   // 3) Invariant Mass of reco decay e- e+ pairs (**USES TRUTH INFO**)
0269   // 4) Invariant Mass of reco scattered e- e+ pairs (**USES TRUTH INFO**)
0270   // 5) Invariant Mass of truth decay e- e+ pairs
0271   // 6) Invariant Mass of truth scattered e- e+ pairs
0272
0273   _vect1 = dvmp->calculateInvariantMass_1();
0274   _vect2 = dvmp->calculateInvariantMass_2();
0275   _vect3 = dvmp->calculateInvariantMass_3();
0276   _vect4 = dvmp->calculateInvariantMass_4();
0277   _vect5 = dvmp->calculateInvariantMass_5();
0278   _vect6 = dvmp->calculateInvariantMass_6();
0279
0280   _tree_event_reco->Fill();
0281   _tree_event_truth->Fill();
0282   _tree_invariant_mass->Fill();
0283   return 0;
0284 }
0285
0286 int
0287 ExclusiveReco::End(PHCompositeNode *topNode)
0288 {
0289   _tfile->cd();
0290
0291   if ( _tree_invariant_mass )
0292     {
0293       _tree_event_reco->Write();
0294       _tree_event_truth->Write();
0295       _tree_invariant_mass->Write();
0296     }
0297
0298   _tfile->Close();
0299
0300   return 0;
0301 }