diff --git a/analysis/app/analyzer_histograms.C b/analysis/app/analyzer_histograms.C index 9dcb19b..bc90cff 100644 --- a/analysis/app/analyzer_histograms.C +++ b/analysis/app/analyzer_histograms.C @@ -58,11 +58,17 @@ void analyzer_histograms::InitHistos(){ h_nDTRechits [i] = InitTH1F("h_nDTRechits", "h_nDTRechits", 300, 0, 300); + + h_dtRechitClusterX [i] = InitTH1F("h_dtRechitClusterX", "h_dtRechitClusterX", 50, -800, 800); + h_dtRechitClusterY [i] = InitTH1F("h_dtRechitClusterY", "h_dtRechitClusterY", 50, -800, 800); + h_dtRechitClusterZ [i] = InitTH1F("h_dtRechitClusterZ", "h_dtRechitClusterZ", 50, -700, 700); + h_dtRechitClusterSize [i] = InitTH1F("h_dtRechitClusterSize", "h_dtRechitClusterSize", 450, 50, 500); h_dtRechitClusterSize_FailPass [i] = InitTH1F("h_dtRechitClusterSize_FailPass", "h_dtRechitClusterSize_FailPass", 2, -0.5, 1.5); h_dtRechitClusterSize_FailPass_uw [i] = InitTH1F("h_dtRechitClusterSize_FailPass_uw", "h_dtRechitClusterSize_FailPass_uw", 2, -0.5, 1.5); h_dtRechitClusterSize_v [i] = InitTH1F("h_dtRechitClusterSize_v", "h_dtRechitClusterSize_v", n_b-1, x_bins); h_dtRechitClusterSize_v2 [i] = InitTH1F("h_dtRechitClusterSize_v2", "h_dtRechitClusterSize_v2", n_b2-1, x_bins2); + h_dtRechitClusterPhi [i] = InitTH1F("h_dtRechitClusterPhi", "h_dtRechitClusterPhi", 40, -3.5, 3.5); h_dtRechitClusterEta [i] = InitTH1F("h_dtRechitClusterEta", "h_dtRechitClusterEta", 40, -5., 5); h_dtRechitClusterMuonVetoPt [i] = InitTH1F("h_dtRechitClusterMuonVetoPt", "h_dtRechitClusterMuonVetoPt", 100, 0.,200.); @@ -114,6 +120,13 @@ void analyzer_histograms::FillHistos(int selbin, Float_t ew){ int d = DtClusterPassSel_all[selbin][i]; double dPhi = -999; if(muon_list.size()>0) dPhi = DeltaPhi(lepPhi[muon_list[0]], dtRechitClusterPhi[d]); + + h_dtRechitClusterDPhiLeadMuon [selbin]->Fill(dPhi); + h_dtRechitCluster_match_RPCBx_dPhi0p5 [selbin]->Fill(dtRechitCluster_match_RPCBx_dPhi0p5 [d]); + h_dtRechitClusterX [selbin]->Fill(dtRechitClusterX [d]); + h_dtRechitClusterY [selbin]->Fill(dtRechitClusterY [d]); + h_dtRechitClusterZ [selbin]->Fill(dtRechitClusterZ [d]); + h_dtRechitClusterDPhiLeadMuon [selbin]->Fill(dPhi, ew); h_dtRechitCluster_match_RPCBx_dPhi0p5 [selbin]->Fill(dtRechitCluster_match_RPCBx_dPhi0p5 [d], ew); h_dtRechitClusterSize [selbin]->Fill(dtRechitClusterSize [d], ew); @@ -135,6 +148,7 @@ void analyzer_histograms::FillHistos(int selbin, Float_t ew){ h_dtRechitCluster_match_RPCTime_sameStation_dR0p4 [selbin]->Fill(dtRechitCluster_match_RPCTime_sameStation_dR0p4 [d], ew); h_dtRechitCluster_match_RPCTimeSpread_sameStation_dR0p4 [selbin]->Fill(dtRechitCluster_match_RPCTimeSpread_sameStation_dR0p4[d], ew); h_dtRechitClusterMaxStation [selbin]->Fill(dtRechitClusterMaxStation [d], ew); + } } @@ -163,6 +177,9 @@ void analyzer_histograms::WriteHistos(int selbin){ h_dtRechitClusterDPhiLeadMuon [selbin]->Write(); h_dtRechitCluster_match_RPCBx_dPhi0p5 [selbin]->Write(); h_nDTRechits [selbin]->Write(); + h_dtRechitClusterX [selbin]->Write(); + h_dtRechitClusterY [selbin]->Write(); + h_dtRechitClusterZ [selbin]->Write(); h_dtRechitClusterSize [selbin]->Write(); h_dtRechitClusterSize_FailPass [selbin]->Write(); h_dtRechitClusterSize_FailPass_uw [selbin]->Write(); @@ -206,6 +223,9 @@ void analyzer_histograms::DeleteHistos(int selbin){ h_dtRechitClusterDPhiLeadMuon [selbin]->Delete(); h_dtRechitCluster_match_RPCBx_dPhi0p5 [selbin]->Delete(); h_nDTRechits [selbin]->Delete(); + h_dtRechitClusterX [selbin]->Delete(); + h_dtRechitClusterY [selbin]->Delete(); + h_dtRechitClusterZ [selbin]->Delete(); h_dtRechitClusterSize [selbin]->Delete(); h_dtRechitClusterSize_FailPass [selbin]->Delete(); h_dtRechitClusterSize_FailPass_uw [selbin]->Delete(); diff --git a/analysis/include/analyzer_histograms.h b/analysis/include/analyzer_histograms.h index ec214b2..907c48f 100644 --- a/analysis/include/analyzer_histograms.h +++ b/analysis/include/analyzer_histograms.h @@ -38,6 +38,9 @@ class analyzer_histograms : public analyzer_tree{ TH1F* h_dtRechitClusterDPhiLeadMuon [SELBINNAMESIZE]; TH1F* h_dtRechitCluster_match_RPCBx_dPhi0p5 [SELBINNAMESIZE]; TH1F* h_nDTRechits [SELBINNAMESIZE]; + TH1F* h_dtRechitClusterX [SELBINNAMESIZE]; + TH1F* h_dtRechitClusterY [SELBINNAMESIZE]; + TH1F* h_dtRechitClusterZ [SELBINNAMESIZE]; TH1F* h_dtRechitClusterSize [SELBINNAMESIZE]; TH1F* h_dtRechitClusterSize_FailPass [SELBINNAMESIZE]; TH1F* h_dtRechitClusterSize_FailPass_uw [SELBINNAMESIZE]; diff --git a/read_roots/README.md b/read_roots/README.md new file mode 100644 index 0000000..e302282 --- /dev/null +++ b/read_roots/README.md @@ -0,0 +1,5 @@ +The read.C file is use to read the g_llp ctau data in the root files, +find the intersection of two different decay length (ctau300 and ctau 1000) +Then can determine a re-weighting factor to extrapolate the behaviour of LLPs at the different decay lengths. + +Currently it calculates the intersection of ctau300 and ctau1000 for Pi0Pi0 category diff --git a/read_roots/read.C b/read_roots/read.C new file mode 100644 index 0000000..dab3604 --- /dev/null +++ b/read_roots/read.C @@ -0,0 +1,155 @@ +#include +#include +#include "TFile.h" +#include "TH1F.h" +#include "TTree.h" +#include "TCanvas.h" + +Float_t ctau_reweighter(Float_t t, Float_t tau0, Float_t tau1) +{ + Float_t numerator = (1.0f / tau1) * expf(-t / tau1); + Float_t denominator = (1.0f / tau0) * expf(-t / tau0); + return numerator / denominator; + +} + +void read() { + // Open the ROOT file + TFile *file1 = TFile::Open("BToKPhi_MuonLLPDecayGenFilter_PhiToPi0Pi0_mPhi1p0_ctau300.root"); + TFile *file2 = TFile::Open("BToKPhi_MuonLLPDecayGenFilter_PhiToPi0Pi0_mPhi1p0_ctau1000.root"); + + + if (!file1 || file1->IsZombie()) { + std::cerr << "Error opening file" << std::endl; + return; + } + if (!file2 || file2->IsZombie()) { + std::cerr << "Error opening file" << std::endl; + return; + } + + // Get the TTree from the file + TTree *tree_300 = nullptr; + file1->GetObject("MuonSystem", tree_300); + if (!tree_300) { + std::cerr << "Error retrieving MuonSystem TTree" << std::endl; + file1->Close(); + return; + } + TTree *tree_1000 = nullptr; + file2->GetObject("MuonSystem", tree_1000); + if (!tree_1000) { + std::cerr << "Error retrieving MuonSystem TTree" << std::endl; + file2->Close(); + return; + } + + // set your binning + TH1F *h_gLLP_ctau_300 = new TH1F("h_gLLP_ctau_300", "h_gLLP_ctau_300", 1000, 0, 1000); + TH1F *h_gLLP_ctau_1000 = new TH1F("h_gLLP_ctau_1000", "h_gLLP_ctau_1000", 1000, 0, 1000); + //try more bins + + + // Set branch address + float gLLP_ctau_300; + float gLLP_ctau_1000; + + tree_300->SetBranchAddress("gLLP_ctau", &gLLP_ctau_300); + tree_1000->SetBranchAddress("gLLP_ctau", &gLLP_ctau_1000); + + // Loop over the tree entries and fill the histogram + // Long64_t nEntries_300 = tree_300->GetEntries(); + // float w = 1.0;//weight + // for (Long64_t i = 0; i < nEntries_300; i++) { + // tree_300->GetEntry(i); + // h_gLLP_ctau_300->Fill(gLLP_ctau_300, w); + // } + //h_gLLP_ctau_300->Scale(1./h_gLLP_ctau_300->Integral(0,-1)); + Long64_t nEntries_300 = tree_300->GetEntries(); + + for (Long64_t i = 0; i < nEntries_300; i++) { + tree_300->GetEntry(i); + float w = 1.0; + float reweight_factor = ctau_reweighter(gLLP_ctau_300, 30, 20); + w *= reweight_factor; + h_gLLP_ctau_300->Fill(gLLP_ctau_300, w); + } + h_gLLP_ctau_300->Scale(1./h_gLLP_ctau_300->Integral(0,-1)); + + + Long64_t nEntries_1000 = tree_1000->GetEntries(); + for (Long64_t i = 0; i < nEntries_1000; i++) { + tree_1000->GetEntry(i); + float w = 1.0; + float reweight_factor = ctau_reweighter(gLLP_ctau_1000, 100, 90); + h_gLLP_ctau_1000->Fill(gLLP_ctau_1000, w); + + } + h_gLLP_ctau_1000->Scale(1./h_gLLP_ctau_1000->Integral(0,-1)); + + + + + float ratio_closest_to_1 = h_gLLP_ctau_300->GetBinContent(1)/h_gLLP_ctau_1000->GetBinContent(1); + int NumBin = 1; + int binNum_300 = h_gLLP_ctau_300->GetNbinsX(); + for (int i = 1; i<= binNum_300; i++) + { + if(h_gLLP_ctau_1000->GetBinContent(i) != 0) + { + if (abs(1-(h_gLLP_ctau_300->GetBinContent(i)/h_gLLP_ctau_1000->GetBinContent(i))) < abs(1-ratio_closest_to_1)) + { + ratio_closest_to_1 = h_gLLP_ctau_300->GetBinContent(i)/h_gLLP_ctau_1000->GetBinContent(i); + NumBin = i; + } + } + + } + float x_cor = h_gLLP_ctau_300->GetBinCenter(NumBin); + std::cout<<"x coordinate of the intersection is: " << x_cor << " with ratio: " << ratio_closest_to_1<< std::endl; + std::cout<<"The bin number is: " << NumBin << endl; + + + + // For Plotting + //combined graph + // TCanvas *c1 = new TCanvas("c1", "gLLP_ctau Distribution", 800, 600); + // h_gLLP_ctau_300->SetLineColor(kBlue); + // h_gLLP_ctau_1000->SetLineColor(kRed); + // h_gLLP_ctau_300->Draw(); + // h_gLLP_ctau_1000->Draw("same"); + + // TLatex *text = new TLatex(); + // text->SetNDC(); + // text->SetTextSize(0.04); + // text->SetTextAlign(22); + + // Convert float to string + // std::ostringstream oss_x; + // std::ostringstream oss_ratio; + // oss_x << x_cor; + // oss_ratio << ratio_closest_to_1; + // std::string x_cor_str = oss_x.str(); + // std::string ratio_str = oss_ratio.str(); + // std::string annotation = "x-intersection " + x_cor_str + " ratio: " + ratio_str; + // text->DrawLatex(0.5, 0.8, annotation.c_str()); + // c1->SaveAs("h_gLLP_PhiToPi0Pi0_ctau_combined.png"); + + //seperate histograms + TCanvas *c2 = new TCanvas("c2", "gLLP_ctau Distribution", 800, 600); + h_gLLP_ctau_300->Draw(); + + TCanvas *c3 = new TCanvas("c3", "gLLP_ctau Distribution", 800, 600); + h_gLLP_ctau_1000->Draw(); + + + c2->SaveAs("reweight200_h_gLLP_PhiToPi0Pi0_ctau300.png"); + c3->SaveAs("reweight900_h_gLLP_PhiToPi0Pi0_ctau1000.png"); + + // c3->SaveAs("h_gLLP_PhiToPi0Pi0_ctau1000.png"); + + // Close the ROOT file + file1->Close(); + file2->Close(); +} +