C++ Blob::InParticle方法代码示例

double Mixing_Handler::DetermineMixingTime(Particle* decayer,
                                           bool checkforpartstatus) const
{
  double time = decayer->Time();

  Blob* motherblob = decayer->ProductionBlob();
  if(motherblob->Type()==btp::Hadron_Mixing) {
    decayer = motherblob->InParticle(0);
    motherblob = decayer->ProductionBlob();
  }
  Particle* sister = NULL;
  if(motherblob->Type()!=btp::Fragmentation ||
     (motherblob->NInP()==2 && motherblob->NOutP()==2 &&
      motherblob->InParticle(0)->Flav()==motherblob->OutParticle(0)->Flav() &&
      motherblob->InParticle(1)->Flav()==motherblob->OutParticle(1)->Flav())) {
    // check if particle was produced coherently
    Particle_Vector sisters = motherblob->GetOutParticles();
    for(Particle_Vector::const_iterator it=sisters.begin(); it!=sisters.end(); it++) {
      if((*it)!=decayer && decayer->Flav()==(*it)->Flav().Bar()) {
        sister = (*it);
        break;
      }
    }
  }
  // in coherent production, special rules apply:
  if(sister) {
    if((checkforpartstatus && sister->Status()==part_status::decayed) ||
       (!checkforpartstatus && sister->DecayBlob()))
    {
      time = time - sister->Time();
    }
    else time = 0.0;
  }
  return time;
}

void Hard_Decay_Handler::CreateDecayBlob(ATOOLS::Particle* inpart)
{
  DEBUG_FUNC(inpart->Flav());
  if(inpart->DecayBlob()) THROW(fatal_error,"Decay blob already exists.");
  if(!Decays(inpart->Flav())) return;
  Blob* blob = p_bloblist->AddBlob(btp::Hard_Decay);
  blob->SetStatus(blob_status::needs_showers);
  blob->AddStatus(blob_status::needs_extraQED);
  blob->AddToInParticles(inpart);
  blob->SetTypeSpec("Sherpa");
  Decay_Table* table=p_decaymap->FindDecay(blob->InParticle(0)->Flav());
  if (table==NULL) {
    msg_Error()<<METHOD<<" decay table not found, retrying event."<<endl
               <<*blob<<endl;
    throw Return_Value::Retry_Event;
  }
  blob->AddData("dc",new Blob_Data<Decay_Channel*>(table->Select()));

  DEBUG_INFO("p_onshell="<<inpart->Momentum());
  blob->AddData("p_onshell",new Blob_Data<Vec4D>(inpart->Momentum()));
  DEBUG_INFO("succeeded.");
}

void AnalyseEvent(Blob_List* blobs)
{
#ifdef USING__ROOT
//   int outgoing = 1;
//   int incoming = -1;
//   Particle_List outparts = blobs->ExtractParticles(part_status::active, outgoing);

  ///////////////////////////////////////////////////////////////////////////////////
  // analyse primary decay blob, ignore subsequent decays                          //
  ///////////////////////////////////////////////////////////////////////////////////
  Blob * primarydecayblob = blobs->FindFirst(btp::Hadron_Decay);
//   msg_Out()<<"primary decay blob:"<<endl<<*primarydecayblob<<endl;
  // photon multiplicity and decay frame radiated energy (total)
  unsigned int photmult = 0;
  double       photener = 0.;
  for (int i=0; i<primarydecayblob->NOutP(); i++) {
    if ((primarydecayblob->OutParticle(i)->Flav().IsPhoton() == true) && 
        (primarydecayblob->OutParticle(i)->Info() == 'S')) {
      photmult++;
      photener = photener + primarydecayblob->OutParticle(i)->Momentum()[0];
    }
  }
  photonmultiplicity->Fill(photmult);
  if (photener != 0.)   decayframeenergy->Fill(photener);
  // multipole rest frame angles
  Vec4D multipolesum  = Vec4D(0.,0.,0.,0.);
  Vec4D axis          = Vec4D(0.,0.,0.,1.);
  std::list<Vec4D> multipole;
  std::list<Vec4D> newphot;
  for (int i=0; i<primarydecayblob->NOutP(); i++) {
    if (primarydecayblob->OutParticle(i)->Flav().Charge() != 0.) {
      multipolesum = multipolesum + primarydecayblob->OutParticle(i)->Momentum();
      multipole.push_back(primarydecayblob->OutParticle(i)->Momentum());
    }
  }
  if (primarydecayblob->InParticle(0)->Flav().Charge() != 0) {
    multipolesum = multipolesum + primarydecayblob->InParticle(0)->Momentum();
    multipole.push_front(primarydecayblob->InParticle(0)->Momentum());
  }
  Poincare boost(multipolesum);
  Poincare rotate;
  // charged initial state: rotate such that initial state at theta = 0
  if (mother_flav.Charge() != 0.) {
    Vec4D inmom = *multipole.begin();
    boost.Boost(inmom);
    rotate = Poincare(inmom,axis);
  }
  // neutral initial state: rotate such that heaviest charged final state at theta = 0
  else {
    std::list<Vec4D>::iterator heaviest = multipole.begin();
    for (std::list<Vec4D>::iterator iter=multipole.begin(); iter!=multipole.end(); iter++) {
      if (abs((iter->Abs2() - heaviest->Abs2())/(iter->Abs2() + heaviest->Abs2())) > 1E-6) {
        heaviest = iter;
      }
    }
    boost.Boost(*heaviest);
    rotate = Poincare(*heaviest,axis);
  }
  for (int i=0; i<primarydecayblob->NOutP(); i++) {
    if (primarydecayblob->OutParticle(i)->Flav().IsPhoton() == true) {
      Vec4D mom = primarydecayblob->OutParticle(i)->Momentum();
      boost.Boost(mom);
      rotate.Rotate(mom);
      double theta = acos((Vec3D(mom)*Vec3D(axis))/(Vec3D(mom).Abs()*Vec3D(axis).Abs()));
      multipoleframeangles->Fill(theta);
    }
  }  

  ///////////////////////////////////////////////////////////////////////////////////
  // inclusive analysis of whole decay chain                                       //
  ///////////////////////////////////////////////////////////////////////////////////

  // to be done ..
#endif
}