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

void Soft_Cluster_Handler::FixHHDecay(Cluster * cluster,Blob * blob,
				      const Flavour had1,const Flavour had2,
				      const bool & constrained)
{
  double M       = cluster->Mass(), M2 = M*M;
  double m12     = sqr(had1.HadMass()), m22 = sqr(had2.HadMass());

  cluster->BoostInCMSAndRotateOnZ();

  double E1((M2+m12-m22)/(2.*M)), pl2(sqr(E1)-m12);
  bool isbeam(false);
  double stheta, pt2;
  double masscor((cluster->GetTrip()->m_flav.HadMass() *
		  cluster->GetAnti()->m_flav.HadMass())/m_pt02);
  do { 
    stheta = 1.-2.*ran->Get(); 
    pt2    = pl2*sqr(stheta);
  } while (pt2>m_pt2max*m_pt2maxfac || 
	   sqr((*p_as)(pt2,false)/p_as->MaxValue())<ran->Get());
  double pt     = sqrt(pt2);
  int sign      = cluster->GetTrip()->m_mom[3]<0?-1:1;
  double pl1    = sign*sqrt(sqr(E1)-sqr(pt)-m12);
  double cosphi = cos(2.*M_PI*ran->Get()), sinphi = sqrt(1.-cosphi*cosphi);
  Vec4D  p1     = Vec4D(E1,pt*cosphi,pt*sinphi,pl1);
  Vec4D  p2     = cluster->Momentum()-p1;

  if (p1[0]<0. || p2[0]<0.) throw Return_Value::Retry_Event;

  cluster->RotateAndBoostBack(p1);
  cluster->RotateAndBoostBack(p2);
  cluster->RotateAndBoostBack();

  Particle * left(new Particle(-1,had1,p1));
  left->SetNumber();
  left->SetInfo('P');
  left->SetFinalMass(had1.HadMass());
  Particle * right(new Particle(-1,had2,p2));
  right->SetNumber();
  right->SetInfo('P');
  right->SetFinalMass(had2.HadMass());
  control::s_AHAparticles+=2;

  if (blob!=NULL) {
    blob->AddToOutParticles(left);
    blob->AddToOutParticles(right);
  }
  // if (cluster->GetTrip()->m_info=='B' || cluster->GetAnti()->m_info=='B') {
  //   msg_Out()<<"==========================================================\n"
  // 	     <<"Cluster decay (pt = "<<pt<<") for cluster \n"<<(*cluster)<<"\n"
  // 	     <<"==> "<<left->Momentum()<<" + "<<right->Momentum()<<" for "
  // 	     <<left->Flav()<<" + "<<right->Flav()<<"\n"
  // 	     <<"==========================================================\n";
  // }
  if (m_ana) {
    Histogram* histo((m_histograms.find(std::string("PT_HH")))->second);
    histo->Insert(pt);
    Histogram* histo2((m_histograms.find(std::string("PT2_HH")))->second);
    histo2->Insert(pt*pt);
  }
}

