All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Groups
TruthVertexFinder.cc
Go to the documentation of this file.
1 #include "TruthVertexFinder.hh"
2 using std::string;
3 using std::vector;
4 using std::map;
5 namespace TTbarAnalysis
6 {
7  TruthVertexFinder aTruthVertexFinder ;
8  TruthVertexFinder::TruthVertexFinder() : Processor("TruthVertexFinder")
9  {
10 
11  _description = "TruthVertexFinder extructs a secondary vertex from generator collection." ;
12 
13 
14  registerInputCollection( LCIO::MCPARTICLE,
15  "CollectionName" ,
16  "Name of the MCParticle collection" ,
17  _colName ,
18  std::string("MCParticleSkimmed")
19  );
20  registerOutputCollection( LCIO::VERTEX,
21  "OutputCollectionName" ,
22  "Name of the Vertex collection" ,
23  _outputcolName ,
24  std::string("MCVertex")
25  );
26  registerOutputCollection( LCIO::MCPARTICLE,
27  "OutputBStarName" ,
28  "Name of the Vertex collection" ,
29  _outputBStarName ,
30  std::string("BStar")
31  );
32  registerOutputCollection( LCIO::MCPARTICLE,
33  "OutputProngsName" ,
34  "Name of the Prongs collection" ,
35  _outputProngsName ,
36  std::string("EGProngs")
37  );
38  /*registerOutputCollection( LCIO::MCPARTICLE,
39  "OutputBStarName" ,
40  "Name of the Vertex collection" ,
41  _outputBStarName ,
42  std::string("BStar")
43  );*/
44  registerOutputCollection(LCIO::MCPARTICLE,
45  "QuarkCollectionName" ,
46  "Name of the b-quark collection" ,
47  _outputquarkcolName,
48  std::string("MCbquarks")
49  );
50  registerInputCollection( LCIO::RECONSTRUCTEDPARTICLE,
51  "RelCollectionName",
52  "Name of the Jet relation collection",
53  _colRelName,
54  std::string("RecoMCTruthLink")
55  );
56  registerProcessorParameter("DecayChainPDGs" ,
57  "PDGs of desired decays" ,
58  inputPdg,
59  vector<int>()
60  );
61  _tagParameter = 6;
62  registerProcessorParameter("tagPDG" ,
63  "PDG of desired particle" ,
64  _tagParameter,
65  _tagParameter
66  );
67  _writeROOTparameter = 0;
68  registerProcessorParameter("writeROOT" ,
69  "write ROOT file?" ,
70  _writeROOTparameter,
71  _writeROOTparameter
72  );
73  _initialQuarkPDGparameter = 5;
74  registerProcessorParameter("initialQuarkPDG" ,
75  "PDG of initial particle" ,
76  _initialQuarkPDGparameter,
77  _initialQuarkPDGparameter
78  );
79  _aParameter = 0.005;
80  registerProcessorParameter("a" ,
81  "a parameter of accuracy in mm" ,
82  _aParameter,
83  _aParameter
84  );
85  _bParameter = 0.01;
86  registerProcessorParameter("b" ,
87  "b parameter of accuracy in mm" ,
88  _bParameter,
89  _bParameter
90  );
91  _writeBonlyParameter = 1;
92  registerProcessorParameter("writeBonly" ,
93  "b parameter" ,
94  _writeBonlyParameter,
95  _writeBonlyParameter
96  );
97  _hfilename = "GenVertexTest.root";
98  registerProcessorParameter("ROOTFileName" ,
99  "ROOT File Name" ,
100  _hfilename,
101  _hfilename
102  );
103  //_pdgs.push_back(BOTTOM_HADRONS);
104  //_pdgs.push_back(CHARMED_HADRONS);
105  /*_pdgs.push_back(BOTTOM_MESONS);
106  _pdgs.push_back(CHARMED_MESONS);
107  _pdgs.push_back(EXCEPTIONAL_PDGTYPE);*/
108  ip[0] = 0.0;
109  ip[1] = 0.0;
110  ip[2] = 0.0;
111  ClearVariables();
112  }
113 
114 
115 
116  void TruthVertexFinder::init()
117  {
118  streamlog_out(MESSAGE) << " init called " << std::endl;
119  printParameters() ;
120  if (inputPdg.size() < 2)
121  {
122  streamlog_out(ERROR) << " ERROR: size " << inputPdg.size() << " - Running on default settings" << std::endl;
123  inputPdg.push_back(500);
124  inputPdg.push_back(400);
125  inputPdg.push_back(0);
126  }
127  for (unsigned int i = 0; i < inputPdg.size(); i++)
128  {
129  _pdgs.push_back((PDGTYPE)inputPdg[i]);
130  }
131  _nRun = 0 ;
132  _nEvt = 0 ;
133  if (_writeROOTparameter == 0)
134  {
135  return;
136  }
137  _hfile = new TFile( _hfilename.c_str(), "RECREATE", _hfilename.c_str() ) ;
138  _hBStarTree = new TTree( "BStar", "My test tree!" );
139  _hBStarTree->Branch("bstarnumber", &_bstarnumber, "bstarnumber/I");
140  _hBStarTree->Branch("bstarmomentum", _bstarmomentum, "bstarmomentum[bstarnumber]/F");
141  _hBStarTree->Branch("bstaroffset", _bstaroffset, "bstaroffset[bstarnumber]/F");
142  _hVertexTree = new TTree( "Vertices", "My test tree!" );
143  _hVertexTree->Branch("numberOfVertices", &_numberOfVertexes, "numberOfVertexes/I");
144  _hVertexTree->Branch("distance", _distanceFromIP, "distance[numberOfVertexes]/F");
145  _hVertexTree->Branch("coordinates", _coordinates, "coordinates[numberOfVertexes][3]/F");
146  _hVertexTree->Branch("PDG", _PDG, "PDG[numberOfVertexes]/I");
147  _hVertexTree->Branch("charge", _charge, "charge[numberOfVertexes]/I");
148  _hVertexTree->Branch("generation", _generation, "generation[numberOfVertexes]/I");
149  _hVertexTree->Branch("numberOfParticles", _numberOfParticles, "numberOfParticles[numberOfVertexes]/I");
150  //_hVertexTree->Branch("energyOfParticles", _energyOfParticles, "energyOfParticles[numberOfVertexes][15]/F");
151  _hVertexTree->Branch("momentumOfParticles", _momentumOfParticles, "momentumOfParticles[numberOfVertexes][15]/F");
152  _hVertexTree->Branch("massOfParticles", _massOfParticles, "massOfParticles[numberOfVertexes][15]/F");
153  _hVertexTree->Branch("interactionOfParticles", _interactionOfParticles, "interactionOfParticles[numberOfVertexes][15]/I");
154  _hTree = new TTree( "Stats", "My test tree!" );
155  _hTree->Branch("tag", &_tag, "tag/I");
156  _hTree->Branch("btotalnumber", &_btotalnumber, "btotalnumber/I");
157  _hTree->Branch("bbartotalnumber", &_bbartotalnumber, "bbartotalnumber/I");
158  if (_writeROOTparameter < 2)
159  {
160  return;
161  }
162  _hTree->Branch("cosquark", &_cosquark, "cosquark/F");
163  _hTree->Branch("cosantiquark", &_cosantiquark, "cosantiquark/F");
164  _hTree->Branch("totalBcharge", &_totalBcharge, "totalBcharge/I");
165  _hTree->Branch("ccharge", &_ccharge, "ccharge/I");
166  _hTree->Branch("cbarcharge", &_cbarcharge, "cbarcharge/I");
167  _hTree->Branch("bcharge", &_bcharge, "bcharge/I");
168  _hTree->Branch("bbarcharge", &_bbarcharge, "bbarcharge/I");
169  _hTree->Branch("baccuracy", &_baccuracy, "baccuracy/F");
170  _hTree->Branch("bbaraccuracy", &_bbaraccuracy, "bbaraccuracy/F");
171  _hTree->Branch("bIPdistance", &_bIPdistance, "bIPdistance/F");
172  _hTree->Branch("bbarIPdistance", &_bbarIPdistance, "bbarIPdistance/F");
173  _hTree->Branch("bdistance", &_bdistance, "bdistance/F");
174  _hTree->Branch("bbardistance", &_bbardistance, "bbardistance/F");
175  _hTree->Branch("bmomentum", &_bmomentum, "bmomentum/F");
176  _hTree->Branch("bbarmomentum", &_bbarmomentum, "bbarmomentum/F");
177  _hTree->Branch("cmomentum", &_cmomentum, "cmomentum/F");
178  _hTree->Branch("cbarmomentum", &_cbarmomentum, "cbarmomentum/F");
179  _hTree->Branch("caccuracy", &_caccuracy, "caccuracy/F");
180  _hTree->Branch("cbaraccuracy", &_cbaraccuracy, "cbaraccuracy/F");
181  _hTree->Branch("bnumber", &_bnumber, "bnumber/I");
182  _hTree->Branch("bbarnumber", &_bbarnumber, "bbarnumber/I");
183  _hTree->Branch("cnumber", &_cnumber, "cnumber/I");
184  _hTree->Branch("cbarnumber", &_cbarnumber, "cbarnumber/I");
185  _hTree->Branch("bptmiss", &_bptmiss, "bptmiss/F");
186  _hTree->Branch("bbarptmiss", &_bbarptmiss, "bbarptmiss/F");
187  _hTree->Branch("bnumber_f", &_bnumber_f, "bnumber_f/I");
188  _hTree->Branch("bbarnumber_f", &_bbarnumber_f, "bbarnumber_f/I");
189  _hTree->Branch("cnumber_f", &_cnumber_f, "cnumber_f/I");
190  _hTree->Branch("cbarnumber_f", &_cbarnumber_f, "cbarnumber_f/I");
191  //_hTree->Branch("firstVertexDistance", _firstVertexDistance, "firstVertexDistance[numberOfB0]/F");
192  //_hTree->Branch("secondVertexDistance", _secondVertexDistance, "secondVertexDistance[numberOfB0]/F");
193  //******************************************************************************************************//
194  _hTrackTree = new TTree( "Tracks", "My test tree!" );
195  _hTrackTree->Branch("bnumber", &_bnumber, "bnumber/I");
196  _hTrackTree->Branch("bbarnumber", &_bbarnumber, "bbarnumber/I");
197  _hTrackTree->Branch("cnumber", &_cnumber, "cnumber/I");
198  _hTrackTree->Branch("cbarnumber", &_cbarnumber, "cbarnumber/I");
199 
200  _hTrackTree->Branch("bptrack", _bptrack, "bptrack[bnumber]/F");
201  _hTrackTree->Branch("cptrack", _cptrack, "cptrack[cnumber]/F");
202  _hTrackTree->Branch("bbarptrack", _bbarptrack, "bbarptrack[bbarnumber]/F");
203  _hTrackTree->Branch("cbarptrack", _cbarptrack, "cbarptrack[cbarnumber]/F");
204 
205  _hTrackTree->Branch("betatrack", _betatrack, "betatrack[bnumber]/F");
206  _hTrackTree->Branch("cetatrack", _cetatrack, "cetatrack[cnumber]/F");
207  _hTrackTree->Branch("bbaretatrack", _bbaretatrack, "bbaretatrack[bbarnumber]/F");
208  _hTrackTree->Branch("cbaretatrack", _cbaretatrack, "cbaretatrack[cbarnumber]/F");
209 
210  _hTrackTree->Branch("boffsettrack", _boffsettrack, "boffsettrack[bnumber]/F");
211  _hTrackTree->Branch("coffsettrack", _coffsettrack, "coffsettrack[cnumber]/F");
212  _hTrackTree->Branch("bbaroffsettrack", _bbaroffsettrack, "bbaroffsettrack[bbarnumber]/F");
213  _hTrackTree->Branch("cbaroffsettrack", _cbaroffsettrack, "cbaroffsetttrack[cbarnumber]/F");
214 
215  _hMisRecoTree = new TTree( "Misreco", "My test tree!" );
216  _hMisRecoTree->Branch("misreconumber", &_misreconumber, "misreconumber/I");
217  _hMisRecoTree->Branch("misrecotheta", _misrecotheta, "misrecotheta[misreconumber]/F");
218  _hMisRecoTree->Branch("misrecocostheta", _misrecocostheta, "misrecocostheta[misreconumber]/F");
219  _hMisRecoTree->Branch("misrecomomentum", _misrecomomentum, "misrecomomentum[misreconumber]/F");
220  _hMisRecoTree->Branch("misrecopt", _misrecopt, "misrecopt[misreconumber]/F");
221 
222  }
223 
224 
225  void TruthVertexFinder::processRunHeader( LCRunHeader* /*run*/)
226  {
227  _nRun++ ;
228  }
229 
230  void TruthVertexFinder::PrintParticle(MCParticle * particle)
231  {
232  if (!particle)
233  {
234  return;
235  }
236  streamlog_out(MESSAGE) << std::fixed << std::setw( 6 ) << std::setprecision( 3 ) << std::setfill( ' ' );
237  streamlog_out(MESSAGE) <<"|"<<particle->getPDG() <<"\t\t|"<<particle->getMass()<<"\t\t|"<<particle->getCharge() <<"\t\t|"<<particle->getEnergy()<<"\t\t|"<<particle->getVertex()[0]<<"\t\t|"<<particle->getVertex()[1]<<"\t\t|"<<particle->getVertex()[2] <<"\t\t|\n";
238 
239  }
240  void TruthVertexFinder::Write(vector< Vertex * > * vertices, int & number)
241  {
242  if (!vertices || vertices->size() < _pdgs.size()-1)
243  {
244  return;
245  }
246  for (unsigned int i = 0; i < vertices->size(); i++)
247  {
248  Vertex * vertex = vertices->at(i);
249  _PDG[_numberOfVertexes] = vertex->getParameters()[1];
250  _generation[_numberOfVertexes] = vertex->getParameters()[2];
251  ReconstructedParticle * particle = vertex->getAssociatedParticle();
252  _coordinates[_numberOfVertexes][0] = vertex->getPosition()[0];
253  _coordinates[_numberOfVertexes][1] = vertex->getPosition()[1];
254  _coordinates[_numberOfVertexes][2] = vertex->getPosition()[2];
255  _charge[_numberOfVertexes] = particle->getCharge();
256  _numberOfParticles[_numberOfVertexes] = particle->getParticles().size();
257  _distanceFromIP[_numberOfVertexes] = vertex->getParameters()[0];
258  MyVertex * myvertex = static_cast< MyVertex * >(vertex);
259  for (unsigned int j = 0; j < particle->getParticles().size(); j++)
260  {
261  _momentumOfParticles[_numberOfVertexes][j] = MathOperator::getModule(particle->getParticles()[j]->getMomentum());
262  _massOfParticles[_numberOfVertexes][j] = particle->getParticles()[j]->getMass();
263  _interactionOfParticles[_numberOfVertexes][j] = myvertex->__GetMCParticles()[j]->isDecayedInCalorimeter();
264  }
265  _numberOfVertexes++;
266  }
267  }
268  void TruthVertexFinder::PrintChain(vector< MCParticle * > * chain)
269  {
270  if (!chain)
271  {
272  return;
273  }
274  for (unsigned int i = 0; i < chain->size(); i++)
275  {
276  PrintParticle(chain->at(i));
277  }
278  }
279  void TruthVertexFinder::WriteQuarksCollection(LCEvent * evt, std::vector< MCParticle * > & quarks)
280  {
281  IMPL::LCCollectionVec * mc = new IMPL::LCCollectionVec ( LCIO::MCPARTICLE ) ;
282  if (quarks.size() == 2)
283  {
284  //streamlog_out(MESSAGE) << "PDG: " << quarks[0]->getPDG() << '\n';
285  //streamlog_out(MESSAGE) << "PDG: " << quarks[1]->getPDG() << '\n';
286  mc->addElement(quarks[0]);
287  mc->addElement(quarks[1]);
288  }
289  evt->addCollection( mc , _outputquarkcolName ) ;
290  }
291 
292  void TruthVertexFinder::AddProngs( VertexMCOperator & vertexOperator, MCOperator & opera, DecayChain * chain, vector< Vertex * > * verticies, std::vector<int> & parameters, IMPL::LCCollectionVec * col)
293  {
294  if (!verticies || !chain)
295  {
296  return;
297  }
298  int vsize = verticies->size();
299  bool useRelation = false;
300  if (vsize > 1)
301  {
302  vector< MCParticle * > bdaughters = opera.SelectStableCloseDaughters(chain->Get(0), chain->Get(1)->getPDG());
303  vertexOperator.AddProngs(verticies->at(0), bdaughters, useRelation);
304  vector< MCParticle * > cdaughters = opera.SelectStableCloseDaughters(chain->Get(1));
305  vertexOperator.AddProngs(verticies->at(1), cdaughters, useRelation);
306  if (col)
307  {
308  for (unsigned int i = 0; i < bdaughters.size(); i++)
309  {
310  parameters.push_back( (float) chain->GetParentPDG() / 5.0 * 2);
311  col->addElement(bdaughters[i]);
312  }
313  for (unsigned int i = 0; i < cdaughters.size(); i++)
314  {
315  parameters.push_back( (float) chain->GetParentPDG() / 5.0 * 3);
316  col->addElement(cdaughters[i]);
317  }
318  }
319  }
320  else
321  {
322  _bbartotalnumber = (_bbartotalnumber > 0)? _bbartotalnumber : 0;
323  _btotalnumber = (_btotalnumber > 0)? _btotalnumber: 0;
324  PrintParticle(chain->Get(0));
325  vector< MCParticle * > bdaughters = opera.SelectStableCloseDaughters(chain->Get(0));
326  streamlog_out(MESSAGE) << "Prongs for quark " << chain->GetParentPDG() << ": \n";
327  for (unsigned int i = 0; i < bdaughters.size(); i++)
328  {
329  PrintParticle(bdaughters[i]);
330  }
331  vertexOperator.AddProngs(verticies->at(0), bdaughters, useRelation);
332  if (col)
333  {
334  for (unsigned int i = 0; i < bdaughters.size(); i++)
335  {
336  col->addElement(bdaughters[i]);
337  if (chain->GetParentPDG() > 0)
338  {
339  _btotalnumber++;
340  }
341  else
342  {
343  _bbartotalnumber++;
344  }
345  }
346  }
347  }
348  }
349 
350  void TruthVertexFinder::processEvent( LCEvent * evt )
351  {
352  try
353  {
354  LCCollection* col = evt->getCollection( _colName );
355  streamlog_out(MESSAGE) << "***********TruthVertexFinder*"<<_nEvt<<"***************\n";
356  LCCollection* rel = evt->getCollection(_colRelName);
357  MCOperator opera(col,rel);
358  VertexMCOperator vertexOperator(rel);
359 
360  _tag = (opera.CheckProcessForPair(_tagParameter))? 1 : 0;
361  vector< MCParticle * > bquarks = opera.GetPairParticles(_pdgs[0]);
362  //GetAsymmetry(bquarks);
363  _nEvt ++ ;
364  int initQuark = abs(_initialQuarkPDGparameter);
365  streamlog_out(MESSAGE) <<"\t|PDG\t\t|Mass\t\t|Charge\t\t|Energy\t\t|Vtx X\t\t|Vtx Y\t\t|Vtx Z\t\t|\n";
366  DecayChain * bChainRaw = opera.Construct(string("b-quark decay chain"), initQuark, _pdgs);
367  DecayChain * bChain = opera.RefineDecayChain(bChainRaw, _pdgs);
368  IMPL::LCCollectionVec * mc = new IMPL::LCCollectionVec ( LCIO::MCPARTICLE ) ;
369  mc->setSubset();
370  if (bChain)
371  {
372  vector< MCParticle * > daughters = opera.SelectStableCloseDaughters(bChainRaw->Get(0), bChain->Get(0)->getPDG());
373  streamlog_out(MESSAGE) <<"Additional B particles: \n";
374  for (unsigned int i = 0; i < daughters.size(); i++)
375  {
376  vector< float > direction = MathOperator::getDirection(daughters[i]->getMomentum());
377  _bstaroffset[_bstarnumber] = MathOperator::getDistanceTo(ip, direction, bChain->Get(1)->getVertex());
378  _bstarmomentum[_bstarnumber++] = MathOperator::getModule(daughters[i]->getMomentum());
379  mc->addElement(daughters[i]);
380  PrintParticle(daughters[i]);
381  }
382  }
383  DecayChain * bbarChainRaw = opera.Construct(string("bbar-quark decay chain"), 0-initQuark, _pdgs);
384  DecayChain * bbarChain = opera.RefineDecayChain(bbarChainRaw, _pdgs);
385  if (bbarChain)
386  {
387  vector< MCParticle * > daughters = opera.SelectStableCloseDaughters(bbarChainRaw->Get(0), bbarChain->Get(0)->getPDG());
388  streamlog_out(MESSAGE) <<"Additional Bbar particles: \n";
389  for (unsigned int i = 0; i < daughters.size(); i++)
390  {
391  vector< float > direction = MathOperator::getDirection(daughters[i]->getMomentum());
392  _bstaroffset[_bstarnumber] = MathOperator::getDistanceTo(ip, direction, bbarChain->Get(1)->getVertex());
393  _bstarmomentum[_bstarnumber++] = MathOperator::getModule(daughters[i]->getMomentum());
394  mc->addElement(daughters[i]);
395  PrintParticle(daughters[i]);
396  }
397  }
398  evt->addCollection( mc , _outputBStarName ) ;
399 
400  vector< Vertex * > * bverticies = vertexOperator.Construct(bChain);
401  vector< Vertex * > * bbarverticies = vertexOperator.Construct(bbarChain);
402  IMPL::LCCollectionVec * prongs = new IMPL::LCCollectionVec ( LCIO::MCPARTICLE ) ;
403  vector <int> parameters;
404  prongs->setSubset ();
405  //prongs->setFlag (16);
406  AddProngs(vertexOperator, opera, bChain, bverticies,parameters, prongs);
407  AddProngs(vertexOperator, opera, bbarChain, bbarverticies,parameters, prongs);
408 
409 
410 
411  WriteQuarksCollection(evt, bquarks);
412  WriteVertexCollection(evt, bverticies, bbarverticies);
413 
414  evt->addCollection( prongs , _outputProngsName);
415  prongs->parameters().setValues("trackIDs", parameters);
416  prongs->parameters().setValue("GenVertexTag", _tag);
417  //int number = 0;
418  if (_writeROOTparameter > 0)
419  {
420  _numberOfVertexes = 0;
421  Write(bverticies,_numberOfVertexes);
422  Write(bbarverticies, _numberOfVertexes);
423  if (_writeROOTparameter > 1)
424  {
425  Write(opera, bChain,bverticies);
426  Write(opera, bbarChain,bbarverticies);
427  _bnumber = (_bnumber < 0)? 0: _bnumber;
428  _bbarnumber = (_bbarnumber < 0)? 0: _bbarnumber;
429  _cnumber = (_cnumber < 0)? 0: _cnumber;
430  _cbarnumber = (_cbarnumber < 0)? 0: _cbarnumber;
431  _hTrackTree->Fill();
432 
433  }
434  _hTree->Fill();
435  _hBStarTree->Fill();
436  _hVertexTree->Fill();
437  ClearVariables();
438  }
439  //streamlog_out(MESSAGE) <<"B cos: " << _cosquark << '\n';
440  //streamlog_out(MESSAGE) <<"Bbar cos: " <<_cosantiquark << '\n';
441  streamlog_out(MESSAGE) << "b number: " << _btotalnumber << "\n";
442  streamlog_out(MESSAGE) << "bbar number: " << _bbartotalnumber << "\n";
443 
444  }
445  catch( DataNotAvailableException &e)
446  {
447  streamlog_out(MESSAGE) << "No collection!" << std::endl ;
448  }
449  }
450  void TruthVertexFinder::GetAsymmetry(std::vector< MCParticle * > & particles)
451  {
452  /*if (!particles)
453  {
454  return;
455  }*/
456  if (particles.size() < 2)
457  {
458  return;
459  }
460  vector< float > direction = MathOperator::getDirection(particles.at(0)->getEndpoint());
461  vector< float > antidirection = MathOperator::getDirection(particles.at(1)->getEndpoint());
462  _cosquark = std::cos(MathOperator::getAngles(direction)[1]);
463  _cosantiquark = std::cos(MathOperator::getAngles(antidirection)[1]);
464  //delete particles;
465  }
466  void TruthVertexFinder::Write(MCOperator & opera, DecayChain * chain, vector< Vertex * > * verticies)
467  {
468  if (!chain || !chain->Get(0))
469  {
470  return;
471  }
472  for (int i = 0; i < chain->GetSize(); i++)
473  {
474  PrintParticle(chain->Get(i));
475  }
476  //vector< MCParticle * > * misreco = new vector< MCParticle * >();
477  vector< float > direction = MathOperator::getDirection(chain->Get(0)->getMomentum());
478  if (chain->GetParentPDG() > 0 && chain->GetSize() > 1)
479  {
480 
481  const vector< MCParticle * > daughters = static_cast< MyVertex * >(verticies->at(0))->__GetMCParticles(); // opera.SelectStableCloseDaughters(chain->Get(0), chain->Get(1)->getPDG()); //opera.ScanForVertexParticles(bverticies->at(0)->getPosition(), 1e-3);
482  _bcharge = (int)chain->Get(0)->getCharge();
483  _bnumber = daughters.size();
484  _bnumber_f = opera.SelectStableCloseDaughters(chain->Get(0), chain->Get(1)->getPDG(),true).size();//opera.SelectStableCloseDaughters(chain->Get(0), chain->Get(1)->getPDG(),true).size();//opera.CheckDaughterVisibility(daughters).size();
485  _bIPdistance = verticies->at(0)->getParameters()[0];
486  Write(daughters, 1);
487  streamlog_out(MESSAGE) <<"Vertex b-quark: " << verticies->at(0)->getParameters()[0]<< " n-tracks: " << _bnumber << '\n';
488 
489  const vector< MCParticle * > cdaughters =static_cast< MyVertex * >(verticies->at(1))->__GetMCParticles(); //opera.SelectStableCloseDaughters(chain->Get(1)); //opera.ScanForVertexParticles(bverticies->at(1)->getPosition(), 1e-3);
490  _cnumber = cdaughters.size();
491  _cnumber_f = opera.SelectStableCloseDaughters(chain->Get(1),0,true).size();//opera.CheckDaughterVisibility(cdaughters).size();
492 
493  opera.CheckDaughterVisibility(daughters);
494  Write(cdaughters, 2);
495 
496  streamlog_out(MESSAGE) <<"Vertex c-quark: " << verticies->at(1)->getParameters()[0] <<" n-tracks: " << _cnumber << '\n';
497  opera.CheckDaughterVisibility(cdaughters);
498  _bdistance = MathOperator::getDistance(verticies->at(1)->getPosition(), verticies->at(0)->getPosition());
499  _btotalnumber = _cnumber + _bnumber;
500  streamlog_out(MESSAGE) <<"Checking b-quark meson...\n";
501  bool compatible = opera.CheckCompatibility(daughters, chain->Get(0), chain->Get(1)->getCharge());
502 
503  streamlog_out(MESSAGE) <<"Checking c-quark meson...\n";
504  compatible = opera.CheckCompatibility(cdaughters, chain->Get(1));
505  _bptmiss = getMissingPt(daughters, cdaughters, verticies->at(0));
506  streamlog_out(MESSAGE) <<"Missing pt for b-quark hadron: " << _bptmiss << "\n";
507  _ccharge = (int)chain->Get(1)->getCharge();
508  _bmomentum = MathOperator::getModule(chain->Get(0)->getMomentum());
509  _baccuracy = opera.GetAccuracy(chain->Get(0), _aParameter, _bParameter);
510  _cmomentum = MathOperator::getModule(chain->Get(1)->getMomentum());
511  _caccuracy = opera.GetAccuracy(chain->Get(1), _aParameter, _bParameter);
512  _cosquark = std::cos(MathOperator::getAngles(direction)[1]);
513  }
514  if (chain->GetParentPDG() < 0 && chain->GetSize() > 1)
515  {
516  const vector< MCParticle * > daughters = static_cast< MyVertex * >(verticies->at(0))->__GetMCParticles();// opera.SelectStableCloseDaughters(chain->Get(0), chain->Get(1)->getPDG()); // opera.ScanForVertexParticles(bbarverticies->at(0)->getPosition(), 1e-3);
517  Write(daughters, -1);
518  _bbarcharge = (int) chain->Get(0)->getCharge();
519  _bbarnumber = daughters.size();
520  _bbarIPdistance = verticies->at(0)->getParameters()[0];
521  _bbarnumber_f = opera.SelectStableCloseDaughters(chain->Get(0), chain->Get(1)->getPDG(),true).size();//opera.CheckDaughterVisibility(daughters).size();
522  streamlog_out(MESSAGE) <<"Vertex bbar-quark"<< 0 <<": " << verticies->at(0)->getParameters()[0] <<" n-tracks: " << _bbarnumber << '\n';
523  const vector< MCParticle * > cdaughters= static_cast< MyVertex * >(verticies->at(1))->__GetMCParticles(); //opera.SelectStableCloseDaughters(chain->Get(1)); //opera.ScanForVertexParticles(bbarverticies->at(1)->getPosition(), 1e-3);
524  _cbarnumber = cdaughters.size();
525  opera.CheckDaughterVisibility(daughters);
526  streamlog_out(MESSAGE) <<"Vertex cbar-quark: " << verticies->at(1)->getParameters()[0] <<" n-tracks: " <<_cbarnumber << '\n';
527  opera.CheckDaughterVisibility(cdaughters);
528  _cbarnumber_f = opera.SelectStableCloseDaughters(chain->Get(1),0,true).size();//opera.CheckDaughterVisibility(cdaughters).size();
529  _bbardistance = MathOperator::getDistance(verticies->at(1)->getPosition(), verticies->at(0)->getPosition());
530  _bbartotalnumber = _cbarnumber + _bbarnumber;
531  _cbarcharge = (int) chain->Get(1)->getCharge();
532  Write(cdaughters, -2);
533  _bbarptmiss = getMissingPt(daughters, cdaughters, verticies->at(0));
534  streamlog_out(MESSAGE) <<"Missing pt for bbar-quark hadron: " << _bbarptmiss << "\n";
535  _bbarmomentum = MathOperator::getModule(chain->Get(0)->getMomentum());
536  _bbaraccuracy = opera.GetAccuracy(chain->Get(0), _aParameter, _bParameter);
537  _cbarmomentum = MathOperator::getModule(chain->Get(1)->getMomentum());
538  _cbaraccuracy = opera.GetAccuracy(chain->Get(1), _aParameter, _bParameter);
539  _cosantiquark = std::cos(MathOperator::getAngles(direction)[1]);
540 
541  }
542  //WriteMisReco(misreco);
543  }
544  void TruthVertexFinder::WriteMisReco(vector< MCParticle * > * particles)
545  {
546 
547  streamlog_out(MESSAGE) << "Misreco: " << _misreconumber << " particles: " << particles->size() << '\n';
548  for (unsigned int i = _misreconumber; i < _misreconumber + particles->size(); i++)
549  {
550  vector< float > direction = MathOperator::getDirection(particles->at(i-_misreconumber)->getMomentum());
551  _misrecotheta[i] = MathOperator::getAngles(direction)[1];
552  _misrecocostheta[i] = std::cos(MathOperator::getAngles(direction)[1]);
553  _misrecomomentum[i] = MathOperator::getModule(particles->at(i-_misreconumber)->getMomentum());
554  _misrecopt[i] = MathOperator::getPt(particles->at(i-_misreconumber)->getMomentum());
555  }
556  _misreconumber += particles->size();
557  }
558  double TruthVertexFinder::getMissingPt(const vector< MCParticle * > & bdaugthers, const vector< MCParticle * > & cdaughters, Vertex * vertex)
559  {
560  const float * position = vertex->getPosition();
561  for (int i = 0; i < 3; i++)
562  {
563  streamlog_out(MESSAGE) << i << ": " << position[i];
564  }
565  streamlog_out(MESSAGE) << '\n';
566  vector< const double * > vectors;
567  for (unsigned int i = 0; i < bdaugthers.size(); i++)
568  {
569  vectors.push_back(bdaugthers[i]->getMomentum());
570  }
571  for (unsigned int i = 0; i < cdaughters.size(); i++)
572  {
573  vectors.push_back(cdaughters[i]->getMomentum());
574  }
575  double missing = MathOperator::getMissingPt(vectors, position);
576  return missing;
577  }
578  void TruthVertexFinder::Write(const vector< MCParticle * > daughters, int v)
579  {
580  float * offset = NULL;
581  float * p = NULL;
582  float * eta = NULL;
583  switch(v)
584  {
585  case 1:
586  offset = _boffsettrack;
587  p = _bptrack;
588  eta = _betatrack;
589  break;
590  case 2:
591  offset = _coffsettrack;
592  p = _cptrack;
593  eta = _cetatrack;
594  break;
595  case -2:
596  offset = _cbaroffsettrack;
597  p = _cbarptrack;
598  eta = _cbaretatrack;
599  break;
600  case -1:
601  offset = _bbaroffsettrack;
602  p = _bbarptrack;
603  eta = _bbaretatrack;
604  break;
605 
606  }
607  int size = daughters.size();
608  for (int i = 0; i < size; i++)
609  {
610  MCParticle * daughter = daughters[i];
611  vector< float > direction = MathOperator::getDirection(daughter->getMomentum());
612  offset[i] = MathOperator::getDistanceTo(ip, direction, daughter->getVertex());
613  //float pt[i] = MathOperator::getPt(daughter->getMomentum());
614  p[i] = MathOperator::getModule(daughter->getMomentum());
615  eta[i] = MathOperator::getAngles(direction)[1];
616  }
617  }
618  void TruthVertexFinder::WriteVertexCollection(LCEvent * evt, vector< Vertex * > * bvertexes, vector< Vertex * > * bbarvertexes)
619  {
620  IMPL::LCCollectionVec * mc = new IMPL::LCCollectionVec ( EVENT::LCIO::VERTEX ) ;
621  if (bvertexes)
622  {
623  for (unsigned int i = 0; i < bvertexes->size(); i++)
624  {
625 
626  mc->addElement(bvertexes->at(i));
627  }
628  }
629  if (bbarvertexes)
630  {
631  for (unsigned int i = 0; i < bbarvertexes->size(); i++)
632  {
633  mc->addElement(bbarvertexes->at(i));
634  }
635  }
636  mc->parameters().setValue("GenVertexTag", _tag);
637  evt->addCollection( mc , _outputcolName ) ;
638  }
639 
640  void TruthVertexFinder::ClearVariables()
641  {
642  _misreconumber = 0;
643  _tag = false;
644  _bstarnumber = 0;
645  _bptmiss = -1.0;
646  _bbarptmiss = -1.0;
647  _cosquark = -2.0;
648  _cosantiquark = -2.0;
649  _totalBcharge = -3.0;
650  _ccharge = -3.0;
651  _cbarcharge = -3.0;
652  _bcharge = -3.0;
653  _bbarcharge = -3.0;
654  _bdistance = -1.0;
655  _bbardistance = -1.0;
656  _bmomentum = -1.0;
657  _bbarmomentum = -1.0;
658  _bbaraccuracy = -1.0;
659  _baccuracy = -1.0;
660  _cmomentum = -1.0;
661  _cbarmomentum = -1.0;
662  _cbaraccuracy = -1.0;
663  _caccuracy = -1.0;
664  _bnumber = -1;
665  _bbarnumber = -1;
666  _cnumber = -1;
667  _cbarnumber = -1;
668  _bbartotalnumber = -1;
669  _btotalnumber = -1;
670  _bnumber_f = -1;
671  _bbarnumber_f = -1;
672  _cnumber_f = -1;
673  _cbarnumber_f = -1;
674  for (int i = 0; i < 2; i++)
675  {
676  _firstVertexDistance[i] = 0.0;
677  _secondVertexDistance[i] = 0.0;
678  _numberOfB0 = 0;
679  }
680  for (int i = 0; i < MAXV; i++)
681  {
682  _bstaroffset[i] = -1.0;
683  _bstarmomentum[i] = -1.0;
684  _bptrack[i] = -1.0;
685  _betatrack[i] = -1.0;
686  _boffsettrack[i] = -1.0;
687  _bbarptrack[i] = -1.0;
688  _bbaretatrack[i] = -1.0;
689  _bbaroffsettrack[i] = -1.0;
690  _cptrack[i] = -1.0;
691  _cetatrack[i] = -1.0;
692  _coffsettrack[i] = -1.0;
693  _cbarptrack[i] = -1.0;
694  _cbaretatrack[i] = -1.0;
695  _cbaroffsettrack[i] = -1.0;
696  _numberOfParticles[i] = -1;
697  _charge[i] = 0;
698  _PDG[i] = 0;
699  _generation[i] = 0;
700  _numberOfParticles[i] = 0;
701  for (int j = 0; j < 3; j++)
702  {
703  _coordinates[i][j] = -100.0;
704  }
705  for (int j = 0; j < MAXV; j++)
706  {
707  _energyOfParticles[i][j] = -1.0;
708  _momentumOfParticles[i][j] = -1.0;
709  _massOfParticles[i][j] = -1.0;
710  _interactionOfParticles[i][j] = -1;
711  }
712  }
713  }
714 
715 
716  void TruthVertexFinder::check( LCEvent * /*evt*/ )
717  {
718  // nothing to check here - could be used to fill checkplots in reconstruction processor
719 
720  }
721 
722 
723  void TruthVertexFinder::end()
724  {
725  if (_writeROOTparameter == 0)
726  {
727  return;
728  }
729  _hfile->cd();
730  _hfile->Write();
731  _hfile->Close();
732 
733  // streamlog_out(MESSAGE) << "TruthVertexFinder::end() " << name()
734  // << " processed " << _nEvt << " events in " << _nRun << " runs "
735  // << std::endl ;
736 
737  }
738 } /* TTbarAnalysis */
TTbarAnalysis::MCOperator::Construct
DecayChain * Construct(std::string name, int pdg, std::vector< PDGTYPE > typeOfProducts)
Definition: MCOperator.cc:49
TTbarAnalysis::TruthVertexFinder::Write
void Write(std::vector< EVENT::Vertex * > *vertices, int &number)
TTbarAnalysis::TruthVertexFinder::AddProngs
void AddProngs(VertexMCOperator &vertexOperator, MCOperator &opera, DecayChain *chain, std::vector< Vertex * > *vertices, std::vector< int > &parameters, IMPL::LCCollectionVec *col=NULL)
Definition: TruthVertexFinder.cc:292
TTbarAnalysis::MCOperator::CheckCompatibility
bool CheckCompatibility(const std::vector< EVENT::MCParticle * > &daughters, EVENT::MCParticle *parent, int plusCharge=0)
Definition: MCOperator.cc:594
TTbarAnalysis::MathOperator::getDirection
static std::vector< float > getDirection(std::vector< int > &vectorPoint1, std::vector< int > &vectorPoint2)
Definition: MathOperator.cc:133
TTbarAnalysis::DecayChain::GetParentPDG
int GetParentPDG() const
Definition: DecayChain.cc:50
TTbarAnalysis::TruthVertexFinder::PrintChain
void PrintChain(std::vector< MCParticle * > *chain)
Definition: TruthVertexFinder.cc:268
TTbarAnalysis::TruthVertexFinder::getMissingPt
double getMissingPt(const std::vector< MCParticle * > &bdaugthers, const std::vector< MCParticle * > &cdaughters, Vertex *vertex)
Definition: TruthVertexFinder.cc:558
TTbarAnalysis::MCOperator::CheckDaughterVisibility
std::vector< EVENT::MCParticle * > CheckDaughterVisibility(const std::vector< EVENT::MCParticle * > &daughters)
Definition: MCOperator.cc:616
TTbarAnalysis::MathOperator::getModule
static float getModule(const std::vector< int > &v)
Definition: MathOperator.cc:52
TTbarAnalysis::TruthVertexFinder::processRunHeader
virtual void processRunHeader(LCRunHeader *run)
Definition: TruthVertexFinder.cc:225
TTbarAnalysis::MCOperator::GetPairParticles
std::vector< EVENT::MCParticle * > GetPairParticles(int pdg)
Definition: MCOperator.cc:417
TTbarAnalysis::aTruthVertexFinder
TruthVertexFinder aTruthVertexFinder
Definition: TruthVertexFinder.cc:7
TTbarAnalysis::TruthVertexFinder::PrintParticle
void PrintParticle(MCParticle *particle)
Definition: TruthVertexFinder.cc:230
TTbarAnalysis::TruthVertexFinder::WriteQuarksCollection
void WriteQuarksCollection(LCEvent *evt, std::vector< MCParticle * > &quarks)
Definition: TruthVertexFinder.cc:279
TTbarAnalysis::TruthVertexFinder::WriteMisReco
void WriteMisReco(std::vector< MCParticle * > *particles)
Definition: TruthVertexFinder.cc:544
TTbarAnalysis::MathOperator::getPt
static float getPt(const double *momentum)
Definition: MathOperator.cc:312
TTbarAnalysis::MCOperator::CheckProcessForPair
std::vector< EVENT::MCParticle * > * CheckProcessForPair(int pdg)
Definition: MCOperator.cc:496
TTbarAnalysis::MathOperator::getAngles
static std::vector< float > getAngles(std::vector< float > &direction)
Definition: MathOperator.cc:97
sistrip::e
static const float e
Definition: PhysicalConstants.h:26
TTbarAnalysis::MCOperator::GetAccuracy
float GetAccuracy(EVENT::MCParticle *particle, float a, float b)
Definition: MCOperator.cc:584
TTbarAnalysis::TruthVertexFinder::processEvent
virtual void processEvent(LCEvent *evt)
Definition: TruthVertexFinder.cc:350
TTbarAnalysis::TruthVertexFinder::WriteVertexCollection
void WriteVertexCollection(LCEvent *evt, std::vector< Vertex * > *bvertexes, std::vector< Vertex * > *bbarvertexes)
Definition: TruthVertexFinder.cc:618
TTbarAnalysis::MCOperator::RefineDecayChain
DecayChain * RefineDecayChain(DecayChain *initial, std::vector< PDGTYPE > typeOfProducts)
Definition: MCOperator.cc:19
sistrip::LCCollectionVec
std::vector< LCCollection * > LCCollectionVec
Definition: SiStripClus.h:55
TTbarAnalysis::MathOperator::getMissingPt
static double getMissingPt(std::vector< const double * > &vectors, const float *target)
Definition: MathOperator.cc:349
TTbarAnalysis::TruthVertexFinder::GetAsymmetry
void GetAsymmetry(std::vector< MCParticle * > &particles)
Definition: TruthVertexFinder.cc:450
TTbarAnalysis::VertexMCOperator::Construct
std::vector< EVENT::Vertex * > * Construct(DecayChain *chain)
Definition: VertexMCOperator.cc:17
TTbarAnalysis::TruthVertexFinder
This processor is designed to extruct a secondary/// vertices from collection of generated particles...
Definition: TruthVertexFinder.hh:48
TTbarAnalysis::MathOperator::getDistance
static float getDistance(const double *start, const double *end)
Definition: MathOperator.cc:22
TTbarAnalysis::MathOperator::getDistanceTo
static float getDistanceTo(const std::vector< int > &vectorPoint1, const std::vector< float > &vector1, const std::vector< int > *pointOfLine)
Definition: MathOperator.cc:188
TTbarAnalysis::DecayChain::Get
EVENT::MCParticle * Get(int i)
Definition: DecayChain.cc:23
TTbarAnalysis::TruthVertexFinder::check
virtual void check(LCEvent *evt)
Definition: TruthVertexFinder.cc:716
TTbarAnalysis::VertexMCOperator::AddProngs
void AddProngs(EVENT::Vertex *vertex, std::vector< EVENT::MCParticle * > &particles, bool usingRelation=false)
Definition: VertexMCOperator.cc:55
TTbarAnalysis::MCOperator::SelectStableCloseDaughters
std::vector< EVENT::MCParticle * > SelectStableCloseDaughters(EVENT::MCParticle *parent, int excludePDG=0, bool selectReco=false, std::vector< EVENT::MCParticle * > *misReconstructed=NULL)
Definition: MCOperator.cc:529