add read.C file

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();

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];

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

read_roots/read.C

@@ -0,0 +1,155 @@
+#include <iostream>
+#include <algorithm>
+#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();
+}
+