Fork me on GitHub

source: git/examples/ExampleKsVtxFit.C

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

Add feature to force starting radius in vertex fit
  • Property mode set to 100644
File size: 6.8 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#include <TVectorD.h>
15#include <TMath.h>
16
17#endif
18
19
20
21//------------------------------------------------------------------------------
22
23void ExampleKsVtxFit(const char* inputFile, Int_t Nevent = 5)
24{
25 // Create chain of root trees
26 TChain chain("Delphes");
27 chain.Add(inputFile);
28
29 // Create object of class ExRootTreeReader
30 ExRootTreeReader* treeReader = new ExRootTreeReader(&chain);
31 Long64_t numberOfEntries = treeReader->GetEntries();
32
33 // Get pointers to branches used in this analysis
34 TClonesArray* branchGenPart = treeReader->UseBranch("Particle");
35 TClonesArray* branchTrack = treeReader->UseBranch("Track");
36
37 // Book histograms
38 Int_t Nbin = 100;
39 TH1D* hXpull = new TH1D("hXpull", "Pull X vertex component", Nbin, -10., 10.);
40 TH1D* hYpull = new TH1D("hYpull", "Pull Y vertex component", Nbin, -10., 10.);
41 TH1D* hZpull = new TH1D("hZpull", "Pull Z vertex component", Nbin, -10., 10.);
42 TH1D* hChi2 = new TH1D("hChi2", "Vertex #chi^{2}/N_{dof}", Nbin, 0., 10.);
43 //
44 // Loop over all events
45 Int_t Nev = TMath::Min(Nevent, (Int_t)numberOfEntries);
46 for (Int_t entry = 0; entry < Nev; ++entry)
47 {
48 // Load selected branches with data from specified event
49 treeReader->ReadEntry(entry);
50 //
51 Int_t NtrG = branchTrack->GetEntries();
52 const Int_t NtFit = 2;
53 TVectorD** pr = new TVectorD * [NtFit];
54 TMatrixDSym** cv = new TMatrixDSym * [NtFit];
55 std::cout<<"Start of event "<<entry<<std::endl;
56 //
57 // Search for generated Ks--> pi+ pi-
58 Int_t KsID = 310;
59 Int_t NKs = 0;
60 std::vector<GenParticle*> gpi1;
61 std::vector<GenParticle*> gpi2;
62 std::vector<Track*> gti1;
63 std::vector<Track*> gti2;
64 std::vector<TVectorD*> gKsVtx;
65 std::vector<TVectorD*> rKsVtx;
66 //
67 if (branchGenPart->GetEntries() > 1)
68 {
69 for(Int_t ig=0; ig<branchGenPart->GetEntries(); ig++)
70 {
71 GenParticle* gp = (GenParticle*)branchGenPart->At(ig);
72 if(gp->PID == KsID) // Found Ks
73 {
74 // Store daughters
75 gpi1.push_back( (GenParticle*)branchGenPart->At(gp->D1));
76 gpi2.push_back( (GenParticle*)branchGenPart->At(gp->D2));
77 Track* Tnull = 0;
78 gti1.push_back(Tnull);
79 gti2.push_back(Tnull);
80 // Store vertex
81 Double_t Vnull[3] = {0.,0.,0.};
82 TVectorD rVnull(3,Vnull);
83 Double_t Vgen[3] = {gpi1[NKs]->X,gpi1[NKs]->Y,gpi1[NKs]->Z};
84 TVectorD gVert(3,Vgen);
85 gKsVtx.push_back(new TVectorD(gVert));
86 rKsVtx.push_back(new TVectorD(rVnull));
87 //
88 NKs++;
89 }
90 }
91 } // Done searching
92 std::cout<<"Done searching. Found "<<NKs<<" K0s"<<std::endl;
93 //
94 // If event contains at least 1 Ks and 1 track
95 //
96 if (branchTrack->GetEntries() > 1 && NKs > 0)
97 {
98 // Loop on tracks
99 for (Int_t it = 0; it < branchTrack->GetEntries(); it++)
100 {
101 //std::cout<<"Start track loop. it ="<<it<<std::endl;
102 Track* trk = (Track*)branchTrack->At(it);
103 //
104 // Get associated generated particle
105 GenParticle* gt = (GenParticle*)trk->Particle.GetObject();
106 Int_t mp = gt->M1; // Mother
107 GenParticle* gm = (GenParticle*)branchGenPart->At(mp);
108 if(gm->PID == KsID) // Mother is Ks
109 {
110 for(Int_t k=0; k<NKs; k++)
111 {
112 if(gt == gpi1[k]) gti1[k] = trk;
113 if(gt == gpi2[k]) gti2[k] = trk;
114 }
115 }
116 //std::cout<<"End track loop"<<std::endl;
117 } // End loop on tracks
118 //
119 // Look for K0s to fit
120 //
121 for(Int_t k=0; k<NKs; k++) //Scan Ks
122 {
123 Bool_t found = (gti1[k] != 0) && (gti2[k] != 0);
124 //std::cout<<"Prepare Ks # "<<k<<" for fitting. Found = "<<found<<std::endl;
125 TVectorD gVtx = *gKsVtx[k];
126 if(gti1[k] != 0 && gti2[k] != 0) // if both tracks found
127 {
128 //std::cout<<"Both tracks found"<<std::endl;
129 // First leg
130 Double_t obsD0 = gti1[k]->D0;
131 Double_t obsPhi = gti1[k]->Phi;
132 Double_t obsC = gti1[k]->C;
133 Double_t obsZ0 = gti1[k]->DZ;
134 Double_t obsCtg = gti1[k]->CtgTheta;
135 //
136 Double_t oPar1[5] = { obsD0, obsPhi, obsC, obsZ0, obsCtg };
137 TVectorD obsPar1(5, oPar1); // Fill observed parameters
138 pr[0] = new TVectorD(obsPar1);
139 cv[0] = new TMatrixDSym(gti1[k]->CovarianceMatrix());
140 // Second leg
141 obsD0 = gti2[k]->D0;
142 obsPhi = gti2[k]->Phi;
143 obsC = gti2[k]->C;
144 obsZ0 = gti2[k]->DZ;
145 obsCtg = gti2[k]->CtgTheta;
146 //
147 Double_t oPar2[5] = { obsD0, obsPhi, obsC, obsZ0, obsCtg };
148 TVectorD obsPar2(5, oPar2); // Fill observed parameters
149 pr[1] = new TVectorD(obsPar2);
150 cv[1] = new TMatrixDSym(gti2[k]->CovarianceMatrix());
151 //
152 // Fit vertex
153 //
154 Double_t Rmin = 17.; // Lowest layer
155 Double_t Rvtx = TMath::Sqrt(gVtx(0)*gVtx(0)+gVtx(0)*gVtx(0));
156 VertexFit* Vtx = new VertexFit(NtFit, pr, cv);
157 if(Rvtx > Rmin)Vtx->SetStartR(Rvtx+10.);
158 //std::cout<<"BEFORE VERTEX FIT. Rvtx = "<<Rvtx<<std::endl;
159 TVectorD xvtx = Vtx->GetVtx();
160 TMatrixDSym covX = Vtx->GetVtxCov();
161 Double_t Chi2 = Vtx->GetVtxChi2();
162 //
163 // Fill plots
164 //
165 Double_t PullX = (xvtx(0)-gVtx(0)) / TMath::Sqrt(covX(0, 0));
166 Double_t PullY = (xvtx(1)-gVtx(1)) / TMath::Sqrt(covX(1, 1));
167 Double_t PullZ = (xvtx(2)-gVtx(2)) / TMath::Sqrt(covX(2, 2));
168 hXpull->Fill(PullX);
169 hYpull->Fill(PullY);
170 hZpull->Fill(PullZ);
171 hChi2->Fill(Chi2);
172 std::cout<<"xg: "<<gVtx(0)<<", yg: "<<gVtx(1)<<", zg: "<<gVtx(2)<<std::endl;
173 std::cout<<"xr: "<<xvtx(0)<<", yr: "<<xvtx(1)<<", zr: "<<xvtx(2)<<std::endl;
174 // Cleanup
175 for (Int_t i = 0; i < NtFit; i++) delete pr[i];
176 for (Int_t i = 0; i < NtFit; i++) delete cv[i];
177 }
178 } // End loop on Ks
179 //std::cout<<"Joust out of Ks loop"<<std::endl;
180 //
181 // Clean
182 delete[] pr;
183 delete[] cv;
184 for(Int_t k=0; k<NKs; k++)
185 {
186 delete gKsVtx[k];
187 delete rKsVtx[k];
188 }
189 gpi1.clear();;
190 gpi2.clear();
191 gti1.clear();
192 gti2.clear();
193 gKsVtx.clear();
194 rKsVtx.clear();
195 //std::cout<<"Finished cleanup"<<std::endl;
196 } // End if on Ks present
197 }
198 //
199 // Show resulting histograms
200 //
201 TCanvas* Cnv = new TCanvas("Cnv", "Delphes Ks vertex pulls", 50, 50, 900, 500);
202 Cnv->Divide(2, 2);
203 Cnv->cd(1); gPad->SetLogy(1); gStyle->SetOptFit(1111);
204 hXpull->Fit("gaus"); hXpull->Draw();
205 Cnv->cd(2); gPad->SetLogy(1); gStyle->SetOptFit(1111);
206 hYpull->Fit("gaus"); hYpull->Draw();
207 Cnv->cd(3); gPad->SetLogy(1); gStyle->SetOptFit(1111);
208 hZpull->Fit("gaus"); hZpull->Draw();
209 Cnv->cd(4); hChi2->Draw();
210}
Note: See TracBrowser for help on using the repository browser.