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

void AnalyseEvent(Blob_List* blobs)
{
  Blob* blob = blobs->FindFirst(btp::Hadron_Decay);
  if(blob==NULL) return;
#ifdef USING__ROOT
  Hadron_Decay_Channel* hdc=(*blob)["hdc"]->Get<Hadron_Decay_Channel*>();

  int currenthist=0;
  for(int i=0; i<blob->NOutP(); i++) {
    double costheta = blob->OutParticle(i)->Momentum().CosTheta();
    ThetaHists[hdc][currenthist]->Fill(costheta);
    currenthist++;
  }
  
  currenthist=0;
  for(int i=0; i<blob->NOutP(); i++) {
    double energy = blob->OutParticle(i)->Momentum()[0];
    EnergyHists[hdc][currenthist]->Fill(energy);
    currenthist++;
  }

  currenthist=0;
  for(int i=0; i<blob->NOutP(); i++) {
    if(blob->NOutP()<3) break;
    for(int j=i+1; j<blob->NOutP(); j++) {
      Vec4D mom1 = blob->OutParticle(i)->Momentum();
      Vec4D mom2 = blob->OutParticle(j)->Momentum();
      double q2 = (mom1+mom2).Abs2();
      TwoInvMassHists[hdc][currenthist]->Fill(q2);
      currenthist++;
    }
  }

  currenthist=0;
  for(int i=0; i<blob->NOutP(); i++) {
    if(blob->NOutP()<4) break;
    for(int j=i+1; j<blob->NOutP(); j++) {
      for(int k=j+1; k<blob->NOutP(); k++) {
        Vec4D mom1 = blob->OutParticle(i)->Momentum();
        Vec4D mom2 = blob->OutParticle(j)->Momentum();
        Vec4D mom3 = blob->OutParticle(k)->Momentum();
        double q2 = (mom1+mom2+mom3).Abs2();
        ThreeInvMassHists[hdc][currenthist]->Fill(q2);
        currenthist++;
      }
    }
  }
#endif
}

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 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
}