Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
plot_residuals.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file plot_residuals.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 void plot_residuals(const std::string inputFile = "residuals_G4sPHENIX.root", const std::string outputFile = "residuals_plots.root")
9 {
10  bool plotBool = false; // Determines whether plots are displayed or not
11  bool savePlot = true; // Determines whether histograms are saved to output file or not
12 
13  TFile fin(inputFile.c_str());
14 
15  TNtuple *ntuple;
16  fin.GetObject("ntp_residuals",ntuple);
17 
18  float seed_id;
19  float layer;
20  float dphi;
21  float dz;
22  float x;
23  float y;
24  float z;
25  float pt;
26  float px;
27  float py;
28  float pz;
29  float crossing;
30  float isSilicon;
31  float isTpc;
32 
33  ntuple->SetBranchAddress("seed_id",&seed_id);
34  ntuple->SetBranchAddress("layer",&layer);
35  ntuple->SetBranchAddress("dphi",&dphi);
36  ntuple->SetBranchAddress("dz",&dz);
37  ntuple->SetBranchAddress("x",&x);
38  ntuple->SetBranchAddress("y",&y);
39  ntuple->SetBranchAddress("z",&z);
40  ntuple->SetBranchAddress("pt",&pt);
41  ntuple->SetBranchAddress("px",&px);
42  ntuple->SetBranchAddress("py",&py);
43  ntuple->SetBranchAddress("pz",&pz);
44  ntuple->SetBranchAddress("crossing",&crossing);
45  ntuple->SetBranchAddress("isSilicon",&isSilicon);
46  ntuple->SetBranchAddress("isTpc",&isTpc);
47 
48 
49  int entries = ntuple->GetEntries();
50 
51  TH2D *xy = new TH2D("xy","cluster x vs. y",5000,-80,80,5000,-80,80);
52  TH2D *tpc_xy = new TH2D("tpc_xy","tpc cluster x vs. y",5000,-80,80,5000,-80,80);
53  TH2D *si_xy = new TH2D("si_xy","si cluster x vs. y",5000,-10,10,5000,-10,10);
54  TH2D *dphiLayer = new TH2D("dphiLayer","dphi vs. layer",5000,-0.02,0.02,5000,0,56);
55  TH2D *dzLayer = new TH2D("dzLayer","dz vs. layer",5000,-1,1,5000,0,56);
56  TH2D *resLayer = new TH2D("resLayer","residual vs. layer",5000,0,1,5000,0,56);
57  TH2D *dphiPt = new TH2D("dphiPt","dphi vs. pt",5000,-0.02,0.02,5000,0,7);
58  TH2D *dphiPz = new TH2D("dphiPz","dphi vs. pz",5000,-0.02,0.02,5000,0,5);
59 
60 
61  //take seedid count entries() whihc is nclusters) with seedid and switches
62  //count how many entries in ntuple have same seedid possibly thousands
63  //make seedid array
64  //loop over entries
65  // if seedid is allready in array add 1
66  // if seed id not in array add to array then add 1
67  // deltaphi v nhits
68 
69  //make multimap with key as seedid and add
70 
71  for(int i=0; i<entries; ++i)
72  {
73  ntuple->GetEntry(i);
74 
75  float r = sqrt(pow(x,2)+pow(y,2)+pow(z,2));
76  float residual = sqrt(pow(r*dphi,2)+pow(dz,2));
77 
78  resLayer->Fill(residual,layer);
79  dphiLayer->Fill(dphi,layer);
80  dzLayer->Fill(dz,layer);
81  xy->Fill(x,y);
82  dphiPt->Fill(dphi,pt);
83  dphiPz->Fill(dphi,pz);
84 
85  if(isSilicon==0)
86  {
87  tpc_xy->Fill(x,y);
88  }
89  if(isTpc==0)
90  {
91  si_xy->Fill(x,y);
92  }
93  }
94 
95  if(plotBool)
96  {
97  TCanvas *c1 = new TCanvas("c1","",10,10,800,800);
98  xy->DrawCopy();
99 
100  TCanvas *c2 = new TCanvas("c2","",10,10,800,800);
101  tpc_xy->DrawCopy();
102 
103  TCanvas *c3 = new TCanvas("c3","",10,10,600,600);
104  si_xy->DrawCopy();
105 
106  TCanvas *c4 = new TCanvas("c4","",10,10,800,800);
107  dphiLayer->DrawCopy();
108 
109  TCanvas *c5 = new TCanvas("c5","",10,10,800,800);
110  dzLayer->DrawCopy();
111 
112  TCanvas *c6 = new TCanvas("c6","",10,10,800,800);
113  resLayer->DrawCopy();
114 
115  TCanvas *c7 = new TCanvas("c7","",10,10,800,800);
116  dphiPt->DrawCopy();
117 
118  TCanvas *c8 = new TCanvas("c8","",10,10,800,800);
119  dphiPz->DrawCopy();
120  }
121 
122  if(savePlot)
123  {
124  std::unique_ptr<TFile> residualsFile(TFile::Open(outputFile.c_str(), "recreate"));
125  residualsFile->WriteObject(xy,"cluster_xvy");
126  residualsFile->WriteObject(tpc_xy,"tpc_cluster_xvy ");
127  residualsFile->WriteObject(si_xy,"si_cluster_xvy");
128  residualsFile->WriteObject(dphiLayer,"dphi_v_layer");
129  residualsFile->WriteObject(dzLayer,"dz_v_layer");
130  residualsFile->WriteObject(resLayer,"residual_v_layer");
131  residualsFile->WriteObject(dphiPt,"dphi_v_pt");
132  residualsFile->WriteObject(dphiPz,"dphi_v_pz");
133  }
134 
135 
136 }