Fork me on GitHub

source: git/examples/ExamplePVtxFitEC.C@ d612dec

Last change on this file since d612dec was 46d3442, checked in by Franco BEDESCHI <bed@…>, 4 years ago

Vertex contrained fit and primary vertex finder added
  • Property mode set to 100644
File size: 5.6 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//------------------------------------------------------------------------------
20
21void ExamplePVtxFitEC(const char* inputFile, Int_t Nevent = 5)
22{
23 // Create chain of root trees
24 TChain chain("Delphes");
25 chain.Add(inputFile);
26
27 // Create object of class ExRootTreeReader
28 ExRootTreeReader* treeReader = new ExRootTreeReader(&chain);
29 Long64_t numberOfEntries = treeReader->GetEntries();
30
31 // Get pointers to branches used in this analysis
32 TClonesArray* branchGenPart = treeReader->UseBranch("Particle");
33 TClonesArray* branchTrack = treeReader->UseBranch("Track");
34
35 // Book histograms
36 Int_t Nbin = 100;
37 TH1D* hXpull = new TH1D("hXpull", "Pull X vertex component", Nbin, -10., 10.);
38 TH1D* hYpull = new TH1D("hYpull", "Pull Y vertex component", Nbin, -10., 10.);
39 TH1D* hZpull = new TH1D("hZpull", "Pull Z vertex component", Nbin, -10., 10.);
40 TH1D* hChi2 = new TH1D("hChi2", "Vertex #chi^{2}/N_{dof}", Nbin, 0., 10.);
41 //
42 // Loop over all events
43 Int_t Nev = TMath::Min(Nevent, (Int_t)numberOfEntries);
44 for (Int_t entry = 0; entry < Nev; ++entry)
45 {
46 // Load selected branches with data from specified event
47 treeReader->ReadEntry(entry);
48 //
49 Int_t Ntr = 0; // # of tracks from primary vertex
50 Int_t NtrG = branchTrack->GetEntries();
51 TVectorD** pr = new TVectorD * [NtrG];
52 TMatrixDSym** cv = new TMatrixDSym * [NtrG];
53 //
54 // True vertex
55 Double_t xpv, ypv, zpv;
56 // If event contains at least 1 track
57 //
58 if (branchTrack->GetEntries() > 0)
59 {
60 // Loop on tracks
61 for (Int_t it = 0; it < branchTrack->GetEntries(); it++)
62 {
63 Track* trk = (Track*)branchTrack->At(it);
64 //
65 // Get associated generated particle
66 GenParticle* gp = (GenParticle*)trk->Particle.GetObject();
67 TVector3 xg(1.e-3*gp->X,1.e-3*gp->Y,1.e-3*gp->Z);
68 TVector3 pg(gp->Px,gp->Py,gp->Pz);
69 Double_t Q = (Double_t)gp->Charge;
70 Double_t Bz = 2.0;
71 TVectorD genParM =TrkUtil:: XPtoPar(xg, pg, Q, Bz);
72 TVectorD genPar = TrkUtil::ParToMm(genParM);
73 //
74 // Position of origin in mm
75 Double_t x = gp->X;
76 Double_t y = gp->Y;
77 Double_t z = gp->Z;
78 Bool_t prim = kTRUE; // Is primary?
79 Int_t mp = gp->M1; // Mother
80 while(mp>0){
81 GenParticle* gm =
82 (GenParticle*)branchGenPart->At(mp);
83 Double_t xm = gm->X;
84 Double_t ym = gm->Y;
85 Double_t zm = gm->Z;
86 if(x!=xm || y!=ym || z!=zm){
87 prim = kFALSE;
88 break;
89 }else mp = gm->M1;
90 }
91 //
92 // group tracks originating from the primary vertex
93 if (prim)
94 {
95 //std::cout<<"Event: "<<entry<<", track: "<<it;
96 //std::cout<<", x = "<<x<<", y = "<<y<<", z = "<<z<<std::endl;
97 xpv = x;
98 ypv = y;
99 zpv = z;
100 //
101 // Reconstructed track parameters
102 Double_t obsD0 = trk->D0;
103 Double_t obsPhi = trk->Phi;
104 Double_t obsC = trk->C;
105 Double_t obsZ0 = trk->DZ;
106 //std::cout<<"Z0 track = "<< obsZ0
107 //<<", gen Z0 = "<<genPar(3)
108 //<<", gen cot = "<<genPar(4)<<std::endl;
109 Double_t obsCtg = trk->CtgTheta;
110 Double_t oPar[5] = { obsD0, obsPhi, obsC, obsZ0, obsCtg };
111 TVectorD obsPar(5, oPar); // Fill observed parameters
112 //
113 pr[Ntr] = new TVectorD(obsPar);
114 //std::cout<<"Cov Matrix:"<<std::endl;
115 //trk->CovarianceMatrix().Print();
116 cv[Ntr] = new TMatrixDSym(trk->CovarianceMatrix());
117 Ntr++;
118 }
119 } // End loop on tracks
120 }
121 //
122 // Fit primary vertex with beam constraint
123 //
124 //Beam constraint
125 TVectorD xpvc(3);
126 xpvc(0) = 1.0;
127 xpvc(1) = -2.0;
128 xpvc(2) = 10.0;
129 TMatrixDSym covpvc(3); covpvc.Zero();
130 covpvc(0, 0) = 0.0097 * 0.0097;
131 covpvc(1, 1) = 2.55e-05 * 2.55e-05;
132 covpvc(2, 2) = 0.64 * 0.64;
133 Int_t MinTrk = 2; // Minumum # tracks for vertex fit
134 if (Ntr >= MinTrk) {
135 VertexFit* Vtx = new VertexFit(Ntr, pr, cv);
136 Vtx->AddVtxConstraint(xpvc, covpvc);
137 TVectorD xvtx = Vtx->GetVtx();
138 TMatrixDSym covX = Vtx->GetVtxCov();
139 Double_t Chi2 = Vtx->GetVtxChi2();
140 Double_t Ndof = 2 * (Double_t)Ntr;
141 Double_t PullX = (xvtx(0)-xpv) / TMath::Sqrt(covX(0, 0));
142 Double_t PullY = (xvtx(1)-ypv) / TMath::Sqrt(covX(1, 1));
143 Double_t PullZ = (xvtx(2)-zpv) / TMath::Sqrt(covX(2, 2));
144 //std::cout<<"**** True vertex (x, y, z) "<<xpv<<", "
145 //<<ypv<<", "<<zpv<<std::endl;
146 //std::cout<<"**** Found vertex (x, y, z) "<<xvtx(0)<<", "
147 //<<xvtx(1)<<", "<<xvtx(2)<<std::endl;
148 //
149 // Fill histograms
150 hXpull->Fill(PullX);
151 hYpull->Fill(PullY);
152 hZpull->Fill(PullZ);
153 hChi2->Fill(Chi2 / Ndof);
154 }
155
156 //std::cout << "Vertex chi2/Ndof = " << Chi2 / Ndof << std::endl;
157 //
158 // Cleanup
159 for (Int_t i = 0; i < Ntr; i++) delete pr[i];
160 for (Int_t i = 0; i < Ntr; i++) delete cv[i];
161 delete[] pr;
162 delete[] cv;
163 }
164 //
165 // Show resulting histograms
166 //
167 TCanvas* Cnv = new TCanvas("Cnv", "Delphes primary vertex pulls", 50, 50, 900, 500);
168 Cnv->Divide(2, 2);
169 Cnv->cd(1); gPad->SetLogy(1); gStyle->SetOptFit(1111);
170 hXpull->Fit("gaus"); hXpull->Draw();
171 Cnv->cd(2); gPad->SetLogy(1); gStyle->SetOptFit(1111);
172 hYpull->Fit("gaus"); hYpull->Draw();
173 Cnv->cd(3); gPad->SetLogy(1); gStyle->SetOptFit(1111);
174 hZpull->Fit("gaus"); hZpull->Draw();
175 Cnv->cd(4); hChi2->Draw();
176}
Note: See TracBrowser for help on using the repository browser.