tutorial/ExampleAnalysisModule/ExampleAnalysisModule.C

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