NewCode2024/TruthHadronCheck/quick_check_truedNdEta.C

0001 #include "/sphenix/user/ChengWei/INTT/INTT_dNdeta_repo/NewCode2024/Constants.cpp"
0002
0003 int quick_check_truedNdEta(
0004     int process_id,
0005     string input_directory,
0006     string input_filename,
0007     string output_directory
0008 )
0009 {
0010     TFile * file_in = TFile::Open(Form("%s/%s", input_directory.c_str(), input_filename.c_str()));
0011     TTree * tree_in = (TTree*)file_in->Get("EventTree");
0012     tree_in -> SetBranchStatus("*", 0);
0013
0014     tree_in -> SetBranchStatus("MBD_z_vtx", 1);
0015     tree_in -> SetBranchStatus("is_min_bias", 1);
0016     tree_in -> SetBranchStatus("MBD_centrality", 1);
0017     tree_in -> SetBranchStatus("MBD_south_charge_sum", 1);
0018     tree_in -> SetBranchStatus("MBD_north_charge_sum", 1);
0019     tree_in -> SetBranchStatus("MBD_charge_sum", 1);
0020     tree_in -> SetBranchStatus("MBD_charge_asymm", 1);
0021     // tree_in -> SetBranchStatus("InttBcoFullDiff_next", 1);
0022     // tree_in -> SetBranchStatus("MBDNSg2", 1);
0023
0024     // tree_in -> SetBranchStatus("ClusLayer", 1);
0025     // tree_in -> SetBranchStatus("ClusEta_INTTz", 1);
0026
0027     tree_in -> SetBranchStatus("INTTvtxZ", 1);
0028     tree_in -> SetBranchStatus("INTTvtxZError", 1);
0029     tree_in -> SetBranchStatus("NgroupTrapezoidal", 1);
0030     tree_in -> SetBranchStatus("NgroupCoarse", 1);
0031     tree_in -> SetBranchStatus("TrapezoidalFitWidth", 1);
0032     tree_in -> SetBranchStatus("TrapezoidalFWHM", 1);
0033
0034     tree_in -> SetBranchStatus("NTruthVtx", 1);
0035     tree_in -> SetBranchStatus("PrimaryG4P_Pt", 1);
0036     tree_in -> SetBranchStatus("PrimaryG4P_Eta", 1);
0037     tree_in -> SetBranchStatus("PrimaryG4P_Phi", 1);
0038     tree_in -> SetBranchStatus("PrimaryG4P_E", 1);
0039     tree_in -> SetBranchStatus("PrimaryG4P_PID", 1);
0040     tree_in -> SetBranchStatus("PrimaryG4P_isChargeHadron", 1);
0041
0042     tree_in -> SetBranchStatus("TruthPV_trig_z", 1);
0043
0044
0045
0046     float MBD_z_vtx;
0047     bool is_min_bias;
0048     float MBD_centrality;
0049     float MBD_south_charge_sum;
0050     float MBD_north_charge_sum;
0051     float MBD_charge_sum;
0052     float MBD_charge_asymm;
0053     // int InttBcoFullDiff_next;
0054
0055     // note : trigger tag
0056     // int MBDNSg2;
0057
0058     std::vector<int> *ClusLayer = 0;
0059     std::vector<double> *ClusEta_INTTz = 0;
0060
0061     // note : INTT vertex Z
0062     double INTTvtxZ;
0063     double INTTvtxZError;
0064     double NgroupTrapezoidal;
0065     double NgroupCoarse;
0066     double TrapezoidalFitWidth;
0067     double TrapezoidalFWHM;
0068
0069     int NTruthVtx;
0070     float TruthPV_trig_z;
0071     std::vector<float> *PrimaryG4P_Pt = 0;
0072     std::vector<float> *PrimaryG4P_Eta = 0;
0073     std::vector<float> *PrimaryG4P_Phi = 0;
0074     std::vector<float> *PrimaryG4P_E = 0;
0075     std::vector<float> *PrimaryG4P_PID = 0;
0076     std::vector<int> *PrimaryG4P_isChargeHadron = 0;
0077
0078     tree_in -> SetBranchAddress("MBD_z_vtx", &MBD_z_vtx);
0079     tree_in -> SetBranchAddress("is_min_bias", &is_min_bias);
0080     tree_in -> SetBranchAddress("MBD_centrality", &MBD_centrality);
0081     tree_in -> SetBranchAddress("MBD_south_charge_sum", &MBD_south_charge_sum);
0082     tree_in -> SetBranchAddress("MBD_north_charge_sum", &MBD_north_charge_sum);
0083     tree_in -> SetBranchAddress("MBD_charge_sum", &MBD_charge_sum);
0084     tree_in -> SetBranchAddress("MBD_charge_asymm", &MBD_charge_asymm);
0085     // tree_in -> SetBranchAddress("InttBcoFullDiff_next", &InttBcoFullDiff_next);
0086
0087     // tree_in -> SetBranchAddress("MBDNSg2", &MBDNSg2);
0088
0089     // tree_in -> SetBranchAddress("ClusLayer", &ClusLayer);
0090     // tree_in -> SetBranchAddress("ClusEta_INTTz", &ClusEta_INTTz);
0091
0092     tree_in -> SetBranchAddress("INTTvtxZ", &INTTvtxZ);
0093     tree_in -> SetBranchAddress("INTTvtxZError", &INTTvtxZError);
0094     tree_in -> SetBranchAddress("NgroupTrapezoidal", &NgroupTrapezoidal);
0095     tree_in -> SetBranchAddress("NgroupCoarse", &NgroupCoarse);
0096     tree_in -> SetBranchAddress("TrapezoidalFitWidth", &TrapezoidalFitWidth);
0097     tree_in -> SetBranchAddress("TrapezoidalFWHM", &TrapezoidalFWHM);
0098
0099     tree_in -> SetBranchAddress("NTruthVtx", &NTruthVtx);
0100     tree_in -> SetBranchAddress("PrimaryG4P_Pt", &PrimaryG4P_Pt);
0101     tree_in -> SetBranchAddress("PrimaryG4P_Eta", &PrimaryG4P_Eta);
0102     tree_in -> SetBranchAddress("PrimaryG4P_Phi", &PrimaryG4P_Phi);
0103     tree_in -> SetBranchAddress("PrimaryG4P_E", &PrimaryG4P_E);
0104     tree_in -> SetBranchAddress("PrimaryG4P_PID", &PrimaryG4P_PID);
0105     tree_in -> SetBranchAddress("PrimaryG4P_isChargeHadron", &PrimaryG4P_isChargeHadron);
0106     tree_in -> SetBranchAddress("TruthPV_trig_z", &TruthPV_trig_z);
0107
0108
0109     // Division : ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
0110
0111     std::string job_index = std::to_string( process_id );
0112     int job_index_len = 5;
0113     job_index.insert(0, job_index_len - job_index.size(), '0');
0114
0115     std::string output_filename = Form("INTTz_New_QuickCheck_%s.root", job_index.c_str());
0116
0117     system(Form("if [ -f %s/completed/%s ]; then rm %s/completed/%s; fi", output_directory.c_str(), output_filename.c_str(), output_directory.c_str(), output_filename.c_str()));
0118
0119     TFile * file_out = new TFile(Form("%s/%s", output_directory.c_str(), output_filename.c_str()), "recreate");
0120
0121     TH1D * h1D_NHadronCount = new TH1D("h1D_NHadronCount","h1D_NHadronCount;#eta;Entries",27,-2.7,2.7);
0122     TH1D * h1D_TrueEvtCount = new TH1D("h1D_TrueEvtCount","h1D_TrueEvtCount;#True vtx Z;Entries",20,-10,10);
0123     TH1D * h1D_TruedNdEta = new TH1D("h1D_TruedNdEta","h1D_TruedNdEta;#eta;Entries",27,-2.7,2.7);
0124
0125     TH1D * h1D_NHadronCount_NTrueVtxCut = new TH1D("h1D_NHadronCount_NTrueVtxCut","h1D_NHadronCount_NTrueVtxCut;#eta;Entries",27,-2.7,2.7);
0126     TH1D * h1D_TrueEvtCount_NTrueVtxCut = new TH1D("h1D_TrueEvtCount_NTrueVtxCut","h1D_TrueEvtCount_NTrueVtxCut;#True vtx Z;Entries",20,-10,10);
0127     TH1D * h1D_TruedNdEta_NTrueVtxCut = new TH1D("h1D_TruedNdEta_NTrueVtxCut","h1D_TruedNdEta_NTrueVtxCut;#eta;Entries",27,-2.7,2.7);
0128
0129     for (int i = 0; i < tree_in -> GetEntries(); i++)
0130     {
0131         tree_in -> GetEntry(i);
0132
0133         // if (i % 10 == 0)
0134         // {
0135         //     std::cout<<"i: "<<i<<std::endl;
0136         // }
0137
0138         int nHadron_count = 0;
0139
0140         // note : for MC
0141         // if (NTruthVtx != 1) {continue;}
0142
0143         // if (is_min_bias != 1) {continue;}
0144         // if (INTTvtxZ != INTTvtxZ) {continue;}
0145         // if (INTTvtxZ < -10 || INTTvtxZ >= 10) {continue;}
0146         // if (MBD_centrality < 0) {continue;}
0147         // if (MBD_centrality > 3) {continue;} // note : [0-3)% centrality interval
0148
0149
0150          // note : both data and MC
0151         if (is_min_bias != 1) {continue;}
0152         if (MBD_z_vtx != MBD_z_vtx) {continue;}
0153         if (MBD_centrality != MBD_centrality) {continue;}
0154         if (MBD_centrality < 0 || MBD_centrality > 100) {continue;}
0155         if (INTTvtxZ != INTTvtxZ) {continue;}
0156         // =======================================================================================================================================================
0157         // note : optional cut
0158         if ((MBD_z_vtx - INTTvtxZ < Constants::cut_vtxZDiff.first || MBD_z_vtx - INTTvtxZ > Constants::cut_vtxZDiff.second) ) {continue;}
0159         if ((TrapezoidalFitWidth < Constants::cut_TrapezoidalFitWidth.first || TrapezoidalFitWidth > Constants::cut_TrapezoidalFitWidth.second)){continue;}
0160         if ((TrapezoidalFWHM < Constants::cut_TrapezoidalFWHM.first || TrapezoidalFWHM > Constants::cut_TrapezoidalFWHM.second)){continue;}
0161         if ((INTTvtxZError < Constants::cut_INTTvtxZError.first || INTTvtxZError > Constants::cut_INTTvtxZError.second)){continue;}
0162         // =======================================================================================================================================================
0163
0164         if (INTTvtxZ < -10 || INTTvtxZ >= 10) {continue;}
0165         if (MBD_centrality < 0) {continue;}
0166         if (MBD_centrality > 3) {continue;} // note : [0-3)% centrality interval
0167
0168
0169         // if (MBD_z_vtx < -10 || MBD_z_vtx >= 10) {continue;}
0170
0171         // note : both data and MC
0172         // if (MBD_z_vtx != MBD_z_vtx) {continue;}
0173         // if (MBD_centrality != MBD_centrality) {continue;}
0174         // if (MBD_centrality < 0 || MBD_centrality > 100) {continue;}
0175         // if (INTTvtxZ != INTTvtxZ) {continue;}
0176         // if (MBD_z_vtx < -60 || MBD_z_vtx > 60) {continue;} // todo: the hard cut 60 cm
0177
0178
0179         // if ((MBD_z_vtx - INTTvtxZ < Constants::cut_vtxZDiff.first || MBD_z_vtx - INTTvtxZ > Constants::cut_vtxZDiff.second) ) {continue;}
0180         // if ((TrapezoidalFitWidth < Constants::cut_TrapezoidalFitWidth.first || TrapezoidalFitWidth > Constants::cut_TrapezoidalFitWidth.second)){continue;}
0181         // if ((TrapezoidalFWHM < Constants::cut_TrapezoidalFWHM.first || TrapezoidalFWHM > Constants::cut_TrapezoidalFWHM.second)){continue;}
0182         // if ((INTTvtxZError < Constants::cut_INTTvtxZError.first || INTTvtxZError > Constants::cut_INTTvtxZError.second)){continue;}
0183
0184         // if (TruthPV_trig_z < -10 || TruthPV_trig_z >= 10) {continue;}
0185         // if (MBD_centrality > 0.7) {continue;}
0186
0187         // if (MBD_centrality > 3) {continue;} // note : [0-3)% centrality interval
0188
0189         for (int true_i = 0; true_i < PrimaryG4P_isChargeHadron -> size(); true_i++)
0190         {
0191             if (PrimaryG4P_isChargeHadron -> at(true_i) == 1)
0192             {
0193                 h1D_NHadronCount -> Fill(PrimaryG4P_Eta -> at(true_i));
0194
0195                 if (fabs(PrimaryG4P_Eta->at(true_i)) < 4){
0196                     nHadron_count++;
0197                 }
0198             }
0199         }
0200
0201         std::cout<<"eID: "<<i<<", nHadron_count: "<<nHadron_count<<", TruthPV_trig_z: "<<TruthPV_trig_z<<std::endl;
0202
0203         h1D_TrueEvtCount -> Fill(INTTvtxZ);
0204
0205
0206         if (NTruthVtx != 1) {continue;}
0207
0208
0209         for (int true_i = 0; true_i < PrimaryG4P_isChargeHadron -> size(); true_i++)
0210         {
0211             if (PrimaryG4P_isChargeHadron -> at(true_i) == 1)
0212             {
0213                 h1D_NHadronCount_NTrueVtxCut -> Fill(PrimaryG4P_Eta -> at(true_i));
0214
0215                 // if (fabs(PrimaryG4P_Eta->at(true_i)) < 4){
0216                 //     nHadron_count++;
0217                 // }
0218             }
0219         }
0220
0221         h1D_TrueEvtCount_NTrueVtxCut -> Fill(INTTvtxZ);
0222     }
0223
0224     h1D_TruedNdEta = (TH1D*) h1D_NHadronCount -> Clone("h1D_TruedNdEta");
0225     h1D_TruedNdEta -> Sumw2(true);
0226     h1D_TruedNdEta -> Scale(1./h1D_TrueEvtCount->Integral());
0227     h1D_TruedNdEta -> Scale(1./ h1D_TruedNdEta->GetBinWidth(1));
0228
0229     file_out -> cd();
0230     h1D_NHadronCount -> Write();
0231     h1D_TrueEvtCount -> Write();
0232     h1D_TruedNdEta -> Write();
0233
0234     h1D_NHadronCount_NTrueVtxCut -> Write();
0235     h1D_TrueEvtCount_NTrueVtxCut -> Write();
0236
0237     file_out -> Close();
0238
0239     system(Form("mv %s/%s %s/completed/%s", output_directory.c_str(), output_filename.c_str(), output_directory.c_str(), output_filename.c_str()));
0240
0241     return 0;
0242
0243 }