NewCode2024/BkgEstimation/BkgEstimation.C

0001 int BkgEstimation()
0002 {
0003     TFile * file_in = TFile::Open("/sphenix/tg/tg01/commissioning/INTT/work/cwshih/seflgendata/run_54280/bco_diff/completed/BcoDiffNtuple_54280_clonehitremoveBCO_hotchannelremove_hitQA_checkclonehit_merged.root");
0004
0005     std::pair<int, int> dead_region = {27, 54};
0006     long long  trigger_peak = 56;
0007
0008     TTree * tree = (TTree *)file_in->Get("tree");
0009     int nentries = tree->GetEntries();
0010     std::cout << "nentries = " << nentries << std::endl;
0011
0012     tree->SetBranchStatus("*", 0);
0013     tree->SetBranchStatus("MBDNS_tight_vtx30cm", 1);
0014     tree->SetBranchStatus("felix0_bcodiff", 1);
0015     tree->SetBranchStatus("felix1_bcodiff", 1);
0016     tree->SetBranchStatus("felix2_bcodiff", 1);
0017     tree->SetBranchStatus("felix3_bcodiff", 1);
0018     tree->SetBranchStatus("felix4_bcodiff", 1);
0019     tree->SetBranchStatus("felix5_bcodiff", 1);
0020     tree->SetBranchStatus("felix6_bcodiff", 1);
0021     tree->SetBranchStatus("felix7_bcodiff", 1);
0022
0023     int in_MBDNS_tight_vtx30cm;
0024     std::vector<int> * in_felix0_bcodiff = 0;
0025     std::vector<int> * in_felix1_bcodiff = 0;
0026     std::vector<int> * in_felix2_bcodiff = 0;
0027     std::vector<int> * in_felix3_bcodiff = 0;
0028     std::vector<int> * in_felix4_bcodiff = 0;
0029     std::vector<int> * in_felix5_bcodiff = 0;
0030     std::vector<int> * in_felix6_bcodiff = 0;
0031     std::vector<int> * in_felix7_bcodiff = 0;
0032
0033     tree->SetBranchAddress("MBDNS_tight_vtx30cm", &in_MBDNS_tight_vtx30cm);
0034     tree->SetBranchAddress("felix0_bcodiff", &in_felix0_bcodiff);
0035     tree->SetBranchAddress("felix1_bcodiff", &in_felix1_bcodiff);
0036     tree->SetBranchAddress("felix2_bcodiff", &in_felix2_bcodiff);
0037     tree->SetBranchAddress("felix3_bcodiff", &in_felix3_bcodiff);
0038     tree->SetBranchAddress("felix4_bcodiff", &in_felix4_bcodiff);
0039     tree->SetBranchAddress("felix5_bcodiff", &in_felix5_bcodiff);
0040     tree->SetBranchAddress("felix6_bcodiff", &in_felix6_bcodiff);
0041     tree->SetBranchAddress("felix7_bcodiff", &in_felix7_bcodiff);
0042
0043
0044     long long selected_event = 0;
0045     std::vector<TH1D *> all_felix_BcoDiffHist;
0046
0047     for (int i = 0; i < 8; i++)
0048     {
0049         // all_felix_BcoDiffHist.push_back((TH1D *)file_in->Get(Form("all_felix_BcoDiffHist_%d", i)));
0050         all_felix_BcoDiffHist.push_back(new TH1D(Form("all_felix_BcoDiffHist_%d", i), Form("all_felix_BcoDiffHist_%d;Time_bucket;Entries", i), 128,0,128));
0051     }
0052
0053     for (long long i = 0; i < nentries; i++){
0054         tree -> GetEntry(i);
0055
0056         if (i % 10000 == 0) {std::cout << "i = " << i << std::endl;}
0057
0058         if (in_MBDNS_tight_vtx30cm != 1) {continue;}
0059
0060         selected_event += 1;
0061
0062         for (int j = 0; j < 128; j++){
0063             all_felix_BcoDiffHist[0]->Fill(j, in_felix0_bcodiff->at(j));
0064             all_felix_BcoDiffHist[1]->Fill(j, in_felix1_bcodiff->at(j));
0065             all_felix_BcoDiffHist[2]->Fill(j, in_felix2_bcodiff->at(j));
0066             all_felix_BcoDiffHist[3]->Fill(j, in_felix3_bcodiff->at(j));
0067             all_felix_BcoDiffHist[4]->Fill(j, in_felix4_bcodiff->at(j));
0068             all_felix_BcoDiffHist[5]->Fill(j, in_felix5_bcodiff->at(j));
0069             all_felix_BcoDiffHist[6]->Fill(j, in_felix6_bcodiff->at(j));
0070             all_felix_BcoDiffHist[7]->Fill(j, in_felix7_bcodiff->at(j));
0071         }
0072     }
0073
0074     std::cout << "selected_event = " << selected_event << std::endl;
0075
0076     TFile * file_out = TFile::Open("BcoDiff_Hist_MBDNS30cm.root", "recreate");
0077     for (int i = 0; i < all_felix_BcoDiffHist.size(); i++)
0078     {
0079         all_felix_BcoDiffHist[i]->Write();
0080     }
0081     file_out->Close();
0082
0083     // std::vector<TH1D *> selected_felix_BcoDiffHist;
0084     // std::vector<TH1D * > bkg_count_vec; bkg_count_vec.clear();
0085
0086     // std::vector<TH1D *> selected_felix_BcoDiffHist_ratio;
0087     // std::vector<TH1D * > bkg_count_vec_ratio; bkg_count_vec_ratio.clear();
0088
0089     // for (TH1D * hist : all_felix_BcoDiffHist)
0090     // {
0091     //     selected_felix_BcoDiffHist.push_back((TH1D *)hist->Clone(Form("selected_%s", hist->GetName())));
0092
0093     //     bkg_count_vec.push_back(new TH1D(Form("BkgCount_%s",hist->GetName()),Form("BkgCount_%s;NInttRawHits (bkg);Entries",hist->GetName()),200,0.1 * pow(10,6), 5 * pow(10,6)));
0094     //     bkg_count_vec_ratio.push_back(new TH1D(Form("BkgCount_%s_ratio",hist->GetName()),Form("BkgCount_%s_ratio;NInttRawHits (per Bkg. Evt.);Entries",hist->GetName()),200,0.1 * pow(10,6) / double (nentries), 5 * pow(10,6) / double (nentries)));
0095
0096     //     for (int i = 0; i < selected_felix_BcoDiffHist.back()->GetNbinsX(); i++)
0097     //     {
0098     //         if (i+1 >= dead_region.first && i+1 <= dead_region.second) // note : dead region
0099     //         {
0100     //             selected_felix_BcoDiffHist.back()->SetBinContent(i + 1, 0);
0101     //         }
0102
0103     //         if (i+1 >= trigger_peak - 1 && i+1 <= trigger_peak + 1) // note : trigger peak
0104     //         {
0105     //             selected_felix_BcoDiffHist.back()->SetBinContent(i + 1, 0);
0106     //         }
0107
0108     //         if ( (i + 1) % 2 != 0)
0109     //         {
0110     //             selected_felix_BcoDiffHist.back()->SetBinContent(i + 1, 0);
0111     //         }
0112
0113     //         if (selected_felix_BcoDiffHist.back()->GetBinContent(i + 1) > 0)
0114     //         {
0115     //             bkg_count_vec.back()->Fill(selected_felix_BcoDiffHist.back()->GetBinContent(i + 1));
0116     //         }
0117
0118     //     }
0119
0120     //     selected_felix_BcoDiffHist_ratio.push_back((TH1D *)selected_felix_BcoDiffHist.back()->Clone(Form("selected_%s_ratio", hist->GetName())));
0121     //     selected_felix_BcoDiffHist_ratio.back()->Scale(1.0 / double(selected_event));
0122
0123     //     for (int j = 0; j < selected_felix_BcoDiffHist_ratio.back()->GetNbinsX(); j++)
0124     //     {
0125     //         if (j+1 >= dead_region.first && j+1 <= dead_region.second) // note : dead region
0126     //         {
0127     //             selected_felix_BcoDiffHist_ratio.back()->SetBinContent(j + 1, 0);
0128     //         }
0129     //     }
0130     // }
0131
0132
0133
0134
0135     // TFile * file_out = TFile::Open("BkgEstimation.root", "recreate");
0136     // for (int i = 0; i < all_felix_BcoDiffHist.size(); i++)
0137     // {
0138     //     all_felix_BcoDiffHist[i]->Write();
0139     // }
0140
0141     // for (int i = 0; i < selected_felix_BcoDiffHist.size(); i++)
0142     // {
0143     //     selected_felix_BcoDiffHist[i]->Write();
0144     // }
0145
0146     // THStack * hstack1D_total_bkg_count = new THStack("hstack1D_total_bkg_count", "hstack1D_total_bkg_count;NInttRawHits (bkg);Entries");
0147
0148     // for (int i = 0; i < bkg_count_vec.size(); i++)
0149     // {
0150     //     bkg_count_vec[i]->Write();
0151     //     bkg_count_vec[i]->SetFillColor(i + 1);
0152     //     hstack1D_total_bkg_count->Add(bkg_count_vec[i]);
0153     // }
0154
0155     // hstack1D_total_bkg_count->Write();
0156     // ((TH1D*)hstack1D_total_bkg_count->GetStack()->Last())->Write("h1D_total_bkg_count");
0157
0158
0159
0160     // file_out->Close();
0161
0162     return 0;
0163 }