AntiSigma/src/TowerTiming.cc

0001 #include "TowerTiming.h"
0002
0003 #include <calobase/RawTowerGeomContainer.h>
0004 #include <calobase/RawTowerContainer.h>
0005 #include <calobase/RawTower.h>
0006
0007 #include <g4main/PHG4TruthInfoContainer.h>
0008 #include <g4main/PHG4Particle.h>
0009 #include <g4main/PHG4VtxPoint.h>
0010
0011 #include <fun4all/Fun4AllHistoManager.h>
0012
0013 #include <phool/getClass.h>
0014
0015 #include <TFile.h>
0016 #include <TH1.h>
0017 #include <TH2.h>
0018 #include <TNtuple.h>
0019
0020 #include<sstream>
0021
0022 using namespace std;
0023
0024 TowerTiming::TowerTiming(const std::string &name, const std::string &filename):
0025   SubsysReco( name ),
0026   hm(nullptr),
0027   _filename(filename),
0028   ntuptwr(nullptr),
0029   ntupe(nullptr),
0030   outfile(nullptr)
0031 {}
0032
0033 TowerTiming::~TowerTiming()
0034 {
0035   //  delete ntup;
0036   delete hm;
0037 }
0038
0039
0040 int
0041 TowerTiming::Init( PHCompositeNode* )
0042 {
0043   ostringstream hname, htit;
0044   hm = new  Fun4AllHistoManager(Name());
0045   outfile = new TFile(_filename.c_str(), "RECREATE");
0046   ntuptwr = new TNtuple("twr","Towers","detid:phi:eta:edep:t");
0047   ntupe = new TNtuple("truth", "The Absolute Truth", "phi:theta:eta:e:p");
0048   return 0;
0049 }
0050
0051 int
0052 TowerTiming::process_event( PHCompositeNode* topNode )
0053 {
0054   map<int, PHG4Particle*>::const_iterator particle_iter;
0055   PHG4TruthInfoContainer *_truth_container = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
0056
0057   PHG4TruthInfoContainer::ConstRange primary_range =
0058     _truth_container->GetPrimaryParticleRange();
0059   float ntvars[5] = {0};
0060   for (PHG4TruthInfoContainer::ConstIterator particle_iter = primary_range.first;
0061        particle_iter != primary_range.second; ++particle_iter)
0062   {
0063     PHG4Particle *particle = particle_iter->second;
0064     ntvars[0] = atan2(particle->get_py(), particle->get_px());
0065     ntvars[1] = atan(sqrt(particle->get_py()*particle->get_py() +
0066               particle->get_px()*particle->get_px()) /
0067              particle->get_pz());
0068     ntvars[2] = 0.5*log((particle->get_e()+particle->get_pz())/
0069             (particle->get_e()-particle->get_pz()));
0070     ntvars[3] = particle->get_e();
0071     ntvars[4] = particle->get_pz();
0072   }
0073   ntupe->Fill(ntvars);
0074   RawTowerContainer  *towers = nullptr;
0075   // double emc[10] = {0};
0076   // double ih[10] = {0};
0077   // double oh[10] = {0};
0078   // double mag[10] = {0};
0079   // double bh[10] = {0};
0080   // double eall[5] = {0};
0081   string nodename[3] = {"TOWER_SIM_CEMC", "TOWER_SIM_HCALIN", "TOWER_SIM_HCALOUT"};
0082   string geonodename[3] = {"TOWERGEOM_CEMC", "TOWERGEOM_HCALIN", "TOWERGEOM_HCALOUT"};
0083   for (int j=0; j<3;j++)
0084   {
0085     towers = findNode::getClass<RawTowerContainer>(topNode,nodename[j]);
0086     if (towers)
0087     {
0088       RawTowerGeomContainer* towergeom = findNode::getClass<RawTowerGeomContainer>(topNode,geonodename[j]);
0089       if (!towergeom)
0090     {
0091       cout << "no geometry node " << geonodename << endl;
0092       continue;
0093     }
0094       RawTowerContainer::ConstRange twr_range = towers->getTowers();
0095       for ( RawTowerContainer::ConstIterator twr_iter = twr_range.first ; twr_iter !=  twr_range.second; ++twr_iter )
0096       {
0097         double energy = twr_iter->second->get_energy();
0098           int phibin = twr_iter->second->get_binphi();
0099           int etabin = twr_iter->second->get_bineta();
0100           double phi = towergeom->get_phicenter(phibin);
0101           double eta = towergeom->get_etacenter(etabin);
0102                ntuptwr->Fill(j,
0103               phi,
0104               eta,
0105                           energy,
0106                  twr_iter->second->get_time());
0107     // double radius = sqrt(hit_iter->second->get_avg_x() * hit_iter->second->get_avg_x() +
0108     //           hit_iter->second->get_avg_y() * hit_iter->second->get_avg_y() +
0109     //           hit_iter->second->get_avg_z() * hit_iter->second->get_avg_z());
0110 //  double phi =  atan2(hit_iter->second->get_avg_y(),hit_iter->second->get_avg_x());
0111 //  double theta = atan(sqrt(hit_iter->second->get_avg_x() * hit_iter->second->get_avg_x() +
0112 //               hit_iter->second->get_avg_y() * hit_iter->second->get_avg_y()) /
0113 //              hit_iter->second->get_avg_z());
0114 // // handle rollover from pi to -pi
0115 //  double diffphi = phi-ntvars[0];
0116 //  if (diffphi > M_PI)
0117 //  {
0118 //    diffphi -= 2*M_PI;
0119 //  }
0120 //  else if (diffphi < - M_PI)
0121 //  {
0122 //    diffphi += 2*M_PI;
0123 //  }
0124
0125 //  double deltasqrt = sqrt(diffphi*diffphi+(theta-ntvars[1])*(theta-ntvars[1]));
0126 //  double edep = hit_iter->second->get_edep();
0127 //  eall[0] += edep;
0128 //  for (int i=0; i<10; i++)
0129 //  {
0130 //    if (deltasqrt < i*0.1)
0131 //    {
0132 //      emc[i]+=edep;
0133 //    }
0134 //  }
0135 //       }
0136       }
0137     }
0138   }
0139   return 0;
0140 }
0141
0142 int
0143 TowerTiming::End(PHCompositeNode * topNode)
0144 {
0145   outfile->cd();
0146   //  ntup->Write();
0147   ntuptwr->Write();
0148   ntupe->Write();
0149   outfile->Write();
0150   outfile->Close();
0151   delete outfile;
0152   hm->dumpHistos(_filename, "UPDATE");
0153   return 0;
0154 }
0155
0156 void
0157 TowerTiming::AddNode(const std::string &name, const int detid)
0158 {
0159   _node_postfix.insert(name);
0160   _detid[name] = detid;
0161   return;
0162 }