tutorial/ExampleAnalysisModule/ExampleAnalysisModule.C

0001 #include "ExampleAnalysisModule.h"
0002
0003 #include <calobase/RawTowerContainer.h>
0004 #include <ffaobjects/EventHeader.h>
0005 #include <pdbcalbase/PdbParameterMap.h>
0006 #include <phparameter/PHParameters.h>
0007 #include <prototype4/RawTower_Prototype4.h>
0008 #include <prototype4/RawTower_Temperature.h>
0009
0010 #include <Event/EventTypes.h>
0011 #include <fun4all/Fun4AllHistoManager.h>
0012 #include <fun4all/Fun4AllReturnCodes.h>
0013 #include <fun4all/Fun4AllServer.h>
0014 #include <fun4all/PHTFileServer.h>
0015 #include <fun4all/SubsysReco.h>
0016 #include <phool/PHCompositeNode.h>
0017 #include <phool/getClass.h>
0018
0019 #include <phool/PHCompositeNode.h>
0020
0021 #include <g4main/PHG4Particle.h>
0022 #include <g4main/PHG4TruthInfoContainer.h>
0023 #include <g4main/PHG4VtxPoint.h>
0024
0025 #include <TFile.h>
0026 #include <TH1F.h>
0027 #include <TH2F.h>
0028 #include <TLorentzVector.h>
0029 #include <TNtuple.h>
0030 #include <TVector3.h>
0031
0032 #include <algorithm>
0033 #include <cassert>
0034 #include <cmath>
0035 #include <exception>
0036 #include <iostream>
0037 #include <set>
0038 #include <sstream>
0039 #include <stdexcept>
0040 #include <vector>
0041
0042 using namespace std;
0043
0044 ClassImp(ExampleAnalysisModule::Eval_Cluster);
0045 ClassImp(ExampleAnalysisModule::Eval_Run);
0046
0047 ExampleAnalysisModule::ExampleAnalysisModule(const std::string &filename)
0048   : SubsysReco("ExampleAnalysisModule")
0049   , _is_sim(false)
0050   , _filename(filename)
0051   , _ievent(
0052         0)
0053 {
0054   verbosity = 1;
0055
0056   _eval_run.reset();
0057   _eval_5x5_CEMC.reset();
0058 }
0059
0060 ExampleAnalysisModule::~ExampleAnalysisModule()
0061 {
0062 }
0063
0064 Fun4AllHistoManager *
0065 ExampleAnalysisModule::get_HistoManager()
0066 {
0067   Fun4AllServer *se = Fun4AllServer::instance();
0068   Fun4AllHistoManager *hm = se->getHistoManager("ExampleAnalysisModule_HISTOS");
0069
0070   if (not hm)
0071   {
0072     cout
0073         << "ExampleAnalysisModule::get_HistoManager - Making Fun4AllHistoManager ExampleAnalysisModule_HISTOS"
0074         << endl;
0075     hm = new Fun4AllHistoManager("ExampleAnalysisModule_HISTOS");
0076     se->registerHistoManager(hm);
0077   }
0078
0079   assert(hm);
0080
0081   return hm;
0082 }
0083
0084 int ExampleAnalysisModule::InitRun(PHCompositeNode *topNode)
0085 {
0086   if (verbosity)
0087     cout << "ExampleAnalysisModule::InitRun" << endl;
0088
0089   _ievent = 0;
0090
0091   PHNodeIterator iter(topNode);
0092   PHCompositeNode *dstNode = static_cast<PHCompositeNode *>(iter.findFirst(
0093       "PHCompositeNode", "DST"));
0094   if (!dstNode)
0095   {
0096     std::cerr << PHWHERE << "DST Node missing, doing nothing." << std::endl;
0097     throw runtime_error(
0098         "Failed to find DST node in EmcRawTowerBuilder::CreateNodes");
0099   }
0100
0101   return Fun4AllReturnCodes::EVENT_OK;
0102 }
0103
0104 int ExampleAnalysisModule::End(PHCompositeNode *topNode)
0105 {
0106   cout << "ExampleAnalysisModule::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   //  if (_T_EMCalTrk)
0115   //    _T_EMCalTrk->Write();
0116
0117   return Fun4AllReturnCodes::EVENT_OK;
0118 }
0119
0120 int ExampleAnalysisModule::Init(PHCompositeNode *topNode)
0121 {
0122   _ievent = 0;
0123
0124   cout << "ExampleAnalysisModule::get_HistoManager - Making PHTFileServer "
0125        << _filename << endl;
0126   PHTFileServer::get().open(_filename, "RECREATE");
0127
0128   Fun4AllHistoManager *hm = get_HistoManager();
0129   assert(hm);
0130
0131   //! Histogram of Cherenkov counters
0132   TH2F *hCheck_Cherenkov = new TH2F("hCheck_Cherenkov", "hCheck_Cherenkov",
0133                                     110, -200, 20000, 5, .5, 5.5);
0134   hCheck_Cherenkov->GetYaxis()->SetBinLabel(1, "C1");
0135   hCheck_Cherenkov->GetYaxis()->SetBinLabel(2, "C2 in");
0136   hCheck_Cherenkov->GetYaxis()->SetBinLabel(3, "C2 out");
0137   hCheck_Cherenkov->GetYaxis()->SetBinLabel(4, "C2 sum");
0138   hCheck_Cherenkov->GetXaxis()->SetTitle("Amplitude");
0139   hCheck_Cherenkov->GetYaxis()->SetTitle("Cherenkov Channel");
0140   hm->registerHisto(hCheck_Cherenkov);
0141
0142   //! Envent nomalization
0143   TH1F *hNormalization = new TH1F("hNormalization", "hNormalization", 10, .5,
0144                                   10.5);
0145   hNormalization->GetXaxis()->SetBinLabel(1, "ALL");
0146   hNormalization->GetXaxis()->SetBinLabel(2, "C2-e");
0147   hNormalization->GetXaxis()->SetBinLabel(3, "trigger_veto_pass");
0148   hNormalization->GetXaxis()->SetBinLabel(4, "valid_hodo_h");
0149   hNormalization->GetXaxis()->SetBinLabel(5, "valid_hodo_v");
0150   hNormalization->GetXaxis()->SetBinLabel(6, "good_e");
0151   hNormalization->GetXaxis()->SetBinLabel(6, "good_anti_e");
0152   hNormalization->GetXaxis()->SetTitle("Cuts");
0153   hNormalization->GetYaxis()->SetTitle("Event count");
0154   hm->registerHisto(hNormalization);
0155
0156   hm->registerHisto(new TH1F("hCheck_Veto", "hCheck_Veto", 1000, -500, 1500));
0157   hm->registerHisto(
0158       new TH1F("hCheck_Hodo_H", "hCheck_Hodo_H", 1000, -500, 1500));
0159   hm->registerHisto(
0160       new TH1F("hCheck_Hodo_V", "hCheck_Hodo_V", 1000, -500, 1500));
0161
0162   hm->registerHisto(new TH1F("hBeam_Mom", "hBeam_Mom", 1200, -120, 120));
0163
0164   // help index files with TChain
0165   TTree *T = new TTree("T", "T");
0166   assert(T);
0167   hm->registerHisto(T);
0168
0169   T->Branch("info", &_eval_run);
0170   T->Branch("clus_5x5_CEMC", &_eval_5x5_CEMC);
0171
0172   return Fun4AllReturnCodes::EVENT_OK;
0173 }
0174
0175 int ExampleAnalysisModule::process_event(PHCompositeNode *topNode)
0176 {
0177   if (verbosity > 2)
0178     cout << "ExampleAnalysisModule::process_event() entered" << endl;
0179
0180   // init eval objects
0181   _eval_run.reset();
0182   _eval_5x5_CEMC.reset();
0183
0184   Fun4AllHistoManager *hm = get_HistoManager();
0185   assert(hm);
0186
0187   if (not _is_sim)
0188   {
0189     PdbParameterMap *info = findNode::getClass<PdbParameterMap>(topNode,
0190                                                                 "RUN_INFO");
0191
0192     assert(info);
0193
0194     PHParameters run_info_copy("RunInfo");
0195     run_info_copy.FillFrom(info);
0196
0197     _eval_run.beam_mom = run_info_copy.get_double_param("beam_MTNRG_GeV");
0198
0199     TH1F *hBeam_Mom = dynamic_cast<TH1F *>(hm->getHisto("hBeam_Mom"));
0200     assert(hBeam_Mom);
0201
0202     hBeam_Mom->Fill(_eval_run.beam_mom);
0203
0204     _eval_run.beam_2CH_mm = run_info_copy.get_double_param("beam_2CH_mm");
0205     _eval_run.beam_2CV_mm = run_info_copy.get_double_param("beam_2CV_mm");
0206   }
0207
0208   EventHeader *eventheader = findNode::getClass<EventHeader>(topNode,
0209                                                              "EventHeader");
0210
0211   if (eventheader->get_EvtType() != DATAEVENT)
0212   {
0213     cout << __PRETTY_FUNCTION__
0214          << " disgard this event: " << endl;
0215
0216     eventheader->identify();
0217
0218     return Fun4AllReturnCodes::DISCARDEVENT;
0219   }
0220
0221   if (not _is_sim)
0222   {
0223     assert(eventheader);
0224
0225     _eval_run.run = eventheader->get_RunNumber();
0226     if (verbosity > 4)
0227       cout << __PRETTY_FUNCTION__ << _eval_run.run << endl;
0228
0229     _eval_run.event = eventheader->get_EvtSequence();
0230   }
0231
0232   if (_is_sim)
0233   {
0234     PHG4TruthInfoContainer *truthInfoList = findNode::getClass<
0235         PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
0236
0237     assert(truthInfoList);
0238
0239     _eval_run.run = -1;
0240
0241     const PHG4Particle *p =
0242         truthInfoList->GetPrimaryParticleRange().first->second;
0243     assert(p);
0244
0245     const PHG4VtxPoint *v = truthInfoList->GetVtx(p->get_vtx_id());
0246     assert(v);
0247
0248     _eval_run.beam_mom = sqrt(
0249         p->get_px() * p->get_px() + p->get_py() * p->get_py() + p->get_pz() * p->get_pz());
0250     _eval_run.truth_y = v->get_y();
0251     _eval_run.truth_z = v->get_z();
0252   }
0253
0254   // normalization
0255   TH1F *hNormalization = dynamic_cast<TH1F *>(hm->getHisto("hNormalization"));
0256   assert(hNormalization);
0257
0258   hNormalization->Fill("ALL", 1);
0259
0260   // get DST nodes
0261
0262   RawTowerContainer *TOWER_RAW_CEMC = findNode::getClass<RawTowerContainer>(
0263       topNode, _is_sim ? "TOWER_RAW_LG_CEMC" : "TOWER_RAW_CEMC");
0264   assert(TOWER_RAW_CEMC);
0265
0266   RawTowerContainer *TOWER_CALIB_CEMC = findNode::getClass<RawTowerContainer>(
0267       topNode, _is_sim ? "TOWER_CALIB_LG_CEMC" : "TOWER_CALIB_CEMC");
0268   assert(TOWER_CALIB_CEMC);
0269
0270   RawTowerContainer *TOWER_CALIB_LG_HCALIN = findNode::getClass<
0271       RawTowerContainer>(topNode, "TOWER_CALIB_LG_HCALIN");
0272   assert(TOWER_CALIB_LG_HCALIN);
0273
0274   RawTowerContainer *TOWER_CALIB_LG_HCALOUT = findNode::getClass<
0275       RawTowerContainer>(topNode, "TOWER_CALIB_LG_HCALOUT");
0276   assert(TOWER_CALIB_LG_HCALOUT);
0277
0278   // other nodes
0279   RawTowerContainer *TOWER_CALIB_TRIGGER_VETO = findNode::getClass<
0280       RawTowerContainer>(topNode, "TOWER_CALIB_TRIGGER_VETO");
0281   if (not _is_sim)
0282   {
0283     assert(TOWER_CALIB_TRIGGER_VETO);
0284   }
0285
0286   RawTowerContainer *TOWER_CALIB_HODO_HORIZONTAL = findNode::getClass<
0287       RawTowerContainer>(topNode, "TOWER_CALIB_HODO_HORIZONTAL");
0288   if (not _is_sim)
0289   {
0290     assert(TOWER_CALIB_HODO_HORIZONTAL);
0291   }
0292   RawTowerContainer *TOWER_CALIB_HODO_VERTICAL = findNode::getClass<
0293       RawTowerContainer>(topNode, "TOWER_CALIB_HODO_VERTICAL");
0294   if (not _is_sim)
0295   {
0296     assert(TOWER_CALIB_HODO_VERTICAL);
0297   }
0298
0299   RawTowerContainer *TOWER_CALIB_C1 = findNode::getClass<RawTowerContainer>(
0300       topNode, "TOWER_CALIB_C1");
0301   if (not _is_sim)
0302   {
0303     assert(TOWER_CALIB_C1);
0304   }
0305   RawTowerContainer *TOWER_CALIB_C2 = findNode::getClass<RawTowerContainer>(
0306       topNode, "TOWER_CALIB_C2");
0307   if (not _is_sim)
0308   {
0309     assert(TOWER_CALIB_C2);
0310   }
0311
0312   // Cherenkov
0313   bool cherekov_e = false;
0314   bool cherekov_anti_e = false;
0315   if (not _is_sim)
0316   {
0317     RawTower *t_c2_in = NULL;
0318     RawTower *t_c2_out = NULL;
0319
0320     t_c2_in = TOWER_CALIB_C2->getTower(0);
0321     t_c2_out = TOWER_CALIB_C2->getTower(1);
0322
0323     assert(t_c2_in);
0324     assert(t_c2_out);
0325
0326     const double c2_in = t_c2_in->get_energy();
0327     const double c2_out = t_c2_out->get_energy();
0328     const double c1 = TOWER_CALIB_C1->getTower(0)->get_energy();
0329
0330     _eval_run.C2_sum = c2_in + c2_out;
0331     _eval_run.C1 = c1;
0332
0333     cherekov_e = (_eval_run.C2_sum) > (abs(_eval_run.beam_mom) >= 10 ? 1000 : 1500);
0334     cherekov_anti_e = (_eval_run.C2_sum) < 100;
0335     hNormalization->Fill("C2-e", cherekov_e);
0336
0337     TH2F *hCheck_Cherenkov = dynamic_cast<TH2F *>(hm->getHisto(
0338         "hCheck_Cherenkov"));
0339     assert(hCheck_Cherenkov);
0340     hCheck_Cherenkov->Fill(c1, "C1", 1);
0341     hCheck_Cherenkov->Fill(c2_in, "C2 in", 1);
0342     hCheck_Cherenkov->Fill(c2_out, "C2 out", 1);
0343     hCheck_Cherenkov->Fill(c2_in + c2_out, "C2 sum", 1);
0344   }
0345
0346   // veto
0347   TH1F *hCheck_Veto = dynamic_cast<TH1F *>(hm->getHisto("hCheck_Veto"));
0348   assert(hCheck_Veto);
0349   bool trigger_veto_pass = true;
0350   if (not _is_sim)
0351   {
0352     auto range = TOWER_CALIB_TRIGGER_VETO->getTowers();
0353     for (auto it = range.first; it != range.second; ++it)
0354     {
0355       RawTower *tower = it->second;
0356       assert(tower);
0357
0358       hCheck_Veto->Fill(tower->get_energy());
0359
0360       if (abs(tower->get_energy()) > 0.2)
0361         trigger_veto_pass = false;
0362     }
0363   }
0364   hNormalization->Fill("trigger_veto_pass", trigger_veto_pass);
0365   _eval_run.trigger_veto_pass = trigger_veto_pass;
0366
0367   // hodoscope
0368   TH1F *hCheck_Hodo_H = dynamic_cast<TH1F *>(hm->getHisto("hCheck_Hodo_H"));
0369   assert(hCheck_Hodo_H);
0370   int hodo_h_count = 0;
0371   if (not _is_sim)
0372   {
0373     auto range = TOWER_CALIB_HODO_HORIZONTAL->getTowers();
0374     for (auto it = range.first; it != range.second; ++it)
0375     {
0376       RawTower *tower = it->second;
0377       assert(tower);
0378
0379       hCheck_Hodo_H->Fill(tower->get_energy());
0380
0381       if (abs(tower->get_energy()) >  0.5)
0382       {
0383         hodo_h_count++;
0384         _eval_run.hodo_h = tower->get_id();
0385       }
0386     }
0387   }
0388   const bool valid_hodo_h = hodo_h_count == 1;
0389   hNormalization->Fill("valid_hodo_h", valid_hodo_h);
0390   _eval_run.valid_hodo_h = valid_hodo_h;
0391
0392   TH1F *hCheck_Hodo_V = dynamic_cast<TH1F *>(hm->getHisto("hCheck_Hodo_V"));
0393   assert(hCheck_Hodo_V);
0394   int hodo_v_count = 0;
0395   if (not _is_sim)
0396   {
0397     auto range = TOWER_CALIB_HODO_VERTICAL->getTowers();
0398     for (auto it = range.first; it != range.second; ++it)
0399     {
0400       RawTower *tower = it->second;
0401       assert(tower);
0402
0403       hCheck_Hodo_V->Fill(tower->get_energy());
0404
0405       if (abs(tower->get_energy()) >  0.5)
0406       {
0407         hodo_v_count++;
0408         _eval_run.hodo_v = tower->get_id();
0409       }
0410     }
0411   }
0412   const bool valid_hodo_v = hodo_v_count == 1;
0413   _eval_run.valid_hodo_v = valid_hodo_v;
0414   hNormalization->Fill("valid_hodo_v", valid_hodo_v);
0415
0416   const bool good_e = (valid_hodo_v and valid_hodo_h and cherekov_e and trigger_veto_pass) and (not _is_sim);
0417   const bool good_anti_e = (valid_hodo_v and valid_hodo_h and cherekov_anti_e and trigger_veto_pass) and (not _is_sim);
0418   hNormalization->Fill("good_e", good_e);
0419   hNormalization->Fill("good_anti_e", good_anti_e);
0420
0421   // simple clustering by matching to cluster of max energy
0422   pair<int, int> max_tower = find_max(TOWER_CALIB_CEMC, 1);
0423
0424   _eval_5x5_CEMC.max_col = max_tower.first;
0425   _eval_5x5_CEMC.max_row = max_tower.second;
0426
0427   // process EMCals
0428   {
0429     double sum_energy_calib = 0;
0430
0431     auto range = TOWER_CALIB_CEMC->getTowers();
0432     for (auto it = range.first; it != range.second; ++it)
0433     {
0434       RawTower *tower = it->second;
0435       assert(tower);
0436
0437       const int col = tower->get_bineta();
0438       const int row = tower->get_binphi();
0439
0440       if (col < 0 or col >= 8)
0441         continue;
0442       if (row < 0 or row >= 8)
0443         continue;
0444
0445       const double energy_calib = tower->get_energy();
0446       sum_energy_calib += energy_calib;
0447
0448       // cluster 5x5
0449       if (col >= max_tower.first - 2 and col <= max_tower.first + 2)
0450         if (row >= max_tower.second - 2 and row <= max_tower.second + 2)
0451         {
0452           // in cluster
0453
0454           _eval_5x5_CEMC.average_col += energy_calib * col;
0455           _eval_5x5_CEMC.average_row += energy_calib * row;
0456           _eval_5x5_CEMC.sum_E += energy_calib;
0457         }
0458     }  //       for (auto it = range.first; it != range.second; ++it)
0459     _eval_5x5_CEMC.reweight_clus_pol();
0460
0461     _eval_run.sum_E_CEMC = sum_energy_calib;
0462   }  // process EMCals
0463
0464   // process inner HCAL
0465   {
0466     double sum_energy_calib = 0;
0467
0468     auto range = TOWER_CALIB_LG_HCALIN->getTowers();
0469     for (auto it = range.first; it != range.second; ++it)
0470     {
0471       RawTower *tower = it->second;
0472       assert(tower);
0473
0474       const int col = tower->get_bineta();
0475       const int row = tower->get_binphi();
0476
0477       if (col < 0 or col >= 4)
0478         continue;
0479       if (row < 0 or row >= 4)
0480         continue;
0481
0482       const double energy_calib = tower->get_energy();
0483       sum_energy_calib += energy_calib;
0484
0485     }  //       for (auto it = range.first; it != range.second; ++it)
0486     _eval_run.sum_E_HCAL_IN = sum_energy_calib;
0487   }  // process inner HCAL
0488
0489   // process outer HCAL
0490   {
0491     double sum_energy_calib = 0;
0492
0493     auto range = TOWER_CALIB_LG_HCALOUT->getTowers();
0494     for (auto it = range.first; it != range.second; ++it)
0495     {
0496       RawTower *tower = it->second;
0497       assert(tower);
0498
0499       const int col = tower->get_bineta();
0500       const int row = tower->get_binphi();
0501
0502       if (col < 0 or col >= 4)
0503         continue;
0504       if (row < 0 or row >= 4)
0505         continue;
0506
0507       const double energy_calib = tower->get_energy();
0508       sum_energy_calib += energy_calib;
0509
0510     }  //       for (auto it = range.first; it != range.second; ++it)
0511     _eval_run.sum_E_HCAL_OUT = sum_energy_calib;
0512   }  // process outer HCAL
0513
0514   _eval_run.good_e = good_e;
0515   _eval_run.good_anti_e = good_anti_e;
0516
0517   TTree *T = dynamic_cast<TTree *>(hm->getHisto("T"));
0518   assert(T);
0519   T->Fill();
0520
0521   return Fun4AllReturnCodes::EVENT_OK;
0522 }
0523
0524 pair<int, int>
0525 ExampleAnalysisModule::find_max(RawTowerContainer *towers, int cluster_size)
0526 {
0527   const int clus_edge_size = (cluster_size - 1) / 2;
0528   assert(clus_edge_size >= 0);
0529
0530   pair<int, int> max(-1000, -1000);
0531   double max_e = 0;
0532
0533   for (int col = 0; col < n_size; ++col)
0534     for (int row = 0; row < n_size; ++row)
0535     {
0536       double energy = 0;
0537
0538       for (int dcol = col - clus_edge_size; dcol <= col + clus_edge_size;
0539            ++dcol)
0540         for (int drow = row - clus_edge_size; drow <= row + clus_edge_size;
0541              ++drow)
0542         {
0543           if (dcol < 0 or drow < 0)
0544             continue;
0545
0546           RawTower *t = towers->getTower(dcol, drow);
0547           if (t)
0548             energy += t->get_energy();
0549         }
0550
0551       if (energy > max_e)
0552       {
0553         max = make_pair(col, row);
0554         max_e = energy;
0555       }
0556     }
0557
0558   return max;
0559 }