calotriggeremulator/src/DijetEventDisplay.cc

0001 #include "DijetEventDisplay.h"
0002 #include <phool/PHNode.h>
0003 #include <phool/getClass.h>
0004 #include <ffaobjects/EventHeader.h>
0005 #include <ffaobjects/EventHeaderv1.h>
0006 #include <ffaobjects/EventHeaderv2.h>
0007 DijetEventDisplay::DijetEventDisplay(const std::string f)
0008   : file{f}
0009 {
0010   hm = new Fun4AllHistoManager("DIJETEVENTDISPLAY");
0011 }
0012 int DijetEventDisplay::FillEvent(PHCompositeNode *topNode, float aj, float dphi)
0013 {
0014   EventHeader *eventheader = findNode::getClass<EventHeader>(topNode, "EventHeader");
0015   RawTowerGeomContainer *tower_geomIH = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALIN");
0016   RawTowerGeomContainer *tower_geomOH = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALOUT");
0017   TowerInfoContainer *hcalin_towers = findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_CALIB_HCALIN");
0018   TowerInfoContainer *hcalout_towers = findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_CALIB_HCALOUT");
0019   TowerInfoContainer *emcalre_towers = findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_CALIB_CEMC_RETOWER");
0020
0021   JetContainer *jets_4 = findNode::getClass<JetContainer>(topNode, "AntiKt_Tower_r04_Sub1");
0022   int eventnumber = eventheader->get_EvtSequence();
0023
0024   std::string aj_name = "h_aj_" + std::to_string(eventnumber);
0025   std::string dphi_name = "h_dphi_" + std::to_string(eventnumber);
0026   TH1F *haj = new TH1F(aj_name.c_str(), "", 1, -1, 1);
0027   haj->SetBinContent(1, aj);
0028   hm->registerHisto(haj);
0029   TH1F *hdphi = new TH1F(dphi_name.c_str(), "", 1, -1, 1);
0030   hdphi->SetBinContent(1, dphi);
0031   hm->registerHisto(hdphi);
0032   float binedges_eta_in[25];
0033   float binedges_eta_out[25];
0034   float binedges_phi_in[65];
0035   float binedges_phi_out[65];
0036   for (int ieta = 0; ieta < 24; ieta++)
0037     {
0038       const RawTowerDefs::keytype key = RawTowerDefs::encode_towerid(RawTowerDefs::CalorimeterId::HCALIN, ieta, 0);
0039       float tower_eta_in = tower_geomIH->get_tower_geometry(key)->get_eta();
0040       const RawTowerDefs::keytype keyo = RawTowerDefs::encode_towerid(RawTowerDefs::CalorimeterId::HCALOUT, ieta, 0);
0041       float tower_eta_out = tower_geomOH->get_tower_geometry(keyo)->get_eta();
0042       binedges_eta_in[ieta] = tower_eta_in;
0043       binedges_eta_out[ieta] = tower_eta_out;
0044     }
0045
0046   binedges_eta_in[24] = binedges_eta_in[23] + (binedges_eta_in[23] - binedges_eta_in[22]);
0047   binedges_eta_out[24] = binedges_eta_out[23] + (binedges_eta_out[23] - binedges_eta_out[22]);
0048   for (int ieta = 0; ieta < 25; ieta++)
0049     {
0050       binedges_eta_in[ieta] = binedges_eta_in[ieta] - (binedges_eta_in[1] - binedges_eta_in[0])/2.;
0051       binedges_eta_out[ieta] = binedges_eta_out[ieta] - (binedges_eta_out[1] - binedges_eta_out[0])/2.;
0052     }
0053   std::vector<float> binedge_phi_in;
0054   std::vector<float> binedge_phi_out;
0055   for (int iphi = 0; iphi < 64; iphi++)
0056     {
0057       const RawTowerDefs::keytype key = RawTowerDefs::encode_towerid(RawTowerDefs::CalorimeterId::HCALIN, 0, iphi);
0058       float tower_phi_in = tower_geomIH->get_tower_geometry(key)->get_phi();
0059       const RawTowerDefs::keytype keyo = RawTowerDefs::encode_towerid(RawTowerDefs::CalorimeterId::HCALOUT, 0, iphi);
0060       float tower_phi_out = tower_geomOH->get_tower_geometry(keyo)->get_phi();
0061
0062       binedge_phi_in.push_back(tower_phi_in);
0063       binedge_phi_out.push_back(tower_phi_out);
0064
0065     }
0066   std::sort(binedge_phi_in.begin(), binedge_phi_in.end());
0067   std::sort(binedge_phi_out.begin(), binedge_phi_out.end());
0068   for (int iphi = 0; iphi < 64; iphi++)
0069     {
0070       binedges_phi_in[iphi] = binedge_phi_in.at(iphi);
0071       binedges_phi_out[iphi] = binedge_phi_out.at(iphi);
0072     }
0073   binedges_phi_in[64] = binedges_phi_in[63] + (binedges_phi_in[63] - binedges_phi_in[62]);
0074   binedges_phi_out[64] = binedges_phi_out[63] + (binedges_phi_out[63] - binedges_phi_out[62]);
0075   for (int iphi = 0; iphi < 65; iphi++)
0076     {
0077
0078       binedges_phi_in[iphi] = binedges_phi_in[iphi] - (binedges_phi_in[1] - binedges_phi_in[0])/2.;
0079       binedges_phi_out[iphi] = binedges_phi_out[iphi] - (binedges_phi_out[1] - binedges_phi_out[0])/2.;
0080     }
0081   int ijet = 0;
0082
0083   for (auto jet : *jets_4)
0084     {
0085       if (jet->get_pt() < 1) continue;
0086
0087       std::string emcal_name = "h_emcal_dijet_" + std::to_string(eventnumber) + "_jet" + std::to_string(ijet);
0088       std::string hcalin_name = "h_hcalin_dijet_" + std::to_string(eventnumber) + "_jet" + std::to_string(ijet);
0089       std::string hcalout_name = "h_hcalout_dijet_" + std::to_string(eventnumber) + "_jet" + std::to_string(ijet);
0090       std::string all_name = "h_all_dijet_" + std::to_string(eventnumber) + "_jet" + std::to_string(ijet);
0091
0092       ijet++;
0093
0094       h_emcal = new TH2D(emcal_name.c_str(), ";#eta;#phi", 24, binedges_eta_in, 64, binedges_phi_in);
0095       h_hcalin = new TH2D(hcalin_name.c_str(), ";#eta;#phi", 24, binedges_eta_in, 64, binedges_phi_in);
0096       h_hcalout = new TH2D(hcalout_name.c_str(), ";#eta;#phi", 24, binedges_eta_out, 64, binedges_phi_out);
0097       h_all = new TH2D(all_name.c_str(), ";#eta;#phi", 24, binedges_eta_out, 64, binedges_phi_out);
0098
0099       hm->registerHisto(h_all);
0100       hm->registerHisto(h_emcal);
0101       hm->registerHisto(h_hcalin);
0102       hm->registerHisto(h_hcalout);
0103
0104       for (auto comp : jet->get_comp_vec())
0105     {
0106       unsigned int channel = comp.second;
0107       TowerInfo *tower;
0108       float tower_eT = 0;
0109       if (comp.first == 26 || comp.first == 30)
0110         {  // IHcal
0111           tower = hcalin_towers->get_tower_at_channel(channel);
0112           if (!tower || !tower_geomIH)
0113         {
0114           continue;
0115         }
0116           if(!tower->get_isGood()) continue;
0117
0118           unsigned int calokey = hcalin_towers->encode_key(channel);
0119
0120           int ieta = hcalin_towers->getTowerEtaBin(calokey);
0121           int iphi = hcalin_towers->getTowerPhiBin(calokey);
0122           const RawTowerDefs::keytype key = RawTowerDefs::encode_towerid(RawTowerDefs::CalorimeterId::HCALIN, ieta, iphi);
0123           float tower_phi = tower_geomIH->get_tower_geometry(key)->get_phi();
0124           float tower_eta = tower_geomIH->get_tower_geometry(key)->get_eta();
0125           tower_eT = tower->get_energy() / std::cosh(tower_eta);
0126
0127           h_hcalin->Fill(tower_eta, tower_phi, tower_eT);
0128           h_all->Fill(tower_eta, tower_phi, tower_eT);
0129         }
0130       else if (comp.first == 27 || comp.first == 31)
0131         {  // IHcal
0132           tower = hcalout_towers->get_tower_at_channel(channel);
0133
0134           if (!tower || !tower_geomOH)
0135         {
0136           continue;
0137         }
0138
0139           unsigned int calokey = hcalout_towers->encode_key(channel);
0140           int ieta = hcalout_towers->getTowerEtaBin(calokey);
0141           int iphi = hcalout_towers->getTowerPhiBin(calokey);
0142           const RawTowerDefs::keytype key = RawTowerDefs::encode_towerid(RawTowerDefs::CalorimeterId::HCALOUT, ieta, iphi);
0143           float tower_phi = tower_geomOH->get_tower_geometry(key)->get_phi();
0144           float tower_eta = tower_geomOH->get_tower_geometry(key)->get_eta();
0145           tower_eT = tower->get_energy() / std::cosh(tower_eta);
0146
0147           h_hcalout->Fill(tower_eta, tower_phi, tower_eT);
0148           h_all->Fill(tower_eta, tower_phi, tower_eT);
0149         }
0150       else if (comp.first == 28 || comp.first == 29)
0151         {  // IHcal
0152           tower = emcalre_towers->get_tower_at_channel(channel);
0153
0154           if (!tower || !tower_geomIH)
0155         {
0156           continue;
0157         }
0158
0159           unsigned int calokey = emcalre_towers->encode_key(channel);
0160           int ieta = emcalre_towers->getTowerEtaBin(calokey);
0161           int iphi = emcalre_towers->getTowerPhiBin(calokey);
0162           const RawTowerDefs::keytype key = RawTowerDefs::encode_towerid(RawTowerDefs::CalorimeterId::HCALIN, ieta, iphi);
0163           float tower_phi = tower_geomIH->get_tower_geometry(key)->get_phi();
0164           float tower_eta = tower_geomIH->get_tower_geometry(key)->get_eta();
0165           tower_eT = tower->get_energy() / std::cosh(tower_eta);
0166
0167           h_emcal->Fill(tower_eta, tower_phi, tower_eT);
0168           h_all->Fill(tower_eta, tower_phi, tower_eT);
0169         }
0170     }
0171     }
0172   std::string emcal_namefull = "h_emcal_full_" + std::to_string(eventnumber);
0173   std::string hcalin_namefull = "h_hcalin_full_" + std::to_string(eventnumber);
0174   std::string hcalout_namefull = "h_hcalout_full_" + std::to_string(eventnumber);
0175   std::string all_namefull = "h_all_full_" + std::to_string(eventnumber);
0176
0177   ijet++;
0178
0179   h_emcal = new TH2D(emcal_namefull.c_str(), ";#eta;#phi", 24, binedges_eta_in, 64, binedges_phi_in);
0180   h_hcalin = new TH2D(hcalin_namefull.c_str(), ";#eta;#phi", 24, binedges_eta_in, 64, binedges_phi_in);
0181   h_hcalout = new TH2D(hcalout_namefull.c_str(), ";#eta;#phi", 24, binedges_eta_out, 64, binedges_phi_out);
0182   h_all = new TH2D(all_namefull.c_str(), ";#eta;#phi", 24, binedges_eta_out, 64, binedges_phi_out);
0183
0184   hm->registerHisto(h_all);
0185   hm->registerHisto(h_emcal);
0186   hm->registerHisto(h_hcalin);
0187   hm->registerHisto(h_hcalout);
0188
0189   int size;
0190   size = emcalre_towers->size();
0191   for (int channel = 0; channel < size; channel++)
0192     {
0193       TowerInfo *tower = emcalre_towers->get_tower_at_channel(channel);
0194
0195       unsigned int calokey = emcalre_towers->encode_key(channel);
0196       int ieta = emcalre_towers->getTowerEtaBin(calokey);
0197       int iphi = emcalre_towers->getTowerPhiBin(calokey);
0198       const RawTowerDefs::keytype key = RawTowerDefs::encode_towerid(RawTowerDefs::CalorimeterId::HCALIN, ieta, iphi);
0199       float tower_phi = tower_geomIH->get_tower_geometry(key)->get_phi();
0200       float tower_eta = tower_geomIH->get_tower_geometry(key)->get_eta();
0201       float tower_eT = tower->get_energy() / std::cosh(tower_eta);
0202
0203       h_emcal->Fill(tower_eta, tower_phi, tower_eT);
0204       h_all->Fill(tower_eta, tower_phi, tower_eT);
0205     }
0206   size = hcalin_towers->size();
0207   for (int channel = 0; channel < size; channel++)
0208     {
0209       TowerInfo *tower = hcalin_towers->get_tower_at_channel(channel);
0210
0211       unsigned int calokey = hcalin_towers->encode_key(channel);
0212       int ieta = hcalin_towers->getTowerEtaBin(calokey);
0213       int iphi = hcalin_towers->getTowerPhiBin(calokey);
0214       const RawTowerDefs::keytype key = RawTowerDefs::encode_towerid(RawTowerDefs::CalorimeterId::HCALIN, ieta, iphi);
0215       float tower_phi = tower_geomIH->get_tower_geometry(key)->get_phi();
0216       float tower_eta = tower_geomIH->get_tower_geometry(key)->get_eta();
0217       float tower_eT = tower->get_energy() / std::cosh(tower_eta);
0218
0219       h_hcalin->Fill(tower_eta, tower_phi, tower_eT);
0220       h_all->Fill(tower_eta, tower_phi, tower_eT);
0221     }
0222   size = hcalout_towers->size();
0223   for (int channel = 0; channel < size; channel++)
0224     {
0225       TowerInfo *tower = hcalout_towers->get_tower_at_channel(channel);
0226
0227       unsigned int calokey = hcalout_towers->encode_key(channel);
0228       int ieta = hcalout_towers->getTowerEtaBin(calokey);
0229       int iphi = hcalout_towers->getTowerPhiBin(calokey);
0230       const RawTowerDefs::keytype key = RawTowerDefs::encode_towerid(RawTowerDefs::CalorimeterId::HCALOUT, ieta, iphi);
0231       float tower_phi = tower_geomOH->get_tower_geometry(key)->get_phi();
0232       float tower_eta = tower_geomOH->get_tower_geometry(key)->get_eta();
0233       float tower_eT = tower->get_energy() / std::cosh(tower_eta);
0234
0235       h_hcalout->Fill(tower_eta, tower_phi, tower_eT);
0236       h_all->Fill(tower_eta, tower_phi, tower_eT);
0237     }
0238   return 0;
0239 }
0240
0241 void DijetEventDisplay::Dump()
0242 {
0243   hm->dumpHistos(file);
0244 }