Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
plot_kshort.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file plot_kshort.C
1 #include <cmath>
2 #include <TFile.h>
3 #include <TNtuple.h>
4 #include <TH2D.h>
5 #include <TCut.h>
6 #include <Eigen/Dense>
7 
8 
9 void plot_kshort(std::string inputFile = "./kshort_G4sPHENIX.root", const std::string outputFile = "./kshort_plots.root")
10 {
11 
12 
13  bool plotBool = true; // Determines whether plots are displayed or not
14  bool savePlot = true; // Determines whether histograms are saved to output file or not
15  bool localVerbosity = true;
16  bool analyzeMassSpectrum = true;
17 
18  TFile fin(inputFile.c_str());
19 
20  TNtuple *ntuple;
21  fin.GetObject("ntp_reco_info",ntuple);
22 
23 
24  /* Below are values in massRecoAnalysis ntuple 3/6/23
25 "ntp_reco_info","decay_pairs","x1:y1:z1:px1:py1:pz1:dca3dxy1:dca3dz1:phi1:pca_rel1_x:pca_rel1_y:pca_rel1_z:eta1:charge1:tpcClusters_1:x2:y2:z2:px2:py2:pz2:dca3dxy2:dca3dz2:phi2:pca_rel2_x:pca_rel2_y:pca_rel2_z:eta2:charge2:tpcClusters_2:vertex_x:vertex_y:vertex_z:pair_dca:invariant_mass:invariant_pt:pathlength_x:pathlength_y:pathlength_z:pathlength:rapidity:pseudorapidity:projected_pos1_x:projected_pos1_y:projected_pos1_z:projected_pos2_x:projected_pos2_y:projected_pos2_z:projected_mom1_x:projected_mom1_y:projected_mom1_z:projected_mom2_x:projected_mom2_y:projected_mom2_z:projected_pca_rel1_x:projected_pca_rel1_y:projected_pca_rel1_z:projected_pca_rel2_x:projected_pca_rel2_y:projected_pca_rel2_z:projected_pair_dca:projected_pathlength_x:projected_pathlength_y:projected_pathlength_z:projected_pathlength:quality1:quality2:cosThetaReco"
26  */
27 
28 
29  float invariant_mass;
30  float projected_pathlength_x;
31  float projected_pathlength_y;
32  float projected_pathlength_z;
33  float projected_pathlength;
34  float vertex_x;
35  float vertex_y;
36  float vertex_z;
37  float projected_mom1_x;
38  float projected_mom1_y;
39  float projected_mom1_z;
40  float projected_mom2_x;
41  float projected_mom2_y;
42  float projected_mom2_z;
43  float projected_pair_dca;
44  float invariant_pt;
45  float quality1;
46  float quality2;
47  float dca3dxy1;
48  float dca3dxy2;
49  float dca3dz1;
50  float dca3dz2;
51 
52 
53  ntuple->SetBranchAddress("invariant_mass",&invariant_mass);
54  ntuple->SetBranchAddress("projected_pathlength_x",&projected_pathlength_x);
55  ntuple->SetBranchAddress("projected_pathlength_y",&projected_pathlength_y);
56  ntuple->SetBranchAddress("projected_pathlength_z",&projected_pathlength_z);
57  ntuple->SetBranchAddress("projected_pathlength",&projected_pathlength);
58  ntuple->SetBranchAddress("vertex_x",&vertex_x);
59  ntuple->SetBranchAddress("vertex_y",&vertex_y);
60  ntuple->SetBranchAddress("vertex_z",&vertex_z);
61  ntuple->SetBranchAddress("projected_mom1_x",&projected_mom1_x);
62  ntuple->SetBranchAddress("projected_mom1_y",&projected_mom1_y);
63  ntuple->SetBranchAddress("projected_mom1_z",&projected_mom1_z);
64  ntuple->SetBranchAddress("projected_mom2_x",&projected_mom2_x);
65  ntuple->SetBranchAddress("projected_mom2_y",&projected_mom2_y);
66  ntuple->SetBranchAddress("projected_mom2_z",&projected_mom2_z);
67  ntuple->SetBranchAddress("projected_pair_dca",&projected_pair_dca);
68  ntuple->SetBranchAddress("invariant_pt",&invariant_pt);
69  ntuple->SetBranchAddress("quality1",&quality1);
70  ntuple->SetBranchAddress("quality2",&quality2);
71  ntuple->SetBranchAddress("dca3dxy1",&dca3dxy1);
72  ntuple->SetBranchAddress("dca3dxy2",&dca3dxy2);
73  ntuple->SetBranchAddress("dca3dz1",&dca3dz1);
74  ntuple->SetBranchAddress("dca3dz2",&dca3dz2);
75 
76 
77  int entries = ntuple->GetEntries();
78 
79  TH1D *invariant_mass_hist = new TH1D("invariant_mass_hist","Invariant Mass",1000,0.35,0.65);
80  invariant_mass_hist->SetMinimum(0);
81  TH2D *pathlengthpathlengthz = new TH2D("pathlengthpathlengthz","path v. pathlength_z",5000,-20.0,20.0,5000,0.0,10.0);
82  TH2D *pathMass = new TH2D("pathMass","Invariant Mass vs. Path Length",5000,0.0,1.0,5000,0.0,10.0);
83 
84 
85  int qual_cut = 10;
86  float pathlength_cut = 0.09;
87  double track_dca_cut = 0.0105;
88 
89  for(int i=0; i<entries; ++i)
90  {
91  ntuple->GetEntry(i);
92 
93  Eigen::Vector3f mom1(projected_mom1_x,projected_mom1_y,projected_mom1_z);
94  Eigen::Vector3f mom2(projected_mom2_x,projected_mom2_y,projected_mom2_z);
95  Eigen::Vector3f projected_momentum = mom1 + mom2;
96  Eigen::Vector3f pathLength (projected_pathlength_x,projected_pathlength_y,projected_pathlength_z);
97  float pathLengthMag = pathLength.norm();
98 
99  Eigen::Vector3f vertex(vertex_x, vertex_y, vertex_z);
100  float costheta = pathLength.dot(projected_momentum)/(projected_momentum.norm()*pathLength.norm());
101  //float radius = sqrt(pow(projected_pathlength_x,2) + pow(projected_pathlength_y,2))
102  /*
103  //Aligned optimized cuts
104  if(costheta < 0.9995) continue;
105  //if(abs(projected_pathlength_x) < 0.1 || abs(projected_pathlength_y) < 0.1) continue;
106  if(abs(projected_pathlength_x) < pathlength_cut || abs(projected_pathlength_y) < pathlength_cut) continue;
107  if(abs(projected_pair_dca) > 0.035) continue;
108  if(invariant_pt < 0.1) continue;
109  if(quality1 > qual_cut || quality2 > qual_cut) continue;
110  if(dca3dxy1 < track_dca_cut || dca3dxy2 < track_dca_cut || dca3dz1 < track_dca_cut || dca3dz2 < track_dca_cut) continue;
111  */
112  if(costheta < 0.99) continue;
113  if(abs(projected_pathlength_x) < pathlength_cut || abs(projected_pathlength_y) < pathlength_cut) continue;
114  if(abs(projected_pair_dca) > 0.055) continue;
115  if(invariant_pt < 0.1) continue;
116  if(quality1 > qual_cut || quality2 > qual_cut) continue;
117  if(dca3dxy1 < track_dca_cut || dca3dxy2 < track_dca_cut || dca3dz1 < track_dca_cut || dca3dz2 < track_dca_cut) continue;
118 
119  invariant_mass_hist->Fill(invariant_mass);
120  pathlengthpathlengthz->Fill(projected_pathlength_z,pathLengthMag);
121  pathMass->Fill(invariant_mass,pathLengthMag);
122  }
123 
124  TF1* f = new TF1("f","gaus(0)+[3]+[4]*x",0.45,0.55);
125 
126  if(plotBool)
127  {
128  TCanvas *c1 = new TCanvas("c1","",10,10,800,800);
129 
130  f->SetParameter(0,20); //100 for aligned data
131  f->SetParameter(1,0.497);
132  f->SetParameter(2,0.020); //0.010 for aligned data
133  f->SetParameter(3,0.0);
134  f->SetParameter(4,0.0);
135 
136  invariant_mass_hist->Fit("f","R L");//R intially for perfect geometry
137  invariant_mass_hist->DrawCopy();
138  }
139 
140 
141  int binLow = invariant_mass_hist->FindBin(0.48);
142  int binHigh = invariant_mass_hist->FindBin(0.52);
143  double integral = invariant_mass_hist->Integral(binLow,binHigh);
144 
145  TF1* fbackground = new TF1("fbackground","[0]+[1]*x");
146  fbackground->SetParameter(0,f->GetParameter(3));
147  fbackground->SetParameter(1,f->GetParameter(4));
148 
149  double backgroundIntegral = fbackground->Integral(0.485,0.51);
150  double foregroundIntegral = f->Integral(0.485,0.51);
151  double signal = (foregroundIntegral - backgroundIntegral);
152  double signalBackgroundRatio = signal/backgroundIntegral;
153 
154 
155  TF1* fsignal = new TF1("fsignal","gaus(0)");
156  fsignal->SetParameter(0,f->GetParameter(0));
157  fsignal->SetParameter(1,f->GetParameter(1));
158  fsignal->SetParameter(2,f->GetParameter(2));
159  double signal_integral = fsignal->Integral(0.35,0.65);
160  double binwidth = invariant_mass_hist->GetBinWidth(1);
161 
162  if(plotBool)
163  {
164  TCanvas *c2 = new TCanvas("c2","",10,10,600,600); // path v path_z
165  pathlengthpathlengthz->DrawCopy();
166 
167  TCanvas *c3 = new TCanvas("c3","",10,10,600,600);
168  pathMass->DrawCopy();
169  }
170 
171  if(localVerbosity)
172  {
173  std::cout << " Integral: "<< integral << std::endl;
174  std::cout<< "signal: " << signal<<" Signal to background ratio: " << signalBackgroundRatio << std::endl;
175  std::cout << " Signal Integral: "<< signal_integral<< " yield: "<< signal_integral/binwidth<<std::endl;
176  }
177 
178  if(savePlot)
179  {
180  std::unique_ptr<TFile> kShortFile(TFile::Open(outputFile.c_str(), "recreate"));
181  kShortFile->WriteObject(invariant_mass_hist,"invariant_mass");
182  kShortFile->WriteObject(pathlengthpathlengthz,"pathlength_v_pathlengthz");
183  kShortFile->WriteObject(pathMass,"pathlength_v_invariantMass");
184  }
185 
186 }