singlePhotonClusterAna/src/singlePhotonClusterAna.cc

0001 #include "singlePhotonClusterAna.h"
0002
0003 //Fun4all stuff
0004 #include <fun4all/Fun4AllReturnCodes.h>
0005 #include <fun4all/Fun4AllServer.h>
0006 #include <fun4all/Fun4AllHistoManager.h>
0007 #include <phool/PHCompositeNode.h>
0008 #include <phool/getClass.h>
0009 #include <phool/phool.h>
0010 #include <ffaobjects/EventHeader.h>
0011
0012 //ROOT stuff
0013 #include <TH1F.h>
0014 #include <TH2F.h>
0015 #include <TH3F.h>
0016 #include <TFile.h>
0017 #include <TLorentzVector.h>
0018 #include <TTree.h>
0019
0020 //for emc clusters
0021 #include <calobase/RawCluster.h>
0022 #include <calobase/RawClusterContainer.h>
0023 #include <calobase/RawClusterUtility.h>
0024 #include <calobase/RawTowerGeomContainer.h>
0025 #include <calobase/RawTower.h>
0026 #include <calobase/RawTowerContainer.h>
0027
0028 //for vetex information
0029 #include <globalvertex/GlobalVertex.h>
0030 #include <globalvertex/GlobalVertexMap.h>
0031
0032 //tracking
0033 #include <trackbase_historic/SvtxTrack.h>
0034 #include <trackbase_historic/SvtxTrackMap.h>
0035 #include <trackbase_historic/SvtxVertex.h>
0036 #include <trackbase_historic/SvtxVertexMap.h>
0037
0038 //truth information
0039 #include <g4main/PHG4TruthInfoContainer.h>
0040 #include <g4main/PHG4Particle.h>
0041 #include <g4main/PHG4VtxPoint.h>
0042 #pragma GCC diagnostic push
0043 #pragma GCC diagnostic ignored "-Wdeprecated-declarations"
0044 #include <HepMC/GenEvent.h>
0045 #include <HepMC/GenParticle.h>
0046 #include <HepMC/GenVertex.h>
0047 #include <HepMC/IteratorRange.h>
0048 #include <HepMC/SimpleVector.h>
0049 #include <HepMC/GenParticle.h>
0050 #pragma GCC diagnostic pop
0051 #include <CLHEP/Geometry/Point3D.h>
0052
0053
0054 //____________________________________________________________________________..
0055 singlePhotonClusterAna::singlePhotonClusterAna(const std::string &name, const std::string &outName = "singlePhotonClusterAnaOut"):
0056 SubsysReco(name)
0057   , clusters_Towers(nullptr)
0058   , truth_photon(nullptr)
0059   , conversion(nullptr)
0060   //, caloevalstack(NULL)
0061   , m_eta_center()
0062   , m_phi_center()
0063   , m_tower_energy()
0064 , m_cluster_eta()
0065   , m_cluster_phi()
0066 , m_cluster_e()
0067 , m_cluster_ecore()
0068   , m_cluster_chi2()
0069   , m_cluster_prob()
0070 , m_cluster_nTowers()
0071   , m_photon_E()
0072 , m_photon_eta()
0073 , m_photon_phi()
0074   , m_electron_E()
0075 , m_electron_eta()
0076 , m_electron_phi()
0077   , m_positron_E()
0078 , m_positron_eta()
0079 , m_positron_phi()
0080   , m_vtx_x()
0081   , m_vtx_y()
0082   , m_vtx_z()
0083   , m_vtx_t()
0084   , m_vtx_r()
0085   , m_isConversionEvent()
0086   , Outfile(outName)
0087
0088 {
0089
0090
0091   std::cout << "singlePhotonClusterAna::singlePhotonClusterAna(const std::string &name) Calling ctor" << std::endl;
0092 }
0093
0094 //____________________________________________________________________________..
0095 singlePhotonClusterAna::~singlePhotonClusterAna()
0096 {
0097   std::cout << "singlePhotonClusterAna::~singlePhotonClusterAna() Calling dtor" << std::endl;
0098 }
0099
0100 //____________________________________________________________________________..
0101 int singlePhotonClusterAna::Init(PHCompositeNode *topNode)
0102 {
0103
0104   clusters_Towers = new TTree("TreeClusterTower","Tree for cluster and tower information");
0105   clusters_Towers -> Branch("towerEtaCEnter",&m_eta_center);
0106   clusters_Towers -> Branch("towerPhiCenter",& m_phi_center);
0107   clusters_Towers -> Branch("towerEnergy",&m_tower_energy);
0108   clusters_Towers -> Branch("clusterEta",&m_cluster_eta);
0109   clusters_Towers -> Branch("clusterPhi", &m_cluster_phi);
0110   clusters_Towers -> Branch("clusterE",&m_cluster_e);
0111   clusters_Towers -> Branch("clusterEcore",&m_cluster_ecore);
0112   clusters_Towers -> Branch("clusterChi2",&m_cluster_chi2);
0113   clusters_Towers -> Branch("clusterProb",&m_cluster_prob);
0114   clusters_Towers -> Branch("clusterNTowers",&m_cluster_nTowers);
0115
0116   truth_photon = new TTree("TreeTruthPhoton", "Tree for Truth Photon");
0117   truth_photon -> Branch("photonE",&m_photon_E);
0118   truth_photon -> Branch("photon_phi",&m_photon_phi);
0119   truth_photon -> Branch("photon_eta",&m_photon_eta);
0120
0121   conversion = new TTree("TreeConversion", "Tree for Photon Conversion Info");
0122   conversion -> Branch("electronE",&m_electron_E);
0123   conversion -> Branch("electron_phi",&m_electron_phi);
0124   conversion -> Branch("electron_eta",&m_electron_eta);
0125   conversion -> Branch("positronE",&m_positron_E);
0126   conversion -> Branch("positron_phi",&m_positron_phi);
0127   conversion -> Branch("positron_eta",&m_positron_eta);
0128   conversion -> Branch("vtx_x",&m_vtx_x);
0129   conversion -> Branch("vtx_y",&m_vtx_y);
0130   conversion -> Branch("vtx_z",&m_vtx_z);
0131   conversion -> Branch("vtx_t",&m_vtx_t);
0132   conversion -> Branch("vtx_r",&m_vtx_r);
0133   conversion -> Branch("isConversionEvent",&m_isConversionEvent);
0134
0135   //so that the histos actually get written out
0136   Fun4AllServer *se = Fun4AllServer::instance();
0137   se -> Print("NODETREE");
0138   hm = new Fun4AllHistoManager("MYHISTOS");
0139
0140   se -> registerHistoManager(hm);
0141
0142   se -> registerHisto(truth_photon -> GetName(), truth_photon);
0143   se -> registerHisto(clusters_Towers -> GetName(), clusters_Towers);
0144
0145   out = new TFile(Outfile.c_str(),"RECREATE");
0146
0147   std::cout << "singlePhotonClusterAna::Init(PHCompositeNode *topNode) Initializing" << std::endl;
0148   return Fun4AllReturnCodes::EVENT_OK;
0149 }
0150
0151 //____________________________________________________________________________..
0152 int singlePhotonClusterAna::InitRun(PHCompositeNode *topNode)
0153 {
0154   std::cout << "singlePhotonClusterAna::InitRun(PHCompositeNode *topNode) Initializing for Run XXX" << std::endl;
0155
0156   return Fun4AllReturnCodes::EVENT_OK;
0157 }
0158
0159 //____________________________________________________________________________..
0160 int singlePhotonClusterAna::process_event(PHCompositeNode *topNode)
0161 {
0162
0163   //Information on clusters
0164   RawClusterContainer *clusterContainer = findNode::getClass<RawClusterContainer>(topNode,"CLUSTER_POS_COR_CEMC");
0165   //RawClusterContainer *clusterContainer = findNode::getClass<RawClusterContainer>(topNode,"CLUSTER_CEMC");
0166   if(!clusterContainer)
0167     {
0168       std::cout << PHWHERE << "singlePhotonClusterAna::process_event - Fatal Error - CLUSTER_CEMC node is missing. " << std::endl;
0169
0170       return 0;
0171     }
0172
0173
0174   /* //Vertex information */
0175   GlobalVertexMap *vtxContainer = findNode::getClass<GlobalVertexMap>(topNode,"GlobalVertexMap");
0176   if (!vtxContainer)
0177     {
0178       std::cout << PHWHERE << "singlePhotonClusterAna::process_event - Fatal Error - GlobalVertexMap node is missing. Please turn on the do_global flag in the main macro in order to reconstruct the global vertex." << std::endl;
0179       assert(vtxContainer);  // force quit
0180
0181       return 0;
0182     }
0183
0184   if (vtxContainer->empty())
0185     {
0186       std::cout << PHWHERE << "singlePhotonClusterAna::process_event - Fatal Error - GlobalVertexMap node is empty. Please turn on the do_global flag in the main macro in order to reconstruct the global vertex." << std::endl;
0187       return 0;
0188     }
0189
0190   //More vertex information
0191   GlobalVertex *vtx = vtxContainer->begin()->second;
0192   if(!vtx)
0193     {
0194
0195       std::cout << PHWHERE << "singlePhotonClusterAna::process_event Could not find vtx from vtxContainer"  << std::endl;
0196       return Fun4AllReturnCodes::ABORTEVENT;
0197     }
0198
0199   //Tower geometry node for location information
0200   RawTowerGeomContainer *towergeom = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_CEMC");
0201   if (!towergeom)
0202     {
0203       std::cout << PHWHERE << "singlePhotonClusterAna::process_event Could not find node TOWERGEOM_CEMC"  << std::endl;
0204       return Fun4AllReturnCodes::ABORTEVENT;
0205     }
0206
0207   //Raw tower information
0208   RawTowerContainer *towerContainer = findNode::getClass<RawTowerContainer>(topNode,"TOWER_CALIB_CEMC");
0209   if(!towerContainer)
0210     {
0211       std::cout << PHWHERE << "singlePhotonClusterAna::process_event Could not find node TOWER_CALIB_CEMC"  << std::endl;
0212       return Fun4AllReturnCodes::ABORTEVENT;
0213     }
0214
0215   //truth particle information
0216   PHG4TruthInfoContainer *truthinfo = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
0217   if(!truthinfo)
0218     {
0219       std::cout << PHWHERE << "singlePhotonClusterAna::process_event Could not find node G4TruthInfo"  << std::endl;
0220       return Fun4AllReturnCodes::ABORTEVENT;
0221     }
0222
0223
0224   float photonEta = -99999;
0225   float photonEtaMax = 1.1;
0226   TLorentzVector photon;
0227   //Now we go and find our truth photon and just take its eta coordinate for now
0228   PHG4TruthInfoContainer::Range truthRange = truthinfo->GetPrimaryParticleRange();
0229   PHG4TruthInfoContainer::ConstIterator truthIter;
0230
0231   PHG4Particle *truthPar = NULL;
0232
0233   for(truthIter = truthRange.first; truthIter != truthRange.second; truthIter++)
0234     {
0235       truthPar = truthIter->second;
0236
0237       // photon pid is 22
0238       if(truthPar -> get_pid() == 22 && truthPar -> get_parent_id() == 0)
0239     {
0240       photonEta = getEta(truthPar);
0241       if(abs(photonEta) >= photonEtaMax)
0242         {
0243           return 0;
0244         }
0245       /* std::cout << "\n\nGreg info:\n\nFound truth photon.\n"; */
0246       /* std::cout << "Truth photon E = " << truthPar->get_e() << "; eta = " << getEta(truthPar) << "\n"; */
0247       /* std::cout << "truthPar = " << truthPar << "\n"; */
0248       /* std::cout << "truthPar->get_primary_id() = " << truthPar->get_primary_id() << "\n"; */
0249       /* std::cout << "truthPar->get_pid() = " << truthPar->get_pid() << "\n"; */
0250       /* std::cout << "truthPar->get_vtx_id() = " << truthPar->get_vtx_id() << "\n"; */
0251       /* std::cout << "truthPar->get_track_id() = " << truthPar->get_track_id() << "\n"; */
0252       photon.SetPxPyPzE(truthPar -> get_px(), truthPar -> get_py(), truthPar -> get_pz(), truthPar -> get_e());
0253     }
0254     }
0255
0256   // Check the daughter particles of truth photon
0257   int decay_vtx_idx = 0;
0258   int decay_vtx_idx1 = 0; // first decay particle vertex
0259   int decay_vtx_idx2 = 0; // second decay particle vertex; SHOULD be the same!
0260   bool first_decay_particle = false;
0261   bool second_decay_particle = false;
0262   bool decay_electron = false;
0263   bool decay_positron = false;
0264   int n_children = 0;
0265   TLorentzVector electron;
0266   TLorentzVector positron;
0267   bool isGoodEvent = true;
0268
0269   PHG4TruthInfoContainer::Range truthRangeDecay1 = truthinfo->GetSecondaryParticleRange();
0270   PHG4TruthInfoContainer::ConstIterator truthIterDecay1;
0271
0272   // First just check that the truth photon has exactly 2 daughters
0273   // and they share the same vertex
0274   for(truthIterDecay1 = truthRangeDecay1.first; truthIterDecay1 != truthRangeDecay1.second; truthIterDecay1++)
0275     {
0276       PHG4Particle *decay = truthIterDecay1 -> second;
0277
0278       int dumparentid = decay -> get_parent_id();
0279       if(dumparentid == 1) {
0280     n_children++;
0281     if (! first_decay_particle) {
0282         first_decay_particle = true;
0283         decay_vtx_idx1 = decay->get_vtx_id();
0284     }
0285     else {
0286         second_decay_particle = true;
0287         decay_vtx_idx2 = decay->get_vtx_id();
0288     }
0289     /* std::cout << "\nDecay particle E = " << decay->get_e() << "; eta = " << getEta(decay) << "\n"; */
0290     /* std::cout << "decay = " << decay << "\n"; */
0291     /* std::cout << "decay->get_primary_id() = " << decay->get_primary_id() << "\n"; */
0292     /* std::cout << "decay->get_pid() = " << decay->get_pid() << "\n"; */
0293     /* std::cout << "decay->get_vtx_id() = " << decay->get_vtx_id() << "\n"; */
0294     /* std::cout << "decay->get_track_id() = " << decay->get_track_id() << "\n"; */
0295       }
0296     }
0297
0298   if (n_children != 2) isGoodEvent = false;
0299   if (!(first_decay_particle && second_decay_particle)) isGoodEvent = false;
0300   if (decay_vtx_idx1 != decay_vtx_idx2) isGoodEvent = false;
0301
0302   // Now get decay electron and positron kinematics
0303   for(truthIterDecay1 = truthRangeDecay1.first; truthIterDecay1 != truthRangeDecay1.second; truthIterDecay1++)
0304     {
0305       PHG4Particle *decay = truthIterDecay1 -> second;
0306
0307       int dumtruthpid = decay -> get_pid();
0308       int dumparentid = decay -> get_parent_id();
0309       if(dumparentid == 1) {
0310     if (dumtruthpid == 11) {
0311         // electron
0312         decay_electron = true;
0313         electron.SetPxPyPzE(decay -> get_px(), decay -> get_py(), decay -> get_pz(), decay -> get_e());
0314     }
0315     if (dumtruthpid == -11) {
0316         // positron
0317         decay_positron = true;
0318         positron.SetPxPyPzE(decay -> get_px(), decay -> get_py(), decay -> get_pz(), decay -> get_e());
0319     }
0320       }
0321     }
0322
0323   if (!(decay_electron && decay_positron)) isGoodEvent = false;
0324
0325   // Conversion vertex information
0326   /* std::cout << "\nVertex information:\n"; */
0327   /* PHG4VtxPoint* photon_vtx = truthinfo->GetVtx(photon_vxt_idx); */
0328   /* std::cout << "\ntruthinfo->is_primary_vtx(photon_vtx) = " << truthinfo->is_primary_vtx(photon_vtx) << "\n"; */
0329   /* std::cout << "photon_vtx->get_id() = " << photon_vtx->get_id() << "\n"; */
0330   /* std::cout << "photon_vtx (x,y,z,t) = (" << photon_vtx->get_x() << ", " << photon_vtx->get_y() << ", " << photon_vtx->get_z() << ", " << photon_vtx->get_t() << ")\n"; */
0331   if (isGoodEvent) {
0332       m_photon_E.push_back(photon.Energy());
0333       m_photon_eta.push_back(photon.PseudoRapidity());
0334       m_photon_phi.push_back(photon.Phi());
0335       m_electron_E.push_back(electron.Energy());
0336       m_electron_eta.push_back(electron.PseudoRapidity());
0337       m_electron_phi.push_back(electron.Phi());
0338       m_positron_E.push_back(positron.Energy());
0339       m_positron_eta.push_back(positron.PseudoRapidity());
0340       m_positron_phi.push_back(positron.Phi());
0341
0342       decay_vtx_idx = decay_vtx_idx1;
0343       PHG4VtxPoint* decay_vtx = truthinfo->GetVtx(decay_vtx_idx);
0344       m_vtx_x.push_back(decay_vtx->get_x());
0345       m_vtx_y.push_back(decay_vtx->get_y());
0346       m_vtx_z.push_back(decay_vtx->get_z());
0347       m_vtx_t.push_back(decay_vtx->get_t());
0348       float vtx_x = decay_vtx->get_x();
0349       float vtx_y = decay_vtx->get_y();
0350       float vtx_r = sqrt(vtx_x*vtx_x + vtx_y*vtx_y);
0351       m_vtx_r.push_back(vtx_r);
0352       bool isConversionEvent = false;
0353       float EMCal_inner_radius = 95; // nominal radius is 97.5cm
0354       if (vtx_r < EMCal_inner_radius) isConversionEvent = true;
0355       m_isConversionEvent.push_back(isConversionEvent);
0356       /* std::cout << "\ntruthinfo->is_primary_vtx(decay_vtx) = " << truthinfo->is_primary_vtx(decay_vtx) << "\n"; */
0357       /* std::cout << "decay_vtx->get_id() = " << decay_vtx->get_id() << "\n"; */
0358       /* std::cout << "decay_vtx (x,y,z,t) = (" << decay_vtx->get_x() << ", " << decay_vtx->get_y() << ", " << decay_vtx->get_z() << ", " << decay_vtx->get_t() << ")\n"; */
0359   }
0360   else {
0361       std::cout << "\nGreg info:\nBad conversion event\n";
0362       std::cout << "n_children = " << n_children << "\n";
0363       std::cout << "first_decay_particle = " << first_decay_particle << "; second_decay_particle = " << second_decay_particle << "\n";
0364       std::cout << "decay_vtx_idx1 = " << decay_vtx_idx1 << "; decay_vtx_idx2 = " << decay_vtx_idx2 << "\n";
0365       std::cout << "decay_electron = " << decay_electron << "; decay_positron = " << decay_positron << "\n";
0366       return 0;
0367   }
0368
0369   // Anthony's code for pi0 decays
0370   /*
0371
0372   int firstphotonflag = 0;
0373   int firstfirstphotonflag = 0;
0374   int secondphotonflag = 0;
0375
0376   //int secondsecondphotonflag = 0;
0377
0378   PHG4TruthInfoContainer::Range truthRangeDecay1 = truthinfo->GetSecondaryParticleRange();
0379   PHG4TruthInfoContainer::ConstIterator truthIterDecay1;
0380
0381
0382   TLorentzVector photon1;
0383   TLorentzVector photon2;
0384   int nParticles = 0;
0385
0386   //loop over all our decay photons.
0387   //Make sure they fall within the desired acceptance
0388   //Toss Dalitz decays
0389   //Retain photon kinematics to determine lead photon for eta binning
0390   for(truthIterDecay1 = truthRangeDecay1.first; truthIterDecay1 != truthRangeDecay1.second; truthIterDecay1++)
0391     {
0392       PHG4Particle *decay = truthIterDecay1 -> second;
0393
0394       int dumtruthpid = decay -> get_pid();
0395       int dumparentid = decay -> get_parent_id();
0396
0397       //if the parent is the pi0 and it's a photon and we haven't marked one yet
0398       if(dumparentid == 1 && dumtruthpid == 22 && !firstphotonflag)
0399     {
0400       if(abs(getEta(decay)) > photonEtaMax)
0401         {
0402           return 0;
0403         }
0404       photon1.SetPxPyPzE(decay -> get_px(), decay -> get_py(), decay -> get_pz(), decay -> get_e());
0405
0406       firstphotonflag = 1;
0407     }
0408
0409       if(dumparentid == 1 && dumtruthpid == 22 && firstphotonflag && firstfirstphotonflag)
0410     {
0411       if(abs(getEta(decay)) > photonEtaMax)
0412         {
0413           return 0;
0414         }
0415       photon2.SetPxPyPzE(decay -> get_px(), decay -> get_py(), decay -> get_pz(), decay -> get_e()) ;
0416
0417       secondphotonflag = 1;
0418     }
0419
0420       //Need this flag to make it skip the first photon slot
0421       if(firstphotonflag) firstfirstphotonflag = 1;
0422       if(dumparentid == 1)nParticles ++;
0423     }
0424
0425   if((!firstphotonflag || !secondphotonflag) && nParticles > 1) //Dalitz
0426     {
0427       return 0;
0428     }
0429   else if((!firstphotonflag || !secondphotonflag) && nParticles == 1) //One photon falls outside simulation acceptance
0430     {
0431       return 0;
0432     }
0433   else if((!firstphotonflag || !secondphotonflag) && nParticles == 0) //Both photons fall outside simulation acceptance
0434     {
0435       return 0;
0436     }
0437   */ // End check on decay photons
0438
0439   //grab all the towers and fill their energies.
0440   RawTowerContainer::ConstRange tower_range = towerContainer -> getTowers();
0441   for(RawTowerContainer::ConstIterator tower_iter = tower_range.first; tower_iter!= tower_range.second; tower_iter++)
0442     {
0443       int phibin = tower_iter -> second -> get_binphi();
0444       int etabin = tower_iter -> second -> get_bineta();
0445       double phi = towergeom -> get_phicenter(phibin);
0446       double eta = towergeom -> get_etacenter(etabin);
0447       double energy = tower_iter -> second -> get_energy();
0448
0449       m_phi_center.push_back(phi);
0450       m_eta_center.push_back(eta);
0451       m_tower_energy.push_back(energy);
0452     }
0453
0454   //This is how we iterate over items in the container.
0455   RawClusterContainer::ConstRange clusterEnd = clusterContainer -> getClusters();
0456   RawClusterContainer::ConstIterator clusterIter;
0457
0458   for(clusterIter = clusterEnd.first; clusterIter != clusterEnd.second; clusterIter++)
0459     {
0460       RawCluster *recoCluster = clusterIter -> second;
0461
0462       CLHEP::Hep3Vector vertex(vtx->get_x(), vtx->get_y(), vtx->get_z());
0463       /* CLHEP::Hep3Vector vertex(0., 0., 0.); */
0464       CLHEP::Hep3Vector E_vec_cluster = RawClusterUtility::GetECoreVec(*recoCluster, vertex);
0465       float clusE = E_vec_cluster.mag();
0466       float clusEcore = recoCluster->get_ecore();
0467       float clus_eta = E_vec_cluster.pseudoRapidity();
0468       float clus_phi = E_vec_cluster.phi();
0469
0470       m_cluster_eta.push_back(clus_eta);
0471       m_cluster_phi.push_back(clus_phi);
0472       m_cluster_e.push_back(clusE);
0473       m_cluster_ecore.push_back(clusEcore);
0474
0475       m_cluster_chi2.push_back(recoCluster -> get_chi2());
0476       m_cluster_prob.push_back(recoCluster -> get_prob());
0477       m_cluster_nTowers.push_back(recoCluster -> getNTowers());
0478     }
0479
0480
0481   clusters_Towers -> Fill();
0482   truth_photon -> Fill();
0483   conversion -> Fill();
0484
0485
0486   return Fun4AllReturnCodes::EVENT_OK;
0487 }
0488
0489 //____________________________________________________________________________..
0490 int singlePhotonClusterAna::ResetEvent(PHCompositeNode *topNode)
0491 {
0492   //std::cout << "singlePhotonClusterAna::ResetEvent(PHCompositeNode *topNode) Resetting internal structures, prepare for next event" << std::endl;
0493
0494   m_eta_center.clear();
0495   m_phi_center.clear();
0496   m_tower_energy.clear();
0497   m_cluster_eta.clear();
0498   m_cluster_phi.clear();
0499   m_cluster_e.clear();
0500   m_cluster_chi2.clear();
0501   m_cluster_prob.clear();
0502   m_cluster_nTowers.clear();
0503   m_photon_E.clear();
0504   m_photon_eta.clear();
0505   m_photon_phi.clear();
0506   m_electron_E.clear();
0507   m_electron_eta.clear();
0508   m_electron_phi.clear();
0509   m_positron_E.clear();
0510   m_positron_eta.clear();
0511   m_positron_phi.clear();
0512   m_vtx_x.clear();
0513   m_vtx_y.clear();
0514   m_vtx_z.clear();
0515   m_vtx_t.clear();
0516   m_vtx_r.clear();
0517   m_isConversionEvent.clear();
0518   return Fun4AllReturnCodes::EVENT_OK;
0519 }
0520
0521 //____________________________________________________________________________..
0522 int singlePhotonClusterAna::EndRun(const int runnumber)
0523 {
0524   std::cout << "singlePhotonClusterAna::EndRun(const int runnumber) Ending Run for Run " << runnumber << std::endl;
0525   return Fun4AllReturnCodes::EVENT_OK;
0526 }
0527
0528 //____________________________________________________________________________..
0529 int singlePhotonClusterAna::End(PHCompositeNode *topNode)
0530 {
0531   std::cout << "singlePhotonClusterAna::End(PHCompositeNode *topNode) This is the End..." << std::endl;
0532   out -> cd();
0533
0534   truth_photon -> Write();
0535   clusters_Towers -> Write();
0536   conversion -> Write();
0537
0538   out -> Close();
0539   delete out;
0540
0541   hm -> dumpHistos(Outfile.c_str(),"UPDATE");
0542
0543
0544   return Fun4AllReturnCodes::EVENT_OK;
0545 }
0546
0547 //____________________________________________________________________________..
0548 int singlePhotonClusterAna::Reset(PHCompositeNode *topNode)
0549 {
0550   std::cout << "singlePhotonClusterAna::Reset(PHCompositeNode *topNode) being Reset" << std::endl;
0551   return Fun4AllReturnCodes::EVENT_OK;
0552 }
0553
0554 //____________________________________________________________________________..
0555 void singlePhotonClusterAna::Print(const std::string &what) const
0556 {
0557   std::cout << "singlePhotonClusterAna::Print(const std::string &what) const Printing info for " << what << std::endl;
0558 }
0559 //____________________________________________________________________________..
0560 float singlePhotonClusterAna::getEta(PHG4Particle *particle)
0561 {
0562   float px = particle -> get_px();
0563   float py = particle -> get_py();
0564   float pz = particle -> get_pz();
0565   float p = sqrt(pow(px,2) + pow(py,2) + pow(pz,2));
0566
0567   return 0.5*log((p+pz)/(p-pz));
0568 }