diff --git a/Hadronization/ClusterHadronizationHandler.cc b/Hadronization/ClusterHadronizationHandler.cc --- a/Hadronization/ClusterHadronizationHandler.cc +++ b/Hadronization/ClusterHadronizationHandler.cc @@ -1,307 +1,320 @@ // -*- C++ -*- // // ClusterHadronizationHandler.cc is a part of Herwig - A multi-purpose Monte Carlo event generator // Copyright (C) 2002-2017 The Herwig Collaboration // // Herwig is licenced under version 3 of the GPL, see COPYING for details. // Please respect the MCnet academic guidelines, see GUIDELINES for details. // // // This is the implementation of the non-inlined, non-templated member // functions of the ClusterHadronizationHandler class. // #include "ClusterHadronizationHandler.h" #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include "Herwig/Utilities/EnumParticles.h" #include "CluHadConfig.h" #include "Cluster.h" #include using namespace Herwig; ClusterHadronizationHandler * ClusterHadronizationHandler::currentHandler_ = 0; DescribeClass describeClusterHadronizationHandler("Herwig::ClusterHadronizationHandler",""); IBPtr ClusterHadronizationHandler::clone() const { return new_ptr(*this); } IBPtr ClusterHadronizationHandler::fullclone() const { return new_ptr(*this); } void ClusterHadronizationHandler::persistentOutput(PersistentOStream & os) const { os << _partonSplitter << _clusterFinder << _colourReconnector << _clusterFissioner << _lightClusterDecayer << _clusterDecayer << ounit(_minVirtuality2,GeV2) << ounit(_maxDisplacement,mm) << _underlyingEventHandler << _reduceToTwoComponents; } void ClusterHadronizationHandler::persistentInput(PersistentIStream & is, int) { is >> _partonSplitter >> _clusterFinder >> _colourReconnector >> _clusterFissioner >> _lightClusterDecayer >> _clusterDecayer >> iunit(_minVirtuality2,GeV2) >> iunit(_maxDisplacement,mm) >> _underlyingEventHandler >> _reduceToTwoComponents; } void ClusterHadronizationHandler::Init() { static ClassDocumentation documentation ("This is the main handler class for the Cluster Hadronization", "The hadronization was performed using the cluster model of \\cite{Webber:1983if}.", "%\\cite{Webber:1983if}\n" "\\bibitem{Webber:1983if}\n" " B.~R.~Webber,\n" " ``A QCD Model For Jet Fragmentation Including Soft Gluon Interference,''\n" " Nucl.\\ Phys.\\ B {\\bf 238}, 492 (1984).\n" " %%CITATION = NUPHA,B238,492;%%\n" // main manual ); static Reference interfacePartonSplitter("PartonSplitter", "A reference to the PartonSplitter object", &Herwig::ClusterHadronizationHandler::_partonSplitter, false, false, true, false); static Reference interfaceClusterFinder("ClusterFinder", "A reference to the ClusterFinder object", &Herwig::ClusterHadronizationHandler::_clusterFinder, false, false, true, false); static Reference interfaceColourReconnector("ColourReconnector", "A reference to the ColourReconnector object", &Herwig::ClusterHadronizationHandler::_colourReconnector, false, false, true, false); static Reference interfaceClusterFissioner("ClusterFissioner", "A reference to the ClusterFissioner object", &Herwig::ClusterHadronizationHandler::_clusterFissioner, false, false, true, false); static Reference interfaceLightClusterDecayer("LightClusterDecayer", "A reference to the LightClusterDecayer object", &Herwig::ClusterHadronizationHandler::_lightClusterDecayer, false, false, true, false); static Reference interfaceClusterDecayer("ClusterDecayer", "A reference to the ClusterDecayer object", &Herwig::ClusterHadronizationHandler::_clusterDecayer, false, false, true, false); static Parameter interfaceMinVirtuality2 ("MinVirtuality2", "Minimum virtuality^2 of partons to use in calculating distances (unit [GeV2]).", &ClusterHadronizationHandler::_minVirtuality2, GeV2, 0.1*GeV2, ZERO, 10.0*GeV2,false,false,false); static Parameter interfaceMaxDisplacement ("MaxDisplacement", "Maximum displacement that is allowed for a particle (unit [millimeter]).", &ClusterHadronizationHandler::_maxDisplacement, mm, 1.0e-10*mm, 0.0*mm, 1.0e-9*mm,false,false,false); static Reference interfaceUnderlyingEventHandler ("UnderlyingEventHandler", "Pointer to the handler for the Underlying Event. " "Set to NULL to disable.", &ClusterHadronizationHandler::_underlyingEventHandler, false, false, true, true, false); static Switch interfaceReduceToTwoComponents ("ReduceToTwoComponents", "Whether or not to reduce three component baryon-number violating clusters to two components before cluster splitting or leave" " this till after the cluster splitting", &ClusterHadronizationHandler::_reduceToTwoComponents, true, false, false); static SwitchOption interfaceReduceToTwoComponentsYes (interfaceReduceToTwoComponents, "BeforeSplitting", "Reduce to two components", true); static SwitchOption interfaceReduceToTwoComponentsNo (interfaceReduceToTwoComponents, "AfterSplitting", "Treat as three components", false); } namespace { void extractChildren(tPPtr p, set & all) { if (p->children().empty()) return; for (PVector::const_iterator child = p->children().begin(); child != p->children().end(); ++child) { all.insert(*child); extractChildren(*child, all); } } } void ClusterHadronizationHandler:: handle(EventHandler & ch, const tPVector & tagged, const Hint &) { useMe(); currentHandler_ = this; PVector currentlist(tagged.begin(),tagged.end()); // set the scale for coloured particles to just above the gluon mass squared // if less than this so they are classed as perturbative Energy2 Q02 = 1.01*sqr(getParticleData(ParticleID::g)->constituentMass()); for(unsigned int ix=0;ixscale()scale(Q02); } // split the gluons _partonSplitter->split(currentlist); // form the clusters ClusterVector clusters = _clusterFinder->formClusters(currentlist); // reduce BV clusters to two components now if needed if(_reduceToTwoComponents) _clusterFinder->reduceToTwoComponents(clusters); // perform colour reconnection if needed and then // decay the clusters into one hadron bool lightOK = false; short tried = 0; const ClusterVector savedclusters = clusters; tPVector finalHadrons; // only needed for partonic decayer while (!lightOK && tried++ < 10) { // no colour reconnection with baryon-number-violating (BV) clusters ClusterVector CRclusters, BVclusters; CRclusters.reserve( clusters.size() ); BVclusters.reserve( clusters.size() ); for (size_t ic = 0; ic < clusters.size(); ++ic) { ClusterPtr cl = clusters.at(ic); bool hasClusterParent = false; for (unsigned int ix=0; ix < cl->parents().size(); ++ix) { if (cl->parents()[ix]->id() == ParticleID::Cluster) { hasClusterParent = true; break; } } if (cl->numComponents() > 2 || hasClusterParent) BVclusters.push_back(cl); else CRclusters.push_back(cl); } // colour reconnection _colourReconnector->rearrange(CRclusters); + // tag new clusters as children of the partons to hadronize _setChildren(CRclusters); - + + // forms diquarks _clusterFinder->reduceToTwoComponents(CRclusters); - + // recombine vectors of (possibly) reconnected and BV clusters clusters.clear(); clusters.insert( clusters.end(), CRclusters.begin(), CRclusters.end() ); clusters.insert( clusters.end(), BVclusters.begin(), BVclusters.end() ); // fission of heavy clusters // NB: during cluster fission, light hadrons might be produced straight away finalHadrons = _clusterFissioner->fission(clusters,isSoftUnderlyingEventON()); + // if clusters not previously reduced to two components do it now if(!_reduceToTwoComponents) _clusterFinder->reduceToTwoComponents(clusters); lightOK = _lightClusterDecayer->decay(clusters,finalHadrons); // if the decay of the light clusters was not successful, undo the cluster // fission and decay steps and revert to the original state of the event // record if (!lightOK) { clusters = savedclusters; for_each(clusters.begin(), clusters.end(), std::mem_fn(&Particle::undecay)); } } if (!lightOK) { - throw Exception("CluHad::handle(): tried LightClusterDecayer 10 times!", + throw Exception( "CluHad::handle(): tried LightClusterDecayer 10 times!", Exception::eventerror); } - // decay the remaining clusters _clusterDecayer->decay(clusters,finalHadrons); // ***************************************** // ***************************************** // ***************************************** + bool finalStateCluster=false; StepPtr pstep = newStep(); set allDecendants; for (tPVector::const_iterator it = tagged.begin(); it != tagged.end(); ++it) { extractChildren(*it, allDecendants); } for(set::const_iterator it = allDecendants.begin(); it != allDecendants.end(); ++it) { // this is a workaround because the set sometimes // re-orders parents after their children - if ((*it)->children().empty()) + if ((*it)->children().empty()){ + // If there is a cluster in the final state throw an event error + if((*it)->id()==81) { + finalStateCluster=true; + } pstep->addDecayProduct(*it); + } else { pstep->addDecayProduct(*it); pstep->addIntermediate(*it); } } + // For very small center of mass energies it might happen that baryonic clusters cannot decay into hadrons + if (finalStateCluster){ + throw Exception( "CluHad::Handle(): Cluster in the final state", + Exception::eventerror); + } // ***************************************** // ***************************************** // ***************************************** // soft underlying event if needed if (isSoftUnderlyingEventON()) { assert(_underlyingEventHandler); ch.performStep(_underlyingEventHandler,Hint::Default()); } } // Sets parent child relationship of all clusters with two components // Relationships for clusters with more than two components are set elsewhere in the Colour Reconnector void ClusterHadronizationHandler::_setChildren(const ClusterVector & clusters) const { // erase existing information about the partons' children tPVector partons; for ( const auto & cl : clusters ) { if ( cl->numComponents() > 2 ) continue; partons.push_back( cl->colParticle() ); partons.push_back( cl->antiColParticle() ); } // erase all previous information about parent child relationship for_each(partons.begin(), partons.end(), std::mem_fn(&Particle::undecay)); // give new parents to the clusters: their constituents for ( const auto & cl : clusters ) { if ( cl->numComponents() > 2 ) continue; cl->colParticle()->addChild(cl); cl->antiColParticle()->addChild(cl); } } diff --git a/Hadronization/ColourReconnector.cc b/Hadronization/ColourReconnector.cc --- a/Hadronization/ColourReconnector.cc +++ b/Hadronization/ColourReconnector.cc @@ -1,823 +1,826 @@ // -*- C++ -*- // // ColourReconnector.cc is a part of Herwig - A multi-purpose Monte Carlo event generator // Copyright (C) 2002-2017 The Herwig Collaboration // // Herwig is licenced under version 3 of the GPL, see COPYING for details. // Please respect the MCnet academic guidelines, see GUIDELINES for details. // // // This is the implementation of the non-inlined, non-templated member // functions of the ColourReconnector class. // #include "ColourReconnector.h" #include "Cluster.h" #include #include #include #include #include #include #include #include "Herwig/Utilities/Maths.h" using namespace Herwig; using CluVecIt = ColourReconnector::CluVecIt; using Constants::pi; using Constants::twopi; DescribeClass describeColourReconnector("Herwig::ColourReconnector",""); IBPtr ColourReconnector::clone() const { return new_ptr(*this); } IBPtr ColourReconnector::fullclone() const { return new_ptr(*this); } void ColourReconnector::rearrange(ClusterVector & clusters) { if (_clreco == 0) return; // need at least two clusters if (clusters.size() < 2) return; // do the colour reconnection switch (_algorithm) { case 0: _doRecoPlain(clusters); break; case 1: _doRecoStatistical(clusters); break; case 2: _doRecoBaryonic(clusters); break; } } Energy2 ColourReconnector::_clusterMassSum(const PVector & q, const PVector & aq) const { const size_t nclusters = q.size(); assert (aq.size() == nclusters); Energy2 sum = ZERO; for (size_t i = 0; i < nclusters; i++) sum += ( q[i]->momentum() + aq[i]->momentum() ).m2(); return sum; } bool ColourReconnector::_containsColour8(const ClusterVector & cv, const vector & P) const { assert (P.size() == cv.size()); for (size_t i = 0; i < cv.size(); i++) { tcPPtr p = cv[i]->colParticle(); tcPPtr q = cv[P[i]]->antiColParticle(); if (_isColour8(p, q)) return true; } return false; } void ColourReconnector::_doRecoStatistical(ClusterVector & cv) const { const size_t nclusters = cv.size(); // initially, enumerate (anti)quarks as given in the cluster vector ParticleVector q, aq; for (size_t i = 0; i < nclusters; i++) { q.push_back( cv[i]->colParticle() ); aq.push_back( cv[i]->antiColParticle() ); } // annealing scheme Energy2 t, delta; Energy2 lambda = _clusterMassSum(q,aq); const unsigned _ntries = _triesPerStepFactor * nclusters; // find appropriate starting temperature by measuring the largest lambda // difference in some dry-run random rearrangements { vector typical; for (int i = 0; i < 10; i++) { const pair toswap = _shuffle(q,aq,5); ParticleVector newaq = aq; swap (newaq[toswap.first], newaq[toswap.second]); Energy2 newlambda = _clusterMassSum(q,newaq); typical.push_back( abs(newlambda - lambda) ); } t = _initTemp * Math::median(typical); } // anneal in up to _annealingSteps temperature steps for (unsigned step = 0; step < _annealingSteps; step++) { // For this temperature step, try to reconnect _ntries times. Stop the // algorithm if no successful reconnection happens. unsigned nSuccess = 0; for (unsigned it = 0; it < _ntries; it++) { // make a random rearrangement const unsigned maxtries = 10; const pair toswap = _shuffle(q,aq,maxtries); const int i = toswap.first; const int j = toswap.second; // stop here if we cannot find any allowed reconfiguration if (i == -1) break; // create a new antiquark vector with the two partons swapped ParticleVector newaq = aq; swap (newaq[i], newaq[j]); // Check if lambda would decrease. If yes, accept the reconnection. If no, // accept it only with a probability given by the current Boltzmann // factor. In the latter case we set p = 0 if the temperature is close to // 0, to avoid division by 0. Energy2 newlambda = _clusterMassSum(q,newaq); delta = newlambda - lambda; double prob = 1.0; if (delta > ZERO) prob = ( abs(t) < 1e-8*MeV2 ) ? 0.0 : exp(-delta/t); if (UseRandom::rnd() < prob) { lambda = newlambda; swap (newaq, aq); nSuccess++; } } if (nSuccess == 0) break; // reduce temperature t *= _annealingFactor; } // construct the new cluster vector ClusterVector newclusters; for (size_t i = 0; i < nclusters; i++) { ClusterPtr cl = new_ptr( Cluster( q[i], aq[i] ) ); newclusters.push_back(cl); } swap(newclusters,cv); return; } void ColourReconnector::_doRecoPlain(ClusterVector & cv) const { ClusterVector newcv = cv; // try to avoid systematic errors by randomising the reconnection order long (*p_irnd)(long) = UseRandom::irnd; random_shuffle( newcv.begin(), newcv.end(), p_irnd ); // iterate over all clusters for (CluVecIt cit = newcv.begin(); cit != newcv.end(); cit++) { // find the cluster which, if reconnected with *cit, would result in the // smallest sum of cluster masses // NB this method returns *cit if no reconnection partner can be found CluVecIt candidate = _findRecoPartner(cit, newcv); // skip this cluster if no possible reshuffling partner can be found if (candidate == cit) continue; // accept the reconnection with probability _preco. if (UseRandom::rnd() < _preco) { pair reconnected = _reconnect(*cit, *candidate); // Replace the clusters in the ClusterVector. The order of the // colour-triplet partons in the cluster vector is retained here. // replace *cit by reconnected.first *cit = reconnected.first; // replace candidate by reconnected.second *candidate = reconnected.second; } } swap(cv,newcv); return; } namespace { inline bool hasDiquark(CluVecIt cit) { for(int i = 0; i<(*cit)->numComponents(); i++) { if (DiquarkMatcher::Check(*((*cit)->particle(i)->dataPtr()))) return true; } return false; } } // Implementation of the baryonic reconnection algorithm void ColourReconnector::_doRecoBaryonic(ClusterVector & cv) const { ClusterVector newcv = cv; ClusterVector deleted; deleted.reserve(cv.size()); // try to avoid systematic errors by randomising the reconnection order long (*p_irnd)(long) = UseRandom::irnd; random_shuffle( newcv.begin(), newcv.end(), p_irnd ); // iterate over all clusters for (CluVecIt cit = newcv.begin(); cit != newcv.end(); ++cit) { //avoid clusters already containing diuarks if (hasDiquark(cit)) continue; //skip the cluster to be deleted later 3->2 cluster if (find(deleted.begin(), deleted.end(), *cit) != deleted.end()) continue; // Skip all found baryonic clusters, this biases the algorithm but implementing // something like re-reconnection is ongoing work if ((*cit)->numComponents()==3) continue; // Find a candidate suitable for reconnection CluVecIt baryonic1, baryonic2; bool isBaryonicCandidate = false; CluVecIt candidate = _findPartnerBaryonic(cit, newcv, isBaryonicCandidate, deleted, baryonic1, baryonic2); // skip this cluster if no possible reconnection partner can be found if ( !isBaryonicCandidate && candidate==cit ) continue; if ( isBaryonicCandidate && UseRandom::rnd() < _precoBaryonic ) { deleted.push_back(*baryonic2); // Function that does the reconnection from 3 -> 2 clusters ClusterPtr b1, b2; _makeBaryonicClusters(*cit,*baryonic1,*baryonic2, b1, b2); *cit = b1; *baryonic1 = b2; // Baryonic2 is easily skipped in the next loop } // Normal 2->2 Colour reconnection if ( !isBaryonicCandidate && UseRandom::rnd() < _preco ) { auto reconnected = _reconnectBaryonic(*cit, *candidate); *cit = reconnected.first; *candidate = reconnected.second; } } // create a new vector of clusters except for the ones which are "deleted" during // baryonic reconnection ClusterVector clustervector; for ( const auto & cluster : newcv ) if ( find(deleted.begin(), deleted.end(), cluster) == deleted.end() ) clustervector.push_back(cluster); swap(cv,clustervector); + + + } namespace { double calculateRapidityRF(const Lorentz5Momentum & q1, const Lorentz5Momentum & p2) { //calculate rapidity wrt the direction of q1 //angle between the particles in the RF of cluster of q1 // calculate the z component of p2 w.r.t the direction of q1 if(q1.rho2()==ZERO) return 0.; const Energy pz = p2.vect() * q1.vect().unit(); if ( pz == ZERO ) return 0.; // Transverse momentum of p2 w.r.t the direction of q1 const Energy pt = sqrt(p2.vect().mag2() - sqr(pz)); // Transverse mass pf p2 w.r.t to the direction of q1 const Energy mtrans = sqrt(p2.mass()*p2.mass() + (pt*pt)); // Correct formula const double y2 = log((p2.t() + abs(pz))/mtrans); return ( pz < ZERO ) ? -y2 : y2; } } CluVecIt ColourReconnector::_findPartnerBaryonic( CluVecIt cl, ClusterVector & cv, bool & baryonicCand, const ClusterVector& deleted, CluVecIt &baryonic1, CluVecIt &baryonic2 ) const { using Constants::pi; using Constants::twopi; // Returns a candidate for possible reconnection CluVecIt candidate = cl; bool bcand = false; double maxrap = 0.0; double minrap = 0.0; double maxrapNormal = 0.0; double minrapNormal = 0.0; double maxsumnormal = 0.0; double maxsum = 0.0; double secondsum = 0.0; // boost into RF of cl Lorentz5Momentum cl1 = (*cl)->momentum(); const Boost boostv(-cl1.boostVector()); cl1.boost(boostv); // boost constituents of cl into RF of cl Lorentz5Momentum p1col = (*cl)->colParticle()->momentum(); Lorentz5Momentum p1anticol = (*cl)->antiColParticle()->momentum(); p1col.boost(boostv); p1anticol.boost(boostv); for (CluVecIt cit=cv.begin(); cit != cv.end(); ++cit) { //avoid looping over clusters containing diquarks if ( hasDiquark(cit) ) continue; if ( (*cit)->numComponents()==3 ) continue; if ( cit==cl ) continue; //skip the cluster to be deleted later 3->2 cluster if ( find(deleted.begin(), deleted.end(), *cit) != deleted.end() ) continue; if ( (*cl)->isBeamCluster() && (*cit)->isBeamCluster() ) continue; // stop it putting far apart clusters together if ( ( (**cl).vertex()-(**cit).vertex() ).m() >_maxDistance ) continue; const bool Colour8 = _isColour8( (*cl)->colParticle(), (*cit)->antiColParticle() ) || _isColour8( (*cit)->colParticle(), (*cl)->antiColParticle() ) ; if ( Colour8 ) continue; // boost constituents of cit into RF of cl Lorentz5Momentum p2col = (*cit)->colParticle()->momentum(); Lorentz5Momentum p2anticol = (*cit)->antiColParticle()->momentum(); p2col.boost(boostv); p2anticol.boost(boostv); // calculate the rapidity of the other constituents of the clusters // w.r.t axis of p1anticol.vect.unit const double rapq = calculateRapidityRF(p1anticol,p2col); const double rapqbar = calculateRapidityRF(p1anticol,p2anticol); // configuration for normal CR if ( rapq > 0.0 && rapqbar < 0.0 && rapq > maxrap && rapqbar < minrap ) { maxrap = rapq; minrap = rapqbar; //sum of rapidities of quarks const double normalsum = abs(rapq) + abs(rapqbar); if ( normalsum > maxsumnormal ) { maxsumnormal = normalsum; maxrapNormal = rapq; minrapNormal = rapqbar; bcand = false; candidate = cit; } } if ( rapq < 0.0 && rapqbar >0.0 && rapqbar > maxrapNormal && rapq < minrapNormal ) { maxrap = rapqbar; minrap = rapq; const double sumrap = abs(rapqbar) + abs(rapq); // first candidate gets here. If second baryonic candidate has higher Ysum than the first // one, the second candidate becomes the first one and the first the second. if (sumrap > maxsum) { if(maxsum != 0){ baryonic2 = baryonic1; baryonic1 = cit; bcand = true; } else { baryonic1 = cit; } maxsum = sumrap; } else { if (sumrap > secondsum && sumrap != maxsum) { secondsum = sumrap; bcand = true; baryonic2 = cit; } } } } if(bcand == true){ baryonicCand = true; } return candidate; } CluVecIt ColourReconnector::_findRecoPartner(CluVecIt cl, ClusterVector & cv) const { CluVecIt candidate = cl; Energy minMass = 1*TeV; for (CluVecIt cit=cv.begin(); cit != cv.end(); ++cit) { // don't even look at original cluster if(cit==cl) continue; // don't allow colour octet clusters if ( _isColour8( (*cl)->colParticle(), (*cit)->antiColParticle() ) || _isColour8( (*cit)->colParticle(), (*cl)->antiColParticle() ) ) { continue; } // stop it putting beam remnants together if((*cl)->isBeamCluster() && (*cit)->isBeamCluster()) continue; // stop it putting far apart clusters together if(((**cl).vertex()-(**cit).vertex()).m()>_maxDistance) continue; // momenta of the old clusters Lorentz5Momentum p1 = (*cl)->colParticle()->momentum() + (*cl)->antiColParticle()->momentum(); Lorentz5Momentum p2 = (*cit)->colParticle()->momentum() + (*cit)->antiColParticle()->momentum(); // momenta of the new clusters Lorentz5Momentum p3 = (*cl)->colParticle()->momentum() + (*cit)->antiColParticle()->momentum(); Lorentz5Momentum p4 = (*cit)->colParticle()->momentum() + (*cl)->antiColParticle()->momentum(); Energy oldMass = abs( p1.m() ) + abs( p2.m() ); Energy newMass = abs( p3.m() ) + abs( p4.m() ); if ( newMass < oldMass && newMass < minMass ) { minMass = newMass; candidate = cit; } } return candidate; } // forms two baryonic clusters from three clusters void ColourReconnector::_makeBaryonicClusters( ClusterPtr &c1, ClusterPtr &c2, ClusterPtr &c3, ClusterPtr &newcluster1, ClusterPtr &newcluster2) const{ //make sure they all have 2 components assert(c1->numComponents()==2); assert(c2->numComponents()==2); assert(c3->numComponents()==2); //abandon children c1->colParticle()->abandonChild(c1); c1->antiColParticle()->abandonChild(c1); c2->colParticle()->abandonChild(c2); c2->antiColParticle()->abandonChild(c2); c3->colParticle()->abandonChild(c3); c3->antiColParticle()->abandonChild(c3); newcluster1 = new_ptr(Cluster(c1->colParticle(),c2->colParticle(), c3->colParticle())); c1->colParticle()->addChild(newcluster1); c2->colParticle()->addChild(newcluster1); c3->colParticle()->addChild(newcluster1); newcluster1->setVertex(LorentzPoint()); newcluster2 = new_ptr(Cluster(c1->antiColParticle(), c2->antiColParticle(), c3->antiColParticle())); c1->antiColParticle()->addChild(newcluster2); c2->antiColParticle()->addChild(newcluster2); c3->antiColParticle()->addChild(newcluster2); newcluster2->setVertex(LorentzPoint()); } pair ColourReconnector::_reconnect(ClusterPtr &c1, ClusterPtr &c2) const { // choose the other possibility to form two clusters from the given // constituents assert(c1->numComponents()==2); assert(c2->numComponents()==2); int c1_col(-1),c1_anti(-1),c2_col(-1),c2_anti(-1); for(unsigned int ix=0;ix<2;++ix) { if (c1->particle(ix)->hasColour(false)) c1_col = ix; else if(c1->particle(ix)->hasColour(true )) c1_anti = ix; if (c2->particle(ix)->hasColour(false)) c2_col = ix; else if(c2->particle(ix)->hasColour(true )) c2_anti = ix; } assert(c1_col>=0&&c2_col>=0&&c1_anti>=0&&c2_anti>=0); ClusterPtr newCluster1 = new_ptr( Cluster( c1->colParticle(), c2->antiColParticle() ) ); newCluster1->setVertex(0.5*( c1->colParticle()->vertex() + c2->antiColParticle()->vertex() )); if(c1->isBeamRemnant(c1_col )) newCluster1->setBeamRemnant(0,true); if(c2->isBeamRemnant(c2_anti)) newCluster1->setBeamRemnant(1,true); ClusterPtr newCluster2 = new_ptr( Cluster( c2->colParticle(), c1->antiColParticle() ) ); newCluster2->setVertex(0.5*( c2->colParticle()->vertex() + c1->antiColParticle()->vertex() )); if(c2->isBeamRemnant(c2_col )) newCluster2->setBeamRemnant(0,true); if(c1->isBeamRemnant(c1_anti)) newCluster2->setBeamRemnant(1,true); return pair (newCluster1, newCluster2); } pair ColourReconnector::_reconnectBaryonic(ClusterPtr &c1, ClusterPtr &c2) const { // choose the other possibility to form two clusters from the given // constituents assert(c1->numComponents()==2); assert(c2->numComponents()==2); int c1_col(-1),c1_anti(-1),c2_col(-1),c2_anti(-1); for(unsigned int ix=0;ix<2;++ix) { if (c1->particle(ix)->hasColour(false)) c1_col = ix; else if(c1->particle(ix)->hasColour(true )) c1_anti = ix; if (c2->particle(ix)->hasColour(false)) c2_col = ix; else if(c2->particle(ix)->hasColour(true )) c2_anti = ix; } assert(c1_col>=0&&c2_col>=0&&c1_anti>=0&&c2_anti>=0); c1->colParticle()->abandonChild(c1); c2->antiColParticle()->abandonChild(c2); ClusterPtr newCluster1 = new_ptr( Cluster( c1->colParticle(), c2->antiColParticle() ) ); c1->colParticle()->addChild(newCluster1); c2->antiColParticle()->addChild(newCluster1); newCluster1->setVertex(0.5*( c1->colParticle()->vertex() + c2->antiColParticle()->vertex() )); if(c1->isBeamRemnant(c1_col )) newCluster1->setBeamRemnant(0,true); if(c2->isBeamRemnant(c2_anti)) newCluster1->setBeamRemnant(1,true); c1->antiColParticle()->abandonChild(c1); c2->colParticle()->abandonChild(c2); ClusterPtr newCluster2 = new_ptr( Cluster( c2->colParticle(), c1->antiColParticle() ) ); c1->antiColParticle()->addChild(newCluster2); c2->colParticle()->addChild(newCluster2); newCluster2->setVertex(0.5*( c2->colParticle()->vertex() + c1->antiColParticle()->vertex() )); if(c2->isBeamRemnant(c2_col )) newCluster2->setBeamRemnant(0,true); if(c1->isBeamRemnant(c1_anti)) newCluster2->setBeamRemnant(1,true); return pair (newCluster1, newCluster2); } pair ColourReconnector::_shuffle (const PVector & q, const PVector & aq, unsigned maxtries) const { const size_t nclusters = q.size(); assert (nclusters > 1); assert (aq.size() == nclusters); int i, j; unsigned tries = 0; bool octet; do { // find two different random integers in the range [0, nclusters) i = UseRandom::irnd( nclusters ); do { j = UseRandom::irnd( nclusters ); } while (i == j); // check if one of the two potential clusters would be a colour octet state octet = _isColour8( q[i], aq[j] ) || _isColour8( q[j], aq[i] ) ; tries++; } while (octet && tries < maxtries); if (octet) i = j = -1; return make_pair(i,j); } bool ColourReconnector::_isColour8(tcPPtr p, tcPPtr q) const { bool octet = false; // make sure we have a triplet and an anti-triplet if ( ( p->hasColour() && q->hasAntiColour() ) || ( p->hasAntiColour() && q->hasColour() ) ) { // true if p and q are originated from a colour octet if ( !p->parents().empty() && !q->parents().empty() ) { octet = ( p->parents()[0] == q->parents()[0] ) && ( p->parents()[0]->data().iColour() == PDT::Colour8 ); } // (Final) option: check if same colour8 parent // or already found an octet. if(_octetOption==0||octet) return octet; // (All) option handling more octets // by browsing particle history/colour lines. tColinePtr cline,aline; // Get colourlines form final states. if(p->hasColour() && q->hasAntiColour()) { cline = p-> colourLine(); aline = q->antiColourLine(); } else { cline = q-> colourLine(); aline = p->antiColourLine(); } // Follow the colourline of p. if ( !p->parents().empty() ) { tPPtr parent = p->parents()[0]; while (parent) { if(parent->data().iColour() == PDT::Colour8) { // Coulour8 particles should have a colour // and an anticolour line. Currently the // remnant has none of those. Since the children // of the remnant are not allowed to emit currently, // the colour octet remnant is handled by the return // statement above. The assert also catches other // colour octets without clines. If the children of // a remnant should be allowed to emit, the remnant // should get appropriate colour lines and // colour states. // See Ticket: #407 // assert(parent->colourLine()&&parent->antiColourLine()); octet = (parent-> colourLine()==cline && parent->antiColourLine()==aline); } if(octet||parent->parents().empty()) break; parent = parent->parents()[0]; } } } return octet; } void ColourReconnector::persistentOutput(PersistentOStream & os) const { os << _clreco << _preco << _precoBaryonic << _algorithm << _initTemp << _annealingFactor << _annealingSteps << _triesPerStepFactor << ounit(_maxDistance,femtometer) << _octetOption; } void ColourReconnector::persistentInput(PersistentIStream & is, int) { is >> _clreco >> _preco >> _precoBaryonic >> _algorithm >> _initTemp >> _annealingFactor >> _annealingSteps >> _triesPerStepFactor >> iunit(_maxDistance,femtometer) >> _octetOption; } void ColourReconnector::Init() { static ClassDocumentation documentation ("This class is responsible of the colour reconnection."); static Switch interfaceColourReconnection ("ColourReconnection", "Colour reconnections", &ColourReconnector::_clreco, 0, true, false); static SwitchOption interfaceColourReconnectionNo (interfaceColourReconnection, "No", "Colour reconnections off", 0); static SwitchOption interfaceColourReconnectionYes (interfaceColourReconnection, "Yes", "Colour reconnections on", 1); static Parameter interfaceMtrpAnnealingFactor ("AnnealingFactor", "The annealing factor is the ratio of the temperatures in two successive " "temperature steps.", &ColourReconnector::_annealingFactor, 0.9, 0.0, 1.0, false, false, Interface::limited); static Parameter interfaceMtrpAnnealingSteps ("AnnealingSteps", "Number of temperature steps in the statistical annealing algorithm", &ColourReconnector::_annealingSteps, 50, 1, 10000, false, false, Interface::limited); static Parameter interfaceMtrpTriesPerStepFactor ("TriesPerStepFactor", "The number of reconnection tries per temperature steps is the number of " "clusters times this factor.", &ColourReconnector::_triesPerStepFactor, 5.0, 0.0, 100.0, false, false, Interface::limited); static Parameter interfaceMtrpInitialTemp ("InitialTemperature", "Factor used to determine the initial temperature from the median of the " "energy change in a few random rearrangements.", &ColourReconnector::_initTemp, 0.1, 0.00001, 100.0, false, false, Interface::limited); static Parameter interfaceRecoProb ("ReconnectionProbability", "Probability that a found reconnection possibility is actually accepted", &ColourReconnector::_preco, 0.5, 0.0, 1.0, false, false, Interface::limited); static Parameter interfaceRecoProbBaryonic ("ReconnectionProbabilityBaryonic", "Probability that a found reconnection possibility is actually accepted", &ColourReconnector::_precoBaryonic, 0.5, 0.0, 1.0, false, false, Interface::limited); static Switch interfaceAlgorithm ("Algorithm", "Specifies the colour reconnection algorithm", &ColourReconnector::_algorithm, 0, true, false); static SwitchOption interfaceAlgorithmPlain (interfaceAlgorithm, "Plain", "Plain colour reconnection as in Herwig 2.5.0", 0); static SwitchOption interfaceAlgorithmStatistical (interfaceAlgorithm, "Statistical", "Statistical colour reconnection using simulated annealing", 1); static SwitchOption interfaceAlgorithmBaryonic (interfaceAlgorithm, "Baryonic", "Baryonic cluster reconnection", 2); static Parameter interfaceMaxDistance ("MaxDistance", "Maximum distance between the clusters at which to consider rearrangement" " to avoid colour reconneections of displaced vertices", &ColourReconnector::_maxDistance, femtometer, 1000.*femtometer, 0.0*femtometer, 1e100*femtometer, false, false, Interface::limited); static Switch interfaceOctetTreatment ("OctetTreatment", "Which octets are not allowed to be reconnected", &ColourReconnector::_octetOption, 0, false, false); static SwitchOption interfaceOctetTreatmentFinal (interfaceOctetTreatment, "Final", "Only prevent for the final (usuaslly non-perturbative) g -> q qbar splitting", 0); static SwitchOption interfaceOctetTreatmentAll (interfaceOctetTreatment, "All", "Prevent for all octets", 1); } diff --git a/Hadronization/LightClusterDecayer.cc b/Hadronization/LightClusterDecayer.cc --- a/Hadronization/LightClusterDecayer.cc +++ b/Hadronization/LightClusterDecayer.cc @@ -1,386 +1,389 @@ // -*- C++ -*- // // LightClusterDecayer.cc is a part of Herwig - A multi-purpose Monte Carlo event generator // Copyright (C) 2002-2017 The Herwig Collaboration // // Herwig is licenced under version 3 of the GPL, see COPYING for details. // Please respect the MCnet academic guidelines, see GUIDELINES for details. // // // This is the implementation of the non-inlined, non-templated member // functions of the LightClusterDecayer class. // #include "LightClusterDecayer.h" #include #include #include #include #include #include #include #include #include "Cluster.h" #include "CheckId.h" #include "Herwig/Utilities/Kinematics.h" #include using namespace Herwig; DescribeClass describeLightClusterDecayer("Herwig::LightClusterDecayer",""); IBPtr LightClusterDecayer::clone() const { return new_ptr(*this); } IBPtr LightClusterDecayer::fullclone() const { return new_ptr(*this); } void LightClusterDecayer::persistentOutput(PersistentOStream & os) const { os << _hadronSelector; } void LightClusterDecayer::persistentInput(PersistentIStream & is, int) { is >> _hadronSelector; } void LightClusterDecayer::Init() { static ClassDocumentation documentation ("There is the class responsible for the one-hadron decay of light clusters"); static Reference interfaceHadronSelector("HadronSelector", "A reference to the HadronSelector object", &Herwig::LightClusterDecayer::_hadronSelector, false, false, true, false); } bool LightClusterDecayer::decay(ClusterVector & clusters, tPVector & finalhadrons) { // Loop over all clusters, and for those that were not heavy enough // to undergo to fission, check if they are below the threshold // for normal two-hadron decays. If this is the case, then the cluster // should be decayed into a single hadron: this can happen only if // it is possible to reshuffle momenta between the cluster and // another one; in the rare occasions in which such exchange of momenta // is not possible (because all of the clusters are too light) then // the event is skipped. // Notice that, differently from what happens in Fortran Herwig, // light (that is below the threshold for the production of the lightest // pair of hadrons with the proper flavours) fission products, produced // by the fission of heavy clusters in class ClusterFissioner // have been already "decayed" into single hadron (the lightest one // with proper flavour) by the same latter class, without requiring // any reshuffling. Therefore the light clusters that are treated in // this LightClusterDecayer class are produced directly // (originally) by the class ClusterFinder. // To preserve all of the information, the cluster partner with which // the light cluster (that decays into a single hadron) exchanges // momentum in the reshuffling procedure is redefined and inserted // in the vector vecNewRedefinedCluPtr. Only at the end, when all // light clusters have been examined, the elements this vector will be // copied in collecCluPtr (the reason is that it is not allowed to // modify a STL container while iterating over it. At the same time, // this ensures that a cluster can be redefined only once, which seems // sensible although not strictly necessary). // Notice that the cluster reshuffling partner is normally redefined // and inserted in the vector vecNewRedefinedCluPtr, but not always: // in the case it is also light, then it is also decayed immediately // into a single hadron, without redefining it (the reason being that, // otherwise, the would-be redefined cluster could have undefined // components). vector redefinedClusters; + for (ClusterVector::const_iterator it = clusters.begin(); it != clusters.end(); ++it) { - + // Skip the clusters that are not available or that are // heavy, intermediate, clusters that have undergone to fission, - if ( ! (*it)->isAvailable() || ! (*it)->isReadyToDecay() ) continue; - + if ( ! (*it)->isAvailable() || ! (*it)->isReadyToDecay() ){ + continue; + } // We need to require (at least at the moment, maybe in the future we // could change it) that the cluster has exactly two components, // because otherwise we don't know how to deal with the kinematics. // If this is not the case, then send a warning because it is not suppose // to happen, and then do nothing with (ignore) such cluster. if ( (*it)->numComponents() != 2 ) { generator()->logWarning( Exception("LightClusterDecayer::decay " "***Still cluster with not exactly" " 2 components*** ", Exception::warning) ); continue; } // select the hadron for single hadron decay tcPDPtr hadron = _hadronSelector->chooseSingleHadron((*it)->particle(0)->dataPtr(), (*it)->particle(1)->dataPtr(), (**it).mass()); // if not single decay continue - if(!hadron) continue; - + if(!hadron){ +continue; + } // We assume that the candidate reshuffling cluster partner, // with whom the light cluster can exchange momenta, // is chosen as the closest in space-time between the available // clusters. Notice that an alternative, sensible approach // could be to consider instead the "closeness" in the colour // structure... // Notice that nor a light cluster (which decays into a single hadron) // neither its cluster reshuffling partner (which either has a // redefined cluster or also decays into a single hadron) can be // a reshuffling partner of another light cluster. // This because we are requiring that the considered candidate cluster // reshuffling partner has the status "isAvailable && isReadyToDecay" true; // furthermore, the new redefined clusters are not added to the collection // of cluster before the end of the entire reshuffling procedure, avoiding // in this way that the redefined cluster of a cluster reshuffling partner // is used again later. Needless to say, this is just an assumption, // although reasonable, but nothing more than that! // Build a multimap of available reshuffling cluster partners, // with key given by the module of the invariant space-time distance // w.r.t. the light cluster, so that this new collection is automatically // ordered in increasing distance values. // We use a multimap, rather than a map, just for precaution against not properly // defined cluster positions which could produce all identical (null) distances. multimap candidates; for ( ClusterVector::iterator jt = clusters.begin(); jt != clusters.end(); ++jt ) { if ((*jt)->isAvailable() && (*jt)->isReadyToDecay() && jt != it) { Length distance = abs (((*it)->vertex() - (*jt)->vertex()).m()); candidates.insert(pair(distance,*jt)); } } // Loop sequentially the multimap. multimap::const_iterator mmapIt = candidates.begin(); bool found = false; while (!found && mmapIt != candidates.end()) { found = reshuffling(hadron, *it, (*mmapIt).second, redefinedClusters, finalhadrons); if (!found) ++mmapIt; } if (!found) return partonicReshuffle(hadron,*it,finalhadrons); } // end loop over collecCluPtr // Add to collecCluPtr all of the redefined new clusters (indeed the // pointers to them are added) contained in vecNewRedefinedCluPtr. for (tClusterVector::const_iterator it = redefinedClusters.begin(); it != redefinedClusters.end(); ++it) { clusters.push_back(*it); } return true; } bool LightClusterDecayer::reshuffling(const tcPDPtr pdata1, tClusterPtr cluPtr1, tClusterPtr cluPtr2, tClusterVector & redefinedClusters, tPVector & finalhadrons) { // don't reshuffle with beam clusters if(cluPtr2->isBeamCluster()) return false; // This method does the reshuffling of momenta between the cluster "1", // that must decay into a single hadron (with id equal to idhad1), and // the candidate cluster "2". It returns true if the reshuffling succeed, // false otherwise. PPtr ptrhad1 = pdata1->produceParticle(); if ( ! ptrhad1 ) { generator()->logWarning( Exception("LightClusterDecayer::reshuffling" "***Cannot create a particle with specified id***", Exception::warning) ); return false; } Energy mhad1 = ptrhad1->mass(); // Let's call "3" and "4" the two constituents of the second cluster tPPtr part3 = cluPtr2->particle(0); tPPtr part4 = cluPtr2->particle(1); // Check if the system of the two clusters can kinematically be replaced by // an hadron of mass mhad1 (which is the lightest single hadron with the // same flavour numbers as the first cluster) and the second cluster. // If not, then try to replace the second cluster with the lightest hadron // with the same flavour numbers; if it still fails, then give up! Lorentz5Momentum pSystem = cluPtr1->momentum() + cluPtr2->momentum(); pSystem.rescaleMass(); // set the mass as the invariant of the quadri-vector Energy mSystem = pSystem.mass(); Energy mclu2 = cluPtr2->mass(); bool singleHadron = false; Energy mLHP2 = _hadronSelector->massLightestHadronPair(part3->dataPtr(),part4->dataPtr()); Energy mLH2 = _hadronSelector->massLightestHadron(part3->dataPtr(),part4->dataPtr()); if(mSystem > mhad1 + mclu2 && mclu2 > mLHP2) { singleHadron = false; } else if(mSystem > mhad1 + mLH2) { singleHadron = true; mclu2 = mLH2; } else return false; // Let's call from now on "Sys" the system of the two clusters, and // had1 (of mass mhad1) the lightest hadron in which the first // cluster decays, and clu2 (of mass mclu2) either the second // cluster or the lightest hadron in which it decays (depending // which one is kinematically allowed, see above). // The idea behind the reshuffling is to replace the system of the // two clusters by the system of the hadron had1 and (cluster or hadron) clu2, // but leaving the overall system unchanged. Furthermore, the motion // of had1 and clu2 in the Sys frame is assumed to be parallel to, respectively, // those of the original cluster1 and cluster2 in the same Sys frame. // Calculate the unit three-vector, in the frame "Sys" along which the // two initial clusters move. Lorentz5Momentum u( cluPtr1->momentum() ); u.boost( - pSystem.boostVector() ); // boost from LAB to Sys // Calculate the momenta of had1 and clu2 in the Sys frame first, // and then boost back in the LAB frame. Lorentz5Momentum phad1, pclu2; if (pSystem.m() < mhad1 + mclu2 ) { throw Exception() << "Impossible Kinematics in LightClusterDecayer::reshuffling()" << Exception::eventerror; } Kinematics::twoBodyDecay(pSystem, mhad1, mclu2, u.vect().unit(), phad1, pclu2); ptrhad1->set5Momentum( phad1 ); // set momentum of first hadron. ptrhad1->setVertex(cluPtr1->vertex()); // set hadron vertex position to the // parent cluster position. cluPtr1->addChild(ptrhad1); finalhadrons.push_back(ptrhad1); cluPtr1->flagAsReshuffled(); cluPtr2->flagAsReshuffled(); if(singleHadron) { // In the case that also the cluster reshuffling partner is light // it is decayed into a single hadron, *without* creating the // redefined cluster (this choice is justified in order to avoid // clusters that could have undefined components). PPtr ptrhad2 = _hadronSelector->lightestHadron(part3->dataPtr(),part4->dataPtr()) ->produceParticle(); ptrhad2->set5Momentum( pclu2 ); ptrhad2->setVertex( cluPtr2->vertex() ); // set hadron vertex position to the // parent cluster position. cluPtr2->addChild(ptrhad2); finalhadrons.push_back(ptrhad2); } else { // Create the new cluster which is the redefinitions of the cluster // partner (cluster "2") used in the reshuffling procedure of the // light cluster (cluster "1"). // The rationale of this is to preserve completely all of the information. ClusterPtr cluPtr2new = ClusterPtr(); if(part3 && part4) cluPtr2new = new_ptr(Cluster(part3,part4)); cluPtr2new->set5Momentum( pclu2 ); cluPtr2new->setVertex( cluPtr2->vertex() ); cluPtr2->addChild( cluPtr2new ); redefinedClusters.push_back( cluPtr2new ); // Set consistently the momenta of the two components of the second cluster // after the reshuffling. To do that we first calculate the momenta of the // constituents in the initial cluster rest frame; then we boost them back // in the lab but using this time the new cluster rest frame. Finally we store // these information in the new cluster. Notice that we do *not* set // consistently also the momenta of the (eventual) particles pointed by the // two components: that's because we do not need to do so, being the momentum // an explicit private member of the class Component (which is set equal // to the momentum of the eventual particle pointed only in the constructor, // but then later should not necessary be the same), and furthermore it allows // us not to loose any information, in the sense that we can always, later on, // to find the original momenta of the two components before the reshuffling. Lorentz5Momentum p3 = part3->momentum(); //p3new->momentum(); p3.boost( - (cluPtr2->momentum()).boostVector() ); // from LAB to clu2 (old) frame p3.boost( pclu2.boostVector() ); // from clu2 (new) to LAB frame Lorentz5Momentum p4 = part4->momentum(); //p4new->momentum(); p4.boost( - (cluPtr2->momentum()).boostVector() ); // from LAB to clu2 (old) frame p4.boost( pclu2.boostVector() ); // from clu2 (new) to LAB frame cluPtr2new->particle(0)->set5Momentum(p3); cluPtr2new->particle(1)->set5Momentum(p4); } // end of if (singleHadron) return true; } bool LightClusterDecayer::partonicReshuffle(const tcPDPtr had, const PPtr cluster, tPVector & finalhadrons) { tPPtr meson(cluster); if(!meson->parents().empty()) meson=meson->parents()[0]; if(!meson->parents().empty()) meson=meson->parents()[0]; // check b/c hadron decay int ptype(abs(meson->id())%10000); bool heavy = (ptype/1000 == 5 || ptype/1000 ==4 ); heavy |= (ptype/100 == 5 || ptype/100 ==4 ); heavy |= (ptype/10 == 5 || ptype/10 ==4 ); if(!heavy) return false; // find the leptons tPVector leptons; for(unsigned int ix=0;ixchildren().size();++ix) { if(!(meson->children()[ix]->dataPtr()->coloured())) { leptons.push_back(meson->children()[ix]); } } if(leptons.size()==1) { tPPtr w=leptons[0]; leptons.pop_back(); for(unsigned int ix=0;ixchildren().size();++ix) { if(!w->children()[ix]->dataPtr()->coloured()) { leptons.push_back(w->children()[ix]); } } } if(leptons.size()!=2) return false; // get momentum of leptonic system and the its minimum possible mass Energy mmin(ZERO); Lorentz5Momentum pw; for(unsigned int ix=0;ixmomentum(); mmin+=leptons[ix]->mass(); } pw.rescaleMass(); // check we can do the reshuffling PPtr ptrhad = had->produceParticle(); // total momentum fo the system Lorentz5Momentum pSystem = pw + cluster->momentum(); pSystem.rescaleMass(); // normal case get additional energy by rescaling momentum in rest frame of // system if(pSystem.mass()>ptrhad->mass()+pw.mass()&&pw.mass()>mmin) { // Calculate the unit three-vector, in the frame "Sys" along which the // two initial clusters move. Lorentz5Momentum u(cluster->momentum()); u.boost( - pSystem.boostVector() ); // Calculate the momenta of had1 and clu2 in the Sys frame first, // and then boost back in the LAB frame. Lorentz5Momentum phad1, pclu2; Kinematics::twoBodyDecay(pSystem, ptrhad->mass(), pw.mass(), u.vect().unit(), phad1, pclu2); // set momentum of first hadron. ptrhad->set5Momentum( phad1 ); // set hadron vertex position to the parent cluster position. ptrhad->setLabVertex(cluster->vertex()); // add hadron cluster->addChild(ptrhad); finalhadrons.push_back(ptrhad); // reshuffle the leptons // boost the leptons to the rest frame of the system Boost boost1(-pw.boostVector()); Boost boost2( pclu2.boostVector()); for(unsigned int ix=0;ixdeepBoost(boost1); leptons[ix]->deepBoost(boost2); } return true; } else { return false; } }