2 #include "DSTPlotProcessor.h"
17 #include "TMultiGraph.h"
19 #endif // of #ifdef USEROOT
21 #include "EVENT/LCCollection.h"
22 #include "UTIL/PIDHandler.h"
23 #include "IMPL/ReconstructedParticleImpl.h"
24 #include "EVENT/LCParameters.h"
25 #include "EVENT/LCFloatVec.h"
26 #include "EVENT/LCIntVec.h"
28 #include "util/inc/vector3.h"
29 #include "util/inc/util.h"
37 using std::stringstream;
48 DSTPlotProcessor::DSTPlotProcessor() : marlin::Processor(
"DSTPlotProcessor")
50 _description =
"Plots various outputs from the flavour tag" ;
54 registerInputCollection( LCIO::RECONSTRUCTEDPARTICLE,
56 "Name of the collection of ReconstructedParticles that is the jet" ,
58 string(
"FTSelectedJets") ) ;
61 registerProcessorParameter(
"OutputFilename" ,
62 "Filename for the output" ,
64 string(
"DSTPlotProcessorOutput") ) ;
66 registerProcessorParameter(
"CheckDSTParameters" ,
67 "Checks the DST parameters against the full parameters - gives a lot of print out" ,
74 registerOptionalParameter(
"FlavourTagCollection" ,
75 "Name of the LCFloatVec Collection containing the flavour tags, only need this if you want to check DST parameters" ,
76 _FlavourTagCollectionName,
77 std::string(
"FlavourTag") ) ;
79 registerOptionalParameter(
"FlavourTagInputsCollection" ,
80 "Name of the LCFloatVec Collection that contains the flavour tag inputs, only need this if you want to check DST parameters" ,
81 _FlavourTagInputsCollectionName,
82 std::string(
"FlavourTagInputs") ) ;
86 registerOptionalParameter(
"TrueJetFlavourCollection" ,
87 "Name of the LCIntVec collection containing the true flavour of the jets, only need this if you want to check DST parameters" ,
88 _TrueJetFlavourColName ,
89 std::string(
"TrueJetFlavour") ) ;
94 DSTPlotProcessor::~DSTPlotProcessor()
98 void DSTPlotProcessor::init()
101 cout << _description << endl
102 <<
"-------------------------------------------------" << endl
110 void DSTPlotProcessor::processRunHeader( LCRunHeader* pRun )
117 std::vector<std::string> VarNames;
118 (pRun->parameters()).getStringVals(_FlavourTagCollectionName,VarNames);
119 for (
size_t i = 0;i < VarNames.size();++i)
120 _IndexOfForEachTag[VarNames[i]] = i;
124 std::vector<std::string> trueJetFlavourVarNames;
126 for (
size_t i = 0;i < trueJetFlavourVarNames.size();++i)
127 _FlavourIndex[trueJetFlavourVarNames[i]] = i;
131 std::vector<std::string> InputVarNames;
132 (pRun->parameters()).getStringVals(_FlavourTagInputsCollectionName,InputVarNames);
133 for (
size_t i = 0;i < InputVarNames.size();++i)
134 _InputsIndex[InputVarNames[i]] = i;
140 void DSTPlotProcessor::processEvent( LCEvent* pEvent )
144 _checkDSTParameters( pEvent);
151 if( pJetCollection->getTypeName()!=LCIO::RECONSTRUCTEDPARTICLE )
153 stringstream message;
155 <<
"########################################################################################\n"
156 <<
"# DSTPlotProcessor - #\n"
157 <<
"# The jet collection requested (\"" <<
_JetCollectionName <<
"\") is not of the type \"" << LCIO::RECONSTRUCTEDPARTICLE <<
"\" #\n"
158 <<
"########################################################################################" << endl;
159 throw EventException( message.str() );
165 ReconstructedParticle* pJet;
167 for(
int a=0; a<pJetCollection->getNumberOfElements(); ++a )
171 pJet=
dynamic_cast<ReconstructedParticle*
>( pJetCollection->getElementAt(a) );
184 void DSTPlotProcessor::end()
188 cout <<
"The largest jet energy was " <<
_jetEMax << endl;
218 double CThJ_upper=0.95;
222 const double* mom=pJet->getMomentum();
225 jetMomentum.makeUnit();
226 double cosTheta=jetMomentum.dot( zAxis );
227 if( fabs(cosTheta)<=CThJ_lower || fabs(cosTheta)>=CThJ_upper )
return false;
241 ReconstructedParticle* pJet=
dynamic_cast<ReconstructedParticle*
>( pJetCollection->getElementAt(jet) );
244 PIDHandler pidh( pJetCollection ) ;
247 int alid = pidh.getAlgorithmID(
"LCFIFlavourTag");
250 const ParticleID& pid = pidh.getParticleID(pJet,alid);
253 FloatVec params = pid.getParameters();
256 double bTag= params[pidh.getParameterIndex(alid,
"BTag")];
257 double cTag= params[pidh.getParameterIndex(alid,
"CTag")];
258 double cTagBBack= params[pidh.getParameterIndex(alid,
"BCTag")];
261 int mcid = pidh.getAlgorithmID(
"MCTruth");
264 const ParticleID& mcpid = pidh.getParticleID(pJet,mcid);
267 FloatVec mcparams = mcpid.getParameters();
270 float jetType = mcparams[pidh.getParameterIndex(mcid,
"TrueJetFlavour")];
277 if( bTag > -10 ) _BTagEfficiencyPurity.add_signal( bTag );
278 if( cTag > -10 ) _CTagEfficiencyPurity.add_background( cTag );
279 if( cTagBBack > -10 ) _BCTagEfficiencyPurity.add_background( cTagBBack );
281 else if( jetType==
C_JET )
283 if( bTag > -10 ) _BTagEfficiencyPurity.add_background( bTag );
284 if( cTag > -10 ) _CTagEfficiencyPurity.add_signal( cTag );
285 if( cTagBBack > -10 ) _BCTagEfficiencyPurity.add_signal( cTagBBack );
289 if( bTag > -10 ) _BTagEfficiencyPurity.add_background( bTag );
290 if( cTag > -10 ) _CTagEfficiencyPurity.add_background( cTag );
299 stringstream filenameStream;
300 filenameStream << filename <<
".root";
302 TFile rootFile( filenameStream.str().c_str(),
"RECREATE",
"Various plots from the plot processor" );
308 int numberOfPoints=100;
311 vector< pair<double,double> > BTagefficiencyPurityPairs=_BTagEfficiencyPurity.eff_pur( numberOfPoints );
312 vector< pair<double,double> > CTagefficiencyPurityPairs=_CTagEfficiencyPurity.eff_pur( numberOfPoints );
313 vector< pair<double,double> > BCTagefficiencyPurityPairs=_BCTagEfficiencyPurity.eff_pur( numberOfPoints );
317 double BNetEfficiency[100];
318 double BNetPurity[100];
319 double CNetEfficiency[100];
320 double CNetPurity[100];
321 double BCNetEfficiency[100];
322 double BCNetPurity[100];
324 for(
int a=0; a<numberOfPoints; ++a )
326 BNetEfficiency[a]=BTagefficiencyPurityPairs[a].first;
327 BNetPurity[a]=BTagefficiencyPurityPairs[a].second;
328 CNetEfficiency[a]=CTagefficiencyPurityPairs[a].first;
329 CNetPurity[a]=CTagefficiencyPurityPairs[a].second;
330 BCNetEfficiency[a]=BCTagefficiencyPurityPairs[a].first;
331 BCNetPurity[a]=BCTagefficiencyPurityPairs[a].second;
334 TGraph BNetGraph( 99, BNetEfficiency, BNetPurity );
335 BNetGraph.SetName( (
" B Tag") );
336 TGraph CNetGraph( 99, CNetEfficiency, CNetPurity );
337 CNetGraph.SetName( (
" C Tag") );
338 TGraph BCNetGraph( 99, BCNetEfficiency, BCNetPurity );
339 BCNetGraph.SetName( (
" BC Tag") );
341 TMultiGraph TagGraphs;
342 TagGraphs.SetName(
"TagGraphs" );
343 TagGraphs.Add( &BNetGraph );
344 TagGraphs.Add( &CNetGraph );
345 TagGraphs.Add( &BCNetGraph );
348 TagGraphs.SetTitle( (
string(
"Efficiency-purity for ") ).c_str() );
351 pAxis=TagGraphs.GetXaxis();
352 if(pAxis) pAxis->SetTitle(
"Efficiency" );
354 pAxis=TagGraphs.GetYaxis();
355 if(pAxis) pAxis->SetTitle(
"Purity" );
357 rootFile.Add( &TagGraphs );
363 TH1F jetEnergyHistogram(
"jetEnergyHistogram",
"Jet energies", 200, 0, 120 );
364 const vector<double> jetEnergyData=
_jetEnergy.sorted_data();
365 for( vector<double>::const_iterator i=jetEnergyData.begin(); i<jetEnergyData.end(); ++i ) jetEnergyHistogram.Fill( (*i) );
375 const vector<string>* pCollectionNames=pEvent->getCollectionNames();
377 cout <<
"The available collections are: (name - type)" << endl;
378 for( vector<string>::const_iterator i=pCollectionNames->begin(); i<pCollectionNames->end(); ++i )
380 LCCollection* pCollection=pEvent->getCollection( (*i) );
381 const string typeName=pCollection->getTypeName();
382 cout <<
" " << (*i) <<
" - " << typeName << endl;
388 void DSTPlotProcessor::_checkDSTParameters( LCEvent* pEvent)
397 LCCollection* flavourTagInputsCollection = pEvent->getCollection( _FlavourTagInputsCollectionName);
400 LCCollection* TagCollection=pEvent->getCollection( _FlavourTagCollectionName );
404 for(
int jet=0; jet<JetCollection->getNumberOfElements(); jet++ )
409 ReconstructedParticle* pJet =
dynamic_cast<ReconstructedParticle*
>( JetCollection->getElementAt(jet) );
412 LCFloatVec* flavourTags =
dynamic_cast<LCFloatVec*
>(TagCollection->getElementAt( jet ));
413 float fullBTag = (*flavourTags)[_IndexOfForEachTag[
"BTag"]];
414 float fullCTag = (*flavourTags)[_IndexOfForEachTag[
"CTag"]];
415 float fullBCTag =(*flavourTags)[_IndexOfForEachTag[
"BCTag"]];
417 LCFloatVec* pJetFlavour =
dynamic_cast<LCFloatVec*
>(trueJetFlavourCollection->getElementAt( jet ));
418 float fullTruePDGCode = (*pJetFlavour)[_FlavourIndex[
"TruePDGCode"]];
419 float fullTruePartonCharge =(*pJetFlavour)[_FlavourIndex[
"TruePartonCharge"]];
420 float fullTrueHadronCharge= (*pJetFlavour)[_FlavourIndex[
"TrueHadronCharge"]];
421 float fullTrueJetFlavour = (*pJetFlavour)[_FlavourIndex[
"TrueJetFlavour"]];
423 LCFloatVec* tagInputs =
dynamic_cast<LCFloatVec*
>(flavourTagInputsCollection->getElementAt( jet ));
424 int fullNumVertices = int((*tagInputs)[_InputsIndex[
"NumVertices"]]);
426 float fullJointProbRPhi=999;
427 float fullJointProbZ=999;
428 float fullNumTracksInVertices=999;
429 float fullDecayLength=999;
430 float fullDecayLengthSignificance=999;
431 float fullRawMomentum=999;
432 float fullPTCorrectedMass=999;
433 float fullSecondaryVertexProbability=999;
435 if(fullNumVertices > 1)
437 fullJointProbRPhi=(*tagInputs)[_InputsIndex[
"JointProbRPhi"]];
438 fullJointProbZ=(*tagInputs)[_InputsIndex[
"JointProbZ"]];
439 fullNumTracksInVertices=(*tagInputs)[_InputsIndex[
"NumTracksInVertices"]];
440 fullDecayLength=(*tagInputs)[_InputsIndex[
"DecayLength"]];
441 fullDecayLengthSignificance=(*tagInputs)[_InputsIndex[
"DecayLengthSignificance"]];
442 fullRawMomentum=(*tagInputs)[_InputsIndex[
"RawMomentum"]];
443 fullPTCorrectedMass=(*tagInputs)[_InputsIndex[
"PTCorrectedMass"]];
444 fullSecondaryVertexProbability=(*tagInputs)[_InputsIndex[
"SecondaryVertexProbability"]];
449 PIDHandler pidh( JetCollection ) ;
450 int alid = pidh.getAlgorithmID(
"LCFIFlavourTag");
451 const ParticleID& pid = pidh.getParticleID(pJet,alid);
452 FloatVec params = pid.getParameters();
454 float dstBTag = params[pidh.getParameterIndex(alid,
"BTag")];
455 float dstCTag = params[pidh.getParameterIndex(alid,
"CTag")];
456 float dstBCTag =params[pidh.getParameterIndex(alid,
"BCTag")];
458 int mcid = pidh.getAlgorithmID(
"MCTruth");
459 const ParticleID& mcpid = pidh.getParticleID(pJet,mcid);
460 FloatVec mcparams = mcpid.getParameters();
462 float dstTruePDGCode = mcparams[pidh.getParameterIndex(mcid,
"TruePDGCode")];
463 float dstTruePartonCharge =mcparams[pidh.getParameterIndex(mcid,
"TruePartonCharge")];
464 float dstTrueHadronChange= mcparams[pidh.getParameterIndex(mcid,
"TrueHadronCharge")];
465 float dstTrueJetFlavour = mcparams[pidh.getParameterIndex(mcid,
"TrueJetFlavour")];
467 int dstNumVertices = int(params[pidh.getParameterIndex(alid,
"NumVertices")]);
469 float dstJointProbRPhi=999;
470 float dstJointProbZ=999;
471 float dstNumTracksInVertices=999;
472 float dstDecayLength=999;
473 float dstDecayLengthSignificance=999;
474 float dstRawMomentum=999;
475 float dstPTCorrectedMass=999;
476 float dstSecondaryVertexProbability=999;
478 if(dstNumVertices > 1)
480 dstJointProbRPhi=params[pidh.getParameterIndex(alid,
"JointProbRPhi")];
481 dstJointProbZ=params[pidh.getParameterIndex(alid,
"JointProbZ")];
482 dstNumTracksInVertices=params[pidh.getParameterIndex(alid,
"NumTracksInVertices")];
483 dstDecayLength=params[pidh.getParameterIndex(alid,
"DecayLength")];
484 dstDecayLengthSignificance=params[pidh.getParameterIndex(alid,
"DecayLengthSignificance")];
485 dstRawMomentum=params[pidh.getParameterIndex(alid,
"RawMomentum")];
486 dstPTCorrectedMass=params[pidh.getParameterIndex(alid,
"PTCorrectedMass")];
487 dstSecondaryVertexProbability=params[pidh.getParameterIndex(alid,
"SecondaryVertexProbability")];
492 cout<< dstBTag <<
" "<<fullBTag<<endl;;
493 cout<< dstCTag <<
" "<<fullCTag<<endl;;
494 cout<< dstBCTag <<
" "<<fullBCTag<<endl;;
497 cout<< dstTruePDGCode <<
" "<<fullTruePDGCode<<endl;;
498 cout<< dstTruePartonCharge <<
" "<<fullTruePartonCharge<<endl;;
499 cout<< dstTrueHadronChange<<
" "<<fullTrueHadronCharge<<endl;;
500 cout<< dstTrueJetFlavour <<
" "<<fullTrueJetFlavour<<endl;;
501 cout<< dstNumVertices <<
" "<<fullNumVertices<<endl;;
503 if(dstNumVertices > 1)
505 cout<< dstJointProbRPhi <<
" "<<fullJointProbRPhi<<endl;;
506 cout<< dstJointProbZ <<
" "<<fullJointProbZ<<endl;;
507 cout<< dstNumTracksInVertices <<
" "<<fullNumTracksInVertices<<endl;;
508 cout<< dstDecayLength <<
" "<<fullDecayLength<<endl;;
509 cout<< dstDecayLengthSignificance <<
" "<<fullDecayLengthSignificance<<endl;;
510 cout<< dstRawMomentum <<
" "<<fullRawMomentum<<endl;;
511 cout<< dstPTCorrectedMass <<
" "<<fullPTCorrectedMass<<endl;;
512 cout<< dstSecondaryVertexProbability <<
" "<<fullSecondaryVertexProbability<<endl;;
517 if(fabs( dstBTag -fullBTag) > ep)
518 cout<<
"B Parameters not equal"<<endl;
519 if(fabs( dstCTag -fullCTag) > ep)
520 cout<<
"C Parameters not equal"<<endl;
521 if(fabs( dstBCTag -fullBCTag) > ep)
522 cout<<
"BC Parameters not equal"<<endl;
525 if(fabs( dstTruePDGCode -fullTruePDGCode) > ep)
526 cout<<
"TPDG Parameters not equal"<<endl;
527 if(fabs( dstTruePartonCharge -fullTruePartonCharge) > ep)
528 cout<<
"TP Parameters not equal"<<endl;
529 if(fabs( dstTrueHadronChange-fullTrueHadronCharge) > ep)
530 cout<<
"TH Parameters not equal"<<endl;
531 if(fabs( dstTrueJetFlavour -fullTrueJetFlavour) > ep)
532 cout<<
"TJ Parameters not equal"<<endl;
533 if(fabs( dstNumVertices -fullNumVertices) > ep)
534 cout<<
"NV Parameters not equal"<<endl;
536 if(dstNumVertices > 1)
538 if(fabs( dstJointProbRPhi -fullJointProbRPhi) > ep)
539 cout<<
"JPRP Parameters not equal"<<endl;
540 if(fabs( dstJointProbZ -fullJointProbZ) > ep)
541 cout<<
"JPZ Parameters not equal"<<endl;
542 if(fabs( dstNumTracksInVertices -fullNumTracksInVertices) > ep)
543 cout<<
"NT Parameters not equal"<<endl;
544 if(fabs( dstDecayLength -fullDecayLength) > ep)
545 cout<<
"DL Parameters not equal"<<endl;
546 if(fabs( dstDecayLengthSignificance -fullDecayLengthSignificance) > ep)
547 cout<<
"DLS Parameters not equal"<< dstDecayLengthSignificance -fullDecayLengthSignificance<<endl;
548 if(fabs( dstRawMomentum -fullRawMomentum) > ep)
549 cout<<
"raw p Parameters not equal"<<endl;
550 if(fabs( dstPTCorrectedMass -fullPTCorrectedMass) > ep)
551 cout<<
"pt cor Parameters not equal"<<endl;
552 if(fabs( dstSecondaryVertexProbability -fullSecondaryVertexProbability) > ep)
553 cout<<
"sv Parameters not equal "<<endl;
562 catch(DataNotAvailableException &e){
566 catch(DataNotAvailableException &e){
570 catch(DataNotAvailableException &e){
571 std::cout <<
"Collection " << _FlavourTagInputsCollectionName <<
" is unavailable " << std::endl;
574 catch(DataNotAvailableException &e){
575 std::cout <<
"Collection " <<_FlavourTagCollectionName <<
" is unavailable " << std::endl;
void _fillPlots(LCEvent *pEvent, unsigned int jet)
bool _passesJetCuts(lcio::ReconstructedParticle *pJet)
std::string _JetCollectionName
std::string _TrueJetFlavourColName
std::string _OutputFilename
bool _passesEventCuts(lcio::LCEvent *pEvent)
histogram_data< double > _jetEnergy
void _outputDataToFile(std::string filename)
void _displayCollectionNames(lcio::LCEvent *pEvent)
Creates some sample plots from the data calculated by the LCFI vertex package.