2 #ifndef MARLIN_USE_AIDA
4 #warning "-------------------------------------------------------------------------------"
5 #warning "- DSTAIDAPlotProcessor requires MARLIN_USE_AIDA to be defined. Did you enable -"
6 #warning "- AIDA when compiling Marlin? DSTAIDAPlotProcessor will not be compiled. -"
7 #warning "-------------------------------------------------------------------------------"
12 #include "DSTAIDAPlotProcessor.h"
23 #include <marlin/AIDAProcessor.h>
24 #include <AIDA/IMeasurement.h>
25 #include <AIDA/IHistogramFactory.h>
26 #include <AIDA/IHistogram1D.h>
27 #include <AIDA/IHistogram2D.h>
28 #include <AIDA/IAxis.h>
29 #include <AIDA/ITree.h>
30 #include <AIDA/ITupleFactory.h>
31 #include <AIDA/ICloud2D.h>
32 #include <UTIL/LCRelationNavigator.h>
39 #warning "USING_RAIDA defined"
42 #warning "USING_JAIDA defined"
46 #ifdef USING_JAIDA//Data point sets aren't implemented in RAIDA - which is a shame as they have functionality not given by histograms
48 #include <AIDA/IDataPointSet.h>
49 #include <AIDA/IDataPointSetFactory.h>
50 #include <AIDA/IDataPoint.h>
55 #include "EVENT/LCCollection.h"
56 #include "UTIL/PIDHandler.h"
57 #include "IMPL/ReconstructedParticleImpl.h"
58 #include "EVENT/LCParameters.h"
59 #include "EVENT/LCFloatVec.h"
60 #include "EVENT/LCIntVec.h"
62 #include "util/inc/vector3.h"
63 #include "util/inc/util.h"
71 using std::stringstream;
79 #include "TypesafeCollection.h"
84 DSTAIDAPlotProcessor::DSTAIDAPlotProcessor() : marlin::Processor(
"DSTAIDAPlotProcessor")
86 _description =
"Plots various outputs from the flavour tag" ;
90 registerInputCollection( LCIO::RECONSTRUCTEDPARTICLE,
92 "Name of the collection of ReconstructedParticles that is the jet" ,
94 string(
"FTSelectedJets") ) ;
99 DSTAIDAPlotProcessor::~DSTAIDAPlotProcessor()
103 void DSTAIDAPlotProcessor::init()
106 cout << _description << endl
107 <<
"-------------------------------------------------" << endl
113 _VertexCatNames.resize(N_VERTEX_CATEGORIES+1);
114 _VertexCatNames[0]=
"AnyNumberOfVertices";
115 _VertexCatNames[1]=
"OneVertex";
116 _VertexCatNames[2]=
"TwoVertices";
117 _VertexCatNames[3]=
"ThreeOrMoreVertices";
119 _NumVertexCatDir.resize(N_VERTEX_CATEGORIES+1);
120 _NumVertexCatDir[1]=
"OneVertex";
121 _NumVertexCatDir[2]=
"TwoVertices";
122 _NumVertexCatDir[3]=
"ThreeOrMoreVertices";
123 _NumVertexCatDir[0]=
"AnyNumberOfVertices";
127 AIDA::IHistogramFactory* pHistogramFactory=marlin::AIDAProcessor::histogramFactory(
this );
128 AIDA::ITree* pTree=marlin::AIDAProcessor::tree(
this );
130 if( pHistogramFactory!=0 )
132 if (!(pTree->cd(
"/" + name() +
"/"))) {
133 pTree->mkdir(
"/" + name() +
"/" );
134 pTree->cd(
"/" + name() +
"/");
140 CreateFlavourTagTuple();
144 std::cerr <<
"### " << __FILE__ <<
"(" << __LINE__ <<
"): Unable to get the histogram factory! No histograms will be made."<< std::endl;
151 void DSTAIDAPlotProcessor::processRunHeader( LCRunHeader* pRun )
155 CreateFlavourTagInputPlots(pRun);
159 void DSTAIDAPlotProcessor::processEvent( LCEvent* pEvent )
163 LCCollection* jetCollection;
165 jetCollection=pEvent->getCollection( _JetCollectionName );
167 catch(DataNotAvailableException &e){
168 std::cout <<
"Collection " <<_JetCollectionName <<
" is unavailable, can't plot anything " << std::endl;
172 PIDHandler pidh( jetCollection ) ;
175 pidh.getAlgorithmID(
"MCTruth");
177 catch(UnknownAlgorithm &e)
179 std::cout <<
"MCTruth Algorithm is unavailable, can't plot anything " << std::endl;
184 pidh.getAlgorithmID(
"LCFIFlavourTag");
186 catch(UnknownAlgorithm &e){
187 std::cout <<
"LCFIFlavourTag Algorithm is unavailable, can't plot anything " << std::endl;
193 if( _passesEventCuts(pEvent) )
196 ReconstructedParticle* pJet;
198 for(
int jetNumber=0; jetNumber<jetCollection->getNumberOfElements(); ++jetNumber )
200 pJet=
dynamic_cast<ReconstructedParticle*
>( jetCollection->getElementAt(jetNumber) );
203 if( _passesJetCuts(pJet) )
205 FillTagPlots( pEvent, jetNumber );
206 FillInputsPlots( pEvent, jetNumber );
218 void DSTAIDAPlotProcessor::end()
221 CalculateIntegralAndBackgroundPlots();
222 CalculateEfficiencyPurityPlots();
250 double CThJ_upper=0.95;
254 const double* mom=pJet->getMomentum();
257 jetMomentum.makeUnit();
258 double cosTheta=jetMomentum.dot( zAxis );
259 if( fabs(cosTheta)<=CThJ_lower || fabs(cosTheta)>=CThJ_upper )
return false;
269 const vector<string>* pCollectionNames=pEvent->getCollectionNames();
271 cout <<
"The available collections are: (name - type)" << endl;
272 for( vector<string>::const_iterator i=pCollectionNames->begin(); i<pCollectionNames->end(); ++i )
274 LCCollection* pCollection=pEvent->getCollection( (*i) );
275 const string typeName=pCollection->getTypeName();
276 cout <<
" " << (*i) <<
" - " << typeName << endl;
282 void DSTAIDAPlotProcessor::CreateTagPlots()
284 AIDA::IHistogramFactory* pHistogramFactory=marlin::AIDAProcessor::histogramFactory(
this );
285 AIDA::ITree* pTree=marlin::AIDAProcessor::tree(
this );
289 for (
unsigned int iVertexCat=0; iVertexCat < N_VERTEX_CATEGORIES+1; ++iVertexCat ){
290 _pLightJetBTag[_VertexCatNames[iVertexCat]]=0;
291 _pLightJetCTag[_VertexCatNames[iVertexCat]]=0;
292 _pLightJetBCTag[_VertexCatNames[iVertexCat]]=0;
294 _pBJetBTag[_VertexCatNames[iVertexCat]]=0;
295 _pBJetCTag[_VertexCatNames[iVertexCat]]=0;
296 _pBJetBCTag[_VertexCatNames[iVertexCat]]=0;
298 _pCJetBTag[_VertexCatNames[iVertexCat]]=0;
299 _pCJetCTag[_VertexCatNames[iVertexCat]]=0;
300 _pCJetBCTag[_VertexCatNames[iVertexCat]]=0;
302 _pBTagBackgroundValues[_VertexCatNames[iVertexCat]]=0;
303 _pCTagBackgroundValues[_VertexCatNames[iVertexCat]]=0;
304 _pBCTagBackgroundValues[_VertexCatNames[iVertexCat]]=0;
311 for (
unsigned int iVertexCat=0; iVertexCat <= N_VERTEX_CATEGORIES; ++iVertexCat ){
313 std::string nvname = _VertexCatNames[iVertexCat];
316 pTree->cd(
"/" + name() +
"/" +
"/");
318 if (!pTree->cd( _NumVertexCatDir[iVertexCat]+
"/")) {
319 pTree->mkdir( _NumVertexCatDir[iVertexCat]+
"/");
320 pTree->cd( _NumVertexCatDir[iVertexCat]+
"/");
323 _pLightJetBTag[nvname] = pHistogramFactory->createHistogram1D(
"Numbers of light jets by B-tag NN value. ("+ nvname +
")",_numberOfPoints , 0, 1.0 );
324 _pLightJetCTag[nvname] = pHistogramFactory->createHistogram1D(
"Numbers of light jets by C-tag NN value. ("+ nvname +
")", _numberOfPoints, 0, 1.0 );
325 _pLightJetBCTag[nvname] = pHistogramFactory->createHistogram1D(
"Numbers of light jets by BC-tag NN value. ("+ nvname +
")", _numberOfPoints, 0, 1.0 );
326 _pBJetBTag[nvname] = pHistogramFactory->createHistogram1D(
"Numbers of B jets by B-tag NN value. ("+ nvname +
")", _numberOfPoints, 0, 1.0 );
327 _pBJetCTag[nvname] = pHistogramFactory->createHistogram1D(
"Numbers of B jets by C-tag NN value. ("+ nvname +
")", _numberOfPoints, 0, 1.0 );
328 _pBJetBCTag[nvname] = pHistogramFactory->createHistogram1D(
"Numbers of B jets by BC-tag NN value. ("+ nvname +
")", _numberOfPoints, 0, 1.0 );
329 _pCJetBTag[nvname] = pHistogramFactory->createHistogram1D(
"Numbers of C jets by B-tag NN value. ("+ nvname +
")", _numberOfPoints, 0, 1.0 );
330 _pCJetCTag[nvname] = pHistogramFactory->createHistogram1D(
"Numbers of C jets by C-tag NN value. ("+ nvname +
")", _numberOfPoints, 0, 1.0 );
331 _pCJetBCTag[nvname] = pHistogramFactory->createHistogram1D(
"Numbers of C jets by BC-tag NN value. ("+ nvname +
")", _numberOfPoints, 0, 1.0 );
337 void DSTAIDAPlotProcessor::CreateFlavourTagTuple()
341 AIDA::ITree* pTree=marlin::AIDAProcessor::tree(
this );
346 AIDA::ITupleFactory* pTupleFactory=marlin::AIDAProcessor::tupleFactory(
this );
350 pTree->cd(
"/" + name());
354 std::string columnNames=
"int TrueJetFlavour=-1, int NumberOfVertices=-1, int NumberOfTracksInVertices=-1, float DecayLength = -999., float DecayLengthSignificance= -999., float JointProbRPhi= -999., float JointProbZ= -999., float PTCorrectedMass= -999., float RawMomentum= -999., float SecondaryVertexProbability= -999. ";
357 if (!pTree->cd(
"/" + name() +
"/" )) {
358 pTree->cd(
"/" + name() +
"/");
359 if (!pTree->cd(
"/")) {
365 if (!pTree->cd(
"TupleDir/")) {
366 pTree->mkdir(
"TupleDir/");
367 pTree->cd(
"TupleDir/");
370 _pMyTuple=pTupleFactory->create(
"FlavourTagInputsTuple",
"FlavourTagInputsTuple", columnNames);
381 void DSTAIDAPlotProcessor::CreateFlavourTagInputPlots(LCRunHeader* pRun )
385 AIDA::IHistogramFactory* pHistogramFactory=marlin::AIDAProcessor::histogramFactory(
this );
386 AIDA::ITree* pTree=marlin::AIDAProcessor::tree(
this );
387 std::vector<std::string> VarNames;
389 VarNames.push_back(
"NumVertices");
390 VarNames.push_back(
"JointProbRPhi");
391 VarNames.push_back(
"JointProbZ");
392 VarNames.push_back(
"NumTracksInVertices");
393 VarNames.push_back(
"DecayLength");
394 VarNames.push_back(
"DecayLengthSignificance");
395 VarNames.push_back(
"RawMomentum");
396 VarNames.push_back(
"PTCorrectedMass");
397 VarNames.push_back(
"SecondaryVertexProbability");
399 for (
size_t i = 0;i < VarNames.size();++i)
403 if( _inputsHistogramsBJets[VarNames[i]]==0 )
405 pTree->cd(
"/" + name() +
"/");
406 if( !pTree->cd(
"NNInputs" ) )
408 pTree->mkdir(
"NNInputs" ) ;
409 pTree->cd(
"NNInputs" ) ;
412 int numberOfPoints=_numberOfPoints/4;
417 if( VarNames[i]==
"NumVertices" )
423 else if( VarNames[i]==
"NumTracksInVertices" )
432 else if (VarNames[i]==
"DecayLengthSignificance")
438 else if (VarNames[i]==
"DecayLength")
444 else if (VarNames[i]==
"JointProbRPhi" || VarNames[i]==
"JointProbZ"|| VarNames[i]==
"SecondaryVertexProbability")
450 else if (VarNames[i]==
"RawMomentum" )
456 else if (VarNames[i]==
"PTCorrectedMass" )
463 pTree->cd(
"/" + name() +
"/" +
"NNInputs");
464 if( !pTree->cd(
"bJets" ) )
466 pTree->mkdir(
"bJets" ) ;
467 pTree->cd(
"bJets" ) ;
470 _inputsHistogramsBJets[VarNames[i]]=pHistogramFactory->createHistogram1D( VarNames[i], numberOfPoints, lowerBin, higerBin );
472 pTree->cd(
"/" + name() +
"/" +
"NNInputs");
473 if( !pTree->cd(
"cJets" ) )
475 pTree->mkdir(
"cJets" ) ;
476 pTree->cd(
"cJets" ) ;
479 _inputsHistogramsCJets[VarNames[i]]=pHistogramFactory->createHistogram1D( VarNames[i], numberOfPoints, lowerBin, higerBin );
481 pTree->cd(
"/" + name() +
"/" +
"NNInputs");
482 if( !pTree->cd(
"udsJets/" ) )
484 pTree->mkdir(
"udsJets/" ) ;
485 pTree->cd(
"udsJets/" ) ;
488 _inputsHistogramsUDSJets[VarNames[i]]=pHistogramFactory->createHistogram1D( VarNames[i], numberOfPoints, lowerBin, higerBin );
500 void DSTAIDAPlotProcessor::FillTagPlots( LCEvent* pEvent,
unsigned int jetNumber)
506 LCCollection* jetCollection=pEvent->getCollection( _JetCollectionName );
510 ReconstructedParticle* pJet=
dynamic_cast<ReconstructedParticle*
>( jetCollection->getElementAt(jetNumber) );
513 PIDHandler pidh( jetCollection ) ;
518 int mcid = pidh.getAlgorithmID(
"MCTruth");
521 const ParticleID& mcpid = pidh.getParticleID(pJet,mcid);
524 FloatVec mcparams = mcpid.getParameters();
527 int jetType = int(mcparams[pidh.getParameterIndex(mcid,
"TrueJetFlavour")]);
530 int alid = pidh.getAlgorithmID(
"LCFIFlavourTag");
533 const ParticleID& pid = pidh.getParticleID(pJet,alid);
536 FloatVec params = pid.getParameters();
539 double bTag= params[pidh.getParameterIndex(alid,
"BTag")];
540 double cTag= params[pidh.getParameterIndex(alid,
"CTag")];
541 double cTagBBack= params[pidh.getParameterIndex(alid,
"BCTag")];
543 unsigned int NumVertices = int(params[pidh.getParameterIndex(alid,
"NumVertices")]);
548 std::string nvname = _VertexCatNames[ (NumVertices>=N_VERTEX_CATEGORIES) ? (N_VERTEX_CATEGORIES) : (NumVertices)];
552 if( jetType==B_JET ) {
554 if( bTag<=1 && bTag>=0 )
556 _pBJetBTag[nvname]->fill( bTag );
560 _pBJetBTag[nvname]->fill( -0.005 );
563 if( cTag<=1 && cTag>=0 )
565 _pBJetCTag[nvname]->fill( cTag );
569 _pBJetCTag[nvname]->fill( -0.005 );
571 if( cTagBBack<=1 && cTagBBack>=0 )
573 _pBJetBCTag[nvname]->fill( cTagBBack );
577 _pBJetBCTag[nvname]->fill( -0.005 );
580 }
else if( jetType==C_JET ) {
582 if( bTag<=1 && bTag>=0 )
584 _pCJetBTag[nvname]->fill( bTag );
588 _pCJetBTag[nvname]->fill( -0.005 );
591 if( cTag<=1 && cTag>=0 )
593 _pCJetCTag[nvname]->fill( cTag );
597 _pCJetCTag[nvname]->fill( -0.005 );
600 if( cTagBBack<=1 && cTagBBack>=0 )
602 _pCJetBCTag[nvname]->fill( cTagBBack );
606 _pCJetBCTag[nvname]->fill( -0.005 );
609 if( bTag<=1 && bTag>=0 )
611 _pLightJetBTag[nvname]->fill( bTag );
615 _pLightJetBTag[nvname]->fill( -0.005 );
618 if( cTag<=1 && cTag>=0 )
620 _pLightJetCTag[nvname]->fill( cTag );
624 _pLightJetCTag[nvname]->fill( -0.005 );
627 if( cTagBBack<=1 && cTagBBack>=0 )
629 _pLightJetBCTag[nvname]->fill( cTagBBack );
633 _pLightJetBCTag[nvname]->fill( -0.005 );
639 void DSTAIDAPlotProcessor::FillInputsPlots( LCEvent* pEvent,
unsigned int jetNumber )
645 LCCollection* jetCollection=pEvent->getCollection( _JetCollectionName );
649 ReconstructedParticle* pJet=
dynamic_cast<ReconstructedParticle*
>( jetCollection->getElementAt(jetNumber) );
652 PIDHandler pidh( jetCollection ) ;
656 int mcid = pidh.getAlgorithmID(
"MCTruth");
659 const ParticleID& mcpid = pidh.getParticleID(pJet,mcid);
662 FloatVec mcparams = mcpid.getParameters();
665 int jetType = int(mcparams[pidh.getParameterIndex(mcid,
"TrueJetFlavour")]);
668 int alid = pidh.getAlgorithmID(
"LCFIFlavourTag");
671 const ParticleID& pid = pidh.getParticleID(pJet,alid);
674 FloatVec params = pid.getParameters();
679 int NumVertices = int(params[pidh.getParameterIndex(alid,
"NumVertices")]);
683 float JointProbRPhi=-999.;
684 float JointProbZ=-999.;
685 int NumTracksInVertices=-999;
686 float DecayLength=-999.;
687 float DecayLengthSignificance=-999.;
688 float RawMomentum=-999.;
689 float PTCorrectedMass=-999.;
690 float SecondaryVertexProbability=-999.;
694 JointProbRPhi= params[pidh.getParameterIndex(alid,
"JointProbRPhi")];
695 JointProbZ= params[pidh.getParameterIndex(alid,
"JointProbZ")];
696 NumTracksInVertices= int(params[pidh.getParameterIndex(alid,
"NumTracksInVertices")]);
697 DecayLength=params[pidh.getParameterIndex(alid,
"DecayLength")];
698 DecayLengthSignificance= params[pidh.getParameterIndex(alid,
"DecayLengthSignificance")];
699 RawMomentum= params[pidh.getParameterIndex(alid,
"RawMomentum")];
700 PTCorrectedMass= params[pidh.getParameterIndex(alid,
"PTCorrectedMass")];
701 SecondaryVertexProbability=params[pidh.getParameterIndex(alid,
"SecondaryVertexProbability")];
706 _pMyTuple->fill( 0, jetType );
707 _pMyTuple->fill( 1, NumVertices );
708 if(NumTracksInVertices > 1)
710 _pMyTuple->fill( 2, NumTracksInVertices );
711 _pMyTuple->fill( 3, DecayLength);
712 _pMyTuple->fill( 4, DecayLengthSignificance);
713 _pMyTuple->fill( 5, JointProbRPhi);
714 _pMyTuple->fill( 6, JointProbZ);
715 _pMyTuple->fill( 7, PTCorrectedMass);
716 _pMyTuple->fill( 8, RawMomentum);
717 _pMyTuple->fill( 9, SecondaryVertexProbability);
724 if( jetType==B_JET ) _inputsHistogramsBJets[
"NumVertices"]->fill(NumVertices);
725 else if( jetType==C_JET ) _inputsHistogramsCJets[
"NumVertices"]->fill(NumVertices);
726 else _inputsHistogramsUDSJets[
"NumVertices"]->fill(NumVertices);
731 if( jetType==B_JET ) _inputsHistogramsBJets[
"NumTracksInVertices"]->fill(NumTracksInVertices);
732 else if( jetType==C_JET ) _inputsHistogramsCJets[
"NumTracksInVertices"]->fill(NumTracksInVertices);
733 else _inputsHistogramsUDSJets[
"NumTracksInVertices"]->fill(NumTracksInVertices);
735 if( jetType==B_JET ) _inputsHistogramsBJets[
"DecayLength"]->fill(DecayLength);
736 else if( jetType==C_JET ) _inputsHistogramsCJets[
"DecayLength"]->fill(DecayLength);
737 else _inputsHistogramsUDSJets[
"DecayLength"]->fill(DecayLength);
739 if( jetType==B_JET ) _inputsHistogramsBJets[
"DecayLengthSignificance"]->fill(DecayLengthSignificance);
740 else if( jetType==C_JET ) _inputsHistogramsCJets[
"DecayLengthSignificance"]->fill(DecayLengthSignificance);
741 else _inputsHistogramsUDSJets[
"DecayLengthSignificance"]->fill(DecayLengthSignificance);
743 if( jetType==B_JET ) _inputsHistogramsBJets[
"JointProbRPhi"]->fill(JointProbRPhi);
744 else if( jetType==C_JET ) _inputsHistogramsCJets[
"JointProbRPhi"]->fill(JointProbRPhi);
745 else _inputsHistogramsUDSJets[
"JointProbRPhi"]->fill(JointProbRPhi);
747 if( jetType==B_JET ) _inputsHistogramsBJets[
"JointProbZ"]->fill(JointProbZ);
748 else if( jetType==C_JET ) _inputsHistogramsCJets[
"JointProbZ"]->fill(JointProbZ);
749 else _inputsHistogramsUDSJets[
"JointProbZ"]->fill(JointProbZ);
751 if( jetType==B_JET ) _inputsHistogramsBJets[
"PTCorrectedMass"]->fill(PTCorrectedMass);
752 else if( jetType==C_JET ) _inputsHistogramsCJets[
"PTCorrectedMass"]->fill(PTCorrectedMass);
753 else _inputsHistogramsUDSJets[
"PTCorrectedMass"]->fill(PTCorrectedMass);
755 if( jetType==B_JET ) _inputsHistogramsBJets[
"RawMomentum"]->fill(RawMomentum);
756 else if( jetType==C_JET ) _inputsHistogramsCJets[
"RawMomentum"]->fill(RawMomentum);
757 else _inputsHistogramsUDSJets[
"RawMomentum"]->fill(RawMomentum);
759 if( jetType==B_JET ) _inputsHistogramsBJets[
"SecondaryVertexProbability"]->fill(SecondaryVertexProbability);
760 else if( jetType==C_JET ) _inputsHistogramsCJets[
"SecondaryVertexProbability"]->fill(SecondaryVertexProbability);
761 else _inputsHistogramsUDSJets[
"SecondaryVertexProbability"]->fill(SecondaryVertexProbability);
773 void DSTAIDAPlotProcessor::CalculateIntegralAndBackgroundPlots() {
775 AIDA::IHistogramFactory* pHistogramFactory=marlin::AIDAProcessor::histogramFactory(
this );
776 AIDA::ITree* pTree=marlin::AIDAProcessor::tree(
this );
779 for (
unsigned int iVertexCat=1; iVertexCat <= N_VERTEX_CATEGORIES; ++iVertexCat ) {
782 _pBJetBTag[_VertexCatNames[0]] -> add(*_pBJetBTag[_VertexCatNames[iVertexCat]]);
783 _pBJetCTag[_VertexCatNames[0]] -> add(*_pBJetCTag[_VertexCatNames[iVertexCat]]);
784 _pBJetBCTag[_VertexCatNames[0]] -> add(*_pBJetBCTag[_VertexCatNames[iVertexCat]]);
785 _pCJetBTag[_VertexCatNames[0]] -> add(*_pCJetBTag[_VertexCatNames[iVertexCat]]);
786 _pCJetCTag[_VertexCatNames[0]] -> add(*_pCJetCTag[_VertexCatNames[iVertexCat]]);
787 _pCJetBCTag[_VertexCatNames[0]] -> add(*_pCJetBCTag[_VertexCatNames[iVertexCat]]);
788 _pLightJetBTag[_VertexCatNames[0]] -> add(*_pLightJetBTag[_VertexCatNames[iVertexCat]]);
789 _pLightJetCTag[_VertexCatNames[0]] -> add(*_pLightJetCTag[_VertexCatNames[iVertexCat]]);
790 _pLightJetBCTag[_VertexCatNames[0]] -> add(*_pLightJetBCTag[_VertexCatNames[iVertexCat]]);
793 for (
unsigned int iVertexCat=0; iVertexCat <= N_VERTEX_CATEGORIES; ++iVertexCat ) {
796 pTree->cd(
"/" + name() +
"/" +_NumVertexCatDir[iVertexCat]);
798 _pBTagBackgroundValues[_VertexCatNames[iVertexCat]] = pHistogramFactory->add(
"Numbers of non-B jets by B-tag NN value. ("+ _VertexCatNames[iVertexCat]+
")",*_pLightJetBTag[_VertexCatNames[iVertexCat]],*_pCJetBTag[_VertexCatNames[iVertexCat]]);
799 _pCTagBackgroundValues[_VertexCatNames[iVertexCat]] = pHistogramFactory->add(
"Numbers of non-C jets by C-tag NN value. ("+ _VertexCatNames[iVertexCat]+
")",*_pLightJetCTag[_VertexCatNames[iVertexCat]],*_pBJetCTag[_VertexCatNames[iVertexCat]]);
800 _pBCTagBackgroundValues[_VertexCatNames[iVertexCat]] = pHistogramFactory->add(
"Numbers of non-C jets by BC-tag NN value. ("+ _VertexCatNames[iVertexCat]+
")",*_pLightJetBCTag[_VertexCatNames[iVertexCat]],*_pBJetBCTag[_VertexCatNames[iVertexCat]]);
806 void DSTAIDAPlotProcessor::CalculateEfficiencyPurityPlots()
808 AIDA::IHistogramFactory* pHistogramFactory=marlin::AIDAProcessor::histogramFactory(
this );
809 AIDA::ITree* pTree=marlin::AIDAProcessor::tree(
this );
812 AIDA::IDataPointSetFactory* pDataPointSetFactory=marlin::AIDAProcessor::dataPointSetFactory(
this);
816 for (
unsigned int iVertexCat=0; iVertexCat <= N_VERTEX_CATEGORIES; ++iVertexCat ) {
818 pTree->cd(
"/" + name() +
"/" +_NumVertexCatDir[iVertexCat] );
820 std::string nvname = _VertexCatNames[iVertexCat];
823 AIDA::IDataPointSet* _pBJetBTagEfficiency = CreateEfficiencyPlot( _pBJetBTag[nvname] , pDataPointSetFactory->create(
"B-Tag efficiency ("+ nvname +
")",2));
824 AIDA::IDataPointSet* _pCJetCTagEfficiency = CreateEfficiencyPlot( _pCJetCTag[nvname] , pDataPointSetFactory->create(
"C-Tag efficiency ("+ nvname +
")",2));
825 AIDA::IDataPointSet* _pCJetBCTagEfficiency = CreateEfficiencyPlot( _pCJetBCTag[nvname] , pDataPointSetFactory->create(
"BC-Tag efficiency ("+ nvname +
")",2));
829 _pBJetBTagIntegral[nvname] =
830 CreateIntegralHistogram( _pBJetBTag[nvname],
831 pHistogramFactory->createHistogram1D(
"B-Jets: Numbers of events passing B-Tag NN Cut ("+ nvname +
") (DON'T TRUST ERRORS!)",
832 _pBJetBTag[nvname]->axis().bins(),_pBJetBTag[nvname]->axis().lowerEdge(),_pBJetBTag[nvname]->axis().upperEdge()));
834 _pCJetBTagIntegral[nvname] =
835 CreateIntegralHistogram( _pCJetBTag[nvname],
836 pHistogramFactory->createHistogram1D(
"C-Jets: Numbers of events passing B-Tag NN Cut ("+ nvname +
") (DON'T TRUST ERRORS!)",
837 _pCJetBTag[nvname]->axis().bins(),_pCJetBTag[nvname]->axis().lowerEdge(),_pCJetBTag[nvname]->axis().upperEdge()));
839 _pLightJetBTagIntegral[nvname] =
840 CreateIntegralHistogram( _pLightJetBTag[nvname],
841 pHistogramFactory->createHistogram1D(
"Light-Jets: Numbers of events passing B-Tag NN Cut ("+ nvname +
") (DON'T TRUST ERRORS!)",
842 _pLightJetBTag[nvname]->axis().bins(),_pLightJetBTag[nvname]->axis().lowerEdge(),_pLightJetBTag[nvname]->axis().upperEdge()));
844 _pBJetCTagIntegral[nvname] =
845 CreateIntegralHistogram( _pBJetCTag[nvname],
846 pHistogramFactory->createHistogram1D(
"B-Jets: Numbers of events passing C-Tag NN Cut ("+ nvname +
") (DON'T TRUST ERRORS!)",
847 _pBJetCTag[nvname]->axis().bins(),_pBJetCTag[nvname]->axis().lowerEdge(),_pBJetCTag[nvname]->axis().upperEdge()));
849 _pCJetCTagIntegral[nvname] =
850 CreateIntegralHistogram( _pCJetCTag[nvname],
851 pHistogramFactory->createHistogram1D(
"C-Jets: Numbers of events passing C-Tag NN Cut ("+ nvname +
") (DON'T TRUST ERRORS!)",
852 _pCJetCTag[nvname]->axis().bins(),_pCJetCTag[nvname]->axis().lowerEdge(),_pCJetCTag[nvname]->axis().upperEdge()));
854 _pLightJetCTagIntegral[nvname] =
855 CreateIntegralHistogram( _pLightJetCTag[nvname],
856 pHistogramFactory->createHistogram1D(
"Light-Jets: Numbers of events passing C-Tag NN Cut ("+ nvname +
") (DON'T TRUST ERRORS!)",
857 _pLightJetCTag[nvname]->axis().bins(),_pLightJetCTag[nvname]->axis().lowerEdge(),_pLightJetCTag[nvname]->axis().upperEdge()));
859 _pBJetBCTagIntegral[nvname] =
860 CreateIntegralHistogram( _pBJetBCTag[nvname],
861 pHistogramFactory->createHistogram1D(
"B-Jets: Numbers of events passing BC-Tag NN Cut ("+ nvname +
") (DON'T TRUST ERRORS!)",
862 _pBJetBCTag[nvname]->axis().bins(),_pBJetBCTag[nvname]->axis().lowerEdge(),_pBJetBCTag[nvname]->axis().upperEdge()));
864 _pCJetBCTagIntegral[nvname] =
865 CreateIntegralHistogram( _pCJetBCTag[nvname],
866 pHistogramFactory->createHistogram1D(
"C-Jets: Numbers of events passing BC-Tag NN Cut ("+ nvname +
") (DON'T TRUST ERRORS!)",
867 _pCJetBCTag[nvname]->axis().bins(),_pCJetBCTag[nvname]->axis().lowerEdge(),_pCJetBCTag[nvname]->axis().upperEdge()));
869 _pLightJetBCTagIntegral[nvname] =
870 CreateIntegralHistogram( _pLightJetBCTag[nvname],
871 pHistogramFactory->createHistogram1D(
"Light-Jets: Numbers of events passing BC-Tag NN Cut ("+ nvname +
")",
872 _pLightJetBCTag[nvname]->axis().bins(),_pLightJetBCTag[nvname]->axis().lowerEdge(),_pLightJetBCTag[nvname]->axis().upperEdge()));
877 AIDA::IDataPointSet* _pBJetBTagPurity = CreatePurityPlot( _pBJetBTag[nvname], _pBTagBackgroundValues[nvname] , pDataPointSetFactory->create(
"B-Jet purity for B-Tag ("+ nvname +
")",2));
878 AIDA::IDataPointSet* _pCJetCTagPurity = CreatePurityPlot( _pCJetCTag[nvname], _pCTagBackgroundValues[nvname] , pDataPointSetFactory->create(
"C-Jet purity for C-Tag ("+ nvname +
")",2));
879 AIDA::IDataPointSet* _pCJetBCTagPurity = CreatePurityPlot( _pCJetBCTag[nvname], _pBJetBCTag[nvname], pDataPointSetFactory->create(
"C-Jet purity for BC-Tag ("+ nvname +
")",2));
881 AIDA::IDataPointSet* _pCJetBTagLeakage = CreateLeakageRatePlot( _pCJetBTag[nvname], pDataPointSetFactory->create(
"C-Jets: Leakage Rate into B-Tag Sample ("+ nvname +
")",2));
882 AIDA::IDataPointSet* _pLightJetBTagLeakage = CreateLeakageRatePlot( _pLightJetBTag[nvname], pDataPointSetFactory->create(
"Light-Jets: Leakage Rate into B-Tag Sample ("+ nvname +
")",2));
883 AIDA::IDataPointSet* _pBJetCTagLeakage = CreateLeakageRatePlot( _pBJetCTag[nvname], pDataPointSetFactory->create(
"B-Jets: Leakage Rate into C-Tag Sample ("+ nvname +
")",2));
884 AIDA::IDataPointSet* _pLightJetCTagLeakage = CreateLeakageRatePlot( _pLightJetCTag[nvname], pDataPointSetFactory->create(
"Light-Jets: Leakage Rate into C-Tag Sample ("+ nvname +
")",2));
885 AIDA::IDataPointSet* _pBJetBCTagLeakage = CreateLeakageRatePlot( _pBJetBCTag[nvname], pDataPointSetFactory->create(
"B-Jets: Leakage Rate into BC-Tag Sample ("+ nvname +
")",2));
886 AIDA::IDataPointSet* _pLightJetBCTagLeakage = CreateLeakageRatePlot( _pLightJetBCTag[nvname], pDataPointSetFactory->create(
"Light-Jets: Leakage Rate into BC-Tag Sample ("+ nvname +
")",2));
887 AIDA::IDataPointSet* _pNonBJetBTagLeakage = CreateLeakageRatePlot( _pBTagBackgroundValues[nvname], pDataPointSetFactory->create(
"C-Jets: Leakage Rate into B-Tag Sample ("+ nvname +
")",2));
888 AIDA::IDataPointSet* _pNonCJetCTagLeakage = CreateLeakageRatePlot( _pCTagBackgroundValues[nvname], pDataPointSetFactory->create(
"Light-Jets: Leakage Rate into B-Tag Sample ("+ nvname +
")",2));
889 AIDA::IDataPointSet* _pNonCJetBCTagLeakage = CreateLeakageRatePlot( _pBCTagBackgroundValues[nvname], pDataPointSetFactory->create(
"B-Jets: Leakage Rate into C-Tag Sample ("+ nvname +
")",2));
891 CreateXYPlot(_pBJetBTagEfficiency, _pBJetBTagPurity, pDataPointSetFactory->create(
"Purity-Efficiency for B-Tag ("+ nvname +
")",2), 1, 1);
892 CreateXYPlot(_pCJetCTagEfficiency, _pCJetCTagPurity, pDataPointSetFactory->create(
"Purity-Efficiency for C-Tag ("+ nvname +
")",2), 1, 1);
893 CreateXYPlot(_pCJetBCTagEfficiency, _pCJetBCTagPurity, pDataPointSetFactory->create(
"Purity-Efficiency for BC-Tag ("+ nvname +
")",2), 1, 1);
895 CreateXYPlot(_pBJetBTagEfficiency, _pNonBJetBTagLeakage, pDataPointSetFactory->create(
"Leakage Rate-Efficiency for B-Tag ("+ nvname +
")",2), 1, 1);
896 CreateXYPlot(_pCJetCTagEfficiency, _pNonCJetCTagLeakage, pDataPointSetFactory->create(
"Leakage Rate-Efficiency for C-Tag ("+ nvname +
")",2), 1, 1);
897 CreateXYPlot(_pCJetBCTagEfficiency,_pNonCJetBCTagLeakage, pDataPointSetFactory->create(
"Leakage Rate-Efficiency for BC-Tag ("+ nvname +
")",2), 1, 1);
899 CreateXYPlot(_pBJetBTagEfficiency, _pCJetBTagLeakage, pDataPointSetFactory->create(
"Leakage Rate (C into B) vs Efficiency for B-Tag ("+ nvname +
")",2), 1, 1);
900 CreateXYPlot(_pBJetBTagEfficiency, _pLightJetBTagLeakage, pDataPointSetFactory->create(
"Leakage Rate (Light into B) Efficiency for B-Tag ("+ nvname +
")",2), 1, 1);
901 CreateXYPlot(_pCJetCTagEfficiency, _pBJetCTagLeakage, pDataPointSetFactory->create(
"Leakage Rate (B into C) vs Efficiency for C-Tag ("+ nvname +
")",2), 1, 1);
902 CreateXYPlot(_pCJetCTagEfficiency, _pLightJetCTagLeakage, pDataPointSetFactory->create(
"Leakage Rate (Light into C) Efficiency for C-Tag ("+ nvname +
")",2), 1, 1);
903 CreateXYPlot(_pCJetBCTagEfficiency, _pBJetBCTagLeakage, pDataPointSetFactory->create(
"Leakage Rate (B into BC) vs Efficiency for BC-Tag ("+ nvname +
")",2), 1, 1);
904 CreateXYPlot(_pCJetCTagEfficiency, _pLightJetBCTagLeakage, pDataPointSetFactory->create(
"Leakage Rate (Light into BC) Efficiency for BC-Tag ("+ nvname +
")",2), 1, 1);
914 AIDA::IHistogram1D* DSTAIDAPlotProcessor::CreateIntegralHistogram(
const AIDA::IHistogram1D* pNN, AIDA::IHistogram1D* pIntegral)
918 const int numberOfBins=pNN->axis().bins();
920 double integral=pNN->binHeight(AIDA::IAxis::OVERFLOW_BIN);
921 pIntegral->fill(pNN->axis().binLowerEdge(AIDA::IAxis::OVERFLOW_BIN)+pNN->axis().binWidth(numberOfBins-1)/2.,integral);
923 for(
int binNumber=numberOfBins-1; binNumber>=0; --binNumber )
925 integral+= pNN->binHeight( binNumber );
926 pIntegral->fill( pNN->axis().binLowerEdge(binNumber)+pNN->axis().binWidth(binNumber)/2.,integral);
929 integral+= pNN->binHeight(AIDA::IAxis::UNDERFLOW_BIN);
930 pIntegral->fill(pNN->axis().binUpperEdge(AIDA::IAxis::UNDERFLOW_BIN)-pNN->axis().binWidth(0)/2.,integral);
938 AIDA::IDataPointSet* DSTAIDAPlotProcessor::CreateEfficiencyPlot(
const AIDA::IHistogram1D* pSignal, AIDA::IDataPointSet* pDataPointSet)
941 double totalSignal=pSignal->sumBinHeights();
942 double signalPassedCut=0;
944 const int numberOfBins=pSignal->axis().bins();
947 for(
int binNumber=numberOfBins-1; binNumber>=0; --binNumber, iPoint++ )
949 signalPassedCut+=pSignal->binHeight( binNumber );
951 double efficiency = signalPassedCut/totalSignal;
953 double efficiencyError = efficiency * (1. - efficiency) / totalSignal;
954 if (efficiencyError>0) efficiencyError = sqrt(efficiencyError);
957 pDataPointSet->addPoint();
958 pDataPointSet->point(iPoint)->coordinate(0)->setValue(pSignal->axis().binLowerEdge(binNumber)+pSignal->axis().binWidth(binNumber)/2.);
959 pDataPointSet->point(iPoint)->coordinate(1)->setValue( efficiency );
960 pDataPointSet->point(iPoint)->coordinate(1)->setErrorPlus(efficiencyError);
961 pDataPointSet->point(iPoint)->coordinate(1)->setErrorMinus(efficiencyError);
964 return pDataPointSet;
967 AIDA::IDataPointSet* DSTAIDAPlotProcessor::CreatePurityPlot(
const AIDA::IHistogram1D* pSignal,
const AIDA::IHistogram1D* pBackground, AIDA::IDataPointSet* pDataPointSet)
969 const int numberOfBins=pSignal->axis().bins();
972 double signalPassedCut=0;
973 double backgroundPassedCut=0;
975 for (
int binNumber = numberOfBins-1; binNumber >= 0 ; --binNumber, iPoint++ )
978 signalPassedCut+=pSignal->binHeight( binNumber );
979 backgroundPassedCut+=pBackground->binHeight( binNumber );
981 double purity = signalPassedCut/(signalPassedCut+backgroundPassedCut);
982 double purityError = purity * (1. - purity) / (signalPassedCut+backgroundPassedCut);
983 if (purityError>0) purityError = sqrt(purityError);
986 pDataPointSet->addPoint();
987 pDataPointSet->point(iPoint)->coordinate(0)->setValue(pSignal->axis().binLowerEdge(binNumber)+pSignal->axis().binWidth(binNumber)/2.);
988 pDataPointSet->point(iPoint)->coordinate(1)->setValue(purity);
989 pDataPointSet->point(iPoint)->coordinate(1)->setErrorPlus(purityError);
990 pDataPointSet->point(iPoint)->coordinate(1)->setErrorMinus(purityError);
995 return pDataPointSet;
997 AIDA::IDataPointSet* DSTAIDAPlotProcessor::CreateLeakageRatePlot(
const AIDA::IHistogram1D* pBackground, AIDA::IDataPointSet* pDataPointSet)
1000 double totalBackground = pBackground->sumBinHeights();
1001 double backgroundPassedCut=0;
1003 const int numberOfBins=pBackground->axis().bins();
1006 for(
int binNumber=numberOfBins-1; binNumber>=0; --binNumber , iPoint++ )
1009 backgroundPassedCut+=pBackground->binHeight( binNumber );
1011 double leakageRate = backgroundPassedCut/totalBackground;
1013 double leakageRateError = leakageRate * (1. - leakageRate) / totalBackground;
1014 if (leakageRateError>0) leakageRateError = sqrt(leakageRateError);
1017 pDataPointSet->addPoint();
1018 pDataPointSet->point(iPoint)->coordinate(0)->setValue(pBackground->axis().binLowerEdge(binNumber)+pBackground->axis().binWidth(binNumber)/2.);
1019 pDataPointSet->point(iPoint)->coordinate(1)->setValue(leakageRate);
1020 pDataPointSet->point(iPoint)->coordinate(1)->setErrorPlus(leakageRateError);
1021 pDataPointSet->point(iPoint)->coordinate(1)->setErrorMinus(leakageRateError);
1024 return pDataPointSet;
1028 AIDA::IDataPointSet* DSTAIDAPlotProcessor::CreateXYPlot(
const AIDA::IDataPointSet* pDataPointSet0,
const AIDA::IDataPointSet* pDataPointSet1, AIDA::IDataPointSet* xyPointSet,
const int dim0,
const int dim1 )
1032 if (pDataPointSet0->size() == pDataPointSet1->size()) {
1034 for (
int iPoint = 0 ; iPoint != pDataPointSet1->size(); iPoint++)
1036 xyPointSet->addPoint();
1037 xyPointSet->point(iPoint)->coordinate(0)->setValue(pDataPointSet0->point(iPoint)->coordinate(dim0)->value());
1038 xyPointSet->point(iPoint)->coordinate(1)->setValue(pDataPointSet1->point(iPoint)->coordinate(dim1)->value());
1039 xyPointSet->point(iPoint)->coordinate(0)->setErrorPlus(pDataPointSet0->point(iPoint)->coordinate(dim0)->errorPlus());
1040 xyPointSet->point(iPoint)->coordinate(1)->setErrorPlus(pDataPointSet1->point(iPoint)->coordinate(dim1)->errorPlus());
1041 xyPointSet->point(iPoint)->coordinate(0)->setErrorMinus(pDataPointSet0->point(iPoint)->coordinate(dim0)->errorMinus());
1042 xyPointSet->point(iPoint)->coordinate(1)->setErrorMinus(pDataPointSet1->point(iPoint)->coordinate(dim1)->errorMinus());
1054 #endif // end of "#ifndef MARLIN_USE_AIDA"
void _displayCollectionNames(lcio::LCEvent *pEvent)
bool _passesEventCuts(lcio::LCEvent *pEvent)
bool _passesJetCuts(lcio::ReconstructedParticle *pJet)