LCFIVertex  0.7.2
RPCutProcessor.cc
1 #include "RPCutProcessor.h"
2 #include <iostream>
3 
4 #include <EVENT/LCCollection.h>
5 #include <IMPL/LCCollectionVec.h>
6 #include <EVENT/ReconstructedParticle.h>
7 #include <IMPL/ReconstructedParticleImpl.h>
8 #include <EVENT/Track.h>
9 #include <EVENT/MCParticle.h>
10 #include <UTIL/LCRelationNavigator.h>
11 #include <EVENT/LCObject.h>
12 #include "../vertex_lcfi/inc/lciointerface.h"
13 #include "util/inc/memorymanager.h"
14 
15 #include <marlin/Global.h>
16 #include <gear/GEAR.h>
17 #include <gearimpl/Vector3D.h>
18 
19 #include <vector>
20 #include <string>
21 
22 using namespace marlin ;
23 using namespace lcio;
24 
25 RPCutProcessor aRPCutProcessor ;
26 
27 RPCutProcessor::RPCutProcessor() : Processor("RPCutProcessor") {
28 
29  // modify processor description
30  _description = "RPCutProcessor - cuts RPs based on several criteria removing those that have no tracks" ;
31 
32 
33  // register steering parameters: name, description, class-variable, default value
34 
35  registerInputCollection( lcio::LCIO::RECONSTRUCTEDPARTICLE,
36  "InputRCPCollection" ,
37  "Name of the ReconstructedParticle collection which will be cut" ,
38  _InRCPColName ,
39  std::string("Jets") ) ;
40  registerOutputCollection( lcio::LCIO::RECONSTRUCTEDPARTICLE,
41  "OutputRCPCollection" ,
42  "Name of the output ReconstructedParticle collection when WriteNewCollection is true" ,
43  _OutRCPColName ,
44  std::string("CutJets") ) ;
45  registerProcessorParameter( "WriteNewCollection" ,
46  "If true writes the cut jet to a new collection (OutputJetRCPCollection) if false overwrites input" ,
47  _WriteNewCollection,
48  bool(1) ) ;
49  registerProcessorParameter( "SubParticleLists" ,
50  "If true cuts tracks from the particle lists of particles in InputRCPCollection, if false just cuts particles from InputRCPCollection" ,
51  _SubParticleLists,
52  bool(1) ) ;
53 
54  registerOptionalParameter( "a1_Chi2OverDOFEnable" ,
55  "Enable a cut on the value of each tracks chi squared over degrees of freedom" ,
56  _Chi2OverDOFEnable,
57  bool(0) ) ;
58  registerOptionalParameter( "a2_Chi2OverDOFCutLowerThan" ,
59  "If true values lower than the cut value will be cut, if false values higher will be cut" ,
60  _Chi2OverDOFCutLowerThan,
61  bool(0) ) ;
62  registerOptionalParameter( "a3_Chi2OverDOFCutValue" ,
63  "Cut Value" ,
64  _Chi2OverDOFCutValue,
65  float(10.0) ) ;
66 
67  registerOptionalParameter( "b1_D0Enable" ,
68  "Enable a cut on the value of each tracks d0 (no correction for ref point position)" ,
69  _D0Enable,
70  bool(0) ) ;
71  registerOptionalParameter( "b2_D0CutLowerThan" ,
72  "If true values lower than the cut value will be cut, if false values higher will be cut" ,
73  _D0CutLowerThan,
74  bool(0) ) ;
75  registerOptionalParameter( "b3_D0CutValue" ,
76  "Cut Value" ,
77  _D0CutValue,
78  float(20.0) ) ;
79 
80  registerOptionalParameter( "c1_D0ErrEnable" ,
81  "Enable a cut on the value of each tracks d0 Error (sqrt(covariance(d0,d0))" ,
82  _D0ErrEnable,
83  bool(0) ) ;
84  registerOptionalParameter( "c2_D0ErrCutLowerThan" ,
85  "If true values lower than the cut value will be cut, if false values higher will be cut" ,
86  _D0ErrCutLowerThan,
87  bool(0) ) ;
88  registerOptionalParameter( "c3_D0ErrCutValue" ,
89  "Cut Value" ,
90  _D0ErrCutValue,
91  float(250.0/1000.0) ) ;
92 
93  registerOptionalParameter( "d1_Z0Enable" ,
94  "Enable a cut on the value of each tracks Z0 (no correction for ref point position)" ,
95  _Z0Enable,
96  bool(0) ) ;
97  registerOptionalParameter( "d2_Z0CutLowerThan" ,
98  "If true values lower than the cut value will be cut, if false values higher will be cut" ,
99  _Z0CutLowerThan,
100  bool(0) ) ;
101  registerOptionalParameter( "d3_Z0CutValue" ,
102  "Cut Value" ,
103  _Z0CutValue,
104  float(20.0) ) ;
105 
106  registerOptionalParameter( "e1_Z0ErrEnable" ,
107  "Enable a cut on the value of each tracks z0 Error (sqrt(covariance(z0,z0))" ,
108  _Z0ErrEnable,
109  bool(0) ) ;
110  registerOptionalParameter( "e2_Z0ErrCutLowerThan" ,
111  "If true values lower than the cut value will be cut, if false values higher will be cut" ,
112  _Z0ErrCutLowerThan,
113  bool(0) ) ;
114  registerOptionalParameter( "e3_Z0ErrCutValue" ,
115  "Cut Value" ,
116  _Z0ErrCutValue,
117  float(250.0/1000.0) ) ;
118 
119  registerOptionalParameter( "f1_PTEnable" ,
120  "Enable a cut on the value of each tracks PT (radial magnitude of ReconstructedParticle->momentum())" ,
121  _PTEnable,
122  bool(0) ) ;
123  registerOptionalParameter( "f2_PTCutLowerThan" ,
124  "If true values lower than the cut value will be cut, if false values higher will be cut" ,
125  _PTCutLowerThan,
126  bool(1) ) ;
127  registerOptionalParameter( "f3_PTCutValue" ,
128  "Cut Value" ,
129  _PTCutValue,
130  float(0.1) ) ;
131 
132  registerOptionalParameter( "g1_DetectorHitsEnable" ,
133  "Enable a cut on the number seen hits in sub detectors - for more details see documentation" ,
134  _DetectorHitsEnable,
135  bool(0) ) ;
136  std::vector<std::string> SubDetectorNames;
137  SubDetectorNames.push_back("VTX");
138  SubDetectorNames.push_back("FTD");
139  SubDetectorNames.push_back("SIT");
140  SubDetectorNames.push_back("TPC");
141  registerOptionalParameter( "g2_SubDetectorNames" ,
142  "Sub detector names in same order as result of Track->getSubdetectorHitNumbers()" ,
143  _DetectorNames,
144  SubDetectorNames ) ;
145  std::vector<std::string> BoundaryDetectorNames;
146  BoundaryDetectorNames.push_back("TPC");
147  BoundaryDetectorNames.push_back("FTD");
148  registerOptionalParameter( "g3_DetectorHitsBoundaryDetectorNames" ,
149  "Sub detector names of detectors defining the boundary between region 1 and region 2" ,
150  _DetectorHitsBoundaryDetectorNames,
151  BoundaryDetectorNames ) ;
152  std::vector<int> BoundaryCuts;
153  BoundaryCuts.push_back(20);
154  BoundaryCuts.push_back(3);
155  registerOptionalParameter( "g4_DetectorHitsBoundaryCuts" ,
156  "Corresponding to the order of DetectorHitsBoundaryDetectorNames the max number of hits for each detector for region 1, if any of the sub detectors has a higher number of hits then region 2 is used" ,
157  _DetectorHitsBoundaryCuts,
158  BoundaryCuts ) ;
159 
160  std::vector<std::string> Region1DetectorNames;
161  Region1DetectorNames.push_back("VTX");
162  registerOptionalParameter( "g5_DetectorHitsRegion1DetectorNames" ,
163  "Sub detector names of detectors used for cutting in region 1" ,
164  _DetectorHitsRegion1DetectorNames,
165  Region1DetectorNames ) ;
166  std::vector<int> Region1Cuts;
167  Region1Cuts.push_back(3);
168  registerOptionalParameter( "g6_DetectorHitsRegion1Cuts" ,
169  "Corresponding to the order of DetectorHitsRegion1DetectorNames the minimum number of hits for region 1 for a track to pass" ,
170  _DetectorHitsRegion1Cuts,
171  Region1Cuts ) ;
172 
173  std::vector<std::string> Region2DetectorNames;
174  Region2DetectorNames.push_back("VTX");
175  registerOptionalParameter( "g7_DetectorHitsRegion2DetectorNames" ,
176  "Sub detector names of detectors used for cutting in region 2" ,
177  _DetectorHitsRegion2DetectorNames,
178  Region2DetectorNames );
179  std::vector<int> Region2Cuts;
180  Region2Cuts.push_back(0);
181  registerOptionalParameter( "g8_DetectorHitsRegion2Cuts" ,
182  "Corresponding to the order of DetectorHitsRegion2DetectorNames the minimum number of hits for region 2 for a track to pass" ,
183  _DetectorHitsRegion2Cuts,
184  Region2Cuts );
185 
186  // The Processor parameters that will cut on Monte Carlo PDG code of a particles parent. A boolean to enable this cut and an integer
187  // vector of the PDG codes to cut on. Also need the name of the LCRelation collection which refers back to the
188  // Monte Carlo data.
189  registerOptionalParameter( "h1_MCPIDEnable" ,
190  "Enable a cut on the PDG code of the parent of Monte Carlo particle associated with the track. Set the codes to cut on in CutPIDS and the LCRelation collection name to the MC data in MonteCarloLCRelationCollection" ,
191  _MonteCarloPIDEnable,
192  bool(0) ) ;
193  std::vector<int> PDGCodeCuts; //just set up a blank vector for the default input
194  PDGCodeCuts.push_back(0);
195  registerOptionalParameter( "h2_CutPIDS" ,
196  "A list of all the PDG codes of the parent Monte Carlo particle to cut" ,
197  _MonteCarloPIDsToCut,
198  PDGCodeCuts );
199  registerInputCollection( lcio::LCIO::LCRELATION,
200  "h3_MonteCarloLCRelationCollection" ,
201  "Name of the LCRelation collection which links InputRCPCollection to the Monte Carlo data. This can be either the relation collection between the Track objects and the MC, or the ReconstructedParticle objects and the MC. Required only if MCPIDEnable is true" ,
202  _MonteCarloRelationColName ,
203  std::string("Relations") ) ;
204  registerOptionalParameter( "i1_BadParametersEnable" ,
205  "Enable a cut on tracks with NaN parameters or covariances" ,
206  _BadParametersEnable,
207  bool(0) ) ;
208  registerOptionalParameter( "j1_MCVertexCut" ,
209  "Enable a cut on tracks with MC Production Vertices in material" ,
210  _MCVertexEnable,
211  bool(0) ) ;
212 
213 }
214 
215 
216 void RPCutProcessor::init() {
217 
218  // usually a good idea to
219  printParameters() ;
220 
221  _nRun = 0 ;
222  _nEvt = 0 ;
223 
224  if (_MCVertexEnable) {
225  _VxdPar = &(Global::GEAR->getVXDParameters());
226  const gear::GearParameters& BeamPipePar
227  = Global::GEAR->getGearParameters("BeamPipe");
228  _BeamPipeInnerRadius=BeamPipePar.getDoubleVal("BeamPipeRadius");
229  _BeamPipeOuterRadius=_BeamPipeInnerRadius
230  +BeamPipePar.getDoubleVal("BeamPipeThickness");
231  _BeamPipeHalfZ=BeamPipePar.getDoubleVal("BeamPipeHalfZ");
232 
233  // tracks originating from anywhere closer than the following value
234  // to a ladder or to the beam pipe will be removed if the user selects
235  // hadronic interaction suppression. The unit is mm.
236  // I chose a distance of 20 micron based on a histogram of the distribution
237  // of track origin distances from ladder material, which shows a clear
238  // excess between 0 and 20 micron, probably due to the GEANT4 step size.
239  _CutDistance=0.02;
240 
241  // widen beam pipe size parameters accordingly
242  _BeamPipeInnerRadius-=_CutDistance;
243  _BeamPipeOuterRadius+=_CutDistance;
244  _BeamPipeHalfZ+=_CutDistance;
245 
246 #ifdef MCFAIL_DIAGNOSTICS
247  _diaghist_vxmat_xy = new TH2F("vxmathist_xy","xy distribution of track origins supposedly in VXD material",1000,-60,60,1000,-60,60);
248  _diaghist_bpmat_xy = new TH2F("bpmathist_xy","xy distribution of track origins supposedly in beam pipe material",1000,-60,60,1000,-60,60);
249  _diaghist_nomat_xy = new TH2F("nomathist_xy","xy distribution of track origins supposedly outside of material",1000,-60,60,1000,-60,60);
250  _diaghist_vxmat_rz = new TH2F("vxmathist_rz","rz distribution of track origins supposedly in VXD material",1000,-300,300,1000,0,60);
251  _diaghist_bpmat_rz = new TH2F("bpmathist_rz","rz distribution of track origins supposedly in beam pipe material",1000,-300,300,1000,0,60);
252  _diaghist_nomat_rz = new TH2F("nomathist_rz","rz distribution of track origins supposedly outside of material",1000,-300,300,1000,0,60);
253  _diaghist_dist = new TH1F("disthist","distance of track origins from nearest ladder",100,0.0,2);
254  _diaghist_dist_vxmat = new TH1F("disthist_vxmat","distance of track origins from nearest ladder, inside VXD material",100,0.0,2);
255  _diaghist_dist_nomat = new TH1F("disthist_nomat","distance of track origins from nearest ladder, outside VXD material",100,0.0,0.2);
256 #endif
257  }
258 }
259 
260 void RPCutProcessor::processRunHeader( LCRunHeader* run) {
261  _nRun++ ;
262 }
263 
264 void RPCutProcessor::processEvent( LCEvent * evt ) {
265 
266  // this gets called for every event
267  // usually the working horse ...
268 
269  LCCollection* InCol = evt->getCollection( _InRCPColName );
270 
271  //Get the collection of associated Monte Carlo particles if the cut on MC PDG code is enabled.
272  //Use a pointer for the LCRelationNavigator, because I only want it created under certain conditions but it needs to be refered to later.
273  UTIL::LCRelationNavigator* pMCRelationNavigator=0;
274  if(_MonteCarloPIDEnable)
275  {
276  lcio::LCCollection* RelCol=0;
277  try
278  {
279  RelCol = evt->getCollection( _MonteCarloRelationColName );
280  }
281  catch( lcio::Exception exception )
282  {
283  // Don't want to quit, so don't do much here. Just ignore MC PID cuts and apply the other cuts.
284  // Print a warning though.
285  std::cerr << "Unable to get the Monte Carlo data for event " << _nEvt << ". The requested cuts on the true PID will not be applied." << std::endl;
286  RelCol=0;
287  }
288  if(RelCol!=0) // the above "try" was successful
289  {
290  pMCRelationNavigator=new UTIL::LCRelationNavigator(RelCol);
291 
292  // Check to make sure the relation is between the correct types. It has to be "to" MCParticle but can be
293  // "from" either Track or ReconstructedParticle.
294  if( pMCRelationNavigator->getToType()!=lcio::LCIO::MCPARTICLE || (pMCRelationNavigator->getFromType()!=lcio::LCIO::TRACK && pMCRelationNavigator->getFromType()!=lcio::LCIO::RECONSTRUCTEDPARTICLE) )
295  {
296  message<marlin::ERROR4>( "The relation collection provided for the PID cuts is not the correct type (it should be from either Track or ReconstructedParticle, and to MCParticle)." );
297  delete pMCRelationNavigator;
298  pMCRelationNavigator=0;
299  }
300  else vertex_lcfi::MemoryManager<UTIL::LCRelationNavigator>::Event()->registerObject(pMCRelationNavigator);
301  }
302  else // the above "try" was unsuccessful
303  {
304  pMCRelationNavigator=0;
305  }
306  }
307 
308  LCCollection* OutRPCollection;
309  if (_WriteNewCollection)
310  {
311  std::vector<std::string>::const_iterator it = find(evt->getCollectionNames()->begin(),evt->getCollectionNames()->end(),_OutRCPColName);
312  if (it == evt->getCollectionNames()->end())
313  {
314  //Not found do add - TODO do these need memory managment?
315  LCCollection* MyCollection = new LCCollectionVec("ReconstructedParticle");
316  evt->addCollection(MyCollection,_OutRCPColName);
317  }
318 
319  }
320  OutRPCollection = evt->getCollection(_OutRCPColName);
321  if (_WriteNewCollection && !_SubParticleLists)
322  {
323  dynamic_cast<LCCollectionVec*>(OutRPCollection)->setSubset(true);
324  }
325 
326 
327  //TODO Only do this if cut enables and check for data existance
328  //Retrive the list of detector names and create a map so that we know at what index to find them
329  std::map<std::string,int> SubdetectorIndex;
330  if (_DetectorHitsEnable)
331  {
332  if (_DetectorNames.empty())
333  {
334  std::cerr << "Subdetector Names not set" << std::endl;
335  }
336  int i=0;
337  for (StringVec::iterator iDetector = _DetectorNames.begin(); iDetector != _DetectorNames.end();++iDetector)
338  {
339  SubdetectorIndex[*iDetector] = i;
340  ++i;
341  }
342  }
343  //RP Loop
344  int nRCP = InCol->getNumberOfElements() ;
345  for(int i=0; i< nRCP ; i++)
346  {
347  ReconstructedParticle* InputRP = dynamic_cast<ReconstructedParticle*>( InCol->getElementAt( i ) );
348 
349  if (_SubParticleLists)
350  {
351  //We are cutting particles listed in the InputRP
352  ReconstructedParticle* OutputRP;
353  if (_WriteNewCollection)
354  {
355  //Get a copy of the input to cut RPs from
356  OutputRP = new ReconstructedParticleImpl(*dynamic_cast<ReconstructedParticleImpl*>(InputRP));
357  //To do a proper copy we need to copy the PID objects seperatly and remove the read only
358  //TODO Get around nasty cast
359  ((ReconstructedParticleLCFI*)OutputRP)->makeWritable();
360  ((ReconstructedParticleLCFI*)OutputRP)->wipePIDs();
361  ((ReconstructedParticleLCFI*)OutputRP)->copyPIDsFrom(InputRP);
362  }
363  else
364  {
365  OutputRP = InputRP;
366  }
367 
368  ReconstructedParticleVec RPTracks = OutputRP->getParticles();
369  for (ReconstructedParticleVec::iterator iRPTrack = RPTracks.begin();iRPTrack != RPTracks.end();++iRPTrack)
370  {
371  if ((*iRPTrack)->getTracks().size() == 0 || ((*iRPTrack)->getCharge() == 0 && ! _MonteCarloPIDEnable)
372  ||((_D0Enable && _D0Fail(*iRPTrack)) ||
373  (_D0ErrEnable && _D0ErrFail(*iRPTrack)) ||
374  (_Z0Enable && _Z0Fail(*iRPTrack)) ||
375  (_Z0ErrEnable && _Z0ErrFail(*iRPTrack)) ||
376  (_PTEnable && _PTFail(*iRPTrack)) ||
377  (_Chi2OverDOFEnable && _Chi2OverDOFFail(*iRPTrack)) ||
378  (_DetectorHitsEnable && _DetectorHitsFail(*iRPTrack,SubdetectorIndex)) ||
379  (_MonteCarloPIDEnable && _MCPIDFail(*iRPTrack,pMCRelationNavigator)) ||
380  (_BadParametersEnable && _BadParametersFail(*iRPTrack)) ||
381  (_MCVertexEnable &&_MCVertexFail((*iRPTrack),pMCRelationNavigator)) ))
382  {
383  //TODO Remove nasty cast
384  ((ReconstructedParticleLCFI*)OutputRP)->removeParticle(*iRPTrack);
385  }
386  }
387  if (_WriteNewCollection)
388  {
389  OutRPCollection->addElement(OutputRP);
390  }
391  else
392  {
393  //NO OP//
394  }
395  }
396  else //NOT SubParticleLists
397  {
398  //We are deciding weather to cut InputRP itself
399  if ((InputRP)->getTracks().size() == 0 || ((InputRP)->getCharge() == 0 && ! _MonteCarloPIDEnable)
400  ||((_D0Enable && _D0Fail(InputRP)) ||
401  (_D0ErrEnable && _D0ErrFail(InputRP)) ||
402  (_Z0Enable && _Z0Fail(InputRP)) ||
403  (_Z0ErrEnable && _Z0ErrFail(InputRP)) ||
404  (_PTEnable && _PTFail(InputRP)) ||
405  (_Chi2OverDOFEnable && _Chi2OverDOFFail(InputRP)) ||
406  (_DetectorHitsEnable && _DetectorHitsFail(InputRP,SubdetectorIndex)) ||
407  (_MonteCarloPIDEnable && _MCPIDFail(InputRP,pMCRelationNavigator)) ||
408  (_BadParametersEnable && _BadParametersFail(InputRP)) ||
409  (_MCVertexEnable &&_MCVertexFail(InputRP,pMCRelationNavigator)) ))
410  {
411  //Cut InputRP either by removing it from this collection or not adding it to a new one
412  if (_WriteNewCollection)
413  {
414  }
415  else
416  {
417  InCol->removeElementAt(i);
418  }
419  }
420  else
421  {
422  //Keep InputRP either by not removing it from this collection or adding it to a new one
423  if (_WriteNewCollection)
424  {
425  OutRPCollection->addElement(InputRP);
426  }
427  else
428  {
429  }
430  }
431  }
432  }
433  //Clear all objects created
435  _nEvt ++ ;
436 }
437 
438 
439 
440 void RPCutProcessor::check( LCEvent * evt ) {
441  // nothing to check here - could be used to fill checkplots in reconstruction processor
442 }
443 
444 
445 void RPCutProcessor::end(){
446 
447 #ifdef MCFAIL_DIAGNOSTICS
448  if (_MCVertexEnable) {
449  TFile dumpfile((name()+std::string(".root")).c_str(),"RECREATE");
450  _diaghist_vxmat_xy->Write();
451  _diaghist_bpmat_xy->Write();
452  _diaghist_nomat_xy->Write();
453  _diaghist_vxmat_rz->Write();
454  _diaghist_bpmat_rz->Write();
455  _diaghist_nomat_rz->Write();
456  _diaghist_dist->Write();
457  _diaghist_dist_vxmat->Write();
458  _diaghist_dist_nomat->Write();
459  dumpfile.Close();
460  }
461 #endif
462 
463  std::cout << "RPCutProcessor::end() " << name()
464  << " processed " << _nEvt << " events in " << _nRun << " runs "
465  << std::endl ;
466 
467 }
468 
469 bool RPCutProcessor::_Chi2OverDOFFail(ReconstructedParticle* RPTrack)
470 {
471  lcio::Track* Track = RPTrack->getTracks()[0];
472  if (_Chi2OverDOFCutLowerThan)
473  return (Track->getChi2()/(float)Track->getNdf() < _Chi2OverDOFCutValue);
474  else
475  return (Track->getChi2()/(float)Track->getNdf() > _Chi2OverDOFCutValue);
476 }
477 bool RPCutProcessor::_D0Fail(ReconstructedParticle* RPTrack)
478 {
479  lcio::Track* Track = RPTrack->getTracks()[0];
480  if (_D0CutLowerThan)
481  return (fabs(Track->getD0()) < _D0CutValue);
482  else
483  return (fabs(Track->getD0()) > _D0CutValue);
484 }
485 bool RPCutProcessor::_D0ErrFail(ReconstructedParticle* RPTrack)
486 {
487  lcio::Track* Track = RPTrack->getTracks()[0];
488  if (_D0ErrCutLowerThan)
489  return (Track->getCovMatrix()[0] < _D0ErrCutValue);
490  else
491  return (Track->getCovMatrix()[0] > _D0ErrCutValue);
492 }
493 bool RPCutProcessor::_Z0Fail(ReconstructedParticle* RPTrack)
494 {
495  lcio::Track* Track = RPTrack->getTracks()[0];
496  if (_Z0CutLowerThan)
497  return (fabs(Track->getZ0()) < _Z0CutValue);
498  else
499  return (fabs(Track->getZ0()) > _Z0CutValue);
500 }
501 bool RPCutProcessor::_Z0ErrFail(ReconstructedParticle* RPTrack)
502 {
503  lcio::Track* Track = RPTrack->getTracks()[0];
504  if (_Z0ErrCutLowerThan)
505  return (Track->getCovMatrix()[9] < _Z0ErrCutValue);
506  else
507  return (Track->getCovMatrix()[9] > _Z0ErrCutValue);
508 }
509 bool RPCutProcessor::_PTFail(ReconstructedParticle* RPTrack)
510 {
511  const double* p = RPTrack->getMomentum();
512  double pt = sqrt(p[0]*p[0] + p[1]*p[1]);
513  if (_PTCutLowerThan)
514  return (pt < _PTCutValue);
515  else
516  return (pt > _PTCutValue);
517 }
518 bool RPCutProcessor::_DetectorHitsFail(ReconstructedParticle* RPTrack, std::map<std::string,int> SubdetectorIndex)
519 {
520  //First find out if this track is in region 1 or 2
521  bool Region2 = 0;
522  std::vector<int>::const_iterator iCut = _DetectorHitsBoundaryCuts.begin();
523  for (StringVec::iterator iDetector = _DetectorHitsBoundaryDetectorNames.begin(); iDetector != _DetectorHitsBoundaryDetectorNames.end();++iDetector)
524  {
525  //TODO Check for exisance of data
526  if (RPTrack->getTracks()[0]->getSubdetectorHitNumbers()[SubdetectorIndex[*iDetector]] >= (*iCut))
527  {
528  Region2 = 1;
529  break;
530  }
531  ++iCut;
532  }
533  //Cut accordingly 1 = fail
534  if (Region2)
535  {
536  iCut = _DetectorHitsRegion2Cuts.begin();
537  for (StringVec::iterator iDetector = _DetectorHitsRegion2DetectorNames.begin(); iDetector != _DetectorHitsRegion2DetectorNames.end();++iDetector)
538  {
539  if (RPTrack->getTracks()[0]->getSubdetectorHitNumbers()[SubdetectorIndex[*iDetector]] < (*iCut))
540  {
541  return 1;
542  }
543  ++iCut;
544  }
545  }
546  else
547  {
548  iCut = _DetectorHitsRegion1Cuts.begin();
549  for (StringVec::iterator iDetector = _DetectorHitsRegion1DetectorNames.begin(); iDetector != _DetectorHitsRegion1DetectorNames.end();++iDetector)
550  {
551  if (RPTrack->getTracks()[0]->getSubdetectorHitNumbers()[SubdetectorIndex[*iDetector]] < (*iCut))
552  {
553  return 1;
554  }
555  ++iCut;
556  }
557  }
558  return 0;
559 }
560 
561 bool RPCutProcessor::_MCPIDFail( lcio::ReconstructedParticle* RPTrack, UTIL::LCRelationNavigator* pMCRelationNavigator )
562 {
563  if( pMCRelationNavigator==0 )
564  {
565  // Unable to get the MC relation data earlier, so just return false because MC PID cuts are
566  // being ignored. A warning to cerr will have already been printed so no need to do so here.
567  return false;
568  }
569 
570  lcio::Track* Track = RPTrack->getTracks()[0];
571  std::vector<lcio::LCObject*> RelatedMCParticles;
572 
573  // Decide whether to use the Track object or the ReconstructedParticle object to get the MCParticle
574  if( pMCRelationNavigator->getFromType()==lcio::LCIO::TRACK ) RelatedMCParticles = pMCRelationNavigator->getRelatedToObjects(Track);
575  else RelatedMCParticles = pMCRelationNavigator->getRelatedToObjects(RPTrack);
576 
577  //Not sure what to do if it has more than one related MC particle. For now just print a warning and ignore.
578  if( RelatedMCParticles.size()!=1 )
579  {
580  //vjm 03/07/08 comment this out to stop continued warnings
581  //std::cerr << RelatedMCParticles.size() << " MCParticles related to this track! No cuts performed on MC PDG code." << std::endl;
582  return false;
583  }
584  else
585  {
586  std::vector<lcio::MCParticle*> Parents = dynamic_cast<lcio::MCParticle*>(RelatedMCParticles[0])->getParents();
587 
588  if (Parents.empty())
589  return false;
590  else
591  {
592  int truePDGCode=Parents[0]->getPDG();
593  //search for this code in the integer vector of codes to cut on
594  std::vector<int>::iterator iSearchResult=std::find( _MonteCarloPIDsToCut.begin(), _MonteCarloPIDsToCut.end(), truePDGCode );
595  if( iSearchResult==_MonteCarloPIDsToCut.end() )
596  return false; // the PDG code was not found in the vector
597  else
598  {
599  return true; // the PDG code was found
600  }
601  }
602  }
603 }
604 
605 bool RPCutProcessor::_BadParametersFail(lcio::ReconstructedParticle* RPTrack)
606 {
607  lcio::Track* Track = RPTrack->getTracks()[0];
608  bool fail=false;
609  for (short i=0;!fail && i<15;++i)
610  {
611  fail = std::isnan(Track->getCovMatrix()[i]);
612  }
613  return (fail
614  || std::isnan(Track->getD0())
615  || std::isnan(Track->getPhi())
616  || std::isnan(Track->getZ0())
617  || std::isnan(Track->getOmega())
618  || std::isnan(Track->getTanLambda()) );
619 }
620 
621 bool RPCutProcessor::_MCVertexFail(lcio::ReconstructedParticle* RPTrack, UTIL::LCRelationNavigator* pMCRelationNavigator )
622 {
623  if( pMCRelationNavigator==0 )
624  {
625  // Unable to get the MC relation data earlier, so just return false because MC PID cuts are
626  // being ignored. A warning to cerr will have already been printed so no need to do so here.
627  return false;
628  }
629 
630 
631  LCObjectVec objectVec;
632 
633  // Decide whether to use the Track object or the ReconstructedParticle object to get the MCParticle
634  if( pMCRelationNavigator->getFromType()==lcio::LCIO::TRACK ) objectVec = pMCRelationNavigator->getRelatedToObjects(RPTrack->getTracks()[0]);
635  else objectVec = pMCRelationNavigator->getRelatedToObjects(RPTrack);
636 
637  if (objectVec.size() > 0)
638  {
639  MCParticle* mcp = dynamic_cast<MCParticle*> (objectVec[0]);
640  const double * Vert = mcp->getVertex();
641  gear::Vector3D VertVec(Vert[0],Vert[1],Vert[2]);
642 
643  double radius = sqrt((Vert[0] * Vert[0])+(Vert[1] * Vert[1]));
644  if (radius <_BeamPipeInnerRadius) {
645  // this track originates from inside the beam pipe.
646  // we assume that there is no material interaction there.
647 #ifdef MCFAIL_DIAGNOSTICS
648  _diaghist_nomat_xy->Fill(Vert[0],Vert[1]);
649  _diaghist_nomat_rz->Fill(Vert[2],radius);
650 #endif
651  return false;
652  }else if (radius >=_BeamPipeInnerRadius
653  && radius <=_BeamPipeOuterRadius
654  && fabs(Vert[2])<=_BeamPipeHalfZ) {
655  // this track originates from within the beam pipe wall.
656  // thus it is likely a hadronic interaction product
657 #ifdef MCFAIL_DIAGNOSTICS
658  _diaghist_bpmat_xy->Fill(Vert[0],Vert[1]);
659  _diaghist_bpmat_rz->Fill(Vert[2],radius);
660 #endif
661  return true;
662  } else {
663  // next check: are we in VXD ladder material?
664  double dist=_VxdPar->distanceToNearestLadder(VertVec).r();
665 #ifdef MCFAIL_DIAGNOSTICS
666  _diaghist_dist->Fill(dist);
667 #endif
668  if (dist<=_CutDistance) {
669 #ifdef MCFAIL_DIAGNOSTICS
670  _diaghist_dist_vxmat->Fill(dist);
671  _diaghist_vxmat_xy->Fill(Vert[0],Vert[1]);
672  _diaghist_vxmat_rz->Fill(Vert[2],radius);
673 #endif
674  return true;
675  } else {
676 #ifdef MCFAIL_DIAGNOSTICS
677  _diaghist_dist_nomat->Fill(dist);
678  _diaghist_nomat_xy->Fill(Vert[0],Vert[1]);
679  _diaghist_nomat_rz->Fill(Vert[2],radius);
680 #endif
681  return false;
682  }
683  }
684  }
685  else
686  std::cerr << "Warning RPCutProcessor::637 No MCParticle information" << std::endl;
687  return false;
688 }
689 
static MetaMemoryManager * Event()
Returns the Event duration singleton instance of the controller.
void delAllObjects()
Delete all objects of all types held by this instance.
static MemoryManager< T > * Event()
Returns the Event duration singleton instance of the MemoryManager for type T.
Unique Track representation.
Cuts ReconstuctedParticles(RPs) from a collection (or from a list of RPs held by another RP) based on...