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

double Continued_PDF::XPDF(const Flavour & flav,const bool & defmax) {
  if (flav.IsDiQuark()) {
    if (m_bunch==Flavour(kf_p_plus) && !flav.IsAnti()) {
      return XPDF(Flavour(kf_u));
    }
    if (m_bunch==Flavour(kf_p_plus).Bar() && flav.IsAnti()) {
      return XPDF(Flavour(kf_u).Bar());
    }
    return 0.;
  }
  if (m_x<m_xmin) return 0.;
  if (m_Q2>m_Q02) return p_pdf->GetXPDF(flav);
  double seapart(0.), valpart(0.);
  if (m_bunch==Flavour(kf_p_plus)) {
    if (flav==Flavour(kf_u) || flav==Flavour(kf_d)) {
      seapart = p_pdf->GetXPDF(flav.Bar())*(m_Q2/m_Q02); 
      valpart = p_pdf->GetXPDF(flav)-p_pdf->GetXPDF(flav.Bar());
    }
    else if (flav==Flavour(kf_gluon)) {
      seapart = p_pdf->GetXPDF(flav)*(m_Q2/m_Q02); 
      valpart = 1./m_gnorm * pow(1.-m_x,m_geta) * pow(m_x,m_glambda) *
	m_Snorm * (1.-m_Q2/m_Q02);
      /*
	((p_pdf->GetXPDF(Flavour(kf_u))-p_pdf->GetXPDF(Flavour(kf_u).Bar()) +
	p_pdf->GetXPDF(Flavour(kf_d))-p_pdf->GetXPDF(Flavour(kf_d).Bar()))/  
	m_Vnorm) *
      */ 
    }
    else
      seapart = p_pdf->GetXPDF(flav)*(m_Q2/m_Q02);
  }
  else if (m_bunch==Flavour(kf_p_plus).Bar()) {
    if (flav==Flavour(kf_u).Bar() || flav==Flavour(kf_d).Bar()) {
      seapart = p_pdf->GetXPDF(flav.Bar())*(m_Q2/m_Q02); 
      valpart = p_pdf->GetXPDF(flav)-p_pdf->GetXPDF(flav.Bar());
    }
    else if (flav==Flavour(kf_gluon)) {
      seapart = p_pdf->GetXPDF(flav)*(m_Q2/m_Q02); 
      valpart = 
	(p_pdf->GetXPDF(Flavour(kf_u).Bar())-p_pdf->GetXPDF(Flavour(kf_u)) +
	 p_pdf->GetXPDF(Flavour(kf_d).Bar())-p_pdf->GetXPDF(Flavour(kf_d))) *  
	m_Snorm/m_Vnorm * (1.-m_Q2/m_Q02);
    }
    else
      seapart = p_pdf->GetXPDF(flav)*(m_Q2/m_Q02);
  }
  double total = seapart+valpart;
  if (defmax && total>m_xpdfmax[flav]) {
    m_xmaxpdf[flav] = m_x;
    m_xpdfmax[flav] = total;
  }
  return total;
}

double Hadron_Remnant::MinimalEnergy(const ATOOLS::Flavour &flavour) 
{
  if (!m_initialized) {
    if (!flavour.Strong()) {
      return p_beam->Beam().HadMass();
    }
    bool found(false);
    kf_code di[3];
    if (flavour.IsQuark()) {
      short unsigned int j=0;
      for (Flavour_Vector::const_iterator flit(m_constit.begin());
	   flit!=m_constit.end();++flit) {
	if (found || flavour!=*flit) di[j++]=flit->Kfcode();
	else found=true;
      }
    }
    Flavour difl;
    if (!found || flavour.IsGluon()) {
      int single=-1;
      for (size_t i=0;i<m_constit.size();++i) {
	for (size_t j=i+1;j<m_constit.size();++j) {
	  if (m_constit[i]==m_constit[j]) { 
	    single=j; 
	    break; 
	  }
	}
	if (single>0) break;
      }
      Flavour fl(m_constit[single]);
      for (short unsigned int j=0, i=0;i<3;i++) 
	if (i!=single) di[j++]=m_constit[i].Kfcode();
      if (di[0]>di[1]) difl=Flavour((kf_code)(di[0]*1000+di[1]*100+1));
      else if (di[1]>di[0]) difl=Flavour((kf_code)(di[1]*1000+di[0]*100+1));
      else difl=Flavour((kf_code)(di[0]*1100+3));
      if (m_constit[single].IsAnti()) difl=difl.Bar();
      return difl.Mass()+fl.Mass()+flavour.Bar().Mass();
    }
    if (di[0]>di[1]) difl=Flavour((kf_code)(di[0]*1000+di[1]*100+1));
    else if (di[1]>di[0]) difl=Flavour((kf_code)(di[1]*1000+di[0]*100+1));
    else difl=Flavour((kf_code)(di[0]*1100+3));
    if (m_constit[0].IsAnti()) difl=difl.Bar();
    return difl.Mass();
  }
  else {
    if (flavour.IsQuark()) return flavour.Bar().Mass();
  }
  return 0.;
}

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

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

