C++ Flavour::IDName方法代码示例

Hadron_Decay_Channel* Mixing_Handler::Select(Particle* decayer,
                                             const Hadron_Decay_Table& ot) const
{
  Flavour flav = decayer->Flav();
  string tag = flav.IsAnti() ? flav.Bar().IDName() : flav.IDName();
  if(m_model("Interference_"+tag,0.0)!=0.0) {
    double lifetime = DetermineMixingTime(decayer,false);
    bool anti_at_t0 = decayer->Flav().IsAnti();
    if(decayer->ProductionBlob()->Type()==btp::Hadron_Mixing) anti_at_t0 = !anti_at_t0;
    if(lifetime!=0.0) {
      Hadron_Decay_Table table(ot);
      double cos_term = cos(flav.DeltaM()/rpa->hBar()*lifetime);
      double sin_term = sin(flav.DeltaM()/rpa->hBar()*lifetime);
      double GX, GR, asymmetry, a;
      for(size_t i=0; i<table.size(); i++) {
        Hadron_Decay_Channel* hdc = table.at(i);
        if(hdc->CPAsymmetryS()==0.0 && hdc->CPAsymmetryC()==0.0) continue;
        if(flav.DeltaGamma()==0.0) {
          asymmetry = hdc->CPAsymmetryS()*sin_term - hdc->CPAsymmetryC()*cos_term;
        }
        else {
          Complex lambda = hdc->CPAsymmetryLambda();
          double l2 = sqr(abs(lambda));
          double dG = flav.DeltaGamma()/rpa->hBar();
          asymmetry = (2.0*lambda.imag()/(1.0+l2)*sin_term - (1.0-l2)/(1.0+l2)*cos_term)/
              (cosh(dG*lifetime/2.0) - 2.0*lambda.real()/(1.0+l2)*sinh(dG*lifetime/2.0));
        }
        GX = hdc->Width(); // partial width of this DC
        GR = table.TotalWidth()-GX; // partial width of other DCs
        if(asymmetry>0.0)
          a = -1.0*GR/2.0/GX/asymmetry+sqrt(sqr(GR)/4.0/sqr(GX)/sqr(asymmetry)+(GR+GX)/GX);
        else if(asymmetry<0.0)
          a = -1.0*GR/2.0/GX/asymmetry-sqrt(sqr(GR)/4.0/sqr(GX)/sqr(asymmetry)+(GR+GX)/GX);
        else
          a = 0.0;
        if(anti_at_t0) table.UpdateWidth(hdc, (1.0+a)*GX);
        else           table.UpdateWidth(hdc, (1.0-a)*GX);
      }
      return table.Select();
    }
  }
  return ot.Select();
}

ATOOLS::Blob* Mixing_Handler::PerformMixing(Particle* decayer) const
{
  // explicit mixing in event record
  Flavour flav = decayer->Flav();
  string tag = flav.IsAnti() ? flav.Bar().IDName() : flav.IDName();
  if(m_model("Mixing_"+tag,0.0)!=0.0 && decayer->Info()!=char('M')) {
    double t = DetermineMixingTime(decayer,true)/rpa->hBar();
    if(t==0.0) return NULL;
    double factor = decayer->Flav().QOverP2();
    if(decayer->Flav().IsAnti()) factor = 1.0/factor;
    double dG = decayer->Flav().DeltaGamma()*t/4.0;
    double dm = decayer->Flav().DeltaM()*t/2.0;
    Complex i(0.0,1.0);
    double prob_not_mix = sqr(abs(exp(i*dm)*exp(dG)+exp(-i*dm)*exp(-dG)));
    double prob_mix = factor*sqr(abs(exp(i*dm)*exp(dG)-exp(-i*dm)*exp(-dG)));
    if(prob_mix > ran->Get()*(prob_mix+prob_not_mix)) {
      if(decayer->DecayBlob()) {
        decayer->DecayBlob()->RemoveOwnedParticles();
        delete decayer->DecayBlob();
      }
      decayer->SetStatus(part_status::decayed);
      decayer->SetInfo('m');
      Particle* mixed_part = new Particle(0, decayer->Flav().Bar(),
                                          decayer->Momentum(), 'M');
      mixed_part->SetFinalMass(decayer->FinalMass());
      mixed_part->SetStatus(part_status::active);
      mixed_part->SetTime(decayer->Time());
      Blob* mixingblob = new Blob();
      mixingblob->SetType(btp::Hadron_Mixing);
      mixingblob->SetId();
      mixingblob->SetStatus(blob_status::inactive);
      mixingblob->SetTypeSpec("HADRONS");
      mixingblob->AddToInParticles(decayer);
      mixingblob->AddToOutParticles(mixed_part);
      mixingblob->SetPosition(decayer->ProductionBlob()->Position());
      mixingblob->SetStatus(blob_status::needs_hadrondecays);
      return mixingblob;
    }
  }
  return NULL;
}

