NewCode2024/TruthHadronCheck/quick_check_truedNdEta_chain.C

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