int Lund_Interface::PrepareFragmentationBlob(Blob * blob) 
{
  int nhep = 0;
  hepevt.idhep[nhep]=SherpaToIdhep(Flavour(kf_photon));
  for (short int j=1;j<4;++j) hepevt.phep[nhep][j-1]=blob->CMS()[j];
  hepevt.phep[nhep][3]=blob->CMS()[0];
  double pabs=(blob->CMS()).Abs2();
  if (pabs<0) hepevt.phep[nhep][4]=0.0;
  else hepevt.phep[nhep][4]=sqrt(pabs);
  for (short int j=0;j<4;++j) hepevt.vhep[nhep][j]=0.0;
  hepevt.isthep[nhep]=1;
  hepevt.jmohep[nhep][0]=0;
  hepevt.jmohep[nhep][1]=0;
  hepevt.jdahep[nhep][0]=0;
  hepevt.jdahep[nhep][1]=0;
  
  // gluon splittings
  for (int i(0);i<blob->NInP();++i) {
    Particle * part = blob->InParticle(i);
  if (part->GetFlow(1)!=0 && part->GetFlow(2)!=0) {
    Flavour            flav = Flavour(kf_d);
    if (ran->Get()<0.5) flav = Flavour(kf_u);
      Particle *help1(new Particle(-1,flav,0.5*part->Momentum()));
      Particle *help2(new Particle(-1,flav.Bar(),help1->Momentum()));
    help1->SetStatus(part_status::active);
    help2->SetStatus(part_status::active);
    AddPartonToString(help1,nhep);
    delete help1;
      unsigned int lastc(part->GetFlow(2));
      for (++i;i<blob->NInP();++i) {
      part = blob->InParticle(i);
      AddPartonToString(part,nhep);
	if (part->GetFlow(1)==lastc) {
	  lastc=0;
	  break;
	}
    }      
      if (lastc!=0)
	msg_Error()<<METHOD<<"(): Error. Open color string."<<std::endl;
    AddPartonToString(help2,nhep);
    delete help2;
      lastc=0;
  }
  else {
      for (;i<blob->NInP();i++) {
      part = blob->InParticle(i);
      AddPartonToString(part,nhep);
	if (part->GetFlow(1)==0) break;
      }  
    }  
  }
  return nhep;
}

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

