PrepareAllDist/vtxZDist/vtxZDist.cpp

0001 #include "vtxZDist.h"
0002
0003 vtxZDist::vtxZDist(
0004     int process_id_in,
0005     int runnumber_in,
0006     int run_nEvents_in,
0007     std::string input_directory_in,
0008     std::string input_file_name_in,
0009     std::string output_directory_in,
0010
0011     std::string output_file_name_suffix_in,
0012
0013     bool Apply_cut_in,
0014     bool ApplyVtxZReWeighting_in,
0015     std::pair<bool, int> ApplyEvtBcoFullDiffCut_in
0016 ):
0017     process_id(process_id_in),
0018     runnumber(runnumber_in),
0019     run_nEvents(run_nEvents_in),
0020     input_directory(input_directory_in),
0021     input_file_name(input_file_name_in),
0022     output_directory(output_directory_in),
0023     output_file_name_suffix(output_file_name_suffix_in),
0024     Apply_cut(Apply_cut_in),
0025     ApplyVtxZReWeighting(ApplyVtxZReWeighting_in),
0026     ApplyEvtBcoFullDiffCut(ApplyEvtBcoFullDiffCut_in)
0027 {
0028     PrepareOutPutFileName();
0029     PrepareInputFile();
0030
0031     run_nEvents = (run_nEvents == -1) ? tree_in->GetEntries() : run_nEvents;
0032     run_nEvents = (run_nEvents > tree_in->GetEntries()) ? tree_in->GetEntries() : run_nEvents;
0033     nCentrality_bin = centrality_edges.size() - 1;
0034
0035     PrepareOutputFile();
0036     PrepareHist();
0037
0038 }
0039
0040 void vtxZDist::PrepareOutPutFileName()
0041 {
0042     std::string job_index = std::to_string( process_id );
0043     int job_index_len = 5;
0044     job_index.insert(0, job_index_len - job_index.size(), '0');
0045
0046     std::string runnumber_str = std::to_string( runnumber );
0047     if (runnumber != -1){
0048         int runnumber_str_len = 8;
0049         runnumber_str.insert(0, runnumber_str_len - runnumber_str.size(), '0');
0050     }
0051
0052     if (output_file_name_suffix.size() > 0 && output_file_name_suffix[0] != '_') {
0053         output_file_name_suffix = "_" + output_file_name_suffix;
0054     }
0055
0056     output_filename = "vtxZDist";
0057     output_filename = (runnumber != -1) ? "Data_" + output_filename : "MC_" + output_filename;
0058
0059     output_filename += (Apply_cut) ? "_VtxZQA" : "_NoVtxZQA";
0060     output_filename += (ApplyVtxZReWeighting) ? "_VtxZReWeighting" : "";
0061     output_filename += (runnumber != -1 && ApplyEvtBcoFullDiffCut.first) ? Form("_EvtBcoFullDiffCut%d", ApplyEvtBcoFullDiffCut.second) : "";
0062
0063     output_filename += output_file_name_suffix;
0064     output_filename += (runnumber != -1) ? Form("_%s_%s.root",runnumber_str.c_str(),job_index.c_str()) : Form("_%s.root",job_index.c_str());
0065 }
0066
0067 std::map<std::string, int> vtxZDist::GetInputTreeBranches(TTree * m_tree_in)
0068 {
0069     std::map<std::string, int> branch_map;
0070     TObjArray * branch_list = m_tree_in -> GetListOfBranches();
0071     for (int i = 0; i < branch_list -> GetEntries(); i++)
0072     {
0073         TBranch * branch = dynamic_cast<TBranch*>(branch_list->At(i));
0074         branch_map[branch -> GetName()] = 1;
0075     }
0076     return branch_map;
0077 }
0078
0079 void vtxZDist::PrepareInputFile()
0080 {
0081     file_in = TFile::Open(Form("%s/%s", input_directory.c_str(), input_file_name.c_str()));
0082     tree_in = (TTree*)file_in->Get(tree_name.c_str());
0083
0084     std::map<std::string, int> branch_map = GetInputTreeBranches(tree_in);
0085
0086     tree_in -> SetBranchStatus("*", 0);
0087
0088     tree_in -> SetBranchStatus("is_min_bias", 1);
0089     tree_in -> SetBranchStatus("MBD_centrality", 1);
0090     tree_in -> SetBranchStatus("MBD_z_vtx", 1);
0091     tree_in -> SetBranchStatus("MBD_charge_sum", 1);
0092
0093     tree_in -> SetBranchStatus("INTTvtxZ", 1);
0094     tree_in -> SetBranchStatus("INTTvtxZError", 1);
0095     tree_in -> SetBranchStatus("NgroupTrapezoidal", 1);
0096     tree_in -> SetBranchStatus("NgroupCoarse", 1);
0097     tree_in -> SetBranchStatus("TrapezoidalFitWidth", 1);
0098     tree_in -> SetBranchStatus("TrapezoidalFWHM", 1);
0099
0100     tree_in -> SetBranchStatus("NClus", 1);
0101     tree_in -> SetBranchStatus("NClus_Layer1", 1);
0102
0103     // note : for data
0104     if (branch_map.find("MBDNSg2") != branch_map.end()) {
0105         tree_in -> SetBranchStatus("MBDNSg2", 1);
0106         m_withTrig = true;
0107     }
0108     if (branch_map.find("MBDNSg2_vtxZ10cm") != branch_map.end()) {tree_in -> SetBranchStatus("MBDNSg2_vtxZ10cm", 1);}
0109     if (branch_map.find("MBDNSg2_vtxZ30cm") != branch_map.end()) {tree_in -> SetBranchStatus("MBDNSg2_vtxZ30cm", 1);}
0110
0111     if (branch_map.find("InttBcoFullDiff_next") != branch_map.end()) {tree_in -> SetBranchStatus("InttBcoFullDiff_next", 1); }
0112
0113     // note :for MC
0114     if (branch_map.find("NTruthVtx") != branch_map.end()) {tree_in -> SetBranchStatus("NTruthVtx", 1);}
0115     if (branch_map.find("TruthPV_trig_z") != branch_map.end()) {tree_in -> SetBranchStatus("TruthPV_trig_z", 1);}
0116
0117     // Division : ---SetBranchAddress-----------------------------------------------------------------------------------------------
0118     tree_in -> SetBranchAddress("is_min_bias", &is_min_bias);
0119     tree_in -> SetBranchAddress("MBD_centrality", &MBD_centrality);
0120     tree_in -> SetBranchAddress("MBD_z_vtx", &MBD_z_vtx);
0121     tree_in -> SetBranchAddress("MBD_charge_sum", &MBD_charge_sum);
0122
0123     tree_in -> SetBranchAddress("INTTvtxZ", &INTTvtxZ);
0124     tree_in -> SetBranchAddress("INTTvtxZError", &INTTvtxZError);
0125     tree_in -> SetBranchAddress("NgroupTrapezoidal", &NgroupTrapezoidal);
0126     tree_in -> SetBranchAddress("NgroupCoarse", &NgroupCoarse);
0127     tree_in -> SetBranchAddress("TrapezoidalFitWidth", &TrapezoidalFitWidth);
0128     tree_in -> SetBranchAddress("TrapezoidalFWHM", &TrapezoidalFWHM);
0129
0130     tree_in -> SetBranchAddress("NClus", &NClus);
0131     tree_in -> SetBranchAddress("NClus_Layer1", &NClus_Layer1);
0132
0133     // note : for data
0134     if (branch_map.find("MBDNSg2") != branch_map.end()) {tree_in -> SetBranchAddress("MBDNSg2", &MBDNSg2);}
0135     if (branch_map.find("MBDNSg2_vtxZ10cm") != branch_map.end()) {tree_in -> SetBranchAddress("MBDNSg2_vtxZ10cm", &MBDNSg2_vtxZ10cm);}
0136     if (branch_map.find("MBDNSg2_vtxZ30cm") != branch_map.end()) {tree_in -> SetBranchAddress("MBDNSg2_vtxZ30cm", &MBDNSg2_vtxZ30cm);}
0137
0138     if (branch_map.find("InttBcoFullDiff_next") != branch_map.end()) {tree_in -> SetBranchAddress("InttBcoFullDiff_next", &InttBcoFullDiff_next); }
0139
0140     // note : for MC
0141     if (branch_map.find("NTruthVtx") != branch_map.end()) {tree_in -> SetBranchAddress("NTruthVtx", &NTruthVtx);}
0142     if (branch_map.find("TruthPV_trig_z") != branch_map.end()) {tree_in -> SetBranchAddress("TruthPV_trig_z", &TruthPV_trig_z);}
0143 }
0144
0145 void vtxZDist::PrepareOutputFile()
0146 {
0147     file_out = new TFile(Form("%s/%s", output_directory.c_str(), output_filename.c_str()), "RECREATE");
0148 }
0149
0150 void vtxZDist::PrepareHist()
0151 {
0152     h1D_centrality_bin = new TH1D("h1D_centrality_bin","h1D_centrality_bin;Centrality [%];Entries",nCentrality_bin,&centrality_edges[0]); // note : the 0-5%
0153
0154     h1D_map.clear();
0155
0156     h1D_map.insert(
0157         std::make_pair(
0158             "h1D_INTTz_Inclusive100",
0159             new TH1D("h1D_INTTz_Inclusive100","h1D_INTTz_Inclusive100;INTT vtxZ [cm];Entries", nVtxZ_bin, vtxZ_range.first, vtxZ_range.second)
0160         )
0161     );
0162
0163     h1D_map.insert(
0164         std::make_pair(
0165             "h1D_INTTz_Inclusive70",
0166             new TH1D("h1D_INTTz_Inclusive70","h1D_INTTz_Inclusive70;INTT vtxZ [cm];Entries", nVtxZ_bin, vtxZ_range.first, vtxZ_range.second)
0167         )
0168     );
0169
0170     h1D_map.insert(
0171         std::make_pair(
0172             "h1D_INTTz_Inclusive100_tight",
0173             new TH1D("h1D_INTTz_Inclusive100_tight","h1D_INTTz_Inclusive100_tight;INTT vtxZ [cm];Entries", nVtxZ_bin, vtxZ_range.first, vtxZ_range.second)
0174         )
0175     );
0176
0177     h1D_map.insert(
0178         std::make_pair(
0179             "h1D_INTTz_Inclusive70_tight",
0180             new TH1D("h1D_INTTz_Inclusive70_tight","h1D_INTTz_Inclusive70_tight;INTT vtxZ [cm];Entries", nVtxZ_bin, vtxZ_range.first, vtxZ_range.second)
0181         )
0182     );
0183
0184     for (int i = 0; i < nCentrality_bin; i++)
0185     {
0186         h1D_map.insert(
0187             std::make_pair(
0188                 Form("h1D_INTTz_Mbin%d", i),
0189                 new TH1D(Form("h1D_INTTz_Mbin%d", i), Form("h1D_INTTz_Mbin%d;INTT vtxZ [cm];Entries", i), nVtxZ_bin, vtxZ_range.first, vtxZ_range.second)
0190             )
0191         );
0192     }
0193
0194
0195     // h1D_map.insert(
0196     //     std::make_pair(
0197     //         "h1D_INTTzNarrow_Inclusive100",
0198     //         new TH1D("h1D_INTTzNarrow_Inclusive100","h1D_INTTzNarrow_Inclusive100;INTT vtxZ [cm];Entries", nVtxZ_bin_narrow, vtxZ_range_narrow.first, vtxZ_range_narrow.second)
0199     //     )
0200     // );
0201
0202     // h1D_map.insert(
0203     //     std::make_pair(
0204     //         "h1D_INTTzNarrow_Inclusive70",
0205     //         new TH1D("h1D_INTTzNarrow_Inclusive70","h1D_INTTzNarrow_Inclusive70;INTT vtxZ [cm];Entries", nVtxZ_bin_narrow, vtxZ_range_narrow.first, vtxZ_range_narrow.second)
0206     //     )
0207     // );
0208
0209     // for (int i = 0; i < nCentrality_bin; i++)
0210     // {
0211     //     h1D_map.insert(
0212     //         std::make_pair(
0213     //             Form("h1D_INTTzNarrow_Mbin%d", i),
0214     //             new TH1D(Form("h1D_INTTzNarrow_Mbin%d", i), Form("h1D_INTTzNarrow_Mbin%d;INTT vtxZ [cm];Entries", i), nVtxZ_bin_narrow, vtxZ_range_narrow.first, vtxZ_range_narrow.second)
0215     //         )
0216     //     );
0217     // }
0218
0219
0220     // note : MBD-INTT vtxZ Diff
0221     h1D_map.insert(
0222         std::make_pair(
0223             "h1D_VtxZDiff_MBD_INTT_Inclusive100",
0224             new TH1D("h1D_VtxZDiff_MBD_INTT_Inclusive100","h1D_VtxZDiff_MBD_INTT_Inclusive100;MBD vtxZ - INTT vtxZ [cm];Entries", 100, -15, 15)
0225         )
0226     );
0227
0228     h1D_map.insert(
0229         std::make_pair(
0230             "h1D_VtxZDiff_MBD_INTT_Inclusive70",
0231             new TH1D("h1D_VtxZDiff_MBD_INTT_Inclusive70","h1D_VtxZDiff_MBD_INTT_Inclusive70;MBD vtxZ - INTT vtxZ [cm];Entries", 100, -15, 15)
0232         )
0233     );
0234
0235     h1D_map.insert(
0236         std::make_pair(
0237             "h1D_VtxZDiff_MBD_INTT_HighNClus",
0238             new TH1D("h1D_VtxZDiff_MBD_INTT_HighNClus","h1D_VtxZDiff_MBD_INTT_HighNClus;MBD vtxZ - INTT vtxZ [cm];Entries", 100, -15, 15)
0239         )
0240     );
0241
0242
0243
0244     // note : for the INTT vtxZ profile
0245     h1D_map.insert(
0246         std::make_pair(
0247             "h1D_TrapezoidalFitWidth",
0248             new TH1D("h1D_TrapezoidalFitWidth","h1D_TrapezoidalFitWidth;INTT vtxZ Dist Fit Width [cm];Entries", 100, 0, 15)
0249         )
0250     );
0251
0252     h1D_map.insert(
0253         std::make_pair(
0254             "h1D_TrapezoidalFWHM",
0255             new TH1D("h1D_TrapezoidalFWHM","h1D_TrapezoidalFWHM;INTT vtxZ Dist FWHM [cm];Entries", 100, 0, 15)
0256         )
0257     );
0258
0259     h1D_map.insert(
0260         std::make_pair(
0261             "h1D_INTTvtxZError",
0262             new TH1D("h1D_INTTvtxZError","h1D_INTTvtxZError;INTT vtxZ StdDev [cm];Entries", 100, 0, 15)
0263         )
0264     );
0265
0266     h1D_map.insert(
0267         std::make_pair(
0268             "h1D_NgroupTrapezoidal",
0269             new TH1D("h1D_NgroupTrapezoidal","h1D_NgroupTrapezoidal;N group (Trapezoidal);Entries", 15, 0, 15)
0270         )
0271     );
0272
0273     h1D_map.insert(
0274         std::make_pair(
0275             "h1D_NgroupCoarse",
0276             new TH1D("h1D_NgroupCoarse","h1D_NgroupCoarse;N group (Coarse);Entries", 15, 0, 15)
0277         )
0278     );
0279
0280     if (runnumber == -1)
0281     {
0282         h1D_map.insert(
0283             std::make_pair(
0284                 "h1D_INTTvtxZ_resolution_Inclusive100",
0285                 new TH1D("h1D_INTTvtxZ_resolution_Inclusive100", "h1D_INTTvtxZ_resolution_Inclusive100;INTTvtxZ - Truth_{z} [cm];Entries (/0.1)",100, -2.5,2.5)
0286             )
0287         );
0288         h1D_map.insert(
0289             std::make_pair(
0290                 "h1D_INTTvtxZ_resolution_Inclusive70",
0291                 new TH1D("h1D_INTTvtxZ_resolution_Inclusive70", "h1D_INTTvtxZ_resolution_Inclusive70;INTTvtxZ - Truth_{z} [cm];Entries (/0.1)",100, -2.5,2.5)
0292             )
0293         );
0294         h1D_map.insert(
0295             std::make_pair(
0296                 "h1D_INTTvtxZ_resolution_HighNClus",
0297             new TH1D("h1D_INTTvtxZ_resolution_HighNClus", "h1D_INTTvtxZ_resolution_HighNClus;INTTvtxZ - Truth_{z} [cm];Entries (/0.1)",100, -2.5, 2.5)
0298         )
0299         );
0300     }
0301
0302     for (auto &hist : h1D_map)
0303     {
0304         std::string YTitle = hist.second -> GetYaxis() -> GetTitle();
0305         std::string XTitle = hist.second -> GetXaxis() -> GetTitle();
0306
0307         std::string YTitle_post;
0308
0309         if (YTitle.find("Entries") != std::string::npos) // note : found the (Entries)
0310         {
0311             YTitle_post = Form("Entries  (/%.2f)", hist.second -> GetBinWidth(1));
0312             hist.second -> GetYaxis() -> SetTitle(YTitle_post.c_str());
0313         }
0314     }
0315
0316     // Division: ---For TH2D---------------------------------------------------------------------------------------
0317     h2D_map.clear();
0318
0319     h2D_map.insert(
0320         std::make_pair(
0321             "h2D_INTTz_MBDz_Inclusive100",
0322             new TH2D("h2D_INTTz_MBDz_Inclusive100","h2D_INTTz_MBDz_Inclusive100;INTT vtxZ [cm]; MBD vtxZ [cm]",200, vtxZ_range.first, vtxZ_range.second, 200, vtxZ_range.first, vtxZ_range.second)
0323         )
0324     );
0325
0326     h2D_map.insert(
0327         std::make_pair(
0328             "h2D_INTTz_MBDz_Inclusive70",
0329             new TH2D("h2D_INTTz_MBDz_Inclusive70","h2D_INTTz_MBDz_Inclusive70;INTT vtxZ [cm]; MBD vtxZ [cm]",200, vtxZ_range.first, vtxZ_range.second, 200, vtxZ_range.first, vtxZ_range.second)
0330         )
0331     );
0332
0333     h2D_map.insert(
0334         std::make_pair(
0335             "h2D_INTTz_MBDz_HighNClus",
0336             new TH2D("h2D_INTTz_MBDz_HighNClus","h2D_INTTz_MBDz_HighNClus;INTT vtxZ [cm]; MBD vtxZ [cm]",200, vtxZ_range.first, vtxZ_range.second, 200, vtxZ_range.first, vtxZ_range.second)
0337         )
0338     );
0339
0340     h2D_map.insert(
0341         std::make_pair(
0342             "h2D_MBDvtxZ_MBDINTTvtxZDiff_HighNClus",
0343             new TH2D("h2D_MBDvtxZ_MBDINTTvtxZDiff_HighNClus","h2D_MBDvtxZ_MBDINTTvtxZDiff_HighNClus;MBD vtxZ [cm]; MBDz - INTTz [cm]",200, vtxZ_range.first, vtxZ_range.second, 200, -15, 15)
0344         )
0345     );
0346
0347     h2D_map.insert(
0348         std::make_pair(
0349             "h2D_NClus_INTTvtxZ",
0350             new TH2D("h2D_NClus_INTTvtxZ","h2D_NClus_INTTvtxZ;NClus (INTT);INTT vtxZ [cm]",200, 0, 10000, 200, vtxZ_range.first, vtxZ_range.second)
0351         )
0352     );
0353
0354     h2D_map.insert(
0355         std::make_pair(
0356             "h2D_NClus_TrapezoidalFitWidth",
0357             new TH2D("h2D_NClus_TrapezoidalFitWidth","h2D_NClus_TrapezoidalFitWidth;NClus;INTT vtxZ Dist Fit Width [cm]", 100, 0, 10000, 100, 0, 15)
0358         )
0359     );
0360
0361     h2D_map.insert(
0362         std::make_pair(
0363             "h2D_NClus_TrapezoidalFWHM",
0364             new TH2D("h2D_NClus_TrapezoidalFWHM","h2D_NClus_TrapezoidalFWHM;NClus;INTT vtxZ Dist FWHM [cm]", 100, 0, 10000, 100, 0, 15)
0365         )
0366     );
0367
0368     h2D_map.insert(
0369         std::make_pair(
0370             "h2D_NClus_INTTvtxZError",
0371             new TH2D("h2D_NClus_INTTvtxZError","h2D_NClus_INTTvtxZError;NClus;INTT vtxZ StdDev [cm]", 100, 0, 10000, 100, 0, 15)
0372         )
0373     );
0374
0375     h2D_map.insert(
0376         std::make_pair(
0377             "h2D_confirm_InttNClus_MbdChargeSum",
0378             new TH2D("h2D_confirm_InttNClus_MbdChargeSum","h2D_confirm_InttNClus_MbdChargeSum;INTT NClus;MBD charge sum", 200, 0, 10000, 200, 0, 2500)
0379         )
0380     );
0381
0382     h2D_map.insert(
0383         std::make_pair(
0384             "h2D_confirm_InttInnerOuterNClus",
0385             new TH2D("h2D_confirm_InttInnerOuterNClus","h2D_confirm_InttInnerOuterNClus;NClus (inner);NClus (outer)", 200, 0, 5500, 200, 0, 5500)
0386         )
0387     );
0388
0389
0390     // note : for MC
0391     if (runnumber == -1)
0392     {
0393         h2D_map.insert(
0394             std::make_pair(
0395                 "h2D_NClus_Zresolution",
0396                 new TH2D("h2D_NClus_Zresolution","h2D_NClus_Zresolution;NClus;INTTvtxZ - Truth_{z} [cm]",100,0,10000,200,-4,4)
0397             )
0398         );
0399
0400         h2D_map.insert(
0401             std::make_pair(
0402                 "h2D_truthZ_Zresolution_Inclusive100",
0403                 new TH2D("h2D_truthZ_Zresolution_Inclusive100","h2D_truthZ_Zresolution_Inclusive100;Truth_{z} [cm];INTTvtxZ - Truth_{z} [cm]",100,-50,50,200,-4,4)
0404             )
0405         );
0406
0407         h2D_map.insert(
0408             std::make_pair(
0409                 "h2D_truthZ_Zresolution_HighNClus",
0410                 new TH2D("h2D_truthZ_Zresolution_HighNClus","h2D_truthZ_Zresolution_HighNClus;Truth_{z} [cm];INTTvtxZ - Truth_{z} [cm]",100,-50,50,200,-4,4)
0411             )
0412         );
0413
0414         h2D_map.insert(
0415             std::make_pair(
0416                 "h2D_GoodRecoZ_TruthZ_Centrality",
0417                 new TH2D("h2D_GoodRecoZ_TruthZ_Centrality","h2D_GoodRecoZ_TruthZ_Centrality;Truth_{z} [cm];Centrality [%]",nVtxZBin, VtxZEdge_min, VtxZEdge_max, nCentrality_bin, &centrality_edges[0])
0418             )
0419         );
0420
0421         h2D_map.insert(
0422             std::make_pair(
0423                 "h2D_TruthCount_TruthZ_Centrality",
0424                 new TH2D("h2D_TruthCount_TruthZ_Centrality","h2D_TruthCount_TruthZ_Centrality;Truth_{z} [cm];Centrality [%]",nVtxZBin, VtxZEdge_min, VtxZEdge_max, nCentrality_bin, &centrality_edges[0])
0425             )
0426         );
0427
0428         h2D_map.insert(
0429             std::make_pair(
0430                 "h2D_RecoZEffi_TruthZ_Centrality",
0431                 new TH2D("h2D_RecoZEffi_TruthZ_Centrality","h2D_RecoZEffi_TruthZ_Centrality;Truth_{z} [cm];Centrality [%]",nVtxZBin, VtxZEdge_min, VtxZEdge_max, nCentrality_bin, &centrality_edges[0])
0432             )
0433         );
0434
0435     }
0436
0437 }
0438
0439 void vtxZDist::PrepareEvent()
0440 {
0441     for (int i = 0; i < run_nEvents; i++)
0442     {
0443         tree_in -> GetEntry(i);
0444
0445         if (i % 500 == 0) {std::cout << "Processing Event " << i << " / " << run_nEvents << std::endl;}
0446
0447         // =======================================================================================================================================================
0448         // note : hard cut
0449
0450         // note : for data
0451         if (runnumber != -1 && ApplyEvtBcoFullDiffCut.first && InttBcoFullDiff_next <= ApplyEvtBcoFullDiffCut.second) {continue;}
0452         if (runnumber != -1 && MBDNSg2 != 1) {continue;} // todo: assume MC no trigger
0453
0454         // note : for MC
0455         // if (runnumber == -1 && NTruthVtx != 1) {continue;}
0456
0457         // note : both data and MC
0458         if (is_min_bias != 1) {continue;}
0459         if (MBD_z_vtx != MBD_z_vtx) {continue;}
0460         if (MBD_centrality != MBD_centrality) {continue;}
0461         if (MBD_centrality < 0 || MBD_centrality > 100) {continue;}
0462         if (INTTvtxZ != INTTvtxZ) {continue;}
0463         if (MBD_z_vtx < cut_GlobalMBDvtxZ.first || MBD_z_vtx > cut_GlobalMBDvtxZ.second) {continue;} // todo: the hard cut 60 cm
0464
0465         // =======================================================================================================================================================
0466         // note : optional cut
0467         if (Apply_cut && (MBD_z_vtx - INTTvtxZ < cut_vtxZDiff.first || MBD_z_vtx - INTTvtxZ > cut_vtxZDiff.second) ) {continue;}
0468         if (Apply_cut && (TrapezoidalFitWidth < cut_TrapezoidalFitWidth.first || TrapezoidalFitWidth > cut_TrapezoidalFitWidth.second)){continue;}
0469         if (Apply_cut && (TrapezoidalFWHM < cut_TrapezoidalFWHM.first || TrapezoidalFWHM > cut_TrapezoidalFWHM.second)){continue;}
0470         if (Apply_cut && (INTTvtxZError < cut_INTTvtxZError.first || INTTvtxZError > cut_INTTvtxZError.second)){continue;}
0471
0472         // =======================================================================================================================================================
0473         double INTTvtxZWeighting;
0474         if (ApplyVtxZReWeighting && h1D_INTT_vtxZ_reweighting != nullptr){
0475             INTTvtxZWeighting = h1D_INTT_vtxZ_reweighting -> GetBinContent(h1D_INTT_vtxZ_reweighting -> FindBin(INTTvtxZ));
0476         }
0477         else if (ApplyVtxZReWeighting && h1D_INTT_vtxZ_reweighting == nullptr){
0478             std::cout << "ApplyVtxZReWeighting is true, but h1D_INTT_vtxZ_reweighting is nullptr" << std::endl;
0479             exit(1);
0480         }
0481         else {
0482             INTTvtxZWeighting = 1.0;
0483         }
0484
0485         int Mbin = h1D_centrality_bin -> Fill(MBD_centrality, INTTvtxZWeighting);
0486         Mbin = (Mbin == -1) ? -1 : Mbin - 1;
0487         if (Mbin == -1) {
0488             std::cout << "Mbin == -1, MBD_centrality = " << MBD_centrality << std::endl;
0489             continue;
0490         }
0491
0492         if (i % 500 == 0) {std::cout << "MBD_centrality: "<< MBD_centrality <<", Mbin: " << Mbin << ", nCentrality_bin: " << nCentrality_bin << std::endl;}
0493
0494
0495         // note : inclusive100
0496         h2D_map["h2D_confirm_InttNClus_MbdChargeSum"] -> Fill(NClus, MBD_charge_sum, INTTvtxZWeighting);
0497         h2D_map["h2D_confirm_InttInnerOuterNClus"] -> Fill(NClus_Layer1, NClus - NClus_Layer1, INTTvtxZWeighting);
0498
0499         h1D_map["h1D_INTTz_Inclusive100"] -> Fill(INTTvtxZ, INTTvtxZWeighting);
0500         h1D_map["h1D_VtxZDiff_MBD_INTT_Inclusive100"] -> Fill(MBD_z_vtx - INTTvtxZ, INTTvtxZWeighting);
0501         h2D_map["h2D_INTTz_MBDz_Inclusive100"] -> Fill(INTTvtxZ, MBD_z_vtx, INTTvtxZWeighting);
0502
0503         h1D_map["h1D_TrapezoidalFitWidth"] -> Fill(TrapezoidalFitWidth, INTTvtxZWeighting);
0504         h1D_map["h1D_TrapezoidalFWHM"] -> Fill(TrapezoidalFWHM, INTTvtxZWeighting);
0505         h1D_map["h1D_INTTvtxZError"] -> Fill(INTTvtxZError, INTTvtxZWeighting);
0506         h1D_map["h1D_NgroupTrapezoidal"] -> Fill(NgroupTrapezoidal, INTTvtxZWeighting);
0507         h1D_map["h1D_NgroupCoarse"] -> Fill(NgroupCoarse, INTTvtxZWeighting);
0508
0509         h2D_map["h2D_NClus_INTTvtxZ"] -> Fill(NClus, INTTvtxZ, INTTvtxZWeighting);
0510         h2D_map["h2D_NClus_TrapezoidalFitWidth"] -> Fill(NClus, TrapezoidalFitWidth, INTTvtxZWeighting);
0511         h2D_map["h2D_NClus_TrapezoidalFWHM"] -> Fill(NClus, TrapezoidalFWHM, INTTvtxZWeighting);
0512         h2D_map["h2D_NClus_INTTvtxZError"] -> Fill(NClus, INTTvtxZError, INTTvtxZWeighting);
0513
0514
0515         h1D_map[Form("h1D_INTTz_Mbin%d", Mbin)] -> Fill(INTTvtxZ, INTTvtxZWeighting);
0516
0517
0518         // todo: the Semi_inclusive_bin
0519         if (Mbin <= Semi_inclusive_bin){
0520             h1D_map["h1D_INTTz_Inclusive70"] -> Fill(INTTvtxZ, INTTvtxZWeighting);
0521             h1D_map["h1D_VtxZDiff_MBD_INTT_Inclusive70"] -> Fill(MBD_z_vtx - INTTvtxZ, INTTvtxZWeighting);
0522             h2D_map["h2D_INTTz_MBDz_Inclusive70"] -> Fill(INTTvtxZ, MBD_z_vtx, INTTvtxZWeighting);
0523         }
0524
0525         if (NClus > HighNClus){
0526             h1D_map["h1D_VtxZDiff_MBD_INTT_HighNClus"]->Fill(MBD_z_vtx - INTTvtxZ, INTTvtxZWeighting);
0527             h2D_map["h2D_INTTz_MBDz_HighNClus"]->Fill(INTTvtxZ, MBD_z_vtx, INTTvtxZWeighting);
0528             h2D_map["h2D_MBDvtxZ_MBDINTTvtxZDiff_HighNClus"]->Fill(MBD_z_vtx, MBD_z_vtx - INTTvtxZ, INTTvtxZWeighting);
0529         }
0530
0531
0532         // note : for MC
0533         // todo : no vtxZ reweighting for the resolution study
0534         if (runnumber == -1){
0535             h1D_map["h1D_INTTvtxZ_resolution_Inclusive100"] -> Fill(INTTvtxZ - TruthPV_trig_z);
0536             h2D_map["h2D_truthZ_Zresolution_Inclusive100"] -> Fill(TruthPV_trig_z, INTTvtxZ - TruthPV_trig_z);
0537             h2D_map["h2D_NClus_Zresolution"]  -> Fill(NClus, INTTvtxZ - TruthPV_trig_z);
0538
0539             if (NClus > HighNClus){
0540                 h1D_map["h1D_INTTvtxZ_resolution_HighNClus"] -> Fill(INTTvtxZ - TruthPV_trig_z);
0541                 h2D_map["h2D_truthZ_Zresolution_HighNClus"] -> Fill(TruthPV_trig_z, INTTvtxZ - TruthPV_trig_z);
0542             }
0543
0544             if (Mbin <= Semi_inclusive_bin){
0545                 h1D_map["h1D_INTTvtxZ_resolution_Inclusive70"] -> Fill(INTTvtxZ - TruthPV_trig_z);
0546             }
0547
0548             // todo : no weight
0549             if (fabs(INTTvtxZ - TruthPV_trig_z) <= cut_GoodRecoVtxZ) {h2D_map["h2D_GoodRecoZ_TruthZ_Centrality"] -> Fill(TruthPV_trig_z, MBD_centrality);}
0550             h2D_map["h2D_TruthCount_TruthZ_Centrality"] -> Fill(TruthPV_trig_z, MBD_centrality);
0551         }
0552
0553         // =======================================================================================================================================================
0554         // todo: the vtxZ range cut for the analysis is based on INTT only
0555         if (INTTvtxZ < cut_AnaVtxZ.first || INTTvtxZ > cut_AnaVtxZ.second) {continue;}
0556         // if (MBD_z_vtx < cut_AnaVtxZ.first || MBD_z_vtx > cut_AnaVtxZ.second) {continue;}
0557
0558         h1D_map["h1D_INTTz_Inclusive100_tight"] -> Fill(INTTvtxZ, INTTvtxZWeighting);
0559
0560         if (Mbin <= Semi_inclusive_bin){
0561             h1D_map["h1D_INTTz_Inclusive70_tight"] -> Fill(INTTvtxZ, INTTvtxZWeighting);
0562         }
0563
0564     } // note : end of the event loop
0565 }
0566
0567 void vtxZDist::EndRun()
0568 {
0569     file_out -> cd();
0570     for (auto &hist : h1D_map)
0571     {
0572         hist.second -> Write();
0573     }
0574
0575     if (runnumber == -1) {h2D_map["h2D_RecoZEffi_TruthZ_Centrality"] -> Divide(h2D_map["h2D_GoodRecoZ_TruthZ_Centrality"], h2D_map["h2D_TruthCount_TruthZ_Centrality"]);}
0576
0577     for (auto &hist : h2D_map)
0578     {
0579         hist.second -> Write();
0580     }
0581
0582     if (h1D_INTT_vtxZ_reweighting != nullptr){
0583         h1D_INTT_vtxZ_reweighting -> Write();
0584     }
0585
0586     h1D_centrality_bin -> Write();
0587
0588     file_out -> Close();
0589 }