All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Groups
TrueJet.cc
Go to the documentation of this file.
1 #include "TrueJet.h"
2 #include <stdlib.h>
3 #include <math.h>
4 #include <iostream>
5 #include <iomanip>
6 
7 // ----- include for verbosity dependend logging ---------
8 #include "marlin/VerbosityLevels.h"
9 
10 #ifdef MARLIN_USE_AIDA
11 #include <marlin/AIDAProcessor.h>
12 #include <AIDA/IHistogramFactory.h>
13 #include <AIDA/ICloud1D.h>
14 //#include <AIDA/IHistogram1D.h>
15 #endif // MARLIN_USE_AIDA
16 
17 
18 using namespace lcio ;
19 using namespace marlin ;
20 
21 struct MCPpyjet : LCIntExtension<MCPpyjet> {} ;
22 struct TJindex : LCIntExtension<TJindex> {} ;
23 LCRelationNavigator* reltrue =0;
24 
25 TrueJet aTrueJet ;
26 
27 TrueJet::TrueJet() : Processor("TrueJet") {
28 
29  // modify processor description
30  _description = "TrueJet does whatever it does ..." ;
31 
32  // register steering parameters: name, description, class-variable, default value
33 
34 
35  // Inputs: MC-particles, Reco-particles, the link between the two
36 
37  registerInputCollection( LCIO::MCPARTICLE,
38  "MCParticleCollection" ,
39  "Name of the MCParticle collection" ,
40  _MCParticleColllectionName ,
41  std::string("MCParticlesSkimmed") );
42 
43  registerInputCollection( LCIO::RECONSTRUCTEDPARTICLE,
44  "RecoParticleCollection" ,
45  "Name of the ReconstructedParticles input collection" ,
46  _recoParticleCollectionName ,
47  std::string("PandoraPFOs") ) ;
48 
49  registerInputCollection( LCIO::LCRELATION,
50  "RecoMCTruthLink",
51  "Name of the RecoMCTruthLink input collection" ,
52  _recoMCTruthLink,
53  std::string("RecoMCTruthLink") ) ;
54 
55  // Outputs: True jets (as a recoparticle, will be the sum of the _reconstructed particles_
56  // created by the true particles in each true jet, in the RecoMCTruthLink sense.
57  // link jet-to-reco particles, link jet-to-MC-particles.
58 
59  registerOutputCollection( LCIO::RECONSTRUCTEDPARTICLE,
60  "TrueJets" ,
61  "Name of the TrueJetCollection output collection" ,
62  _trueJetCollectionName ,
63  std::string("TrueJets") ) ;
64 
65  registerOutputCollection( LCIO::RECONSTRUCTEDPARTICLE,
66  "FinalColourNeutrals" ,
67  "Name of the FinalColourNeutralCollection output collection" ,
68  _finalColourNeutralCollectionName ,
69  std::string("FinalColourNeutrals") ) ;
70 
71  registerOutputCollection( LCIO::RECONSTRUCTEDPARTICLE,
72  "InitialColourNeutrals" ,
73  "Name of the InitialColourNeutralCollection output collection" ,
74  _initialColourNeutralCollectionName ,
75  std::string("InitialColourNeutrals") ) ;
76 
77 
78  registerOutputCollection( LCIO::LCRELATION,
79  "TrueJetPFOLink" ,
80  "Name of the TrueJetPFOLink output collection" ,
81  _trueJetPFOLink,
82  std::string("TrueJetPFOLink") ) ;
83 
84  registerOutputCollection( LCIO::LCRELATION,
85  "TrueJetMCParticleLink" ,
86  "Name of the TrueJetMCParticleLink output collection" ,
87  _trueJetMCParticleLink,
88  std::string("TrueJetMCParticleLink") ) ;
89 
90  registerOutputCollection( LCIO::LCRELATION,
91  "FinalElementonLink" ,
92  "Name of the FinalElementonLink output collection" ,
93  _finalElementonLink,
94  std::string("FinalElementonLink") ) ;
95 
96  registerOutputCollection( LCIO::LCRELATION,
97  "InitialElementonLink" ,
98  "Name of the InitialElementonLink output collection" ,
99  _initialElementonLink,
100  std::string("InitialElementonLink") ) ;
101 
102  registerOutputCollection( LCIO::LCRELATION,
103  "FinalColourNeutralLink" ,
104  "Name of the FinalColourNeutralLink output collection" ,
105  _finalColourNeutralLink,
106  std::string("FinalColourNeutralLink") ) ;
107 
108  registerOutputCollection( LCIO::LCRELATION,
109  "InitialColourNeutralLink" ,
110  "Name of the InitialColourNeutralLink output collection" ,
111  _initialColourNeutralLink,
112  std::string("InitialColourNeutralLink") ) ;
113  // steering
114 
115  registerProcessorParameter("Whizard1",
116  "true: Input has MCParticles in Whizard1 convention ",
117  _whiz1,
118  bool(false)
119  );
120  }
121 
122 
123 
124 void TrueJet::init() {
125 
126  streamlog_out(DEBUG7) << " init called " << std::endl ;
127 
128  // usually a good idea to
129  printParameters() ;
130 
131  _nRun = 0 ;
132  _nEvt = 0 ;
133 
134 }
135 
136 
137 void TrueJet::processRunHeader( LCRunHeader* /*run*/) {
138 
139  _nRun++ ;
140 }
141 
142 
143 
144 void TrueJet::processEvent( LCEvent * event ) {
145 
146 
147  // For the concepts this method works with, see the comments in TrueJet.h
148 
149  evt=event;
150  streamlog_out(WARNING) << " processing event: " << evt->getEventNumber()
151  << " in run: " << evt->getRunNumber() << std::endl ;
152  streamlog_out(MESSAGE4) << " =====================================" << std::endl;
153 
154 
155  LCCollection* mcpcol = NULL;
156  try{
157  mcpcol = evt->getCollection( _MCParticleColllectionName );
158  }
159  catch( lcio::DataNotAvailableException& e )
160  {
161  streamlog_out(WARNING) << _MCParticleColllectionName << " collection not available" << std::endl;
162  mcpcol = NULL;
163  return;
164  }
165  LCCollection* rmclcol = NULL;
166  try{
167  rmclcol = evt->getCollection( _recoMCTruthLink );
168  }
169  catch( lcio::DataNotAvailableException& e )
170  {
171  streamlog_out(MESSAGE4) << _recoMCTruthLink << " collection not available" << std::endl;
172  rmclcol = NULL;
173  }
174 
175  // TODO: in w2 : get polarization1 &2 ; count sum of W:s and L/R : if this is odd, then an odd number
176  // of leptons are expected ( one beam-remnant + zero or more pairs ). For now, just skip the
177  // warning about odd number of leptons.
178 
179  if( mcpcol != NULL ){
180 
181  if ( rmclcol != NULL ) { reltrue = new LCRelationNavigator( rmclcol ); }
182 
183  // recast the MCParticles into a PYJETS look-alike, needed for the logic to work:
184 
185 
186  mcp_pyjets.reserve(4000);
187  getPyjets(mcpcol );
188 
189  if ( _higgs_to_glue_glue ) {
190  streamlog_out(DEBUG4) << " Higgs-> gluon gluon event : TrueJet will not work correctly " << std::endl;
191  if ( reltrue ) delete reltrue ;
192  _nEvt++ ;
193  LCCollectionVec* jet_vec = new LCCollectionVec(LCIO::RECONSTRUCTEDPARTICLE ) ;
194  evt->addCollection( jet_vec, _trueJetCollectionName);
195  return;
196  }
197 
198 
199 
200  njet=0 ; current_jet=0 ; nstr=0 ;
201  for ( int kk=1 ; kk <= 25 ; kk++ ) {
202  boson[kk]=0;
203  elementon[kk]=0;
204  fafp_last[kk]=0;
205  fafp_boson[kk]=0;
206  fafp[kk]=0;
207  nfsr[kk]=0 ;
208  }
209  bool seen[4001];
210  for ( int kk=1 ; kk <= 4000 ; kk++ ) {
211  jet[kk]=0;
212  companion[kk] = 0;
213  seen[kk] = false;
214  }
215 
216  // the actual work happens in the following routines. They are each responible to
217  // find jets stemming from different sources, as indicated by their names.
218 
219  cluster() ;
220 
221  true_lepton() ;
222 
223  photon();
224 
225  isr() ;
226 
227  nstr=0;
228  for ( int k_py=1 ; k_py <= nlund ; k_py++ ) {
229  streamlog_out(DEBUG0) << " Checking string at " << k_py << " k[k_py][2], jet[k_py] and k[k_py][1] : " <<
230  " " << k[k_py][2] << " , " << jet[k_py] << " , " << k[k_py][1] << std::endl;
231  if ( k[k_py][2] == 92 && jet[k_py] == 0 && k[k_py][1] < 30 ) {
232  streamlog_out(DEBUG4) << " string found at " << k_py <<std::endl;
233  nstr++;
234  }
235  }
236 
237  njet = 2*nstr+n_hard_lepton+2*nclu+nisr+nphot;
238 
239  if ( nstr != 0 ) {
240  string();
241  }
242 
243  n_beam_jet=0;
244  for (int k_py=1 ; k_py<=nlund ; k_py++) {
245  if ( k[k_py][1] >= 30 ) {
246  if ( n_beam_jet == 0 ) {
247  n_beam_jet=1;
248  }
249  jet[k_py]=njet+n_beam_jet;
250  }
251  }
252 
253  njet=njet+n_beam_jet;
254 
255  grouping();
256 
257  // At this point, the job is done. All the rest is printouts, putting things into collections and setting up navigators.
258  // Non-debuggung code is beween ******:s and =========:s
259 
260  // FIXME : determine form the event header is an odd or even number of leptons are expected
261  // if ( ! ((n_mixed==0) && njet==2*nstr+n_hard_lepton+2*nclu+nisr+n_beam_jet &&
262  // n_jetless==0 && n_hard_lepton%2 == 0 )) {
263 
264  if ( ! ((n_mixed==0) && njet==2*nstr+n_hard_lepton+2*nclu+nisr+n_beam_jet+nphot &&
265  n_jetless==0 )) {
266  if ( ! _higgs_to_glue_glue ) {
267  streamlog_out(ERROR) << " inconsiency in jet finding : " << std::endl;
268  streamlog_out(ERROR) << " n_mixed/ njet/ 2*nstr+n_hard_lepton+2*nclu+nisr+n_beam_jet / n_jetless / n_hard_lepton: " << std::endl;
269  streamlog_out(ERROR) << n_mixed <<" / " << njet<< " / " <<2*nstr+n_hard_lepton+2*nclu+nisr+n_beam_jet<< " / "
270  << n_jetless << " / " << n_hard_lepton << std::endl;
271  streamlog_out(ERROR) << " Event: " << evt->getEventNumber() << ", run: " << evt->getRunNumber() << std::endl ;
272  streamlog_out(ERROR) << std::endl ;
273  }
274  }
275 
276  streamlog_out(DEBUG3) << " HEPEVT relation table up with jet assignment " <<std::endl;
277  streamlog_out(DEBUG3) << " line status pdg parent first last px py "
278  << " pz E M jet companion type" << std::endl;
279  streamlog_out(DEBUG3) << " daughter "
280  << " jet" << std::endl;
281 
282  for ( int k_py=1 ; k_py <= nlund ; k_py++){
283  streamlog_out(DEBUG3) << std::setw(7) << k_py << std::setw(7) << k[k_py][1] <<
284  std::setw(12) << k[k_py][2] << std::setw(7) << k[k_py][3] <<
285  std::setw(7) << k[k_py][4] << std::setw(7) << k[k_py][5] <<
286  std::setw(12) << p[k_py][1] <<
287  std::setw(12) << p[k_py][2] << std::setw(12) << p[k_py][3] <<
288  std::setw(12) << p[k_py][4] << std::setw(12) << p[k_py][5] <<
289 
290  std::setw(7) << jet[k_py] <<
291  std::setw(7) << companion[abs(jet[k_py])] << std::setw(7) << type[abs(jet[k_py])]<<std::endl;
292  }
293  streamlog_out(DEBUG3) << " jet fafp-beg qrk/lept fafp-end fafp-boson type dijet-beg dijet-end boson" << std::endl;
294  for ( int k_py=1 ; k_py <= njet ; k_py++){
295  streamlog_out(DEBUG3) << std::setw(10) << k_py << std::setw(10) << fafp[k_py] <<
296  std::setw(10) << elementon[k_py] << std::setw(10) << fafp_last[k_py] <<
297  std::setw(10) << fafp_boson[k_py] << std::setw(10) << type[k_py] << std::setw(10) << dijet_begining[k_py] <<
298  std::setw(10) << dijet_end[k_py] << std::setw(10) << boson[k_py]<<std::endl;
299  }
300 
301 
302  // now for py-partic i I know which true jet it belongs to.
303  // mcp_pyjets tells me which MCParticle corresponds to each pyjets one.
304  // I then loop all jets creating the true_jet object for each one, loop
305  // all mcparticles to find which ones
306  // belong to the current jet, find if they contribute to a reco-particle, and
307  // if so, add the reco energy etc. to the true_jet-object. While loopoing, the
308  // links are set up, as well.
309 
310  //*****************************
311  // create the navigators:
312 
313  LCCollection* tjrcol = 0;
314  LCCollection* tjtcol = 0;
315  LCRelationNavigator truejet_pfo_Nav(LCIO::RECONSTRUCTEDPARTICLE , LCIO::RECONSTRUCTEDPARTICLE ) ;
316  LCRelationNavigator truejet_truepart_Nav(LCIO::RECONSTRUCTEDPARTICLE , LCIO::MCPARTICLE ) ;
317 
318 
319 
320  // create the collection-vector that will contain the true-jet objects:
321 
322  LCCollectionVec* jet_vec = new LCCollectionVec(LCIO::RECONSTRUCTEDPARTICLE ) ;
323 
324  //============================
325 
326  double str_tmom[3]={0}, str_mom[3]={0}, str_tmomS[3]={0} , str_tE=0 , str_E=0 , str_tES=0;
327 
328  streamlog_out(DEBUG8) << " Number of jets found : " << njet << std::endl;
329 
330  double tmomS[25][3]={0.} ;
331  double tES[25]={0.};
332  int pid_type[25]={0} ;
333  for ( int i_jet=1; i_jet<=njet ; i_jet++ ) { // jet-loop
334 
335  //*****************************
336  ReconstructedParticleImpl* true_jet = new ReconstructedParticleImpl;
337 
338  pid_type[i_jet] = -type[i_jet] ;
339 
340  jet_vec->addElement(true_jet);
341  //============================
342 
343  streamlog_out(DEBUG1) << std::endl;
344  streamlog_out(DEBUG1) << " Following jet " << i_jet << " ( "<< true_jet << ")" << std::endl;
345  streamlog_out(DEBUG1) << " ============================= " << std::endl;
346  streamlog_out(DEBUG1) << " HEPEVT relation table with jet assignment and pfo(s), if any " <<std::endl;
347  streamlog_out(DEBUG1) << " line mcp status pdg parent first last jet pfo(s)" << std::endl;
348  streamlog_out(DEBUG1) << " daughter daugther " << std::endl;
349  }
350 
351  ReconstructedParticleImpl* true_jet ;
352  for ( int k_py=1 ; k_py<=nlund ; k_py++){ // pyjets loop
353  int i_jet = abs(jet[k_py]);
354  if ( i_jet > 0 ) {
355 
356  true_jet = dynamic_cast<ReconstructedParticleImpl*>(jet_vec->getElementAt(i_jet-1));
357 
358  // OK, add jet-to-true link
359 
360  streamlog_out(DEBUG1) << std::setw(10) << k_py << " (" <<mcp_pyjets[k_py] << ")" << std::setw(10) << k[k_py][1]%30 <<
361  std::setw(10) << k[k_py][2] << std::setw(10) << k[k_py][3] <<
362  std::setw(10) << k[k_py][4] << std::setw(10) << k[k_py][5] << std::setw(10) << jet[k_py] ;
363 
364 
365  //*****************************
366  if ( k[k_py][1] == 1 && k[k_py][2] == 22 && jet[k[k_py][3]] < 0 && k[k[k_py][3]][2] != 22 ) {
367  // to be able to find FSRs, set weight to -ve for them
368  truejet_truepart_Nav.addRelation( true_jet, mcp_pyjets[k_py], ( jet[k_py]>0 ? -(k[k_py][1])%30 : 0.0 )) ;
369  } else {
370  truejet_truepart_Nav.addRelation( true_jet, mcp_pyjets[k_py], ( jet[k_py]>0 ? k[k_py][1]%30 : 0.0 )) ;
371  }
372  LCObjectVec recovec;
373  if ( rmclcol != NULL ) {
374  recovec = reltrue->getRelatedFromObjects( mcp_pyjets[k_py]);
375  }
376  //============================
377 
378  // (maybe this will be needed: If the current assumption that also MCPartiles created in
379  // simulation can in a consitent way be put into pyjets turns out to be wrong, then
380  // getPyjets should only load gen. stat. /= 0 particles, and then code below
381  // would be needed to assign seen gen. stat. = 0 particles to jets based on their
382  // gen. stat. /= 0 ancestor ?!)
383  // if ( recovec.size() == 0 && k[k_py][1] == 1 ) { // not reconstructed stable particle
384  // if ( abs(k[k_py][2]) != 12 && abs(k[k_py][2]) != 14 && abs(k[k_py][2]) != 16 ) { // not neutrino
385  // // mcp_pyjets[k_py].getDaughter ... etc . If any decendant is seen, make association
386  // }
387  // }
388 
389  if ( recovec.size() > 0 ) { // if reconstructed
390  double mom2[3] ;
391  double E,q ;
392 
393  // true quanaties for this jet (all true seen)
394 
395  if ( ! seen[k[k_py][3]] ) {
396  // add to true-of-seen only if no ancestor of this true has already been counted
397  tmomS[i_jet][0]+=p[k_py][1] ;
398  tmomS[i_jet][1]+=p[k_py][2] ;
399  tmomS[i_jet][2]+=p[k_py][3] ;
400  tES[i_jet]+=p[k_py][4] ;
401  } else {
402  streamlog_out(DEBUG3) << " Ancestor of true particle " << k_py << " already seen, not added again to true-of-seen " << std::endl;
403  }
404  seen[k_py] = true ;
405  streamlog_out(DEBUG1) << " recovec size is " << recovec.size() ;
406  int last_winner = i_jet ;
407  for ( unsigned i_reco=0 ; i_reco<recovec.size() ; i_reco++ ) { // reco-of-this-true loop
408 
409  // add things of this reco-particle to the current jet:
410 
411  //*****************************
412  ReconstructedParticle* reco_part = dynamic_cast<ReconstructedParticle*>(recovec[i_reco]);
413 
414  if (truejet_pfo_Nav.getRelatedFromObjects(reco_part).size() == 0 ) { // only if not yet used
415  streamlog_out(DEBUG3) << " recopart " << i_reco ;
416  bool split_between_jets = false;
417  int winner=i_jet , wgt_trk[26]={0} , wgt_clu[26]={0} ;
418 
419  LCObjectVec recomctrues = reltrue->getRelatedToObjects(reco_part);
420  static FloatVec recomctrueweights;
421  recomctrueweights = reltrue->getRelatedToWeights(reco_part);
422  streamlog_out(DEBUG3) << " mctrues of this reco " << recomctrues.size() << std::endl ;
423  for ( unsigned k_mcp_of_reco=0 ; k_mcp_of_reco<recomctrues.size() ; k_mcp_of_reco++ ) {
424  MCParticle* an_mcp = dynamic_cast<MCParticle*>(recomctrues[k_mcp_of_reco]);
425  int jetoftrue=jet[an_mcp->ext<MCPpyjet>()];
426  if ( jetoftrue != i_jet ) split_between_jets = true;
427  wgt_trk[jetoftrue]+=int(recomctrueweights[k_mcp_of_reco])%10000 ;
428  wgt_clu[jetoftrue]+=int(recomctrueweights[k_mcp_of_reco])/10000 ;
429  streamlog_out(DEBUG1) << " weights for jet " << jetoftrue << " : " << int(recomctrueweights[k_mcp_of_reco])
430  << " " << int(recomctrueweights[k_mcp_of_reco])%10000
431  << " " << int(recomctrueweights[k_mcp_of_reco])/10000 << std::endl ;
432  }
433  if ( split_between_jets ) {
434  double wgt_trk_max=0, wgt_clu_max=0 ;
435  int imax_trk=0, imax_clu=0;
436 
437  for ( int j_jet=1 ; j_jet <= njet ; j_jet++ ) {
438  streamlog_out(DEBUG3) << " jetweights for " << j_jet << " " << wgt_trk[j_jet] << " " << wgt_clu[j_jet] ;
439  if ( wgt_trk[j_jet] > wgt_trk_max ) {
440  imax_trk= j_jet ; wgt_trk_max = wgt_trk[j_jet] ;
441  }
442  if ( wgt_clu[j_jet] > wgt_clu_max ) {
443  imax_clu= j_jet ; wgt_clu_max = wgt_clu[j_jet] ;
444  }
445  }
446  streamlog_out(DEBUG3) << " " << imax_clu << " " << imax_trk << " " << wgt_clu_max << " " << wgt_trk_max ;
447  if ( imax_clu == imax_trk || wgt_trk_max >= 750 ) {
448  winner = imax_trk;
449  } else if ( wgt_trk_max == 0 ) {
450  winner = imax_clu;
451  } else {
452  streamlog_out(DEBUG3) << " was ? " << imax_clu << " " << imax_trk << " " << wgt_clu_max << " " << wgt_trk_max << std::endl;
453  for ( int j_jet=1 ; j_jet <= njet ; j_jet++ ) {
454  streamlog_out(DEBUG3) << " was jetweights for " << j_jet << " " << wgt_trk[j_jet] << " " << wgt_clu[j_jet] << std::endl;
455  }
456  // probably the cluster is somewhat more reliable in this case. Anyways, this happens very rarely (for 3 PFOs in 6400 4f_had...)
457  if ( imax_clu > 0 ) {
458  winner = imax_clu;
459  } else {
460  // ... as a last resort ...
461  winner = imax_trk ;
462  }
463  }
464  }
465 
466  if ( winner > 0 && winner != last_winner ) true_jet = dynamic_cast<ReconstructedParticleImpl*>(jet_vec->getElementAt(winner-1));
467 
468  streamlog_out(DEBUG3) << " and the winner is " << winner << " ( i_jet is " << i_jet << " ) " ;
469  if ( winner > 0 ) {
470  last_winner = winner ;
471 
472  const double* mom = reco_part->getMomentum();
473  double psq=0.;
474  for ( int xyz=0 ; xyz<3 ; xyz++) {
475  mom2[xyz]=mom[xyz]+true_jet->getMomentum()[xyz] ;
476  psq+= mom2[xyz]*mom2[xyz];
477  }
478 
479  E = reco_part->getEnergy()+true_jet->getEnergy();
480  q = reco_part->getCharge()+true_jet->getCharge();
481  true_jet->setCharge(q);
482  true_jet->setEnergy(E);
483  true_jet->setMass(sqrt(E*E-psq));
484  true_jet->setMomentum(mom2);
485  streamlog_out(DEBUG1) << " " << reco_part ;
486  pid_type[winner] = type[winner] ;
487 
488  true_jet->addParticle(reco_part);
489 
490  // add jet-to-reco (and v.v.) link (good to have both ways easily, reco-particle member of the
491  // true jet object just gives jet->particle, not particle->jet!)
492 
493  truejet_pfo_Nav.addRelation( true_jet, reco_part , 1.0 );
494  //============================
495 
496  } else {
497  // FIXME: a follow on error for h->glueglue: If a particle from a non-higgs jet gets reconstructed into a PFO where
498  // the majority contribution is from a higgs jet (which are not found in h->glueglue, and hence has jet# 0 ).
499  if ( ! _higgs_to_glue_glue ) {
500  streamlog_out(WARNING) << " reco particle " << i_reco << " of pythia-particle " << k_py
501  << " has weights = 0 to ALL MCPs ??? " << std::endl;
502  streamlog_out(WARNING) << " Event: " << evt->getEventNumber() << ", run: " << evt->getRunNumber() << std::endl ;
503  }
504  }
505  } // end if not yet used
506  else {
507  streamlog_out(DEBUG1) << " " << reco_part << " seen more than once ! " << std::endl;
508  }
509  } // end reco-of-this-true loop
510 
511  } else { // end if reconstructed
512  if ( seen[k[k_py][3]] ) {
513  seen[k_py] = true ; // this makes sure that in any decay-chain, only the first seen particle will includes in true-of-seen
514  streamlog_out(DEBUG2) << " Ancestor of true particle " << k_py << " was seen, so particle " << k_py
515  << " is marked as seen (even if it wasn't seen itself " << std::endl;
516  }
517  }
518  streamlog_out(DEBUG1) << std::endl;
519 
520  } // end if i_jet > 0
521 
522  } // end pyjets loop
523 
524  for ( int i_jet=1; i_jet<=njet ; i_jet++ ) { // jet-loop
525 
526  true_jet = dynamic_cast<ReconstructedParticleImpl*>(jet_vec->getElementAt(i_jet-1));
527  ParticleIDImpl* pid = new ParticleIDImpl;
528  pid->setPDG(k[elementon[i_jet]][2]);
529  pid->setType(pid_type[i_jet]);
530  true_jet->addParticleID(pid);
531  true_jet->ext<TJindex>()=i_jet;
532 
533  }
534 
535  for ( int i_jet=1; i_jet<=njet ; i_jet++ ) { // jet-loop
536 
537  true_jet = dynamic_cast<ReconstructedParticleImpl*>(jet_vec->getElementAt(i_jet-1));
538 
539  const double* mom = true_jet->getMomentum();
540 
541  streamlog_out(DEBUG5) << std::endl;
542  streamlog_out(DEBUG9) << " summary " << i_jet << " " << type[i_jet]%100 << " " << true_jet->getEnergy() << " "
543  << tE[i_jet] << " " << tES[i_jet] << std::endl ;
544  streamlog_out(DEBUG5) << " jet " <<std::setw(4)<< i_jet << " (seen) p : "
545  << mom[0] << " " << mom[1] << " " << mom[2]
546  << " " << " E " << true_jet->getEnergy()<< std::endl;
547  streamlog_out(DEBUG5) << " " << " (true) p : "
548  << tmom[i_jet][0] << " " << tmom[i_jet][1] << " " << tmom[i_jet][2]
549  << " " << " E " << tE[i_jet] << std::endl;
550  streamlog_out(DEBUG5) << " " << " (true of seen) p : "
551  << tmomS[i_jet][0] << " " << tmomS[i_jet][1] << " " << tmomS[i_jet][2]
552  << " " << " E " << tES[i_jet] << std::endl;
553  streamlog_out(DEBUG5) << " ancestor 1" << std::setw(4)<<elementon[i_jet] << " p : "
554  << p[elementon[i_jet]][1] << " " << p[elementon[i_jet]][2] << " " << p[elementon[i_jet]][3]
555  << " " << " E " << p[elementon[i_jet]][4]<< std::endl;
556  streamlog_out(DEBUG5) << " ancestor 2" << std::setw(4)<<fafp_last[i_jet] << " p : "
557  << p[fafp_last[i_jet]][1] << " " << p[fafp_last[i_jet]][2] << " " << p[fafp_last[i_jet]][3]
558  << " " << " E " << p[fafp_last[i_jet]][4]<< std::endl;
559  streamlog_out(DEBUG5) << " ancestor 3" << std::setw(4)<<fafp[i_jet] << " p : "
560  << p[fafp[i_jet]][1] << " " << p[fafp[i_jet]][2] << " " << p[fafp[i_jet]][3]
561  << " " << " E " << p[fafp[i_jet]][4]<< std::endl;
562  streamlog_out(DEBUG5) << std::endl;
563 
564  streamlog_out(DEBUG5) << " Character of jet " << i_jet << " : "<<type[i_jet]%100 <<" (1=string, 2=lepton , "
565  <<"3=cluster, 4=isr, 5=overlay, 6=M.E. photon)." << std::endl;
566  streamlog_out(DEBUG5) << " Jet from boson ? " << (type[i_jet]>=100 ? type[i_jet]/100 : 0)
567  << " (0 not from boson, !=0 from boson radiated of that jet)." << std::endl;
568  streamlog_out(DEBUG5) << " No. of FSRs : " << nfsr[i_jet] << std::endl;
569  streamlog_out(DEBUG5) << " Jet has no energy ? " << (tE[i_jet]<0.001) << std::endl;
570 
571  streamlog_out(DEBUG5) << std::endl;
572 
573  str_E+= true_jet->getEnergy(); str_tE+= tE[i_jet]; str_tES+= tES[i_jet];
574  for ( int i=0 ; i<3 ; i++ ) {
575  str_mom[i]+=mom[i];
576  str_tmom[i]+=tmom[i_jet][i];
577  str_tmomS[i]+=tmomS[i_jet][i];
578  }
579 
580  int odd_even = 0 ;
581  if ( type[i_jet]%100 == 3 && n_hard_lepton%2 == 1 ) odd_even = 1;
582 
583  if ( i_jet%2 == odd_even && type[i_jet]%100 <= 3 ) {
584  streamlog_out(DEBUG6) << " Mass of jet " << i_jet-1 << " and " << i_jet << std::endl;
585  double masssq = str_E*str_E - str_mom[0]*str_mom[0] - str_mom[1]*str_mom[1] - str_mom[2]*str_mom[2] ;
586  if ( masssq > 0. ) {
587  streamlog_out(DEBUG6) << " Seen " << sqrt(masssq) << std::endl;
588  }
589  double masssq_t = str_tE*str_tE - str_tmom[0]*str_tmom[0] - str_tmom[1]*str_tmom[1] - str_tmom[2]*str_tmom[2] ;
590  if ( masssq_t > 0. ) {
591  streamlog_out(DEBUG6) << " true " << sqrt(masssq_t) << std::endl;
592  }
593  masssq = str_tES*str_tES - str_tmomS[0]*str_tmomS[0] - str_tmomS[1]*str_tmomS[1] - str_tmomS[2]*str_tmomS[2] ;
594  if ( masssq > 0. ) {
595  streamlog_out(DEBUG6) << " true of seen " << sqrt(masssq) << std::endl;
596  }
597  if ( type[i_jet]%100 == 3 ) {
598  // cluster - the 91 object does NOT have the mass of the following physical states,
599  // so we sum up the descendants - in the vast majority there is only one
600  double clu_mass=0 , clu_E= 0., clu_P[4]= {0,0,0};
601  for ( int j_py=k[k[elementon[i_jet]][4]][4] ; j_py<=k[k[elementon[i_jet]][4]][5] ; j_py++ ) {
602  if ( k[j_py][3] == k[elementon[i_jet]][4] ) {
603  clu_E += p[j_py][4] ; clu_P[1] += p[j_py][1] ; clu_P[2] += p[j_py][2] ; clu_P[3] += p[j_py][3] ;
604  }
605  }
606  clu_mass=sqrt(clu_E*clu_E-clu_P[1]*clu_P[1]-clu_P[2]*clu_P[2]-clu_P[3]*clu_P[3] );
607  streamlog_out(DEBUG6) << " string mass : " << clu_mass << std::endl;
608  } else if ( type[i_jet]%100 == 2 ) {
609  double dilept_mass=0 , dilept_E= 0., dilept_P[4]= {0,0,0};
610  for ( int j_jet=i_jet-1 ; j_jet<=i_jet ; j_jet++ ) {
611  dilept_E += p[elementon[j_jet]][4] ;
612  dilept_P[1] += p[elementon[j_jet]][1] ;
613  dilept_P[2] += p[elementon[j_jet]][2] ;
614  dilept_P[3] += p[elementon[j_jet]][3] ;
615  }
616  dilept_mass=sqrt(dilept_E*dilept_E-dilept_P[1]*dilept_P[1]-dilept_P[2]*dilept_P[2]-dilept_P[3]*dilept_P[3] );
617  streamlog_out(DEBUG6) << " string mass : " << dilept_mass << std::endl;
618  } else {
619  // string, just take the number from the 92-object
620  streamlog_out(DEBUG6) << " string mass : " << p[k[elementon[i_jet]][4]][5] << std::endl;
621  }
622 
623  streamlog_out(DEBUG6) << std::endl;
624  str_tE=0; str_tES=0; str_E=0;
625  for ( int i=0 ; i<3 ; i++ ) {
626  str_tmom[i]=0. ; str_tmomS[i]=0. ; str_mom[i]=0. ;
627  }
628  if ( masssq_t > 0. ) {
629  if ( abs( sqrt(masssq_t)- p[k[elementon[i_jet]][4]][5] )/ p[k[elementon[i_jet]][4]][5] > 0.007 ) {
630  if ( type[i_jet]%100 == 1 && type[i_jet-1]%100 == 1 && nfsr[i_jet]+nfsr[i_jet-1]==0) {
631  if ( abs( tE[i_jet]+tE[i_jet-1]- p[k[elementon[i_jet]][4]][4] )/ p[k[elementon[i_jet]][4]][4] > 0.001 ) {
632  streamlog_out(ERROR) << " bad match M (sum/initial) " << " " << sqrt(masssq_t)
633  << " / " << p[k[elementon[i_jet]][4]][5] << std::endl;
634  streamlog_out(ERROR) << " E (sum/initial) " << " " << tE[i_jet]+tE[i_jet-1]
635  << " / " << p[k[elementon[i_jet]][4]][4] << std::endl;
636 
637  streamlog_out(DEBUG9) << "list HEPEVT relation table up with jet assignment " <<std::endl;
638  streamlog_out(DEBUG9) << "list line status pdg parent first last px py pz "
639  <<" E M jet companion type" << std::endl;
640  streamlog_out(DEBUG9) << "list daughter "
641  <<" jet" << std::endl;
642  for ( int i_py=1 ; i_py <= nlund ; i_py++){
643  streamlog_out(DEBUG9) <<"list "<< std::setw(7) << i_py << std::setw(7) << k[i_py][1] <<
644  std::setw(12) << k[i_py][2] << std::setw(7) << k[i_py][3] <<
645  std::setw(7) << k[i_py][4] << std::setw(7) << k[i_py][5] <<
646  std::setw(12) << p[i_py][1] <<
647  std::setw(12) << p[i_py][2] << std::setw(12) << p[i_py][3] <<
648  std::setw(12) << p[i_py][4] << std::setw(12) << p[i_py][5] <<
649 
650  std::setw(7) << jet[i_py] <<
651  std::setw(7) << companion[abs(jet[i_py])] << std::setw(7) << type[abs(jet[i_py])]<<std::endl;
652  }
653  streamlog_out(DEBUG9) << "list of individual particles with bad parent/kid energy " <<std::endl;
654  for ( int i_py=1 ; i_py<=nlund ; i_py++ ) {
655  if ( jet[i_py] == i_jet || jet[i_py] == i_jet -1 ) {
656  if ( k[i_py][1] == 11 ) {
657  double e_kid=0;
658  for ( int jj= k[i_py][4] ; jj <= k[i_py][5] ; jj++ ) {
659  e_kid+=p[jj][4];
660  }
661  if ( (abs(e_kid - p[i_py][4])/ p[i_py][4] ) > 0.001 ) {
662  streamlog_out(DEBUG9) << i_py << " " << k[i_py][4] << " " << k[i_py][5] << " "
663  << e_kid << " " << p[i_py][4] << std::endl ;
664  }
665  }
666  }
667  }
668  streamlog_out(WARNING) << " jet fafp-beg qrk/lept fafp-end fafp-boson "
669  <<"type dijet-beg dijet-end boson" << std::endl;
670  for ( int j_jet=i_jet-1 ; j_jet <= i_jet ; j_jet++){
671  streamlog_out(WARNING) << std::setw(10) << i_jet << std::setw(10) << fafp[j_jet] <<
672  std::setw(10) << elementon[j_jet] << std::setw(10) << fafp_last[j_jet] <<
673  std::setw(10) << fafp_boson[4] << std::setw(10) << type[j_jet] << std::setw(10) << dijet_begining[j_jet] <<
674  std::setw(10) << dijet_end[j_jet] << std::setw(10) << boson[j_jet]<<std::endl;
675  }
676  streamlog_out(ERROR) << " Event: " << evt->getEventNumber() << ", run: " << evt->getRunNumber() << std::endl ;
677  streamlog_out(ERROR) << std::endl ;
678  } else {
679  streamlog_out(DEBUG8) << " bad match M (sum/initial) " << " " << sqrt(masssq_t)
680  << " / " << p[k[elementon[i_jet]][4]][5] << std::endl;
681  streamlog_out(DEBUG8) << " E (sum/initial) " << " " << tE[i_jet]+tE[i_jet-1]
682  << " / " << p[k[elementon[i_jet]][4]][4] << std::endl;
683  streamlog_out(DEBUG8) << " (As the energy matches well, this is probably due to the B-field) " <<std::endl;
684  streamlog_out(DEBUG8) << " Event: " << evt->getEventNumber() << ", run: " << evt->getRunNumber() << std::endl ;
685  streamlog_out(DEBUG8) << std::endl;
686  }
687  }
688  }
689  }
690  }
691  } // end jet-loop
692 
693 
694 
695  //*****************************
696 
697  // ! fill the two color-singlet blocks
698  // post-PS part
699 
700  LCCollectionVec* fafpf_vec = new LCCollectionVec(LCIO::RECONSTRUCTEDPARTICLE ) ;
701 
702  LCRelationNavigator FinalColourNeutral_Nav(LCIO::RECONSTRUCTEDPARTICLE , LCIO::RECONSTRUCTEDPARTICLE ) ;
703  LCRelationNavigator FinalElementon_Nav(LCIO::RECONSTRUCTEDPARTICLE , LCIO::MCPARTICLE ) ;
704 
705  for ( int k_dj_end=1; k_dj_end<=n_dje ; k_dj_end++ ) {
706  double E=0, M=0 , mom[3]={} ; int pdg[26]={};
707  if ( type[jets_end[1][k_dj_end]]%100 == 5 ) {
708  // beam-jet: No ancestor - just use sum of true-stable
709  E=tE[jets_end[1][k_dj_end]];
710  double psqrs=0;
711  for (int jj=0 ; jj<3 ; jj++ ) {
712  mom[jj] = tmom[jets_end[1][k_dj_end]][jj];
713  psqrs += mom[jj]*mom[jj];
714  }
715  M = sqrt(E*E-psqrs);
716 
717  } else {
718 
719  // all other cases: sum of (unique) direct decendants of the
720  // di-jet generator. For strings, clusters and ISR, there is only one
721  // (a 92, a 91 or a 22, respectively), so in these cases the sum is
722  // actualy just copying the descendant. For leptons, however, there are
723  // two (the next generation leptons), so the sum is non-trivial.
724  // (actually, sometimes also for clustes 0 see below)
725 
726  int first=nlund , last=0;
727  for ( int j_jet_end=1 ; j_jet_end<=jets_end[0][k_dj_end] ; j_jet_end++ ) {
728 
729  // decide which line should be used - this (leptons) - child (strings,gammas) - grandchild (cluster)
730  // NB. In the case of leptons, it might be that the first and last leptons are not next to
731  // eachother.
732 
733  int this_first=0, this_last=nlund;
734  if ( type[jets_end[1][k_dj_end]]%100 == 2 || type[jets_end[1][k_dj_end]]%100 == 4 || type[jets_end[1][k_dj_end]]%100 == 6 ) {
735 
736  // lepton : need to go one step less as the di-jet generator might
737  // be stable, and hence have no descendants.
738 
739  this_first=elementon[jets_end[j_jet_end][k_dj_end]];
740  this_last=elementon[jets_end[j_jet_end][k_dj_end]] ;
741 
742  } else if ( type[jets_end[1][k_dj_end]]%100 == 3 ) {
743 
744  // cluster: need to go one step further (the 91 object doesn't always
745  // have the mass of the particle it goes to). Sometimes first and
746  // last will actually become different - the cluster goes to more than one
747  // hadron.
748 
749  this_first = k[k[elementon[jets_end[j_jet_end][k_dj_end]]][4]][4];
750  this_last = k[k[elementon[jets_end[j_jet_end][k_dj_end]]][5]][5];
751 
752  } else {
753 
754  if (k[elementon[jets_end[j_jet_end][k_dj_end]]][4] != 0 ) {
755  this_first = k[elementon[jets_end[j_jet_end][k_dj_end]]][4] ;
756  this_last = k[elementon[jets_end[j_jet_end][k_dj_end]]][5];
757  } else {
758  this_first = elementon[jets_end[j_jet_end][k_dj_end]] ;
759  this_last = elementon[jets_end[j_jet_end][k_dj_end]];
760  }
761 
762  }
763 
764  if ( this_first < first ) {
765  first= this_first;
766  }
767  if ( this_last > last ) {
768  last= this_last;
769  }
770 
771  pdg[j_jet_end]=k[elementon[jets_end[j_jet_end][k_dj_end]]][2];
772  if ( pdg[0] == 0 ) {
773  if ( k[elementon[jets_end[j_jet_end][k_dj_end]]][4] != 0 ) {
774  pdg[0]=k[k[elementon[jets_end[j_jet_end][k_dj_end]]][4]][2];
775  } else {
776  pdg[0]=k[elementon[jets_end[j_jet_end][k_dj_end]]][2];
777  }
778  }
779  }
780  if ( last == first ) {
781  E= p[first][4];
782  for ( int ll=1 ; ll<=3 ; ll++ ) {
783  mom[ll-1] = p[first][ll];
784  }
785  } else {
786  // need double loop for non-consecutive case ( if this could only be F95 !)
787  for ( int j_jet_end=1 ; j_jet_end<=jets_end[0][k_dj_end] ; j_jet_end++ ) {
788  for ( int k_py=first ; k_py<=last ; k_py++ ) {
789 
790  if ( abs(jet[k_py]) == jets_end[j_jet_end][k_dj_end] ) {
791  E+= p[k_py][4];
792  for ( int ll=1 ; ll<=3 ; ll++ ) {
793  mom[ll-1] += p[k_py][ll];
794  }
795  if ( type[jets_end[j_jet_end][k_dj_end]]%100 != 3 ) {break;}
796  }
797 
798  }
799  }
800  }
801  double psqrs=0 ;
802  for (int ii=0 ; ii<3 ; ii++ ) {
803  psqrs += mom[ii]*mom[ii];
804  }
805  M = sqrt(E*E-psqrs);
806  }
807 
808  ReconstructedParticleImpl* fafpf = new ReconstructedParticleImpl;
809  fafpf->setCharge(0);
810  fafpf->setEnergy(E);
811  fafpf->setMass(M);
812  fafpf->setMomentum(mom);
813  for(int j_jet_end=1 ; j_jet_end<=jets_end[0][k_dj_end] ; j_jet_end++) {
814  fafpf->addParticle(dynamic_cast<ReconstructedParticle*>(jet_vec->getElementAt(jets_end[j_jet_end][k_dj_end]-1)) );
815  }
816  ParticleIDImpl* pid[26];
817  pid[0] = new ParticleIDImpl;
818  pid[0]->setPDG(pdg[0]);
819  pid[0]->setType(type[jets_end[1][k_dj_end]]%100); // maybe flag from boson? could be 0,1, or 2 jets in the fafp that's from boson ..
820  fafpf->addParticleID(pid[0]);
821  for ( int j_jet_end=1 ; j_jet_end<=jets_end[0][k_dj_end] ; j_jet_end++ ) {
822  pid[j_jet_end]= new ParticleIDImpl;
823  pid[j_jet_end]->setPDG(pdg[j_jet_end]);
824  pid[j_jet_end]->setType(type[jets_end[j_jet_end][k_dj_end]]);
825  fafpf->addParticleID(pid[j_jet_end]);
826  if ( elementon[jets_end[j_jet_end][k_dj_end]] > 0 ) {
827  FinalElementon_Nav.addRelation( jet_vec->getElementAt(jets_end[j_jet_end][k_dj_end]-1) ,
828  mcp_pyjets[elementon[jets_end[j_jet_end][k_dj_end]]], 1.0 );
829  }
830  FinalColourNeutral_Nav.addRelation( jet_vec->getElementAt(jets_end[j_jet_end][k_dj_end]-1) , fafpf);
831  }
832  fafpf_vec->addElement(fafpf);
833 
834  } // di-jet(end) loop
835 
836 
837  // pre-PS part
838 
839  LCCollectionVec* fafpi_vec = new LCCollectionVec(LCIO::RECONSTRUCTEDPARTICLE ) ;
840 
841  LCRelationNavigator InitialElementon_Nav(LCIO::RECONSTRUCTEDPARTICLE , LCIO::MCPARTICLE ) ;
842  LCRelationNavigator InitialColourNeutral_Nav(LCIO::RECONSTRUCTEDPARTICLE , LCIO::RECONSTRUCTEDPARTICLE ) ;
843 
844  for (int k_dj_begin=1 ; k_dj_begin<=n_djb ; k_dj_begin++ ) {
845  double E=0, M=0 , mom[3]={0} ;
846  int pdg[26]= {0};
847 
848  int bosonid=23 ;
849 
850  // E and P is sum of (unique) direct decendants of the
851  // initial ffbar. Normally, there is only one
852  // (a 94), so the sum is actualy just copying the descendant.
853  // For leptons, however, there is no 94 - they go directly to
854  // the next generation leptons), so the sum is non-trivial.
855  // Sometimes this also happens for quraks, ie. they go directly to
856  // the genertor quarks of the string.
857 
858  int first=nlund , last=0;
859  if ( type[jets_begin[1][k_dj_begin]] == 4 ) { // i.e. ISR
860 
861  pdg[1]=k[elementon[jets_begin[1][k_dj_begin]]][2];
862  first= fafp[jets_begin[1][k_dj_begin]] ;
863  last= fafp[jets_begin[1][k_dj_begin]] ;
864  bosonid=22;
865 
866  } else {
867 
868  for ( int j_jet_begin=1 ; j_jet_begin<=jets_begin[0][k_dj_begin] ; j_jet_begin++ ) {
869 
870 
871  int gen_part = fafp [jets_begin[j_jet_begin][k_dj_begin]];
872  bool take_genpart = false;
873  if ( k[gen_part][4] == 0 ) {
874  take_genpart = true;
875  } else if ( k [k [gen_part] [4]][2] != k[gen_part][2] || k [k [gen_part] [4]][1] == 0 ) {
876  take_genpart = true;
877  }
878  if ( take_genpart ) {
879  if ( gen_part < first ) {
880  first=gen_part;
881  }
882  if ( gen_part > last ) {
883  last=gen_part;
884  }
885  } else {
886  if ( k[gen_part][4] < first ) {
887  first= k[gen_part][4] ;
888  }
889  if ( k[gen_part][5] > last ) {
890  last= k[gen_part][5] ;
891  }
892 
893  }
894 
895  pdg[j_jet_begin]=k[elementon[jets_begin[j_jet_begin][k_dj_begin]]][2];
896  int last_backtrack = fafp[jets_begin[j_jet_begin][k_dj_begin]] ;
897  int this_pdg = k[last_backtrack][2] ;
898 
899  while ( k[last_backtrack][2] == this_pdg ) {
900  last_backtrack = k[last_backtrack][3];
901  if ( last_backtrack == 0 ) break;
902  }
903 
904  if ( last_backtrack > 0) {
905  if ( abs(k[last_backtrack][2]) > 21 && abs(k[last_backtrack][2]) <= 39 ) {
906  // ie. any IVB except the gluon
907  bosonid = abs(k[last_backtrack][2]) ;
908  streamlog_out(DEBUG3) << " elementon " << fafp[jets_begin[j_jet_begin][k_dj_begin]]
909  << " has explicit mother with PDG = " << bosonid << std::endl;
910  }
911  }
912 
913  if ( pdg[0] == 0 ) {
914 
915  pdg[0]=abs(k[fafp[jets_begin[j_jet_begin][k_dj_begin]]][2]);
916 
917  }
918  else if ( bosonid==23 && abs(k[fafp[jets_begin[j_jet_begin][k_dj_begin]]][2]) != pdg[0] ) {
919 
920  bosonid=24;
921 
922  }
923  }
924  }
925  if ( k[last][2] == 92 ) { first = last ; } // IS THIS RIGHT ? If there is a gluon radiation ?!
926  // should it rather be the other way last="first"
927  pdg[0]=bosonid;
928 
929  if ( last == first ) {
930  E= p[first][4];
931  for ( int ll=1 ; ll<=3 ; ll++ ) {
932  mom[ll-1] = p[first][ll];
933  }
934  } else {
935 
936  // if several, sum up but make sure that there are no intruders between first and last.
937  for ( int j_jets_begin=1 ; j_jets_begin<=jets_begin[0][k_dj_begin] ; j_jets_begin++ ) {
938  for ( int k_py=first ; k_py<=last ; k_py++ ) {
939  if ( abs(jet[k_py]) == jets_begin[j_jets_begin][k_dj_begin] && k[k_py][3] < first ) {
940  E+= p[k_py][4];
941  for ( int ll=1 ; ll<=3 ; ll++ ) {
942  mom[ll-1] += p[k_py][ll];
943  }
944  }
945 
946  }
947  }
948  }
949 
950  double psqrs=0;
951  for ( int ii=1 ; ii<=3 ; ii++ ) {
952  psqrs += mom[ii-1]*mom[ii-1];
953  }
954 
955  M = sqrt(E*E-psqrs);
956 
957 
958  ReconstructedParticleImpl* fafpi = new ReconstructedParticleImpl;
959  fafpi->setCharge(0);
960  fafpi->setEnergy(E);
961  fafpi->setMass(M);
962  fafpi->setMomentum(mom);
963  ParticleIDImpl* pid[26] ;
964  pid[0]= new ParticleIDImpl;
965  pid[0]->setPDG(pdg[0]);
966  pid[0]->setType(type[jets_begin[1][k_dj_begin]]%100);
967  fafpi->addParticleID(pid[0]);
968  for(int j_jet_begin=1 ; j_jet_begin<=jets_begin[0][k_dj_begin] ; j_jet_begin++) {
969  pid[j_jet_begin]= new ParticleIDImpl;
970  pid[j_jet_begin]->setPDG(pdg[j_jet_begin]);
971  pid[j_jet_begin]->setType(type[jets_begin[j_jet_begin][k_dj_begin]]);
972  fafpi->addParticleID(pid[j_jet_begin]);
973  fafpi->addParticle(dynamic_cast<ReconstructedParticle*>(jet_vec->getElementAt(jets_begin[j_jet_begin][k_dj_begin]-1)));
974  }
975  fafpi_vec->addElement(fafpi);
976 
977  for ( int j_jet_begin=1 ; j_jet_begin<=jets_begin[0][k_dj_begin] ; j_jet_begin++ ) {
978  InitialElementon_Nav.addRelation( jet_vec->getElementAt(jets_begin[j_jet_begin][k_dj_begin]-1),
979  mcp_pyjets[fafp[jets_begin[j_jet_begin][k_dj_begin]]], 1.0 );
980  InitialColourNeutral_Nav.addRelation( jet_vec->getElementAt(jets_begin[j_jet_begin][k_dj_begin]-1) , fafpi);
981  }
982 
983  } // di-jet (begining) loop
984 
985  //============================
986 
987 
988  //*****************************
989  evt->addCollection( jet_vec, _trueJetCollectionName);
990  evt->addCollection( fafpf_vec, _finalColourNeutralCollectionName);
991  evt->addCollection( fafpi_vec, _initialColourNeutralCollectionName);
992 
993  LCCollection* tjfpcol = 0;
994  LCCollection* tjlpcol = 0;
995  LCCollection* tjfccol = 0;
996  LCCollection* tjlccol = 0;
997  tjrcol = truejet_pfo_Nav.createLCCollection() ;
998  tjtcol = truejet_truepart_Nav.createLCCollection() ;
999  tjfpcol = InitialElementon_Nav.createLCCollection() ;
1000  tjlpcol = FinalElementon_Nav.createLCCollection() ;
1001  tjfccol = InitialColourNeutral_Nav.createLCCollection() ;
1002  tjlccol = FinalColourNeutral_Nav.createLCCollection() ;
1003  evt->addCollection( tjrcol , _trueJetPFOLink);
1004  evt->addCollection( tjtcol , _trueJetMCParticleLink);
1005  evt->addCollection( tjfpcol , _initialElementonLink);
1006  evt->addCollection( tjlpcol , _finalElementonLink);
1007  evt->addCollection( tjfccol , _initialColourNeutralLink);
1008  evt->addCollection( tjlccol , _finalColourNeutralLink);
1009  //============================
1010 
1011 
1012 
1013  streamlog_out(DEBUG4) << " HEPEVT relation table up with jet assignment, extended status codes, and PFOs w. energy (if any) " <<std::endl;
1014  streamlog_out(DEBUG4) << " line status pdg first last first last colour anti- px py pz "
1015  <<" E M jet companion type PFO/Energy" << std::endl;
1016  streamlog_out(DEBUG4) << " init ext'd parent daughter colour "
1017  <<" jet" << std::endl;
1018 
1019  // for ( int i_py=1 ; i_py <= nlund ; i_py++){
1020  for ( int i_py=1 ; i_py <= std::min(nlund,200) ; i_py++){
1021 
1022  int istat=0;
1023  if ( jet[i_py] != 0 ) {
1024  static FloatVec www;
1025  www = truejet_truepart_Nav.getRelatedFromWeights(mcp_pyjets[i_py]);
1026  istat=int(www[0]);
1027  }
1028 
1029  if ( k[i_py][1] > 0 && k[i_py][1] < 30 ) {
1030 
1031  streamlog_out(DEBUG4) << std::setw(4) << i_py << std::setw(7) << k[i_py][1] << std::setw(5) << istat <<
1032  std::setw(7) << k[i_py][2] << std::setw(7) << k[i_py][3] << std::setw(5) << k[i_py][6] <<
1033  std::setw(7) << k[i_py][4] << std::setw(5) << k[i_py][5] <<
1034  std::setw(7) << k[i_py][7] << std::setw(7) << k[i_py][8] <<
1035  std::setw(12) << p[i_py][1] <<
1036  std::setw(12) << p[i_py][2] << std::setw(12) << p[i_py][3] <<
1037  std::setw(12) << p[i_py][4] << std::setw(12) << p[i_py][5] <<
1038 
1039  std::setw(7) << jet[i_py] <<
1040  std::setw(7) << companion[abs(jet[i_py])] << std::setw(7) << type[abs(jet[i_py])];
1041  } else {
1042  streamlog_out(DEBUG3) << std::setw(4) << i_py << std::setw(7) << k[i_py][1] << std::setw(5) << istat <<
1043  std::setw(5) << k[i_py][2] << std::setw(7) << k[i_py][3] << std::setw(5) << k[i_py][6] <<
1044  std::setw(7) << k[i_py][4] << std::setw(5) << k[i_py][5] <<
1045  std::setw(7) << k[i_py][7] << std::setw(7) << k[i_py][8] <<
1046  std::setw(12) << p[i_py][1] <<
1047  std::setw(12) << p[i_py][2] << std::setw(12) << p[i_py][3] <<
1048  std::setw(12) << p[i_py][4] << std::setw(12) << p[i_py][5] <<
1049 
1050  std::setw(7) << jet[i_py] <<
1051  std::setw(7) << companion[abs(jet[i_py])] << std::setw(7) << type[abs(jet[i_py])];
1052  }
1053 
1054  if ( reltrue ) {
1055  LCObjectVec recovec = reltrue->getRelatedFromObjects( mcp_pyjets[i_py]);
1056  if ( recovec.size() > 0 ) { // if reconstructed
1057  for ( unsigned i_reco=0 ; i_reco<recovec.size() ; i_reco++ ) { // reco-of-this-true loop
1058  if ( k[i_py][1] > 0 && k[i_py][1] < 30 ) {
1059 
1060  ReconstructedParticle* reco_part = dynamic_cast<ReconstructedParticle*>(recovec[i_reco]);
1061  streamlog_out(DEBUG4) << " [" << std::setw(8) << reco_part->getEnergy() <<"] ";
1062  } else {
1063  ReconstructedParticle* reco_part = dynamic_cast<ReconstructedParticle*>(recovec[i_reco]);
1064  streamlog_out(DEBUG3) << " [" << std::setw(10) << reco_part <<" /"<< std::setw(8) << reco_part->getEnergy() <<"] ";
1065  }
1066  }
1067  } else if ( istat == 1 ) {
1068  if ( k[i_py][1] > 0 && k[i_py][1] < 30 ) {
1069  streamlog_out(DEBUG4) << " [" << std::setw(10) << "N.A "<<" /"<< std::setw(8) << "N.A "<<"] ";
1070  } else {
1071  streamlog_out(DEBUG3) << " [" << std::setw(10) << "N.A "<<" /"<< std::setw(8) << "N.A "<<"] ";
1072  }
1073  }
1074  }
1075  if ( k[i_py][1] > 0 && k[i_py][1] < 30 ) {
1076  streamlog_out(DEBUG4) << std::endl;
1077  } else {
1078  streamlog_out(DEBUG3) << std::endl;
1079  }
1080  }
1081  streamlog_out(DEBUG6) << std::endl;
1082 
1083  streamlog_out(DEBUG6) << " jet fafp-beg qrk/lept fafp-end fafp-boson type "
1084  <<"dijet-beg dijet-end boson" << std::endl;
1085  for ( int i_jet=1 ; i_jet <= njet ; i_jet++){
1086  streamlog_out(DEBUG6) << std::setw(10) << i_jet << std::setw(10) << fafp[i_jet] <<
1087  std::setw(10) << elementon[i_jet] << std::setw(10) << fafp_last[i_jet] <<
1088  std::setw(10) << fafp_boson[i_jet] << std::setw(10) << type[i_jet] <<
1089  std::setw(10) << dijet_begining[i_jet] <<
1090  std::setw(10) << dijet_end[i_jet] << std::setw(10) << boson[i_jet]<<std::endl;
1091  }
1092  streamlog_out(DEBUG6) << std::endl;
1093  streamlog_out(DEBUG6) << " Colour-neutral objects at the beginning of the parton shower " << std::endl;
1094  streamlog_out(DEBUG6) << " number px py pz E M "
1095  <<" type PDGs jets " << std::endl;
1096  int glurad=0;
1097  for ( unsigned i_cn_b=0 ; i_cn_b<fafpi_vec->size() ; i_cn_b++ ) {
1098 
1099  ReconstructedParticle* fafpi = dynamic_cast<ReconstructedParticle*>(fafpi_vec->at(i_cn_b));
1100  LCObjectVec jts = InitialColourNeutral_Nav.getRelatedFromObjects(fafpi);
1101 
1102  const double* mom=fafpi->getMomentum();
1103  streamlog_out(DEBUG6) << std::setw(7) << i_cn_b+1 << std::setw(12) << mom[0] <<
1104  std::setw(12) << mom[1] << std::setw(12) << mom[2] <<
1105  std::setw(12) <<fafpi->getEnergy() << std::setw(12) << fafpi->getMass() <<
1106  std::setw(7) << fafpi->getParticleIDs()[0]->getType() << " " ;
1107  for ( unsigned i_pid=0 ; i_pid < fafpi->getParticleIDs().size() ; i_pid++ ) {
1108  if (fafpi->getParticleIDs()[i_pid]->getType()>=100 ) glurad=1;
1109  if (i_pid < fafpi->getParticleIDs().size()-1 ) {
1110  streamlog_out(DEBUG6) << std::setw(3)<< fafpi->getParticleIDs()[i_pid]->getPDG() << "," ;
1111  } else {
1112  streamlog_out(DEBUG6) << std::setw(3)<< fafpi->getParticleIDs()[i_pid]->getPDG() << " " ;
1113  }
1114  }
1115  for ( unsigned j_jet_begin=0 ; j_jet_begin<jts.size() ; j_jet_begin++ ) {
1116  int i_jet=jts[j_jet_begin]->ext<TJindex>();
1117  if ( j_jet_begin < jts.size()-1 ) {
1118  streamlog_out(DEBUG6) << std::setw(3) << i_jet << "," ;
1119  } else {
1120  streamlog_out(DEBUG6) << std::setw(3) << i_jet ;
1121  }
1122  }
1123  streamlog_out(DEBUG6) << std::endl;
1124  }
1125  if ( glurad == 1 ) {
1126  streamlog_out(DEBUG6) << std::endl;
1127  streamlog_out(DEBUG6) << " gluon radiation in event : "<< std::endl;
1128  for ( unsigned j_cn_b=0 ; j_cn_b<fafpi_vec->size() ; j_cn_b++ ) {
1129 
1130  ReconstructedParticle* fafpi = dynamic_cast<ReconstructedParticle*>(fafpi_vec->at(j_cn_b));
1131  LCObjectVec jts = InitialColourNeutral_Nav.getRelatedFromObjects(fafpi);
1132 
1133  for ( unsigned j_pid_cn_b=1 ; j_pid_cn_b < fafpi->getParticleIDs().size() ; j_pid_cn_b++ ) {
1134  if (fafpi->getParticleIDs()[j_pid_cn_b]->getType()>=100 ) {
1135  int i_jet=jts[j_pid_cn_b-1]->ext<TJindex>();
1136  streamlog_out(DEBUG6) << std::setw(18)<< i_jet << " is from "
1137  << fafpi->getParticleIDs()[j_pid_cn_b]->getType()/100 << std::endl;
1138  }
1139  }
1140  }
1141  }
1142 
1143  streamlog_out(DEBUG6) << std::endl;
1144  streamlog_out(DEBUG6) << " Colour-neutral objects at the end of the parton shower " << std::endl;
1145  streamlog_out(DEBUG6) << " number px py pz E M "
1146  <<" type PDGs jets " << std::endl;
1147  for ( unsigned i_cn_e=0 ; i_cn_e<fafpf_vec->size() ; i_cn_e++ ) {
1148 
1149  ReconstructedParticle* fafpf = dynamic_cast<ReconstructedParticle*>(fafpf_vec->at(i_cn_e));
1150  LCObjectVec jts = FinalColourNeutral_Nav.getRelatedFromObjects(fafpf);
1151 
1152  const double* mom=fafpf->getMomentum();
1153  streamlog_out(DEBUG6) << std::setw(7) << i_cn_e+1 << std::setw(12) << mom[0] <<
1154  std::setw(12) << mom[1] << std::setw(12) << mom[2] <<
1155  std::setw(12) <<fafpf->getEnergy() << std::setw(12) << fafpf->getMass() <<
1156  std::setw(7) << fafpf->getParticleIDs()[0]->getType() << " ";
1157  if ( fafpf->getParticleIDs().size() < 3 ) streamlog_out(DEBUG6) << std::setw(3) <<" ";
1158  for ( unsigned i_pid=0 ; i_pid < fafpf->getParticleIDs().size() ; i_pid++ ) {
1159  if (i_pid < fafpf->getParticleIDs().size()-1 ) {
1160  streamlog_out(DEBUG6) << std::setw(3)<< fafpf->getParticleIDs()[i_pid]->getPDG() << "," ;
1161  } else {
1162  streamlog_out(DEBUG6) << std::setw(3)<< fafpf->getParticleIDs()[i_pid]->getPDG() << " " ;
1163  }
1164  }
1165  if ( jts.size() < 2 ) streamlog_out(DEBUG6) << std::setw(3) <<" ";
1166  for ( unsigned j_jet_end=0 ; j_jet_end<jts.size() ; j_jet_end++ ) {
1167  int i_jet=jts[j_jet_end]->ext<TJindex>();
1168  if ( j_jet_end < jts.size()-1 ) {
1169  streamlog_out(DEBUG6) << std::setw(3) << i_jet << "," ;
1170  } else {
1171  streamlog_out(DEBUG6) << std::setw(3) << i_jet ;
1172  }
1173  }
1174  streamlog_out(DEBUG6) << std::endl;
1175  }
1176  streamlog_out(DEBUG6) << std::endl;
1177 
1178 
1179  } // endif( mcpcol != NULL )
1180 
1181 
1182  //*****************************
1183  // (this guy is not deleted by lcio, so ...
1184  if ( reltrue ) delete reltrue ;
1185  // ... to avoid a memory leak!)
1186  _nEvt ++ ;
1187  //============================
1188 }
1189 
1190 
1191 
1192 void TrueJet::getPyjets(LCCollection* mcpcol )
1193 {
1194  const double* ptemp;
1195  _higgs_to_glue_glue = false;
1196  bool _generator_input_format = false;
1197  int higgs_line = 0 ;
1198  int k_py=0;
1199  int nMCP = mcpcol->getNumberOfElements() ;
1200  _top_event = false ;
1201  first_line=1;
1202 
1203  if ( nMCP > 4000 ) {
1204  streamlog_out(WARNING) << " More than 4000 MCParticles in event " << evt->getEventNumber()
1205  << ", run " << evt->getRunNumber() << std::endl ;
1206  }
1207 
1208  for(int j_mcp=0; j_mcp< std::min(nMCP,4000) ; j_mcp++){
1209 
1210  MCParticle* mcp = dynamic_cast<MCParticle*>( mcpcol->getElementAt( j_mcp ) ) ;
1211  // (assume this is not needed, ie. that the k-vector can be set up correctly
1212  // also for created-in-simulation particles. The issue is if decay-products
1213  // of a given particle will be consecutive).
1214  // if (mcp->getGeneratorStatus() == 1 || mcp->getGeneratorStatus() == 2 ) {
1215 
1216  if ( j_mcp == 0 && mcp->getSimulatorStatus()==0 ) { _generator_input_format = true ; }
1217  k_py++;
1218  mcp_pyjets[k_py]=mcp;
1219  mcp->ext<MCPpyjet>()=k_py;
1220  ptemp=mcp->getMomentum();
1221  p[k_py][1]=ptemp[0];
1222  p[k_py][2]=ptemp[1];
1223  p[k_py][3]=ptemp[2];
1224  p[k_py][4]=mcp->getEnergy();
1225  p[k_py][5]=mcp->getMass();
1226  k[k_py][1]=mcp->getGeneratorStatus();
1227  if ( k[k_py][1] == 2 ) { k[k_py][1]=11 ;}
1228  if ( k[k_py][1] == 3 ) { k[k_py][1]=21 ;}
1229  if ( k[k_py][1] == 4 ) { k[k_py][1]=22 ;}
1230  if (mcp->isOverlay() ) { k[k_py][1]=k[k_py][1]+30 ; }
1231  k[k_py][2]=mcp->getPDG();
1232  if ( abs(k[k_py][2]) == 6 ) { _top_event = true ; }
1233  if ( k[k_py][2] == 25 ) {
1234  higgs_line = k_py ;
1235  }
1236 
1237  // }
1238  }
1239  nlund=k_py;
1240 
1241  for (int j_py=1; j_py <= nlund ; j_py++ ) {
1242  k[j_py][3] = 0;
1243  k[j_py][4] = 0;
1244  k[j_py][5] = 0;
1245  k[j_py][7] = 0;
1246  k[j_py][8] = 0;
1247  }
1248  streamlog_out(DEBUG1) << " before Motherless check: first_line = " << first_line << std::endl ;
1249 
1250  for(int j_mcp=0; j_mcp< nMCP ; j_mcp++){
1251 
1252  MCParticle* mcp = dynamic_cast<MCParticle*>( mcpcol->getElementAt( j_mcp ) ) ;
1253  int i_py=0;
1254  if ( mcp->ext<MCPpyjet>() != 0 ) {
1255 
1256  i_py= mcp->ext<MCPpyjet>();
1257  if ( mcp->getParents().size() != 0 ) {
1258  k[i_py][3]=mcp->getParents()[0]->ext<MCPpyjet>();
1259  if ( k[i_py][3] == higgs_line ) {
1260  if ( k[i_py][2] == 21 ) {
1261  _higgs_to_glue_glue = true;
1262  }
1263  }
1264  } else {
1265  if ( i_py > first_line+1 && ! mcp->isOverlay() && k[i_py][3] == 0 && abs(mcp->getPDG()) != 6 && abs(mcp->getPDG()) != 24 ) {
1266  streamlog_out(MESSAGE4) << " Motherless generator particle " << i_py << ". pdg: " << mcp->getPDG() << std::endl ;
1267  streamlog_out(MESSAGE4) << " Event: " << evt->getEventNumber() << ", run: " << evt->getRunNumber() << std::endl ;
1268  k[i_py][3]=-1;
1269  }
1270 
1271  }
1272 
1273  k[i_py][7]=mcp->getColorFlow()[0];
1274  k[i_py][8]=mcp->getColorFlow()[1];
1275 
1276  if ( mcp->getParents().size() > 1 ) {
1277  k[i_py][6]=mcp->getParents()[mcp->getParents().size()-1]->ext<MCPpyjet>();
1278  } else {
1279  k[i_py][6]=0;
1280  }
1281 
1282  if ( mcp->getDaughters().size() != 0 ) {
1283 
1284  // check daughter-status. If any of the daughters are created in simulation,
1285  // or has simulator-status 0 (indicating that it was ignored by Geant),
1286  // set k[][1] to stable of the mother, and change it to 0 for direct descendants.
1287  // if the current particle has k[][1] == 0, propagate to its descendants, that way
1288  // all generations of a k[][1] == 0 ancestor will have k[][1] == 0 once
1289  // the loop is done.
1290 
1291  // (A further case is if all daughters have genstat 1, simstat !=0, but the sum of their energies is
1292  // different from that of the mother (indicating that Geant dis something un-documented with them)
1293  // This case can't be fixed, yet as the check should only be applied to lines after the parton-shower,
1294  // which we can't determine yet. Instead it is taken care of in the gouping method.)
1295 
1296  if ( ! _generator_input_format ) { // none of the below can happen if we are looking at the generator input MCparticle collection
1297  if ( k[i_py][1]!=0 ) {
1298  int has_genstat0_daughter=0;
1299  for ( unsigned i_dau=0 ; i_dau < mcp->getDaughters().size() ; i_dau++ ) {
1300  int j_dau=mcp->getDaughters()[i_dau]->ext<MCPpyjet>();
1301  if ( mcp->getDaughters()[i_dau]->getGeneratorStatus()==0 || (mcp->getDaughters()[i_dau]->getSimulatorStatus()==0 && k[j_dau][1]==1 )) {
1302  has_genstat0_daughter=1;
1303  }
1304  }
1305  if ( has_genstat0_daughter==1 ) {
1306 
1307  if ( k[i_py][1] < 30 ) {
1308  k[i_py][1]=1;
1309  }
1310 
1311  for ( unsigned i_dau=0 ; i_dau < mcp->getDaughters().size() ; i_dau++ ) {
1312  int j_dau=mcp->getDaughters()[i_dau]->ext<MCPpyjet>();
1313  if ( k[j_dau][1] < 30 ) {
1314  k[j_dau][1]=0;
1315  }
1316  }
1317 
1318  }
1319 
1320  } else {
1321 
1322  for ( unsigned i_dau=0 ; i_dau < mcp->getDaughters().size() ; i_dau++ ) {
1323  int j_dau=mcp->getDaughters()[i_dau]->ext<MCPpyjet>();
1324  if ( k[j_dau][1] < 30 ) {
1325  k[j_dau][1]=0;
1326  }
1327  }
1328  }
1329  }
1330 
1331  // fill in k[][4] and [5]. Note that it is possible that one gets
1332  // k[][1]=1 but k[][4] and [5] != 0. This means that a generator
1333  // particle end in a simulation particle, which is OK for the
1334  // logic elsewhere !
1335 
1336  for ( unsigned i_dau=0 ; i_dau < mcp->getDaughters().size() ; i_dau++ ) {
1337  int j_dau=mcp->getDaughters()[i_dau]->ext<MCPpyjet>();
1338  if ( j_dau < k[i_py][4] || k[i_py][4] == 0 ) {
1339  k[i_py][4]=j_dau;
1340  }
1341  if ( j_dau > k[i_py][5] ) {
1342  k[i_py][5]=j_dau;
1343  }
1344  if ( abs(k[i_py][2]) == 6 ) { // fix for top quark - possibly not needed with whiz2 ?
1345  if (k[k[i_py][4]][2] != 94 ) {
1346  k[i_py][5]= k[i_py][4]+1;
1347  k[k[i_py][5]][3] = i_py;
1348  }
1349  }
1350 
1351  }
1352  }
1353  } // endif ( mcp->ext<MCPpyjet>() != 0 )
1354  } // end MCP loop
1355 
1356 
1357  // At this point basically 99% of the job is done. However, there is still the CMShower lines (pdg=94).
1358  // These have several mothers and several daughters. The 4-momentum is conserved in total, but not for
1359  // the individual components one-by-one. Here we will short-circuit them, i.e. the history codes will
1360  // connect mother and daughters directly. This is no problem if all pdg:s are different in the group
1361  // of particles going in and out of the 94. The problem arrises when this is not the case (e.g. e+e-e+e-).
1362  // Here we need to look at the momenta of the particles and match which dauther is closest to which
1363  // mother. Remember, they will *not* be identical, so there is a bit of guessing.
1364  // Note that the mothers of a 94 are not necessarily adjecent. They *are* in the range between k[line94][3]
1365  // and k[line94][6], but there might be "intruders" between these lines.
1366 
1367 
1368  streamlog_out(DEBUG4) << " HEPEVT relation table before 94-fixing" <<std::endl;
1369  streamlog_out(DEBUG4) << " line status pdg parent first last second anti " << std::endl;
1370  streamlog_out(DEBUG4) << " daughter daugther parent colour colour " << std::endl;
1371  for ( int j_py=1 ; j_py<=nlund ; j_py++ ) {
1372  if ( k[j_py][1] != 0 && k[j_py][1] < 30 ) {
1373  streamlog_out(DEBUG4) << std::setw(10) << j_py << std::setw(10) << k[j_py][1] <<
1374  std::setw(10) << k[j_py][2] << std::setw(10) << k[j_py][3] <<
1375  std::setw(10) << k[j_py][4] << std::setw(10) << k[j_py][5] <<
1376  std::setw(10) << k[j_py][6] << std::setw(10)<< k[j_py][7] << std::setw(10) << k[j_py][8] << std::endl;
1377  } else {
1378  streamlog_out(DEBUG3) << std::setw(10) << j_py << std::setw(10) << k[j_py][1] <<
1379  std::setw(10) << k[j_py][2] << std::setw(10) << k[j_py][3] <<
1380  std::setw(10) << k[j_py][4] << std::setw(10) << k[j_py][5] <<
1381  std::setw(10) << k[j_py][6] << std::setw(10)<< k[j_py][7] << std::setw(10) << k[j_py][8] << std::endl;
1382  }
1383  }
1384 
1385  // 94 fix.
1386 
1387  for ( int i_py=1 ; i_py<=nlund ; i_py++ ) {
1388 
1389  // if the daughters of the parents to a 94 is not (only) the 94,
1390  // check if between the 94 line and the second daughter of the
1391  // parent of the 94 actually has the particle *going in* to the 94 (not the
1392  // 94) as parent. If so, change the dautghter lines of the incomming particle
1393  // to these lines (and ignore the 94). If not, move those that are not daughters
1394  // to the 94 instead, and set also the second daughter of the parent to the
1395  // 94.
1396 
1397  if ( k[i_py][2] == 94 && k[i_py][1] < 30 ) {
1398  // A 94: treat that:
1399  // case 1: 2->2 , clear-cut connect in to out dirctly
1400  // case 2: 1->1+gamma : FSR, chnage 94 to 95, else unchanged.
1401  // case 3: 4->4, unique flavours, incomming has only 94 as daughter, outgoind only 94 as mother: : as case 1
1402  // case 4: as case 3, but flavours not unique: Check directions etc, then connect directly
1403  // case 5: 4->4, some weidness in parent-daughter relations: to study
1404  // case 6: n->m, n != m: to study
1405  // case 7: odd->m : Does that ever happen, except for the FSR case?
1406  // case 8: Anything else : Does that ever happen
1407  // case X: any of the above, but with intruders between first and last mother and/or daughter: to study
1408  //
1409 
1410  int mothers[10] ={0};
1411  int daughters[10][10] ={0};
1412 
1413  int line94 = i_py ;
1414  if ( _whiz1 ) { stdhep_reader_bug_workaround(line94) ; }
1415 
1416  int first_94_mother= k[line94][3];
1417  int last_94_mother= k[line94][6];
1418  int first_94_daughter = k[line94][4];
1419  int last_94_daughter = k[line94][5];
1420 
1421  int i_mo=0;
1422  for ( int mother = first_94_mother ; mother <= last_94_mother; mother++ ) {
1423  if ( (abs(k[mother][2]) == 24 && abs(k[k[mother][3]][2]) !=6 ) || k[mother][2]==23 || k[mother][2]==25 ) continue;
1424  if ( abs(k[mother][4]) != line94 ) continue; // intruder, go to next line
1425 
1426  i_mo++;
1427  mothers[0]=i_mo;
1428  mothers[i_mo]=mother;
1429  int i_dau=0;
1430  daughters[i_mo][0]=0;
1431  int pdg_mother=k[mother][2];
1432  float p_mother[6];
1433  for ( int ip=1; ip <=5 ; ip++ ) {
1434  p_mother[ip]=p[mother][ip];
1435  }
1436 
1437  if ( k[mother][4] == line94 && k[mother][5] == line94 ) {
1438 
1439  // normal case
1440 
1441  if ( _whiz1 ) { // Maybe useful not only for whiz1 ? If so, also should check k[][7:8] == 0
1442  int n_94_mother = 0 ;
1443  for ( int j_py = first_94_mother ; j_py <= last_94_mother ; j_py++ ) {
1444  if ( k[j_py][4] == line94 ) { n_94_mother++ ; }
1445  }
1446  if ( n_94_mother == 2 ) {
1447  if ( k[mother][2] > 0 ) {
1448  k[mother][7] = i_py ; k[mother][8] = 0 ;
1449  k[k[mother][3]][7] = i_py ; k[k[mother][3]][8] = 0 ;
1450  } else {
1451  k[mother][7] = 0 ; k[mother][8] = i_py ;
1452  k[k[mother][3]][7] = 0 ; k[k[mother][3]][8] = i_py ;
1453  }
1454  int oma=mother;
1455  while ( oma > 3 ) {
1456  k[oma][7] = k[mother][7] ; k[oma][8]=k[mother][8] ;
1457  oma=k[oma][3];
1458  }
1459  }
1460  } // endif ( _whiz1 )
1461 
1462  // find the daughter that best matches the current mother
1463 
1464  int best_match = 0;
1465  double min_dist = 1.0e20;
1466  for ( int daughter=first_94_daughter ; daughter <= last_94_daughter ; daughter++ ) {
1467  if ( k[daughter][2] == pdg_mother ) {
1468  double dist=((p_mother[1]-p[daughter][1])*(p_mother[1]-p[daughter][1]) +
1469  (p_mother[2]-p[daughter][2])*(p_mother[2]-p[daughter][2]) +
1470  (p_mother[3]-p[daughter][3])*(p_mother[3]-p[daughter][3]));
1471  if ( dist < min_dist ) {
1472  best_match=daughter;
1473  min_dist=dist;
1474  }
1475  }
1476  }
1477 
1478  if ( best_match == 0 ) {
1479  streamlog_out(ERROR) << " No best match ... " << std::endl;
1480  } else {
1481 
1482  i_dau++;
1483  daughters[i_mo][0]=i_dau;
1484  daughters[i_mo][i_dau]=best_match;
1485 
1486  }
1487 
1488  } else { // mother has several daughters, find the daughter that has the same pdg as the mother (actually an intruder)
1489 
1490  bool found_daughter = false ;
1491  for ( int daughter=k[mother][4] ; daughter <= k[mother][5] ; daughter++ ) {
1492  if ( k[daughter][3] == mother && daughter != line94 ) {
1493  found_daughter = true;
1494  i_dau++;
1495  daughters[i_mo][0]=i_dau;
1496  daughters[i_mo][i_dau]=daughter;
1497  }
1498  }
1499 
1500  if ( ! found_daughter ) {
1501  streamlog_out(ERROR) << " Only daughter of the mother of a 94 is neither the 94, "
1502  <<"nor a particle with same pdg as the mother " << std::endl;
1503  }
1504  }
1505 
1506  if (daughters[i_mo][0] == 0 ) {
1507  streamlog_out(ERROR) << " Did not find any daughters to the mother of a 94" << std::endl;
1508  }
1509  }
1510 
1511  // Here all is sorted out, update the found history relations.
1512 
1513  for ( int j_mo=1 ; j_mo <= mothers[0] ; j_mo++ ) {
1514 
1515  int mother=mothers[j_mo];
1516 
1517  for ( int j_dau=1 ; j_dau <= daughters[j_mo][0] ; j_dau++ ) {
1518  int daughter=daughters[j_mo][j_dau];
1519  if ( j_dau == 1 ) k[mother][4]=daughter;
1520  k[mother][5]=daughter;
1521  k[daughter][3]=mother;
1522  }
1523 
1524  }
1525  } // endif this is a 94 line
1526  } // for all lines
1527 
1528  streamlog_out(DEBUG3) << " HEPEVT relation table, after 94 fixing" <<std::endl;
1529  streamlog_out(DEBUG3) << " line status pdg parent first last second anti " << std::endl;
1530  streamlog_out(DEBUG3) << " daughter daugther parent colour coulour " << std::endl;
1531 
1532  for ( int j_py=1 ; j_py<=nlund ; j_py++ ) {
1533  streamlog_out(DEBUG3) << std::setw(10) << j_py << std::setw(10) << k[j_py][1] <<
1534  std::setw(10) << k[j_py][2] << std::setw(10) << k[j_py][3] <<
1535  std::setw(10) << k[j_py][4] << std::setw(10) << k[j_py][5] << std::setw(10) << k[j_py][6] <<
1536  std::setw(10) << k[j_py][7] << std::setw(10) << k[j_py][8] << std::endl;
1537  }
1538 }
1539 
1540 void TrueJet::stdhep_reader_bug_workaround( int line94 )
1541 {
1542  int first_94_mother= k[line94][3];
1543  int last_94_mother= k[line94][6];
1544  int first_94_daughter = k[line94][4];
1545  int last_94_daughter = k[line94][5];
1546 
1547 
1548 
1549  bool inconsitent_pid =
1550  (( k[first_94_daughter][2] != k[first_94_mother][2] && k[first_94_daughter][2] != k[last_94_mother][2] ) || // inconsistent pdgs
1551  ( k[ last_94_daughter][2] != k[first_94_mother][2] && k[ last_94_daughter][2] != k[last_94_mother][2] )) ;
1552 
1553  bool inconsitent_E =
1554  ( abs((p[last_94_mother][4]+p[first_94_mother][4]) -p[line94][4])/p[line94][4] > 0.001 ) ; // inconsitent E
1555 
1556  bool gluon_daughters = ( k[first_94_daughter][2] == 21 || k[last_94_daughter][2] == 21 ) ; // 94->gluons is OK as is.
1557 
1558  if ( ( inconsitent_pid || inconsitent_E ) && ( ! gluon_daughters )) {
1559  int right_mother = 0 ;
1560  int sought_mother = 0 ;
1561  int new_mother = 0;
1562  bool first = true;
1563  while ( 1 ) {
1564 
1565 
1566  if ( first ) {
1567  right_mother = first_94_mother;
1568  } else {
1569  right_mother = last_94_mother;
1570  }
1571 
1572  if ( first ) {
1573  if ( inconsitent_pid ) {
1574  if ( k[first_94_daughter][2]*k[last_94_daughter][2] > 0 ) {
1575 
1576  streamlog_out(ERROR) << " 94 DAUGHTERS wrong of 94 on line " << line94
1577  << " : both fermions or both anti-fermions " <<std::endl;
1578  // Probably does not happen ...
1579 
1580  } else if ( k[first_94_mother][2]*k[last_94_mother][2] > 0 ) {
1581 
1582  streamlog_out(DEBUG4) << " 94 MOTHERS wrong of 94 on line " << line94
1583  << " : both fermions or both anti-fermions " <<std::endl;
1584 
1585  // In this case, already the particle or antiparticle nature of the daughter tells
1586  // us the mother of which particle we need
1587  // to find in the record to correct things.
1588 
1589  } else {
1590 
1591  streamlog_out(ERROR) << " 94 on line " << line94 << " is wrong : differnt PDGs, but both mothers and daughters are "
1592  << "fermions/anti-fermion pairs ??? " <<std::endl;
1593  // Probably does not happen ...
1594 
1595  }
1596  }
1597  }
1598 
1599 
1600  if ( k[first_94_daughter][2] == k[right_mother][2] ) {
1601  sought_mother = last_94_daughter;
1602  } else {
1603  sought_mother = first_94_daughter;
1604  }
1605 
1606  new_mother = 0 ;
1607  double target[5]={0};
1608  for ( int jjj=1 ; jjj<=4 ; jjj++ ){target[jjj]=p[line94][jjj]-p[right_mother][jjj];}
1609  for ( int search_line = line94-1 ; search_line >= 1 ; search_line-- ) {
1610  if ( k[search_line][2] == k[sought_mother][2] ) {
1611  if ( abs(p[search_line][1] - target[1])/ abs(target[1]) < 0.05 &&
1612  abs(p[search_line][2] - target[2])/ abs(target[2]) < 0.05 &&
1613  abs(p[search_line][3] - target[3])/ abs(target[3]) < 0.05 ) {
1614  new_mother=search_line;
1615  break ;
1616  }
1617  }
1618  }
1619  if ( new_mother != 0 ) { break ; }
1620  if ( ! first ) { break ; }
1621  first = false;
1622  } // while ( 1 )
1623 
1624 
1625  if ( new_mother == 0 ) {
1626  streamlog_out(ERROR) << " 94 on line " << line94 << " is wrong, but could not figure out un what way. " ;
1627  } else {
1628  if ( new_mother < right_mother ) {
1629  first_94_mother = new_mother ; last_94_mother = right_mother ;
1630  } else {
1631  first_94_mother = right_mother ; last_94_mother = new_mother ;
1632  }
1633  }
1634 
1635  k[line94][3] = first_94_mother ; k[line94][6] = last_94_mother ;
1636  k[first_94_mother][4] = line94 ;k[first_94_mother][5] = line94 ;
1637  k[last_94_mother][4] = line94 ;k[last_94_mother][5] = line94 ;
1638 
1639  } // if(inconsitent_pid || inconsitent_E )
1640 }
1641 void TrueJet::true_lepton()
1642 {
1643  int ihard_lepton_1[4011]={0} ;
1644  int ihard_lepton[4011]={0};
1645  int lept=0;
1646 
1647  //! count number of hard leptons. The condition is: it should be a lepton lepton (pdg 11 to 16).
1648  //! In whiz2, it is hard if the status of its parent is 21, or it comes from a W or a Z.
1649  //! In whiz1, it is hard either if it is comming from line 3, or if comes from line 1, goes to
1650  //! a single daughter which is an electron. (The second part is needed for beam-remenant
1651  //! electrons in odd-f events)
1652 
1653  n_hard_lepton = 0;
1654  int start_loop=1;
1655 
1656  if ( _top_event ) { // start looking for leptons after the last top - there are un-connected leptons before that to avoid!
1657  for ( int k_py=nlund ; k_py>= 1 ; k_py-- ) {
1658  if ( abs(k[k_py][2]) == 6 ) {
1659  start_loop = k_py ;
1660  break ;
1661  }
1662  }
1663  }
1664 
1665  for ( int k_py= start_loop ; k_py<= nlund ; k_py++ ) {
1666  if ( ( abs(k[k_py][2]) >= 11 && abs(k[k_py][2]) <= 16 ) ) {
1667  streamlog_out(DEBUG1) << " found lepton, k_py = " << k_py << ", k[k_py][3] = "
1668  << k[k_py][3] << ", k[3][3] = " << k[3][3] << std::endl ;
1669  if (
1670  (( ! _whiz1 ) && ((k [k[k_py][3]][1] == 21 ) || ( (abs(k [k[k_py][3]][2]) == 23 ||
1671  abs(k [k[k_py][3]][2]) == 24 ) && k[k [k[k_py][3]][3]][1] == 21 )) ) ||
1672  ( _whiz1 && ((k[k_py][3] == 3 && k[3][3] > 0) || (k[k_py][3] == -1) ||
1673  (k[k_py][3] == 1 && k[k_py][4] ==k[k_py][5] && abs(k[k[k_py][5]][2])==11 ))) ) {
1674 
1675  n_hard_lepton++;
1676  ihard_lepton_1[n_hard_lepton]=k_py;
1677  ihard_lepton[n_hard_lepton]=0;
1678  streamlog_out(DEBUG4) << " lepton on line " << k_py <<" hard according to condition 1 " << std::endl;
1679  }
1680  else if (abs(k[k[k_py][3]][2]) == 25 ||
1681  ((abs(k[k[k_py][3]][2]) == 24 || abs(k[k[k_py][3]][2]) == 23) &&
1682  abs(k[k[k[k_py][3]][3]][2]) == 25)) { // lepton from higgs, or from H->ZZ*/WW*
1683  n_hard_lepton++;
1684  ihard_lepton_1[n_hard_lepton]=k_py;
1685  ihard_lepton[n_hard_lepton]=0;
1686  streamlog_out(DEBUG4) << " lepton on line " << k_py <<" hard according to condition 2 " << std::endl;
1687  }
1688  else if ( k[k[k[k_py][3]][3]][1] == 22 && k[k_py][1] < 20 ) { // grandmother has code 22 (= incomming),
1689  // lepton has code < 20
1690  n_hard_lepton++;
1691  ihard_lepton_1[n_hard_lepton]=k_py;
1692  ihard_lepton[n_hard_lepton]=0;
1693  streamlog_out(DEBUG4) << " lepton on line " << k_py <<" hard according to condition 3 " << std::endl;
1694  }
1695  else if ( _top_event && ( abs(k[k[k_py][3]][2]) == 24) ) { // ultimate ( via only W:s ) ancestor is a top quark
1696  int oma = k[k_py][3] ;
1697  while ( 1 ) {
1698  if ( oma < 3 ) break ;
1699  if ( abs(k[k[oma][3]][2]) == 24 ) {
1700  oma=k[oma][3] ;
1701  } else if ( abs(k[k[oma][3]][2]) == 6 ) {
1702  n_hard_lepton++;
1703  ihard_lepton_1[n_hard_lepton]=k_py;
1704  ihard_lepton[n_hard_lepton]=0;
1705  streamlog_out(DEBUG4) << " lepton on line " << k_py <<" hard according to condition 4 " << std::endl;
1706  break;
1707  } else {
1708  break ;
1709  }
1710  }
1711  }
1712  }
1713  }
1714 
1715  // here n_hard_leptons is the number of hard leptons (i.e. directly from the initial state of from a boson).
1716  // and ihard_lepton_1[..] contains the list of pyjets lines of these
1717 
1718  if ( n_hard_lepton == 0 ) return ;
1719 
1720  //! Sort the leptons in "color singlets", ie. in groups of flavour/anti-flavour. In most events this is
1721  //! the order they come in, but in the case of all flavour being the same, it isn't
1722 
1723  // Outcome of this loop is the list of hard leptons in ihard_lepton which is ordered in pairs
1724  // from the same boson, or if that cant't be figured out, by flavour - anti-flavour pairs.
1725 
1726  for (int j_lept= 1; j_lept<=n_hard_lepton; j_lept++ ) {
1727  if ( ihard_lepton[j_lept] == 0 ) {
1728  for ( int l_lept=1; l_lept<=n_hard_lepton ; l_lept++ ) {
1729  if ( ihard_lepton_1[l_lept] != 0 ) {
1730  ihard_lepton[j_lept]= ihard_lepton_1[l_lept];
1731  ihard_lepton_1[l_lept]= 0;
1732  break ;
1733  }
1734  }
1735 
1736  int k_py = ihard_lepton[j_lept] ;
1737  int n_hint = 0;
1738  int n_boson = 0;
1739  int hint_companion = 0 ;
1740  int boson_companion = 0 ;
1741  int fallback_companion = 0;
1742  for ( int l_lept=1; l_lept <= n_hard_lepton ; l_lept++ ) {
1743 
1744  int l_py = ihard_lepton_1[l_lept] ;
1745 
1746  if ( flavour(k[l_py][2]) == flavour(k[k_py][2]) && k[l_py][2]*k[k_py][2] <0 ) { // consider only
1747  // lepton-flavour/anti-lepton flavour groupings
1748  if ( ( k[l_py][2] > 0 && k[l_py][7] == k[k_py][8] && k[l_py][7] != 0 ) ||
1749  ( k[l_py][2] < 0 && k[l_py][8] == k[k_py][7] && k[l_py][8] != 0 ) ){// a hint on grouping (leptons from the same 94)
1750  streamlog_out(DEBUG4) << " hint match for lepton on lines " << l_py << " and " << k_py << std::endl;
1751  hint_companion = l_lept;
1752  n_hint++ ;
1753  }
1754  if ( k[l_py][3] == k[k_py][3] ) {
1755  if ( k[k[l_py][3]][2] == 23 || k[k[l_py][3]][2] == 24 || k[k[l_py][3]][2] == 25 ) { // leptons from the same boson
1756  streamlog_out(DEBUG4) << " boson match for lepton on lines " << l_py << " and " << k_py << std::endl;
1757  boson_companion = l_lept;
1758  n_boson++;
1759  }
1760  if ( k[l_py][7] == 0 && k[l_py][8]==0 && k[k_py][7] == 0 && k[k_py][8]==0 ) {// if no hint nor boson parent found,
1761  // use this: at least it has the right flavour.
1762  fallback_companion = l_lept ;
1763  }
1764  }
1765  }
1766  }
1767 
1768  if ( n_hint == 1 ) {
1769 
1770  ihard_lepton[j_lept+1]= ihard_lepton_1[hint_companion] ;
1771  ihard_lepton_1[hint_companion]=0;
1772 
1773  } else if (n_boson == 1 ) {
1774 
1775  ihard_lepton[j_lept+1]=ihard_lepton_1[boson_companion] ;
1776  ihard_lepton_1[boson_companion]=0;
1777 
1778  } else {
1779 
1780  if ( n_hint > 1 || n_boson > 1 ) {
1781  streamlog_out(DEBUG4)<< " Problem finding companion lepton to lepton on line " << k_py
1782  << " : n_boson = " << n_boson << " , n_hint = " << n_hint << std::endl;
1783  }
1784 
1785  if (fallback_companion != 0 ) { // no hint nor boson parent: take any lepton with the right
1786  // lepton flavour as companion
1787 
1788  ihard_lepton[j_lept+1]=ihard_lepton_1[fallback_companion];
1789  ihard_lepton_1[fallback_companion]=0;
1790 
1791  }
1792  }
1793  }
1794  }
1795 
1796  // To summerise: at this point we have in ihard_lepton[1 - n_hard_lepton ] the pyjets line numbers of
1797  // all found hard leptons. These can come directly from the in-state or from a boson. If they are involved in
1798  // a CM-shower, the ones *going in* into it are the ones in ihard_lepton.
1799 
1800  //! assign to jets. Basically the jet is the index in the list, -ve if it later on goes into a
1801  //! 94, +ve otherwise.
1802 
1803  // outcome of this loop is the -(jet number) for each hard lepton in jet[ (each hard-lepton pyjets-line) ] and the elementon
1804  // of the corresponding jet elementon[jet] = hard-lepton pyjets-line. lept on exit will be the first post-CM shower
1805  // incarnation of the lepton.
1806 
1807  if ( n_hard_lepton%2 != 0 ) {
1808  if ( abs(k[ihard_lepton[n_hard_lepton]][2]) != 22 ) {
1809 
1810  // If there is an odd number of initial leptons or photons, possibly we are looking at an
1811  // odd-fermion sample. In that case, assume - resonably - that the beam-remnant is the first
1812  // hard lepton seen. Make it the last instead. That way it will be the one not grouped.
1813  //
1814  // TODO: Keep an eye on this : gamma in ME (not tested) ?
1815  // Also: better condition for this case (best would be process-type from the run header)
1816 
1817  int beamrem=ihard_lepton[1];
1818  for ( int j_lept=2; j_lept<=n_hard_lepton ; j_lept++ ) {
1819  ihard_lepton[j_lept-1]=ihard_lepton[j_lept];
1820  }
1821  ihard_lepton[n_hard_lepton]=beamrem;
1822  }
1823  }
1824  for ( int j_lept=1; j_lept<=n_hard_lepton ; j_lept++ ) {
1825 
1826  lept = ihard_lepton[j_lept];
1827  streamlog_out(DEBUG4) << " assigning leptons to jets: ihard_lepton[" << j_lept << "] = " << ihard_lepton[j_lept]
1828  << " jet : " << current_jet+j_lept <<std::endl ;
1829  }
1830  for ( int j_lept=1; j_lept<=n_hard_lepton ; j_lept++ ) {
1831 
1832  lept = ihard_lepton[j_lept];
1833  streamlog_out(DEBUG4) << " assigning leptons to jets: ihard_lepton[" << j_lept << "] = " << ihard_lepton[j_lept]
1834  << " and has PDG " << k[lept][2] << " jet : " << current_jet+j_lept<< std::endl ;
1835 
1836  // ! initially, set to -j_lept
1837 
1838  jet[lept]=-(current_jet+j_lept);
1839 
1840  fafp_last[current_jet+j_lept]=lept; // so it goes when too specific variable names were choosen ...
1841 
1842 
1843  if ( _top_event ) { // we want to back-track the lepton to the top, so that jet grouping can group all top daughters into on "c.n."
1844  int oma = lept ;
1845  int fafp_top = 0;
1846  while ( 1 ) {
1847  if ( oma < 3 ) break ;
1848  if ( abs(k[k[oma][3]][2]) == 24 ) {
1849  oma=k[oma][3] ;
1850  } else if ( abs(k[k[oma][3]][2]) == 6 ) {
1851  fafp_top = k[oma][3];
1852  oma=k[oma][3] ;
1853  } else {
1854  break ;
1855  }
1856  }
1857  if ( fafp_top != 0 ) { fafp_last[current_jet+j_lept]=fafp_top; }
1858  }
1859 
1860  // ! loop until either stable descendent found
1861 
1862  while ( k[lept][2] != 94 && k[lept][1] != 1 && k[lept][4] == k[lept][5] ) {
1863  lept = k[lept][4] ;
1864  }
1865 
1866  streamlog_out(DEBUG2)<< " after looping to stable : " << lept << " " << k[lept][2] << std::endl;
1867 
1868  while ( true ) {
1869  int ida=0;
1870  bool is_elementon = true ;
1871  if ( k[lept][1] != 1 && ( k[lept][4] == k[lept][5] && k[lept][4] != 0 )) {
1872  ida=k[lept][4] ;
1873  if ( k[ida][2] == k[lept][2] ) {
1874  is_elementon = false;
1875  }
1876  }
1877 
1878  if ( is_elementon ) {
1879 
1880  jet[lept]=current_jet+j_lept ;
1881  elementon[current_jet+j_lept]=lept ;
1882  break;
1883  } else {
1884 
1885  lept=ida ;
1886 
1887  }
1888 
1889  }
1890 
1891  }
1892 
1893  // ! set jet for all further descendants of hard leptons . Don't change already set assignments. Never assign a jet to a 94.
1894 
1895  for ( int k_py=1; k_py<=nlund ; k_py++ ) {
1896  if ( jet[k_py] == 0 && k[k_py][2] != 94 && abs(jet[k[k_py][3]]) > current_jet && abs(jet[k[k_py][3]]) <= current_jet+n_hard_lepton ) {
1897 
1898  // ! this is indeed a descendant of a hard lepton
1899  if ( k[k_py][1] == 1 ) {
1900  jet[k_py] = abs(jet[k[k_py][3]]);
1901  if ( k[k_py][2] == 22 && k[k[k_py][3] ][2] == k[elementon[jet[k_py]]][2] ) {
1902  nfsr[jet[k_py]]++;
1903  }
1904  } else {
1905  jet[k_py] = jet[k[k_py][3]];
1906  }
1907  }
1908  }
1909  current_jet=current_jet+n_hard_lepton ;
1910 
1911 }
1912 
1913 void TrueJet::cluster()
1914 {
1915  int clus[4011]={0};
1916  int clu=0 , jet1=0, jet2=0;
1917 
1918  nclu=0;
1919  for (int k_py=1; k_py<=nlund ; k_py++ ) {
1920  streamlog_out(DEBUG0) << " Checking cluster at " << k_py << " k[k_py][2], and k[k_py][1] : " <<
1921  " " << k[k_py][2] << " , " << k[k_py][1] << std::endl;
1922  if ( k[k_py][2] == 91 && k[k_py][1] < 30 ) {
1923  streamlog_out(DEBUG4) << " cluster found at " << k_py <<std::endl;
1924  nclu++;
1925  clus[nclu]=k_py;
1926  }
1927  }
1928  if ( nclu == 0 ) return;
1929 
1930  for ( int i_clu=1 ; i_clu<=nclu ; i_clu++ ) {
1931  clu=clus[i_clu];
1932 
1933  //! Jet-assignment: Decide if the particle should be
1934  //! assigned to the first or last quark end of the cluster, by checking
1935  //! the sign of the projection of the hadron momentum on the difference
1936  //! between the two quark directions.
1937 
1938  //! First of all: find first and second quark going into the cluster, normally the first is
1939  //! the mother of the cluster, and the second is the line after, or possibly
1940  //! a few linese later.
1941 
1942  jet1=current_jet+2*(i_clu-1)+1 ; jet2= current_jet+2*(i_clu-1)+2;
1943  elementon[jet1]=k[clu][3];
1944  elementon[jet2]=elementon[jet1];
1945  while ( k[elementon[jet2]+1][4] == clu ) {
1946  elementon[jet2]=elementon[jet2]+1;
1947  }
1948  streamlog_out(DEBUG4) << " Assigning jets " << jet1 << " or " << jet2
1949  << " to CLUSTER at " << clu << std::endl;
1950  assign_jet(jet1,jet2,clu) ;
1951 
1952  }
1953  current_jet=current_jet+ 2*nclu;
1954 
1955 }
1956 
1957 void TrueJet::photon()
1958 {
1959  int phots[4011]={0};
1960 
1961  nphot=0;
1962  for (int k_py=1; k_py<=nlund ; k_py++ ) {
1963  streamlog_out(DEBUG0) << " Checking photon at " << k_py << " k[k_py][2], and k[k_py][1] : " <<
1964  " " << k[k_py][2] << " , " << k[k_py][1] << std::endl;
1965  if ( k[k_py][2] == 22 && k[k_py][1] < 30 ) {
1966  if ( ( ! _whiz1 ) && (k [k[k_py][3]][1] == 21 ) ) {
1967  nphot++;
1968  phots[nphot]=k_py;
1969  streamlog_out(DEBUG4) << " photon on line " << k_py <<" hard according to condition 1 " << std::endl;
1970  }
1971  else if (abs(k[k[k_py][3]][2]) == 25 ) {
1972  nphot++;
1973  phots[nphot]=k_py;
1974  streamlog_out(DEBUG4) << " photon on line " << k_py <<" hard according to condition 2 " << std::endl;
1975  }
1976  }
1977  }
1978 
1979  if ( nphot == 0 ) return;
1980 
1981  for ( int i_phot=1 ; i_phot<=nphot ; i_phot++ ) {
1982 
1983  // ! jet assignment - index in the list. -ve if "decayimg", which might be the case
1984  // ! for explicitly requested gammas.
1985 
1986  elementon[current_jet+i_phot] = phots[i_phot];
1987  fafp_last[current_jet+i_phot] = phots[i_phot];
1988  if ( k[ phots[i_phot]][1] == 1 ) {
1989  jet[ phots[i_phot]]=current_jet+i_phot;
1990  } else {
1991  jet[ phots[i_phot]]=-(current_jet+i_phot);
1992  }
1993  }
1994  for ( int k_py=1; k_py<=nlund ; k_py++ ) {
1995  if ( jet[k_py] == 0 && abs(jet[k[k_py][3]]) > current_jet && abs(jet[k[k_py][3]]) <= current_jet+nphot ) {
1996 
1997  // ! this is indeed a descendant of a hard photon
1998 
1999  jet[k_py] = abs(jet[k[k_py][3]]);
2000  }
2001  }
2002 
2003  current_jet=current_jet+ nphot;
2004 
2005 }
2006 
2007 void TrueJet::isr()
2008 {
2009  int iisr[4011]={0};
2010 
2011  nisr=0;
2012  for (int k_py=1; k_py<=nlund ; k_py++ ) {
2013  if ( (( ! _whiz1 ) && ( (k[k_py][3] == 3 || k[k_py][3] == 4 ) && ( k[k_py][2]==22 ) && ( k[k_py][1] < 20 ) )) ||
2014  (( _whiz1 ) && ( (k[k_py][3] == 1 || k[k_py][3] == 3 || k[k_py][3] == 0 ) && ( k[k_py][2]==22 ) &&
2015  (k[k_py][4]==k[k_py][5]) && (k[k[k_py][4]][2] ==22) && (k[k[k_py][4]][1] != 0) ) )){
2016  nisr++;
2017  iisr[nisr] = k_py ;
2018  }
2019  }
2020  if ( nisr > 0 ) {
2021  // ! jet assignment - index in the list. -ve if "decayimg", which might be the case
2022  // ! for explicitly requested gammas.
2023 
2024  for ( int j_isr=1 ; j_isr<=nisr ; j_isr++ ) {
2025  elementon[current_jet+j_isr] = iisr[j_isr];
2026  fafp_last[current_jet+j_isr] = iisr[j_isr];
2027  if ( k[ iisr[j_isr]][1] == 1 ) {
2028  jet[ iisr[j_isr]]=current_jet+j_isr;
2029  } else {
2030  jet[ iisr[j_isr]]=-(current_jet+j_isr);
2031  }
2032  }
2033 
2034  for ( int k_py=1; k_py<=nlund ; k_py++ ) {
2035  if ( jet[k_py] == 0 && abs(jet[k[k_py][3]]) > current_jet && abs(jet[k[k_py][3]]) <= current_jet+nisr ) {
2036 
2037  // ! this is indeed a descendant of an isr
2038 
2039  jet[k_py] = abs(jet[k[k_py][3]]);
2040  }
2041  }
2042  current_jet=current_jet+nisr;
2043  }
2044 }
2045 
2046 void TrueJet::string()
2047 {
2048  int n_left=0, sstr=0,lquark=0,lstr=0, istr=0,str_nb=0,jet1=0,jet2=0,istr_previous=0,iback=0;
2049  int str_index[4011]={0}, rev_index[4011]={0};
2050 
2051  //!! Move jets away to make room for string-induced ones at the lowest indicies
2052 
2053  for ( int k_py=1 ; k_py<=nlund ; k_py++ ) {
2054  rev_index[k_py]=0;
2055  if ( jet[k_py] != 0 ) {
2056  jet[k_py]=(jet[k_py]>0 ? jet[k_py]+2*nstr : jet[k_py]-2*nstr);
2057  }
2058  }
2059 
2060  if ( nstr > 0 ) {
2061  for (int j_jet=njet-2*nstr ; j_jet>=1 ; j_jet-- ) {
2062  elementon[j_jet+2*nstr]=elementon[j_jet] ;
2063  elementon[j_jet]=0;
2064  nfsr[j_jet+2*nstr]=nfsr[j_jet] ;
2065  nfsr[j_jet]=0;
2066  fafp_last[j_jet+2*nstr]=fafp_last[j_jet];
2067  fafp_last[j_jet]=0;
2068  boson[j_jet+2*nstr]=boson[j_jet];
2069  boson[j_jet]=0;
2070  fafp_boson[j_jet+2*nstr]=fafp_boson[j_jet];
2071  fafp_boson[j_jet]=0;
2072  }
2073  }
2074 
2075  current_jet=0 ;
2076 
2077  //! Decode the event record. First find the line of the first string,
2078  //! the last quark of the last string (lquark, at the line before the first string),
2079  //! the last string (lstr, the daughter of the lquark quark), and the number of strings (nstr).
2080 
2081 
2082  for ( int k_py=1 ; k_py<= nlund ; k_py++ ) {
2083  streamlog_out(DEBUG0) << " particle " << k_py << " with PDG " << k[k_py][2]
2084  << " is assigned to jet " << jet[k_py] << std::endl;
2085 
2086  if (jet[k_py]==0 && k[k_py][2] != 91 && k[k_py][2] != 94 && k[k_py][2] != 21 && k[k_py][2] != 23 &&
2087  k[k_py][2] != 24 && k[k_py][2] != 25 && k[k_py][1] < 30 ){ // only look for fermions or gammas (not already
2088  // assigned to jets because they are leptons, from
2089  // clusters, isr, whatever...) in the P.S.
2090  n_left++;
2091  str_index[n_left]=k_py;
2092  rev_index[k_py]=n_left;
2093  if ( sstr == 0 && k[k_py][2] == 92 ) {
2094  sstr=n_left;
2095  }
2096  }
2097  }
2098 
2099 
2100  // sstr is the index in str_index of the first 92 found, and lquark is the quark on the line before that.
2101  // This is the last quark of the last string in the event.
2102 
2103  lquark=str_index[sstr-1];
2104 
2105  if (k[k[lquark][4]][2] != 92 ) { // there might be an intruder between the lquark and the first string. This is taken care of here
2106  if ( ! _higgs_to_glue_glue ) {
2107  int lll = lquark ;
2108  while ( lll >= k[str_index[sstr]][3] ) {
2109  if ( k[lll][4] == str_index[sstr] ) {
2110  break ;
2111  } else {
2112  lll-- ;
2113  }
2114  }
2115  if ( lll > k[str_index[sstr]][3] ) {
2116  lquark = lll ;
2117  } else {
2118  // FIXME : this happens in H->glueglue, Need radical re-thinking for that.
2119  streamlog_out(ERROR) << " INFO: Non-contigous string ancestors (1)" << std::endl;
2120  streamlog_out(ERROR) << " Event: " << evt->getEventNumber() << ", run: " << evt->getRunNumber() << std::endl ;
2121  streamlog_out(ERROR) << std::endl ;
2122  }
2123  }
2124  }
2125 
2126 
2127  // lstr is the last string and is the daughter of the last quark of the last string
2128 
2129  lstr=k[lquark][4] ; istr=lstr ; str_nb = nstr ;
2130  streamlog_out(DEBUG3) << " sstr, lstr, lquark, current_jet, str_nb " << sstr << " " << lstr << " "
2131  << lquark << " " << current_jet << " " << str_nb << std::endl;
2132 
2133  //! now the jet-assignment. Start by looking at the last string (lstr),
2134  //! stop once the first has been reached. decide if the particle should be
2135  //! assigned to the first or last quark end of the string, by checking
2136  //! the sign of the projection of the hadron momentum on the difference
2137  //! between the two quark directions.
2138 
2139  while( istr >= str_index[sstr]) {
2140 
2141  streamlog_out(DEBUG1) << " calculating jet1/2 from current_jet = " << current_jet << ", str_nb = " << str_nb << std::endl;
2142  streamlog_out(DEBUG1) << " istr = " << istr << ", k[istr][3] = " << k[istr][3] << ", lquark = " << lquark << std::endl;
2143 
2144  jet1= current_jet+2*(str_nb-1)+1 ; jet2=current_jet+2*(str_nb-1)+2 ;
2145 
2146  elementon[jet1]=k[istr][3] ; elementon[jet2]=lquark ; // elementons: the last quarks before hadronisation, the ones
2147  // jet-numbering pivots around.
2148 
2149  streamlog_out(DEBUG4) << " Assigning jets " << jet1 << " or " << jet2 << " to STRING at " << istr << std::endl;
2150 
2151  assign_jet(jet1,jet2,istr) ; // assigns jets *both* to ancestors and decendants of the elementons
2152 
2153  streamlog_out(DEBUG4) << " Assigned jets " << std::endl;
2154 
2155  //! the previous string is the mother of its hadrons, eg. the last one, which is
2156  //! on the line before the current string:
2157 
2158  //! new string:
2159 
2160  istr_previous=istr ;
2161  iback=str_index[rev_index[istr]-1] ; istr = k[iback][3] ;
2162  if ( istr < str_index[sstr] ) break ; // before the py-line of the first string -> done
2163 
2164  //! back-up: end-quark of the previous string is the line before the start-quark
2165  //! of the present string:
2166 
2167  int lquark_start = str_index[rev_index[k[istr_previous][3]]-1] ;
2168  lquark=str_index[rev_index[k[istr_previous][3]]-1] ;
2169 
2170  if ( k[lquark][4] != istr && lquark > 0 ) {
2171 
2172  //!! Once again the very unusual case: need to search for the end of the new string backwards
2173  //!! print a message and the event for now. This one actually happens now and then: If there
2174  //!! is a cluster in the event, there might be a part of the parton shower after the direct parents
2175  //!! of one string that eventually ends up in a later _string_ (not a cluster, that would have
2176  //!! been taken care of),
2177  streamlog_out(DEBUG3) << " Non-contigous string ancestors (2) " << std:: endl;
2178  while ( lquark>-1 && k[lquark][4] != istr ) {
2179  lquark=lquark-1;
2180  }
2181 
2182  if ( lquark <= 0 ) {
2183  lquark = lquark_start ;
2184  while ( lquark<nlund && ( k[lquark][4] != istr || lquark == k[istr][3] )) {
2185  lquark=lquark+1;
2186  }
2187  if ( lquark >= nlund ) {
2188  streamlog_out(ERROR) << " ERROR: Quack ?! " << std::endl ;
2189  streamlog_out(ERROR) << " Event: " << evt->getEventNumber() << ", run: " << evt->getRunNumber() << std::endl ;
2190  streamlog_out(ERROR) << std::endl ;
2191  }
2192  }
2193  }
2194 
2195  //! decrement the string-counter
2196 
2197  str_nb--;
2198 
2199  }
2200 }
2201 
2202 void TrueJet::assign_jet(int jet1,int jet2,int this_string)
2203 {
2204  double dir_diff[4]={0};
2205  int first_p1=0, first_p2=0;
2206  streamlog_out(DEBUG4) << " assign_jet: elementon[" << jet1 << "] = " << elementon[jet1]
2207  << ", elementon[" << jet2 << "] = " << elementon[jet2] << std::endl;
2208 
2209  // elementon is the last quark before hadronisation (i.e. either of the two going into to 92 objects).
2210  // jet1/2 is the jet-number of that one. The output is first_p1/2 which is the first quark in the
2211  // chain leading up to the elementon. fafp_last is the last quark before any 94 or boson (last in the sense
2212  // of going back the history!) - often the same as first_p, while nfsr is the number of FSRs in the jet,
2213  // boson is the line of any boson in the P.S. yieding a qqbar pair (typically a gluon), and fafp_boson is
2214  // the ultimate ancestor of the boson.
2215 
2216  first_parton( elementon[jet1], jet1 , first_p1, fafp_last[jet1], nfsr[jet1], boson[jet1], fafp_boson[jet1]);
2217  streamlog_out(DEBUG4) << " back in assign_jets for jet " << jet1 << " fafp_last = " << fafp_last[jet1]
2218  << " nfsr = " << nfsr[jet1]<< " boson = " << boson[jet1]
2219  << " fafp_boson = " << fafp_boson[jet1] << std::endl;
2220 
2221  first_parton( elementon[jet2], jet2 , first_p2, fafp_last[jet2], nfsr[jet2], boson[jet2], fafp_boson[jet2]);
2222  streamlog_out(DEBUG4) << " back in assign_jets for jet " << jet2 << " fafp_last = " << fafp_last[jet2]
2223  << " nfsr = " << nfsr[jet2]<< " boson = " << boson[jet2]
2224  << " fafp_boson = " << fafp_boson[jet2] << std::endl;
2225 
2226  // ! difference between the two quark directions: one is at elementon(jet1), the
2227  // ! other is elementon(jet2) = three lines of fortran, 20 lines of C++ ...
2228  // ! The closest match will determine which elementon the final hadrons are assigned to
2229  // ! and hence the jet-numbering.
2230 
2231 
2232  double absp_1=0;
2233  double absp_2=0;
2234 
2235  for ( int kk=1; kk<=3 ; kk++ ) {
2236  absp_1+=p[elementon[jet1]][kk]*p[elementon[jet1]][kk] ;
2237  absp_2+=p[elementon[jet2]][kk]*p[elementon[jet2]][kk] ;
2238  }
2239 
2240  absp_1=sqrt(absp_1);
2241  absp_2=sqrt(absp_2);
2242  for ( int kk=1; kk<=3 ; kk++ ) {
2243  dir_diff[kk] =
2244  p[elementon[jet2]][kk]/absp_2- p[elementon[jet1]][kk]/absp_1;
2245  }
2246 
2247  for (int k_py_had=k[this_string][4]; k_py_had<=k[this_string][5] ; k_py_had++ ) { // this loops over the first generation hadrons
2248  if ( k[k_py_had][3] == this_string ) {
2249  double dot=0;
2250  for ( int jj=1 ; jj<=3 ; jj++ ) {
2251  dot+=dir_diff[jj]*p[k_py_had][jj] ;
2252  }
2253  if ( dot <= 0. ) { // by projcting on the *difference* between the two elelementon directions,
2254  // one just need to check the sign to know which is closest!
2255  jet[k_py_had] = jet1 ;
2256  } else {
2257  jet[k_py_had] = jet2 ;
2258  }
2259  }
2260  }
2261 
2262  for ( int k_py=k[elementon[jet2]][5]+1; k_py <=nlund ; k_py++ ) { // k[elementon[jet2]][5]+1 is the first
2263  // hadron of the string (actually any of
2264  // elementon[jet1/jet2][4/5] would do...)
2265 
2266  if ( jet[k_py] == 0 && k[k_py][1] < 30 ) { // if not already assigned, and not overlay. NB that the first
2267  // generation already was assigned above.
2268  jet[k_py] = abs(jet[k[k_py][3]]); // i.e. jet is the same as it's mother - fair enough !
2269  if ( k[k_py][2] == 21 || abs( k[k_py][2])<=6 ) {
2270  // can happen if a paricle decay is done by gluons (Ypsilon etc.)
2271  jet[k_py] = -jet[k_py];
2272  }
2273  if ( jet[k_py] != 0 ) streamlog_out(DEBUG1) << " Particle " << k_py
2274  << " assigned to jet " << abs(jet[k[k_py][3]]) << std::endl;
2275  }
2276  }
2277 
2278 }
2279 
2280 void TrueJet::first_parton(int this_partic,int this_jet,int& first_partic,int& last_94_parent,int& nfsr_here,int& info,int& info2)
2281 {
2282  int this_quark=0,quark_pdg=0,boson_ancestor=0,istat=0,boson_last_94_parent=0,istat2=0;
2283 
2284  this_quark=this_partic;
2285  jet[this_quark]=-this_jet;
2286  quark_pdg=k[this_quark][2];
2287 
2288  streamlog_out(DEBUG4) << " first_parton: this_quark = " << this_quark
2289  << ", quark_pdg = " << quark_pdg << std::endl;
2290 
2291 
2292  last_94_parent = 0;
2293 
2294  //! back-track quarks to the begining of the event record, or until the quark-chain is broken
2295  //! Beginning of record is hit when the parent of a quark is an electron, photon or a whizard inserted boson
2296  //! and the chain is brocken if the parent of a quark is a gluon (21) or a W (24, non-resonanc-insertion).
2297 
2298  //! start at this_quark, which might on entry actually be something else (a gluon or a W)
2299 
2300  while ( k[this_quark][3] > 0 &&
2301  abs(k[k[this_quark][3]][2]) != 11 && abs(k[k[this_quark][3]][2]) != 21 && abs(k[k[this_quark][3]][2]) != 22
2302  && abs(k[k[this_quark][3]][2]) != 23 && abs(k[k[this_quark][3]][2]) != 24 && abs(k[k[this_quark][3]][2]) != 25 ) {
2303  // (this end-condition works both for whiz1 and whiz2)
2304 
2305  streamlog_out(DEBUG3) << " start of while loop: this_quark = " << this_quark
2306  << ", its PDG = " << k[this_quark][2]
2307  << ", its parent = " << k[this_quark][3]
2308  << ", parent's PDG = " << abs(k[k[this_quark][3]][2])
2309  << ", quark_pdg = " << quark_pdg << std::endl;
2310 
2311  //! step back in the history
2312 
2313  this_quark=k[this_quark][3];
2314  if ( quark_pdg == 21 || abs(quark_pdg) == 24 ) {
2315  //! looking at a gluon or a W. We've alredy steped to the
2316  //! parent of the boson (which should be a quark), so now we might
2317  //! be loking at a different flavour:
2318  //!
2319 
2320  quark_pdg=k[this_quark][2];
2321 
2322  }
2323 
2324  //! jet-assignment if requested. Set to -ve to indicate that we are in the
2325  //! parton shower,
2326 
2327  if ( this_jet != 0 ) {jet[this_quark]=-this_jet;}
2328 
2329  //! Internal brems ? If the line after this quark is a photon, and this quark and
2330  //! the photon has the same mother, it is an FSR
2331 
2332  if ( k[this_quark+1][2] == 22 && k[this_quark+1][3] == k[this_quark][3] ) {
2333 
2334  streamlog_out(DEBUG5) << " FSR from " << this_quark << std::endl;
2335  if ( this_jet != 0 ) {
2336 
2337  //! assign the FST to the jet of the quark. It's a real particle, so it
2338  //! assigned a +ve jet number
2339 
2340  jet[this_quark+1]=this_jet;
2341  nfsr_here=nfsr_here+1;
2342  }
2343  }
2344 
2345  streamlog_out(DEBUG3) << " end of while loop: this_quark = " << this_quark
2346  << ", its PDG = " << k[this_quark][2]
2347  << ", its parent = " << k[this_quark][3]
2348  << ", parent's PDG = " << abs(k[k[this_quark][3]][2])
2349  << ", quark_pdg = " << quark_pdg << std::endl;
2350  }
2351 
2352  // check if a boson parent is a "real" boson or a resonance insertion, which is flagged by the
2353  // fact that the parent of the boson is a beam-particle (electron or photon), with status code 21,
2354  // Resonance insertion were never done in Whiz1, so it is no problem that statuscode 21 is not
2355  // assigned for Whiz1 - the condition should always evaluate to false, anyhow!
2356  //
2357  // We also stop in the parent of the boson is a higgs, both for Whiz1 and 2.
2358  int resonance_insertion = 0;
2359  streamlog_out(DEBUG3) << " resonance_insertion check : " << abs(k[k[this_quark][3]][2]) << std::endl;
2360  if ( abs(k[k[this_quark][3]][2]) == 23 || abs(k[k[this_quark][3]][2]) == 24 ) {
2361  int lanc = k[k[this_quark][3]][3];
2362  streamlog_out(DEBUG1) << " lanc = " << lanc << " " << k[lanc][2] << " " << k[lanc][1] << std::endl ;
2363  if ( ( ( abs(k[lanc][2]) == 11 || abs(k[lanc][2]) == 22 ) && k[lanc][1] == 21 ) ||
2364  ( abs(k[lanc][2]) == 25 ) ){ // This is not a "real" boson, but a technical resonance insert
2365  resonance_insertion = 1;
2366  }
2367  }
2368 
2369  streamlog_out(DEBUG3) << " resonance_insertion = " << resonance_insertion << std::endl;
2370 
2371  //! we are here because we've reached the end of the chain. If this happened
2372  //! becasue the parent was a boson radiated during the PS, we recurse. (Else we arrived at the
2373  //! initial state and are done.)
2374 
2375  if ( resonance_insertion == 0 && ( abs(k[k[this_quark][3]][2]) == 21 || abs(k[k[this_quark][3]][2]) == 24 )) {
2376  //! recurse.
2377 
2378  //! jet-number set to 0 => do not assign jet numbers anymore. The untimate first
2379  //! partic will go into boson_ancestor, parent of the last 94 into boson_last_94_parent.
2380 
2381  streamlog_out(DEBUG3) << " in while loop calling first_parton with: k[this_quark][3] = " << k[this_quark][3] << std::endl;
2382  first_parton(k[this_quark][3],0 ,boson_ancestor, boson_last_94_parent, nfsr_here, istat,istat2);
2383 
2384  //! the way this bottoms-out: second to last argument (istat in the call) is set to 0 if
2385  //! we didn't find a boson.. In that case ...
2386 
2387  if ( istat == 0 ) {
2388 
2389  //! ... we end up here. Here we set the argument to something different from 0
2390  //! (namely boson_ancestor), so whoever called me in this step will ...
2391 
2392  info= boson_ancestor;
2393  info2= boson_last_94_parent;
2394 
2395  } else {
2396 
2397  //! ... end up here, where the same information is further transmitted up
2398  //! the stack.
2399 
2400  info = istat;
2401  info2 = istat2;
2402 
2403  }
2404  //! ... and so on until we are back at the top of the stack.
2405 
2406  } else {
2407 
2408  //! bottom-out
2409 
2410  info=0;
2411  info2=0;
2412 
2413  }
2414 
2415  if ( last_94_parent == 0 ) {last_94_parent=this_quark;}
2416 
2417  first_partic = this_quark;
2418  streamlog_out(DEBUG4) << " first_parton, return: = first_partic " << first_partic
2419  << ", last_94_parent = " << last_94_parent << std::endl;
2420 
2421 
2422 }
2423 
2424 int TrueJet::flavour(int k2)
2425 {
2426  // note that flavour is a bit ambigous: for quarks it is the type of the quark (6 values in the SM),
2427  // for leptons it is the family of the lepton (3 valuse in the SM). Here we mean "family" in both cases.
2428  int k2l=0;
2429 
2430  if ( abs(k2)>16) return 0 ;
2431  if ( abs(k2) > 10 ) {
2432  k2l=abs(k2)-11;
2433  } else {
2434  k2l=abs(k2)-1;
2435  }
2436  return k2 > 0 ? k2l/2 : k2l/2 ;
2437  // sign(k2l/2,k2);
2438 }
2439 
2440 void TrueJet::grouping(){
2441  int ijet=0 ;
2442  for ( int k_jet=1 ; k_jet<=njet ; k_jet++ ) {
2443  tE[k_jet]=0;
2444  for ( int i=0 ; i<3 ; i++ ) {
2445  tmom[k_jet][i]=0. ;
2446  }
2447  }
2448  for ( int i_py=1 ; i_py<=nlund ; i_py++){
2449  if ( jet[i_py]>0 ) {
2450  if ( k[i_py][1] == 11 ) {
2451  double Ekid=0;
2452  for ( int jj=k[i_py][4] ; jj <= k[i_py][5] ; jj++ ) {
2453  Ekid+= p[jj][4];
2454  }
2455  if ( abs(Ekid-p[i_py][4])/p[i_py][4] > 0.001 ) {
2456  streamlog_out(DEBUG3) << " Particle " << i_py << " has energy " << p[i_py][4] <<
2457  " , but the sum of its genstat1 kids is " << Ekid << std::endl;
2458  streamlog_out(DEBUG3) << " indicating that Geant did something fishy and un-documented in MCParticle "
2459  << std::endl ;
2460  streamlog_out(DEBUG3) << " (Counter-meassures taken, so it should be OK.) " << std::endl ;
2461  streamlog_out(DEBUG3) << " Event: " << evt->getEventNumber() << ", run: " << evt->getRunNumber() << std::endl;
2462  streamlog_out(DEBUG3) << std::endl ;
2463  const double* v = mcp_pyjets[k[i_py][4]]->getVertex();
2464  if ( v[0]*v[0] + v[1]*v[1] > 10. ) {
2465  if ( k[i_py][1] < 30 ) {
2466  k[i_py][1]=1;
2467  }
2468  for ( int j_dau=k[i_py][4] ; j_dau <= k[i_py][5] ; j_dau++ ) {
2469 
2470  if ( k[j_dau][1] < 30 ) {
2471  k[j_dau][1]=0;
2472  }
2473  }
2474  }
2475  }
2476 
2477  } else if ( k[i_py][1] == 0 ) {
2478  for ( unsigned i_dau=0 ; i_dau < mcp_pyjets[i_py]->getDaughters().size() ; i_dau++ ) {
2479  int j_dau=mcp_pyjets[i_py]->getDaughters()[i_dau]->ext<MCPpyjet>();
2480  if ( k[j_dau][1] < 30 ) {
2481  k[j_dau][1]=0;
2482  }
2483  }
2484 
2485  }
2486  }
2487  }
2488  for ( int i_py=1 ; i_py<=nlund ; i_py++){
2489  ijet=abs(jet[i_py]);
2490 
2491  if ( k[i_py][1]%30 == 1 ) { // true quanaties for this jet (all true stable) NB. Checks HEPEVT code,
2492  // not wether there are daughters. This is right, since Mokka/DDsim does not
2493  // change the HEPEVT status in MCParticle, even if an interaction with
2494  // seeable daughters takes place. To check: "un-used daughters", ie.
2495  // the opposite case, where the generator has decayed a particle,
2496  // but the decay happens after Mokka/DDsim has destroyed the parent.
2497  // Also: generator-decays of longlived charged-particles: probaly
2498  // (chek that) the MCParticle info is after application of the
2499  // B-field, so the direction of the momentum has changed.
2500  // [Finally: the crossing-angle. To what has it been applied?
2501  // everything it seems, so that should be OK]
2502 
2503  tmom[ijet][0]+=p[i_py][1] ;
2504  tmom[ijet][1]+=p[i_py][2] ;
2505  tmom[ijet][2]+=p[i_py][3] ;
2506  tE[ijet]+=p[i_py][4] ;
2507 
2508  }
2509  }
2510  //!! Add info on jet-type to return (1=from string,2=lepton,3=from cluster,4=isr,5=overlay,6=M.E. photon.
2511  //!! (+ x00=comes from boson, boson from jet x) and indication of the companion.
2512 
2513 
2514  for ( int k_jet=1 ; k_jet <= njet ; k_jet++ ) {
2515  if ( k_jet <= 2*nstr ){
2516  type[k_jet]=1;
2517  } else if ( k_jet<= 2*nstr+2*nclu ){
2518  type[k_jet]=3;
2519  } else if ( k_jet<= 2*nstr+2*nclu+n_hard_lepton ){
2520  type[k_jet]=2;
2521  } else if ( k_jet<= 2*nstr+2*nclu+n_hard_lepton+nphot ){
2522  type[k_jet]=6;
2523  } else if ( k_jet<= 2*nstr+n_hard_lepton+2*nclu+nphot+nisr ){
2524  type[k_jet]=4;
2525  } else if ( k_jet<= njet ){
2526  type[k_jet]=5;
2527  }
2528  if ( boson[k_jet] != 0 ) {
2529  type[k_jet]=abs(jet[boson[k_jet]])*100+type[k_jet];
2530  }
2531  if ( fafp_boson[k_jet] != 0 ) {
2532  fafp[k_jet]=fafp_boson[k_jet];
2533  } else {
2534  fafp[k_jet]=fafp_last[k_jet];
2535  }
2536  }
2537 
2538  for (int k_jet=1 ; k_jet<=njet-nisr-n_beam_jet-nphot-n_hard_lepton%2 ; k_jet++) {
2539  if ( k_jet%2 == 0 ) {
2540  companion[k_jet]=k_jet-1;
2541  } else {
2542  companion[k_jet]=k_jet+1;
2543  }
2544  }
2545 
2546 
2547  // TODO: the following would be much simpler if the groups would be sets, not lists: then the
2548  // order of the entries would not be important, and one would not need to explictly add
2549  // k_jet to it's own group. In fact, except for top-events , any jets is always grouped with
2550  // it's companion, so all what the following code is doing is to check is the initial c.n.
2551  // jets connect to a jet which is not already it's companion...
2552  // [In top events, where the initial b:s from the t decay are each-others companions, but should
2553  // not be grouped - rather we want to group jets from the same top together into an "initial
2554  // colour-neutral" - a misnomer in this case.]
2555 
2556  for ( int k_jet=1 ; k_jet<= njet ; k_jet++ ) {
2557  group[0][k_jet]=0;
2558  }
2559 
2560  for ( int k_jet=1 ; k_jet<= njet ; k_jet++ ) {
2561  if ( type[k_jet] == 4 || type[k_jet] == 6 ) { // trivial grouping of photons
2562  group[0][k_jet]=1 ; group[1][k_jet]=k_jet ;
2563  }
2564 
2565  if ( type[k_jet]%100 <= 3 ) { // only group strings, clusters and hard leptons
2566 
2567  bool k_lept=(abs(k[fafp[k_jet]][2])>=11 && abs(k[fafp[k_jet]][2])<=16);
2568 
2569  streamlog_out(DEBUG4) << std::endl;
2570 
2571  for ( int j_jet=1 ; j_jet<= njet ; j_jet++ ) {
2572  if ( type[j_jet]%100 <= 3 ) {
2573 
2574  bool j_lept=(abs(k[fafp[j_jet]][2])>=11 && abs(k[fafp[j_jet]][2])<=16);
2575 
2576  // the reasons to group k_jet with j_jet:
2577 
2578  bool colour_connected = (( k[fafp[j_jet]][2]<0 && k[fafp[k_jet]][2]>0 &&
2579  k[fafp[j_jet]][8] == k[fafp[k_jet]][7] && k[fafp[j_jet]][8] != 0 ) ||
2580  ( k[fafp[j_jet]][2]>0 && k[fafp[k_jet]][2]<0 &&
2581  k[fafp[j_jet]][7] == k[fafp[k_jet]][8] && k[fafp[j_jet]][7] != 0 ));
2582  bool no_colour = ( k[fafp[j_jet]][7]+ k[fafp[j_jet]][8]+ k[fafp[k_jet]][7]+ k[fafp[k_jet]][8] == 0 );
2583 
2584  bool boson_connected = (k[fafp[j_jet]][3]==k[fafp[k_jet]][3] &&
2585  ( abs(k[k[fafp[k_jet]][3]][2]) == 6 ||
2586  k[k[fafp[k_jet]][3]][2] == 23 ||
2587  abs(k[k[fafp[k_jet]][3]][2]) == 24 ||
2588  k[k[fafp[k_jet]][3]][2] == 25 ));
2589  bool same_parent = ( k[fafp[j_jet]][3]==k[fafp[k_jet]][3] && k[fafp[k_jet]][3] != 0 );
2590  bool same_family = ( flavour( k[fafp[k_jet]][2]) == flavour( k[fafp[j_jet]][2])) ;
2591 
2592  if (k_lept == j_lept || _top_event ) { // only group quarks with quarks, leptons with leptons, except in top events,
2593  // where we find it more useful to group all decay-products of the tops together,
2594  // which might imply grouping quarks with leptons. Note that the concept of
2595  // "initial colour neutral" is therefore not really apt in this case!
2596 
2597  if (( fafp[j_jet] == fafp[k_jet] ) || // the fafp:s are the same or have the same parent. NB will group k_jet with itself!
2598  colour_connected ||
2599  boson_connected ||
2600  (no_colour && (same_parent && same_family)) ||
2601  (companion[k_jet] == j_jet && abs( k[fafp[k_jet]][2]) != 6) ) {
2602 
2603  group[0][k_jet]++;
2604  group[ group[0][k_jet] ][k_jet] = j_jet;
2605 
2606  if ( k_jet != j_jet ) {
2607  streamlog_out(DEBUG4) << " Jet " << k_jet << " was grouped with jet " << j_jet <<" .";
2608  streamlog_out(DEBUG4) << " fafp of the jets are " << fafp[k_jet] << " (pdg " << k[fafp[k_jet]][2]
2609  << ") and " << fafp[j_jet] <<" (pdg " << k[fafp[j_jet]][2] <<" ) ." ;
2610  if ( k[fafp[k_jet]][3] != 0 ) {
2611  streamlog_out(DEBUG4) << " parents of the fafp of the jets are " << k[fafp[k_jet]][3]<< " and " << k[fafp[j_jet]][3] <<" .";
2612  }
2613  streamlog_out(DEBUG4) << " colour-connection : " << colour_connected<< std::endl;
2614  }
2615  }
2616  }
2617  }
2618  }
2619  streamlog_out(DEBUG4) << " group of jet " << k_jet << " : " ;
2620  bool compthere=false;
2621  for ( int i_grp=1 ; i_grp <= group[0][k_jet] ; i_grp++ ) {
2622  streamlog_out(DEBUG4) <<group[i_grp][k_jet] << " " ;
2623  if (group[i_grp][k_jet] == companion[k_jet] || companion[k_jet] == 0 ) { compthere = true ; }
2624  }
2625  if ( ! compthere && not _top_event ) {
2626  streamlog_out(ERROR) << " jet " << k_jet << " : companion " << companion[k_jet] << " not in group " << std::endl;
2627  }
2628  streamlog_out(DEBUG4) << std::endl;
2629  if ( group[0][k_jet] <= 1 && ( type[k_jet] != 2 || n_hard_lepton%2 == 0 ) ) {
2630  streamlog_out(ERROR) << " jet " << k_jet << " ungrouped " << std::endl;
2631  }
2632  }
2633  }
2634 
2635  n_djb = 0 ;
2636  n_dje = 0;
2637 
2638  //! figure out the different jet combinations. (If there are no bosons, it's straight forward:
2639  //! it's just the odd+even combinations. Also if one only cares about the final singlet
2640  //! the same aplies, boson or not)
2641 
2642  // We want to find how many colour-neutrals we have before the P.S., and which jets they
2643  // go into. dijet_begining[k_jet] will be the number of that c.n., and conversely, jets_begin[initial cn#]
2644  // is the list of the jets each initial cn goes to.
2645  // Likewise (but trivially) dijet_end[k_jet] is the number of the final cn the jet, belongs to -
2646  // jets_end[final cn#] is the list (which always, by construction, contains two adjecent jets).
2647 
2648  // TODO: the below would be simpler if group would be sets. Then just set-equality would be needed to check
2649  // The logic below needs identical groups not only to contain the same numbers, but also that they are in the same
2650  // order!
2651 
2652  for ( int k_jet=1; k_jet<=njet ; k_jet++ ) {
2653 
2654  dijet_begining[k_jet] = 0;
2655 
2656  if ( group[0][k_jet] > 0 ) {
2657  for ( int i_djb=1; i_djb<=n_djb ; i_djb++ ) {
2658  int gotit=1;
2659  for ( int i_grpmbr=1 ; i_grpmbr<= group[0][k_jet] ; i_grpmbr++ ) {
2660  if ( ! (group[i_grpmbr][k_jet]==jets_begin[i_grpmbr][i_djb] )) {
2661  gotit=0;
2662  }
2663  }
2664  if ( gotit ) {
2665  dijet_begining[k_jet] = i_djb;
2666  }
2667  }
2668  }
2669 
2670  if ( dijet_begining[k_jet] == 0 && group[0][k_jet] > 0 ) {
2671  n_djb++;
2672  dijet_begining[k_jet] = n_djb;
2673  for ( int i_grpmbr=0 ; i_grpmbr<= group[0][k_jet] ; i_grpmbr++ ){
2674  jets_begin[i_grpmbr][n_djb] = group[i_grpmbr][k_jet];
2675  }
2676  }
2677 
2678  dijet_end[k_jet] = 0;
2679  for ( int i_dje=1 ; i_dje<=n_dje ; i_dje++ ) {
2680  if ( std::min(k_jet,companion[k_jet])==jets_end[1][i_dje] &&
2681  std::max(k_jet,companion[k_jet])==jets_end[2][i_dje] ) {
2682  dijet_end[k_jet] = i_dje;
2683  }
2684  }
2685 
2686  if ( dijet_end[k_jet] == 0 ) {
2687  n_dje++;
2688  dijet_end[k_jet] = n_dje;
2689  if ( companion[k_jet] > 0 ) {
2690  jets_end[1][n_dje] = std::min(k_jet,companion[k_jet]);
2691  jets_end[2][n_dje] = std::max(k_jet,companion[k_jet]);
2692  jets_end[0][n_dje] = 2 ;
2693  } else {
2694  jets_end[1][n_dje] = k_jet;
2695  jets_end[0][n_dje] = 1;
2696  }
2697  }
2698 
2699  }
2700 
2701 
2702  // collect some topology and consitency information:
2703 
2704  nboson=0;n_jetless=0; n_0_E_jets=0; n_2jet_clu=0;
2705  for ( int k_jet=1 ; k_jet<=njet ; k_jet++ ) {
2706  if (type[k_jet]/100 != 0 ){
2707  nboson++;
2708  }
2709  if (tE[k_jet]<0.000001 ) {
2710  n_0_E_jets++;
2711  }
2712  if ( jet[k_jet]==0 && k[k_jet][1]==1 ) {
2713  n_jetless++;
2714  }
2715  if ( (type[k_jet]%100== 3) && ( tE[k_jet]> 0.0000001 ) &&
2716  ( tE[companion[k_jet]] > 0.0000001)) {
2717  n_2jet_clu++;
2718  }
2719  }
2720 
2721  n_mixed=0;
2722  for ( int k_jet=1 ; k_jet<=njet-1 ; k_jet+=2 ) {
2723  if ((type[k_jet]%100==1) &&
2724  (k[fafp[k_jet]][4] != k[fafp[k_jet+1]][4]) &&
2725  ( !( ( k[fafp[k_jet]][2] > 0 && k[fafp[k_jet]][7] == k[fafp[k_jet+1]][8]) ||
2726  ( k[fafp[k_jet]][2] < 0 && k[fafp[k_jet]][8] == k[fafp[k_jet+1]][7] )) ) &&
2727  (k[fafp[k_jet]][3] != k[fafp[k_jet+1]][3]) ) {
2728 
2729  int ii = fafp[k_jet];
2730  while ( ii <= nlund ) {
2731  if ( ii == elementon[k_jet] ) {
2732  break;
2733  } else {
2734  ii = k[ii][4] ;
2735  if ( ii > elementon[k_jet] ) {
2736  n_mixed++;
2737  break ;
2738  }
2739  }
2740  }
2741  }
2742  }
2743  // jsum=jets_summary_t(njet,n_djb,n_dje,2*nstr, n_hard_lepton, nboson, 2*nclu,nisr , n_mixed, n_jetless, n_2jet_clu, n_0_E_jets,jt(1:njet),dijb(1:n_djb),dije(1:n_dje))
2744 
2745 }
2746 
2747 void TrueJet::check( LCEvent * /*evt*/ ) {
2748  // nothing to check here - could be used to fill checkplots in reconstruction processor
2749 }
2750 
2751 
2752 void TrueJet::end(){
2753 
2754  // std::cout << "TrueJet::end() " << name()
2755  // << " processed " << _nEvt << " events in " << _nRun << " runs "
2756  // << std::endl ;
2757 
2758 }
TrueJet::_recoParticleCollectionName
std::string _recoParticleCollectionName
Definition: TrueJet.h:79
TrueJet::grouping
void grouping()
Definition: TrueJet.cc:2440
TrueJet::_nRun
int _nRun
Definition: TrueJet.h:91
TrueJet::evt
LCEvent * evt
Definition: TrueJet.h:108
TrueJet::flavour
int flavour(int k2)
Definition: TrueJet.cc:2424
TrueJet::_finalColourNeutralCollectionName
std::string _finalColourNeutralCollectionName
Definition: TrueJet.h:83
TrueJet::fafp
int fafp[26]
Definition: TrueJet.h:190
TrueJet
Example processor for marlin.
Definition: TrueJet.h:38
TrueJet::n_jetless
int n_jetless
Definition: TrueJet.h:173
TrueJet::n_0_E_jets
int n_0_E_jets
Definition: TrueJet.h:175
TrueJet::init
virtual void init()
Called at the begin of the job before anything is read.
Definition: TrueJet.cc:124
TrueJet::p
double p[4001][6]
Definition: TrueJet.h:118
TrueJet::end
virtual void end()
Called after data processing for clean up.
Definition: TrueJet.cc:2752
TrueJet::processRunHeader
virtual void processRunHeader(LCRunHeader *run)
Called for every run.
Definition: TrueJet.cc:137
TrueJet::jet
int jet[4001]
Definition: TrueJet.h:151
TrueJet::_whiz1
bool _whiz1
Definition: TrueJet.h:206
TrueJet::_top_event
bool _top_event
Definition: TrueJet.h:207
TrueJet::_trueJetPFOLink
std::string _trueJetPFOLink
Definition: TrueJet.h:85
TrueJet::dijet_begining
int dijet_begining[26]
Definition: TrueJet.h:194
TrueJet::boson
int boson[26]
Definition: TrueJet.h:188
TrueJet::companion
int companion[4001]
Definition: TrueJet.h:152
TrueJet::nisr
int nisr
Definition: TrueJet.h:170
TrueJet::_finalColourNeutralLink
std::string _finalColourNeutralLink
Definition: TrueJet.h:89
TrueJet::cluster
void cluster()
Definition: TrueJet.cc:1913
TrueJet::k
int k[4001][9]
Definition: TrueJet.h:119
TrueJet::njet
int njet
Definition: TrueJet.h:160
TrueJet::dijet_end
int dijet_end[26]
Definition: TrueJet.h:196
TrueJet::elementon
int elementon[26]
Definition: TrueJet.h:184
TrueJet::nfsr
int nfsr[26]
Definition: TrueJet.h:191
TrueJet::_finalElementonLink
std::string _finalElementonLink
Definition: TrueJet.h:87
TrueJet::processEvent
virtual void processEvent(LCEvent *evt)
Called for every event - the working horse.
Definition: TrueJet.cc:144
TrueJet::nstr
int nstr
Definition: TrueJet.h:164
reltrue
LCRelationNavigator * reltrue
Definition: TrueJet.cc:23
TrueJet::check
virtual void check(LCEvent *evt)
Definition: TrueJet.cc:2747
TrueJet::n_dje
int n_dje
Definition: TrueJet.h:162
TrueJet::n_2jet_clu
int n_2jet_clu
Definition: TrueJet.h:174
TrueJet::nphot
int nphot
Definition: TrueJet.h:171
TrueJet::_recoMCTruthLink
std::string _recoMCTruthLink
Definition: TrueJet.h:80
TrueJet::photon
void photon()
Definition: TrueJet.cc:1957
TrueJet::stdhep_reader_bug_workaround
void stdhep_reader_bug_workaround(int line94)
Definition: TrueJet.cc:1540
TrueJet::_initialColourNeutralLink
std::string _initialColourNeutralLink
Definition: TrueJet.h:90
TrueJet::_initialElementonLink
std::string _initialElementonLink
Definition: TrueJet.h:88
TrueJet::isr
void isr()
Definition: TrueJet.cc:2007
TrueJet::type
int type[26]
Definition: TrueJet.h:192
TrueJet::_trueJetCollectionName
std::string _trueJetCollectionName
ouput collection names.
Definition: TrueJet.h:82
TrueJet::group
int group[26][26]
Definition: TrueJet.h:193
TrueJet::current_jet
int current_jet
Definition: TrueJet.h:155
TrueJet::true_lepton
void true_lepton()
Definition: TrueJet.cc:1641
TrueJet::n_hard_lepton
int n_hard_lepton
Definition: TrueJet.h:165
sistrip::e
static const float e
Definition: PhysicalConstants.h:26
TrueJet::_higgs_to_glue_glue
bool _higgs_to_glue_glue
Definition: TrueJet.h:208
TrueJet::nlund
int nlund
Definition: TrueJet.h:146
TrueJet::string
void string()
Definition: TrueJet.cc:2046
TrueJet::_trueJetMCParticleLink
std::string _trueJetMCParticleLink
Definition: TrueJet.h:86
TrueJet::jets_end
int jets_end[26][26]
Definition: TrueJet.h:204
TrueJet::getPyjets
void getPyjets(LCCollection *mcpcol)
Definition: TrueJet.cc:1192
sistrip::LCCollectionVec
std::vector< LCCollection * > LCCollectionVec
Definition: SiStripClus.h:55
TrueJet::_initialColourNeutralCollectionName
std::string _initialColourNeutralCollectionName
Definition: TrueJet.h:84
TrueJet::mcp_pyjets
MCParticleVec mcp_pyjets
Definition: TrueJet.h:147
TrueJet::tE
double tE[26]
Definition: TrueJet.h:198
TrueJet::n_djb
int n_djb
Definition: TrueJet.h:161
TrueJet::_MCParticleColllectionName
std::string _MCParticleColllectionName
Input collection name.
Definition: TrueJet.h:78
TrueJet::_nEvt
int _nEvt
Definition: TrueJet.h:92
TrueJet::first_line
int first_line
Definition: TrueJet.h:156
TrueJet::jets_begin
int jets_begin[26][26]
Definition: TrueJet.h:201
TrueJet::first_parton
void first_parton(int this_partic, int this_jet, int &first_partic, int &last_94_parent, int &nfsr, int &info, int &info2)
Definition: TrueJet.cc:2280
TrueJet::nclu
int nclu
Definition: TrueJet.h:167
TrueJet::fafp_boson
int fafp_boson[26]
Definition: TrueJet.h:189
TrueJet::nboson
int nboson
Definition: TrueJet.h:166
TrueJet::TrueJet
TrueJet()
Definition: TrueJet.cc:27
TrueJet::fafp_last
int fafp_last[26]
Definition: TrueJet.h:182
TrueJet::n_beam_jet
int n_beam_jet
Definition: TrueJet.h:176
TrueJet::n_mixed
int n_mixed
Definition: TrueJet.h:172
aTrueJet
TrueJet aTrueJet
Definition: TrueJet.cc:25
TrueJet::assign_jet
void assign_jet(int jet1, int jet2, int this_fafp)
Definition: TrueJet.cc:2202
TrueJet::tmom
double tmom[26][3]
Definition: TrueJet.h:197