From 1a252280445c6f4dfb35510024e5d556466b92a6 Mon Sep 17 00:00:00 2001
From: Ziyue Zhang <zzhan70@uic.edu>
Date: Tue, 23 Mar 2021 20:56:01 +0000
Subject: [PATCH] Dvmp jpsi mass resolution
---
benchmarks/dvmp/analysis/vm_mass.cxx | 48 +++++++++++++++++-----------
include/util.h | 2 ++
options/tracker_reconstruction.py | 2 +-
3 files changed, 33 insertions(+), 19 deletions(-)
diff --git a/benchmarks/dvmp/analysis/vm_mass.cxx b/benchmarks/dvmp/analysis/vm_mass.cxx
index 9d9b32b1..577e4a17 100644
--- a/benchmarks/dvmp/analysis/vm_mass.cxx
+++ b/benchmarks/dvmp/analysis/vm_mass.cxx
@@ -68,8 +68,8 @@ int vm_mass(const std::string& config_name)
{"description", "Invariant Mass Resolution calculated from raw "
"tracking data using a Gaussian fit."},
{"quantity", "resolution"},
- {"target", ".1"}}};
-
+ {"target", ".2"}}}; //these 2 need to be consistent
+ double width_target = 0.2; //going to find a way to use the same variable
// Run this in multi-threaded mode if desired
ROOT::EnableImplicitMT(kNumThreads);
@@ -112,8 +112,8 @@ int vm_mass(const std::string& config_name)
// Define output histograms
//auto h_im_rec = d_im.Histo1D({"h_im_rec", ";m_{ll'} (GeV/c^{2});#", (int)(vm_mass+0.5)*2*100, 0., 2.*(int)(vm_mass+0.5)}, "mass_rec"); //real rec
- auto h_im_rec = d_im.Histo1D({"h_im_rec", ";m_{ll'} (GeV/c^{2});#", 50, 1., 5.}, "mass_rec");//for dummy_rec
- auto h_im_sim = d_im.Histo1D({"h_im_sim", ";m_{ll'} (GeV/c^{2});#", 50, 1., 5.}, "mass_sim");
+ auto h_im_rec = d_im.Histo1D({"h_im_rec", ";m_{ll'} (GeV/c^{2});#", 30, 1.5, 4.5}, "mass_rec");//for dummy_rec
+ auto h_im_sim = d_im.Histo1D({"h_im_sim", ";m_{ll'} (GeV/c^{2});#", 30, 1.5, 4.5}, "mass_sim");
auto h_pt_rec = d_im.Histo1D({"h_pt_rec", ";p_{T} (GeV/c);#", 50, 0., 10.}, "pt_rec");
auto h_pt_sim = d_im.Histo1D({"h_pt_sim", ";p_{T} (GeV/c);#", 50, 0., 10.}, "pt_sim");
@@ -127,7 +127,7 @@ int vm_mass(const std::string& config_name)
// Plot our histograms.
// TODO: to start I'm explicitly plotting the histograms, but want to
// factorize out the plotting code moving forward.
- {
+ //{
TCanvas c{"canvas", "canvas", 1200, 1200};
c.Divide(2, 2, 0.0001, 0.0001);
// pad 1 mass
@@ -147,8 +147,22 @@ int vm_mass(const std::string& config_name)
// draw everything
h11.DrawClone("hist");
h12.DrawClone("hist same");
+
+ //Fit
+ TF1* mfMass = new TF1("mfMass", "[2]*TMath::Gaus(x, [0], [1], kFALSE)", 1.5, 4.5);
+ mfMass->SetParameters(3.096, 0.1, 100.);
+ mfMass->SetParLimits(0, 3.0, 3.2);
+ mfMass->SetParLimits(1, 0., 10.);
+ mfMass->SetParLimits(2, 0., 1000.);
+ mfMass->SetNpx(1000);
+ mfMass->SetLineColor(2);
+ mfMass->SetLineStyle(7);
+
+ TFitResultPtr myFitPtr = h12.Fit(mfMass, "S 0", "", 1.5, 4.5);
+ mfMass->Draw("same");
+
// FIXME hardcoded beam configuration
- plot::draw_label(10, 100, detector);
+ plot::draw_label(10, 100, detector);
TText* tptr1;
auto t1 = new TPaveText(.6, .8417, .9, .925, "NB NDC");
t1->SetFillColorAlpha(kWhite, 0);
@@ -175,18 +189,8 @@ int vm_mass(const std::string& config_name)
h21.GetXaxis()->CenterTitle();
h21.GetYaxis()->CenterTitle();
// draw everything
- //TF1* mfPt = new TF1("mfPt", "[0]*exp(-[1]*x)", 1.2, 10.);
- //mfPt->SetParameters(5., 1.);
- //mfPt->SetParLimits(0, 0., 1000000.);
- //mfPt->SetParLimits(1, 0., 1000000.);
- //mfPt->SetNpx(1000);
- //mfPt->SetLineColor(2);
- //mfPt->SetLineStyle(7);
- //h21.Fit(mfPt, "", "", 1.2, 10.);
h21.DrawClone("hist");
h22.DrawClone("hist same");
- //mfPt->Draw("same");
-
// FIXME hardcoded beam configuration
plot::draw_label(10, 100, detector);
@@ -262,11 +266,19 @@ int vm_mass(const std::string& config_name)
t4->Draw();
c.Print(fmt::format("{}vm_mass_pt_phi_rapidity.png", output_prefix).c_str());
- }
+ //}
// TODO we're not actually doing an IM fit yet, so for now just return an
// error for the test result
- mass_resolution_test.error(-1);
+ double width = mfMass->GetParameter(1);
+ if(myFitPtr->Status()!=0){
+ mass_resolution_test.error(-1);
+ }else if(width > width_target){
+ mass_resolution_test.fail(width);
+ }else{
+ mass_resolution_test.pass(width);
+ }
+
// write out our test data
eic::util::write_test(mass_resolution_test, fmt::format("{}mass.json", output_prefix));
diff --git a/include/util.h b/include/util.h
index eddf1393..8083d38c 100644
--- a/include/util.h
+++ b/include/util.h
@@ -11,6 +11,8 @@
#include <string>
#include <vector>
#include "TF1.h"
+#include "TFitResult.h"
+#include "TFitResultPtr.h"
#include <Math/Vector4D.h>
diff --git a/options/tracker_reconstruction.py b/options/tracker_reconstruction.py
index 6bef5d7b..e303b38e 100644
--- a/options/tracker_reconstruction.py
+++ b/options/tracker_reconstruction.py
@@ -57,7 +57,7 @@ podioinput = PodioInput("PodioReader",
dummy = MC2DummyParticle("MC2Dummy",
inputCollection="mcparticles",
outputCollection="DummyReconstructedParticles",
- smearing = 0.0)
+ smearing = 0.1)
## copiers to get around input --> output copy bug. Note the "2" appended to the output collection.
copier = MCCopier("MCCopier",
--
GitLab