SiCalo/anamacro/analyze_SiSeedPair.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 const float electron_mass = 0.000511; // GeV/c^2
0014 const float pion_mass     = 0.1396; // GeV/c^2
0015
0016 void analyze_SiSeedPair(
0017      //const std::string& filename="/sphenix/user/jaein213/tracking/SiliconSeeding/MC/macro/ana/jpsi/ana_all.root",
0018 //     const std::string& filename="/sphenix/tg/tg01/commissioning/INTT/work/mahiro/SIliconCalo/MC/ana_jpsi/merged_200k.root",
0019 //     //const std::string& filename="/sphenix/user/jaein213/tracking/SiliconSeeding/MC/macro/ana/e+/inner_ana_0601_all.root",
0020 //     const std::string& filename="/sphenix/tg/tg01/commissioning/INTT/work/mahiro/SIliconCalo/MC/ana_e-/merged_10k.root",
0021 //     const std::string& filename="/sphenix/tg/tg01/commissioning/INTT/work/mahiro/SIliconCalo/MC/ana/ana_0_1kevt.root",
0022 //     const std::string& filename="/sphenix/tg/tg01/commissioning/INTT/work/mahiro/SIliconCalo/MC/ana/merged_100k.root", // pythia
0023      const std::string& filename="/sphenix/tg/tg01/commissioning/INTT/work/hachiya/SiCaloTrack/data/ana_em/all_em.root", // pythia
0024 //     const std::string& filename="/sphenix/tg/tg01/commissioning/INTT/work/hachiya/SiCaloTrack/data/ana_ep/all_ep.root", // pythia
0025 //     const std::string& filename="/sphenix/tg/tg01/commissioning/INTT/work/hachiya/SiCaloTrack/data/ana_jpsi/all_jpsi.root", // pythia
0026 //     const std::string& filename="/sphenix/user/hachiya/INTT/INTT/general_codes/hachiya/SiliconSeeding/macro/data/test_ana2k.root", // pythia
0027      //const std::string& filename="/sphenix/user/hachiya/INTT/INTT/general_codes/hachiya/SiliconSeeding/macro/data/test_ana5k.root", // pythia
0028      //const std::string& filename="/sphenix/user/hachiya/INTT/INTT/general_codes/hachiya/SiliconSeeding/macro/data/test_ana10k.root", // pythia
0029      //const std::string& filename="/sphenix/user/hachiya/INTT/INTT/general_codes/hachiya/SiliconSeeding/macro/data/test_ana50k.root", // pythia
0030      //const std::string& filename="/sphenix/user/hachiya/INTT/INTT/general_codes/hachiya/SiliconSeeding/macro/data/test_ana.root", // pythia
0031      //const std::string& filename="/sphenix/tg/tg01/commissioning/INTT/work/hachiya/SiCaloTrack/data/ana_pythia/ana_addall.root",
0032      //const std::string& filename="/sphenix/tg/tg01/commissioning/INTT/work/hachiya/SiCaloTrack/data/ana/all_pythia.root",
0033      float phi_threshold = 0.05
0034                        )
0035 {
0036   TFile *f = TFile::Open(filename.c_str(), "READ");
0037
0038   TTree *trackTree = (TTree *)f->Get("trackTree");
0039   TTree *caloTree  = (TTree *)f->Get("caloTree");
0040   TTree *siClusTree  = (TTree *)f->Get("SiClusTree");
0041   TTree *truthTree = (TTree *)f->Get("truthTree");
0042
0043   int evt;
0044   std::vector<float> *track_phi    = 0, *track_pt = 0, *track_pz = 0, *track_eta = 0, *track_z = 0;
0045   std::vector<float> *track_chi2ndf= 0;
0046   std::vector<int>   *track_charge = 0, *track_nmaps = 0, *track_nintt = 0;
0047   std::vector<float> *track_phi_emc= 0, *track_z_emc = 0;
0048
0049   int calo_evt;
0050   std::vector<float> *calo_phi = 0, *calo_energy = 0;
0051   std::vector<float> *calo_x = 0, *calo_y = 0, *calo_z = 0;
0052   std::vector<float> *calo_chi2 = 0, *calo_prob = 0;
0053
0054
0055   std::vector<float>* truth_pt=0, *truth_eta=0, *truth_phi=0;
0056   std::vector<float> *truth_pz = 0;
0057   //std::vector<float>* truth_px, *truth_py, *truth_pz, *truth_e;
0058
0059   // SiClus
0060   std::vector<int> *SiClus_trackid=0;
0061   std::vector<int> *SiClus_layer=0;
0062   std::vector<float> *SiClus_x=0;
0063   std::vector<float> *SiClus_y=0;
0064   std::vector<float> *SiClus_z=0;
0065
0066
0067
0068   trackTree->SetBranchAddress("evt",     &evt);
0069   trackTree->SetBranchAddress("chi2ndf", &track_chi2ndf);
0070   trackTree->SetBranchAddress("charge",  &track_charge);
0071   trackTree->SetBranchAddress("nmaps",   &track_nmaps);
0072   trackTree->SetBranchAddress("nintt",   &track_nintt);
0073   trackTree->SetBranchAddress("pt0",     &track_pt);
0074   trackTree->SetBranchAddress("pz0",     &track_pz);
0075   trackTree->SetBranchAddress("eta0",    &track_eta);
0076   trackTree->SetBranchAddress("phi0",    &track_phi);
0077   trackTree->SetBranchAddress("z0",      &track_z);
0078   trackTree->SetBranchAddress("phi_proj_emc", &track_phi_emc);
0079   trackTree->SetBranchAddress("z_proj_emc",   &track_z_emc);
0080 /*
0081   trackTree->SetBranchAddress("pt0",      &track_pt);
0082   trackTree->SetBranchAddress("eta0",     &track_eta);
0083   trackTree->SetBranchAddress("phi0",     &track_phi);
0084   trackTree->SetBranchAddress("z0",       &track_z);
0085   trackTree->SetBranchAddress("phi_proj_emc", &track_phi_emc);
0086   trackTree->SetBranchAddress("z_proj_emc",   &track_z_emc);
0087 */
0088
0089 ///////////
0090    // Truth tree and branches
0091   if(truthTree!=nullptr){
0092     truthTree->SetBranchAddress("truth_pt",  &truth_pt);
0093     truthTree->SetBranchAddress("truth_pz",  &truth_pz);
0094     truthTree->SetBranchAddress("truth_eta", &truth_eta);
0095     truthTree->SetBranchAddress("truth_phi", &truth_phi);
0096   }
0097
0098   //--truthTree->SetBranchAddress("truth_pid", &truth_pid);
0099   //--truthTree->SetBranchAddress("truth_px",  &truth_px);
0100   //--truthTree->SetBranchAddress("truth_py",  &truth_py);
0101   //--truthTree->SetBranchAddress("truth_e",   &truth_e);
0102   //--truthTree->SetBranchAddress("truth_pt",  &truth_pt);
0103
0104 ///////////
0105   siClusTree->SetBranchAddress("Siclus_trackid", &SiClus_trackid);
0106   siClusTree->SetBranchAddress("Siclus_layer",   &SiClus_layer);
0107   siClusTree->SetBranchAddress("Siclus_x",       &SiClus_x);
0108   siClusTree->SetBranchAddress("Siclus_y",       &SiClus_y);
0109   siClusTree->SetBranchAddress("Siclus_z",       &SiClus_z);
0110
0111
0112 ///////////
0113   caloTree->SetBranchAddress("calo_evt", &calo_evt);
0114   caloTree->SetBranchAddress("phi",      &calo_phi);
0115   caloTree->SetBranchAddress("x",        &calo_x);
0116   caloTree->SetBranchAddress("y",        &calo_y);
0117   caloTree->SetBranchAddress("z",        &calo_z);
0118   caloTree->SetBranchAddress("energy",   &calo_energy);
0119   caloTree->SetBranchAddress("chi2",     &calo_chi2);
0120   caloTree->SetBranchAddress("prob",     &calo_prob);
0121
0122 ///////////
0123
0124   //TFile *outFile = new TFile("dphi_distribution_em.root", "RECREATE");
0125   TFile *outFile = new TFile("dphi_distribution.root", "RECREATE");
0126   TH1F *h_dphi = new TH1F("h_dphi", "Track - Calo #Delta#phi;#Delta#phi;Counts", 200, -0.3, 0.3);
0127   TH1F *h_dphi_emc = new TH1F("h_dphi (EMCal Proj)", "Track - Calo #Delta#phi;#Delta#phi;Counts", 200, -0.3, 0.3);
0128   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);
0129   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);
0130   TH1F *h_EoverP_all = new TH1F("h_EoverP_all", "E/p of matched track-calo;E/p;Counts", 100, 0, 5);
0131   TH1F *h_EoverP_cut = new TH1F("h_EoverP_cut", "E/p of matched track-calo phi cut;E/p;Counts", 100, 0, 5);
0132   TH1F *h_dz = new TH1F("h_dz", "Track - Calo #Delta z;#Delta z (cm);Counts", 200, -50, 50);
0133   TH1F *h_dz_emc = new TH1F("h_dz_emc", "Track(EMCal Proj) - Calo #Delta z;#Delta z (cm);Counts", 200, -50, 50);
0134
0135   TH1F *h_mass     = new TH1F("h_mass", "Invariant mass of matched track pairs (e^{+}e^{-});Mass (GeV/c^{2});Counts", 200, 0, 5);
0136   // Invariant mass histograms for different J/psi pT bins
0137   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);
0138   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);
0139   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);
0140   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);
0141   TH1F *h_mass_pt4up = new TH1F("h_mass_pt4up", "Invariant mass (pT >4 GeV/c);Mass (GeV/c^{2});Counts", 200, 0, 5);
0142   TH1F *h_track_chi2ndf_matched = new TH1F("h_track_chi2ndf_matched", "Chi2/NDF of matched tracks;#chi^{2}/ndf;Counts", 100, 0, 10);
0143
0144   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);
0145
0146   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_intt:pt_tr:pz_tr:phi_tr:phi_emc:z_emc:e_emc:chi2_emc:phi_calo:pt_calo");
0147
0148   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_intt:pt_tr:pz_tr:phi_tr:phi_emc:z_emc:e_emc:chi2_emc:phi_calo:pt_calo");
0149
0150   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");
0151
0152
0153   // Define a struct to hold matched track info
0154   struct MatchedTrack {
0155     size_t track_idx;
0156     size_t calo_idx;
0157     float  min_dphi;
0158     float  min_dphi_emc;
0159     float  min_dz;
0160     float  min_dz_emc;
0161     float  eop;
0162     float  pt_calo;
0163     float  eop_calo;
0164     int    nmaps;
0165     int    nintt;
0166     int    charge;
0167     float  chi2ndf;
0168   };
0169
0170   struct stCluster {
0171     int   m_clsTrkid;
0172     int   m_clslyr  ;
0173     float m_clsx    ;
0174     float m_clsy    ;
0175     float m_clsz    ;
0176   };
0177   ////////////////////////////////////////////////////
0178
0179   Long64_t nentries = trackTree->GetEntries();
0180
0181   int nskip=0;
0182
0183   for (Long64_t i = 0; i < nentries; ++i) // Event Loop
0184   {
0185
0186     trackTree->GetEntry(i);
0187     caloTree->GetEntry(i);
0188     siClusTree->GetEntry(i);
0189     if(truthTree!=nullptr) truthTree->GetEntry(i);
0190
0191     if (evt != calo_evt)
0192     {
0193       std::cerr << "Warning: evt mismatch at entry " << i << ": track evt = " << evt << ", calo evt = " << calo_evt << std::endl;
0194       continue;
0195     }
0196
0197     if((evt%1000)==0) {
0198        std::cout << "Matching evt = " << evt << std::endl;
0199     }
0200
0201     std::vector<stCluster> vClus;
0202
0203
0204
0205     int nclusSize = SiClus_trackid->size();
0206
0207     cout<<evt <<" : ntrk, nemc : "<<track_phi->size()<<", "<< calo_phi->size()<< " "<<nclusSize<<endl;
0208     std::vector<MatchedTrack> matched_tracks;
0209
0210     int itotalClus=0;
0211     ///////////////////////////
0212     // First pass: find matched tracks and store info
0213     for (size_t it = 0; it < track_phi->size(); ++it)
0214     {
0215       int nmaps = (*track_nmaps)[it];
0216       int nintt = (*track_nintt)[it];
0217
0218       float& pt     = ((*track_pt)[it]);
0219       float& pz     = ((*track_pz)[it]);
0220       float& eta    = ((*track_eta)[it]);
0221       int&   charge = ((*track_charge)[it]);
0222
0223       int nclus = nmaps+nintt;
0224       //std::cout << "charge : " << charge<<" "<<pt<<" "<<nmaps<<" "<<nintt<<" "<<calo_phi->size()<<" "<<nclus<<std::endl;
0225
0226       // cluster:
0227       int hitbit=0;
0228       stCluster iCl, oCl;
0229       for(int ic=0; ic<nclus; ic++) {
0230         int   clsid    = ic+itotalClus;
0231
0232         int layer =  ((*SiClus_layer)[clsid]);
0233         stCluster Cl {
0234           ((*SiClus_trackid)[clsid]),
0235           layer,
0236           ((*SiClus_x)[clsid]),
0237           ((*SiClus_y)[clsid]),
0238           ((*SiClus_z)[clsid])
0239         };
0240         vClus.push_back(Cl);
0241
0242         if(0<=Cl.m_clslyr&&Cl.m_clslyr<5) iCl = Cl;
0243         else                              oCl = Cl;
0244
0245         hitbit |= (1<<layer);
0246
0247         //cout<<"cls : "<<clsid<<" "<<Cl.m_clsTrkid<<" "<<Cl.m_clslyr<<" "<<Cl.m_clsx<<" "<<Cl.m_clsy<<" "<<Cl.m_clsz<<endl;
0248       }
0249       itotalClus+=nclus;
0250
0251       // truth track
0252       float pt_tr=0, pz_tr=0, phi_tr=0, eta_tr=0;
0253
0254       if(truthTree!=nullptr){
0255         int ntruth = truth_pt->size();
0256         for (size_t itr = 0; itr < (ntruth>0?1:0); ++itr)
0257         {
0258           pt_tr  = ((*truth_pt)[itr]);
0259           pz_tr  = ((*truth_pz)[itr]);
0260           phi_tr = ((*truth_phi)[itr]);
0261           eta_tr = ((*truth_eta)[itr]);
0262           cout<<"tru : "<<itr<<" "<<pt_tr<<" "<<pz_tr<<" "<<phi_tr<<" "<<eta_tr<<endl;
0263         }
0264       }
0265
0266       if (pt < 0.3) continue;
0267
0268       if (nmaps < 1 || nintt < 1)
0269       {
0270         if(nskip<1000) {
0271           std::cout << "Skipping track with nmaps = " << nmaps << " and nintt = " << nintt << std::endl;
0272         } else if(nskip==1000) {
0273           std::cout << "exceed nskip. comment suppressed" << std::endl;
0274         }
0275         nskip++;
0276
0277         continue;
0278       }
0279
0280       //cout<<"cls-in : "<<iCl.m_clsTrkid<<" "<<iCl.m_clslyr<<" "<<iCl.m_clsx<<" "<<iCl.m_clsy<<" "<<iCl.m_clsz<<endl;
0281       //cout<<"cls-out: "<<oCl.m_clsTrkid<<" "<<oCl.m_clslyr<<" "<<oCl.m_clsx<<" "<<oCl.m_clsy<<" "<<oCl.m_clsz<<endl;
0282       //std::cout << "charge : " << charge<<" "<<pt<<std::endl;
0283
0284       float phi_intt = atan2(oCl.m_clsy-iCl.m_clsy, oCl.m_clsx-iCl.m_clsx);
0285
0286
0287       // Find calo cluster with minimum dphi_emc for this track
0288       float ntp_value[30] = {
0289                 pt,                   //0
0290                 pz,                   //1
0291                 (float)(charge),      //2
0292                 (*track_chi2ndf)[it], //3
0293                 (float)nintt,         //4
0294                 (float)nmaps,         //5
0295                 (float)hitbit,        //6
0296                 (*track_eta)[it],     //7
0297                 (float)(2.* std::atan( std::exp(-eta) )),//8
0298                 (*track_phi)[it],     //9
0299                 (*track_z)[it],       //10
0300                 (*track_phi_emc)[it], //11
0301                 (*track_z_emc)[it],   //12
0302                 phi_intt,             //13
0303                 pt_tr,                //14
0304                 pz_tr,                //15
0305                 phi_tr                //16
0306               };
0307
0308
0309       const float par[2] = {0.21, -0.986}; // par for pT_calo calculation
0310
0311       float min_dr = 1e9;
0312       size_t min_ic = calo_phi->size();
0313       float min_dphi = 0, min_dphi_emc = 0, min_dz = 0, min_dz_emc = 0;
0314       float match_pt_calo=0;
0315       //cout<<"aaaa"<<endl;
0316       for (size_t ic = 0; ic < calo_phi->size(); ++ic)
0317       {
0318       //--cout<<"bbb"<<endl;
0319
0320         float dphi_emc = (*track_phi_emc)[it] - (*calo_phi)[ic];
0321         if (dphi_emc >  M_PI) dphi_emc -= 2 * M_PI;
0322         if (dphi_emc < -M_PI) dphi_emc += 2 * M_PI;
0323
0324         const float p0 = -0.181669;
0325         const float p1 =  0.00389827;
0326         //float dphi_emc_corr = dphi_emc - charge*(p0/pt + p1);
0327         float dphi_emc_corr = charge*dphi_emc - 0.18*std::pow(pt, -0.986);
0328
0329         float x_calo = (*calo_x)[ic];
0330         float y_calo = (*calo_y)[ic];
0331         float z_calo = (*calo_z)[ic];
0332
0333         float dz_emc = (*track_z_emc)[it] - z_calo;
0334
0335
0336         float phi_calo = atan2(y_calo - oCl.m_clsy,  x_calo - oCl.m_clsx);
0337
0338         float dphi = phi_calo - phi_intt;
0339         float pt_calo = par[0]*pow(-charge*dphi, par[1]);//cal_CaloPt(dphi);
0340
0341         ntp_value[17]  = (*calo_phi)[ic];
0342         ntp_value[18] = z_calo;
0343         ntp_value[19] = (*calo_energy)[ic];
0344         ntp_value[20] = (*calo_chi2)[ic];
0345         ntp_value[21] = phi_calo;
0346         ntp_value[22] = pt_calo;
0347
0348         h_ntp_sicalo->Fill(ntp_value);
0349         //--std::cout << "   energy : " << ntp_value[10]<<std::endl;
0350
0351
0352         if (fabs(dphi_emc) < min_dr)
0353         //if (fabs(dphi_emc_corr) < min_dr)
0354         //if (fabs(dz_emc) < min_dr)
0355         {
0356           //min_dr = fabs(dphi_emc);
0357           //min_dr = fabs(dphi_emc_corr);
0358
0359           min_dr = fabs(dz_emc);
0360           min_ic = ic;
0361
0362           //min_dphi = (*track_phi)[it] - (*calo_phi)[ic];
0363           //if (min_dphi >  M_PI) min_dphi -= 2 * M_PI;
0364           //if (min_dphi < -M_PI) min_dphi += 2 * M_PI;
0365
0366           //min_dphi_emc = dphi_emc_corr;
0367           //min_dz = (*track_z)[it] - (*calo_z)[ic];
0368           //min_dz_emc = (*track_z_emc)[it] - (*calo_z)[ic];
0369
0370           //match_pt_calo = pt_calo;
0371         }
0372       }
0373
0374       if (min_ic != calo_phi->size())
0375       { // for matched track
0376         float x_calo = (*calo_x)[min_ic];
0377         float y_calo = (*calo_y)[min_ic];
0378         float z_calo = (*calo_z)[min_ic];
0379
0380         float phi_calo = atan2(y_calo - oCl.m_clsy,  x_calo - oCl.m_clsx);
0381
0382         float dphi = phi_calo - phi_intt;
0383         float pt_calo = par[0]*pow(-charge*dphi, par[1]);//cal_CaloPt(dphi);
0384
0385         ntp_value[17]  = (*calo_phi)[min_ic];
0386         ntp_value[18] = z_calo;
0387         ntp_value[19] = (*calo_energy)[min_ic];
0388         ntp_value[20] = (*calo_chi2)[min_ic];
0389         ntp_value[21] = phi_calo;
0390         ntp_value[22] = pt_calo;
0391       }
0392       else {
0393         ntp_value[17] = -9999.;
0394         ntp_value[18] = -9999.;
0395         ntp_value[19] = -9999.;
0396         ntp_value[20] = -9999.;
0397         ntp_value[21] = -9999.;
0398         ntp_value[22] = -9999.;
0399       }
0400       h_ntp_sicalotrk->Fill(ntp_value);
0401
0402       if (min_ic == calo_phi->size()) {
0403         continue; // no match
0404       }
0405
0406       float p = (*track_pt)[it] * std::cosh((*track_eta)[it]);
0407       float E = (*calo_energy)[min_ic];
0408       float eop = -999;
0409       if (p > 0)
0410       {
0411         eop = E / p;
0412         h_EoverP_all->Fill(eop);
0413       }
0414       float p_calo = match_pt_calo * std::cosh((*track_eta)[it]);
0415       float eop_calo = E/p_calo;
0416
0417       h_dphi->Fill(min_dphi);
0418       h_dphi_emc->Fill(min_dphi_emc);
0419       h_dphi_emc_pt->Fill((*track_pt)[it], min_dphi_emc);
0420       //h_dphi_emc_pt_truth->Fill((*truth_pt)[0], min_dphi_emc);
0421       h_dz->Fill(min_dz);
0422       h_dz_emc->Fill(min_dz_emc);
0423      // std::cout << "  ➤ Δφ = " << min_dphi << ", Δz = " << min_dz << std::endl;
0424
0425       h_track_chi2ndf_matched->Fill((*track_chi2ndf)[it]);
0426       if ((*track_chi2ndf)[it] > 5)
0427         continue; // Skip tracks with high chi2/ndf
0428       if (min_dz_emc > 4)
0429         continue;
0430      // h_reco_vs_truth_pt->Fill((*truth_pt)[0], (*track_pt)[it]);
0431
0432       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]});
0433     }
0434
0435     //cout<<"nclus-end : "<<itotalClus<<endl;
0436
0437     /////////////////////////
0438     //--// 1.5 pass: split positive and negative to different array;
0439     //--std::vector<MatchedTrack> v_trkp, v_trkm;
0440     //--for (size_t idx = 0; idx < matched_tracks.size(); ++idx)
0441     //--{
0442     //--  const auto& mt = matched_tracks[idx];
0443     //--  if(mt.charge>0) v_trkp.push_back(mt);
0444     //--  else            v_trkm.push_back(mt);
0445     //--}
0446
0447     /////////////////////////
0448     // Second pass: loop through matched tracks and fill histograms and compute invariant mass
0449     cout<<"matched track : "<<matched_tracks.size()<<endl;
0450     float v_pair[20];
0451     for (size_t idx = 0; idx < matched_tracks.size(); ++idx)
0452     {
0453       const auto& mt = matched_tracks[idx];
0454
0455       if (mt.nmaps > 2 && mt.nintt > 1 && std::abs(mt.min_dz_emc) < 4)
0456       {
0457         h_EoverP_cut->Fill(mt.eop);
0458
0459         if (!(mt.eop > 0.6 && mt.eop < 2)) continue;
0460
0461         for (size_t jdx = idx+1; jdx < matched_tracks.size(); ++jdx)
0462         {
0463           const auto& mt2 = matched_tracks[jdx];
0464
0465           //if (mt2.nmaps < 2 || mt2.nintt < 1) continue;
0466           //if (std::abs(mt2.min_dz_emc) > 4) continue;
0467
0468           if (mt2.nmaps > 2 && mt2.nintt > 1 && std::abs(mt2.min_dz_emc) < 4) {
0469
0470             if (mt.charge * mt2.charge >= 0) continue; // opposite sign only
0471
0472             //float eop2 = -999;
0473             //float p2 = (*track_pt)[mt2.track_idx] * std::cosh((*track_eta)[mt2.track_idx]);
0474             //float E2_raw = (*calo_energy)[mt2.calo_idx];
0475             //if (p2 > 0) eop2 = E2_raw / p2;
0476             //if (!(eop2 > 0.8 && eop2 < 2)) continue;
0477
0478             if (!(mt2.eop > 0.6 && mt2.eop < 2)) continue;
0479
0480             float pt1 = (*track_pt)[mt.track_idx],  eta1 = (*track_eta)[mt.track_idx],  phi1 = (*track_phi)[mt.track_idx];
0481             float pt2 = (*track_pt)[mt2.track_idx], eta2 = (*track_eta)[mt2.track_idx], phi2 = (*track_phi)[mt2.track_idx];
0482             TLorentzVector v1, v2;
0483             v1.SetPtEtaPhiM(pt1, eta1, phi1, electron_mass);
0484             v2.SetPtEtaPhiM(pt2, eta2, phi2, electron_mass);
0485             //v1.SetPtEtaPhiM(pt1, eta1, phi1, pion_mass);
0486             //v2.SetPtEtaPhiM(pt2, eta2, phi2, pion_mass);
0487
0488             TLorentzVector vpair = v1 + v2;
0489
0490             //float mass = (v1 + v2).M();
0491             //float total_pt = pt1 + pt2;
0492             float mass     = vpair.M();
0493             float total_pt = vpair.Pt();
0494
0495             // mass with pt_calo
0496             float ptc1 = mt.pt_calo;
0497             float ptc2 = mt2.pt_calo;
0498             TLorentzVector vc1, vc2;
0499             vc1.SetPtEtaPhiM(ptc1, eta1, phi1, electron_mass);
0500             vc2.SetPtEtaPhiM(ptc2, eta2, phi2, electron_mass);
0501             //vc1.SetPtEtaPhiM(ptc1, eta1, phi1, pion_mass);
0502             //vc2.SetPtEtaPhiM(ptc2, eta2, phi2, pion_mass);
0503             TLorentzVector vpairc = vc1 + vc2;
0504
0505             float massc     = vpairc.M();
0506             float total_ptc = vpairc.Pt();
0507
0508             //if ((0.8 < mt.eop  && mt.eop  < 2) &&
0509             //    (0.8 < mt2.eop && mt2.eop < 2))
0510             {
0511               h_mass->Fill(mass);
0512
0513               // Fill the appropriate mass histogram based on total pT
0514               if (total_pt < 1.0) {
0515                 h_mass_pt0_1->Fill(mass);
0516               } else if (total_pt < 2.0) {
0517                 h_mass_pt1_2->Fill(mass);
0518               } else if (total_pt < 3.0) {
0519                 h_mass_pt2_3->Fill(mass);
0520               } else if (total_pt < 4.0) {
0521                 h_mass_pt3_4->Fill(mass);
0522               } else {
0523                 h_mass_pt4up->Fill(mass);
0524               }
0525             }
0526
0527     //TNtuple *h_ntp_pair = new TNtuple("ntp_pair", "SiSeed + Calo combination",
0528     //"mass:pt:phi:eta:ptp:pzp:phip:dpp:dzp:ptm:pzm:phim:dpm:dzm");
0529
0530             v_pair[ 0] = mass;
0531             v_pair[ 1] = total_pt;
0532             v_pair[ 2] = vpair.Phi();
0533             v_pair[ 3] = vpair.Eta();
0534             v_pair[ 4] = massc; // mass with pt_calo
0535             v_pair[ 5] = total_ptc; // pair-pt with pt_calo
0536             v_pair[ 6] = pt1;
0537             v_pair[ 7] = v1.Pz();
0538             v_pair[ 8] = mt.eop_calo;
0539             v_pair[ 9] = mt.min_dphi_emc;
0540             v_pair[10] = mt.min_dz_emc;
0541             v_pair[11] = pt2;
0542             v_pair[12] = v2.Pz();
0543             v_pair[13] = mt2.eop_calo;
0544             v_pair[15] = mt2.min_dphi_emc;
0545             v_pair[15] = mt2.min_dz_emc;
0546
0547             h_ntp_pair->Fill(v_pair);
0548           }
0549         }
0550       }
0551     }
0552   }
0553
0554   std::cout<<"total nskip : "<< nskip <<std::endl;
0555
0556  // TFile *outFile = new TFile("dphi_distribution_e-.root", "RECREATE");
0557   //TFile *outFile = new TFile("dphi_distribution_PYTHIA_temp.root", "RECREATE");
0558   std::cout << "Saving histogram to dphi_distribution.root" << std::endl;
0559   h_dphi->Write();
0560   h_dphi_emc->Write();
0561   h_dphi_emc_pt->Write();
0562   h_dphi_emc_pt_truth->Write();
0563   h_EoverP_all->Write();
0564   h_EoverP_cut->Write();
0565   h_dz->Write();
0566   h_dz_emc->Write();
0567   h_mass->Write();
0568   h_mass_pt0_1->Write();
0569   h_mass_pt1_2->Write();
0570   h_mass_pt2_3->Write();
0571   h_mass_pt3_4->Write();
0572   h_mass_pt4up->Write();
0573   h_track_chi2ndf_matched->Write();
0574   h_reco_vs_truth_pt->Write();
0575   h_ntp_sicalo->Write();
0576   h_ntp_sicalotrk->Write();
0577   h_ntp_pair->Write();
0578
0579   outFile->Close();
0580
0581   std::cout << "Δφ histogram saved to 'dphi_distribution.root'" << std::endl;
0582   f->Close();
0583 }