Fork me on GitHub

source: git/examples/ExamplePVtxFit.C@ 7bca620

Last change on this file since 7bca620 was a47edcc, checked in by Franco BEDESCHI <bed@…>, 3 years ago

Fixed issue with backward generated tracks
  • Property mode set to 100644
File size: 5.3 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 ExamplePVtxFit(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
123 Int_t MinTrk = 2; // Minumum # tracks for vertex fit
124 if (Ntr >= MinTrk) {
125 VertexFit* Vtx = new VertexFit(Ntr, pr, cv);
126 TVectorD xvtx = Vtx->GetVtx();
127 TMatrixDSym covX = Vtx->GetVtxCov();
128 Double_t Chi2 = Vtx->GetVtxChi2();
129 Double_t Ndof = 2 * (Double_t)Ntr - 3;
130 Double_t PullX = (xvtx(0)-xpv) / TMath::Sqrt(covX(0, 0));
131 Double_t PullY = (xvtx(1)-ypv) / TMath::Sqrt(covX(1, 1));
132 Double_t PullZ = (xvtx(2)-zpv) / TMath::Sqrt(covX(2, 2));
133 //std::cout<<"**** True vertex (x, y, z) "<<xpv<<", "
134 //<<ypv<<", "<<zpv<<std::endl;
135 //std::cout<<"**** Found vertex (x, y, z) "<<xvtx(0)<<", "
136 //<<xvtx(1)<<", "<<xvtx(2)<<std::endl;
137 //
138 // Fill histograms
139 hXpull->Fill(PullX);
140 hYpull->Fill(PullY);
141 hZpull->Fill(PullZ);
142 hChi2->Fill(Chi2 / Ndof);
143 }
144
145 //std::cout << "Vertex chi2/Ndof = " << Chi2 / Ndof << std::endl;
146 //
147 // Cleanup
148 for (Int_t i = 0; i < Ntr; i++) delete pr[i];
149 for (Int_t i = 0; i < Ntr; i++) delete cv[i];
150 delete[] pr;
151 delete[] cv;
152 }
153 //
154 // Show resulting histograms
155 //
156 TCanvas* Cnv = new TCanvas("Cnv", "Delphes primary vertex pulls", 50, 50, 900, 500);
157 Cnv->Divide(2, 2);
158 Cnv->cd(1); gPad->SetLogy(1); gStyle->SetOptFit(1111);
159 hXpull->Fit("gaus"); hXpull->Draw();
160 Cnv->cd(2); gPad->SetLogy(1); gStyle->SetOptFit(1111);
161 hYpull->Fit("gaus"); hYpull->Draw();
162 Cnv->cd(3); gPad->SetLogy(1); gStyle->SetOptFit(1111);
163 hZpull->Fit("gaus"); hZpull->Draw();
164 Cnv->cd(4); hChi2->Draw();
165}
Note: See TracBrowser for help on using the repository browser.