EMCal-analysis/EMCalAna/EMCalLikelihood.C

0001 #include "EMCalLikelihood.h"
0002 #include "EMCalTrk.h"
0003
0004 #include <phool/getClass.h>
0005 #include <fun4all/SubsysReco.h>
0006 #include <fun4all/Fun4AllServer.h>
0007 #include <fun4all/PHTFileServer.h>
0008 #include <phool/PHCompositeNode.h>
0009 #include <fun4all/Fun4AllReturnCodes.h>
0010 #include <fun4all/Fun4AllHistoManager.h>
0011
0012 #include <phool/PHCompositeNode.h>
0013
0014 #include <g4main/PHG4TruthInfoContainer.h>
0015 #include <g4main/PHG4VtxPoint.h>
0016 #include <g4main/PHG4Particle.h>
0017
0018 #include <g4hough/SvtxVertexMap.h>
0019
0020 #include <calobase/RawTowerContainer.h>
0021 #include <calobase/RawTowerGeom.h>
0022 #include <calobase/RawTower.h>
0023 #include <calobase/RawClusterContainer.h>
0024 #include <calobase/RawCluster.h>
0025
0026 #include <g4eval/CaloEvalStack.h>
0027 #include <g4eval/CaloRawClusterEval.h>
0028 #include <g4eval/CaloRawTowerEval.h>
0029 #include <g4eval/CaloTruthEval.h>
0030 #include <g4eval/SvtxEvalStack.h>
0031
0032 #include <TNtuple.h>
0033 #include <TFile.h>
0034 #include <TH1F.h>
0035 #include <TH2F.h>
0036 #include <TVector3.h>
0037 #include <TLorentzVector.h>
0038
0039 #include <exception>
0040 #include <stdexcept>
0041 #include <iostream>
0042 #include <vector>
0043 #include <set>
0044 #include <algorithm>
0045 #include <cassert>
0046 #include <cmath>
0047
0048 using namespace std;
0049
0050 EMCalLikelihood::EMCalLikelihood(const std::string &filename) :
0051     SubsysReco("EMCalLikelihood"), _filename(filename), _ievent(0), _trk(NULL) //
0052         , center_cemc_iphi(NAN), center_cemc_ieta(NAN), //
0053     center_hcalin_iphi(NAN), center_hcalin_ieta(NAN), //
0054     width_cemc_iphi(NAN), width_cemc_ieta(NAN), //
0055     width_hcalin_iphi(NAN), width_hcalin_ieta(NAN), //
0056     h2_Edep_Distribution_e(NULL), h2_Edep_Distribution_pi(NULL), //
0057     h1_ep_Distribution_e(NULL), h1_ep_Distribution_pi(NULL), //
0058     _do_ganging(false), _ganging_size(1, 1)
0059 {
0060
0061 }
0062
0063 EMCalLikelihood::~EMCalLikelihood()
0064 {
0065
0066 }
0067
0068 int
0069 EMCalLikelihood::InitRun(PHCompositeNode *topNode)
0070 {
0071   _ievent = 0;
0072
0073   PHNodeIterator iter(topNode);
0074   PHCompositeNode *dstNode = static_cast<PHCompositeNode*>(iter.findFirst(
0075       "PHCompositeNode", "DST"));
0076   if (!dstNode)
0077     {
0078       std::cerr << PHWHERE << "DST Node missing, doing nothing." << std::endl;
0079       throw runtime_error(
0080           "Failed to find DST node in EmcRawTowerBuilder::CreateNodes");
0081     }
0082
0083     {
0084       _trk = findNode::getClass<EMCalTrk>(dstNode, "EMCalTrk");
0085
0086       if (not _trk)
0087         {
0088
0089           cout << __PRETTY_FUNCTION__ << "::" << Name() << " - Error - "
0090               << "Can not find the EMCalTrk node" << endl;
0091
0092           return Fun4AllReturnCodes::ABORTRUN;
0093
0094         }
0095
0096       assert(_trk);
0097     }
0098
0099   return Fun4AllReturnCodes::EVENT_OK;
0100 }
0101
0102 int
0103 EMCalLikelihood::End(PHCompositeNode *topNode)
0104 {
0105
0106   cout << "EMCalLikelihood::End - write to " << _filename << endl;
0107   PHTFileServer::get().cd(_filename);
0108
0109   Fun4AllHistoManager *hm = get_HistoManager();
0110   assert(hm);
0111   for (unsigned int i = 0; i < hm->nHistos(); i++)
0112     hm->getHisto(i)->Write();
0113
0114   return Fun4AllReturnCodes::EVENT_OK;
0115 }
0116
0117 int
0118 EMCalLikelihood::Init(PHCompositeNode *topNode)
0119 {
0120
0121   cout << "EMCalLikelihood::Init - Making PHTFileServer " << _filename << endl;
0122   PHTFileServer::get().open(_filename, "RECREATE");
0123
0124   Fun4AllHistoManager *hm = get_HistoManager();
0125   assert(hm);
0126
0127   TH2F * h_merge_direction = new TH2F("h_merge_direction",
0128       "h_merge_direction;merge eta shift;merge phi shift", 10, -.5, 9.5, 10,
0129       -.5, 9.5);
0130   hm->registerHisto(h_merge_direction);
0131
0132   if (h2_Edep_Distribution_e)
0133     {
0134       h2_Edep_Distribution_e =
0135           dynamic_cast<TH2 *>(h2_Edep_Distribution_e->Clone(
0136               "h2_Edep_Distribution_e"));
0137       hm->registerHisto(h2_Edep_Distribution_e);
0138
0139       h1_ep_Distribution_e = h2_Edep_Distribution_e->ProjectionX(
0140           "h1_ep_Distribution_e");
0141       assert(h1_ep_Distribution_e);
0142       hm->registerHisto(h1_ep_Distribution_e);
0143
0144       // regulate minimal probability value
0145       double min_prob = 1;
0146       for (int x = 1; x <= h2_Edep_Distribution_e->GetNbinsX(); ++x)
0147         for (int y = 1; y <= h2_Edep_Distribution_e->GetNbinsX(); ++y)
0148           {
0149             const double binc = h2_Edep_Distribution_e->GetBinContent(x, y);
0150             assert(binc >= 0);
0151
0152             if (binc > 0 and binc < min_prob)
0153               min_prob = binc;
0154           }
0155       assert(min_prob > 0);
0156
0157       for (int x = 0; x <= h2_Edep_Distribution_e->GetNbinsX() + 1; ++x)
0158         for (int y = 0; y <= h2_Edep_Distribution_e->GetNbinsX() + 1; ++y)
0159           {
0160             const double binc = h2_Edep_Distribution_e->GetBinContent(x, y);
0161
0162             if (binc == 0)
0163               h2_Edep_Distribution_e->SetBinContent(x, y, min_prob);
0164           }
0165
0166       for (int x = 0; x <= h1_ep_Distribution_e->GetNbinsX() + 1; ++x)
0167         {
0168           const double binc = h1_ep_Distribution_e->GetBinContent(x);
0169
0170           if (binc == 0)
0171             h1_ep_Distribution_e->SetBinContent(x, min_prob);
0172         }
0173     }
0174   if (h2_Edep_Distribution_pi)
0175     {
0176       h2_Edep_Distribution_pi =
0177           dynamic_cast<TH2 *>(h2_Edep_Distribution_pi->Clone(
0178               "h2_Edep_Distribution_pi"));
0179       hm->registerHisto(h2_Edep_Distribution_pi);
0180
0181       h1_ep_Distribution_pi = h2_Edep_Distribution_pi->ProjectionX(
0182           "h1_ep_Distribution_pi");
0183       assert(h1_ep_Distribution_pi);
0184       hm->registerHisto(h1_ep_Distribution_pi);
0185
0186       // regulate minimal probability value
0187       double min_prob = 1;
0188       for (int x = 1; x <= h2_Edep_Distribution_pi->GetNbinsX(); ++x)
0189         for (int y = 1; y <= h2_Edep_Distribution_pi->GetNbinsX(); ++y)
0190           {
0191             const double binc = h2_Edep_Distribution_pi->GetBinContent(x, y);
0192             assert(binc >= 0);
0193
0194             if (binc > 0 and binc < min_prob)
0195               min_prob = binc;
0196           }
0197       assert(min_prob > 0);
0198
0199       for (int x = 0; x <= h2_Edep_Distribution_pi->GetNbinsX() + 1; ++x)
0200         for (int y = 0; y <= h2_Edep_Distribution_pi->GetNbinsX() + 1; ++y)
0201           {
0202             const double binc = h2_Edep_Distribution_pi->GetBinContent(x, y);
0203
0204             if (binc == 0)
0205               h2_Edep_Distribution_pi->SetBinContent(x, y, min_prob);
0206           }
0207
0208       for (int x = 0; x <= h1_ep_Distribution_pi->GetNbinsX() + 1; ++x)
0209         {
0210           const double binc = h1_ep_Distribution_pi->GetBinContent(x);
0211
0212           if (binc == 0)
0213             h1_ep_Distribution_pi->SetBinContent(x, min_prob);
0214         }
0215     }
0216   return Fun4AllReturnCodes::EVENT_OK;
0217 }
0218
0219 int
0220 EMCalLikelihood::process_event(PHCompositeNode *topNode)
0221 {
0222 //  const double significand = _ievent / TMath::Power(10, (int) (log10(_ievent)));
0223
0224 //  if (fmod(significand, 1.0) == 0 && significand <= 10)
0225 //    cout << "EMCalLikelihood::process_event - " << _ievent << endl;
0226
0227   _ievent++;
0228
0229   if (_do_ganging)
0230     ApplyEMCalGanging(_trk);
0231
0232   UpdateEnergyDeposition(_trk);
0233
0234   if (h2_Edep_Distribution_e and h2_Edep_Distribution_pi)
0235     UpdateEnergyDepositionLikelihood(_trk);
0236
0237   return Fun4AllReturnCodes::EVENT_OK;
0238 }
0239
0240 Fun4AllHistoManager *
0241 EMCalLikelihood::get_HistoManager()
0242 {
0243
0244   Fun4AllServer *se = Fun4AllServer::instance();
0245   Fun4AllHistoManager *hm = se->getHistoManager("EMCalLikelihood_HISTOS");
0246
0247   if (not hm)
0248     {
0249       cout
0250           << "EMCalLikelihood::get_HistoManager - Making Fun4AllHistoManager EMCalLikelihood_HISTOS"
0251           << endl;
0252       hm = new Fun4AllHistoManager("EMCalLikelihood_HISTOS");
0253       se->registerHistoManager(hm);
0254     }
0255
0256   assert(hm);
0257
0258   return hm;
0259 }
0260
0261 void
0262 EMCalLikelihood::ApplyEMCalGanging(EMCalTrk * trk)
0263 {
0264   assert(trk);
0265   assert(_ganging_size.first > 1 or _ganging_size.second > 1);
0266   assert(_ganging_size.first > 0 and _ganging_size.second > 0);
0267
0268   static bool once = true;
0269   if (once)
0270     {
0271       once = false;
0272
0273       cout << "EMCalLikelihood::ApplyEMCalGanging - apply EMCal ganging "
0274           << _ganging_size.first << "x" << _ganging_size.second << endl;
0275     }
0276
0277   // estimate an eta-phi bin
0278   TVector3 vertex(trk->gvx, trk->gvy, trk->gvz);
0279   TVector3 momentum(trk->gpx, trk->gpy, trk->gpz);
0280   const double radius = 100; // center of CEMC.
0281   const double eta_bin_cm = 2.5;
0282   const double phi_bin_rad = 0.025;
0283
0284   assert(momentum.Pt() > 0);
0285
0286   TVector3 proj = vertex + radius / momentum.Pt() * momentum;
0287   const int eta_bin = floor(proj.Z() / eta_bin_cm);
0288   const int phi_bin = floor(proj.Phi() / phi_bin_rad);
0289
0290   const int merge_plus_eta = (eta_bin % _ganging_size.first);
0291   const bool merge_plus_phi = (phi_bin % _ganging_size.second);
0292
0293   Fun4AllHistoManager *hm = get_HistoManager();
0294   assert(hm);
0295
0296   TH2F * h_merge_direction = (TH2F *) hm->getHisto("h_merge_direction");
0297   assert(h_merge_direction);
0298
0299   h_merge_direction->Fill(merge_plus_eta, merge_plus_phi);
0300
0301   for (int ieta = 0; ieta < trk->Max_N_Tower; ++ieta)
0302     {
0303       for (int iphi = 0; iphi < trk->Max_N_Tower; ++iphi)
0304         {
0305           bool out_of_edge = false;
0306
0307           double cemc_ieta = 0;
0308           double cemc_iphi = 0;
0309           double cemc_energy = 0;
0310
0311           for (unsigned int dieta = 0; dieta < _ganging_size.first; ++dieta)
0312             {
0313               const unsigned int fetch_eta = _ganging_size.first * ieta + dieta + merge_plus_eta;
0314
0315               if ( fetch_eta>= trk->Max_N_Tower or out_of_edge)
0316                 {
0317                   out_of_edge = true;
0318                   break;
0319                 }
0320
0321               for (unsigned int diphi = 0; diphi < _ganging_size.second; ++diphi)
0322                 {
0323                   const unsigned int fetch_phi = _ganging_size.first * iphi + diphi + merge_plus_phi;
0324
0325                   if (fetch_phi>= trk->Max_N_Tower or out_of_edge)
0326                     {
0327                       out_of_edge = true;
0328                       break;
0329                     }
0330
0331                   cemc_ieta += trk->cemc_ieta[fetch_eta][fetch_phi];
0332                   cemc_iphi += trk->cemc_iphi[fetch_eta][fetch_phi];
0333                   cemc_energy += trk->cemc_energy[fetch_eta][fetch_phi];
0334                 }
0335             }
0336
0337           if (out_of_edge)
0338             {
0339               trk->cemc_ieta[ieta][iphi] = -9999;
0340               trk->cemc_iphi[ieta][iphi] = -9999;
0341               trk->cemc_energy[ieta][iphi] = -0;
0342             }
0343           else
0344             {
0345               trk->cemc_ieta[ieta][iphi] = cemc_ieta / _ganging_size.first
0346                   / _ganging_size.second;
0347               trk->cemc_iphi[ieta][iphi] = cemc_iphi / _ganging_size.first
0348                   / _ganging_size.second;
0349               trk->cemc_energy[ieta][iphi] = cemc_energy;
0350             }
0351         }
0352     }
0353 }
0354
0355 void
0356 EMCalLikelihood::UpdateEnergyDeposition(EMCalTrk * trk)
0357 {
0358   assert(trk);
0359
0360   trk->cemc_sum_energy = 0;
0361   trk->cemc_sum_n_tower = 0;
0362   trk->hcalin_sum_energy = 0;
0363   trk->hcalin_sum_n_tower = 0;
0364
0365   for (int ieta = 0; ieta < trk->Max_N_Tower; ++ieta)
0366     {
0367       for (int iphi = 0; iphi < trk->Max_N_Tower; ++iphi)
0368         {
0369
0370           //cemc
0371           const double cemc_radius2 = pow(
0372               (trk->cemc_ieta[ieta][iphi] - center_cemc_ieta) / width_cemc_ieta,
0373               2.)
0374               + pow(
0375                   (trk->cemc_iphi[ieta][iphi] - center_cemc_iphi)
0376                       / width_cemc_iphi, 2.);
0377           trk->cemc_radius[ieta][iphi] = sqrt(cemc_radius2);
0378
0379           if (cemc_radius2 <= 1.0)
0380             {
0381               trk->cemc_sum_energy += trk->cemc_energy[ieta][iphi];
0382               trk->cemc_sum_n_tower++;
0383             }
0384
0385           // hcalin
0386           const double hcalin_radius2 = pow(
0387               (trk->hcalin_ieta[ieta][iphi] - center_hcalin_ieta)
0388                   / width_hcalin_ieta, 2.)
0389               + pow(
0390                   (trk->hcalin_iphi[ieta][iphi] - center_hcalin_iphi)
0391                       / width_hcalin_iphi, 2.);
0392           trk->hcalin_radius[ieta][iphi] = sqrt(hcalin_radius2);
0393
0394           if (hcalin_radius2 <= 1.0)
0395             {
0396               trk->hcalin_sum_energy += trk->hcalin_energy[ieta][iphi];
0397               trk->hcalin_sum_n_tower++;
0398             }
0399         }
0400     }
0401 }
0402
0403 void
0404 EMCalLikelihood::UpdateEnergyDepositionLikelihood(EMCalTrk * trk)
0405 {
0406   assert(trk);
0407   assert(h2_Edep_Distribution_e);
0408   assert(h2_Edep_Distribution_pi);
0409   assert(h1_ep_Distribution_e);
0410   assert(h1_ep_Distribution_pi);
0411
0412   assert(
0413       h2_Edep_Distribution_e->GetNbinsX() == h2_Edep_Distribution_pi->GetNbinsX());
0414   assert(
0415       h2_Edep_Distribution_e->GetNbinsY() == h2_Edep_Distribution_pi->GetNbinsY());
0416   assert(
0417       h2_Edep_Distribution_e->GetNbinsX() == h1_ep_Distribution_e->GetNbinsX());
0418   assert(
0419       h1_ep_Distribution_pi->GetNbinsX() == h2_Edep_Distribution_pi->GetNbinsX());
0420
0421   const int binx = h2_Edep_Distribution_e->GetXaxis()->FindBin(trk->get_ep());
0422   const int biny = h2_Edep_Distribution_e->GetYaxis()->FindBin(
0423       trk->hcalin_sum_energy);
0424
0425     {
0426       double prob_e = h2_Edep_Distribution_e->GetBinContent(binx, biny);
0427       if (!prob_e > 0)
0428         {
0429           cout << __PRETTY_FUNCTION__ << Name()
0430               << " - Error - invalid likelihood value prob_e = " << prob_e
0431               << " @ bin " << binx << ", " << biny << endl;
0432         }
0433       assert(prob_e > 0);
0434
0435       double prob_pi = h2_Edep_Distribution_pi->GetBinContent(binx, biny);
0436       if (!prob_pi > 0)
0437         {
0438           cout << __PRETTY_FUNCTION__ << Name()
0439               << " - Error - invalid likelihood value prob_pi = " << prob_pi
0440               << " @ bin " << binx << ", " << biny << endl;
0441         }
0442       assert(prob_pi > 0);
0443
0444       trk->ll_edep_e = log(prob_e);
0445       trk->ll_edep_h = log(prob_pi);
0446     }
0447
0448     {
0449       double prob_e = h1_ep_Distribution_e->GetBinContent(binx);
0450       assert(prob_e > 0);
0451
0452       double prob_pi = h1_ep_Distribution_pi->GetBinContent(binx);
0453       assert(prob_pi > 0);
0454
0455       trk->ll_ep_e = log(prob_e);
0456       trk->ll_ep_h = log(prob_pi);
0457     }
0458 }