HCal/SampleFrac/Proto4SampleFrac.C

0001 #include "Proto4SampleFrac.h"
0002
0003 #include <g4main/PHG4Hit.h>
0004 #include <g4main/PHG4HitContainer.h>
0005 #include <g4main/PHG4Particle.h>
0006 #include <g4main/PHG4TruthInfoContainer.h>
0007 #include <g4main/PHG4VtxPoint.h>
0008
0009 #include <calobase/RawTowerContainer.h>
0010 #include <calobase/RawTowerGeomContainer.h>
0011 #include <calobase/RawTower.h>
0012 #include <calobase/RawClusterContainer.h>
0013 #include <calobase/RawCluster.h>
0014
0015 #include <g4eval/CaloEvalStack.h>
0016 #include <g4eval/CaloRawClusterEval.h>
0017 #include <g4eval/CaloRawTowerEval.h>
0018 #include <g4eval/CaloTruthEval.h>
0019 #include <g4eval/SvtxEvalStack.h>
0020
0021 #include <fun4all/SubsysReco.h>
0022 #include <fun4all/Fun4AllServer.h>
0023 #include <fun4all/PHTFileServer.h>
0024 #include <fun4all/Fun4AllReturnCodes.h>
0025
0026 #include <phool/getClass.h>
0027 #include <phool/PHCompositeNode.h>
0028
0029 #include <TString.h>
0030 #include <TFile.h>
0031 #include <TH1F.h>
0032 #include <TH2F.h>
0033 #include <TH3F.h>
0034 #include <TVector3.h>
0035 #include <TLorentzVector.h>
0036 #include <TMath.h>
0037
0038 #include <exception>
0039 #include <stdexcept>
0040 #include <iostream>
0041 #include <vector>
0042 #include <set>
0043 #include <algorithm>
0044 #include <cassert>
0045 #include <cmath>
0046
0047 using namespace std;
0048
0049 Proto4SampleFrac::Proto4SampleFrac(const std::string &calo_name, const std::string &filename)
0050   : SubsysReco("Proto4SampleFrac_" + calo_name),
0051   _is_sim(true),
0052   _filename(filename),
0053   _calo_name(calo_name),
0054   _calo_hit_container(nullptr), _calo_abs_hit_container(nullptr), _truth_container(nullptr)
0055 {
0056   Verbosity(1);
0057 }
0058
0059 Proto4SampleFrac::~Proto4SampleFrac()
0060 {
0061 }
0062
0063 Fun4AllHistoManager* Proto4SampleFrac::get_HistoManager()
0064 {
0065   Fun4AllServer *se = Fun4AllServer::instance();
0066   Fun4AllHistoManager *hm = se->getHistoManager("Proto4SampleFrac_HISTOS");
0067
0068   if (not hm)
0069   {
0070     cout
0071         << "Proto4SampleFrac::get_HistoManager - Making Fun4AllHistoManager Proto4SampleFrac_HISTOS"
0072         << endl;
0073     hm = new Fun4AllHistoManager("Proto4SampleFrac_HISTOS");
0074     se->registerHistoManager(hm);
0075   }
0076
0077   assert(hm);
0078
0079   return hm;
0080 }
0081
0082
0083 string Proto4SampleFrac::get_histo_prefix()
0084 {
0085   return "h_QAG4Sim_" + string(_calo_name);
0086 }
0087
0088 int Proto4SampleFrac::InitRun(PHCompositeNode *topNode)
0089 {
0090   if (Verbosity())
0091     cout << "Proto4SampleFrac::InitRun" << endl;
0092
0093   PHNodeIterator iter(topNode);
0094   PHCompositeNode *dstNode = static_cast<PHCompositeNode *>(iter.findFirst(
0095       "PHCompositeNode", "DST"));
0096   assert(dstNode);
0097
0098   if (!dstNode)
0099   {
0100     std::cerr << PHWHERE << "DST Node missing, doing nothing." << std::endl;
0101     throw runtime_error(
0102         "Failed to find DST node in EmcRawTowerBuilder::CreateNodes");
0103   }
0104
0105   _calo_hit_container = findNode::getClass<PHG4HitContainer>(topNode,
0106       "G4HIT_" + _calo_name);
0107   if (!_calo_hit_container)
0108   {
0109     cout << "Proto4SampleFrac::InitRun - Fatal Error - "
0110       << "unable to find DST node " << "G4HIT_" + _calo_name << endl;
0111     assert(_calo_hit_container);
0112   }
0113
0114   _calo_abs_hit_container = findNode::getClass<PHG4HitContainer>(topNode,
0115       "G4HIT_ABSORBER_" + _calo_name);
0116   if (!_calo_abs_hit_container)
0117   {
0118     cout << "Proto4SampleFrac::InitRun - Fatal Error - "
0119       << "unable to find DST node " << "G4HIT_ABSORBER_" + _calo_name
0120       << endl;
0121     assert(_calo_abs_hit_container);
0122   }
0123
0124   _truth_container = findNode::getClass<PHG4TruthInfoContainer>(topNode,
0125       "G4TruthInfo");
0126   if (!_truth_container)
0127   {
0128     cout << "Proto4SampleFrac::InitRun - Fatal Error - "
0129       << "unable to find DST node " << "G4TruthInfo" << endl;
0130     assert(_truth_container);
0131   }
0132
0133   return Fun4AllReturnCodes::EVENT_OK;
0134 }
0135
0136 int Proto4SampleFrac::End(PHCompositeNode *topNode)
0137 {
0138   cout << "Proto4SampleFrac::End - write to " << _filename << endl;
0139   PHTFileServer::get().cd(_filename);
0140
0141   Fun4AllHistoManager *hm = get_HistoManager();
0142   assert(hm);
0143   for (unsigned int i = 0; i < hm->nHistos(); i++)
0144     hm->getHisto(i)->Write();
0145
0146   //  if (_T_EMCalTrk)
0147   //    _T_EMCalTrk->Write();
0148
0149   return Fun4AllReturnCodes::EVENT_OK;
0150 }
0151
0152 int Proto4SampleFrac::Init(PHCompositeNode *topNode)
0153 {
0154   cout << "Proto4SampleFrac::get_HistoManager - Making PHTFileServer "
0155        << _filename << endl;
0156   PHTFileServer::get().open(_filename, "RECREATE");
0157
0158   Fun4AllHistoManager *hm = get_HistoManager();
0159   assert(hm);
0160
0161   TH1D * h_info = new TH1D(TString(get_histo_prefix()) + "_Normalization", //
0162       TString(_calo_name) + " Normalization;Items;Count", 10, .5, 10.5);
0163   int i = 1;
0164   h_info->GetXaxis()->SetBinLabel(i++, "Event");
0165   h_info->GetXaxis()->SetBinLabel(i++, "G4Hit Active");
0166   h_info->GetXaxis()->SetBinLabel(i++, "G4Hit Absor.");
0167   h_info->GetXaxis()->LabelsOption("v");
0168   hm->registerHisto(h_info);
0169
0170   hm->registerHisto(
0171       new TH2F(TString(get_histo_prefix()) + "_G4Hit_RZ", //
0172     TString(_calo_name) + " RZ projection;Z (cm);R (cm)", 1200, -300, 300,
0173     600, 0, 300));
0174
0175   hm->registerHisto(
0176       new TH2F(TString(get_histo_prefix()) + "_G4Hit_XY_cal", //
0177     TString(_calo_name) + " XY projection;X (cm);Y (cm)", 3000, 50, 350,
0178     2000, -100, 100));
0179     // TString(_calo_name) + " XY projection;X (cm);Y (cm)", 1200, -300, 300,
0180     // 1200, -300, 300));
0181
0182   hm->registerHisto(
0183       new TH2F(TString(get_histo_prefix()) + "_G4Hit_XY_abs", //
0184     TString(_calo_name) + " XY projection;X (cm);Y (cm)", 3000, 50, 350,
0185     2000, -100, 100));
0186
0187   hm->registerHisto(
0188       new TH2F(TString(get_histo_prefix()) + "_G4Hit_Mat_XY_cal", //
0189     TString(_calo_name) + " XY projection;X (cm);Y (cm)", 3000, 50, 350,
0190     2000, -100, 100));
0191
0192   hm->registerHisto(
0193       new TH2F(TString(get_histo_prefix()) + "_G4Hit_Mat_XY_abs", //
0194     TString(_calo_name) + " XY projection;X (cm);Y (cm)", 3000, 50, 350,
0195     2000, -100, 100));
0196
0197   hm->registerHisto(
0198       new TH2F(TString(get_histo_prefix()) + "_G4Hit_LateralTruthProjection", //
0199     TString(_calo_name)
0200     + " shower lateral projection (last primary);Polar direction (cm);Azimuthal direction (cm)",
0201     200, -30, 30, 200, -30, 30));
0202
0203   hm->registerHisto(new TH1F(TString(get_histo_prefix()) + "_G4Hit_SF", //
0204     TString(_calo_name) + " sampling fraction;Sampling fraction", 1000, 0, 1.0));
0205
0206   hm->registerHisto(
0207       new TH1F(TString(get_histo_prefix()) + "_G4Hit_VSF", //
0208     TString(_calo_name)
0209     + " visible sampling fraction;Visible sampling fraction", 1000, 0, 1.0));
0210
0211   TH1F * h =
0212     new TH1F(TString(get_histo_prefix()) + "_G4Hit_HitTime", //
0213     TString(_calo_name)
0214     + " hit time (edep weighting);Hit time - T0 (ns);Geant4 energy density",
0215     1000, 0.5, 10000);
0216   hm->registerHisto(h);
0217
0218   hm->registerHisto(
0219       new TH1F(TString(get_histo_prefix()) + "_G4Hit_FractionTruthEnergy", //
0220     TString(_calo_name)
0221     + " fraction truth energy ;G4 energy (active + absorber) / total truth energy",
0222     1000, 0, 1));
0223
0224   hm->registerHisto(
0225       new TH1F(TString(get_histo_prefix()) + "_G4Hit_FractionEMVisibleEnergy", //
0226     TString(_calo_name)
0227     + " fraction visible energy from EM; visible energy from e^{#pm} / total visible energy",
0228     100, 0, 1));
0229
0230   return Fun4AllReturnCodes::EVENT_OK;
0231 }
0232
0233 int Proto4SampleFrac::process_event(PHCompositeNode *topNode)
0234 {
0235   if (Verbosity() > 2) cout << "Proto4SampleFrac::process_event() entered" << endl;
0236
0237   TH1F *h = nullptr;
0238
0239   Fun4AllHistoManager *hm = get_HistoManager();
0240   assert(hm);
0241
0242   TH1D* h_info = dynamic_cast<TH1D*>(hm->getHisto(get_histo_prefix() + "_Normalization"));
0243   assert(h_info);
0244
0245   // get primary
0246   assert(_truth_container);
0247   PHG4TruthInfoContainer::ConstRange primary_range =
0248     _truth_container->GetPrimaryParticleRange();
0249   double total_primary_energy = 1e-9; //make it zero energy epsilon samll so it can be used for denominator
0250   for (PHG4TruthInfoContainer::ConstIterator particle_iter = primary_range.first;
0251       particle_iter != primary_range.second; ++particle_iter)
0252   {
0253     const PHG4Particle *particle = particle_iter->second;
0254     assert(particle);
0255     total_primary_energy += particle->get_e();
0256   }
0257
0258   assert(not _truth_container->GetMap().empty());
0259   const PHG4Particle * last_primary = _truth_container->GetMap().rbegin()->second;
0260   assert(last_primary);
0261
0262   if (Verbosity() > 2)
0263   {
0264     cout
0265       << "QAG4SimulationCalorimeter::process_event_G4Hit() handle this truth particle"
0266       << endl;
0267     last_primary->identify();
0268   }
0269   const PHG4VtxPoint* primary_vtx = //
0270     _truth_container->GetPrimaryVtx(last_primary->get_vtx_id());
0271   assert(primary_vtx);
0272   if (Verbosity() > 2)
0273   {
0274     cout
0275       << "QAG4SimulationCalorimeter::process_event_G4Hit() handle this vertex"
0276       << endl;
0277     primary_vtx->identify();
0278   }
0279
0280   const double t0 = primary_vtx->get_t();
0281   const TVector3 vertex(primary_vtx->get_x(), primary_vtx->get_y(),
0282       primary_vtx->get_z());
0283
0284   // projection axis
0285   TVector3 axis_proj(last_primary->get_px(), last_primary->get_py(),
0286       last_primary->get_pz());
0287   if (axis_proj.Mag() == 0)
0288     axis_proj.SetXYZ(0, 0, 1);
0289   axis_proj = axis_proj.Unit();
0290
0291   // azimuthal direction axis
0292   TVector3 axis_azimuth = axis_proj.Cross(TVector3(0, 0, 1));
0293   if (axis_azimuth.Mag() == 0)
0294     axis_azimuth.SetXYZ(1, 0, 0);
0295   axis_azimuth = axis_azimuth.Unit();
0296
0297   // polar direction axis
0298   TVector3 axis_polar = axis_proj.Cross(axis_azimuth);
0299   assert(axis_polar.Mag() > 0);
0300   axis_polar = axis_polar.Unit();
0301
0302   double e_calo = 0.0; // active energy deposition
0303   double ev_calo = 0.0; // visible energy
0304   double ea_calo = 0.0; // absorber energy
0305   double ev_calo_em = 0.0; // EM visible energy
0306
0307   if (_calo_hit_container)
0308   {
0309     TH2F * hrz = dynamic_cast<TH2F*>(hm->getHisto(
0310       get_histo_prefix() + "_G4Hit_RZ"));
0311     assert(hrz);
0312
0313     TH2F * hxy_cal = dynamic_cast<TH2F*>(hm->getHisto(
0314       get_histo_prefix() + "_G4Hit_XY_cal"));
0315     assert(hxy_cal);
0316
0317     TH2F * hmat_xy_cal = dynamic_cast<TH2F*>(hm->getHisto(
0318       get_histo_prefix() + "_G4Hit_Mat_XY_cal"));
0319     assert(hmat_xy_cal);
0320
0321     TH1F * ht = dynamic_cast<TH1F*>(hm->getHisto(
0322       get_histo_prefix() + "_G4Hit_HitTime"));
0323     assert(ht);
0324
0325     TH2F * hlat = dynamic_cast<TH2F*>(hm->getHisto(
0326       get_histo_prefix() + "_G4Hit_LateralTruthProjection"));
0327     assert(hlat);
0328
0329     h_info->Fill("G4Hit Active", _calo_hit_container->size());
0330     PHG4HitContainer::ConstRange calo_hit_range =
0331       _calo_hit_container->getHits();
0332     for (PHG4HitContainer::ConstIterator hit_iter = calo_hit_range.first;
0333     hit_iter != calo_hit_range.second; hit_iter++)
0334     {
0335
0336       PHG4Hit *this_hit = hit_iter->second;
0337       assert(this_hit);
0338
0339       e_calo += this_hit->get_edep();
0340       ev_calo += this_hit->get_light_yield();
0341
0342       // EM visible energy that is only associated with electron energy deposition
0343       PHG4Particle* particle = _truth_container->GetParticle(
0344       this_hit->get_trkid());
0345       if (!particle)
0346       {
0347     cout <<__PRETTY_FUNCTION__<<" - Error - this PHG4hit missing particle: "; this_hit -> identify();
0348       }
0349       assert(particle);
0350       if (abs(particle->get_pid()) == 11)
0351     ev_calo_em += this_hit->get_light_yield();
0352
0353       const TVector3 hit(this_hit->get_avg_x(), this_hit->get_avg_y(),
0354       this_hit->get_avg_z());
0355
0356       hrz->Fill(hit.Z(), hit.Perp(), this_hit->get_edep());
0357       hxy_cal->Fill(hit.X(), hit.Y(), this_hit->get_edep());
0358       ht->Fill(this_hit->get_avg_t() - t0, this_hit->get_edep());
0359
0360       const double hit_azimuth = axis_azimuth.Dot(hit - vertex);
0361       const double hit_polar = axis_polar.Dot(hit - vertex);
0362       hlat->Fill(hit_polar, hit_azimuth, this_hit->get_edep());
0363
0364       const TVector3 hit_entry(this_hit->get_x(0), this_hit->get_y(0), this_hit->get_z(0));
0365       const TVector3 hit_exit(this_hit->get_x(1), this_hit->get_y(1), this_hit->get_z(1));
0366       hmat_xy_cal->Fill(hit_entry.X(),hit_entry.Y());
0367       hmat_xy_cal->Fill(hit_exit.X(),hit_exit.Y());
0368     }
0369   }
0370
0371   if (_calo_abs_hit_container)
0372   {
0373     h_info->Fill("G4Hit Absor.", _calo_abs_hit_container->size());
0374
0375     TH2F * hxy_abs = dynamic_cast<TH2F*>(hm->getHisto(
0376       get_histo_prefix() + "_G4Hit_XY_abs"));
0377     assert(hxy_abs);
0378
0379     TH2F * hmat_xy_abs = dynamic_cast<TH2F*>(hm->getHisto(
0380       get_histo_prefix() + "_G4Hit_Mat_XY_abs"));
0381     assert(hmat_xy_abs);
0382
0383     PHG4HitContainer::ConstRange calo_abs_hit_range =
0384       _calo_abs_hit_container->getHits();
0385     for (PHG4HitContainer::ConstIterator hit_iter = calo_abs_hit_range.first;
0386     hit_iter != calo_abs_hit_range.second; hit_iter++)
0387     {
0388
0389       PHG4Hit *this_hit = hit_iter->second;
0390       assert(this_hit);
0391
0392       ea_calo += this_hit->get_edep();
0393
0394       const TVector3 hit(this_hit->get_avg_x(), this_hit->get_avg_y(),
0395       this_hit->get_avg_z());
0396
0397       hxy_abs->Fill(hit.X(), hit.Y(), this_hit->get_edep());
0398
0399       const TVector3 hit_entry(this_hit->get_x(0), this_hit->get_y(0), this_hit->get_z(0));
0400       const TVector3 hit_exit(this_hit->get_x(1), this_hit->get_y(1), this_hit->get_z(1));
0401       hmat_xy_abs->Fill(hit_entry.X(),hit_entry.Y());
0402       hmat_xy_abs->Fill(hit_exit.X(),hit_exit.Y());
0403     }
0404   }
0405
0406   if (Verbosity() > 3)
0407     cout << "QAG4SimulationCalorimeter::process_event_G4Hit::" << _calo_name
0408       << " - SF = " << e_calo / (e_calo + ea_calo + 1e-9) << ", VSF = "
0409       << ev_calo / (e_calo + ea_calo + 1e-9) << endl;
0410
0411   if (e_calo + ea_calo > 0)
0412   {
0413     h = dynamic_cast<TH1F*>(hm->getHisto(get_histo_prefix() + "_G4Hit_SF"));
0414     assert(h);
0415     h->Fill(e_calo / (e_calo + ea_calo));
0416
0417     h = dynamic_cast<TH1F*>(hm->getHisto(get_histo_prefix() + "_G4Hit_VSF"));
0418     assert(h);
0419     h->Fill(ev_calo / (e_calo + ea_calo));
0420   }
0421
0422   h = dynamic_cast<TH1F*>(hm->getHisto(
0423     get_histo_prefix() + "_G4Hit_FractionTruthEnergy"));
0424   assert(h);
0425   h->Fill((e_calo + ea_calo) / total_primary_energy);
0426
0427   if (ev_calo > 0)
0428   {
0429     h = dynamic_cast<TH1F*>(hm->getHisto(
0430       get_histo_prefix() + "_G4Hit_FractionEMVisibleEnergy"));
0431     assert(h);
0432     h->Fill(ev_calo_em / (ev_calo));
0433   }
0434
0435   if (Verbosity() > 3)
0436     cout << "QAG4SimulationCalorimeter::process_event_G4Hit::" << _calo_name
0437       << " - histogram " << h->GetName() << " Get Sum = " << h->GetSum()
0438       << endl;
0439
0440
0441   // at the end, count success events
0442   h_info->Fill("Event", 1);
0443
0444   return Fun4AllReturnCodes::EVENT_OK;
0445 }
0446
0447 pair<int, int>
0448 Proto4SampleFrac::find_max(RawTowerContainer *towers, int cluster_size)
0449 {
0450   const int clus_edge_size = (cluster_size - 1) / 2;
0451   assert(clus_edge_size >= 0);
0452
0453   pair<int, int> max(-1000, -1000);
0454   double max_e = 0;
0455
0456   for (int col = 0; col < n_size; ++col)
0457     for (int row = 0; row < n_size; ++row)
0458     {
0459       double energy = 0;
0460
0461       for (int dcol = col - clus_edge_size; dcol <= col + clus_edge_size;
0462            ++dcol)
0463         for (int drow = row - clus_edge_size; drow <= row + clus_edge_size;
0464              ++drow)
0465         {
0466           if (dcol < 0 or drow < 0)
0467             continue;
0468
0469           RawTower *t = towers->getTower(dcol, drow);
0470           if (t)
0471             energy += t->get_energy();
0472         }
0473
0474       if (energy > max_e)
0475       {
0476         max = make_pair(col, row);
0477         max_e = energy;
0478       }
0479     }
0480
0481   return max;
0482 }