Skip to content

Commit

Permalink
Merge pull request #118 from fbedesch/master
Browse files Browse the repository at this point in the history
Fixed error handling with external vertex constraint and test example
  • Loading branch information
selvaggi committed Mar 14, 2023
2 parents 68eead1 + 225a839 commit 7b7cc48
Show file tree
Hide file tree
Showing 5 changed files with 285 additions and 66 deletions.
83 changes: 31 additions & 52 deletions examples/ExamplePVtxFitEC.C
Original file line number Diff line number Diff line change
Expand Up @@ -11,6 +11,7 @@ R__LOAD_LIBRARY(libDelphes)
#include "modules/TrackCovariance.h"
#include "external/TrackCovariance/TrkUtil.h"
#include "external/TrackCovariance/VertexFit.h"
#include "VtxPulls.h"

#endif

Expand Down Expand Up @@ -49,17 +50,7 @@ void ExamplePVtxFitEC(const char* inputFile, Int_t Nevent = 5)
TClonesArray* branchTrack = treeReader->UseBranch("Track");

// Book histograms
Int_t Nbin = 100;
// Vertex fit pulls
TH1D* hXpull = new TH1D("hXpull", "Pull X vertex component", Nbin, -10., 10.);
TH1D* hYpull = new TH1D("hYpull", "Pull Y vertex component", Nbin, -10., 10.);
TH1D* hZpull = new TH1D("hZpull", "Pull Z vertex component", Nbin, -10., 10.);
TH1D* hChi2 = new TH1D("hChi2", "Vertex #chi^{2}/N_{dof}", Nbin, 0., 10.);
// Generation check
TH1D* hXvpull = new TH1D("hXvpull", "Pull X generated vertex", Nbin, -10., 10.);
TH1D* hYvpull = new TH1D("hYvpull", "Pull Y generated vertex", Nbin, -10., 10.);
TH1D* hZvpull = new TH1D("hZvpull", "Pull Z generated vertex", Nbin, -10., 10.);

VtxPulls* Pulls = new VtxPulls();
//
// Loop over all events
Int_t Nev = TMath::Min(Nevent, (Int_t)numberOfEntries);
Expand All @@ -72,6 +63,15 @@ void ExamplePVtxFitEC(const char* inputFile, Int_t Nevent = 5)
Int_t NtrG = branchTrack->GetEntries();
TVectorD** pr = new TVectorD * [NtrG];
TMatrixDSym** cv = new TMatrixDSym * [NtrG];
TVectorD vTrue(3);
std::vector<TVectorD> pGen; // True momenta
//
// Print event numbers
if(entry % 500 == 0){
std::cout<<std::endl;
std::cout <<"Event " << entry << std::endl;
std::cout<<"Nr. of tracks " << NtrG << std::endl;
}
//
// True vertex
Double_t xpv, ypv, zpv;
Expand Down Expand Up @@ -122,6 +122,9 @@ void ExamplePVtxFitEC(const char* inputFile, Int_t Nevent = 5)
xpv = x;
ypv = y;
zpv = z;
vTrue(0) = x;
vTrue(1) = y;
vTrue(2) = z;
//
// Reconstructed track parameters
Double_t obsD0 = trk->D0;
Expand All @@ -135,40 +138,32 @@ void ExamplePVtxFitEC(const char* inputFile, Int_t Nevent = 5)
pr[Ntr] = new TVectorD(obsPar);
cv[Ntr] = new TMatrixDSym(trk->CovarianceMatrix());
Ntr++;
// Store true momenta
GenParticle* gen = (GenParticle*) trk->Particle.GetObject();
TVectorD GenP(3);
GenP(0) = gen->Px;
GenP(1) = gen->Py;
GenP(2) = gen->Pz;
pGen.push_back(GenP);
}
} // End loop on tracks
}
//
// Fit primary vertex with beam constraint
//
Int_t MinTrk = 3; // Minumum # tracks for vertex fit
Int_t MinTrk = 1; // Minumum # tracks for vertex fit
if (Ntr >= MinTrk) {
VertexFit* Vtx = new VertexFit(Ntr, pr, cv);
Vtx->AddVtxConstraint(xpvc, covpvc);
TVectorD xvtx = Vtx->GetVtx();
TMatrixDSym covX = Vtx->GetVtxCov();
Double_t Chi2 = Vtx->GetVtxChi2();
Double_t Ndof = 2 * (Double_t)Ntr;
//std::cout<<"N tracks = "<<Ntr<<", X,Y,Z in: "<<xpvc(0)<<", "<<xpvc(1)<<", "<<xpvc(2)<<std::endl;
//std::cout<<"X, Y, Z out: "<<xvtx(0)<<", "<<xvtx(1)<<", "<<xvtx(2)<<std::endl;
for(Int_t i=0; i<Ntr; i++){
TVectorD pMom = pGen[i];
//pMom.Print();
}
Pulls->Fill(vTrue, pGen, Vtx);
delete Vtx;
//
Double_t PullX = (xvtx(0)-xpv) / TMath::Sqrt(covX(0, 0));
Double_t PullY = (xvtx(1)-ypv) / TMath::Sqrt(covX(1, 1));
Double_t PullZ = (xvtx(2)-zpv) / TMath::Sqrt(covX(2, 2));
//
Double_t PullXv = (xpvc(0)-xpv) / TMath::Sqrt(covpvc(0, 0));
Double_t PullYv = (xpvc(1)-ypv) / TMath::Sqrt(covpvc(1, 1));
Double_t PullZv = (xpvc(2)-zpv) / TMath::Sqrt(covpvc(2, 2));
//
//
// Fill histograms
hXpull->Fill(PullX);
hYpull->Fill(PullY);
hZpull->Fill(PullZ);
hChi2->Fill(Chi2 / Ndof);
//
hXvpull->Fill(PullXv);
hYvpull->Fill(PullYv);
hZvpull->Fill(PullZv);
}

