NewCode2024/lleeCheck/GetCorrelation.cpp

0001 #include "/sphenix/user/ChengWei/INTT/INTT_dNdeta_repo/NewCode2024/Constants.cpp"
0002
0003 class GetCorrelation
0004 {
0005     public:
0006         GetCorrelation(
0007             string input_directory_in,
0008             string input_filename_in,
0009             string output_directory_in,
0010             string output_filename_in,
0011             int Nfiles_in,
0012             string MBD_z_vtx_directory_in,
0013
0014             bool DoSelfCentrality_in,
0015             bool DoVtxZReWeight_in,
0016             bool ReadRecoEvtFile_in = false
0017         );
0018
0019
0020     private:
0021         string input_directory;
0022         string input_filename;
0023         string output_directory;
0024         string output_filename;
0025         int Nfiles;
0026         string MBD_z_vtx_directory;
0027         bool DoSelfCentrality;
0028         bool DoVtxZReWeight;
0029         bool ReadRecoEvtFile;
0030
0031         // Division : ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
0032         TH1D * h1D_reweight;
0033         void GetVtxZReWeight();
0034
0035         // Division : ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
0036         TFile * file_in;
0037         TChain * tree_in;
0038         float MBD_z_vtx;
0039         bool is_min_bias;
0040         float MBD_centrality;
0041         float MBD_south_charge_sum;
0042         float MBD_north_charge_sum;
0043         float MBD_charge_sum;
0044         float MBD_charge_asymm;
0045         int InttBcoFullDiff_next;
0046         int NClus;
0047         std::vector<int> *ClusLayer;
0048         std::vector<int> *ClusAdc;
0049         std::vector<float> *ClusPhiSize;
0050         void PrepareTree();
0051
0052         // Division : ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
0053         int self_MBD_centrality;
0054         std::vector<double> MBD_charge_sum_vec_percentage_cut;
0055         void GetSelfCentralitySet();
0056         int GetSelfCentrality(double MBD_charge_sum_in);
0057
0058         // Division : ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
0059         TFile * file_out;
0060         TH2D * h2D_Mbin_NClus;
0061         TH2D * h2D_Mbin_NClus_Self; // note : self-centrality
0062         TH1D * h1D_mbd_z_vtx_narrow;
0063         TH2D * h2D_Centrality_correlation;
0064         void PrepareOut();
0065
0066         // Division : ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
0067         void mainAN();
0068
0069         // Division : ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
0070         void EndRun();
0071 };
0072
0073 GetCorrelation::GetCorrelation(
0074     string input_directory_in,
0075     string input_filename_in,
0076     string output_directory_in,
0077     string output_filename_in,
0078     int Nfiles_in,
0079     string MBD_z_vtx_directory_in,
0080
0081     bool DoSelfCentrality_in,
0082     bool DoVtxZReWeight_in,
0083     bool ReadRecoEvtFile_in
0084 ):
0085     input_directory(input_directory_in),
0086     input_filename(input_filename_in),
0087     output_directory(output_directory_in),
0088     output_filename(output_filename_in),
0089     Nfiles(Nfiles_in),
0090     MBD_z_vtx_directory(MBD_z_vtx_directory_in),
0091
0092     DoSelfCentrality(DoSelfCentrality_in),
0093     DoVtxZReWeight(DoVtxZReWeight_in),
0094     ReadRecoEvtFile(ReadRecoEvtFile_in)
0095 {
0096     if (DoVtxZReWeight){GetVtxZReWeight();}
0097     PrepareTree();
0098     PrepareOut();
0099     if (DoSelfCentrality){GetSelfCentralitySet();}
0100
0101     mainAN();
0102     EndRun();
0103 }
0104
0105 void GetCorrelation::GetVtxZReWeight()
0106 {
0107     std::string Run54280_file_directory = "/sphenix/tg/tg01/commissioning/INTT/work/cwshih/seflgendata/run_54280_HR_Feb102025/lltoee_check/MBD_z_vtx.root";
0108     std::string hist_name = "h1D_mbd_z_vtx";
0109
0110     TFile * file_in = TFile::Open(Run54280_file_directory.c_str());
0111     TH1D * hist = (TH1D*)file_in->Get(hist_name.c_str());
0112     hist->Scale(1./hist->Integral());
0113     hist->Sumw2(true);
0114
0115     TFile * file_in_2 = TFile::Open(MBD_z_vtx_directory.c_str());
0116     if (!file_in_2) {std::cout<<"file_in_2 not found"<<std::endl; exit(1);}
0117
0118     TH1D * hist_2 = (TH1D*)file_in_2->Get(hist_name.c_str());
0119     hist_2->Scale(1./hist_2->Integral());
0120
0121     h1D_reweight = (TH1D*)hist->Clone("h1D_reweight");
0122     h1D_reweight->Reset("ICESM");
0123     h1D_reweight->Divide(hist, hist_2);
0124 }
0125
0126 void GetCorrelation::PrepareTree()
0127 {
0128     tree_in = new TChain("EventTree");
0129
0130     for (int i = 0; i < Nfiles; i++)
0131     {
0132         std::string job_index = std::to_string( i );
0133         int job_index_len = 5;
0134         job_index.insert(0, job_index_len - job_index.size(), '0');
0135
0136         std::string input_filename_full = Form("%s/%s%s.root", input_directory.c_str(), input_filename.c_str(), job_index.c_str());
0137         std::cout<<"input_filename_full: "<<input_filename_full<<std::endl;
0138         tree_in -> Add(input_filename_full.c_str());
0139     }
0140
0141     std::cout<<"tree_in -> GetEntries(): "<<tree_in -> GetEntries()<<std::endl;
0142
0143     tree_in -> SetBranchStatus("*", 0);
0144
0145     tree_in -> SetBranchStatus("MBD_z_vtx", 1);
0146     tree_in -> SetBranchStatus("is_min_bias", 1);
0147     tree_in -> SetBranchStatus("MBD_centrality", 1);
0148     tree_in -> SetBranchStatus("MBD_south_charge_sum", 1);
0149     tree_in -> SetBranchStatus("MBD_north_charge_sum", 1);
0150     tree_in -> SetBranchStatus("MBD_charge_sum", 1);
0151     tree_in -> SetBranchStatus("MBD_charge_asymm", 1);
0152
0153     tree_in -> SetBranchStatus("NClus", 1);
0154     tree_in -> SetBranchStatus("ClusAdc", 1);
0155     tree_in -> SetBranchStatus("ClusPhiSize", 1);
0156
0157     if (ReadRecoEvtFile){tree_in -> SetBranchStatus("InttBcoFullDiff_next", 1);}
0158
0159
0160     ClusLayer = 0;
0161     ClusAdc = 0;
0162     ClusPhiSize = 0;
0163
0164     tree_in -> SetBranchAddress("MBD_z_vtx", &MBD_z_vtx);
0165     tree_in -> SetBranchAddress("is_min_bias", &is_min_bias);
0166     tree_in -> SetBranchAddress("MBD_centrality", &MBD_centrality);
0167     tree_in -> SetBranchAddress("MBD_south_charge_sum", &MBD_south_charge_sum);
0168     tree_in -> SetBranchAddress("MBD_north_charge_sum", &MBD_north_charge_sum);
0169     tree_in -> SetBranchAddress("MBD_charge_sum", &MBD_charge_sum);
0170     tree_in -> SetBranchAddress("MBD_charge_asymm", &MBD_charge_asymm);
0171
0172     tree_in -> SetBranchAddress("NClus",&NClus);
0173     tree_in -> SetBranchAddress("ClusAdc",&ClusAdc);
0174     tree_in -> SetBranchAddress("ClusPhiSize",&ClusPhiSize);
0175
0176     if (ReadRecoEvtFile){tree_in -> SetBranchAddress("InttBcoFullDiff_next", &InttBcoFullDiff_next);}
0177 }
0178
0179 void GetCorrelation::PrepareOut()
0180 {
0181     vector<double> centrality_edges = Constants::centrality_edges;
0182
0183     file_out = new TFile(Form("%s/%s", output_directory.c_str(), output_filename.c_str()), "recreate");
0184     h2D_Mbin_NClus = new TH2D("h2D_Mbin_NClus","h2D_Mbin_NClus;Centrality;NClus",centrality_edges.size() - 1, &centrality_edges[0], 500, 0, 10000);
0185     h2D_Mbin_NClus_Self = new TH2D("h2D_Mbin_NClus_Self","h2D_Mbin_NClus_Self;Centrality;NClus",centrality_edges.size() - 1, &centrality_edges[0], 500, 0, 10000);
0186     h1D_mbd_z_vtx_narrow = new TH1D("h1D_mbd_z_vtx_narrow", "h1D_mbd_z_vtx_narrow;MBD_z_vtx [cm];Counts", 120, -60, 60);
0187     h2D_Centrality_correlation = new TH2D("h2D_Centrality_correlation","h2D_Centrality_correlation;Official Centrality;Self Centrality", 101, -0.5, 100.5, 101, -0.5, 100.5);
0188 }
0189
0190 void GetCorrelation::GetSelfCentralitySet()
0191 {
0192
0193     std::vector<int> MBD_charge_sum_vec; MBD_charge_sum_vec.clear();
0194
0195     for (int i = 0; i < tree_in -> GetEntries(); i++)
0196     {
0197         tree_in -> GetEntry(i);
0198
0199         if (i % 10000 == 0)
0200         {
0201             std::cout<<"self-centrality, i: "<<i<<std::endl;
0202         }
0203
0204         // note : both data and MC
0205         if (is_min_bias != 1) {continue;}
0206         if (MBD_z_vtx != MBD_z_vtx) {continue;}
0207         if (MBD_centrality != MBD_centrality) {continue;}
0208         if (MBD_centrality < 0 || MBD_centrality > 100) {continue;}
0209         if (NClus <= 5) {continue;}
0210         if (MBD_z_vtx < -10 || MBD_z_vtx >= 10) {continue;}
0211
0212         if (ReadRecoEvtFile && InttBcoFullDiff_next <= 61) {continue;}
0213
0214         if (MBD_centrality > 92) {continue;}
0215
0216         MBD_charge_sum_vec.push_back(MBD_charge_sum);
0217     }
0218
0219     std::sort(MBD_charge_sum_vec.begin(), MBD_charge_sum_vec.end(),std::greater<int>());
0220
0221     MBD_charge_sum_vec_percentage_cut.clear();
0222
0223
0224     for (int i = 1; i <= 92; i++)
0225     {
0226         int index = (int)(MBD_charge_sum_vec.size() * i / 100);
0227         MBD_charge_sum_vec_percentage_cut.push_back(MBD_charge_sum_vec[index]);
0228
0229         std::cout<<"i: "<<i<<", MBD_charge_sum_vec[index]: "<<MBD_charge_sum_vec[index]<<std::endl;
0230     }
0231 }
0232
0233 int GetCorrelation::GetSelfCentrality(double MBD_charge_sum_in)
0234 {
0235     for (int i = 0; i < MBD_charge_sum_vec_percentage_cut.size(); i++)
0236     {
0237         if (MBD_charge_sum_in >= MBD_charge_sum_vec_percentage_cut[i])
0238         {
0239             return i + 1;
0240         }
0241     }
0242
0243     return 0;
0244 }
0245
0246 void GetCorrelation::mainAN()
0247 {
0248     for (int i = 0; i < tree_in -> GetEntries(); i++)
0249     {
0250         tree_in -> GetEntry(i);
0251
0252         if (i % 10000 == 0)
0253         {
0254             std::cout<<"i: "<<i<<std::endl;
0255         }
0256
0257         int nHadron_count = 0;
0258
0259         // note : both data and MC
0260         if (is_min_bias != 1) {continue;}
0261         if (MBD_z_vtx != MBD_z_vtx) {continue;}
0262         if (MBD_centrality != MBD_centrality) {continue;}
0263         if (MBD_centrality < 0 || MBD_centrality > 100) {continue;}
0264         if (NClus <= 5) {continue;}
0265         if (MBD_z_vtx < -10 || MBD_z_vtx >= 10) {continue;}
0266
0267         if (ReadRecoEvtFile && InttBcoFullDiff_next <= 61) {continue;}
0268
0269         if (MBD_centrality > 92) {continue;}
0270
0271         double INTTvtxZWeighting;
0272         INTTvtxZWeighting = (DoVtxZReWeight) ? h1D_reweight -> GetBinContent(h1D_reweight -> FindBin(MBD_z_vtx)) : 1.;
0273
0274
0275         int NClus_count = 0;
0276         for (int clu_i = 0; clu_i < NClus; clu_i++)
0277         {
0278             if (ClusAdc -> at(clu_i) <= 35) {continue;}
0279             if (ClusPhiSize -> at(clu_i) > 40) {continue;}
0280             NClus_count++;
0281         }
0282
0283         h2D_Mbin_NClus -> Fill(MBD_centrality, NClus_count, INTTvtxZWeighting);
0284         h1D_mbd_z_vtx_narrow -> Fill(MBD_z_vtx, INTTvtxZWeighting);
0285
0286         if (DoSelfCentrality)
0287         {
0288             self_MBD_centrality = GetSelfCentrality(MBD_charge_sum);
0289             h2D_Mbin_NClus_Self -> Fill(self_MBD_centrality, NClus_count, INTTvtxZWeighting);
0290             h2D_Centrality_correlation -> Fill(MBD_centrality, self_MBD_centrality, INTTvtxZWeighting);
0291         }
0292     }
0293 }
0294
0295 void GetCorrelation::EndRun()
0296 {
0297     file_out -> cd();
0298     h2D_Mbin_NClus -> Write();
0299     h1D_mbd_z_vtx_narrow -> Write();
0300     if (DoVtxZReWeight) {h1D_reweight -> Write("h1D_reweight");}
0301     if (DoSelfCentrality){h2D_Centrality_correlation -> Write();}
0302     if (DoSelfCentrality){h2D_Mbin_NClus_Self -> Write();}
0303
0304     file_out -> Close();
0305 }
0306
0307 void GetCorrelation_main()
0308 {
0309     string input_directory = "/sphenix/tg/tg01/commissioning/INTT/work/cwshih/seflgendata/run_54280_HR_Feb102025/Run6_EvtZFitWidthChange/EvtVtxZ/completed";
0310     string input_filename = "Data_EvtVtxZProtoTracklet_FieldOff_BcoFullDiff_VtxZReco_00054280_";
0311
0312     string output_directory = "/sphenix/tg/tg01/commissioning/INTT/work/cwshih/seflgendata/run_54280_HR_Feb102025/lltoee_check";
0313     string output_filename = Form("Data_MBin92_NClus.root");
0314
0315     int Nfiles = 1000;
0316
0317     string MBD_z_vtx_directory = "no_map";
0318
0319     bool DoSelfCentrality = true;
0320     bool DoVtxZReWeight = false;
0321     bool ReadRecoEvtFile = true;
0322
0323     GetCorrelation * GetCorrelation_ = new GetCorrelation(
0324         input_directory,
0325         input_filename,
0326         output_directory,
0327         output_filename,
0328         Nfiles,
0329         MBD_z_vtx_directory,
0330
0331         DoSelfCentrality,
0332         DoVtxZReWeight,
0333         ReadRecoEvtFile
0334     );
0335 }
0336
0337
0338 void GetCorrelation_MCRun54280()
0339 {
0340     string input_directory = "/sphenix/tg/tg01/commissioning/INTT/work/cwshih/sPH_dNdeta/Run24AuAuMC/Sim_HIJING_MDC2_ana472_20250307/per5k";
0341     string input_filename = "ntuple_per5k_";
0342
0343     string output_directory = "/sphenix/tg/tg01/commissioning/INTT/work/cwshih/sPH_dNdeta/Run24AuAuMC/Sim_HIJING_MDC2_ana472_20250307/lltoee_check";
0344     string output_filename = Form("MC_MBin92_NClus.root");
0345
0346     int Nfiles = 146;
0347     string MBD_z_vtx_directory = "/sphenix/tg/tg01/commissioning/INTT/work/cwshih/sPH_dNdeta/Run24AuAuMC/Sim_HIJING_MDC2_ana472_20250307/lltoee_check/MBD_z_vtx.root";
0348
0349     bool DoSelfCentrality = true;
0350     bool DoVtxZReWeight = true;
0351     bool ReadRecoEvtFile = false;
0352
0353     GetCorrelation * GetCorrelation_ = new GetCorrelation(
0354         input_directory,
0355         input_filename,
0356         output_directory,
0357         output_filename,
0358         Nfiles,
0359         MBD_z_vtx_directory,
0360
0361         DoSelfCentrality,
0362         DoVtxZReWeight,
0363         ReadRecoEvtFile
0364     );
0365 }
0366
0367
0368 void GetCorrelation_MCFieldON()
0369 {
0370     string input_directory = "/sphenix/tg/tg01/commissioning/INTT/work/cwshih/sPH_dNdeta/Run24AuAuMC/WrongStrangeIncrease_FieldOn_Sim_HIJING_strangeness_MDC2_ana467_20250226/per5k";
0371     string input_filename = "ntuple_per5k_";
0372
0373     string output_directory = "/sphenix/tg/tg01/commissioning/INTT/work/cwshih/sPH_dNdeta/Run24AuAuMC/WrongStrangeIncrease_FieldOn_Sim_HIJING_strangeness_MDC2_ana467_20250226/lltoee_check";
0374     string output_filename = Form("MC_MBin92_NClus.root");
0375
0376     int Nfiles = 196;
0377     string MBD_z_vtx_directory = "/sphenix/tg/tg01/commissioning/INTT/work/cwshih/sPH_dNdeta/Run24AuAuMC/WrongStrangeIncrease_FieldOn_Sim_HIJING_strangeness_MDC2_ana467_20250226/lltoee_check/MBD_z_vtx.root";
0378
0379     bool DoSelfCentrality = true;
0380     bool DoVtxZReWeight = true;
0381     bool ReadRecoEvtFile = false;
0382
0383     GetCorrelation * GetCorrelation_ = new GetCorrelation(
0384         input_directory,
0385         input_filename,
0386         output_directory,
0387         output_filename,
0388         Nfiles,
0389         MBD_z_vtx_directory,
0390
0391         DoSelfCentrality,
0392         DoVtxZReWeight,
0393         ReadRecoEvtFile
0394     );
0395 }
0396
0397
0398 void GetCorrelation_DataRun54912()
0399 {
0400     string input_directory = "/sphenix/tg/tg01/commissioning/INTT/work/cwshih/seflgendata/run_54912/completed";
0401     string input_filename = "testNtuple_";
0402
0403     string output_directory = "/sphenix/tg/tg01/commissioning/INTT/work/cwshih/seflgendata/run_54912/lltoee_check";
0404     string output_filename = Form("Data_MBin92_NClus.root");
0405
0406     int Nfiles = 1;
0407     string MBD_z_vtx_directory = "/sphenix/tg/tg01/commissioning/INTT/work/cwshih/seflgendata/run_54912/lltoee_check/MBD_z_vtx.root";
0408
0409     bool DoSelfCentrality = true;
0410     bool DoVtxZReWeight = true;
0411     bool ReadRecoEvtFile = false;
0412
0413     GetCorrelation * GetCorrelation_ = new GetCorrelation(
0414         input_directory,
0415         input_filename,
0416         output_directory,
0417         output_filename,
0418         Nfiles,
0419         MBD_z_vtx_directory,
0420
0421         DoSelfCentrality,
0422         DoVtxZReWeight,
0423         ReadRecoEvtFile
0424     );
0425 }