File indexing completed on 2025-08-05 08:11:12
0001 #include "DijetHistoMaker.h"
0002
0003 #include <TH1.h>
0004 #include <TH2.h>
0005 #include <TH2D.h>
0006 #include <TFile.h>
0007 #include <TF1.h>
0008 #include <TH1D.h>
0009 #include <g4main/PHG4TruthInfoContainer.h>
0010 #include <g4main/PHG4Particle.h>
0011 #include <TMath.h>
0012 #include <calotrigger/TriggerRunInfov1.h>
0013
0014 #include <calobase/RawTowerGeom.h>
0015 #include <jetbackground/TowerBackground.h>
0016 #include <jetbase/Jet.h>
0017 #include <jetbase/JetContainer.h>
0018 #include <jetbase/JetContainerv1.h>
0019 #include <jetbase/Jetv2.h>
0020 #include <calobase/TowerInfoContainer.h>
0021 #include <calobase/TowerInfoContainerv1.h>
0022 #include <calobase/TowerInfoContainerv2.h>
0023 #include <calobase/TowerInfoContainerv3.h>
0024 #include <calobase/TowerInfo.h>
0025 #include <calobase/TowerInfov1.h>
0026 #include <calobase/TowerInfov2.h>
0027 #include <calobase/TowerInfov3.h>
0028
0029 #include <calobase/RawCluster.h>
0030 #include <calobase/RawClusterContainer.h>
0031 #include <calobase/RawClusterUtility.h>
0032 #include <calobase/RawTowerGeomContainer.h>
0033 #include <calobase/RawTower.h>
0034 #include <calobase/RawTowerContainer.h>
0035 #include <fun4all/Fun4AllHistoManager.h>
0036 #include <HepMC/SimpleVector.h>
0037
0038 #include <globalvertex/GlobalVertex.h>
0039 #include <globalvertex/GlobalVertexMap.h>
0040 #include "TrashInfov1.h"
0041 #include "DijetEventDisplay.h"
0042 #include <vector>
0043
0044 #include <fun4all/Fun4AllReturnCodes.h>
0045 #include <phool/PHCompositeNode.h>
0046 #include <phool/PHIODataNode.h>
0047 #include <phool/PHNode.h>
0048 #include <phool/PHNodeIterator.h>
0049 #include <phool/PHObject.h>
0050 #include <phool/getClass.h>
0051
0052
0053 #include <iostream>
0054
0055 #include <map>
0056
0057
0058 DijetHistoMaker::DijetHistoMaker(const std::string &name, const std::string &outfilename):
0059 SubsysReco(name)
0060
0061 {
0062 isSim = false;
0063 _foutname = outfilename;
0064 m_calo_nodename = "TOWERINFO_CALIB";
0065 m_jet_triggernames = {"Jet 6 GeV + MBD NS >= 1",
0066 "Jet 8 GeV + MBD NS >= 1",
0067 "Jet 10 GeV + MBD NS >= 1",
0068 "Jet 12 GeV + MBD NS >= 1",
0069 "Jet 6 GeV",
0070 "Jet 8 GeV",
0071 "Jet 10 GeV",
0072 "Jet 12 GeV"};
0073 }
0074
0075
0076 DijetHistoMaker::~DijetHistoMaker()
0077 {
0078
0079
0080
0081
0082
0083
0084
0085
0086 }
0087
0088
0089 int DijetHistoMaker::Init(PHCompositeNode *topNode)
0090 {
0091
0092
0093
0094
0095
0096
0097
0098 std::string outfile_display = _foutname;
0099 auto it = outfile_display.find("DIJET");
0100 outfile_display.replace(it, 5, "DIJETEVT");
0101 if (!isSim)
0102 {
0103 if (drawDijets) dijeteventdisplay = new DijetEventDisplay(outfile_display.c_str());
0104 triggeranalyzer = new TriggerAnalyzer();
0105 }
0106 hm = new Fun4AllHistoManager("DIJET_HISTOGRAMS");
0107
0108 TH1D *h;
0109 TH2D *h2;
0110
0111 h2 = new TH2D("h_pt1_pt2",";p_{T1};p_{T2}", 50, 0, 50, 50, 0, 50);
0112 hm->registerHisto(h2);
0113 h2 = new TH2D("h_et1_et2",";e_{T1};e_{T2}", 50, 0, 50, 50, 0, 50);
0114 hm->registerHisto(h2);
0115
0116 if (!isSim)
0117 {
0118 h2 = new TH2D("h_pt1_pt2_mbd",";p_{T1};p_{T2}", 50, 0, 50, 50, 0, 50);
0119 hm->registerHisto(h2);
0120 h2 = new TH2D("h_et1_et2_mbd",";e_{T1};e_{T2}", 50, 0, 50, 50, 0, 50);
0121 hm->registerHisto(h2);
0122
0123 for (int j = 0; j < 8; j++)
0124 {
0125 h2 = new TH2D(Form("h_pt1_pt2_gl1_%d", j),";p_{T1};p_{T2}", 50, 0, 50, 50, 0, 50);
0126 hm->registerHisto(h2);
0127 h2 = new TH2D(Form("h_et1_et2_gl1_%d", j),";e_{T1};e_{T2}", 50, 0, 50, 50, 0, 50);
0128 hm->registerHisto(h2);
0129 }
0130 }
0131 else
0132 {
0133 h2 = new TH2D("h_truth_pt1_pt2",";p_{T1};p_{T2}", 50, 0, 50, 50, 0, 50);
0134 hm->registerHisto(h2);
0135 h2 = new TH2D("h_truth_smear_pt1_pt2",";p_{T1};p_{T2}", 50, 0, 50, 50, 0, 50);
0136 hm->registerHisto(h2);
0137
0138 }
0139 for (int i = 0; i < 3; i++)
0140 {
0141
0142 h = new TH1D(Form("h_dphi_et_%d", i), ";#Delta #phi;N", 64, 0, TMath::Pi());
0143 hm->registerHisto(h);
0144 h = new TH1D(Form("h_dphi_%d", i), ";#Delta #phi;N", 64, 0, TMath::Pi());
0145 hm->registerHisto(h);
0146 if (!isSim)
0147 {
0148 h = new TH1D(Form("h_dphi_mb_%d", i), ";#Delta #phi;N", 64, 0, TMath::Pi());
0149 hm->registerHisto(h);
0150 h = new TH1D(Form("h_dphi_et_mb_%d", i), ";#Delta #phi;N", 64, 0, TMath::Pi());
0151 hm->registerHisto(h);
0152
0153 }
0154 else
0155 {
0156 h = new TH1D(Form("h_truth_dphi_%d", i), ";#Delta #phi;N", 64, 0, TMath::Pi());
0157 hm->registerHisto(h);
0158 h = new TH1D(Form("h_truth_smear_dphi_%d", i), ";#Delta #phi;N", 64, 0, TMath::Pi());
0159 hm->registerHisto(h);
0160
0161 }
0162 for (int ip = 0; ip < 3; ip++)
0163 {
0164
0165 h = new TH1D(Form("h_Aj_et_%d_dp%d", i, ip), ";A_j;N", 20, 0, 1);
0166 hm->registerHisto(h);
0167 h = new TH1D(Form("h_xj_et_%d_dp%d", i, ip), ";x_j;N", 20, 0, 1);
0168 hm->registerHisto(h);
0169
0170 h = new TH1D(Form("h_Aj_%d_dp%d", i, ip), ";A_j;N", 20, 0, 1);
0171 hm->registerHisto(h);
0172 h = new TH1D(Form("h_xj_%d_dp%d", i, ip), ";x_j;N", 20, 0, 1);
0173 hm->registerHisto(h);
0174 for (int j = 0; j < 2; j++)
0175 {
0176 h = new TH1D(Form("h_Aj_et_iso%d_%d_dp%d", j, i, ip), ";A_j;N", 20, 0, 1);
0177 hm->registerHisto(h);
0178 h = new TH1D(Form("h_xj_et_iso%d_%d_dp%d", j, i, ip), ";x_j;N", 20, 0, 1);
0179 hm->registerHisto(h);
0180
0181 h = new TH1D(Form("h_Aj_iso%d_%d_dp%d", j, i, ip), ";A_j;N", 20, 0, 1);
0182 hm->registerHisto(h);
0183 h = new TH1D(Form("h_xj_iso%d_%d_dp%d", j, i, ip), ";x_j;N", 20, 0, 1);
0184 hm->registerHisto(h);
0185 }
0186
0187 if (!isSim)
0188 {
0189 h = new TH1D(Form("h_Aj_et_mb_%d_dp%d", i, ip), ";A_j;N", 20, 0, 1);
0190 hm->registerHisto(h);
0191 h = new TH1D(Form("h_xj_et_mb_%d_dp%d", i, ip), ";x_j;N", 20, 0, 1);
0192 hm->registerHisto(h);
0193 h = new TH1D(Form("h_Aj_mb_%d_dp%d", i, ip), ";A_j;N", 20, 0, 1);
0194 hm->registerHisto(h);
0195 h = new TH1D(Form("h_xj_mb_%d_dp%d", i, ip), ";x_j;N", 20, 0, 1);
0196 hm->registerHisto(h);
0197
0198
0199 for (int j = 0; j < 8; j++)
0200 {
0201 h = new TH1D(Form("h_Aj_et_gl1_%d_%d_dp%d", j, i, ip), ";A_j;N", 20, 0, 1);
0202 hm->registerHisto(h);
0203 h = new TH1D(Form("h_xj_et_gl1_%d_%d_dp%d", j, i, ip), ";x_j;N", 20, 0, 1);
0204 hm->registerHisto(h);
0205 h = new TH1D(Form("h_Aj_gl1_%d_%d_dp%d", j, i, ip), ";A_j;N", 20, 0, 1);
0206 hm->registerHisto(h);
0207 h = new TH1D(Form("h_xj_gl1_%d_%d_dp%d", j, i, ip), ";x_j;N", 20, 0, 1);
0208 hm->registerHisto(h);
0209 if (ip == 0)
0210 {
0211 h = new TH1D(Form("h_dphi_gl1_%d_%d", j, i), ";#Delta #phi;N", 64, 0, TMath::Pi());
0212 hm->registerHisto(h);
0213 h = new TH1D(Form("h_dphi_et_gl1_%d_%d", j, i), ";#Delta #phi;N", 64, 0, TMath::Pi());
0214 hm->registerHisto(h);
0215 }
0216 }
0217 }
0218 else
0219 {
0220
0221 h = new TH1D(Form("h_truth_Aj_%d_dp%d", i, ip), ";A_j;N", 20, 0, 1);
0222 hm->registerHisto(h);
0223 h = new TH1D(Form("h_truth_xj_%d_dp%d", i, ip), ";x_j;N", 20, 0, 1);
0224 hm->registerHisto(h);
0225 h = new TH1D(Form("h_truth_smear_Aj_%d_dp%d", i, ip), ";A_j;N", 20, 0, 1);
0226 hm->registerHisto(h);
0227 h = new TH1D(Form("h_truth_smear_xj_%d_dp%d", i, ip), ";x_j;N", 20, 0, 1);
0228 hm->registerHisto(h);
0229
0230 }
0231 }
0232 }
0233
0234 if (!isSim)
0235 {
0236 h2 = new TH2D("h_trash_pt1_pt2",";p_{T1};p_{T2}", 50, 0, 50, 50, 0, 50);
0237 hm->registerHisto(h2);
0238 h2 = new TH2D("h_trash_et1_et2",";e_{T1};e_{T2}", 50, 0, 50, 50, 0, 50);
0239 hm->registerHisto(h2);
0240
0241 h2 = new TH2D("h_trash_pt1_pt2_mb",";p_{T1};p_{T2}", 50, 0, 50, 50, 0, 50);
0242 hm->registerHisto(h2);
0243 h2 = new TH2D("h_trash_et1_et2_mb",";e_{T1};e_{T2}", 50, 0, 50, 50, 0, 50);
0244 hm->registerHisto(h2);
0245
0246 for (int j = 0; j < 8; j++)
0247 {
0248 h2 = new TH2D(Form("h_trash_pt1_pt2_gl1_%d", j),";p_{T1};p_{T2}", 50, 0, 50, 50, 0, 50);
0249 hm->registerHisto(h2);
0250 h2 = new TH2D(Form("h_trash_et1_et2_gl1_%d",j),";e_{T1};e_{T2}", 50, 0, 50, 50, 0, 50);
0251 hm->registerHisto(h2);
0252
0253 }
0254 for (int i = 0; i < 3; i++)
0255 {
0256
0257 h = new TH1D(Form("h_trash_dphi_et_%d", i), ";#Delta #phi;N", 64, 0, TMath::Pi());
0258 hm->registerHisto(h);
0259 h = new TH1D(Form("h_trash_dphi_%d", i), ";#Delta #phi;N", 64, 0, TMath::Pi());
0260 hm->registerHisto(h);
0261
0262 h = new TH1D(Form("h_trash_dphi_mb_%d", i), ";#Delta #phi;N", 64, 0, TMath::Pi());
0263 hm->registerHisto(h);
0264 h = new TH1D(Form("h_trash_dphi_et_mb_%d", i), ";#Delta #phi;N", 64, 0, TMath::Pi());
0265 hm->registerHisto(h);
0266
0267 for (int ip = 0; ip < 3; ip++)
0268 {
0269
0270 h = new TH1D(Form("h_trash_Aj_et_%d_dp%d", i, ip), ";A_j;N", 20, 0, 1);
0271 hm->registerHisto(h);
0272 h = new TH1D(Form("h_trash_xj_et_%d_dp%d", i, ip), ";x_j;N", 20, 0, 1);
0273 hm->registerHisto(h);
0274 h = new TH1D(Form("h_trash_Aj_%d_dp%d", i, ip), ";A_j;N", 20, 0, 1);
0275 hm->registerHisto(h);
0276 h = new TH1D(Form("h_trash_xj_%d_dp%d", i, ip), ";x_j;N", 20, 0, 1);
0277 hm->registerHisto(h);
0278
0279 h = new TH1D(Form("h_trash_Aj_et_mb_%d_dp%d", i, ip), ";A_j;N", 20, 0, 1);
0280 hm->registerHisto(h);
0281 h = new TH1D(Form("h_trash_xj_et_mb_%d_dp%d", i, ip), ";x_j;N", 20, 0, 1);
0282 hm->registerHisto(h);
0283 h = new TH1D(Form("h_trash_Aj_mb_%d_dp%d", i, ip), ";A_j;N", 20, 0, 1);
0284 hm->registerHisto(h);
0285 h = new TH1D(Form("h_trash_xj_mb_%d_dp%d", i, ip), ";x_j;N", 20, 0, 1);
0286 hm->registerHisto(h);
0287
0288 for (int j = 0; j < 8; j++)
0289 {
0290
0291 h = new TH1D(Form("h_trash_Aj_et_gl1_%d_%d_dp%d", j, i, ip), ";A_j;N", 20, 0, 1);
0292 hm->registerHisto(h);
0293 h = new TH1D(Form("h_trash_xj_et_gl1_%d_%d_dp%d", j, i, ip), ";x_j;N", 20, 0, 1);
0294 hm->registerHisto(h);
0295 h = new TH1D(Form("h_trash_Aj_gl1_%d_%d_dp%d", j, i, ip), ";A_j;N", 20, 0, 1);
0296 hm->registerHisto(h);
0297 h = new TH1D(Form("h_trash_xj_gl1_%d_%d_dp%d", j, i, ip), ";x_j;N", 20, 0, 1);
0298 hm->registerHisto(h);
0299 if (ip == 0)
0300 {
0301 h = new TH1D(Form("h_trash_dphi_gl1_%d_%d", j, i), ";#Delta #phi;N", 64, 0, TMath::Pi());
0302 hm->registerHisto(h);
0303 h = new TH1D(Form("h_trash_dphi_et_gl1_%d_%d", j, i), ";#Delta #phi;N", 64, 0, TMath::Pi());
0304 hm->registerHisto(h);
0305 }
0306 }
0307 }
0308 }
0309 }
0310
0311 return Fun4AllReturnCodes::EVENT_OK;
0312 }
0313
0314
0315 int DijetHistoMaker::InitRun(PHCompositeNode *topNode)
0316 {
0317 return Fun4AllReturnCodes::EVENT_OK;
0318 }
0319
0320
0321 void DijetHistoMaker::process_truth_jets(JetContainer *jets)
0322 {
0323
0324 fsmear = new TF1("fg","gaus", 0, 2);
0325 float jes = 1;
0326 float jer = 0.18;
0327 fsmear->SetParameters(1, jes, jer);
0328 std::vector<jetty> truth_jets{};
0329 std::vector<jetty> truth_smear_jets{};
0330 if (Verbosity()) jets->identify(std::cout);
0331 for (auto jet : *jets)
0332 {
0333 if (Verbosity()) std::cout << __LINE__ << std::endl;
0334 struct jetty myjet;
0335 myjet.pt = jet->get_pt();
0336 myjet.eta = jet->get_eta();
0337 myjet.phi = jet->get_phi();
0338 truth_jets.push_back(myjet);
0339
0340 float pt = jet->get_pt() * fsmear->GetRandom();
0341
0342 struct jetty myjetsmear;
0343 myjetsmear.pt = pt;
0344 myjetsmear.eta = jet->get_eta();
0345 myjetsmear.phi = jet->get_phi();
0346 truth_smear_jets.push_back(myjetsmear);
0347
0348 }
0349 if ((int) truth_jets.size() >= 2)
0350 {
0351 std::sort(truth_jets.begin(), truth_jets.end(), [] (auto a, auto b){ return a.pt > b.pt; });
0352 float eta1 = truth_jets.at(0).eta;
0353 float phi1 = truth_jets.at(0).phi;
0354 float eta2 = truth_jets.at(1).eta;
0355 float phi2 = truth_jets.at(1).phi;
0356 if (fabs(eta1) < 0.7 && fabs(eta2) < 0.7)
0357 {
0358
0359 float dphi = phi1 - phi2;
0360
0361 if (dphi > TMath::Pi())
0362 {
0363 dphi = dphi - 2.*TMath::Pi();
0364 }
0365 if (dphi <= -1*TMath::Pi())
0366 {
0367 dphi = dphi + 2. * TMath::Pi();
0368 }
0369
0370 dphi = fabs(dphi);
0371
0372 if (dphi > TMath::Pi())
0373 return;
0374
0375 if (Verbosity())
0376 {
0377 std::cout << "Truth Dijet" << std::endl;
0378 std::cout << "Leading Jet" << std::endl;
0379 int jl = 0;
0380 std::cout << " pt = " << truth_jets.at(jl).pt << std::endl;
0381 std::cout << " eta = " << truth_jets.at(jl).eta << std::endl;
0382 std::cout << " phi = " << truth_jets.at(jl).phi << std::endl;
0383 std::cout << "Subleading Jet" << std::endl;
0384 jl++;
0385 std::cout << " pt = " << truth_jets.at(jl).pt << std::endl;
0386 std::cout << " eta = " << truth_jets.at(jl).eta << std::endl;
0387 std::cout << " phi = " << truth_jets.at(jl).phi << std::endl;
0388 std::cout << " " << std::endl;
0389 }
0390
0391
0392 if (truth_jets.at(0).pt > pt_cut1[0] && truth_jets.at(1).pt > pt_cut2[0])
0393 {
0394 auto h_pt1_pt2 = dynamic_cast<TH2*>(hm->getHisto("h_truth_pt1_pt2"));
0395 h_pt1_pt2->Fill(truth_jets.at(0).pt, truth_jets.at(1).pt);
0396 h_pt1_pt2->Fill(truth_jets.at(1).pt, truth_jets.at(0).pt);
0397 }
0398
0399 float Aj_pt = (truth_jets.at(0).pt - truth_jets.at(1).pt)/(truth_jets.at(0).pt + truth_jets.at(1).pt);
0400 float xj_pt = (truth_jets.at(1).pt / truth_jets.at(0).pt);
0401
0402 if (Verbosity())
0403 {
0404 std::cout << "---- DIJET ---- " << std::endl;
0405 std::cout << "Calculated: " << std::endl;
0406 std::cout << " A_j = " << Aj_pt << std::endl;
0407 std::cout << " x_j = " << xj_pt << std::endl;
0408 std::cout << " dPhi = " << dphi << std::endl;
0409 std::cout << " ------------------------ " << std::endl;
0410 }
0411
0412 for (int i = 0; i < 3; i++)
0413 {
0414 if (truth_jets.at(0).pt > pt_cut1[i] && truth_jets.at(1).pt > pt_cut2[i])
0415 {
0416 auto h_dphi = dynamic_cast<TH1*>(hm->getHisto(Form("h_truth_dphi_%d", i)));
0417 h_dphi->Fill(dphi);
0418
0419 for (int ip = 0; ip < 3; ip++)
0420 {
0421
0422 if (dphi <= dphi_cut[ip])
0423 break;
0424
0425 auto h_Aj = dynamic_cast<TH1*>(hm->getHisto(Form("h_truth_Aj_%d_dp%d", i, ip)));
0426 h_Aj->Fill(Aj_pt);
0427 auto h_xj = dynamic_cast<TH1*>(hm->getHisto(Form("h_truth_xj_%d_dp%d", i, ip)));
0428 h_xj->Fill(xj_pt);
0429
0430 }
0431 }
0432 }
0433 }
0434 }
0435 if ((int) truth_smear_jets.size() >= 2)
0436 {
0437 std::sort(truth_smear_jets.begin(), truth_smear_jets.end(), [] (auto a, auto b){ return a.pt > b.pt; });
0438 float eta1 = truth_smear_jets.at(0).eta;
0439 float phi1 = truth_smear_jets.at(0).phi;
0440 float eta2 = truth_smear_jets.at(1).eta;
0441 float phi2 = truth_smear_jets.at(1).phi;
0442 if (fabs(eta1) < 0.7 && fabs(eta2) < 0.7)
0443 {
0444
0445 float dphi = phi1 - phi2;
0446
0447 if (dphi > TMath::Pi())
0448 {
0449 dphi = dphi - 2.*TMath::Pi();
0450 }
0451 if (dphi <= -1*TMath::Pi())
0452 {
0453 dphi = dphi + 2. * TMath::Pi();
0454 }
0455
0456 dphi = fabs(dphi);
0457
0458 if (dphi > TMath::Pi())
0459 return;
0460
0461 if (Verbosity())
0462 {
0463 std::cout << "Truth Dijet" << std::endl;
0464 std::cout << "Leading Jet" << std::endl;
0465 int jl = 0;
0466 std::cout << " pt = " << truth_smear_jets.at(jl).pt << std::endl;
0467 std::cout << " eta = " << truth_smear_jets.at(jl).eta << std::endl;
0468 std::cout << " phi = " << truth_smear_jets.at(jl).phi << std::endl;
0469 std::cout << "Subleading Jet" << std::endl;
0470 jl++;
0471 std::cout << " pt = " << truth_smear_jets.at(jl).pt << std::endl;
0472 std::cout << " eta = " << truth_smear_jets.at(jl).eta << std::endl;
0473 std::cout << " phi = " << truth_smear_jets.at(jl).phi << std::endl;
0474 std::cout << " " << std::endl;
0475 }
0476
0477
0478 if (truth_smear_jets.at(0).pt > pt_cut1[0] && truth_smear_jets.at(1).pt > pt_cut2[0])
0479 {
0480 auto h_pt1_pt2 = dynamic_cast<TH2*>(hm->getHisto("h_truth_smear_pt1_pt2"));
0481 h_pt1_pt2->Fill(truth_smear_jets.at(0).pt, truth_smear_jets.at(1).pt);
0482 h_pt1_pt2->Fill(truth_smear_jets.at(1).pt, truth_smear_jets.at(0).pt);
0483 }
0484
0485 float Aj_pt = (truth_smear_jets.at(0).pt - truth_smear_jets.at(1).pt)/(truth_smear_jets.at(0).pt + truth_smear_jets.at(1).pt);
0486 float xj_pt = (truth_smear_jets.at(1).pt / truth_smear_jets.at(0).pt);
0487
0488 if (Verbosity())
0489 {
0490 std::cout << "---- DIJET ---- " << std::endl;
0491 std::cout << "Calculated: " << std::endl;
0492 std::cout << " A_j = " << Aj_pt << std::endl;
0493 std::cout << " x_j = " << xj_pt << std::endl;
0494 std::cout << " dPhi = " << dphi << std::endl;
0495 std::cout << " ------------------------ " << std::endl;
0496 }
0497
0498 for (int i = 0; i < 3; i++)
0499 {
0500 if (truth_smear_jets.at(0).pt > pt_cut1[i] && truth_smear_jets.at(1).pt > pt_cut2[i])
0501 {
0502 auto h_dphi = dynamic_cast<TH1*>(hm->getHisto(Form("h_truth_smear_dphi_%d", i)));
0503 h_dphi->Fill(dphi);
0504
0505 for (int ip = 0; ip < 3; ip++)
0506 {
0507
0508 if (dphi <= dphi_cut[ip])
0509 break;
0510 auto h_Aj = dynamic_cast<TH1*>(hm->getHisto(Form("h_truth_smear_Aj_%d_dp%d", i, ip)));
0511 h_Aj->Fill(Aj_pt);
0512 auto h_xj = dynamic_cast<TH1*>(hm->getHisto(Form("h_truth_smear_xj_%d_dp%d", i, ip)));
0513 h_xj->Fill(xj_pt);
0514
0515 }
0516 }
0517 }
0518 }
0519 }
0520 return;
0521 }
0522 int DijetHistoMaker::process_event(PHCompositeNode *topNode)
0523 {
0524
0525
0526 _i_event++;
0527
0528 if (!isSim)
0529 {
0530 triggeranalyzer->decodeTriggers(topNode);
0531 std::cout << triggeranalyzer->getTriggerName(17) << " : " << triggeranalyzer->getTriggerScalers(m_jet_triggernames[1]) << std::endl;
0532 }
0533 RawTowerGeomContainer *tower_geomIH = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALIN");
0534 RawTowerGeomContainer *tower_geomOH = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALOUT");
0535
0536
0537 std::string recoJetName = "AntiKt_Tower_r04_Sub1";
0538
0539 JetContainer *jets_4 = findNode::getClass<JetContainer>(topNode, recoJetName);
0540 recoJetName = "AntiKt_Truth_r04";
0541 JetContainer *jets_truth_4 = findNode::getClass<JetContainer>(topNode, recoJetName);
0542
0543 if (jets_truth_4)
0544 {
0545 if (Verbosity()) std::cout << __LINE__ << std::endl;
0546 process_truth_jets(jets_truth_4);
0547 if (Verbosity()) std::cout << __LINE__ << std::endl;
0548
0549 }
0550
0551 TrashInfo *trashinfo = findNode::getClass<TrashInfo>(topNode, "TrashInfo");
0552 if (!jets_4)
0553 {
0554 std::cout << " Nothing in the jets"<<std::endl;
0555 return Fun4AllReturnCodes::ABORTRUN;
0556 }
0557 int trash = 0;
0558 if (trashinfo)
0559 {
0560 trash = (trashinfo->isTrash()? 1: 0);
0561 }
0562
0563 TowerInfoContainer *hcalin_towers = findNode::getClass<TowerInfoContainer>(topNode, m_calo_nodename + "_HCALIN");
0564 TowerInfoContainer *hcalout_towers = findNode::getClass<TowerInfoContainer>(topNode, m_calo_nodename + "_HCALOUT");
0565 TowerInfoContainer *emcalre_towers = findNode::getClass<TowerInfoContainer>(topNode, m_calo_nodename + "_CEMC_RETOWER");
0566
0567 float background_v2 = 0;
0568 float background_Psi2 = 0;
0569 bool has_tower_background = false;
0570 TowerBackground *towBack = findNode::getClass<TowerBackground>(topNode, "TowerInfoBackground_Sub2");
0571 if (!towBack)
0572 {
0573
0574 }
0575 else
0576 {
0577 has_tower_background = true;
0578 background_v2 = towBack->get_v2();
0579 background_Psi2 = towBack->get_Psi2();
0580 }
0581
0582 if (!jets_4)
0583 {
0584 return Fun4AllReturnCodes::ABORTRUN;
0585 }
0586
0587 std::vector<jetty> jets{};
0588 for (auto jet : *jets_4)
0589 {
0590 float jet_total_eT = 0;
0591 float eFrac_ihcal = 0;
0592 float eFrac_ohcal = 0;
0593 float eFrac_emcal = 0;
0594
0595 for (auto comp : jet->get_comp_vec())
0596 {
0597
0598 unsigned int channel = comp.second;
0599 TowerInfo *tower;
0600 float tower_eT = 0;
0601 if (comp.first == 26 || comp.first == 30)
0602 {
0603 tower = hcalin_towers->get_tower_at_channel(channel);
0604 if (!tower || !tower_geomIH)
0605 {
0606 continue;
0607 }
0608 if(!tower->get_isGood()) continue;
0609
0610 unsigned int calokey = hcalin_towers->encode_key(channel);
0611
0612 int ieta = hcalin_towers->getTowerEtaBin(calokey);
0613 int iphi = hcalin_towers->getTowerPhiBin(calokey);
0614 const RawTowerDefs::keytype key = RawTowerDefs::encode_towerid(RawTowerDefs::CalorimeterId::HCALIN, ieta, iphi);
0615 float tower_phi = tower_geomIH->get_tower_geometry(key)->get_phi();
0616 float tower_eta = tower_geomIH->get_tower_geometry(key)->get_eta();
0617 tower_eT = tower->get_energy() / std::cosh(tower_eta);
0618
0619 if (comp.first == 30)
0620 {
0621 if (has_tower_background)
0622 {
0623 float UE = towBack->get_UE(1).at(ieta);
0624 float tower_UE = UE * (1 + 2 * background_v2 * std::cos(2 * (tower_phi - background_Psi2)));
0625 tower_eT = (tower->get_energy() - tower_UE) / std::cosh(tower_eta);
0626 }
0627 }
0628
0629 eFrac_ihcal += tower_eT;
0630 jet_total_eT += tower_eT;
0631 }
0632 else if (comp.first == 27 || comp.first == 31)
0633 {
0634 tower = hcalout_towers->get_tower_at_channel(channel);
0635
0636 if (!tower || !tower_geomOH)
0637 {
0638 continue;
0639 }
0640
0641 unsigned int calokey = hcalout_towers->encode_key(channel);
0642 int ieta = hcalout_towers->getTowerEtaBin(calokey);
0643 int iphi = hcalout_towers->getTowerPhiBin(calokey);
0644 const RawTowerDefs::keytype key = RawTowerDefs::encode_towerid(RawTowerDefs::CalorimeterId::HCALOUT, ieta, iphi);
0645 float tower_phi = tower_geomOH->get_tower_geometry(key)->get_phi();
0646 float tower_eta = tower_geomOH->get_tower_geometry(key)->get_eta();
0647 tower_eT = tower->get_energy() / std::cosh(tower_eta);
0648
0649 if (comp.first == 31)
0650 {
0651 if (has_tower_background)
0652 {
0653 float UE = towBack->get_UE(2).at(ieta);
0654 float tower_UE = UE * (1 + 2 * background_v2 * std::cos(2 * (tower_phi - background_Psi2)));
0655 tower_eT = (tower->get_energy() - tower_UE) / std::cosh(tower_eta);
0656 }
0657 }
0658
0659 eFrac_ohcal += tower_eT;
0660 jet_total_eT += tower_eT;
0661 }
0662 else if (comp.first == 28 || comp.first == 29)
0663 {
0664 tower = emcalre_towers->get_tower_at_channel(channel);
0665
0666 if (!tower || !tower_geomIH)
0667 {
0668 continue;
0669 }
0670
0671 unsigned int calokey = emcalre_towers->encode_key(channel);
0672 int ieta = emcalre_towers->getTowerEtaBin(calokey);
0673 int iphi = emcalre_towers->getTowerPhiBin(calokey);
0674 const RawTowerDefs::keytype key = RawTowerDefs::encode_towerid(RawTowerDefs::CalorimeterId::HCALIN, ieta, iphi);
0675 float tower_phi = tower_geomIH->get_tower_geometry(key)->get_phi();
0676 float tower_eta = tower_geomIH->get_tower_geometry(key)->get_eta();
0677 tower_eT = tower->get_energy() / std::cosh(tower_eta);
0678
0679 if (comp.first == 29)
0680 {
0681 if (has_tower_background)
0682 {
0683 float UE = towBack->get_UE(0).at(ieta);
0684 float tower_UE = UE * (1 + 2 * background_v2 * std::cos(2 * (tower_phi - background_Psi2)));
0685 tower_eT = (tower->get_energy() - tower_UE) / std::cosh(tower_eta);
0686 }
0687 }
0688
0689
0690 eFrac_emcal += tower_eT;
0691 jet_total_eT += tower_eT;
0692 }
0693 }
0694 struct jetty myjet;
0695 myjet.pt = jet->get_pt();
0696 myjet.eta = jet->get_eta();
0697 myjet.phi = jet->get_phi();
0698 myjet.et = jet_total_eT;
0699 myjet.emcal = eFrac_emcal;
0700 myjet.ihcal = eFrac_ihcal;
0701 myjet.ohcal = eFrac_ohcal;
0702 jets.push_back(myjet);
0703 }
0704
0705 if (jets.size() < 2)
0706 return Fun4AllReturnCodes::EVENT_OK;
0707
0708
0709 auto pt_ordered = jets;
0710 std::sort(jets.begin(), jets.end(), [](auto a, auto b) { return a.et > b.et; });
0711 std::sort(pt_ordered.begin(), pt_ordered.end(), [](auto a, auto b) { return a.pt > b.pt; });
0712 if (pt_ordered.at(0).pt < 10 && jets.at(0).et < 10) return Fun4AllReturnCodes::EVENT_OK;
0713
0714 const RawTowerDefs::keytype bottomkey = RawTowerDefs::encode_towerid(RawTowerDefs::CalorimeterId::HCALIN, 0, 0);
0715 const RawTowerDefs::keytype topkey = RawTowerDefs::encode_towerid(RawTowerDefs::CalorimeterId::HCALIN, 23, 0);
0716 float tower_eta_low = tower_geomIH->get_tower_geometry(bottomkey)->get_eta() + 0.4;
0717 float tower_eta_high = tower_geomIH->get_tower_geometry(topkey)->get_eta() - 0.4;
0718
0719
0720 for (int ij = 0; ij < 2; ij++)
0721 {
0722 if (jets.at(ij).eta > tower_eta_high || jets.at(ij).eta < tower_eta_low)
0723 return Fun4AllReturnCodes::EVENT_OK;
0724 }
0725
0726 float dphi = jets.at(0).phi - jets.at(1).phi;
0727 if (dphi > TMath::Pi())
0728 {
0729 dphi = dphi - 2.*TMath::Pi();
0730 }
0731 if (dphi <= -1*TMath::Pi())
0732 {
0733 dphi = dphi + 2. * TMath::Pi();
0734 }
0735
0736 dphi = fabs(dphi);
0737
0738
0739
0740 bool isopt = true;
0741 bool isoet = true;
0742 if (jets.size() > 2)
0743 {
0744 for (int ijet = 0; ijet < (int) jets.size() - 2; ijet++)
0745 {
0746 if (jets.at(ijet+2).pt < 2 ) continue;
0747
0748 for (int ij = 0; ij < 2; ij++)
0749 {
0750 double dR = get_Dr(jets.at(ij), jets.at(ijet+2));
0751 isoet &= (dR > isocut);
0752 }
0753 }
0754 }
0755
0756 if (pt_ordered.size() > 2)
0757 {
0758 for (int ijet = 0; ijet < (int) pt_ordered.size() - 2; ijet++)
0759 {
0760 if (pt_ordered.at(ijet+2).pt < 2 ) continue;
0761
0762 for (int ij = 0; ij < 2; ij++)
0763 {
0764 double dR = get_Dr(pt_ordered.at(ij), pt_ordered.at(ijet+2));
0765 isopt &= (dR > isocut);
0766 }
0767 }
0768 }
0769 if (dphi > TMath::Pi())
0770 return Fun4AllReturnCodes::EVENT_OK;
0771
0772 if (Verbosity())
0773 {
0774
0775 std::cout << "Leading Jet" << std::endl;
0776 int jl = 0;
0777 std::cout << " pt = " << jets.at(jl).pt << std::endl;
0778 std::cout << " et = " << jets.at(jl).et << std::endl;
0779 std::cout << " eta = " << jets.at(jl).eta << std::endl;
0780 std::cout << " phi = " << jets.at(jl).phi << std::endl;
0781 std::cout << "Subleading Jet" << std::endl;
0782 jl++;
0783 std::cout << " pt = " << jets.at(jl).pt << std::endl;
0784 std::cout << " et = " << jets.at(jl).et << std::endl;
0785 std::cout << " eta = " << jets.at(jl).eta << std::endl;
0786 std::cout << " phi = " << jets.at(jl).phi << std::endl;
0787 std::cout << " " << std::endl;
0788
0789
0790
0791
0792
0793 }
0794
0795 if (trash)
0796 {
0797 if (pt_ordered.at(0).pt > pt_cut1[0] && pt_ordered.at(1).pt > pt_cut2[0])
0798 {
0799 auto h_trash_pt1_pt2 = dynamic_cast<TH2*>(hm->getHisto("h_trash_pt1_pt2"));
0800 h_trash_pt1_pt2->Fill(pt_ordered.at(0).pt, pt_ordered.at(1).pt);
0801 h_trash_pt1_pt2->Fill(pt_ordered.at(1).pt, pt_ordered.at(0).pt);
0802 }
0803 if (jets.at(0).et > pt_cut1[0] && jets.at(1).et > pt_cut2[0])
0804 {
0805 auto h_trash_et1_et2 = dynamic_cast<TH2*>(hm->getHisto("h_trash_et1_et2"));
0806 h_trash_et1_et2->Fill(jets.at(0).et, jets.at(1).et);
0807 h_trash_et1_et2->Fill(jets.at(1).et, jets.at(0).et);
0808 }
0809 if (triggeranalyzer->didTriggerFire("MBD N&S >= 1"))
0810 {
0811 if (pt_ordered.at(0).pt > pt_cut1[0] && pt_ordered.at(1).pt > pt_cut2[0])
0812 {
0813 auto h_trash_pt1_pt2_mbd = dynamic_cast<TH2*>(hm->getHisto("h_trash_pt1_pt2_mbd"));
0814 h_trash_pt1_pt2_mbd->Fill(pt_ordered.at(0).pt, pt_ordered.at(1).pt);
0815 h_trash_pt1_pt2_mbd->Fill(pt_ordered.at(1).pt, pt_ordered.at(0).pt);
0816 }
0817 if (jets.at(0).et > pt_cut1[0] && jets.at(1).et > pt_cut2[0])
0818 {
0819 auto h_trash_et1_et2_mbd = dynamic_cast<TH2*>(hm->getHisto("h_trash_et1_et2_mbd"));
0820 h_trash_et1_et2_mbd->Fill(jets.at(0).et, jets.at(1).et);
0821 h_trash_et1_et2_mbd->Fill(jets.at(1).et, jets.at(0).et);
0822 }
0823 }
0824 for (int i = 0; i < 8; i++)
0825 {
0826 if (triggeranalyzer->didTriggerFire(m_jet_triggernames[i]))
0827 {
0828 if (pt_ordered.at(0).pt > pt_cut1[0] && pt_ordered.at(1).pt > pt_cut2[0])
0829 {
0830 auto h_trash_pt1_pt2_gl1 = dynamic_cast<TH2*>(hm->getHisto(Form("h_trash_pt1_pt2_gl1_%d", i)));
0831 h_trash_pt1_pt2_gl1->Fill(pt_ordered.at(0).pt, pt_ordered.at(1).pt);
0832 h_trash_pt1_pt2_gl1->Fill(pt_ordered.at(1).pt, pt_ordered.at(0).pt);
0833 }
0834 if (jets.at(0).et > pt_cut1[0] && jets.at(1).et > pt_cut2[0])
0835 {
0836 auto h_trash_et1_et2_gl1 = dynamic_cast<TH2*>(hm->getHisto(Form("h_trash_et1_et2_gl1_%d", i)));
0837 h_trash_et1_et2_gl1->Fill(jets.at(0).et, jets.at(1).et);
0838 h_trash_et1_et2_gl1->Fill(jets.at(1).et, jets.at(0).et);
0839 }
0840 }
0841 }
0842
0843 float Aj_pt = (pt_ordered.at(0).pt - pt_ordered.at(1).pt)/(pt_ordered.at(0).pt + pt_ordered.at(1).pt);
0844 float Aj_et = (jets.at(0).et - jets.at(1).et)/(jets.at(0).et + jets.at(1).et);
0845 float xj_pt = (pt_ordered.at(1).pt / pt_ordered.at(0).pt);
0846 float xj_et = (jets.at(1).et / jets.at(0).et);
0847
0848 if (Verbosity())
0849 {
0850 std::cout << "---- TRASH ---- " << std::endl;
0851 std::cout << "Calculated: " << std::endl;
0852 std::cout << " A_j = " << Aj_pt << " ( " << Aj_et << " ) "<< std::endl;
0853 std::cout << " x_j = " << xj_pt << " ( " << xj_et << " ) "<< std::endl;
0854 std::cout << " dPhi = " << dphi << std::endl;
0855 std::cout << " ------------------------ " << std::endl;
0856 }
0857 bool trashalreadyfilled = false;
0858 for (int i = 0; i < 3; i++)
0859 {
0860 if (pt_ordered.at(0).pt > pt_cut1[i] && pt_ordered.at(1).pt > pt_cut2[i])
0861 {
0862 auto h_trash_dphi = dynamic_cast<TH1*>(hm->getHisto(Form("h_trash_dphi_%d", i)));
0863 h_trash_dphi->Fill(dphi);
0864 if (!isSim)
0865 {
0866 if (triggeranalyzer->didTriggerFire("MBD N&S >= 1"))
0867 {
0868 auto h_trash_dphi_mb = dynamic_cast<TH1*>(hm->getHisto(Form("h_trash_dphi_mb_%d", i)));
0869 h_trash_dphi_mb->Fill(dphi);
0870 }
0871
0872 for (int j = 0; j < 8; j++)
0873 {
0874 if (triggeranalyzer->didTriggerFire(m_jet_triggernames[j]))
0875 {
0876 auto h_trash_dphi_gl1 = dynamic_cast<TH1*>(hm->getHisto(Form("h_trash_dphi_gl1_%d_%d", j, i)));
0877 h_trash_dphi_gl1->Fill(dphi);
0878 }
0879 }
0880 }
0881
0882 if (i == 0 && drawDijets)
0883 {
0884 trashalreadyfilled = true;
0885 std::cout << " FOUND " << std::endl;
0886 dijeteventdisplay->FillEvent(topNode, Aj_et, dphi);
0887 }
0888 for (int ip = 0; ip < 3; ip++)
0889 {
0890 if (dphi <= dphi_cut[ip])
0891 break;
0892 auto h_trash_Aj = dynamic_cast<TH1*>(hm->getHisto(Form("h_trash_Aj_%d_dp%d", i, ip)));
0893 h_trash_Aj->Fill(Aj_pt);
0894 auto h_trash_xj = dynamic_cast<TH1*>(hm->getHisto(Form("h_trash_xj_%d_dp%d", i, ip)));
0895 h_trash_xj->Fill(xj_pt);
0896
0897 if (!isSim)
0898 {
0899 if (triggeranalyzer->didTriggerFire("MBD N&S >= 1"))
0900 {
0901 auto h_trash_Aj_mb = dynamic_cast<TH1*>(hm->getHisto(Form("h_trash_Aj_mb_%d_dp%d", i, ip)));
0902 h_trash_Aj_mb->Fill(Aj_pt);
0903 auto h_trash_xj_mb = dynamic_cast<TH1*>(hm->getHisto(Form("h_trash_xj_mb_%d_dp%d", i, ip)));
0904 h_trash_xj_mb->Fill(xj_pt);
0905 }
0906 }
0907 for (int j = 0; j < 8; j++)
0908 {
0909 if (triggeranalyzer->didTriggerFire(m_jet_triggernames[j]))
0910 {
0911 auto h_trash_Aj_gl1 = dynamic_cast<TH1*>(hm->getHisto(Form("h_trash_Aj_gl1_%d_%d_dp%d", j, i, ip)));
0912 h_trash_Aj_gl1->Fill(Aj_pt);
0913 auto h_trash_xj_gl1 = dynamic_cast<TH1*>(hm->getHisto(Form("h_trash_xj_gl1_%d_%d_dp%d", j, i, ip)));
0914 h_trash_xj_gl1->Fill(xj_pt);
0915 }
0916 }
0917
0918 }
0919 }
0920 if (jets.at(0).et > pt_cut1[i] && jets.at(1).et > pt_cut2[i])
0921 {
0922 auto h_trash_dphi = dynamic_cast<TH1*>(hm->getHisto(Form("h_trash_dphi_et_%d", i)));
0923 h_trash_dphi->Fill(dphi);
0924
0925 if (!isSim)
0926 {
0927 if (triggeranalyzer->didTriggerFire("MBD N&S >= 1"))
0928 {
0929
0930 auto h_trash_dphi_et_mb = dynamic_cast<TH1*>(hm->getHisto(Form("h_trash_dphi_et_mb_%d", i)));
0931 h_trash_dphi_et_mb->Fill(dphi);
0932 }
0933 for (int j = 0; j < 8; j++)
0934 {
0935 if (triggeranalyzer->didTriggerFire(m_jet_triggernames[j]))
0936 {
0937 auto h_trash_dphi_et_gl1 = dynamic_cast<TH1*>(hm->getHisto(Form("h_trash_dphi_et_gl1_%d_%d", j, i)));
0938 h_trash_dphi_et_gl1->Fill(dphi);
0939 }
0940 }
0941 }
0942
0943 if (i == 0 && !trashalreadyfilled && drawDijets)
0944 {
0945 std::cout << " FOUND " << std::endl;
0946 dijeteventdisplay->FillEvent(topNode, Aj_et, dphi);
0947 }
0948 for (int ip = 0 ; ip < 3; ip++)
0949 {
0950 if (dphi <= dphi_cut[ip]) break;
0951 auto h_trash_Aj = dynamic_cast<TH1*>(hm->getHisto(Form("h_trash_Aj_et_%d_dp%d", i, ip)));
0952 h_trash_Aj->Fill(Aj_et);
0953 auto h_trash_xj = dynamic_cast<TH1*>(hm->getHisto(Form("h_trash_xj_et_%d_dp%d", i, ip)));
0954 h_trash_xj->Fill(xj_et);
0955
0956 if (!isSim)
0957 {
0958 if (triggeranalyzer->didTriggerFire("MBD N&S >= 1"))
0959 {
0960 auto h_trash_Aj_et_mb = dynamic_cast<TH1*>(hm->getHisto(Form("h_trash_Aj_et_mb_%d_dp%d", i, ip)));
0961 h_trash_Aj_et_mb->Fill(Aj_et);
0962 auto h_trash_xj_et_mb = dynamic_cast<TH1*>(hm->getHisto(Form("h_trash_xj_et_mb_%d_dp%d", i, ip)));
0963 h_trash_xj_et_mb->Fill(xj_et);
0964 }
0965 for (int j = 0; j < 8; j++)
0966 {
0967 if (triggeranalyzer->didTriggerFire(m_jet_triggernames[j]))
0968 {
0969 auto h_trash_Aj_et_gl1 = dynamic_cast<TH1*>(hm->getHisto(Form("h_trash_Aj_et_gl1_%d_%d_dp%d", j, i, ip)));
0970 h_trash_Aj_et_gl1->Fill(Aj_et);
0971 auto h_trash_xj_et_gl1 = dynamic_cast<TH1*>(hm->getHisto(Form("h_trash_xj_et_gl1_%d_%d_dp%d", j, i, ip)));
0972 h_trash_xj_et_gl1->Fill(xj_et);
0973 }
0974 }
0975 }
0976 }
0977 }
0978 }
0979 return Fun4AllReturnCodes::EVENT_OK;
0980 }
0981
0982 if (pt_ordered.at(0).pt > pt_cut1[0] && pt_ordered.at(1).pt > pt_cut2[0])
0983 {
0984 auto h_pt1_pt2 = dynamic_cast<TH2*>(hm->getHisto("h_pt1_pt2"));
0985 h_pt1_pt2->Fill(pt_ordered.at(0).pt, pt_ordered.at(1).pt);
0986 h_pt1_pt2->Fill(pt_ordered.at(1).pt, pt_ordered.at(0).pt);
0987 }
0988 if (jets.at(0).et > pt_cut1[0] && jets.at(1).et > pt_cut2[0])
0989 {
0990 auto h_et1_et2 = dynamic_cast<TH2*>(hm->getHisto("h_et1_et2"));
0991 h_et1_et2->Fill(jets.at(0).et, jets.at(1).et);
0992 h_et1_et2->Fill(jets.at(1).et, jets.at(0).et);
0993 }
0994 if (!isSim)
0995 {
0996 if (triggeranalyzer->didTriggerFire("MBD N&S >= 1"))
0997 {
0998 if (pt_ordered.at(0).pt > pt_cut1[0] && pt_ordered.at(1).pt > pt_cut2[0])
0999 {
1000 auto h_pt1_pt2_mbd = dynamic_cast<TH2*>(hm->getHisto("h_pt1_pt2_mbd"));
1001 h_pt1_pt2_mbd->Fill(pt_ordered.at(0).pt, pt_ordered.at(1).pt);
1002 h_pt1_pt2_mbd->Fill(pt_ordered.at(1).pt, pt_ordered.at(0).pt);
1003 }
1004 if (jets.at(0).et > pt_cut1[0] && jets.at(1).et > pt_cut2[0])
1005 {
1006 auto h_et1_et2_mbd = dynamic_cast<TH2*>(hm->getHisto("h_et1_et2_mbd"));
1007 h_et1_et2_mbd->Fill(jets.at(0).et, jets.at(1).et);
1008 h_et1_et2_mbd->Fill(jets.at(1).et, jets.at(0).et);
1009 }
1010 }
1011 for (int i = 0; i < 8; i++)
1012 {
1013 if (triggeranalyzer->didTriggerFire(m_jet_triggernames[i]))
1014 {
1015 if (pt_ordered.at(0).pt > pt_cut1[0] && pt_ordered.at(1).pt > pt_cut2[0])
1016 {
1017 auto h_pt1_pt2_gl1 = dynamic_cast<TH2*>(hm->getHisto(Form("h_pt1_pt2_gl1_%d", i)));
1018 h_pt1_pt2_gl1->Fill(pt_ordered.at(0).pt, pt_ordered.at(1).pt);
1019 h_pt1_pt2_gl1->Fill(pt_ordered.at(1).pt, pt_ordered.at(0).pt);
1020 }
1021 if (jets.at(0).et > pt_cut1[0] && jets.at(1).et > pt_cut2[0])
1022 {
1023 auto h_et1_et2_gl1 = dynamic_cast<TH2*>(hm->getHisto(Form("h_et1_et2_gl1_%d", i)));
1024 h_et1_et2_gl1->Fill(jets.at(0).et, jets.at(1).et);
1025 h_et1_et2_gl1->Fill(jets.at(1).et, jets.at(0).et);
1026 }
1027 }
1028 }
1029 }
1030
1031 float Aj_pt = (pt_ordered.at(0).pt - pt_ordered.at(1).pt)/(pt_ordered.at(0).pt + pt_ordered.at(1).pt);
1032 float Aj_et = (jets.at(0).et - jets.at(1).et)/(jets.at(0).et + jets.at(1).et);
1033 float xj_pt = (pt_ordered.at(1).pt / pt_ordered.at(0).pt);
1034 float xj_et = (jets.at(1).et / jets.at(0).et);
1035
1036 if (Verbosity())
1037 {
1038 std::cout << "---- DIJET ---- " << std::endl;
1039 std::cout << "Calculated: " << std::endl;
1040 std::cout << " A_j = " << Aj_pt << " ( " << Aj_et << " ) "<< std::endl;
1041 std::cout << " x_j = " << xj_pt << " ( " << xj_et << " ) "<< std::endl;
1042 std::cout << " dPhi = " << dphi << std::endl;
1043 std::cout << " ------------------------ " << std::endl;
1044 }
1045 bool alreadyfilled = false;
1046 for (int i = 0; i < 3; i++)
1047 {
1048 if (pt_ordered.at(0).pt > pt_cut1[i] && pt_ordered.at(1).pt > pt_cut2[i])
1049 {
1050 auto h_dphi = dynamic_cast<TH1*>(hm->getHisto(Form("h_dphi_%d", i)));
1051 h_dphi->Fill(dphi);
1052
1053 if (!isSim)
1054 {
1055 if (i == 0 && drawDijets)
1056 {
1057 alreadyfilled = true;
1058 std::cout << " FOUND " << std::endl;
1059 dijeteventdisplay->FillEvent(topNode, Aj_et, dphi);
1060 }
1061 if (triggeranalyzer->didTriggerFire("MBD N&S >= 1"))
1062 {
1063 auto h_dphi_mb = dynamic_cast<TH1*>(hm->getHisto(Form("h_dphi_mb_%d", i)));
1064 h_dphi_mb->Fill(dphi);
1065 }
1066
1067 for (int j = 0; j < 8; j++)
1068 {
1069 if (triggeranalyzer->didTriggerFire(m_jet_triggernames[j]))
1070 {
1071 auto h_dphi_gl1 = dynamic_cast<TH1*>(hm->getHisto(Form("h_dphi_gl1_%d_%d", j, i)));
1072 h_dphi_gl1->Fill(dphi);
1073 }
1074 }
1075 }
1076 for (int ip = 0; ip < 3; ip++)
1077 {
1078 if (dphi <= dphi_cut[ip]) break;
1079
1080 auto h_Aj = dynamic_cast<TH1*>(hm->getHisto(Form("h_Aj_%d_dp%d", i, ip)));
1081 h_Aj->Fill(Aj_pt);
1082 auto h_xj = dynamic_cast<TH1*>(hm->getHisto(Form("h_xj_%d_dp%d", i, ip)));
1083 h_xj->Fill(xj_pt);
1084
1085 auto h_Ajiso = dynamic_cast<TH1*>(hm->getHisto(Form("h_Aj_iso%d_%d_dp%d", (isopt?1:0), i, ip)));
1086 h_Ajiso->Fill(Aj_pt);
1087 auto h_xjiso = dynamic_cast<TH1*>(hm->getHisto(Form("h_xj_iso%d_%d_dp%d", (isopt?1:0), i, ip)));
1088 h_xjiso->Fill(xj_pt);
1089
1090 if (!isSim)
1091 {
1092 if (i == 0 && drawDijets)
1093 {
1094 alreadyfilled = true;
1095 std::cout << " FOUND " << std::endl;
1096 dijeteventdisplay->FillEvent(topNode, Aj_et, dphi);
1097 }
1098 if (triggeranalyzer->didTriggerFire("MBD N&S >= 1"))
1099 {
1100 auto h_Aj_mb = dynamic_cast<TH1*>(hm->getHisto(Form("h_Aj_mb_%d_dp%d", i, ip)));
1101 h_Aj_mb->Fill(Aj_pt);
1102 auto h_xj_mb = dynamic_cast<TH1*>(hm->getHisto(Form("h_xj_mb_%d_dp%d", i, ip)));
1103 h_xj_mb->Fill(xj_pt);
1104 }
1105
1106 for (int j = 0; j < 8; j++)
1107 {
1108 if (triggeranalyzer->didTriggerFire(m_jet_triggernames[j]))
1109 {
1110 auto h_Aj_gl1 = dynamic_cast<TH1*>(hm->getHisto(Form("h_Aj_gl1_%d_%d_dp%d", j, i, ip)));
1111 h_Aj_gl1->Fill(Aj_pt);
1112 auto h_xj_gl1 = dynamic_cast<TH1*>(hm->getHisto(Form("h_xj_gl1_%d_%d_dp%d", j, i, ip)));
1113 h_xj_gl1->Fill(xj_pt);
1114 }
1115 }
1116 }
1117 }
1118 }
1119 if (jets.at(0).et > pt_cut1[i] && jets.at(1).et > pt_cut2[i])
1120 {
1121 auto h_dphi = dynamic_cast<TH1*>(hm->getHisto(Form("h_dphi_et_%d", i)));
1122 h_dphi->Fill(dphi);
1123 if (!isSim)
1124 {
1125 auto h_dphi_et_mb = dynamic_cast<TH1*>(hm->getHisto(Form("h_dphi_et_mb_%d", i)));
1126 h_dphi_et_mb->Fill(dphi);
1127 for (int j = 0; j < 8; j++)
1128 {
1129 if (triggeranalyzer->didTriggerFire(m_jet_triggernames[j]))
1130 {
1131 auto h_dphi_et_gl1 = dynamic_cast<TH1*>(hm->getHisto(Form("h_dphi_et_gl1_%d_%d", j, i)));
1132 h_dphi_et_gl1->Fill(dphi);
1133 }
1134 }
1135
1136 }
1137 for (int ip = 0 ; ip < 3; ip++)
1138 {
1139 if (dphi <= dphi_cut[ip]) break;
1140
1141 auto h_Aj = dynamic_cast<TH1*>(hm->getHisto(Form("h_Aj_et_%d_dp%d", i, ip)));
1142 h_Aj->Fill(Aj_et);
1143 auto h_xj = dynamic_cast<TH1*>(hm->getHisto(Form("h_xj_et_%d_dp%d", i, ip)));
1144 h_xj->Fill(xj_et);
1145 auto h_Ajiso = dynamic_cast<TH1*>(hm->getHisto(Form("h_Aj_iso%d_%d_dp%d", (isoet?1:0), i, ip)));
1146 h_Ajiso->Fill(Aj_et);
1147 auto h_xjiso = dynamic_cast<TH1*>(hm->getHisto(Form("h_xj_iso%d_%d_dp%d", (isoet?1:0), i, ip)));
1148 h_xjiso->Fill(xj_et);
1149
1150 if (!isSim)
1151 {
1152 if (i == 0 && !alreadyfilled && drawDijets)
1153 {
1154 std::cout << " FOUND " << std::endl;
1155 dijeteventdisplay->FillEvent(topNode, Aj_et, dphi);
1156 }
1157 if (triggeranalyzer->didTriggerFire("MBD N&S >= 1"))
1158 {
1159 auto h_Aj_et_mb = dynamic_cast<TH1*>(hm->getHisto(Form("h_Aj_et_mb_%d_dp%d", i, ip)));
1160 h_Aj_et_mb->Fill(Aj_et);
1161 auto h_xj_et_mb = dynamic_cast<TH1*>(hm->getHisto(Form("h_xj_et_mb_%d_dp%d", i, ip)));
1162 h_xj_et_mb->Fill(xj_et);
1163 }
1164
1165 for (int j = 0; j < 8; j++)
1166 {
1167 if (triggeranalyzer->didTriggerFire(m_jet_triggernames[j]))
1168 {
1169 auto h_Aj_et_gl1 = dynamic_cast<TH1*>(hm->getHisto(Form("h_Aj_et_gl1_%d_%d_dp%d", j, i, ip)));
1170 h_Aj_et_gl1->Fill(Aj_et);
1171 auto h_xj_et_gl1 = dynamic_cast<TH1*>(hm->getHisto(Form("h_xj_et_gl1_%d_%d_dp%d", j, i, ip)));
1172 h_xj_et_gl1->Fill(xj_et);
1173 }
1174 }
1175 }
1176 }
1177 }
1178 }
1179 return Fun4AllReturnCodes::EVENT_OK;
1180 }
1181
1182
1183
1184 void DijetHistoMaker::GetNodes(PHCompositeNode* topNode)
1185 {
1186
1187
1188 }
1189
1190 int DijetHistoMaker::ResetEvent(PHCompositeNode *topNode)
1191 {
1192
1193 return Fun4AllReturnCodes::EVENT_OK;
1194 }
1195
1196 double DijetHistoMaker::get_Dr(struct jetty jet1, struct jetty jet2)
1197 {
1198 float dphi = jet1.phi - jet2.phi;
1199
1200 if (dphi > TMath::Pi())
1201 {
1202 dphi = dphi - 2.*TMath::Pi();
1203 }
1204 if (dphi <= -1*TMath::Pi())
1205 {
1206 dphi = dphi + 2. * TMath::Pi();
1207 }
1208
1209 return sqrt(TMath::Power(dphi, 2) + TMath::Power(jet1.eta - jet2.eta, 2));
1210 }
1211
1212
1213 int DijetHistoMaker::EndRun(const int runnumber)
1214 {
1215 return Fun4AllReturnCodes::EVENT_OK;
1216 }
1217
1218
1219 int DijetHistoMaker::End(PHCompositeNode *topNode)
1220 {
1221 if (Verbosity() > 0)
1222 {
1223 std::cout << "DijetHistoMaker::End(PHCompositeNode *topNode) This is the End..." << std::endl;
1224 }
1225 std::cout<<"Total events: "<<_i_event<<std::endl;
1226 hm->dumpHistos(_foutname.c_str());
1227 if (!isSim && drawDijets) dijeteventdisplay->Dump();
1228 return Fun4AllReturnCodes::EVENT_OK;
1229 }
1230
1231
1232 int DijetHistoMaker::Reset(PHCompositeNode *topNode)
1233 {
1234 if (Verbosity() > 0)
1235 {
1236 std::cout << "DijetHistoMaker::Reset(PHCompositeNode *topNode) being Reset" << std::endl;
1237 }
1238 return Fun4AllReturnCodes::EVENT_OK;
1239 }
