move old files, extend do_all

do_all.cpp

@@ -35,67 +35,78 @@
 #include <vector>
 
 template <typename T>
-class Cut {
+class Cut
+{
 public:
   virtual bool Check(T y) = 0;
   virtual std::string ToString() = 0;
 };
 
 template <typename T>
-class LT : public Cut<T> {
+class LT : public Cut<T>
+{
 private:
   T value;
 
 public:
-    LT(T val) : value{val} {
-
+  LT(T val) : value{val}
+  {
   }
 
-    bool Check(T y) {
+  bool Check(T y)
+  {
     return y < value;
   }
 
-    std::string ToString() {
+  std::string ToString()
+  {
     return TString::Format("(x < %s)", std::to_string(value).c_str()).Data();
   }
 };
 
 template <typename T>
-class GT : public Cut<T> {
+class GT : public Cut<T>
+{
 private:
   T value;
 
 public:
   GT(T val) : value{val} {}
 
-    bool Check(T y) {
+  bool Check(T y)
+  {
     return y > value;
   }
 
-    std::string ToString() {
+  std::string ToString()
+  {
     return TString::Format("(x > %s)", std::to_string(value).c_str()).Data();
   }
 };
 
 template <typename T>
-class EQ : public Cut<T> {
+class EQ : public Cut<T>
+{
 private:
   T value;
 
 public:
   EQ(T val) : value{val} {}
 
-    bool Check(T y) {
+  bool Check(T y)
+  {
     return y == value;
   }
 
-    std::string ToString() {
+  std::string ToString()
+  {
     return TString::Format("(x == %s)", std::to_string(value).c_str()).Data();
   }
 };
 
 template <typename T>
-class BT : public Cut<T> {
+class BT : public Cut<T>
+{
 private:
   T low;
   T up;
@@ -103,17 +114,20 @@ class BT : public Cut<T> {
 public:
   BT(T lower, T upper) : low{lower}, up{upper} {}
 
-    bool Check(T y) {
+  bool Check(T y)
+  {
     return low < y && y < up;
   }
 
-    std::string ToString() {
+  std::string ToString()
+  {
     return TString::Format("(x in [%s, %s])", std::to_string(low).c_str(), std::to_string(up).c_str()).Data();
   }
 };
 
 template <typename T>
-class Variable {
+class Variable
+{
 private:
   const char *name;
   std::vector<Cut<T> *> cuts;
@@ -121,113 +135,362 @@ class Variable {
   T value;
 
 public:
-
-    Variable() {
+  Variable()
+  {
   }
 
-    Variable(const char* name, std::vector<Cut<T>*> cuts) {
+  Variable(const char *name, std::vector<Cut<T> *> cuts)
+  {
     this->name = name;
     this->cuts = cuts;
   }
 
-    T* GetValuePtr() {
+  T *GetValuePtr()
+  {
     return &value;
   }
 
-    T GetValue() const {
+  T GetValue() const
+  {
     return value;
   }
+
+  bool CheckCuts() const
+  {
+    for (int i = 0; i < cuts.size(); i++)
+    {
+      if (!cuts[i]->Check(value))
+      {
+        return false;
+      }
+    }
+    return true;
+  }
 };
 
-class VariableCollection {
+class VariableCollection
+{
 private:
   std::map<std::string, Variable<Double_t>> doubles;
   std::map<std::string, Variable<Float_t>> floats;
   std::map<std::string, Variable<Int_t>> ints;
   std::map<std::string, Variable<Bool_t>> bools;
+  std::vector<std::vector<std::string>> trigger_cuts;
+
+  TLorentzVector ComputeFourMomentum(const char *name, bool do_mass_sub, double subbed_mass)
+  {
+    std::string pxname = TString::Format("%s_PX", name).Data();
+    std::string pyname = TString::Format("%s_PY", name).Data();
+    std::string pzname = TString::Format("%s_PZ", name).Data();
+    std::string ename = TString::Format("%s_ENERGY", name).Data();
+
+    auto itx = floats.find(pxname);
+    if (itx == floats.end())
+    {
+      return TLorentzVector{};
+    }
+
+    auto ity = floats.find(pyname);
+    if (ity == floats.end())
+    {
+      return TLorentzVector{};
+    }
+
+    auto itz = floats.find(pzname);
+    if (itz == floats.end())
+    {
+      return TLorentzVector{};
+    }
+
+    TVector3 momentum(itx->second.GetValue(), ity->second.GetValue(), itz->second.GetValue());
+    double energy = 0;
+
+    if (do_mass_sub)
+    {
+      energy = TMath::Sqrt(TMath::Sq(subbed_mass) + momentum.Mag2());
+    }
+    else
+    {
+      auto ite = floats.find(ename);
+      if (ite == floats.end())
+      {
+        return TLorentzVector{};
+      }
+      energy = ite->second.GetValue();
+    }
+
+    return TLorentzVector(momentum, energy);
+  }
 
 public:
-    VariableCollection() {
+  VariableCollection()
+  {
     doubles = std::map<std::string, Variable<Double_t>>{};
     floats = std::map<std::string, Variable<Float_t>>{};
     ints = std::map<std::string, Variable<Int_t>>{};
     bools = std::map<std::string, Variable<Bool_t>>{};
+    trigger_cuts = std::vector<std::vector<std::string>>{};
   }
 
-    void Connect(TChain *chain) {
-      for (auto& [key, var] : doubles) {
+  void Connect(TChain *chain)
+  {
+    for (auto &[key, var] : doubles)
+    {
       chain->SetBranchAddress(key.c_str(), var.GetValuePtr());
     }
 
-      for (auto& [key, var] : floats) {
+    for (auto &[key, var] : floats)
+    {
       chain->SetBranchAddress(key.c_str(), var.GetValuePtr());
     }
 
-      for (auto& [key, var] : ints) {
+    for (auto &[key, var] : ints)
+    {
       chain->SetBranchAddress(key.c_str(), var.GetValuePtr());
     }
 
-      for (auto& [key, var] : bools) {
+    for (auto &[key, var] : bools)
+    {
       chain->SetBranchAddress(key.c_str(), var.GetValuePtr());
     }
   }
 
-    void AddDouble(const char* name, const std::vector<Cut<Double_t>*>& cuts) {
+  void AddDouble(const char *name, const std::vector<Cut<Double_t> *> &cuts)
+  {
     doubles[name] = Variable<Double_t>(name, cuts);
   }
 
-    void AddFloat(const char* name, const std::vector<Cut<Float_t>*>& cuts) {
+  void AddFloat(const char *name, const std::vector<Cut<Float_t> *> &cuts)
+  {
     floats[name] = Variable<Float_t>(name, cuts);
   }
 
-    void AddInt(const char* name, const std::vector<Cut<Int_t>*>& cuts) {
+  void AddInt(const char *name, const std::vector<Cut<Int_t> *> &cuts)
+  {
     ints[name] = Variable<Int_t>(name, cuts);
   }
 
-    void AddBool(const char* name, const std::vector<Cut<Bool_t>*>& cuts) {
+  void AddBool(const char *name, const std::vector<Cut<Bool_t> *> &cuts)
+  {
     bools[name] = Variable<Bool_t>(name, cuts);
   }
 
-    void PrintValues() {
-      for (const auto& [key, var] : doubles) {
+  void AddTriggerCut(const std::vector<std::string> &triggers)
+  {
+    trigger_cuts.push_back(triggers);
+  }
+
+  void AddFourMomentumFor(const char *name)
+  {
+    std::string pxname = TString::Format("%s_PX", name).Data();
+    std::string pyname = TString::Format("%s_PY", name).Data();
+    std::string pzname = TString::Format("%s_PZ", name).Data();
+    std::string ename = TString::Format("%s_ENERGY", name).Data();
+
+    floats[pxname] = Variable<Float_t>(pxname.c_str(), {});
+    floats[pyname] = Variable<Float_t>(pyname.c_str(), {});
+    floats[pzname] = Variable<Float_t>(pzname.c_str(), {});
+    floats[ename] = Variable<Float_t>(ename.c_str(), {});
+  }
+
+  Double_t GetDouble(const char *name) const
+  {
+    auto it = doubles.find(name);
+    if (it != doubles.end())
+    {
+      return it->second.GetValue();
+    }
+    return 0;
+  }
+
+  Float_t GetFloat(const char *name) const
+  {
+    auto it = floats.find(name);
+    if (it != floats.end())
+    {
+      return it->second.GetValue();
+    }
+    return 0;
+  }
+
+  Int_t GetInt(const char *name) const
+  {
+    auto it = ints.find(name);
+    if (it != ints.end())
+    {
+      return it->second.GetValue();
+    }
+    return 0;
+  }
+
+  Bool_t GetBool(const char *name) const
+  {
+    auto it = bools.find(name);
+    if (it != bools.end())
+    {
+      return it->second.GetValue();
+    }
+    return 0;
+  }
+
+  TLorentzVector GetFourMomentum(const char *name)
+  {
+    return ComputeFourMomentum(name, false, 0.);
+  }
+
+  TLorentzVector GetFourMomentumWithMassSub(const char *name, double subbed_mass)
+  {
+    return ComputeFourMomentum(name, true, subbed_mass);
+  }
+
+  bool CheckCuts() const
+  {
+    for (const auto &[key, var] : doubles)
+    {
+      if (!var.CheckCuts())
+      {
+        return false;
+      }
+    }
+
+    for (const auto &[key, var] : floats)
+    {
+      if (!var.CheckCuts())
+      {
+        return false;
+      }
+    }
+
+    for (const auto &[key, var] : ints)
+    {
+      if (!var.CheckCuts())
+      {
+        return false;
+      }
+    }
+
+    for (const auto &[key, var] : bools)
+    {
+      if (!var.CheckCuts())
+      {
+        return false;
+      }
+    }
+
+    for (const auto &triggers : trigger_cuts)
+    {
+      bool cut_okay = false;
+      for (const auto &trigger : triggers)
+      {
+        cut_okay = cut_okay | this->GetBool(trigger.c_str());
+      }
+
+      if (!cut_okay)
+      {
+        return false;
+      }
+    }
+
+    return true;
+  }
+
+  void PrintValues() const
+  {
+    for (const auto &[key, var] : doubles)
+    {
       std::cout << key << ": " << var.GetValue() << std::endl;
     }
 
-      for (const auto& [key, var] : floats) {
+    for (const auto &[key, var] : floats)
+    {
       std::cout << key << ": " << var.GetValue() << std::endl;
     }
 
-      for (const auto& [key, var] : ints) {
+    for (const auto &[key, var] : ints)
+    {
       std::cout << key << ": " << var.GetValue() << std::endl;
     }
 
-      for (const auto& [key, var] : bools) {
+    for (const auto &[key, var] : bools)
+    {
       std::cout << key << ": " << var.GetValue() << std::endl;
     }
   }
 };
 
-int do_all() {
+int do_all()
+{
 
-  VariableCollection col{};
-  col.AddDouble("B_M", {});
-  col.AddDouble("Jpsi_M", {});
-  col.AddFloat("Kplus_PT", {});
-  col.AddBool("muminus_ISMUON", {});
-  col.AddInt("muplus_Q", {});
+  TChain *data_chain = new TChain("BuToHpMuMu/DecayTree");
+  data_chain->Add("/auto/data/pfeiffer/inclusive_detached_dilepton/data_samples/spruce_magdown_2023_v0r1_tuple_90000000_2023_v0r0p6288631.root");
+  data_chain->Add("/auto/data/pfeiffer/inclusive_detached_dilepton/data_samples/spruce_magdown_2023_v0_tuple_90000000_v0r0p6288631.root");
 
-  TChain *data_chain = new TChain(TString::Format("%s/DecayTree", "Hlt2RD_BuToKpMuMu").Data());
-  data_chain->Add("/auto/data/pfeiffer/inclusive_detached_dilepton/data_samples/Collision23_Beam6800GeV-VeloClosed-MagDown-Excl-UT_RealData_SprucingPass23r1_94000000_RD.root");
+  VariableCollection col{};
+  col.AddFloat("B_BPVFDCHI2", {new GT<float>(36.)});
+  col.AddDouble("B_CHI2VXNDOF", {new LT(16.)});
+  col.AddFloat("B_BPVIPCHI2", {new LT<float>(25.)});
+  col.AddDouble("Jpsi_MAXDOCACHI2", {new LT(36.)});
+  col.AddFloat("L1_BPVIPCHI2", {new GT<float>(9.)});
+  col.AddFloat("L2_BPVIPCHI2", {new GT<float>(9.)});
+  col.AddDouble("L1_PID_MU", {new GT(-3.)});
+  col.AddDouble("L2_PID_MU", {new GT(-3.)});
+  col.AddBool("L1_ISMUON", {new EQ(true)});
+  col.AddBool("L2_ISMUON", {new EQ(true)});
+  col.AddFloat("L1_PT", {new GT<float>(350.)});
+  col.AddFloat("L2_PT", {new GT<float>(350.)});
+  col.AddDouble("Jpsi_CHI2DOF", {new LT(9.)});
+  col.AddFloat("Jpsi_BPVFDCHI2", {new GT<float>(16.)});
+  col.AddDouble("Jpsi_M", {new LT(5500.), new BT(3096.9 - 100., 3096.9 + 100.)});
+  col.AddFloat("Hp_PT", {new GT<float>(400.)});
+  col.AddFloat("Hp_BPVIPCHI2", {new GT<float>(6.)});
+  col.AddFloat("Hp_P", {new GT<float>(2000.)});
+  col.AddDouble("Hp_PID_K", {new GT(-4.)});
+  col.AddBool("Hlt2_InclDetDiMuon_4BodyDecision", {});
+  col.AddBool("Hlt2_InclDetDiMuon_3BodyDecision", {});
+  col.AddBool("Hlt2_InclDetDiMuonDecision", {});
+  col.AddBool("Hlt2RD_BuToKpMuMuDecision", {});
+  col.AddBool("Hlt1TrackMVADecision", {});
+  col.AddBool("Hlt1TwoTrackMVADecision", {});
+  col.AddFourMomentumFor("L1");
+  col.AddFourMomentumFor("L2");
+  col.AddFourMomentumFor("Hp");
+
+  col.AddTriggerCut({"Hlt2_InclDetDiMuonDecision", "Hlt2_InclDetDiMuon_3BodyDecision", "Hlt2_InclDetDiMuon_4BodyDecision"});
+  col.AddTriggerCut({"Hlt1TrackMVADecision", "Hlt1TwoTrackMVADecision"});
+  col.AddTriggerCut({"Hlt2RD_BuToKpMuMuDecision"});
 
   col.Connect(data_chain);
 
-  data_chain->GetEntry(1);
+  const Double_t MASS_HIST_MIN = 5100.;
+  const Double_t MASS_HIST_MAX = 6000.;
+  const int N_BINS = 70;
+  TH1D *h1_B_M = new TH1D("h1_B_M", "B Mass", N_BINS, MASS_HIST_MIN, MASS_HIST_MAX);
+
+  unsigned int entries = data_chain->GetEntries();
+  for (unsigned int i = 0; i < entries; i++)
+  {
+    data_chain->GetEntry(i);
 
-  col.PrintValues();
+    if (col.CheckCuts())
+    {
+      auto l1_4v = col.GetFourMomentum("L1");
+      auto l2_4v = col.GetFourMomentum("L2");
+      auto K_4v = col.GetFourMomentumWithMassSub("Hp", 493.677);
 
-  data_chain->GetEntry(2);
+      h1_B_M->Fill((l1_4v + l2_4v + K_4v).M());
+    }
+
+    if ((i + 1) % 10000 == 0 || i + 1 == entries)
+    {
+      std::cout << "["
+                << "BuToHpMuMu"
+                << "] Processed event: " << i + 1 << " (" << TString::Format("%.2f", ((double)(i + 1) / (double)entries) * 100.) << "%)" << std::endl;
+    }
+  }
 
-  col.PrintValues();
+  TCanvas *c1 = new TCanvas("c1", "c1", 0, 0, 800, 600);
+  h1_B_M->Draw();
+  c1->Draw();
 
   return 0;
 }

analysis_fullstream_b02hphmmumu.cpp → old_analysis_macros/analysis_fullstream_b02hphmmumu.cpp

@@ -40,8 +40,8 @@ const Double_t MASS_HIST_MAX = 5700.;
 const Double_t MASS_HIST_FIT_MIN = 5100.;
 const Double_t MASS_HIST_FIT_MAX = 5700.;
 
-const char* TITLE = "SpruceRD_B0ToHpHmMuMu (#pi^{+} #rightarrow K^{+})";
-const char* FILE_NAME = "SpruceRD_B0ToHpHmMuMu_Pip2Kp";
+const char* TITLE = "SpruceRD_B0ToHpHmMuMu (#pi^{+} #rightarrow K^{+}, w/o excl Hlt2 decision cut)";
+const char* FILE_NAME = "SpruceRD_B0ToHpHmMuMu_Pip2Kp_wohlt2cut";
 const char* MASS_LITERAL = "m(#pi^{+}_{(#rightarrow K^{+})}#pi^{-}#mu^{+}#mu^{-})";
 
 const Double_t K_MASS = 493.677;
@@ -203,7 +203,7 @@ int analysis_fullstream_b02hphmmumu()
                 (Hp_BPVIPCHI2 > 6) & (Hp_PT > 250) & (Hp_P > 2000) & (Hp_PID_K > 0)
             ) &
             //
-            (Hlt2RD_B0ToKpPimMuMuDecision) & ((Hlt2_InclDetDiMuon_4BodyDecision) | (Hlt2_InclDetDiMuon_3BodyDecision) | (Hlt2_InclDetDiMuonDecision)))
+            /*(Hlt2RD_B0ToKpPimMuMuDecision) &*/ ((Hlt2_InclDetDiMuon_4BodyDecision) | (Hlt2_InclDetDiMuon_3BodyDecision) | (Hlt2_InclDetDiMuonDecision)))
         {
             if ((TMath::Abs(Jpsi_M - 3096.9) < 100) & ((Hlt1TrackMVADecision) | (Hlt1TwoTrackMVADecision)))
             {

analysis_fullstream_bu2hpmumu.cpp → old_analysis_macros/analysis_fullstream_bu2hpmumu.cpp

@@ -43,8 +43,8 @@ const Double_t K_MASS = 493.677;
 
 void CreateRooFitAndDraw(TH1D *hist, int fitting_entries);
 
-const char *TITLE = "SpruceRD_BuToHpMuMu (#pi^{+} #rightarrow K^{+}, Small Cuts 2)";
-const char *FILE_NAME = "SpruceRD_BuToHpMuMu_Pip2Kp_small2";
+const char *TITLE = "SpruceRD_BuToHpMuMu (#pi^{+} #rightarrow K^{+}, w/o excl Hlt2 decision cut)";
+const char *FILE_NAME = "SpruceRD_BuToHpMuMu_Pip2Kp_wohlt2cut";
 const char *MASS_LITERAL = "m(#pi^{+}_{(#rightarrow K^{+})}#mu^{+}#mu^{-})";
 
 int analysis_fullstream_bu2hpmumu()
@@ -145,7 +145,8 @@ int analysis_fullstream_bu2hpmumu()
         TLorentzVector l2_4v(L2_PX, L2_PY, L2_PZ, L2_ENERGY);
         Double_t reconstructed_B_Mass = (K_4v + l1_4v + l2_4v).M();
 
-
+        if ((((B_BPVFDCHI2 > 36) & (B_CHI2VXNDOF < 16) & (B_BPVIPCHI2 < 25) & (Jpsi_MAXDOCACHI2 < 36) & (L1_BPVIPCHI2 > 9) & (L2_BPVIPCHI2 > 9) & (L1_PID_MU > -3) & (L2_PID_MU > -3) & (L1_ISMUON == 1) & (L2_ISMUON == 1) & (L1_PT > 350) & (L2_PT > 350) & (Jpsi_CHI2DOF < 9) & (Jpsi_BPVFDCHI2 > 16) & (Jpsi_M < 5500) & (Hp_PT > 400) & (Hp_BPVIPCHI2 > 6) & (Hp_PT > 400) & (Hp_BPVIPCHI2 > 6) & (Hp_P > 2000) & (Hp_PID_K > -4)) /*& (Hlt2RD_BuToKpMuMuDecision == 1)*/ & (((Hlt2_InclDetDiMuon_4BodyDecision == 1) | (Hlt2_InclDetDiMuon_3BodyDecision == 1) | (Hlt2_InclDetDiMuonDecision == 1)))))
+        {
             if ((TMath::Abs(Jpsi_M - 3096.9) < 100) & ((Hlt1TrackMVADecision == 1) | (Hlt1TwoTrackMVADecision == 1)))
             {
                 h1_B_M->Fill(reconstructed_B_Mass);
@@ -155,11 +156,7 @@ int analysis_fullstream_bu2hpmumu()
                     fitting_entries++;
                 }
             }
-
-        // if ((((B_BPVFDCHI2 > 36) & (B_CHI2VXNDOF < 16) & (B_BPVIPCHI2 < 25) & (Jpsi_MAXDOCACHI2 < 36) & (L1_BPVIPCHI2 > 9) & (L2_BPVIPCHI2 > 9) & (L1_PID_MU > -3) & (L2_PID_MU > -3) & (L1_ISMUON == 1) & (L2_ISMUON == 1) & (L1_PT > 350) & (L2_PT > 350) & (Jpsi_CHI2DOF < 9) & (Jpsi_BPVFDCHI2 > 16) & (Jpsi_M < 5500) & (Hp_PT > 400) & (Hp_BPVIPCHI2 > 6) & (Hp_PT > 400) & (Hp_BPVIPCHI2 > 6) & (Hp_P > 2000) & (Hp_PID_K > -4)) & (Hlt2RD_BuToKpMuMuDecision == 1) & (((Hlt2_InclDetDiMuon_4BodyDecision == 1) | (Hlt2_InclDetDiMuon_3BodyDecision == 1) | (Hlt2_InclDetDiMuonDecision == 1)))))
-        // {
-            
-        // }
+        }
 
         if (i % 10000 == 0)
         {
@@ -190,7 +187,6 @@ void CreateRooFitAndDraw(TH1D *hist, int fitting_entries)
     roohist_B_M.plotOn(roo_frame_mass, RooFit::Binning(nBins), RooFit::Name("B Mass Distribution"));
     roo_frame_mass->GetXaxis()->SetTitle(MASS_LITERAL);
 
-    
     // Crystal Ball for Signal
     RooRealVar var_mass_x0("var_mass_x0", "#mu", 5278., 5170., 5500.);
     RooRealVar var_mass_sigmaLR("var_mass_sigmaLR", "#sigma_{LR}", 16., 5., 40.);
@@ -202,10 +198,10 @@ void CreateRooFitAndDraw(TH1D *hist, int fitting_entries)
     // RooRealVar var_mass_alphaR("var_mass_alphaR", "#alpha_{R}", 2., 0., 4.);
     // RooRealVar var_mass_nR("var_mass_nR", "n_{R}", 5., 0., 15.);
 
-    auto var_mass_alphaL = RooRealConstant::value(1.309);
-    auto var_mass_nL = RooRealConstant::value(0.72);
-    auto var_mass_alphaR = RooRealConstant::value(1.226);
-    auto var_mass_nR = RooRealConstant::value(0.992);
+    auto var_mass_alphaL = RooRealConstant::value(1.894);
+    auto var_mass_nL = RooRealConstant::value(1.120);
+    auto var_mass_alphaR = RooRealConstant::value(2.446);
+    auto var_mass_nR = RooRealConstant::value(1.146);
 
     RooCrystalBall sig_cb("sig_cb", "Signal Crystal Ball", roo_var_mass, var_mass_x0, var_mass_sigmaLR, var_mass_alphaL, var_mass_nL, var_mass_alphaR, var_mass_nR);
 
@@ -221,7 +217,7 @@ void CreateRooFitAndDraw(TH1D *hist, int fitting_entries)
 
     RooFitResult *fitres = sigplusbkg.fitTo(roohist_B_M, RooFit::Save(), RooFit::PrintLevel(1), RooFit::Range("fitting_range"));
 
-    sigplusbkg.plotOn(roo_frame_mass, RooFit::VisualizeError(*fitres, 1), RooFit::FillColor(kOrange + 1), RooFit::FillStyle(3144));
+    // sigplusbkg.plotOn(roo_frame_mass, RooFit::VisualizeError(*fitres, 1), RooFit::FillColor(kOrange + 1), RooFit::FillStyle(3144));
     sigplusbkg.plotOn(roo_frame_mass, RooFit::LineColor(kRed), RooFit::LineStyle(kSolid), RooFit::Range("fitting_range"), RooFit::Name("sigplusbkg"));
     sigplusbkg.plotOn(roo_frame_mass, RooFit::Components(RooArgSet(bkg_exp)), RooFit::LineColor(kBlue - 7), RooFit::LineStyle(kDashed), RooFit::Range("fitting_range"), RooFit::Name("bkg_exp"));
     sigplusbkg.plotOn(roo_frame_mass, RooFit::Components(RooArgSet(sig_cb)), RooFit::LineColor(kRed - 7), RooFit::LineStyle(kDashed), RooFit::Range("fitting_range"), RooFit::Name("sig_cb"));

analysis_turbo.cpp → old_analysis_macros/analysis_turbo.cpp

@@ -34,22 +34,30 @@
 
 const int nBins = 70;
 
-const Double_t MASS_HIST_MIN = 5100.;
+const Double_t MASS_HIST_MIN = 5150.;
 const Double_t MASS_HIST_MAX = 6000.;
 
-const char *B0_DECAY = "Hlt2RD_B0ToKpPimMuMu";
-const char *Bu_DECAY = "Hlt2RD_BuToKpMuMu";
-
-const char *TITLE = B0_DECAY;
-const char *X_AXIS = "m(K^{+}#pi^{-}#mu^{+}#mu^{-})";
-
+const char *B0_DECAY = "Hlt2RD_B0ToKpPimMuMu_woHtl1cut";
+const char *Bu_DECAY = "Hlt2RD_BuToKpMuMu_woHtl1cut";
 
+const char *TITLE = Bu_DECAY;
+const char *HIST_TITLE = "Hlt2RD_BuToKpMuMu w/o Hlt1 Decision Cut";
+const char *X_AXIS = "m(K^{+}#mu^{+}#mu^{-})";
 
 void CreateRooFitAndDraw(TH1D *hist, int fitting_entries);
 
 int analysis_turbo()
 {
-    TChain *data_chain = new TChain(TString::Format("%s/DecayTree", TITLE).Data());
+    std::string tree = "";
+    if (TITLE == B0_DECAY)
+    {
+        tree = "Hlt2RD_B0ToKpPimMuMu";
+    }
+    else
+    {
+        tree = "Hlt2RD_BuToKpMuMu";
+    }
+    TChain *data_chain = new TChain(TString::Format("%s/DecayTree", tree.c_str()).Data());
     data_chain->Add("/auto/data/pfeiffer/inclusive_detached_dilepton/data_samples/Collision23_Beam6800GeV-VeloClosed-MagDown-Excl-UT_RealData_SprucingPass23r1_94000000_RD.root");
 
     Double_t B_M, Jpsi_M, B0_M, Kst_M;
@@ -72,7 +80,7 @@ int analysis_turbo()
         data_chain->SetBranchAddress("Hlt1TwoTrackMVADecision", &Hlt1TwoTrackMVADecision);
     }
 
-    TH1D *h1_B_M = new TH1D("h1_B_M", TITLE, nBins, MASS_HIST_MIN, MASS_HIST_MAX);
+    TH1D *h1_B_M = new TH1D("h1_B_M", HIST_TITLE, nBins, MASS_HIST_MIN, MASS_HIST_MAX);
     h1_B_M->GetXaxis()->SetTitle(X_AXIS);
 
     int fitting_entries = 0;
@@ -83,7 +91,7 @@ int analysis_turbo()
 
         if (TITLE == B0_DECAY)
         {
-            if ((TMath::Abs(Jpsi_M - 3096.9) < 100) && (B0_M > 4500) && (B0_M < 6000) && (TMath::Abs(Kst_M - 895.55) < 50) && ((Hlt1TrackMVADecision) || (Hlt1TwoTrackMVADecision)))
+            if ((TMath::Abs(Jpsi_M - 3096.9) < 100) && (B0_M > 4500) && (B0_M < 6000) && (TMath::Abs(Kst_M - 895.55) < 50) /*&& ((Hlt1TrackMVADecision) || (Hlt1TwoTrackMVADecision))*/)
             {
                 h1_B_M->Fill(B0_M);
                 if (MASS_HIST_MIN <= B0_M && B0_M <= MASS_HIST_MAX)
@@ -94,7 +102,7 @@ int analysis_turbo()
         }
         else
         {
-            if ((TMath::Abs(Jpsi_M - 3096.9) < 100.) && (B_M > 4500.) && (B_M < 6000.) && ((Hlt1TrackMVADecision) || (Hlt1TwoTrackMVADecision)))
+            if ((TMath::Abs(Jpsi_M - 3096.9) < 100.) && (B_M > 4500.) && (B_M < 6000.) /*&& ((Hlt1TrackMVADecision) || (Hlt1TwoTrackMVADecision))*/)
             {
                 h1_B_M->Fill(B_M);
                 if (MASS_HIST_MIN <= B_M && B_M <= MASS_HIST_MAX)
@@ -133,7 +141,7 @@ void CreateRooFitAndDraw(TH1D *hist, int fitting_entries)
 
     RooDataHist roohist_B_M("roohist_B_M", "B Mass Histogram", roo_var_mass, RooFit::Import(*hist));
 
-    RooPlot *roo_frame_mass = roo_var_mass.frame(RooFit::Title(TITLE));
+    RooPlot *roo_frame_mass = roo_var_mass.frame(RooFit::Title(HIST_TITLE));
     roohist_B_M.plotOn(roo_frame_mass, RooFit::Binning(nBins), RooFit::Name("B Mass Distribution"));
     roo_frame_mass->GetXaxis()->SetTitle(X_AXIS);
 
@@ -148,10 +156,17 @@ void CreateRooFitAndDraw(TH1D *hist, int fitting_entries)
     // RooRealVar var_mass_alphaR("var_mass_alphaR", "#alpha_{R}", 2., 0., 4.);
     // RooRealVar var_mass_nR("var_mass_nR", "n_{R}", 5., 0., 15.);
 
-    auto var_mass_alphaL = RooRealConstant::value(1.948);
-    auto var_mass_nL = RooRealConstant::value(0.618);
-    auto var_mass_alphaR = RooRealConstant::value(2.320);
-    auto var_mass_nR = RooRealConstant::value(0.473);
+    // B+/-
+    auto var_mass_alphaL = RooRealConstant::value(1.912);
+    auto var_mass_nL = RooRealConstant::value(1.125);
+    auto var_mass_alphaR = RooRealConstant::value(2.447);
+    auto var_mass_nR = RooRealConstant::value(1.491);
+
+    // B0
+    // auto var_mass_alphaL = RooRealConstant::value(1.948);
+    // auto var_mass_nL = RooRealConstant::value(0.618);
+    // auto var_mass_alphaR = RooRealConstant::value(2.320);
+    // auto var_mass_nR = RooRealConstant::value(0.473);
 
     RooCrystalBall sig_cb("sig_cb", "Signal Crystal Ball", roo_var_mass, var_mass_x0, var_mass_sigmaLR, var_mass_alphaL, var_mass_nL, var_mass_alphaR, var_mass_nR);
 
@@ -167,7 +182,7 @@ void CreateRooFitAndDraw(TH1D *hist, int fitting_entries)
 
     RooFitResult *fitres = sigplusbkg.fitTo(roohist_B_M, RooFit::Save(), RooFit::PrintLevel(1), RooFit::Range("fitting_range"));
 
-    sigplusbkg.plotOn(roo_frame_mass, RooFit::VisualizeError(*fitres, 1), RooFit::FillColor(kOrange + 1), RooFit::FillStyle(3144));
+    sigplusbkg.plotOn(roo_frame_mass, RooFit::VisualizeError(*fitres, 1), RooFit::FillColor(kYellow -7), RooFit::FillStyle(3144));
     sigplusbkg.plotOn(roo_frame_mass, RooFit::LineColor(kRed), RooFit::LineStyle(kSolid), RooFit::Range("fitting_range"), RooFit::Name("sigplusbkg"));
     sigplusbkg.plotOn(roo_frame_mass, RooFit::Components(RooArgSet(bkg_exp)), RooFit::LineColor(kBlue - 7), RooFit::LineStyle(kDashed), RooFit::Range("fitting_range"), RooFit::Name("bkg_exp"));
     sigplusbkg.plotOn(roo_frame_mass, RooFit::Components(RooArgSet(sig_cb)), RooFit::LineColor(kRed - 7), RooFit::LineStyle(kDashed), RooFit::Range("fitting_range"), RooFit::Name("sig_cb"));

simulation_fullstream_B02KpPimMuMu.cpp → old_analysis_macros/simulation_fullstream_B02KpPimMuMu.cpp

@@ -44,7 +44,7 @@ const int PID_KAON = 321;
 const int PID_PION = 211;
 const int PID_MUON = 13;
 
-const char *NAME = "B0ToHpHmMuMu_Fullstream";
+const char *NAME = "B0ToHpHmMuMu_Fullstream_wobkgcat";
 const char *X_AXIS = "m(#pi^{+}_{#rightarrow K^{+}}#pi^{-}#mu^{+}#mu^{-})";
 const bool USE_HYP_REPLACE = false;
 const char *B_NAME = "B0";
@@ -120,7 +120,7 @@ int simulation_fullstream_B02KpPimMuMu()
 
         Double_t reconstructed_B_Mass = 0;
 
-        if (TMath::Abs(L1_TRUEID) == PID_MUON && L2_TRUEID == -L1_TRUEID)
+        if (B_BKGCAT == 30 && TMath::Abs(L1_TRUEID) == PID_MUON && L2_TRUEID == -L1_TRUEID)
         {
             if (TMath::Abs(Hp_TRUEID) == PID_KAON && TMath::Abs(Hm_TRUEID) == PID_PION)
             {