//
Expand All @@ -177,27 +172,11 @@ void ExamplePVtxFitEC(const char* inputFile, Int_t Nevent = 5)
for (Int_t i = 0; i < Ntr; i++) delete cv[i];
delete[] pr;
delete[] cv;
pGen.clear();
}

//
// Show resulting histograms
//
TCanvas* Cnv = new TCanvas("Cnv", "Delphes primary vertex pulls", 50, 50, 900, 500);
Cnv->Divide(2, 2);
Cnv->cd(1); gPad->SetLogy(1); gStyle->SetOptFit(1111);
hXpull->Fit("gaus"); hXpull->Draw();
Cnv->cd(2); gPad->SetLogy(1); gStyle->SetOptFit(1111);
hYpull->Fit("gaus"); hYpull->Draw();
Cnv->cd(3); gPad->SetLogy(1); gStyle->SetOptFit(1111);
hZpull->Fit("gaus"); hZpull->Draw();
Cnv->cd(4); hChi2->Draw();
//
TCanvas* Cnv1 = new TCanvas("Cnv1", "Generated primary vertex pulls", 100, 100, 900, 500);
Cnv1->Divide(3, 1);
Cnv1->cd(1); gPad->SetLogy(1); gStyle->SetOptFit(1111);
hXvpull->Fit("gaus"); hXvpull->Draw();
Cnv1->cd(2); gPad->SetLogy(1); gStyle->SetOptFit(1111);
hYvpull->Fit("gaus"); hYvpull->Draw();
Cnv1->cd(3); gPad->SetLogy(1); gStyle->SetOptFit(1111);
hZvpull->Fit("gaus"); hZvpull->Draw();
Pulls->Print();
}
159 changes: 159 additions & 0 deletions examples/VtxPulls.cc
Original file line number Diff line number Diff line change
@@ -0,0 +1,159 @@
#include "VtxPulls.h"

