LCFIVertex  0.7.2
TrueAngularJetFlavourProcessor.cc
1 #include "TrueAngularJetFlavourProcessor.h"
2 #include <iostream>
3 #include <fstream>
4 #include <sstream>
5 #include <cmath>
6 #include <cstdlib>
7 #include <algorithm>
8 
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"
17 
18 #include <UTIL/LCRelationNavigator.h>
19 #include <IMPL/ParticleIDImpl.h>
20 #include <IMPL/LCCollectionVec.h>
21 
22 #include <vector>
23 #include <string>
24 #include <map>
25 
26 #ifdef __APPLE__
27 #define uint uint32_t
28 #endif
29 
30 using namespace marlin ;
31 using namespace lcio;
32 using std::vector;
33 using std::string;
34 using std::map;
35 using EVENT::Track;
36 
37 TrueAngularJetFlavourProcessor aTrueAngularJetFlavourProcessor ;
38 
39 TrueAngularJetFlavourProcessor::TrueAngularJetFlavourProcessor() : Processor("TrueAngularJetFlavourProcessor") {
40 
41  // modify processor description
42  _description = "TrueAngularJetFlavourProcessor - Determines the true flavour of a jet from the MC Paticles associated to the Jets RP" ;
43 
44 
45  // register steering parameters: name, description, class-variable, default value
46 
47  registerInputCollection( lcio::LCIO::RECONSTRUCTEDPARTICLE,
48  "JetRPCollection" ,
49  "Name of the ReconstructedParticle collection that represents jets" ,
50  _JetRPColName ,
51  std::string("Jets") ) ;
52  registerInputCollection( lcio::LCIO::MCPARTICLE,
53  "MCParticleCollection" ,
54  "Name of the collection that holds all MC particles. " ,
55  _MCParticleColName ,
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" ,
64  _MaximumAngle,
65  double(180));
66  }
67 
68 void TrueAngularJetFlavourProcessor::init()
69 {
70  // usually a good idea to
71  printParameters() ;
72 
73  _nRun = 0 ;
74  _nEvt = 0 ;
75 }
76 
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");
83 
84  run->parameters().setValues(_TrueJetFlavourColName, _TrueJetVariableNames);
85 
86  _nRun++ ;
87 }
88 
89 void TrueAngularJetFlavourProcessor::processEvent( LCEvent * evt ) {
90 
91  LCCollection* JetRPCol = evt->getCollection( _JetRPColName );
92  LCCollection* MCParticleCol = evt->getCollection( _MCParticleColName );
93 
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;
99 
100 
101  LCCollectionVec* OutCollection = new LCCollectionVec("LCFloatVec");
102 
103  evt->addCollection(OutCollection,_TrueJetFlavourColName);
104 
105  int nRCP = JetRPCol->getNumberOfElements() ;
106 
107  if(nRCP ==0 ) std::cerr<<"Warning: TrueAngularFetFlavourProcessor.cc:88 : NO jets present "<<std::endl;
108 
109  for(int i=0; i< nRCP ; i++)
110  {
111  MyJets.push_back(dynamic_cast<ReconstructedParticle*>(JetRPCol->getElementAt(i)));
112  }
113 
114  vector<int> jetflavourdata(MyJets.size(),1);
115  vector<int> jetcodedata(MyJets.size(),0);
116  vector<double> jetangledata(MyJets.size(),1000);
117 
118  // Set default charge for light jets to -100
119  vector<float> hadronchargedata(MyJets.size(),-100);
120  vector<float> partonchargedata(MyJets.size(),-100);
121 
122  std::vector<std::vector<int> > Ntrkjet;
123  bool originalhadron = 0;
124  int code = 0;
125 
126  float charge =-100;
127  float pcharge =-100;
128 
129  std::vector<MCParticle*> HeavyMCs;
130  std::vector<MCParticle*> JetMCs;
131  std::vector<double> TrueFlavour;
132 
133  std::vector<int> typepart;
134  std::vector<MCParticle*> daughtercharge;
135 
136 
137  if( MCParticleCol->getNumberOfElements() == 0 ) std::cerr<<"Warning: TrueAngularFetFlavourProcessor.cc:107 : NO MC data presentjets present "<<std::endl;
138 
139 
140  for(int nMCpart = 0; nMCpart< MCParticleCol->getNumberOfElements();nMCpart++ )
141  {
142 
143  pcharge = -100;
144  MCParticle* MCused = dynamic_cast<MCParticle*> (MCParticleCol->getElementAt(nMCpart));
145 
146  //note this whole process relies on 551/100 = 5 in integer cast!!! (which is true here)
147 
148  code = abs( MCused->getPDG() );
149  charge = MCused->getCharge();
150 
151  //NOTE: this works only with mokka-06-04-p02 or above.
152  // In case you are looking at older versions of mokka please remove the if statement.
153  if(charge == -1000 && code >100)
154  {
155  charge = chargefromPDG(MCused->getPDG());
156  }
157 
158  //codes with six digits or more are for special purposes only
159  if( code > 100000 )
160  {
161  code = 0;
162  }
163 
164  //this little addition will take care of wierder mesons
165  if( code > 10000 )
166  {
167  code = code%1000;
168  }
169 
170  if( code > 1000 )
171  {
172 
173  if( code/ 1000 == 4 )
174  {
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 )
180  { pcharge = 1; }
181  else
182  { pcharge = -1;}
183  MCPCharge.push_back(pcharge);
184  }
185 
186  if( code/ 1000 == 5 )
187  {
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 )
193  {pcharge = 1; }
194  else
195  {pcharge = -1;}
196  MCPCharge.push_back(pcharge);
197  }
198  }
199 
200  if( code > 100 )
201  {
202  if( code / 100 == 4 )
203  {
204  HeavyMCs.push_back(MCused);
205  MCflavour.push_back(4);
206  MCCharge.push_back(charge);
207  MCCode.push_back(MCused->getPDG());
208 
209  if ( (code/10) %11 == 0)
210  {
211  pcharge = 0;
212  }
213  else
214  {
215  if( MCused->getPDG() > 0 )
216  { pcharge = 1; }
217  else
218  { pcharge = -1;}
219  }
220  MCPCharge.push_back(pcharge);
221  }
222 
223  if( code/100 == 5 )
224  {
225  HeavyMCs.push_back(MCused);
226  MCflavour.push_back(5);
227  MCCharge.push_back(charge);
228  MCCode.push_back(MCused->getPDG());
229 
230  if ( code %11 == 0)
231  {
232  pcharge = 0;
233  }
234  else
235  {
236  if( MCused->getPDG() < 0 )
237  { pcharge = 1; }
238  else
239  { pcharge = -1;}
240  }
241  MCPCharge.push_back(pcharge);
242  }
243 
244  }
245  //std::cout<< "hadron, "<<MCused->getPDG() << ", "<< pcharge<<std::endl;
246  }
247 
248 
249  for (std::vector<MCParticle*>::const_iterator iParticle = HeavyMCs.begin(); iParticle != HeavyMCs.end() ;++iParticle)
250  {
251  originalhadron =0;
252  MCParticle* MCused = *iParticle;
253  while(!originalhadron)
254  {
255  // at hadron level particles can have only 1 parent only at parton level
256  // there might be more than 1 parent. if we are at parton level we are too far anyway.
257  // refer to LCIO manual for details.
258  if( MCused->getParents().size() == 1)
259  {
260  // std::cout<<" ParentalPDG "<<MCused->getParents()[0]->getPDG()<<std::endl;
261  if (abs( MCused->getParents()[0]->getPDG()) >100 || (abs(MCused->getParents()[0]->getPDG())>10 && abs(MCused->getParents()[0]->getPDG())<81))
262  {
263  MCused = MCused->getParents()[0];
264  }
265  else
266  {
267  JetMCs.push_back(MCused);
268  originalhadron = 1;
269  }
270  }
271  else
272  {
273  JetMCs.push_back(MCused);
274  originalhadron = 1;
275  }
276  }
277  }
278 
279  // std::cout<<JetMCs.size()<<std::endl;
280 
281  if(MyJets.size()>0)
282  {
283 
284 
285  int MCcounter =0;
286  for(std::vector<MCParticle*>::const_iterator iParticle3 = JetMCs.begin(); iParticle3 != JetMCs.end() ;++iParticle3)
287  {
288  double dist = -999999999;
289  int JetCounter =0;
290  int JetValue =-1;
291  double angle = 0;
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);
299 
300  // std::cout<<"numberUSED "<< (*iParticle3)->getPDG() <<std::endl;
301 
302  for(std::vector<ReconstructedParticle*>::const_iterator iParticle2 = MyJets.begin(); iParticle2 != MyJets.end() ;++iParticle2)
303  {
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);
310 
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])));
314 
315  // std::cout<<"Disttemp "<<JetCounter<< " "<<disttemp<<std::endl;
316 
317 
318 
319  if (disttemp < fabs(dist))
320  {
321  angle = 180*asin( disttemp/2 )*2/3.14159265;
322 
323  if(angle<_MaximumAngle)
324  {
325  dist = disttemp;
326  JetValue = JetCounter;
327  anglestored = angle;
328  // std::cout<<"Inside "<<(*iParticle3)->getPDG()<<" "<<JetValue<<" "<<anglestored<<std::endl;
329  }
330  else
331  {
332  std::cout<< "Heavy MC Particle not assigned due to angle cut "<<std::endl;
333  }
334  }
335 
336  JetCounter++;
337  }
338 
339  if(dist>0 && JetValue >=0)
340  {
341  if(jetflavourdata[JetValue]>3)
342  {
343  if( jetflavourdata[JetValue] == MCflavour[MCcounter] )
344  {
345  if(fabs(jetangledata[JetValue]) > fabs(anglestored))
346  {
347 
348  hadronchargedata[JetValue] = MCCharge[MCcounter];
349  partonchargedata[JetValue] = MCPCharge[MCcounter];
350  jetcodedata[JetValue] = int(MCCode[MCcounter]);
351  jetangledata[JetValue] = anglestored;
352 
353  }
354  }
355  if( jetflavourdata[JetValue] < MCflavour[MCcounter])
356  {
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;
362  }
363  }
364  else
365  {
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;
371  }
372  }
373 
374  MCcounter++;
375 
376  }
377  }
378 
379  for(unsigned int iii=0; iii<jetflavourdata.size(); iii++)
380  {
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]);
386 
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;
393  };
394  //Create the collection to store the results
395 
396  _nEvt ++ ;
397 }
398 
399 
400 
401 void TrueAngularJetFlavourProcessor::check( LCEvent * evt ) {
402  // nothing to check here - could be used to fill checkplots in reconstruction processor
403 }
404 
405 
406 void TrueAngularJetFlavourProcessor::end(){
407 
408  std::cout << "TrueAngularJetFlavourProcessor::end() " << name()
409  << " processed " << _nEvt << " events in " << _nRun << " runs "
410  << std::endl ;
411 
412 }
413 
414 //reconstruct quark content from pdg code
415 //I am surprised I could not find this anywhere already coded
416 //This is needed because sometimes mokka gives us the default -1000 and not the correct charge
417 float TrueAngularJetFlavourProcessor::chargefromPDG(int code)
418 {
419  double first = 0;
420  double second = 0;
421  double third = 0;
422  double charge = -1000;
423 
424  if( abs(code) > 1000000000 )
425  {
426  std::cout<<" Note: TrueAngularFetFlavourProcessor.cc: Hevy nucleus found and ignored. "<<std::endl;
427  return -1000;
428  }
429 
430  if( abs(code) > 10000 )
431  {
432  code = code%1000;
433  }
434  if(fabs(code) >1000)
435  {
436  first = (code/1000) %2;
437  second = (code/100) %2;
438  third = (code/10) %2;
439 
440  if( third == 0)
441  {third = third+2;}
442  else
443  {third = -1; }
444  if (second == 0)
445  {second = second+2;}
446  else
447  {second = -1; }
448  if (first ==0)
449  {first = first+2;}
450  else
451  {first = -1; }
452  charge = (first+second+third)/3;
453  if(charge ==0)
454  {return 0;}
455  else if (code>0)
456  {return charge;}
457  else
458  {return (-charge);}
459  }
460  else if (abs(code)>100)
461  {
462  first = (code/100) % 2;
463  second = (code/10) %2;
464  if (second ==0)
465  {second = second+2;}
466  if (first ==0)
467  {first = first+2;}
468  charge = (first-second)/3;
469  if (charge != 0)
470  {
471  if(code>0)
472  {return 1;}
473  if(code<0)
474  {return -1;}
475  }
476  else
477  {return 0;}
478  }
479  else
480  {
481  std::cerr<<"Warning: TrueAngularFetFlavourProcessor.cc: A non hadronic particle has been choosen for vertex calculations "<<std::endl;
482  return -1000;
483  }
484  // if we are here something probably went wrong.
485  return -1000;
486 }
TrueAngularJetFlavourProcessor
Determine MC Jet Flavour by angular matching of heavy quarks to jets, also determine hadronic and par...
Definition: TrueAngularJetFlavourProcessor.h:53
Generated on Fri Dec 2 2016 12:36:49 for LCFIVertex by   doxygen 1.8.6