void InitialiseAnalysis()
{
#ifdef USING__ROOT
  map<string, string> tags=Read_Write_Base::GlobalTags();
  string filepiece;
  map<string, string>::const_iterator it=tags.find("TAG_FILE_PIECE");
  if(it!=tags.end())
    filepiece="."+it->second;
  Hadron_Decay_Table* table=hadrons->DecayMap()->FindDecay(mother_flav);
  for(int i(0);i<table->size();++i) {
    Hadron_Decay_Channel* hdc=table->at(i);
    string fname = hdc->FileName();
    rootfiles[hdc]=
      new TFile(("Analysis/"+fname.erase(fname.length()-4)+filepiece+".root").c_str(),"RECREATE");

    TH1D* hist;
    int currenthist=0;
    int noutp(hdc->NOut());
    for(int i=0; i<noutp; i++) {
      Flavour flav = hdc->GetDecayProduct(i);
      string name = "costheta_"+flav.IDName()+"_"+ToString(i);
      string xtitle = "cos(#Theta) of "+flav.RootName();
      string ytitle = "#frac{1}{#Gamma} #frac{d#Gamma}{dcos(#Theta)}";
      hist = makeTH1D(name, "", 50, -1.0, 1.0, xtitle, ytitle);
      ThetaHists[hdc].push_back(hist);
      currenthist++;
    }
  
    currenthist=0;
    for(int i=0; i<noutp; i++) {
      Flavour flav = hdc->GetDecayProduct(i);
      string name = "E_"+flav.IDName()+"_"+ToString(i);
      string xtitle = "E of "+flav.RootName()+" [GeV]";
      string ytitle = "#frac{1}{#Gamma} #frac{d#Gamma}{dE} [GeV^{-1}]";
      double M = mother_flav.HadMass();
      double othermass = 0.0;
      for(int k=0; k<noutp; k++) {
	if(k!=i) othermass+=hdc->GetDecayProduct(k).HadMass();
      }
      double Emin = 0.9*flav.HadMass();
      double Emax = 1.1/2.0/M*(sqr(M)+sqr(flav.HadMass())-othermass);
      hist = makeTH1D(name, "", 50, Emin, Emax, xtitle, ytitle);
      EnergyHists[hdc].push_back(hist);
      currenthist++;
    }

    currenthist=0;
    for(int i=0; i<noutp; i++) {
      if(noutp<3) break;
      for(int j=i+1; j<noutp; j++) {
	Flavour flav1 = hdc->GetDecayProduct(i);
	Flavour flav2 = hdc->GetDecayProduct(j);
	string name = "q2_"+flav1.IDName()+"_"+flav2.IDName()+"_"+ToString(i)+ToString(j);
	string xtitle = "Inv. Mass q^{2}=(p_{"+flav1.RootName()+"}+p_{"
	  +flav2.RootName()+"})^{2} [GeV^{2}]";
	string ytitle = "#frac{1}{#Gamma} #frac{d#Gamma}{dq^{2}} [GeV^{-2}]";
	double othermass = 0.0;
	for(int k=0; k<noutp; k++) {
	  if(k==i || k==j) continue;
	  othermass+=hdc->GetDecayProduct(k).HadMass();
	}
	double q2min = 0.0;
	double q2max = 1.1*sqr(hdc->GetDecaying().HadMass()-othermass);
	hist = makeTH1D(name.c_str(), "", 50, q2min, q2max, xtitle, ytitle);
	TwoInvMassHists[hdc].push_back(hist);
	currenthist++;
      }
    }

    currenthist=0;
    for(int i=0; i<noutp; i++) {
      if(noutp<4) break;
      for(int j=i+1; j<noutp; j++) {
	for(int k=j+1; k<noutp; k++) {
	  Flavour flav1 = hdc->GetDecayProduct(i);
	  Flavour flav2 = hdc->GetDecayProduct(j);
	  Flavour flav3 = hdc->GetDecayProduct(k);
	  string name = "q2_"+flav1.IDName()+"_"+flav2.IDName()+"_"+flav3.IDName()
	    +"_"+ToString(i)+ToString(j)+ToString(k);
	  string xtitle = "Inv. Mass q^{2}=(p_{"+flav1.RootName()+
	    "}+p_{"+flav2.RootName()+"}+p_{"+flav3.RootName()+"})^{2} [GeV^{2}]";
	  string ytitle = "#frac{1}{#Gamma} #frac{d#Gamma}{dq^{2}} [GeV^{-2}]";
	  double othermass = 0.0;
	  for(int l=0; l<noutp; l++) {
	    if(l==i || l==j || l==k) continue;
	    othermass+=hdc->GetDecayProduct(k).HadMass();
	  }
	  double q2min = 0.0;
	  double q2max = 1.1*sqr(hdc->GetDecaying().HadMass());
	  hist = makeTH1D(name.c_str(), "", 50, q2min, q2max, xtitle, ytitle);
	  ThreeInvMassHists[hdc].push_back(hist);
	  currenthist++;
	}
      }
    }
  }
#endif
}