InttSeedTrackPerformance-Jingyu/macro/calc_pt_fromFunc.C

0001 #include <TFile.h>
0002 #include <TTree.h>
0003 #include <TH3D.h>
0004 #include <TCanvas.h>
0005 #include <TVector2.h>
0006 #include <vector>
0007 #include <cmath>
0008
0009 void calc_pt_fromFunc(
0010     // const char* filename="/mnt/e/sphenix/INTT-EMCAL/InttSeedingTrackDev/ML4Reco/dataset/merged_last500kNew.root",
0011     const char* filename="/mnt/e/sphenix/INTT-EMCAL/InttSeedingTrackDev/ML4Reco/dataset/positron_last500kNew.root",
0012     // const char* filename="/mnt/e/sphenix/INTT-EMCAL/InttSeedingTrackDev/ParticleGen/output/electron_1M.root",
0013     const char* treename="tree",
0014
0015     // const char* ptfuncfilename="/mnt/e/sphenix/INTT-EMCAL/InttSeedingTrackDev/InttSeedTrackPerformance/output/Ptfunc_fitCofEta_positron.root",
0016     const char* ptfuncfilename="/mnt/e/sphenix/INTT-EMCAL/InttSeedingTrackDev/InttSeedTrackPerformance/output/Ptfunc_fitCofEta_takashi_p.root",
0017     // const char* ptfuncfilename="/mnt/e/sphenix/INTT-EMCAL/InttSeedingTrackDev/InttSeedTrackPerformance/output/Ptfunc_fitCofEta.root",
0018     // const char* ptfuncfilename="/mnt/e/sphenix/INTT-EMCAL/InttSeedingTrackDev/InttSeedTrackPerformance/output/Ptfunc_fitCofEta_takashi_e.root",
0019
0020     const char* outputname="/mnt/e/sphenix/INTT-EMCAL/InttSeedingTrackDev/InttSeedTrackPerformance/output/Ptfunc_Performance_positron_wTa.root"
0021     // const char* outputname="/mnt/e/sphenix/INTT-EMCAL/InttSeedingTrackDev/InttSeedTrackPerformance/output/Ptfunc_Performance_electron_wTa.root"
0022 )
0023 {
0024     TFile outfile(outputname, "RECREATE");
0025     TH2D *h2_ptreso = new TH2D("h2_ptreso", "h2_ptreso; p_{T}^{reco} [GeV]; (p_{T}^{reco} - p_{T}^{truth}) / p_{T}^{truth}", 110, 0, 11, 200, -2, 2);
0026
0027     TFile *file_in = TFile::Open(filename);
0028     if (!file_in || file_in->IsZombie()) {
0029         std::cerr<<"Error: cannot open "<<filename<<std::endl;
0030         return;
0031     }
0032
0033     TTree *tree_data = (TTree*)file_in->Get(treename);
0034     if (!tree_data) {
0035         std::cerr<<"Error: cannot find tree "<<treename<<std::endl;
0036         return;
0037     }
0038
0039     TFile* ptfuncfile = TFile::Open(ptfuncfilename, "READ");
0040     if (!ptfuncfile || ptfuncfile->IsZombie()) {
0041         std::cerr << "Error: cannot open function file " << ptfuncfilename << std::endl;
0042         return;
0043     }
0044     TF1* ptfuncTP = nullptr;
0045     TF1* ptfuncTG = nullptr;
0046     ptfuncfile->GetObject("poly_func_profile", ptfuncTP);
0047     ptfuncfile->GetObject("poly_func_graph", ptfuncTG);
0048
0049     TFile* file_C = new TFile("/mnt/e/sphenix/INTT-EMCAL/InttSeedingTrackDev/ParticleGen/output/calo_positron_dphi_ptbin_woC.root", "READ");
0050     // TFile* file_C = new TFile("/mnt/e/sphenix/INTT-EMCAL/InttSeedingTrackDev/ParticleGen/output/calo_electron_dphi_ptbin_woC.root", "READ");
0051     TF1* f1_dphi_C = (TF1*)file_C->Get("fitCurve");
0052     TGraph* g1_dphi_C = (TGraph*)file_C->Get("grPeakVsX");
0053
0054     // 定义变量
0055     int NClus;
0056     std::vector <int> *clus_system = nullptr;
0057     std::vector <int> *clus_layer = nullptr;
0058     std::vector <int> *clus_adc = nullptr;
0059     std::vector <double> *clus_X = nullptr;
0060     std::vector <double> *clus_Y = nullptr;
0061     std::vector <double> *clus_Z = nullptr;
0062
0063     std::vector<double> *CEMC_Pr_Hit_x = nullptr;
0064     std::vector<double> *CEMC_Pr_Hit_y = nullptr;
0065     std::vector<double> *CEMC_Pr_Hit_z = nullptr;
0066     std::vector<double> *CEMC_Pr_Hit_R = nullptr;
0067
0068     std::vector<double> *PrimaryG4P_Pt_ = nullptr;
0069     std::vector<double> *PrimaryG4P_Eta_ = nullptr;
0070
0071     std::vector<double> *caloClus_innr_x  = nullptr;
0072     std::vector<double> *caloClus_innr_y = nullptr;
0073     std::vector<double> *caloClus_innr_z = nullptr;
0074     std::vector<double> *caloClus_innr_phi  = nullptr;
0075     std::vector<double> *caloClus_innr_edep = nullptr;
0076
0077     tree_data->SetBranchAddress("trk_NClus", & NClus);
0078     tree_data->SetBranchAddress("trk_system", & clus_system);
0079     tree_data->SetBranchAddress("trk_layer", & clus_layer);
0080     tree_data->SetBranchAddress("trk_X", & clus_X);
0081     tree_data->SetBranchAddress("trk_Y", & clus_Y);
0082     tree_data->SetBranchAddress("trk_Z", & clus_Z);
0083
0084     tree_data->SetBranchAddress("CEMC_Pr_Hit_x",  &CEMC_Pr_Hit_x);
0085     tree_data->SetBranchAddress("CEMC_Pr_Hit_y",  &CEMC_Pr_Hit_y);
0086     tree_data->SetBranchAddress("CEMC_Pr_Hit_z",  &CEMC_Pr_Hit_z);
0087     tree_data->SetBranchAddress("CEMC_Pr_Hit_R",  &CEMC_Pr_Hit_R);
0088
0089     tree_data->SetBranchAddress("PrimaryG4P_Pt", &PrimaryG4P_Pt_);
0090     tree_data->SetBranchAddress("PrimaryG4P_Eta", &PrimaryG4P_Eta_);
0091
0092     tree_data->SetBranchAddress("caloClus_innr_X", &caloClus_innr_x);
0093     tree_data->SetBranchAddress("caloClus_innr_Y", &caloClus_innr_y);
0094     tree_data->SetBranchAddress("caloClus_innr_Z", &caloClus_innr_z);
0095     tree_data->SetBranchAddress("caloClus_innr_Phi", &caloClus_innr_phi);
0096     tree_data->SetBranchAddress("caloClus_innr_edep", &caloClus_innr_edep);
0097
0098     Long64_t nentries = tree_data->GetEntries();
0099     for (Long64_t i=0; i<nentries; ++i)
0100     {
0101         tree_data->GetEntry(i);
0102
0103         // 选事件：层4/5 恰有1个、层6/7 恰有1个
0104         int idx34 = -1, idx56 = -1;
0105         int cnt34 = 0, cnt56 = 0;
0106         for (int j=0; j<(int)clus_layer->size(); ++j) {
0107             int L = (*clus_layer)[j];
0108             if (L==3 || L==4) { ++cnt34; idx34=j; }
0109             if (L==5 || L==6) { ++cnt56; idx56=j; }
0110         }
0111         if (cnt34!=1 || cnt56!=1) continue;
0112
0113         // 要求 CEMC 记录为1
0114         if (CEMC_Pr_Hit_x->size()!=1) continue;
0115
0116         // 要求 prim particle 为1
0117         if (PrimaryG4P_Pt_->size()!=1) continue;
0118
0119         // 取出所有量
0120         double trk34_X = (*clus_X)[idx34];
0121         double trk34_Y = (*clus_Y)[idx34];
0122         double trk56_X = (*clus_X)[idx56];
0123         double trk56_Y = (*clus_Y)[idx56];
0124         double cemc_Pr_X  = (*CEMC_Pr_Hit_x)[0];
0125         double cemc_Pr_Y  = (*CEMC_Pr_Hit_y)[0];
0126         double Pr_pt = (*PrimaryG4P_Pt_)[0];
0127         double Pr_eta = (*PrimaryG4P_Eta_)[0];
0128
0129         double cemc_Pr_R = std::sqrt(cemc_Pr_X * cemc_Pr_X + cemc_Pr_Y * cemc_Pr_Y);
0130         // if (cemc_Pr_R > 102) continue; // 限制 CEMC 半径范围
0131
0132         // 计算 Δφ
0133         double phi1 = std::atan2(trk56_Y - trk34_Y, trk56_X - trk34_X);
0134         double phi2 = std::atan2(cemc_Pr_Y - trk56_Y, cemc_Pr_X - trk56_X);
0135         double dphi = TVector2::Phi_mpi_pi(phi2 - phi1);
0136
0137         for (size_t j = 0; j < caloClus_innr_x->size(); j++)
0138         {
0139             double cluster_innr_x = caloClus_innr_x->at(j);
0140             double cluster_innr_y = caloClus_innr_y->at(j);
0141             double cluster_innr_z = caloClus_innr_z->at(j);
0142             double cluster_innr_phi = caloClus_innr_phi->at(j);
0143             double cluster_innr_edep = caloClus_innr_edep->at(j);
0144             double cluster_innr_R = std::sqrt(cluster_innr_x*cluster_innr_x + cluster_innr_y*cluster_innr_y);
0145             double cluster_innr_eta = std::asinh(cluster_innr_z / cluster_innr_R);
0146
0147             double Ecalo_threshold = 0.5;
0148             if((cluster_innr_edep < Ecalo_threshold)&&(cluster_innr_R>140)) continue;
0149
0150             // clus innr position correction
0151             // double dphi_C = f1_dphi_C->Eval(cluster_innr_edep);
0152             double dphi_C = g1_dphi_C->Eval(cluster_innr_edep);
0153
0154             if (cluster_innr_phi < 0)
0155             {
0156                 cluster_innr_phi += 2 * TMath::Pi();
0157             }
0158             cluster_innr_phi = cluster_innr_phi + dphi_C;
0159             cluster_innr_x = cluster_innr_R*cos(cluster_innr_phi); // update position xy
0160             cluster_innr_y = cluster_innr_R*sin(cluster_innr_phi); // update position xy
0161
0162             double phi2_reco = std::atan2(cluster_innr_y - trk56_Y, cluster_innr_x - trk56_X);
0163             double dphi_reco = TVector2::Phi_mpi_pi(phi2_reco - phi1);
0164
0165             dphi_reco = std::abs(dphi_reco); // 确保 Δφ 为正值
0166
0167             // double Cval_eval = ptfuncTP->Eval(cluster_innr_eta);
0168             double Cval_eval = ptfuncTG->Eval(cluster_innr_eta);
0169             // double pt_reco = Cval_eval / dphi_reco;
0170             double pt_reco = Cval_eval / (std::pow(dphi_reco, 0.986));
0171
0172             h2_ptreso->Fill(Pr_pt, (pt_reco-Pr_pt)/Pr_pt);
0173         }
0174     }
0175
0176     gStyle->SetOptFit(1111);
0177     TH1D* projY = h2_ptreso->ProjectionY("projY", 10, 11, "e");
0178     // Fit with Gaussian
0179     TF1* gausFit = new TF1("gausFit", "gaus", -1, 1);
0180     projY->Fit(gausFit, "R", "", -0.02, 0.1);
0181
0182     TCanvas* c1 = new TCanvas("c1","Pt reso about 1 GeV",800,600);
0183     projY->Draw();
0184     gausFit->SetRange(-0.02, 0.1);
0185     gausFit->Draw("same");
0186     c1->SaveAs("Pt_reso_fromFunc_1GeV.pdf");
0187
0188     outfile.cd();
0189     h2_ptreso->Write();
0190     projY->Write();
0191     gausFit->Write();
0192 }
0193
0194