1 #include "TrueAngularJetFlavourProcessor.h"
9 #include <EVENT/LCParameters.h>
10 #include <EVENT/LCCollection.h>
11 #include <EVENT/ReconstructedParticle.h>
12 #include <EVENT/Track.h>
13 #include <EVENT/LCRelation.h>
14 #include <EVENT/MCParticle.h>
15 #include <EVENT/LCIntVec.h>
16 #include "EVENT/LCFloatVec.h"
18 #include <UTIL/LCRelationNavigator.h>
19 #include <IMPL/ParticleIDImpl.h>
20 #include <IMPL/LCCollectionVec.h>
30 using namespace marlin ;
39 TrueAngularJetFlavourProcessor::TrueAngularJetFlavourProcessor() : Processor(
"TrueAngularJetFlavourProcessor") {
42 _description =
"TrueAngularJetFlavourProcessor - Determines the true flavour of a jet from the MC Paticles associated to the Jets RP" ;
47 registerInputCollection( lcio::LCIO::RECONSTRUCTEDPARTICLE,
49 "Name of the ReconstructedParticle collection that represents jets" ,
51 std::string(
"Jets") ) ;
52 registerInputCollection( lcio::LCIO::MCPARTICLE,
53 "MCParticleCollection" ,
54 "Name of the collection that holds all MC particles. " ,
56 std::string(
"MCParticle") ) ;
57 registerOutputCollection( lcio::LCIO::LCFLOATVEC,
58 "TrueJetFlavourCollection" ,
59 "Name of the output collection of LCIntVec (same order as jets)" ,
60 _TrueJetFlavourColName ,
61 std::string(
"TrueJetFlavour") ) ;
62 registerOptionalParameter(
"MaximumAngle" ,
63 "Maximum value allowed between MCParticle and jet momentum expressed in degrees" ,
68 void TrueAngularJetFlavourProcessor::init()
77 void TrueAngularJetFlavourProcessor::processRunHeader( LCRunHeader* run) {
78 _TrueJetVariableNames.clear();
79 _TrueJetVariableNames.push_back(
"TrueJetFlavour");
80 _TrueJetVariableNames.push_back(
"TruePartonCharge");
81 _TrueJetVariableNames.push_back(
"TrueHadronCharge");
82 _TrueJetVariableNames.push_back(
"TruePDGCode");
84 run->parameters().setValues(_TrueJetFlavourColName, _TrueJetVariableNames);
89 void TrueAngularJetFlavourProcessor::processEvent( LCEvent * evt ) {
91 LCCollection* JetRPCol = evt->getCollection( _JetRPColName );
92 LCCollection* MCParticleCol = evt->getCollection( _MCParticleColName );
94 std::vector<ReconstructedParticle*> MyJets;
95 std::vector<int> MCflavour;
96 std::vector<float> MCCharge;
97 std::vector<float> MCPCharge;
98 std::vector<int> MCCode;
101 LCCollectionVec* OutCollection =
new LCCollectionVec(
"LCFloatVec");
103 evt->addCollection(OutCollection,_TrueJetFlavourColName);
105 int nRCP = JetRPCol->getNumberOfElements() ;
107 if(nRCP ==0 ) std::cerr<<
"Warning: TrueAngularFetFlavourProcessor.cc:88 : NO jets present "<<std::endl;
109 for(
int i=0; i< nRCP ; i++)
111 MyJets.push_back(dynamic_cast<ReconstructedParticle*>(JetRPCol->getElementAt(i)));
114 vector<int> jetflavourdata(MyJets.size(),1);
115 vector<int> jetcodedata(MyJets.size(),0);
116 vector<double> jetangledata(MyJets.size(),1000);
119 vector<float> hadronchargedata(MyJets.size(),-100);
120 vector<float> partonchargedata(MyJets.size(),-100);
122 std::vector<std::vector<int> > Ntrkjet;
123 bool originalhadron = 0;
129 std::vector<MCParticle*> HeavyMCs;
130 std::vector<MCParticle*> JetMCs;
131 std::vector<double> TrueFlavour;
133 std::vector<int> typepart;
134 std::vector<MCParticle*> daughtercharge;
137 if( MCParticleCol->getNumberOfElements() == 0 ) std::cerr<<
"Warning: TrueAngularFetFlavourProcessor.cc:107 : NO MC data presentjets present "<<std::endl;
140 for(
int nMCpart = 0; nMCpart< MCParticleCol->getNumberOfElements();nMCpart++ )
144 MCParticle* MCused =
dynamic_cast<MCParticle*
> (MCParticleCol->getElementAt(nMCpart));
148 code = abs( MCused->getPDG() );
149 charge = MCused->getCharge();
153 if(charge == -1000 && code >100)
155 charge = chargefromPDG(MCused->getPDG());
173 if( code/ 1000 == 4 )
175 HeavyMCs.push_back(MCused);
176 MCflavour.push_back(4);
177 MCCharge.push_back(charge);
178 MCCode.push_back(MCused->getPDG());
179 if( MCused->getPDG() > 0 )
183 MCPCharge.push_back(pcharge);
186 if( code/ 1000 == 5 )
188 HeavyMCs.push_back(MCused);
189 MCflavour.push_back(5);
190 MCCharge.push_back(charge);
191 MCCode.push_back(MCused->getPDG());
192 if( MCused->getPDG() > 0 )
196 MCPCharge.push_back(pcharge);
202 if( code / 100 == 4 )
204 HeavyMCs.push_back(MCused);
205 MCflavour.push_back(4);
206 MCCharge.push_back(charge);
207 MCCode.push_back(MCused->getPDG());
209 if ( (code/10) %11 == 0)
215 if( MCused->getPDG() > 0 )
220 MCPCharge.push_back(pcharge);
225 HeavyMCs.push_back(MCused);
226 MCflavour.push_back(5);
227 MCCharge.push_back(charge);
228 MCCode.push_back(MCused->getPDG());
236 if( MCused->getPDG() < 0 )
241 MCPCharge.push_back(pcharge);
249 for (std::vector<MCParticle*>::const_iterator iParticle = HeavyMCs.begin(); iParticle != HeavyMCs.end() ;++iParticle)
252 MCParticle* MCused = *iParticle;
253 while(!originalhadron)
258 if( MCused->getParents().size() == 1)
261 if (abs( MCused->getParents()[0]->getPDG()) >100 || (abs(MCused->getParents()[0]->getPDG())>10 && abs(MCused->getParents()[0]->getPDG())<81))
263 MCused = MCused->getParents()[0];
267 JetMCs.push_back(MCused);
273 JetMCs.push_back(MCused);
286 for(std::vector<MCParticle*>::const_iterator iParticle3 = JetMCs.begin(); iParticle3 != JetMCs.end() ;++iParticle3)
288 double dist = -999999999;
292 double anglestored = 100000;
293 vector<double> normMomMC;
294 const double* MomentumMC = (*iParticle3)->getMomentum();
295 double length = sqrt(MomentumMC[0]*MomentumMC[0]+MomentumMC[1]*MomentumMC[1]+MomentumMC[2]*MomentumMC[2]);
296 normMomMC.push_back(MomentumMC[0]/length);
297 normMomMC.push_back(MomentumMC[1]/length);
298 normMomMC.push_back(MomentumMC[2]/length);
302 for(std::vector<ReconstructedParticle*>::const_iterator iParticle2 = MyJets.begin(); iParticle2 != MyJets.end() ;++iParticle2)
304 vector<double> normMom;
305 const double* Momentum = (*iParticle2)->getMomentum();
306 double length = sqrt(Momentum[0]*Momentum[0]+Momentum[1]*Momentum[1]+Momentum[2]*Momentum[2]);
307 normMom.push_back(Momentum[0]/length);
308 normMom.push_back(Momentum[1]/length);
309 normMom.push_back(Momentum[2]/length);
311 double disttemp= sqrt(((normMomMC[0]-normMom[0])*(normMomMC[0]-normMom[0]))+
312 ((normMomMC[1]-normMom[1])*(normMomMC[1]-normMom[1]))+
313 ((normMomMC[2]-normMom[2])*(normMomMC[2]-normMom[2])));
319 if (disttemp < fabs(dist))
321 angle = 180*asin( disttemp/2 )*2/3.14159265;
323 if(angle<_MaximumAngle)
326 JetValue = JetCounter;
332 std::cout<<
"Heavy MC Particle not assigned due to angle cut "<<std::endl;
339 if(dist>0 && JetValue >=0)
341 if(jetflavourdata[JetValue]>3)
343 if( jetflavourdata[JetValue] == MCflavour[MCcounter] )
345 if(fabs(jetangledata[JetValue]) > fabs(anglestored))
348 hadronchargedata[JetValue] = MCCharge[MCcounter];
349 partonchargedata[JetValue] = MCPCharge[MCcounter];
350 jetcodedata[JetValue] = int(MCCode[MCcounter]);
351 jetangledata[JetValue] = anglestored;
355 if( jetflavourdata[JetValue] < MCflavour[MCcounter])
357 jetflavourdata[JetValue] = MCflavour[MCcounter];
358 hadronchargedata[JetValue] = MCCharge[MCcounter];
359 partonchargedata[JetValue] = MCPCharge[MCcounter];
360 jetcodedata[JetValue] = int(MCCode[MCcounter]);
361 jetangledata[JetValue] = anglestored;
366 jetflavourdata[JetValue] = MCflavour[MCcounter];
367 hadronchargedata[JetValue] = MCCharge[MCcounter];
368 partonchargedata[JetValue] = MCPCharge[MCcounter];
369 jetcodedata[JetValue] = int(MCCode[MCcounter]);
370 jetangledata[JetValue] = anglestored;
379 for(
unsigned int iii=0; iii<jetflavourdata.size(); iii++)
381 LCFloatVec *OutVec =
new LCFloatVec();
382 OutVec->push_back(jetflavourdata[iii]);
383 OutVec->push_back(partonchargedata[iii]);
384 OutVec->push_back(hadronchargedata[iii]);
385 OutVec->push_back(jetcodedata[iii]);
387 OutCollection->addElement(OutVec);
388 std::cout <<
"TJ: Jet # "<< iii<<
" PDG: ";
389 std::cout <<jetcodedata[iii]<<
" Flavour: ";
390 std::cout <<jetflavourdata[iii]<<
" Parton Charge: ";
391 std::cout <<partonchargedata[iii] <<
" Hadron Charge: ";
392 std::cout <<hadronchargedata[iii] <<std::endl;
401 void TrueAngularJetFlavourProcessor::check( LCEvent * evt ) {
406 void TrueAngularJetFlavourProcessor::end(){
408 std::cout <<
"TrueAngularJetFlavourProcessor::end() " << name()
409 <<
" processed " << _nEvt <<
" events in " << _nRun <<
" runs "
417 float TrueAngularJetFlavourProcessor::chargefromPDG(
int code)
422 double charge = -1000;
424 if( abs(code) > 1000000000 )
426 std::cout<<
" Note: TrueAngularFetFlavourProcessor.cc: Hevy nucleus found and ignored. "<<std::endl;
430 if( abs(code) > 10000 )
436 first = (code/1000) %2;
437 second = (code/100) %2;
438 third = (code/10) %2;
452 charge = (first+second+third)/3;
460 else if (abs(code)>100)
462 first = (code/100) % 2;
463 second = (code/10) %2;
468 charge = (first-second)/3;
481 std::cerr<<
"Warning: TrueAngularFetFlavourProcessor.cc: A non hadronic particle has been choosen for vertex calculations "<<std::endl;
Determine MC Jet Flavour by angular matching of heavy quarks to jets, also determine hadronic and par...