C++ Particle_Vector::push_back方法代码示例

void Lund_Interface::FillOutgoingParticlesInBlob(Blob *blob)
{
  Flavour    flav;
  Vec4D      momentum;
  Particle * particle;

  int n_q(0), n_g(0); Particle_Vector partons;
  for (int j=hepevt.jdahep[0][0]-1;j<hepevt.jdahep[0][1];j++) {
    int idhep = hepevt.idhep[j];
    flav=IdhepToSherpa(idhep);
    momentum=Vec4D(hepevt.phep[j][3],hepevt.phep[j][0],
		   hepevt.phep[j][1],hepevt.phep[j][2]);
    // don't fill blob position vector here, because decay is in CMS until boosting back
    particle = new Particle(-1,flav,momentum);
    particle->SetNumber(0);
    particle->SetStatus(part_status::active);
    particle->SetInfo('D');
    particle->SetFinalMass(hepevt.phep[j][4]);
    blob->AddToOutParticles(particle);
    
    if(abs(idhep)>0 && abs(idhep)<7) {
      n_q++;
      partons.push_back(particle);
    }
    else if(abs(idhep)==21) {
      n_g++;
      partons.push_back(particle);
    }
  }
  
  size_t n=partons.size();
  if(n>0) blob->SetStatus(blob_status::needs_showers);
  if(n_q==2 && n_g==0 && n==2) {
    if(partons[0]->Flav().IsAnti()) {
      partons[0]->SetFlow(2,-1);
      partons[1]->SetFlow(1,partons[0]->GetFlow(2));
    }
    else {
      partons[0]->SetFlow(1,-1);
      partons[1]->SetFlow(2,partons[0]->GetFlow(1));
    }
  }
  else if(n_q==0 && n_g==2 && n==2) {
    partons[0]->SetFlow(2,-1);
    partons[0]->SetFlow(1,-1);
    partons[1]->SetFlow(2,partons[0]->GetFlow(1));
    partons[1]->SetFlow(1,partons[0]->GetFlow(2));
  }
  else if(n_q==0 && n_g==3 && n==3) {
    partons[0]->SetFlow(2,-1);
    partons[0]->SetFlow(1,-1);
    partons[1]->SetFlow(2,partons[0]->GetFlow(1));
    partons[1]->SetFlow(1,-1);
    partons[2]->SetFlow(2,partons[1]->GetFlow(1));
    partons[2]->SetFlow(1,partons[0]->GetFlow(2));
  }
  else if(n>0) {
    msg_Error()<<METHOD<<" wasn't able to set the color flow for"<<endl<<*blob<<endl;
  }
}

void Signal_Process_FS_QED_Correction::FillBlob
(Blob * blob, const Subprocess_Info& spi, Particle_Vector& pv)
{
  // find the leptons owned by this subprocess
  Particle_Vector localpv;
  bool onlyleptons(true);
  for (size_t i=0;i<spi.m_ps.size();++i) {
    if (spi.m_ps[i].m_fl.Strong()) onlyleptons=false;
    for (Particle_Vector::iterator it=pv.begin();it!=pv.end();) {
      if ((*it)->Flav()==spi.m_ps[i].m_fl) {
        localpv.push_back(*it);
        pv.erase(it);
      }
      else ++it;
    }
  }
  if (onlyleptons) FillBlob(blob,spi.m_fl,localpv);
  else FillBlob(blob,DetermineGenericResonance(localpv),localpv);
}

void Hard_Decay_Handler::AssignPhotons(const Particle_Vector& daughters,
                                       ParticlePair_Vector& photons)
{
  // for every photon, find charged particle that's closest
  // ignore radiation off charged resonance for now
  if (photons.size()) {
    Particle_Vector cdaughters;
    for (size_t i(0);i<daughters.size();++i)
      if (daughters[i]->Flav().Charge()) cdaughters.push_back(daughters[i]);
    if (cdaughters.size()==1) {
      for (size_t i(0);i<photons.size();++i)
        photons[i].second=cdaughters[0];
    }
    else {
      Vec4D cmom(0.,0.,0.,0.);
      Vec4D_Vector cmoms;
      for (size_t i(0);i<cdaughters.size();++i) {
        cmoms.push_back(cdaughters[i]->Momentum());
        cmom+=cmoms[i];
      }
      Poincare ccms(cmom);
      for (size_t i(0);i<cdaughters.size();++i) ccms.Boost(cmoms[i]);
      for (size_t i(0);i<photons.size();++i){
        Vec4D pmom(photons[i].first->Momentum());
        ccms.Boost(pmom);
        size_t id(0);
        double dR(pmom.DR(cmoms[0]));
        for (size_t j(1);j<cmoms.size();++j) {
          double dRj(pmom.DR(cmoms[j]));
          if (dRj<dR) { id=j; dR=dRj; }
        }
        photons[i].second=cdaughters[id];
      }
    }
    for (size_t i(0);i<photons.size();++i) {
      if (photons[i].first==photons[i].second)
        THROW(fatal_error,"Photon has not been assigned.");
      msg_Debugging()<<photons[i].first->Flav()<<" "
                     <<photons[i].first->Momentum()
                     <<" assigned to "<<photons[i].second->Flav()<<std::endl;
    }
  }
}

