NewCode2024/lleeCheck/make_MBD_dist.C

0001 #include "/sphenix/user/ChengWei/INTT/INTT_dNdeta_repo/NewCode2024/Constants.cpp"
0002
0003 int make_MBD_dist(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("MBD_z_vtx.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     TH1D * h1D_mbd_z_vtx = new TH1D("h1D_mbd_z_vtx", "h1D_mbd_z_vtx;MBD_z_vtx [cm];Counts", 120, -60, 60);
0188
0189     TFile * file_out = new TFile(Form("%s/%s", output_directory.c_str(), output_filename.c_str()), "recreate");
0190
0191     std::cout<<222<<std::endl;
0192
0193     for (int i = 0; i < tree_in -> GetEntries(); i++)
0194     {
0195         // std::cout<<333<<std::endl;
0196         tree_in -> GetEntry(i);
0197         // std::cout<<444<<std::endl;
0198
0199         if (i % 1000 == 0)
0200         {
0201             std::cout<<"i: "<<i<<std::endl;
0202         }
0203
0204         int nHadron_count = 0;
0205
0206         // note : for MC
0207         // if (NTruthVtx != 1) {continue;}
0208
0209         // if (is_min_bias != 1) {continue;}
0210         // if (INTTvtxZ != INTTvtxZ) {continue;}
0211         // if (INTTvtxZ < -10 || INTTvtxZ >= 10) {continue;}
0212         // if (MBD_centrality < 0) {continue;}
0213         // if (MBD_centrality > 3) {continue;} // note : [0-3)% centrality interval
0214
0215
0216          // note : both data and MC
0217         if (is_min_bias != 1) {continue;}
0218         if (MBD_z_vtx != MBD_z_vtx) {continue;}
0219         if (MBD_centrality != MBD_centrality) {continue;}
0220         if (MBD_centrality < 0 || MBD_centrality > 100) {continue;}
0221         if (isData && InttBcoFullDiff_next <= 61) {continue;}
0222
0223         if (MBD_centrality > 70) {continue;}
0224
0225         h1D_mbd_z_vtx -> Fill(MBD_z_vtx);
0226     }
0227
0228     file_out -> cd();
0229     h1D_mbd_z_vtx -> Write();
0230
0231     file_out -> Close();
0232
0233     // system(Form("mv %s/%s %s/completed/%s", output_directory.c_str(), output_filename.c_str(), output_directory.c_str(), output_filename.c_str()));
0234
0235     return 0;
0236
0237 }
0238
0239
0240 int make_MBD_dist_2(bool isData = true) // note : for field ON
0241 {
0242     string input_directory;
0243     string input_filename;
0244     string output_directory;
0245     string output_filename;
0246     int Nfiles;
0247
0248     if (isData == true){
0249         input_directory = "/sphenix/tg/tg01/commissioning/INTT/work/cwshih/seflgendata/run_54912/completed";
0250         input_filename = "testNtuple_"; // note : testNtuple_00000.root
0251
0252         output_directory = "/sphenix/tg/tg01/commissioning/INTT/work/cwshih/seflgendata/run_54912/lltoee_check";
0253         output_filename = Form("MBD_z_vtx.root");
0254
0255         Nfiles = 1;
0256     }
0257     else {
0258         input_directory = "/sphenix/user/ChengWei/sPH_dNdeta/Run24AuAuMC/WrongStrangeIncrease_FieldOn_Sim_HIJING_strangeness_MDC2_ana467_20250226/per5k";
0259         input_filename = "ntuple_per5k_";
0260
0261         output_directory = input_directory;
0262         output_filename = Form("MBD_z_vtx.root");
0263
0264         Nfiles = 196;
0265     }
0266
0267
0268     // TFile * file_in = TFile::Open(Form("%s/%s", input_directory.c_str(), input_filename.c_str()));
0269     // TTree * tree_in = (TTree*)file_in->Get("EventTree");
0270
0271     TChain * tree_in = new TChain("EventTree");
0272     // tree_in -> Add(Form("%s/%s000{00..24}.root", input_directory.c_str(), input_filename.c_str()));
0273
0274     for (int i = 0; i < Nfiles; i++)
0275     {
0276         std::string job_index = std::to_string( i );
0277         int job_index_len = 5;
0278         job_index.insert(0, job_index_len - job_index.size(), '0');
0279
0280         std::string input_filename_full = Form("%s/%s%s.root", input_directory.c_str(), input_filename.c_str(), job_index.c_str());
0281         std::cout<<"input_filename_full: "<<input_filename_full<<std::endl;
0282         tree_in -> Add(input_filename_full.c_str());
0283     }
0284
0285     // std::cout<<"aaaa"<<std::endl;
0286     // tree_in->GetEntry(0);
0287     // std::cout<<"aaaa"<<std::endl;
0288
0289     std::cout<<"tree_in -> GetEntries(): "<<tree_in -> GetEntries()<<std::endl;
0290
0291     tree_in -> SetBranchStatus("*", 0);
0292
0293     tree_in -> SetBranchStatus("MBD_z_vtx", 1);
0294     tree_in -> SetBranchStatus("is_min_bias", 1);
0295     tree_in -> SetBranchStatus("MBD_centrality", 1);
0296     tree_in -> SetBranchStatus("MBD_south_charge_sum", 1);
0297     tree_in -> SetBranchStatus("MBD_north_charge_sum", 1);
0298     tree_in -> SetBranchStatus("MBD_charge_sum", 1);
0299     tree_in -> SetBranchStatus("MBD_charge_asymm", 1);
0300
0301     float MBD_z_vtx;
0302     bool is_min_bias;
0303     float MBD_centrality;
0304     float MBD_south_charge_sum;
0305     float MBD_north_charge_sum;
0306     float MBD_charge_sum;
0307     float MBD_charge_asymm;
0308     int InttBcoFullDiff_next;
0309
0310
0311     tree_in -> SetBranchAddress("MBD_z_vtx", &MBD_z_vtx);
0312     tree_in -> SetBranchAddress("is_min_bias", &is_min_bias);
0313     tree_in -> SetBranchAddress("MBD_centrality", &MBD_centrality);
0314     tree_in -> SetBranchAddress("MBD_south_charge_sum", &MBD_south_charge_sum);
0315     tree_in -> SetBranchAddress("MBD_north_charge_sum", &MBD_north_charge_sum);
0316     tree_in -> SetBranchAddress("MBD_charge_sum", &MBD_charge_sum);
0317     tree_in -> SetBranchAddress("MBD_charge_asymm", &MBD_charge_asymm);
0318
0319
0320
0321
0322     // Division : ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
0323
0324     vector<double> centrality_edges = Constants::centrality_edges;
0325     // TH2D * h2D_Mbin_NClus = new TH2D("h2D_Mbin_NClus","h2D_Mbin_NClus;Centrality;NClus",centrality_edges.size() - 1, &centrality_edges[0], 500, 0, 10000);
0326
0327     TH1D * h1D_mbd_z_vtx = new TH1D("h1D_mbd_z_vtx", "h1D_mbd_z_vtx;MBD_z_vtx [cm];Counts", 120, -60, 60);
0328
0329     TFile * file_out = new TFile(Form("%s/%s", output_directory.c_str(), output_filename.c_str()), "recreate");
0330
0331     std::cout<<222<<std::endl;
0332
0333     for (int i = 0; i < tree_in -> GetEntries(); i++)
0334     {
0335         // std::cout<<333<<std::endl;
0336         tree_in -> GetEntry(i);
0337         // std::cout<<444<<std::endl;
0338
0339         if (i % 1000 == 0)
0340         {
0341             std::cout<<"i: "<<i<<std::endl;
0342         }
0343
0344         // int nHadron_count = 0;
0345
0346         // note : for MC
0347         // if (NTruthVtx != 1) {continue;}
0348
0349         // if (is_min_bias != 1) {continue;}
0350         // if (INTTvtxZ != INTTvtxZ) {continue;}
0351         // if (INTTvtxZ < -10 || INTTvtxZ >= 10) {continue;}
0352         // if (MBD_centrality < 0) {continue;}
0353         // if (MBD_centrality > 3) {continue;} // note : [0-3)% centrality interval
0354
0355
0356          // note : both data and MC
0357         if (is_min_bias != 1) {continue;}
0358         if (MBD_z_vtx != MBD_z_vtx) {continue;}
0359         if (MBD_centrality != MBD_centrality) {continue;}
0360         if (MBD_centrality < 0 || MBD_centrality > 100) {continue;}
0361
0362         if (MBD_centrality > 70) {continue;}
0363
0364         h1D_mbd_z_vtx -> Fill(MBD_z_vtx);
0365     }
0366
0367     file_out -> cd();
0368     h1D_mbd_z_vtx -> Write();
0369
0370     file_out -> Close();
0371
0372     // system(Form("mv %s/%s %s/completed/%s", output_directory.c_str(), output_filename.c_str(), output_directory.c_str(), output_filename.c_str()));
0373
0374     return 0;
0375
0376 }