//
// Final state configuration description
//
// Constructor
VtxPulls::VtxPulls()
{
//
// Pulls
// Ds vertex
h_Xv = new TH1D("h_Xv","X-vertex pull",100,-10.,10.);
h_Yv = new TH1D("h_Yv","Y-vertex pull",100,-10.,10.);
h_Zv = new TH1D("h_Zv","Z-vertex pull",100,-10.,10.);
h_Chi2 = new TH1D("h_Chi2", "Vertex #chi^{2}/N_{dof}", 100, 0., 10.);
// Updated track parameter pulls after vertex fit
h_D = new TH1D("h_D","D impact parameter pull",100,-10.,10.);
h_P0 = new TH1D("h_P0","Phi0 parameter pull",100,-10.,10.);
h_C = new TH1D("h_C","Half curvature pull",100,-10.,10.);
h_Z0 = new TH1D("h_Z0","Z0 pull",100,-10.,10.);
h_Ct = new TH1D("h_Ct","Cot theta pull",100,-10.,10.);
// track momenta at vertex
h_Px = new TH1D("h_Px","Momentum x pull",100,-10.,10.);
h_Py = new TH1D("h_Py","Momentum y pull",100,-10.,10.);
h_Pz = new TH1D("h_Pz","Momentum z pull",100,-10.,10.);
}

//
// Destructor
VtxPulls::~VtxPulls()
{
}

void VtxPulls::Fill(TVectorD vTrue, std::vector<TVectorD> pTrue, VertexFit *Vfit)
{
// Global variables
Double_t Bz = 2.; // Magnetic field
//
// Vertex pulls
Double_t VxPull[3];
TVectorD XvFit = Vfit->GetVtx();
TMatrixDSym XvCov = Vfit->GetVtxCov();
Int_t Ntr = Vfit->GetNtrk();
Double_t Chi2 = Vfit->GetVtxChi2();
for(Int_t i=0; i<3; i++) VxPull[i] = (XvFit(i)-vTrue[i])/TMath::Sqrt(XvCov(i,i));
h_Xv->Fill(VxPull[0]);
h_Yv->Fill(VxPull[1]);
h_Zv->Fill(VxPull[2]);
//std::cout<<"Pulls: Vx = "<<VxPull[0]<<", "<<VxPull[1]<<", "<<VxPull[2]<<std::endl;
// Chi/Ndof
Double_t Ndof = 2 * (Double_t)Ntr;
h_Chi2->Fill(Chi2/Ndof);
//
// Updated track parameter pulls
Int_t Ntst = (Int_t) pTrue.size();
if(Ntr != Ntst) {
std::cout<<"VtxPulls:: mismatch on input track lists. Abort."<<std::endl;
exit(-1);
}else{
Bool_t Mm = kTRUE;
VertexMore* Vmore = new VertexMore(Vfit, Mm);
// Track loop
for(Int_t i=0; i<Ntr; i++){
TVector3 pTrk(pTrue[i](0),pTrue[i](1),pTrue[i](2));
TVector3 Xvm(vTrue[0]*1.e-3,vTrue[1]*1.e-3,vTrue[2]*1.e-3);
Double_t Q = Vmore->GetCharge(i);
TVectorD gPar = TrkUtil:: XPtoPar(Xvm, pTrk, Q, Bz);
TVectorD GenPar = TrkUtil::ParToMm(gPar); // Back to mm
TVectorD NewPar = Vfit->GetNewPar(i);
TMatrixDSym ParCov = Vfit->GetNewCov(i);
//
// Track Parameter pulls
Double_t ParPull[5];
for(Int_t k=0; k<5; k++)ParPull[k] = (NewPar(k)-GenPar(k))/TMath::Sqrt(ParCov(k,k));
h_D ->Fill(ParPull[0]);
h_P0->Fill(ParPull[1]);
h_C ->Fill(ParPull[2]);
h_Z0->Fill(ParPull[3]);
h_Ct->Fill(ParPull[4]);
//
// Get fit momenta momenta
//
TVector3 pRec = Vmore->GetMomentum(i);
TMatrixDSym pCov = Vmore->GetMomentumC(i);
Double_t MomPull[3];
for(Int_t k=0; k<3; k++)MomPull[k] = (pRec(k)-pTrk(k))/TMath::Sqrt(pCov(k,k));
h_Px->Fill(MomPull[0]);
h_Py->Fill(MomPull[1]);
h_Pz->Fill(MomPull[2]);
}
delete Vmore;
}

}