bool Signal_Process_FS_QED_Correction::FindResonances
(Particle_Vector& pv,std::vector<Particle_Vector>& rpvs,Flavour_Vector& rfl,
 const Vertex_List& vlist)
{
  if (vlist.empty()) return false;
  DEBUG_FUNC("find resonances in "<<pv.size()<<" particles");
  // find a combination in pv for which a vertex exists such that the
  // IS flavour is on-shell within m_resdist times its width
  // book-keep first to later disentangle competing resonances
  // need to book-keep i,j,k,abs(mij-mk)/wk
  std::map<double,std::vector<size_t> > restab;
  for (size_t i(0);i<pv.size();++i) {
    for (size_t j(i+1);j<pv.size();++j) {
      for (size_t k(0);k<vlist.size();++k) {
        double mdist(abs((pv[i]->Momentum()+pv[j]->Momentum()).Mass()
                         -vlist[k]->in[0].Mass())/vlist[k]->in[0].Width());
        if (vlist[k]->nleg==3 &&
            ((pv[i]->Flav()==vlist[k]->in[1] &&
              pv[j]->Flav()==vlist[k]->in[2]) ||
             (pv[i]->Flav()==vlist[k]->in[2] &&
              pv[j]->Flav()==vlist[k]->in[1])) &&
            mdist<m_resdist) {
          size_t ida[3]={i,j,k};
          restab[mdist]=std::vector<size_t>(ida,ida+3);
        }
      }
    }
  }
  if (restab.empty()) {
    msg_Debugging()<<"no resonances found"<<std::endl;
    return false;
  }
  if (msg_LevelIsDebugging()) {
    msg_Debugging()<<"resonances found:\n";
    for (std::map<double,std::vector<size_t> >::const_iterator
         it=restab.begin();it!=restab.end();++it)
      msg_Debugging()<<it->second[0]<<it->second[1]<<it->second[2]<<": "
                     <<vlist[it->second[2]]->in[0]<<" -> "
                     <<vlist[it->second[2]]->in[1]<<" "
                     <<vlist[it->second[2]]->in[2]
                     <<", |m-M|/W="<<it->first<<std::endl;
  }
  Particle_Vector usedparts;
  for (std::map<double,std::vector<size_t> >::const_iterator it=restab.begin();
       it!=restab.end();++it) {
    bool valid(true);
    for (size_t i(0);i<usedparts.size();++i)
      if (pv[it->second[0]]==usedparts[i] ||
          pv[it->second[1]]==usedparts[i]) { valid=false; break; }
    if (!valid) continue;
    usedparts.push_back(pv[it->second[0]]);
    usedparts.push_back(pv[it->second[1]]);
    msg_Debugging()<<"constructing decay: "<<vlist[it->second[2]]->in[0]<<" -> "
                                           <<vlist[it->second[2]]->in[1]<<" "
                                           <<vlist[it->second[2]]->in[2]<<"\n";
    rfl.push_back(vlist[it->second[2]]->in[0]);
    Particle_Vector parts;
    parts.push_back(pv[it->second[0]]);
    parts.push_back(pv[it->second[1]]);
    rpvs.push_back(parts);
  }
  for (Particle_Vector::iterator it=usedparts.begin();it!=usedparts.end();++it)
    for (Particle_Vector::iterator pit=pv.begin();pit!=pv.end();++pit)
      if (*it==*pit) { pv.erase(pit); break; }
  return true;
}

