Report problems to ATLAS LXR Team (with time and IP address indicated)

The LXR Cross Referencer
[ source navigation ]
[ diff markup ]
[ identifier search ]
[ general search ]
 
 
Architecture: [ linux ]
Version: [ head ] [ nightly ] [ GaudiDev ]
  Links to LXR source navigation pages for stable releases [ 12.*.* ]   [ 13.*.* ]   [ 14.*.* ]   [ 15.*.* ]  
001 #include "JetTagTools/SharedHitMapper.h"
002 #include "Particle/TrackParticle.h"
003 #include "Particle/TrackParticleContainer.h"
004 #include "StoreGate/StoreGateSvc.h"
005 #include "TrkTrack/Track.h"
006 #include <string>
007 #include "InDetIdentifier/PixelID.h"
008 #include "InDetIdentifier/SCT_ID.h"
009 #include "TrkRIO_OnTrack/RIO_OnTrack.h"
010 #include "TrkTrackSummary/TrackSummary.h"
011 #include <vector>
012 #include <utility>
013 #include <iostream>
014 
015 typedef std::pair<const Rec::TrackParticle*, const Trk::RIO_OnTrack*> pairTrkRio;
016 typedef std::vector<pairTrkRio> vecpairTrkRio;
017 
018 
019 SharedHitTrackAssoc::SharedHitTrackAssoc() {
020 }
021 SharedHitTrackAssoc::~SharedHitTrackAssoc() {
022 }
023 void SharedHitTrackAssoc::add(const Rec::TrackParticle* const trk, int shPattern) {
024   m_assocs.insert(std::make_pair(trk,shPattern));
025 }
026 void SharedHitTrackAssoc::add(const Rec::TrackParticle* const trk, int shB, int shP, int shS) {
027   int shPattern = 100000*shB + 1000*shS + 100*shP;
028   m_assocs.insert(std::make_pair(trk,shPattern));
029 }
030 int SharedHitTrackAssoc::numberSharedBLayer(const Rec::TrackParticle* const trk) const {
031   int tempSharedPattern(0);
032   AssocIter aE = m_assocs.end();
033   AssocIter aI = m_assocs.find(trk);
034   if(aI!=aE) tempSharedPattern = aI->second;
035   return tempSharedPattern/100000;
036 }
037 int SharedHitTrackAssoc::numberSharedPix(const Rec::TrackParticle* const trk) const {
038   int tempSharedPattern(0);
039   AssocIter aE = m_assocs.end();
040   AssocIter aI = m_assocs.find(trk);
041   if(aI!=aE) tempSharedPattern = aI->second;
042   return ((tempSharedPattern%100000)%1000)/100;
043 }
044 int SharedHitTrackAssoc::numberSharedSct(const Rec::TrackParticle* const trk) const {
045   int tempSharedPattern(0);
046   AssocIter aE = m_assocs.end();
047   AssocIter aI = m_assocs.find(trk);
048   if(aI!=aE) tempSharedPattern = aI->second;
049   return (tempSharedPattern%100000)/1000;
050 }
051 
052 void SharedHitTrackAssoc::dump() const {
053   std::cout<<"SharedHitTrackAssoc has "<<this->size()<<" elements:"<<std::endl;
054   AssocIter aI = m_assocs.begin();
055   AssocIter aE = m_assocs.end();
056   for(; aI!=aE; ++aI) {
057     std::cout << "--> track: "
058               << " Eta= " << aI->first->eta()
059               << " Phi= " << aI->first->phi()
060               << " pT= " << aI->first->pt()
061               << " Shared pattern= " << aI->second 
062       // what follows is not nice...
063               << " shB= " << this->numberSharedBLayer(aI->first)
064               << " shP= " << this->numberSharedPix(aI->first)
065               << " shS= " << this->numberSharedSct(aI->first)
066               << std::endl;
067   }
068 }
069 
070 
071 SharedHitMapper::SharedHitMapper(const std::string& name, ISvcLocator* pSvcLocator) 
072   : Algorithm(name, pSvcLocator) {
073   declareProperty("inputTrackCollection", m_inputTrackCollection = "TrackParticleCandidate");
074   declareProperty("sharedHitMapLocation", m_sharedHitMapLocation = "SharedHitMap");
075   declareProperty("deltaRCut", m_deltaRCut = 1.0);
076   declareProperty("QualityOrdering", m_qualOrder = false);
077   declareProperty("useTrackSummaryShared", m_useTrackSummaryShared = true);
078 }
079 
080 SharedHitMapper::~SharedHitMapper() {
081 }
082 
083 StatusCode SharedHitMapper::initialize() {
084   MsgStream mlog(msgSvc(), name());
085   StatusCode sc = service("StoreGateSvc",m_StoreGate);
086   if (sc.isFailure()) {
087     mlog << MSG::FATAL << "StoreGate service not found !" << endreq;
088     return StatusCode::FAILURE;
089   }
090   if(!m_useTrackSummaryShared) {
091     StoreGateSvc* detStore = 0;
092     sc = service( "DetectorStore", detStore );
093     if (sc.isFailure()) {
094       mlog << MSG::FATAL << "Could not get DetectorStore" << endreq;
095       return sc;
096     }
097     sc = detStore->retrieve(m_pixelId, "PixelID");
098     if (sc.isFailure()) {
099       mlog << MSG::ERROR << "Could not get PixelID helper !" << endreq;
100     }
101     sc = detStore->retrieve(m_sctId, "SCT_ID");
102     if (sc.isFailure()) {
103       mlog << MSG::ERROR << "Could not get SCT_ID helper !" << endreq;
104     }
105   }
106   std::cout<<name()<<" with input tracks "<<m_inputTrackCollection<<std::endl;
107   std::cout<<name()<<" initialized"<<std::endl;
108   return StatusCode::SUCCESS;
109 }
110 
111 StatusCode SharedHitMapper::execute() {
112 
113   MsgStream mlog(msgSvc(), name());
114 
115   /** retrieve input tracks: */
116   const Rec::TrackParticleContainer* inputTracks(0);
117   StatusCode sc = m_StoreGate->retrieve(inputTracks, m_inputTrackCollection);
118   if (sc.isFailure()) {
119     mlog << MSG::ERROR << "TrackParticles " << m_inputTrackCollection << " not found in StoreGate." << endreq;
120     return StatusCode::SUCCESS;
121   }
122   mlog << MSG::VERBOSE << "TrackParticleContainer " << m_inputTrackCollection 
123        << " found." << endreq; 
124 
125   SharedHitTrackAssoc* assoc = new SharedHitTrackAssoc();
126   mlog << MSG::INFO << "Number of initial TrackParticles: " << inputTracks->size() << endreq;
127 
128   if(m_useTrackSummaryShared) {
129     // --- Use shared hit information from TrackSummary:
130     // loop on first track:
131     Rec::TrackParticleContainer::const_iterator trk1I(inputTracks->begin());
132     Rec::TrackParticleContainer::const_iterator trkE(inputTracks->end());
133     for (; trk1I!=trkE; ++trk1I) {
134       const Trk::TrackSummary* trks  = (*trk1I)->trackSummary();
135       int nbs = trks->get(Trk::numberOfBLayerSharedHits); if(nbs < 0) nbs = 0;
136       int nps = trks->get(Trk::numberOfPixelSharedHits); if(nps < 0) nps = 0;
137       int nss = trks->get(Trk::numberOfSCTSharedHits); if(nss < 0) nss = 0;
138       if(nbs+nps+nss>0) assoc->add( (*trk1I), nbs, nps, nss);
139     }
140   } else {
141     // --- Recompute shared hit information from TrkTracks:
142 
143     vecpairTrkRio myvecTrkRio(0);
144 
145     // loop on first track:
146     Rec::TrackParticleContainer::const_iterator trk1I(inputTracks->begin());
147     Rec::TrackParticleContainer::const_iterator trkE(inputTracks->end());
148     for (; trk1I!=trkE; ++trk1I) {
149       const Trk::Track* trk1 = (*trk1I)->originalTrack();
150       if(trk1==0) {
151         mlog << MSG::ERROR << "Can't get the original track ! Please run on ESD" << endreq;
152         return StatusCode::SUCCESS;
153       }
154       //
155       // For the quality computation
156       int q1 = 0 ,q2 = 0;
157       int np1 = 0, np2 = 0, ns1 = 0, ns2 = 0;
158       int nhp1 = 0, nhp2 = 0, nhs1 = 0, nhs2 = 0, nh1 = 0, nh2 = 0;
159       int nprec1 = 0, nprec2 = 0;
160       double chi21 = 0., chi22 = 0.;
161       if (m_qualOrder) { 
162         chi21 = (*trk1I)->fitQuality()->chiSquared();
163         const Trk::TrackSummary* trks  = (*trk1I)->trackSummary();
164         np1    = trks->get(Trk::numberOfPixelHits); if(np1 < 0) np1 = 0;
165         ns1    = trks->get(Trk::numberOfSCTHits);   if(ns1 < 0) ns1 = 0;
166         nprec1 = np1 + ns1;
167         nhp1   = trks->get(Trk::numberOfPixelHoles); if(nhp1 < 0) nhp1 = 0;
168         nhs1   = trks->get(Trk::numberOfSCTHoles);   if(nhs1 < 0) nhs1 = 0;
169         nh1    = nhp1 + nhs1;
170         q1     = 10000*(2*nprec1-nh1) - int(20.*chi21/float(nprec1));
171       }
172       //
173       HepLorentzVector t1p4 = (*trk1I)->hlv();
174       const DataVector<const Trk::MeasurementBase>* mb = trk1->measurementsOnTrack();
175       DataVector<const Trk::MeasurementBase>::const_iterator nextMb;
176       DataVector<const Trk::MeasurementBase>::const_iterator lastMb(mb->end());
177       mlog << MSG::VERBOSE << "First track (pT="<<t1p4.perp()<<", eta="<<t1p4.pseudoRapidity ()<<") has " << mb->size() << " RIOs on track" << endreq;
178       // init arrays for pixels:
179       bool pl[m_nbpl]; // layers
180       for(int i=0;i<m_nbpl;i++) m_npl[i] = 0;
181       bool pd[m_nbpd]; // disks
182       for(int i=0;i<m_nbpd;i++) m_npd[i] = 0;
183       // init arrays for sct:
184       bool sl[m_nbsl];
185       for(int i=0;i<m_nbsl;i++) m_nsl[i] = 0;
186       bool sd[m_nbsd];
187       for(int i=0;i<m_nbsd;i++) m_nsd[i] = 0;
188       // loop on second track:
189       Rec::TrackParticleContainer::const_iterator trk2I(inputTracks->begin());
190       for (; trk2I!=trkE; ++trk2I) {
191         if(trk2I==trk1I) continue;
192         const Trk::Track* trk2 = (*trk2I)->originalTrack();
193         if(trk2==0) {
194           mlog << MSG::ERROR << "Can't get the original track ! Please run on ESD" << endreq;
195           return StatusCode::SUCCESS;
196         }
197         HepLorentzVector t2p4 = (*trk2I)->hlv();
198         double deltaR = t1p4.deltaR(t2p4);
199         if(deltaR<m_deltaRCut) {
200           // Assign the shared hit to the worst quality track 
201           // If a hit is shared by more than 2 tracks, this will not be perfect...
202           // From IG : Quality = 10000*(2*Nhit-Nhole) - int(20.*Xi2/float(Nhit))
203           if (m_qualOrder) { 
204             chi22 = (*trk2I)->fitQuality()->chiSquared();
205             const Trk::TrackSummary* trks2 = (*trk2I)->trackSummary();
206             np2    = trks2->get(Trk::numberOfPixelHits); if(np2 < 0) np2 = 0;
207             ns2    = trks2->get(Trk::numberOfSCTHits);   if(ns2 < 0) ns2 = 0;
208             nprec2 = np2 + ns2;
209             nhp2   = trks2->get(Trk::numberOfPixelHoles); if(nhp2 < 0) nhp2 = 0;
210             nhs2   = trks2->get(Trk::numberOfSCTHoles);   if(nhs2 < 0) nhs2 = 0;        
211             nh2    = nhp2 + nhs2;
212             q2     = 10000*(2*nprec2-nh2) - int(20.*chi22/float(nprec2));
213             if (q1 > q2) continue;
214           }
215           // loop on 1st track clusters:
216           nextMb = mb->begin(); 
217           for(; nextMb != lastMb; ++nextMb) {
218             const Trk::RIO_OnTrack* nextRio = dynamic_cast<const Trk::RIO_OnTrack*>(*nextMb);
219             bool pixl = m_pixelId->is_pixel(nextRio->identify());
220             bool ssct = m_sctId->is_sct(nextRio->identify());
221             if( (!pixl) && (!ssct) ) continue;
222             // pixels:
223             for(int i=0;i<m_nbpl;i++) pl[i] = false;
224             for(int i=0;i<m_nbpd;i++) pd[i] = false;
225             if(pixl) {
226               pl[0] = m_pixelId->is_blayer(nextRio->identify());
227               if(m_pixelId->is_barrel(nextRio->identify())) {
228                 int player = m_pixelId->layer_disk(nextRio->identify());
229                 for(int i=1;i<m_nbpl;i++) pl[i] = ( player == i );
230               } else {
231                 int pdisk = m_pixelId->layer_disk(nextRio->identify());
232                 for(int i=0;i<m_nbpd;i++) pd[i] = ( pdisk == i );
233               }
234             }
235             // sct:
236             for(int i=0;i<m_nbsl;i++) sl[i] = false;
237             for(int i=0;i<m_nbsd;i++) sd[i] = false;
238             if(ssct) {
239               if(m_sctId->is_barrel(nextRio->identify())) {
240                 int slayer = m_sctId->layer_disk(nextRio->identify());
241                 for(int i=0;i<m_nbsl;i++) sl[i] = ( slayer == i );
242               } else {
243                 int sdisk = m_sctId->layer_disk(nextRio->identify());
244                 for(int i=0;i<m_nbsd;i++) sd[i] = ( sdisk == i );
245               }
246             }
247             // loop on 2nd track clusters:
248             const DataVector<const Trk::MeasurementBase>* mb2 = trk2->measurementsOnTrack();
249             DataVector<const Trk::MeasurementBase>::const_iterator nextMb2(mb2->begin());
250             DataVector<const Trk::MeasurementBase>::const_iterator lastMb2(mb2->end());
251             if (nextMb == mb->begin()) mlog << MSG::VERBOSE << "Second track (pT="<<t2p4.perp()<<") has " << mb2->size() << " RIOs on track" << endreq;
252             int nriosi = 0;
253             for(; nextMb2 != lastMb2; ++nextMb2) {
254               const Trk::RIO_OnTrack* nextRio2 = dynamic_cast<const Trk::RIO_OnTrack*>(*nextMb2);
255               bool pixl2 = m_pixelId->is_pixel(nextRio2->identify());
256               bool ssct2 = m_sctId->is_sct(nextRio2->identify());
257               if( (!pixl2) && (!ssct2) ) continue;
258               nriosi++;
259               if( nextRio2->identify() == nextRio->identify() ) {
260                 mlog << MSG::DEBUG << "Shared hit ! Track1 Pt = " << t1p4.perp() <<" Track2 Pt = " << t2p4.perp() << endreq;
261                 if (m_qualOrder) {
262                   mlog << MSG::VERBOSE << "Track1 chi2 = " << chi21 
263                        << " nprec = " << nprec1 << " (np,ns) = (" << np1 << " , " << ns1 << ") "
264                        << " nhole = " << nh1    << " (nhp,nhs) = (" << nhp1 << " , " << nhs1 << ") "
265                        << " quality = " << q1 << endreq;
266                   mlog << MSG::VERBOSE << "Track2 chi2 = " << chi22 
267                        << " nprec = " << nprec2 << " (np,ns) = (" << np2 << " , " << ns2 << ") "
268                        << " nhole = " << nh2    << " (nhp,nhs) = (" << nhp2 << " , " << nhs2 << ") "
269                        << " quality = " << q2 << endreq;
270                 }
271                 mlog << MSG::DEBUG << pl[0] << " "<< pl[1] << " " << pl[2] << endreq;
272                 mlog << MSG::DEBUG << pd[0] << " "<< pd[1] << " " << pd[2] << endreq;
273                 mlog << MSG::DEBUG << sl[0] << " "<< sl[1] << " " << sl[2] << " " << sl[3] << endreq;
274                 mlog << MSG::DEBUG << sd[0] << " "<< sd[1] << " " << sd[2] << " " << sd[3] << " " << sd[4] << " "<< sd[5] << " " << sd[6] << " " << sd[7] << " " << sd[8] << endreq;
275                 pairTrkRio p = std::make_pair(*trk1I,nextRio);
276                 for(int i=0;i<m_nbpl;i++) {
277                   if (pl[i] && std::find(myvecTrkRio.begin(), myvecTrkRio.end(), p) == myvecTrkRio.end()) {m_npl[i]++;myvecTrkRio.push_back(p);}
278                 }
279                 for(int i=0;i<m_nbpd;i++) {
280                   if (pd[i] && std::find(myvecTrkRio.begin(), myvecTrkRio.end(), p) == myvecTrkRio.end()) {m_npd[i]++;myvecTrkRio.push_back(p);}
281                 }
282                 for(int i=0;i<m_nbsl;i++) {
283                   if (sl[i] && std::find(myvecTrkRio.begin(), myvecTrkRio.end(), p) == myvecTrkRio.end()) {m_nsl[i]++;myvecTrkRio.push_back(p);}
284                 }
285                 for(int i=0;i<m_nbsd;i++) {
286                   if (sd[i] && std::find(myvecTrkRio.begin(), myvecTrkRio.end(), p) == myvecTrkRio.end()) {m_nsd[i]++;myvecTrkRio.push_back(p);}
287                 }
288               }
289             } // end loop on rio2
290           } // end loop on rio1
291         }
292       } // end loop on trk2
293         // now store the result:
294       if(this->numberSharedBLayer()+this->numberSharedPix()+this->numberSharedSct()>0) {
295         assoc->add( (*trk1I), this->numberSharedBLayer(), 
296                     this->numberSharedPix(),
297                     this->numberSharedSct() );
298       }
299     } // end loop on trk1
300   }
301 
302   mlog << MSG::INFO << "Size of SharedHitMap: " << assoc->size() << endreq;
303 
304   if (mlog.level() <= MSG::DEBUG) assoc->dump();
305 
306   if(m_StoreGate->record(assoc,m_sharedHitMapLocation,false).isFailure()) {
307     mlog << MSG::ERROR << "recording of SharedHitAssoc " << m_sharedHitMapLocation 
308          << " failed." << endreq;
309   } else {
310     mlog << MSG::VERBOSE << "Recording of SharedHitAssoc " << m_sharedHitMapLocation 
311          << " succeeded, coll is now const." << endreq;
312   }
313 
314 
315   return StatusCode::SUCCESS;
316 }
317 
318 StatusCode SharedHitMapper::finalize() {
319   return StatusCode::SUCCESS;
320 }
321 
322 int SharedHitMapper::numberSharedBLayer() const {
323   return m_npl[0];
324 }
325 
326 int SharedHitMapper::numberSharedPix() const {
327   int nbs = 0;
328   for(int i=0;i<m_nbpl;i++) nbs += m_npl[i];
329   for(int i=0;i<m_nbpd;i++) nbs += m_npd[i];
330   return nbs;
331 }
332 
333 int SharedHitMapper::numberSharedSct() const {
334   int nbs = 0;
335   for(int i=0;i<m_nbsl;i++) nbs += m_nsl[i];
336   for(int i=0;i<m_nbsd;i++) nbs += m_nsd[i];
337   return nbs;
338 }
339 
340
[ source navigation ] [ diff markup ] [ identifier search ] [ general search ]
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. Valid HTML 4.01!