ATLAS Offline Software
Loading...
Searching...
No Matches
EFTrackingSmearingAlg.cxx
Go to the documentation of this file.
1/*
2 Copyright (C) 2002-2024 CERN for the benefit of the ATLAS collaboration
3*/
4
5
9
10#include "TH1.h"
11
12#ifndef XAOD_ANALYSIS
15#else
16
17namespace Amg
18{
19 // expression template to evaluate the required size of the compressed matrix at compile time
20 template<int N>
22 static const int value = N + CalculateCompressedSize<N - 1>::value;
23 };
24
25 template<>
26 struct CalculateCompressedSize<1> {
27 static const int value = 1;
28 };
29
30 template<int N>
31 inline void compress(const AmgSymMatrix(N)& covMatrix, std::vector<float>& vec ) {
33 for (unsigned int i = 0; i < N ; ++i)
34 for (unsigned int j = 0; j <= i; ++j)
35 vec.push_back(covMatrix(i,j));
36 }
37 inline void compress(const MatrixX& covMatrix, std::vector<float>& vec) {
38 int rows = covMatrix.rows();
39 for (int i = 0; i < rows; ++i) {
40 for (int j = 0; j <= i; ++j) {
41 vec.push_back(covMatrix(i, j));
42 }
43 }
44 }
45}
46
47namespace Trk
48{
49 enum ParamDefs {
50 // Enums for LocalParameters - LocalPosition/
51 loc1=0, loc2=1,
52 locX=0, locY=1 ,
53 locRPhi=0, locPhiR=0, locZ=1,
54 locR=0, locPhi=1,
55 iPhi=0, iEta=1,
56 distPhi=0, distEta=1,
57 driftRadius=0,
58 // Enums for const Amg::Vector3D & GlobalMomentum /
59 x=0, y=1, z=2,
60 px=0, py=1, pz=2,
61 // Enums for Perigee //
62 d0=0, z0=1, phi0=2, theta=3, qOverP=4,
63 /* Enums for TrackState on Surfaces
64 The first two enums in the TrackParameters refer to the local Frame, i.e.
65 - LocalCartesian for AtanArbitraryPlane
66 - LocalCylindrical for AtaCylinder (includes line)
67 - LocalPolar for AtaDisc
68 The other three enums are standard \f$(\phi, \eta, \frac{q}{p_{T}})\f$
69 */
70 phi=2,
72 u=0, v=1,
73
75 trkMass=5
76
77 };
78
79}
80#endif
81
82
83// from Tracking/TrkEventCnv/TrkEventTPCnv/TrkEventTPCnv/helpers/EigenHelpers.h
84namespace EigenHelpers
85{
86 template <class VECTOR, class COVARIANCE>
87 inline static void eigenMatrixToVector(VECTOR& vec, COVARIANCE& cov, const char* ) {
88 Amg::compress(cov, vec);
89 }
90}
91
92
93
94EFTrackingSmearingAlg::EFTrackingSmearingAlg( const std::string& name, ISvcLocator* pSvcLocator )
95: AthHistogramAlgorithm( name, pSvcLocator ){}
96
97
99
100
102 ATH_MSG_INFO ("Initializing " << name() << "...");
103
104 ATH_CHECK( m_inputTrackParticleKey.initialize() );
105 ATH_CHECK( m_outputTrackParticleKey.initialize() );
106
107 ATH_CHECK( m_inputTruthParticleKey.initialize() );
108 ATH_CHECK( m_outputTruthParticleKey.initialize() );
109
110 // Decoration keys
114
115 ATH_CHECK(m_d0DecoratorKey.initialize());
116 ATH_CHECK(m_z0DecoratorKey.initialize());
117 ATH_CHECK(m_ptDecoratorKey.initialize());
118
119 ATH_MSG_INFO("########## EFTrackingSmearingAlg Configurations are ########## ");
120 ATH_MSG_INFO("------- InputTrackParticleKey: "<<m_inputTrackParticleKey.key());
121 ATH_MSG_INFO("------- OutputTrackParticleKey: "<<m_outputTrackParticleKey.key());
122 ATH_MSG_INFO("------- InputTruthParticleKey: " <<m_inputTruthParticleKey.key());
123 ATH_MSG_INFO("------- OutputTruthParticleKey: "<<m_outputTruthParticleKey.key());
124 ATH_MSG_INFO("------- inputTracksPtCut [GeV]: "<<m_inputTracksPtCut);
125 ATH_MSG_INFO("------- outputTracksPtCut [GeV]: "<<m_outputTracksPtCut);
126 ATH_MSG_INFO("------- SmearingSigma: "<<m_SigmaScaleFactor);
127 ATH_MSG_INFO("------- trackEfficiency: "<<m_smearedTrackEfficiency);
128 ATH_MSG_INFO("------- parameterizedEfficiency: "<<m_parameterizedTrackEfficiency);
129 ATH_MSG_INFO("------- parameterizedEfficiency LRT: "<<m_parameterizedTrackEfficiency_LRT);
130 ATH_MSG_INFO("------- parameterizedEfficiency LRT high d0 cut: "<<m_smearedTrackEfficiency_d0high_LRT );
131 ATH_MSG_INFO("------- parameterizedEfficiency LRT low d0 cut: "<<m_smearedTrackEfficiency_d0low_LRT );
132 ATH_MSG_INFO("------- UseResolutionPtCutOff: "<<m_UseResolutionPtCutOff);
133 ATH_MSG_INFO("------- SetResolutionPtCutOff: "<<m_SetResolutionPtCutOff);
134 ATH_MSG_INFO("------- EnableMonitoring:" <<m_enableMonitoring);
135 ATH_MSG_INFO("------- SmearTruthParticle:"<< m_smearTruthParticle);
136
137 ATH_MSG_INFO("------- IncludeDuplicatesAndFakes: "<<m_EnableFakes);
138 ATH_MSG_INFO("------- RandomSeed: " <<m_RandomSeed);
139 ATH_MSG_INFO("------- UseCoinToss: " <<m_UseCoinToss);
140 ATH_MSG_INFO("------- FakeKillerEnable: "<<m_FakeKillerEnable);
141 ATH_MSG_INFO("------- IncludeFakesInResolutionCalculation: "<<m_IncludeFakesInResolutionCalculation);
142 ATH_MSG_INFO("########## EFTrackingSmearingAlg Configurations: That's it. ########## ");
143
144 std::string smearerName;
146 ATH_MSG_INFO("Will output new truth container with name "<<m_outputTruthParticleKey.key());
147 smearerName = m_outputTruthParticleKey.key()+"_smearer";
148 }
149 else {
150 ATH_MSG_INFO("Will output new track container with name "<<m_outputTrackParticleKey.key());
151 smearerName = m_outputTrackParticleKey.key()+"_smearer";
152 }
153
154 // configure the Smearer
155 m_mySmearer = std::make_unique<FakeTrackSmearer>(smearerName, m_RandomSeed, msgLvl (MSG::DEBUG));
156 m_mySmearer->SetInputTracksPtCut(m_inputTracksPtCut);
157 m_mySmearer->SetOutputTracksPtCut(m_outputTracksPtCut);
158 m_mySmearer->SetTrackingEfficiency(m_smearedTrackEfficiency);
159 m_mySmearer->SetParameterizedEfficiency(m_parameterizedTrackEfficiency);
160 m_mySmearer->SetParameterizedEfficiency_LRT(m_parameterizedTrackEfficiency_LRT);
161 m_mySmearer->SetParameterizedEfficiency_highd0_LRT(m_smearedTrackEfficiency_d0high_LRT);
162 m_mySmearer->SetParameterizedEfficiency_lowd0_LRT(m_smearedTrackEfficiency_d0low_LRT);
163
164 m_mySmearer->SetSigmaScaleFactor(m_SigmaScaleFactor.value());
165 m_mySmearer->UseResolutionPtCutOff(m_UseResolutionPtCutOff.value());
166 m_mySmearer->SetResolutionPtCutOff(m_SetResolutionPtCutOff.value());
167
168 m_mySmearer->EnableFakes(m_EnableFakes.value());
169 m_mySmearer->UseCoinToss(m_UseCoinToss.value());
170 m_mySmearer->FakeKillerEnable(m_FakeKillerEnable.value());
171 m_mySmearer->IncludeFakesInResolutionCalculation(m_IncludeFakesInResolutionCalculation.value());
172
173 if (m_enableMonitoring) {
174 // store the smearing functions
175 TF1 *d0res_eta = m_mySmearer->d0res_eta;
176 TF1 *z0res_eta = m_mySmearer->z0res_eta;
177 TF1 *curvres_eta = m_mySmearer->curvres_eta;
178 TF1 *d0res_pt = m_mySmearer->d0res_pt;
179 TF1 *z0res_pt = m_mySmearer->z0res_pt;
180 TF1 *curvres_pt = m_mySmearer->curvres_pt;
181 TF1 *effLRT_d0 = m_mySmearer->effLRT_d0;
182
183 CHECK(book(new TH1F("d0res_function_vs_eta","#eta of track (p_{T}=10GeV);#eta",100, 0.0,4.0)));
184 CHECK(book(new TH1F("z0res_function_vs_eta","#eta of track (p_{T}=10GeV);#eta",100, 0.0,4.0)));
185 CHECK(book(new TH1F("curvres_function_vs_eta","#eta of track (p_{T}=10GeV);#eta",100, 0.0,4.0)));
186 CHECK(book(new TH1F("d0res_function_vs_pt","p_{T} of track (#eta=1);p_{T} [GeV]",100, 1.0,200.0)));
187 CHECK(book(new TH1F("z0res_function_vs_pt","p_{T} of track (#eta=1);p_{T} [GeV]",100, 1.0,200.0)));
188 CHECK(book(new TH1F("curvres_function_vs_pt","p_{T} of track (#eta=1);p_{T} [GeV]",100, 1.0,200.0)));
189 CHECK(book(new TH1F("effLRT_function_vs_d0","d_{0} of track;d_{0} [mm]",100, 0.001,600.0)));
190 hist("d0res_function_vs_eta")->Add(d0res_eta);
191 hist("z0res_function_vs_eta")->Add(z0res_eta);
192 hist("curvres_function_vs_eta")->Add(curvres_eta);
193 hist("d0res_function_vs_pt")->Add(d0res_pt);
194 hist("z0res_function_vs_pt")->Add(z0res_pt);
195 hist("curvres_function_vs_pt")->Add(curvres_pt);
196 hist("effLRT_function_vs_d0")->Add(effLRT_d0);
197 // book historgams
199 }
200 return StatusCode::SUCCESS;
201}
202
203
204StatusCode EFTrackingSmearingAlg::smearTruthParticles(const EventContext& ctx) {
205
207 const xAOD::TruthParticleContainer* inputTruth = inputTruth_handle.cptr();
208 if (not inputTruth) {
209 ATH_MSG_FATAL("Unable to retrieve input truth particle");
210 return StatusCode::FAILURE;
211 }
213 ATH_CHECK( outputTruth_handle.record( std::make_unique<xAOD::TruthParticleContainer>(), std::make_unique<xAOD::TruthParticleAuxContainer>() ) );
214 auto outputTruth = outputTruth_handle.ptr();
215
216 // create decorators
220
221
222 // clear the smearear
223 m_mySmearer->Clear();
224
225 static const SG::ConstAccessor<float> ptAcc("pt");
226 static const SG::ConstAccessor<float> d0Acc("d0");
227 static const SG::ConstAccessor<float> z0Acc("z0");
228 static const SG::ConstAccessor<float> thetaAcc("theta");
229
230 int n_input_tracks=0;
231 int n_output_tracks=0;
232 int n_output_broad_tracks=0;
233 int n_output_narrow_tracks=0;
234 ATH_MSG_DEBUG ("Found "<<inputTruth->size()<< " input truth particles");
235 for ( const auto* part : *inputTruth )
236 {
237 double pt = part->pt();// MeV
238 float theta = thetaAcc(*part);
239 float z0 = z0Acc(*part);
240 float d0 = d0Acc(*part);
241 float eta = part->eta();
242 float phi = part->phi();
243 if (part->isNeutral()) continue;
244 if (pt <=0.) continue;
245
246 ATH_MSG_DEBUG ("===> New Truth: "
247 <<" curv=" << 1./part->pt()
248 <<" phi=" << part->phi()
249 <<" eta=" << part->eta()
250 <<" d0=" << d0Acc(*part)
251 <<" z0=" << z0Acc(*part)
252 <<" pT=" << part->pt()
253 <<" PDGID=" << part->pdgId()
254 <<" status=" << part->status()
255 );
256 if (part->parent(0)) ATH_MSG_DEBUG (" parent status=" << part->parent(0)->pdgId());
257
258 if (std::abs(pt)/1000. > m_inputTracksPtCut) //GeV cut
259 {
260 n_input_tracks++;
261 if (m_enableMonitoring) {
262 hist("track_input_eta")->Fill(eta);
263 hist("track_input_theta")->Fill(theta);
264 hist("track_input_pt" )->Fill(pt/1000.);
265 hist("track_input_phi")->Fill(phi);
266 hist("track_input_z0" )->Fill(z0);
267 hist("track_input_d0" )->Fill(d0);
268 }
269 double qoverPt = part->charge()*1000./pt; //this must be in GeV
270 m_mySmearer->AddTrack(d0,z0,qoverPt,eta,phi); // smearing here
271 n_output_tracks += m_mySmearer->GetNTracks();
272
273 ATH_MSG_DEBUG ("Looping on output tracks #"<< m_mySmearer->GetNTracks());
274 for (const auto& otrack : m_mySmearer->Tracks)
275 {
276 xAOD::TruthParticle * newtrk = new xAOD::TruthParticle(*part);
277 outputTruth->push_back(newtrk);
278 *newtrk = *part;
279 auto newpt = part->pt();
280 if (m_SigmaScaleFactor !=0) {
281 // set the decorators
282 d0Decorator(*newtrk) = otrack.d0();
283 z0Decorator(*newtrk) = otrack.z0();
284 ptDecorator(*newtrk) = otrack.pt()*1000.; //MeV
285 //TrackParticle has already ::pt(), so the smeared value is in the decorator
286 // and can be accessed by
287 newpt = ptAcc(*newtrk);
288 }
289 if (newpt==0.) continue;
290 ATH_MSG_DEBUG ("Smeared Truth: "
291 <<" curv=" << 1./newpt
292 <<" phi=" << newtrk->phi()
293 <<" eta=" << newtrk->eta()
294 <<" d0=" << d0Acc(*newtrk)
295 <<" z0=" << z0Acc(*newtrk)
296 <<" pT=" << newpt
297 <<" PDGID=" << newtrk->pdgId()
298 <<" status=" << newtrk->status()
299 );
300 if (newtrk->parent(0)) ATH_MSG_DEBUG (" parent status=" << newtrk->parent(0)->pdgId());
301
302 if (m_enableMonitoring) {
303 hist("track_output_eta")->Fill(otrack.eta());
304 hist("track_output_theta")->Fill(otrack.theta());
305 hist("track_output_pt" )->Fill(part->charge()*otrack.pt() );
306 hist("track_output_phi")->Fill(otrack.phi());
307 hist("track_output_z0" )->Fill(otrack.z0() );
308 hist("track_output_d0" )->Fill(otrack.d0() );
309
310 hist("track_outputcoll_eta")->Fill(newtrk->eta());
311 hist("track_outputcoll_theta")->Fill(thetaAcc(*newtrk));
312 hist("track_outputcoll_pt" )->Fill(part->charge()* newpt/1000.);
313 hist("track_outputcoll_phi")->Fill(newtrk->phi());
314 hist("track_outputcoll_z0" )->Fill(z0Acc(*newtrk));
315 hist("track_outputcoll_d0" )->Fill(d0Acc(*newtrk));
316
317 hist("track_delta_eta")->Fill(newtrk->eta() - part->eta());
318 hist("track_delta_pt") ->Fill((newpt - part->pt())/1000.);
319 hist("track_delta_crv")->Fill(newtrk->charge()*1000./newpt - ((part->charge()*1000./part->pt())));
320 hist("track_delta_phi")->Fill(newtrk->phi() - part->phi());
321 hist("track_delta_z0" )->Fill(z0Acc(*newtrk) - z0Acc(*part));
322 hist("track_delta_d0" )->Fill(d0Acc(*newtrk) - d0Acc(*part));
323 }
324 } // end of loop
325 }
326 m_mySmearer->Clear(); // clear the smearer after each input track
327 }
328
329 ATH_MSG_DEBUG ("End of loop track #"<<n_input_tracks<<" ---> "<<" "<< n_output_tracks
330 <<" "<<n_output_narrow_tracks<<" "<<n_output_broad_tracks);
331 if (m_enableMonitoring) {
332 hist("n_input_tracks")->Fill(n_input_tracks);
333 hist("n_output_tracks")->Fill(n_output_tracks);
334 hist("n_output_narrow_tracks")->Fill(n_output_narrow_tracks);
335 hist("n_output_broad_tracks")->Fill(n_output_broad_tracks);
336 }
337 return StatusCode::SUCCESS;
338}
339
340
341
342StatusCode EFTrackingSmearingAlg::execute(const EventContext& ctx) {
343 ATH_MSG_DEBUG ("Executing " << name() << "...");
344
346 return smearTruthParticles(ctx);
347
349 const xAOD::TrackParticleContainer* inputTracks = inputTracks_handle.cptr();
350 if (not inputTracks) {
351 ATH_MSG_FATAL("Unable to retrieve input ID tacks");
352 return StatusCode::FAILURE;
353 }
354
356 ATH_CHECK( outputTracks_handle.record( std::make_unique<xAOD::TrackParticleContainer>(), std::make_unique<xAOD::TrackParticleAuxContainer>() ) );
357 auto outputTracks = outputTracks_handle.ptr();
358
359 // clear the smearear
360 m_mySmearer->Clear();
361
362 //int trackno=0;
363 int n_input_tracks=0;
364 int n_output_tracks=0;
365 int n_output_broad_tracks=0;
366 int n_output_narrow_tracks=0;
367 ATH_MSG_DEBUG ("Found "<<inputTracks->size()<< " input tracks");
368 for ( const auto* trk : *inputTracks )
369 {
370 // get Cov matrix of input track
371 xAOD::ParametersCovMatrix_t trkcov = trk->definingParametersCovMatrix();
372 auto trkcovvec = trk->definingParametersCovMatrixVec();
373 double theta=trk->theta();
374 double pt = std::sin(theta)/std::abs(trk->qOverP());// MeV
375
376 ATH_MSG_DEBUG ("===> New Track: "
377 <<" curv=" << 1./trk->pt()
378 <<" phi=" << trk->phi0()
379 <<" eta=" << trk->eta()
380 <<" d0=" << trk->d0()
381 <<" z0=" << trk->z0()
382 <<" pT=" << trk->pt()
383 <<" cov_d0=" << trkcov(Trk::d0,Trk::d0)
384 <<" cov_z0=" << trkcov(Trk::z0,Trk::z0)
385 <<" sigma_d0=" << std::sqrt(std::abs(trkcov(Trk::d0,Trk::d0)))
386 <<" sigma_z0=" << std::sqrt(std::abs(trkcov(Trk::z0,Trk::z0))) );
387
388
389
390 if (std::abs(pt) > m_inputTracksPtCut)
391 {
392 n_input_tracks++;
393 if (m_enableMonitoring) {
394 hist("track_input_eta")->Fill(trk->eta());
395 hist("track_input_theta")->Fill(trk->theta());
396 hist("track_input_pt" )->Fill(pt/1000.);
397 hist("track_input_phi")->Fill(trk->phi0());
398 hist("track_input_z0" )->Fill(trk->z0());
399 hist("track_input_d0" )->Fill(trk->d0());
400
401 hist("track_input_sigma_theta") ->Fill(std::sqrt(trkcov(Trk::theta,Trk::theta)));
402 hist("track_input_sigma_qOverP")->Fill(std::sqrt(trkcov(Trk::qOverP,Trk::qOverP)));
403 hist("track_input_sigma_phi") ->Fill(std::sqrt(trkcov(Trk::phi0,Trk::phi0)));
404 hist("track_input_sigma_z0") ->Fill(std::sqrt(trkcov(Trk::z0,Trk::z0)));
405 hist("track_input_sigma_d0") ->Fill(std::sqrt(trkcov(Trk::d0,Trk::d0)));
406 }
407
408 // get Cov matrix of input track
409 auto trkcovvec = trk->definingParametersCovMatrixVec();
410 double qoverPt = trk->charge()*1000./pt; //this must be in GeV
411 m_mySmearer->AddTrack(trk->d0(),trk->z0(),qoverPt,trk->eta(),trk->phi0());
412 n_output_tracks += m_mySmearer->GetNTracks();
413
414 ATH_MSG_DEBUG ("Looping on output tracks #"<< m_mySmearer->GetNTracks());
415 for (const auto& otrack : m_mySmearer->Tracks)
416 {
417 xAOD::TrackParticle * newtrk = new xAOD::TrackParticle(*trk);
418 outputTracks->push_back(newtrk);
419 *newtrk = *trk;
420
421 double sintheta=std::sin(otrack.theta());
422 trkcov = trk->definingParametersCovMatrix();
423 auto newtrkcov = trkcov;
424 if (m_SigmaScaleFactor !=0) {
425 // modify the cov matrix
426 ATH_MSG_DEBUG ("Setting parameters in covariance vector");
427 for (int ii=0;ii<5;ii++) for (int jj=0;jj<5;jj++) {
428 // correct the CM to be consistent with the smearing ofthe parameters
429 newtrkcov(ii,jj)=trkcov(ii,jj) * (std::pow(m_SigmaScaleFactor,2) + 1);
430 }
431 trkcovvec.clear();
432 EigenHelpers::eigenMatrixToVector(trkcovvec,newtrkcov,"");
433 newtrk->setDefiningParametersCovMatrixVec(trkcovvec);
434 ATH_MSG_DEBUG ("Setting parameters covariance");
435 newtrk->setDefiningParametersCovMatrix(newtrkcov);
436
437 /* old method stays for reference, since it's not yet validated
438 for (int ii=0;ii<5;ii++) for (int jj=0;jj<5;jj++) if (ii!=jj) trkcov(ii,jj)=0.0;
439 trkcov(Trk::d0 ,Trk::d0 ) = otrack.sigma_d0() * otrack.sigma_d0() ;
440 trkcov(Trk::z0 ,Trk::z0 ) = otrack.sigma_z0() * otrack.sigma_z0() ;
441 trkcov(Trk::theta ,Trk::theta ) = otrack.sigma_theta()* otrack.sigma_theta();
442 trkcov(Trk::phi0 ,Trk::phi0 ) = otrack.sigma_phi() * otrack.sigma_phi() ;
443 trkcov(Trk::qOverP,Trk::qOverP) = otrack.sigma_curv()/1000. * otrack.sigma_curv()/1000. *sintheta*sintheta;
444
445 trkcovvec.clear();
446 EigenHelpers::eigenMatrixToVector(trkcovvec,trkcov,"");
447 newtrk->setDefiningParametersCovMatrixVec(trkcovvec);
448 ATH_MSG_DEBUG ("Setting parameters covariance");
449 newtrk->setDefiningParametersCovMatrix(trkcov);
450 */
451
452 //(float d0, float z0, float phi0, float theta, float qOverP)
453 newtrk->setDefiningParameters(
454 otrack.d0(),
455 otrack.z0(),
456 otrack.phi(),
457 otrack.theta(),
458 (otrack.curv()*sintheta/1000.) //oneOverp in MeV
459 );
460 }
461 xAOD::ParametersCovMatrix_t trkcov_out = newtrk->definingParametersCovMatrix();
462
463
464 ATH_MSG_DEBUG ("Smeared Track: "
465 <<" curv=" << 1./newtrk->pt()
466 <<" phi=" << newtrk->phi()
467 <<" eta=" << newtrk->eta()
468 <<" d0=" << newtrk->d0()
469 <<" z0=" << newtrk->z0()
470 <<" pT=" << newtrk->pt()
471 <<" cov_d0=" << trkcov_out(Trk::d0,Trk::d0)
472 <<" cov_z0=" << trkcov_out(Trk::z0,Trk::z0)
473 <<" sigma_d0=" << std::sqrt(std::abs(trkcov_out(Trk::d0,Trk::d0)))
474 <<" sigma_z0=" << std::sqrt(std::abs(trkcov_out(Trk::z0,Trk::z0))) );
475
476
477 if (m_enableMonitoring) {
478 hist("track_output_eta")->Fill(otrack.eta());
479 hist("track_output_theta")->Fill(otrack.theta());
480 hist("track_output_pt" )->Fill(trk->charge()*otrack.pt() );
481 hist("track_output_phi")->Fill(otrack.phi());
482 hist("track_output_z0" )->Fill(otrack.z0() );
483 hist("track_output_d0" )->Fill(otrack.d0() );
484
485 hist("track_outputcoll_eta")->Fill(newtrk->eta());
486 hist("track_outputcoll_theta")->Fill(newtrk->theta());
487 hist("track_outputcoll_pt" )->Fill(trk->charge()*std::sin(newtrk->theta())/(1000.0*newtrk->qOverP()));
488 hist("track_outputcoll_phi")->Fill(newtrk->phi0());
489 hist("track_outputcoll_z0" )->Fill(newtrk->z0());
490 hist("track_outputcoll_d0" )->Fill(newtrk->d0());
491
492 hist("track_delta_eta")->Fill(newtrk->eta() - trk->eta());
493 hist("track_delta_pt")->Fill((newtrk->pt() - trk->pt())/1000.);
494 hist("track_delta_crv")->Fill((1000.0*newtrk->qOverP()/std::sin(theta))-((1000.0*trk->qOverP())/std::sin(theta)));
495 hist("track_delta_phi")->Fill(newtrk->phi() - trk->phi0());
496 hist("track_delta_z0" )->Fill(newtrk->z0() - trk->z0());
497 hist("track_delta_d0" )->Fill(newtrk->d0() - trk->d0());
498
499 auto trkcov_original = trk->definingParametersCovMatrix();
500 hist("track_delta_sigma_theta") ->Fill(std::sqrt(trkcov_out(Trk::theta,Trk::theta)) - std::sqrt(trkcov_original(Trk::theta,Trk::theta)));
501 hist("track_delta_sigma_qOverP")->Fill(std::sqrt(trkcov_out(Trk::qOverP,Trk::qOverP)) - std::sqrt(trkcov_original(Trk::qOverP,Trk::qOverP)));
502 hist("track_delta_sigma_phi") ->Fill(std::sqrt(trkcov_out(Trk::phi0,Trk::phi0)) - std::sqrt(trkcov_original(Trk::phi0,Trk::phi0)));
503 hist("track_delta_sigma_z0") ->Fill(std::sqrt(trkcov_out(Trk::z0,Trk::z0)) - std::sqrt(trkcov_original(Trk::z0,Trk::z0)));
504 hist("track_delta_sigma_d0") ->Fill(std::sqrt(trkcov_out(Trk::d0,Trk::d0)) - std::sqrt(trkcov_original(Trk::d0,Trk::d0)));
505
506
507 hist("track_outputcoll_sigma_theta") ->Fill(std::sqrt(trkcov(Trk::theta,Trk::theta)));
508 hist("track_outputcoll_sigma_qOverP")->Fill(std::sqrt(trkcov(Trk::qOverP,Trk::qOverP)));
509 hist("track_outputcoll_sigma_phi") ->Fill(std::sqrt(trkcov(Trk::phi0,Trk::phi0)));
510 hist("track_outputcoll_sigma_z0") ->Fill(std::sqrt(trkcov(Trk::z0,Trk::z0)));
511 hist("track_outputcoll_sigma_d0") ->Fill(std::sqrt(trkcov(Trk::d0,Trk::d0)));
512
513 }
514 // do we need this hack?
515 // hack!!!
516 // Chi2: 1 -> "Broad Fake"
517 // 0 -> "core tracks"
518 //
519 // NDF: index of parent track in input track collection
520 //newtrk->setFitQuality(otrack.FakeFlags(),trackno);
521
522
523 }
524
525 }
526 m_mySmearer->Clear(); // clear teh smearer after each input track
527 //trackno++;
528 }
529
530 ATH_MSG_DEBUG ("End of loop track #"<<n_input_tracks<<" ---> "<<" "<< n_output_tracks
531 <<" "<<n_output_narrow_tracks<<" "<<n_output_broad_tracks);
532 if (m_enableMonitoring) {
533 hist("n_input_tracks")->Fill(n_input_tracks);
534 hist("n_output_tracks")->Fill(n_output_tracks);
535 hist("n_output_narrow_tracks")->Fill(n_output_narrow_tracks);
536 hist("n_output_broad_tracks")->Fill(n_output_broad_tracks);
537 }
538 return StatusCode::SUCCESS;
539}
540
541
542
543
545{
546 double b_eta[3]={50,-5,5};
547 CHECK(book(new TH1F("track_input_eta","#eta of input tracks",b_eta[0],b_eta[1],b_eta[2])));
548 CHECK(book(new TH1F("track_input_theta","#Theta of input tracks",50,0.0,4.0)));
549 CHECK(book(new TH1F("track_input_pt" ,"p_T of input tracks",50.0,-10.0,10.0)));
550 CHECK(book(new TH1F("track_input_phi","#phi of input tracks",50,-6.28,6.28)));
551 CHECK(book(new TH1F("track_input_z0" ,"z_0 of input tracks",100,-50,50)));
552 CHECK(book(new TH1F("track_input_d0" ,"d_0 of input tracks",50,-5,5)));
553
554 CHECK(book(new TH1F("track_input_sigma_theta","#sigma_{#Theta} of input tracks" ,50,0.0,0.001)));
555 CHECK(book(new TH1F("track_input_sigma_qOverP" ,"#sigma_{q/P} of input tracks" ,50.0,0.0,0.0001)));
556 CHECK(book(new TH1F("track_input_sigma_phi" ,"#sigma_{#phi} of input tracks" ,50,0.0,0.002)));
557 CHECK(book(new TH1F("track_input_sigma_z0" ,"#sigma_{z_0} of input tracks" ,100,0.0,5)));
558 CHECK(book(new TH1F("track_input_sigma_d0" ,"#sigma_{d_0} of input tracks" ,50,0.0,5)));
559
560 int maxtracks=500;
561 CHECK(book(new TH1F("n_input_tracks" ,"Number of input tracks" ,maxtracks,0,maxtracks)));
562 CHECK(book(new TH1F("n_output_tracks" ,"Number of output tracks" ,maxtracks,0,maxtracks)));
563 CHECK(book(new TH1F("n_output_narrow_tracks","Number of output tracks (narrow)" ,maxtracks,0,maxtracks)));
564 CHECK(book(new TH1F("n_output_broad_tracks" ,"Number of output tracks (broad)" ,maxtracks,0,maxtracks)));
565
566 CHECK(book(new TH1F("track_output_eta","#eta of output tracks",50,-5,5)));
567 CHECK(book(new TH1F("track_output_theta","#Theta of output tracks",50,0.0,4.0)));
568 CHECK(book(new TH1F("track_output_pt" ,"p_T of output tracks [GeV]",50,-10,10)));
569 CHECK(book(new TH1F("track_output_phi","#phi of output tracks",50,-6.28,6.28)));
570 CHECK(book(new TH1F("track_output_z0" ,"z_0 of output tracks",100,-50,50)));
571 CHECK(book(new TH1F("track_output_d0" ,"d_0 of output tracks",50,-5,5)));
572
573 CHECK(book(new TH1F("track_outputcoll_eta","#eta of output tracks collection",50,-5,5)));
574 CHECK(book(new TH1F("track_outputcoll_theta","#Theta of output tracks collection",50,0.0,4.0)));
575 CHECK(book(new TH1F("track_outputcoll_pt" ,"p_T of output tracks collection [GeV]",50,-10,10)));
576 CHECK(book(new TH1F("track_outputcoll_phi","#phi of output tracks collection",50,-6.28,6.28)));
577 CHECK(book(new TH1F("track_outputcoll_z0" ,"z_0 of output tracks collection",100,-50,50)));
578 CHECK(book(new TH1F("track_outputcoll_d0" ,"d_0 of output tracks collection",50,-5,5)));
579
580 CHECK(book(new TH1F("track_outputcoll_sigma_theta","#sigma_{#Theta} of output tracks collection" ,50,0.0,0.001)));
581 CHECK(book(new TH1F("track_outputcoll_sigma_qOverP" ,"#sigma_{q/P} of output tracks collection",50.0,0.0,0.0001)));
582 CHECK(book(new TH1F("track_outputcoll_sigma_phi" ,"#sigma_{#phi} of output tracks collection" ,50,0.0,0.002)));
583 CHECK(book(new TH1F("track_outputcoll_sigma_z0" ,"#sigma_{z_0} of output tracks collection" ,100,0.0,5)));
584 CHECK(book(new TH1F("track_outputcoll_sigma_d0" ,"#sigma_{d_0} of output tracks collection" ,50,0.0,5)));
585
586 CHECK(book(new TH1F("track_delta_sigma_theta", "#sigma_{#Theta} of output tracks collection",100,-0.001,0.001)));
587 CHECK(book(new TH1F("track_delta_sigma_qOverP", "#sigma_{q/P} of output tracks collection",100,-0.0001,0.0001)));
588 CHECK(book(new TH1F("track_delta_sigma_phi", "#sigma_{#phi} of output tracks collection" ,100,-0.002,0.002)));
589 CHECK(book(new TH1F("track_delta_sigma_z0", "#sigma_{z_0} of output tracks collection" ,100,-1.,1.)));
590 CHECK(book(new TH1F("track_delta_sigma_d0", "#sigma_{d_0} of output tracks collection" ,100,-0.5,0.5)));
591
592 CHECK(book(new TH1F("track_delta_eta","tracks #Delta #eta",50,-1,1)));
593 CHECK(book(new TH1F("track_delta_pt","tracks #Delta #pt [GeV]",50,-2.,2.)));
594 CHECK(book(new TH1F("track_delta_crv" ,"tracks #Delta crv",50,-1,1)));
595 CHECK(book(new TH1F("track_delta_phi","tracks #Delta #phi",50,-0.5,0.5)));
596 CHECK(book(new TH1F("track_delta_z0" ,"tracks #Delta z_0 ",100,-10,10)));
597 CHECK(book(new TH1F("track_delta_d0" ,"tracks #Delta d_0 ",50,-2,2)));
598 return StatusCode::SUCCESS;
599}
600
601
602
Scalar eta() const
pseudorapidity method
Scalar phi() const
phi method
Scalar theta() const
theta method
#define ATH_CHECK
Evaluate an expression and check for errors.
#define ATH_MSG_FATAL(x)
#define ATH_MSG_INFO(x)
#define ATH_MSG_DEBUG(x)
std::vector< size_t > vec
Helper class to provide constant type-safe access to aux data.
#define CHECK(...)
Evaluate an expression and check for errors.
#define AmgSymMatrix(dim)
bool msgLvl(const MSG::Level lvl) const
AthHistogramAlgorithm(const std::string &name, ISvcLocator *pSvcLocator)
Constructor with parameters:
StatusCode book(const TH1 &hist, const std::string &tDir="", const std::string &stream="")
Simplify the booking and registering (into THistSvc) of histograms.
TH1 * hist(const std::string &histName, const std::string &tDir="", const std::string &stream="")
Simplify the retrieval of registered histograms of any type.
size_type size() const noexcept
Returns the number of elements in the collection.
EFTrackingSmearingAlg(const std::string &name, ISvcLocator *pSvcLocator)
BooleanProperty m_IncludeFakesInResolutionCalculation
DoubleProperty m_SetResolutionPtCutOff
DoubleProperty m_smearedTrackEfficiency
virtual StatusCode execute(const EventContext &ctx) override
Execute method.
SG::WriteDecorHandleKey< xAOD::TruthParticleContainer > m_ptDecoratorKey
SG::WriteDecorHandleKey< xAOD::TruthParticleContainer > m_z0DecoratorKey
BooleanProperty m_parameterizedTrackEfficiency
BooleanProperty m_UseResolutionPtCutOff
BooleanProperty m_enableMonitoring
SG::WriteDecorHandleKey< xAOD::TruthParticleContainer > m_d0DecoratorKey
DoubleProperty m_smearedTrackEfficiency_d0low_LRT
SG::ReadHandleKey< xAOD::TrackParticleContainer > m_inputTrackParticleKey
BooleanProperty m_FakeKillerEnable
SG::ReadHandleKey< xAOD::TruthParticleContainer > m_inputTruthParticleKey
virtual StatusCode initialize() override
BooleanProperty m_parameterizedTrackEfficiency_LRT
StatusCode smearTruthParticles(const EventContext &ctx)
std::unique_ptr< FakeTrackSmearer > m_mySmearer
DoubleProperty m_smearedTrackEfficiency_d0high_LRT
BooleanProperty m_smearTruthParticle
SG::WriteHandleKey< xAOD::TrackParticleContainer > m_outputTrackParticleKey
SG::WriteHandleKey< xAOD::TruthParticleContainer > m_outputTruthParticleKey
Helper class to provide constant type-safe access to aux data.
const_pointer_type cptr()
Dereference the pointer.
Handle class for adding a decoration to an object.
StatusCode record(std::unique_ptr< T > data)
Record a const object to the store.
pointer_type ptr()
Dereference the pointer.
void setDefiningParametersCovMatrix(const ParametersCovMatrix_t &cov)
Set the defining parameters covariance matrix.
float z0() const
Returns the parameter.
float theta() const
Returns the parameter, which has range 0 to .
virtual double phi() const override final
The azimuthal angle ( ) of the particle (has range to .).
const ParametersCovMatrix_t definingParametersCovMatrix() const
Returns the 5x5 symmetric matrix containing the defining parameters covariance matrix.
float d0() const
Returns the parameter.
void setDefiningParameters(float d0, float z0, float phi0, float theta, float qOverP)
Set the defining parameters.
float qOverP() const
Returns the parameter.
virtual double pt() const override final
The transverse momentum ( ) of the particle.
virtual double eta() const override final
The pseudorapidity ( ) of the particle.
float phi0() const
Returns the parameter, which has range to .
void setDefiningParametersCovMatrixVec(const std::vector< float > &cov)
int status() const
Status code.
int pdgId() const
PDG ID code.
const TruthParticle_v1 * parent(size_t i) const
Retrieve the i-th mother (TruthParticle) of this TruthParticle.
double charge() const
Physical charge.
virtual double eta() const override final
The pseudorapidity ( ) of the particle.
virtual double phi() const override final
The azimuthal angle ( ) of the particle.
Definition of ATLAS Math & Geometry primitives (Amg).
Eigen::Matrix< double, Eigen::Dynamic, Eigen::Dynamic > MatrixX
Dynamic Matrix - dynamic allocation.
void compress(const AmgSymMatrix(N) &covMatrix, std::vector< float > &vec)
constexpr int CalculateCompressedSize(int n)
static void eigenMatrixToVector(VECTOR &vec, COVARIANCE &cov, const char *)
constexpr std::size_t N
float j(const xAOD::IParticle &, const xAOD::TrackMeasurementValidation &hit, const Eigen::Matrix3d &jab_inv)
Ensure that the ATLAS eigen extensions are properly loaded.
ParamDefs
This file defines the parameter enums in the Trk namespace.
Definition ParamDefs.h:32
@ driftRadius
trt, straws
Definition ParamDefs.h:53
@ iEta
(old readout) will be skipped
Definition ParamDefs.h:48
@ locY
local cartesian
Definition ParamDefs.h:38
@ x
Definition ParamDefs.h:55
@ locX
Definition ParamDefs.h:37
@ z
global position (cartesian)
Definition ParamDefs.h:57
@ pz
global momentum (cartesian)
Definition ParamDefs.h:61
@ iPhi
Definition ParamDefs.h:47
@ locR
Definition ParamDefs.h:44
@ v
Definition ParamDefs.h:78
@ locRPhi
Definition ParamDefs.h:40
@ u
Enums for curvilinear frames.
Definition ParamDefs.h:77
@ distEta
readout for silicon
Definition ParamDefs.h:51
@ phi0
Definition ParamDefs.h:65
@ distPhi
Definition ParamDefs.h:50
@ locPhiR
Definition ParamDefs.h:41
@ theta
Definition ParamDefs.h:66
@ qOverP
perigee
Definition ParamDefs.h:67
@ y
Definition ParamDefs.h:56
@ loc2
generic first and second local coordinate
Definition ParamDefs.h:35
@ phi
Definition ParamDefs.h:75
@ loc1
Definition ParamDefs.h:34
@ locZ
local cylindrical
Definition ParamDefs.h:42
@ d0
Definition ParamDefs.h:63
@ px
Definition ParamDefs.h:59
@ z0
Definition ParamDefs.h:64
@ trkMass
Extended perigee: mass.
Definition ParamDefs.h:81
@ py
Definition ParamDefs.h:60
@ locPhi
local polar
Definition ParamDefs.h:45
TrackParticle_v1 TrackParticle
Reference the current persistent version:
TruthParticle_v1 TruthParticle
Typedef to implementation.
TrackParticleContainer_v1 TrackParticleContainer
Definition of the current "TrackParticle container version".
TruthParticleContainer_v1 TruthParticleContainer
Declare the latest version of the truth particle container.
TVectorD VECTOR