Fork me on GitHub

source: git/examples/ExamplePVtxFitEC.C@ 4564cba

Last change on this file since 4564cba was b750b0a, checked in by Franco BEDESCHI <bed@…>, 3 years ago

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