NewCode2024/lleeCheck/quick_check_llee_chain.C

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