001
002
003
004
005
006 #include "GaudiKernel/SmartDataPtr.h"
007 #include "GeneratorObjects/McEventCollection.h"
008 #include "InDetAlignTrkInfo/AlignTrkContainer.h"
009 #include "TrkTrack/Track.h"
010 #include "TrkTrack/TrackCollection.h"
011 #include "TrkTruthData/TrackTruthCollection.h"
012 #include "TrkMaterialOnTrack/MaterialEffectsOnTrack.h"
013 #include "TrkMaterialOnTrack/ScatteringAngles.h"
014 #include "TrkParameters/MeasuredPerigee.h"
015 #include "TrkRIO_OnTrack/RIO_OnTrack.h"
016 #include "InDetRIO_OnTrack/SiClusterOnTrack.h"
017 #include "InDetRIO_OnTrack/TRT_DriftCircleOnTrack.h"
018 #include "TrkEventPrimitives/ParamDefs.h"
019 #include "TrkExInterfaces/IExtrapolator.h"
020 #include "TrkEventPrimitives/PropDirection.h"
021 #include "TrkPrepRawData/PrepRawData.h"
022 #include "InDetPrepRawData/TRT_DriftCircle.h"
023 #include "InDetIdentifier/PixelID.h"
024 #include "InDetIdentifier/SCT_ID.h"
025 #include "InDetIdentifier/TRT_ID.h"
026 #include "InDetAlignGenAlgs/AlignSiModuleList.h"
027 #include "InDetAlignGenTools/InDetAlignTrackSelTool.h"
028 #include "InDetAlignGenAlgs/InDetAlignNt.h"
029 #include "TrkPseudoMeasurementOnTrack/PseudoMeasurementOnTrack.h"
030
031
032
033
034
035 #define MAXHITS 60
036
037 void InDetAlignNt::FillFromTrkTrack(AlignTrkContainer* tracks, int& ntrk) {
038
039 const DataVector<Trk::Track>* trks;
040 if (StatusCode::SUCCESS!=evtStore()->retrieve(trks,par_newtrkcol)) {
041 msg (MSG::ERROR) << "Cannot find " << par_newtrkcol << endreq;
042
043 } else {
044 ATH_MSG_DEBUG ( "New track filling called found " << trks->size() );
045
046
047
048
049 for (DataVector<Trk::Track>::const_iterator t=trks->begin();
050 t!=trks->end();++t) {
051 ++ntrk;
052 bool status=false;
053 if (iTrkTrkSel(*t)) {
054 AlignTrk* atp=new AlignTrk(TrkTrkProc(trks,*t,ntrk,status));
055 if (status) {
056 if (par_alignnt) FillTrkNtuple(*atp);
057 if (par_overlap) {
058 OverlapProc(*atp);
059 FillOverlap(*atp);
060 }
061 tracks->push_back(atp);
062 }
063 }
064 ATH_MSG_DEBUG( "Exec loop for track " << ntrk );
065 }
066 }
067 msg (MSG::INFO)<<"Created a total of " << tracks->size() << " aligntrks" << endreq;
068 }
069
070 bool InDetAlignNt::iTrkTrkSel(const Trk::Track* pnt) {
071
072 if (par_trksel) {
073 return (m_trkseltool->getStatus(*pnt)!=0);
074 } else {
075
076
077 bool momok;
078
079 const Trk::MeasuredPerigee *mesp = 0;
080 DataVector<const Trk::TrackStateOnSurface>::const_iterator tsos_it =
081 pnt->trackStateOnSurfaces()->begin();
082 DataVector<const Trk::TrackStateOnSurface>::const_iterator tsos_end =
083 pnt->trackStateOnSurfaces()->end();
084 for (; tsos_it != tsos_end; ++tsos_it) {
085 if ((**tsos_it).type(Trk::TrackStateOnSurface::Perigee))
086 mesp = dynamic_cast<const Trk::MeasuredPerigee *>((**tsos_it).trackParameters());
087 }
088
089 if (mesp!=0) {
090 momok=(fabs(mesp->parameters()[Trk::qOverP])/
091 sin(mesp->parameters()[Trk::theta]) < 1/par_ptmin);
092 } else {
093 momok=false;
094 }
095 return momok;
096 }
097 }
098
099 AlignTrk InDetAlignNt::TrkTrkProc(const DataVector<Trk::Track>* trks,
100 const Trk::Track* pnt,const int nt,bool& status) {
101 AlignTrk aligntrk;
102 SetCommon(aligntrk,nt);
103 ATH_MSG_DEBUG( "TrkTrkProc for track " << nt );
104 ATH_MSG_VERBOSE(*pnt );
105
106
107
108
109 float dx1=1.;
110 float dx2=1.;
111 const Trk::MeasuredPerigee *mesp = 0;
112 DataVector<const Trk::TrackStateOnSurface>::const_iterator tsos_it =
113 pnt->trackStateOnSurfaces()->begin();
114 DataVector<const Trk::TrackStateOnSurface>::const_iterator tsos_end =
115 pnt->trackStateOnSurfaces()->end();
116 for (; tsos_it != tsos_end; ++tsos_it) {
117 if ((**tsos_it).type(Trk::TrackStateOnSurface::Perigee))
118 mesp = dynamic_cast<const Trk::MeasuredPerigee *>((**tsos_it).trackParameters());
119 }
120
121 if (mesp!=NULL) {
122
123 double d0 = mesp->parameters()[Trk::d0];
124 double z0 = mesp->parameters()[Trk::z0];
125 double phi0 = mesp->parameters()[Trk::phi0];
126 double theta = mesp->parameters()[Trk::theta];
127 double qOverP = mesp->parameters()[Trk::qOverP];
128
129 ATH_MSG_DEBUG( nt << ". Track Parameters (d0, z0, phi0, theta, q/p)" );
130 ATH_MSG_DEBUG( " " << d0 << ", " << z0 << ", "
131 << phi0 << ", " << theta << ", " << qOverP );
132
133
134 dx1=1./(pow(cos(mesp->parameters()[Trk::theta]),2));
135 dx2=1./sin(mesp->parameters()[Trk::theta]);
136 aligntrk.set_trkpar(mesp->parameters()[Trk::d0],1);
137 aligntrk.set_trkpar(mesp->parameters()[Trk::z0],2);
138 aligntrk.set_trkpar(mesp->parameters()[Trk::phi],3);
139 aligntrk.set_trkpar(1/tan(mesp->parameters()[Trk::theta]),4);
140 aligntrk.set_trkpar(dx2*(mesp->parameters()[Trk::qOverP]),5);
141 const Trk::ErrorMatrix& ermatrix=mesp->localErrorMatrix();
142 aligntrk.set_trkcov(ermatrix.covValue(Trk::d0,Trk::d0),1);
143 aligntrk.set_trkcov(ermatrix.covValue(Trk::z0,Trk::d0),2);
144 aligntrk.set_trkcov(ermatrix.covValue(Trk::z0,Trk::z0),3);
145 aligntrk.set_trkcov(ermatrix.covValue(Trk::phi,Trk::d0),4);
146 aligntrk.set_trkcov(ermatrix.covValue(Trk::phi,Trk::z0),5);
147 aligntrk.set_trkcov(ermatrix.covValue(Trk::phi,Trk::phi),6);
148 aligntrk.set_trkcov(dx1*ermatrix.covValue(Trk::theta,Trk::d0),7);
149 aligntrk.set_trkcov(dx1*ermatrix.covValue(Trk::theta,Trk::z0),8);
150 aligntrk.set_trkcov(dx1*ermatrix.covValue(Trk::theta,Trk::phi),9);
151 aligntrk.set_trkcov(dx1*dx1*ermatrix.covValue(Trk::theta,Trk::theta),10);
152 aligntrk.set_trkcov(dx2*ermatrix.covValue(Trk::qOverP,Trk::d0),11);
153 aligntrk.set_trkcov(dx2*ermatrix.covValue(Trk::qOverP,Trk::z0),12);
154 aligntrk.set_trkcov(dx2*ermatrix.covValue(Trk::qOverP,Trk::phi),13);
155 aligntrk.set_trkcov(dx1*dx2*ermatrix.covValue(Trk::qOverP,Trk::theta),14);
156 aligntrk.set_trkcov(dx2*dx2*ermatrix.covValue(Trk::qOverP,Trk::qOverP),15);
157 } else {
158 msg (MSG::ERROR) << "Could not get Trk::MeasuredPerigee" << endreq;
159 for (int i=1;i<=5;++i) aligntrk.set_trkpar(0.,i);
160 for (int i=1;i<=15;++i) aligntrk.set_trkcov(0.,i);
161 }
162
163
164 if (par_truth) TrkSetTruth(trks,pnt,aligntrk);
165 bool outofslice=false;
166
167
168
169
170 int nscat=0;
171 int nhitscat[MAXHITS];
172 HepPoint3D posscat[MAXHITS];
173
174 aligntrk.set_trtcons(0., 0.);
175 const Trk::TrackParameters* tparpscat[MAXHITS];
176 const Trk::ScatteringAngles* saospscat[MAXHITS];
177
178 for (std::vector<const Trk::TrackStateOnSurface*>::const_iterator tsos=pnt->trackStateOnSurfaces()->begin();
179 tsos!=pnt->trackStateOnSurfaces()->end(); ++tsos) {
180
181 const Trk::TrackParameters* tparp=((*tsos)->trackParameters());
182
183 const Trk::MeasurementBase* mesb=(*tsos)->measurementOnTrack();
184 const Trk::PseudoMeasurementOnTrack* trtm=dynamic_cast<const Trk::PseudoMeasurementOnTrack*>(mesb);
185 if (trtm!=0) {
186
187 ATH_MSG_DEBUG( "The TRT pseudo-measurement has been found! " );
188 ATH_MSG_DEBUG( "Q/p_t(perigee), Q/p_t(TRT), sigma(1,1) = " <<
189 dx2*(mesp->parameters()[Trk::qOverP]) << ", " <<
190 dx2*(trtm->localParameters()[Trk::qOverP]) << ", " <<
191 dx2*getErrorQonP(trtm) );
192 aligntrk.set_trtcons(dx2*(trtm->localParameters()[Trk::qOverP]), dx2*getErrorQonP(trtm));
193 }
194
195 if ((*tsos)->type(Trk::TrackStateOnSurface::Scatterer)) {
196 const Trk::GlobalPosition& globp=tparp->position();
197 posscat[nscat]=globp;
198 tparpscat[nscat]=tparp;
199 const Trk::MaterialEffectsBase* meb = (*tsos)->materialEffectsOnTrack();
200 const Trk::MaterialEffectsOnTrack* meot = dynamic_cast<const Trk::MaterialEffectsOnTrack*>(meb);
201 if (meot!=NULL) {
202 saospscat[nscat]=(meot)->scatteringAngles();
203 ATH_MSG_VERBOSE( "Scattering angle (" << nscat << "): " << saospscat[nscat]);
204
205 }
206 nhitscat[nscat]=0;
207 ++nscat;
208 }
209
210 }
211
212 ATH_MSG_DEBUG( "Number of Scatterers found: " << nscat );
213
214
215 int nhits=0;
216 bool seenhit=false;
217 Trk::GlobalPosition lastpos,pos;
218 ATH_MSG_DEBUG( "Begin loop over " << pnt->trackStateOnSurfaces()->size() << " TSOS" );
219 for (std::vector<const Trk::TrackStateOnSurface*>::const_iterator tsos=
220 pnt->trackStateOnSurfaces()->begin();
221 tsos!=pnt->trackStateOnSurfaces()->end(); ++tsos) {
222
223 const Trk::MeasurementBase* mesb=(*tsos)->measurementOnTrack();
224 const Trk::RIO_OnTrack* rotp=dynamic_cast<const Trk::RIO_OnTrack*>(mesb);
225 const Trk::TrackParameters* tparp=((*tsos)->trackParameters());
226 const Trk::ScatteringAngles* saosp;
227 const Trk::MaterialEffectsBase* meb = (*tsos)->materialEffectsOnTrack();
228 const Trk::MaterialEffectsOnTrack* meot = dynamic_cast<const Trk::MaterialEffectsOnTrack*>(meb);
229 if (meot!=NULL) {
230 saosp=(meot)->scatteringAngles();
231 ATH_MSG_DEBUG( "Scattering angles del " << saosp->deltaTheta()
232 << "," << saosp->deltaPhi() << "," << saosp->sigmaDeltaTheta() << "," << saosp->sigmaDeltaPhi() );
233 }
234 else saosp = NULL;
235 ATH_MSG_VERBOSE( "Si " << nhits
236 << " pointers ROT/TP/SO " << rotp
237 << "," << tparp << "," << saosp);
238
239
240
241 if (tparp!=0) {
242 const Trk::GlobalPosition& globp=tparp->position();
243 ATH_MSG_DEBUG( "Track par position [" << globp.x() << "," <<
244 globp.y() << "," << globp.z() << "]" );
245 }
246
247 bool sharedscat=false;
248 if (rotp!=0 && (*tsos)->type(Trk::TrackStateOnSurface::Measurement)) {
249 const Trk::PrepRawData* rio=rotp->prepRawData();
250 ATH_MSG_VERBOSE( "RIO PRD pointer " << rio );
251
252 Identifier ident=rotp->identify();
253 Identifier mident;
254 bool sihit=false;
255 bool pix=false;
256 bool sctec=false;
257
258
259 if (m_pixid->is_pixel(ident) || (par_layout==11 &&
260 (!(m_sctid->is_sct(ident) || m_trtid->is_trt(ident))))) {
261 mident=m_pixid->wafer_id(ident);
262 sihit=true;
263 pix=true;
264 } else if (m_sctid->is_sct(ident)) {
265
266 mident=m_sctid->wafer_id(m_sctid->barrel_ec(ident),
267 m_sctid->layer_disk(ident),m_sctid->phi_module(ident),
268 m_sctid->eta_module(ident),m_sctid->side(ident));
269 sihit=true;
270 }
271 if (sihit) {
272 AlignSiHit newhit;
273 nhits++;
274
275 const InDet::SiClusterOnTrack* siclus=
276 dynamic_cast<const InDet::SiClusterOnTrack*>(rotp);
277 if (siclus!=0) {
278 pos=siclus->globalPosition();
279 } else {
280 msg (MSG::ERROR) << "Could not determine hit position" << endreq;
281 }
282 newhit.set_sinlocal(0.0);
283 AlignSiModuleList::ModuleIDMap::const_iterator
284 itmod=p_modlist->idmap.find(mident);
285 if (itmod != p_modlist->idmap.end()) {
286 newhit.set_index(itmod->second);
287 if (itmod->second!=0) {
288 if (!pix) newhit.set_sinlocal((m_detelep[(itmod->second)-1])->
289 sinStereoLocal(pos));
290 const AlignSiModule& simod=p_modlist->vec[(itmod->second)-1];
291 simod.inchits();
292
293
294 if (!pix) {
295 const InDetDD::SiDetectorElement*
296 detelm=m_detelep[(itmod->second)-1];
297 if (detelm!=0 && detelm->design().shape()==
298 InDetDD::Trapezoid) {
299
300 Hep3Vector globy=
301 ((m_global2local[itmod->second-1]).inverse()*HepYHat);
302 float ster=atan2(simod.pos(1),simod.pos(0))-
303 atan2(globy.y(),globy.x());
304 if (ster<-M_PI) ster+=2*M_PI;
305 if (ster>M_PI) ster-=2*M_PI;
306 newhit.set_stereo(ster);
307 sctec=true;
308 } else {
309
310 Hep3Vector globx=
311 ((m_global2local[itmod->second-1]).inverse()*HepXHat);
312 newhit.set_stereo(asin(globx.z()));
313 }
314 }
315 } else {
316 outofslice=true;
317 }
318 } else {
319 msg (MSG::ERROR) << "Module ident " << mident.getString() <<
320 "not found in moduleidmap" << endreq;
321 newhit.set_index(-1);
322 }
323
324
325
326
327
328
329 if (siclus!=0) {
330 pos=siclus->globalPosition();
331 newhit.set_r(pos.perp());
332 newhit.set_phi(pos.phi());
333 newhit.set_z(pos.z());
334 if (seenhit) {
335 newhit.set_trklen((pos-lastpos).mag());
336 newhit.set_status(2);
337 } else {
338 newhit.set_trklen(0);
339 newhit.set_status(1);
340 seenhit=true;
341 }
342 } else {
343 msg (MSG::ERROR) << "RIO cast to SiClusterOnTrack failed" << endreq;
344 }
345
346
347 const Trk::ErrorMatrix& locale=rotp->localErrorMatrix();
348 if (pix) {
349 newhit.set_erphi_resid(locale.error(Trk::locX));
350 newhit.set_ez_resid(locale.error(Trk::locY));
351 } else {
352
353
354 if (sctec) {
355 newhit.set_erphi_resid(locale.error(Trk::locX)*sqrt
356 (1.-locale.covValue(Trk::locX,Trk::locY)*
357 locale.covValue(Trk::locX,Trk::locY)/
358 (locale.covValue(Trk::locX)*locale.covValue(Trk::locY))));
359 } else {
360 newhit.set_erphi_resid(locale.error(Trk::locX));
361 }
362 newhit.set_ez_resid(0.);
363 }
364
365
366 const Trk::MeasuredTrackParameters* mtparp = dynamic_cast<const Trk::MeasuredTrackParameters*>(tparp);
367 if (mtparp) {
368 const Trk::ErrorMatrix& extrapolationError=mtparp->errorMatrixInMeasurementFrame();
369 newhit.set_trkExtrapLocalXError(extrapolationError.error(Trk::loc1));
370 if(pix)
371 newhit.set_trkExtrapLocalYError(extrapolationError.error(Trk::loc2));
372
373 }
374
375
376
377
378 bool tplook=false;
379 if (tparp==0 || saosp==0) {
380 float d2min=1.E10;
381 int jhit=-1;
382 for (int j=0;j<nscat;++j) {
383 float d2=pow(pos.x()-posscat[j].x(),2.)+
384 pow(pos.y()-posscat[j].y(),2.)+
385 pow(pos.z()-posscat[j].z(),2.);
386 if (d2<d2min) {
387 d2min=d2;
388 jhit=j;
389 }
390 }
391 if (jhit>=0) {
392
393 ATH_MSG_VERBOSE( "Associated scatterer " << jhit <<
394 " at distance " << sqrt(d2min) << " hit at r/z " << pos.perp() <<
395 "/" << pos.z() );
396
397 if (tparp==0) {
398 tparp=tparpscat[jhit];
399 tplook=true;
400 }
401 saosp=saospscat[jhit];
402
403 if (nhitscat[jhit]>0) sharedscat=true;
404 ++nhitscat[jhit];
405 }
406 }
407 if (tparp!=0) {
408
409 const HepVector& tparv=tparp->parameters();
410 if (!tplook) {
411
412
413 Trk::LocalParameters localp=rotp->localParameters();
414
415 if (pix) {
416 newhit.set_rphi_resid(localp[Trk::locX]-tparv[Trk::locX]);
417 newhit.set_z_resid(localp[Trk::locY]-tparv[Trk::locY]);
418
419
420 } else {
421
422
423
424 float trapcor=0.0;
425 if ( sctec ) {
426 const InDetDD::SiDetectorElement* detelm=
427 m_detelep[newhit.index()-1];
428 if (localp.contains(Trk::locY) && detelm!=0) {
429 trapcor = detelm->sinStereoLocal(Trk::LocalPosition(localp[Trk::locX],localp[Trk::locY]))
430 *(tparv[Trk::locY]-localp[Trk::locY]);
431 newhit.set_sinlocal(detelm->sinStereoLocal(Trk::LocalPosition(localp[Trk::locX],localp[Trk::locY])));
432 } else {
433 trapcor = newhit.sinlocal()*tparv[Trk::locY];
434 }
435 if (newhit.sinlocal()!=0 && !localp.contains(Trk::locY)) msg(MSG::WARNING) <<
436 "No hit locY information for an SCT EC module" << endreq;
437 if (detelm==0) msg (MSG::ERROR) <<
438 "Zero detelm pointer - no trapezoid correction" << endreq;
439 }
440
441
442 newhit.set_rphi_resid((localp[Trk::locX]-trapcor
443 -tparv[Trk::locX])*sqrt(1-newhit.sinlocal()*newhit.sinlocal()));
444
445
446
447
448 }
449 if (localp.contains(Trk::locY)) {
450 newhit.set_delzr(tparv[Trk::locY]-localp[Trk::locY]);
451 } else {
452 newhit.set_delzr(tparv[Trk::locY]);
453 }
454
455
456 } else {
457
458 FillResid(rotp,tparp,newhit);
459 }
460
461 if (par_numdiv<31)
462 FillDeriv(newhit,aligntrk.trkpar(1),aligntrk.trkpar(3),
463 aligntrk.trkpar(5),pos,1/tan(tparv[Trk::theta]),tparv[Trk::phi]);
464
465
466 if (par_numdiv>0) FillNumDeriv(aligntrk,tparp,newhit);
467 }
468
469 if (pix) {
470 newhit.set_chan1(m_pixid->phi_index(ident));
471 newhit.set_chan2(m_pixid->eta_index(ident));
472 } else {
473 newhit.set_z_resid(0.);
474 newhit.set_ez_resid(0.);
475 newhit.set_chan1(m_sctid->strip(ident));
476 newhit.set_chan2(0);
477 }
478
479
480
481
482
483 newhit.set_locdel(0.,0);
484 newhit.set_locdel(0.,1);
485 newhit.set_locedel(0.,0);
486 newhit.set_locedel(0.,1);
487 newhit.set_globdel(0.,0);
488 newhit.set_globdel(0.,1);
489 newhit.set_globedel(0.,0);
490 newhit.set_globedel(0.,1);
491
492 if (sharedscat) saosp=0;
493 if (tparp!=0) {
494 const HepVector& tparv=tparp->parameters();
495 FillAngles(newhit,tparv,saosp);
496
497 lastpos=pos;
498 }
499
500 aligntrk.add_hit(newhit);
501 }
502 }
503 }
504 aligntrk.set_nscat(nscat);
505 ATH_MSG_DEBUG( "Track has " << nhits << " Si hits and " << nscat <<
506 " scatterers" );
507
508
509 if (par_trt!=0) {
510 int ntrthits=0;
511
512
513 for (std::vector<const Trk::TrackStateOnSurface*>::const_iterator tsos=
514 pnt->trackStateOnSurfaces()->begin();
515 tsos!=pnt->trackStateOnSurfaces()->end();++tsos) {
516 const Trk::MeasurementBase* mesb=(*tsos)->measurementOnTrack();
517 const Trk::RIO_OnTrack* rotp=dynamic_cast<const Trk::RIO_OnTrack*>(mesb);
518 const Trk::TrackParameters* tparp=((*tsos)->trackParameters());
519 ATH_MSG_VERBOSE( "TRT " << ntrthits <<
520 " pointers ROT/TP " << rotp << "," << tparp );
521 if (rotp!=0 && (*tsos)->type(Trk::TrackStateOnSurface::Measurement)) {
522 const Trk::PrepRawData* rio=rotp->prepRawData();
523 Identifier ident=rotp->identify();
524
525 if (m_trtid->is_trt(ident)) {
526 AlignTRTHit newhit;
527 newhit.set_index(m_trtid->barrel_ec(ident),
528 m_trtid->layer_or_wheel(ident),m_trtid->phi_module(ident),
529 m_trtid->straw_layer(ident),m_trtid->straw(ident));
530
531 const InDet::TRT_DriftCircleOnTrack* trtcirc=
532 dynamic_cast<const InDet::TRT_DriftCircleOnTrack*>(rotp);
533 if (trtcirc!=0) {
534 Trk::GlobalPosition tpos=trtcirc->globalPosition();
535 newhit.set_r(tpos.perp());
536 newhit.set_phi(tpos.phi());
537 newhit.set_z(tpos.z());
538 newhit.set_highth(trtcirc->highLevel());
539
540
541 float xdrift=rotp->localParameters()[Trk::locX];
542 newhit.set_xdrift(xdrift);
543 bool isvalid=false;
544
545 const InDet::TRT_DriftCircle* rdc=dynamic_cast<const
546 InDet::TRT_DriftCircle*>(rio);
547 if (rdc) {
548 float tdrift=rdc->driftTime(isvalid);
549 if (!isvalid) tdrift=-1.;
550 newhit.set_tdrift(tdrift);
551 }
552
553 if (tparp!=0) {
554 const HepVector& tparv=tparp->parameters();
555 newhit.set_res(xdrift-tparv[0]);
556 } else {
557
558 float d2min=1.E10;
559 int jhit=-1;
560 for (int j=0;j<nscat;++j) {
561 float d2=pow(tpos.x()-posscat[j].x(),2.)+
562 pow(tpos.y()-posscat[j].y(),2.)+
563 pow(tpos.z()-posscat[j].z(),2.);
564 if (d2<d2min) {
565 d2min=d2;
566 jhit=j;
567 }
568 }
569 if (jhit>=0) {
570
571
572 }
573 }
574 ntrthits++;
575 aligntrk.add_trthit(newhit);
576 } else {
577 msg (MSG::ERROR) << "TRT drift RIO cast failed - no hit stored"
578 << endreq;
579 }
580 }
581 }
582 }
583 ATH_MSG_DEBUG( "Track has " << ntrthits << " TRT hits " );
584 }
585 status=true;
586 return aligntrk;
587 }
588
589 void InDetAlignNt::FillAngles(AlignSiHit& ahit,const HepVector& tparv,
590 const Trk::ScatteringAngles* saosp) {
591
592
593 HepRotation& global2local=m_global2local[ahit.index()-1];
594
595 float gtheta=tparv[Trk::theta];
596 float gphi=tparv[Trk::phi];
597
598 HepPoint3D local_dir=global2local*
599 HepPoint3D(sin(gtheta)*cos(gphi),sin(gtheta)*sin(gphi),cos(gtheta));
600
601 ahit.set_locang(local_dir.theta(),0);
602 ahit.set_locang(local_dir.phi(),1);
603 ahit.set_globang(1/tan(gtheta),0);
604 ahit.set_globang(gphi,1);
605
606
607 if (saosp!=NULL) {
608
609 ATH_MSG_DEBUG ( "-> gtheta " << gtheta << ", gphi " << gphi );
610 ATH_MSG_DEBUG ( "Filled Scattering angles del " << saosp->deltaTheta()
611 << "," << saosp->deltaPhi() << "," << saosp->sigmaDeltaTheta() << "," <<
612 saosp->sigmaDeltaPhi() );
613
614 float ngtheta=gtheta+saosp->deltaTheta();
615 float ngphi=gphi+saosp->deltaPhi();
616 float rgtheta=gtheta+saosp->sigmaDeltaTheta();
617 float rgphi=gphi+saosp->sigmaDeltaPhi();
618
619 HepPoint3D new_local_dir=global2local*
620 HepPoint3D(sin(ngtheta)*cos(ngphi),sin(ngtheta)*sin(ngphi),cos(ngtheta));
621 HepPoint3D rms_local_dir=global2local*
622 HepPoint3D(sin(rgtheta)*cos(rgphi),sin(rgtheta)*sin(rgphi),cos(rgtheta));
623
624 ahit.set_locdel(new_local_dir.theta()-local_dir.theta(),0);
625 ahit.set_locdel(new_local_dir.phi()-local_dir.phi(),1);
626 ahit.set_locedel(fabs(rms_local_dir.theta()-local_dir.theta()),0);
627 ahit.set_locedel(fabs(rms_local_dir.phi()-local_dir.phi()),1);
628
629 if (m_ipattrk) {
630 ahit.set_locedel(sqrt(ahit.locedel(0)),0);
631 ahit.set_locedel(sqrt(ahit.locedel(1)),1);
632 }
633 ahit.set_globdel(saosp->deltaTheta(),0);
634 ahit.set_globdel(saosp->deltaPhi(),1);
635 ahit.set_globedel(saosp->sigmaDeltaTheta(),0);
636 ahit.set_globedel(saosp->sigmaDeltaPhi(),1);
637 }
638 }
639
640 void InDetAlignNt::FillResid(const Trk::RIO_OnTrack* rotp,
641 const Trk::TrackParameters* tparp, AlignSiHit& newhit) {
642 Trk::LocalParameters localp=rotp->localParameters();
643
644
645
646
647
648
649
650 newhit.set_rphi_resid(-1.);
651 newhit.set_z_resid(-1.);
652
653 if (newhit.index()>0) {
654
655
656
657
658
659 const InDetDD::SiDetectorElement* detelm=m_detelep[newhit.index()-1];
660 const Trk::TrackParameters* paratsurf=
661 m_extraptool->extrapolate(*tparp,
662 detelm->surface(),Trk::anyDirection,false);
663 const HepVector& tpars=paratsurf->parameters();
664 float xresid;
665 if (detelm->isSCT() && detelm->design().shape()==InDetDD::Trapezoid) {
666
667 const HepPoint3D& rtrk=paratsurf->position();
668 float rcor=newhit.r()/sqrt(rtrk.x()*rtrk.x()+rtrk.y()*rtrk.y());
669 xresid=localp[Trk::locX]-rcor*tpars[Trk::locX];
670 xresid*=sqrt(1-newhit.sinlocal()*newhit.sinlocal());
671 } else {
672
673 xresid=localp[Trk::locX]-tpars[Trk::locX];
674 }
675
676
677
678
679 newhit.set_rphi_resid(xresid);
680 if (detelm->isPixel()) {
681
682 newhit.set_z_resid(localp[Trk::locY]-tpars[Trk::locY]);
683 }
684 newhit.set_delzr(tpars[Trk::locY]-localp[Trk::locY]);
685 } else {
686 msg (MSG::ERROR) << "Cannot extrapolate as hit index zero " << endreq;
687 }
688 }
689
690 void InDetAlignNt::TrkSetTruth(const DataVector<Trk::Track>* trks,
691 const Trk::Track* ptrk,AlignTrk& aligntrk) {
692
693
694
695
696 const TrackTruthCollection* trucp;
697 if (StatusCode::SUCCESS==evtStore()->retrieve(trucp,par_trktruthcol)) {
698 ElementLink<TrackCollection> tlink;
699 tlink.setElement(const_cast<Trk::Track*>(ptrk));
700 tlink.setStorableObject(*trks);
701 TrackTruthCollection::const_iterator itmp=
702 trucp->find(Trk::TrackTruthKey(tlink));
703 if (itmp!=trucp->end()) {
704 ATH_MSG_DEBUG( "Track has truth barcode " <<
705 itmp->second.barcode());
706 if ((itmp->second.barcode()!=0) && m_mccolln!=0)
707 FillTruth(aligntrk,(*m_mccolln->begin())
708 ->barcode_to_particle(itmp->second.barcode()));
709 } else {
710 ATH_MSG_DEBUG( "Cannot find truth in map!" );
711 }
712 } else {
713 msg (MSG::ERROR) << "Cannot find " << par_trktruthcol << endreq;
714 }
715 }
716
717 void InDetAlignNt::FillNumDeriv(const AlignTrk& altrk,
718 const Trk::TrackParameters* tparp,
719 AlignSiHit& newhit) {
720
721
722
723
724
725
726
727 double peritrk[5];
728 const Trk::PerigeeSurface perisurf;
729 const Trk::TrackParameters* pperitrk=
730 m_extraptool->extrapolate(*tparp,perisurf,
731 Trk::anyDirection,false);
732 if (pperitrk!=0) {
733 const HepVector& peripars=pperitrk->parameters();
734 peritrk[0]=peripars[Trk::d0];
735 peritrk[1]=peripars[Trk::z0];
736 peritrk[2]=peripars[Trk::phi0];
737 peritrk[3]=peripars[Trk::theta];
738 peritrk[4]=peripars[Trk::qOverP];
739 double sinth=sin(peritrk[3]);
740 double jac[5];
741 jac[0]=1.;
742 jac[1]=1.;
743 jac[2]=1.;
744 jac[3]=-sinth*sinth;
745 jac[4]=sinth;
746
747
748 const Trk::Surface& detsurf=(m_detelep[newhit.index()-1])->surface();
749 const Trk::TrackParameters* pbasetrk=
750 m_extraptool->extrapolate(*pperitrk,detsurf,
751 Trk::alongMomentum,false);
752 if (pbasetrk!=0) {
753 const HepVector& basepars=pbasetrk->parameters();
754
755
756
757 double varsig=0.2;
758 double varpar[5],deltrkp[5];
759 varpar[0]=varsig*sqrt(altrk.trkcov(1));
760 varpar[1]=varsig*sqrt(altrk.trkcov(3));
761 varpar[2]=varsig*sqrt(altrk.trkcov(6));
762 varpar[3]=varsig*sqrt(altrk.trkcov(10))*jac[3];
763 varpar[4]=varsig*sqrt(altrk.trkcov(15))*jac[4];
764
765 for (int ipar=0;ipar<5;++ipar) {
766
767 if (par_numdiv & (1 << ipar)) {
768
769 for (int j=0;j<5;++j) deltrkp[j]=peritrk[j];
770 deltrkp[ipar]=deltrkp[ipar]+varpar[ipar];
771
772 if (ipar==3) {
773 deltrkp[4]=deltrkp[4]-peritrk[4]*varpar[3]/tan(peritrk[3]);
774 }
775 const Trk::TrackParameters* pperivar=new Trk::Perigee(deltrkp[0],
776 deltrkp[1],
777 deltrkp[2],
778 deltrkp[3],
779 deltrkp[4],
780 perisurf);
781
782
783 const Trk::TrackParameters* pdeltrk= m_extraptool->extrapolate(*pperivar,detsurf,
784 Trk::alongMomentum,false);
785 if (pdeltrk!=0) {
786 const HepVector& delpars=pdeltrk->parameters();
787
788 if(std::isnan(delpars[Trk::locX])) msg (MSG::WARNING) << " delpars[Trk::locX] is nan" << endreq;
789 if(std::isnan(delpars[Trk::locY])) msg (MSG::WARNING) << " delpars[Trk::locY] is nan" << endreq;
790
791
792 if (std::isnan(delpars[Trk::locX]) || std::isnan(delpars[Trk::locY])) {
793 msg (MSG::WARNING )<< " -> skip extrapolation" << endreq;
794 continue;
795 }
796
797 if(std::isnan(delpars[Trk::locX])) msg (MSG::ERROR) << " delpars[Trk::locX] is nan" << endreq;
798 if(std::isnan(delpars[Trk::locY])) msg (MSG::ERROR) << " delpars[Trk::locY] is nan" << endreq;
799
800 ATH_MSG_VERBOSE( "Extrapolation " << ipar << " new local " <<
801 delpars[Trk::locX] << "," << delpars[Trk::locY] );
802
803 ATH_MSG_VERBOSE( "Shifts are " << ipar << " by " <<
804 varpar[ipar] << " : " << delpars[Trk::locX]-basepars[Trk::locX] <<
805 "," << delpars[Trk::locY]-basepars[Trk::locY]);
806
807
808 double sl(newhit.sinlocal());
809 double cl(sqrt(1.0-newhit.sinlocal()*newhit.sinlocal()));
810 double dlocx(cl*(basepars[Trk::locX]-delpars[Trk::locX])
811 + sl*(basepars[Trk::locY]-delpars[Trk::locY]));
812 double dlocy(cl*(basepars[Trk::locY]-delpars[Trk::locY])
813 - sl*(basepars[Trk::locX]-delpars[Trk::locX]));
814
815
816 newhit.set_dtpar1(jac[ipar]*dlocx/varpar[ipar],ipar);
817 newhit.set_dtpar2(jac[ipar]*dlocy/varpar[ipar],ipar);
818 delete pdeltrk;
819 }
820 else {
821 msg (MSG::ERROR) << "FillNumDeriv Detla extraplation " << ipar << "failed" << endreq;
822 }
823 delete pperivar;
824 }
825 }
826 delete pbasetrk;
827 } else {
828 msg (MSG::ERROR) << "FillNumDeriv base extrapolation failed" << endreq;
829 }
830 delete pperitrk;
831 } else {
832 msg (MSG::ERROR) << "FillNumDeriv backextrapolation failed from track: "
833 << *tparp << endreq;
834 }
835 }
836
837 void InDetAlignNt::FillNumDeriviPat(const AlignTrk& aligntrk,
838 AlignSiHit& newhit) {
839
840
841
842 const Trk::PerigeeSurface perisurf;
843 float theta=atan(1./aligntrk.trkpar(4));
844 if (theta<0) theta+=3.1415926;
845 const Trk::TrackParameters* tparp=new Trk::Perigee(aligntrk.trkpar(1),
846 aligntrk.trkpar(2),
847 aligntrk.trkpar(3),
848 theta,
849 sin(theta)*aligntrk.trkpar(5),
850 perisurf);
851 FillNumDeriv(aligntrk,tparp,newhit);
852 delete tparp;
853 }
854
855
856 double InDetAlignNt::getErrorQonP( const Trk::PseudoMeasurementOnTrack* pmot ) const{
857
858
859
860 int nmeas(0);
861 if (pmot->localParameters().contains(Trk::locRPhi)){
862 nmeas++;
863 }
864 if (pmot->localParameters().contains(Trk::locZ)){
865 nmeas++;
866 }
867 if (pmot->localParameters().contains(Trk::phi)){
868 nmeas++;
869 }
870 if (pmot->localParameters().contains(Trk::theta)){
871 nmeas++;
872 }
873 if (pmot->localParameters().contains(Trk::qOverP)){
874 nmeas++;
875 }
876 if (nmeas == 1){
877 return pmot->localErrorMatrix().error(Trk::locRPhi);
878 } else if (nmeas == 2){
879 return pmot->localErrorMatrix().error(Trk::locZ);
880 } else if (nmeas == 3){
881 return pmot->localErrorMatrix().error(Trk::phi);
882 } else if (nmeas == 4){
883 return pmot->localErrorMatrix().error(Trk::theta);
884 } else if (nmeas == 5){
885 return pmot->localErrorMatrix().error(Trk::qOverP);
886 }
887
888 return 0.;
889
890 }
891
Due to the LXR bug, the updates fail sometimes to remove references to deleted files. The Saturday's full rebuilds fix these problems |
This page was automatically generated by the
LXR engine.
|
|