LCFIVertex  0.7.2
lciointerface.cpp
1 #include <inc/lciointerface.h>
2 #include <EVENT/ReconstructedParticle.h>
3 #include <EVENT/Track.h>
4 #include <EVENT/LCCollection.h>
5 #include <EVENT/LCFloatVec.h>
6 #include <IMPL/LCCollectionVec.h>
7 #include <EVENT/LCEvent.h>
8 #include <IMPL/VertexImpl.h>
9 #include <IMPL/ParticleIDImpl.h>
10 #include <IMPL/ReconstructedParticleImpl.h>
11 #include <inc/track.h>
12 #include <inc/trackstate.h>
13 #include <inc/event.h>
14 #include <inc/vertex.h>
15 #include <util/inc/memorymanager.h>
16 #include <inc/decaychain.h>
17 
18 #include <algorithm>
19 #include <vector>
20 
21 using namespace lcio;
22 using namespace vertex_lcfi;
23 using std::vector;
24 using std::map;
25 namespace vertex_lcfi{
26 
27 vertex_lcfi::Track* trackFromLCIORP(Event* MyEvent, lcio::ReconstructedParticle* RP)
28 {
29  //Get the track from the RP
30  lcio::Track* RPTrack = *(RP->getTracks().begin());
31 
32  HelixRep H;
33  H.d0()=RPTrack->getD0();
34  H.z0()=RPTrack->getZ0();
35  H.phi()=RPTrack->getPhi();
36  H.invR()=RPTrack->getOmega();
37  H.tanLambda()=RPTrack->getTanLambda();
38  //std::cout << "H= " << H << std::endl;
39  Vector3 Mom;
40  Mom.x() = RP->getMomentum()[0];
41  Mom.y() = RP->getMomentum()[1];
42  Mom.z() = RP->getMomentum()[2];
43 
44  SymMatrix5x5 Cov;
45 
46  // Classic order
47  Cov(0,0)=RPTrack->getCovMatrix()[0]; // d0d0
48  Cov(3,0)=RPTrack->getCovMatrix()[6]; // d0z0
49  Cov(3,3)=RPTrack->getCovMatrix()[9]; // z0z0
50 
51  Cov(1,0)=RPTrack->getCovMatrix()[1];
52  Cov(1,1)=RPTrack->getCovMatrix()[2];
53  Cov(0,2)=RPTrack->getCovMatrix()[3];
54  Cov(1,2)=RPTrack->getCovMatrix()[4];
55  Cov(2,2)=RPTrack->getCovMatrix()[5];
56  Cov(1,3)=RPTrack->getCovMatrix()[7];
57 
58  Cov(2,3)=RPTrack->getCovMatrix()[8];
59  Cov(0,4)=RPTrack->getCovMatrix()[10];
60  Cov(1,4)=RPTrack->getCovMatrix()[11];
61  Cov(2,4)=RPTrack->getCovMatrix()[12];
62  Cov(3,4)=RPTrack->getCovMatrix()[13];
63  Cov(4,4)=RPTrack->getCovMatrix()[14];
64 
65  /*
66  FIXED BY ALEXEI IN CODE OF 070416
67  Email from Alexei:
68  Track parameterization follows ILC convention
69  described in Thomas Kraemer's LC-note LC-DET-2006-004. Only the
70  lower triangle of the cov. matrix is stored and the order of the
71  elements is:
72 
73  (omega,omega),
74  (omega,tanl), (tanl,tanl),
75  (omega,phi0), (tanl,phi0), (phi0,phi0),
76  (omega,d0), (tanl,d0), (phi0,d0), (d0,d0)
77  (omega,z0), (tanl,z0), (phi0,z0), (d0,z0), (z0,z0)
78 
79  Unfortunately, I've just noticed that the sequence of cov.
80  matrix elements is different from that proposed by Thomas in his LC-Note.
81  //Alexei's Order
82  Cov(0,0)=RPTrack->getCovMatrix()[9];
83  Cov(3,0)=RPTrack->getCovMatrix()[13];
84  Cov(3,3)=RPTrack->getCovMatrix()[14];
85  */
86  /*
87  Unneeded Fix for small errors
88  double errorvar = 1.0/(Mom.mag() * pow(sin(H.theta()),(3.0/2.0)));
89  double minerror = 3.0 + 6.5 * errorvar;
90  minerror = pow((minerror/1000.0),2.0);
91  if (Cov(0,0) < minerror)
92  Cov(0,0) = minerror;
93  */
94 
95  //std::cout << "Track cov:" << sqrt(Cov(0,0))*1000.0 << " " << sqrt(Cov(3,0))*1000.0 << " " << sqrt(Cov(3,3))*1000.0 << std::endl;
96  Track* MyTrack = new vertex_lcfi::Track(MyEvent,
97  H,
98  Mom,
99  RP->getCharge(),
100  Cov,
101  RPTrack->getSubdetectorHitNumbers(),
102  (void *)RP);
104 
105  //Commented Out as unneeded.
106  /*//LCIO Tracks have non origin PCA, correct for this
107  Vector3 refPoint(RPTrack->getReferencePoint()[0],
108  RPTrack->getReferencePoint()[1],
109  RPTrack->getReferencePoint()[2]);
110  if(RPTrack->isReferencePointPCA() && refPoint.distanceTo(Vector3(0,0,0)) > 0.0)
111  {
112  //Shift frame of reference so that the refpoint is the origin
113  TrackState* TS = MyTrack->makeState();
114  TS->swimToStateNearestXY(-refPoint);
115 
116  //Calculate the new track parameters at this point
117  double NewZ0 = TS->position().z();
118  double NewPhi = TS->phi();
119  double NewD0 = TS->position().mag(RPhi);
120  //We need to sign d0 appropriatly - same as invR if origin in circle, opposite if outside
121  //Find the center of the circle
122  double invR = MyTrack->helixRep().invR();
123  Vector3 center = Vector3(TS->position().x() + (sin(NewPhi)/invR),
124  TS->position().y() - (cos(NewPhi)/invR),
125  0.0);
126  //Test if the orgin is inside the circle
127  if(pow(Vector3(0.0,0.0,0.0).distanceTo(center),2.0) - pow(1.0/invR,2.0) > 0.0)
128  {
129  //Outside the circle - opposite sign to invR
130  NewD0 = (invR > 0.0 ? -NewD0 : NewD0);
131  }
132  else
133  {
134  //Inside the circle - same sign as invR
135  NewD0 = (invR > 0.0 ? NewD0 : -NewD0);
136  }
137 
138  //Set new parameters
139  MyTrack->helixRep().d0() = NewD0;
140  MyTrack->helixRep().z0() = NewZ0;
141  MyTrack->helixRep().phi() = NewPhi;
142  TS->swimToStateNearestXY(refPoint);
143  }
144  */
145  return MyTrack;
146 }
147 
148 vertex_lcfi::Jet* jetFromLCIORP(Event* MyEvent,lcio::ReconstructedParticle* RP)
149 {
150  //Make Jet
151  Jet* MyJet = new Jet(MyEvent, vector<vertex_lcfi::Track*>(),RP->getEnergy(),Vector3(RP->getMomentum()[0],RP->getMomentum()[1],RP->getMomentum()[2]),(void*)RP); //Empty Jet
152  MemoryManager<Jet>::Event()->registerObject(MyJet);
153  MyEvent->addJet(MyJet);
154 
155  //RPs that represent tracks Loop
156  vector<ReconstructedParticle*> LCIOJetRPs = RP->getParticles();
157  for (vector<ReconstructedParticle*>::const_iterator iRP = LCIOJetRPs.begin();iRP!=LCIOJetRPs.end();++iRP)
158  {
159  vector<lcio::Track*> RPsTracks = (*iRP)->getTracks();
160  if( ! RPsTracks.empty() ) {
161  vertex_lcfi::Track* MyTrack = trackFromLCIORP(MyEvent,*iRP);
162  MyEvent->addTrack(MyTrack);
163  MyJet->addTrack(MyTrack);
164  } else {
165  std::cout << "Warning lciointerface.cpp:31 RP with no Track - excluding" << std::endl;
166  }
167  }
168  return MyJet;
169 }
170 
171 lcio::Vertex* vertexFromLCFIVertex(vertex_lcfi::Vertex* MyLCFIVertex)
172 {
173  //Make the vertex
174  lcio::VertexImpl* MyLCIOVertex = new VertexImpl(); //Empty Vertex
175  //MemoryManager<lcio::VertexImpl>::Event()->registerObject(MyLCIOVertex);
176  MyLCIOVertex->setChi2(MyLCFIVertex->chi2()) ;
177  MyLCIOVertex->setProbability(MyLCFIVertex->probability()) ;
178  MyLCIOVertex->setPosition(MyLCFIVertex->position().x(),MyLCFIVertex->position().y(),MyLCFIVertex->position().z()) ;
179  EVENT::FloatVec Error;
180  Error.push_back(MyLCFIVertex->positionError()(0,0));
181  Error.push_back(MyLCFIVertex->positionError()(1,0));
182  Error.push_back(MyLCFIVertex->positionError()(1,1));
183  Error.push_back(MyLCFIVertex->positionError()(2,0));
184  Error.push_back(MyLCFIVertex->positionError()(2,1));
185  Error.push_back(MyLCFIVertex->positionError()(2,2));
186  MyLCIOVertex->setCovMatrix(Error);
187  MyLCIOVertex->setPrimary(MyLCFIVertex->isPrimary());
188  return MyLCIOVertex;
189 }
190 
191 vertex_lcfi::Vertex* vertexFromLCIOVertex(lcio::Vertex* LCIOVertex, Event* MyEvent)
192 {
193  Vector3 Pos(LCIOVertex->getPosition()[0],LCIOVertex->getPosition()[1],LCIOVertex->getPosition()[2]);
194  SymMatrix3x3 PosErr;
195  PosErr(0,0)=LCIOVertex->getCovMatrix()[0];
196  PosErr(1,0)=LCIOVertex->getCovMatrix()[1];
197  PosErr(1,1)=LCIOVertex->getCovMatrix()[2];
198  PosErr(2,0)=LCIOVertex->getCovMatrix()[3];
199  PosErr(2,1)=LCIOVertex->getCovMatrix()[4];
200  PosErr(2,2)=LCIOVertex->getCovMatrix()[5];
201 
202  vertex_lcfi::Vertex* LCFIVertex = new vertex_lcfi::Vertex(MyEvent, vector<Track*>(), Pos, PosErr, LCIOVertex->isPrimary(), LCIOVertex->getChi2(), LCIOVertex->getProbability());
204 
205  return LCFIVertex;
206 }
207 
208 //TODO NB - Does nto support a deacy chain that has tracks not in the jet
209 DecayChain* decayChainFromLCIORP(Jet* LCFIJet, ReconstructedParticle* DecayChainRP)
210 {
211  //First make a map associating LCFI tracks to LCIO Reconstructed particles
212  map<ReconstructedParticle*,Track*> LCFITrack;
213  vector<Track*> LCFITracks = LCFIJet->tracks();
214  for (vector<Track*>::const_iterator iTrack = LCFITracks.begin();iTrack < LCFITracks.end();++iTrack)
215  {
216  //std::cout << "Adding " << (ReconstructedParticle*)(*iTrack)->trackingNum() << " =L " << *iTrack << std::endl;
217  LCFITrack[(ReconstructedParticle*)(*iTrack)->trackingNum()] = *iTrack;
218  }
219 
220  vector<lcio::Vertex*> LCIOVertices;
221  //Add the primary first
222  LCIOVertices.push_back(DecayChainRP->getStartVertex());
223 
224  //Loop over RPs to find all the vertices
225  vector<ReconstructedParticle*> RPs = DecayChainRP->getParticles();
226 
227  for (vector<ReconstructedParticle*>::const_iterator iRP = RPs.begin();iRP < RPs.end();++iRP)
228  {
229  lcio::Vertex* MyVertex = (*iRP)->getStartVertex();
230  if(MyVertex)
231  {
232  vector<lcio::Vertex*>::const_iterator it = find(LCIOVertices.begin(),LCIOVertices.end(),MyVertex);
233  if (it == LCIOVertices.end())
234  {
235  LCIOVertices.push_back(MyVertex);
236  }
237  }
238  }
239 
240  DecayChain* NewDecayChain = new DecayChain(LCFIJet,vector<Track*>(),vector<vertex_lcfi::Vertex*>());
242 
243  vector<vertex_lcfi::Vertex*> LCFIVertices;
244  map<lcio::Vertex*,vertex_lcfi::Vertex*> LCFIVertex;
245  for (vector<lcio::Vertex*>::const_iterator iVertex = LCIOVertices.begin();iVertex < LCIOVertices.end();++iVertex)
246  {
247  vertex_lcfi::Vertex* NewVertex = vertexFromLCIOVertex((*iVertex), LCFIJet->event());
248  LCFIVertices.push_back(NewVertex);
249  LCFIVertex[(*iVertex)] = NewVertex;
250  NewDecayChain->addVertex(NewVertex);
251  }
252 
253  for (vector<ReconstructedParticle*>::const_iterator iRP = RPs.begin();iRP < RPs.end();++iRP)
254  {
255  if (!(*iRP)->getParticles().empty())
256  {
257  //Only add if there is a corresponding LCFI track
258  map<ReconstructedParticle*,Track*>::const_iterator iLCFITrack= LCFITrack.find((*iRP)->getParticles()[0]);
259  if (iLCFITrack != LCFITrack.end())
260  {
261  if ((*iRP)->getStartVertex())
262  {
263  LCFIVertex[(*iRP)->getStartVertex()]->addTrack(LCFITrack[(*iRP)->getParticles()[0]]);
264  }
265  else
266  {
267  NewDecayChain->addTrack(LCFITrack[(*iRP)->getParticles()[0]]);
268  }
269  }
270  }
271  }
272  return NewDecayChain;
273 }
274 
275 
276 ReconstructedParticle* addDecayChainToLCIOEvent(LCEvent* MyLCIOEvent, DecayChain* MyDecayChain, std::string VertexCollectionName, std::string TrackRPCollectionName, bool StoreTrackChiSquareds)
277 {
278  //Check for cols and add if needed
279  vector<std::string>::const_iterator it = find(MyLCIOEvent->getCollectionNames()->begin(),MyLCIOEvent->getCollectionNames()->end(),VertexCollectionName);
280  if (it == MyLCIOEvent->getCollectionNames()->end())
281  {
282  //Not found do add - TODO do these need memory managment?
283  LCCollection* MyCollection = new LCCollectionVec("Vertex");
284  MyLCIOEvent->addCollection(MyCollection,VertexCollectionName);
285  }
286  it = find(MyLCIOEvent->getCollectionNames()->begin(),MyLCIOEvent->getCollectionNames()->end(),TrackRPCollectionName);
287  if (it == MyLCIOEvent->getCollectionNames()->end())
288  {
289  //Not found do add
290  LCCollection* MyCollection = new LCCollectionVec("ReconstructedParticle");
291  MyLCIOEvent->addCollection(MyCollection,TrackRPCollectionName);
292  }
293  //Get the collections
294  LCCollection* VertexCollection = MyLCIOEvent->getCollection(VertexCollectionName);
295  LCCollection* TrackRPCollection = MyLCIOEvent->getCollection(TrackRPCollectionName);
296 
297  //Check that the decay chain has vertices
298  if (MyDecayChain->vertices().empty())
299  {
300  //We have nothing to do so return a null RP
301  return 0;
302  }
303 
304  //Not needed as the first one is used even if not primary
306  //if (!(*MyDecayChain->vertices().begin())->isPrimary())
307  //{
308  // //Throw something
309  //}
310 
311  //The first vertex is always assumed to be the IP
312  //We need to check to see if the primary vertex is already in the LCIO record
313  lcio::Vertex* PrimaryVertex = 0;
314  int nVerts = VertexCollection->getNumberOfElements() ;
315  for(int i=0; i< nVerts && !PrimaryVertex ; i++)
316  {
317  if(dynamic_cast<lcio::Vertex*>(VertexCollection->getElementAt(i))->isPrimary())
318  PrimaryVertex = dynamic_cast<lcio::Vertex*>(VertexCollection->getElementAt(i));
319  }
320 
321  //If no primary found make one and add it to the vertex collection
322  if (!PrimaryVertex)
323  {
324  //Take the first vertex and convert it
325  PrimaryVertex = vertexFromLCFIVertex(*(MyDecayChain->vertices().begin()));
326  //Add to the vertex collection
327  VertexCollection->addElement(PrimaryVertex);
328  }
329 
330  //If we are storing chi squareds add a collection for them and make a map of tracks and their chi squareds
331  //Note we only currently cope with a track being in one vertex (ie external legs)
332  map<Track*,double> ChiSquaredOf;
333  if (StoreTrackChiSquareds)
334  {
335  LCCollection* MyCollection = new LCCollectionVec("LCFloatVec");
336  MyLCIOEvent->addCollection(MyCollection,TrackRPCollectionName+"TrackChiSquareds");
337 
338  //Loop over vertices adding chi squareds of each track to the map
339  for (vector<vertex_lcfi::Vertex*>::const_iterator iVertex = MyDecayChain->vertices().begin();iVertex < MyDecayChain->vertices().end();++iVertex)
340  {
341  ChiSquaredOf.insert((*iVertex)->chi2OfTracks().begin(),(*iVertex)->chi2OfTracks().end());
342  }
343  }
344 
345  //We now take a copy of all tracks in all vertices, linking the new tracks to the original tracks by addParticle
346  //We keep track of correspondance between old and new with this map
347  map<ReconstructedParticle*,ReconstructedParticleImpl*> NewTrackRP;
348 
349  vector<Track*> AllTracks = MyDecayChain->allTracks();
350  for (vector<Track*>::const_iterator iTrack = AllTracks.begin();iTrack < AllTracks.end();++iTrack)
351  {
352  //Make an RP
353  ReconstructedParticleImpl* NewRP = new ReconstructedParticleImpl();
354  //MemoryManager<ReconstructedParticleImpl>::Event()->registerObject(NewRP);
355  ReconstructedParticle* OriginalRP = (ReconstructedParticle*)(*iTrack)->trackingNum();
356  //Link the RP to the orginal RP
357  NewRP->addParticle(OriginalRP);
358  NewTrackRP[OriginalRP] = NewRP;
359  TrackRPCollection->addElement(NewRP);
360 
361  //We now optionally store the chi squareds of this track in each vertex it is found in
362  //These are stored as a colletion of LCFloatVecs in the same order as TrackRPCollection
363  //The first element is the chi squared in the start vertex
364  //Eventually the second will be the chi squared in the second vertex (if any) but
365  //as this code does not yet cope with seen decaying particles that will be a future
366  //upgrade
367  if (StoreTrackChiSquareds)
368  {
369  LCFloatVec* ChiSquared = new LCFloatVec();
370  if(ChiSquaredOf.find(*iTrack) == ChiSquaredOf.end())
371  {
372  //Not found in any vertex - shouldn't happen store -1
373  ChiSquared->push_back(-1);
374  }
375  else
376  {
377  ChiSquared->push_back(ChiSquaredOf[*iTrack]);
378  }
379 
380  MyLCIOEvent->getCollection(TrackRPCollectionName+"TrackChiSquareds")->addElement(ChiSquared);
381  }
382  }
383  //We now associate tracks to the primary vertex
384  //We assume that the tracking pointer points to the ReconstructedParticle that the track was
385  //formed of, in future we may want to make tracks that don't have RP's already
386  vector<Track*> PrimaryTracks = (*MyDecayChain->vertices().begin())->tracks();
387  for (vector<Track*>::const_iterator iTrack = PrimaryTracks.begin();iTrack < PrimaryTracks.end();++iTrack)
388  {
389  NewTrackRP[(ReconstructedParticle*)(*iTrack)->trackingNum()]->setStartVertex(PrimaryVertex);
390  }
391 
392  //Remember the RPs we make
393  vector<ReconstructedParticle*> CreatedRPs;
394 
395  //Keep track of the last vertex so we can link them
396  lcio::Vertex* PreviousVertex = PrimaryVertex;
397 
398  //Now iterate over the other verts!
399  for (vector<vertex_lcfi::Vertex*>::const_iterator iVertex = ++(MyDecayChain->vertices().begin());iVertex < MyDecayChain->vertices().end();++iVertex)
400  {
401  //Make the LCIO vertex
402  VertexImpl* NewVertex = dynamic_cast<VertexImpl*>(vertexFromLCFIVertex(*iVertex));
403  //Make the RP that represents the decayed particle
404  ReconstructedParticleImpl* NewRP = new ReconstructedParticleImpl();
405  //MemoryManager<ReconstructedParticleImpl>::Event()->registerObject(NewRP);
406  CreatedRPs.push_back(NewRP);
407  //Set the parameters of this particle
408  double mom[3];
409  mom[0] = (*iVertex)->momentum().x();
410  mom[1] = (*iVertex)->momentum().y();
411  mom[2] = (*iVertex)->momentum().z();
412  NewRP->setMomentum(mom);
413  NewRP->setCharge((*iVertex)->charge());
414  //Link this IP to the previous vertex - currently only sorted in distance from IP by ZVTOP
415  NewRP->setStartVertex(PreviousVertex);
416  //Now link the particle to the vertex
417  NewVertex->setAssociatedParticle(NewRP);
418 
419  //And link in the tracks contained in this vertex
420  for (vector<Track*>::const_iterator iTrack = (*iVertex)->tracks().begin();iTrack < (*iVertex)->tracks().end();++iTrack)
421  {
422  ReconstructedParticleImpl* iRP = NewTrackRP[(ReconstructedParticle*)(*iTrack)->trackingNum()];
423  NewRP->addParticle(iRP);
424  iRP->setStartVertex(NewVertex);
425  }
426 
427  PreviousVertex = NewVertex;
428  //Store them in the LCIO collections
429  VertexCollection->addElement(NewVertex);
430  TrackRPCollection->addElement(NewRP);
431  }
432 
433  //Then make a a big RP of all tracks in the decay chain and set the start vertex
434  ReconstructedParticleImpl* DecayChainRP = new ReconstructedParticleImpl();
435  //MemoryManager<ReconstructedParticleImpl>::Event()->registerObject(DecayChainRP);
436  for (vector<Track*>::const_iterator iTrack = MyDecayChain->allTracks().begin();iTrack < MyDecayChain->allTracks().end();++iTrack)
437  {
438  DecayChainRP->addParticle(NewTrackRP[(ReconstructedParticle*)(*iTrack)->trackingNum()]);
439  }
440 
441  //Add the RPs we created as decayed Particles
442  for (vector<ReconstructedParticle*>::const_iterator iRP = CreatedRPs.begin();iRP < CreatedRPs.end();++iRP)
443  {
444  DecayChainRP->addParticle(*iRP);
445  }
446  DecayChainRP->setStartVertex(PrimaryVertex);
447 
448  return DecayChainRP;
449 }
450 
451  void ReconstructedParticleLCFI::removeParticle(EVENT::ReconstructedParticle* particle)
452  {
453  //checkAccess("ParticleIDLCFI::removeParticle") ;
454  ReconstructedParticleVec::iterator it = find(_particles.begin(),_particles.end(),particle);
455  if (it != _particles.end())
456  {
457  _particles.erase(it);
458  }
459  }
460  void ReconstructedParticleLCFI::wipePIDs()
461  {
462  _pid = EVENT::ParticleIDVec();
463  }
464 
465  void ReconstructedParticleLCFI::copyPIDsFrom(ReconstructedParticle* InputRP)
466  {
467  for(EVENT::ParticleIDVec::const_iterator iPID = InputRP->getParticleIDs().begin();iPID != InputRP->getParticleIDs().end();++iPID)
468  {
469  this->addParticleID(new ParticleIDImpl(*dynamic_cast<ParticleIDImpl*>(*iPID)));
470  }
471 
472  }
473 
474  void ReconstructedParticleLCFI::makeWritable()
475  {
476  this->setReadOnly(false);
477  }
478 }
void registerObject(T *pointer)
Register an object for memory management.
void addTrack(Track *Track)
Add Track.
Definition: event.cpp:45
const std::vector< Track * > & allTracks() const
All tracks contained in DecayChain.
Definition: decaychain.cpp:52
const std::vector< Track * > & tracks() const
Tracks.
Definition: jet.cpp:23
bool isPrimary() const
Is this vertex primary.
void addVertex(Vertex *Vertex)
Add Vertex.
Definition: decaychain.cpp:180
Event * event() const
Event.
Definition: jet.cpp:18
void addTrack(Track *Track)
Add Track.
Definition: decaychain.cpp:133
const std::vector< Vertex * > & vertices() const
Vertices contained in DecayChain.
Definition: decaychain.cpp:80
const Vector3 & position() const
Position.
double chi2() const
Chi Squared.
const SymMatrix3x3 & positionError() const
Position.
Unique Track representation.
double probability() const
Probability.