int Vertex::SetVertex(Single_Vertex& orig, Single_Vertex& probe, int i0, int i1, int i2, int i3,int mode)
{
  probe = orig;

  if ((orig.dec>0 && orig.dec&4) &&
      !((i0==1 && i1==2 && i2==3) || 
	(i0==-2 && i1==3 && i2==-1) ||
	(i0==-3 && i1==-1 && i2==2))) return 0;
  if ((orig.dec>0 && orig.dec&2) && !(i0==1 && i1==2 && i2==3)) return 0;
  
  if (i0<0) probe.in[0] = orig.in[-i0-1].Bar();
       else probe.in[0] = orig.in[i0-1];
  if (i1<0) probe.in[1] = orig.in[-i1-1].Bar();
       else probe.in[1] = orig.in[i1-1];
  if (i2<0) probe.in[2] = orig.in[-i2-1].Bar();
       else probe.in[2] = orig.in[i2-1];
  if (orig.nleg==4) {
    if (i3<0) probe.in[3] = orig.in[-i3-1].Bar();
         else if (i3<99) probe.in[3] = orig.in[i3-1];
  }
  if (CheckExistence(probe)==0) return 0;
  if (probe.nleg==3) {
    if (FermionRule(probe)==0) return 0;}
  int hc = 0;

  int cnt = 0;
  for(int i=0;i<orig.nleg;i++) {
    if(orig.in[i]!=orig.in[i].Bar()) cnt++;
  }
  if(cnt>0){
    Flavour *flavlist= new Flavour[cnt];
    int *flaglist= new int[cnt];
    cnt = 0;
    for(int i=0;i<orig.nleg;i++){
      if(orig.in[i]!=orig.in[i].Bar()) {
	flavlist[cnt] = orig.in[i];
	if (i==0) flavlist[cnt] = flavlist[cnt].Bar();
	flaglist[cnt] = 0;
	cnt++;
      }
    }
    for (int i=0;i<cnt;i++) {
      for (int j=i+1;j<cnt;j++) {
	if (flavlist[i]==flavlist[j].Bar()) flaglist[i]=flaglist[j]=1;
      }
    }
    for (int i=0;i<cnt;i++) {
      if (flaglist[i]==0) hc = 1;
    }
    delete[] flavlist;
    delete[] flaglist;
  }
  
  int probehc = 0;
  if (hc) {
    // probe = h.c. ???
    for (short int i=0;i<orig.nleg;i++) { 
      Flavour flav = orig.in[i];
      if (flav!=flav.Bar()) {
	if (i==0) flav = flav.Bar();
	for (short int j=0;j<orig.nleg;j++) {
	  Flavour flav2 = probe.in[j];
	  if (j==0) flav2 = flav2.Bar();
	  if (flav2!=flav2.Bar()) {
	    if (flav==flav2.Bar()) {
	      probehc = 1;
	      break;
	    }
	  }
	}
	if (probehc) break;
      }
    }
    if (probehc) {
      int conjugate = 1;

      for (short int i=0;i<orig.nleg;i++) {
	//pseudoscalar
	if (orig.in[i]==Flavour(kf_A0)) conjugate *= -1;
      }
      
      if (orig.Lorentz.front()->Type()=="SSV" ||
	  orig.Lorentz.front()->Type()=="VVV") conjugate *= -1;
      
      if (conjugate==-1) {
	for (short int i=0;i<4;i++) probe.cpl[i] = -probe.cpl[i];
      }

      probe.Color.front().Conjugate();

       if (probe.Lorentz.front()->String()=="1") {
	//exchange left and right
	Kabbala help = probe.cpl[0];
	probe.cpl[0] = probe.cpl[1];
	probe.cpl[1] = help;
       }
    }
  } 
  if (orig.nleg==3 && (orig.Lorentz.front()->Type()=="FFV")) {
    //exchange left and right for 'barred' FFV vertices
//.........这里部分代码省略.........

Comix1to3::Comix1to3(const vector<Flavour>& flavs, const Flavour& prop,
                     size_t nonprop, size_t propi, size_t propj) :
  Spin_Amplitudes(flavs,Complex(0.0,0.0)), m_cur(4), m_anticur(4), m_nhel(4),
  m_prop(prop)
{
  DEBUG_FUNC(flavs<<" with prop "<<prop<<" in "<<propi<<","<<propj);
  assert(nonprop>0 && propi>0 && propj>0);
  if (flavs.size()!=4) THROW(fatal_error,"Internal error.");
  Vec4D k(1.0,0.0,1.0,0.0);

  for (size_t i(0);i<4;++i) {
    Current_Key ckey(i==0?flavs[i].Bar():flavs[i],MODEL::s_model,1);
    m_cur[i] = Current_Getter::GetObject("D"+ckey.Type(),ckey);
    if (m_cur[i]==NULL) THROW(fatal_error, "current not found");
    m_cur[i]->SetDirection(i==0?1:-1);
    m_cur[i]->SetId(std::vector<int>(1,i));
    m_cur[i]->InitPols(std::vector<int>(1,m_spins[i]));
    m_cur[i]->SetKey(i);
    m_cur[i]->SetGauge(k);
    m_nhel[i]=NHel(flavs[i]);
  }
  // s-channel for prop (i,j)
  Current_Key ckey(prop,MODEL::s_model,2);
  m_scur = Current_Getter::GetObject("D"+ckey.Type(),ckey);
  METOOLS::Int_Vector isfs(2), ids(2), pols(2);
  isfs[0]=flavs[propi].IsFermion();
  isfs[1]=flavs[propj].IsFermion();
  pols[0]=m_spins[ids[0]=propi];
  pols[1]=m_spins[ids[1]=propj];
  m_scur->SetId(ids);
  m_scur->SetFId(isfs);
  m_scur->FindPermutations();
  // final current (1,2,3)
  ckey=Current_Key(flavs[0],MODEL::s_model,1);
  m_fcur = Current_Getter::GetObject("D"+ckey.Type(),ckey);
  METOOLS::Int_Vector isfs2(3), ids2(3), pols2(3);
  isfs2[0]=flavs[1].IsFermion();
  isfs2[1]=flavs[2].IsFermion();
  isfs2[2]=flavs[3].IsFermion();
  pols2[0]=m_spins[ids2[0]=1];
  pols2[1]=m_spins[ids2[1]=2];
  pols2[2]=m_spins[ids2[2]=3];
  m_fcur->SetId(ids2);
  m_fcur->SetFId(isfs2);
  m_fcur->FindPermutations();
  // connect (2) & (3) into (2,3)
  m_v1=ConstructVertices(m_cur[propi], m_cur[propj], m_scur);
  DEBUG_VAR(m_v1.size());
  // connect (1) & (2,3) into (1,2,3)
  m_v2=ConstructVertices(m_cur[nonprop],m_scur,m_fcur);
  DEBUG_VAR(m_v2.size());
  m_scur->Print();
  m_fcur->Print();
  m_scur->InitPols(pols);
  m_fcur->InitPols(pols2);
  m_fcur->HM().resize(m_n);
  for (size_t i(0);i<m_n;++i) m_fcur->HM()[i]=i;

  for (size_t i(0);i<4;++i) {
    ckey=Current_Key(i==0?flavs[i]:flavs[i].Bar(),MODEL::s_model,1);
    m_anticur[i] = Current_Getter::GetObject("D"+ckey.Type(),ckey);
    if (m_anticur[i]==NULL) THROW(fatal_error, "current not found");
    m_anticur[i]->SetDirection(i==0?1:-1);
    m_anticur[i]->SetId(std::vector<int>(1,i));
    m_anticur[i]->InitPols(std::vector<int>(1,m_spins[i]));
    m_anticur[i]->SetKey(i);
    m_anticur[i]->SetGauge(k);
  }
  // s-channel for prop (2,3)
  ckey=Current_Key(prop.Bar(),MODEL::s_model,2);
  m_antiscur = Current_Getter::GetObject("D"+ckey.Type(),ckey);
  m_antiscur->SetId(ids);
  m_antiscur->SetFId(isfs);
  m_antiscur->FindPermutations();
  // final current (1,2,3)
  ckey=Current_Key(flavs[0].Bar(),MODEL::s_model,1);
  m_antifcur = Current_Getter::GetObject("D"+ckey.Type(),ckey);
  m_antifcur->SetId(ids2);
  m_antifcur->SetFId(isfs2);
  m_antifcur->FindPermutations();
  // connect (2) & (3) into (2,3)
  m_antiv1=ConstructVertices(m_anticur[propi], m_anticur[propj], m_antiscur);
  DEBUG_VAR(m_antiv1.size());
  // connect (1) & (2,3) into (1,2,3)
  m_antiv2=ConstructVertices(m_anticur[nonprop],m_antiscur,m_antifcur);
  DEBUG_VAR(m_antiv2.size());
  m_antiscur->Print();
  m_antifcur->Print();
  m_antiscur->InitPols(pols);
  m_antifcur->InitPols(pols2);
  m_antifcur->HM().resize(m_n);
  for (size_t i(0);i<m_n;++i) m_antifcur->HM()[i]=i;

  p_ci=new Color_Integrator();
  Idx_Vector cids(4,0);
  METOOLS::Int_Vector acts(4,0), types(4,0);
  for (size_t i(0);i<flavs.size();++i) {
    cids[i]=i;
    acts[i]=flavs[i].Strong();
    if (acts[i]) {
//.........这里部分代码省略.........