void SM_U1_B::FillSpectrum(const PDF::ISR_Handler_Map& isr) {
  p_dataread->RereadInFile();
  p_constants->insert(make_pair(string("M_Z'"),    
				p_dataread->GetValue<double>("M_Z'",-1.)));
  p_constants->insert(make_pair(string("g'_1"),    
				p_dataread->GetValue<double>("g'_1",0.)));

  double MZprime(ScalarConstant(string("M_Z'"))),GZprime(0.),inc;
  double g1prime(sqr(ScalarConstant(std::string("g'_1"))));
  FixMix();
  msg_Out()<<"Calculate width for Z' with mass = "<<MZprime<<" GeV.\n";

  for (short int i=1;i<=6;i++) {
    Flavour quark = Flavour((kf_code)(i));
    double  massq = quark.HadMass(); 
    if (!quark.IsOn() || !quark.Strong()) continue;
    if (i%2) {
      if (2.*massq>MZprime) continue;
      GZprime += inc = 
	3.*MZprime*g1prime/(24.*M_PI)*pow(1.-sqr(2.*massq/MZprime),3./2.);
      msg_Out()<<"   Add Z' --> "<<quark<<"+"<<quark.Bar()<<" ("<<i<<"), "
	       <<" add "<<inc<<" GeV.\n";
    }
    else {
      double M2 = MZprime*MZprime, mq2 = massq*massq;
      for (short int j=2;j<=7;j+=2) {
	Flavour anti = Flavour((kf_code)(j)).Bar();
	if (!anti.IsOn() || !anti.Strong()) continue;
	double massa  = anti.HadMass(), ma2 = massa*massa; 
	if (massq+massa>MZprime) continue;
	Complex mix   = ComplexMatrixElement(std::string("UpMix"),i/2-1,j/2-1);
	double absmix = ATOOLS::sqr(mix.real())+ATOOLS::sqr(mix.imag());
	if (ATOOLS::IsZero(absmix)) continue;
	GZprime += inc = 
	  3.*g1prime*absmix/(24.*M_PI*M2)*
	  ((2.*M2-mq2-ma2)/2.-3.*massq*massa-ATOOLS::sqr(ma2-mq2)/(2.*M2)) *
	  sqrt(ATOOLS::sqr(M2-mq2-ma2)-4.*ma2*mq2)/(2.*MZprime);
	msg_Out()<<"   Add Z' --> "<<quark<<"+"<<anti<<" ("<<i<<"), "
		 <<" add "<<inc<<" GeV.\n";
      }
    }
  }

  Flavour flav;
  flav = Flavour(kf_Z0_2);
  flav.SetMass(ScalarConstant(string("M_Z'")));
  flav.SetHadMass(MZprime);
  flav.SetMassOn(true);
  flav.SetWidth(GZprime);

  msg_Out()<<METHOD<<" initializes Z' boson with \n"
	   <<"    mass = "<<MZprime<<" GeV and width = "<<GZprime<<" GeV\n"
	   <<"    for g'_1 = "<<ScalarConstant(std::string("g'_1"))<<".\n";
}

void InitialiseGenerator(int argc, char *argv[])
{
  if(argc<2) {
    cout<<"Usage: ./SingleDecay <PDG_CODE>"<<endl;
    THROW(normal_exit,"you didn't specify the decaying particle by PDG code.");
  }

  small_sherpa_init(argc, argv);

  hadrons = new SHERPA::Hadron_Decay_Handler(".", "Fragmentation.dat");

  mother_flav = Flavour( (kf_code) abs(ToType<int>(argv[1])) );
  mother_flav.SetStable(false);
  if(ToType<int>(argv[1])<0) mother_flav=mother_flav.Bar();
  if(hadrons->DecayMap()->FindDecay(mother_flav)==NULL)
    THROW(fatal_error, "Didn't find "+ToString<Flavour>(mother_flav)+
	  " in HadronDecays.dat.");
  
  // set all decay channel BR's equal, such that we generate the same amount of
  // each decay channel to be tested
  PHASIC::Decay_Table* table=hadrons->DecayMap()->FindDecay(mother_flav);
  for(size_t i(0);i<table->size();++i)
    table->UpdateWidth(table->at(i), 1.0);

  rpa->gen.SetEcms(mother_flav.HadMass());
  msg_Info()<<"Welcome. I am decaying a "<<mother_flav<<endl;
}

double Soft_Cluster_Handler::
TransformKin(const double MC,const Flavour & flav,const bool & enforce) {
  double mass2(sqr(flav.HadMass()));
  double width2(sqr(Max(flav.Width(),1.e-6)));
  if (!enforce && sqr(MC*MC-mass2)>10.*mass2*width2) return 0.;
  return
    pow(sqr(mass2)/(sqr(MC*MC-mass2) + mass2*width2),m_kappa) * 
    pow(mass2*width2/(sqr(MC*MC-mass2) + mass2*width2),m_lambda);
}