Return_Value::code Signal_Process_FS_QED_Correction::Treat
(Blob_List * bloblist, double & weight)
{
  if (!m_on) return Return_Value::Nothing;
  if (bloblist->empty()) {
    msg_Error()<<"Signal_Process_FS_QED_Correction::Treat"
	       <<"("<<bloblist<<","<<weight<<"): "<<endl
               <<"   Blob list contains "<<bloblist->size()<<" entries."<<endl
               <<"   Continue and hope for the best."<<endl;
    return Return_Value::Error;
  }
  // look for QCD corrected hard process in need for QED
  Blob * sigblob(bloblist->FindLast(btp::Shower));
  if (!sigblob) return Return_Value::Nothing;
  // if already treated -> nothing to do
  if (sigblob->TypeSpec()=="YFS-type_QED_Corrections_to_ME")
    return Return_Value::Nothing;
  if (sigblob->TypeSpec()=="setting_leptons_on-shell")
    return Return_Value::Nothing;
  // extract FS leptons
  // two vectors -> the ones from the blob and the ones to be massive
  DEBUG_FUNC(m_qed);
  Particle_Vector fslep(sigblob->GetOutParticles());
  Particle_Vector mfslep;
  for (Particle_Vector::iterator it=fslep.begin();it!=fslep.end();) {
    if ((*it)->Flav().Strong() || (*it)->Flav().IsDiQuark() || 
	(*it)->DecayBlob()!=NULL) {
      fslep.erase(it);
    }
    else {
      mfslep.push_back(new Particle(-1,(*it)->Flav(),(*it)->Momentum(),'F'));
      (*mfslep.rbegin())->SetNumber(0);
      (*mfslep.rbegin())->SetOriginalPart(*it);
      (*mfslep.rbegin())->SetFinalMass((*it)->FinalMass());
      ++it;
    }
  }
  // if no leptons, nothing to do
  // if only one lepton, cannot do anything
  if (fslep.size()<2) {
    sigblob->UnsetStatus(blob_status::needs_extraQED);
    for (Particle_Vector::iterator it=mfslep.begin();it!=mfslep.end();++it)
      delete *it;
    return Return_Value::Nothing;
  }
  // if switched off or no need for QED stop here and build a blob
  if (!m_qed || !sigblob->Has(blob_status::needs_extraQED)) {
    Blob * onshellblob = bloblist->AddBlob(btp::QED_Radiation);
    onshellblob->SetTypeSpec("setting_leptons_on-shell");
    if (sigblob->Has(blob_status::needs_extraQED))
      sigblob->UnsetStatus(blob_status::needs_extraQED);
    for (Particle_Vector::iterator it=fslep.begin();it!=fslep.end();++it) {
      (*it)->SetInfo('H');
      (*it)->SetStatus(part_status::decayed);
      onshellblob->AddToInParticles(*it);
    }
    for (Particle_Vector::iterator it=mfslep.begin();it!=mfslep.end();++it) {
      onshellblob->AddToOutParticles(*it);
    }
    onshellblob->SetStatus(blob_status::needs_hadronization);
    return Return_Value::Success;
  }
  // put them on-shell (spoils consistency of pertubative calculation,
  // but necessary for YFS)
  if (!PutOnMassShell(mfslep)) {
    msg_Error()<<"Signal_Process_FS_QED_Correction::Treat("
	       <<bloblist<<","<<weight<<"): \n"
               <<"  Leptons could not be put on their mass shell.\n"
               <<"  Trying new event.\n"
               <<"  The event contained a ";
    for (Particle_Vector::iterator it=mfslep.begin();it!=mfslep.end();++it)
       msg_Error()<<(*it)->Flav().ShellName()<<"-";
    if (mfslep.size()==2) msg_Error()<<"pair";
    else                  msg_Error()<<"set";
    msg_Error()<<" of too little invariant mass to be put\n"
	       <<"  on their mass shell. If you are sensitive to this specific"
	       <<" signature consider\n  to set the respective particles"
	       <<" massive in the perturbative calculation using\n"
	       <<"  'MASSIVE[<id>]=1' to avoid this problem.\n";
    for (Particle_Vector::iterator it=mfslep.begin();it!=mfslep.end();++it)
      delete *it;
    return Return_Value::New_Event;
  }
  // build effective verteces for resonant production
  // use subprocess infos if possible
  Blob_Vector blobs = BuildResonantBlobs(mfslep);
  // add radiation
  for (Blob_Vector::iterator it=blobs.begin();it!=blobs.end();++it) {
    // do nothing if no resonance determined
    if ((*it)->InParticle(0)->Flav().Kfcode()!=kf_none) {
      (*it)->SetStatus(blob_status::needs_extraQED);
      if (!p_sphotons->AddRadiation(*it)) {
        msg_Error()<<"Signal_Process_FS_QED_Correction::Treat("<<bloblist
                   <<","<<weight<<"): "<<endl
                   <<"  Higher order QED corrections failed."<<endl
                   <<"  Retrying event."<<endl;
        for (Particle_Vector::iterator it=mfslep.begin();it!=mfslep.end();++it)
          delete *it;
        for (Blob_Vector::iterator it=blobs.begin();it!=blobs.end();++it)
          delete *it;
//.........这里部分代码省略.........

void Hard_Decay_Handler::AddDecayClustering(ATOOLS::Cluster_Amplitude*& ampl,
                                            Blob* blob,
                                            size_t& imax,
                                            size_t idmother)
{
  DEBUG_FUNC("blob->Id()="<<blob->Id()<<" idmother="<<ID(idmother));
  DEBUG_VAR(*blob);
  Particle_Vector daughters;
  ParticlePair_Vector photons;
  for (size_t i(0);i<blob->GetOutParticles().size();++i) {
    Particle * p(blob->OutParticle(i));
    if (p->Info()=='S') photons.push_back(make_pair(p,p));
    else                daughters.push_back(p);
  }
  std::sort(photons.begin(),photons.end(),ParticlePairFirstEnergySort());
  msg_Debugging()<<"daughters: ";
  for (size_t i(0);i<daughters.size();++i)
    msg_Debugging()<<daughters[i]->Flav().IDName()<<" ";
  msg_Debugging()<<" +  "<<photons.size()<<" soft photon(s)"<<std::endl;
  AssignPhotons(daughters,photons);
  if (daughters.size()==2) {
    msg_Debugging()<<"1 to 2 case"<<std::endl;
    Cluster_Amplitude* copy=ampl->InitPrev();
    copy->CopyFrom(ampl);
    copy->SetFlag(1);
    copy->SetMS(ampl->MS());
    Cluster_Leg *lij(ampl->IdLeg(idmother));
    copy->SetKT2(lij->Mom().Abs2());
    for (size_t i=0; i<ampl->Legs().size(); ++i)
      ampl->Leg(i)->SetStat(ampl->Leg(i)->Stat()|1);
    lij->SetStat(1|2|4);
    size_t idk(0);
    for (size_t i=0; i<copy->Legs().size(); ++i) {
      copy->Leg(i)->SetK(0);
      if (copy->Leg(i)->Id()!=idmother &&
          (copy->Leg(i)->Col().m_i==lij->Col().m_j ||
           copy->Leg(i)->Col().m_j==lij->Col().m_i))
        idk=copy->Leg(i)->Id();
    }
    if (lij->Col().m_i==0 && lij->Col().m_j==0) {
      // Ad hoc EW partner
      size_t ampl_nout=ampl->Legs().size()-ampl->NIn();
      if (ampl_nout==1) idk=ampl->Leg(0)->Id();
      else {
        size_t select=ampl->Legs().size();
        do {
          select=ampl->NIn()+floor(ran->Get()*ampl_nout);
        } while (ampl->Leg(select)->Id()&idmother ||
                 select>ampl->Legs().size()-1);
        idk=ampl->Leg(select)->Id();
      }
    }
    if (idk==0) THROW(fatal_error,"Colour partner not found");
    lij->SetK(idk);
    Cluster_Leg *d1(copy->IdLeg(idmother));
    size_t stat1(0), stat2(0);
    d1->SetMom(RecombinedMomentum(daughters[0],photons,stat1));
    d1->SetStat(stat1);
    d1->SetFlav(daughters[0]->Flav());
    copy->CreateLeg(RecombinedMomentum(daughters[1],photons,stat2),
                    daughters[1]->RefFlav());
    size_t idnew=1<<(++imax);
    copy->Legs().back()->SetId(idnew);
    copy->Legs().back()->SetStat(stat2);
    Cluster_Amplitude::SetColours(ampl->IdLeg(idmother),
                                  copy->IdLeg(idmother),
                                  copy->Legs().back());
    DEBUG_VAR(*copy);
    Cluster_Amplitude* tmp=copy;
    while (tmp->Next()) {
      tmp=tmp->Next();
      for (size_t i=0; i<tmp->Legs().size(); ++i) {
        if (tmp->Leg(i)->Id()&idmother) {
          tmp->Leg(i)->SetId(tmp->Leg(i)->Id()|idnew);
	}
        if (tmp->Leg(i)->K()&idmother) {
          tmp->Leg(i)->SetK(tmp->Leg(i)->K()|idnew);
        }
      }
      DEBUG_VAR(*tmp);
    }
    std::vector<size_t> ids;
    ids.push_back(idmother);
    ids.push_back(idnew);
    while (photons.size())
      AddPhotonsClustering(copy, daughters, photons, imax, ids);
    if (daughters[0]->DecayBlob())
      AddDecayClustering(copy, daughters[0]->DecayBlob(), imax, idmother);
    if (daughters[1]->DecayBlob())
      AddDecayClustering(copy, daughters[1]->DecayBlob(), imax, idnew);
    ampl=copy;
  }
  else if (daughters.size()==3) {
    msg_Debugging()<<"1 to 3 case"<<std::endl;
    // structure m -> 0 P[->1 2]
    // propagator always combines daughters 1+2
    Cluster_Amplitude* step1=ampl->InitPrev();
    step1->CopyFrom(ampl);
    step1->SetFlag(1);
    step1->SetMS(ampl->MS());
//.........这里部分代码省略.........