1 #include "RPCutProcessor.h"
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"
15 #include <marlin/Global.h>
16 #include <gear/GEAR.h>
17 #include <gearimpl/Vector3D.h>
22 using namespace marlin ;
27 RPCutProcessor::RPCutProcessor() : Processor(
"RPCutProcessor") {
30 _description =
"RPCutProcessor - cuts RPs based on several criteria removing those that have no tracks" ;
35 registerInputCollection( lcio::LCIO::RECONSTRUCTEDPARTICLE,
36 "InputRCPCollection" ,
37 "Name of the ReconstructedParticle collection which will be cut" ,
39 std::string(
"Jets") ) ;
40 registerOutputCollection( lcio::LCIO::RECONSTRUCTEDPARTICLE,
41 "OutputRCPCollection" ,
42 "Name of the output ReconstructedParticle collection when WriteNewCollection is true" ,
44 std::string(
"CutJets") ) ;
45 registerProcessorParameter(
"WriteNewCollection" ,
46 "If true writes the cut jet to a new collection (OutputJetRCPCollection) if false overwrites input" ,
49 registerProcessorParameter(
"SubParticleLists" ,
50 "If true cuts tracks from the particle lists of particles in InputRCPCollection, if false just cuts particles from InputRCPCollection" ,
54 registerOptionalParameter(
"a1_Chi2OverDOFEnable" ,
55 "Enable a cut on the value of each tracks chi squared over degrees of freedom" ,
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,
62 registerOptionalParameter(
"a3_Chi2OverDOFCutValue" ,
67 registerOptionalParameter(
"b1_D0Enable" ,
68 "Enable a cut on the value of each tracks d0 (no correction for ref point position)" ,
71 registerOptionalParameter(
"b2_D0CutLowerThan" ,
72 "If true values lower than the cut value will be cut, if false values higher will be cut" ,
75 registerOptionalParameter(
"b3_D0CutValue" ,
80 registerOptionalParameter(
"c1_D0ErrEnable" ,
81 "Enable a cut on the value of each tracks d0 Error (sqrt(covariance(d0,d0))" ,
84 registerOptionalParameter(
"c2_D0ErrCutLowerThan" ,
85 "If true values lower than the cut value will be cut, if false values higher will be cut" ,
88 registerOptionalParameter(
"c3_D0ErrCutValue" ,
91 float(250.0/1000.0) ) ;
93 registerOptionalParameter(
"d1_Z0Enable" ,
94 "Enable a cut on the value of each tracks Z0 (no correction for ref point position)" ,
97 registerOptionalParameter(
"d2_Z0CutLowerThan" ,
98 "If true values lower than the cut value will be cut, if false values higher will be cut" ,
101 registerOptionalParameter(
"d3_Z0CutValue" ,
106 registerOptionalParameter(
"e1_Z0ErrEnable" ,
107 "Enable a cut on the value of each tracks z0 Error (sqrt(covariance(z0,z0))" ,
110 registerOptionalParameter(
"e2_Z0ErrCutLowerThan" ,
111 "If true values lower than the cut value will be cut, if false values higher will be cut" ,
114 registerOptionalParameter(
"e3_Z0ErrCutValue" ,
117 float(250.0/1000.0) ) ;
119 registerOptionalParameter(
"f1_PTEnable" ,
120 "Enable a cut on the value of each tracks PT (radial magnitude of ReconstructedParticle->momentum())" ,
123 registerOptionalParameter(
"f2_PTCutLowerThan" ,
124 "If true values lower than the cut value will be cut, if false values higher will be cut" ,
127 registerOptionalParameter(
"f3_PTCutValue" ,
132 registerOptionalParameter(
"g1_DetectorHitsEnable" ,
133 "Enable a cut on the number seen hits in sub detectors - for more details see documentation" ,
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()" ,
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,
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,
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,
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,
193 std::vector<int> PDGCodeCuts;
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,
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,
208 registerOptionalParameter(
"j1_MCVertexCut" ,
209 "Enable a cut on tracks with MC Production Vertices in material" ,
216 void RPCutProcessor::init() {
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");
242 _BeamPipeInnerRadius-=_CutDistance;
243 _BeamPipeOuterRadius+=_CutDistance;
244 _BeamPipeHalfZ+=_CutDistance;
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);
260 void RPCutProcessor::processRunHeader( LCRunHeader* run) {
264 void RPCutProcessor::processEvent( LCEvent * evt ) {
269 LCCollection* InCol = evt->getCollection( _InRCPColName );
273 UTIL::LCRelationNavigator* pMCRelationNavigator=0;
274 if(_MonteCarloPIDEnable)
276 lcio::LCCollection* RelCol=0;
279 RelCol = evt->getCollection( _MonteCarloRelationColName );
281 catch( lcio::Exception exception )
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;
290 pMCRelationNavigator=
new UTIL::LCRelationNavigator(RelCol);
294 if( pMCRelationNavigator->getToType()!=lcio::LCIO::MCPARTICLE || (pMCRelationNavigator->getFromType()!=lcio::LCIO::TRACK && pMCRelationNavigator->getFromType()!=lcio::LCIO::RECONSTRUCTEDPARTICLE) )
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;
304 pMCRelationNavigator=0;
308 LCCollection* OutRPCollection;
309 if (_WriteNewCollection)
311 std::vector<std::string>::const_iterator it = find(evt->getCollectionNames()->begin(),evt->getCollectionNames()->end(),_OutRCPColName);
312 if (it == evt->getCollectionNames()->end())
315 LCCollection* MyCollection =
new LCCollectionVec(
"ReconstructedParticle");
316 evt->addCollection(MyCollection,_OutRCPColName);
320 OutRPCollection = evt->getCollection(_OutRCPColName);
321 if (_WriteNewCollection && !_SubParticleLists)
323 dynamic_cast<LCCollectionVec*
>(OutRPCollection)->setSubset(
true);
329 std::map<std::string,int> SubdetectorIndex;
330 if (_DetectorHitsEnable)
332 if (_DetectorNames.empty())
334 std::cerr <<
"Subdetector Names not set" << std::endl;
337 for (StringVec::iterator iDetector = _DetectorNames.begin(); iDetector != _DetectorNames.end();++iDetector)
339 SubdetectorIndex[*iDetector] = i;
344 int nRCP = InCol->getNumberOfElements() ;
345 for(
int i=0; i< nRCP ; i++)
347 ReconstructedParticle* InputRP =
dynamic_cast<ReconstructedParticle*
>( InCol->getElementAt( i ) );
349 if (_SubParticleLists)
352 ReconstructedParticle* OutputRP;
353 if (_WriteNewCollection)
356 OutputRP =
new ReconstructedParticleImpl(*dynamic_cast<ReconstructedParticleImpl*>(InputRP));
368 ReconstructedParticleVec RPTracks = OutputRP->getParticles();
369 for (ReconstructedParticleVec::iterator iRPTrack = RPTracks.begin();iRPTrack != RPTracks.end();++iRPTrack)
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)) ))
387 if (_WriteNewCollection)
389 OutRPCollection->addElement(OutputRP);
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)) ))
412 if (_WriteNewCollection)
417 InCol->removeElementAt(i);
423 if (_WriteNewCollection)
425 OutRPCollection->addElement(InputRP);
440 void RPCutProcessor::check( LCEvent * evt ) {
445 void RPCutProcessor::end(){
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();
463 std::cout <<
"RPCutProcessor::end() " << name()
464 <<
" processed " << _nEvt <<
" events in " << _nRun <<
" runs "
469 bool RPCutProcessor::_Chi2OverDOFFail(ReconstructedParticle* RPTrack)
471 lcio::Track*
Track = RPTrack->getTracks()[0];
472 if (_Chi2OverDOFCutLowerThan)
473 return (Track->getChi2()/(float)Track->getNdf() < _Chi2OverDOFCutValue);
475 return (Track->getChi2()/(float)Track->getNdf() > _Chi2OverDOFCutValue);
477 bool RPCutProcessor::_D0Fail(ReconstructedParticle* RPTrack)
479 lcio::Track* Track = RPTrack->getTracks()[0];
481 return (fabs(Track->getD0()) < _D0CutValue);
483 return (fabs(Track->getD0()) > _D0CutValue);
485 bool RPCutProcessor::_D0ErrFail(ReconstructedParticle* RPTrack)
487 lcio::Track* Track = RPTrack->getTracks()[0];
488 if (_D0ErrCutLowerThan)
489 return (Track->getCovMatrix()[0] < _D0ErrCutValue);
491 return (Track->getCovMatrix()[0] > _D0ErrCutValue);
493 bool RPCutProcessor::_Z0Fail(ReconstructedParticle* RPTrack)
495 lcio::Track* Track = RPTrack->getTracks()[0];
497 return (fabs(Track->getZ0()) < _Z0CutValue);
499 return (fabs(Track->getZ0()) > _Z0CutValue);
501 bool RPCutProcessor::_Z0ErrFail(ReconstructedParticle* RPTrack)
503 lcio::Track* Track = RPTrack->getTracks()[0];
504 if (_Z0ErrCutLowerThan)
505 return (Track->getCovMatrix()[9] < _Z0ErrCutValue);
507 return (Track->getCovMatrix()[9] > _Z0ErrCutValue);
509 bool RPCutProcessor::_PTFail(ReconstructedParticle* RPTrack)
511 const double* p = RPTrack->getMomentum();
512 double pt = sqrt(p[0]*p[0] + p[1]*p[1]);
514 return (pt < _PTCutValue);
516 return (pt > _PTCutValue);
518 bool RPCutProcessor::_DetectorHitsFail(ReconstructedParticle* RPTrack, std::map<std::string,int> SubdetectorIndex)
522 std::vector<int>::const_iterator iCut = _DetectorHitsBoundaryCuts.begin();
523 for (StringVec::iterator iDetector = _DetectorHitsBoundaryDetectorNames.begin(); iDetector != _DetectorHitsBoundaryDetectorNames.end();++iDetector)
526 if (RPTrack->getTracks()[0]->getSubdetectorHitNumbers()[SubdetectorIndex[*iDetector]] >= (*iCut))
536 iCut = _DetectorHitsRegion2Cuts.begin();
537 for (StringVec::iterator iDetector = _DetectorHitsRegion2DetectorNames.begin(); iDetector != _DetectorHitsRegion2DetectorNames.end();++iDetector)
539 if (RPTrack->getTracks()[0]->getSubdetectorHitNumbers()[SubdetectorIndex[*iDetector]] < (*iCut))
548 iCut = _DetectorHitsRegion1Cuts.begin();
549 for (StringVec::iterator iDetector = _DetectorHitsRegion1DetectorNames.begin(); iDetector != _DetectorHitsRegion1DetectorNames.end();++iDetector)
551 if (RPTrack->getTracks()[0]->getSubdetectorHitNumbers()[SubdetectorIndex[*iDetector]] < (*iCut))
561 bool RPCutProcessor::_MCPIDFail( lcio::ReconstructedParticle* RPTrack, UTIL::LCRelationNavigator* pMCRelationNavigator )
563 if( pMCRelationNavigator==0 )
570 lcio::Track* Track = RPTrack->getTracks()[0];
571 std::vector<lcio::LCObject*> RelatedMCParticles;
574 if( pMCRelationNavigator->getFromType()==lcio::LCIO::TRACK ) RelatedMCParticles = pMCRelationNavigator->getRelatedToObjects(Track);
575 else RelatedMCParticles = pMCRelationNavigator->getRelatedToObjects(RPTrack);
578 if( RelatedMCParticles.size()!=1 )
586 std::vector<lcio::MCParticle*> Parents =
dynamic_cast<lcio::MCParticle*
>(RelatedMCParticles[0])->getParents();
592 int truePDGCode=Parents[0]->getPDG();
594 std::vector<int>::iterator iSearchResult=std::find( _MonteCarloPIDsToCut.begin(), _MonteCarloPIDsToCut.end(), truePDGCode );
595 if( iSearchResult==_MonteCarloPIDsToCut.end() )
605 bool RPCutProcessor::_BadParametersFail(lcio::ReconstructedParticle* RPTrack)
607 lcio::Track* Track = RPTrack->getTracks()[0];
609 for (
short i=0;!fail && i<15;++i)
611 fail = std::isnan(Track->getCovMatrix()[i]);
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()) );
621 bool RPCutProcessor::_MCVertexFail(lcio::ReconstructedParticle* RPTrack, UTIL::LCRelationNavigator* pMCRelationNavigator )
623 if( pMCRelationNavigator==0 )
631 LCObjectVec objectVec;
634 if( pMCRelationNavigator->getFromType()==lcio::LCIO::TRACK ) objectVec = pMCRelationNavigator->getRelatedToObjects(RPTrack->getTracks()[0]);
635 else objectVec = pMCRelationNavigator->getRelatedToObjects(RPTrack);
637 if (objectVec.size() > 0)
639 MCParticle* mcp =
dynamic_cast<MCParticle*
> (objectVec[0]);
640 const double * Vert = mcp->getVertex();
641 gear::Vector3D VertVec(Vert[0],Vert[1],Vert[2]);
643 double radius = sqrt((Vert[0] * Vert[0])+(Vert[1] * Vert[1]));
644 if (radius <_BeamPipeInnerRadius) {
647 #ifdef MCFAIL_DIAGNOSTICS
648 _diaghist_nomat_xy->Fill(Vert[0],Vert[1]);
649 _diaghist_nomat_rz->Fill(Vert[2],radius);
652 }
else if (radius >=_BeamPipeInnerRadius
653 && radius <=_BeamPipeOuterRadius
654 && fabs(Vert[2])<=_BeamPipeHalfZ) {
657 #ifdef MCFAIL_DIAGNOSTICS
658 _diaghist_bpmat_xy->Fill(Vert[0],Vert[1]);
659 _diaghist_bpmat_rz->Fill(Vert[2],radius);
664 double dist=_VxdPar->distanceToNearestLadder(VertVec).r();
665 #ifdef MCFAIL_DIAGNOSTICS
666 _diaghist_dist->Fill(dist);
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);
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);
686 std::cerr <<
"Warning RPCutProcessor::637 No MCParticle information" << std::endl;
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...