InttSeedTrackPerformance-Jingyu/macro/CaloInfo.C

0001 #include <iostream>
0002 #include <string>
0003 #include <vector>
0004 #include <cmath>
0005 #include <TChain.h>
0006 #include <TTree.h>
0007 #include <TH1F.h>
0008 #include <TFile.h>
0009 #include <TMath.h>
0010 #include <TProfile.h>
0011
0012 double Cal_Distance_R(const double xyz_1[3], const double xyz_2[3]);
0013 double Cal_dR(const double xyz_1[3], const double xyz_2[3]);
0014 double Cal_dphiR(const double xyz_1[3], const double xyz_2[3]);
0015 TGraph* Cor_var_Z(TH2D* h2, double xrange_min = -0.05, double xrange_max = 0.05);
0016 Double_t crystalball(Double_t* x, Double_t* par);
0017
0018 TGraph* FitPeakVsX(TH2D* h2_XY, TF1* projectionFit = nullptr, TF1* peakVsXFit = nullptr, const char* graphName = "grPeakVsX");
0019
0020
0021 void CaloInfo()
0022 {
0023     // 工作目录
0024     std::string fDir = "/mnt/e/sphenix/INTT-EMCAL/InttSeedingTrackDev/ParticleGen/output";
0025
0026     TFile* fin2M = new TFile("/mnt/e/sphenix/INTT-EMCAL/InttSeedingTrackDev/ParticleGen/output/Var_Z_electron.root", "READ");
0027     TGraph* g_dphiM_z = (TGraph*)fin2M->Get("g_cor_z_dphi_distri_test");
0028     TGraph* g_dRadM_z = (TGraph*)fin2M->Get("g_cor_z_dRad_distri_test");
0029
0030     TFile* file_C = new TFile("/mnt/e/sphenix/INTT-EMCAL/InttSeedingTrackDev/ParticleGen/output/calo_positron_dphi_ptbin_woC.root", "READ");
0031     // TFile* file_C = new TFile("/mnt/e/sphenix/INTT-EMCAL/InttSeedingTrackDev/ParticleGen/output/calo_electron_dphi_ptbin_woC.root", "READ");
0032     TF1* f1_dphi_C = (TF1*)file_C->Get("fitCurve");
0033     TGraph* g1_dphi_C = (TGraph*)file_C->Get("grPeakVsX");
0034
0035     // 输入文件
0036     TChain *tc = new TChain("tree");
0037     // std::string fInputName = fDir + "/ana457_Electron_5GeV_10k.root";
0038     // std::string fInputName = fDir + "/ana457_Positron_5GeV_10k.root";
0039     // std::string fInputName = "/mnt/e/sphenix/INTT-EMCAL/InttSeedingTrackDev/ML4Reco/dataset/positron_last500kNew.root";
0040     // std::string fInputName = "/mnt/e/sphenix/INTT-EMCAL/InttSeedingTrackDev/ML4Reco/dataset/merged_last500kNew.root";
0041     std::string fInputName = "/mnt/e/sphenix/INTT-EMCAL/InttSeedingTrackDev/ParticleGen/output/electron_1M.root";
0042
0043     tc->Add(fInputName.c_str());
0044
0045     // 输出文件
0046     // std::string fOutputName = fDir + "/Var_Z_electron.root";
0047     // std::string fOutputName = fDir + "/calo_positron_dphi_ptbin_wC.root";
0048     std::string fOutputName = fDir + "/calo_electron_dphi_ptbin_woC.root";
0049
0050     TTree *caloinfotree = (TTree*)tc;
0051     if (!caloinfotree)
0052     {
0053         std::cerr << "Error: Input tree is null!" << std::endl;
0054         return;
0055     }
0056
0057     // 定义变量
0058     std::vector<double> *CEMC_Pr_Hit_x = nullptr;
0059     std::vector<double> *CEMC_Pr_Hit_y = nullptr;
0060     std::vector<double> *CEMC_Pr_Hit_z = nullptr;
0061     std::vector<double> *CEMC_Pr_Hit_R = nullptr;
0062
0063     std::vector<double> *CEMC_Hit_x = nullptr;
0064     std::vector<double> *CEMC_Hit_y = nullptr;
0065     std::vector<double> *CEMC_Hit_z = nullptr;
0066     std::vector<double> *CEMC_Hit_R = nullptr;
0067     std::vector<double> *CEMC_Hit_Edep = nullptr;
0068
0069     std::vector<double> *tower_innr_x = nullptr;
0070     std::vector<double> *tower_innr_y = nullptr;
0071     std::vector<double> *tower_innr_z = nullptr;
0072     std::vector<double> *tower_innr_R = nullptr;
0073     std::vector<double> *tower_innr_phi = nullptr;
0074     std::vector<double> *tower_innr_eta = nullptr;
0075     std::vector<double> *tower_innr_edep = nullptr;
0076
0077     std::vector<double> *tower_system = nullptr;
0078     std::vector<double> *tower_x = nullptr;
0079     std::vector<double> *tower_y = nullptr;
0080     std::vector<double> *tower_z = nullptr;
0081     std::vector<double> *tower_R = nullptr;
0082     std::vector<double> *tower_phi = nullptr;
0083     std::vector<double> *tower_eta = nullptr;
0084     std::vector<double> *tower_edep = nullptr;
0085
0086     std::vector<double> *caloClus_innr_x  = nullptr;
0087     std::vector<double> *caloClus_innr_y  = nullptr;
0088     std::vector<double> *caloClus_innr_z  = nullptr;
0089     std::vector<double> *caloClus_innr_R  = nullptr;
0090     std::vector<double> *caloClus_innr_phi  = nullptr;
0091     std::vector<double> *caloClus_innr_edep = nullptr;
0092
0093     std::vector<double> *caloClus_system = nullptr;
0094     std::vector<double> *caloClus_x  = nullptr;
0095     std::vector<double> *caloClus_y  = nullptr;
0096     std::vector<double> *caloClus_z  = nullptr;
0097     std::vector<double> *caloClus_R  = nullptr;
0098     std::vector<double> *caloClus_phi  = nullptr;
0099     std::vector<double> *caloClus_edep = nullptr;
0100
0101     std::vector<double> *PrimaryG4P_Pt = nullptr;
0102     std::vector<double> *PrimaryG4P_E  = nullptr;
0103
0104     // 定义直方图
0105     auto h_g4hitphi = new TH1F("h_g4hitphi", "Azimuthal angle;Phi (rad);Counts", 100, -TMath::Pi(), TMath::Pi());
0106     auto h_towerphi = new TH1F("h_towerphi", "Azimuthal angle;Phi (rad);Counts", 100, -TMath::Pi(), TMath::Pi());
0107     auto h_clustphi = new TH1F("h_clustphi", "Azimuthal angle;Phi (rad);Counts", 100, -TMath::Pi(), TMath::Pi());
0108     auto h_innr_towerphi = new TH1F("h_innr_towerphi", "Azimuthal angle;Phi (rad);Counts", 100, -TMath::Pi(), TMath::Pi());
0109     auto h_innr_clustphi = new TH1F("h_innr_clustphi", "Azimuthal angle;Phi (rad);Counts", 100, -TMath::Pi(), TMath::Pi());
0110
0111     auto h1_DR_pr1_clustinnr   = new TH1F("h1_DR_pr1_clustinnr", "h_DR_pr1_clustinnr;R(cm);Counts", 100, 0, 5);
0112     auto h1_DR_pr1_towerinnr   = new TH1F("h1_DR_pr1_towerinnr", "h_DR_pr1_towerinnr;R(cm);Counts", 100, 0, 20);
0113     auto h1_DR_g4hit_clustgeom = new TH1F("h1_DR_g4hit_clustgeom", "h_DR_g4hit_clustgeom;R(cm);Counts", 200, -20, 20);
0114     auto h1_DR_g4hit_towergeom = new TH1F("h1_DR_g4hit_towergeom", "h_DR_g4hit_towergeom;R(cm);Counts", 100, 0, 20);
0115
0116     auto h1_dRadius_pr1_clustinnr  = new TH1F("h1_dRadius_pr1_clustinnr", "h1_dRadius_pr1_clustinnr;R(cm);Counts", 50, -2, 2);
0117     auto h1_dRadius_g4hit_clustgeom= new TH1F("h1_dRadius_g4hit_clustgeom", "h1_dRadius_g4hit_clustgeom;R(cm);Counts", 50, -20, 20);
0118     auto h1_dphiR_pr1_clustinnr    = new TH1F("h1_dphiR_pr1_clustinnr", "h1_dphiR_pr1_clustinnr;R(cm);Counts", 50, -2, 2);
0119     auto h1_dphiR_g4hit_clustgeom  = new TH1F("h1_dphiR_g4hit_clustgeom", "h1_dphiR_g4hit_clustgeom;R(cm);Counts", 50, -20, 20);
0120
0121     auto h2_DR_pr1_clustinnr_pt   = new TH2F("h2_DR_pr1_clustinnr_pt", "h2_DR_pr1_clustinnr_pt;pt(GeV);R(cm);Counts", 11, -0.5, 10.5, 100, 0, 20);
0122     auto h2_DR_g4hit_clustgeom_pt = new TH2F("h2_DR_g4hit_clustgeom_pt", "h2_DR_g4hit_clustgeom_pt;pt(GeV);R(cm);Counts", 11, -0.5, 10.5, 100, 0, 20);
0123     auto h2_DR_pr1_clustinnr_E    = new TH2F("h2_DR_pr1_clustinnr_E", "h2_DR_pr1_clustinnr_E;E(GeV);R(cm);Counts", 15, -0.5, 14.5, 100, 0, 20);
0124     auto h2_DR_g4hit_clustgeom_E  = new TH2F("h2_DR_g4hit_clustgeom_E", "h2_DR_g4hit_clustgeom_E;E(GeV);R(cm);Counts", 15, -0.5, 14.5, 100, 0, 20);
0125
0126     auto h2_zr_g4hit = new TH2D("h2_zr_g4hit", "g4hit energy-weighted RZ distribution;Z (cm);R (cm);Weighted Counts", 100, -150, 150, 200, 0, 200);
0127     auto h2_zr_tower = new TH2D("h2_zr_tower", "tower energy-weighted RZ distribution;Z (cm);R (cm);Weighted Counts", 100, -150, 150, 200, 0, 200);
0128     auto h2_zr_clust = new TH2D("h2_zr_clust", "clust energy-weighted RZ distribution;Z (cm);R (cm);Weighted Counts", 100, -150, 150, 200, 0, 200);
0129     auto h2_zr_tower_innr = new TH2D("h2_zr_tower_innr", "tower energy-weighted RZ distribution;Z (cm);R (cm);Weighted Counts", 100, -150, 150, 200, 0, 200);
0130     auto h2_zr_clust_innr = new TH2D("h2_zr_clust_innr", "clust energy-weighted RZ distribution;Z (cm);R (cm);Weighted Counts", 100, -150, 150, 200, 0, 200);
0131
0132     auto h2_corr_dR_dphiR = new TH2D("h2_corr_dR_dphiR", "h2_corr_dR_dphiR; dR; dphiR; Weighted Counts", 30, -150, 150, 100, -5, 5);
0133
0134     auto h2_weird_XY = new TH2D("h2_weird_XY", "h2_weird; X; Y; Weighted Counts", 300, -150, 150, 300, -150, 150);
0135     auto h2_weird_pt_E = new TH2D("h2_weird_pt_E", "h2_weird_pt_E; pt(GeV); E(GeV); Weighted Counts", 10, 0, 10, 20, 0, 20);
0136     auto h1_weird_pt = new TH1D("h1_weird_pt", "h1_weird_pt; pt(GeV); Weighted Counts", 50, 0, 10);
0137     auto h1_weird_E  = new TH1D("h1_weird_E", "h1_weird_E; E(GeV); Weighted Counts", 100, 0, 20);
0138
0139     // delta phi distribution
0140     auto h1_dphi_distri_pr1_clustinnr = new TH1D(" h1_dphi_distri_pr1_clustinnr", "h1_dphi_distri_pr1_clustinnr; dphi; Counts", 50, -0.02, 0.02);
0141     auto h1_dphi_distri_pr1_clustgeom = new TH1D(" h1_dphi_distri_pr1_clustgeom", "h1_dphi_distri_pr1_clustgeom; dphi; Counts", 50, -0.02, 0.02);
0142     auto h1_dphi_distri_g4h_clustinnr = new TH1D(" h1_dphi_distri_g4h_clustinnr", "h1_dphi_distri_g4h_clustinnr; dphi; Counts", 50, -0.02, 0.02);
0143     auto h1_dphi_distri_g4h_clustgeom = new TH1D(" h1_dphi_distri_g4h_clustgeom", "h1_dphi_distri_g4h_clustgeom; dphi; Counts", 50, -0.02, 0.02);
0144
0145     auto h2_dphi_distri_test = new TH2D("h2_dphi_distri_test", "h2_dphi_distri_test; pt(GeV); dphi", 55, 0, 11, 40, -0.02, 0.02);
0146     auto h2_dphi_Pr_innr_Etruth = new TH2D("h2_dphi_Pr_innr_Etruth", "h2_dphi_Pr_innr_Etruth; E(GeV); dphi", 100, 0, 20, 60, -0.03, 0.03);
0147     auto h2_dphi_Pr_geom_Etruth = new TH2D("h2_dphi_Pr_geom_Etruth", "h2_dphi_Pr_geom_Etruth; E(GeV); dphi", 100, 0, 20, 60, -0.03, 0.03);
0148     auto h2_dphi_Pr_innr_Ereco = new TH2D("h2_dphi_Pr_innr_Ereco", "h2_dphi_Pr_innr_Ereco; E(GeV); dphi", 100, 0, 20, 60, -0.03, 0.03);
0149     auto h2_dphi_Pr_geom_Ereco = new TH2D("h2_dphi_Pr_geom_Ereco", "h2_dphi_Pr_geom_Ereco; E(GeV); dphi", 100, 0, 20, 60, -0.03, 0.03);
0150
0151     auto h2_dphi_Pr_innr_eta = new TH2D("h2_dphi_Pr_innr_eta", "h2_dphi_Pr_innr_eta; #eta; dphi", 120, -1.2, 1.2, 60, -0.03, 0.03);
0152     auto h2_dphi_Pr_geom_eta = new TH2D("h2_dphi_Pr_geom_eta", "h2_dphi_Pr_geom_eta; #eta; dphi", 120, -1.2, 1.2, 60, -0.03, 0.03);
0153
0154     // delta R distribution
0155     auto h2_dRadius_Pr_innr_Ereco = new TH2D("h2_dRadius_Pr_innr_Ereco", "h2_dRadius_Pr_innr_Ereco; E(GeV); dRadius", 100, 0, 20, 200, -10, 10);
0156     auto h2_dRadius_Pr_geom_Ereco = new TH2D("h2_dRadius_Pr_geom_Ereco", "h2_dRadius_Pr_geom_Ereco; E(GeV); dRadius", 100, 0, 20, 200, -15, 5);
0157
0158     // TProfile setting
0159     auto TP_DR_pr1_clustinnr_pt   = new TProfile("TP_DR_pr1_clustinnr_pt", "TP_DR_pr1_clustinnr_pt;pt(GeV);R(cm)", 11, 0, 11, 0, 20);
0160     auto TP_DR_g4hit_clustgeom_pt = new TProfile("TP_DR_g4hit_clustgeom_pt", "TP_DR_g4hit_clustgeom_pt;pt(GeV);R(cm)", 11, -0.5, 10.5, 0, 20);
0161     auto TP_DR_pr1_clustinnr_E    = new TProfile("TP_DR_pr1_clustinnr_E", "TP_DR_pr1_clustinnr_E;E(GeV);R(cm)", 15, -0.5, 14.5, 0, 20);
0162     auto TP_DR_g4hit_clustgeom_E  = new TProfile("TP_DR_g4hit_clustgeom_E", "TP_DR_g4hit_clustgeom_E;E(GeV);R(cm)", 15, -0.5, 14.5, 0, 20);
0163
0164     // check eta/z dependence
0165     auto h2_z_dphi_distri_test = new TH2D("h2_z_dphi_distri_test", "h2_z_dphi_distri_test; Z(cm); dphi(rad)",30, -150, 150, 40, -0.02, 0.02);
0166     auto h2_z_dRad_distri_test = new TH2D("h2_z_dRad_distri_test", "h2_z_dRad_distri_test; Z(cm); dRad(cm)", 30, -150, 150, 50, -5, 5);
0167     auto h2_z_dDis_distri_test = new TH2D("h2_z_dDis_distri_test", "h2_z_dDis_distri_test; Z(cm); dDis(cm)", 30, -150, 150, 50, -5, 5);
0168
0169     TGraph* g_cor_z_dphi_distri_test = nullptr;
0170     TGraph* g_cor_z_dRad_distri_test = nullptr;
0171     TGraph* g_cor_z_dDis_distri_test = nullptr;
0172
0173     // 定义 TGraph
0174     TGraph *g_pr1hit = new TGraph();
0175     g_pr1hit->SetMarkerStyle(20);  // 实心圆点
0176     g_pr1hit->SetMarkerColor(kRed); // 红色
0177     g_pr1hit->SetMarkerSize(0.5);  // 设置红色点大小
0178
0179     TGraph *g_clustinnr = new TGraph();
0180     g_clustinnr->SetMarkerStyle(21);  // 实心方点
0181     g_clustinnr->SetMarkerColor(kBlue); // 蓝色
0182     g_clustinnr->SetMarkerSize(0.5);  // 设置蓝色点大小
0183
0184     int pointIndex = 0;
0185
0186     // 设置分支地址
0187     caloinfotree->SetBranchAddress("CEMC_Pr_Hit_x", &CEMC_Pr_Hit_x);
0188     caloinfotree->SetBranchAddress("CEMC_Pr_Hit_y", &CEMC_Pr_Hit_y);
0189     caloinfotree->SetBranchAddress("CEMC_Pr_Hit_z", &CEMC_Pr_Hit_z);
0190     caloinfotree->SetBranchAddress("CEMC_Pr_Hit_R", &CEMC_Pr_Hit_R);
0191
0192     caloinfotree->SetBranchAddress("CEMC_Hit_x", &CEMC_Hit_x);
0193     caloinfotree->SetBranchAddress("CEMC_Hit_y", &CEMC_Hit_y);
0194     caloinfotree->SetBranchAddress("CEMC_Hit_z", &CEMC_Hit_z);
0195     // caloinfotree->SetBranchAddress("CEMC_Hit_R", &CEMC_Hit_R);
0196     caloinfotree->SetBranchAddress("CEMC_Hit_Edep", &CEMC_Hit_Edep);
0197
0198     caloinfotree->SetBranchAddress("tower_system", &tower_system);
0199     caloinfotree->SetBranchAddress("tower_X", &tower_x);
0200     caloinfotree->SetBranchAddress("tower_Y", &tower_y);
0201     caloinfotree->SetBranchAddress("tower_Z", &tower_z);
0202     caloinfotree->SetBranchAddress("tower_R", &tower_R);
0203     caloinfotree->SetBranchAddress("tower_Phi", &tower_phi);
0204     caloinfotree->SetBranchAddress("tower_Eta", &tower_eta);
0205     caloinfotree->SetBranchAddress("tower_edep", &tower_edep);
0206
0207     caloinfotree->SetBranchAddress("tower_int_X", &tower_innr_x);
0208     caloinfotree->SetBranchAddress("tower_int_Y", &tower_innr_y);
0209     caloinfotree->SetBranchAddress("tower_int_Z", &tower_innr_z);
0210     caloinfotree->SetBranchAddress("tower_int_R", &tower_innr_R);
0211     // caloinfotree->SetBranchAddress("tower_int_Phi", &tower_innr_phi);
0212     // caloinfotree->SetBranchAddress("tower_int_Eta", &tower_innr_eta);
0213     // caloinfotree->SetBranchAddress("tower_int_edep", &tower_innr_edep);
0214
0215     caloinfotree->SetBranchAddress("caloClus_system", &caloClus_system);
0216     caloinfotree->SetBranchAddress("caloClus_X", &caloClus_x);
0217     caloinfotree->SetBranchAddress("caloClus_Y", &caloClus_y);
0218     caloinfotree->SetBranchAddress("caloClus_Z", &caloClus_z);
0219     caloinfotree->SetBranchAddress("caloClus_R", &caloClus_R);
0220     caloinfotree->SetBranchAddress("caloClus_Phi", &caloClus_phi);
0221     caloinfotree->SetBranchAddress("caloClus_edep", &caloClus_edep);
0222
0223     caloinfotree->SetBranchAddress("caloClus_innr_X", &caloClus_innr_x);
0224     caloinfotree->SetBranchAddress("caloClus_innr_Y", &caloClus_innr_y);
0225     caloinfotree->SetBranchAddress("caloClus_innr_Z", &caloClus_innr_z);
0226     caloinfotree->SetBranchAddress("caloClus_innr_R", &caloClus_innr_R);
0227     caloinfotree->SetBranchAddress("caloClus_innr_Phi", &caloClus_innr_phi);
0228     caloinfotree->SetBranchAddress("caloClus_innr_edep", &caloClus_innr_edep);
0229
0230     // moment pt and energy E
0231     caloinfotree->SetBranchAddress("PrimaryG4P_Pt", &PrimaryG4P_Pt);
0232     caloinfotree->SetBranchAddress("PrimaryG4P_E", &PrimaryG4P_E);
0233
0234     int num_of_clus = 0;
0235
0236     // 遍历树的条目
0237     Long64_t nentries = caloinfotree->GetEntries();
0238     // for (Long64_t i = 0; i < nentries; i++)
0239     for (Long64_t i = 0; i < nentries; i++)
0240     {
0241         caloinfotree->GetEntry(i);
0242
0243         num_of_clus = 0;
0244
0245         // CEMC_Pr_Hit ---------------------------------------------------------------------------------------
0246         if(CEMC_Pr_Hit_x->size() != 1)
0247         {
0248             // std::cerr << "Error: CEMC_Pr_Hit_x size is: "<< CEMC_Pr_Hit_x->size() << std::endl;
0249             continue;
0250         }
0251         double PrHit_x = CEMC_Pr_Hit_x->at(0);
0252         double PrHit_y = CEMC_Pr_Hit_y->at(0);
0253         double PrHit_z = CEMC_Pr_Hit_z->at(0);
0254         double PrHit_R = CEMC_Pr_Hit_R->at(0);
0255         h2_zr_g4hit->Fill(PrHit_z, PrHit_R);
0256         double pr1cemc_xyz[3] = {PrHit_x, PrHit_y, PrHit_z};
0257         double PrHit_phi = TMath::ATan2(PrHit_y, PrHit_x);
0258         if (PrHit_phi < 0)
0259         {
0260             PrHit_phi += 2 * TMath::Pi();
0261         }
0262
0263         // 遍历 CEMC_Hit vector 中的每个元素 --------------------------------------------------------------------
0264         Double_t g4hit_TotEMCalE = 0.;
0265         Double_t g4hit_ModifEMCal_x = 0.;
0266         Double_t g4hit_ModifEMCal_y = 0.;
0267         Double_t g4hit_ModifEMCal_z = 0.;
0268
0269         for (size_t j = 0; j < CEMC_Hit_x->size(); j++)
0270         {
0271             double hitval_x = CEMC_Hit_x->at(j);
0272             double hitval_y = CEMC_Hit_y->at(j);
0273             double hitval_z = CEMC_Hit_z->at(j);
0274             double hitval_R = TMath::Sqrt(hitval_x * hitval_x + hitval_y * hitval_y);
0275             double hitval_phi = TMath::ATan2(hitval_y, hitval_x);
0276             double hitval_eta = TMath::ATan2(hitval_R, hitval_z);
0277             double hitval_edep = CEMC_Hit_Edep->at(j);
0278
0279             if(hitval_R<200)
0280             {
0281                 g4hit_TotEMCalE += hitval_edep;
0282
0283                 g4hit_ModifEMCal_x += hitval_edep * hitval_x;
0284                 g4hit_ModifEMCal_y += hitval_edep * hitval_y;
0285                 g4hit_ModifEMCal_z += hitval_edep * hitval_z;
0286
0287                 // h2_zr_g4hit->Fill(hitval_z, hitval_R, hitval_edep);
0288             }
0289         }
0290         g4hit_ModifEMCal_x /= g4hit_TotEMCalE;
0291         g4hit_ModifEMCal_y /= g4hit_TotEMCalE;
0292         g4hit_ModifEMCal_z /= g4hit_TotEMCalE;
0293         double g4hit_xyz[3] = {g4hit_ModifEMCal_x, g4hit_ModifEMCal_y, g4hit_ModifEMCal_z};
0294
0295         double g4hit_ModifEMCal_phi = TMath::ATan2(g4hit_ModifEMCal_y, g4hit_ModifEMCal_x);
0296         h_g4hitphi->Fill(g4hit_ModifEMCal_phi);
0297         if (g4hit_ModifEMCal_phi < 0)
0298         {
0299             g4hit_ModifEMCal_phi += 2 * TMath::Pi();
0300         }
0301
0302         // 遍历 tower vector 中的每个元素 -----------------------------------------------------------------------------------
0303         Double_t tower_TotEMCalE = 0.;
0304         Double_t tower_ModifEMCal_x = 0.;
0305         Double_t tower_ModifEMCal_y = 0.;
0306         Double_t tower_ModifEMCal_z = 0.;
0307
0308         Double_t tower_innr_ModifEMCal_x = 0.;
0309         Double_t tower_innr_ModifEMCal_y = 0.;
0310         Double_t tower_innr_ModifEMCal_z = 0.;
0311         int j_innr = 0;
0312
0313         for (size_t j = 0; j < tower_x->size(); j++)
0314         {
0315             if ((tower_system->at(j)) != 0)
0316             {
0317                 j_innr = j_innr - 1;
0318                 continue;
0319             }
0320             double towerval_x = tower_x->at(j);
0321             double towerval_y = tower_y->at(j);
0322             double towerval_z = tower_z->at(j);
0323             double towerval_R = tower_R->at(j);
0324             double towerval_phi = tower_phi->at(j);
0325             double towerval_eta = tower_eta->at(j);
0326             double towerval_edep = tower_edep->at(j);
0327
0328             double towerval_innr_x = tower_innr_x->at(j_innr);
0329             double towerval_innr_y = tower_innr_y->at(j_innr);
0330             double towerval_innr_z = tower_innr_z->at(j_innr);
0331             j_innr += 1;
0332             double towerval_innr_R = sqrt(towerval_innr_x*towerval_innr_x + towerval_innr_y*towerval_innr_y);
0333
0334             double Ecalo_threshold = 0.07;
0335
0336             if(towerval_edep > Ecalo_threshold)
0337             {
0338                 tower_TotEMCalE += towerval_edep;
0339
0340                 // tower geom center
0341                 tower_ModifEMCal_x += towerval_edep * towerval_x;
0342                 tower_ModifEMCal_y += towerval_edep * towerval_y;
0343                 tower_ModifEMCal_z += towerval_edep * towerval_z;
0344
0345                 h2_zr_tower->Fill(towerval_z, towerval_R, towerval_edep);
0346
0347                 // innr face center
0348                 tower_innr_ModifEMCal_x += towerval_edep * towerval_innr_x;
0349                 tower_innr_ModifEMCal_y += towerval_edep * towerval_innr_y;
0350                 tower_innr_ModifEMCal_z += towerval_edep * towerval_innr_z;
0351
0352                 // cout<<"towerval_x is: "<<towerval_x<<", "<<towerval_y<<", "<<towerval_z<<endl;
0353                 // cout<<"towerval_innr_x is: "<<towerval_innr_x<<", "<<towerval_innr_y<<", "<<towerval_innr_z<<endl;
0354
0355                 h2_zr_tower_innr->Fill(towerval_innr_z, towerval_innr_R, towerval_edep);
0356             }
0357         }
0358         // geom center
0359         tower_ModifEMCal_x /= tower_TotEMCalE;
0360         tower_ModifEMCal_y /= tower_TotEMCalE;
0361         tower_ModifEMCal_z /= tower_TotEMCalE;
0362         double towergeom_xyz[3] = {tower_ModifEMCal_x, tower_ModifEMCal_y, tower_ModifEMCal_z};
0363
0364         double tower_ModifEMCal_phi = TMath::ATan2(tower_ModifEMCal_y, tower_ModifEMCal_x);
0365         h_towerphi->Fill(tower_ModifEMCal_phi);
0366
0367         // innr center
0368         tower_innr_ModifEMCal_x /= tower_TotEMCalE;
0369         tower_innr_ModifEMCal_y /= tower_TotEMCalE;
0370         tower_innr_ModifEMCal_z /= tower_TotEMCalE;
0371         double towerinnr_xyz[3] = {tower_innr_ModifEMCal_x, tower_innr_ModifEMCal_y, tower_innr_ModifEMCal_z};
0372
0373         double tower_innr_ModifEMCal_phi = TMath::ATan2(tower_innr_ModifEMCal_y, tower_innr_ModifEMCal_x);
0374         h_innr_towerphi->Fill(tower_innr_ModifEMCal_phi);
0375
0376         // 遍历 cluster innersurface vector 中的每个元素 --------------------------------------------------------------------
0377         Double_t cluster_innr_TotEMCalE = 0.;
0378         Double_t cluster_innr_ModifEMCal_x = 0.;
0379         Double_t cluster_innr_ModifEMCal_y = 0.;
0380         Double_t cluster_innr_ModifEMCal_z = 0.;
0381
0382         for (size_t j = 0; j < caloClus_innr_x->size(); j++)
0383         {
0384             double clusterval_innr_x = caloClus_innr_x->at(j);
0385             double clusterval_innr_y = caloClus_innr_y->at(j);
0386             double clusterval_innr_z = caloClus_innr_z->at(j);
0387             double clusterval_innr_R = caloClus_innr_R->at(j);
0388             double clusterval_innr_phi = caloClus_innr_phi->at(j);
0389             double clusterval_innr_edep = caloClus_innr_edep->at(j);
0390
0391             double Ecalo_threshold = 0.5;
0392             if((clusterval_innr_edep > Ecalo_threshold)&&(clusterval_innr_R<140))
0393             {
0394                 cluster_innr_ModifEMCal_x += clusterval_innr_x;
0395                 cluster_innr_ModifEMCal_y += clusterval_innr_y;
0396                 cluster_innr_ModifEMCal_z += clusterval_innr_z;
0397                 cluster_innr_TotEMCalE += clusterval_innr_edep;
0398
0399                 num_of_clus += 1;
0400                 // cout<<"cluster_innr_xyz: "<<clusterval_innr_x<<" "<<clusterval_innr_y<<" "<<clusterval_innr_z<<endl;
0401
0402                 h2_zr_clust_innr->Fill(clusterval_innr_z, clusterval_innr_R, clusterval_innr_edep);
0403             }
0404         }
0405         cluster_innr_ModifEMCal_x /= cluster_innr_TotEMCalE;
0406         cluster_innr_ModifEMCal_y /= cluster_innr_TotEMCalE;
0407         cluster_innr_ModifEMCal_z /= cluster_innr_TotEMCalE;
0408         double cluster_innr_phi = TMath::ATan2(cluster_innr_ModifEMCal_y, cluster_innr_ModifEMCal_x);
0409
0410         // double dphi_C = g_dphiM_z->Eval(cluster_innr_ModifEMCal_z);
0411         // double dRad_C = g_dRadM_z->Eval(cluster_innr_ModifEMCal_z);
0412         // double dphi_C = f1_dphi_C->Eval(cluster_innr_TotEMCalE);
0413         double dphi_C = g1_dphi_C->Eval(cluster_innr_TotEMCalE);
0414
0415         if (cluster_innr_phi < 0)
0416         {
0417             cluster_innr_phi += 2 * TMath::Pi();
0418         }
0419         // cluster_innr_phi = cluster_innr_phi + dphi_C; // ave = - 0.0083
0420
0421         double R_temp = sqrt(cluster_innr_ModifEMCal_x*cluster_innr_ModifEMCal_x + cluster_innr_ModifEMCal_y*cluster_innr_ModifEMCal_y);
0422         // R_temp = R_temp + dRad_C; // ave = 0.4
0423         cluster_innr_ModifEMCal_x = R_temp*cos(cluster_innr_phi);
0424         cluster_innr_ModifEMCal_y = R_temp*sin(cluster_innr_phi);
0425
0426         double clustinnr_xyz[3] = {cluster_innr_ModifEMCal_x, cluster_innr_ModifEMCal_y, cluster_innr_ModifEMCal_z};
0427
0428         // 遍历 cluster vector 中的每个元素  --------------------------------------------------------------------
0429         Double_t cluster_TotEMCalE = 0.;
0430         Double_t cluster_ModifEMCal_x = 0.;
0431         Double_t cluster_ModifEMCal_y = 0.;
0432         Double_t cluster_ModifEMCal_z = 0.;
0433
0434         for (size_t j = 0; j < caloClus_x->size(); j++)
0435         {
0436             if ((caloClus_system->at(j)) != 0) continue;
0437
0438             double clusterval_x = caloClus_x->at(j);
0439             double clusterval_y = caloClus_y->at(j);
0440             double clusterval_z = caloClus_z->at(j);
0441             double clusterval_R = sqrt(clusterval_x*clusterval_x + clusterval_y*clusterval_y);
0442             double clusterval_phi = caloClus_phi->at(j);
0443             double clusterval_edep = caloClus_edep->at(j);
0444
0445             double Ecalo_threshold = 0.3;
0446             if((clusterval_edep > Ecalo_threshold)&&(clusterval_R<200))
0447             {
0448                 h_clustphi->Fill(clusterval_phi);
0449                 cluster_ModifEMCal_x += clusterval_x;
0450                 cluster_ModifEMCal_y += clusterval_y;
0451                 cluster_ModifEMCal_z += clusterval_z;
0452                 cluster_TotEMCalE += 1;
0453
0454                 // cout<<"cluster_geom_xyz: "<<clusterval_x<<" "<<clusterval_y<<" "<<clusterval_z<<endl;
0455
0456                 h2_zr_clust->Fill(clusterval_z, clusterval_R, clusterval_edep);
0457             }
0458         }
0459         cluster_ModifEMCal_x /= cluster_TotEMCalE;
0460         cluster_ModifEMCal_y /= cluster_TotEMCalE;
0461         cluster_ModifEMCal_z /= cluster_TotEMCalE;
0462         double clustgeom_xyz[3]  = {cluster_ModifEMCal_x, cluster_ModifEMCal_y, cluster_ModifEMCal_z};
0463
0464         double cluster_geom_phi = TMath::ATan2(cluster_ModifEMCal_y, cluster_ModifEMCal_x);
0465         if (cluster_geom_phi < 0)
0466         {
0467             cluster_geom_phi += 2 * TMath::Pi();
0468         }
0469
0470         // calculate delta truth tower cluster  --------------------------------------------------------------------
0471         double DR_pr1_clustinnr = Cal_Distance_R(pr1cemc_xyz, clustinnr_xyz);
0472         double DR_pr1_towerinnr = Cal_Distance_R(pr1cemc_xyz, towerinnr_xyz);
0473         double DR_g4hit_clustgeom = Cal_Distance_R(g4hit_xyz, clustgeom_xyz);
0474         double DR_g4hit_towergeom = Cal_Distance_R(g4hit_xyz, towergeom_xyz);
0475
0476         double dRadius_pr1_clustinnr   = Cal_dR(pr1cemc_xyz, clustinnr_xyz);
0477         double dRadius_g4hit_clustgeom = Cal_dR(g4hit_xyz, clustgeom_xyz);
0478
0479         double dphiR_pr1_clustinnr   = Cal_dphiR(pr1cemc_xyz, clustinnr_xyz);
0480         double dphiR_g4hit_clustgeom = Cal_dphiR(g4hit_xyz, clustgeom_xyz);
0481
0482         // distance between reco and truth
0483         h1_DR_pr1_clustinnr->Fill(DR_pr1_clustinnr);
0484         h1_DR_pr1_towerinnr->Fill(DR_pr1_towerinnr);
0485         h1_DR_g4hit_clustgeom->Fill(DR_g4hit_clustgeom);
0486         h1_DR_g4hit_towergeom->Fill(DR_g4hit_towergeom);
0487
0488         // num of clus checking
0489         if(DR_pr1_clustinnr>5)
0490         {
0491             // cout<<"DR_pr1_clustinnr is: "<< DR_pr1_clustinnr <<", number of clus is: "<< num_of_clus <<endl;
0492             if(num_of_clus==1)
0493             {
0494                 double clus_R_tem = sqrt(clustinnr_xyz[0]*clustinnr_xyz[0] + clustinnr_xyz[1]*clustinnr_xyz[1]);
0495                 h2_weird_XY->Fill(clustinnr_xyz[2], clus_R_tem);
0496                 h2_weird_pt_E->Fill(PrimaryG4P_Pt->at(0), PrimaryG4P_E->at(0));
0497                 h1_weird_pt->Fill(PrimaryG4P_Pt->at(0));
0498                 h1_weird_E ->Fill(PrimaryG4P_E->at(0));
0499
0500                 // 添加点
0501                 g_pr1hit->SetPoint(pointIndex, pr1cemc_xyz[0], pr1cemc_xyz[1]);
0502                 g_clustinnr->SetPoint(pointIndex, clustinnr_xyz[0], clustinnr_xyz[1]);
0503                 pointIndex++;
0504             }
0505         }
0506
0507         h1_dRadius_pr1_clustinnr  ->Fill(dRadius_pr1_clustinnr);
0508         h1_dRadius_g4hit_clustgeom->Fill(dRadius_g4hit_clustgeom);
0509
0510         PrHit_phi - cluster_innr_phi > 0 ? h1_dphiR_pr1_clustinnr->Fill(dphiR_pr1_clustinnr) : h1_dphiR_pr1_clustinnr->Fill(-dphiR_pr1_clustinnr);
0511         // h1_dphiR_pr1_clustinnr    ->Fill(dphiR_pr1_clustinnr);
0512         h1_dphiR_g4hit_clustgeom  ->Fill(dphiR_g4hit_clustgeom);
0513
0514         PrHit_phi - cluster_innr_phi > 0 ? h2_corr_dR_dphiR->Fill(dRadius_pr1_clustinnr,dphiR_pr1_clustinnr) : h2_corr_dR_dphiR->Fill(dRadius_pr1_clustinnr,-dphiR_pr1_clustinnr);
0515
0516         // right or left?
0517         // PrHit_phi - cluster_innr_phi > 0 ? h1_DR_pr1_clustinnr->Fill(DR_pr1_clustinnr) : h1_DR_pr1_clustinnr->Fill(-DR_pr1_clustinnr);
0518         // g4hit_ModifEMCal_phi - cluster_innr_phi > 0 ? h1_DR_g4hit_clustgeom->Fill(DR_g4hit_clustgeom) : h1_DR_g4hit_clustgeom->Fill(-DR_g4hit_clustgeom);
0519         h1_dphi_distri_pr1_clustinnr->Fill(PrHit_phi - cluster_innr_phi);
0520         h1_dphi_distri_pr1_clustgeom->Fill(PrHit_phi - cluster_geom_phi);
0521         h1_dphi_distri_g4h_clustinnr->Fill(g4hit_ModifEMCal_phi - cluster_innr_phi);
0522         h1_dphi_distri_g4h_clustgeom->Fill(g4hit_ModifEMCal_phi - cluster_geom_phi);
0523
0524         // h2_dphi_distri_test->Fill(abs(cluster_innr_ModifEMCal_z),(PrHit_phi - cluster_innr_phi));
0525
0526         // 2d map
0527         double pr_electron_pt = PrimaryG4P_Pt->at(0);
0528         double pr_electron_E  = PrimaryG4P_E->at(0);
0529
0530         h2_DR_pr1_clustinnr_pt  ->Fill(pr_electron_pt, DR_pr1_clustinnr);
0531         h2_DR_g4hit_clustgeom_pt->Fill(pr_electron_pt, DR_g4hit_clustgeom);
0532         h2_DR_pr1_clustinnr_E   ->Fill(pr_electron_E, DR_pr1_clustinnr);
0533         h2_DR_g4hit_clustgeom_E ->Fill(pr_electron_E, DR_g4hit_clustgeom);
0534
0535         TP_DR_pr1_clustinnr_pt  ->Fill(pr_electron_pt, DR_pr1_clustinnr);
0536         TP_DR_g4hit_clustgeom_pt->Fill(pr_electron_pt, DR_g4hit_clustgeom);
0537         TP_DR_pr1_clustinnr_E   ->Fill(pr_electron_E, DR_pr1_clustinnr);
0538         TP_DR_g4hit_clustgeom_E ->Fill(pr_electron_E, DR_g4hit_clustgeom);
0539
0540         // 2d map for dphi
0541         h2_dphi_distri_test->Fill(pr_electron_pt,(PrHit_phi - cluster_innr_phi));
0542         h2_dphi_Pr_innr_Etruth->Fill(pr_electron_E,(PrHit_phi - cluster_innr_phi));
0543         h2_dphi_Pr_geom_Etruth->Fill(pr_electron_E,(PrHit_phi - cluster_geom_phi));
0544
0545         if (caloClus_innr_edep->size() == 1)
0546         {
0547             // 计算 inner‐surface cluster 的 pseudorapidity and 计算 geometric cluster 的 pseudorapidity
0548             double clust_innr_x = caloClus_innr_x->at(0);
0549             double clust_innr_y = caloClus_innr_y->at(0);
0550             double clust_innr_z = caloClus_innr_z->at(0);
0551             double clust_innr_R = std::sqrt(clust_innr_x*clust_innr_x + clust_innr_y*clust_innr_y);
0552             double clust_innr_eta = std::asinh(clust_innr_z / clust_innr_R);
0553
0554             double clust_geo_x = caloClus_x->at(0);
0555             double clust_geo_y = caloClus_y->at(0);
0556             double clust_geo_z = caloClus_z->at(0);
0557             double clust_geo_R = std::sqrt(clust_geo_x*clust_geo_x + clust_geo_y*clust_geo_y);
0558             double clust_geo_eta  = std::asinh(clust_geo_z  / clust_geo_R);
0559
0560             double reco_electron_E = caloClus_innr_edep->at(0);
0561             h2_dphi_Pr_innr_Ereco->Fill(reco_electron_E,(PrHit_phi - cluster_innr_phi));
0562             h2_dphi_Pr_geom_Ereco->Fill(reco_electron_E,(PrHit_phi - cluster_geom_phi));
0563
0564             h2_dphi_Pr_innr_eta->Fill(clust_innr_eta,(PrHit_phi - cluster_innr_phi));
0565             h2_dphi_Pr_geom_eta->Fill(clust_geo_eta,(PrHit_phi - cluster_geom_phi));
0566
0567             h2_dRadius_Pr_innr_Ereco->Fill(reco_electron_E,(PrHit_R - clust_innr_R));
0568             h2_dRadius_Pr_geom_Ereco->Fill(reco_electron_E,(PrHit_R - clust_geo_R));
0569         }
0570
0571         // z dependent check
0572         h2_z_dphi_distri_test->Fill(cluster_innr_ModifEMCal_z,(PrHit_phi - cluster_innr_phi));
0573         h2_z_dRad_distri_test->Fill(cluster_innr_ModifEMCal_z,dRadius_pr1_clustinnr);
0574         h2_z_dDis_distri_test->Fill(cluster_innr_ModifEMCal_z,DR_pr1_clustinnr);
0575
0576         // g_cor_z_dphi_distri_test = Cor_var_Z(h2_z_dphi_distri_test, -0.012, -0.004);
0577         // g_cor_z_dRad_distri_test = Cor_var_Z(h2_z_dRad_distri_test, -1, 2);
0578         // g_cor_z_dDis_distri_test = Cor_var_Z(h2_z_dDis_distri_test);
0579     }
0580     cout<<"pointIndex is: "<<pointIndex<<endl;
0581
0582     // 定义一个用于二次曲线拟合的 TF1
0583     TF1* fitCurve_dphi = new TF1("fitCurve_dphi", "[0] + [1]*x + [2]*x*x + [3]*x*x*x + [4]*x*x*x*x + [5]*x*x*x*x", 0.5, 16);
0584     TF1* projectionFit_dphi = new TF1("projectionFit_dphi","gaus", -1, 1);
0585     TGraph* g1fitPeakVsX_dphi = FitPeakVsX(h2_dphi_Pr_innr_Ereco, projectionFit_dphi, fitCurve_dphi, "g1fitPeakVsX_dphi");
0586
0587     TF1* fitCurve_dRadius = new TF1("fitCurve_dRadius", "[0] + [1]*x + [2]*x*x + [3]*x*x*x + [4]*x*x*x*x + [5]*x*x*x*x", 0.5, 16);
0588     TF1* projectionFit_dRadius = new TF1("projectionFit_dRadius","gaus", -1, 1);
0589     TGraph* g1fitPeakVsX_dRadius = FitPeakVsX(h2_dRadius_Pr_innr_Ereco, projectionFit_dRadius, fitCurve_dRadius, "g1fitPeakVsX_dRadius");
0590
0591     // 保存直方图到文件
0592     TFile outfile(fOutputName.c_str(), "RECREATE");
0593     outfile.cd();
0594
0595     h_g4hitphi->Write();
0596     h_towerphi->Write();
0597     h_innr_towerphi->Write();
0598
0599     h1_DR_pr1_clustinnr->Write();
0600     h1_DR_pr1_towerinnr->Write();
0601     h1_DR_g4hit_clustgeom->Write();
0602     h1_DR_g4hit_towergeom->Write();
0603
0604     h1_dRadius_pr1_clustinnr  ->Write();
0605     h1_dRadius_g4hit_clustgeom->Write();
0606     h1_dphiR_pr1_clustinnr    ->Write();
0607     h1_dphiR_g4hit_clustgeom  ->Write();
0608
0609     h2_zr_g4hit->Write();
0610     TH1D* h1_zr_g4Pr = h2_zr_g4hit->ProjectionY("h1_zr_g4Pr");
0611     h1_zr_g4Pr->Write();
0612     h2_zr_tower->Write();
0613     h2_zr_clust->Write();
0614
0615     h2_zr_tower_innr->Write();
0616     h2_zr_clust_innr->Write();
0617
0618     h2_DR_pr1_clustinnr_pt->Write();
0619     h2_DR_g4hit_clustgeom_pt->Write();
0620     h2_DR_pr1_clustinnr_E->Write();
0621     h2_DR_g4hit_clustgeom_E->Write();
0622
0623     // profile write
0624     TP_DR_pr1_clustinnr_pt  ->Write();
0625     TP_DR_g4hit_clustgeom_pt->Write();
0626     TP_DR_pr1_clustinnr_E   ->Write();
0627     TP_DR_g4hit_clustgeom_E ->Write();
0628
0629     // weird hit pos
0630     h2_weird_XY->Write();
0631     h2_weird_pt_E->Write();
0632     h1_weird_pt->Write();
0633     h1_weird_E->Write();
0634
0635     h1_dphi_distri_pr1_clustinnr->Write();
0636     h1_dphi_distri_pr1_clustgeom->Write();
0637     h1_dphi_distri_g4h_clustinnr->Write();
0638     h1_dphi_distri_g4h_clustgeom->Write();
0639
0640     h2_dphi_distri_test->Write();
0641     h2_dphi_Pr_innr_Etruth->Write();
0642     h2_dphi_Pr_geom_Etruth->Write();
0643     h2_dphi_Pr_innr_Ereco->Write();
0644     h2_dphi_Pr_geom_Ereco->Write();
0645
0646     h2_dphi_Pr_innr_eta->Write();
0647     h2_dphi_Pr_geom_eta->Write();
0648
0649     h2_dRadius_Pr_innr_Ereco->Write();
0650     h2_dRadius_Pr_geom_Ereco->Write();
0651
0652     h2_corr_dR_dphiR->Write();
0653
0654     h2_z_dphi_distri_test->Write();
0655     h2_z_dRad_distri_test->Write();
0656     h2_z_dDis_distri_test->Write();
0657
0658     g1fitPeakVsX_dphi->Write();
0659     fitCurve_dphi->Write();
0660     g1fitPeakVsX_dRadius->Write();
0661     fitCurve_dRadius->Write();
0662
0663     // g_cor_z_dphi_distri_test->Write("g_cor_z_dphi_distri_test");
0664     // g_cor_z_dRad_distri_test->Write("g_cor_z_dRad_distri_test");
0665     // g_cor_z_dDis_distri_test->Write("g_cor_z_dDis_distri_test");
0666
0667     // Draw the results
0668     TCanvas *c1 = new TCanvas("c1", "Hit Position", 800, 600);
0669     g_pr1hit->Draw("AP");
0670     g_clustinnr->Draw("P SAME"); // 在同一图层绘制
0671     c1->SaveAs("pr1hit_vs_clustinnr_large.pdf");
0672
0673     std::cout << "Histograms saved to " << fOutputName << std::endl;
0674 }
0675
0676 double Cal_Distance_R(const double xyz_1[3], const double xyz_2[3])
0677 {
0678     double delta_x = xyz_1[0] - xyz_2[0];
0679     double delta_y = xyz_1[1] - xyz_2[1];
0680     double delta_z = xyz_1[2] - xyz_2[2];
0681
0682     double distance_R = sqrt( delta_x*delta_x + delta_y*delta_y);
0683     // double distance_R = sqrt( delta_x*delta_x + delta_y*delta_y + delta_z*delta_z);
0684
0685     return distance_R;
0686 }
0687
0688 double Cal_dR(const double xyz_1[3], const double xyz_2[3])
0689 {
0690     double R_1 = sqrt(xyz_1[0]*xyz_1[0] + xyz_1[1]*xyz_1[1]);
0691     double R_2 = sqrt(xyz_2[0]*xyz_2[0] + xyz_2[1]*xyz_2[1]);
0692
0693     double dR = R_1 - R_2;
0694
0695     return dR;
0696 }
0697
0698 double Cal_dphiR(const double xyz_1[3], const double xyz_2[3])
0699 {
0700     double R_1 = sqrt(xyz_1[0]*xyz_1[0] + xyz_1[1]*xyz_1[1]);
0701     double R_2 = sqrt(xyz_2[0]*xyz_2[0] + xyz_2[1]*xyz_2[1]);
0702     double R_ratio = R_2 / R_1;
0703
0704     double x1 = xyz_1[0] * R_ratio;
0705     double y1 = xyz_1[1] * R_ratio;
0706     double x2 = xyz_2[0];
0707     double y2 = xyz_2[1];
0708
0709     double delta_x = x1 - x2;
0710     double delta_y = y1 - y2;
0711
0712     double dphiR = sqrt(delta_x*delta_x + delta_y*delta_y);
0713     // double dphiR = sqrt(x1*x1 + y1*y1) - sqrt(x2*x2 + y2*y2);
0714     // double dphiR = R_2 / R_1;
0715
0716
0717     return dphiR;
0718 }
0719
0720 TGraph* Cor_var_Z(TH2D* h2, double xrange_min = -0.05, double xrange_max = 0.05)
0721 {
0722     const int nBinsX = h2->GetNbinsX();  // Z 的 bin 数
0723     std::vector<double> z_vals;
0724     std::vector<double> dphi_means;
0725
0726     int binstep = 3; // bin 步长
0727     for (int i = 1; i <= nBinsX; i+=binstep)
0728     {
0729         double z1 = h2->GetXaxis()->GetBinCenter(i);
0730         double z2 = h2->GetXaxis()->GetBinCenter(i + 1);
0731         double z3 = h2->GetXaxis()->GetBinCenter(i + 2);
0732         double z_center = (z1 + z2 + z3)/binstep;
0733
0734         // 在每个 Z bin 上做投影（对 dphi 轴）
0735         TH1D* projY = h2->ProjectionY(Form("projY_%d", i), i, i + binstep - 1);
0736         if (projY->GetEntries() < 10)
0737         {
0738             delete projY;
0739             continue; // 跳过统计太少的 bin
0740         }
0741
0742         // TF1* cbfit = new TF1("cbfit", crystalball, xrange_min, xrange_max, 5);
0743         // // cbfit->SetParameters(1, 1.5, 2, 0, 0.01);
0744         // projY->Fit(cbfit, "RQ");
0745         // double dphi_peak = cbfit->GetParameter(3);
0746
0747         TF1* gausfit = new TF1("gausfit", "gaus", xrange_min, xrange_max);
0748         // gausfit->SetParameters(projY->GetMaximum(), projY->GetMean(), projY->GetRMS());
0749         projY->Fit(gausfit, "RQ");
0750         double dphi_peak = gausfit->GetParameter(1);
0751
0752         double dphi_mean = projY->GetMean();
0753
0754         z_vals.push_back(z_center);
0755         dphi_means.push_back(dphi_peak);
0756         delete projY; // 清理内存
0757     }
0758
0759     // 用 TGraph 存储 correction 曲线
0760     TGraph* gr_correction = new TGraph(z_vals.size(), &z_vals[0], &dphi_means[0]);
0761     gr_correction->SetTitle("Z-dependent var correction;Z (cm)");
0762     gr_correction->SetMarkerStyle(20);
0763
0764     return gr_correction;
0765 }
0766
0767 Double_t crystalball(Double_t* x, Double_t* par)
0768 {
0769     // par: [0]=norm, [1]=alpha, [2]=n, [3]=mean, [4]=sigma
0770     Double_t t = (x[0]-par[3])/par[4];
0771     if (par[1] < 0) t = -t;
0772     Double_t abs_alpha = fabs(par[1]);
0773     if (t >= -abs_alpha)
0774         return par[0]*exp(-0.5*t*t);
0775     else {
0776         Double_t a = pow(par[2]/abs_alpha, par[2]) * exp(-0.5*abs_alpha*abs_alpha);
0777         Double_t b = par[2]/abs_alpha - abs_alpha;
0778         return par[0]*a / pow(b - t, par[2]);
0779     }
0780 }
0781
0782
0783 TGraph* FitPeakVsX(TH2D* h2_XY,
0784                    TF1* projectionFit = nullptr, // 可选：对每个 projY 拟合用，若不用可设为 nullptr
0785                    TF1* peakVsXFit = nullptr,    // 用来拟合峰值随 x 变化的函数
0786                    const char* graphName = "grPeakVsX")
0787 {
0788     std::vector<double> vx, vy;
0789     int nBinsX = h2_XY->GetNbinsX();
0790
0791     for (int ix = 1; ix <= nBinsX; ++ix) {
0792         // 投影第 ix 个 bin 到 Y 轴
0793         TH1D* projY = h2_XY->ProjectionY(
0794             Form("projY_%d", ix),
0795             ix, ix, "e"
0796         );
0797         if (projY->GetEntries() < 500) {
0798             delete projY;
0799             continue;
0800         }
0801
0802         double yPeak = 0;
0803         int maxBin = projY->GetMaximumBin();
0804         yPeak = projY->GetBinCenter(maxBin);
0805
0806         double fitMin = yPeak - 0.004;
0807         double fitMax = yPeak + 0.004;
0808         projectionFit->SetRange(fitMin, fitMax);
0809         projY->Fit(projectionFit, "RQ");
0810
0811         // yPeak update!
0812         yPeak = projectionFit->GetParameter(1);
0813
0814         double xCenter = h2_XY->GetXaxis()->GetBinCenter(ix);
0815         vx.push_back(xCenter);
0816         vy.push_back(yPeak);
0817
0818         delete projY;
0819     }
0820
0821     // 构造 TGraph 并拟合
0822     int n = vx.size();
0823     TGraph* gr = new TGraph(n, vx.data(), vy.data());
0824     gr->SetName(graphName);
0825     gr->SetTitle("Peak position vs X;X;Y_{peak}");
0826     gr->Fit(peakVsXFit, "RQ");
0827
0828     return gr;
0829 }