Implement/testcode/analyze_SiSeedPair_pT.C

0001 #include <TFile.h>
0002 #include <TTree.h>
0003 #include <TH1.h>
0004 #include <TH2.h>
0005 #include <TNtuple.h>
0006
0007 #include <TLorentzVector.h>
0008
0009 #include <cmath>
0010 #include <string>
0011 #include <vector>
0012
0013 #include "PtCalculator.h"  // SiCaloPt::PtCalculator & friends
0014 R__LOAD_LIBRARY(libPtCalc.so)
0015
0016
0017 const float electron_mass = 0.000511; // GeV/c^2
0018 const float pion_mass     = 0.1396; // GeV/c^2
0019
0020 SiCaloPt::PtCalculator pTCalc;
0021
0022 // ---- Weights(onnx) and Scalers(json) Path ---------------------------
0023 struct DemoPaths
0024 {
0025   std::string file_dir             = "/sphenix/user/jzhang1/testcode4all/INTT-EMCAL/InttSeedingTrackDev/ML4Reco/Implement/ML_Weight_Scaler/";
0026   std::string emd_onnx             = file_dir + "model_MLEMD.onnx";
0027   std::string emd_scaler_json      = "";
0028
0029   std::string eproj_onnx           = file_dir + "model_MLEproj.onnx";
0030   std::string eproj_scaler_json    = file_dir + "scaler_MLEproj.json";
0031
0032   std::string combined_onnx        = file_dir + "model_MLCombined.onnx";
0033   std::string combined_scaler_json = "";
0034
0035   void print(){
0036     std::cout<<"emd_onnx             : "<<emd_onnx            <<std::endl;
0037     std::cout<<"emd_scaler_json      : "<<emd_scaler_json     <<std::endl;
0038     std::cout<<"eproj_onnx           : "<<eproj_onnx          <<std::endl;
0039     std::cout<<"eproj_scaler_json    : "<<eproj_scaler_json   <<std::endl;
0040     std::cout<<"combined_onnx        : "<<combined_onnx       <<std::endl;
0041     std::cout<<"combined_scaler_json : "<<combined_scaler_json<<std::endl;
0042   };
0043 };
0044
0045 // ---- turn string into optional<string> for Config ------------------------
0046 template<typename Opt>
0047 Opt make_opt(const std::string& s)
0048 {
0049   if (s.empty()) return std::nullopt;
0050   return s;
0051 }
0052
0053 struct PtInput {
0054   float SiIn[3];
0055   float SiOut[3];
0056   float Calo[3];
0057   float CaloE;
0058 };
0059
0060 float calcPt(PtInput& in)
0061 {
0062   // for EMD and ML-EMD
0063   float phi_si   = atan2(in.SiOut[1]-in.SiIn[1], in.SiOut[0]-in.SiIn[0]);
0064   float phi_calo = atan2(in.Calo[1]-in.SiOut[1], in.Calo[0]-in.SiOut[0]);
0065   float dphi = phi_calo - phi_si;
0066
0067   // for EMD
0068   SiCaloPt::InputEMD inEMD;
0069   inEMD.EMD_Angle  = dphi; // delta_phi - EM Deflection angle in rad
0070   inEMD.EMD_Eta    = 0.00; // EMCal Cluster eta
0071   inEMD.EMD_Radius = 93.5; // EMCal Cluster radius in cm
0072
0073   auto r = pTCalc.ComputePt(SiCaloPt::Method::MethodEMD, SiCaloPt::AnyInput{inEMD});
0074   std::cout << "[EMD-analytic] ok=" << r.ok
0075             << "  pt=" << r.pt_reco
0076             << "  err=\"" << r.err << "\"\n";
0077
0078   float pT = r.pt_reco;
0079
0080   // for ML-EMD
0081   // 2-d input features:{ dphi_EMD, eta_track }
0082   std::vector<float> featsMLEMD = {dphi, 0};
0083   SiCaloPt::InputMLEMD inMLEMD{featsMLEMD};
0084
0085   r = pTCalc.ComputePt(SiCaloPt::Method::MethodMLEMD, SiCaloPt::AnyInput{inMLEMD});
0086   std::cout << "[MLEMD-2D] ok=" << r.ok
0087           << "  pt=" << r.pt_reco
0088           << "  err=\"" << r.err << "\"\n";
0089
0090
0091   // for ML-Eproj
0092   float r_siin  = sqrt(pow(in.SiIn[0], 2) + pow(in.SiIn[1], 2));
0093   float r_siout = sqrt(pow(in.SiOut[0], 2)+ pow(in.SiOut[1], 2));
0094   float r_calo  = sqrt(pow(in.Calo[0], 2)+ pow(in.Calo[1], 2));
0095
0096
0097   // 7-d input features:{ INTT 3/4 layer R, INTT 3/4 layer Z, INTT 5/4 layer R, INTT 5/4 layer Z, Calo R, Calo Z, Calo Energy }
0098   std::vector<float> featsMLEproj = { r_siin,  in.SiIn[2],   // INTT 3/4 layer
0099                                       r_siout, in.SiOut[2],  // INTT 5/4 layer
0100                                       r_calo,  in.Calo[2], in.CaloE}; // Calo R,Z,Energy
0101
0102   SiCaloPt::InputMLEproj inMLEproj{featsMLEproj};
0103   r = pTCalc.ComputePt(SiCaloPt::Method::MethodMLEproj, SiCaloPt::AnyInput{inMLEproj});
0104   std::cout << "[MLEproj-7D] ok=" << r.ok
0105           << "  pt=" << r.pt_reco
0106           << "  err=\"" << r.err << "\"\n";
0107
0108
0109
0110   return pT;
0111 }
0112
0113 float demoPt(){
0114
0115   // ======================== EMD Formula ==========================
0116   {
0117     SiCaloPt::InputEMD in;
0118     in.EMD_Angle  = 0.025;  // delta_phi - EM Deflection angle in rad
0119     in.EMD_Eta    = 0.00;   // EMCal Cluster eta
0120     in.EMD_Radius = 93.5;   // EMCal Cluster radius in cm
0121
0122     pTCalc.setParCeta(0.2);   // set C_eta parameter if needed
0123     pTCalc.setParPower(-1.0); // set Power parameter if needed
0124
0125     auto r = pTCalc.ComputePt(SiCaloPt::Method::MethodEMD, SiCaloPt::AnyInput{in});
0126     std::cout << "[EMD-analytic] ok=" << r.ok
0127               << "  pt=" << r.pt_reco
0128               << "  err=\"" << r.err << "\"\n";
0129
0130     in.EMD_Angle  = -0.025;  // delta_phi - EM Deflection angle in rad
0131     r = pTCalc.ComputePt(SiCaloPt::Method::MethodEMD, SiCaloPt::AnyInput{in});
0132     std::cout << "[EMD-analytic] ok=" << r.ok
0133               << "  pt=" << r.pt_reco
0134               << "  err=\"" << r.err << "\"\n";
0135   }
0136
0137   // ======================== Eproj Formula ==========================
0138   {
0139     SiCaloPt::InputEproj in;
0140     in.Energy_Calo   = 1.8;   // EMCal Cluster energy in GeV
0141     in.Radius_Calo   = 93.5;  // EMCal Cluster radius in cm
0142     in.Z_Calo        = 0.0;   // EMCal Cluster z in cm
0143     in.Radius_vertex = 0.0;   // Vertex radius in cm
0144     in.Z_vertex      = 0.0;   // Vertex z in cm
0145
0146     auto r = pTCalc.ComputePt(SiCaloPt::Method::MethodEproj, SiCaloPt::AnyInput{in});
0147     std::cout << "[Eproj-analytic] ok=" << r.ok
0148       << "  pt=" << r.pt_reco
0149       << "  err=\"" << r.err << "\"\n";
0150   }
0151
0152
0153   // ============= ML：MLEMD (2-d input: 1/dphi_EMD, eta_track = 0 now) ===============
0154   {
0155     // 2-d input features:{ dphi_EMD, eta_track }
0156     std::vector<float> featsMLEMD = {15, 0};
0157
0158     SiCaloPt::InputMLEMD in{featsMLEMD};
0159     auto r = pTCalc.ComputePt(SiCaloPt::Method::MethodMLEMD, SiCaloPt::AnyInput{in});
0160     std::cout << "[MLEMD-2D] ok=" << r.ok
0161       << "  pt=" << r.pt_reco
0162       << "  err=\"" << r.err << "\"\n";
0163   }
0164
0165   // ====== ML：MLEproj (7-d input: INTT 3/4 layer R,Z, INTT 5/4 layer R,Z, Calo R,Z,Energy) ======
0166   {
0167     // 7-d input features:{ INTT 3/4 layer R, INTT 3/4 layer Z, INTT 5/4 layer R, INTT 5/4 layer Z, Calo R, Calo Z, Calo Energy }
0168     std::vector<float> featsMLEproj = {
0169       10.0,  5.0,   // INTT 3/4 layer
0170       15.0,  7.5,   // INTT 5/4 layer
0171       100.0, 50.0, 8.0 }; // Calo R,Z,Energy
0172
0173     SiCaloPt::InputMLEproj in{featsMLEproj};
0174     auto r = pTCalc.ComputePt(SiCaloPt::Method::MethodMLEproj, SiCaloPt::AnyInput{in});
0175     std::cout << "[MLEproj-7D] ok=" << r.ok
0176               << "  pt=" << r.pt_reco
0177               << "  err=\"" << r.err << "\"\n";
0178   }
0179
0180   // ============ ML：Combined/Gate (2-d input: pt_from_MLEMD, pt_from_MLEproj) ===================
0181   {
0182     // 2-d input features:{ pt_from_MLEMD, pt_from_MLEproj }
0183     std::vector<float> featsMLCombined = {8.0, 9.5};
0184
0185     SiCaloPt::InputMLCombined in{featsMLCombined};
0186     auto r = pTCalc.ComputePt(SiCaloPt::Method::MethodMLCombined, SiCaloPt::AnyInput{in});
0187     std::cout << "[MLCombined] ok=" << r.ok
0188               << "  pt=" << r.pt_reco
0189               << "  err=\"" << r.err << "\"\n";
0190   }
0191
0192
0193
0194
0195
0196   return 0;
0197 };
0198
0199
0200 void analyze_SiSeedPair_pT(
0201      //const std::string& filename="/sphenix/user/jaein213/tracking/SiliconSeeding/MC/macro/ana/jpsi/ana_all.root",
0202 //     const std::string& filename="/sphenix/tg/tg01/commissioning/INTT/work/mahiro/SIliconCalo/MC/ana_jpsi/merged_200k.root",
0203 //     //const std::string& filename="/sphenix/user/jaein213/tracking/SiliconSeeding/MC/macro/ana/e+/inner_ana_0601_all.root",
0204 //     const std::string& filename="/sphenix/tg/tg01/commissioning/INTT/work/mahiro/SIliconCalo/MC/ana_e-/merged_10k.root",
0205 //     const std::string& filename="/sphenix/tg/tg01/commissioning/INTT/work/mahiro/SIliconCalo/MC/ana/ana_0_1kevt.root",
0206 //     const std::string& filename="/sphenix/tg/tg01/commissioning/INTT/work/mahiro/SIliconCalo/MC/ana/merged_100k.root", // pythia
0207 //     const std::string& filename="/sphenix/tg/tg01/commissioning/INTT/work/hachiya/SiCaloTrack/data/ana_em/all_em.root", // pythia
0208 //     const std::string& filename="/sphenix/tg/tg01/commissioning/INTT/work/hachiya/SiCaloTrack/data/ana_em/all_em1.root", // pythia
0209 //     const std::string& filename="/sphenix/tg/tg01/commissioning/INTT/work/hachiya/SiCaloTrack/data/ana_ep/all_ep.root", // pythia
0210      const std::string& filename="/sphenix/tg/tg01/commissioning/INTT/work/hachiya/SiCaloTrack/data/ana_jpsi/all_jpsi.root", // pythia
0211 //     const std::string& filename="/sphenix/tg/tg01/commissioning/INTT/work/hachiya/SiCaloTrack/data/ana_jpsi/all_jpsi1.root", // pythia
0212 //     const std::string& filename="/sphenix/tg/tg01/commissioning/INTT/work/hachiya/SiCaloTrack/data/ana_jpsi/ana_jpsi_0.root", // pythia
0213 //     const std::string& filename="/sphenix/tg/tg01/commissioning/INTT/work/hachiya/SiCaloTrack/data/ana_em/ana_em_0.root", // pythia
0214 //     const std::string& filename="/sphenix/user/hachiya/INTT/INTT/general_codes/hachiya/SiliconSeeding/macro/data/test_ana2k.root", // pythia
0215      //const std::string& filename="/sphenix/user/hachiya/INTT/INTT/general_codes/hachiya/SiliconSeeding/macro/data/test_ana5k.root", // pythia
0216      //const std::string& filename="/sphenix/user/hachiya/INTT/INTT/general_codes/hachiya/SiliconSeeding/macro/data/test_ana10k.root", // pythia
0217      //const std::string& filename="/sphenix/user/hachiya/INTT/INTT/general_codes/hachiya/SiliconSeeding/macro/data/test_ana50k.root", // pythia
0218      //const std::string& filename="/sphenix/user/hachiya/INTT/INTT/general_codes/hachiya/SiliconSeeding/macro/data/test_ana.root", // pythia
0219      //const std::string& filename="/sphenix/tg/tg01/commissioning/INTT/work/hachiya/SiCaloTrack/data/ana_pythia/ana_addall.root",
0220 //     const std::string& filename="/sphenix/tg/tg01/commissioning/INTT/work/hachiya/SiCaloTrack/data/ana/all_pythia.root",
0221 //     const std::string& filename="/gpfs/mnt/gpfs02/sphenix/user/hachiya/INTT/INTT/general_codes/hachiya/SiliconSeeding/macro/chargerecalc/ana_test_53879_10evt.root",
0222      //const std::string& filename="/sphenix/tg/tg01/commissioning/INTT/work/hachiya/SiCaloTrack/data/ana/ana10_pythia.root",
0223      float phi_threshold = 0.05
0224 )
0225 {
0226
0227   /////////////////////
0228   // for PtCalculator
0229   DemoPaths WS_Path;
0230   WS_Path.print();
0231
0232   //SiCaloPt::PtCalculatorConfig setting(if you want to use the ML models, Formula-method donnot need these)
0233   SiCaloPt::PtCalculatorConfig cfg;
0234   cfg.mlEMD_model_path        = make_opt<decltype(cfg.mlEMD_model_path)      >(WS_Path.emd_onnx);
0235   cfg.mlEMD_scaler_json       = make_opt<decltype(cfg.mlEMD_scaler_json)     >(WS_Path.emd_scaler_json);
0236   cfg.mlEproj_model_path      = make_opt<decltype(cfg.mlEproj_model_path)    >(WS_Path.eproj_onnx);
0237   cfg.mlEproj_scaler_json     = make_opt<decltype(cfg.mlEproj_scaler_json)   >(WS_Path.eproj_scaler_json);
0238   cfg.mlCombined_model_path   = make_opt<decltype(cfg.mlCombined_model_path) >(WS_Path.combined_onnx);
0239   cfg.mlCombined_scaler_json  = make_opt<decltype(cfg.mlCombined_scaler_json)>(WS_Path.combined_scaler_json);
0240
0241   // PtCalculator instance
0242   //SiCaloPt::PtCalculator calcTutorial(cfg);
0243   pTCalc.setConfig(cfg);
0244
0245   // initialize (load models and scalers for ML methods)
0246   std::string err;
0247   if (!pTCalc.init(&err))
0248   {
0249     std::cout << "[init] failed: " << err << std::endl;
0250     return;
0251   }
0252   std::cout << "[init] OK\n";
0253
0254
0255   demoPt();
0256   //calcPt();
0257   //return;
0258   /////////////////////
0259
0260
0261
0262
0263
0264
0265
0266
0267
0268   TFile *f = TFile::Open(filename.c_str(), "READ");
0269
0270   TTree *trackTree = (TTree *)f->Get("trackTree");
0271   TTree *caloTree  = (TTree *)f->Get("caloTree");
0272   TTree *siClusTree  = (TTree *)f->Get("SiClusTree");
0273   TTree *truthTree = (TTree *)f->Get("truthTree");
0274
0275   int evt;
0276   std::vector<float> *track_phi    = 0, *track_pt = 0, *track_pz = 0, *track_px = 0, *track_py = 0, *track_eta = 0, *track_z = 0;
0277   std::vector<float> *track_chi2ndf= 0;
0278   std::vector<int>   *track_charge = 0, *track_nmaps = 0, *track_nintt = 0;
0279   std::vector<float> *track_phi_emc= 0, *track_x_emc = 0, *track_y_emc = 0, *track_z_emc = 0;
0280   std::vector<float> *track_x_oemc = 0, *track_y_oemc = 0, *track_z_oemc = 0;
0281   std::vector<float> *track_rv_x_emc = 0, *track_rv_y_emc = 0, *track_rv_z_emc = 0;
0282
0283   int calo_evt;
0284   std::vector<float> *calo_phi = 0, *calo_energy = 0;
0285   std::vector<float> *calo_x = 0, *calo_y = 0, *calo_z = 0;
0286   std::vector<float> *calo_chi2 = 0, *calo_prob = 0;
0287
0288
0289   std::vector<float>* truth_pt=0, *truth_eta=0, *truth_phi=0;
0290   std::vector<float> *truth_px = 0, *truth_py = 0, *truth_pz = 0;
0291   std::vector<int>   *truth_pid=0, *truth_id=0;
0292   //std::vector<float>* truth_px, *truth_py, *truth_pz, *truth_e;
0293
0294   // SiClus
0295   std::vector<int> *SiClus_trackid=0;
0296   std::vector<int> *SiClus_layer=0;
0297   std::vector<float> *SiClus_x=0;
0298   std::vector<float> *SiClus_y=0;
0299   std::vector<float> *SiClus_z=0;
0300
0301
0302
0303   trackTree->SetBranchAddress("evt",     &evt);
0304   trackTree->SetBranchAddress("chi2ndf", &track_chi2ndf);
0305   trackTree->SetBranchAddress("charge",  &track_charge);
0306   trackTree->SetBranchAddress("nmaps",   &track_nmaps);
0307   trackTree->SetBranchAddress("nintt",   &track_nintt);
0308   trackTree->SetBranchAddress("pt0",     &track_pt);
0309   trackTree->SetBranchAddress("px0",     &track_px);
0310   trackTree->SetBranchAddress("py0",     &track_py);
0311   trackTree->SetBranchAddress("pz0",     &track_pz);
0312   trackTree->SetBranchAddress("eta0",    &track_eta);
0313   trackTree->SetBranchAddress("phi0",    &track_phi);
0314   trackTree->SetBranchAddress("z0",      &track_z);
0315   trackTree->SetBranchAddress("phi_proj_emc", &track_phi_emc);
0316   trackTree->SetBranchAddress("x_proj_emc",   &track_x_emc);
0317   trackTree->SetBranchAddress("y_proj_emc",   &track_y_emc);
0318   trackTree->SetBranchAddress("z_proj_emc",   &track_z_emc);
0319   trackTree->SetBranchAddress("x_proj_oemc",  &track_x_oemc);
0320   trackTree->SetBranchAddress("y_proj_oemc",  &track_y_oemc);
0321   trackTree->SetBranchAddress("z_proj_oemc",  &track_z_oemc);
0322   trackTree->SetBranchAddress("x_rv_proj_emc",  &track_rv_x_emc);
0323   trackTree->SetBranchAddress("y_rv_proj_emc",  &track_rv_y_emc);
0324   trackTree->SetBranchAddress("z_rv_proj_emc",  &track_rv_z_emc);
0325 /*
0326   trackTree->SetBranchAddress("pt0",      &track_pt);
0327   trackTree->SetBranchAddress("eta0",     &track_eta);
0328   trackTree->SetBranchAddress("phi0",     &track_phi);
0329   trackTree->SetBranchAddress("z0",       &track_z);
0330   trackTree->SetBranchAddress("phi_proj_emc", &track_phi_emc);
0331   trackTree->SetBranchAddress("z_proj_emc",   &track_z_emc);
0332 */
0333
0334 ///////////
0335    // Truth tree and branches
0336   if(truthTree!=nullptr){
0337     truthTree->SetBranchAddress("truth_pt",  &truth_pt);
0338     truthTree->SetBranchAddress("truth_pz",  &truth_pz);
0339     truthTree->SetBranchAddress("truth_eta", &truth_eta);
0340     truthTree->SetBranchAddress("truth_phi", &truth_phi);
0341     truthTree->SetBranchAddress("truth_px",  &truth_px);
0342     truthTree->SetBranchAddress("truth_py",  &truth_py);
0343     truthTree->SetBranchAddress("truth_pid", &truth_pid);
0344     truthTree->SetBranchAddress("truth_id",  &truth_id);
0345   }
0346
0347   //--truthTree->SetBranchAddress("truth_pid", &truth_pid);
0348   //--truthTree->SetBranchAddress("truth_px",  &truth_px);
0349   //--truthTree->SetBranchAddress("truth_py",  &truth_py);
0350   //--truthTree->SetBranchAddress("truth_e",   &truth_e);
0351   //--truthTree->SetBranchAddress("truth_pt",  &truth_pt);
0352
0353 ///////////
0354   siClusTree->SetBranchAddress("Siclus_trackid", &SiClus_trackid);
0355   siClusTree->SetBranchAddress("Siclus_layer",   &SiClus_layer);
0356   siClusTree->SetBranchAddress("Siclus_x",       &SiClus_x);
0357   siClusTree->SetBranchAddress("Siclus_y",       &SiClus_y);
0358   siClusTree->SetBranchAddress("Siclus_z",       &SiClus_z);
0359
0360
0361 ///////////
0362   caloTree->SetBranchAddress("calo_evt", &calo_evt);
0363   caloTree->SetBranchAddress("phi",      &calo_phi);
0364   caloTree->SetBranchAddress("x",        &calo_x);
0365   caloTree->SetBranchAddress("y",        &calo_y);
0366   caloTree->SetBranchAddress("z",        &calo_z);
0367   caloTree->SetBranchAddress("energy",   &calo_energy);
0368   caloTree->SetBranchAddress("chi2",     &calo_chi2);
0369   caloTree->SetBranchAddress("prob",     &calo_prob);
0370
0371 ///////////
0372
0373   //TFile *outFile = new TFile("dphi_distribution_em.root", "RECREATE");
0374   TFile *outFile = new TFile("dphi_distribution.root", "RECREATE");
0375   TH1F *h_dphi = new TH1F("h_dphi", "Track - Calo #Delta#phi;#Delta#phi;Counts", 200, -0.3, 0.3);
0376   TH1F *h_dphi_emc = new TH1F("h_dphi (EMCal Proj)", "Track - Calo #Delta#phi;#Delta#phi;Counts", 200, -0.3, 0.3);
0377   TH2F *h_dphi_emc_pt = new TH2F("h_dphi_pt", "Track - Calo #Delta#phi;pT;#Delta#phi", 200, 0, 20, 200, -0.3, 0.3);
0378   TH2F *h_dphi_emc_pt_truth = new TH2F("h_dphi_pt_truth", "Track - Calo #Delta#phi;pT{true};#Delta#phi", 200, 0, 20, 200, -0.3, 0.3);
0379   TH1F *h_EoverP_all = new TH1F("h_EoverP_all", "E/p of matched track-calo;E/p;Counts", 100, 0, 5);
0380   TH1F *h_EoverP_cut = new TH1F("h_EoverP_cut", "E/p of matched track-calo phi cut;E/p;Counts", 100, 0, 5);
0381   TH1F *h_dz = new TH1F("h_dz", "Track - Calo #Delta z;#Delta z (cm);Counts", 200, -50, 50);
0382   TH1F *h_dz_emc = new TH1F("h_dz_emc", "Track(EMCal Proj) - Calo #Delta z;#Delta z (cm);Counts", 200, -50, 50);
0383
0384   TH1F *h_mass     = new TH1F("h_mass", "Invariant mass of matched track pairs (e^{+}e^{-});Mass (GeV/c^{2});Counts", 200, 0, 5);
0385   // Invariant mass histograms for different J/psi pT bins
0386   TH1F *h_mass_pt0_1 = new TH1F("h_mass_pt0_1", "Invariant mass (pT 0-1 GeV/c);Mass (GeV/c^{2});Counts", 200, 0, 5);
0387   TH1F *h_mass_pt1_2 = new TH1F("h_mass_pt1_2", "Invariant mass (pT 1-2 GeV/c);Mass (GeV/c^{2});Counts", 200, 0, 5);
0388   TH1F *h_mass_pt2_3 = new TH1F("h_mass_pt2_3", "Invariant mass (pT 2-3 GeV/c);Mass (GeV/c^{2});Counts", 200, 0, 5);
0389   TH1F *h_mass_pt3_4 = new TH1F("h_mass_pt3_4", "Invariant mass (pT 3-4 GeV/c);Mass (GeV/c^{2});Counts", 200, 0, 5);
0390   TH1F *h_mass_pt4up = new TH1F("h_mass_pt4up", "Invariant mass (pT >4 GeV/c);Mass (GeV/c^{2});Counts", 200, 0, 5);
0391   TH1F *h_track_chi2ndf_matched = new TH1F("h_track_chi2ndf_matched", "Chi2/NDF of matched tracks;#chi^{2}/ndf;Counts", 100, 0, 10);
0392
0393   TH2F *h_reco_vs_truth_pt = new TH2F("h_reco_vs_truth_pt", "Reconstructed pT vs Truth pT;Truth pT (GeV/c);Reconstructed pT (GeV/c)", 100, 0, 20, 100, 0, 20);
0394
0395   TNtuple *h_ntp_sicalo = new TNtuple("ntp_sicalo", "SiSeed + Calo combination",       "pt:pz:c:chi2:nintt:nmaps:hitbit:eta:the:phi:z:phi_pemc:z_pemc:phi_poemc:z_poemc:phi_rvpemc:z_rvpemc:phi_intt:pt_tr:pz_tr:phi_tr:pid_tr:nemc:phi_emc:z_emc:e_emc:chi2_emc:phi_calo:pt_calo");
0396
0397   TNtuple *h_ntp_sicalotrk = new TNtuple("ntp_sicalotrk", "SiSeed + Calo combination", "pt:pz:c:chi2:nintt:nmaps:hitbit:eta:the:phi:z:phi_pemc:z_pemc:phi_poemc:z_poemc:phi_rvpemc:z_rvpemc:phi_intt:pt_tr:pz_tr:phi_tr:pid_tr:nemc:phi_emc:z_emc:e_emc:chi2_emc:phi_calo:pt_calo");
0398
0399   TNtuple *h_ntp_pair = new TNtuple("ntp_pair", "pair", "mass:pt:phi:eta:massc:ptc:ptp:pzp:eopp:dpp:dzp:ptm:pzm:eopm:dpm:dzm");
0400
0401
0402   // Define a struct to hold matched track info
0403   struct MatchedTrack {
0404     size_t track_idx;
0405     size_t calo_idx;
0406     float  min_dphi;
0407     float  min_dphi_emc;
0408     float  min_dz;
0409     float  min_dz_emc;
0410     float  eop;
0411     float  pt_calo;
0412     float  eop_calo;
0413     int    nmaps;
0414     int    nintt;
0415     int    charge;
0416     float  chi2ndf;
0417   };
0418
0419   struct stCluster {
0420     int   m_clsTrkid;
0421     int   m_clslyr  ;
0422     float m_clsx    ;
0423     float m_clsy    ;
0424     float m_clsz    ;
0425   };
0426   ////////////////////////////////////////////////////
0427
0428   Long64_t nentries = trackTree->GetEntries();
0429
0430   int nskip=0;
0431
0432   for (Long64_t i = 0; i < nentries; ++i) // Event Loop
0433   {
0434
0435     trackTree->GetEntry(i);
0436     caloTree->GetEntry(i);
0437     siClusTree->GetEntry(i);
0438     if(truthTree!=nullptr) truthTree->GetEntry(i);
0439
0440     if (evt != calo_evt)
0441     {
0442       std::cerr << "Warning: evt mismatch at entry " << i << ": track evt = " << evt << ", calo evt = " << calo_evt << std::endl;
0443       continue;
0444     }
0445
0446     if((evt%1000)==0) {
0447        std::cout << "Matching evt = " << evt << std::endl;
0448     }
0449
0450     std::vector<stCluster> vClus;
0451
0452
0453
0454     int nclusSize = SiClus_trackid->size();
0455
0456     cout<<evt <<" : ntrk, nemc : "<<track_phi->size()<<", "<< calo_phi->size()<< " "<<nclusSize<<endl;
0457     std::vector<MatchedTrack> matched_tracks;
0458
0459     int itotalClus=0;
0460     ///////////////////////////
0461     // First pass: find matched tracks and store info
0462     for (size_t it = 0; it < track_phi->size(); ++it)
0463     {
0464       int nmaps = (*track_nmaps)[it];
0465       int nintt = (*track_nintt)[it];
0466
0467       float& pt     = ((*track_pt)[it]);
0468       float& px     = ((*track_px)[it]);
0469       float& py     = ((*track_py)[it]);
0470       float& pz     = ((*track_pz)[it]);
0471       float& eta    = ((*track_eta)[it]);
0472       int&   charge = ((*track_charge)[it]);
0473
0474       int nclus = nmaps+nintt;
0475       std::cout << "charge : " << charge<<" "<<pt<<" "<<nmaps<<" "<<nintt<<" "<<calo_phi->size()<<" "<<nclus<<std::endl;
0476
0477       // cluster:
0478       int hitbit=0;
0479       stCluster iCl, oCl;
0480       for(int ic=0; ic<nclus; ic++) {
0481         int   clsid    = ic+itotalClus;
0482
0483         int layer =  ((*SiClus_layer)[clsid]);
0484         stCluster Cl {
0485           ((*SiClus_trackid)[clsid]),
0486           layer,
0487           ((*SiClus_x)[clsid]),
0488           ((*SiClus_y)[clsid]),
0489           ((*SiClus_z)[clsid])
0490         };
0491         vClus.push_back(Cl);
0492
0493         if(0<=Cl.m_clslyr&&Cl.m_clslyr<5) iCl = Cl;
0494         else                              oCl = Cl;
0495
0496         hitbit |= (1<<layer);
0497
0498         //cout<<"cls : "<<clsid<<" "<<Cl.m_clsTrkid<<" "<<Cl.m_clslyr<<" "<<Cl.m_clsx<<" "<<Cl.m_clsy<<" "<<Cl.m_clsz<<endl;
0499       }
0500       itotalClus+=nclus;
0501
0502       // truth track
0503       float pt_tr=0, pz_tr=0, phi_tr=0, eta_tr=0;
0504       int pid_tr=0;
0505
0506       if(truthTree!=nullptr){
0507         int ntruth = truth_pt->size();
0508         float min_dp = 99999;
0509         size_t min_itr=1000000;
0510         for (size_t itr = 0; itr < ntruth; ++itr)
0511         {
0512           int pid = ((*truth_pid)[itr]);
0513           //if( !(abs(pid)==11 // e-
0514           //  || abs(pid)==13 // mu-
0515           //  || abs(pid)==211 // pi+
0516           //  || abs(pid)==321 // K+
0517           //  || abs(pid)==2212 // p
0518           //  )) continue;
0519
0520
0521           float dpx = (px - ((*truth_px)[itr]));
0522           float dpy = (py - ((*truth_py)[itr]));
0523           float dpz = (pz - ((*truth_pz)[itr]));
0524           float dp  = sqrt(dpx*dpx+dpy*dpy+dpz*dpz);
0525           if(dp<min_dp) {
0526             min_dp = dp;
0527             min_itr = itr;
0528           }
0529
0530           //--cout<<"tru itr: "<<itr<<" "<<pid<<" "<<dp<<endl;
0531           //--pt_tr  = ((*truth_pt)[itr]);
0532           //--pz_tr  = ((*truth_pz)[itr]);
0533           //--phi_tr = ((*truth_phi)[itr]);
0534           //--eta_tr = ((*truth_eta)[itr]);
0535           //--px_tr  = ((*truth_px)[itr]);
0536           //--py_tr  = ((*truth_py)[itr]);
0537         }
0538
0539         if(min_itr<ntruth){
0540           pt_tr  = ((*truth_pt)[min_itr]);
0541           pz_tr  = ((*truth_pz)[min_itr]);
0542           phi_tr = ((*truth_phi)[min_itr]);
0543           eta_tr = ((*truth_eta)[min_itr]);
0544           pid_tr = ((*truth_pid)[min_itr]);
0545         }
0546         cout<<"tru : "<<ntruth<<" "<<min_itr<<" "<<pt_tr<<" "<<pz_tr<<" "<<phi_tr<<" "<<eta_tr<<" "<<pid_tr<<endl;
0547       }
0548
0549       if (pt < 0.3) continue;
0550
0551       if (nmaps < 1 || nintt < 1)
0552       {
0553         if(nskip<1000) {
0554           std::cout << "Skipping track with nmaps = " << nmaps << " and nintt = " << nintt << std::endl;
0555         } else if(nskip==1000) {
0556           std::cout << "exceed nskip. comment suppressed" << std::endl;
0557         }
0558         nskip++;
0559
0560         continue;
0561       }
0562
0563       //cout<<"cls-in : "<<iCl.m_clsTrkid<<" "<<iCl.m_clslyr<<" "<<iCl.m_clsx<<" "<<iCl.m_clsy<<" "<<iCl.m_clsz<<endl;
0564       //cout<<"cls-out: "<<oCl.m_clsTrkid<<" "<<oCl.m_clslyr<<" "<<oCl.m_clsx<<" "<<oCl.m_clsy<<" "<<oCl.m_clsz<<endl;
0565       //std::cout << "charge : " << charge<<" "<<pt<<std::endl;
0566
0567       float phi_intt = atan2(oCl.m_clsy-iCl.m_clsy, oCl.m_clsx-iCl.m_clsx);
0568
0569       float& x_emc = (*track_x_emc)[it];   //12
0570       float& y_emc = (*track_y_emc)[it];   //12
0571       float phi_proj = atan2(y_emc, x_emc);
0572
0573       float& x_oemc = (*track_x_oemc)[it];   //12
0574       float& y_oemc = (*track_y_oemc)[it];   //12
0575       float phi_proj_o = atan2(y_oemc, x_oemc);
0576
0577       float& x_rv_emc = (*track_rv_x_emc)[it];   //12
0578       float& y_rv_emc = (*track_rv_y_emc)[it];   //12
0579       float phi_proj_rv = atan2(y_rv_emc, x_rv_emc);
0580
0581       // Find calo cluster with minimum dphi_emc for this track
0582       float ntp_value[30] = {
0583                 pt,                   //0
0584                 pz,                   //1
0585                 (float)(charge),      //2
0586                 (*track_chi2ndf)[it], //3
0587                 (float)nintt,         //4
0588                 (float)nmaps,         //5
0589                 (float)hitbit,        //6
0590                 (*track_eta)[it],     //7
0591                 (float)(2.* std::atan( std::exp(-eta) )),//8
0592                 (*track_phi)[it],     //9
0593                 (*track_z)[it],       //10
0594                 phi_proj, //(*track_phi_emc)[it], //11
0595                 (*track_z_emc)[it],   //12
0596                 phi_proj_o, //(*track_phi_emc)[it], //13
0597                 (*track_z_oemc)[it],   //14
0598                 phi_proj_rv, //(*track_phi_emc)[it], //15
0599                 (*track_rv_z_emc)[it],   //16
0600                 phi_intt,             //17
0601                 pt_tr,                //18
0602                 pz_tr,                //19
0603                 phi_tr,               //20
0604                 (float)pid_tr,        //21
0605                 (float)calo_phi->size() //22
0606               };
0607
0608       PtInput input;
0609       input.SiIn[0] = iCl.m_clsx;
0610       input.SiIn[1] = iCl.m_clsy;
0611       input.SiIn[2] = iCl.m_clsz;
0612       input.SiOut[0] = oCl.m_clsx;
0613       input.SiOut[1] = oCl.m_clsy;
0614       input.SiOut[2] = oCl.m_clsz;
0615
0616       //const float par[2] = {0.21, -0.986}; // par for pT_calo calculation
0617       const float par[2] = {0.2, -1}; // par for pT_calo calculation
0618
0619       float min_dr = 1e9;
0620       size_t min_ic = calo_phi->size();
0621       float min_dphi = 0, min_dphi_emc = 0, min_dz = 0, min_dz_emc = 0;
0622       float match_pt_calo=0;
0623       //cout<<"aaaa"<<endl;
0624       for (size_t ic = 0; ic < calo_phi->size(); ++ic)
0625       {
0626       //--cout<<"bbb"<<endl;
0627
0628         //float dphi_emc = (*track_phi_emc)[it] - (*calo_phi)[ic];
0629         float dphi_emc = (*calo_phi)[ic] - phi_proj;
0630         if (dphi_emc >  M_PI) dphi_emc -= 2 * M_PI;
0631         if (dphi_emc < -M_PI) dphi_emc += 2 * M_PI;
0632
0633         const float p0 = -0.181669;
0634         const float p1 =  0.00389827;
0635         //float dphi_emc_corr = dphi_emc - charge*(p0/pt + p1);
0636         float dphi_emc_corr = charge*dphi_emc - 0.18*std::pow(pt, -0.986);
0637
0638         float x_calo = (*calo_x)[ic];
0639         float y_calo = (*calo_y)[ic];
0640         float z_calo = (*calo_z)[ic];
0641
0642         float dz_emc = z_calo - (*track_z_emc)[it];
0643
0644
0645         float phi_calo = atan2(y_calo - oCl.m_clsy,  x_calo - oCl.m_clsx);
0646
0647         float dphi = phi_calo - phi_intt;
0648         //float pt_calo = par[0]*pow(-charge*dphi, par[1]);//cal_CaloPt(dphi);
0649
0650         float pt_inv = charge*0.02+4.9*(-charge*dphi)-0.6*pow(-charge*dphi, 2);
0651         float pt_calo = 1./pt_inv;
0652
0653         input.Calo[0] = x_calo;
0654         input.Calo[1] = y_calo;
0655         input.Calo[2] = z_calo;
0656         input.CaloE   = (*calo_energy)[ic];
0657
0658
0659         //pt_calo = calcPt(input);
0660
0661         ntp_value[23] = (*calo_phi)[ic];
0662         ntp_value[24] = z_calo;
0663         ntp_value[25] = (*calo_energy)[ic];
0664         ntp_value[26] = (*calo_chi2)[ic];
0665         ntp_value[27] = phi_calo;
0666         ntp_value[28] = pt_calo;
0667
0668         h_ntp_sicalo->Fill(ntp_value);
0669         //--std::cout << "   energy : " << ntp_value[10]<<std::endl;
0670
0671         float dphi_emc_cm = dphi_emc*93.5/3.; // normalized by sigma=3
0672         float dr = sqrt(dphi_emc_cm*dphi_emc_cm + dz_emc*dz_emc);
0673
0674         //if (fabs(dphi_emc) < min_dr)
0675         //if (fabs(dphi_emc_corr) < min_dr)
0676         //if (fabs(dz_emc) < min_dr)
0677         if (fabs(dr) < min_dr)
0678         {
0679           //min_dr = fabs(dphi_emc);
0680           //min_dr = fabs(dphi_emc_corr);
0681
0682           //min_dr = fabs(dz_emc);
0683           min_dr = fabs(dr);
0684           min_ic = ic;
0685
0686           //min_dphi = (*track_phi)[it] - (*calo_phi)[ic];
0687           //if (min_dphi >  M_PI) min_dphi -= 2 * M_PI;
0688           //if (min_dphi < -M_PI) min_dphi += 2 * M_PI;
0689
0690           //min_dphi_emc = dphi_emc_corr;
0691           //min_dz = (*track_z)[it] - (*calo_z)[ic];
0692           //min_dz_emc = (*track_z_emc)[it] - (*calo_z)[ic];
0693           min_dphi_emc = dphi_emc;
0694           min_dz_emc   = dz_emc;
0695
0696           match_pt_calo = pt_calo;
0697         }
0698       }
0699
0700       if (min_ic != calo_phi->size())
0701       { // for matched track
0702         float x_calo = (*calo_x)[min_ic];
0703         float y_calo = (*calo_y)[min_ic];
0704         float z_calo = (*calo_z)[min_ic];
0705
0706         float phi_calo = atan2(y_calo - oCl.m_clsy,  x_calo - oCl.m_clsx);
0707
0708         float dphi = phi_calo - phi_intt;
0709         //float pt_calo = par[0]*pow(-charge*dphi, par[1]);//cal_CaloPt(dphi);
0710         float pt_inv = charge*0.02+4.9*(-charge*dphi)-0.6*pow(-charge*dphi, 2);
0711         float pt_calo = 1./pt_inv;
0712
0713         cout<<"new pT : "<<pt_calo<<endl;
0714
0715         input.Calo[0] = x_calo;
0716         input.Calo[1] = y_calo;
0717         input.Calo[2] = z_calo;
0718         input.CaloE   = (*calo_energy)[min_ic];
0719         calcPt(input);
0720         //pt_calo = calcPt(input);
0721
0722         ntp_value[23] = (*calo_phi)[min_ic];
0723         ntp_value[24] = z_calo;
0724         ntp_value[25] = (*calo_energy)[min_ic];
0725         ntp_value[26] = (*calo_chi2)[min_ic];
0726         ntp_value[27] = phi_calo;
0727         ntp_value[28] = pt_calo;
0728       }
0729       else {
0730         ntp_value[23] = -9999.;
0731         ntp_value[24] = -9999.;
0732         ntp_value[25] = -9999.;
0733         ntp_value[26] = -9999.;
0734         ntp_value[27] = -9999.;
0735         ntp_value[28] = -9999.;
0736       }
0737       h_ntp_sicalotrk->Fill(ntp_value);
0738
0739       if (min_ic == calo_phi->size()) {
0740         continue; // no match
0741       }
0742
0743       float p = (*track_pt)[it] * std::cosh((*track_eta)[it]);
0744       float E = (*calo_energy)[min_ic];
0745       float eop = -999;
0746       if (p > 0)
0747       {
0748         eop = E / p;
0749         h_EoverP_all->Fill(eop);
0750       }
0751       float p_calo = match_pt_calo * std::cosh((*track_eta)[it]);
0752       float eop_calo = E/p_calo;
0753
0754       h_dphi->Fill(min_dphi);
0755       h_dphi_emc->Fill(min_dphi_emc);
0756       h_dphi_emc_pt->Fill((*track_pt)[it], min_dphi_emc);
0757       //h_dphi_emc_pt_truth->Fill((*truth_pt)[0], min_dphi_emc);
0758       h_dz->Fill(min_dz);
0759       h_dz_emc->Fill(min_dz_emc);
0760      // std::cout << "  ➤ Δφ = " << min_dphi << ", Δz = " << min_dz << std::endl;
0761
0762       h_track_chi2ndf_matched->Fill((*track_chi2ndf)[it]);
0763       //--if ((*track_chi2ndf)[it] > 5)
0764       //--  continue; // Skip tracks with high chi2/ndf
0765       if (min_dz_emc > 4)
0766         continue;
0767      // h_reco_vs_truth_pt->Fill((*truth_pt)[0], (*track_pt)[it]);
0768
0769       matched_tracks.push_back({it, min_ic, min_dphi, min_dphi_emc, min_dz, min_dz_emc, eop, match_pt_calo, eop_calo, nmaps, nintt, (*track_charge)[it], (*track_chi2ndf)[it]});
0770     }
0771
0772     //cout<<"nclus-end : "<<itotalClus<<endl;
0773
0774     /////////////////////////
0775     //--// 1.5 pass: split positive and negative to different array;
0776     //--std::vector<MatchedTrack> v_trkp, v_trkm;
0777     //--for (size_t idx = 0; idx < matched_tracks.size(); ++idx)
0778     //--{
0779     //--  const auto& mt = matched_tracks[idx];
0780     //--  if(mt.charge>0) v_trkp.push_back(mt);
0781     //--  else            v_trkm.push_back(mt);
0782     //--}
0783
0784     /////////////////////////
0785     // Second pass: loop through matched tracks and fill histograms and compute invariant mass
0786     cout<<"matched track : "<<matched_tracks.size()<<endl;
0787     float v_pair[20];
0788     for (size_t idx = 0; idx < matched_tracks.size(); ++idx)
0789     {
0790       const auto& mt = matched_tracks[idx];
0791
0792       if (mt.nmaps > 2 && mt.nintt > 1 && std::abs(mt.min_dz_emc) < 4)
0793       {
0794         h_EoverP_cut->Fill(mt.eop);
0795
0796         //--if (!(mt.eop > 0.01 && mt.eop < 2)) continue;
0797
0798         for (size_t jdx = idx+1; jdx < matched_tracks.size(); ++jdx)
0799         {
0800           const auto& mt2 = matched_tracks[jdx];
0801
0802           //if (mt2.nmaps < 2 || mt2.nintt < 1) continue;
0803           //if (std::abs(mt2.min_dz_emc) > 4) continue;
0804
0805           if (mt2.nmaps > 2 && mt2.nintt > 1 && std::abs(mt2.min_dz_emc) < 4) {
0806
0807             if (mt.charge * mt2.charge >= 0) continue; // opposite sign only
0808
0809             //float eop2 = -999;
0810             //float p2 = (*track_pt)[mt2.track_idx] * std::cosh((*track_eta)[mt2.track_idx]);
0811             //float E2_raw = (*calo_energy)[mt2.calo_idx];
0812             //if (p2 > 0) eop2 = E2_raw / p2;
0813             //if (!(eop2 > 0.8 && eop2 < 2)) continue;
0814
0815             //--if (!(mt2.eop > 0.6 && mt2.eop < 2)) continue;
0816
0817             float pt1 = (*track_pt)[mt.track_idx],  eta1 = (*track_eta)[mt.track_idx],  phi1 = (*track_phi)[mt.track_idx];
0818             float pt2 = (*track_pt)[mt2.track_idx], eta2 = (*track_eta)[mt2.track_idx], phi2 = (*track_phi)[mt2.track_idx];
0819             TLorentzVector v1, v2;
0820             v1.SetPtEtaPhiM(pt1, eta1, phi1, electron_mass);
0821             v2.SetPtEtaPhiM(pt2, eta2, phi2, electron_mass);
0822             //v1.SetPtEtaPhiM(pt1, eta1, phi1, pion_mass);
0823             //v2.SetPtEtaPhiM(pt2, eta2, phi2, pion_mass);
0824
0825             TLorentzVector vpair = v1 + v2;
0826
0827             //float mass = (v1 + v2).M();
0828             //float total_pt = pt1 + pt2;
0829             float mass     = vpair.M();
0830             float total_pt = vpair.Pt();
0831
0832             // mass with pt_calo
0833             float ptc1 = mt.pt_calo;
0834             float ptc2 = mt2.pt_calo;
0835             TLorentzVector vc1, vc2;
0836             vc1.SetPtEtaPhiM(ptc1, eta1, phi1, electron_mass);
0837             vc2.SetPtEtaPhiM(ptc2, eta2, phi2, electron_mass);
0838             //vc1.SetPtEtaPhiM(ptc1, eta1, phi1, pion_mass);
0839             //vc2.SetPtEtaPhiM(ptc2, eta2, phi2, pion_mass);
0840             TLorentzVector vpairc = vc1 + vc2;
0841
0842             float massc     = vpairc.M();
0843             float total_ptc = vpairc.Pt();
0844
0845             //if ((0.8 < mt.eop  && mt.eop  < 2) &&
0846             //    (0.8 < mt2.eop && mt2.eop < 2))
0847             {
0848               h_mass->Fill(mass);
0849
0850               // Fill the appropriate mass histogram based on total pT
0851               if (total_pt < 1.0) {
0852                 h_mass_pt0_1->Fill(mass);
0853               } else if (total_pt < 2.0) {
0854                 h_mass_pt1_2->Fill(mass);
0855               } else if (total_pt < 3.0) {
0856                 h_mass_pt2_3->Fill(mass);
0857               } else if (total_pt < 4.0) {
0858                 h_mass_pt3_4->Fill(mass);
0859               } else {
0860                 h_mass_pt4up->Fill(mass);
0861               }
0862             }
0863
0864     //TNtuple *h_ntp_pair = new TNtuple("ntp_pair", "SiSeed + Calo combination",
0865     //"mass:pt:phi:eta:ptp:pzp:phip:dpp:dzp:ptm:pzm:phim:dpm:dzm");
0866
0867             v_pair[ 0] = mass;
0868             v_pair[ 1] = total_pt;
0869             v_pair[ 2] = vpair.Phi();
0870             v_pair[ 3] = vpair.Eta();
0871             v_pair[ 4] = massc; // mass with pt_calo
0872             v_pair[ 5] = total_ptc; // pair-pt with pt_calo
0873             v_pair[ 6] = pt1;
0874             v_pair[ 7] = v1.Pz();
0875             v_pair[ 8] = mt.eop; //mt.eop_calo;
0876             v_pair[ 9] = mt.min_dphi_emc;
0877             v_pair[10] = mt.min_dz_emc;
0878             v_pair[11] = pt2;
0879             v_pair[12] = v2.Pz();
0880             v_pair[13] = mt2.eop; //mt2.eop_calo;
0881             v_pair[14] = mt2.min_dphi_emc;
0882             v_pair[15] = mt2.min_dz_emc;
0883
0884             h_ntp_pair->Fill(v_pair);
0885           }
0886         }
0887       }
0888     }
0889   }
0890
0891   std::cout<<"total nskip : "<< nskip <<std::endl;
0892
0893  // TFile *outFile = new TFile("dphi_distribution_e-.root", "RECREATE");
0894   //TFile *outFile = new TFile("dphi_distribution_PYTHIA_temp.root", "RECREATE");
0895   std::cout << "Saving histogram to dphi_distribution.root" << std::endl;
0896   h_dphi->Write();
0897   h_dphi_emc->Write();
0898   h_dphi_emc_pt->Write();
0899   h_dphi_emc_pt_truth->Write();
0900   h_EoverP_all->Write();
0901   h_EoverP_cut->Write();
0902   h_dz->Write();
0903   h_dz_emc->Write();
0904   h_mass->Write();
0905   h_mass_pt0_1->Write();
0906   h_mass_pt1_2->Write();
0907   h_mass_pt2_3->Write();
0908   h_mass_pt3_4->Write();
0909   h_mass_pt4up->Write();
0910   h_track_chi2ndf_matched->Write();
0911   h_reco_vs_truth_pt->Write();
0912   h_ntp_sicalo->Write();
0913   h_ntp_sicalotrk->Write();
0914   h_ntp_pair->Write();
0915
0916   outFile->Close();
0917
0918   std::cout << "Δφ histogram saved to 'dphi_distribution.root'" << std::endl;
0919   f->Close();
0920 }