analysis/runTrackingDiagnostics/plot_residuals.C

0001 #include <cmath>
0002 #include <TFile.h>
0003 #include <TNtuple.h>
0004 #include <TH2D.h>
0005 #include <TCut.h>
0006 #include <Eigen/Dense>
0007
0008 void plot_residuals(const std::string inputFile = "residuals_G4sPHENIX.root", const std::string outputFile = "residuals_plots.root")
0009 {
0010   bool plotBool = false; // Determines whether plots are displayed or not
0011   bool savePlot = true; // Determines whether histograms are saved to output file or not
0012
0013   TFile fin(inputFile.c_str());
0014
0015   TNtuple *ntuple;
0016   fin.GetObject("ntp_residuals",ntuple);
0017
0018   float seed_id;
0019   float layer;
0020   float dphi;
0021   float dz;
0022   float x;
0023   float y;
0024   float z;
0025   float pt;
0026   float px;
0027   float py;
0028   float pz;
0029   float crossing;
0030   float isSilicon;
0031   float isTpc;
0032
0033   ntuple->SetBranchAddress("seed_id",&seed_id);
0034   ntuple->SetBranchAddress("layer",&layer);
0035   ntuple->SetBranchAddress("dphi",&dphi);
0036   ntuple->SetBranchAddress("dz",&dz);
0037   ntuple->SetBranchAddress("x",&x);
0038   ntuple->SetBranchAddress("y",&y);
0039   ntuple->SetBranchAddress("z",&z);
0040   ntuple->SetBranchAddress("pt",&pt);
0041   ntuple->SetBranchAddress("px",&px);
0042   ntuple->SetBranchAddress("py",&py);
0043   ntuple->SetBranchAddress("pz",&pz);
0044   ntuple->SetBranchAddress("crossing",&crossing);
0045   ntuple->SetBranchAddress("isSilicon",&isSilicon);
0046   ntuple->SetBranchAddress("isTpc",&isTpc);
0047
0048
0049   int entries = ntuple->GetEntries();
0050
0051   TH2D *xy         = new TH2D("xy","cluster x vs. y",5000,-80,80,5000,-80,80);
0052   TH2D *tpc_xy     = new TH2D("tpc_xy","tpc cluster x vs. y",5000,-80,80,5000,-80,80);
0053   TH2D *si_xy      = new TH2D("si_xy","si cluster x vs. y",5000,-10,10,5000,-10,10);
0054   TH2D *dphiLayer  = new TH2D("dphiLayer","dphi vs. layer",5000,-0.02,0.02,5000,0,56);
0055   TH2D *dzLayer    = new TH2D("dzLayer","dz vs. layer",5000,-1,1,5000,0,56);
0056   TH2D *resLayer   = new TH2D("resLayer","residual vs. layer",5000,0,1,5000,0,56);
0057   TH2D *dphiPt     = new TH2D("dphiPt","dphi vs. pt",5000,-0.02,0.02,5000,0,7);
0058   TH2D *dphiPz     = new TH2D("dphiPz","dphi vs. pz",5000,-0.02,0.02,5000,0,5);
0059
0060
0061   //take seedid count entries() whihc is nclusters) with seedid and switches
0062   //count how many entries in ntuple have same seedid possibly thousands
0063   //make seedid array
0064   //loop over entries
0065   // if seedid is allready in array add 1
0066   // if seed id not in array add to array then add 1
0067   // deltaphi v nhits
0068
0069   //make multimap with key as seedid and add
0070
0071   for(int i=0; i<entries; ++i)
0072     {
0073       ntuple->GetEntry(i);
0074
0075       float r = sqrt(pow(x,2)+pow(y,2)+pow(z,2));
0076       float residual = sqrt(pow(r*dphi,2)+pow(dz,2));
0077
0078       resLayer->Fill(residual,layer);
0079       dphiLayer->Fill(dphi,layer);
0080       dzLayer->Fill(dz,layer);
0081       xy->Fill(x,y);
0082       dphiPt->Fill(dphi,pt);
0083       dphiPz->Fill(dphi,pz);
0084
0085       if(isSilicon==0)
0086     {
0087       tpc_xy->Fill(x,y);
0088     }
0089       if(isTpc==0)
0090     {
0091       si_xy->Fill(x,y);
0092     }
0093     }
0094
0095   if(plotBool)
0096     {
0097       TCanvas *c1 = new TCanvas("c1","",10,10,800,800);
0098       xy->DrawCopy();
0099
0100       TCanvas *c2 = new TCanvas("c2","",10,10,800,800);
0101       tpc_xy->DrawCopy();
0102
0103       TCanvas *c3 = new TCanvas("c3","",10,10,600,600);
0104       si_xy->DrawCopy();
0105
0106       TCanvas *c4 = new TCanvas("c4","",10,10,800,800);
0107       dphiLayer->DrawCopy();
0108
0109       TCanvas *c5 = new TCanvas("c5","",10,10,800,800);
0110       dzLayer->DrawCopy();
0111
0112       TCanvas *c6 = new TCanvas("c6","",10,10,800,800);
0113       resLayer->DrawCopy();
0114
0115       TCanvas *c7 = new TCanvas("c7","",10,10,800,800);
0116       dphiPt->DrawCopy();
0117
0118       TCanvas *c8 = new TCanvas("c8","",10,10,800,800);
0119       dphiPz->DrawCopy();
0120     }
0121
0122   if(savePlot)
0123     {
0124       std::unique_ptr<TFile> residualsFile(TFile::Open(outputFile.c_str(), "recreate"));
0125       residualsFile->WriteObject(xy,"cluster_xvy");
0126       residualsFile->WriteObject(tpc_xy,"tpc_cluster_xvy ");
0127       residualsFile->WriteObject(si_xy,"si_cluster_xvy");
0128       residualsFile->WriteObject(dphiLayer,"dphi_v_layer");
0129       residualsFile->WriteObject(dzLayer,"dz_v_layer");
0130       residualsFile->WriteObject(resLayer,"residual_v_layer");
0131       residualsFile->WriteObject(dphiPt,"dphi_v_pt");
0132       residualsFile->WriteObject(dphiPz,"dphi_v_pz");
0133     }
0134
0135
0136 }