Fork me on GitHub

source: git/modules/HighMassVertexRecover.cc@ 584e2e9

Timing
Last change on this file since 584e2e9 was 0976f6a, checked in by Kaan Yüksel Oyulmaz <kaanyukseloyulmaz@…>, 5 years ago

Card updated and new module added to modules directory as you said also I updated a propagator card for Ozgun to test his analysis including High vertex recover

  • Property mode set to 100755
File size: 8.3 KB
Line 
1/** \class HighMassVertexRecover
2 *
3 * Try to assign a vertex also to tracks which have
4 * not been clusterized changing the mass hypotesis.
5 *
6 * \author Olmo Cerri, Caltech
7 *
8 */
9
10#include "modules/HighMassVertexRecover.h"
11
12#include "classes/DelphesClasses.h"
13#include "classes/DelphesFactory.h"
14#include "classes/DelphesFormula.h"
15
16#include "ExRootAnalysis/ExRootResult.h"
17#include "ExRootAnalysis/ExRootFilter.h"
18#include "ExRootAnalysis/ExRootClassifier.h"
19
20#include "TMath.h"
21#include "TString.h"
22#include "TObjArray.h"
23#include "TLorentzVector.h"
24
25#include <iostream>
26#include <sstream>
27
28
29using namespace std;
30namespace vtx_DAZT
31{
32 static const Double_t c_light = 2.99792458e+8; // [m/s]
33}
34using namespace vtx_DAZT;
35
36//------------------------------------------------------------------------------
37
38HighMassVertexRecover::HighMassVertexRecover()
39{
40}
41
42//------------------------------------------------------------------------------
43
44HighMassVertexRecover::~HighMassVertexRecover()
45{
46}
47
48//------------------------------------------------------------------------------
49
50void HighMassVertexRecover::Init()
51{
52 // Get parameters from the card
53 fVerbose = GetInt("Verbose", 0);
54
55 fSigmaCompatibility = GetDouble("SigmaCompatibility", 2.0);
56
57 ExRootConfParam param = GetParam("MassHypot");
58 UInt_t size = param.GetSize();
59 if(size > 0)
60 {
61 for(UInt_t i = 0; i < size; ++i)
62 {
63 Double_t mass = GetDouble(param[i].GetString(), 1.);
64 fMassList.push_back(mass);
65 }
66 }
67 else
68 {
69 // Ordered by production probability in the SM
70 fMassList.push_back(0.13957); //pion (+/-)
71 fMassList.push_back(0.49367); // k (+/-)
72 fMassList.push_back(0.93827); // proton
73 fMassList.push_back(0.00051); //electron
74 fMassList.push_back(0.10558); //muon
75 }
76
77 if(fVerbose)
78 {
79 for(unsigned int i = 0; i < fMassList.size(); i++)
80 {
81 cout << "mass: " << fMassList[i] << endl;
82 }
83 }
84
85
86 // import input array
87 fTrackInputArray = ImportArray(GetString("TrackInputArray", "VertexFinderDAClusterizerZT/tracks"));
88 fItTrackInputArray = fTrackInputArray->MakeIterator();
89
90 fVertexInputArray = ImportArray(GetString("VertexInputArray", "VertexFinderDAClusterizerZT/vertices"));
91 fItVertexInputArray = fVertexInputArray->MakeIterator();
92
93 // create output array
94 fTrackOutputArray = ExportArray(GetString("TrackOutputArray", "tracks"));
95 fVertexOutputArray = ExportArray(GetString("VertexOutputArray", "vertices"));
96}
97
98//------------------------------------------------------------------------------
99
100void HighMassVertexRecover::Finish()
101{
102 delete fItTrackInputArray;
103 delete fItVertexInputArray;
104 delete fItVertexOutputArray;
105}
106
107//------------------------------------------------------------------------------
108
109void HighMassVertexRecover::Process()
110{
111 if(fVerbose > 1)
112 {
113 cout << endl << endl << endl << "------------------ EVENT ----------------------------------------------" << endl;
114 }
115 // Make a copy of input VerticesfItVertexInputArray->Reset();
116 Candidate * vertex;
117 fItVertexInputArray->Reset();
118 while((vertex = static_cast<Candidate*>(fItVertexInputArray->Next())))
119 {
120 fVertexOutputArray->Add(static_cast<Candidate*>(vertex->Clone()));
121 }
122 fItVertexOutputArray = fVertexOutputArray->MakeIterator();
123
124 fItTrackInputArray->Reset();
125 Candidate * track;
126 while((track = static_cast<Candidate*>(fItTrackInputArray->Next())))
127 {
128 Candidate* mother = track;
129 track = static_cast<Candidate*>(track->Clone());
130 track->AddCandidate(mother);
131
132 if(track->ClusterIndex == -1)
133 {
134 if(fVerbose > 1)
135 {
136 cout << endl << "PID: " << track->PID << endl;
137 cout << "Pt: " << track->Momentum.Pt() << endl;
138 cout << Form("Zd: %.3f +/- %.3f", track->Zd, track->ErrorDZ) << endl;
139 }
140 vector<Candidate*> vtx; //vector on vertices compatible in z
141 // Fill it with onyl the vertexes compatible
142 fItVertexOutputArray->Reset();
143 Candidate * vertex;
144 while((vertex = static_cast<Candidate*>(fItVertexOutputArray->Next())))
145 {
146 double dz = vertex->Position.Z() - track->Zd;
147 dz /= TMath::Hypot(track->ErrorDZ, vertex->PositionError.Z());
148
149 Double_t sigd0 = fabs(track->D0)/track->ErrorD0;
150 if(fabs(dz) < fSigmaCompatibility && sigd0 < 3) //Should be done better...is assuming v_xy = 0 and the stat threatment is poor/wrong
151 {
152 UInt_t i = 0;
153 UInt_t nv = vtx.size();
154 while(i<nv)
155 {
156 if(vtx[i]->SumPT2 > vertex->SumPT2) i++;
157 else break;
158 }
159 vtx.insert(vtx.begin()+i, vertex);
160 }
161 }
162
163 if(fVerbose > 2)
164 {
165 cout << "Compatible vertexes: " << vtx.size() << endl;
166 for(unsigned int i = 0; i<vtx.size(); i++)
167 {
168 cout << Form("SumPT2: %.0f GeV^2 Z: %.3f +/- %.3f mm T: %.3f +/- %.3f ps", vtx[i]->SumPT2, vtx[i]->Position.Z(), vtx[i]->PositionError.Z(), vtx[i]->Position.T(), vtx[i]->PositionError.T()) << endl;
169 }
170 }
171
172 if(vtx.size() > 0)
173 {
174 // Try at first to see if the mass hypothesys can work
175 // The vertexes are order by SumPT2 and the mass hypotesys for relevance
176 UInt_t i = 0, match = 0;
177 while(i<vtx.size() && match == 0)
178 {
179 UInt_t j = 0;
180 while(j<fMassList.size() && match == 0)
181 {
182 auto Tpair = ComputeCATime(track, fMassList[j]);
183 double dt = vtx[i]->Position.T()*1.E9/c_light - Tpair.first; // [ps]
184 // dt /= Tpair.second;
185 dt /= TMath::Hypot(Tpair.second, vtx[i]->PositionError.T()*1.E9/c_light);
186
187 if(fabs(dt)<fSigmaCompatibility)
188 {
189 match = 1;
190 track->Mass = fMassList[j];
191 track->ClusterIndex = vtx[i]->ClusterIndex;
192 track->InitialPosition.SetT(vtx[i]->Position.T());
193 track->InitialPosition.SetZ(vtx[i]->Position.Z());
194 track->Td = Tpair.first * 1E-9 * c_light;
195
196 vtx[i]->SumPT2 += track->Momentum.Pt()*track->Momentum.Pt();
197 vtx[i]->SumPt += track->Momentum.Pt();
198 if(fVerbose>5)
199 {
200 cout << "SM VTX found at " << i << ". dt = " << dt << "sigma" << endl;
201 cout << "vtx=> SumPT2 " << vtx[i]->SumPT2 << " - DOF: " << vtx[i]->ClusterNDF << endl;
202 cout << "Mass: " << fMassList[j] << endl;
203 }
204 }
205
206 j++;
207 }
208
209 i++;
210 }
211
212 if(match == 0)
213 {
214 Double_t p = track->Momentum.Pt() * sqrt(1 + track->CtgTheta*track->CtgTheta);
215 // Only Z
216 // Double_t beta_z = vtx[0]->Position.Z() - track->Position.Z();
217 // beta_z /= vtx[0]->Position.T() - track->Position.T();
218 // Double_t e = track->Momentum.Pt() * track->CtgTheta / beta_z;
219
220 //Full path length
221 Double_t beta = track->L / (vtx[0]->Position.T() - track->Position.T());
222 Double_t e = p / beta;
223
224 track->Mass = sqrt(e*e - p*p);
225 track->ClusterIndex = vtx[0]->ClusterIndex;
226 track->InitialPosition.SetT(vtx[0]->Position.T());
227 track->InitialPosition.SetZ(vtx[0]->Position.Z());
228 track->Td = vtx[0]->Position.T();
229
230 vtx[0]->SumPT2 += track->Momentum.Pt()*track->Momentum.Pt();
231 vtx[0]->SumPt += track->Momentum.Pt();
232 if(fVerbose>5)
233 {
234 cout << "BSM VTX fitted:" << endl;
235 cout << "vtx=> SumPT2 " << vtx[0]->SumPT2 << " - DOF: " << vtx[0]->ClusterNDF << endl;
236 cout << "Mass: " << track->Mass << endl;
237 }
238 }
239
240 }
241 }
242
243 fTrackOutputArray->Add(track);
244 }
245}
246
247//------------------------------------------------------------------------------
248// Auxiliary function to compute the CA time and erro given the mass
249pair<Double_t, Double_t> HighMassVertexRecover::ComputeCATime(Candidate * tk, Double_t m)
250{
251 Double_t p = tk->Momentum.Pt() * sqrt(1 + tk->CtgTheta*tk->CtgTheta);
252 Double_t e = sqrt(p*p + m*m);
253
254 Double_t t = tk->Position.T()*1.E9/c_light; // from [mm] to [ps]
255 //Full path length
256 t -= tk->L*1E9/(c_light*p/e);
257
258 // Only Z
259 // Double_t bz = tk->Momentum.Pt() * tk->CtgTheta/e;
260 // t += (tk->Zd - tk->Position.Z())*1E9/(c_light*bz);
261
262 pair<Double_t, Double_t> out;
263 out.first = t;
264 out.second = tk->ErrorT*1.E9/c_light; // [ps]
265
266 return out;
267}
Note: See TracBrowser for help on using the repository browser.