File indexing completed on 2026-04-07 08:08:32
0001 #include "AnNeutralMeson_micro_dst.h"
0002 #include "ClusterSmallInfoContainer.h"
0003 #include "ClusterSmallInfo.h"
0004 #include <fun4all/Fun4AllReturnCodes.h>
0005 #include <phool/getClass.h>
0006 #include <globalvertex/GlobalVertex.h>
0007
0008
0009 #include <odbc++/connection.h>
0010 #include <odbc++/drivermanager.h>
0011 #include <odbc++/resultset.h>
0012 #include <odbc++/resultsetmetadata.h>
0013 #include <uspin/SpinDBContent.h>
0014 #include <uspin/SpinDBContentv1.h>
0015 #include <uspin/SpinDBOutput.h>
0016
0017 #include <jetbase/JetContainer.h>
0018 #include <jetbase/Jet.h>
0019
0020 #include <format>
0021 #include <fstream>
0022 #include <iomanip>
0023 #include <iostream>
0024 #include <sstream>
0025 #include <numeric> // for std::iota
0026 #include <cstdlib> // for rand
0027 #include <cmath> // for fmod
0028
0029 #include <Math/Vector4D.h>
0030 #include <Math/Vector3D.h>
0031 #include <Math/GenVector/VectorUtil.h> // For DeltaR
0032 #include <Math/LorentzRotation.h>
0033 #include <Math/Rotation3D.h>
0034 #include <Math/AxisAngle.h>
0035 #include <TFile.h>
0036 #include <TH1F.h>
0037 #include <TH2F.h>
0038 #include <TTree.h>
0039 #include <bitset>
0040
0041
0042 #include <Riostream.h>
0043 #include <TSystem.h>
0044 #include <malloc.h>
0045
0046 void monitorMemoryUsage(const std::string &label)
0047 {
0048 ProcInfo_t procInfo;
0049 gSystem->GetProcInfo(&procInfo);
0050 std::cout << label << " - Memory usage: "
0051 << "Resident = " << procInfo.fMemResident << " kB, "
0052 << "Virtual = " << procInfo.fMemVirtual << " kB" << std::endl;
0053 }
0054
0055 AnNeutralMeson_micro_dst::AnNeutralMeson_micro_dst(const std::string &name, const int runnb,
0056 const std::string &outputfilename,
0057 const std::string &outputfiletreename)
0058 : SubsysReco(name)
0059 , runnumber(runnb)
0060 , outfilename(outputfilename)
0061 , outtreename(outputfiletreename)
0062 {
0063 }
0064
0065 AnNeutralMeson_micro_dst::~AnNeutralMeson_micro_dst()
0066 {
0067 delete outfile;
0068 delete outfile_tree;
0069 }
0070
0071 int AnNeutralMeson_micro_dst::Init(PHCompositeNode *)
0072 {
0073
0074 if (store_tree)
0075 {
0076 outfile_tree = new TFile(outtreename.c_str(),
0077 "RECREATE");
0078 output_tree = new TTree(
0079 "tree_diphoton_compact",
0080 "Minimal diphoton info");
0081 output_tree->Branch(
0082 "diphoton_vertex_z", &diphoton_vertex_z,
0083 "diphoton_vertex_z/F");
0084 output_tree->Branch(
0085 "diphoton_bunchnumber", &diphoton_bunchnumber,
0086 "diphoton_bunchnumber/I");
0087 output_tree->Branch(
0088 "diphoton_mass", &diphoton_mass,
0089 "diphoton_mass/F");
0090 output_tree->Branch(
0091 "diphoton_eta", &diphoton_eta,
0092 "diphoton_eta/F");
0093 output_tree->Branch(
0094 "diphoton_phi", &diphoton_phi,
0095 "diphoton_phi/F");
0096 output_tree->Branch(
0097 "diphoton_pt", &diphoton_pt,
0098 "diphoton_pt/F");
0099 output_tree->Branch(
0100 "diphoton_xf", &diphoton_xf,
0101 "diphoton_xf/F");
0102 }
0103
0104
0105 outfile = new TFile(outfilename.c_str(),
0106 "RECREATE");
0107
0108 if (require_emulator_matching) {
0109 h_emulator_selection = new TH1I(
0110 "h_emulator_selection",
0111 ";Trigger Tile Count; Count",
0112 10,0,10);
0113 }
0114
0115 if (require_efficiency_matching && require_emulator_matching)
0116 {
0117 h_efficiency_x_matching = new TH1I(
0118 "h_efficiency_x_matching",
0119 ";case;Counts",
0120 4, 1, 4);
0121
0122 for (int iPt = 0; iPt < nPtBins; iPt++) {
0123 std::stringstream h_efficiency_x_matching_name;
0124 h_efficiency_x_matching_name << "h_efficiency_x_matching_pt_" << iPt;
0125
0126 h_efficiency_x_matching_pt[iPt] = new TH1I(h_efficiency_x_matching_name.str().c_str(),
0127 ";case;Counts",
0128 4, 1, 4);
0129 }
0130 }
0131
0132 if (store_qa)
0133 {
0134
0135 h_multiplicity_efficiency = new TH1F(
0136 "h_multiplicity_efficiency",
0137 ";multiplicity; counts",
0138 10, 0, 10);
0139
0140 h_multiplicity_emulator = new TH1F(
0141 "h_multiplicity_emulator",
0142 ";multiplicity; counts",
0143 10, 0, 10);
0144
0145 h_event_zvtx = new TH1F(
0146 "h_event_zvtx",
0147 ";vertex z [cm]; counts",
0148 200, -200, 200);
0149
0150 h_triggered_event_zvtx = new TH1F(
0151 "h_triggered_event_zvtx",
0152 ";vertex z [cm]; counts",
0153 200, -200, 200);
0154
0155 h_triggered_event_photons = new TH1F(
0156 "h_triggered_event_photons",
0157 ";# photons; counts",
0158 20, 0, 10);
0159
0160 h_analysis_event_zvtx = new TH1F(
0161 "h_analysis_event_zvtx",
0162 ";vertex z [cm]; counts",
0163 200, -200, 200);
0164
0165
0166 h_trigger_live =
0167 new TH1F(
0168 "h_trigger_live",
0169 "Trigger Live;;Counts", 32, -0.5, 32 - 0.5);
0170 h_trigger_scaled = new TH1F(
0171 "h_trigger_scaled",
0172 "Trigger Scaled;;Counts", 32,
0173 -0.5, 32 - 0.5);
0174
0175 for (const auto &[trigger_number,
0176 trigger_name] : trigger_names)
0177 {
0178 h_trigger_live->GetXaxis()->SetBinLabel(trigger_number + 1,
0179 trigger_name.c_str());
0180 h_trigger_scaled->GetXaxis()->SetBinLabel(trigger_number + 1,
0181 trigger_name.c_str());
0182 }
0183
0184 h_count_diphoton_nomatch = new TH1I(
0185 "h_count_diphoton_nomatch",
0186 ";>= 1 diphoton; Count",
0187 1,0,1);
0188
0189
0190 h_current_E = new TH1F(
0191 "h_current_E",
0192 ";E [GeV]; counts",
0193 200, 0, 20);
0194 h_current_eta = new TH1F(
0195 "h_current_eta",
0196 ";#eta [rad];counts",
0197 200, -2.0, 2.0);
0198 h_current_phi = new TH1F(
0199 "h_current_phi",
0200 ";#phi [rad];counts",
0201 128, -M_PI, M_PI);
0202 h_pool_E = new TH1F(
0203 "h_pool_E",
0204 ";E [GeV]; counts",
0205 200, 0, 20);
0206 h_pool_eta = new TH1F(
0207 "h_pool_eta",
0208 ";#eta [rad];counts",
0209 200, -2.0, 2.0);
0210 h_pool_phi = new TH1F(
0211 "h_pool_phi",
0212 ";#phi [rad];counts",
0213 128, -M_PI, M_PI);
0214
0215 h_same_delta_eta = new TH1F(
0216 "h_same_delta_eta",
0217 ";#Delta #eta [rad];counts",
0218 500, 0, 10);
0219 h_same_delta_phi = new TH1F(
0220 "h_same_delta_phi",
0221 ";#Delta #phi [rad];counts",
0222 500, 0, 10);
0223 h_same_delta_R = new TH1F(
0224 "h_same_delta_R",
0225 ";#Delta R [rad];counts",
0226 500, 0, 10);
0227 h_same_alpha = new TH1F(
0228 "h_same_alpha",
0229 ";#alpha [rad];counts",
0230 500, 0, 1);
0231 h_mixed_delta_eta = new TH1F(
0232 "h_mixed_delta_eta",
0233 ";#Delta #eta [rad];counts",
0234 500, 0, 10);
0235 h_mixed_delta_phi = new TH1F(
0236 "h_mixed_delta_phi",
0237 ";#Delta #phi [rad];counts",
0238 500, 0, 10);
0239 h_mixed_delta_R = new TH1F(
0240 "h_mixed_delta_R",
0241 ";#Delta R [rad];counts",
0242 500, 0, 10);
0243 h_mixed_alpha = new TH1F(
0244 "h_mixed_alpha",
0245 ";#alpha [rad];counts",
0246 500, 0, 1);
0247
0248
0249 h_photon_eta = new TH1F(
0250 "h_photon_eta",
0251 ";#eta [rad];counts",
0252 200, -2.0, 2.0);
0253 h_photon_phi = new TH1F(
0254 "h_photon_phi",
0255 ";#phi [rad];counts",
0256 128, -M_PI, M_PI);
0257 h_photon_pt = new TH1F(
0258 "h_photon_pt",
0259 ";p_{T} [GeV];counts",
0260 200, 0, 20);
0261 h_photon_zvtx = new TH1F(
0262 "h_photon_zvtx",
0263 ";z vertex [cm];counts",
0264 200, -200, 200);
0265
0266 h_photon_eta_phi = new TH2F(
0267 "h_photon_eta_phi",
0268 ";#eta [rad];#phi [rad]",
0269 200, -2.0, 2.0, 128, -M_PI, M_PI);
0270 h_photon_eta_pt = new TH2F(
0271 "h_photon_eta_pt",
0272 ";#eta [rad];p_{T} [GeV]",
0273 200, -2.0, 2.0, 200, 0, 20);
0274 h_photon_eta_zvtx = new TH2F(
0275 "h_photon_eta_zvtx",
0276 ";#eta [rad];vertex z [cm]; counts",
0277 200, -2, 2, 200, -200, 200);
0278 h_photon_phi_pt = new TH2F(
0279 "h_photon_phi_pt",
0280 ";#phi [rad];p_{T} [GeV]; counts",
0281 128, -M_PI, M_PI, 200, 0, 20);
0282 h_photon_phi_zvtx = new TH2F(
0283 "h_photon_phi_zvtx",
0284 ";#phi [rad];vertex z [cm]; counts",
0285 128, -M_PI, M_PI, 200, -200, 200);
0286 h_photon_pt_zvtx = new TH2F(
0287 "h_photon_pt_zvtx",
0288 ";p_{T} [GeV];vertex z [cm]; counts",
0289 200, 0, 20, 200, -200, 200);
0290
0291 h_selected_photon_eta = new TH1F(
0292 "h_selected_photon_eta",
0293 ";#eta [rad];counts",
0294 200, -2.0, 2.0);
0295 h_selected_photon_phi = new TH1F(
0296 "h_selected_photon_phi",
0297 ";#phi [rad];counts",
0298 128, -M_PI, M_PI);
0299 h_selected_photon_pt = new TH1F(
0300 "h_selected_photon_pt",
0301 ";p_{T} [GeV];counts",
0302 200, 0, 20);
0303 h_selected_photon_zvtx = new TH1F(
0304 "h_selected_photon_zvtx",
0305 ";z vertex [cm];counts",
0306 200, -200, 200);
0307
0308 h_selected_photon_eta_phi = new TH2F(
0309 "h_selected_photon_eta_phi",
0310 ";#eta [rad];#phi [rad]",
0311 200, -2.0, 2.0, 128, -M_PI, M_PI);
0312 h_selected_photon_eta_pt = new TH2F(
0313 "h_selected_photon_eta_pt",
0314 ";#eta [rad];p_{T} [GeV]",
0315 200, -2.0, 2.0, 200, 0, 20);
0316 h_selected_photon_eta_zvtx = new TH2F(
0317 "h_selected_photon_eta_zvtx",
0318 ";#eta [rad];vertex z [cm]; counts",
0319 200, -2, 2, 200, -200, 200);
0320 h_selected_photon_phi_pt = new TH2F(
0321 "h_selected_photon_phi_pt",
0322 ";#phi [rad];p_{T} [GeV]; counts",
0323 128, -M_PI, M_PI, 200, 0, 20);
0324 h_selected_photon_phi_zvtx = new TH2F(
0325 "h_selected_photon_phi_zvtx",
0326 ";#phi [rad];vertex z [cm]; counts",
0327 128, -M_PI, M_PI, 200, -200, 200);
0328 h_selected_photon_pt_zvtx = new TH2F(
0329 "h_selected_photon_pt_zvtx",
0330 ";p_{T} [GeV];vertex z [cm]; counts",
0331 200, 0, 20, 200, -200, 200);
0332
0333
0334 h_pair_eta = new TH1F(
0335 "h_pair_eta",
0336 ";#eta [rad];counts",
0337 200, -2.0, 2.0);
0338 h_pair_phi = new TH1F(
0339 "h_pair_phi",
0340 ";#phi [rad];counts",
0341 128, -M_PI, M_PI);
0342 h_pair_pt = new TH1F(
0343 "h_pair_pt",
0344 ";p_{T} [GeV];counts",
0345 200, 0, 20);
0346 h_pair_zvtx = new TH1F(
0347 "h_pair_zvtx",
0348 ";z vertex [cm];counts",
0349 200, -200, 200);
0350 h_pair_DeltaR = new TH1F(
0351 "h_pair_DeltaR",
0352 ";#Delta R [rad];counts",
0353 500, 0, 5);
0354 h_pair_alpha = new TH1F(
0355 "h_pair_alpha",
0356 ";#alpha [rad];counts",
0357 200, 0, 2);
0358 h_pair_alpha_pt = new TH2F(
0359 "h_pair_alpha_pt",
0360 ";#alpha; p_{T}; counts",
0361 200, 0, 2, 200, 0, 20);
0362 h_pair_eta_phi = new TH2F(
0363 "h_pair_eta_phi",
0364 ";#eta [rad];#phi [rad]",
0365 200, -2.0, 2.0, 128, -M_PI, M_PI);
0366 h_pair_eta_pt = new TH2F(
0367 "h_pair_eta_pt",
0368 ";#eta [rad];p_{T} [GeV]",
0369 200, -2.0, 2.0, 200, 0, 20);
0370 h_pair_eta_zvtx = new TH2F(
0371 "h_pair_eta_zvtx",
0372 ";#eta [rad];vertex z [cm]; counts",
0373 200, -2, 2, 200, -200, 200);
0374 h_pair_phi_pt = new TH2F(
0375 "h_pair_phi_pt",
0376 ";#phi [rad];p_{T} [GeV]; counts",
0377 128, -M_PI, M_PI, 200, 0, 20);
0378 h_pair_phi_zvtx = new TH2F(
0379 "h_pair_phi_zvtx",
0380 ";#phi [rad];vertex z [cm]; counts",
0381 128, -M_PI, M_PI, 200, -200, 200);
0382 h_pair_pt_zvtx = new TH2F(
0383 "h_pair_pt_zvtx",
0384 ";p_{T} [GeV];vertex z [cm]; counts",
0385 200, 0, 20, 200, -200, 200);
0386
0387
0388 h_pair_pi0_eta_pt_1 = new TH1F(
0389 "h_pair_pi0_eta_pt_1",
0390 ";#eta [rad]; Counts / [20 mrad]",
0391 200,-2.0,2.0);
0392
0393 h_pair_pi0_eta_pt_2 = new TH1F(
0394 "h_pair_pi0_eta_pt_2",
0395 ";#eta [rad]; Counts / [20 mrad]",
0396 200,-2.0,2.0);
0397
0398 h_pair_pi0_eta_pt_3 = new TH1F(
0399 "h_pair_pi0_eta_pt_3",
0400 ";#eta [rad]; Counts / [20 mrad]",
0401 200,-2.0,2.0);
0402
0403 h_pair_pi0_eta_pt_4 = new TH1F(
0404 "h_pair_pi0_eta_pt_4",
0405 ";#eta [rad]; Counts / [20 mrad]",
0406 200,-2.0,2.0);
0407
0408 h_pair_pi0_eta_pt_5 = new TH1F(
0409 "h_pair_pi0_eta_pt_5",
0410 ";#eta [rad]; Counts / [20 mrad]",
0411 200,-2.0,2.0);
0412
0413 h_pair_pi0_eta_pt_6 = new TH1F(
0414 "h_pair_pi0_eta_pt_6",
0415 ";#eta [rad]; Counts / [20 mrad]",
0416 200,-2.0,2.0);
0417
0418 h_pair_pi0_eta_pt_7 = new TH1F(
0419 "h_pair_pi0_eta_pt_7",
0420 ";#eta [rad]; Counts / [20 mrad]",
0421 200,-2.0,2.0);
0422
0423 h_pair_pi0_eta_pt_8 = new TH1F(
0424 "h_pair_pi0_eta_pt_8",
0425 ";#eta [rad]; Counts / [20 mrad]",
0426 200,-2.0,2.0);
0427
0428 h_pair_pi0_eta_pt_9 = new TH1F(
0429 "h_pair_pi0_eta_pt_9",
0430 ";#eta [rad]; Counts / [20 mrad]",
0431 200,-2.0,2.0);
0432
0433 h_pair_eta_eta_pt_1 = new TH1F(
0434 "h_pair_eta_eta_pt_1",
0435 ";#eta [rad]; Counts / [20 mrad]",
0436 200,-2.0,2.0);
0437
0438 h_pair_eta_eta_pt_2 = new TH1F(
0439 "h_pair_eta_eta_pt_2",
0440 ";#eta [rad]; Counts / [20 mrad]",
0441 200,-2.0,2.0);
0442
0443 h_pair_eta_eta_pt_3 = new TH1F(
0444 "h_pair_eta_eta_pt_3",
0445 ";#eta [rad]; Counts / [20 mrad]",
0446 200,-2.0,2.0);
0447
0448 h_pair_eta_eta_pt_4 = new TH1F(
0449 "h_pair_eta_eta_pt_4",
0450 ";#eta [rad]; Counts / [20 mrad]",
0451 200,-2.0,2.0);
0452
0453 h_pair_eta_eta_pt_5 = new TH1F(
0454 "h_pair_eta_eta_pt_5",
0455 ";#eta [rad]; Counts / [20 mrad]",
0456 200,-2.0,2.0);
0457
0458 h_pair_eta_eta_pt_6 = new TH1F(
0459 "h_pair_eta_eta_pt_6",
0460 ";#eta [rad]; Counts / [20 mrad]",
0461 200,-2.0,2.0);
0462
0463 h_pair_eta_eta_pt_7 = new TH1F(
0464 "h_pair_eta_eta_pt_7",
0465 ";#eta [rad]; Counts / [20 mrad]",
0466 200,-2.0,2.0);
0467
0468 h_pair_eta_eta_pt_8 = new TH1F(
0469 "h_pair_eta_eta_pt_8",
0470 ";#eta [rad]; Counts / [20 mrad]",
0471 200,-2.0,2.0);
0472
0473 h_pair_eta_eta_pt_9 = new TH1F(
0474 "h_pair_eta_eta_pt_9",
0475 ";#eta [rad]; Counts / [20 mrad]",
0476 200,-2.0,2.0);
0477
0478
0479
0480 h_pair_pi0_xf_pt_1 = new TH1F(
0481 "h_pair_pi0_xf_pt_1",
0482 ";x_{F} [rad]; Counts / [2 mrad]",
0483 200,-0.2,0.2);
0484
0485 h_pair_pi0_xf_pt_2 = new TH1F(
0486 "h_pair_pi0_xf_pt_2",
0487 ";x_{F} [rad]; Counts / [2 mrad]",
0488 200,-0.2,0.2);
0489
0490 h_pair_pi0_xf_pt_3 = new TH1F(
0491 "h_pair_pi0_xf_pt_3",
0492 ";x_{F} [rad]; Counts / [2 mrad]",
0493 200,-0.2,0.2);
0494
0495 h_pair_pi0_xf_pt_4 = new TH1F(
0496 "h_pair_pi0_xf_pt_4",
0497 ";x_{F} [rad]; Counts / [2 mrad]",
0498 200,-0.2,0.2);
0499
0500 h_pair_pi0_xf_pt_5 = new TH1F(
0501 "h_pair_pi0_xf_pt_5",
0502 ";x_{F} [rad]; Counts / [2 mrad]",
0503 200,-0.2,0.2);
0504
0505 h_pair_pi0_xf_pt_6 = new TH1F(
0506 "h_pair_pi0_xf_pt_6",
0507 ";x_{F} [rad]; Counts / [2 mrad]",
0508 200,-0.2,0.2);
0509
0510 h_pair_pi0_xf_pt_7 = new TH1F(
0511 "h_pair_pi0_xf_pt_7",
0512 ";x_{F} [rad]; Counts / [2 mrad]",
0513 200,-0.2,0.2);
0514
0515 h_pair_pi0_xf_pt_8 = new TH1F(
0516 "h_pair_pi0_xf_pt_8",
0517 ";x_{F} [rad]; Counts / [2 mrad]",
0518 200,-0.2,0.2);
0519
0520 h_pair_pi0_xf_pt_9 = new TH1F(
0521 "h_pair_pi0_xf_pt_9",
0522 ";x_{F} [rad]; Counts / [2 mrad]",
0523 200,-0.2,0.2);
0524
0525 h_pair_eta_xf_pt_1 = new TH1F(
0526 "h_pair_eta_xf_pt_1",
0527 ";x_{F} [rad]; Counts / [2 mrad]",
0528 200,-0.2,0.2);
0529
0530 h_pair_eta_xf_pt_2 = new TH1F(
0531 "h_pair_eta_xf_pt_2",
0532 ";x_{F} [rad]; Counts / [2 mrad]",
0533 200,-0.2,0.2);
0534
0535 h_pair_eta_xf_pt_3 = new TH1F(
0536 "h_pair_eta_xf_pt_3",
0537 ";x_{F} [rad]; Counts / [2 mrad]",
0538 200,-0.2,0.2);
0539
0540 h_pair_eta_xf_pt_4 = new TH1F(
0541 "h_pair_eta_xf_pt_4",
0542 ";x_{F} [rad]; Counts / [2 mrad]",
0543 200,-0.2,0.2);
0544
0545 h_pair_eta_xf_pt_5 = new TH1F(
0546 "h_pair_eta_xf_pt_5",
0547 ";x_{F} [rad]; Counts / [2 mrad]",
0548 200,-0.2,0.2);
0549
0550 h_pair_eta_xf_pt_6 = new TH1F(
0551 "h_pair_eta_xf_pt_6",
0552 ";x_{F} [rad]; Counts / [2 mrad]",
0553 200,-0.2,0.2);
0554
0555 h_pair_eta_xf_pt_7 = new TH1F(
0556 "h_pair_eta_xf_pt_7",
0557 ";x_{F} [rad]; Counts / [2 mrad]",
0558 200,-0.2,0.2);
0559
0560 h_pair_eta_xf_pt_8 = new TH1F(
0561 "h_pair_eta_xf_pt_8",
0562 ";x_{F} [rad]; Counts / [2 mrad]",
0563 200,-0.2,0.2);
0564
0565 h_pair_eta_xf_pt_9 = new TH1F(
0566 "h_pair_eta_xf_pt_9",
0567 ";x_{F} [rad]; Counts / [2 mrad]",
0568 200,-0.2,0.2);
0569
0570 h_meson_pi0_pt = new TH1F(
0571 "h_meson_pi0_pt",
0572 ";p_{T} [GeV];counts",
0573 200, 0, 20);
0574
0575 h_meson_pi0_E = new TH1F(
0576 "h_meson_pi0_E",
0577 ";E [GeV];counts",
0578 200, 0, 20);
0579
0580 h_meson_eta_pt = new TH1F(
0581 "h_meson_eta_pt",
0582 ";p_{T} [GeV];counts",
0583 200, 0, 20);
0584
0585 h_meson_eta_E = new TH1F(
0586 "h_meson_eta_E",
0587 ";E [GeV];counts",
0588 200, 0, 20);
0589 }
0590
0591
0592 h_pair_mass = new TH1F(
0593 "h_pair_mass",
0594 ";M_{#gamma} [GeV];counts", 500, 0, 1);
0595
0596 h_pair_mass_mixing = new TH1F(
0597 "h_pair_mass_mixing",
0598 ";M_{#gamma #gamma} [GeV];counts", 500, 0, 1);
0599
0600 for (int iPt = 0; iPt < nPtBins; iPt++)
0601 {
0602 std::stringstream h_pair_mass_name;
0603 h_pair_mass_name << std::fixed << std::setprecision(0) << "h_pair_mass_pt_"
0604 << iPt;
0605 std::stringstream h_pair_mass_title;
0606 h_pair_mass_title << std::fixed << std::setprecision(1) << "diphoton mass ["
0607 << pTBins[iPt] << " < p_{T} [GeV/c] < "
0608 << pTBins[iPt + 1]
0609 << "];M_{#gamma#gamma} [GeV/c^{2}];counts";
0610 h_pair_mass_pt[iPt] = new TH1F(h_pair_mass_name.str().c_str(),
0611 h_pair_mass_title.str().c_str(), 500, 0, 1);
0612 h_pair_mass_name.str("");
0613 h_pair_mass_name << std::fixed << std::setprecision(0) << "h_pair_mass_mixing_pt_"
0614 << iPt;
0615 h_pair_mass_mixing_pt[iPt] = new TH1F(h_pair_mass_name.str().c_str(),
0616 h_pair_mass_title.str().c_str(), 500, 0, 1);
0617 }
0618
0619 for (int izvtx = 0; izvtx < nZvtxBins; izvtx++)
0620 {
0621 std::stringstream h_pair_mass_name;
0622 h_pair_mass_name << std::fixed << std::setprecision(0) << "h_pair_mass_zvtx_"
0623 << izvtx;
0624 std::stringstream h_pair_mass_title;
0625 h_pair_mass_title << std::fixed << std::setprecision(1) << "diphoton mass ["
0626 << - zvtxBins[izvtx + 1] << " < z_{vtx} [cm] < "
0627 << zvtxBins[izvtx + 1]
0628 << "];M_{#gamma#gamma} [GeV/c^{2}];counts";
0629 h_pair_mass_zvtx[izvtx] = new TH1F(h_pair_mass_name.str().c_str(),
0630 h_pair_mass_title.str().c_str(), 500, 0, 1);
0631 }
0632
0633 for (int ietaBin = 0; ietaBin < nEtaBins; ietaBin++)
0634 {
0635 std::stringstream h_pair_mass_name;
0636 h_pair_mass_name << std::fixed << std::setprecision(0) << "h_pair_mass_eta_"
0637 << ietaBin;
0638 std::stringstream h_pair_mass_title;
0639 h_pair_mass_title << std::fixed << std::setprecision(2) << "diphoton mass ["
0640 << etaBins[ietaBin] << " < #eta < "
0641 << etaBins[ietaBin + 1] << "];M_{#gamma#gamma};counts";
0642 h_pair_mass_eta[ietaBin] =
0643 new TH1F(h_pair_mass_name.str().c_str(),
0644 h_pair_mass_title.str().c_str(), 500, 0, 1);
0645 }
0646
0647 for (int ixfBin = 0; ixfBin < nXfBins; ixfBin++)
0648 {
0649 std::stringstream h_pair_mass_name;
0650 h_pair_mass_name << std::fixed << std::setprecision(0) << "h_pair_mass_xf_"
0651 << ixfBin;
0652 std::stringstream h_pair_mass_title;
0653 h_pair_mass_title << std::fixed << std::setprecision(2) << "diphoton mass ["
0654 << xfBins[ixfBin] << " < x_{F} < " << xfBins[ixfBin + 1]
0655 << "];M_{#gamma#gamma};counts";
0656 h_pair_mass_xf[ixfBin] =
0657 new TH1F(h_pair_mass_name.str().c_str(),
0658 h_pair_mass_title.str().c_str(), 500, 0, 1);
0659 }
0660
0661
0662 h_average_pt[0] = new TH1F(
0663 "h_average_pt_pi0",
0664 ";iPt; pT average [GeV/c]", nPtBins,
0665 0, nPtBins);
0666 h_average_eta[0] = new TH1F(
0667 "h_average_eta_pi0",
0668 ";ieta; #eta average", nEtaBins, 0,
0669 nEtaBins);
0670 h_average_xf[0] =
0671 new TH1F(
0672 "h_average_xf_pi0",
0673 ";ixf; pT average", nXfBins, 0, nXfBins);
0674 h_norm_pt[0] = new TH1F(
0675 "h_norm_pt_pi0",
0676 ";iPt; pT norm [GeV/c]", nPtBins,
0677 0, nPtBins);
0678 h_norm_eta[0] = new TH1F(
0679 "h_norm_eta_pi0",
0680 ";ieta; #eta norm", nEtaBins, 0,
0681 nEtaBins);
0682 h_norm_xf[0] =
0683 new TH1F(
0684 "h_norm_xf_pi0",
0685 ";ixf; pT norm", nXfBins, 0, nXfBins);
0686 h_average_pt[1] = new TH1F(
0687 "h_average_pt_eta",
0688 ";iPt; pT average [GeV/c]", nPtBins,
0689 0, nPtBins);
0690 h_average_eta[1] = new TH1F(
0691 "h_average_eta_eta",
0692 ";ieta; #eta average", nEtaBins, 0,
0693 nEtaBins);
0694 h_average_xf[1] =
0695 new TH1F(
0696 "h_average_xf_eta",
0697 ";ixf; pT average", nXfBins, 0, nXfBins);
0698 h_norm_pt[1] = new TH1F(
0699 "h_norm_pt_eta",
0700 ";iPt; pT norm [GeV/c]", nPtBins,
0701 0, nPtBins);
0702 h_norm_eta[1] = new TH1F(
0703 "h_norm_eta_eta",
0704 ";ieta; #eta norm", nEtaBins, 0,
0705 nEtaBins);
0706 h_norm_xf[1] =
0707 new TH1F(
0708 "h_norm_xf_eta",
0709 ";ixf; pT norm", nXfBins, 0, nXfBins);
0710
0711
0712 for (int iB = 0; iB < nBeams; iB++)
0713 {
0714 for (int iP = 0; iP < nParticles; iP++)
0715 {
0716 for (int iR = 0; iR < nRegions; iR++)
0717 {
0718 for (int iS = 0; iS < nSpins; iS++)
0719 {
0720 for (int iPt = 0; iPt < nPtBins; iPt++)
0721 {
0722 for (int ieta = 0; ieta < nEtaRegions; ieta++)
0723 {
0724
0725 std::stringstream h_yield_name;
0726 h_yield_name << "h_yield_" << beams[iB] << "_" << particle[iP]
0727 << "_" << regions[iR] << "_"
0728 << "pT_" << iPt << "_"
0729 << etaRegions[ieta]
0730 << "_" << spins[iS];
0731 std::stringstream h_yield_title;
0732 h_yield_title << std::fixed << std::setprecision(1)
0733 << "P_{T} #in [" << pTBins[iPt] << ", "
0734 << pTBins[iPt + 1] << "] "
0735 << (ieta == 0 ? "|#eta| < 0.35"
0736 : "|#eta| > 0.35")
0737 << " " << (iP == 0 ? "(#pi^{0} " : "(#eta ")
0738 << regions[iR] << " range);"
0739 << "#phi_{" << (iB == 0 ? "yellow" : "blue")
0740 << "};counts";
0741 h_yield_1[iB][iP][iR][iPt][ieta][iS] =
0742 new TH1F(h_yield_name.str().c_str(),
0743 h_yield_title.str().c_str(), 12, -M_PI, M_PI);
0744 for (int iDir = 0; iDir < nDirections; iDir++)
0745 {
0746 if (ieta == 0)
0747 {
0748
0749 h_yield_name.str(
0750 "");
0751 h_yield_name << "h_yield_" << beams[iB] << "_" << particle[iP]
0752 << "_" << regions[iR] << "_"
0753 << "pT_" << iPt << "_"
0754 << directions[iDir]
0755 << "_" << spins[iS];
0756 h_yield_title.str(
0757 "");
0758 h_yield_title << std::fixed << std::setprecision(1)
0759 << "P_{T} #in [" << pTBins[iPt] << ", "
0760 << pTBins[iPt + 1] << "] " << directions[iDir]
0761 << " " << (iP == 0 ? "(#pi^{0} " : "(#eta ")
0762 << regions[iR] << " range);"
0763 << "#phi_{" << (iB == 0 ? "yellow" : "blue")
0764 << "};counts";
0765 h_yield_2[iB][iP][iR][iPt][iDir][iS] =
0766 new TH1F(h_yield_name.str().c_str(),
0767 h_yield_title.str().c_str(), 12, -M_PI, M_PI);
0768 }
0769
0770 h_yield_name.str(
0771 "");
0772 h_yield_name << "h_yield_" << beams[iB] << "_" << particle[iP]
0773 << "_" << regions[iR] << "_"
0774 << "pT_" << iPt << "_"
0775 << etaRegions[ieta]
0776 << "_" << directions[iDir] << "_" << spins[iS];
0777 h_yield_title.str(
0778 "");
0779 h_yield_title << std::fixed << std::setprecision(1)
0780 << "P_{T} #in [" << pTBins[iPt] << ", "
0781 << pTBins[iPt + 1] << "] "
0782 << (ieta == 0 ? "|#eta| < 0.35"
0783 : "|#eta| > 0.35")
0784 << " " << directions[iDir] << " "
0785 << (iP == 0 ? "(#pi^{0} "
0786 : "(#eta ")
0787 << regions[iR] << " range);"
0788 << "#phi_{" << (iB == 0 ? "yellow" : "blue")
0789 << "};counts";
0790 h_yield_3[iB][iP][iR][iPt][ieta][iDir][iS] =
0791 new TH1F(h_yield_name.str().c_str(),
0792 h_yield_title.str().c_str(), 12, -M_PI, M_PI);
0793 }
0794 }
0795 }
0796 }
0797 }
0798 }
0799 }
0800
0801
0802 for (int iB = 0; iB < nBeams; iB++)
0803 {
0804 for (int iP = 0; iP < nParticles; iP++)
0805 {
0806 for (int iR = 0; iR < nRegions; iR++)
0807 {
0808 for (int iS = 0; iS < nSpins; iS++)
0809 {
0810 for (int ietaBin = 0; ietaBin < nEtaBins; ietaBin++)
0811 {
0812 std::stringstream h_yield_name;
0813 h_yield_name << "h_yield_" << beams[iB] << "_" << particle[iP]
0814 << "_" << regions[iR] << "_"
0815 << "eta_" << (int) ietaBin << "_" << spins[iS];
0816 std::stringstream h_yield_title;
0817 h_yield_title << std::fixed << std::setprecision(1) << "#eta #in ["
0818 << etaBins[ietaBin] << ", " << etaBins[ietaBin]
0819 << "] " << (iP == 0 ? "(#pi^{0} " : "(#eta ")
0820 << regions[iR] << " range);"
0821 << "#phi_{" << (iB == 0 ? "yellow" : "blue")
0822 << "};counts";
0823 h_yield_eta[iB][iP][iR][ietaBin][iS] =
0824 new TH1F(h_yield_name.str().c_str(),
0825 h_yield_title.str().c_str(), 12, -M_PI, M_PI);
0826 }
0827 }
0828 }
0829 }
0830 }
0831
0832
0833 for (int iB = 0; iB < nBeams; iB++)
0834 {
0835 for (int iP = 0; iP < nParticles; iP++)
0836 {
0837 for (int iR = 0; iR < nRegions; iR++)
0838 {
0839 for (int iS = 0; iS < nSpins; iS++)
0840 {
0841 for (int ixfBin = 0; ixfBin < nXfBins; ixfBin++)
0842 {
0843 std::stringstream h_yield_name;
0844 h_yield_name << "h_yield_" << beams[iB] << "_" << particle[iP]
0845 << "_" << regions[iR] << "_"
0846 << "xf_" << (int) ixfBin << "_" << spins[iS];
0847 std::stringstream h_yield_title;
0848 h_yield_title << std::fixed << std::setprecision(1) << "x_{F} #in ["
0849 << xfBins[ixfBin] << ", " << xfBins[ixfBin + 1]
0850 << "] " << (iP == 0 ? "(#pi^{0} " : "(#eta ")
0851 << regions[iR] << " range);"
0852 << "#phi_{" << (iB == 0 ? "yellow" : "blue")
0853 << "};counts";
0854 h_yield_xf[iB][iP][iR][ixfBin][iS] =
0855 new TH1F(h_yield_name.str().c_str(),
0856 h_yield_title.str().c_str(), 12, -M_PI, M_PI);
0857 }
0858 }
0859 }
0860 }
0861 }
0862
0863
0864 chi2_cuts_labels.resize(n_chi2_cuts);
0865 ecore_cuts_labels.resize(n_ecore_cuts);
0866
0867 for (int i = 0; i < n_chi2_cuts; i++)
0868 {
0869 std::stringstream label;
0870 label << std::fixed << std::setprecision(1) << chi2_cuts[i];
0871 chi2_cuts_labels[i] = label.str();
0872 }
0873
0874 for (int i = 0; i < n_ecore_cuts; i++)
0875 {
0876 std::stringstream label;
0877 label << std::fixed << std::setprecision(2) << ecore_cuts[i] << " GeV";
0878 ecore_cuts_labels[i] = label.str();
0879 }
0880
0881
0882
0883 if (n_chi2_cuts > 1 || n_ecore_cuts > 1)
0884 {
0885 h_cluster_level_cuts_total = new TH2F(
0886 "h_cluster_level_cuts_total",
0887 "Total number of Diphotons;Cluster #chi^{2} Cut;Min Cluster Ecore Cut",
0888 n_chi2_cuts, 0, (float) n_chi2_cuts, n_ecore_cuts, 0, (float) n_ecore_cuts);
0889 for (int il = 0; il < n_chi2_cuts; il++)
0890 h_cluster_level_cuts_total->GetXaxis()->SetBinLabel(
0891 il + 1, chi2_cuts_labels[il].c_str());
0892 for (int il = 0; il < n_ecore_cuts; il++)
0893 h_cluster_level_cuts_total->GetYaxis()->SetBinLabel(
0894 il + 1, ecore_cuts_labels[il].c_str());
0895
0896 for (int iP = 0; iP < nParticles; iP++)
0897 {
0898 for (int iR = 0; iR < nRegions; iR++)
0899 {
0900 std::stringstream h_name;
0901 h_name << "h_cluster_level_cuts_" << particle[iP] << "_" << regions[iR];
0902 std::stringstream h_title;
0903 h_title << "Total Number of Diphotons (" << (iP == 0 ? "#pi^{0}" : "#eta")
0904 << " " << regions[iR]
0905 << ");Cluster #chi^{2} Cut;Min Cluster Energy Cut";
0906 h_cluster_level_cuts_particle[iP][iR] =
0907 new TH2F(h_name.str().c_str(), h_title.str().c_str(), n_chi2_cuts, 0,
0908 (float) n_chi2_cuts, n_ecore_cuts, 0, (float) n_ecore_cuts);
0909 for (int il = 0; il < n_chi2_cuts; il++)
0910 h_cluster_level_cuts_particle[iP][iR]->GetXaxis()->SetBinLabel(
0911 il + 1, chi2_cuts_labels[il].c_str());
0912 for (int il = 0; il < n_ecore_cuts; il++)
0913 h_cluster_level_cuts_particle[iP][iR]->GetYaxis()->SetBinLabel(
0914 il + 1, ecore_cuts_labels[il].c_str());
0915
0916 for (int iPt = 0; iPt < nPtBins; iPt++)
0917 {
0918 if ((iP == 0) && (iR == 0))
0919 {
0920 h_name.str("");
0921 h_name << "h_cluster_level_cuts_total_pt_" << iPt;
0922 h_title.str("");
0923 h_title << "Total Number of Diphotons (p_{T} #in [" << pTBins[iPt]
0924 << "," << pTBins[iPt + 1]
0925 << "] GeV/c);Cluster #chi^{2} Cut;Min Cluster Energy Cut";
0926 h_cluster_level_cuts_total_pt[iPt] = new TH2F(
0927 h_name.str().c_str(), h_title.str().c_str(), n_chi2_cuts, 0,
0928 (float) n_chi2_cuts, n_ecore_cuts, 0, (float) n_ecore_cuts);
0929 for (int il = 0; il < n_chi2_cuts; il++)
0930 h_cluster_level_cuts_total_pt[iPt]->GetXaxis()->SetBinLabel(
0931 il + 1, chi2_cuts_labels[il].c_str());
0932 for (int il = 0; il < n_ecore_cuts; il++)
0933 h_cluster_level_cuts_total_pt[iPt]->GetYaxis()->SetBinLabel(
0934 il + 1, ecore_cuts_labels[il].c_str());
0935 }
0936 h_name.str(
0937 "");
0938 h_name << "h_cluster_level_cuts_" << particle[iP] << "_" << regions[iR]
0939 << "_pt_" << iPt;
0940 h_title.str(
0941 "");
0942 h_title << "Total Number of Diphotons ("
0943 << (iP == 0 ? "#pi^{0}" : "#eta")
0944 << " " << regions[iR]
0945 << ", p_{T} #in [" << pTBins[iPt] << "," << pTBins[iPt + 1]
0946 << "] GeV/c);Cluster #chi^{2} Cut;Min Cluster Energy Cut";
0947 h_cluster_level_cuts_particle_pt[iP][iR][iPt] = new TH2F(
0948 h_name.str().c_str(), h_title.str().c_str(), n_chi2_cuts, 0,
0949 (float) n_chi2_cuts, n_ecore_cuts, 0, (float) n_ecore_cuts);
0950 for (int il = 0; il < n_chi2_cuts; il++)
0951 h_cluster_level_cuts_particle_pt[iP][iR][iPt]
0952 ->GetXaxis()
0953 ->SetBinLabel(il + 1, chi2_cuts_labels[il].c_str());
0954 for (int il = 0; il < n_ecore_cuts; il++)
0955 h_cluster_level_cuts_particle_pt[iP][iR][iPt]
0956 ->GetYaxis()
0957 ->SetBinLabel(il + 1, ecore_cuts_labels[il].c_str());
0958 }
0959 }
0960 }
0961 }
0962
0963 return 0;
0964 }
0965
0966 int AnNeutralMeson_micro_dst::InitRun(PHCompositeNode *topNode)
0967 {
0968 syncobject = findNode::getClass<SyncObject>(topNode, syncdefs::SYNCNODENAME);
0969 if (!syncobject)
0970 {
0971 std::cout << "syncobject not found" << std::endl;
0972 }
0973
0974 towergeom =
0975 findNode::getClass<RawTowerGeomContainer>(topNode,
0976 "TOWERGEOM_CEMC");
0977
0978 if (require_emulator_matching)
0979 {
0980 emcaltiles =
0981 findNode::getClass<TriggerTile>(topNode,
0982 "TILES_EMCAL");
0983 if (!emcaltiles)
0984 {
0985 std::cerr << "AnNeutralMeson_micro_dst TILES_EMCAL node is missing for emulator matching" << "\n";
0986 std::cerr << "Either load DST_TRIGGER_EMULATOR files with Fun4AllDstInputManager" << "\n";
0987 std::cerr << "Or call method set_photon_trigger_emulator_matching_requirement(false)" << std::endl;
0988 }
0989 }
0990
0991 vertexmap =
0992 findNode::getClass<GlobalVertexMap>(topNode,
0993 "GlobalVertexMap");
0994 if (!vertexmap)
0995 {
0996 std::cerr << "AnNeutralMeson_micro_dst GlobalVertexMap node is missing"
0997 << std::endl;
0998 return Fun4AllReturnCodes::ABORTRUN;
0999 }
1000
1001 _smallclusters = findNode::getClass<ClusterSmallInfoContainer>(topNode, "CLUSTER_SMALLINFO_CEMC");
1002 if (!_smallclusters)
1003 {
1004 std::cerr << "AnNeutralMeson_micro_dst CLUSTER_SMALLINFO_CEMC node is missing"
1005 << std::endl;
1006 return Fun4AllReturnCodes::ABORTRUN;
1007 }
1008
1009
1010 if (n_chi2_cuts == 0 || n_ecore_cuts == 0)
1011 {
1012 std::cerr << "No chi2 or no ecore cuts defined. Abort run" << std::endl;
1013 return Fun4AllReturnCodes::ABORTRUN;
1014 }
1015
1016
1017 odbc::Connection *con = nullptr;
1018 try
1019 {
1020 con = odbc::DriverManager::getConnection(db_name, user_name, "");
1021 }
1022 catch (odbc::SQLException &e)
1023 {
1024 std::cout << std::format("Error: {}.", e.getMessage()) << std::endl;
1025 return Fun4AllReturnCodes::ABORTRUN;
1026 }
1027
1028 std::stringstream cmd;
1029 cmd << "SELECT * FROM " << table_name << " WHERE runnumber=" << runnumber << ";";
1030 odbc::Statement *stmt = con->createStatement();
1031 odbc::ResultSet *rs = nullptr;
1032 try
1033 {
1034 rs = stmt->executeQuery(cmd.str());
1035 }
1036 catch (odbc::SQLException &e)
1037 {
1038 std::cout << std::format("Error: {}.", e.getMessage()) << std::endl;
1039 delete rs;
1040 delete stmt;
1041 delete con;
1042 return Fun4AllReturnCodes::ABORTRUN;
1043 }
1044
1045 if (rs->next() == 0)
1046 {
1047 std::cout << std::format("Error: Can't find GL1 scaledown data for run {}", runnumber) << std::endl;
1048 delete rs;
1049 delete stmt;
1050 delete con;
1051 return Fun4AllReturnCodes::ABORTRUN;
1052 }
1053
1054 int ncols = rs->getMetaData()->getColumnCount();
1055 for (int icol = 0; icol < ncols; icol++)
1056 {
1057 std::string col_name = rs->getMetaData()->getColumnName(icol + 1);
1058
1059 if (col_name.starts_with("scaledown"))
1060 {
1061 int scaledown_index = std::stoi(col_name.substr(9,2));
1062 std::string cont = rs->getString(col_name);
1063 scaledown[scaledown_index] = std::stoi(cont);
1064 }
1065 }
1066
1067
1068
1069 for (int index = 0; index < 64; index++)
1070 {
1071 cmd.str("");
1072 cmd << "SELECT raw, live FROM gl1_scalers WHERE runnumber=" << runnumber << " AND index=" << index << ";";
1073 stmt = nullptr; stmt = con->createStatement();
1074 rs = nullptr;
1075 try
1076 {
1077 rs = stmt->executeQuery(cmd.str());
1078 }
1079 catch (odbc::SQLException &e)
1080 {
1081 std::cout << std::format("Error: {}.", e.getMessage()) << std::endl;
1082 delete rs;
1083 delete stmt;
1084 delete con;
1085 return Fun4AllReturnCodes::ABORTRUN;
1086 }
1087
1088 if (rs->next() == 0)
1089 {
1090 delete rs;
1091 delete stmt;
1092 livetime[index] = 0;
1093 continue;
1094 }
1095
1096 double total_live = std::stod(rs->getString("live"));
1097 double total_raw = std::stod(rs->getString("raw"));
1098 livetime[index] = total_live / total_raw;
1099 }
1100
1101
1102
1103
1104
1105 SpinDBOutput spin_out("phnxrc");
1106 SpinDBContent *spin_cont = new SpinDBContentv1();
1107 spin_out.StoreDBContent(runnumber, runnumber);
1108 spin_out.GetDBContentStore(spin_cont, runnumber);
1109
1110 crossingshift = spin_cont->GetCrossingShift();
1111
1112 for (int i = 0; i < nBunches; i++)
1113 {
1114 beamspinpat[0][i] = -1 * spin_cont->GetSpinPatternYellow(i);
1115 beamspinpat[1][i] = -1 * spin_cont->GetSpinPatternBlue(i);
1116
1117
1118
1119 }
1120
1121 _eventcounter = -1;
1122 return Fun4AllReturnCodes::EVENT_OK;
1123 }
1124
1125 int AnNeutralMeson_micro_dst::process_event(PHCompositeNode *topNode)
1126 {
1127 _eventcounter++;
1128
1129 if (_eventcounter % 100'000 == 0)
1130 {
1131 monitorMemoryUsage("check");
1132 std::cout << "event: " << _eventcounter << std::endl;
1133 }
1134
1135 // Store event-level entries
1136 live_trigger = _smallclusters->get_live_trigger();
1137 scaled_trigger = _smallclusters->get_scaled_trigger();
1138 diphoton_bunchnumber = _smallclusters->get_bunch_number();
1139 cluster_number = _smallclusters->size();
1140
1141 if (vertexmap && !vertexmap->empty())
1142 {
1143 GlobalVertex *vtx = vertexmap->begin()->second;
1144 if (vtx)
1145 {
1146 // vtx_x and vtx_y always estimated to 0
1147 vertex_z = vtx->get_z();
1148 }
1149 else
1150 {
1151 return Fun4AllReturnCodes::ABORTEVENT;
1152 }
1153 }
1154 else
1155 {
1156 return Fun4AllReturnCodes::ABORTEVENT;
1157 }
1158
1159 if (store_qa)
1160 {
1161 // Trigger count per bit (QA)
1162 for (const auto &[trigger_number,
1163 trigger_name] : trigger_names)
1164 {
1165 if ((live_trigger >> trigger_number & 0x1U) == 0x1U)
1166 {
1167 h_trigger_live->Fill(trigger_number + 1);
1168 }
1169 if ((scaled_trigger >> trigger_number & 0x1U) == 0x1U)
1170 {
1171 h_trigger_scaled->Fill(trigger_number + 1);
1172 }
1173 }
1174 }
1175
1176 // Vertex cut
1177 if (std::abs(vertex_z) > vertex_max)
1178 {
1179 return Fun4AllReturnCodes::ABORTEVENT;
1180 }
1181 diphoton_vertex_z = vertex_z;
1182
1183 if (store_qa) h_event_zvtx->Fill(vertex_z);
1184
1185 // Additional trigger selection
1186 if (require_mbd_trigger_bit && !require_photon_trigger_bit)
1187 {
1188 // Minimum Bias Event
1189 if (!mbd_trigger_bit())
1190 {
1191 return Fun4AllReturnCodes::ABORTEVENT;
1192 }
1193 mbd_trigger_bit_event = true;
1194 }
1195 else if (require_photon_trigger_bit && !require_mbd_trigger_bit)
1196 {
1197 // Photon Trigger Event
1198 if (!photon_trigger_bit())
1199 {
1200 return Fun4AllReturnCodes::ABORTEVENT;
1201 }
1202 photon_trigger_bit_event = true;
1203 }
1204 else { // If both triggers are used, apply a pT threshold between them
1205 mbd_trigger_bit_event = mbd_trigger_bit();
1206 photon_trigger_bit_event = photon_trigger_bit();
1207 }
1208
1209 if (store_qa) {
1210 h_triggered_event_photons->Fill(cluster_number);
1211
1212 // Fill vertex info, for QA
1213 h_triggered_event_zvtx->Fill(vertex_z);
1214 }
1215
1216 // If more than one cut for ecore and chi2 is given, apply the following method for cluster number comparison (QA only)
1217 if (n_chi2_cuts > 1 || n_ecore_cuts > 1)
1218 {
1219 cluster_cuts();
1220 return Fun4AllReturnCodes::EVENT_OK;
1221 }
1222
1223 // Select the "good" photon clusters
1224 good_photons.clear();
1225 num_photons = 0;
1226 for (int iclus = 0; iclus < cluster_number; iclus++)
1227 {
1228 ClusterSmallInfo *smallcluster = _smallclusters->get_cluster_at(iclus);
1229 float eta = smallcluster->get_eta();
1230 float phi = smallcluster->get_phi();
1231 float chi2 = smallcluster->get_chi2();
1232 float ecore = smallcluster->get_ecore();
1233
1234 if (chi2 <= chi2_cuts[0] &&
1235 ecore >= ecore_cuts[0])
1236 {
1237 ROOT::Math::PtEtaPhiMVector photon(ecore / std::cosh(eta), eta, phi, 0);
1238 if (store_qa)
1239 {
1240 h_photon_eta->Fill(photon.Eta());
1241 h_photon_phi->Fill(photon.Phi());
1242 h_photon_pt->Fill(photon.Pt());
1243 h_photon_zvtx->Fill(vertex_z);
1244 h_photon_eta_phi->Fill(photon.Eta(), photon.Phi());
1245 h_photon_eta_pt->Fill(photon.Eta(), photon.Pt());
1246 h_photon_eta_zvtx->Fill(photon.Eta(), vertex_z);
1247 h_photon_phi_pt->Fill(photon.Phi(), photon.Pt());
1248 h_photon_phi_zvtx->Fill(photon.Phi(), vertex_z);
1249 h_photon_pt_zvtx->Fill(photon.Pt(), vertex_z);
1250 }
1251 Cluster cluster(photon, false); // by default, a cluster does not activate the trigger
1252 good_photons.push_back(cluster);
1253 num_photons++;
1254 }
1255 }
1256
1257 // Trigger emulator -> determine what cluster(s) fired the trigger (cluster.isTrigger = true)
1258 if (require_emulator_matching)
1259 {
1260 trigger_emulator_input();
1261 }
1262
1263 // Select the photon pairs
1264 int multiplicity_efficiency = 0;
1265 int multiplicity_emulator = 0;
1266 bool first_diphoton_nomatch = true;
1267 for (int iclus1 = 0; iclus1 < num_photons - 1; iclus1++)
1268 {
1269 ROOT::Math::PtEtaPhiMVector photon1 = good_photons[iclus1].p4;
1270 for (int iclus2 = iclus1 + 1; iclus2 < num_photons; iclus2++)
1271 {
1272 ROOT::Math::PtEtaPhiMVector photon2 = good_photons[iclus2].p4;
1273
1274 ROOT::Math::PtEtaPhiMVector diphoton = photon1 + photon2;
1275
1276 diphoton_pt = diphoton.Pt();
1277
1278 float eta_diff = photon1.Eta() - photon2.Eta();
1279 float phi_diff = WrapAngleDifference(photon1.Phi(), photon2.Phi());
1280 if (store_qa) h_pair_DeltaR->Fill(std::sqrt(std::pow(eta_diff, 2) + std::pow(phi_diff, 2)));
1281
1282 // General diphoton cut
1283 if (diphoton_cut(photon1, photon2, diphoton)) continue;
1284
1285 if (store_qa && diphoton.M() > 0.25 && diphoton.M() < 0.38)
1286 {
1287 h_same_delta_eta->Fill(std::abs(photon1.Eta() - photon2.Eta()));
1288 h_same_delta_phi->Fill(std::abs(ROOT::Math::VectorUtil::DeltaPhi(photon1, photon2)));
1289 h_same_delta_R->Fill(ROOT::Math::VectorUtil::DeltaR(photon1, photon2));
1290 h_same_alpha->Fill(std::abs(photon1.E() - photon2.E())/(photon1.E() + photon2.E()));
1291 }
1292
1293 if (store_qa && first_diphoton_nomatch)
1294 {
1295 h_count_diphoton_nomatch->Fill(1);
1296 first_diphoton_nomatch = false;
1297 }
1298
1299 // Distinct diphoton selection between MBD and Photon Trigger
1300
1301 // For a photon-triggered event
1302 // Keep candidates at high-pT which satisfy the trigger matching criterion
1303
1304 // Trigger emulator matching cut
1305 if (require_emulator_matching && !require_efficiency_matching)
1306 {
1307 emulator_match = (good_photons[iclus1].isTrigger || good_photons[iclus2].isTrigger);
1308 if (!(emulator_match))
1309 {
1310 continue;
1311 }
1312 multiplicity_emulator++;
1313 }
1314
1315 // Trigger efficiency matching cut ("Proxy method")
1316 if (require_efficiency_matching && !require_emulator_matching)
1317 {
1318 efficiency_match = trigger_efficiency_matching(photon1, photon2, diphoton);
1319 if (!(efficiency_match))
1320 {
1321 continue;
1322 }
1323 }
1324
1325 // Emulator vs proxy comparison
1326 if (require_efficiency_matching && require_emulator_matching)
1327 {
1328 int iPt = FindBinBinary(diphoton_pt, pTBins, nPtBins + 1);
1329 if (iPt < 0 || iPt >= nPtBins) continue;
1330 emulator_match = (good_photons[iclus1].isTrigger || good_photons[iclus2].isTrigger);
1331 efficiency_match = trigger_efficiency_matching(photon1, photon2, diphoton);
1332 if (emulator_match && efficiency_match)
1333 {
1334 h_efficiency_x_matching_pt[iPt]->AddBinContent(1);
1335 h_efficiency_x_matching->AddBinContent(1);
1336 }
1337 if (!emulator_match && !efficiency_match)
1338 {
1339 h_efficiency_x_matching_pt[iPt]->AddBinContent(2);
1340 h_efficiency_x_matching->AddBinContent(2);
1341 }
1342 if (emulator_match && !efficiency_match)
1343 {
1344 h_efficiency_x_matching_pt[iPt]->AddBinContent(3);
1345 h_efficiency_x_matching->AddBinContent(3);
1346 }
1347 if (!emulator_match && efficiency_match) {
1348 h_efficiency_x_matching_pt[iPt]->AddBinContent(4);
1349 h_efficiency_x_matching->AddBinContent(4);
1350 }
1351 continue;
1352 }
1353
1354 if (store_qa)
1355 {
1356 h_pair_alpha->Fill(std::abs(photon1.E() - photon2.E())/(photon1.E() + photon2.E()));
1357 h_pair_alpha_pt->Fill(std::abs(photon1.E() - photon2.E())/(photon1.E() + photon2.E()), diphoton_pt);
1358 }
1359
1360 if (multiplicity_efficiency == 1) {
1361 h_analysis_event_zvtx->Fill(vertex_z);
1362 }
1363
1364 // Store diphoton properties
1365 diphoton_mass = diphoton.mag();
1366 diphoton_eta = diphoton.Eta();
1367 diphoton_phi = diphoton.Phi();
1368 diphoton_xf = 2 * diphoton.Pz() / Vs;
1369
1370 if (store_qa)
1371 {
1372 h_selected_photon_eta->Fill(photon1.Eta());
1373 h_selected_photon_phi->Fill(photon1.Phi());
1374 h_selected_photon_pt->Fill(photon1.Pt());
1375 h_selected_photon_zvtx->Fill(vertex_z);
1376 h_selected_photon_eta_phi->Fill(photon1.Eta(), photon1.Phi());
1377 h_selected_photon_eta_pt->Fill(photon1.Eta(), photon1.Pt());
1378 h_selected_photon_eta_zvtx->Fill(photon1.Eta(), vertex_z);
1379 h_selected_photon_phi_pt->Fill(photon1.Phi(), photon1.Pt());
1380 h_selected_photon_phi_zvtx->Fill(photon1.Phi(), vertex_z);
1381 h_selected_photon_pt_zvtx->Fill(photon1.Pt(), vertex_z);
1382 h_selected_photon_eta->Fill(photon2.Eta());
1383 h_selected_photon_phi->Fill(photon2.Phi());
1384 h_selected_photon_pt->Fill(photon2.Pt());
1385 h_selected_photon_zvtx->Fill(vertex_z);
1386 h_selected_photon_eta_phi->Fill(photon2.Eta(), photon2.Phi());
1387 h_selected_photon_eta_pt->Fill(photon2.Eta(), photon2.Pt());
1388 h_selected_photon_eta_zvtx->Fill(photon2.Eta(), vertex_z);
1389 h_selected_photon_phi_pt->Fill(photon2.Phi(), photon2.Pt());
1390 h_selected_photon_phi_zvtx->Fill(photon2.Phi(), vertex_z);
1391 h_selected_photon_pt_zvtx->Fill(photon2.Pt(), vertex_z);
1392
1393 h_pair_eta->Fill(diphoton_eta);
1394 h_pair_phi->Fill(diphoton_phi);
1395 h_pair_pt->Fill(diphoton_pt);
1396 h_pair_zvtx->Fill(vertex_z);
1397 h_pair_eta_phi->Fill(diphoton_eta, diphoton_phi);
1398 h_pair_eta_pt->Fill(diphoton_eta, diphoton_pt);
1399 h_pair_eta_zvtx->Fill(diphoton_eta, vertex_z);
1400 h_pair_phi_pt->Fill(diphoton_phi, diphoton_pt);
1401 h_pair_phi_zvtx->Fill(diphoton_phi, vertex_z);
1402 h_pair_pt_zvtx->Fill(diphoton_pt, vertex_z);
1403 }
1404
1405 // Select the pt bin:
1406 int iPt = FindBinBinary(diphoton_pt, pTBins, nPtBins + 1);
1407
1408 // Select the zvtx bin:
1409 int izvtx = FindBinBinary(std::abs(vertex_z), zvtxBins, nZvtxBins + 1);
1410
1411 // Select the eta bin:
1412 int ietaBin = FindBinBinary(diphoton_eta, etaBins, nEtaBins + 1);
1413
1414 // Select the xf bin:
1415 int ixfBin = FindBinBinary(diphoton_xf, xfBins, nXfBins + 1);
1416
1417 // Select the eta region (small eta vs high eta)
1418 int ietaRegion = (std::abs(diphoton_eta) < etaThreshold ? 0 : 1);
1419
1420 // Invariant mass distributions
1421 h_pair_mass->Fill(diphoton_mass);
1422 if (izvtx >= 0 && izvtx < nZvtxBins)
1423 for (int izvtxtmp = izvtx; izvtxtmp < nZvtxBins; izvtxtmp++)
1424 h_pair_mass_zvtx[izvtxtmp]->Fill(diphoton_mass);
1425 if (iPt >= 0 && iPt < nPtBins)
1426 h_pair_mass_pt[iPt]->Fill(diphoton_mass);
1427 if (ietaBin >= 0 && ietaBin < nEtaBins)
1428 h_pair_mass_eta[ietaBin]->Fill(diphoton_mass);
1429 if (ixfBin >= 0 && ixfBin < nXfBins)
1430 h_pair_mass_xf[ixfBin]->Fill(diphoton_mass);
1431
1432 // Fill minimal output ttree
1433 // Can be used for future short analysis (see AnNeutralMeson_nano)
1434 if (store_tree)
1435 {
1436 output_tree->Fill();
1437 }
1438
1439 // Select particle and region index;
1440 int iP = -1;
1441 int iR = -1;
1442 int ival = FindBinBinary(diphoton_mass, band_limits,
1443 nParticles * (nRegions + 1) * 2);
1444 // Ignore any diphoton which is not within the side band or the peak band
1445 if (ival >= 0 && ival % 2 == 0)
1446 {
1447 iP = ival / 2 / (nRegions + 1); // particle index
1448 iR = (ival / 2 % (nRegions + 1) + 1) % 2; // region index
1449 }
1450
1451 // Store pT, eta and xF values in order to compute the average value of
1452 // each bin
1453 if (iP != -1)
1454 {
1455 if (iPt >= 0 && iPt < nPtBins)
1456 {
1457 h_average_pt[iP]->Fill(iPt, diphoton_pt);
1458 h_norm_pt[iP]->Fill(iPt, 1);
1459 }
1460 if (ietaBin >= 0 && ietaBin < nEtaBins)
1461 {
1462 h_average_eta[iP]->Fill(ietaBin, diphoton_eta);
1463 h_norm_eta[iP]->Fill(ietaBin, 1);
1464 }
1465 if (ixfBin >= 0 && ixfBin < nXfBins)
1466 {
1467 h_average_xf[iP]->Fill(ixfBin, diphoton_xf);
1468 h_norm_xf[iP]->Fill(ixfBin, 1);
1469 }
1470 }
1471
1472 // Get Kinematic distributions:
1473 if (store_qa)
1474 {
1475 if (iP == 0 && iR == 0)
1476 {
1477 h_meson_pi0_pt->Fill(diphoton_pt);
1478 h_meson_pi0_E->Fill(diphoton.E());
1479 if (1 < diphoton_pt && diphoton_pt < 2)
1480 {
1481 h_pair_pi0_eta_pt_1->Fill(diphoton_eta);
1482 h_pair_pi0_xf_pt_1->Fill(diphoton_xf);
1483 }
1484 else if (2 < diphoton_pt && diphoton_pt < 3)
1485 {
1486 h_pair_pi0_eta_pt_2->Fill(diphoton_eta);
1487 h_pair_pi0_xf_pt_2->Fill(diphoton_xf);
1488 }
1489 else if (3 < diphoton_pt && diphoton_pt < 4)
1490 {
1491 h_pair_pi0_eta_pt_3->Fill(diphoton_eta);
1492 h_pair_pi0_xf_pt_3->Fill(diphoton_xf);
1493 }
1494 else if (4 < diphoton_pt && diphoton_pt < 5)
1495 {
1496 h_pair_pi0_eta_pt_4->Fill(diphoton_eta);
1497 h_pair_pi0_xf_pt_4->Fill(diphoton_xf);
1498 }
1499 else if (5 < diphoton_pt && diphoton_pt < 6)
1500 {
1501 h_pair_pi0_eta_pt_5->Fill(diphoton_eta);
1502 h_pair_pi0_xf_pt_5->Fill(diphoton_xf);
1503 }
1504 else if (6 < diphoton_pt && diphoton_pt < 7)
1505 {
1506 h_pair_pi0_eta_pt_6->Fill(diphoton_eta);
1507 h_pair_pi0_xf_pt_6->Fill(diphoton_xf);
1508 }
1509 else if (7 < diphoton_pt && diphoton_pt < 8)
1510 {
1511 h_pair_pi0_eta_pt_7->Fill(diphoton_eta);
1512 h_pair_pi0_xf_pt_7->Fill(diphoton_xf);
1513 }
1514 else if (8 < diphoton_pt && diphoton_pt < 10)
1515 {
1516 h_pair_pi0_eta_pt_8->Fill(diphoton_eta);
1517 h_pair_pi0_xf_pt_8->Fill(diphoton_xf);
1518 }
1519 else if (10 < diphoton_pt && diphoton_pt < 20)
1520 {
1521 h_pair_pi0_eta_pt_9->Fill(diphoton_eta);
1522 h_pair_pi0_xf_pt_9->Fill(diphoton_xf);
1523 }
1524 }
1525
1526 if (iP == 1 && iR == 0)
1527 {
1528 h_meson_eta_pt->Fill(diphoton_pt);
1529 h_meson_eta_E->Fill(diphoton.E());
1530 if (1 < diphoton_pt && diphoton_pt < 2)
1531 {
1532 h_pair_eta_eta_pt_1->Fill(diphoton_eta);
1533 h_pair_eta_xf_pt_1->Fill(diphoton_xf);
1534 }
1535 else if (2 < diphoton_pt && diphoton_pt < 3)
1536 {
1537 h_pair_eta_eta_pt_2->Fill(diphoton_eta);
1538 h_pair_eta_xf_pt_2->Fill(diphoton_xf);
1539 }
1540 else if (3 < diphoton_pt && diphoton_pt < 4)
1541 {
1542 h_pair_eta_eta_pt_3->Fill(diphoton_eta);
1543 h_pair_eta_xf_pt_3->Fill(diphoton_xf);
1544 }
1545 else if (4 < diphoton_pt && diphoton_pt < 5)
1546 {
1547 h_pair_eta_eta_pt_4->Fill(diphoton_eta);
1548 h_pair_eta_xf_pt_4->Fill(diphoton_xf);
1549 }
1550 else if (5 < diphoton_pt && diphoton_pt < 6)
1551 {
1552 h_pair_eta_eta_pt_5->Fill(diphoton_eta);
1553 h_pair_eta_xf_pt_5->Fill(diphoton_xf);
1554 }
1555 else if (6 < diphoton_pt && diphoton_pt < 7)
1556 {
1557 h_pair_eta_eta_pt_6->Fill(diphoton_eta);
1558 h_pair_eta_xf_pt_6->Fill(diphoton_xf);
1559 }
1560 else if (7 < diphoton_pt && diphoton_pt < 8)
1561 {
1562 h_pair_eta_eta_pt_7->Fill(diphoton_eta);
1563 h_pair_eta_xf_pt_7->Fill(diphoton_xf);
1564 }
1565 else if (8 < diphoton_pt && diphoton_pt < 10)
1566 {
1567 h_pair_eta_eta_pt_8->Fill(diphoton_eta);
1568 h_pair_eta_xf_pt_8->Fill(diphoton_xf);
1569 }
1570 else if (10 < diphoton_pt && diphoton_pt < 20)
1571 {
1572 h_pair_eta_eta_pt_9->Fill(diphoton_eta);
1573 h_pair_eta_xf_pt_9->Fill(diphoton_xf);
1574 }
1575 }
1576 }
1577
1578 // Compute yield
1579 float phi_shift = 0;
1580 float phi_beam = 0;
1581 for (int iB = 0; iB < nBeams; iB++)
1582 {
1583 int iDir = (diphoton_eta * beamDirection[iB] > 0 ? 0 : 1);
1584 int ietaBinRelative =
1585 (beamDirection[iB] > 0 ? ietaBin : nEtaBins - 1 - ietaBin);
1586 int ixfBinRelative =
1587 (beamDirection[iB] > 0 ? ixfBin : nXfBins - 1 - ixfBin);
1588
1589 // sanity check (eta, xf > 0 -> forward)
1590 if (((ietaBin >= 0 && ietaBin < nEtaBins) && (ixfBin >= 0 && ixfBin < nXfBins)) &&
1591 (((iDir == 0) && ((etaBins[ietaBinRelative + 1] <= 0) ||
1592 (xfBins[ixfBinRelative + 1] <= 0))) ||
1593 ((iDir == 1) && ((etaBins[ietaBinRelative] >= 0) ||
1594 (xfBins[ixfBinRelative] >= 0)))))
1595 {
1596 std::cerr << "Error: Bad definition of relative eta, xf bins.\n";
1597 std::cerr << (iDir == 0 ? "forward"
1598 : "backward")
1599 << " -> ixf = " << ixfBinRelative << " ("
1600 << xfBins[ixfBinRelative] << ", "
1601 << xfBins[ixfBinRelative + 1]
1602 << "), ieta = " << ietaBinRelative << " ("
1603 << etaBins[ietaBinRelative] << ", "
1604 << etaBins[ietaBinRelative + 1] << ")\n";
1605 }
1606
1607 phi_shift = (iB == 0 ? M_PI / 2 : -M_PI / 2);
1608 phi_beam = diphoton_phi + phi_shift; // phi global to phi yellow/blue
1609 if (phi_beam < -M_PI)
1610 phi_beam += 2 * M_PI;
1611 else if (phi_beam > M_PI)
1612 phi_beam -= 2 * M_PI;
1613 if (iR != -1 && iP != -1)
1614 {
1615 if (beamspinpat[iB][(crossingshift + diphoton_bunchnumber) % nBunches] == 1)
1616 {
1617 if (iPt >= 0 && iPt < nPtBins)
1618 {
1619 h_yield_1[iB][iP][iR][iPt][ietaRegion][0]->Fill(phi_beam);
1620 h_yield_2[iB][iP][iR][iPt][iDir][0]->Fill(phi_beam);
1621 h_yield_3[iB][iP][iR][iPt][ietaRegion][iDir][0]->Fill(phi_beam);
1622 }
1623 if (ietaBin >= 0 && ietaBin < nEtaBins)
1624 {
1625 h_yield_eta[iB][iP][iR < 2 ? iR : 1][ietaBinRelative][0]->Fill(phi_beam);
1626 }
1627 if (ixfBin >= 0 && ixfBin < nXfBins)
1628 {
1629 h_yield_xf[iB][iP][iR < 2 ? iR : 1][ixfBinRelative][0]->Fill(phi_beam);
1630 }
1631 }
1632 else if (beamspinpat[iB][(crossingshift + diphoton_bunchnumber) % nBunches] == -1)
1633 {
1634 if (iPt >= 0 && iPt < nPtBins)
1635 {
1636 h_yield_1[iB][iP][iR < 2 ? iR : 1][iPt][ietaRegion][1]->Fill(phi_beam);
1637 h_yield_2[iB][iP][iR < 2 ? iR : 1][iPt][iDir][1]->Fill(phi_beam);
1638 h_yield_3[iB][iP][iR < 2 ? iR : 1][iPt][ietaRegion][iDir][1]->Fill(phi_beam);
1639 }
1640 if (ietaBin >= 0 && ietaBin < nEtaBins)
1641 {
1642 h_yield_eta[iB][iP][iR < 2 ? iR : 1][ietaBinRelative][1]->Fill(phi_beam);
1643 }
1644 if (ixfBin >= 0 && ixfBin < nXfBins)
1645 {
1646 h_yield_xf[iB][iP][iR < 2 ? iR : 1][ixfBinRelative][1]->Fill(phi_beam);
1647 }
1648 }
1649 }
1650 }
1651 }
1652 }
1653
1654 if (store_qa)
1655 {
1656 h_multiplicity_efficiency->Fill(multiplicity_efficiency);
1657 h_multiplicity_emulator->Fill(multiplicity_emulator);
1658 }
1659
1660 // Event Mixing: If there is at least one good pair, store the event in the pool
1661 if (use_event_mixing)
1662 {
1663 if (require_mbd_trigger_bit)
1664 {
1665 event_mixing_mbd(topNode);
1666 }
1667 else if (require_photon_trigger_bit)
1668 {
1669 event_mixing_photon();
1670 }
1671 }
1672
1673 return Fun4AllReturnCodes::EVENT_OK;
1674 }
1675
1676 int AnNeutralMeson_micro_dst::End(PHCompositeNode *)
1677 {
1678 outfile->cd();
1679 outfile->Write();
1680 outfile->Close();
1681 delete outfile;
1682 outfile = nullptr; // To prevent double deletion with the destructor
1683 if (store_tree)
1684 {
1685 outfile_tree->cd();
1686 outfile_tree->Write();
1687 outfile_tree->Close();
1688 delete outfile_tree;
1689 outfile_tree = nullptr;
1690 }
1691 return 0;
1692 }
1693
1694 void AnNeutralMeson_micro_dst::cluster_cuts()
1695 {
1696 // Define the good clusters with respect to the chi2 and energy (corrected) cuts
1697 // set before the execution
1698 std::vector<std::vector<ROOT::Math::PtEtaPhiMVector>> good_photons_by_cut(n_chi2_cuts * n_ecore_cuts);
1699 for (int iclus = 0; iclus < cluster_number; iclus++)
1700 {
1701 const ClusterSmallInfo *smallcluster = _smallclusters->get_cluster_at(iclus);
1702 float eta = smallcluster->get_eta();
1703 float phi = smallcluster->get_phi();
1704 float chi2 = smallcluster->get_chi2();
1705 float ecore = smallcluster->get_ecore();
1706 int ichi2cut_max = FindBinBinary(chi2, chi2_cuts.data(), n_chi2_cuts);
1707 int iecorecut_max = FindBinBinary(ecore, ecore_cuts.data(), n_ecore_cuts);
1708 ROOT::Math::PtEtaPhiMVector photon(ecore / eta, eta, phi, 0);
1709 for (int ichi2cut = 0; ichi2cut <= ichi2cut_max; ichi2cut++)
1710 {
1711 for (int iecorecut = 0; iecorecut <= iecorecut_max; iecorecut++)
1712 {
1713 good_photons_by_cut[ichi2cut * n_ecore_cuts + iecorecut].push_back(photon);
1714 }
1715 }
1716 }
1717
1718 for (int ichi2cut = 0; ichi2cut < n_chi2_cuts; ichi2cut++)
1719 {
1720 for (int iecorecut = 0; iecorecut < n_ecore_cuts; iecorecut++)
1721 {
1722 num_photons = good_photons_by_cut[ichi2cut * n_ecore_cuts + iecorecut].size();
1723 for (int iclus1 = 0; iclus1 < num_photons; iclus1++)
1724 {
1725 ROOT::Math::PtEtaPhiMVector photon1 = good_photons_by_cut[ichi2cut * n_ecore_cuts + iecorecut][iclus1];
1726 for (int iclus2 = iclus1 + 1; iclus2 < num_photons; iclus2++)
1727 {
1728 ROOT::Math::PtEtaPhiMVector photon2 = good_photons_by_cut[ichi2cut * n_ecore_cuts + iecorecut][iclus2];
1729 ROOT::Math::PtEtaPhiMVector diphoton = photon1 + photon2;
1730
1731 if (diphoton_cut(photon1, photon2, diphoton)) continue;
1732
1733 diphoton_mass = diphoton.mag();
1734
1735 // Select particle and region index;
1736 int iP = -1;
1737 int iR = -1;
1738 int ival = FindBinBinary(diphoton_mass, band_limits,
1739 nParticles * (nRegions + 1) * 2);
1740 // Ignore any diphoton which is not within the side band or the peak band
1741 if (ival >= 0 && ival % 2 == 0)
1742 {
1743 iP = ival / 2 / (nRegions + 1); // particle index
1744 iR = (ival / 2 % (nRegions + 1) + 1) % 2; // region index
1745 }
1746
1747 int iPt = FindBinBinary(diphoton_pt, pTBins, nPtBins + 1);
1748
1749 if (iPt >= 0 && iPt < nPtBins)
1750 h_cluster_level_cuts_total_pt[iPt]
1751 ->Fill((float) ichi2cut + 0.5, (float) iecorecut + 0.5);
1752 // Check within pi0 and eta mass range
1753 if (iP != -1 && iR != -1)
1754 {
1755 // OK
1756 h_cluster_level_cuts_particle[iP][iR]
1757 ->Fill((float) ichi2cut + 0.5, (float) iecorecut + 0.5);
1758 if (iPt >= 0 && iPt < nPtBins)
1759 h_cluster_level_cuts_particle_pt[iP][iR][iPt]
1760 ->Fill((float) ichi2cut + 0.5, (float) iecorecut + 0.5);
1761 }
1762 }
1763 }
1764 }
1765 }
1766 }
1767
1768 bool AnNeutralMeson_micro_dst::diphoton_cut(ROOT::Math::PtEtaPhiMVector p1,
1769 ROOT::Math::PtEtaPhiMVector p2,
1770 ROOT::Math::PtEtaPhiMVector ppair)
1771 {
1772 float alpha = std::abs(p1.E() - p2.E()) / (p1.E() + p2.E());
1773 float pt = ppair.Pt();
1774 return (alpha > alphaCut || pt < pTCutMin || pt > pTCutMax);
1775 }
1776
1777 bool AnNeutralMeson_micro_dst::mbd_trigger_bit()
1778 {
1779 trigger_mbd_any_vtx = ((scaled_trigger >> 10 & 0x1U) == 0x1U);
1780 trigger_mbd_vtx_10 = ((scaled_trigger >> 12 & 0x1U) == 0x1U);
1781 trigger_mbd = (trigger_mbd_any_vtx || trigger_mbd_vtx_10);
1782
1783 return trigger_mbd;
1784 }
1785
1786 bool AnNeutralMeson_micro_dst::photon_trigger_bit()
1787 {
1788 trigger_mbd_photon_3 = false;
1789 trigger_mbd_photon_4 = false;
1790
1791 // If (photon 3 GeV + MBD NS >= 1) is enabled for recording
1792 if ((scaledown[25] != -1) &&
1793 ((scaled_trigger >> 25 & 0x1U) == 0x1U))
1794 {
1795 trigger_mbd_photon_3 = true;
1796 }
1797 // Else if (photon 3 GeV + MBD NS >=1, vtx < 10) is enabled for recording
1798 if ((scaledown[36] != -1) &&
1799 ((scaled_trigger >> 36 & 0x1U) == 0x1U))
1800 {
1801 trigger_mbd_photon_3 = true;
1802 }
1803 // Else if (photon 3 GeV) is enabled for recording with no prescale
1804 if ((scaledown[29] == 0) &&
1805 (((live_trigger >> 25 & 0x1U) == 0x1U) ||
1806 ((live_trigger >> 36 & 0x1U) == 0x1U)))
1807 {
1808 if ((live_trigger >> 25 & 0x1U) == 0x1U)
1809 {
1810 trigger_mbd_photon_3 = true;
1811 }
1812 if ((live_trigger >> 36 & 0x1U) == 0x1U)
1813 {
1814 trigger_mbd_photon_3 = true;
1815 }
1816 }
1817
1818 // If (photon 4 GeV + MBD NS >= 1) is enabled for recording
1819 if ((scaledown[26] != -1) &&
1820 ((scaled_trigger >> 26 & 0x1U) == 0x1U))
1821 {
1822 trigger_mbd_photon_4 = true;
1823 }
1824 // Else if (photon 4 GeV + MBD NS >=1, vtx < 10) is enabled for recording
1825 if ((scaledown[37] != -1) &&
1826 ((scaled_trigger >> 37 & 0x1U) == 0x1U))
1827 {
1828 trigger_mbd_photon_4 = true;
1829 }
1830 // Else if (photon 4 GeV) is enabled for recording with no prescale
1831 if ((scaledown[30] == 0) &&
1832 (((live_trigger >> 26 & 0x1U) == 0x1U) ||
1833 ((live_trigger >> 37 & 0x1U) == 0x1U)))
1834 {
1835 if ((live_trigger >> 26 & 0x1U) == 0x1U)
1836 {
1837 trigger_mbd_photon_4 = true;
1838 }
1839 if ((live_trigger >> 37 & 0x1U) == 0x1U)
1840 {
1841 trigger_mbd_photon_4 = true;
1842 }
1843 }
1844 return (trigger_mbd_photon_3 || trigger_mbd_photon_4);
1845 }
1846
1847 void AnNeutralMeson_micro_dst::trigger_emulator_input()
1848 {
1849 // 1) Identify the tile(s) (8x8 tower windows) which fired the trigger
1850 // And if so, to what energy: 3 GeV or 4 GeV
1851 // In most cases, only one tile fired
1852 emulator_selection = 0;
1853 fired_indices.clear();
1854
1855 tileNumber = emcaltiles->get_tile_number();
1856 if (trigger_mbd_photon_3)
1857 {
1858 for (int itile = 0; itile < tileNumber; itile++)
1859 {
1860 if (emcaltiles->get_tile_energy_adc(itile) >= adc_threshold_3)
1861 {
1862 fired_indices.push_back(itile);
1863 emulator_selection++;
1864 }
1865 }
1866 }
1867 else if (trigger_mbd_photon_4)
1868 {
1869 for (int itile = 0; itile < tileNumber; itile++)
1870 {
1871 if (emcaltiles->get_tile_energy_adc(itile) >= adc_threshold_4)
1872 {
1873 fired_indices.push_back(itile);
1874 emulator_selection++;
1875 }
1876 }
1877 }
1878 h_emulator_selection->Fill(emulator_selection);
1879
1880 // Identify all the good photon clusters
1881 // located in the firing tile
1882 // And store them by decreasing energy
1883 emulator_energies.clear();
1884 emulator_clusters.clear();
1885 if (emulator_selection) // i.e. emulator detects at least one firing tile
1886 {
1887 for (int i = 0; i < emulator_selection; i++)
1888 {
1889 // Tile dimensions
1890 int itile = fired_indices[i];
1891 int ieta_min = emcaltiles->get_tile_eta(itile) * 8;
1892 int ieta_max = (emcaltiles->get_tile_eta(itile) + 1) * 8 - 1;
1893 int iphi_min = emcaltiles->get_tile_phi(itile) * 8;
1894 int iphi_max = (emcaltiles->get_tile_phi(itile) + 1) * 8 - 1;
1895 double eta_min = towergeom->get_etabounds(std::max(0,ieta_min)).first;
1896 double eta_max = towergeom->get_etabounds(std::min(ieta_max,95)).second;
1897 double z_min = radius * std::sinh(eta_min);
1898 double z_max = radius * std::sinh(eta_max);
1899 double phi_min = towergeom->get_phibounds(std::max(0,iphi_min)).first;
1900 double phi_max = towergeom->get_phibounds(std::min(255,iphi_max)).second;
1901 phi_min = (phi_min >= 0 ? phi_min : phi_min + 2 * M_PI);
1902 phi_max = (phi_max >= 0 ? phi_max : phi_max + 2 * M_PI);
1903 for (int iclus1 = 0; iclus1 < num_photons; iclus1++)
1904 {
1905 ROOT::Math::PtEtaPhiMVector photon1 = good_photons[iclus1].p4;
1906 double photon1_z = vertex_z + radius * std::sinh(photon1.Eta());
1907 double photon1_phi = (photon1.Phi() >= 0 ? photon1.Phi() : photon1.Phi() + 2 * M_PI);
1908
1909 if ((photon1_z >= z_min && photon1_z <= z_max) &&
1910 ((phi_max >= phi_min && photon1_phi >= phi_min && photon1_phi <= phi_max) ||
1911 (phi_max <= phi_min && (photon1_phi <= phi_max || photon1_phi >= phi_min)))
1912 )
1913 {
1914 emulator_energies[itile].push_back(photon1.E());
1915 emulator_clusters[itile].push_back(iclus1);
1916 }
1917 }
1918 if (auto it = emulator_energies.find(itile); it != emulator_energies.end() && !it->second.empty())
1919 {
1920 emulator_multiplicities[itile] = emulator_energies[itile].size();
1921
1922 // Sort the clusters by decreasing energy
1923 std::vector<int> indices(emulator_energies[itile].size());
1924 std::iota(indices.begin(), indices.end(), 0);
1925 std::sort(indices.begin(), indices.end(), [&](int a, int b) -> bool { return emulator_energies[itile][a] > emulator_energies[itile][b]; });
1926 std::vector<int> vec_temp_clusters(emulator_clusters[itile]);
1927 std::vector<float> vec_temp_energies(emulator_energies[itile]);
1928 for (int iclus = 0; iclus < (int) emulator_energies[itile].size(); iclus++)
1929 {
1930 emulator_clusters[itile][iclus] = vec_temp_clusters[indices[iclus]];
1931 emulator_energies[itile][iclus] = vec_temp_energies[indices[iclus]];
1932 }
1933 vec_temp_clusters.clear();
1934 vec_temp_energies.clear();
1935 }
1936 }
1937 }
1938
1939 // Identify the trigger firing clusters:
1940 // Clusters which account for more than 20% of the tile energy
1941 // Usually, there is only 1 per event
1942 std::map<int, std::vector<int>>::const_iterator it_trigger;
1943 std::map<int, std::vector<float>>::const_iterator it_energy;
1944 it_trigger = emulator_clusters.begin();
1945 it_energy = emulator_energies.begin();
1946 for (; it_trigger != emulator_clusters.end(); it_trigger++, it_energy++)
1947 {
1948 std::vector<float>::const_iterator it_local_energy;
1949 float trigger_window_sum = std::accumulate(it_energy->second.begin(), it_energy->second.end(), 0.0);
1950 float emulator_iclus = 0;
1951 float clus_energy = trigger_window_sum;
1952 for (size_t local_index = 0; local_index < it_trigger->second.size(); local_index++)
1953 {
1954 emulator_iclus = it_trigger->second[local_index];
1955 clus_energy = it_energy->second[local_index];
1956 if (clus_energy < 0.2 * trigger_window_sum)
1957 {
1958 break;
1959 }
1960 good_photons[emulator_iclus].isTrigger = true;
1961 }
1962 }
1963 }
1964
1965 bool AnNeutralMeson_micro_dst::trigger_efficiency_matching(const ROOT::Math::PtEtaPhiMVector& photon1,
1966 const ROOT::Math::PtEtaPhiMVector& photon2,
1967 const ROOT::Math::PtEtaPhiMVector& diphoton)
1968
1969 {
1970 float energy_threshold = 0;
1971
1972 if (trigger_mbd_photon_3) energy_threshold = energy_threshold_3; // Photon 3 GeV trigger
1973 else if (trigger_mbd_photon_4) energy_threshold = energy_threshold_4; // Photon 4 GeV trigger
1974
1975 else
1976 {
1977 std::cerr << "Error: either photon 3 GeV or photon 4 GeV scaled trigger bit should be fired at this point.\n";
1978 exit(1);
1979 }
1980
1981 // If one of the two clusters has sufficient energy to fire the specific trigger
1982 // that led to the event written on disk, keep the event
1983 if (!(photon1.E() >= energy_threshold || photon2.E() >= energy_threshold))
1984 {
1985 // If no, check if:
1986 // - distance between two clusters is small enough for them to have been in the same trigger tile
1987 // - if their energy sum is sufficient to fire the specific trigger.
1988 float delta_eta = std::abs(photon1.Eta() - photon2.Eta());
1989 float delta_phi = WrapAngleDifference(photon1.Phi(), photon2.Phi());
1990 if (!((diphoton.E() > energy_threshold) &&
1991 (delta_eta < delta_eta_threshold) &&
1992 (delta_phi < delta_phi_threshold)))
1993 {
1994 // Otherwise, discard the diphoton
1995 return false;
1996 }
1997 }
1998 return true;
1999 }
2000
2001 void AnNeutralMeson_micro_dst::event_mixing_mbd(PHCompositeNode *topNode)
2002 {
2003
2004 EventInPool current_event;
2005 current_event.zvtx = vertex_z;
2006 current_event.clusters = good_photons;
2007
2008 // Associate cluster from current event to
2009 // clusters in pooled events
2010 int iPoolZvtxBin = FindBinBinary(vertex_z, poolZvtxBins, nPoolZvtxBins + 1);
2011 int iMultiplicityBin = FindBinBinary(num_photons, multiplicityBins, nMultiplicityBins + 1);
2012
2013 // Determine the angle of the leading jet
2014 std::vector<ROOT::Math::PtEtaPhiMVector> photon_vectors;
2015 for (int iclus = 0; iclus < cluster_number; iclus++)
2016 {
2017 ClusterSmallInfo *smallcluster = _smallclusters->get_cluster_at(iclus);
2018 float eta = smallcluster->get_eta();
2019 float phi = smallcluster->get_phi();
2020 float ecore = smallcluster->get_ecore();
2021 ROOT::Math::PtEtaPhiMVector photon(ecore / std::cosh(eta), eta, phi, 0);
2022 photon_vectors.push_back(photon);
2023 }
2024
2025 // Determine the leading jet position
2026 JetContainer *jet_cont = findNode::getClass<JetContainer>(topNode, "AntiKt_Tower_r04");
2027
2028 if (!jet_cont)
2029 {
2030 std::cout << "Jet Container is empty" << std::endl;
2031 return;
2032 }
2033
2034 Jet *leading_jet = nullptr;
2035 float leading_pt = 0;
2036 float sum_relative_pt = 0;
2037 float sum_absolute_pt = 0;
2038 for (Jet *jet : *jet_cont)
2039 {
2040 std::cout << "(" << jet->get_pt() << ", " << jet->get_eta() << ", " << jet->get_phi() << ")" << std::endl;
2041 sum_relative_pt += jet->get_pt();
2042 sum_absolute_pt += std::abs(jet->get_pt());
2043 if (jet->get_pt() > leading_pt)
2044 {
2045 leading_jet = jet;
2046 }
2047 }
2048 std::cout << std::endl;
2049
2050 // Sharpness (to define "jetty" events):
2051 float S = std::abs(sum_relative_pt) / sum_absolute_pt;
2052 bool is_jet_event = false;
2053 if (S > 0.3) is_jet_event = true;
2054
2055 std::cout << "Sharpness = " << S << std::endl;
2056
2057 if (is_jet_event)
2058 {
2059 std::cout << "Jetty event (" << leading_jet->get_eta() << ", " << leading_jet->get_phi() << ")" << std::endl;
2060 current_event.eta_lead = leading_jet->get_eta();
2061 current_event.phi_lead = leading_jet->get_phi();
2062 }
2063 else
2064 {
2065 std::cout << "Not a jetty event" << std::endl;
2066 }
2067
2068 int iPoolEtaBin = FindBinBinary(current_event.eta_lead, poolEtaBins, nPoolEtaBins + 1);
2069
2070 if ((iPoolZvtxBin < 0 || iPoolZvtxBin >= nPoolZvtxBins) ||
2071 (iMultiplicityBin < 0 || iMultiplicityBin >= nMultiplicityBins) ||
2072 (iPoolEtaBin < 0 || iPoolEtaBin >= nPoolEtaBins))
2073 {
2074 return;
2075 }
2076
2077 int CurrentNmix = std::min(int(pool[iPoolZvtxBin][iMultiplicityBin][iPoolEtaBin].size()), Nmix);
2078 float sum_multiplicities_pool = 0;
2079 for (int ipool = 0; ipool < CurrentNmix; ipool++)
2080 {
2081 EventInPool pool_event = pool[iPoolZvtxBin][iMultiplicityBin][iPoolEtaBin][ipool];
2082 sum_multiplicities_pool += (float) pool_event.clusters.size();
2083 }
2084
2085 float correction_weight = (sum_multiplicities_pool > 0.5 ? ((float) current_event.clusters.size() - 1)/ sum_multiplicities_pool : 0);
2086
2087 for (auto current_cluster: current_event.clusters)
2088 {
2089 for (int ipool = 0; ipool < CurrentNmix; ipool++)
2090 {
2091 EventInPool pool_event = pool[iPoolZvtxBin][iMultiplicityBin][iPoolEtaBin][ipool];
2092 for (auto pool_cluster: pool_event.clusters)
2093 {
2094 // Adapt pseudorapidity direction to the current vertex position
2095 float z_pool = pool_event.zvtx + radius * std::sinh(pool_cluster.p4.Eta());
2096 float eta_pool = std::asinh((z_pool - current_event.zvtx) / radius);
2097
2098 ROOT::Math::PtEtaPhiMVector pool_photon(0, 0, 0, 0);
2099 if (is_jet_event) // rotation
2100 {
2101 // Adapt pseudorapidity direction to the current vertex position
2102 float z_lead_pool = pool_event.zvtx + radius * std::sinh(pool_event.eta_lead);
2103 float eta_lead_pool = std::asinh((z_lead_pool - current_event.zvtx) / radius);
2104
2105 // Rotate the event to align the jet axes
2106 float eta_rot = eta_pool + current_event.eta_lead - eta_lead_pool;
2107 float phi_rot = pool_cluster.p4.Phi() + current_event.phi_lead - pool_event.phi_lead;
2108 pool_photon.SetCoordinates(pool_cluster.p4.Pt(), eta_rot, phi_rot, 0);
2109 }
2110 else // Naïve event mixing
2111 {
2112 pool_photon.SetCoordinates(pool_cluster.p4.Pt(), eta_pool, pool_cluster.p4.Phi(), 0);
2113 }
2114
2115 ROOT::Math::PtEtaPhiMVector mixed_pair = current_cluster.p4 + pool_photon;
2116 if (false)
2117 mixed_pair = current_cluster.p4 + pool_photon;
2118
2119 if (store_qa)
2120 {
2121 h_current_E->Fill(current_cluster.p4.E());
2122 h_current_eta->Fill(current_cluster.p4.Eta());
2123 h_current_phi->Fill(current_cluster.p4.Phi());
2124 h_pool_E->Fill(pool_photon.E());
2125 h_pool_eta->Fill(pool_photon.Eta());
2126 h_pool_phi->Fill(pool_photon.Phi());
2127 }
2128
2129 // The mixed event pair should satisfy the main cut as main analysis
2130 if (diphoton_cut(current_cluster.p4, pool_photon, mixed_pair)) continue;
2131 // In addition, require a minimum separation between two clusters
2132 // In reconstruction, two clusters cannot be too close
2133 float dR = ROOT::Math::VectorUtil::DeltaR(current_cluster.p4, pool_photon);
2134 if (dR < dR_min) continue;
2135 if (dR < 0.05) continue;
2136
2137 if (store_qa && mixed_pair.M() > 0.25 && mixed_pair.M() < 0.38)
2138 {
2139 h_mixed_delta_eta->Fill(std::abs(current_cluster.p4.Eta() - pool_photon.Eta()), correction_weight);
2140 h_mixed_delta_phi->Fill(std::abs(ROOT::Math::VectorUtil::DeltaPhi(current_cluster.p4, pool_photon)), correction_weight);
2141 h_mixed_delta_R->Fill(dR, correction_weight);
2142 h_mixed_alpha->Fill(std::abs(current_cluster.p4.E() - pool_photon.E()) / (current_cluster.p4.E() + pool_photon.E()), correction_weight);
2143 }
2144
2145 int iPt = FindBinBinary(mixed_pair.Pt(), pTBins, nPtBins + 1);
2146
2147 float mixed_mass = mixed_pair.M();
2148
2149 h_pair_mass_mixing->Fill(mixed_mass, correction_weight);
2150 h_pair_mass_mixing_pt[iPt]->Fill(mixed_mass, correction_weight);
2151 }
2152 }
2153 }
2154 pool[iPoolZvtxBin][iMultiplicityBin][iPoolEtaBin].push_front(current_event);
2155 if (CurrentNmix == Nmix)
2156 {
2157 pool[iPoolZvtxBin][iMultiplicityBin][iPoolEtaBin].pop_back();
2158 }
2159 }
2160
2161 void AnNeutralMeson_micro_dst::event_mixing_photon()
2162 {
2163 // if ((std::abs(vertex_z) >= 50) || (cluster_number > 11))
2164 // return;
2165
2166 // EventInPool current_event;
2167 // current_event.zvtx = vertex_z;
2168 // current_event.clusters = good_photons;
2169
2170 // // Associate trigger cluster from current event to
2171 // // non-trigger clusters in pooled events
2172 // int iPoolZvtxBin = FindBinBinary(vertex_z, poolZvtxBins, nPoolZvtxBins + 1);
2173 // int iMultiplicityBin = FindBinBinary(cluster_number, multiplicityBins, nMultiplicityBins + 1);
2174 // int CurrentNmix = std::min(int(pool[iPoolZvtxBin][iMultiplicityBin].size()), Nmix);
2175 // for (auto current_cluster: current_event.clusters)
2176 // {
2177 // if (! current_cluster.isTrigger) continue;
2178 // for (int ipool = 0; ipool < CurrentNmix; ipool++)
2179 // {
2180 // EventInPool pool_event = pool[iPoolZvtxBin][iMultiplicityBin][ipool];
2181 // for (auto pool_cluster: pool_event.clusters)
2182 // {
2183 // if (pool_cluster.isTrigger) continue;
2184
2185 // ROOT::Math::PtEtaPhiMVector mixed_pair = current_cluster.p4 + pool_cluster.p4;
2186 // if (diphoton_cut(current_cluster.p4, pool_cluster.p4, mixed_pair)) continue;
2187 // int iPt = FindBinBinary(mixed_pair.Pt(), pTBins, nPtBins + 1);
2188
2189 // float mixed_mass = mixed_pair.M();
2190
2191 // h_pair_mass_mixing->Fill(mixed_mass);
2192 // h_pair_mass_mixing_pt[iPt]->Fill(mixed_mass);
2193 // }
2194 // }
2195 // }
2196 // pool[iPoolZvtxBin][iMultiplicityBin].push_front(current_event);
2197 // if (CurrentNmix == Nmix)
2198 // {
2199 // pool[iPoolZvtxBin][iMultiplicityBin].pop_back();
2200 // }
2201 }
2202
2203 bool AnNeutralMeson_micro_dst::startswith(const std::string& str, const std::string& cmp)
2204 {
2205 return str.compare(0, cmp.length(), cmp) == 0;
2206 }
2207
2208 float AnNeutralMeson_micro_dst::WrapAngleDifference(const float& phi1, const float& phi2)
2209 {
2210 float difference = std::fmod(phi1 - phi2, 2 * M_PI);
2211 if (difference < 0)
2212 {
2213 difference += 2 * M_PI;
2214 }
2215 return std::min(difference, (float)(2 * M_PI - difference));
2216 }
2217
2218 ROOT::Math::XYZVector AnNeutralMeson_micro_dst::EnergyWeightedAverageP3(
2219 const std::vector<ROOT::Math::PtEtaPhiMVector>& v4s)
2220 {
2221 ROOT::Math::XYZVector sum_vector(0.,0.,0.);
2222 double sumE = 0.0;
2223
2224 for (const auto& v : v4s) {
2225 const double E = v.E();
2226 sum_vector += E * ROOT::Math::XYZVector(v.Px(), v.Py(), v.Pz());
2227 sumE += E;
2228 }
2229
2230 return (sumE > 0.0) ? (sum_vector * (1.0 / sumE)) : ROOT::Math::XYZVector(0.,0.,0.);
2231 }
