Fork me on GitHub

source: git/examples/ExamplePVtxFit.C@ bd33fce

Last change on this file since bd33fce was 91ef0b8, checked in by Franco BEDESCHI <bed@…>, 4 years ago

Fixed conflicts and Reginv3

  • Property mode set to 100644
File size: 4.1 KB
Line 
1/*
2Example of using vertex fitter class to fit primary vertex
3assumed to be generated in (0,0,0)
4*/
5
6#ifdef __CLING__
7R__LOAD_LIBRARY(libDelphes)
8
9#include "classes/DelphesClasses.h"
10#include "external/ExRootAnalysis/ExRootTreeReader.h"
11#include "modules/TrackCovariance.h"
12#include "external/TrackCovariance/TrkUtil.h"
13#include "external/TrackCovariance/VertexFit.h"
14
15#endif
16
17
18//------------------------------------------------------------------------------
19
20void ExamplePVtxFit(const char* inputFile, Int_t Nevent = 5)
21{
22 // Create chain of root trees
23 TChain chain("Delphes");
24 chain.Add(inputFile);
25
26 // Create object of class ExRootTreeReader
27 ExRootTreeReader* treeReader = new ExRootTreeReader(&chain);
28 Long64_t numberOfEntries = treeReader->GetEntries();
29
30 // Get pointers to branches used in this analysis
31 TClonesArray* branchGenPart = treeReader->UseBranch("Particle");
32 TClonesArray* branchTrack = treeReader->UseBranch("Track");
33
34 // Book histograms
35 Int_t Nbin = 100;
36 TH1D* hXpull = new TH1D("hXpull", "Pull X vertex component", Nbin, -10., 10.);
37 TH1D* hYpull = new TH1D("hYpull", "Pull Y vertex component", Nbin, -10., 10.);
38 TH1D* hZpull = new TH1D("hZpull", "Pull Z vertex component", Nbin, -10., 10.);
39 TH1D* hChi2 = new TH1D("hChi2", "Vertex #chi^{2}/N_{dof}", Nbin, 0., 10.);
40 //
41 // Loop over all events
42 Int_t Nev = TMath::Min(Nevent, (Int_t)numberOfEntries);
43 for (Int_t entry = 0; entry < Nev; ++entry)
44 {
45 // Load selected branches with data from specified event
46 treeReader->ReadEntry(entry);
47 Int_t Ntr = 0; // # of tracks from primary vertex
48 Int_t NtrG = branchTrack->GetEntries();
49 TVectorD** pr = new TVectorD * [NtrG];
50 TMatrixDSym** cv = new TMatrixDSym * [NtrG];
51 // If event contains at least 1 track
52 //
53 if (branchTrack->GetEntries() > 0)
54 {
55 // Loop on tracks
56 for (Int_t it = 0; it < branchTrack->GetEntries(); it++)
57 {
58 Track* trk = (Track*)branchTrack->At(it);
59 //
60 // Get associated generated particle
61 GenParticle* gp = (GenParticle*)trk->Particle.GetObject();
62 //
63 // Position of origin in mm
64 Double_t x = gp->X;
65 Double_t y = gp->Y;
66 Double_t z = gp->Z;
67 //
68 // group tracks originating from the primary vertex
69 if (x == 0.0 && y == 0.0)
70 {
71 //
72 // Reconstructed track parameters
73 Double_t obsD0 = trk->D0;
74 Double_t obsPhi = trk->Phi;
75 Double_t obsC = trk->C;
76 Double_t obsZ0 = trk->DZ;
77 Double_t obsCtg = trk->CtgTheta;
78 Double_t oPar[5] = { obsD0, obsPhi, obsC, obsZ0, obsCtg };
79 TVectorD obsPar(5, oPar); // Fill observed parameters
80 //
81 pr[Ntr] = new TVectorD(obsPar);
82 //std::cout<<"Cov Matrix:"<<std::endl;
83 //trk->CovarianceMatrix().Print();
84 cv[Ntr] = new TMatrixDSym(trk->CovarianceMatrix());
85 Ntr++;
86 }
87 } // End loop on tracks
88 }
89 //
90 // Fit primary vertex
91 Int_t MinTrk = 2; // Minumum # tracks for vertex fit
92 if (Ntr >= MinTrk) {
93 VertexFit* Vtx = new VertexFit(Ntr, pr, cv);
94 TVectorD xvtx = Vtx->GetVtx();
95 TMatrixDSym covX = Vtx->GetVtxCov();
96 Double_t Chi2 = Vtx->GetVtxChi2();
97 Double_t Ndof = 2 * (Double_t)Ntr - 3;
98 Double_t PullX = xvtx(0) / TMath::Sqrt(covX(0, 0));
99 Double_t PullY = xvtx(1) / TMath::Sqrt(covX(1, 1));
100 Double_t PullZ = xvtx(2) / TMath::Sqrt(covX(2, 2));
101 //
102 // Fill histograms
103 hXpull->Fill(PullX);
104 hYpull->Fill(PullY);
105 hZpull->Fill(PullZ);
106 hChi2->Fill(Chi2 / Ndof);
107 }
108
109 //std::cout << "Vertex chi2/Ndof = " << Chi2 / Ndof << std::endl;
110 //
111 // Cleanup
112 for (Int_t i = 0; i < Ntr; i++) delete pr[i];
113 for (Int_t i = 0; i < Ntr; i++) delete cv[i];
114 delete[] pr;
115 delete[] cv;
116 }
117 //
118 // Show resulting histograms
119 //
120 TCanvas* Cnv = new TCanvas("Cnv", "Delphes primary vertex pulls", 50, 50, 900, 500);
121 Cnv->Divide(2, 2);
122 Cnv->cd(1); gPad->SetLogy(1); gStyle->SetOptFit(1111);
123 hXpull->Fit("gaus"); hXpull->Draw();
124 Cnv->cd(2); gPad->SetLogy(1); gStyle->SetOptFit(1111);
125 hYpull->Fit("gaus"); hYpull->Draw();
126 Cnv->cd(3); gPad->SetLogy(1); gStyle->SetOptFit(1111);
127 hZpull->Fit("gaus"); hZpull->Draw();
128 Cnv->cd(4); hChi2->Draw();
129}
Note: See TracBrowser for help on using the repository browser.