// Print
void VtxPulls::Print()
{
//
// Vertex position pulls
TCanvas * cnv = new TCanvas("cnv","Vertex position pulls",10,10, 600,600);
cnv->Divide(2,2);
// X
cnv->cd(1); gPad->SetLogy(1);
gStyle->SetOptStat(111111); gStyle->SetOptFit(1111);
h_Xv->Fit("gaus"); h_Xv->Draw();
// Y
cnv->cd(2); gPad->SetLogy(1);
gStyle->SetOptStat(111111); gStyle->SetOptFit(1111);
h_Yv->Fit("gaus"); h_Yv->Draw();
// Z
cnv->cd(3); gPad->SetLogy(1);
gStyle->SetOptStat(111111); gStyle->SetOptFit(1111);
h_Zv->Fit("gaus"); h_Zv->Draw();
// Chi2/Ndof
cnv->cd(4);
gStyle->SetOptStat(111111); gStyle->SetOptFit(1111);
h_Chi2->Draw();
//
// Updated track parameter pulls
TCanvas * cnv1 = new TCanvas("cnv1","Updated parameter pulls",50,50, 900,600);
cnv1->Divide(3,2);
// D
cnv1->cd(1); gPad->SetLogy(1);
gStyle->SetOptStat(111111); gStyle->SetOptFit(1111);
h_D->Fit("gaus"); h_D->Draw();
// Phi0
cnv1->cd(2); gPad->SetLogy(1);
gStyle->SetOptStat(111111); gStyle->SetOptFit(1111);
h_P0->Fit("gaus"); h_P0->Draw();
// C
cnv1->cd(3); gPad->SetLogy(1);
gStyle->SetOptStat(111111); gStyle->SetOptFit(1111);
h_C->Fit("gaus"); h_C->Draw();
// Z0
cnv1->cd(4); gPad->SetLogy(1);
gStyle->SetOptStat(111111); gStyle->SetOptFit(1111);
h_Z0->Fit("gaus"); h_Z0->Draw();
// Cotg Theta
cnv1->cd(5); gPad->SetLogy(1);
gStyle->SetOptStat(111111); gStyle->SetOptFit(1111);
h_Ct->Fit("gaus"); h_Ct->Draw();
//
// Momenta pulls
TCanvas * cnv2 = new TCanvas("cnv2","Momenta pulls",100,100, 900,600);
cnv2->Divide(3,1);
// Px
cnv2->cd(1); gPad->SetLogy(1);
gStyle->SetOptStat(111111); gStyle->SetOptFit(1111);
h_Px->Fit("gaus"); h_Px->Draw();
// Py
cnv2->cd(2); gPad->SetLogy(1);
gStyle->SetOptStat(111111); gStyle->SetOptFit(1111);
h_Py->Fit("gaus"); h_Py->Draw();
// Pz
cnv2->cd(3); gPad->SetLogy(1);
gStyle->SetOptStat(111111); gStyle->SetOptFit(1111);
h_Pz->Fit("gaus"); h_Pz->Draw();
}
51 changes: 51 additions & 0 deletions examples/VtxPulls.h
Original file line number Diff line number Diff line change
@@ -0,0 +1,51 @@
//
#ifndef G__VTXPULLS_H
#define G__VTXPULLS_H
//
#include <iostream>
#include <TH1.h>
#include <TCanvas.h>
#include <TStyle.h>
#include "external/TrackCovariance/VertexFit.h"
#include "external/TrackCovariance/VertexMore.h"

//
// Pulls of tracks parameters and vertex after vertex fit

class VtxPulls{
//
//
// Author: F. Bedeschi, INFN-Pisa, Italy
// March 6, 2023
//
private:
// x vertex
TH1D* h_Xv;
TH1D* h_Yv;
TH1D* h_Zv;
TH1D* h_Chi2;
// Track parameters after fit
TH1D* h_D;
TH1D* h_P0;
TH1D* h_C;
TH1D* h_Z0;
TH1D* h_Ct;
// Track momenta after fit
TH1D* h_Px;
TH1D* h_Py;
TH1D* h_Pz;
public:
//
// Constructors
VtxPulls();
// Destructor
~VtxPulls();
//
//
// Methods
//
void Fill(TVectorD vTrue, std::vector<TVectorD> pTrue, VertexFit *Vfit);
void Print();
};

#endif
Loading

0 comments on commit 7b7cc48

Please sign in to comment.