Fork me on GitHub

source: git/examples/vertexAnalyzer.C

Last change on this file was 23389ff, checked in by Michele Selvaggi <michele.selvaggi@…>, 8 years ago

first commit of examples for testing vertexing algos

  • Property mode set to 100644
File size: 4.7 KB
Line 
1/*
2Simple macro showing how to access vertex information, and tracks associated to vertices.
3Computes vertex reconstruction efficiencies, merging, fake and duplicate rates.
4The input file needs to be produced with cards/CMS_PhaseII/testVertexing.tcl
5
6root -l examples/vertexAnalyzer.C'("delphes_output.root")'
7*/
8
9#ifdef __CLING__
10R__LOAD_LIBRARY(libDelphes)
11#include "classes/DelphesClasses.h"
12#include "external/ExRootAnalysis/ExRootTreeReader.h"
13#endif
14
15//------------------------------------------------------------------------------
16
17void vertexAnalyzer(const char *inputFile)
18{
19 gSystem->Load("libDelphes");
20
21 // Create chain of root trees
22 TChain chain("Delphes");
23 chain.Add(inputFile);
24
25 // Create object of class ExRootTreeReader
26 ExRootTreeReader *treeReader = new ExRootTreeReader(&chain);
27 Long64_t numberOfEntries = treeReader->GetEntries();
28
29 // Get pointers to branches used in this analysis
30 TClonesArray *branchParticle = treeReader->UseBranch("Particle");
31 TClonesArray *branchGenVtx = treeReader->UseBranch("GenVertex");
32 TClonesArray *branchVtx = treeReader->UseBranch("Vertex");
33 TClonesArray *branchTrack = treeReader->UseBranch("Track");
34
35 Vertex *bestgenvtx, *bestvtx, *vtx, *recovtx;
36
37 Track *track;
38 GenParticle *particle;
39
40 TObject *object;
41
42 Double_t sumpt2;
43
44 Int_t nGen = 0, nReco = 0, ndenReco = 0, nMerge = 0;
45 Int_t nGenTot = 0, nRecoTot = 0, nFake = 0, nDuplicate = 0, nOverall = 0;
46 Int_t nmatch;
47
48 Double_t distance;
49 // Loop over all events
50 for(Int_t entry = 0; entry < numberOfEntries; ++entry)
51 {
52 //cout<<"------------- New event -------------------- "<<endl;
53 //cout<<" "<<endl;
54 //cout<<" Reconstructed vertices: "<<branchGenVtx->GetEntriesFast()<<endl;
55
56 // Load selected branches with data from specified event
57 treeReader->ReadEntry(entry);
58
59 sumpt2 = - 999;
60
61 // If event contains at least 1 genvtx
62 for(Int_t i = 0; i < branchGenVtx->GetEntriesFast(); ++i)
63 {
64
65 nGenTot++;
66 // Take vtx
67 vtx = (Vertex*) branchGenVtx->At(i);
68 //cout<<"gen vertex, Sumpt2: "<<vtx->SumPT2<<", Z: "<<vtx->Z<<", T:"<<vtx->T<<endl;
69
70 // Loop over gen vertex attached tracks
71 for(Int_t j = 0; j < vtx->Constituents.GetEntriesFast(); ++j)
72 {
73 object = vtx->Constituents.At(j);
74
75 // Check if the constituent is accessible
76 if(object == 0) continue;
77
78 if(object->IsA() == GenParticle::Class())
79 {
80 particle = (GenParticle*) object;
81 //cout << " GenPart pt: " << particle->PT << ", z: " << particle->Z << ", t: " << particle->T <<", PID: "<<particle->PID<<", Status: "<<particle->Status<<", PU: "<<particle->IsPU<<endl;
82
83 if(particle->IsPU == 0) bestgenvtx = (Vertex*) branchGenVtx->At(i);
84 }
85
86
87 }
88
89 nmatch = 0;
90 // loop of reco vertices
91 for(Int_t j = 0; j < branchVtx->GetEntriesFast(); ++j)
92 {
93
94 recovtx = (Vertex*) branchVtx->At(j);
95
96 distance = TMath::Sqrt( (vtx->Z - recovtx->Z)*(vtx->Z - recovtx->Z)/(recovtx->ErrorZ*recovtx->ErrorZ) +
97 (vtx->T - recovtx->T)*(vtx->T - recovtx->T)/(recovtx->ErrorT*recovtx->ErrorT) );
98
99 if(distance < 3.0)
100 {
101 nmatch++;
102 }
103
104 }
105
106 if(nmatch > 1) nDuplicate++;
107 if(nmatch > 0) nOverall++;
108
109 }
110
111 nGen++;
112
113 Bool_t found = false;
114
115 for(Int_t i = 0; i < branchVtx->GetEntriesFast(); ++i)
116 {
117
118 nRecoTot++;
119 recovtx = (Vertex*) branchVtx->At(i);
120
121 distance = TMath::Sqrt( (bestgenvtx->Z - recovtx->Z)*(bestgenvtx->Z - recovtx->Z)/(recovtx->ErrorZ*recovtx->ErrorZ) +
122 (bestgenvtx->T - recovtx->T)*(bestgenvtx->T - recovtx->T)/(recovtx->ErrorT*recovtx->ErrorT) );
123
124
125 if(distance < 3.0 && !found)
126 {
127 nReco++;
128 found = true;
129 }
130
131
132 nmatch = 0;
133 // start loop over gen vertices to find out merging rate
134 for(Int_t i = 0; i < branchGenVtx->GetEntriesFast(); ++i)
135 {
136
137 // Take vtx
138 vtx = (Vertex*) branchGenVtx->At(i);
139 // Loop over gen vertex attached tracks
140 distance = TMath::Sqrt( (vtx->Z - recovtx->Z)*(vtx->Z - recovtx->Z)/(recovtx->ErrorZ*recovtx->ErrorZ) +
141 (vtx->T - recovtx->T)*(vtx->T - recovtx->T)/(recovtx->ErrorT*recovtx->ErrorT) );
142
143 if(distance < 3.0)
144 {
145 nmatch++;
146 }
147
148
149 }
150
151 if(nmatch > 1) nMerge++;
152 if(nmatch == 0) nFake++;
153
154
155 }
156
157 } // end event loop
158
159 cout<<"hard vertex reco eff: "<<double(nReco)/nGen<<endl;
160 cout<<"all vertex reco eff: "<<double(nOverall)/nGenTot<<endl;
161 cout<<"merge rate: "<<double(nMerge)/nRecoTot<<endl;
162 cout<<"fake rate: "<<double(nFake)/nRecoTot<<endl;
163 cout<<"duplicate rate: "<<double(nDuplicate)/nGenTot<<endl;
164}
165
Note: See TracBrowser for help on using the repository browser.