5 #ifndef MARLIN_USE_AIDA
7 #warning "--------------------------------------------------------------------------------"
8 #warning "- LCFIAIDAPlotProcessor requires MARLIN_USE_AIDA to be defined. Did you enable -"
9 #warning "- AIDA when compiling Marlin? LCFIAIDAPlotProcessor will not be compiled. -"
10 #warning "--------------------------------------------------------------------------------"
15 #include "LCFIAIDAPlotProcessor.h"
33 #include "EVENT/LCCollection.h"
34 #include "EVENT/LCParameters.h"
35 #include "EVENT/LCIntVec.h"
36 #include "EVENT/LCFloatVec.h"
37 #include "EVENT/Vertex.h"
38 #include "EVENT/MCParticle.h"
39 #include "EVENT/Track.h"
43 #include <marlin/Exceptions.h>
46 #include <marlin/AIDAProcessor.h>
47 #include <AIDA/IMeasurement.h>
48 #include <AIDA/IHistogramFactory.h>
49 #include <AIDA/IHistogram1D.h>
50 #include <AIDA/IHistogram2D.h>
51 #include <AIDA/IAxis.h>
52 #include <AIDA/ITree.h>
53 #include <AIDA/ITupleFactory.h>
54 #include <AIDA/ICloud2D.h>
55 #include <UTIL/LCRelationNavigator.h>
61 #warning "USING_RAIDA defined"
64 #warning "USING_JAIDA defined"
68 #ifdef USING_JAIDA//Data point sets aren't implemented in RAIDA - which is a shame as they have functionality not given by histograms
70 #include <AIDA/IDataPointSet.h>
71 #include <AIDA/IDataPointSetFactory.h>
72 #include <AIDA/IDataPoint.h>
75 #include "TypesafeCollection.h"
81 LCFIAIDAPlotProcessor::LCFIAIDAPlotProcessor() : marlin::Processor(
"LCFIAIDAPlotProcessor" )
84 _description=
"Creates an AIDA plot of the LCFIVertex tagging efficiency-purity values and various other things. Make sure that MarlinAIDAProcessor is run before this.";
86 registerInputCollection( LCIO::RECONSTRUCTEDPARTICLE,
88 "Name of the collection of ReconstructedParticles that is the jet" ,
90 std::string(
"FTSelectedJets") );
92 std::vector<std::string> FlavourTagCollectionNamesDefault;
93 FlavourTagCollectionNamesDefault.push_back(
"FlavourTag");
94 registerProcessorParameter(
"FlavourTagCollections" ,
95 "Names of the LCFloatVec Collections that contain the flavour tags (one purity efficiency plot per tag) (in same order as jet collection)" ,
96 _FlavourTagCollectionNames,
97 FlavourTagCollectionNamesDefault) ;
99 registerOptionalParameter(
"TrueJetFlavourCollection" ,
100 "Name of the LCFloatVec collection containing the true flavour of the jets (same order as jets)" ,
101 _TrueJetFlavourColName ,
102 std::string(
"TrueJetFlavour") ) ;
104 registerInputCollection( lcio::LCIO::MCPARTICLE,
105 "MCParticleCollection" ,
106 "Name of the collection that holds all MC particles. " ,
108 std::string(
"MCParticle") ) ;
110 registerOptionalParameter(
"VertexCollection",
111 "Name of the collection that holds the Vertices",
113 std::string(
"ZVRESVertices") ) ;
115 registerOptionalParameter(
"CVertexChargeCollection",
116 "Name of collection containing the vertex charge of the jets, assuming they are C-jets",
117 _CVertexChargeCollection,
118 std::string(
"CCharge") );
120 registerOptionalParameter(
"BVertexChargeCollection",
121 "Name of collection containing the vertex charge of the jets, assuming they are B-jets",
122 _BVertexChargeCollection,
123 std::string(
"BCharge") ) ;
125 registerOptionalParameter(
"TrueTracksToMCPCollection",
126 "Name of collection linking the tracks and the Monte Carlo Particles",
127 _TrueTracksToMCPCollection,
128 std::string(
"LDCTracksMCP") ) ;
130 registerOptionalParameter(
"ZVRESDecayChainCollection" ,
131 "Name of the ZVRES DecayChain collection" ,
132 _ZVRESDecayChainCollection ,
133 std::string(
"ZVRESDecayChains") );
135 registerOptionalParameter(
"ZVRESSelectedJetsCollection" ,
136 "Name of the ZVRES Selected Jets collection" ,
137 _ZVRESSelectedJetsCollection ,
138 std::string(
"ZVRESSelectedJets") );
140 registerOptionalParameter(
"ZVRESDecayChainTrackCollection" ,
141 "Name of the ZVRES Decay Chain Tracks Collection" ,
142 _ZVRESDecayChainRPTracksCollection ,
143 std::string(
"ZVRESDecayChainRPTracks") );
145 FlavourTagCollectionNamesDefault.clear();
146 FlavourTagCollectionNamesDefault.push_back(
"FlavourTagInputs");
147 registerProcessorParameter(
"TagInputsCollections" ,
148 "Names of the LCFloatVec Collections that contain the flavour tag inputs (in same order as jet collection)" ,
149 _FlavourTagInputsCollectionNames,
150 FlavourTagCollectionNamesDefault) ;
152 registerOptionalParameter(
"CosThetaJetMax",
153 "Cut determining the maximum cos(theta) of the jet. Default: |cos(theta)|<0.9" ,
157 registerOptionalParameter(
"CosThetaJetMin",
158 "Cut determining the minimum cos(theta) of the jet. Default: no lower cut." ,
162 registerOptionalParameter(
"PJetMax",
163 "Cut determining the maximum momentum of the jet. Default: 10000 GeV/c" ,
167 registerOptionalParameter(
"PJetMin",
168 "Cut determining the minimum momentum of the jet. Default: no lower cut." ,
172 registerOptionalParameter(
"PrintTrackVertexOutput",
173 "Set true if you want a print-out of the track-vertex association purity",
174 _PrintTrackVertexOutput,
177 registerOptionalParameter(
"MakePurityEfficiencyPlots",
178 "Set true if you want to make the purity-efficiency plots, and leakage rates plots for the various flavour tags",
179 _MakePurityEfficiencyPlots,
182 registerOptionalParameter(
"PrintPurityEfficiencyValues",
183 "Set true if you want a print-out of the purity-efficiency for the various flavour tags",
184 _PrintPurityEfficiencyValues,
187 registerOptionalParameter(
"MakeAdditionalPlots",
188 "Set true if you want to make all other plots (i.e. the non purity-efficiency plots) and the other functionality provided by LCFIAIDAPlotProcessor",
189 _MakeAdditionalPlots,
192 registerOptionalParameter(
"PurityEfficiencyOutputFile" ,
193 "Output filename for the Purity-Efficiency values. Only used if PrintPurityEfficiencyValues parameter is true. If left blank, output will be directed to standard out.",
194 _PurityEfficiencyOutputFile,
195 std::string(
"PurityEfficiencyOutput.txt") ) ;
197 registerOptionalParameter(
"TrackVertexOutputFile" ,
198 "Output filename for the table of the Track-Vertex association. Only used if PrintTrackVertexOutput parameter is true. If left blank, output will be directed to standard out.",
199 _TrackVertexOutputFile,
200 std::string(
"TrackVertexOutput.txt") ) ;
202 registerOptionalParameter(
"MakeTuple",
203 "Set true to make a tuple of the flavour tag input variables. Default is false (only works with jaida).",
207 registerOptionalParameter(
"CTagNNCut",
208 "Cut determining the Neural Net cut used to select C-Jets",
212 registerOptionalParameter(
"BTagNNCut",
213 "Cut determining the Neural Net cut used to select B-Jets",
217 registerOptionalParameter(
"UseFlavourTagCollectionForVertexCharge",
218 "Integer parameter determing which FlavourTag Collection to use the determine C-Jets and B-Jets in Vertex Charge Plots",
219 _iVertexChargeTagCollection,
224 LCFIAIDAPlotProcessor::~LCFIAIDAPlotProcessor()
228 void LCFIAIDAPlotProcessor::init()
231 if ((_iVertexChargeTagCollection >=
int(_FlavourTagCollectionNames.size()) || _iVertexChargeTagCollection < 0) && _FlavourTagCollectionNames.size()!=0) {
232 std::cerr <<
" In " << __FILE__ <<
"(" << __LINE__ <<
"): Invalid parameter for UseFlavourTagCollectionForVertexCharge. Setting to 0." << std::endl;
233 _myVertexChargeTagCollection = 0;
234 }
else if (_FlavourTagCollectionNames.size()==0) {
235 _myVertexChargeTagCollection = 0;
237 _myVertexChargeTagCollection = unsigned(_iVertexChargeTagCollection);
240 _ZoomedVarNames.push_back(
"D0Significance1");
241 _ZoomedVarNames.push_back(
"D0Significance2");
242 _ZoomedVarNames.push_back(
"Z0Significance1");
243 _ZoomedVarNames.push_back(
"Z0Significance2");
245 _VertexCatNames.resize(N_VERTEX_CATEGORIES+1);
246 _VertexCatNames[0]=
"AnyNumberOfVertices";
247 _VertexCatNames[1]=
"OneVertex";
248 _VertexCatNames[2]=
"TwoVertices";
249 _VertexCatNames[3]=
"ThreeOrMoreVertices";
252 _NumVertexCatDir.resize(N_VERTEX_CATEGORIES+1);
253 _NumVertexCatDir[1]=
"OneVertex";
254 _NumVertexCatDir[2]=
"TwoVertices";
255 _NumVertexCatDir[3]=
"ThreeOrMoreVertices";
256 _NumVertexCatDir[0]=
"AnyNumberOfVertices";
263 AIDA::IHistogramFactory* pHistogramFactory=marlin::AIDAProcessor::histogramFactory(
this );
264 AIDA::ITree* pTree=marlin::AIDAProcessor::tree(
this );
266 if( pHistogramFactory!=0 )
268 if (!(pTree->cd(
"/" + name() +
"/"))) {
269 pTree->mkdir(
"/" + name() +
"/" );
270 pTree->cd(
"/" + name() +
"/");
274 if (_MakePurityEfficiencyPlots) CreateTagPlots();
275 if (_MakeAdditionalPlots) CreateAdditionalPlots();
276 if (_MakeTuple) CreateFlavourTagTuple();
280 std::cerr <<
"### " << __FILE__ <<
"(" << __LINE__ <<
"): Unable to get the histogram factory! No histograms will be made."<< std::endl;
283 _lastRunHeaderProcessed=-1;
284 _suppressOutputForRun=-1;
287 InternalVectorInitialisation();
291 void LCFIAIDAPlotProcessor::processRunHeader( LCRunHeader* pRun )
299 _lastRunHeaderProcessed=pRun->getRunNumber();
301 InitialiseFlavourTagInputs(pRun);
307 for (
unsigned int iTag=0; iTag < _FlavourTagCollectionNames.size(); ++iTag)
309 std::vector<std::string> VarNames;
310 (pRun->parameters()).getStringVals(_FlavourTagCollectionNames[iTag],VarNames);
314 std::set<std::string> AvailableNames;
315 std::map<std::string,unsigned int> IndexOf;
317 for (
size_t i = 0;i < VarNames.size();++i)
319 AvailableNames.insert(VarNames[i]);
320 IndexOf[VarNames[i]] = i;
324 _IndexOfForEachTag.push_back(IndexOf);
327 std::set<std::string> RequiredNames;
328 RequiredNames.insert(
"BTag");
329 RequiredNames.insert(
"CTag");
330 RequiredNames.insert(
"BCTag");
332 if (!std::includes(AvailableNames.begin(),AvailableNames.end(),RequiredNames.begin(),RequiredNames.end()))
334 std::cerr << __FILE__ <<
"(" << __LINE__ <<
"): The collection \"" << _FlavourTagCollectionNames[iTag]
335 <<
"\" (if it exists) does not contain the tag values required by " << type() <<
"." << std::endl;
336 std::cerr << __FILE__ <<
"(" << __LINE__ <<
"): The collection \"" << _FlavourTagCollectionNames[iTag]
337 <<
"\" (if it exists) does not contain the tag values required by " << type() <<
"." << std::endl;
341 if (_MakeTuple) CreateFlavourTagInputPlots(pRun);
345 void LCFIAIDAPlotProcessor::CreateFlavourTagTuple()
349 AIDA::ITree* pTree=marlin::AIDAProcessor::tree(
this );
354 AIDA::ITupleFactory* pTupleFactory=marlin::AIDAProcessor::tupleFactory(
this );
358 pTree->cd(
"/" + name());
362 std::string columnNames=
"int TrueJetFlavour=-1, int NumberOfVertices=-1, int NumberOfTracksInVertices=-1, float D0Significance1 = -999., float D0Significance2 = -999., float DecayLength = -999., float DecayLength_SeedToIP= -999., float DecayLengthSignificance= -999., float JointProbRPhi= -999., float JointProbZ= -999., float Momentum1= -999.,float Momentum2= -999., float PTCorrectedMass= -999., float RawMomentum= -999., float SecondaryVertexProbability= -999., float Z0Significance1= -999., float Z0Significance2= -999., int BQVtx=-10, int CQVtx=-10";
365 if (!pTree->cd(
"/" + name() +
"/" +_FlavourTagCollectionNames[_myVertexChargeTagCollection] +
"/")) {
366 pTree->cd(
"/" + name() +
"/");
367 if (!pTree->cd( _FlavourTagCollectionNames[_myVertexChargeTagCollection])) {
368 pTree->mkdir( _FlavourTagCollectionNames[_myVertexChargeTagCollection] +
"/");
369 pTree->cd( _FlavourTagCollectionNames[_myVertexChargeTagCollection]+
"/");
373 if (!pTree->cd(
"TupleDir/")) {
374 pTree->mkdir(
"TupleDir/");
375 pTree->cd(
"TupleDir/");
378 _pMyTuple=pTupleFactory->create(
"FlavourTagInputsTuple",
"FlavourTagInputsTuple", columnNames);
385 void LCFIAIDAPlotProcessor::InitialiseFlavourTagInputs(LCRunHeader* pRun )
387 for (
unsigned int iInputCollection=0; iInputCollection < _FlavourTagInputsCollectionNames.size(); ++iInputCollection)
390 std::vector<std::string> VarNames;
391 (pRun->parameters()).getStringVals(_FlavourTagInputsCollectionNames[iInputCollection],VarNames);
393 std::map<std::string,unsigned int> IndexOf;
394 for (
size_t i = 0;i < VarNames.size();++i)
396 IndexOf[VarNames[i]] = i;
399 _InputsIndex.push_back(IndexOf);
402 std::vector<std::string> trueJetFlavourVarNames;
403 (pRun->parameters()).getStringVals(_TrueJetFlavourColName,trueJetFlavourVarNames);
405 for (
size_t i = 0;i < trueJetFlavourVarNames.size();++i)
407 _FlavourIndex[trueJetFlavourVarNames[i]] = i;
414 void LCFIAIDAPlotProcessor::CreateFlavourTagInputPlots(LCRunHeader* pRun )
417 _inputsHistogramsBJets.resize( _FlavourTagInputsCollectionNames.size() );
418 _inputsHistogramsCJets.resize( _FlavourTagInputsCollectionNames.size() );
419 _inputsHistogramsUDSJets.resize( _FlavourTagInputsCollectionNames.size() );
421 _zoomedInputsHistogramsBJets.resize( _FlavourTagInputsCollectionNames.size() );
422 _zoomedInputsHistogramsCJets.resize( _FlavourTagInputsCollectionNames.size() );
423 _zoomedInputsHistogramsUDSJets.resize( _FlavourTagInputsCollectionNames.size() );
425 AIDA::IHistogramFactory* pHistogramFactory=marlin::AIDAProcessor::histogramFactory(
this );
426 AIDA::ITree* pTree=marlin::AIDAProcessor::tree(
this );
428 for (
unsigned int iInputCollection=0; iInputCollection < _FlavourTagInputsCollectionNames.size(); ++iInputCollection)
431 std::vector<std::string> VarNames;
432 (pRun->parameters()).getStringVals(_FlavourTagInputsCollectionNames[iInputCollection],VarNames);
435 for (
size_t i = 0;i < VarNames.size();++i)
439 if( _inputsHistogramsBJets[iInputCollection][VarNames[i]]==0 )
441 pTree->cd(
"/" + name() +
"/");
442 if( !pTree->cd(_FlavourTagInputsCollectionNames[iInputCollection] ) )
444 pTree->mkdir( _FlavourTagInputsCollectionNames[iInputCollection] ) ;
445 pTree->cd(_FlavourTagInputsCollectionNames[iInputCollection] ) ;
448 int numberOfPoints=_numberOfPoints/4;
453 if( VarNames[i]==
"BQVtx" || VarNames[i]==
"CQVtx" )
459 else if( VarNames[i]==
"NumVertices" )
465 else if( VarNames[i]==
"NumTracksInVertices" )
471 else if ( VarNames[i]==
"D0Significance1" || VarNames[i]==
"D0Significance2" )
477 else if ( VarNames[i]==
"Z0Significance1" || VarNames[i]==
"Z0Significance2")
483 else if (VarNames[i]==
"DecayLengthSignificance")
489 else if (VarNames[i]==
"DecayLength" || VarNames[i]==
"DecayLength(SeedToIP)" )
495 else if (VarNames[i]==
"JointProbRPhi" || VarNames[i]==
"JointProbZ"|| VarNames[i]==
"SecondaryVertexProbability")
501 else if (VarNames[i]==
"Momentum1" || VarNames[i]==
"Momentum2" || VarNames[i]==
"RawMomentum" )
507 else if (VarNames[i]==
"PTCorrectedMass" )
514 pTree->cd(
"/" + name() +
"/" + _FlavourTagInputsCollectionNames[iInputCollection]);
515 if( !pTree->cd(
"bJets" ) )
517 pTree->mkdir(
"bJets" ) ;
518 pTree->cd(
"bJets" ) ;
521 _inputsHistogramsBJets[iInputCollection][VarNames[i]]=pHistogramFactory->createHistogram1D( VarNames[i], numberOfPoints, lowerBin, higerBin );
523 pTree->cd(
"/" + name() +
"/" + _FlavourTagInputsCollectionNames[iInputCollection]);
524 if( !pTree->cd(
"cJets" ) )
526 pTree->mkdir(
"cJets" ) ;
527 pTree->cd(
"cJets" ) ;
530 _inputsHistogramsCJets[iInputCollection][VarNames[i]]=pHistogramFactory->createHistogram1D( VarNames[i], numberOfPoints, lowerBin, higerBin );
532 pTree->cd(
"/" + name() +
"/" + _FlavourTagInputsCollectionNames[iInputCollection]);
533 if( !pTree->cd(
"udsJets/" ) )
535 pTree->mkdir(
"udsJets/" ) ;
536 pTree->cd(
"udsJets/" ) ;
539 _inputsHistogramsUDSJets[iInputCollection][VarNames[i]]=pHistogramFactory->createHistogram1D( VarNames[i], numberOfPoints, lowerBin, higerBin );
544 if (isFirstEvent()) {
546 for (
size_t i = 0;i < _ZoomedVarNames.size();++i) {
548 std::string zoomed_name = _ZoomedVarNames[i] +
" (zoomed)";
550 pTree->cd(
"/" + name() +
"/" + _FlavourTagInputsCollectionNames[iInputCollection]);
551 if( !pTree->cd(
"bJets" ) )
553 pTree->mkdir(
"bJets" ) ;
554 pTree->cd(
"bJets" ) ;
557 _zoomedInputsHistogramsBJets[iInputCollection][zoomed_name] = pHistogramFactory->createHistogram1D( zoomed_name, 100, -10., 20.);
559 pTree->cd(
"/" + name() +
"/" + _FlavourTagInputsCollectionNames[iInputCollection]);
560 if( !pTree->cd(
"cJets" ) )
562 pTree->mkdir(
"cJets" ) ;
563 pTree->cd(
"cJets" ) ;
566 _zoomedInputsHistogramsCJets[iInputCollection][zoomed_name] = pHistogramFactory->createHistogram1D( zoomed_name, 100, -10., 20.);
569 pTree->cd(
"/" + name() +
"/" + _FlavourTagInputsCollectionNames[iInputCollection]);
570 if( !pTree->cd(
"udsJets" ) )
572 pTree->mkdir(
"udsJets" ) ;
573 pTree->cd(
"udsJets" ) ;
576 _zoomedInputsHistogramsUDSJets[iInputCollection][zoomed_name] = pHistogramFactory->createHistogram1D( zoomed_name, 100, -10., 20.);
583 void LCFIAIDAPlotProcessor::processEvent( LCEvent* pEvent )
587 if( (_lastRunHeaderProcessed != pEvent->getRunNumber()) && (_suppressOutputForRun != pEvent->getRunNumber()) )
589 std::cerr << __FILE__ <<
"(" << __LINE__ <<
"): processRunHeader() was not called for run " << pEvent->getRunNumber()
590 <<
" (did you use \"SkipNEvents\"?). The order of the information in the flavour tag collection(s) is going to be guessed." << std::endl;
593 _suppressOutputForRun=pEvent->getRunNumber();
596 std::map<std::string,unsigned int> guessedOrder;
597 guessedOrder[
"BTag"]=0; guessedOrder[
"CTag"]=1; guessedOrder[
"BCTag"]=2;
598 _IndexOfForEachTag.clear();
599 for (
unsigned int iTag=0; iTag < _FlavourTagCollectionNames.size(); ++iTag) _IndexOfForEachTag.push_back( guessedOrder );
613 if( PassesEventCuts(pEvent) )
616 ReconstructedParticle* pJet;
618 for(
int jetNumber=0; jetNumber<jetCollection.getNumberOfElements(); ++jetNumber )
620 pJet=jetCollection.getElementAt(jetNumber);
623 if( PassesJetCuts(pJet) )
626 if (_MakePurityEfficiencyPlots || _PrintPurityEfficiencyValues) FillTagPlots( pEvent, jetNumber );
627 if (_MakeAdditionalPlots) FillVertexChargePlots( pEvent, jetNumber );
628 if (_MakeTuple) FillInputsPlots( pEvent, jetNumber );
629 if (_MakeAdditionalPlots) FillVertexPlots( pEvent, jetNumber );
631 if (_PrintTrackVertexOutput) FillZVRESTable(pEvent);
640 void LCFIAIDAPlotProcessor::check( LCEvent* pEvent )
646 void LCFIAIDAPlotProcessor::end()
648 if (_MakePurityEfficiencyPlots||_PrintPurityEfficiencyValues) CalculateIntegralAndBackgroundPlots();
649 if (_MakePurityEfficiencyPlots) CalculateEfficiencyPurityPlots();
650 if (_MakeAdditionalPlots) CalculateAdditionalPlots();
651 if (_PrintPurityEfficiencyValues) PrintNNOutput();
652 if (_PrintTrackVertexOutput) PrintZVRESTable();
658 bool LCFIAIDAPlotProcessor::PassesEventCuts( LCEvent* pEvent )
669 bool LCFIAIDAPlotProcessor::PassesJetCuts( ReconstructedParticle* pJet )
683 const double* jetMomentum=pJet->getMomentum();
685 double momentumMagnitude = sqrt(pow(jetMomentum[0],2)+pow(jetMomentum[1],2)+pow(jetMomentum[2],2));
687 double cosTheta = jetMomentum[2] / momentumMagnitude;
688 if( fabs(cosTheta) < _CosThetaJetMin || fabs(cosTheta) > _CosThetaJetMax )
return false;
690 if (momentumMagnitude > _PJetMax || momentumMagnitude < _PJetMin)
return false;
698 void LCFIAIDAPlotProcessor::CalculateIntegralAndBackgroundPlots() {
700 AIDA::IHistogramFactory* pHistogramFactory=marlin::AIDAProcessor::histogramFactory(
this );
701 AIDA::ITree* pTree=marlin::AIDAProcessor::tree(
this );
703 _pBJetBTagIntegral.resize( _FlavourTagCollectionNames.size() );
704 _pCJetBTagIntegral.resize( _FlavourTagCollectionNames.size() );
705 _pLightJetBTagIntegral.resize( _FlavourTagCollectionNames.size() );
706 _pBJetCTagIntegral.resize( _FlavourTagCollectionNames.size() );
707 _pCJetCTagIntegral.resize( _FlavourTagCollectionNames.size() );
708 _pLightJetCTagIntegral.resize( _FlavourTagCollectionNames.size() );
709 _pBJetBCTagIntegral.resize( _FlavourTagCollectionNames.size() );
710 _pCJetBCTagIntegral.resize( _FlavourTagCollectionNames.size() );
711 _pLightJetBCTagIntegral.resize( _FlavourTagCollectionNames.size() );
714 for (
unsigned int iTagCollection=0; iTagCollection < _FlavourTagCollectionNames.size(); ++iTagCollection )
717 for (
unsigned int iVertexCat=1; iVertexCat <= N_VERTEX_CATEGORIES; ++iVertexCat ) {
720 _pBJetBTag[iTagCollection][_VertexCatNames[0]] -> add(*_pBJetBTag[iTagCollection][_VertexCatNames[iVertexCat]]);
721 _pBJetCTag[iTagCollection][_VertexCatNames[0]] -> add(*_pBJetCTag[iTagCollection][_VertexCatNames[iVertexCat]]);
722 _pBJetBCTag[iTagCollection][_VertexCatNames[0]] -> add(*_pBJetBCTag[iTagCollection][_VertexCatNames[iVertexCat]]);
723 _pCJetBTag[iTagCollection][_VertexCatNames[0]] -> add(*_pCJetBTag[iTagCollection][_VertexCatNames[iVertexCat]]);
724 _pCJetCTag[iTagCollection][_VertexCatNames[0]] -> add(*_pCJetCTag[iTagCollection][_VertexCatNames[iVertexCat]]);
725 _pCJetBCTag[iTagCollection][_VertexCatNames[0]] -> add(*_pCJetBCTag[iTagCollection][_VertexCatNames[iVertexCat]]);
726 _pLightJetBTag[iTagCollection][_VertexCatNames[0]] -> add(*_pLightJetBTag[iTagCollection][_VertexCatNames[iVertexCat]]);
727 _pLightJetCTag[iTagCollection][_VertexCatNames[0]] -> add(*_pLightJetCTag[iTagCollection][_VertexCatNames[iVertexCat]]);
728 _pLightJetBCTag[iTagCollection][_VertexCatNames[0]] -> add(*_pLightJetBCTag[iTagCollection][_VertexCatNames[iVertexCat]]);
731 for (
unsigned int iVertexCat=0; iVertexCat <= N_VERTEX_CATEGORIES; ++iVertexCat ) {
734 pTree->cd(
"/" + name() +
"/" + _FlavourTagCollectionNames[iTagCollection]+
"/" +_NumVertexCatDir[iVertexCat]);
736 _pBTagBackgroundValues[iTagCollection][_VertexCatNames[iVertexCat]] = pHistogramFactory->add(
"Numbers of non-B jets by B-tag NN value. ("+ _VertexCatNames[iVertexCat]+
")",*_pLightJetBTag[iTagCollection][_VertexCatNames[iVertexCat]],*_pCJetBTag[iTagCollection][_VertexCatNames[iVertexCat]]);
737 _pCTagBackgroundValues[iTagCollection][_VertexCatNames[iVertexCat]] = pHistogramFactory->add(
"Numbers of non-C jets by C-tag NN value. ("+ _VertexCatNames[iVertexCat]+
")",*_pLightJetCTag[iTagCollection][_VertexCatNames[iVertexCat]],*_pBJetCTag[iTagCollection][_VertexCatNames[iVertexCat]]);
738 _pBCTagBackgroundValues[iTagCollection][_VertexCatNames[iVertexCat]] = pHistogramFactory->add(
"Numbers of non-C jets by BC-tag NN value. ("+ _VertexCatNames[iVertexCat]+
")",*_pLightJetBCTag[iTagCollection][_VertexCatNames[iVertexCat]],*_pBJetBCTag[iTagCollection][_VertexCatNames[iVertexCat]]);
747 void LCFIAIDAPlotProcessor::CalculateEfficiencyPurityPlots()
749 AIDA::IHistogramFactory* pHistogramFactory=marlin::AIDAProcessor::histogramFactory(
this );
750 AIDA::ITree* pTree=marlin::AIDAProcessor::tree(
this );
753 AIDA::IDataPointSetFactory* pDataPointSetFactory=marlin::AIDAProcessor::dataPointSetFactory(
this);
758 for (
unsigned int iTagCollection=0; iTagCollection < _FlavourTagCollectionNames.size(); ++iTagCollection )
761 for (
unsigned int iVertexCat=0; iVertexCat <= N_VERTEX_CATEGORIES; ++iVertexCat ) {
763 pTree->cd(
"/" + name() +
"/" + _FlavourTagCollectionNames[iTagCollection]+
"/" +_NumVertexCatDir[iVertexCat] );
765 std::string nvname = _VertexCatNames[iVertexCat];
768 AIDA::IDataPointSet* _pBJetBTagEfficiency = CreateEfficiencyPlot( _pBJetBTag[iTagCollection][nvname] , pDataPointSetFactory->create(
"B-Tag efficiency ("+ nvname +
")",2));
769 AIDA::IDataPointSet* _pCJetCTagEfficiency = CreateEfficiencyPlot( _pCJetCTag[iTagCollection][nvname] , pDataPointSetFactory->create(
"C-Tag efficiency ("+ nvname +
")",2));
770 AIDA::IDataPointSet* _pCJetBCTagEfficiency = CreateEfficiencyPlot( _pCJetBCTag[iTagCollection][nvname] , pDataPointSetFactory->create(
"BC-Tag efficiency ("+ nvname +
")",2));
774 _pBJetBTagIntegral[iTagCollection][nvname] =
775 CreateIntegralHistogram( _pBJetBTag[iTagCollection][nvname],
776 pHistogramFactory->createHistogram1D(
"B-Jets: Numbers of events passing B-Tag NN Cut ("+ nvname +
") (DON'T TRUST ERRORS!)",
777 _pBJetBTag[iTagCollection][nvname]->axis().bins(),_pBJetBTag[iTagCollection][nvname]->axis().lowerEdge(),_pBJetBTag[iTagCollection][nvname]->axis().upperEdge()));
779 _pCJetBTagIntegral[iTagCollection][nvname] =
780 CreateIntegralHistogram( _pCJetBTag[iTagCollection][nvname],
781 pHistogramFactory->createHistogram1D(
"C-Jets: Numbers of events passing B-Tag NN Cut ("+ nvname +
") (DON'T TRUST ERRORS!)",
782 _pCJetBTag[iTagCollection][nvname]->axis().bins(),_pCJetBTag[iTagCollection][nvname]->axis().lowerEdge(),_pCJetBTag[iTagCollection][nvname]->axis().upperEdge()));
784 _pLightJetBTagIntegral[iTagCollection][nvname] =
785 CreateIntegralHistogram( _pLightJetBTag[iTagCollection][nvname],
786 pHistogramFactory->createHistogram1D(
"Light-Jets: Numbers of events passing B-Tag NN Cut ("+ nvname +
") (DON'T TRUST ERRORS!)",
787 _pLightJetBTag[iTagCollection][nvname]->axis().bins(),_pLightJetBTag[iTagCollection][nvname]->axis().lowerEdge(),_pLightJetBTag[iTagCollection][nvname]->axis().upperEdge()));
789 _pBJetCTagIntegral[iTagCollection][nvname] =
790 CreateIntegralHistogram( _pBJetCTag[iTagCollection][nvname],
791 pHistogramFactory->createHistogram1D(
"B-Jets: Numbers of events passing C-Tag NN Cut ("+ nvname +
") (DON'T TRUST ERRORS!)",
792 _pBJetCTag[iTagCollection][nvname]->axis().bins(),_pBJetCTag[iTagCollection][nvname]->axis().lowerEdge(),_pBJetCTag[iTagCollection][nvname]->axis().upperEdge()));
794 _pCJetCTagIntegral[iTagCollection][nvname] =
795 CreateIntegralHistogram( _pCJetCTag[iTagCollection][nvname],
796 pHistogramFactory->createHistogram1D(
"C-Jets: Numbers of events passing C-Tag NN Cut ("+ nvname +
") (DON'T TRUST ERRORS!)",
797 _pCJetCTag[iTagCollection][nvname]->axis().bins(),_pCJetCTag[iTagCollection][nvname]->axis().lowerEdge(),_pCJetCTag[iTagCollection][nvname]->axis().upperEdge()));
799 _pLightJetCTagIntegral[iTagCollection][nvname] =
800 CreateIntegralHistogram( _pLightJetCTag[iTagCollection][nvname],
801 pHistogramFactory->createHistogram1D(
"Light-Jets: Numbers of events passing C-Tag NN Cut ("+ nvname +
") (DON'T TRUST ERRORS!)",
802 _pLightJetCTag[iTagCollection][nvname]->axis().bins(),_pLightJetCTag[iTagCollection][nvname]->axis().lowerEdge(),_pLightJetCTag[iTagCollection][nvname]->axis().upperEdge()));
804 _pBJetBCTagIntegral[iTagCollection][nvname] =
805 CreateIntegralHistogram( _pBJetBCTag[iTagCollection][nvname],
806 pHistogramFactory->createHistogram1D(
"B-Jets: Numbers of events passing BC-Tag NN Cut ("+ nvname +
") (DON'T TRUST ERRORS!)",
807 _pBJetBCTag[iTagCollection][nvname]->axis().bins(),_pBJetBCTag[iTagCollection][nvname]->axis().lowerEdge(),_pBJetBCTag[iTagCollection][nvname]->axis().upperEdge()));
809 _pCJetBCTagIntegral[iTagCollection][nvname] =
810 CreateIntegralHistogram( _pCJetBCTag[iTagCollection][nvname],
811 pHistogramFactory->createHistogram1D(
"C-Jets: Numbers of events passing BC-Tag NN Cut ("+ nvname +
") (DON'T TRUST ERRORS!)",
812 _pCJetBCTag[iTagCollection][nvname]->axis().bins(),_pCJetBCTag[iTagCollection][nvname]->axis().lowerEdge(),_pCJetBCTag[iTagCollection][nvname]->axis().upperEdge()));
814 _pLightJetBCTagIntegral[iTagCollection][nvname] =
815 CreateIntegralHistogram( _pLightJetBCTag[iTagCollection][nvname],
816 pHistogramFactory->createHistogram1D(
"Light-Jets: Numbers of events passing BC-Tag NN Cut ("+ nvname +
")",
817 _pLightJetBCTag[iTagCollection][nvname]->axis().bins(),_pLightJetBCTag[iTagCollection][nvname]->axis().lowerEdge(),_pLightJetBCTag[iTagCollection][nvname]->axis().upperEdge()));
841 AIDA::IDataPointSet* _pBJetBTagPurity = CreatePurityPlot( _pBJetBTag[iTagCollection][nvname], _pBTagBackgroundValues[iTagCollection][nvname] , pDataPointSetFactory->create(
"B-Jet purity for B-Tag ("+ nvname +
")",2));
842 AIDA::IDataPointSet* _pCJetCTagPurity = CreatePurityPlot( _pCJetCTag[iTagCollection][nvname], _pCTagBackgroundValues[iTagCollection][nvname] , pDataPointSetFactory->create(
"C-Jet purity for C-Tag ("+ nvname +
")",2));
843 AIDA::IDataPointSet* _pCJetBCTagPurity = CreatePurityPlot( _pCJetBCTag[iTagCollection][nvname], _pBJetBCTag[iTagCollection][nvname], pDataPointSetFactory->create(
"C-Jet purity for BC-Tag ("+ nvname +
")",2));
845 AIDA::IDataPointSet* _pCJetBTagLeakage = CreateLeakageRatePlot( _pCJetBTag[iTagCollection][nvname], pDataPointSetFactory->create(
"C-Jets: Leakage Rate into B-Tag Sample ("+ nvname +
")",2));
846 AIDA::IDataPointSet* _pLightJetBTagLeakage = CreateLeakageRatePlot( _pLightJetBTag[iTagCollection][nvname], pDataPointSetFactory->create(
"Light-Jets: Leakage Rate into B-Tag Sample ("+ nvname +
")",2));
847 AIDA::IDataPointSet* _pBJetCTagLeakage = CreateLeakageRatePlot( _pBJetCTag[iTagCollection][nvname], pDataPointSetFactory->create(
"B-Jets: Leakage Rate into C-Tag Sample ("+ nvname +
")",2));
848 AIDA::IDataPointSet* _pLightJetCTagLeakage = CreateLeakageRatePlot( _pLightJetCTag[iTagCollection][nvname], pDataPointSetFactory->create(
"Light-Jets: Leakage Rate into C-Tag Sample ("+ nvname +
")",2));
849 AIDA::IDataPointSet* _pBJetBCTagLeakage = CreateLeakageRatePlot( _pBJetBCTag[iTagCollection][nvname], pDataPointSetFactory->create(
"B-Jets: Leakage Rate into BC-Tag Sample ("+ nvname +
")",2));
850 AIDA::IDataPointSet* _pLightJetBCTagLeakage = CreateLeakageRatePlot( _pLightJetBCTag[iTagCollection][nvname], pDataPointSetFactory->create(
"Light-Jets: Leakage Rate into BC-Tag Sample ("+ nvname +
")",2));
851 AIDA::IDataPointSet* _pNonBJetBTagLeakage = CreateLeakageRatePlot( _pBTagBackgroundValues[iTagCollection][nvname], pDataPointSetFactory->create(
"C-Jets: Leakage Rate into B-Tag Sample ("+ nvname +
")",2));
852 AIDA::IDataPointSet* _pNonCJetCTagLeakage = CreateLeakageRatePlot( _pCTagBackgroundValues[iTagCollection][nvname], pDataPointSetFactory->create(
"Light-Jets: Leakage Rate into B-Tag Sample ("+ nvname +
")",2));
853 AIDA::IDataPointSet* _pNonCJetBCTagLeakage = CreateLeakageRatePlot( _pBCTagBackgroundValues[iTagCollection][nvname], pDataPointSetFactory->create(
"B-Jets: Leakage Rate into C-Tag Sample ("+ nvname +
")",2));
855 CreateXYPlot(_pBJetBTagEfficiency, _pBJetBTagPurity, pDataPointSetFactory->create(
"Purity-Efficiency for B-Tag ("+ nvname +
")",2), 1, 1);
856 CreateXYPlot(_pCJetCTagEfficiency, _pCJetCTagPurity, pDataPointSetFactory->create(
"Purity-Efficiency for C-Tag ("+ nvname +
")",2), 1, 1);
857 CreateXYPlot(_pCJetBCTagEfficiency, _pCJetBCTagPurity, pDataPointSetFactory->create(
"Purity-Efficiency for BC-Tag ("+ nvname +
")",2), 1, 1);
859 CreateXYPlot(_pBJetBTagEfficiency, _pNonBJetBTagLeakage, pDataPointSetFactory->create(
"Leakage Rate-Efficiency for B-Tag ("+ nvname +
")",2), 1, 1);
860 CreateXYPlot(_pCJetCTagEfficiency, _pNonCJetCTagLeakage, pDataPointSetFactory->create(
"Leakage Rate-Efficiency for C-Tag ("+ nvname +
")",2), 1, 1);
861 CreateXYPlot(_pCJetBCTagEfficiency,_pNonCJetBCTagLeakage, pDataPointSetFactory->create(
"Leakage Rate-Efficiency for BC-Tag ("+ nvname +
")",2), 1, 1);
863 CreateXYPlot(_pBJetBTagEfficiency, _pCJetBTagLeakage, pDataPointSetFactory->create(
"Leakage Rate (C into B) vs Efficiency for B-Tag ("+ nvname +
")",2), 1, 1);
864 CreateXYPlot(_pBJetBTagEfficiency, _pLightJetBTagLeakage, pDataPointSetFactory->create(
"Leakage Rate (Light into B) Efficiency for B-Tag ("+ nvname +
")",2), 1, 1);
865 CreateXYPlot(_pCJetCTagEfficiency, _pBJetCTagLeakage, pDataPointSetFactory->create(
"Leakage Rate (B into C) vs Efficiency for C-Tag ("+ nvname +
")",2), 1, 1);
866 CreateXYPlot(_pCJetCTagEfficiency, _pLightJetCTagLeakage, pDataPointSetFactory->create(
"Leakage Rate (Light into C) Efficiency for C-Tag ("+ nvname +
")",2), 1, 1);
867 CreateXYPlot(_pCJetBCTagEfficiency, _pBJetBCTagLeakage, pDataPointSetFactory->create(
"Leakage Rate (B into BC) vs Efficiency for BC-Tag ("+ nvname +
")",2), 1, 1);
868 CreateXYPlot(_pCJetCTagEfficiency, _pLightJetBCTagLeakage, pDataPointSetFactory->create(
"Leakage Rate (Light into BC) Efficiency for BC-Tag ("+ nvname +
")",2), 1, 1);
879 void LCFIAIDAPlotProcessor::CalculateAdditionalPlots()
883 AIDA::ITree* pTree=marlin::AIDAProcessor::tree(
this );
886 AIDA::IDataPointSetFactory* pDataPointSetFactory=marlin::AIDAProcessor::dataPointSetFactory(
this);
891 for (
unsigned int iTagCollection=0; iTagCollection < _FlavourTagCollectionNames.size(); ++iTagCollection )
894 for (
unsigned int iVertexCat=0; iVertexCat <= N_VERTEX_CATEGORIES; ++iVertexCat ) {
896 pTree->cd(
"/" + name() +
"/" + _FlavourTagCollectionNames[iTagCollection]+
"/" +_NumVertexCatDir[iVertexCat] );
898 std::string nvname = _VertexCatNames[iVertexCat];
903 pTree->cd(
"/" + name() +
"/" + _FlavourTagCollectionNames[iTagCollection]+
"/");
906 CreateEfficiencyPlot2( _pCDecayLengthAll[iTagCollection], _pCDecayLengthTwoVertices[iTagCollection], pDataPointSetFactory->create(
"Efficiency to reconstruct secondary vertex vs true decay length for true C-events",2));
907 CreateEfficiencyPlot2( _pBDecayLengthAll[iTagCollection], _pBDecayLengthTwoVertices[iTagCollection], pDataPointSetFactory->create(
"Efficiency to reconstruct secondary vertex vs true decay length for true B-events",2));
913 pTree->cd(
"/" + name() +
"/VertexChargePlots/");
917 AIDA::IDataPointSet* pBJetVtxChargeDPS = pDataPointSetFactory->create(
"B-Jets: Vertex Charge Leakage",2);
918 AIDA::IDataPointSet* pCJetVtxChargeDPS = pDataPointSetFactory->create(
"C-Jets: Vertex Charge Leakage",2);
919 CreateVertexChargeLeakagePlot(pBJetVtxChargeDPS, pCJetVtxChargeDPS);
925 void LCFIAIDAPlotProcessor::FillInputsPlots( LCEvent* pEvent,
unsigned int jetNumber )
927 int jetType=FindTrueJetType( pEvent, jetNumber );
928 if( jetType==0 )
return;
932 for (
unsigned int iInputsCollection=0; iInputsCollection < _FlavourTagInputsCollectionNames.size(); ++iInputsCollection )
935 if( !inputsCollection.is_valid() )
937 std::cerr << _FlavourTagCollectionNames[iInputsCollection] <<
"is not valid, not filling flavour tag input plots" << std::endl;
938 std::cerr <<
"set MakeAdditionalPlots to false to remove this warning" << std::endl;
943 lcio::LCFloatVec* pInputs=inputsCollection.getElementAt( jetNumber );
953 int CQVtx = FindCQVtx(pEvent, jetNumber);
954 int BQVtx = FindBQVtx(pEvent, jetNumber);
956 if (_MakeTuple && iInputsCollection==0) {
960 int NumVertices = int((*pInputs)[_InputsIndex[iInputsCollection][
"NumVertices"]]);
961 int NumTracksInVertices = int((*pInputs)[_InputsIndex[iInputsCollection] [
"NumTracksInVertices"]]);
962 float D0Significance1=(*pInputs)[_InputsIndex[iInputsCollection] [
"D0Significance1"]];
963 float D0Significance2=(*pInputs)[_InputsIndex[iInputsCollection] [
"D0Significance2"]];
964 float DecayLength=(*pInputs)[_InputsIndex[iInputsCollection] [
"DecayLength"]];
965 float DecayLength_SeedToIP=(*pInputs)[_InputsIndex[iInputsCollection][
"DecayLength(SeedToIP)"]];
966 float DecayLengthSignificance=(*pInputs)[_InputsIndex[iInputsCollection] [
"DecayLengthSignificance"]];
967 float JointProbRPhi=(*pInputs)[_InputsIndex[iInputsCollection] [
"JointProbRPhi"]];
968 float JointProbZ=(*pInputs)[_InputsIndex[iInputsCollection] [
"JointProbZ"]];
969 float Momentum1=(*pInputs)[_InputsIndex[iInputsCollection] [
"Momentum1"]];
970 float Momentum2=(*pInputs)[_InputsIndex[iInputsCollection] [
"Momentum2"]];
971 float PTCorrectedMass=(*pInputs)[_InputsIndex[iInputsCollection] [
"PTCorrectedMass"]];
972 float RawMomentum=(*pInputs)[_InputsIndex[iInputsCollection] [
"RawMomentum"]];
973 float SecondaryVertexProbability=(*pInputs)[_InputsIndex[iInputsCollection][
"SecondaryVertexProbability"]];
974 float Z0Significance1=(*pInputs)[_InputsIndex[iInputsCollection] [
"Z0Significance1"]];
975 float Z0Significance2=(*pInputs)[_InputsIndex[iInputsCollection] [
"Z0Significance2"]];
979 _pMyTuple->fill( 0, jetType );
980 _pMyTuple->fill( 1, NumVertices );
981 _pMyTuple->fill( 2, NumTracksInVertices );
982 _pMyTuple->fill( 3, D0Significance1);
983 _pMyTuple->fill( 4, D0Significance2);
984 _pMyTuple->fill( 5, DecayLength);
985 _pMyTuple->fill( 6, DecayLength_SeedToIP);
986 _pMyTuple->fill( 7, DecayLengthSignificance);
987 _pMyTuple->fill( 8, JointProbRPhi);
988 _pMyTuple->fill( 9, JointProbZ);
989 _pMyTuple->fill( 10, Momentum1);
990 _pMyTuple->fill( 11, Momentum2);
991 _pMyTuple->fill( 12, PTCorrectedMass);
992 _pMyTuple->fill( 13, RawMomentum);
993 _pMyTuple->fill( 14, SecondaryVertexProbability);
994 _pMyTuple->fill( 15, Z0Significance1);
995 _pMyTuple->fill( 16 ,Z0Significance2);
997 _pMyTuple->fill( 17, BQVtx );
998 _pMyTuple->fill( 18, CQVtx );
1000 _pMyTuple->addRow();
1006 for( std::map<std::string,unsigned int>::iterator iTagNames=_InputsIndex[iInputsCollection].begin();
1007 iTagNames!=_InputsIndex[iInputsCollection].end(); ++iTagNames ) {
1009 double input=(*pInputs)[(*iTagNames).second];
1012 if (! ((*pInputs)[_InputsIndex[iInputsCollection][
"NumVertices"]] < 2 &&
1013 ((*iTagNames).first ==
"DecayLength" || (*iTagNames).first ==
"RawMomentum" ||
1014 (*iTagNames).first ==
"SecondaryVertexProbability" || (*iTagNames).first ==
"PTCorrectedMass" ||
1015 (*iTagNames).first ==
"DecayLength(SeedToIP)" || (*iTagNames).first ==
"DecayLengthSignificance") )) {
1017 if( jetType==B_JET ) _inputsHistogramsBJets[iInputsCollection][(*iTagNames).first]->fill(input);
1018 else if( jetType==C_JET ) _inputsHistogramsCJets[iInputsCollection][(*iTagNames).first]->fill(input);
1019 else _inputsHistogramsUDSJets[iInputsCollection][(*iTagNames).first]->fill(input);
1023 for (std::vector<std::string>::const_iterator iter = _ZoomedVarNames.begin() ; iter != _ZoomedVarNames.end(); iter++){
1024 if ((*iTagNames).first == *iter) {
1025 std::string zoomed_name = (*iTagNames).first +
" (zoomed)";
1026 if( jetType==B_JET ) _zoomedInputsHistogramsBJets[iInputsCollection][zoomed_name]->fill(input);
1027 else if( jetType==C_JET ) _zoomedInputsHistogramsCJets[iInputsCollection][zoomed_name]->fill(input);
1028 else _zoomedInputsHistogramsUDSJets[iInputsCollection][zoomed_name]->fill(input);
1042 void LCFIAIDAPlotProcessor::FillTagPlots( LCEvent* pEvent,
unsigned int jetNumber)
1044 int jetType=FindTrueJetType( pEvent, jetNumber );
1045 if( jetType==0 )
return;
1049 ReconstructedParticle* pJet;
1050 pJet=jetCollection.getElementAt(jetNumber);
1052 for (
unsigned int iTagCollection=0; iTagCollection < _FlavourTagCollectionNames.size(); ++iTagCollection)
1055 if( !tagCollection.is_valid() )
1057 std::cerr << _FlavourTagCollectionNames[iTagCollection] <<
"is not valid, not filling tag plots" << std::endl;
1058 std::cerr <<
"set MakeAdditionalPlots to false to remove this warning" << std::endl;
1063 lcio::LCFloatVec* pJetFlavourTags=tagCollection.getElementAt( jetNumber );
1064 if( !pJetFlavourTags )
1073 double bTag= (*pJetFlavourTags)[_IndexOfForEachTag[iTagCollection][
"BTag"]];
1074 double cTag= (*pJetFlavourTags)[_IndexOfForEachTag[iTagCollection][
"CTag"]];
1075 double cTagBBack= (*pJetFlavourTags)[_IndexOfForEachTag[iTagCollection][
"BCTag"]];
1076 unsigned int NumVertices = FindNumVertex(pEvent, jetNumber, iTagCollection);
1081 std::string nvname = _VertexCatNames[ (NumVertices>=N_VERTEX_CATEGORIES) ? (N_VERTEX_CATEGORIES) : (NumVertices)];
1083 if( jetType==B_JET ) {
1085 if( bTag<=1 && bTag>=0 )
1087 _pBJetBTag[iTagCollection][nvname]->fill( bTag );
1091 _pBJetBTag[iTagCollection][nvname]->fill( -0.005 );
1094 if( cTag<=1 && cTag>=0 )
1096 _pBJetCTag[iTagCollection][nvname]->fill( cTag );
1100 _pBJetCTag[iTagCollection][nvname]->fill( -0.005 );
1102 if( cTagBBack<=1 && cTagBBack>=0 )
1104 _pBJetBCTag[iTagCollection][nvname]->fill( cTagBBack );
1108 _pBJetBCTag[iTagCollection][nvname]->fill( -0.005 );
1111 }
else if( jetType==C_JET ) {
1113 if( bTag<=1 && bTag>=0 )
1115 _pCJetBTag[iTagCollection][nvname]->fill( bTag );
1119 _pCJetBTag[iTagCollection][nvname]->fill( -0.005 );
1122 if( cTag<=1 && cTag>=0 )
1124 _pCJetCTag[iTagCollection][nvname]->fill( cTag );
1128 _pCJetCTag[iTagCollection][nvname]->fill( -0.005 );
1131 if( cTagBBack<=1 && cTagBBack>=0 )
1133 _pCJetBCTag[iTagCollection][nvname]->fill( cTagBBack );
1137 _pCJetBCTag[iTagCollection][nvname]->fill( -0.005 );
1140 if( bTag<=1 && bTag>=0 )
1142 _pLightJetBTag[iTagCollection][nvname]->fill( bTag );
1146 _pLightJetBTag[iTagCollection][nvname]->fill( -0.005 );
1149 if( cTag<=1 && cTag>=0 )
1151 _pLightJetCTag[iTagCollection][nvname]->fill( cTag );
1155 _pLightJetCTag[iTagCollection][nvname]->fill( -0.005 );
1158 if( cTagBBack<=1 && cTagBBack>=0 )
1160 _pLightJetBCTag[iTagCollection][nvname]->fill( cTagBBack );
1164 _pLightJetBCTag[iTagCollection][nvname]->fill( -0.005 );
1172 void LCFIAIDAPlotProcessor::CreateAdditionalPlots()
1174 _pBJetCharge.resize( _FlavourTagCollectionNames.size() );
1175 _pCJetCharge.resize( _FlavourTagCollectionNames.size() );
1176 _pCDecayLengthAll.resize( _FlavourTagCollectionNames.size() );
1177 _pBDecayLengthAll.resize( _FlavourTagCollectionNames.size() );
1178 _pCDecayLengthTwoVertices.resize( _FlavourTagCollectionNames.size() );
1179 _pBDecayLengthTwoVertices.resize( _FlavourTagCollectionNames.size() );
1182 AIDA::IHistogramFactory* pHistogramFactory=marlin::AIDAProcessor::histogramFactory(
this );
1183 AIDA::ITree* pTree=marlin::AIDAProcessor::tree(
this );
1185 if (!pTree->cd(
"/" + name() +
"/ZVTOPPlots/")) {
1186 pTree->cd(
"/" + name() +
"/");
1187 pTree->mkdir(
"ZVTOPPlots/");
1188 pTree->cd(
"ZVTOPPlots/");
1193 _decayLengthBJet2D = pHistogramFactory->createHistogram2D(
"B jets: Reconstructed secondary decay length vs MC B decay length (mm)",50,0.,20.,50,0.,20.);
1194 _decayLengthCJet2D = pHistogramFactory->createHistogram2D(
"B jets: Reconstructed sec-ter decay length vs MC D decay length (mm)",50,0.,20.,50,0.,20.);
1195 _decayLengthBJetCloud2D = pHistogramFactory->createCloud2D(
"B jets: Reconstructed secondary decay length vs MC B decay length (mm) (cloud)");
1196 _decayLengthCJetCloud2D = pHistogramFactory->createCloud2D(
"B jets: Reconstructed sec-ter decay length vs MC D decay length (mm) (cloud)");
1198 _reconstructedSecondaryDecayLength = pHistogramFactory->createHistogram1D(
"All jets: Reconstructed secondary decay length (mm)",100,0.,20.);
1199 _reconstructedSecTerDecayLength = pHistogramFactory->createHistogram1D(
"All jets: Reconstructed secondary-tertiary decay length",100,0.,20.);
1201 _recoDecayLengthBJet = pHistogramFactory->createHistogram1D(
"B jets: Reconstructed secondary decay length (mm)",50,0.,20.);
1202 _recoDecayLengthCJet = pHistogramFactory->createHistogram1D(
"C jets: Reconstructed secondary decay length (mm)",50,0.,20.);
1203 _recoDecayLengthBCJet = pHistogramFactory->createHistogram1D(
"B jets: Reconstructed secondary-tertiary decay length (mm)",50,0.,20.);
1204 _recoDecayLengthLightJet = pHistogramFactory->createHistogram1D(
"Light jets: Reconstructed secondary decay length (mm)",50,0.,10.);
1206 _nVerticesBJet = pHistogramFactory->createHistogram1D(
"B jets: Number of reconstructed vertices",10,-0.5,9.5);
1207 _nVerticesCJet = pHistogramFactory->createHistogram1D(
"C jets: Number of reconstructed vertices",10,-0.5,9.5);
1208 _nVerticesLightJet = pHistogramFactory->createHistogram1D(
"Light jets: Number of reconstructed vertices",10,-0.5,9.5);
1210 _decayLengthBJetTrue = pHistogramFactory->createHistogram1D(
"B jets: MC secondary decay length (mm)",50,0.,20.);;
1211 _decayLengthBCJetTrue = pHistogramFactory->createHistogram1D(
"B jets: MC secondary-tertiary decay length (mm)",50,0.,20.);;
1212 _decayLengthBJetTrue = pHistogramFactory->createHistogram1D(
"C jets: MC secondary decay length (mm)",50,0.,20.);;
1216 if (!pTree->cd(
"/" + name() +
"/VertexChargePlots/")) {
1217 pTree->cd(
"/" + name() +
"/");
1218 pTree->mkdir(
"VertexChargePlots/");
1219 pTree->cd(
"VertexChargePlots/");
1222 _pBJetCharge2D = pHistogramFactory->createHistogram2D(
"B Jets: Reconstructed Vertex Charge vs True Jet Charge",7,-3.5,+3.5,7,-3.5,+3.5);
1223 _pCJetCharge2D = pHistogramFactory->createHistogram2D(
"C Jets: Reconstructed Vertex Charge vs True Jet Charge",7,-3.5,+3.5,7,-3.5,+3.5);
1225 _pBJetVertexCharge = pHistogramFactory->createHistogram1D(
"B Jets: Reconstructed Vertex Charge",9,-4.5,+4.5);
1226 _pCJetVertexCharge = pHistogramFactory->createHistogram1D(
"C Jets: Reconstructed Vertex Charge",9,-4.5,+4.5);
1228 _pCJetLeakageRate = pHistogramFactory->createHistogram1D(
"C Jets: Charged Leakage Rate (DON'T TRUST ERRORS)", N_JETANGLE_BINS,0.,1.);
1229 _pBJetLeakageRate = pHistogramFactory->createHistogram1D(
"B Jets: Charged Leakage Rate (DON'T TRUST ERRORS)", N_JETANGLE_BINS,0.,1.);
1233 for (
unsigned int iTagCollection=0; iTagCollection < _FlavourTagCollectionNames.size(); ++iTagCollection )
1236 if (!pTree->cd(
"/" + name() +
"/" +_FlavourTagCollectionNames[iTagCollection] +
"/")) {
1237 pTree->cd(
"/" + name() +
"/");
1238 if (!pTree->cd( _FlavourTagCollectionNames[iTagCollection])) {
1239 pTree->mkdir( _FlavourTagCollectionNames[iTagCollection] +
"/");
1240 pTree->cd( _FlavourTagCollectionNames[iTagCollection]+
"/");
1244 if (!pTree->cd(
"VertexChargePlots/")) {
1245 pTree->mkdir(
"VertexChargePlots/");
1246 pTree->cd(
"VertexChargePlots/");
1249 _pCDecayLengthAll[iTagCollection] = pHistogramFactory->createHistogram1D(
"True MC decay length of all reconstructed C-jets", 25, 0, 10.0 );
1250 _pBDecayLengthAll[iTagCollection] = pHistogramFactory->createHistogram1D(
"True MC decay length of all reconstructed B-jets", 25, 0, 10.5 );
1251 _pCDecayLengthTwoVertices[iTagCollection] = pHistogramFactory->createHistogram1D(
"True MC decay length of C-jets with more than one vertex", 25, 0, 10.0 );
1252 _pBDecayLengthTwoVertices[iTagCollection] = pHistogramFactory->createHistogram1D(
"True MC decay length of B-jets with more than one vertex", 25, 0, 10.5 );
1254 _pBJetCharge[iTagCollection] = pHistogramFactory->createHistogram2D(
"B Jets: Reconstructed Jet Charge vs True Jet Charge",7,-3.5,+3.5,7,-3.5,+3.5);
1255 _pCJetCharge[iTagCollection] = pHistogramFactory->createHistogram2D(
"C Jets: Reconstructed Jet Charge vs True Jet Charge",7,-3.5,+3.5,7,-3.5,+3.5);
1259 pTree->cd(
"/" + name());
1261 if (!pTree->cd(
"VertexPlots/")) {
1262 pTree->mkdir(
"VertexPlots/");
1263 pTree->cd(
"VertexPlots/");
1266 _pVertexDistanceFromIP = pHistogramFactory->createHistogram1D(
"Reconstructed Vertex distance from IP",100, 0., 10.);
1267 _pVertexPositionX = pHistogramFactory->createHistogram1D(
"Non-primary vertex: x-position", 100, -10., 10.) ;
1268 _pVertexPositionY = pHistogramFactory->createHistogram1D(
"Non-primary vertex: y-position", 100, -10., 10.);
1269 _pVertexPositionZ = pHistogramFactory->createHistogram1D(
"Non-primary vertex: z-position", 100, -10., 10.);
1271 _pPrimaryVertexPullX = pHistogramFactory->createHistogram1D(
"Non-primary vertex: x-pull", 100, -10., 10.);
1272 _pPrimaryVertexPullY = pHistogramFactory->createHistogram1D(
"Non-primary vertex: y-pull", 100, -10., 10.);
1273 _pPrimaryVertexPullZ = pHistogramFactory->createHistogram1D(
"Non-primary vertex: z-pull", 100, -10., 10.);
1274 _pPrimaryVertexPositionX = pHistogramFactory->createHistogram1D(
"Primary vertex: x-position", 100, -10., 10.);
1275 _pPrimaryVertexPositionY = pHistogramFactory->createHistogram1D(
"Primary vertex: y-position", 100, -10., 10.);
1276 _pPrimaryVertexPositionZ = pHistogramFactory->createHistogram1D(
"Primary vertex: z-position", 100, -10., 10.);
1282 void LCFIAIDAPlotProcessor::CreateTagPlots()
1284 AIDA::IHistogramFactory* pHistogramFactory=marlin::AIDAProcessor::histogramFactory(
this );
1285 AIDA::ITree* pTree=marlin::AIDAProcessor::tree(
this );
1287 _pBJetBTag.resize( _FlavourTagCollectionNames.size() );
1288 _pBJetCTag.resize( _FlavourTagCollectionNames.size() );
1289 _pBJetBCTag.resize( _FlavourTagCollectionNames.size() );
1291 _pCJetCTag.resize( _FlavourTagCollectionNames.size() );
1292 _pCJetBTag.resize( _FlavourTagCollectionNames.size() );
1293 _pCJetBCTag.resize( _FlavourTagCollectionNames.size() );
1295 _pLightJetBTag.resize( _FlavourTagCollectionNames.size() );
1296 _pLightJetCTag.resize( _FlavourTagCollectionNames.size() );
1297 _pLightJetBCTag.resize( _FlavourTagCollectionNames.size() );
1299 _pBTagBackgroundValues.resize( _FlavourTagCollectionNames.size() );
1300 _pCTagBackgroundValues.resize( _FlavourTagCollectionNames.size() );
1301 _pBCTagBackgroundValues.resize( _FlavourTagCollectionNames.size() );
1304 for (
unsigned int iTagCollection=0; iTagCollection < _FlavourTagCollectionNames.size(); ++iTagCollection )
1306 for (
unsigned int iVertexCat=0; iVertexCat < N_VERTEX_CATEGORIES+1; ++iVertexCat ){
1307 _pLightJetBTag[iTagCollection][_VertexCatNames[iVertexCat]]=0;
1308 _pLightJetCTag[iTagCollection][_VertexCatNames[iVertexCat]]=0;
1309 _pLightJetBCTag[iTagCollection][_VertexCatNames[iVertexCat]]=0;
1311 _pBJetBTag[iTagCollection][_VertexCatNames[iVertexCat]]=0;
1312 _pBJetCTag[iTagCollection][_VertexCatNames[iVertexCat]]=0;
1313 _pBJetBCTag[iTagCollection][_VertexCatNames[iVertexCat]]=0;
1315 _pCJetBTag[iTagCollection][_VertexCatNames[iVertexCat]]=0;
1316 _pCJetCTag[iTagCollection][_VertexCatNames[iVertexCat]]=0;
1317 _pCJetBCTag[iTagCollection][_VertexCatNames[iVertexCat]]=0;
1319 _pBTagBackgroundValues[iTagCollection][_VertexCatNames[iVertexCat]]=0;
1320 _pCTagBackgroundValues[iTagCollection][_VertexCatNames[iVertexCat]]=0;
1321 _pBCTagBackgroundValues[iTagCollection][_VertexCatNames[iVertexCat]]=0;
1326 for (
unsigned int iTagCollection=0; iTagCollection < _FlavourTagCollectionNames.size(); ++iTagCollection )
1329 if (!pTree->cd(
"/" + name() +
"/" +_FlavourTagCollectionNames[iTagCollection] +
"/")) {
1330 pTree->cd(
"/" + name() +
"/");
1331 if (!pTree->cd( _FlavourTagCollectionNames[iTagCollection])) {
1332 pTree->mkdir( _FlavourTagCollectionNames[iTagCollection] +
"/");
1333 pTree->cd( _FlavourTagCollectionNames[iTagCollection]+
"/");
1337 for (
unsigned int iVertexCat=0; iVertexCat <= N_VERTEX_CATEGORIES; ++iVertexCat ){
1339 std::string nvname = _VertexCatNames[iVertexCat];
1342 pTree->cd(
"/" + name() +
"/" +_FlavourTagCollectionNames[iTagCollection] +
"/");
1344 if (!pTree->cd( _NumVertexCatDir[iVertexCat]+
"/")) {
1345 pTree->mkdir( _NumVertexCatDir[iVertexCat]+
"/");
1346 pTree->cd( _NumVertexCatDir[iVertexCat]+
"/");
1349 _pLightJetBTag[iTagCollection][nvname] = pHistogramFactory->createHistogram1D(
"Numbers of light jets by B-tag NN value. ("+ nvname +
")",_numberOfPoints , 0, 1.0 );
1350 _pLightJetCTag[iTagCollection][nvname] = pHistogramFactory->createHistogram1D(
"Numbers of light jets by C-tag NN value. ("+ nvname +
")", _numberOfPoints, 0, 1.0 );
1351 _pLightJetBCTag[iTagCollection][nvname] = pHistogramFactory->createHistogram1D(
"Numbers of light jets by BC-tag NN value. ("+ nvname +
")", _numberOfPoints, 0, 1.0 );
1352 _pBJetBTag[iTagCollection][nvname] = pHistogramFactory->createHistogram1D(
"Numbers of B jets by B-tag NN value. ("+ nvname +
")", _numberOfPoints, 0, 1.0 );
1353 _pBJetCTag[iTagCollection][nvname] = pHistogramFactory->createHistogram1D(
"Numbers of B jets by C-tag NN value. ("+ nvname +
")", _numberOfPoints, 0, 1.0 );
1354 _pBJetBCTag[iTagCollection][nvname] = pHistogramFactory->createHistogram1D(
"Numbers of B jets by BC-tag NN value. ("+ nvname +
")", _numberOfPoints, 0, 1.0 );
1355 _pCJetBTag[iTagCollection][nvname] = pHistogramFactory->createHistogram1D(
"Numbers of C jets by B-tag NN value. ("+ nvname +
")", _numberOfPoints, 0, 1.0 );
1356 _pCJetCTag[iTagCollection][nvname] = pHistogramFactory->createHistogram1D(
"Numbers of C jets by C-tag NN value. ("+ nvname +
")", _numberOfPoints, 0, 1.0 );
1357 _pCJetBCTag[iTagCollection][nvname] = pHistogramFactory->createHistogram1D(
"Numbers of C jets by BC-tag NN value. ("+ nvname +
")", _numberOfPoints, 0, 1.0 );
1363 void LCFIAIDAPlotProcessor::FillVertexChargePlots(LCEvent* pEvent,
unsigned int jetNumber)
1365 int jetType=FindTrueJetType( pEvent, jetNumber );
1366 if( jetType==0 )
return;
1370 ReconstructedParticle* pJet;
1371 pJet=jetCollection.getElementAt(jetNumber);
1373 for (
unsigned int iTagCollection=0; iTagCollection < _FlavourTagCollectionNames.size(); ++iTagCollection)
1376 if( !tagCollection.is_valid() )
1378 std::cerr << _FlavourTagCollectionNames[iTagCollection] <<
"is not valid, not filling vertex charge plots" << std::endl;
1379 std::cerr <<
"set MakeAdditionalPlots to false to remove this warning" << std::endl;
1383 lcio::LCFloatVec* pJetFlavourTags=tagCollection.getElementAt( jetNumber );
1384 if( !pJetFlavourTags )
1389 const double* jetMomentum = pJet->getMomentum();
1390 double cosTheta = jetMomentum[2] / sqrt(pow(jetMomentum[0],2)+pow(jetMomentum[1],2)+pow(jetMomentum[2],2));
1393 if (iTagCollection == _myVertexChargeTagCollection) {
1396 double bTag= (*pJetFlavourTags)[_IndexOfForEachTag[iTagCollection][
"BTag"]];
1397 double cTag= (*pJetFlavourTags)[_IndexOfForEachTag[iTagCollection][
"CTag"]];
1399 unsigned int NumVertices = FindNumVertex(pEvent, jetNumber, iTagCollection);
1400 int CQVtx = FindCQVtx(pEvent, jetNumber);
1401 int BQVtx = FindBQVtx(pEvent, jetNumber);
1402 int trueJetCharge = int(FindTrueJetHadronCharge(pEvent,jetNumber));
1404 std::string nvname = _VertexCatNames[ (NumVertices>=N_VERTEX_CATEGORIES) ? (N_VERTEX_CATEGORIES) : (NumVertices)];
1407 if( jetType==C_JET && cTag > _CTagNNCut) {
1409 int bin = _pCJetLeakageRate->coordToIndex(fabs(cosTheta));
1411 if (trueJetCharge==+2) _cJet_truePlus2++;
1412 if (trueJetCharge==+1) _cJet_truePlus++;
1413 if (trueJetCharge==0) _cJet_trueNeut++;
1414 if (trueJetCharge==-1) _cJet_trueMinus++;
1415 if (trueJetCharge==-2) _cJet_trueMinus2++;
1417 if (trueJetCharge==+2){
1418 if(CQVtx>0) _cJet_truePlus2_recoPlus++;
1419 if(CQVtx==0) _cJet_truePlus2_recoNeut++;
1420 if(CQVtx<0) _cJet_truePlus2_recoMinus++;
1422 if (trueJetCharge==+1){
1423 if(CQVtx>0) _cJet_truePlus_recoPlus++;
1424 if(CQVtx==0) _cJet_truePlus_recoNeut++;
1425 if(CQVtx<0) _cJet_truePlus_recoMinus++;
1427 if (trueJetCharge==0) {
1428 if(CQVtx>0) _cJet_trueNeut_recoPlus++;
1429 if(CQVtx==0) _cJet_trueNeut_recoNeut++;
1430 if(CQVtx<0) _cJet_trueNeut_recoMinus++;
1432 if (trueJetCharge==-1) {
1433 if(CQVtx>0) _cJet_trueMinus_recoPlus++;
1434 if(CQVtx==0) _cJet_trueMinus_recoNeut++;
1435 if(CQVtx<0) _cJet_trueMinus_recoMinus++;
1437 if (trueJetCharge==-2) {
1438 if(CQVtx>0) _cJet_trueMinus2_recoPlus++;
1439 if(CQVtx==0) _cJet_trueMinus2_recoNeut++;
1440 if(CQVtx<0) _cJet_trueMinus2_recoMinus++;
1443 _pCJetVertexCharge->fill(CQVtx);
1444 _pCJetCharge2D->fill(trueJetCharge,CQVtx);
1446 if (trueJetCharge==+2) _cJet_truePlus2_angle[bin]++;
1447 if (trueJetCharge==+1) _cJet_truePlus_angle[bin]++;
1448 if (trueJetCharge==0) _cJet_trueNeut_angle[bin]++;
1449 if (trueJetCharge==-1) _cJet_trueMinus_angle[bin]++;
1450 if (trueJetCharge==-2) _cJet_trueMinus2_angle[bin]++;
1452 if (trueJetCharge==+2){
1453 if(CQVtx>0) _cJet_truePlus2_recoPlus_angle[bin]++;
1454 if(CQVtx==0) _cJet_truePlus2_recoNeut_angle[bin]++;
1455 if(CQVtx<0) _cJet_truePlus2_recoMinus_angle[bin]++;
1457 if (trueJetCharge==+1){
1458 if(CQVtx>0) _cJet_truePlus_recoPlus_angle[bin]++;
1459 if(CQVtx==0) _cJet_truePlus_recoNeut_angle[bin]++;
1460 if(CQVtx<0) _cJet_truePlus_recoMinus_angle[bin]++;
1462 if (trueJetCharge==0) {
1463 if(CQVtx>0) _cJet_trueNeut_recoPlus_angle[bin]++;
1464 if(CQVtx==0) _cJet_trueNeut_recoNeut_angle[bin]++;
1465 if(CQVtx<0) _cJet_trueNeut_recoMinus_angle[bin]++;
1467 if (trueJetCharge==-1) {
1468 if(CQVtx>0) _cJet_trueMinus_recoPlus_angle[bin]++;
1469 if(CQVtx==0) _cJet_trueMinus_recoNeut_angle[bin]++;
1470 if(CQVtx<0) _cJet_trueMinus_recoMinus_angle[bin]++;
1472 if (trueJetCharge==-2) {
1473 if(CQVtx>0) _cJet_trueMinus2_recoPlus_angle[bin]++;
1474 if(CQVtx==0) _cJet_trueMinus2_recoNeut_angle[bin]++;
1475 if(CQVtx<0) _cJet_trueMinus2_recoMinus_angle[bin]++;
1478 }
else if ( jetType==B_JET && bTag > _BTagNNCut) {
1480 int bin = _pBJetLeakageRate->coordToIndex(fabs(cosTheta));
1482 if (trueJetCharge==+2) _bJet_truePlus2++;
1483 if (trueJetCharge==+1) _bJet_truePlus++;
1484 if (trueJetCharge==0) _bJet_trueNeut++;
1485 if (trueJetCharge==-1) _bJet_trueMinus++;
1486 if (trueJetCharge==-2) _bJet_trueMinus2++;
1488 if (trueJetCharge==+2){
1489 if(BQVtx>0) _bJet_truePlus2_recoPlus++;
1490 if(BQVtx==0) _bJet_truePlus2_recoNeut++;
1491 if(BQVtx<0) _bJet_truePlus2_recoMinus++;
1493 if (trueJetCharge==+1){
1494 if(BQVtx>0) _bJet_truePlus_recoPlus++;
1495 if(BQVtx==0) _bJet_truePlus_recoNeut++;
1496 if(BQVtx<0) _bJet_truePlus_recoMinus++;
1498 if (trueJetCharge==0) {
1499 if(BQVtx>0) _bJet_trueNeut_recoPlus++;
1500 if(BQVtx==0) _bJet_trueNeut_recoNeut++;
1501 if(BQVtx<0) _bJet_trueNeut_recoMinus++;
1503 if (trueJetCharge==-1) {
1504 if(BQVtx>0) _bJet_trueMinus_recoPlus++;
1505 if(BQVtx==0) _bJet_trueMinus_recoNeut++;
1506 if(BQVtx<0) _bJet_trueMinus_recoMinus++;
1508 if (trueJetCharge==-2) {
1509 if(BQVtx>0) _bJet_trueMinus2_recoPlus++;
1510 if(BQVtx==0) _bJet_trueMinus2_recoNeut++;
1511 if(BQVtx<0) _bJet_trueMinus2_recoMinus++;
1515 if (trueJetCharge==+2) _bJet_truePlus2_angle[bin]++;
1516 if (trueJetCharge==+1) _bJet_truePlus_angle[bin]++;
1517 if (trueJetCharge==0) _bJet_trueNeut_angle[bin]++;
1518 if (trueJetCharge==-1) _bJet_trueMinus_angle[bin]++;
1519 if (trueJetCharge==-2) _bJet_trueMinus2_angle[bin]++;
1521 if (trueJetCharge==+2){
1522 if(BQVtx>0) _bJet_truePlus2_recoPlus_angle[bin]++;
1523 if(BQVtx==0) _bJet_truePlus2_recoNeut_angle[bin]++;
1524 if(BQVtx<0) _bJet_truePlus2_recoMinus_angle[bin]++;
1526 if (trueJetCharge==+1){
1527 if(BQVtx>0) _bJet_truePlus_recoPlus_angle[bin]++;
1528 if(BQVtx==0) _bJet_truePlus_recoNeut_angle[bin]++;
1529 if(BQVtx<0) _bJet_truePlus_recoMinus_angle[bin]++;
1531 if (trueJetCharge==0) {
1532 if(BQVtx>0) _bJet_trueNeut_recoPlus_angle[bin]++;
1533 if(BQVtx==0) _bJet_trueNeut_recoNeut_angle[bin]++;
1534 if(BQVtx<0) _bJet_trueNeut_recoMinus_angle[bin]++;
1536 if (trueJetCharge==-1) {
1537 if(BQVtx>0) _bJet_trueMinus_recoPlus_angle[bin]++;
1538 if(BQVtx==0) _bJet_trueMinus_recoNeut_angle[bin]++;
1539 if(BQVtx<0) _bJet_trueMinus_recoMinus_angle[bin]++;
1541 if (trueJetCharge==-2) {
1542 if(BQVtx>0) _bJet_trueMinus2_recoPlus_angle[bin]++;
1543 if(BQVtx==0) _bJet_trueMinus2_recoNeut_angle[bin]++;
1544 if(BQVtx<0) _bJet_trueMinus2_recoMinus_angle[bin]++;
1547 _pBJetVertexCharge->fill(BQVtx);
1548 _pBJetCharge2D->fill(trueJetCharge,BQVtx);
1562 if( !trueJetFlavourCollection.is_valid() )
1564 std::cerr <<
" In " << __FILE__ <<
"(" << __LINE__ <<
"): Collection " << _TrueJetFlavourColName <<
" is not valid " << std::endl;
1568 lcio::LCFloatVec* pJetFlavour = trueJetFlavourCollection.getElementAt( jetNumber );
1572 std::cerr << __FILE__ <<
"(" << __LINE__ <<
"): Unable to get the LCFloatVec for jet " << jetNumber <<
" from the collection " << _TrueJetFlavourColName
1573 <<
" for event " << pEvent->getEventNumber() <<
" in run " << pEvent->getRunNumber() <<
"." << std::endl;
1577 pdgCode = (*pJetFlavour)[_FlavourIndex[
"TrueJetPDGFlavour"]];
1580 return int(pdgCode+0.01);
1590 if( !trueJetFlavourCollection.is_valid() )
1592 std::cerr <<
" In " << __FILE__ <<
"(" << __LINE__ <<
"): Collection " << _TrueJetFlavourColName <<
" is not valid " << std::endl;
1596 lcio::LCFloatVec* pJetFlavour = trueJetFlavourCollection.getElementAt( jetNumber );
1599 std::cerr << __FILE__ <<
"(" << __LINE__ <<
"): Unable to get the LCFloatVec for jet " << jetNumber <<
" from the collection " << _TrueJetFlavourColName
1600 <<
" for event " << pEvent->getEventNumber() <<
" in run " << pEvent->getRunNumber() <<
"." << std::endl;
1604 return (*pJetFlavour)[_FlavourIndex[
"TruePartonCharge"]];
1612 float jetFlavour=0.;
1615 if( !trueJetFlavourCollection.is_valid() )
1617 std::cerr <<
" In " << __FILE__ <<
"(" << __LINE__ <<
"): Collection " << _TrueJetFlavourColName <<
" is not valid " << std::endl;
1621 lcio::LCFloatVec* pJetFlavour = trueJetFlavourCollection.getElementAt( jetNumber );
1624 std::cerr << __FILE__ <<
"(" << __LINE__ <<
"): Unable to get the LCFloatVec for jet " << jetNumber <<
" from the collection " << _TrueJetFlavourColName
1625 <<
" for event " << pEvent->getEventNumber() <<
" in run " << pEvent->getRunNumber() <<
"." << std::endl;
1629 jetFlavour=(*pJetFlavour)[_FlavourIndex[
"TrueJetFlavour"]];
1631 return int(jetFlavour+0.001);
1637 if( !trueJetFlavourCollection.is_valid() )
1639 std::cerr <<
" In " << __FILE__ <<
"(" << __LINE__ <<
"): Collection " << _TrueJetFlavourColName <<
" is not valid " << std::endl;
1643 lcio::LCFloatVec* pJetFlavour = trueJetFlavourCollection.getElementAt( jetNumber );
1646 std::cerr << __FILE__ <<
"(" << __LINE__ <<
"): Unable to get the LCFloatVec for jet " << jetNumber <<
" from the collection " << _TrueJetFlavourColName
1647 <<
" for event " << pEvent->getEventNumber() <<
" in run " << pEvent->getRunNumber() <<
"." << std::endl;
1651 return (*pJetFlavour)[_FlavourIndex[
"TrueJetHadronCharge"]];
1661 return FindTrueJetFlavour(pEvent, jetNumber);
1670 if( !inputsCollection.is_valid() ) {
1672 std::cerr << __FILE__ <<
"(" << __LINE__ <<
"): Unable to get the LCFloatVec from the collection " << _FlavourTagInputsCollectionNames[iInputsCollection]
1673 <<
" for event " << pEvent->getEventNumber() <<
" in run " << pEvent->getRunNumber() <<
"." << std::endl;
1680 lcio::LCFloatVec* pInputs=inputsCollection.getElementAt( jetNumber );
1684 std::cerr << __FILE__ <<
"(" << __LINE__ <<
"): Unable to get the LCFloatVec for jet " << jetNumber <<
" from the collection " << _FlavourTagInputsCollectionNames[iInputsCollection]
1685 <<
" for event " << pEvent->getEventNumber() <<
" in run " << pEvent->getRunNumber() <<
"." << std::endl;
1689 return int((*pInputs)[_InputsIndex[iInputsCollection][
"NumVertices"]]);
1701 if( !inputsCollection.is_valid() ) {
1703 std::cerr <<
"In " << __FILE__ <<
"(" << __LINE__ <<
"): Cannot find collection " << _BVertexChargeCollection <<
" for event " << pEvent->getEventNumber() <<
" in run " << pEvent->getRunNumber() <<
" BQVtx will be invalid" << std::endl;
1709 lcio::LCFloatVec* pBVtxChargeVector =inputsCollection.getElementAt( jetNumber );
1712 if( pBVtxChargeVector && pBVtxChargeVector->size() == 1) {
1714 bqvtx = pBVtxChargeVector->back();
1718 std::cerr <<
"In " << __FILE__ <<
"(" << __LINE__ <<
"): Cannot find collection element in " << _BVertexChargeCollection <<
" for event " << pEvent->getEventNumber() <<
" in run " << pEvent->getRunNumber() <<
" corresponding to jet number: " << jetNumber <<
" BQVtx will be invalid" << std::endl;
1735 if( !inputsCollection.is_valid() ) {
1737 std::cerr <<
"In " << __FILE__ <<
"(" << __LINE__ <<
"): Cannot find collection " << _CVertexChargeCollection <<
" for event " << pEvent->getEventNumber() <<
" in run " << pEvent->getRunNumber() <<
" CQVtx will be invalid" << std::endl;
1743 lcio::LCFloatVec* pCVtxChargeVector =inputsCollection.getElementAt( jetNumber );
1746 if( pCVtxChargeVector && pCVtxChargeVector->size() == 1) {
1748 cqvtx = pCVtxChargeVector->back();
1752 std::cerr <<
"In " << __FILE__ <<
"(" << __LINE__ <<
"): Cannot find collection element in " << _CVertexChargeCollection <<
" for event " << pEvent->getEventNumber() <<
" in run " << pEvent->getRunNumber() <<
" corresponding to jet number: " << jetNumber <<
" CQVtx will be invalid" << std::endl;
1767 const int numberOfBins=pNN->axis().bins();
1769 double integral=pNN->binHeight(AIDA::IAxis::OVERFLOW_BIN);
1770 pIntegral->fill(pNN->axis().binLowerEdge(AIDA::IAxis::OVERFLOW_BIN)+pNN->axis().binWidth(numberOfBins-1)/2.,integral);
1772 for(
int binNumber=numberOfBins-1; binNumber>=0; --binNumber )
1774 integral+= pNN->binHeight( binNumber );
1775 pIntegral->fill( pNN->axis().binLowerEdge(binNumber)+pNN->axis().binWidth(binNumber)/2.,integral);
1778 integral+= pNN->binHeight(AIDA::IAxis::UNDERFLOW_BIN);
1779 pIntegral->fill(pNN->axis().binUpperEdge(AIDA::IAxis::UNDERFLOW_BIN)-pNN->axis().binWidth(0)/2.,integral);
1786 void LCFIAIDAPlotProcessor::CalculateTagEfficiency(
const AIDA::IHistogram1D* pSignal, std::vector<double>& tagEfficiency, std::vector<double>& tagEfficiencyError)
1789 double totalSignal=pSignal->sumBinHeights();
1790 double signalPassedCut=0;
1791 const int numberOfBins=pSignal->axis().bins();
1793 tagEfficiency.clear();
1794 tagEfficiencyError.clear();
1797 for(
int binNumber=numberOfBins-1; binNumber>=0; --binNumber)
1799 signalPassedCut+=pSignal->binHeight( binNumber );
1801 double efficiency = signalPassedCut/totalSignal;
1803 double efficiencyError = efficiency * (1. - efficiency) / totalSignal;
1804 if (efficiencyError>0) efficiencyError = sqrt(efficiencyError);
1806 tagEfficiency.push_back(efficiency);
1807 tagEfficiencyError.push_back(efficiencyError);
1813 void LCFIAIDAPlotProcessor::CalculateTagPurity(
const AIDA::IHistogram1D* pSignal,
const AIDA::IHistogram1D* pBackground, std::vector<double>& tagPurity, std::vector<double>& tagPurityError)
1815 const int numberOfBins=pSignal->axis().bins();
1817 double signalPassedCut=0;
1818 double backgroundPassedCut=0;
1821 tagPurityError.clear();
1824 for (
int binNumber = numberOfBins-1; binNumber >= 0 ; --binNumber)
1827 signalPassedCut+=pSignal->binHeight( binNumber );
1828 backgroundPassedCut+=pBackground->binHeight( binNumber );
1830 double purity = ((signalPassedCut+backgroundPassedCut)>0) ? signalPassedCut/(signalPassedCut+backgroundPassedCut) : 1.;
1831 double purityError = ((signalPassedCut+backgroundPassedCut)>0) ? purity * (1. - purity) / (signalPassedCut+backgroundPassedCut) : 0.;
1832 if (purityError>0) purityError = sqrt(purityError);
1833 else purityError = -sqrt(-purityError);
1835 tagPurity.push_back(purity);
1836 tagPurityError.push_back(purityError);
1849 double totalSignal=pSignal->sumBinHeights();
1850 double signalPassedCut=0;
1852 const int numberOfBins=pSignal->axis().bins();
1855 for(
int binNumber=numberOfBins-1; binNumber>=0; --binNumber, iPoint++ )
1857 signalPassedCut+=pSignal->binHeight( binNumber );
1859 double efficiency = signalPassedCut/totalSignal;
1861 double efficiencyError = efficiency * (1. - efficiency) / totalSignal;
1862 if (efficiencyError>0) efficiencyError = sqrt(efficiencyError);
1865 pDataPointSet->addPoint();
1866 pDataPointSet->point(iPoint)->coordinate(0)->setValue(pSignal->axis().binLowerEdge(binNumber)+pSignal->axis().binWidth(binNumber)/2.);
1867 pDataPointSet->point(iPoint)->coordinate(1)->setValue( efficiency );
1868 pDataPointSet->point(iPoint)->coordinate(1)->setErrorPlus(efficiencyError);
1869 pDataPointSet->point(iPoint)->coordinate(1)->setErrorMinus(efficiencyError);
1872 return pDataPointSet;
1881 const int numberOfBins=pAllEvents->axis().bins();
1882 if (pPassEvents->axis().bins() != numberOfBins)
return pDataPointSet;
1886 for(
int binNumber=numberOfBins-1; binNumber>=0; --binNumber, iPoint++ )
1888 double all = pAllEvents->binHeight(binNumber);
1889 double pass = pPassEvents->binHeight(binNumber);
1891 double efficiency = (all>0) ? pass/all : 0.;
1893 double efficiencyError = (all>0) ? efficiency * (1. - efficiency) / all : 0.;
1894 if (efficiencyError>0.) efficiencyError = sqrt(efficiencyError);
1896 pDataPointSet->addPoint();
1897 pDataPointSet->point(iPoint)->coordinate(0)->setValue(pAllEvents->axis().binLowerEdge(binNumber)+pAllEvents->axis().binWidth(binNumber)/2.);
1898 pDataPointSet->point(iPoint)->coordinate(1)->setValue( efficiency );
1899 pDataPointSet->point(iPoint)->coordinate(1)->setErrorPlus(efficiencyError);
1900 pDataPointSet->point(iPoint)->coordinate(1)->setErrorMinus(efficiencyError);
1903 return pDataPointSet;
1912 const int numberOfBins=pNN->axis().bins();
1916 for (
int binNumber = 0; binNumber < numberOfBins ; binNumber++ )
1919 integral+= pNN->binHeight( binNumber );
1920 pDataPointSet->addPoint();
1921 pDataPointSet->point(binNumber)->coordinate(0)->setValue(pNN->axis().binLowerEdge(binNumber)+pNN->axis().binWidth(binNumber)/2.);
1922 pDataPointSet->point(binNumber)->coordinate(1)->setValue(integral);
1923 pDataPointSet->point(binNumber)->coordinate(1)->setErrorPlus(sqrt(integral));
1924 pDataPointSet->point(binNumber)->coordinate(1)->setErrorMinus(sqrt(integral));
1926 return pDataPointSet;
1931 const int numberOfBins=pSignal->axis().bins();
1934 double signalPassedCut=0;
1935 double backgroundPassedCut=0;
1937 for (
int binNumber = numberOfBins-1; binNumber >= 0 ; --binNumber, iPoint++ )
1940 signalPassedCut+=pSignal->binHeight( binNumber );
1941 backgroundPassedCut+=pBackground->binHeight( binNumber );
1943 double purity = signalPassedCut/(signalPassedCut+backgroundPassedCut);
1944 double purityError = purity * (1. - purity) / (signalPassedCut+backgroundPassedCut);
1945 if (purityError>0) purityError = sqrt(purityError);
1948 pDataPointSet->addPoint();
1949 pDataPointSet->point(iPoint)->coordinate(0)->setValue(pSignal->axis().binLowerEdge(binNumber)+pSignal->axis().binWidth(binNumber)/2.);
1950 pDataPointSet->point(iPoint)->coordinate(1)->setValue(purity);
1951 pDataPointSet->point(iPoint)->coordinate(1)->setErrorPlus(purityError);
1952 pDataPointSet->point(iPoint)->coordinate(1)->setErrorMinus(purityError);
1957 return pDataPointSet;
1963 double totalBackground = pBackground->sumBinHeights();
1964 double backgroundPassedCut=0;
1966 const int numberOfBins=pBackground->axis().bins();
1969 for(
int binNumber=numberOfBins-1; binNumber>=0; --binNumber , iPoint++ )
1972 backgroundPassedCut+=pBackground->binHeight( binNumber );
1974 double leakageRate = backgroundPassedCut/totalBackground;
1976 double leakageRateError = leakageRate * (1. - leakageRate) / totalBackground;
1977 if (leakageRateError>0) leakageRateError = sqrt(leakageRateError);
1980 pDataPointSet->addPoint();
1981 pDataPointSet->point(iPoint)->coordinate(0)->setValue(pBackground->axis().binLowerEdge(binNumber)+pBackground->axis().binWidth(binNumber)/2.);
1982 pDataPointSet->point(iPoint)->coordinate(1)->setValue(leakageRate);
1983 pDataPointSet->point(iPoint)->coordinate(1)->setErrorPlus(leakageRateError);
1984 pDataPointSet->point(iPoint)->coordinate(1)->setErrorMinus(leakageRateError);
1987 return pDataPointSet;
1992 AIDA::IDataPointSet*
LCFIAIDAPlotProcessor::CreateXYPlot(
const AIDA::IDataPointSet* pDataPointSet0,
const AIDA::IDataPointSet* pDataPointSet1, AIDA::IDataPointSet* xyPointSet,
const int dim0,
const int dim1 )
1996 if (pDataPointSet0->size() == pDataPointSet1->size()) {
1998 for (
int iPoint = 0 ; iPoint != pDataPointSet1->size(); iPoint++)
2000 xyPointSet->addPoint();
2001 xyPointSet->point(iPoint)->coordinate(0)->setValue(pDataPointSet0->point(iPoint)->coordinate(dim0)->value());
2002 xyPointSet->point(iPoint)->coordinate(1)->setValue(pDataPointSet1->point(iPoint)->coordinate(dim1)->value());
2003 xyPointSet->point(iPoint)->coordinate(0)->setErrorPlus(pDataPointSet0->point(iPoint)->coordinate(dim0)->errorPlus());
2004 xyPointSet->point(iPoint)->coordinate(1)->setErrorPlus(pDataPointSet1->point(iPoint)->coordinate(dim1)->errorPlus());
2005 xyPointSet->point(iPoint)->coordinate(0)->setErrorMinus(pDataPointSet0->point(iPoint)->coordinate(dim0)->errorMinus());
2006 xyPointSet->point(iPoint)->coordinate(1)->setErrorMinus(pDataPointSet1->point(iPoint)->coordinate(dim1)->errorMinus());
2017 for (
int j = 0 ; j < N_JETANGLE_BINS ; j++) {
2019 double c_numerator = _cJet_truePlus2_recoPlus_angle[j]+_cJet_truePlus_recoPlus_angle[j]
2020 +_cJet_trueMinus2_recoMinus_angle[j]+_cJet_trueMinus_recoMinus_angle[j];
2021 double c_domininator = _cJet_truePlus2_angle[j] +_cJet_truePlus_angle[j]
2022 +_cJet_trueMinus2_angle[j] +_cJet_trueMinus_angle[j];
2024 double b_numerator = _bJet_truePlus2_recoPlus_angle[j]+_bJet_truePlus_recoPlus_angle[j]
2025 +_bJet_trueMinus2_recoMinus_angle[j]+_bJet_trueMinus_recoMinus_angle[j];
2026 double b_domininator = _bJet_truePlus2_angle[j] +_bJet_truePlus_angle[j]
2027 +_bJet_trueMinus2_angle[j] +_bJet_trueMinus_angle[j];
2029 double b_leakage = 1. - b_numerator/b_domininator;
2030 double c_leakage = 1. - c_numerator/c_domininator;
2032 double b_leakage_error = sqrt( b_leakage * (1. - b_leakage) / b_domininator);
2033 double c_leakage_error = sqrt( c_leakage * (1. - c_leakage) / c_domininator);
2036 pBJetVtxChargeDPS->addPoint();
2037 pBJetVtxChargeDPS->point(j)->coordinate(0)->setValue(_pBJetLeakageRate->axis().binLowerEdge(j)+_pBJetLeakageRate->axis().binWidth(j)/2.);
2038 pBJetVtxChargeDPS->point(j)->coordinate(1)->setValue(b_leakage);
2039 pBJetVtxChargeDPS->point(j)->coordinate(1)->setErrorPlus(b_leakage_error);
2040 pBJetVtxChargeDPS->point(j)->coordinate(1)->setErrorMinus(b_leakage_error);
2043 pCJetVtxChargeDPS->addPoint();
2044 pCJetVtxChargeDPS->point(j)->coordinate(0)->setValue(_pCJetLeakageRate->axis().binLowerEdge(j)+_pCJetLeakageRate->axis().binWidth(j)/2.);
2045 pCJetVtxChargeDPS->point(j)->coordinate(1)->setValue(c_leakage);
2046 pCJetVtxChargeDPS->point(j)->coordinate(1)->setErrorPlus(c_leakage_error);
2047 pCJetVtxChargeDPS->point(j)->coordinate(1)->setErrorMinus(c_leakage_error);
2049 _pCJetLeakageRate->fill(_pCJetLeakageRate->axis().binLowerEdge(j)+_pCJetLeakageRate->axis().binWidth(j)/2.,c_leakage);
2050 _pBJetLeakageRate->fill(_pBJetLeakageRate->axis().binLowerEdge(j)+_pBJetLeakageRate->axis().binWidth(j)/2.,b_leakage);
2060 for (
int j = 0 ; j < N_JETANGLE_BINS ; j++) {
2062 double c_numerator = _cJet_truePlus2_recoPlus_angle[j]+_cJet_truePlus_recoPlus_angle[j]
2063 +_cJet_trueMinus2_recoMinus_angle[j]+_cJet_trueMinus_recoMinus_angle[j];
2064 double c_domininator = _cJet_truePlus2_angle[j] +_cJet_truePlus_angle[j]
2065 +_cJet_trueMinus2_angle[j] +_cJet_trueMinus_angle[j];
2067 double b_numerator = _bJet_truePlus2_recoPlus_angle[j]+_bJet_truePlus_recoPlus_angle[j]
2068 +_bJet_trueMinus2_recoMinus_angle[j]+_bJet_trueMinus_recoMinus_angle[j];
2069 double b_domininator = _bJet_truePlus2_angle[j] +_bJet_truePlus_angle[j]
2070 +_bJet_trueMinus2_angle[j] +_bJet_trueMinus_angle[j];
2072 double b_leakage = 1. - b_numerator/b_domininator;
2073 double c_leakage = 1. - c_numerator/c_domininator;
2078 _pCJetLeakageRate->fill(_pCJetLeakageRate->axis().binLowerEdge(j)+_pCJetLeakageRate->axis().binWidth(j)/2.,c_leakage);
2079 _pBJetLeakageRate->fill(_pBJetLeakageRate->axis().binLowerEdge(j)+_pBJetLeakageRate->axis().binWidth(j)/2.,b_leakage);
2085 void LCFIAIDAPlotProcessor::FillVertexPlots(LCEvent* pEvent,
unsigned int jetNumber)
2087 int jetType=FindTrueJetType( pEvent, jetNumber );
2089 if( jetType==0 )
return;
2091 std::vector<double> mcDecayLengthVector;
2092 std::vector<double> bMCDecayLengthVector;
2093 std::vector<double> cMCDecayLengthVector;
2094 double b_mc_decay_length, c_mc_decay_length;
2096 FindTrueJetDecayLength(pEvent, jetNumber, mcDecayLengthVector, bMCDecayLengthVector, cMCDecayLengthVector);
2097 FindTrueJetDecayLength2(pEvent, jetNumber, b_mc_decay_length, c_mc_decay_length);
2100 for (
unsigned int iTagCollection=0; iTagCollection < _FlavourTagCollectionNames.size(); ++iTagCollection) {
2102 unsigned int NumVertices = FindNumVertex(pEvent, jetNumber, iTagCollection);
2104 if( jetType==B_JET ) {
2105 _pBDecayLengthAll[iTagCollection]->fill(b_mc_decay_length);
2106 if (NumVertices > 1) _pBDecayLengthTwoVertices[iTagCollection]->fill(b_mc_decay_length);
2107 }
else if( jetType==C_JET ) {
2108 _pCDecayLengthAll[iTagCollection]->fill(c_mc_decay_length);
2109 if (NumVertices > 1) _pCDecayLengthTwoVertices[iTagCollection]->fill(c_mc_decay_length);
2118 if (!zvresjetCollection.is_valid() || !zvresdcCollection.is_valid() ) {
2120 std::cerr <<
"Input Collections " << _ZVRESSelectedJetsCollection <<
" and/or " << _ZVRESDecayChainCollection <<
" not valid. Will not fill vertex plots." <<std::endl;
2121 std::cerr <<
"Set MakeAdditionalPlots to false to remove this warning" << std::endl;
2130 for (
int ijet = 0 ; ijet < zvresdcCollection.getNumberOfElements(); ijet++) {
2133 int jetType = FindTrueJetType( pEvent, ijet );
2135 lcio::ReconstructedParticle* pJet=zvresdcCollection.getElementAt(ijet);
2137 std::vector<Vertex*> myVertexVector;
2138 std::pair<Vertex*,Vertex*> myVertexPair;
2139 std::vector< std::pair<Vertex*,Vertex*> > myVertexPairVector;
2141 std::vector<lcio::ReconstructedParticle*> decayChainPartColl = pJet->getParticles();
2142 for (
unsigned int ipart = 0 ; ipart < decayChainPartColl.size(); ipart++) {
2144 lcio::ReconstructedParticle* decayChainPart =
dynamic_cast< lcio::ReconstructedParticle*
> (decayChainPartColl[ipart]);
2148 bool alreadyContainsStart =
false;
2149 bool alreadyContainsEnd =
false;
2151 for (std::vector<Vertex*>::const_iterator iVert=myVertexVector.begin(); iVert!=myVertexVector.end(); iVert++) {
2152 if (decayChainPart->getStartVertex() == *iVert) alreadyContainsStart =
true;
2153 if (decayChainPart->getEndVertex() == *iVert) alreadyContainsEnd =
true;
2156 if (!alreadyContainsStart && decayChainPart->getStartVertex()!=0) myVertexVector.push_back(decayChainPart->getStartVertex());
2157 if (!alreadyContainsEnd && decayChainPart->getEndVertex()!=0) myVertexVector.push_back(decayChainPart->getEndVertex());
2160 myVertexPair.first = decayChainPart->getStartVertex();
2161 myVertexPair.second = decayChainPart->getEndVertex();
2162 bool pairAlreadyContained =
false;
2164 for (std::vector< std::pair<Vertex*,Vertex*> >::const_iterator iter = myVertexPairVector.begin(); iter != myVertexPairVector.end(); iter++) {
2165 if ((*iter).first == myVertexPair.first && (*iter).second == myVertexPair.second) pairAlreadyContained =
true;
2168 if (myVertexPair.first == 0 || myVertexPair.second == 0) pairAlreadyContained =
true;
2170 if (!pairAlreadyContained) myVertexPairVector.push_back(myVertexPair);
2175 float reconstructedSecondaryDecayLength = -1.;
2176 float reconstructedSecTerDecayLength = - 1.;
2177 std::vector< Vertex* > secondaryVertexVector;
2178 std::vector< float > secondaryDecayLengthVector;
2180 for (std::vector< std::pair<Vertex*,Vertex*> >::const_iterator iter = myVertexPairVector.begin();
2181 iter != myVertexPairVector.end(); iter++) {
2183 const Vertex* startV = (*iter).first;
2185 if (startV->isPrimary()) {
2187 const Vertex* endV = (*iter).second;
2189 const float* startPos = startV->getPosition();
2190 const float* endPos = endV->getPosition();
2192 secondaryVertexVector.push_back((*iter).second);
2194 reconstructedSecondaryDecayLength = CalculateDistance(startPos,endPos);
2196 _reconstructedSecondaryDecayLength->fill(reconstructedSecondaryDecayLength);
2199 if( jetType==B_JET ) {
2200 _recoDecayLengthBJet->fill(reconstructedSecondaryDecayLength);
2202 }
else if( jetType==C_JET ) {
2203 _recoDecayLengthCJet->fill(reconstructedSecondaryDecayLength);
2206 _recoDecayLengthLightJet->fill(reconstructedSecondaryDecayLength);
2214 for (std::vector< std::pair<Vertex*,Vertex*> >::const_iterator iter = myVertexPairVector.begin();
2215 iter != myVertexPairVector.end(); iter++) {
2217 const Vertex* startV = (*iter).first;
2219 bool isSecondary=
false;
2220 for (std::vector< Vertex* >::const_iterator iter2 = secondaryVertexVector.begin(); iter2 !=secondaryVertexVector.end(); iter2++)
2222 if (startV == (*iter2)) isSecondary=
true;
2226 const Vertex* endV = (*iter).second;
2228 const float* startPos = startV->getPosition();
2229 const float* endPos = endV->getPosition();
2230 reconstructedSecTerDecayLength = CalculateDistance(startPos,endPos);
2231 secondaryDecayLengthVector.push_back(reconstructedSecTerDecayLength);
2233 _reconstructedSecTerDecayLength->fill(reconstructedSecTerDecayLength);
2235 if (jetType==B_JET) _recoDecayLengthBCJet->fill(reconstructedSecTerDecayLength);
2241 int nRecoVertices = myVertexVector.size();
2243 if( jetType==B_JET ) {
2244 _nVerticesBJet->fill(nRecoVertices);
2245 _decayLengthBJet2D->fill(b_mc_decay_length,reconstructedSecondaryDecayLength);
2246 _decayLengthBJetCloud2D->fill(b_mc_decay_length,reconstructedSecondaryDecayLength);
2248 _decayLengthCJet2D->fill(c_mc_decay_length,reconstructedSecTerDecayLength);
2249 _decayLengthCJetCloud2D->fill(c_mc_decay_length,reconstructedSecTerDecayLength);
2251 _decayLengthBJetTrue->fill(b_mc_decay_length);
2252 _decayLengthBCJetTrue->fill(c_mc_decay_length);
2254 }
else if( jetType==C_JET ) {
2255 _nVerticesCJet->fill(nRecoVertices);
2256 _decayLengthBJetTrue->fill(c_mc_decay_length);
2259 _nVerticesLightJet->fill(nRecoVertices);
2269 if (vertexCol.is_valid()) {
2271 float primaryVertexPostion[] = {0.,0.,0.};
2273 for (
int ii=0 ; ii < vertexCol.getNumberOfElements() ; ii++){
2275 Vertex* myVertex =
dynamic_cast<Vertex*
>(vertexCol.getElementAt(ii));
2277 const float* vertexPosition = myVertex->getPosition();
2278 if (myVertex->isPrimary()) {
2279 primaryVertexPostion[0] = vertexPosition[0];
2280 primaryVertexPostion[1] = vertexPosition[1];
2281 primaryVertexPostion[2] = vertexPosition[2];
2285 for (
int ii=0 ; ii < vertexCol.getNumberOfElements() ; ii++){
2287 Vertex* myVertex =
dynamic_cast<Vertex*
>(vertexCol.getElementAt(ii));
2288 const float* vertexPosition = myVertex->getPosition();
2290 double px = double(vertexPosition[0]);
2291 double py = double(vertexPosition[1]);
2292 double pz = double(vertexPosition[2]);
2293 double ex = sqrt((
double)myVertex->getCovMatrix()[0]);
2294 double ey = sqrt((
double)myVertex->getCovMatrix()[2]);
2295 double ez = sqrt((
double)myVertex->getCovMatrix()[5]);
2298 if (! myVertex->isPrimary() ) {
2300 double distanceIP = sqrt(pow(px-primaryVertexPostion[0],2)+pow(py-primaryVertexPostion[1],2)+pow(pz-primaryVertexPostion[1],2));
2302 _pVertexDistanceFromIP->fill(distanceIP);
2303 _pVertexPositionX->fill(px);
2304 _pVertexPositionY->fill(py);
2305 _pVertexPositionZ->fill(pz);
2309 _pPrimaryVertexPositionX->fill(px);
2310 _pPrimaryVertexPositionY->fill(py);
2311 _pPrimaryVertexPositionZ->fill(pz);
2313 _pPrimaryVertexPullX->fill(px/ex);
2314 _pPrimaryVertexPullY->fill(py/ey);
2315 _pPrimaryVertexPullZ->fill(pz/ez);
2320 std::cerr <<
"Input Collections " << _VertexColName <<
" is not valid. Not filling all vertex plots." <<std::endl;
2321 std::cerr <<
"Set MakeAdditionalPlots to false to remove this warning" << std::endl;
2329 void LCFIAIDAPlotProcessor::FillZVRESTable(LCEvent* pEvent)
2336 LCCollection* relationTrackCollection = pEvent->getCollection(_TrueTracksToMCPCollection);
2338 UTIL::LCRelationNavigator* MCRelationNavigator = (relationTrackCollection!=NULL) ?
new LCRelationNavigator(relationTrackCollection) : NULL;
2341 if (MCRelationNavigator==NULL)
2342 std::cerr << __FILE__ <<
"(" << __LINE__ <<
"), run number: " << pEvent->getRunNumber() <<
", event number: " << pEvent->getEventNumber()<<
" LCCollection " << _TrueTracksToMCPCollection <<
" is not valid." << std::endl;
2344 if (!ftjetCollection.is_valid())
2345 std::cerr << __FILE__ <<
"(" << __LINE__ <<
"), run number: " << pEvent->getRunNumber() <<
", event number: " << pEvent->getEventNumber()<<
" ReconstructedParticle " << _JetCollectionName <<
" is not valid." << std::endl;
2346 if (!zvresjetCollection.is_valid())
2347 std::cerr << __FILE__ <<
"(" << __LINE__ <<
"), run number: " << pEvent->getRunNumber() <<
", event number: " << pEvent->getEventNumber()<<
" ReconstructedParticle " << _ZVRESSelectedJetsCollection <<
" is not valid." << std::endl;
2348 if (!zvresdcrptrackCollection.is_valid())
2349 std::cerr << __FILE__ <<
"(" << __LINE__ <<
"), run number: " << pEvent->getRunNumber() <<
", event number: " << pEvent->getEventNumber()<<
" ReconstructedParticle " << _ZVRESDecayChainRPTracksCollection <<
" is not valid." << std::endl;
2350 if (!zvresdcCollection.is_valid())
2351 std::cerr << __FILE__ <<
"(" << __LINE__ <<
"), run number: " << pEvent->getRunNumber() <<
", event number: " << pEvent->getEventNumber()<<
" ReconstructedParticle " << _ZVRESDecayChainCollection <<
" is not valid." << std::endl;
2354 if (!zvresdcCollection.is_valid() || !ftjetCollection.is_valid() || !zvresjetCollection.is_valid() || !zvresdcrptrackCollection.is_valid()) {
2356 std::cerr <<
"Vertex-Track Association cannot be analysed. Set MakeAdditionalPlots to false to remove this warning" << std::endl;
2362 std::vector<Vertex*> primaryVertices;
2363 std::vector<Vertex*> secondaryVertices;
2364 std::vector<Vertex*> tertiaryVertices;
2365 std::vector<Vertex*> allVertices;
2367 std::vector< std::pair<Vertex*,Vertex*> > myVertexPairVector;
2368 primaryVertices.clear();
2369 secondaryVertices.clear();
2370 tertiaryVertices.clear();
2371 allVertices.clear();
2372 myVertexPairVector.clear();
2377 for (
int iDecayChain = 0 ; iDecayChain < zvresdcCollection.getNumberOfElements(); iDecayChain++) {
2379 lcio::ReconstructedParticle* pDecayChain=zvresdcCollection.getElementAt(iDecayChain);
2382 std::pair<Vertex*,Vertex*> myVertexPair;
2383 std::vector< std::pair<Vertex*,Vertex*> > myVertexPairVector;
2385 myVertexPairVector.clear();
2388 std::vector<lcio::ReconstructedParticle*> decayChainPartColl = pDecayChain->getParticles();
2390 for (
unsigned int iPart = 0 ; iPart < decayChainPartColl.size(); iPart++) {
2392 lcio::ReconstructedParticle* decayChainPart =
dynamic_cast< lcio::ReconstructedParticle*
> (decayChainPartColl[iPart]);
2394 myVertexPair.first = decayChainPart->getStartVertex();
2395 myVertexPair.second = decayChainPart->getEndVertex();
2397 bool pairAlreadyContained =
false;
2399 for (std::vector< std::pair<Vertex*,Vertex*> >::const_iterator iter = myVertexPairVector.begin(); iter != myVertexPairVector.end(); iter++) {
2400 if ((*iter).first == myVertexPair.first && (*iter).second == myVertexPair.second) pairAlreadyContained =
true;
2403 if (!pairAlreadyContained) myVertexPairVector.push_back(myVertexPair);
2409 for (std::vector< std::pair<Vertex*,Vertex*> >::const_iterator iter = myVertexPairVector.begin(); iter != myVertexPairVector.end(); iter++) {
2412 if ((*iter).first && (*iter).first->isPrimary())
2414 bool primaryAlreadyContained =
false;
2415 for (std::vector<Vertex*>::const_iterator piter=primaryVertices.begin(); piter != primaryVertices.end(); piter++)
2417 if (*piter == (*iter).first) primaryAlreadyContained =
true;
2419 if (!primaryAlreadyContained) primaryVertices.push_back((*iter).first);
2423 if ((*iter).second != 0)
2425 bool secondaryAlreadyContained =
false;
2426 for (std::vector<Vertex*>::const_iterator siter=secondaryVertices.begin(); siter != secondaryVertices.end(); siter++)
2428 if (*siter == (*iter).second) secondaryAlreadyContained =
true;
2430 if (!secondaryAlreadyContained) secondaryVertices.push_back((*iter).second);
2436 bool firstAlreadyContained=
false, secondAlreadyContained=
false;
2437 for (std::vector<Vertex*>::const_iterator aiter = allVertices.begin(); aiter != allVertices.end(); aiter++)
2439 if ( (*iter).first && *aiter == (*iter).first) firstAlreadyContained=
true;
2440 if ( (*iter).second && *aiter == (*iter).second ) secondAlreadyContained=
true;
2443 if ( !firstAlreadyContained && (*iter).first) allVertices.push_back((*iter).first);
2444 if (!secondAlreadyContained && (*iter).second) allVertices.push_back((*iter).second);
2448 for (std::vector< std::pair<Vertex*,Vertex*> >::const_iterator iter = myVertexPairVector.begin(); iter != myVertexPairVector.end(); iter++)
2450 for (std::vector<Vertex*>::const_iterator siter = secondaryVertices.begin(); siter != secondaryVertices.end(); siter++)
2452 if ((*iter).first && (*iter).first == (*siter) && (*iter).second)
2454 bool tertiaryAlreadyContained =
false;
2455 for (std::vector<Vertex*>::const_iterator titer=tertiaryVertices.begin(); titer != tertiaryVertices.end(); titer++)
2457 if (*titer == (*iter).second) tertiaryAlreadyContained=
false;
2459 if (!tertiaryAlreadyContained) tertiaryVertices.push_back((*iter).second);
2470 for (
int iJet = 0 ; iJet < zvresjetCollection.getNumberOfElements(); iJet++) {
2471 lcio::ReconstructedParticle* pJet=zvresjetCollection.getElementAt(iJet);
2474 std::vector<lcio::ReconstructedParticle*> jetParticles = pJet->getParticles();
2477 int jetType = FindTrueJetType( pEvent, iJet );
2479 if (jetType == B_JET || jetType == C_JET) {
2481 for (
unsigned int iPart = 0 ; iPart < jetParticles.size(); iPart++)
2483 lcio::ReconstructedParticle* myJetParticle =
dynamic_cast< lcio::ReconstructedParticle*
> (jetParticles[iPart]);
2486 bool decayChainMatch =
false;
2488 bool decayChainRPMatch =
false;
2491 lcio::ReconstructedParticle* myDecayChain = 0;
2493 lcio::ReconstructedParticle* myDecayChainRP = 0;
2496 std::vector<Vertex*> myVertexVector;
2497 myVertexVector.clear();
2499 for (
int iDCJet = 0 ; iDCJet < zvresdcrptrackCollection.getNumberOfElements(); iDCJet++) {
2501 lcio::ReconstructedParticle* myDCJetRP = zvresdcrptrackCollection.getElementAt(iDCJet);
2503 lcio::ReconstructedParticle* myDCJetRPTrack = myDCJetRP->getParticles()[0];
2505 if (myDCJetRPTrack==myJetParticle) {
2506 decayChainRPMatch =
true;
2507 myDecayChainRP = myDCJetRP;
2514 if (decayChainRPMatch) {
2517 for (
int iDecayChain = 0 ; iDecayChain < zvresdcCollection.getNumberOfElements(); iDecayChain++) {
2518 lcio::ReconstructedParticle* pDecayChain=zvresdcCollection.getElementAt(iDecayChain);
2519 std::vector<lcio::ReconstructedParticle*> decayChainParticles = pDecayChain->getParticles();
2520 for (
unsigned int iDCPart = 0 ; iDCPart < decayChainParticles.size(); iDCPart++) {
2521 lcio::ReconstructedParticle* myDecayChainParticle =
dynamic_cast< lcio::ReconstructedParticle*
> (decayChainParticles[iDCPart]);
2522 if (myDecayChainParticle == myDecayChainRP) {
2523 decayChainMatch =
true;
2524 myDecayChain = pDecayChain;
2533 lcio::ReconstructedParticle* pDecayChain=zvresdcCollection.getElementAt(iJet);
2534 myDecayChain = pDecayChain;
2535 if (pDecayChain) decayChainMatch =
true;
2538 if (decayChainMatch) {
2540 std::vector<lcio::ReconstructedParticle*> myDecayChainParticles = myDecayChain->getParticles();
2541 for (
unsigned int iDCPart = 0 ; iDCPart < myDecayChainParticles.size(); iDCPart++) {
2543 lcio::ReconstructedParticle* myDecayChainParticle =
dynamic_cast< lcio::ReconstructedParticle*
> (myDecayChainParticles[iDCPart]);
2545 bool startAlreadyContained =
false;
2546 bool endAlreadyContained =
false;
2548 for (std::vector<Vertex*>::const_iterator iter=myVertexVector.begin(); iter!=myVertexVector.end(); iter++)
2550 if (*iter == myDecayChainParticle->getStartVertex()) startAlreadyContained =
true;
2552 if (!startAlreadyContained && myDecayChainParticle->getStartVertex()) myVertexVector.push_back(myDecayChainParticle->getStartVertex());
2555 for (std::vector<Vertex*>::const_iterator iter=myVertexVector.begin(); iter!=myVertexVector.end(); iter++)
2557 if (*iter == myDecayChainParticle->getEndVertex()) endAlreadyContained =
true;
2559 if (!endAlreadyContained && myDecayChainParticle->getEndVertex()) myVertexVector.push_back(myDecayChainParticle->getEndVertex());
2564 int numVertex=myVertexVector.size();
2566 bool trackFromPrimaryLight =
false, trackFromSecondaryC =
false, trackFromSecondaryB =
false, trackFromTertiaryC =
false;
2567 bool trackHasNoAssociatedMCP =
false;
2570 std::vector<Track*> myTracks;
2572 if (decayChainRPMatch) myTracks = myDecayChainRP->getParticles()[0]->getTracks();
2574 myTracks = myJetParticle->getTracks();
2578 if (myTracks.size()==0)
2583 trackHasNoAssociatedMCP =
true;
2587 else if (myTracks.size()>0)
2590 lcio::Track* Track = myTracks[0];
2592 std::vector<lcio::LCObject*> RelatedMCParticles = MCRelationNavigator->getRelatedToObjects(Track);
2599 lcio::MCParticle* MCP =
dynamic_cast<lcio::MCParticle*
>(RelatedMCParticles[0]);
2600 std::vector<lcio::MCParticle*> Parents =
dynamic_cast<lcio::MCParticle*
>(RelatedMCParticles[0])->getParents();
2602 if (Parents.empty())
2603 std::cerr <<
" In " << __FILE__ <<
"(" << __LINE__ <<
"), run number: " << pEvent->getRunNumber() <<
", event number: " << pEvent->getEventNumber()<<
" : MCP has no parents" << std::endl;
2607 int truePDGCode=MCP->getPDG();
2608 int truePDGFlavour = GetPDGFlavour(truePDGCode);
2609 std::vector<int> pdgCodeVector;
2610 std::vector<int> pdgFlavourVector;
2611 std::vector<float> pdgDecayLengths;
2613 pdgCodeVector.push_back(truePDGCode);
2614 pdgFlavourVector.push_back(truePDGFlavour);
2616 std::vector<lcio::MCParticle*> ParentVec0 = MCP->getParents();
2623 if (ParentVec0.size()>0) {
2624 lcio::MCParticle* Parent0 = ParentVec0[0];
2628 pdgCodeVector.push_back(Parent0->getPDG());
2629 pdgFlavourVector.push_back(GetPDGFlavour(Parent0->getPDG()));
2630 std::vector<lcio::MCParticle*> ParentVec1 = Parent0->getParents();
2631 pdgDecayLengths.push_back(CalculateDistance(Parent0->getVertex(),Parent0->getEndpoint()));
2633 if (ParentVec1.size()>0) {
2634 lcio::MCParticle* Parent1 = ParentVec1[0];
2637 pdgCodeVector.push_back(Parent1->getPDG());
2638 pdgFlavourVector.push_back(GetPDGFlavour(Parent1->getPDG()));
2639 std::vector<lcio::MCParticle*> ParentVec2 = Parent1->getParents();
2640 pdgDecayLengths.push_back(CalculateDistance(Parent1->getVertex(),Parent1->getEndpoint()));
2642 if (ParentVec2.size()>0) {
2643 lcio::MCParticle* Parent2 = ParentVec2[0];
2646 pdgCodeVector.push_back(Parent2->getPDG());
2647 pdgFlavourVector.push_back(GetPDGFlavour(Parent2->getPDG()));
2648 std::vector<lcio::MCParticle*> ParentVec3 = Parent2->getParents();
2649 pdgDecayLengths.push_back(CalculateDistance(Parent2->getVertex(),Parent2->getEndpoint()));
2651 if (ParentVec3.size()>0) {
2652 lcio::MCParticle* Parent3 = ParentVec3[0];
2655 pdgCodeVector.push_back(Parent3->getPDG());
2656 pdgFlavourVector.push_back(GetPDGFlavour(Parent3->getPDG()));
2657 std::vector<lcio::MCParticle*> ParentVec4 = Parent3->getParents();
2658 pdgDecayLengths.push_back(CalculateDistance(Parent3->getVertex(),Parent3->getEndpoint()));
2660 if (ParentVec4.size()>0) {
2661 lcio::MCParticle* Parent4 = ParentVec4[0];
2664 pdgCodeVector.push_back(Parent4->getPDG());
2665 pdgFlavourVector.push_back(GetPDGFlavour(Parent4->getPDG()));
2666 std::vector<lcio::MCParticle*> ParentVec5 = Parent4->getParents();
2667 pdgDecayLengths.push_back(CalculateDistance(Parent4->getVertex(),Parent4->getEndpoint()));
2669 if (ParentVec5.size()>0) {
2670 lcio::MCParticle* Parent5 = ParentVec5[0];
2673 pdgCodeVector.push_back(Parent5->getPDG());
2674 pdgFlavourVector.push_back(GetPDGFlavour(Parent5->getPDG()));
2675 std::vector<lcio::MCParticle*> ParentVec6 = Parent5->getParents();
2676 pdgDecayLengths.push_back(CalculateDistance(Parent5->getVertex(),Parent5->getEndpoint()));
2678 if (ParentVec6.size()>0) {
2679 lcio::MCParticle* Parent6 = ParentVec6[0];
2682 pdgCodeVector.push_back(Parent6->getPDG());
2683 pdgFlavourVector.push_back(GetPDGFlavour(Parent6->getPDG()));
2684 pdgDecayLengths.push_back(CalculateDistance(Parent6->getVertex(),Parent6->getEndpoint()));
2703 bool isLight =
false;
2708 for (std::vector<int>::const_iterator iter = pdgFlavourVector.begin();
2709 iter != pdgFlavourVector.end(); iter++) {
2710 if ((*iter) == 1) isLight =
true;
2711 if ((*iter) == C_JET) isC =
true;
2712 if ((*iter) == B_JET) isB =
true;
2720 if (isLight && !isC && !isB) trackFromPrimaryLight =
true;
2721 if (isC && !isB) trackFromSecondaryC =
true;
2722 if (!isC && isB) trackFromSecondaryB =
true;
2723 if (isC && isB) trackFromTertiaryC =
true;
2730 bool isPrimary=
false, isSecondary=
false, isTertiary=
false;
2732 if (decayChainRPMatch) {
2734 lcio::Vertex* partVertex = myDecayChainRP->getStartVertex();
2736 for (std::vector<Vertex*>::const_iterator piter = primaryVertices.begin(); piter != primaryVertices.end(); piter++)
2738 if ((*piter) == partVertex) isPrimary=
true;
2740 for (std::vector<Vertex*>::const_iterator siter = secondaryVertices.begin(); siter != secondaryVertices.end(); siter++)
2742 if ((*siter) == partVertex) isSecondary=
true;
2744 for (std::vector<Vertex*>::const_iterator titer = tertiaryVertices.begin(); titer != tertiaryVertices.end(); titer++)
2746 if ((*titer) == partVertex) isTertiary=
true;
2753 if (trackHasNoAssociatedMCP)
2755 if (isPrimary) _nb_twoVertex_Primary_noMCP++;
2756 else if (isSecondary) _nb_twoVertex_Secondary_noMCP++;
2757 else _nb_twoVertex_Isolated_noMCP++;
2760 if (trackFromSecondaryB && jetType == B_JET )
2762 if (isPrimary) _nb_twoVertex_bTrack_Primary++;
2763 else if (isSecondary) _nb_twoVertex_bTrack_Secondary++;
2765 else _nb_twoVertex_bTrack_Isolated++;
2767 else if ((trackFromSecondaryC || trackFromTertiaryC) && jetType == B_JET )
2769 if (isPrimary) _nb_twoVertex_cTrack_Primary++;
2770 else if (isSecondary) _nb_twoVertex_cTrack_Secondary++;
2772 else _nb_twoVertex_cTrack_Isolated++;
2774 else if (trackFromPrimaryLight && jetType == B_JET)
2776 if (isPrimary) _nb_twoVertex_lTrack_Primary++;
2777 else if (isSecondary) _nb_twoVertex_lTrack_Secondary++;
2779 else _nb_twoVertex_lTrack_Isolated++;
2782 if (trackHasNoAssociatedMCP && jetType == C_JET) {
2784 if (isPrimary) _nc_twoVertex_Primary_noMCP++;
2785 else if (isSecondary) _nc_twoVertex_Secondary_noMCP++;
2786 else _nc_twoVertex_Isolated_noMCP++;
2789 if (trackFromSecondaryB && jetType == C_JET )
2791 if (isPrimary) _nc_twoVertex_bTrack_Primary++;
2792 else if (isSecondary) _nc_twoVertex_bTrack_Secondary++;
2794 else _nc_twoVertex_bTrack_Isolated++;
2796 else if ((trackFromSecondaryC || trackFromTertiaryC) && jetType == C_JET)
2798 if (isPrimary) _nc_twoVertex_cTrack_Primary++;
2799 else if (isSecondary) _nc_twoVertex_cTrack_Secondary++;
2801 else _nc_twoVertex_cTrack_Isolated++;
2803 else if (trackFromPrimaryLight && jetType == C_JET)
2805 if (isPrimary) _nc_twoVertex_lTrack_Primary++;
2806 else if (isSecondary) _nc_twoVertex_lTrack_Secondary++;
2808 else _nc_twoVertex_lTrack_Isolated++;
2811 else if (numVertex == 3 )
2813 if (trackHasNoAssociatedMCP && jetType == B_JET)
2815 if (isPrimary) _nb_twoVertex_Primary_noMCP++;
2816 else if (isSecondary) _nb_twoVertex_Secondary_noMCP++;
2817 else if (isTertiary) _nb_twoVertex_Tertiary_noMCP++;
2818 else _nb_twoVertex_Isolated_noMCP++;
2821 if (trackFromSecondaryB && jetType == B_JET)
2823 if (isPrimary) _nb_threeVertex_bTrack_Primary++;
2824 else if (isSecondary) _nb_threeVertex_bTrack_Secondary++;
2825 else if (isTertiary) _nb_threeVertex_bTrack_Tertiary++;
2826 else _nb_threeVertex_bTrack_Isolated++;
2828 else if ((trackFromSecondaryC || trackFromTertiaryC) && jetType == B_JET)
2830 if (isPrimary) _nb_threeVertex_cTrack_Primary++;
2831 else if (isSecondary) _nb_threeVertex_cTrack_Secondary++;
2832 else if (isTertiary) _nb_threeVertex_cTrack_Tertiary++;
2833 else _nb_threeVertex_cTrack_Isolated++;
2835 else if (trackFromPrimaryLight && jetType == B_JET)
2837 if (isPrimary) _nb_threeVertex_lTrack_Primary++;
2838 else if (isSecondary) _nb_threeVertex_lTrack_Secondary++;
2839 else if (isTertiary) _nb_threeVertex_lTrack_Tertiary++;
2840 else _nb_threeVertex_lTrack_Isolated++;
2842 if (trackHasNoAssociatedMCP && jetType == C_JET)
2844 if (isPrimary) _nc_twoVertex_Primary_noMCP++;
2845 else if (isSecondary) _nc_twoVertex_Secondary_noMCP++;
2846 else if (isTertiary) _nc_twoVertex_Tertiary_noMCP++;
2847 else _nc_twoVertex_Isolated_noMCP++;
2849 else if (trackFromSecondaryB && jetType == C_JET)
2851 if (isPrimary) _nc_threeVertex_bTrack_Primary++;
2852 else if (isSecondary) _nc_threeVertex_bTrack_Secondary++;
2853 else if (isTertiary) _nc_threeVertex_bTrack_Tertiary++;
2854 else _nc_threeVertex_bTrack_Isolated++;
2856 else if ((trackFromSecondaryC || trackFromTertiaryC) && jetType == C_JET)
2858 if (isPrimary) _nc_threeVertex_cTrack_Primary++;
2859 else if (isSecondary) _nc_threeVertex_cTrack_Secondary++;
2860 else if (isTertiary) _nc_threeVertex_cTrack_Tertiary++;
2861 else _nc_threeVertex_cTrack_Isolated++;
2863 else if (trackFromPrimaryLight && jetType == C_JET)
2865 if (isPrimary) _nc_threeVertex_lTrack_Primary++;
2866 else if (isSecondary) _nc_threeVertex_lTrack_Secondary++;
2867 else if (isTertiary) _nc_threeVertex_lTrack_Tertiary++;
2868 else _nc_threeVertex_lTrack_Isolated++;
2877 float LCFIAIDAPlotProcessor::CalculateDistance(
const float* pos1,
const float* pos2){
2878 return sqrt(pow((pos1[0]-pos2[0]),2)+pow((pos1[1]-pos2[1]),2)+pow((pos1[2]-pos2[2]),2));
2881 double LCFIAIDAPlotProcessor::CalculateDistance(
const double* pos1,
const double* pos2){
2882 return sqrt(pow((pos1[0]-pos2[0]),2)+pow((pos1[1]-pos2[1]),2)+pow((pos1[2]-pos2[2]),2));
2885 void LCFIAIDAPlotProcessor::PrintZVRESTable()
2889 std::filebuf* fb =
new std::filebuf;
2891 std::ostream outputFile( (!_TrackVertexOutputFile.empty()) ?
2892 fb->open(_TrackVertexOutputFile.c_str(),
2893 std::ios_base::out|std::ios_base::trunc):
2907 float b_twoVertex_Primary = _nb_twoVertex_bTrack_Primary+_nb_twoVertex_cTrack_Primary+_nb_twoVertex_lTrack_Primary;
2908 float b_twoVertex_Secondary = _nb_twoVertex_bTrack_Secondary+_nb_twoVertex_cTrack_Secondary+_nb_twoVertex_lTrack_Secondary;
2909 float b_twoVertex_Isolated = _nb_twoVertex_bTrack_Isolated+_nb_twoVertex_cTrack_Isolated+_nb_twoVertex_lTrack_Isolated;
2910 float b_twoVertex = b_twoVertex_Primary+b_twoVertex_Secondary+b_twoVertex_Isolated;
2912 float b_threeVertex_Primary = _nb_threeVertex_bTrack_Primary+_nb_threeVertex_cTrack_Primary+_nb_threeVertex_lTrack_Primary;
2913 float b_threeVertex_Secondary = _nb_threeVertex_bTrack_Secondary+_nb_threeVertex_cTrack_Secondary+_nb_threeVertex_lTrack_Secondary;
2914 float b_threeVertex_Tertiary = _nb_threeVertex_bTrack_Tertiary+_nb_threeVertex_cTrack_Tertiary+_nb_threeVertex_lTrack_Tertiary;
2915 float b_threeVertex_Isolated = _nb_threeVertex_bTrack_Isolated+_nb_threeVertex_cTrack_Isolated+_nb_threeVertex_lTrack_Isolated;
2916 float b_threeVertex = b_threeVertex_Primary+b_threeVertex_Secondary+b_threeVertex_Tertiary+b_threeVertex_Isolated;
2920 float frac_b_twoVertex_bTrack_Primary = 100.*float(_nb_twoVertex_bTrack_Primary) / float(b_twoVertex_Primary);
2921 float frac_b_twoVertex_bTrack_Secondary = 100.*float(_nb_twoVertex_bTrack_Secondary) / float(b_twoVertex_Secondary);
2923 float frac_b_twoVertex_bTrack_Isolated = 100.*float(_nb_twoVertex_bTrack_Isolated) / float(b_twoVertex_Isolated);
2925 float frac_b_twoVertex_cTrack_Primary = 100.*float(_nb_twoVertex_cTrack_Primary) / float(b_twoVertex_Primary);
2926 float frac_b_twoVertex_cTrack_Secondary = 100.*float(_nb_twoVertex_cTrack_Secondary) / float(b_twoVertex_Secondary);
2928 float frac_b_twoVertex_cTrack_Isolated = 100.*float(_nb_twoVertex_cTrack_Isolated) / float(b_twoVertex_Isolated);
2930 float frac_b_twoVertex_lTrack_Primary = 100.*float(_nb_twoVertex_lTrack_Primary) / float(b_twoVertex_Primary);
2931 float frac_b_twoVertex_lTrack_Secondary = 100.*float(_nb_twoVertex_lTrack_Secondary) / float(b_twoVertex_Secondary);
2933 float frac_b_twoVertex_lTrack_Isolated = 100.*float(_nb_twoVertex_lTrack_Isolated) / float(b_twoVertex_Isolated);
2935 float frac_b_threeVertex_bTrack_Primary = 100.*float(_nb_threeVertex_bTrack_Primary) / float(b_threeVertex_Primary);
2936 float frac_b_threeVertex_bTrack_Secondary = 100.*float(_nb_threeVertex_bTrack_Secondary) / float(b_threeVertex_Secondary);
2937 float frac_b_threeVertex_bTrack_Tertiary = 100.*float(_nb_threeVertex_bTrack_Tertiary) / float(b_threeVertex_Tertiary);
2938 float frac_b_threeVertex_bTrack_Isolated = 100.*float(_nb_threeVertex_bTrack_Isolated) / float(b_threeVertex_Isolated);
2940 float frac_b_threeVertex_cTrack_Primary = 100.*float(_nb_threeVertex_cTrack_Primary) / float(b_threeVertex_Primary);
2941 float frac_b_threeVertex_cTrack_Secondary = 100.*float(_nb_threeVertex_cTrack_Secondary) / float(b_threeVertex_Secondary);
2942 float frac_b_threeVertex_cTrack_Tertiary = 100.*float(_nb_threeVertex_cTrack_Tertiary) / float(b_threeVertex_Tertiary);
2943 float frac_b_threeVertex_cTrack_Isolated = 100.*float(_nb_threeVertex_cTrack_Isolated) / float(b_threeVertex_Isolated);
2945 float frac_b_threeVertex_lTrack_Primary = 100.*float(_nb_threeVertex_lTrack_Primary) / float(b_threeVertex_Primary);
2946 float frac_b_threeVertex_lTrack_Secondary = 100.*float(_nb_threeVertex_lTrack_Secondary) / float(b_threeVertex_Secondary);
2947 float frac_b_threeVertex_lTrack_Tertiary = 100.*float(_nb_threeVertex_lTrack_Tertiary) / float(b_threeVertex_Tertiary);
2948 float frac_b_threeVertex_lTrack_Isolated = 100.*float(_nb_threeVertex_lTrack_Isolated) / float(b_threeVertex_Isolated);
2950 float frac_b_twoVertex_Primary = 100.*float(b_twoVertex_Primary) / float(b_twoVertex);
2951 float frac_b_twoVertex_Secondary = 100.*float(b_twoVertex_Secondary) / float(b_twoVertex);
2952 float frac_b_twoVertex_Isolated = 100.*float(b_twoVertex_Isolated) / float(b_twoVertex);
2954 float frac_b_threeVertex_Primary = 100.*float(b_threeVertex_Primary) / float(b_threeVertex);
2955 float frac_b_threeVertex_Secondary = 100.*float(b_threeVertex_Secondary) / float(b_threeVertex);
2956 float frac_b_threeVertex_Tertiary = 100.*float(b_threeVertex_Tertiary) / float(b_threeVertex);
2957 float frac_b_threeVertex_Isolated = 100.*float(b_threeVertex_Isolated) / float(b_threeVertex);
2995 float c_twoVertex_Primary = _nc_twoVertex_bTrack_Primary+_nc_twoVertex_cTrack_Primary+_nc_twoVertex_lTrack_Primary;
2996 float c_twoVertex_Secondary = _nc_twoVertex_bTrack_Secondary+_nc_twoVertex_cTrack_Secondary+_nc_twoVertex_lTrack_Secondary;
2997 float c_twoVertex_Isolated = _nc_twoVertex_bTrack_Isolated+_nc_twoVertex_cTrack_Isolated+_nc_twoVertex_lTrack_Isolated;
2998 float c_twoVertex = c_twoVertex_Primary+c_twoVertex_Secondary+c_twoVertex_Isolated;
3000 float c_threeVertex_Primary = _nc_threeVertex_bTrack_Primary+_nc_threeVertex_cTrack_Primary+_nc_threeVertex_lTrack_Primary;
3001 float c_threeVertex_Secondary = _nc_threeVertex_bTrack_Secondary+_nc_threeVertex_cTrack_Secondary+_nc_threeVertex_lTrack_Secondary;
3002 float c_threeVertex_Tertiary = _nc_threeVertex_bTrack_Tertiary+_nc_threeVertex_cTrack_Tertiary+_nc_threeVertex_lTrack_Tertiary;
3003 float c_threeVertex_Isolated = _nc_threeVertex_bTrack_Isolated+_nc_threeVertex_cTrack_Isolated+_nc_threeVertex_lTrack_Isolated;
3004 float c_threeVertex = c_threeVertex_Primary+c_threeVertex_Secondary+c_threeVertex_Tertiary+c_threeVertex_Isolated;
3013 float frac_c_twoVertex_cTrack_Primary = 100.*float(_nc_twoVertex_cTrack_Primary) / float(c_twoVertex_Primary);
3014 float frac_c_twoVertex_cTrack_Secondary = 100.*float(_nc_twoVertex_cTrack_Secondary) / float(c_twoVertex_Secondary);
3016 float frac_c_twoVertex_cTrack_Isolated = 100.*float(_nc_twoVertex_cTrack_Isolated) / float(c_twoVertex_Isolated);
3018 float frac_c_twoVertex_lTrack_Primary = 100.*float(_nc_twoVertex_lTrack_Primary) / float(c_twoVertex_Primary);
3019 float frac_c_twoVertex_lTrack_Secondary = 100.*float(_nc_twoVertex_lTrack_Secondary) / float(c_twoVertex_Secondary);
3021 float frac_c_twoVertex_lTrack_Isolated = 100.*float(_nc_twoVertex_lTrack_Isolated) / float(c_twoVertex_Isolated);
3028 float frac_c_threeVertex_cTrack_Primary = 100.*float(_nc_threeVertex_cTrack_Primary) / float(c_threeVertex_Primary);
3029 float frac_c_threeVertex_cTrack_Secondary = 100.*float(_nc_threeVertex_cTrack_Secondary) / float(c_threeVertex_Secondary);
3030 float frac_c_threeVertex_cTrack_Tertiary = 100.*float(_nc_threeVertex_cTrack_Tertiary) / float(c_threeVertex_Tertiary);
3031 float frac_c_threeVertex_cTrack_Isolated = 100.*float(_nc_threeVertex_cTrack_Isolated) / float(c_threeVertex_Isolated);
3033 float frac_c_threeVertex_lTrack_Primary = 100.*float(_nc_threeVertex_lTrack_Primary) / float(c_threeVertex_Primary);
3034 float frac_c_threeVertex_lTrack_Secondary = 100.*float(_nc_threeVertex_lTrack_Secondary) / float(c_threeVertex_Secondary);
3035 float frac_c_threeVertex_lTrack_Tertiary = 100.*float(_nc_threeVertex_lTrack_Tertiary) / float(c_threeVertex_Tertiary);
3036 float frac_c_threeVertex_lTrack_Isolated = 100.*float(_nc_threeVertex_lTrack_Isolated) / float(c_threeVertex_Isolated);
3038 float frac_c_twoVertex_Primary = 100.*float(c_twoVertex_Primary) / float(c_twoVertex);
3039 float frac_c_twoVertex_Secondary = 100.*float(c_twoVertex_Secondary) / float(c_twoVertex);
3040 float frac_c_twoVertex_Isolated = 100.*float(c_twoVertex_Isolated) / float(c_twoVertex);
3042 float frac_c_threeVertex_Primary = 100.*float(c_threeVertex_Primary) / float(c_threeVertex);
3043 float frac_c_threeVertex_Secondary = 100.*float(c_threeVertex_Secondary) / float(c_threeVertex);
3044 float frac_c_threeVertex_Tertiary = 100.*float(c_threeVertex_Tertiary) / float(c_threeVertex);
3045 float frac_c_threeVertex_Isolated = 100.*float(c_threeVertex_Isolated) / float(c_threeVertex);
3083 outputFile << std::endl;
3084 outputFile <<
" ---------------------------------------------------------------------------" << std::endl;
3085 outputFile <<
" B-JETS " << std::endl;
3086 outputFile <<
" ---------------------------------------------------------------------------" << std::endl;
3087 outputFile <<
" Purity of reconstructed track-vertex association (%) " << std::endl;
3088 outputFile <<
" ---------------------------------------------------------------------------" << std::endl;
3089 outputFile <<
" Monte Carlo Reconstructed track-vertex association " << std::endl;
3090 outputFile <<
" track origin Two-vertex case Three-vertex case " << std::endl;
3091 outputFile <<
" ---------------------------------------------------------------------------" << std::endl;
3092 outputFile <<
" Pri. Sec. Iso. Pri. Sec. Ter. Iso. " << std::endl;
3093 outputFile <<
" ---------------------------------------------------------------------------" << std::endl;
3094 outputFile <<
" Primary ";
3095 outputFile.width(7);
3096 outputFile.precision(3);
3097 outputFile << frac_b_twoVertex_lTrack_Primary <<
" " ;
3098 outputFile.width(7);
3099 outputFile.precision(3);
3100 outputFile << frac_b_twoVertex_lTrack_Secondary <<
" " ;
3101 outputFile.width(7);
3102 outputFile.precision(3);
3103 outputFile << frac_b_twoVertex_lTrack_Isolated <<
" " ;
3104 outputFile.width(7);
3105 outputFile.precision(3);
3106 outputFile << frac_b_threeVertex_lTrack_Primary <<
" " ;
3107 outputFile.width(7);
3108 outputFile.precision(3);
3109 outputFile << frac_b_threeVertex_lTrack_Secondary <<
" " ;
3110 outputFile.width(7);
3111 outputFile.precision(3);
3112 outputFile << frac_b_threeVertex_lTrack_Tertiary <<
" " ;
3113 outputFile.width(7);
3114 outputFile.precision(3);
3115 outputFile << frac_b_threeVertex_lTrack_Isolated <<
" " ;
3116 outputFile << std::endl;
3118 outputFile <<
" B decay ";
3119 outputFile.width(7);
3120 outputFile.precision(3);
3121 outputFile << frac_b_twoVertex_bTrack_Primary <<
" " ;
3122 outputFile.width(7);
3123 outputFile.precision(3);
3124 outputFile << frac_b_twoVertex_bTrack_Secondary <<
" " ;
3125 outputFile.width(7);
3126 outputFile.precision(3);
3127 outputFile << frac_b_twoVertex_bTrack_Isolated <<
" " ;
3128 outputFile.width(7);
3129 outputFile.precision(3);
3130 outputFile << frac_b_threeVertex_bTrack_Primary <<
" " ;
3131 outputFile.width(7);
3132 outputFile.precision(3);
3133 outputFile << frac_b_threeVertex_bTrack_Secondary <<
" " ;
3134 outputFile.width(7);
3135 outputFile.precision(3);
3136 outputFile << frac_b_threeVertex_bTrack_Tertiary <<
" " ;
3137 outputFile.width(7);
3138 outputFile.precision(3);
3139 outputFile << frac_b_threeVertex_bTrack_Isolated <<
" " ;
3140 outputFile << std::endl;
3142 outputFile <<
" D decay ";
3143 outputFile.width(7);
3144 outputFile.precision(3);
3145 outputFile << frac_b_twoVertex_cTrack_Primary <<
" " ;
3146 outputFile.width(7);
3147 outputFile.precision(3);
3148 outputFile << frac_b_twoVertex_cTrack_Secondary <<
" " ;
3149 outputFile.width(7);
3150 outputFile.precision(3);
3151 outputFile << frac_b_twoVertex_cTrack_Isolated <<
" " ;
3152 outputFile.width(7);
3153 outputFile.precision(3);
3154 outputFile << frac_b_threeVertex_cTrack_Primary <<
" " ;
3155 outputFile.width(7);
3156 outputFile.precision(3);
3157 outputFile << frac_b_threeVertex_cTrack_Secondary <<
" " ;
3158 outputFile.width(7);
3159 outputFile.precision(3);
3160 outputFile << frac_b_threeVertex_cTrack_Tertiary <<
" " ;
3161 outputFile.width(7);
3162 outputFile.precision(3);
3163 outputFile << frac_b_threeVertex_cTrack_Isolated <<
" " ;
3164 outputFile << std::endl;
3165 outputFile <<
" ----------------------------------------------------------------------------" << std::endl;
3166 outputFile <<
" All above ";
3167 outputFile.width(7);
3168 outputFile.precision(3);
3169 outputFile << frac_b_twoVertex_Primary <<
" " ;
3170 outputFile.width(7);
3171 outputFile.precision(3);
3172 outputFile << frac_b_twoVertex_Secondary <<
" " ;
3173 outputFile.width(7);
3174 outputFile.precision(3);
3175 outputFile << frac_b_twoVertex_Isolated <<
" " ;
3176 outputFile.width(7);
3177 outputFile.precision(3);
3178 outputFile << frac_b_threeVertex_Primary <<
" " ;
3179 outputFile.width(7);
3180 outputFile.precision(3);
3181 outputFile << frac_b_threeVertex_Secondary <<
" " ;
3182 outputFile.width(7);
3183 outputFile.precision(3);
3184 outputFile << frac_b_threeVertex_Tertiary <<
" " ;
3185 outputFile.width(7);
3186 outputFile.precision(3);
3187 outputFile << frac_b_threeVertex_Isolated <<
" " ;
3188 outputFile << std::endl;
3189 outputFile <<
" ----------------------------------------------------------------------------" << std::endl;
3190 outputFile << std::endl;
3191 outputFile << std::endl;
3194 outputFile << std::endl;
3195 outputFile <<
" ---------------------------------------------------------------------------" << std::endl;
3196 outputFile <<
" C-JETS " << std::endl;
3197 outputFile <<
" ---------------------------------------------------------------------------" << std::endl;
3198 outputFile <<
" Purity of reconstructed track-vertex association (%) " << std::endl;
3199 outputFile <<
" ---------------------------------------------------------------------------" << std::endl;
3200 outputFile <<
" Monte Carlo Reconstructed track-vertex association " << std::endl;
3201 outputFile <<
" track origin Two-vertex case Three-vertex case " << std::endl;
3202 outputFile <<
" ---------------------------------------------------------------------------" << std::endl;
3203 outputFile <<
" Pri. Sec. Iso. Pri. Sec. Ter. Iso. " << std::endl;
3204 outputFile <<
" ---------------------------------------------------------------------------" << std::endl;
3205 outputFile <<
" Primary ";
3206 outputFile.width(7);
3207 outputFile.precision(3);
3208 outputFile << frac_c_twoVertex_lTrack_Primary <<
" " ;
3209 outputFile.width(7);
3210 outputFile.precision(3);
3211 outputFile << frac_c_twoVertex_lTrack_Secondary <<
" " ;
3212 outputFile.width(7);
3213 outputFile.precision(3);
3214 outputFile << frac_c_twoVertex_lTrack_Isolated <<
" " ;
3215 outputFile.width(7);
3216 outputFile.precision(3);
3217 outputFile << frac_c_threeVertex_lTrack_Primary <<
" " ;
3218 outputFile.width(7);
3219 outputFile.precision(3);
3220 outputFile << frac_c_threeVertex_lTrack_Secondary <<
" " ;
3221 outputFile.width(7);
3222 outputFile.precision(3);
3223 outputFile << frac_c_threeVertex_lTrack_Tertiary <<
" " ;
3224 outputFile.width(7);
3225 outputFile.precision(3);
3226 outputFile << frac_c_threeVertex_lTrack_Isolated <<
" " ;
3227 outputFile << std::endl;
3253 outputFile <<
" D decay ";
3254 outputFile.width(7);
3255 outputFile.precision(3);
3256 outputFile << frac_c_twoVertex_cTrack_Primary <<
" " ;
3257 outputFile.width(7);
3258 outputFile.precision(3);
3259 outputFile << frac_c_twoVertex_cTrack_Secondary <<
" " ;
3260 outputFile.width(7);
3261 outputFile.precision(3);
3262 outputFile << frac_c_twoVertex_cTrack_Isolated <<
" " ;
3263 outputFile.width(7);
3264 outputFile.precision(3);
3265 outputFile << frac_c_threeVertex_cTrack_Primary <<
" " ;
3266 outputFile.width(7);
3267 outputFile.precision(3);
3268 outputFile << frac_c_threeVertex_cTrack_Secondary <<
" " ;
3269 outputFile.width(7);
3270 outputFile.precision(3);
3271 outputFile << frac_c_threeVertex_cTrack_Tertiary <<
" " ;
3272 outputFile.width(7);
3273 outputFile.precision(3);
3274 outputFile << frac_c_threeVertex_cTrack_Isolated <<
" " ;
3275 outputFile << std::endl;
3276 outputFile <<
" ----------------------------------------------------------------------------" << std::endl;
3277 outputFile <<
" All above ";
3278 outputFile.width(7);
3279 outputFile.precision(3);
3280 outputFile << frac_c_twoVertex_Primary <<
" " ;
3281 outputFile.width(7);
3282 outputFile.precision(3);
3283 outputFile << frac_c_twoVertex_Secondary <<
" " ;
3284 outputFile.width(7);
3285 outputFile.precision(3);
3286 outputFile << frac_c_twoVertex_Isolated <<
" " ;
3287 outputFile.width(7);
3288 outputFile.precision(3);
3289 outputFile << frac_c_threeVertex_Primary <<
" " ;
3290 outputFile.width(7);
3291 outputFile.precision(3);
3292 outputFile << frac_c_threeVertex_Secondary <<
" " ;
3293 outputFile.width(7);
3294 outputFile.precision(3);
3295 outputFile << frac_c_threeVertex_Tertiary <<
" " ;
3296 outputFile.width(7);
3297 outputFile.precision(3);
3298 outputFile << frac_c_threeVertex_Isolated <<
" " ;
3299 outputFile << std::endl;
3300 outputFile <<
" ----------------------------------------------------------------------------" << std::endl;
3301 outputFile << std::endl;
3302 outputFile << std::endl;
3308 void LCFIAIDAPlotProcessor::PrintNNOutput(){
3311 std::filebuf* fb =
new std::filebuf;
3313 std::ostream outputFile( (!_PurityEfficiencyOutputFile.empty()) ?
3314 fb->open(_PurityEfficiencyOutputFile.c_str(),
3315 std::ios_base::out|std::ios_base::trunc):
3318 if (outputFile.rdbuf() != fb)
3321 std::cerr <<
" In " << __FILE__ <<
"(" << __LINE__ <<
"): Unable to open output file " << _PurityEfficiencyOutputFile <<
"! Redirecting output to standard out." << std::endl;
3322 outputFile << std::endl;
3325 for (
unsigned int iTagCollection=0; iTagCollection < _FlavourTagCollectionNames.size(); ++iTagCollection)
3327 outputFile <<
"\n\nRESULTS for " << iTagCollection <<
"-th Flavour Tag Collection named " << _FlavourTagCollectionNames[iTagCollection] << std::endl;
3328 outputFile <<
" any number of ZVTOP vertices found " ;
3329 outputFile <<
" N(b) = " ;
3330 outputFile.width(10);
3331 outputFile << _pBJetBTagIntegral[iTagCollection][_VertexCatNames[0]]->sumAllBinHeights();
3332 outputFile <<
" N(c) = " ;
3333 outputFile.width(10);
3334 outputFile << _pCJetBTagIntegral[iTagCollection][_VertexCatNames[0]]->sumAllBinHeights();
3335 outputFile <<
" N(light) = ";
3336 outputFile.width(10);
3337 outputFile << _pLightJetBTagIntegral[iTagCollection][_VertexCatNames[0]]->sumAllBinHeights();
3338 outputFile << std::endl << std::endl;
3340 for (
unsigned int iVertexCat=1; iVertexCat <= N_VERTEX_CATEGORIES; ++iVertexCat ) {
3342 std::string nvname = _VertexCatNames[iVertexCat];
3347 outputFile.width(2);
3348 outputFile << iVertexCat;
3350 (iVertexCat == 1) ?
" ZVTOP vertex found " :
" ZVTOP vertices found ") ;
3352 outputFile <<
" N(b) = " ;
3353 outputFile.width(10);
3354 outputFile << _pBJetBTagIntegral[iTagCollection][nvname]->sumAllBinHeights();
3355 outputFile <<
" N(c) = " ;
3356 outputFile.width(10);
3357 outputFile << _pCJetBTagIntegral[iTagCollection][nvname]->sumAllBinHeights();
3358 outputFile <<
" N(light) = ";
3359 outputFile.width(10);
3360 outputFile << _pLightJetBTagIntegral[iTagCollection][nvname]->sumAllBinHeights();
3361 outputFile << std::endl << std::endl;
3367 outputFile.precision(5);
3368 outputFile.setf(std::ios::fixed,std::ios::floatfield);
3371 for (
unsigned int iVertexCat=1; iVertexCat <= N_VERTEX_CATEGORIES; ++iVertexCat ) {
3373 std::string nvname = _VertexCatNames[iVertexCat];
3375 outputFile <<
"numbers of jets in cuts for ";
3376 outputFile.width(2);
3377 outputFile << iVertexCat;
3378 outputFile <<
" ZVTOP vertices found" << std::endl;
3379 outputFile <<
"cut b-tag b b-tag other c-tag c c-tag other c-tagbb c c-tagbb other" << std::endl;
3381 int numberOfBins=_pBJetBTagIntegral[iTagCollection][nvname]->axis().bins();
3383 for (
int binNumber = 0; binNumber < numberOfBins ; binNumber++ ) {
3385 outputFile.width(5);
3386 outputFile.precision(3);
3387 outputFile << _pBJetBTagIntegral[iTagCollection][nvname]->axis().binUpperEdge(binNumber) <<
" ";
3388 outputFile.width(10);
3389 outputFile << int(_pBJetBTagIntegral[iTagCollection][nvname]->binHeight(binNumber)) <<
" ";
3390 outputFile.width(10);
3391 outputFile << int(_pCJetBTagIntegral[iTagCollection][nvname]->binHeight(binNumber)+_pLightJetBTagIntegral[iTagCollection][nvname]->binHeight(binNumber))<<
" ";
3392 outputFile.width(10);
3393 outputFile << int(_pCJetCTagIntegral[iTagCollection][nvname]->binHeight(binNumber)) <<
" ";
3394 outputFile.width(10);
3395 outputFile << int(_pBJetCTagIntegral[iTagCollection][nvname]->binHeight(binNumber)+_pLightJetCTagIntegral[iTagCollection][nvname]->binHeight(binNumber))<<
" ";
3396 outputFile.width(10);
3397 outputFile << int(_pCJetBCTagIntegral[iTagCollection][nvname]->binHeight(binNumber)) <<
" ";
3398 outputFile.width(10);
3399 outputFile << int(_pBJetBCTagIntegral[iTagCollection][nvname]->binHeight(binNumber)+_pLightJetBCTagIntegral[iTagCollection][nvname]->binHeight(binNumber))<<
" ";
3400 outputFile << std::endl;
3403 outputFile << std::endl;
3406 outputFile <<
"numbers of jets in cuts summed (for any number of vertices)" << std::endl;
3407 outputFile <<
"cut b-tag b b-tag other c-tag c c-tag other c-tagbb c c-tagbb other" << std::endl;
3409 int numberOfBins=_pBJetBTagIntegral[iTagCollection][_VertexCatNames[0]]->axis().bins();
3411 for (
int binNumber = 0; binNumber < numberOfBins ; binNumber++ ) {
3413 outputFile.width(5);
3414 outputFile.precision(3);
3415 outputFile << _pBJetBTagIntegral[iTagCollection][_VertexCatNames[0]]->axis().binUpperEdge(binNumber) <<
" ";
3416 outputFile.width(10);
3417 outputFile << int(_pBJetBTagIntegral[iTagCollection][_VertexCatNames[0]]->binHeight(binNumber)) <<
" ";
3418 outputFile.width(10);
3419 outputFile << int(_pCJetBTagIntegral[iTagCollection][_VertexCatNames[0]]->binHeight(binNumber)+_pLightJetBTagIntegral[iTagCollection][_VertexCatNames[0]]->binHeight(binNumber))<<
" ";
3420 outputFile.width(10);
3421 outputFile << int(_pCJetCTagIntegral[iTagCollection][_VertexCatNames[0]]->binHeight(binNumber)) <<
" ";
3422 outputFile.width(10);
3423 outputFile << int(_pBJetCTagIntegral[iTagCollection][_VertexCatNames[0]]->binHeight(binNumber)+_pLightJetCTagIntegral[iTagCollection][_VertexCatNames[0]]->binHeight(binNumber))<<
" ";
3424 outputFile.width(10);
3425 outputFile << int(_pCJetBCTagIntegral[iTagCollection][_VertexCatNames[0]]->binHeight(binNumber)) <<
" ";
3426 outputFile.width(10);
3427 outputFile << int(_pBJetBCTagIntegral[iTagCollection][_VertexCatNames[0]]->binHeight(binNumber)+_pLightJetBCTagIntegral[iTagCollection][_VertexCatNames[0]]->binHeight(binNumber))<<
" ";
3428 outputFile << std::endl;
3430 outputFile << std::endl;
3431 outputFile << std::endl;
3435 for (
unsigned int iTagCollection=0; iTagCollection < _FlavourTagCollectionNames.size(); ++iTagCollection)
3438 int numberOfBins=_pBJetBTagIntegral[iTagCollection][_VertexCatNames[0]]->axis().bins();
3440 outputFile <<
"\n\nRESULTS for " << iTagCollection <<
"-th Flavour Tag Collection named " << _FlavourTagCollectionNames[iTagCollection] << std::endl;
3442 for (
unsigned int iVertexCat=0; iVertexCat <= N_VERTEX_CATEGORIES; ++iVertexCat ) {
3444 std::string nvname = _VertexCatNames[iVertexCat];
3446 std::vector<double> bJetBTagEfficiency;
3447 std::vector<double> bJetBTagEfficiencyError;
3448 std::vector<double> cJetCTagEfficiency;
3449 std::vector<double> cJetCTagEfficiencyError;
3450 std::vector<double> cJetBCTagEfficiency;
3451 std::vector<double> cJetBCTagEfficiencyError;
3453 CalculateTagEfficiency( _pBJetBTag[iTagCollection][nvname], bJetBTagEfficiency, bJetBTagEfficiencyError);
3454 CalculateTagEfficiency( _pCJetCTag[iTagCollection][nvname], cJetCTagEfficiency, cJetCTagEfficiencyError);
3455 CalculateTagEfficiency( _pCJetBCTag[iTagCollection][nvname], cJetBCTagEfficiency, cJetBCTagEfficiencyError);
3457 std::vector<double> bJetBTagPurity;
3458 std::vector<double> bJetBTagPurityError;
3459 std::vector<double> cJetCTagPurity;
3460 std::vector<double> cJetCTagPurityError;
3461 std::vector<double> cJetBCTagPurity;
3462 std::vector<double> cJetBCTagPurityError;
3464 CalculateTagPurity( _pBJetBTag[iTagCollection][nvname], _pBTagBackgroundValues[iTagCollection][nvname], bJetBTagPurity, bJetBTagPurityError);
3465 CalculateTagPurity( _pCJetCTag[iTagCollection][nvname], _pCTagBackgroundValues[iTagCollection][nvname], cJetCTagPurity, cJetCTagPurityError);
3466 CalculateTagPurity( _pCJetBCTag[iTagCollection][nvname], _pBJetBCTag[iTagCollection][nvname], cJetBCTagPurity, cJetBCTagPurityError);
3468 outputFile <<
"\n Purity-Efficiencies Values for ";
3469 if (iVertexCat==0) {
3470 outputFile <<
" any number of ZVTOP vertices found " << std::endl;
3472 outputFile << iVertexCat;
3474 (iVertexCat == 1) ?
" ZVTOP vertex found " :
" ZVTOP vertices found ") ;
3476 outputFile <<
"\n-----------------------------------------------------------------------------------------------" << std::endl;
3477 outputFile <<
"cut value efficiency(b) purity(b) efficiency(c) purity(c) efficiency(bc) purity(bc) " << std::endl;
3479 for (
int binNumber = 0; binNumber < numberOfBins ; binNumber++ ) {
3480 outputFile.width(5);
3481 outputFile.precision(3);
3482 outputFile << _pBJetBTagIntegral[iTagCollection][nvname]->axis().binUpperEdge(binNumber) <<
" ";
3484 outputFile << bJetBTagEfficiency[binNumber];
3485 outputFile <<
" +/- ";
3486 outputFile << bJetBTagEfficiencyError[binNumber];
3488 outputFile << bJetBTagPurity[binNumber];
3489 outputFile <<
" +/- ";
3490 outputFile << bJetBTagPurityError[binNumber];
3493 outputFile << cJetCTagEfficiency[binNumber];
3494 outputFile <<
" +/- ";
3495 outputFile << cJetCTagEfficiencyError[binNumber];
3497 outputFile << cJetCTagPurity[binNumber];
3498 outputFile <<
" +/- ";
3499 outputFile << cJetCTagPurityError[binNumber];
3502 outputFile << cJetBCTagEfficiency[binNumber];
3503 outputFile <<
" +/- ";
3504 outputFile << cJetBCTagEfficiencyError[binNumber];
3506 outputFile << cJetBCTagPurity[binNumber];
3507 outputFile <<
" +/- ";
3508 outputFile << cJetBCTagPurityError[binNumber];
3511 outputFile << std::endl;
3515 outputFile << std::endl;
3516 outputFile << std::endl;
3522 void LCFIAIDAPlotProcessor::InternalVectorInitialisation()
3526 _cJet_truePlus2_angle.resize(N_JETANGLE_BINS);
3527 _cJet_truePlus_angle.resize(N_JETANGLE_BINS);
3528 _cJet_trueNeut_angle.resize(N_JETANGLE_BINS);
3529 _cJet_trueMinus_angle.resize(N_JETANGLE_BINS);
3530 _cJet_trueMinus2_angle.resize(N_JETANGLE_BINS);
3532 _cJet_truePlus2_recoPlus_angle.resize(N_JETANGLE_BINS);
3533 _cJet_truePlus2_recoNeut_angle.resize(N_JETANGLE_BINS);
3534 _cJet_truePlus2_recoMinus_angle.resize(N_JETANGLE_BINS);
3535 _cJet_truePlus_recoPlus_angle.resize(N_JETANGLE_BINS);
3536 _cJet_truePlus_recoNeut_angle.resize(N_JETANGLE_BINS);
3537 _cJet_truePlus_recoMinus_angle.resize(N_JETANGLE_BINS);
3538 _cJet_trueNeut_recoPlus_angle.resize(N_JETANGLE_BINS);
3539 _cJet_trueNeut_recoNeut_angle.resize(N_JETANGLE_BINS);
3540 _cJet_trueNeut_recoMinus_angle.resize(N_JETANGLE_BINS);
3541 _cJet_trueMinus_recoPlus_angle.resize(N_JETANGLE_BINS);
3542 _cJet_trueMinus_recoNeut_angle.resize(N_JETANGLE_BINS);
3543 _cJet_trueMinus_recoMinus_angle.resize(N_JETANGLE_BINS);
3544 _cJet_trueMinus2_recoPlus_angle.resize(N_JETANGLE_BINS);
3545 _cJet_trueMinus2_recoNeut_angle.resize(N_JETANGLE_BINS);
3546 _cJet_trueMinus2_recoMinus_angle.resize(N_JETANGLE_BINS);
3548 _bJet_truePlus2_angle.resize(N_JETANGLE_BINS);
3549 _bJet_truePlus_angle.resize(N_JETANGLE_BINS);
3550 _bJet_trueNeut_angle.resize(N_JETANGLE_BINS);
3551 _bJet_trueMinus_angle.resize(N_JETANGLE_BINS);
3552 _bJet_trueMinus2_angle.resize(N_JETANGLE_BINS);
3553 _bJet_truePlus2_recoPlus_angle.resize(N_JETANGLE_BINS);
3554 _bJet_truePlus2_recoNeut_angle.resize(N_JETANGLE_BINS);
3555 _bJet_truePlus2_recoMinus_angle.resize(N_JETANGLE_BINS);
3556 _bJet_truePlus_recoPlus_angle.resize(N_JETANGLE_BINS);
3557 _bJet_truePlus_recoNeut_angle.resize(N_JETANGLE_BINS);
3558 _bJet_truePlus_recoMinus_angle.resize(N_JETANGLE_BINS);
3559 _bJet_trueNeut_recoPlus_angle.resize(N_JETANGLE_BINS);
3560 _bJet_trueNeut_recoNeut_angle.resize(N_JETANGLE_BINS);
3561 _bJet_trueNeut_recoMinus_angle.resize(N_JETANGLE_BINS);
3562 _bJet_trueMinus_recoPlus_angle.resize(N_JETANGLE_BINS);
3563 _bJet_trueMinus_recoNeut_angle.resize(N_JETANGLE_BINS);
3564 _bJet_trueMinus_recoMinus_angle.resize(N_JETANGLE_BINS);
3565 _bJet_trueMinus2_recoPlus_angle.resize(N_JETANGLE_BINS);
3566 _bJet_trueMinus2_recoNeut_angle.resize(N_JETANGLE_BINS);
3567 _bJet_trueMinus2_recoMinus_angle.resize(N_JETANGLE_BINS);
3570 for (
int j = 0 ; j < N_JETANGLE_BINS ; j++) {
3571 _cJet_truePlus2_recoPlus_angle[j]=0;
3572 _cJet_truePlus2_recoNeut_angle[j]=0;
3573 _cJet_truePlus2_recoMinus_angle[j]=0;
3574 _cJet_truePlus_recoPlus_angle[j]=0;
3575 _cJet_truePlus_recoNeut_angle[j]=0;
3576 _cJet_truePlus_recoMinus_angle[j]=0;
3577 _cJet_trueNeut_recoPlus_angle[j]=0;
3578 _cJet_trueNeut_recoNeut_angle[j]=0;
3579 _cJet_trueNeut_recoMinus_angle[j]=0;
3580 _cJet_trueMinus_recoPlus_angle[j]=0;
3581 _cJet_trueMinus_recoNeut_angle[j]=0;
3582 _cJet_trueMinus_recoMinus_angle[j]=0;
3583 _cJet_trueMinus2_recoPlus_angle[j]=0;
3584 _cJet_trueMinus2_recoNeut_angle[j]=0;
3585 _cJet_trueMinus2_recoMinus_angle[j]=0;
3587 _bJet_truePlus2_angle[j]=0;
3588 _bJet_truePlus_angle[j]=0;
3589 _bJet_trueNeut_angle[j]=0;
3590 _bJet_trueMinus_angle[j]=0;
3591 _bJet_trueMinus2_angle[j]=0;
3592 _bJet_truePlus2_recoPlus_angle[j]=0;
3593 _bJet_truePlus2_recoNeut_angle[j]=0;
3594 _bJet_truePlus2_recoMinus_angle[j]=0;
3595 _bJet_truePlus_recoPlus_angle[j]=0;
3596 _bJet_truePlus_recoNeut_angle[j]=0;
3597 _bJet_truePlus_recoMinus_angle[j]=0;
3598 _bJet_trueNeut_recoPlus_angle[j]=0;
3599 _bJet_trueNeut_recoNeut_angle[j]=0;
3600 _bJet_trueNeut_recoMinus_angle[j]=0;
3601 _bJet_trueMinus_recoPlus_angle[j]=0;
3602 _bJet_trueMinus_recoNeut_angle[j]=0;
3603 _bJet_trueMinus_recoMinus_angle[j]=0;
3604 _bJet_trueMinus2_recoPlus_angle[j]=0;
3605 _bJet_trueMinus2_recoNeut_angle[j]=0;
3606 _bJet_trueMinus2_recoMinus_angle[j]=0;
3615 _cJet_truePlus2_recoPlus=0;
3616 _cJet_truePlus2_recoNeut=0;
3617 _cJet_truePlus2_recoMinus=0;
3618 _cJet_truePlus_recoPlus=0;
3619 _cJet_truePlus_recoNeut=0;
3620 _cJet_truePlus_recoMinus=0;
3621 _cJet_trueNeut_recoPlus=0;
3622 _cJet_trueNeut_recoNeut=0;
3623 _cJet_trueNeut_recoMinus=0;
3624 _cJet_trueMinus_recoPlus=0;
3625 _cJet_trueMinus_recoNeut=0;
3626 _cJet_trueMinus_recoMinus=0;
3627 _cJet_trueMinus2_recoPlus=0;
3628 _cJet_trueMinus2_recoNeut=0;
3629 _cJet_trueMinus2_recoMinus=0;
3636 _bJet_truePlus2_recoPlus=0;
3637 _bJet_truePlus2_recoNeut=0;
3638 _bJet_truePlus2_recoMinus=0;
3639 _bJet_truePlus_recoPlus=0;
3640 _bJet_truePlus_recoNeut=0;
3641 _bJet_truePlus_recoMinus=0;
3642 _bJet_trueNeut_recoPlus=0;
3643 _bJet_trueNeut_recoNeut=0;
3644 _bJet_trueNeut_recoMinus=0;
3645 _bJet_trueMinus_recoPlus=0;
3646 _bJet_trueMinus_recoNeut=0;
3647 _bJet_trueMinus_recoMinus=0;
3648 _bJet_trueMinus2_recoPlus=0;
3649 _bJet_trueMinus2_recoNeut=0;
3650 _bJet_trueMinus2_recoMinus=0;
3652 _nb_twoVertex_bTrack_Primary=0;
3653 _nb_twoVertex_bTrack_Secondary=0;
3654 _nb_twoVertex_bTrack_Tertiary=0;
3655 _nb_twoVertex_bTrack_Isolated=0;
3656 _nb_twoVertex_cTrack_Primary=0;
3657 _nb_twoVertex_cTrack_Secondary=0;
3658 _nb_twoVertex_cTrack_Tertiary=0;
3659 _nb_twoVertex_cTrack_Isolated=0;
3660 _nb_twoVertex_lTrack_Primary=0;
3661 _nb_twoVertex_lTrack_Secondary=0;
3662 _nb_twoVertex_lTrack_Tertiary=0;
3663 _nb_twoVertex_lTrack_Isolated=0;
3664 _nb_threeVertex_bTrack_Primary=0;
3665 _nb_threeVertex_bTrack_Secondary=0;
3666 _nb_threeVertex_bTrack_Tertiary=0;
3667 _nb_threeVertex_bTrack_Isolated=0;
3668 _nb_threeVertex_cTrack_Primary=0;
3669 _nb_threeVertex_cTrack_Secondary=0;
3670 _nb_threeVertex_cTrack_Tertiary=0;
3671 _nb_threeVertex_cTrack_Isolated=0;
3672 _nb_threeVertex_lTrack_Primary=0;
3673 _nb_threeVertex_lTrack_Secondary=0;
3674 _nb_threeVertex_lTrack_Tertiary=0;
3675 _nb_threeVertex_lTrack_Isolated=0;
3677 _nb_threeVertex_Primary_noMCP=0;
3678 _nb_threeVertex_Secondary_noMCP=0;
3679 _nb_threeVertex_Tertiary_noMCP=0;
3680 _nb_threeVertex_Isolated_noMCP=0;
3681 _nb_twoVertex_Primary_noMCP=0;
3682 _nb_twoVertex_Secondary_noMCP=0;
3683 _nb_twoVertex_Tertiary_noMCP=0;
3684 _nb_twoVertex_Isolated_noMCP=0;
3686 _nc_twoVertex_bTrack_Primary=0;
3687 _nc_twoVertex_bTrack_Secondary=0;
3688 _nc_twoVertex_bTrack_Tertiary=0;
3689 _nc_twoVertex_bTrack_Isolated=0;
3690 _nc_twoVertex_cTrack_Primary=0;
3691 _nc_twoVertex_cTrack_Secondary=0;
3692 _nc_twoVertex_cTrack_Tertiary=0;
3693 _nc_twoVertex_cTrack_Isolated=0;
3694 _nc_twoVertex_lTrack_Primary=0;
3695 _nc_twoVertex_lTrack_Secondary=0;
3696 _nc_twoVertex_lTrack_Tertiary=0;
3697 _nc_twoVertex_lTrack_Isolated=0;
3698 _nc_threeVertex_bTrack_Primary=0;
3699 _nc_threeVertex_bTrack_Secondary=0;
3700 _nc_threeVertex_bTrack_Tertiary=0;
3701 _nc_threeVertex_bTrack_Isolated=0;
3702 _nc_threeVertex_cTrack_Primary=0;
3703 _nc_threeVertex_cTrack_Secondary=0;
3704 _nc_threeVertex_cTrack_Tertiary=0;
3705 _nc_threeVertex_cTrack_Isolated=0;
3706 _nc_threeVertex_lTrack_Primary=0;
3707 _nc_threeVertex_lTrack_Secondary=0;
3708 _nc_threeVertex_lTrack_Tertiary=0;
3709 _nc_threeVertex_lTrack_Isolated=0;
3711 _nc_threeVertex_Primary_noMCP=0;
3712 _nc_threeVertex_Secondary_noMCP=0;
3713 _nc_threeVertex_Tertiary_noMCP=0;
3714 _nc_threeVertex_Isolated_noMCP=0;
3715 _nc_twoVertex_Primary_noMCP=0;
3716 _nc_twoVertex_Secondary_noMCP=0;
3717 _nc_twoVertex_Tertiary_noMCP=0;
3718 _nc_twoVertex_Isolated_noMCP=0;
3726 std::vector<double>& bjetdecaylengthvector, std::vector<double>& cjetdecaylengthvector)
3732 int jettype = FindTrueJetType( pEvent, jetNumber );
3734 if (abs(jettype)!= C_JET && abs(jettype)!=B_JET)
3739 int pdgcode = FindTrueJetPDGCode( pEvent, jetNumber );
3743 if (!mcParticleCollection.is_valid())
3745 std::cerr <<
"In " << __FILE__ <<
"(" << __LINE__ <<
"): for event " << pEvent->getEventNumber() <<
" in run " << pEvent->getRunNumber() <<
"mcParticleCollection is not valid" << std::endl;
3750 for (
int imc = 0 ; imc < mcParticleCollection.getNumberOfElements(); imc++ ) {
3752 lcio::MCParticle* myMCParticle = mcParticleCollection.getElementAt( imc );
3754 if (myMCParticle->getPDG() == pdgcode) {
3759 const double* vertex = myMCParticle->getVertex();
3760 const double* endpoint = myMCParticle->getEndpoint();
3761 double decay_length = CalculateDistance(vertex, endpoint);
3765 decaylengthvector.push_back(decay_length);
3766 if (GetPDGFlavour(pdgcode) == C_JET ) cjetdecaylengthvector.push_back(decay_length);
3767 if (GetPDGFlavour(pdgcode) == B_JET ) bjetdecaylengthvector.push_back(decay_length);
3770 const MCParticleVec& daughters = myMCParticle->getDaughters();
3772 for ( MCParticleVec::const_iterator idaughter = daughters.begin(); idaughter!=daughters.end(); idaughter++)
3774 int code = (*idaughter)->getPDG();
3778 int flavour = GetPDGFlavour(code);
3780 if (abs(flavour) == C_JET || abs(flavour) == B_JET) {
3786 const double* vertex = (*idaughter)->getVertex();
3787 const double* endpoint = (*idaughter)->getEndpoint();
3788 double decay_length = CalculateDistance(vertex, endpoint);
3791 decaylengthvector.push_back(decay_length);
3792 if (abs(flavour) == C_JET ) cjetdecaylengthvector.push_back(decay_length);
3793 if (abs(flavour) == B_JET ) bjetdecaylengthvector.push_back(decay_length);
3803 const MCParticleVec& daughters2 = (*idaughter)->getDaughters();
3805 for ( MCParticleVec::const_iterator idaughter2 = daughters2.begin(); idaughter2!=daughters2.end(); idaughter2++)
3807 int code2 = (*idaughter2)->getPDG();
3810 int flavour2 = GetPDGFlavour(code2);
3812 if (abs(flavour2) == C_JET || abs (flavour2) == B_JET) {
3817 const double* vertex = (*idaughter2)->getVertex();
3818 const double* endpoint = (*idaughter2)->getEndpoint();
3819 double decay_length = CalculateDistance(vertex, endpoint);
3822 decaylengthvector.push_back(decay_length);
3823 if (abs(flavour2) == C_JET ) cjetdecaylengthvector.push_back(decay_length);
3824 if (abs(flavour2) == B_JET ) bjetdecaylengthvector.push_back(decay_length);
3827 const MCParticleVec& daughters3 = (*idaughter2)->getDaughters();
3829 for ( MCParticleVec::const_iterator idaughter3 = daughters3.begin(); idaughter3!=daughters3.end(); idaughter3++)
3831 int code3 = (*idaughter3)->getPDG();
3834 int flavour3 = GetPDGFlavour(code3);
3836 if (abs(flavour3)==C_JET || abs(flavour3)==B_JET) {
3841 const double* vertex = (*idaughter3)->getVertex();
3842 const double* endpoint = (*idaughter3)->getEndpoint();
3843 double decay_length = CalculateDistance(vertex, endpoint);
3846 decaylengthvector.push_back(decay_length);
3847 if (abs(flavour3) == C_JET ) cjetdecaylengthvector.push_back(decay_length);
3848 if (abs(flavour3) == B_JET ) bjetdecaylengthvector.push_back(decay_length);
3851 const MCParticleVec& daughters4 = (*idaughter3)->getDaughters();
3853 for ( MCParticleVec::const_iterator idaughter4 = daughters4.begin(); idaughter4!=daughters4.end(); idaughter4++)
3855 int code4 = (*idaughter4)->getPDG();
3858 int flavour4 = GetPDGFlavour(code4);
3860 if (abs(flavour4)==C_JET || abs(flavour4)==B_JET) {
3865 const double* vertex = (*idaughter4)->getVertex();
3866 const double* endpoint = (*idaughter4)->getEndpoint();
3867 double decay_length = CalculateDistance(vertex, endpoint);
3869 decaylengthvector.push_back(decay_length);
3870 if (abs(flavour4) == C_JET ) cjetdecaylengthvector.push_back(decay_length);
3871 if (abs(flavour4) == B_JET ) bjetdecaylengthvector.push_back(decay_length);
3895 void LCFIAIDAPlotProcessor::FindTrueJetDecayLength2( LCEvent* pEvent,
unsigned int jetNumber,
double& bjetdecaylength,
double& cjetdecaylength)
3897 double b_decay_length0 = -1.;
3898 double c_decay_length0 = -1.;
3902 int jettype = FindTrueJetType( pEvent, jetNumber );
3904 if (abs(jettype)!=B_JET && abs(jettype)!=C_JET )
3909 int pdgcode = FindTrueJetPDGCode( pEvent, jetNumber );
3913 if (!mcParticleCollection.is_valid())
3915 std::cerr << __FILE__ <<
"(" << __LINE__ <<
"): MCParticle collection " << _MCParticleColName <<
" is not valid"
3916 <<
" for event " << pEvent->getEventNumber() <<
" in run " << pEvent->getRunNumber() << std::endl;
3921 for (
int imc = 0 ; imc < mcParticleCollection.getNumberOfElements(); imc++ ) {
3923 lcio::MCParticle* myMCParticle = mcParticleCollection.getElementAt( imc );
3924 if (myMCParticle->getPDG() == pdgcode) {
3927 std::vector<double> BJetDecayLengthVector;
3928 std::vector<double> CJetDecayLengthVector;
3929 std::vector<int> BJetDecayVector;
3930 std::vector<int> CJetDecayVector;
3933 int code0 = pdgcode;
3934 int flavour0 = GetPDGFlavour(code0);
3938 const double* vertex0 = myMCParticle->getVertex();
3941 const MCParticleVec& daughters = myMCParticle->getDaughters();
3943 for ( MCParticleVec::const_iterator idaughter = daughters.begin(); idaughter!=daughters.end(); idaughter++)
3945 int code1 = (*idaughter)->getPDG();
3946 int flavour1 = GetPDGFlavour(code1);
3948 if (abs(flavour1) != B_JET && abs(flavour0) == B_JET) {
3950 const double* vertex1 = (*idaughter)->getVertex();
3951 double decay_length = CalculateDistance(vertex0, vertex1);
3952 BJetDecayLengthVector.push_back(decay_length);
3955 if (abs(flavour1) != C_JET && abs(flavour0) == C_JET) {
3956 const double* vertex1 = (*idaughter)->getVertex();
3957 double decay_length = CalculateDistance(vertex0, vertex1);
3958 CJetDecayLengthVector.push_back(decay_length);
3963 const MCParticleVec& daughters2 = (*idaughter)->getDaughters();
3965 for ( MCParticleVec::const_iterator idaughter2 = daughters2.begin(); idaughter2!=daughters2.end(); idaughter2++)
3967 int code2 = (*idaughter2)->getPDG();
3968 int flavour2 = GetPDGFlavour(code2);
3970 if (abs(flavour2) != B_JET && abs(flavour1) == B_JET) {
3971 const double* vertex2 = (*idaughter2)->getVertex();
3972 double decay_length = CalculateDistance(vertex0, vertex2);
3973 BJetDecayLengthVector.push_back(decay_length);
3976 if (abs(flavour2) != C_JET && abs(flavour1) == C_JET) {
3977 const double* vertex2 = (*idaughter2)->getVertex();
3978 double decay_length = CalculateDistance(vertex0, vertex2);
3979 CJetDecayLengthVector.push_back(decay_length);
3982 const MCParticleVec& daughters3 = (*idaughter2)->getDaughters();
3984 for ( MCParticleVec::const_iterator idaughter3 = daughters3.begin(); idaughter3!=daughters3.end(); idaughter3++)
3986 int code3 = (*idaughter3)->getPDG();
3987 int flavour3 = GetPDGFlavour(code3);
3989 if (abs(flavour3) != B_JET && abs(flavour2) == B_JET) {
3990 const double* vertex3 = (*idaughter3)->getVertex();
3991 double decay_length = CalculateDistance(vertex0, vertex3);
3992 BJetDecayLengthVector.push_back(decay_length);
3995 if (abs(flavour3) != C_JET && abs(flavour2) == C_JET) {
3996 const double* vertex3 = (*idaughter3)->getVertex();
3997 double decay_length = CalculateDistance(vertex0, vertex3);
3998 CJetDecayLengthVector.push_back(decay_length);
4001 const MCParticleVec& daughters4 = (*idaughter3)->getDaughters();
4003 for ( MCParticleVec::const_iterator idaughter4 = daughters4.begin(); idaughter4!=daughters4.end(); idaughter4++)
4005 int code4 = (*idaughter4)->getPDG();
4006 int flavour4 = GetPDGFlavour(code4);
4008 if (abs(flavour4) != B_JET && abs(flavour3) == B_JET) {
4009 const double* vertex4 = (*idaughter4)->getVertex();
4010 double decay_length = CalculateDistance(vertex0, vertex4);
4011 BJetDecayLengthVector.push_back(decay_length);
4013 if (abs(flavour4) != C_JET && abs(flavour3) == C_JET) {
4014 const double* vertex4 = (*idaughter4)->getVertex();
4015 double decay_length = CalculateDistance(vertex0, vertex4);
4016 CJetDecayLengthVector.push_back(decay_length);
4020 const MCParticleVec& daughters5 = (*idaughter4)->getDaughters();
4022 for ( MCParticleVec::const_iterator idaughter5 = daughters5.begin(); idaughter5!=daughters5.end(); idaughter5++)
4024 int code5 = (*idaughter5)->getPDG();
4025 int flavour5 = GetPDGFlavour(code5);
4027 if (abs(flavour5) != B_JET && abs(flavour4) == B_JET) {
4028 const double* vertex5 = (*idaughter5)->getVertex();
4029 double decay_length = CalculateDistance(vertex0, vertex5);
4030 BJetDecayLengthVector.push_back(decay_length);
4032 if (abs(flavour5) != C_JET && abs(flavour4) == C_JET) {
4033 const double* vertex5 = (*idaughter5)->getVertex();
4034 double decay_length = CalculateDistance(vertex0, vertex5);
4035 CJetDecayLengthVector.push_back(decay_length);
4046 for (std::vector<double>::const_iterator iter = BJetDecayLengthVector.begin(); iter != BJetDecayLengthVector.end() ; iter++)
4048 double decay_length = *iter;
4049 if (decay_length > b_decay_length0 ) b_decay_length0 = decay_length;
4052 for (std::vector<double>::const_iterator iter = CJetDecayLengthVector.begin(); iter != CJetDecayLengthVector.end() ; iter++)
4054 double decay_length = *iter;
4055 if (decay_length > c_decay_length0 ) c_decay_length0 = decay_length;
4063 bjetdecaylength = b_decay_length0;
4064 cjetdecaylength = c_decay_length0;
4070 int LCFIAIDAPlotProcessor::GetPDGFlavour(
int code){
4078 if (code==411)
return C_JET;
4079 if (code==421)
return C_JET;
4080 if (code==431)
return C_JET;
4082 if (code==511)
return B_JET;
4083 if (code==521)
return B_JET;
4084 if (code==531)
return B_JET;
4085 if (code==541)
return B_JET;
4087 if (code==4122)
return C_JET;
4088 if (code==4222)
return C_JET;
4089 if (code==4112)
return C_JET;
4090 if (code==4212)
return C_JET;
4091 if (code==4232)
return C_JET;
4092 if (code==4132)
return C_JET;
4095 if (code==5122)
return B_JET;
4096 if (code==5222)
return B_JET;
4097 if (code==5112)
return B_JET;
4098 if (code==5212)
return B_JET;
4099 if (code==5132)
return B_JET;
4100 if (code==5232)
return B_JET;
4101 if (code==5332)
return B_JET;
4114 if( code/ 1000 == C_JET )
4119 if( code/ 1000 == B_JET )
4127 if( code / 100 == C_JET )
4133 if( code/100 == B_JET )
4144 #endif // endif of the check to see if MARLIN_USE_AIDA has been defined
int FindCQVtx(LCEvent *pEvent, unsigned int jetNumber)
Finds the vertex charge of the jet - using cuts tuned to find vertex charge for C-jets (from CVertexC...
AIDA::IDataPointSet * CreateIntegralPlot(const AIDA::IHistogram1D *pNN, AIDA::IDataPointSet *pIntegral)
Makes a DataPointSet integrating a histogram from the first bin to the last bin – NOT USED...
void FindTrueJetDecayLength(LCEvent *pEvent, unsigned int jetNumber, std::vector< double > &decaylengthvector, std::vector< double > &bjetdecaylengthvector, std::vector< double > &cjetdecaylengthvector)
Finds the true decay length of the longest b- or c- hadron in a jet.
AIDA::IDataPointSet * CreateXYPlot(const AIDA::IDataPointSet *pDataPointSet0, const AIDA::IDataPointSet *pDataPointSet1, AIDA::IDataPointSet *xyPointSet, const int dim0=0, const int dim1=0)
Plots two DataPointSets against each other.
int FindTrueJetFlavour(LCEvent *pEvent, unsigned int jetNumber)
Finds the true flavour of the jet (uses TrueJetFlavourCollection)
AIDA::IDataPointSet * CreateEfficiencyPlot(const AIDA::IHistogram1D *pSignal, AIDA::IDataPointSet *pDataPointSet)
Makes a DataPointSet of the tag efficiency e.g number of B-jets passing a given B-tag NN cut...
int FindNumVertex(LCEvent *pEvent, unsigned int jetNumber, unsigned int iInputsCollection)
Finds the number of vertices in an event (from the flavour tag inputs)
AIDA::IDataPointSet * CreatePurityPlot(const AIDA::IHistogram1D *pSignal, const AIDA::IHistogram1D *pBackground, AIDA::IDataPointSet *pDataPointSet)
Makes a DataPointSet of the tag purity e.g. N(B-jets passing NN cut)/N(all-jets passing NN cut) for a...
int FindBQVtx(LCEvent *pEvent, unsigned int jetNumber)
Finds the vertex charge of the jet - using cuts tuned to find vertex charge for B-jets (from BVertexC...
float FindTrueJetPartonCharge(LCEvent *pEvent, unsigned int jetNumber)
Finds the true charge of the parton producing a jet (uses TrueJetFlavourCollection) ...
AIDA::IDataPointSet * CreateEfficiencyPlot2(const AIDA::IHistogram1D *pAllEvents, const AIDA::IHistogram1D *pPassEvents, AIDA::IDataPointSet *pDataPointSet)
Makes a DataPointSet of histogram 1 divide by histogram 2 - this is an IDataPointSet as a histrogram ...
AIDA::IDataPointSet * CreateLeakageRatePlot(const AIDA::IHistogram1D *pBackground, AIDA::IDataPointSet *pDataPointSet)
Makes a DataPointSet showing the tagging leakage e.g. the number of non-B-jets passing a given B-tag ...
AIDA::IHistogram1D * CreateIntegralHistogram(const AIDA::IHistogram1D *pNN, AIDA::IHistogram1D *pIntegral)
Makes a histogram integrating a histogram from the first bin to the last bin - THE ERRORS RETURNED AR...
int FindTrueJetPDGCode(LCEvent *pEvent, unsigned int jetNumber)
Finds the PDG code of the hadron producing a jet (uses TrueJetFlavourCollection)
int FindTrueJetType(LCEvent *pEvent, unsigned int jetNumber)
Finds the true flavour of a jet (uses TrueJetFlavourCollection)
LCFIAIDAPlotProcessor Class - make plots of the LCFI flavour tag and vertex charge code...
void CreateVertexChargeLeakagePlot(AIDA::IDataPointSet *pBJetVtxChargeDPS, AIDA::IDataPointSet *pCJetVtxChargeDPS)
Makes DataPointSets for the number of.
float FindTrueJetHadronCharge(LCEvent *pEvent, unsigned int jetNumber)
Finds the true charge of the hadron producing a jet (uses TrueJetFlavourCollection) ...