void InitialiseGenerator(int argc, char *argv[])
{
  if(argc<2) {
    cout<<"Usage: ./SingleDecay <PDG_CODE>"<<endl;
    THROW(normal_exit,"you didn't specify the decaying particle by PDG code.");
  }

  small_sherpa_init(argc, argv);

  hadrons = new SHERPA::Hadron_Decay_Handler(".", "Fragmentation.dat");

  Data_Reader * reader = new Data_Reader(" ",";","!","=");
  reader->AddWordSeparator("\t");
  reader->SetInputPath("./");
  reader->SetInputFile("YFS.dat");
  photons = new PHOTONS::Photons(reader,true);

  mother_flav = Flavour( (kf_code) abs(ToType<int>(argv[1])) );
  mother_flav.SetStable(false);
  if(ToType<int>(argv[1])<0) mother_flav=mother_flav.Bar();

  rpa->gen.SetEcms(mother_flav.HadMass());
  msg_Info()<<"Welcome. I am decaying a "<<mother_flav<<endl;

  Particle* mother_part = new Particle( 1,mother_flav,
                                        Vec4D(mother_flav.HadMass(),0.,0.,0.) );
  mother_part->SetTime();
  mother_part->SetFinalMass(mother_flav.HadMass());
  
  ref_blob = new Blob();
  ref_blob->SetType(btp::Hadron_Decay);
  ref_blob->SetStatus(blob_status::needs_hadrondecays);
  ref_blob->AddToInParticles(mother_part);

  try {
    hadrons->FillOnshellDecay(ref_blob, NULL);
  } catch (Return_Value::code ret) {
    msg_Error()<<METHOD<<" Something went wrong for blob: "<<ref_blob<<endl;
    return;
  }
}

void Lund_Interface::AdjustProperties(Flavour flav)
{
  int kc = pycomp(SherpaToIdhep(flav));
  if(kc>500) return;
  // adjust mass
  double pythiamass = pydat2.pmas[1-1][kc-1];
  double sherpamass = flav.HadMass();
  flav.SetMass(pythiamass);
  if( !(abs(sherpamass-pythiamass)/sherpamass < 1.e-2) ) {
    msg_Info()<<METHOD<<" Adjusted mass of "<<flav<<" ("<<flav.Kfcode()
        <<") from "<<sherpamass<<" to "<<pythiamass<<" to allow Pythia decays."<<endl;
  }
}

Blob_List* GenerateEvent()
{
  Blob_List* blobs = new Blob_List();
  Particle* mother_part = new Particle( 1,mother_flav,Vec4D(mother_flav.HadMass(),0.,0.,0.) );
  mother_part->SetTime();
  mother_part->SetFinalMass(mother_flav.HadMass());
  
  Blob* blob = blobs->AddBlob(btp::Hadron_Decay);
  blob->SetStatus(blob_status::needs_hadrondecays);
  blob->AddToInParticles(mother_part);

  try {
    hadrons->FillOnshellDecay(blob, NULL);
  } catch (Return_Value::code ret) {
    msg_Error()<<METHOD<<" Something went wrong for event: "<<*blobs
        <<endl;
    return blobs;
  }

  hadrons->CleanUp();

  msg_Events()<<*blobs<<std::endl;
  return blobs;
}

void InitialiseAnalysis()
{
#ifdef USING__ROOT
  msg_Out()<<"initialising ROOT analysis ..."<<std::endl;
  ATOOLS::MakeDir("PhotonAnalysisDirectory/SingleDecay/"+mother_flav.ShellName()+"_decays",true,493);
  rootfile = new TFile(string("PhotonAnalysisDirectory/SingleDecay/"+
                              mother_flav.ShellName()+"_decays/"+
                              mother_flav.ShellName()+"__"+
                              "DAUGHTERS"+".root").c_str(), "RECREATE");
  photonmultiplicity   = makeTH1D("photon_multiplicity","",
                                  10, -0.5, 10.5,
                                  Flavour(kf_photon).RootName()+" multiplicity","Events");
  decayframeenergy     = makeTH1D("decayframeenergy","",
                                  1000, 0., mother_flav.HadMass(),
                                  "total energy radiated in decay frame","Events");
  multipoleframeangles = makeTH1D("multipoleframeangles","",
                                  1000, 0., M_PI,
                                  "angular radiation pattern","Events");
#endif
}

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
}