All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Groups
KinkFinder.cc
Go to the documentation of this file.
1 #include "KinkFinder.h"
2 #include <EVENT/LCRelation.h>
3 #include <EVENT/MCParticle.h>
4 #include <UTIL/LCRelationNavigator.h>
5 #include "marlin/Global.h"
6 #include "EVENT/LCCollection.h"
7 #include <ClusterShapes.h>
8 #include "EVENT/Track.h"
9 #include "IMPL/LCCollectionVec.h"
10 #include "IMPL/ReconstructedParticleImpl.h"
11 #include "IMPL/VertexImpl.h"
12 #include <math.h>
13 
14 #include <DD4hep/Detector.h>
15 #include <DD4hep/DD4hepUnits.h>
16 #include <DDRec/DetectorData.h>
17 
18 #include "HelixClass.h"
19 
20 using namespace lcio ;
21 using namespace marlin ;
22 
23 const float mMuon = 0.106; // GeV
24 //const float cTauMuon = 657000; // mm
25 const float mPion = 0.140; // GeV
26 const float cTauPion = 7800; // mm
27 const float mKaon = 0.494; // GeV
28 const float cTauKaon = 3714; // mm
29 const float mSigma = 1.189; // GeV
30 const float cTauSigma = 24.0; // mm
31 const float mLambda = 1.115; // GeV
32 //const float cTauLambda = 78.0; // mm
33 const float mHyperon = 1.314; // GeV
34 const float cTauHyperon= 87.0; // mm
35 const float mNeutron = 0.940; // GeV
36 const float mProton = 0.940; // GeV
37 
38 KinkFinder aKinkFinder ;
39 
40 
41 KinkFinder::KinkFinder() : Processor("KinkFinder") {
42 
43  _description = "Kink Finder Processor : finds kinks, prongs and splits" ;
44 
45  registerInputCollection(LCIO::TRACK,
46  "TrackCollection",
47  "Name of input collection of reconstructed particles",
48  _trackColName,
49  std::string("LDCTracks"));
50 
51  registerOutputCollection(LCIO::RECONSTRUCTEDPARTICLE,
52  "KinkRecoParticleCollection",
53  "Name of output collection of kink reconstructed particles",
54  _kinkRecoPartColName,
55  std::string("KinkRecoParticles"));
56 
57  registerOutputCollection(LCIO::VERTEX,
58  "KinkVertexCollection",
59  "Name of output collection of kink vertices",
60  _kinkVertexColName,
61  std::string("KinkVertices"));
62 
63  registerOutputCollection(LCIO::RECONSTRUCTEDPARTICLE,
64  "ProngRecoParticleCollection",
65  "Name of output collection of prong reconstructed particles",
66  _prongRecoPartColName,
67  std::string("ProngRecoParticles"));
68 
69  registerOutputCollection(LCIO::VERTEX,
70  "ProngVertexCollection",
71  "Name of output collection of prong vertices",
72  _prongVertexColName,
73  std::string("ProngVertices"));
74 
75 
76  registerOutputCollection(LCIO::RECONSTRUCTEDPARTICLE,
77  "SplitRecoParticleCollection",
78  "Name of output collection of split reconstructed particles",
79  _splitRecoPartColName,
80  std::string("SplitRecoParticles"));
81 
82  registerOutputCollection(LCIO::VERTEX,
83  "SplitVertexCollection",
84  "Name of output collection of split vertices",
85  _splitVertexColName,
86  std::string("SplitVertices"));
87 
88 
89  registerProcessorParameter("SplitTrackMaxFracDeltaP",
90  "Maximum fractional difference in deltaP for split tracks",
91  _maxSplitTrackFracDeltaP,
92  float(0.02));
93 
94  registerProcessorParameter("SplitTrackMaxDeltaP",
95  "Maximum difference in deltaP for split tracks",
96  _maxSplitTrackDeltaP,
97  float(0.2));
98 
99  registerProcessorParameter("CutOnRadius",
100  "Cuts on Kink radius",
101  _rKinkCut,
102  float(100.0));
103 
104  registerProcessorParameter("MinimumTrackHits",
105  "cuts on number of track hits",
106  _minTrackHits,
107  int(5));
108 
109  registerProcessorParameter("MaxDeltaTpcLayers",
110  "Cuts on Kink DeltaRxy for tracks in TPC",
111  _maxDeltaTpcLayers,
112  int(10));
113 
114 
115  registerProcessorParameter("PionDecayMassCut",
116  "Cuts on kink mass for pion decay",
117  _pionDecayMassCut,
118  float(0.03));
119 
120  registerProcessorParameter("KaonDecayMassCut",
121  "Cuts on kink mass for kaon decay",
122  _kaonDecayMassCut,
123  float(0.075));
124 
125  registerProcessorParameter("SigmaDecayMassCut",
126  "Cuts on kink mass for sigma decay",
127  _sigmaDecayMassCut,
128  float(0.15));
129 
130  registerProcessorParameter("SigmaTimeCut",
131  "Cuts on lifetime for sigma decay",
132  _sigmaTimeCut,
133  float(6.));
134 
135  registerProcessorParameter("HyperonDecayMassCut",
136  "Cuts on kink mass for hyperon decay",
137  _hyperonDecayMassCut,
138  float(0.15));
139 
140  registerProcessorParameter("HyperonTimeCut",
141  "Cuts on lifetime for hyperon decay",
142  _hyperonTimeCut,
143  float(6.));
144 
145  registerProcessorParameter("KinkProjectionCutTPC",
146  "Cuts on kink distance of closest approach in TPC",
147  _drCutTPC,
148  float(20));
149 
150  registerProcessorParameter("TightKinkProjectionCutTPC",
151  "Tight Cut on kink distance of closest approach in TPC",
152  _tightDrCutTPC,
153  float(5));
154 
155  registerProcessorParameter("VeryTightKinkProjectionCutTPC",
156  "Very Tight Cut on kink distance of closest approach in TPC",
157  _veryTightDrCutTPC,
158  float(1));
159 
160  registerProcessorParameter("KinkProjectionCutSIT",
161  "Cuts on kink distance of closest approach in SIT",
162  _drCutSIT,
163  float(10));
164 
165  registerProcessorParameter("LooseProjectionCutSIT",
166  "Cuts on kink distance of closest approach in SIT",
167  _looseDrCutSIT,
168  float(10));
169 
170  registerProcessorParameter("MinELambda",
171  "Minimum Lambda Energy",
172  _minELambda,
173  float(2.));
174 
175 
176  registerProcessorParameter("DebugPrinting",
177  "Debug level",
178  _debugPrinting,
179  int(0));
180 
181 }
182 
183 void KinkFinder::init() {
184 
185  // Print algorithm parameters
186  printParameters() ;
187 
188 
189  dd4hep::Detector& theDet = dd4hep::Detector::getInstance();
190 
191  double bfieldV[3] ;
192  theDet.field().magneticField( { 0., 0., 0. } , bfieldV ) ;
193  _bField = bfieldV[2]/dd4hep::tesla ;
194 
195  dd4hep::DetElement tpcDE = theDet.detector("TPC") ;
196  const dd4hep::rec::FixedPadSizeTPCData* tpc = tpcDE.extension<dd4hep::rec::FixedPadSizeTPCData>() ;
197 
198  _tpcInnerR = tpc->rMinReadout/dd4hep::mm;
199  _tpcOuterR = tpc->rMaxReadout/dd4hep::mm;
200  _tpcZmax = tpc->driftLength/dd4hep::mm;
201  _tpcMaxRow = tpc->maxRow;
202 
203 
204  dd4hep::DetElement vtxDE = theDet.detector("VXD");
205  dd4hep::rec::ZPlanarData* vtx = vtxDE.extension<dd4hep::rec::ZPlanarData>();
206  _nLayersVTX=vtx->layers.size();
207 
208  for(int iL=0;iL<_nLayersVTX;iL++){
209  _rVTX.push_back( vtx->layers[iL].distanceSensitive/dd4hep::mm );
210  }
211 
212  dd4hep::DetElement sitDE = theDet.detector("SIT");
213  dd4hep::rec::ZPlanarData* sit = sitDE.extension<dd4hep::rec::ZPlanarData>();
214  _nLayersSIT=sit->layers.size();
215 
216  for(int iL=0;iL<_nLayersSIT;iL++){
217  _rSIT.push_back( sit->layers[iL].distanceSensitive/dd4hep::mm );
218  }
219 
220 
221  _nRun = -1;
222  _nEvt = 0;
223 
224 }
225 
226 
227 void KinkFinder::processRunHeader( LCRunHeader* ) {
228 
229  _nRun++ ;
230  _nEvt = 0;
231 
232 }
233 
234 void KinkFinder::processEvent( LCEvent * evt ) {
235 
236  TrackVec tracks;
237 
238 
239  if(_debugPrinting>0)std::cout << " Searching for kinks/prongs/split tracks " << std::endl;
240 
241 
242  LCCollection* relationTrackCollection = NULL;
243  LCRelationNavigator* trackNavigator=NULL;
244  try {
245  relationTrackCollection = evt->getCollection("LDCTracksMCP");
246  trackNavigator = new LCRelationNavigator(relationTrackCollection);
247  }
248  catch(DataNotAvailableException &e){
249  }
250 
251  try {
252  LCCollection * col = evt->getCollection( _trackColName.c_str() );
253 
254  int nelem = col->getNumberOfElements();
255  TrackPairVec trkPairs;
256  trkPairs.clear();
257  std::map<Track*,int> trackUsed;
258 
259  for (int i=0;i<nelem;++i) {
260  Track * trk = dynamic_cast<Track*>(col->getElementAt(i));
261  tracks.push_back(trk);
262  trackUsed[trk] = 0;
263  }
264  }
265  catch(DataNotAvailableException &e) {}
266 
267  std::vector<vec3>rin;
268  std::vector<vec3>rout;
269  std::vector<float>rOuter;
270  std::vector<float>rInner;
271  std::vector<vec3>zin;
272  std::vector<vec3>zout;
273  std::vector<float>momentum;
274  std::vector<float>momentumZ;
275  std::vector<float>momentumS;
276  std::vector<float>momentumE;
277  std::vector<float>charge;
278  std::vector<int>hits;
279  std::vector<MCParticle*>mcParticle;
280  std::vector<HelixClass*> helixEnd( tracks.size() );
281  std::vector<HelixClass*> helixStart(tracks.size() );
282  std::vector< std::vector<twoTrackIntersection_t> > kinkDaughters ( tracks.size() ) ;
283  std::vector< std::vector<twoTrackIntersection_t> > prongDaughters( tracks.size() ) ;
284  std::vector< std::vector<twoTrackIntersection_t> > splitDaughters( tracks.size() ) ;
285  std::vector< float> zAtEnd( tracks.size() );
286  std::vector< float> zAtStart(tracks.size() );
287 
288  for(unsigned int itrack=0;itrack< tracks.size();++itrack){
289  TrackerHitVec hitvec = tracks[itrack]->getTrackerHits();
290  int nhits = (int)hitvec.size();
291  hits.push_back(nhits);
292  int outer = 0;
293  float router =-99999.;
294  int inner = 0;
295  float rinner =99999.;
296  int max = 0;
297  int min = 0;
298  float abszmax =-99999.;
299  float abszmin =99999.;
300  float zmax =-99999.;
301  float zmin =99999.;
302  for(int ih =0;ih<nhits;++ih){
303  float x = (float)hitvec[ih]->getPosition()[0];
304  float y = (float)hitvec[ih]->getPosition()[1];
305  float z = (float)hitvec[ih]->getPosition()[2];
306  float r2 = x*x+y*y;
307  float r = sqrt(r2);
308  if(z<zmin)zmin=z;
309  if(z>zmax)zmax=z;
310  if(fabs(z)<abszmin){
311  min=ih;
312  abszmin = fabs(z);
313  }
314  if(fabs(z)>abszmax){
315  max=ih;
316  abszmax = fabs(z);
317  }
318  if(r<rinner && r>15.){
319  inner=ih;
320  rinner = r;
321  }
322  if(r>router){
323  outer=ih;
324  router = r;
325  }
326  }
327  vec3 ri;
328  ri.x = (float)hitvec[inner]->getPosition()[0];
329  ri.y = (float)hitvec[inner]->getPosition()[1];
330  ri.z = (float)hitvec[inner]->getPosition()[2];
331  rin.push_back(ri);
332  vec3 ro;
333  ro.x = (float)hitvec[outer]->getPosition()[0];
334  ro.y = (float)hitvec[outer]->getPosition()[1];
335  ro.z = (float)hitvec[outer]->getPosition()[2];
336  rout.push_back(ro);
337  vec3 zi;
338  zi.x = (float)hitvec[min]->getPosition()[0];
339  zi.y = (float)hitvec[min]->getPosition()[1];
340  zi.z = (float)hitvec[min]->getPosition()[2];
341  vec3 zo;
342  zo.x = (float)hitvec[max]->getPosition()[0];
343  zo.y = (float)hitvec[max]->getPosition()[1];
344  zo.z = (float)hitvec[max]->getPosition()[2];
345  zin.push_back(zi);
346  zout.push_back(zo);
347  rInner.push_back(rinner);
348  rOuter.push_back(router);
349  float zBegin, zEnd;
350  if (fabs(zmin)<fabs(zmax)) {
351  zBegin = zmin;
352  zEnd = zmax;
353  }else {
354  zBegin = zmax;
355  zEnd = zmin;
356  }
357  float signPz = zEnd - zBegin;
358  zAtEnd[itrack] = zEnd;
359  zAtStart[itrack] = zBegin;
360 
361  int _nHitsInFit = 50;
362  int nhitsFit;
363  if (nhits > _nHitsInFit) {
364  for (int iz = 0 ; iz < nhits-1; iz++)
365  for (int jz = 0; jz < nhits-iz-1; jz++){
366  TrackerHit * one = hitvec[jz];
367  TrackerHit * two = hitvec[jz+1];
368  float dist1 = fabs(float(one->getPosition()[2])-zBegin);
369  float dist2 = fabs(float(two->getPosition()[2])-zBegin);
370  if(dist1 > dist2)
371  {
372  TrackerHit * Temp = hitvec[jz];
373  hitvec[jz] = hitvec[jz+1];
374  hitvec[jz+1] = Temp;
375  }
376  }
377  nhitsFit = _nHitsInFit;
378  }else {
379  nhitsFit = nhits;
380  }
381 
382  HelixClass helix;
383  float phi0 = tracks[itrack]->getPhi();
384  float d0 = tracks[itrack]->getD0();
385  float omega = tracks[itrack]->getOmega();
386  float z0 = tracks[itrack]->getZ0();
387  float tanlambda = tracks[itrack]->getTanLambda();
388  helix.Initialize_Canonical(phi0, d0, z0, omega, tanlambda, _bField);
389  float totMomentum =0;
390  float Mom[3];
391  for (int j=0; j < 3; ++j) {
392  Mom[j] = helix.getMomentum()[j];
393  totMomentum += Mom[j]*Mom[j];
394  }
395  totMomentum = sqrt(totMomentum);
396  momentum.push_back(totMomentum);
397  momentumZ.push_back(Mom[2]);
398  charge.push_back(helix.getCharge());
399 
400  MCParticle* mcpart = NULL;
401  if(trackNavigator!=NULL){
402  LCObjectVec objectVec = trackNavigator->getRelatedToObjects(tracks[itrack]);
403  if (objectVec.size() > 0)mcpart = dynamic_cast<MCParticle*>(objectVec[0]);
404  if (objectVec.size() > 1){
405  float dbest = 999999999.;
406  for(unsigned im=0;im<objectVec.size();im++ ){
407  MCParticle* mcparti = dynamic_cast<MCParticle*>(objectVec[im]);
408  float p[3];
409  if(mcparti!=NULL){
410  p[0] = mcparti->getMomentum()[0];
411  p[1] = mcparti->getMomentum()[1];
412  p[2] = mcparti->getMomentum()[2];
413  float ptot = sqrt(p[0]*p[0]+p[1]*p[1]+p[2]*p[2]);
414  float d = fabs(ptot-totMomentum);
415  if(d<dbest){
416  dbest = d;
417  mcpart = mcparti;
418  }
419  }
420  }
421  }
422  }
423  mcParticle.push_back(mcpart);
424 
425  float * xh = new float[nhitsFit];
426  float * yh = new float[nhitsFit];
427  float * zh = new float[nhitsFit];
428  float * ah = new float[nhitsFit];
429  for (int i=0; i<nhitsFit; ++i) {
430  ah[i] = 0;
431  xh[i] = float(hitvec[i]->getPosition()[0]);
432  yh[i] = float(hitvec[i]->getPosition()[1]);
433  zh[i] = float(hitvec[i]->getPosition()[2]);
434  }
435 
436  ClusterShapes * shapesS = new ClusterShapes(nhitsFit,ah,xh,yh,zh);
437  float par[5];
438  float dpar[5];
439  float chi2;
440  float distmax;
441  shapesS->FitHelix(500, 0, 1, par, dpar, chi2, distmax);
442 
443  delete shapesS;
444  float x0 = par[0];
445  float y0 = par[1];
446  float r0 = par[2];
447  float bz = par[3];
448  float phiH = par[4];
449  helixStart[itrack] = new HelixClass();
450  helixStart[itrack]->Initialize_BZ(x0, y0, r0,
451  bz, phiH, _bField,signPz,
452  zBegin);
453 
454 
455  float seeds[3];
456  float refs[3];
457  refs[0] = helixStart[itrack]->getReferencePoint()[0];
458  refs[1] = helixStart[itrack]->getReferencePoint()[1];
459  refs[2] = helixStart[itrack]->getReferencePoint()[2];
460 
461  helixStart[itrack]->getPointInZ(zAtEnd[itrack], refs, seeds);
462  for (int i=0; i<nhitsFit; ++i) {
463  ah[i] = 0;
464  xh[i] = float(hitvec[nhits-1-i]->getPosition()[0]);
465  yh[i] = float(hitvec[nhits-1-i]->getPosition()[1]);
466  zh[i] = float(hitvec[nhits-1-i]->getPosition()[2]);
467  }
468  ClusterShapes * shapesE = new ClusterShapes(nhitsFit,ah,xh,yh,zh);
469  shapesE->FitHelix(500, 0, 1, par, dpar, chi2, distmax);
470  delete shapesE;
471  x0 = par[0];
472  y0 = par[1];
473  r0 = par[2];
474  bz = par[3];
475  phiH = par[4];
476  helixEnd[itrack] = new HelixClass();
477  helixEnd[itrack]->Initialize_BZ(x0, y0, r0,
478  bz, phiH, _bField,signPz,
479  zEnd);
480 
481  delete[] xh;
482  delete[] yh;
483  delete[] zh;
484  delete[] ah;
485 
486  }
487 
488  for(unsigned int i=0;i< tracks.size();++i){
489 // Track* tracki = tracks[i];
490 // float d0i = tracki->getD0();
491 // float z0i = tracki->getZ0();
492  if(hits[i]>=_minTrackHits){
493  float seedi[3];
494  float refe[3];
495  float z = zAtEnd[i];
496  refe[0] = helixEnd[i]->getReferencePoint()[0];
497  refe[1] = helixEnd[i]->getReferencePoint()[1];
498  refe[2] = helixEnd[i]->getReferencePoint()[2];
499  helixEnd[i]->getPointInZ(z, refe, seedi);
500  float rendi = sqrt(seedi[0]*seedi[0]+seedi[1]*seedi[1]);
501  for(unsigned int j=0;j< tracks.size();++j){
502  if(i!=j && hits[j]>_minTrackHits && momentum[j] < (1.0 + _maxSplitTrackFracDeltaP)*momentum[i]){
503 // Track* trackj = tracks[j];
504 // float d0j = trackj->getD0();
505 // float z0j = trackj->getZ0();
506  if(rInner[i]>_rKinkCut || rInner[j]>_rKinkCut){
507  float seedj[3];
508  float refs[3];
509  float deltaz;
510  float ddx;
511  float ddy;
512  float ddz;
513  float deltar;
514  float deltarxy;
515  int ip;
516  int id;
517  float z1;
518  float z2;
519  bool flipped = false;
520  if(fabs(zAtEnd[i]-zAtStart[j])<=fabs(zAtEnd[i]-zAtEnd[j])){
521  // normal kink where both tracks propagate outwards
522  ip = i;
523  id = j;
524  refs[0] = helixStart[j]->getReferencePoint()[0];
525  refs[1] = helixStart[j]->getReferencePoint()[1];
526  refs[2] = helixStart[j]->getReferencePoint()[2];
527  helixStart[j]->getPointInZ(z, refs, seedj);
528  deltaz = zAtEnd[i] - zAtStart[j];
529  ddx = (zout[i].x-zin[j].x);
530  ddy = (zout[i].y-zin[j].y);
531  ddz = (zout[i].z-zin[j].z);
532  deltar = sqrt(ddx*ddx+ddy*ddy+ddz*ddz);
533  deltarxy = fabs(rInner[j]-rOuter[i]);
534  z1 = zAtEnd[ip];
535  z2 = zAtStart[id];
536  }else{
537  // kink where daughter track propagates inwards
538  // or a V0 !!!
539  // or a backscattered track
540  flipped = true;
541  ip = j;
542  id = i;
543  refs[0] = helixEnd[j]->getReferencePoint()[0];
544  refs[1] = helixEnd[j]->getReferencePoint()[1];
545  refs[2] = helixEnd[j]->getReferencePoint()[2];
546  helixEnd[j]->getPointInZ(z, refs, seedj);
547  deltaz = zAtEnd[i] - zAtEnd[j];
548  ddx = (zout[i].x-zout[j].x);
549  ddy = (zout[i].y-zout[j].y);
550  ddz = (zout[i].z-zout[j].z);
551  deltar = sqrt(ddx*ddx+ddy*ddy+ddz*ddz);
552  deltarxy = fabs(rInner[i]-rOuter[j]);
553  z1 = zAtEnd[ip];
554  z2 = zAtEnd[id];
555  }
556  float dx = seedi[0]-seedj[0];
557  float dy = seedi[1]-seedj[1];
558  float dz = seedi[2]-seedj[2];
559  float dr = sqrt(dx*dx+dy*dy+dz*dz);
560 
561  float xkink = 0;
562  float ykink = 0;
563  float zkink = 0;
564  float rkink = 0;
565 
566 
567  bool mcKink = false;
568  if(mcParticle[j]!=NULL && _debugPrinting>0){
569  EVENT::MCParticleVec parents = mcParticle[j]->getParents();
570  if(parents.size()>0 && mcParticle[i]!=NULL){
571  for(unsigned im = 0; im<parents.size();im++){
572  if(parents[im]==mcParticle[i] && momentum[i]>2.0){
573  // if(abs(mcParticle[i]->getPDG())==211 &&
574  // abs(mcParticle[j]->getPDG())==13)mcKink = true;
575  if(abs(mcParticle[i]->getPDG())==321 &&
576  abs(mcParticle[j]->getPDG())==13)mcKink = true;
577  if(abs(mcParticle[i]->getPDG())==321 &&
578  abs(mcParticle[j]->getPDG())==211)mcKink = true;
579  if(abs(mcParticle[i]->getPDG())==3222 &&
580  abs(mcParticle[j]->getPDG())==211)mcKink = true;
581  if(abs(mcParticle[i]->getPDG())==3212 &&
582  abs(mcParticle[j]->getPDG())==211)mcKink = true;
583  if(abs(mcParticle[i]->getPDG())==3112 &&
584  abs(mcParticle[j]->getPDG())==211)mcKink = true;
585  }
586  }
587  }
588  }
589 
590 
591 
592  if(dr<100 || mcKink){
593  dr = 999999.;
594  float refsp[3];
595  float refsd[3];
596  refsp[0] = helixEnd[ip]->getReferencePoint()[0];
597  refsp[1] = helixEnd[ip]->getReferencePoint()[1];
598  refsp[2] = helixEnd[ip]->getReferencePoint()[2];
599  refsd[0] = helixStart[id]->getReferencePoint()[0];
600  refsd[1] = helixStart[id]->getReferencePoint()[1];
601  refsd[2] = helixStart[id]->getReferencePoint()[2];
602  float seedp[3];
603  float seedd[3];
604  float zs = z1;
605  float ze = z2;
606  if(z2 < z1){
607  zs = z2;
608  ze = z1;
609  }
610  float zstep = (ze-zs)/10.;
611  if(zstep<1)zstep=1;
612 
613  for(float zf = zs; zf < ze; zf+= zstep){
614  helixEnd[ip]->getPointInZ(zf, refsp, seedp);
615  helixStart[id]->getPointInZ(zf, refsd, seedd);
616  float dxf = seedp[0]-seedd[0];
617  float dyf = seedp[1]-seedd[1];
618  float dzf = seedp[2]-seedd[2];
619  float newdr = sqrt(dxf*dxf+dyf*dyf+dzf*dzf);
620  if(newdr<dr){
621  dr=newdr;
622  xkink = (seedp[0]+seedd[0])/2;
623  ykink = (seedp[1]+seedd[1])/2;
624  rkink = sqrt(xkink*xkink+ykink*ykink);
625  zkink = zf;
626  if(rkink<rOuter[i])rkink=rOuter[i];
627  if(rkink>rInner[j])rkink=rInner[j];
628  }
629  }
630  }
631 
632  bool ok = true;
633  if(fabs(deltaz)>200)ok=false;
634  if(fabs(deltaz)>100 && dr > 5.0)ok=false;
635 
636  float deltaRxyCut = -100;
637  float drCut = -100;
638 // bool goodRadialSep = false;
639  // require kink to be outside VTX detector
640  if(rkink > _rVTX[_nLayersVTX-1]){
641  // for TPC use
642  drCut = _drCutTPC;
643  deltaRxyCut= (_tpcOuterR-_tpcInnerR)/_tpcMaxRow*_maxDeltaTpcLayers;
644  if(dr<_tightDrCutTPC)deltaRxyCut=deltaRxyCut*1.5;
645  if(dr<_veryTightDrCutTPC)deltaRxyCut=deltaRxyCut/1.5*2.0;
646 
647  // for SIT
648  if(rkink<_tpcInnerR+deltaRxyCut){
649  drCut = _drCutSIT;
650  int iSitLayer=0;
651  for(int il=0;il<_nLayersSIT;il++){
652  if(rkink>_rSIT[il])iSitLayer=il+1;
653  }
654  // kink between VTX and SIT
655  if(iSitLayer==0)deltaRxyCut = _rSIT[0]-_rVTX[_nLayersVTX-1];
656  // kink between SIT and TPC
657  if(iSitLayer==_nLayersSIT){
658  deltaRxyCut = _tpcInnerR - _rSIT[_nLayersSIT-1];
659  drCut = _drCutTPC;
660  }
661  if(iSitLayer>0 && iSitLayer <_nLayersSIT)deltaRxyCut = _rSIT[iSitLayer] - _rSIT[iSitLayer-1];
662  // add some protection for SIT layers - large gap
663  int ili = -999;
664  int ilo = -999;
665  for(int il=0; il<_nLayersSIT;il++){
666  if(fabs(rOuter[i]-_rSIT[il])<10.)ili=il;
667  if(fabs(rInner[j]-_rSIT[il])<10.)ilo=il;
668  }
669  if(ili>=0&&ilo>=0){
670  float fix = (_rSIT[ilo] - _rSIT[ili]);
671  if(fix<10)fix=10.;
672  if(fix>deltaRxyCut)deltaRxyCut = fix;
673  }
674  // add some slop to account for geometry
675  deltaRxyCut = deltaRxyCut*1.1;
676  }
677  }
678 
679 
680  // looser cuts for debug
681  // std::cout << i << " : " << j << " dr = " << dr << " ( " << drCut << " ) deltaRxy = " << deltarxy << " ( " << deltaRxyCut << " ) " << std::endl;
682  if( (dr<drCut && deltarxy < deltaRxyCut*2) || mcKink){
683  bool possibleSplit = false;
684  bool split = false;
685  rkink = sqrt(xkink*xkink+ykink*ykink);
686 
687  if( (rkink > _rKinkCut && !flipped) || mcKink){
688 
689  float massENu = this->kinkMass(helixEnd[i],helixStart[j],0., 0.);
690  float massMuNu = this->kinkMass(helixEnd[i],helixStart[j],mMuon,0.);
691  float massPiPi = this->kinkMass(helixEnd[i],helixStart[j],mPion,mPion);
692  float massPiN = this->kinkMass(helixEnd[i],helixStart[j],mPion,mNeutron);
693  float massPPi0 = this->kinkMass(helixEnd[i],helixStart[j],mProton,mPion);
694  float massPiL = this->kinkMass(helixEnd[i],helixStart[j],mPion,mLambda);
695  float FmassENu = this->kinkMass(helixEnd[i],helixEnd[j],0., 0.);
696  float FmassMuNu = this->kinkMass(helixEnd[i],helixEnd[j],mMuon,0.);
697  float FmassPiPi = this->kinkMass(helixEnd[i],helixEnd[j],mPion,mPion);
698  float FmassPiN = this->kinkMass(helixEnd[i],helixEnd[j],mPion,mNeutron);
699  float FmassPPi0 = this->kinkMass(helixEnd[i],helixEnd[j],mProton,mPion);
700  float FmassPiL = this->kinkMass(helixEnd[i],helixEnd[j],mPion,mLambda);
701  float tPion = fabs(zkink-zAtStart[i])*mPion/fabs(momentumZ[i])/cTauPion;
702  float tSigma = fabs(zkink-zAtStart[i])*mSigma/fabs(momentumZ[i])/cTauSigma;
703  float tKaon = fabs(zkink-zAtStart[i])*mKaon/fabs(momentumZ[i])/cTauKaon;
704  float tHyperon = fabs(zkink-zAtStart[i])*mHyperon/fabs(momentumZ[i])/cTauHyperon;
705  float probPion = 0;
706  float probKaon = 0;
707  float probSigma = 0;
708  float probHyperon = 0;
709  if(_pionDecayMassCut>0){
710  float deltaPion = fabs(massMuNu-mPion)/_pionDecayMassCut;
711  if(deltaPion<1 && tPion > 0.001)probPion = 3.125*deltaPion*deltaPion+tPion;
712  }
713  if(_kaonDecayMassCut>0){
714  float deltaKaonMuNu = fabs(massMuNu-mKaon)/_kaonDecayMassCut;
715  float deltaKaonPiPi = fabs(massPiPi-mKaon)/_kaonDecayMassCut;
716  float deltaKaon = deltaKaonMuNu;
717  if(deltaKaonPiPi<deltaKaon)deltaKaon = deltaKaonPiPi;
718  if(deltaKaon<1 && tKaon > 0.005)probKaon = 3.125*deltaKaon*deltaKaon+tKaon;
719  }
720  if(_sigmaDecayMassCut>0 && momentum[i] > _minELambda){
721  float deltaSigmaPiN = fabs(massPiN-mSigma)/_sigmaDecayMassCut;
722  float deltaSigmaPPi0 = fabs(massPPi0-mSigma)/_sigmaDecayMassCut;
723  float deltaSigma = deltaSigmaPiN;
724  if(deltaSigmaPPi0<deltaSigma)deltaSigma = deltaSigmaPPi0;
725  if(deltaSigma<1 && tSigma < _sigmaTimeCut)probSigma = 3.125*deltaSigma*deltaSigma+tSigma;
726  }
727  if(_hyperonDecayMassCut>0 && momentum[i] > _minELambda){
728  float deltaHyperon = fabs(massPiL-mHyperon)/_hyperonDecayMassCut;
729  if(deltaHyperon<1 && tHyperon < _hyperonTimeCut)probHyperon = 3.125*deltaHyperon*deltaHyperon+tHyperon;
730  }
731  bool goodKinkMass = false;
732  if(probPion>0.0001||probKaon>0.0001||probSigma>0.0001||probHyperon>0.0001)goodKinkMass=true;
733 
734 
735  float dpop = 2*fabs(momentum[i]-momentum[j])/(momentum[i]+momentum[j]);
736  float dp = fabs(momentum[i]-momentum[j]);
737  if(dr<drCut && dpop < _maxSplitTrackFracDeltaP && dp < _maxSplitTrackDeltaP){
738  possibleSplit = true;
739  }
740 
741  if( (deltarxy<deltaRxyCut && dr<drCut) || possibleSplit ){
742  twoTrackIntersection_t kinkij;
743  kinkij.tracki = i;
744  kinkij.trackj = j;
745  kinkij.vtx[0] = xkink;
746  kinkij.vtx[1] = ykink;
747  kinkij.vtx[2] = zkink;
748  kinkij.p[0] = helixStart[i]->getMomentum()[0];
749  kinkij.p[1] = helixStart[i]->getMomentum()[1];
750  kinkij.p[2] = helixStart[i]->getMomentum()[2];
751  kinkij.mass = 0.;
752  kinkij.distance = dr;
753 
754  // split tracks
755  if(deltarxy<2*deltaRxyCut && dr<drCut*2 ){
756  if(possibleSplit && charge[i]*charge[j]>0 ){
757  TrackerHitVec hitveci= tracks[i]->getTrackerHits();
758  TrackerHitVec hitvecj= tracks[j]->getTrackerHits();
759  int nhitsi = (int)hitveci.size();
760  int nhitsj = (int)hitvecj.size();
761  int ntpci = 0;
762  int ntpcj = 0;
763  HelixClass* helixi = helixStart[j];
764  HelixClass* helixj = helixEnd[i];
765  float hitxyz[3];
766  float dist[3];
767  float maxdisti=0;
768  float maxdistj=0;
769  int nclosei = 0;
770  int nclosej = 0;
771  float zmini = 99999;
772  float zminj = 99999;
773  float zmaxi = -99999;
774  float zmaxj = -99999;
775  for(int ih =0;ih<nhitsi;++ih){
776  float x = (float)hitveci[ih]->getPosition()[0];
777  float y = (float)hitveci[ih]->getPosition()[1];
778  float zz = (float)hitveci[ih]->getPosition()[2];
779  if(fabs(zz)<zmini)zmini=fabs(zz);
780  if(fabs(zz)>zmaxi)zmaxi=fabs(zz);
781  float r2 = x*x+y*y;
782  float r = sqrt(r2);
783  if(r>_tpcInnerR)ntpci++;
784  hitxyz[0]=x;
785  hitxyz[1]=y;
786  hitxyz[2]=zz;
787  helixj->getDistanceToPoint(hitxyz, dist);
788  if(dist[2]>maxdisti)maxdisti=dist[2];
789  if(dist[2]<25.)nclosei++;
790  }
791  for(int ih =0;ih<nhitsj;++ih){
792  float x = (float)hitvecj[ih]->getPosition()[0];
793  float y = (float)hitvecj[ih]->getPosition()[1];
794  float zz = (float)hitvecj[ih]->getPosition()[2];
795  if(fabs(zz)<zminj)zminj=fabs(zz);
796  if(fabs(zz)>zmaxj)zmaxj=fabs(zz);
797  float r2 = x*x+y*y;
798  float r = sqrt(r2);
799  if(r>_tpcInnerR)ntpcj++;
800  hitxyz[0]=x;
801  hitxyz[1]=y;
802  hitxyz[2]=zz;
803  helixi->getDistanceToPoint(hitxyz, dist);
804  if(dist[2]>maxdistj)maxdistj=dist[2];
805  if(dist[2]<25.)nclosej++;
806  }
807  float fclosei = (float)nclosei/(float)nhitsi;
808  float fclosej = (float)nclosej/(float)nhitsj;
809  if(_debugPrinting>0){
810  std::cout << " CAND SPLIT I : " << nhitsi << " ntpc " << ntpci << " nclose " << nclosei << " max " << maxdisti << " fclose : " << fclosei << std::endl;
811  std::cout << " CAND SPLIT J : " << nhitsj << " ntpc " << ntpcj << " nclose " << nclosej << " max " << maxdistj << " fclose : " << fclosej << std::endl;
812  }
813  if(maxdistj<50 && maxdisti < 50 && fclosej > 0.95 && fclosei > 0.95 && ntpcj+ntpci < _tpcMaxRow+10.)split = true;
814  splitDaughters[i].push_back(kinkij);
815  }
816  }
817 
818  // kinks
819  if(deltarxy<deltaRxyCut && dr<drCut ){
820  if(charge[i]*charge[j]>0 && goodKinkMass){
821  if(_debugPrinting>0)std::cout << "Found kink candidate " << i << " " << j << std::endl;
822  if(probPion>0||probKaon>0||probSigma>0||probHyperon>0)goodKinkMass=true;
823  if(probPion > probKaon &&
824  probPion > probSigma &&
825  probPion > probHyperon){
826  kinkij.pdgCode = 211;
827  if(charge[i]<0)kinkij.pdgCode = -211;
828  kinkij.mass = mPion;
829  }
830  if(probKaon > probPion &&
831  probKaon > probSigma &&
832  probKaon > probHyperon){
833  kinkij.mass = mKaon;
834  kinkij.pdgCode = 321;
835  if(charge[i]<0)kinkij.pdgCode = -321;
836  }
837  if(probSigma > probPion &&
838  probSigma > probKaon &&
839  probSigma > probHyperon){
840  kinkij.mass = mSigma;
841  kinkij.pdgCode = 3222;
842  if(charge[i]<0)kinkij.pdgCode = 3222;
843  }
844  if(probHyperon > probPion &&
845  probHyperon > probKaon &&
846  probHyperon > probSigma){
847  kinkij.mass = mHyperon;
848  kinkij.pdgCode = 3312;
849  if(charge[i]>0)kinkij.pdgCode = -3312;
850  }
851  kinkDaughters[i].push_back(kinkij);
852  }
853  }
854 
855  // prongs
856  if(deltarxy<deltaRxyCut && dr<drCut ){
857 
858  if(_debugPrinting>0)std::cout << "Found prong candidate " << i << " " << j << std::endl;
859  prongDaughters[i].push_back(kinkij);
860  kinkij.pdgCode = 211;
861  if(charge[i]<0)kinkij.pdgCode = -211;
862  kinkij.mass = mPion;
863  }
864 
865  }
866 
867  if(_debugPrinting>0){
868  if(mcParticle[j]!=NULL){
869  EVENT::MCParticleVec parents = mcParticle[j]->getParents();
870  if(parents.size()>0 && mcParticle[i]!=NULL){
871  for(unsigned im = 0; im<parents.size();im++){
872  if(parents[im]==mcParticle[i])std::cout << " TRUE KINK . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . " << std::endl;
873  }
874  }
875  }
876 
877 
878 
879  std::cout << " Candidate kink tracks : " << i << "," << j << " p[i] " << momentum[i] << " p[j] " << momentum[j] << std::endl;
880  if(flipped)std::cout << " Flipped Track " << std::endl;;
881  std::cout << " MC : ";
882  if(mcParticle[i]!=NULL)std::cout << mcParticle[i]->getPDG();
883  std::cout << " - ";
884  if(mcParticle[j]!=NULL)std::cout << mcParticle[j]->getPDG();
885  std::cout << std::endl;
886  std::cout << " Mass : " << massENu << " " << massMuNu << " " << massPiPi << " " << massPiN << " " << massPPi0 << " " << massPiL << std::endl;
887  if(flipped)std::cout << " MassF : " << FmassENu << " " << FmassMuNu << " " << FmassPiPi << " " << FmassPiN << " " << FmassPPi0 << FmassPiL << std::endl;
888 
889  std::cout << " ProbPion : " << probPion << " TimePion : " << tPion << " Dm = " << massMuNu-mPion << std::endl;
890  std::cout << " ProbKaon : " << probKaon << " TimeKaon : " << tKaon <<
891  " Dm = " << massMuNu-mKaon <<
892  " Dm = " << massPiPi-mKaon << std::endl;
893  std::cout << " ProbSigma : " << probSigma << " TimeSigma : " << tSigma << " Dm = " << massPiN-mSigma <<
894  " Dm = " << massPPi0-mSigma << std::endl;
895  std::cout << " ProbHyperon : " << probHyperon << " TimeHyperon : " << tHyperon << " Dm = " << massPiL-mHyperon << std::endl;
896  std::cout << " Charge : " << i << "," << j << " q[i] " << charge[i] << " q[j] " << charge[j] << std::endl;
897  std::cout << " Rendi : " << i << "," << j << " rend " << rendi << " dr = " << dr << " (deltaR = " << deltar << ")" << std::endl;
898  std::cout << " zkink : " << zkink << " rzykink " << rkink << std::endl;
899  if(deltarxy<deltaRxyCut)std::cout << " DeltaRxy vs cut : " << deltarxy << " < " << deltaRxyCut << std::endl;
900  if(deltarxy>deltaRxyCut)std::cout << " DeltaRxy vs cut : " << deltarxy << " > " << deltaRxyCut << std::endl;
901  if(dr<drCut)std::cout << " dr vs cut : " << dr << " < " << drCut << std::endl;
902  if(dr>drCut)std::cout << " dr vs cut : " << dr << " > " << drCut << std::endl;
903  std::cout << " zi/zj : " << zAtStart[i] << " - " << zAtEnd[i] << " " << zAtStart[j] << " - " << zAtEnd[j] << std::endl;
904  std::cout << " rinner/router : " << rInner[i] << " - " << rOuter[i] << " " << rInner[j] << " - " << rOuter[j] << std::endl;
905  }
906  }
907  }
908  }
909  }
910  }
911  }
912  }
913 
914  int countk = 0;
915  int countp = 0;
916  for(unsigned int i=0;i< tracks.size();++i){
917 
918 
919  if(kinkDaughters[i].size()>0 || prongDaughters[i].size()>0){
920  if(kinkDaughters[i].size()==1 && prongDaughters[i].size()<2){
921  countk++;
922  }else{
923  countp++;
924  }
925  }
926  }
927  LCCollectionVec * colKinkRecoPart = NULL;
928  LCCollectionVec * colKinkVertex = NULL;
929  LCCollectionVec * colProngRecoPart = NULL;
930  LCCollectionVec * colProngVertex = NULL;
931  LCCollectionVec * colSplitRecoPart = NULL;
932  LCCollectionVec * colSplitVertex = NULL;
933 
934  for(unsigned int i=0;i< tracks.size();++i){
935  if(kinkDaughters[i].size()>0 || prongDaughters[i].size()>0 || splitDaughters[i].size()>0){
936  if(_debugPrinting>0){
937  if(kinkDaughters[i].size()>0){
938  std::cout << " Track " << i << " has " << kinkDaughters[i].size() << " kink daughters : ";
939  for(unsigned j=0;j<kinkDaughters[i].size();j++)std::cout << kinkDaughters[i][j].trackj << " ";
940  std::cout << std::endl;
941  }
942  if(prongDaughters[i].size()>0){
943  std::cout << " Track " << i << " has " << prongDaughters[i].size() << " prong daughters : ";
944  for(unsigned j=0;j<prongDaughters[i].size();j++)std::cout << prongDaughters[i][j].trackj << " ";
945  std::cout << std::endl;
946  }
947  if(splitDaughters[i].size()>0){
948  std::cout << " Track " << i << " has " << splitDaughters[i].size() << " split daughters : ";
949  for(unsigned j=0;j<splitDaughters[i].size();j++)std::cout << splitDaughters[i][j].trackj << " ";
950  std::cout << std::endl;
951  }
952  }
953 
954  // Split track collection
955  if(splitDaughters[i].size()==1){
956  if(_debugPrinting>0)std::cout << " Reconstructed SPLIT " << std::endl;
957  ReconstructedParticleImpl * part = new ReconstructedParticleImpl();
958  VertexImpl * vtx = new VertexImpl();
959  float vertex[3];
960  float mom[3];
961  int code = 211;
962  for (int iC=0;iC<3;++iC) {
963  vertex[iC] = splitDaughters[i][0].vtx[iC];
964  mom[iC] = splitDaughters[i][0].p[iC];
965  }
966  float distance = splitDaughters[i][0].distance;
967  vtx->setPosition( vertex );
968  vtx->addParameter( distance );
969  part->setMomentum( mom );
970  part->setType( code );
971 
972  if(_debugPrinting>0){
973  std::cout << " Code = " << code << " Distance = " << distance << std::endl;
974  std::cout << " Vertex = ("
975  << vertex[0] << ","
976  << vertex[1] << ","
977  << vertex[2] << ")" << std::endl;
978  std::cout << " Momentum = ("
979  << mom[0] << ","
980  << mom[1] << ","
981  << mom[2] << ")" << std::endl;
982  }
983  float mass = splitDaughters[i][0].mass;
984  part->setMass( mass );
985  vtx->setAssociatedParticle( part );
986  part->setStartVertex( vtx );
987  part->addTrack( tracks[splitDaughters[i][0].tracki] );
988  part->addTrack( tracks[splitDaughters[i][0].trackj] );
989  if(colSplitRecoPart==NULL){
990  colSplitRecoPart = new LCCollectionVec(LCIO::RECONSTRUCTEDPARTICLE);
991  colSplitVertex = new LCCollectionVec(LCIO::VERTEX);
992  }
993  colSplitRecoPart->addElement( part );
994  colSplitVertex->addElement( vtx );
995  }
996 
997  // Kink collection
998  if(splitDaughters[i].size()!=1){
999  if(kinkDaughters[i].size()==1 && prongDaughters[i].size()<2){
1000  if(_debugPrinting>0)std::cout << " Reconstructed KINKS " << std::endl;
1001  ReconstructedParticleImpl * part = new ReconstructedParticleImpl();
1002  VertexImpl * vtx = new VertexImpl();
1003  float vertex[3];
1004  float mom[3];
1005  int code = kinkDaughters[i][0].pdgCode;
1006  for (int iC=0;iC<3;++iC) {
1007  vertex[iC] = kinkDaughters[i][0].vtx[iC];
1008  mom[iC] = kinkDaughters[i][0].p[iC];
1009  }
1010  float distance = kinkDaughters[i][0].distance;
1011  vtx->setPosition( vertex );
1012  vtx->addParameter( distance );
1013  part->setMomentum( mom );
1014  part->setType( code );
1015  if(_debugPrinting>0){
1016  std::cout << " Code = " << code << " Distance = " << distance << std::endl;
1017  std::cout << " Vertex = ("
1018  << vertex[0] << ","
1019  << vertex[1] << ","
1020  << vertex[2] << ")" << std::endl;
1021  std::cout << " Momentum = ("
1022  << mom[0] << ","
1023  << mom[1] << ","
1024  << mom[2] << ")" << std::endl;
1025  }
1026  float mass = kinkDaughters[i][0].mass;
1027  part->setMass( mass );
1028  vtx->setAssociatedParticle( part );
1029  part->setStartVertex( vtx );
1030  part->addTrack( tracks[kinkDaughters[i][0].tracki] );
1031  part->addTrack( tracks[kinkDaughters[i][0].trackj] );
1032  if(colKinkRecoPart==NULL){
1033  colKinkRecoPart = new LCCollectionVec(LCIO::RECONSTRUCTEDPARTICLE);
1034  colKinkVertex = new LCCollectionVec(LCIO::VERTEX);
1035  }
1036  colKinkRecoPart->addElement( part );
1037  colKinkVertex->addElement( vtx );
1038  // trackUsed[firstTrack] = 1;
1039  // trackUsed[secondTrack] = 1;
1040  } else if(prongDaughters[i].size() > 0 ) {
1041  // Prong track collection
1042  if(_debugPrinting>0)std::cout << " Reconstructed PRONG " << std::endl;
1043  ReconstructedParticleImpl * part = new ReconstructedParticleImpl();
1044  VertexImpl * vtx = new VertexImpl();
1045  float vertex[3];
1046  float mom[3];
1047  int code = 211;
1048  for (int iC=0;iC<3;++iC) {
1049  vertex[iC] = prongDaughters[i][0].vtx[iC];
1050  mom[iC] = prongDaughters[i][0].p[iC];
1051  }
1052  float distance = prongDaughters[i][0].distance;
1053  vtx->setPosition( vertex );
1054  vtx->addParameter( distance );
1055  part->setMomentum( mom );
1056  part->setType( code );
1057  if(_debugPrinting>0){
1058  std::cout << " Code = " << code << " Distance = " << distance << std::endl;
1059  std::cout << " Vertex = ("
1060  << vertex[0] << ","
1061  << vertex[1] << ","
1062  << vertex[2] << ")" << std::endl;
1063  std::cout << " Momentum = ("
1064  << mom[0] << ","
1065  << mom[1] << ","
1066  << mom[2] << ")" << std::endl;
1067  }
1068  float mass = prongDaughters[i][0].mass;
1069  part->setMass( mass );
1070  vtx->setAssociatedParticle( part );
1071  part->setStartVertex( vtx );
1072  part->addTrack( tracks[prongDaughters[i][0].tracki] );
1073  for(unsigned id=0;id<prongDaughters[i].size();id++)part->addTrack( tracks[prongDaughters[i][id].trackj] );
1074  if(colProngRecoPart==NULL){
1075  colProngRecoPart = new LCCollectionVec(LCIO::RECONSTRUCTEDPARTICLE);
1076  colProngVertex = new LCCollectionVec(LCIO::VERTEX);
1077  }
1078  colProngRecoPart->addElement( part );
1079  colProngVertex->addElement( vtx );
1080  }
1081  // trackUsed[firstTrack] = 1;
1082  // trackUsed[secondTrack] = 1;
1083  }
1084  }
1085 
1086  }
1087 
1088 
1089  if(colKinkRecoPart!=NULL){
1090  evt->addCollection(colKinkRecoPart, _kinkRecoPartColName.c_str() );
1091  evt->addCollection(colKinkVertex, _kinkVertexColName.c_str() );
1092  }
1093  if(colProngRecoPart!=NULL){
1094  evt->addCollection(colProngRecoPart, _prongRecoPartColName.c_str() );
1095  evt->addCollection(colProngVertex, _prongVertexColName.c_str() );
1096  }
1097  if(colSplitRecoPart!=NULL){
1098  evt->addCollection(colSplitRecoPart, _splitRecoPartColName.c_str() );
1099  evt->addCollection(colSplitVertex, _splitVertexColName.c_str() );
1100  }
1101 
1102  for(unsigned int itrack=0;itrack< tracks.size();++itrack){
1103  delete helixEnd[itrack];
1104  delete helixStart[itrack];
1105  }
1106 
1107  _nEvt++;
1108  if(_debugPrinting>0)std::cout << " Done " << std::endl;
1109  if(trackNavigator != NULL) delete trackNavigator;
1110 }
1111 
1112 
1113 void KinkFinder::check( LCEvent * ) { }
1114 
1115 void KinkFinder::end(){ }
1116 
1117 void KinkFinder::Sorting( TrackPairVec & trkPairVec ) {
1118 
1119  int sizeOfVector = int(trkPairVec.size());
1120  TrackPair *one,*two,*Temp;
1121 
1122  for (int i = 0 ; i < sizeOfVector-1; i++)
1123  for (int j = 0; j < sizeOfVector-i-1; j++)
1124  {
1125  one = trkPairVec[j];
1126  two = trkPairVec[j+1];
1127  if( one->getDistance() > two->getDistance() )
1128  {
1129  Temp = trkPairVec[j];
1130  trkPairVec[j] = trkPairVec[j+1];
1131  trkPairVec[j+1] = Temp;
1132  }
1133  }
1134 
1135 }
1136 
1137 float KinkFinder::kinkMass(HelixClass* parent, HelixClass* daughter, float daughterMass, float neutralMass){
1138 
1139  // Calculate the invariant mass for a decaying charged particle
1140 
1141  float parentMom[4];
1142  parentMom[0] = parent->getMomentum()[0];
1143  parentMom[1] = parent->getMomentum()[1];
1144  parentMom[2] = parent->getMomentum()[2];
1145  float daughterMom[4];
1146  daughterMom[0] = daughter->getMomentum()[0];
1147  daughterMom[1] = daughter->getMomentum()[1];
1148  daughterMom[2] = daughter->getMomentum()[2];
1149  float daughterE = sqrt(daughterMom[0]*daughterMom[0] + daughterMom[1]*daughterMom[1]+
1150  daughterMom[2]*daughterMom[2] + daughterMass*daughterMass);
1151 
1152  float pnu[3];
1153  pnu[0] = parentMom[0]-daughterMom[0];
1154  pnu[1] = parentMom[1]-daughterMom[1];
1155  pnu[2] = parentMom[2]-daughterMom[2];
1156 
1157  float Enu = sqrt(pnu[0]*pnu[0]+pnu[1]*pnu[1]+pnu[2]*pnu[2]+neutralMass*neutralMass);
1158  float mx = sqrt((daughterE+Enu)*(daughterE+Enu)-parentMom[0]*parentMom[0]-
1159  parentMom[1]*parentMom[1]-parentMom[2]*parentMom[2]);
1160  return mx;
1161 
1162 }
1163 
KinkFinder::_rVTX
std::vector< float > _rVTX
Definition: KinkFinder.h:180
KinkFinder::init
virtual void init()
Definition: KinkFinder.cc:183
sistrip::mm
static const float mm
Definition: PhysicalConstants.h:37
KinkFinder::_drCutTPC
float _drCutTPC
Definition: KinkFinder.h:164
twoTrackIntersection_t::distance
float distance
Definition: KinkFinder.h:105
KinkFinder
KinkFinder Processor KinkFinder processor identify kinked tracks
Definition: KinkFinder.h:113
KinkFinder::_veryTightDrCutTPC
float _veryTightDrCutTPC
Definition: KinkFinder.h:163
KinkFinder::_kinkVertexColName
std::string _kinkVertexColName
Definition: KinkFinder.h:144
KinkFinder::end
virtual void end()
Definition: KinkFinder.cc:1115
vec3::x
float x
Definition: KinkFinder.h:94
twoTrackIntersection_t::trackj
int trackj
Definition: KinkFinder.h:101
KinkFinder::_looseDrCutSIT
float _looseDrCutSIT
Definition: KinkFinder.h:166
KinkFinder::_tpcMaxRow
int _tpcMaxRow
Definition: KinkFinder.h:176
KinkFinder::_sigmaTimeCut
float _sigmaTimeCut
Definition: KinkFinder.h:160
KinkFinder::_minELambda
float _minELambda
Definition: KinkFinder.h:169
KinkFinder::_tpcInnerR
float _tpcInnerR
Definition: KinkFinder.h:173
KinkFinder::_kaonDecayMassCut
float _kaonDecayMassCut
Definition: KinkFinder.h:156
twoTrackIntersection_t::pdgCode
int pdgCode
Definition: KinkFinder.h:106
KinkFinder::_bField
float _bField
Definition: KinkFinder.h:172
KinkFinder::_tightDrCutTPC
float _tightDrCutTPC
Definition: KinkFinder.h:162
vec3::y
float y
Definition: KinkFinder.h:95
mLambda
const float mLambda
Definition: KinkFinder.cc:31
aKinkFinder
KinkFinder aKinkFinder
Definition: KinkFinder.cc:38
KinkFinder::_nLayersSIT
int _nLayersSIT
Definition: KinkFinder.h:177
cTauSigma
const float cTauSigma
Definition: KinkFinder.cc:30
KinkFinder::_tpcOuterR
float _tpcOuterR
Definition: KinkFinder.h:174
twoTrackIntersection_t
Definition: KinkFinder.h:99
mNeutron
const float mNeutron
Definition: KinkFinder.cc:35
cTauPion
const float cTauPion
Definition: KinkFinder.cc:26
KinkFinder::_nLayersVTX
int _nLayersVTX
Definition: KinkFinder.h:178
mProton
const float mProton
Definition: KinkFinder.cc:36
mMuon
const float mMuon
Definition: KinkFinder.cc:23
KinkFinder::_drCutSIT
float _drCutSIT
Definition: KinkFinder.h:165
KinkFinder::_prongRecoPartColName
std::string _prongRecoPartColName
Definition: KinkFinder.h:149
KinkFinder::Sorting
void Sorting(TrackPairVec &trkPairVec)
Definition: KinkFinder.cc:1117
KinkFinder::_splitRecoPartColName
std::string _splitRecoPartColName
Definition: KinkFinder.h:150
KinkFinder::_nRun
int _nRun
Definition: KinkFinder.h:139
mHyperon
const float mHyperon
Definition: KinkFinder.cc:33
cTauKaon
const float cTauKaon
Definition: KinkFinder.cc:28
KinkFinder::check
virtual void check(LCEvent *evt)
Definition: KinkFinder.cc:1113
twoTrackIntersection_t::p
float p[3]
Definition: KinkFinder.h:103
KinkFinder::_prongVertexColName
std::string _prongVertexColName
Definition: KinkFinder.h:145
KinkFinder::_maxSplitTrackDeltaP
float _maxSplitTrackDeltaP
Definition: KinkFinder.h:168
KinkFinder::_nEvt
int _nEvt
Definition: KinkFinder.h:140
vec3::z
float z
Definition: KinkFinder.h:96
sistrip::e
static const float e
Definition: PhysicalConstants.h:26
KinkFinder::_rKinkCut
float _rKinkCut
Definition: KinkFinder.h:152
cTauHyperon
const float cTauHyperon
Definition: KinkFinder.cc:34
KinkFinder::_splitVertexColName
std::string _splitVertexColName
Definition: KinkFinder.h:146
KinkFinder::_maxSplitTrackFracDeltaP
float _maxSplitTrackFracDeltaP
Definition: KinkFinder.h:167
KinkFinder::_tpcZmax
float _tpcZmax
Definition: KinkFinder.h:175
KinkFinder::_rSIT
std::vector< float > _rSIT
Definition: KinkFinder.h:179
sistrip::LCCollectionVec
std::vector< LCCollection * > LCCollectionVec
Definition: SiStripClus.h:55
KinkFinder::_hyperonDecayMassCut
float _hyperonDecayMassCut
Definition: KinkFinder.h:159
KinkFinder::_debugPrinting
int _debugPrinting
Definition: KinkFinder.h:155
mSigma
const float mSigma
Definition: KinkFinder.cc:29
KinkFinder::_kinkRecoPartColName
std::string _kinkRecoPartColName
Definition: KinkFinder.h:148
KinkFinder::KinkFinder
KinkFinder()
Definition: KinkFinder.cc:41
KinkFinder::_minTrackHits
int _minTrackHits
Definition: KinkFinder.h:153
KinkFinder::_pionDecayMassCut
float _pionDecayMassCut
Definition: KinkFinder.h:157
twoTrackIntersection_t::mass
float mass
Definition: KinkFinder.h:104
KinkFinder::_maxDeltaTpcLayers
int _maxDeltaTpcLayers
Definition: KinkFinder.h:154
KinkFinder::processRunHeader
virtual void processRunHeader(LCRunHeader *run)
Definition: KinkFinder.cc:227
KinkFinder::_sigmaDecayMassCut
float _sigmaDecayMassCut
Definition: KinkFinder.h:158
KinkFinder::_trackColName
std::string _trackColName
Definition: KinkFinder.h:142
mKaon
const float mKaon
Definition: KinkFinder.cc:27
twoTrackIntersection_t::tracki
int tracki
Definition: KinkFinder.h:100
KinkFinder::processEvent
virtual void processEvent(LCEvent *evt)
Definition: KinkFinder.cc:234
KinkFinder::kinkMass
float kinkMass(HelixClass *parent, HelixClass *daughter, float daughterMass, float neutralMass)
Definition: KinkFinder.cc:1137
KinkFinder::_hyperonTimeCut
float _hyperonTimeCut
Definition: KinkFinder.h:161
mPion
const float mPion
Definition: KinkFinder.cc:25
twoTrackIntersection_t::vtx
float vtx[3]
Definition: KinkFinder.h:102