File indexing completed on 2025-08-05 08:11:12
0001 #include "DijetTreeMaker.h"
0002
0003 #include <g4main/PHG4TruthInfoContainer.h>
0004 #include <g4main/PHG4Particle.h>
0005 #include <TMath.h>
0006 #include <calotrigger/TriggerRunInfov1.h>
0007
0008 #include <calobase/RawTowerGeom.h>
0009 #include <jetbackground/TowerBackground.h>
0010 #include <calobase/RawCluster.h>
0011 #include <calobase/RawClusterContainer.h>
0012 #include <calobase/RawClusterUtility.h>
0013 #include <calobase/RawTowerGeomContainer.h>
0014 #include <calobase/RawTower.h>
0015 #include <calobase/RawTowerContainer.h>
0016 #include <HepMC/SimpleVector.h>
0017
0018 #include <globalvertex/GlobalVertex.h>
0019 #include <globalvertex/GlobalVertexMap.h>
0020 #include <vector>
0021
0022 #include <fun4all/Fun4AllReturnCodes.h>
0023 #include <phool/PHCompositeNode.h>
0024 #include <phool/PHIODataNode.h>
0025 #include <phool/PHNode.h>
0026 #include <phool/PHNodeIterator.h>
0027 #include <phool/PHObject.h>
0028 #include <phool/getClass.h>
0029
0030
0031 #include <iostream>
0032
0033 #include <map>
0034
0035
0036 DijetTreeMaker::DijetTreeMaker(const std::string &name, const std::string &outfilename):
0037 SubsysReco(name)
0038
0039 {
0040 isSim = false;
0041 _foutname = outfilename;
0042 m_calo_nodename = "TOWERINFO_CALIB";
0043 }
0044
0045
0046 DijetTreeMaker::~DijetTreeMaker()
0047 {
0048
0049 }
0050
0051
0052 int DijetTreeMaker::Init(PHCompositeNode *topNode)
0053 {
0054
0055
0056 if (isSim)
0057 {
0058 pt_cut = 1.;
0059 }
0060
0061
0062 if (Verbosity()) std::cout << __FUNCTION__ << __LINE__<<std::endl;
0063 _f = new TFile( _foutname.c_str(), "RECREATE");
0064
0065 std::cout << " making a file = " << _foutname.c_str() << " , _f = " << _f << std::endl;
0066
0067 _tree = new TTree("ttree","a persevering date tree");
0068
0069 if (isSim)
0070 {
0071
0072
0073
0074
0075
0076
0077 _tree->Branch("truth_jet_pt_2",&b_truth_jet_pt_2);
0078 _tree->Branch("truth_jet_eta_2",&b_truth_jet_eta_2);
0079 _tree->Branch("truth_jet_phi_2",&b_truth_jet_phi_2);
0080 _tree->Branch("truth_jet_pt_4",&b_truth_jet_pt_4);
0081 _tree->Branch("truth_jet_eta_4",&b_truth_jet_eta_4);
0082 _tree->Branch("truth_jet_phi_4",&b_truth_jet_phi_4);
0083 _tree->Branch("truth_jet_pt_6",&b_truth_jet_pt_6);
0084 _tree->Branch("truth_jet_eta_6",&b_truth_jet_eta_6);
0085 _tree->Branch("truth_jet_phi_6",&b_truth_jet_phi_6);
0086
0087 }
0088
0089 _i_event = 0;
0090 if (save_calo)
0091 {
0092 _tree->Branch("emcal_good",&b_emcal_good);
0093 _tree->Branch("emcal_energy",&b_emcal_energy);
0094 _tree->Branch("emcal_time",&b_emcal_time);
0095 _tree->Branch("emcal_phibin",&b_emcal_phibin);
0096 _tree->Branch("emcal_etabin",&b_emcal_etabin);
0097 _tree->Branch("emcal_phi",&b_emcal_phi);
0098 _tree->Branch("emcal_eta",&b_emcal_eta);
0099
0100 _tree->Branch("hcalin_good",&b_hcalin_good);
0101 _tree->Branch("hcalin_energy",&b_hcalin_energy);
0102 _tree->Branch("hcalin_time",&b_hcalin_time);
0103 _tree->Branch("hcalin_phibin",&b_hcalin_phibin);
0104 _tree->Branch("hcalin_etabin",&b_hcalin_etabin);
0105 _tree->Branch("hcalin_phi",&b_hcalin_phi);
0106 _tree->Branch("hcalin_eta",&b_hcalin_eta);
0107
0108 _tree->Branch("hcalout_good",&b_hcalout_good);
0109 _tree->Branch("hcalout_energy",&b_hcalout_energy);
0110 _tree->Branch("hcalout_time",&b_hcalout_time);
0111 _tree->Branch("hcalout_phibin",&b_hcalout_phibin);
0112 _tree->Branch("hcalout_etabin",&b_hcalout_etabin);
0113 _tree->Branch("hcalout_phi",&b_hcalout_phi);
0114 _tree->Branch("hcalout_eta",&b_hcalout_eta);
0115 }
0116
0117 _tree->Branch("mbd_vertex_z", &b_vertex_z, "mbd_vertex_z/F");
0118
0119 _tree->Branch("njet_2", &b_njet_2, "njet_2/I");
0120 _tree->Branch("jet_pt_2", &b_jet_pt_2);
0121 _tree->Branch("jet_et_2", &b_jet_et_2);
0122 _tree->Branch("jet_eta_2", &b_jet_eta_2);
0123 _tree->Branch("jet_phi_2", &b_jet_phi_2);
0124 _tree->Branch("jet_emcal_2", &b_jet_emcal_2);
0125 _tree->Branch("jet_hcalout_2", &b_jet_hcalout_2);
0126 _tree->Branch("jet_hcalin_2", &b_jet_hcalin_2);
0127 _tree->Branch("jet_eccen_2", &b_jet_eccen_2);
0128
0129 _tree->Branch("njet_4", &b_njet_4, "njet_4/I");
0130 _tree->Branch("jet_pt_4", &b_jet_pt_4);
0131 _tree->Branch("jet_et_4", &b_jet_et_4);
0132 _tree->Branch("jet_eta_4", &b_jet_eta_4);
0133 _tree->Branch("jet_phi_4", &b_jet_phi_4);
0134 _tree->Branch("jet_emcal_4", &b_jet_emcal_4);
0135 _tree->Branch("jet_hcalout_4", &b_jet_hcalout_4);
0136 _tree->Branch("jet_hcalin_4", &b_jet_hcalin_4);
0137 _tree->Branch("jet_eccen_4", &b_jet_eccen_4);
0138
0139 _tree->Branch("njet_6", &b_njet_6, "njet_6/I");
0140 _tree->Branch("jet_pt_6", &b_jet_pt_6);
0141 _tree->Branch("jet_et_6", &b_jet_et_6);
0142 _tree->Branch("jet_eta_6", &b_jet_eta_6);
0143 _tree->Branch("jet_phi_6", &b_jet_phi_6);
0144 _tree->Branch("jet_emcal_6", &b_jet_emcal_6);
0145 _tree->Branch("jet_hcalout_6", &b_jet_hcalout_6);
0146 _tree->Branch("jet_hcalin_6", &b_jet_hcalin_6);
0147 _tree->Branch("jet_eccen_6", &b_jet_eccen_6);
0148
0149 _tree->Branch("ncluster", &b_ncluster);
0150 _tree->Branch("cluster_ecore", &b_cluster_ecore);
0151 _tree->Branch("cluster_e", &b_cluster_e);
0152 _tree->Branch("cluster_pt", &b_cluster_pt);
0153 _tree->Branch("cluster_chi2", &b_cluster_chi2);
0154 _tree->Branch("cluster_prob", &b_cluster_prob);
0155 _tree->Branch("cluster_eta", &b_cluster_eta);
0156 _tree->Branch("cluster_phi", &b_cluster_phi);
0157
0158 std::cout << "Done initing the treemaker"<<std::endl;
0159 return Fun4AllReturnCodes::EVENT_OK;
0160 }
0161
0162
0163 int DijetTreeMaker::InitRun(PHCompositeNode *topNode)
0164 {
0165 if (Verbosity()) std::cout << __FUNCTION__ << __LINE__<<std::endl;
0166 return Fun4AllReturnCodes::EVENT_OK;
0167 }
0168
0169
0170
0171 void DijetTreeMaker::SetVerbosity(int verbo){
0172 _verbosity = verbo;
0173 return;
0174 }
0175
0176 void DijetTreeMaker::reset_tree_vars()
0177 {
0178 if (Verbosity()) std::cout << __FUNCTION__ << __LINE__<<std::endl;
0179
0180 b_ncluster = 0;
0181 b_cluster_e.clear();
0182 b_cluster_ecore.clear();
0183 b_cluster_eta.clear();
0184 b_cluster_phi.clear();
0185 b_cluster_pt.clear();
0186 b_cluster_chi2.clear();
0187 b_cluster_prob.clear();
0188
0189 b_njet_2 = 0;
0190 b_jet_pt_2.clear();
0191 b_jet_et_2.clear();
0192 b_jet_eta_2.clear();
0193 b_jet_phi_2.clear();
0194 b_jet_emcal_2.clear();
0195 b_jet_hcalin_2.clear();
0196 b_jet_hcalout_2.clear();
0197 b_jet_eccen_2.clear();
0198
0199 b_njet_4 = 0;
0200 b_jet_pt_4.clear();
0201 b_jet_et_4.clear();
0202 b_jet_eta_4.clear();
0203 b_jet_phi_4.clear();
0204 b_jet_emcal_4.clear();
0205 b_jet_hcalin_4.clear();
0206 b_jet_hcalout_4.clear();
0207 b_jet_eccen_4.clear();
0208
0209 b_njet_6 = 0;
0210 b_jet_pt_6.clear();
0211 b_jet_et_6.clear();
0212 b_jet_eta_6.clear();
0213 b_jet_phi_6.clear();
0214 b_jet_emcal_6.clear();
0215 b_jet_hcalin_6.clear();
0216 b_jet_hcalout_6.clear();
0217 b_jet_eccen_4.clear();
0218
0219 b_emcal_good.clear();
0220 b_emcal_energy.clear();
0221 b_emcal_time.clear();
0222 b_emcal_phibin.clear();
0223 b_emcal_etabin.clear();
0224 b_emcal_phi.clear();
0225 b_emcal_eta.clear();
0226
0227 b_hcalin_good.clear();
0228 b_hcalin_energy.clear();
0229 b_hcalin_time.clear();
0230 b_hcalin_phibin.clear();
0231 b_hcalin_etabin.clear();
0232 b_hcalin_phi.clear();
0233 b_hcalin_eta.clear();
0234
0235 b_hcalout_good.clear();
0236 b_hcalout_energy.clear();
0237 b_hcalout_time.clear();
0238 b_hcalout_phibin.clear();
0239 b_hcalout_etabin.clear();
0240 b_hcalout_phi.clear();
0241 b_hcalout_eta.clear();
0242
0243 if (isSim)
0244 {
0245 b_ntruth_jet_2 = 0;
0246 b_truth_jet_pt_2.clear();
0247 b_truth_jet_eta_2.clear();
0248 b_truth_jet_phi_2.clear();
0249 b_ntruth_jet_4 = 0;
0250 b_truth_jet_pt_4.clear();
0251 b_truth_jet_eta_4.clear();
0252 b_truth_jet_phi_4.clear();
0253 b_ntruth_jet_6 = 0;
0254 b_truth_jet_pt_6.clear();
0255 b_truth_jet_eta_6.clear();
0256 b_truth_jet_phi_6.clear();
0257 }
0258
0259 return;
0260 }
0261
0262 int DijetTreeMaker::process_event(PHCompositeNode *topNode)
0263 {
0264
0265 if (Verbosity()) std::cout << __FILE__ << " "<< __LINE__<<" "<<std::endl;
0266
0267 int dijet_candidate = false;
0268 _i_event++;
0269
0270 reset_tree_vars();
0271
0272 RawClusterContainer *clusters = findNode::getClass<RawClusterContainer>(topNode, "CLUSTERINFO_CEMC");
0273
0274 RawTowerGeomContainer *tower_geomIH = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALIN");
0275 RawTowerGeomContainer *tower_geomEM = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_CEMC");
0276 RawTowerGeomContainer *tower_geomOH = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALOUT");
0277
0278 float max_secondmax_et[2]={0};
0279 float max_secondmax_pt[2]={0};
0280
0281 if (isSim)
0282 {
0283
0284 std::string truthJetName = "AntiKt_Truth_r02";
0285 JetContainer *jetstruth_2 = findNode::getClass<JetContainer>(topNode, truthJetName);
0286 if (jetstruth_2)
0287 {
0288
0289
0290 for (auto jet : *jetstruth_2)
0291 {
0292 if (jet->get_pt() < pt_cut_truth) continue;
0293
0294 b_ntruth_jet_2++;
0295 b_truth_jet_pt_2.push_back(jet->get_pt());
0296 b_truth_jet_eta_2.push_back(jet->get_eta());
0297 b_truth_jet_phi_2.push_back(jet->get_phi());
0298 for (int im = 0; im < 2; im++)
0299 {
0300 if (jet->get_pt() > max_secondmax_pt[im])
0301 {
0302 if (im == 0) max_secondmax_pt[1] = max_secondmax_pt[0];
0303 max_secondmax_pt[im] = jet->get_pt();
0304 }
0305 }
0306
0307 }
0308
0309 if (max_secondmax_pt[0] > 10 && max_secondmax_pt[1] > 5) dijet_candidate = true;
0310 for (int im = 0; im < 2; im++)
0311 {
0312 max_secondmax_pt[im]=0;
0313 }
0314 }
0315 truthJetName = "AntiKt_Truth_r04";
0316 JetContainer *jetstruth_4 = findNode::getClass<JetContainer>(topNode, truthJetName);
0317 if (jetstruth_4)
0318 {
0319
0320
0321 for (auto jet : *jetstruth_4)
0322 {
0323 if (jet->get_pt() < pt_cut_truth) continue;
0324
0325 b_ntruth_jet_4++;
0326 b_truth_jet_pt_4.push_back(jet->get_pt());
0327 b_truth_jet_eta_4.push_back(jet->get_eta());
0328 b_truth_jet_phi_4.push_back(jet->get_phi());
0329 for (int im = 0; im < 2; im++)
0330 {
0331 if (jet->get_pt() > max_secondmax_pt[im])
0332 {
0333 if (im == 0) max_secondmax_pt[1] = max_secondmax_pt[0];
0334 max_secondmax_pt[im] = jet->get_pt();
0335 }
0336 }
0337
0338 }
0339
0340 if (max_secondmax_pt[0] > 10 && max_secondmax_pt[1] > 5) dijet_candidate = true;
0341 for (int im = 0; im < 2; im++)
0342 {
0343 max_secondmax_pt[im]=0;
0344 }
0345 }
0346 truthJetName = "AntiKt_Truth_r06";
0347 JetContainer *jetstruth_6 = findNode::getClass<JetContainer>(topNode, truthJetName);
0348 if (jetstruth_6)
0349 {
0350
0351
0352 for (auto jet : *jetstruth_6)
0353 {
0354 if (jet->get_pt() < pt_cut_truth) continue;
0355
0356 b_ntruth_jet_6++;
0357 b_truth_jet_pt_6.push_back(jet->get_pt());
0358 b_truth_jet_eta_6.push_back(jet->get_eta());
0359 b_truth_jet_phi_6.push_back(jet->get_phi());
0360 for (int im = 0; im < 2; im++)
0361 {
0362 if (jet->get_pt() > max_secondmax_pt[im])
0363 {
0364 if (im == 0) max_secondmax_pt[1] = max_secondmax_pt[0];
0365 max_secondmax_pt[im] = jet->get_pt();
0366 }
0367 }
0368
0369 }
0370
0371 if (max_secondmax_pt[0] > 10 && max_secondmax_pt[1] > 5) dijet_candidate = true;
0372 for (int im = 0; im < 2; im++)
0373 {
0374 max_secondmax_pt[im]=0;
0375 }
0376 }
0377
0378 }
0379 GlobalVertexMap* vertexmap = findNode::getClass<GlobalVertexMap>(topNode, "RealVertexMap");
0380 if (!vertexmap)
0381 {
0382 vertexmap = findNode::getClass<GlobalVertexMap>(topNode, "GlobalVertexMap");
0383 }
0384
0385 float vtx_z = 0;
0386 b_vertex_z = -999;
0387 if (vertexmap && !vertexmap->empty())
0388 {
0389 GlobalVertex* vtx = vertexmap->begin()->second;
0390 if (vtx)
0391 {
0392 vtx_z = vtx->get_z();
0393 b_vertex_z = vtx_z;
0394 }
0395 }
0396 if (Verbosity() > 5) std::cout << "Mbd Vertex = " << b_vertex_z << std::endl;
0397
0398 TowerInfoContainer *hcalin_towers = findNode::getClass<TowerInfoContainer>(topNode, m_calo_nodename + "_HCALIN");
0399 if (!hcalin_towers)
0400 {
0401 if (Verbosity()) std::cout << "no hcalin towers "<<std::endl;
0402 return Fun4AllReturnCodes::ABORTRUN;
0403 }
0404
0405 TowerInfoContainer *hcalout_towers = findNode::getClass<TowerInfoContainer>(topNode, m_calo_nodename + "_HCALOUT");
0406 if (!hcalout_towers)
0407 {
0408 std::cout << "no hcalout towers "<<std::endl;
0409 return Fun4AllReturnCodes::ABORTRUN;
0410 }
0411 TowerInfoContainer *emcalre_towers = findNode::getClass<TowerInfoContainer>(topNode, m_calo_nodename + "_CEMC_RETOWER");
0412 TowerInfoContainer *emcal_towers = findNode::getClass<TowerInfoContainer>(topNode, m_calo_nodename + "_CEMC");
0413
0414 float background_v2 = 0;
0415 float background_Psi2 = 0;
0416 bool has_tower_background = false;
0417 TowerBackground *towBack = findNode::getClass<TowerBackground>(topNode, "TowerInfoBackground_Sub2");
0418 if (towBack)
0419 {
0420 has_tower_background = true;
0421 background_v2 = towBack->get_v2();
0422 background_Psi2 = towBack->get_Psi2();
0423 }
0424
0425 std::string recoJetName = "AntiKt_Tower_r02_Sub1";
0426
0427
0428 JetContainer *jets_2 = findNode::getClass<JetContainer>(topNode, recoJetName);
0429 if (jets_2)
0430 {
0431
0432
0433 for (auto jet : *jets_2)
0434 {
0435 if (jet->get_pt() < pt_cut) continue;
0436
0437 int n_comp_total = 0;
0438 int n_comp_ihcal = 0;
0439 int n_comp_ohcal = 0;
0440 int n_comp_emcal = 0;
0441
0442 float jet_total_eT = 0;
0443 float eFrac_ihcal = 0;
0444 float eFrac_ohcal = 0;
0445 float eFrac_emcal = 0;
0446
0447 b_njet_2++;
0448 b_jet_pt_2.push_back(jet->get_pt());
0449 b_jet_eta_2.push_back(jet->get_eta());
0450 b_jet_phi_2.push_back(jet->get_phi());
0451
0452 std::vector<TVector3> hcal_constituents;
0453
0454 for (auto comp : jet->get_comp_vec())
0455 {
0456
0457 unsigned int channel = comp.second;
0458 TowerInfo *tower;
0459 float tower_eT = 0;
0460 if (comp.first == 26 || comp.first == 30)
0461 {
0462 tower = hcalin_towers->get_tower_at_channel(channel);
0463 if (!tower || !tower_geomIH)
0464 {
0465 continue;
0466 }
0467 if(!tower->get_isGood()) continue;
0468
0469 unsigned int calokey = hcalin_towers->encode_key(channel);
0470
0471 int ieta = hcalin_towers->getTowerEtaBin(calokey);
0472 int iphi = hcalin_towers->getTowerPhiBin(calokey);
0473 const RawTowerDefs::keytype key = RawTowerDefs::encode_towerid(RawTowerDefs::CalorimeterId::HCALIN, ieta, iphi);
0474 float tower_phi = tower_geomIH->get_tower_geometry(key)->get_phi();
0475 float tower_eta = tower_geomIH->get_tower_geometry(key)->get_eta();
0476 tower_eT = tower->get_energy() / std::cosh(tower_eta);
0477
0478 if (comp.first == 30)
0479 {
0480 if (has_tower_background)
0481 {
0482 float UE = towBack->get_UE(1).at(ieta);
0483 float tower_UE = UE * (1 + 2 * background_v2 * std::cos(2 * (tower_phi - background_Psi2)));
0484 tower_eT = (tower->get_energy() - tower_UE) / std::cosh(tower_eta);
0485 }
0486 }
0487
0488 eFrac_ihcal += tower_eT;
0489 jet_total_eT += tower_eT;
0490 n_comp_ihcal++;
0491 n_comp_total++;
0492 }
0493 else if (comp.first == 27 || comp.first == 31)
0494 {
0495 tower = hcalout_towers->get_tower_at_channel(channel);
0496
0497 if (!tower || !tower_geomOH)
0498 {
0499 continue;
0500 }
0501
0502 unsigned int calokey = hcalout_towers->encode_key(channel);
0503 int ieta = hcalout_towers->getTowerEtaBin(calokey);
0504 int iphi = hcalout_towers->getTowerPhiBin(calokey);
0505 const RawTowerDefs::keytype key = RawTowerDefs::encode_towerid(RawTowerDefs::CalorimeterId::HCALOUT, ieta, iphi);
0506 float tower_phi = tower_geomOH->get_tower_geometry(key)->get_phi();
0507 float tower_eta = tower_geomOH->get_tower_geometry(key)->get_eta();
0508 tower_eT = tower->get_energy() / std::cosh(tower_eta);
0509
0510 if (comp.first == 31)
0511 {
0512 if (has_tower_background)
0513 {
0514 float UE = towBack->get_UE(2).at(ieta);
0515 float tower_UE = UE * (1 + 2 * background_v2 * std::cos(2 * (tower_phi - background_Psi2)));
0516 tower_eT = (tower->get_energy() - tower_UE) / std::cosh(tower_eta);
0517 }
0518 }
0519 TVector3 v;
0520 v.SetPtEtaPhi(tower_eT, tower_eta, tower_phi);
0521 hcal_constituents.push_back(v);
0522 eFrac_ohcal += tower_eT;
0523 jet_total_eT += tower_eT;
0524 n_comp_ohcal++;
0525 n_comp_total++;
0526 }
0527 else if (comp.first == 28 || comp.first == 29)
0528 {
0529 tower = emcalre_towers->get_tower_at_channel(channel);
0530
0531 if (!tower || !tower_geomIH)
0532 {
0533 continue;
0534 }
0535
0536 unsigned int calokey = emcalre_towers->encode_key(channel);
0537 int ieta = emcalre_towers->getTowerEtaBin(calokey);
0538 int iphi = emcalre_towers->getTowerPhiBin(calokey);
0539 const RawTowerDefs::keytype key = RawTowerDefs::encode_towerid(RawTowerDefs::CalorimeterId::HCALIN, ieta, iphi);
0540 float tower_phi = tower_geomIH->get_tower_geometry(key)->get_phi();
0541 float tower_eta = tower_geomIH->get_tower_geometry(key)->get_eta();
0542 tower_eT = tower->get_energy() / std::cosh(tower_eta);
0543
0544 if (comp.first == 29)
0545 {
0546 if (has_tower_background)
0547 {
0548 float UE = towBack->get_UE(0).at(ieta);
0549 float tower_UE = UE * (1 + 2 * background_v2 * std::cos(2 * (tower_phi - background_Psi2)));
0550 tower_eT = (tower->get_energy() - tower_UE) / std::cosh(tower_eta);
0551 }
0552 }
0553
0554 eFrac_emcal += tower_eT;
0555 jet_total_eT += tower_eT;
0556 n_comp_emcal++;
0557 n_comp_total++;
0558 }
0559 }
0560
0561 eFrac_emcal /= jet_total_eT;
0562 eFrac_ihcal /= jet_total_eT;
0563
0564 float eccen = get_eccentricity(hcal_constituents, eFrac_ohcal);
0565 eFrac_ohcal /= jet_total_eT;
0566
0567 b_jet_et_2.push_back(jet_total_eT);
0568 b_jet_emcal_2.push_back(eFrac_emcal);
0569 b_jet_hcalin_2.push_back(eFrac_ihcal);
0570 b_jet_hcalout_2.push_back(eFrac_ohcal);
0571 b_jet_eccen_2.push_back(eccen);
0572 for (int im = 0; im < 2; im++)
0573 {
0574 if (jet_total_eT > max_secondmax_et[im])
0575 {
0576 if (im == 0) max_secondmax_et[1] = max_secondmax_et[0];
0577 max_secondmax_et[im] = jet_total_eT;
0578 }
0579 }
0580 for (int im = 0; im < 2; im++)
0581 {
0582 if (jet->get_pt() > max_secondmax_pt[im])
0583 {
0584 if (im == 0) max_secondmax_pt[1] = max_secondmax_pt[0];
0585 max_secondmax_pt[im] = jet->get_pt();
0586 }
0587 }
0588
0589 }
0590 }
0591 dijet_candidate = false;
0592 if (max_secondmax_et[0] > dijetcut) dijet_candidate = true;
0593 if (max_secondmax_pt[0] > dijetcut) dijet_candidate = true;
0594 for (int im = 0; im < 2; im++)
0595 {
0596 max_secondmax_et[im]=0;
0597 max_secondmax_pt[im]=0;
0598 }
0599 recoJetName = "AntiKt_Tower_r04_Sub1";
0600 JetContainer *jets_4 = findNode::getClass<JetContainer>(topNode, recoJetName);
0601 if (jets_4)
0602 {
0603
0604
0605 for (auto jet : *jets_4)
0606 {
0607 if (jet->get_pt() < pt_cut) continue;
0608
0609 int n_comp_total = 0;
0610 int n_comp_ihcal = 0;
0611 int n_comp_ohcal = 0;
0612 int n_comp_emcal = 0;
0613
0614 float jet_total_eT = 0;
0615 float eFrac_ihcal = 0;
0616 float eFrac_ohcal = 0;
0617 float eFrac_emcal = 0;
0618
0619 b_njet_4++;
0620 b_jet_pt_4.push_back(jet->get_pt());
0621 b_jet_eta_4.push_back(jet->get_eta());
0622 b_jet_phi_4.push_back(jet->get_phi());
0623
0624 std::vector<TVector3> hcal_constituents;
0625 for (auto comp : jet->get_comp_vec())
0626 {
0627
0628 unsigned int channel = comp.second;
0629 TowerInfo *tower;
0630 float tower_eT = 0;
0631 if (comp.first == 26 || comp.first == 30)
0632 {
0633 tower = hcalin_towers->get_tower_at_channel(channel);
0634 if (!tower || !tower_geomIH)
0635 {
0636 continue;
0637 }
0638 if(!tower->get_isGood()) continue;
0639
0640 unsigned int calokey = hcalin_towers->encode_key(channel);
0641
0642 int ieta = hcalin_towers->getTowerEtaBin(calokey);
0643 int iphi = hcalin_towers->getTowerPhiBin(calokey);
0644 const RawTowerDefs::keytype key = RawTowerDefs::encode_towerid(RawTowerDefs::CalorimeterId::HCALIN, ieta, iphi);
0645 float tower_phi = tower_geomIH->get_tower_geometry(key)->get_phi();
0646 float tower_eta = tower_geomIH->get_tower_geometry(key)->get_eta();
0647 tower_eT = tower->get_energy() / std::cosh(tower_eta);
0648
0649 if (comp.first == 30)
0650 {
0651 if (has_tower_background)
0652 {
0653 float UE = towBack->get_UE(1).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_ihcal += tower_eT;
0660 jet_total_eT += tower_eT;
0661 n_comp_ihcal++;
0662 n_comp_total++;
0663 }
0664 else if (comp.first == 27 || comp.first == 31)
0665 {
0666 tower = hcalout_towers->get_tower_at_channel(channel);
0667
0668 if (!tower || !tower_geomOH)
0669 {
0670 continue;
0671 }
0672
0673 unsigned int calokey = hcalout_towers->encode_key(channel);
0674 int ieta = hcalout_towers->getTowerEtaBin(calokey);
0675 int iphi = hcalout_towers->getTowerPhiBin(calokey);
0676 const RawTowerDefs::keytype key = RawTowerDefs::encode_towerid(RawTowerDefs::CalorimeterId::HCALOUT, ieta, iphi);
0677 float tower_phi = tower_geomOH->get_tower_geometry(key)->get_phi();
0678 float tower_eta = tower_geomOH->get_tower_geometry(key)->get_eta();
0679 tower_eT = tower->get_energy() / std::cosh(tower_eta);
0680
0681 if (comp.first == 31)
0682 {
0683 if (has_tower_background)
0684 {
0685 float UE = towBack->get_UE(2).at(ieta);
0686 float tower_UE = UE * (1 + 2 * background_v2 * std::cos(2 * (tower_phi - background_Psi2)));
0687 tower_eT = (tower->get_energy() - tower_UE) / std::cosh(tower_eta);
0688 }
0689 }
0690 TVector3 v;
0691 v.SetPtEtaPhi(tower_eT, tower_eta, tower_phi);
0692 hcal_constituents.push_back(v);
0693
0694 eFrac_ohcal += tower_eT;
0695 jet_total_eT += tower_eT;
0696 n_comp_ohcal++;
0697 n_comp_total++;
0698 }
0699 else if (comp.first == 28 || comp.first == 29)
0700 {
0701 tower = emcalre_towers->get_tower_at_channel(channel);
0702
0703 if (!tower || !tower_geomIH)
0704 {
0705 continue;
0706 }
0707
0708 unsigned int calokey = emcalre_towers->encode_key(channel);
0709 int ieta = emcalre_towers->getTowerEtaBin(calokey);
0710 int iphi = emcalre_towers->getTowerPhiBin(calokey);
0711 const RawTowerDefs::keytype key = RawTowerDefs::encode_towerid(RawTowerDefs::CalorimeterId::HCALIN, ieta, iphi);
0712 float tower_phi = tower_geomIH->get_tower_geometry(key)->get_phi();
0713 float tower_eta = tower_geomIH->get_tower_geometry(key)->get_eta();
0714 tower_eT = tower->get_energy() / std::cosh(tower_eta);
0715
0716 if (comp.first == 29)
0717 {
0718 if (has_tower_background)
0719 {
0720 float UE = towBack->get_UE(0).at(ieta);
0721 float tower_UE = UE * (1 + 2 * background_v2 * std::cos(2 * (tower_phi - background_Psi2)));
0722 tower_eT = (tower->get_energy() - tower_UE) / std::cosh(tower_eta);
0723 }
0724 }
0725
0726 eFrac_emcal += tower_eT;
0727 jet_total_eT += tower_eT;
0728 n_comp_emcal++;
0729 n_comp_total++;
0730 }
0731 }
0732 float eccen = get_eccentricity(hcal_constituents, eFrac_ohcal);
0733 eFrac_ohcal /= jet_total_eT;
0734 eFrac_emcal /= jet_total_eT;
0735 eFrac_ihcal /= jet_total_eT;
0736
0737 b_jet_et_4.push_back(jet_total_eT);
0738 b_jet_emcal_4.push_back(eFrac_emcal);
0739 b_jet_hcalin_4.push_back(eFrac_ihcal);
0740 b_jet_hcalout_4.push_back(eFrac_ohcal);
0741
0742 b_jet_eccen_4.push_back(eccen);
0743 for (int im = 0; im < 2; im++)
0744 {
0745 if (jet_total_eT > max_secondmax_et[im])
0746 {
0747 if (im == 0) max_secondmax_et[1] = max_secondmax_et[0];
0748 max_secondmax_et[im] = jet_total_eT;
0749 }
0750 }
0751 for (int im = 0; im < 2; im++)
0752 {
0753 if (jet->get_pt() > max_secondmax_pt[im])
0754 {
0755 if (im == 0) max_secondmax_pt[1] = max_secondmax_pt[0];
0756 max_secondmax_pt[im] = jet->get_pt();
0757 }
0758 }
0759
0760 }
0761 }
0762
0763 if (max_secondmax_et[0] > dijetcut) dijet_candidate = true;
0764 if (max_secondmax_pt[0] > dijetcut) dijet_candidate = true;
0765 for (int im = 0; im < 2; im++)
0766 {
0767 max_secondmax_et[im]=0;
0768 max_secondmax_pt[im]=0;
0769 }
0770
0771
0772 recoJetName = "AntiKt_Tower_r06_Sub1";
0773 JetContainer *jets_6 = findNode::getClass<JetContainer>(topNode, recoJetName);
0774 if (jets_6)
0775 {
0776
0777
0778 for (auto jet : *jets_6)
0779 {
0780 if (jet->get_pt() < pt_cut) continue;
0781
0782 int n_comp_total = 0;
0783 int n_comp_ihcal = 0;
0784 int n_comp_ohcal = 0;
0785 int n_comp_emcal = 0;
0786
0787 float jet_total_eT = 0;
0788 float eFrac_ihcal = 0;
0789 float eFrac_ohcal = 0;
0790 float eFrac_emcal = 0;
0791
0792 b_njet_6++;
0793 b_jet_pt_6.push_back(jet->get_pt());
0794 b_jet_eta_6.push_back(jet->get_eta());
0795 b_jet_phi_6.push_back(jet->get_phi());
0796
0797 std::vector<TVector3> hcal_constituents;
0798 for (auto comp : jet->get_comp_vec())
0799 {
0800
0801 unsigned int channel = comp.second;
0802 TowerInfo *tower;
0803 float tower_eT = 0;
0804 if (comp.first == 26 || comp.first == 30)
0805 {
0806 tower = hcalin_towers->get_tower_at_channel(channel);
0807 if (!tower || !tower_geomIH)
0808 {
0809 continue;
0810 }
0811 if(!tower->get_isGood()) continue;
0812
0813 unsigned int calokey = hcalin_towers->encode_key(channel);
0814
0815 int ieta = hcalin_towers->getTowerEtaBin(calokey);
0816 int iphi = hcalin_towers->getTowerPhiBin(calokey);
0817 const RawTowerDefs::keytype key = RawTowerDefs::encode_towerid(RawTowerDefs::CalorimeterId::HCALIN, ieta, iphi);
0818 float tower_phi = tower_geomIH->get_tower_geometry(key)->get_phi();
0819 float tower_eta = tower_geomIH->get_tower_geometry(key)->get_eta();
0820 tower_eT = tower->get_energy() / std::cosh(tower_eta);
0821
0822 if (comp.first == 30)
0823 {
0824 if (has_tower_background)
0825 {
0826 float UE = towBack->get_UE(1).at(ieta);
0827 float tower_UE = UE * (1 + 2 * background_v2 * std::cos(2 * (tower_phi - background_Psi2)));
0828 tower_eT = (tower->get_energy() - tower_UE) / std::cosh(tower_eta);
0829 }
0830 }
0831
0832 eFrac_ihcal += tower_eT;
0833 jet_total_eT += tower_eT;
0834 n_comp_ihcal++;
0835 n_comp_total++;
0836 }
0837 else if (comp.first == 27 || comp.first == 31)
0838 {
0839 tower = hcalout_towers->get_tower_at_channel(channel);
0840
0841 if (!tower || !tower_geomOH)
0842 {
0843 continue;
0844 }
0845
0846 unsigned int calokey = hcalout_towers->encode_key(channel);
0847 int ieta = hcalout_towers->getTowerEtaBin(calokey);
0848 int iphi = hcalout_towers->getTowerPhiBin(calokey);
0849 const RawTowerDefs::keytype key = RawTowerDefs::encode_towerid(RawTowerDefs::CalorimeterId::HCALOUT, ieta, iphi);
0850 float tower_phi = tower_geomOH->get_tower_geometry(key)->get_phi();
0851 float tower_eta = tower_geomOH->get_tower_geometry(key)->get_eta();
0852 tower_eT = tower->get_energy() / std::cosh(tower_eta);
0853
0854 if (comp.first == 31)
0855 {
0856 if (has_tower_background)
0857 {
0858 float UE = towBack->get_UE(2).at(ieta);
0859 float tower_UE = UE * (1 + 2 * background_v2 * std::cos(2 * (tower_phi - background_Psi2)));
0860 tower_eT = (tower->get_energy() - tower_UE) / std::cosh(tower_eta);
0861 }
0862 }
0863 TVector3 v;
0864 v.SetPtEtaPhi(tower_eT, tower_eta, tower_phi);
0865 hcal_constituents.push_back(v);
0866
0867 eFrac_ohcal += tower_eT;
0868 jet_total_eT += tower_eT;
0869 n_comp_ohcal++;
0870 n_comp_total++;
0871 }
0872 else if (comp.first == 28 || comp.first == 29)
0873 {
0874 tower = emcalre_towers->get_tower_at_channel(channel);
0875
0876 if (!tower || !tower_geomIH)
0877 {
0878 continue;
0879 }
0880
0881 unsigned int calokey = emcalre_towers->encode_key(channel);
0882 int ieta = emcalre_towers->getTowerEtaBin(calokey);
0883 int iphi = emcalre_towers->getTowerPhiBin(calokey);
0884 const RawTowerDefs::keytype key = RawTowerDefs::encode_towerid(RawTowerDefs::CalorimeterId::HCALIN, ieta, iphi);
0885 float tower_phi = tower_geomIH->get_tower_geometry(key)->get_phi();
0886 float tower_eta = tower_geomIH->get_tower_geometry(key)->get_eta();
0887 tower_eT = tower->get_energy() / std::cosh(tower_eta);
0888
0889 if (comp.first == 29)
0890 {
0891 if (has_tower_background)
0892 {
0893 float UE = towBack->get_UE(0).at(ieta);
0894 float tower_UE = UE * (1 + 2 * background_v2 * std::cos(2 * (tower_phi - background_Psi2)));
0895 tower_eT = (tower->get_energy() - tower_UE) / std::cosh(tower_eta);
0896 }
0897 }
0898
0899 eFrac_emcal += tower_eT;
0900 jet_total_eT += tower_eT;
0901 n_comp_emcal++;
0902 n_comp_total++;
0903 }
0904 }
0905 float eccen = get_eccentricity(hcal_constituents, eFrac_ohcal);
0906 eFrac_ohcal /= jet_total_eT;
0907 eFrac_emcal /= jet_total_eT;
0908 eFrac_ihcal /= jet_total_eT;
0909
0910 b_jet_eccen_6.push_back(eccen);
0911
0912 b_jet_et_6.push_back(jet_total_eT);
0913 b_jet_emcal_6.push_back(eFrac_emcal);
0914 b_jet_hcalin_6.push_back(eFrac_ihcal);
0915 b_jet_hcalout_6.push_back(eFrac_ohcal);
0916 for (int im = 0; im < 2; im++)
0917 {
0918 if (jet_total_eT > max_secondmax_et[im])
0919 {
0920 if (im == 0) max_secondmax_et[1] = max_secondmax_et[0];
0921 max_secondmax_et[im] = jet_total_eT;
0922 }
0923 }
0924 for (int im = 0; im < 2; im++)
0925 {
0926 if (jet->get_pt() > max_secondmax_pt[im])
0927 {
0928 if (im == 0) max_secondmax_pt[1] = max_secondmax_pt[0];
0929 max_secondmax_pt[im] = jet->get_pt();
0930 }
0931 }
0932
0933 }
0934 }
0935
0936 if (max_secondmax_et[0] > dijetcut) dijet_candidate = true;
0937 if (max_secondmax_pt[0] > dijetcut) dijet_candidate = true;
0938
0939 for (int im = 0; im < 2; im++)
0940 {
0941 max_secondmax_et[im]=0;
0942 max_secondmax_pt[im]=0;
0943 }
0944
0945 if (!dijet_candidate)
0946 {
0947 return Fun4AllReturnCodes::EVENT_OK;
0948 }
0949
0950 int size;
0951 if (emcal_towers)
0952 {
0953 size = emcal_towers->size();
0954 for (int channel = 0; channel < size;channel++)
0955 {
0956 _tower = emcal_towers->get_tower_at_channel(channel);
0957 float energy = _tower->get_energy();
0958 float time = _tower->get_time();
0959 unsigned int towerkey = emcal_towers->encode_key(channel);
0960 int ieta = emcal_towers->getTowerEtaBin(towerkey);
0961 int iphi = emcal_towers->getTowerPhiBin(towerkey);
0962 short good = (_tower->get_isGood() ? 1:0);
0963 const RawTowerDefs::keytype key = RawTowerDefs::encode_towerid(RawTowerDefs::CalorimeterId::CEMC, ieta, iphi);
0964 float tower_phi = tower_geomEM->get_tower_geometry(key)->get_phi();
0965 float tower_eta = tower_geomEM->get_tower_geometry(key)->get_eta();
0966
0967 b_emcal_good.push_back(good);
0968 b_emcal_energy.push_back(energy);
0969 b_emcal_time.push_back(time);
0970 b_emcal_etabin.push_back(ieta);
0971 b_emcal_phibin.push_back(iphi);
0972 b_emcal_eta.push_back(tower_eta);
0973 b_emcal_phi.push_back(tower_phi);
0974 }
0975 }
0976
0977 if (hcalin_towers)
0978 {
0979
0980 size = hcalin_towers->size();
0981 for (int channel = 0; channel < size;channel++)
0982 {
0983 _tower = hcalin_towers->get_tower_at_channel(channel);
0984 float energy = _tower->get_energy();
0985 float time = _tower->get_time();
0986 short good = (_tower->get_isGood() ? 1:0);
0987 unsigned int towerkey = hcalin_towers->encode_key(channel);
0988 int ieta = hcalin_towers->getTowerEtaBin(towerkey);
0989 int iphi = hcalin_towers->getTowerPhiBin(towerkey);
0990 const RawTowerDefs::keytype key = RawTowerDefs::encode_towerid(RawTowerDefs::CalorimeterId::HCALIN, ieta, iphi);
0991 float tower_phi = tower_geomIH->get_tower_geometry(key)->get_phi();
0992 float tower_eta = tower_geomIH->get_tower_geometry(key)->get_eta();
0993
0994 b_hcalin_good.push_back(good);
0995 b_hcalin_energy.push_back(energy);
0996 b_hcalin_time.push_back(time);
0997 b_hcalin_etabin.push_back(ieta);
0998 b_hcalin_phibin.push_back(iphi);
0999 b_hcalin_eta.push_back(tower_eta);
1000 b_hcalin_phi.push_back(tower_phi);
1001 }
1002 }
1003 if (hcalout_towers)
1004 {
1005
1006 size = hcalout_towers->size();
1007 for (int channel = 0; channel < size;channel++)
1008 {
1009 _tower = hcalout_towers->get_tower_at_channel(channel);
1010 float energy = _tower->get_energy();
1011 float time = _tower->get_time();
1012 unsigned int towerkey = hcalout_towers->encode_key(channel);
1013 int ieta = hcalout_towers->getTowerEtaBin(towerkey);
1014 int iphi = hcalout_towers->getTowerPhiBin(towerkey);
1015 short good = (_tower->get_isGood() ? 1:0);
1016 const RawTowerDefs::keytype key = RawTowerDefs::encode_towerid(RawTowerDefs::CalorimeterId::HCALOUT, ieta, iphi);
1017 float tower_phi = tower_geomOH->get_tower_geometry(key)->get_phi();
1018 float tower_eta = tower_geomOH->get_tower_geometry(key)->get_eta();
1019
1020 b_hcalout_good.push_back(good);
1021 b_hcalout_energy.push_back(energy);
1022 b_hcalout_time.push_back(time);
1023 b_hcalout_etabin.push_back(ieta);
1024 b_hcalout_phibin.push_back(iphi);
1025 b_hcalout_eta.push_back(tower_eta);
1026 b_hcalout_phi.push_back(tower_phi);
1027 }
1028 }
1029
1030
1031 if (clusters)
1032 {
1033 CLHEP::Hep3Vector vertex(0, 0, b_vertex_z);
1034
1035 auto clusterrange = clusters->getClusters();
1036 for (auto iclus = clusterrange.first; iclus != clusterrange.second; ++iclus)
1037 {
1038
1039 RawCluster *cluster = iclus->second;
1040 CLHEP::Hep3Vector E_vec_cluster = RawClusterUtility::GetECoreVec(*cluster, vertex);
1041
1042 float energy = E_vec_cluster.mag();
1043 float pt = E_vec_cluster.perp();
1044 float clus_eta = E_vec_cluster.pseudoRapidity();
1045 float clus_phi = E_vec_cluster.phi();
1046 float chi2 = cluster->get_chi2();
1047 float prob = cluster->get_prob();
1048 float ecore = cluster->get_ecore();
1049 if (energy < 0.5) continue;
1050 b_ncluster++;
1051 b_cluster_e.push_back(energy);
1052 b_cluster_ecore.push_back(ecore);
1053 b_cluster_eta.push_back(clus_eta);
1054 b_cluster_phi.push_back(clus_phi);
1055 b_cluster_pt.push_back(pt);
1056 b_cluster_chi2.push_back(chi2);
1057 b_cluster_prob.push_back(prob);
1058 }
1059 }
1060
1061
1062
1063
1064
1065
1066
1067
1068
1069
1070
1071
1072
1073
1074
1075
1076
1077
1078
1079
1080
1081
1082
1083
1084
1085
1086
1087
1088
1089
1090
1091
1092
1093
1094
1095
1096
1097
1098
1099
1100
1101
1102
1103
1104
1105
1106
1107
1108
1109
1110
1111
1112
1113
1114
1115
1116
1117
1118
1119
1120
1121
1122
1123
1124
1125
1126
1127
1128
1129
1130
1131
1132
1133
1134
1135
1136
1137
1138
1139
1140
1141
1142
1143
1144
1145
1146
1147
1148
1149
1150
1151
1152
1153
1154
1155
1156
1157
1158
1159
1160
1161
1162
1163
1164
1165
1166
1167
1168
1169
1170
1171
1172
1173
1174
1175
1176
1177
1178
1179
1180
1181
1182
1183
1184
1185
1186
1187
1188
1189
1190
1191
1192
1193
1194
1195
1196
1197
1198
1199
1200
1201
1202
1203
1204
1205
1206 if (Verbosity()) std::cout << __FILE__ << " "<< __LINE__<<" "<<std::endl;
1207
1208 _tree->Fill();
1209
1210 return Fun4AllReturnCodes::EVENT_OK;
1211 }
1212
1213
1214
1215 void DijetTreeMaker::GetNodes(PHCompositeNode* topNode)
1216 {
1217
1218
1219 }
1220
1221 int DijetTreeMaker::ResetEvent(PHCompositeNode *topNode)
1222 {
1223 if (Verbosity() > 0)
1224 {
1225 std::cout << "DijetTreeMaker::ResetEvent(PHCompositeNode *topNode) Resetting internal structures, prepare for next event" << std::endl;
1226 }
1227
1228
1229 return Fun4AllReturnCodes::EVENT_OK;
1230 }
1231
1232
1233 int DijetTreeMaker::EndRun(const int runnumber)
1234 {
1235 if (Verbosity() > 0)
1236 {
1237 std::cout << "DijetTreeMaker::EndRun(const int runnumber) Ending Run for Run " << runnumber << std::endl;
1238 }
1239 return Fun4AllReturnCodes::EVENT_OK;
1240 }
1241
1242
1243 int DijetTreeMaker::End(PHCompositeNode *topNode)
1244 {
1245 if (Verbosity() > 0)
1246 {
1247 std::cout << "DijetTreeMaker::End(PHCompositeNode *topNode) This is the End..." << std::endl;
1248 }
1249 std::cout<<"Total events: "<<_i_event<<std::endl;
1250 _f->Write();
1251 _f->Close();
1252
1253 return Fun4AllReturnCodes::EVENT_OK;
1254 }
1255
1256
1257 int DijetTreeMaker::Reset(PHCompositeNode *topNode)
1258 {
1259 if (Verbosity() > 0)
1260 {
1261 std::cout << "DijetTreeMaker::Reset(PHCompositeNode *topNode) being Reset" << std::endl;
1262 }
1263 return Fun4AllReturnCodes::EVENT_OK;
1264 }
1265
1266 Double_t DijetTreeMaker::getDPHI(Double_t phi1, Double_t phi2) {
1267 Double_t dphi = phi1 - phi2;
1268
1269
1270 if ( dphi > TMath::Pi() )
1271 dphi = dphi - 2. * TMath::Pi();
1272 if ( dphi <= -TMath::Pi() )
1273 dphi = dphi + 2. * TMath::Pi();
1274
1275 if ( TMath::Abs(dphi) > TMath::Pi() ) {
1276 std::cout << " commonUtility::getDPHI error!!! dphi is bigger than TMath::Pi() " << std::endl;
1277 std::cout << " " << phi1 << ", " << phi2 << ", " << dphi << std::endl;
1278 }
1279
1280 return dphi;
1281 }
1282
1283 float DijetTreeMaker::get_eccentricity(std::vector<TVector3> hcaltowers, float oh_sum)
1284 {
1285 float avg[2] = { 0 };
1286 for (auto &constituent : hcaltowers)
1287 {
1288 avg[0]+=constituent.Eta();
1289 avg[1]+=constituent.Phi();
1290 }
1291
1292 avg[0]/=oh_sum;
1293 avg[1]/=oh_sum;
1294
1295 float cov[2][2] = {0};
1296 float cov1[2][2] = {0};
1297 float cov2[2][2] = {0};
1298
1299 for (int i = 0; i < 2; i++)
1300 {
1301 for (int j = 0; j < 2; j++)
1302 {
1303 for (auto &tower : hcaltowers)
1304 {
1305 float diff1 = 0;
1306 float diff2 = 0;
1307
1308 if (i == 0)
1309 diff1 = tower.Eta();
1310 else
1311 diff1 = tower.Phi();
1312 if (j == 0)
1313 diff2 = tower.Eta();
1314 else
1315 diff2 = tower.Phi();
1316
1317 cov[i][j] += tower.Pt()*diff1*diff2/oh_sum;
1318 cov1[i][j] += tower.Pt()*diff1/oh_sum;
1319 cov2[i][j] += tower.Pt()*diff2/oh_sum;
1320
1321 }
1322
1323 cov[i][j] -= cov1[i][j]*cov2[i][j];
1324 }
1325 }
1326
1327 if (fabs(cov[0][0]) < fabs(cov[1][1])) return 0.0;
1328
1329 float maxi = std::max( cov[0][0]*cov[0][0] , cov[1][1]*cov[1][1] );
1330 float mini = std::min( cov[0][0]*cov[0][0] , cov[1][1]*cov[1][1] );
1331 float eccen = sqrt(maxi - mini)/maxi;
1332 return eccen;
1333 }
