GitLab's annual major release is around the corner. Along with a lot of new and exciting features, there will be a few breaking changes. Learn more here.

CS_Shower.C 27.6 KB
Newer Older
1
#include "CSSHOWER++/Main/CS_Shower.H"
2

3
#include "CSSHOWER++/Showers/Splitting_Function_Base.H"
4
#include "PHASIC++/Selectors/Jet_Finder.H"
5
#include "PDF/Main/Jet_Criterion.H"
6 7
#include "ATOOLS/Phys/Cluster_Amplitude.H"
#include "ATOOLS/Org/Run_Parameter.H"
8
#include "ATOOLS/Org/MyStrStream.H"
9
#include "ATOOLS/Org/Exception.H"
10
#include "ATOOLS/Org/My_Limits.H"
11 12

using namespace CSSHOWER;
13
using namespace PHASIC;
14 15 16 17 18 19 20 21
using namespace ATOOLS;
using namespace std;

CS_Shower::CS_Shower(PDF::ISR_Handler *const _isr,MODEL::Model_Base *const model,
		     Data_Reader *const _dataread) : 
  Shower_Base("CSS"), p_isr(_isr), 
  p_shower(NULL), p_cluster(NULL), p_ampl(NULL)
{
22
  rpa->gen.AddCitation
23 24 25 26 27 28 29 30 31
    (1,"The Catani-Seymour subtraction based shower is published under \\cite{Schumann:2007mg}.");
  int maxem=_dataread->GetValue<int>("CSS_MAXEM",-1);
  if (maxem<0) m_maxem=std::numeric_limits<size_t>::max();
  else {
    m_maxem=maxem;
    msg_Info()<<METHOD<<"(): Set max emissions "<<m_maxem<<"\n";
  }
  m_kmode=_dataread->GetValue<int>("CSS_KMODE",1);
  if (m_kmode!=1) msg_Info()<<METHOD<<"(): Set kernel mode "<<m_kmode<<"\n";
32 33
  m_recocheck=_dataread->GetValue<int>("CSS_RECO_CHECK",0);
  if (m_recocheck!=0) msg_Info()<<METHOD<<"(): Set reco check mode "<<m_recocheck<<"\n";
34
  m_cmode=ToType<int>(rpa->gen.Variable("METS_CLUSTER_MODE"));
35 36
  int amode(_dataread->GetValue<int>("EXCLUSIVE_CLUSTER_MODE",0));
  if (amode!=0) msg_Info()<<METHOD<<"(): Set exclusive cluster mode "<<amode<<".\n";
37 38 39
  
  m_weightmode = int(_dataread->GetValue<int>("WEIGHT_MODE",1));
  
40
  int _qed=_dataread->GetValue<int>("CSS_EW_MODE",0);
41
  if (_qed==1) {
42
    s_kftable[kf_photon]->SetResummed();
43
  }
44 45 46
  p_shower = new Shower(_isr,_qed,_dataread);
  
  p_next = new All_Singlets();
47
  p_refs = new All_Singlets();
48 49 50

  p_cluster = new CS_Cluster_Definitions(p_shower,m_kmode);
  p_cluster->SetAMode(amode);
51 52 53 54 55 56 57
}

CS_Shower::~CS_Shower() 
{
  CleanUp();
  if (p_shower)      { delete p_shower; p_shower = NULL; }
  if (p_cluster)     { delete p_cluster; p_cluster = NULL; }
58
  if (p_ampl) p_ampl->Delete();
59
  delete p_next;
60
  delete p_refs;
61 62
}

63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93
void CS_Shower::RefCopy()
{
  p_refs->clear();
  std::map<Parton*,Parton*> pmap;
  std::map<Parton*,Parton*>::iterator pit;
  for (All_Singlets::const_iterator sit(m_allsinglets.begin());
       sit!=m_allsinglets.end();++sit) (*sit)->RefCopy(p_refs,pmap);
  for (size_t i(0);i<m_allsinglets.size();++i) {
    pit=pmap.find(m_allsinglets[i]->GetLeft());
    if (pit!=pmap.end()) (*p_refs)[i]->SetLeft(pit->second);
    pit=pmap.find(m_allsinglets[i]->GetRight());
    if (pit!=pmap.end()) (*p_refs)[i]->SetRight(pit->second);
    pit=pmap.find(m_allsinglets[i]->GetSplit());
    if (pit!=pmap.end()) (*p_refs)[i]->SetSplit(pit->second);
    pit=pmap.find(m_allsinglets[i]->GetSpec());
    if (pit!=pmap.end()) (*p_refs)[i]->SetSpec(pit->second);
    for (Singlet::const_iterator it(m_allsinglets[i]->begin());
	 it!=m_allsinglets[i]->end();++it) {
      Parton *c(pmap[*it]);
      pit=pmap.find((*it)->GetLeft());
      if (pit!=pmap.end()) c->SetLeft(pit->second);
      pit=pmap.find((*it)->GetRight());
      if (pit!=pmap.end()) c->SetRight(pit->second);
      pit=pmap.find((*it)->GetPrev());
      if (pit!=pmap.end()) c->SetPrev(pit->second);
      pit=pmap.find((*it)->GetNext());
      if (pit!=pmap.end()) c->SetNext(pit->second);
    }
  }
}

94 95 96
int CS_Shower::PerformShowers(const size_t &maxem,size_t &nem)
{
  if (!p_shower) return 1;
97
  m_weight=1.0;
98 99 100
  RefCopy();
  for (All_Singlets::const_iterator rit(p_refs->begin()),
	 sit(m_allsinglets.begin());sit!=m_allsinglets.end();++sit,++rit) {
101
    msg_Debugging()<<"before shower step\n";
102 103 104 105 106 107 108 109
    for (Singlet::const_iterator it((*sit)->begin());it!=(*sit)->end();++it)
      if ((*it)->GetPrev()) 
	if((*it)->GetPrev()->GetNext()==*it) {
	  if ((*it)->GetPrev()->TMin()==std::numeric_limits<double>::max())
	    (*it)->SetStart((*it)->GetPrev()->KtNext());
	  else (*it)->SetStart((*it)->GetPrev()->KtStart());
	}
    msg_Debugging()<<**sit;
110
    size_t pem(nem);
111
    if (!p_shower->EvolveShower(*sit,maxem,nem)) return 0;
112
    m_weight*=p_shower->Weight();
113 114 115 116
    if ((*sit)->GetLeft()) {
      p_shower->ReconstructDaughters(*sit,1);
      p_shower->ReconstructDaughters(*rit,1);
    }
117 118
    msg_Debugging()<<"after shower step with "<<nem-pem
		   <<" of "<<nem<<" emission(s)\n";
119
    msg_Debugging()<<**sit;
120 121 122 123 124 125 126
    msg_Debugging()<<"\n";
  }
  return 1;
}

int CS_Shower::PerformShowers() 
{
127
  return PerformShowers(m_maxem,m_nem);
128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160
}

int CS_Shower::PerformDecayShowers() {
  if (!p_shower) return 1;
  size_t nem(0);
  for (All_Singlets::const_iterator 
	 asit(m_allsinglets.begin());asit!=m_allsinglets.end();++asit) {
    if (!p_shower->EvolveShower(*asit,m_maxem,nem)) return 0;
  }
  return 1;
}

bool CS_Shower::ExtractPartons(Blob_List *const blist) {
  
  Blob * psblob(blist->FindLast(btp::Shower));
  if (psblob==NULL) THROW(fatal_error,"No Shower blob");
  psblob->SetTypeSpec("CSSHOWER++1.0");
  for (int i=0;i<psblob->NInP();++i) 
    psblob->InParticle(i)->SetStatus(part_status::decayed);
  for (int i=0;i<psblob->NOutP();++i) 
    psblob->OutParticle(i)->SetStatus(part_status::decayed);
  
  psblob->SetStatus(blob_status::needs_beams |
		    blob_status::needs_hadronization);
  
  for (All_Singlets::const_iterator 
	 sit(m_allsinglets.begin());sit!=m_allsinglets.end();++sit)
      (*sit)->ExtractPartons(psblob,p_ms);
  return true;
}

void CS_Shower::CleanUp()
{
161 162
  m_nem=0;
  m_kt2prev=std::numeric_limits<double>::max();
163 164 165 166 167
  for (All_Singlets::const_iterator 
	 sit(m_allsinglets.begin());sit!=m_allsinglets.end();++sit) {
    if (*sit) delete *sit;
  }
  m_allsinglets.clear();
168
  p_refs->clear();
169 170
}

171
PDF::Cluster_Definitions_Base * CS_Shower::GetClusterDefinitions() 
172 173 174 175
{
  return p_cluster;
}

176 177 178 179 180 181 182
void CS_Shower::GetKT2Min(Cluster_Amplitude *const ampl,const size_t &id,
			  KT2X_Map &kt2xmap,std::set<size_t> &aset)
{
  msg_Indent();
  for (size_t i(0);i<ampl->Legs().size();++i) {
    Cluster_Leg *cl(ampl->Leg(i));
    if ((cl->Id()&id)==0) continue;
183 184
    if (ampl->Prev()) GetKT2Min(ampl->Prev(),cl->Id(),kt2xmap,aset);
    if (cl->Stat()==3 || cl->Stat()==5) {
185 186
      double ckt2(cl->Mom().Abs2());
      kt2xmap[cl->Id()].first=kt2xmap[cl->Id()].second=HardScale(ampl);
187
      for (KT2X_Map::iterator kit(kt2xmap.begin());kit!=kt2xmap.end();++kit)
188
	if (kit->first!=cl->Id() && (kit->first&cl->Id()) &&
189
	    aset.find(kit->first)==aset.end()) {
190 191
	  kit->second.first=ckt2;
	  kit->second.second=ckt2;
192 193 194 195 196 197 198
	  aset.insert(kit->first);
	}
    }
    else if (ampl->Prev()==NULL) {
      kt2xmap[cl->Id()].first=kt2xmap[cl->Id()].second=HardScale(ampl); 
    }
    else {
199
      double ckt2max(HardScale(ampl));
200 201 202 203
      double ckt2min(std::numeric_limits<double>::max());
      for (KT2X_Map::iterator kit(kt2xmap.begin());kit!=kt2xmap.end();++kit)
	if (kit->first&cl->Id()) {
	  ckt2min=kit->second.first;
204 205
	  ckt2max=Max(ckt2max,kit->second.second);
	}
206 207
      kt2xmap[cl->Id()].first=ckt2min;
      kt2xmap[cl->Id()].second=ckt2max; 
208
      for (KT2X_Map::iterator kit(kt2xmap.begin());kit!=kt2xmap.end();++kit)
209 210
	if ((kit->first&cl->Id()) && aset.find(kit->first)==aset.end()) {
	  kit->second.first=ckt2min;
211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236
	  kit->second.second=ckt2max;
	}
    }
  }
}

void CS_Shower::GetKT2Min(Cluster_Amplitude *const ampl,KT2X_Map &kt2xmap)
{
  std::set<size_t> aset;
  Cluster_Amplitude *campl(ampl);
  while (campl->Next()) campl=campl->Next();
  GetKT2Min(campl,(1<<ampl->Legs().size())-1,kt2xmap,aset);
  std::vector<size_t> cns;
  double ckt2min(std::numeric_limits<double>::max()), ckt2max(0.0);
  for (KT2X_Map::iterator kit(kt2xmap.begin());kit!=kt2xmap.end();++kit)
    if (aset.find(kit->first)==aset.end()) {
      ckt2min=Min(ckt2min,kit->second.first);
      ckt2max=Max(ckt2max,kit->second.second);
      bool ins(true);
      for (size_t j(0);j<cns.size();++j)
	if (cns[j]&kit->first) {
	  ins=false;
	  break;
	}
      if (ins) cns.push_back(kit->first);
    }
237
  bool smin(ampl->Legs().size()-ampl->NIn()==campl->Leg(2)->NMax());
238 239 240 241
  for (KT2X_Map::iterator kit(kt2xmap.begin());kit!=kt2xmap.end();++kit)
    if (aset.find(kit->first)==aset.end()) {
      if (smin) kit->second.first=ckt2min;
      else kit->second.first=0.0;
Stefan Hoeche's avatar
Stefan Hoeche committed
242
      kit->second.second=ckt2max;
243 244 245
    }
  msg_Debugging()<<"k_{T,min} / k_{T,max} = {\n";
  for (KT2X_Map::const_iterator
246
  	 kit(kt2xmap.begin());kit!=kt2xmap.end();++kit)
247
    msg_Debugging()<<"  "<<ID(kit->first)
248 249
		  <<" -> "<<sqrt(kit->second.first)
		  <<" / "<<sqrt(kit->second.second)<<"\n";
250 251 252
  msg_Debugging()<<"}\n";
}

253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368
bool CS_Shower::PrepareShower(Cluster_Amplitude *const ampl,const bool & soft)
{
  if (soft) return PrepareShowerFromSoft(ampl);
  return PrepareStandardShower(ampl);
}


void CS_Shower::EstablishRelations(Parton * parton,Cluster_Leg * leg,
				   std::map<Cluster_Leg*,Parton*> & pmap) {
  int no = leg->NumberOfSpectators();
  if (no==0) return;

  /*
    msg_Out()<<"   "<<no<<" spectators for "
    <<"["<<leg->Col().m_i<<", "<<leg->Col().m_j<<"]";
    for (std::set<Cluster_Leg *>::iterator clit=leg->GetSpectators().begin();
    clit!=leg->GetSpectators().end();clit++) {
    msg_Out()<<" --> ["<<(*clit)->Col().m_i<<", "<<(*clit)->Col().m_j<<"]";
    }
    msg_Out()<<".\n";
  */

  size_t col1(parton->GetFlow(1)),col2(parton->GetFlow(2));
  bool connect;
  for (std::set<Cluster_Leg *>::iterator clit=leg->GetSpectators().begin();
       clit!=leg->GetSpectators().end();clit++) {
    connect = false;
    Parton * spect = pmap[(*clit)];
    if (col1!=0 && col1==spect->GetFlow(2)) {
      if ((parton->GetLeft() && parton->GetLeft()!=spect) ||
	  (spect->GetRight() && spect->GetRight()!=parton)) {
	msg_Error()<<"Error in "<<METHOD<<":\n"
		   <<"   have already left colour partner for \n"<<(*parton)
		   <<"  --> ["<<parton->GetLeft()->GetFlow(1)
		   <<", "<<parton->GetLeft()->GetFlow(2)<<"] while trying "
		   <<"["<<spect->GetFlow(1)<<", "<<spect->GetFlow(2)<<"].\n"
		   <<"   or right colour partner for \n"<<(*spect);
      }
      else {
	//msg_Tracking()<<"Set ["<<spect->GetFlow(1)<<", "
	//<<spect->GetFlow(2)<<"] "
	//	 <<" as left of ["<<parton->GetFlow(1)
	//	 <<", "<<parton->GetFlow(2)<<"].\n";
	parton->SetLeft(spect);
	spect->SetRight(parton);
	connect = true;
      }
    }
    if (col2!=0 && col2==spect->GetFlow(1)) {
      if ((parton->GetRight() && parton->GetRight()!=spect) ||
	  (spect->GetLeft() && spect->GetLeft()!=parton)) {
	msg_Error()<<"Error in "<<METHOD<<":\n"
		   <<"   have already right colour partner for \n"<<(*parton)
		   <<"  --> ["<<parton->GetRight()->GetFlow(1)
		   <<", "<<parton->GetRight()->GetFlow(2)<<"] while trying "
		   <<"["<<spect->GetLeft()->GetFlow(1)<<", "
		   <<spect->GetLeft()->GetFlow(2)<<"].\n"
		   <<"   or left colour partner for \n"<<(*spect);
      }
      else {
	//msg_Tracking()<<"Set ["<<spect->GetFlow(1)<<", "
	//<<spect->GetFlow(2)<<"] "
	//	 <<" as right of ["<<parton->GetFlow(1)
	//	 <<", "<<parton->GetFlow(2)<<"].\n";
	parton->SetRight(spect);
	spect->SetLeft(parton);
	connect = true;
      }
    }
    if (!connect) {
      if (spect->GetLeft()==parton) {
	msg_Tracking()<<"Set ["<<spect->GetFlow(1)<<", "
		      <<spect->GetFlow(2)<<"] "
		      <<" as help right of ["<<parton->GetFlow(1)
		      <<", "<<parton->GetFlow(2)<<"].\n";
	parton->SetRight(spect);
      }
      else if (spect->GetRight()==parton) {
	msg_Tracking()<<"Set ["<<spect->GetFlow(1)<<", "
		      <<spect->GetFlow(2)<<"] "
		      <<" as help left of ["<<parton->GetFlow(1)
		      <<", "<<parton->GetFlow(2)<<"].\n";
	parton->SetLeft(spect);
      }
      else {
	if (!parton->GetLeft()) {
	  msg_Tracking()<<"Set ["<<spect->GetFlow(1)<<", "
			<<spect->GetFlow(2)<<"] "
			<<" as help left of ["<<parton->GetFlow(1)
			<<", "<<parton->GetFlow(2)<<"].\n";
	  parton->SetLeft(spect);
	}
	else if (!parton->GetRight()) {
	  msg_Tracking()<<"Set ["<<spect->GetFlow(1)<<", "
			<<spect->GetFlow(2)<<"] "
			<<" as help right of ["<<parton->GetFlow(1)
			<<", "<<parton->GetFlow(2)<<"].\n";
	  parton->SetRight(spect);
	}
      }
    }
  }
  parton->SetStat(leg->Stat());
  msg_Tracking()<<(*parton);
}

bool CS_Shower::PrepareShowerFromSoft(Cluster_Amplitude *const ampl)
{
  CleanUp();
  msg_Tracking()<<"===============================================\n"
  	   <<METHOD<<"(ms = "<<ampl->MS()<<"):\n";
  //msg_Indent();
  p_rampl = ampl;
  p_ms    = ampl->MS();
  p_next->clear();
  //msg_Tracking()<<*ampl<<"\n";
369
  std::map<Parton*,Cluster_Leg*,partcomp> lmap;
370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430 431 432 433 434 435 436 437 438
  std::map<Cluster_Leg*,Parton*> pmap;

  Parton      * parton;
  Cluster_Leg * leg;
  Singlet *singlet(new Singlet());
  singlet->SetMS(p_ms);
  for (size_t i(0);i<ampl->Legs().size();++i) {
    leg = ampl->Leg(i);
    if (leg->Flav().IsHadron() && 
	leg->Flav().Kfcode()!=kf_pomeron &&
	leg->Flav().Kfcode()!=kf_reggeon && 
	leg->Id()&((1<<ampl->NIn())-1)) continue;
    bool is(leg->Id()&((1<<ampl->NIn())-1));
    Particle p(1,is?leg->Flav().Bar():leg->Flav(),is?-leg->Mom():leg->Mom());
    if (is) {
      p.SetFlow(2,leg->Col().m_i);
      p.SetFlow(1,leg->Col().m_j);
    }
    else {
      p.SetFlow(1,leg->Col().m_i);
      p.SetFlow(2,leg->Col().m_j);
    }
    parton       = new Parton(&p,is?pst::IS:pst::FS);
    pmap[leg]    = parton;
    lmap[parton] = leg;
    parton->SetRFlow();
    parton->SetKin(p_shower->KinScheme());
    parton->SetId(leg->Id());
    if (is) {
      if (Vec3D(p.Momentum())*Vec3D(rpa->gen.PBeam(0))>0.) {
	parton->SetBeam(0);
      }
      else { 
	parton->SetBeam(1);
      }
    }
    double kt2max(leg->KTMax()),kt2veto(leg->KTVeto()),kt2start(leg->KTStart());
    parton->SetStart(kt2start);   // start scale of shower
    parton->SetKtMax(kt2max);    // no jet veto below ktmax
    parton->SetKtPrev(kt2veto);   // upper kt - acts as veto
    parton->SetKtNext(0.0);       // lower kt - set to 0.
    parton->SetVeto(kt2veto);     // irrelevant 
    parton->SetConnected(leg->Connected());
    parton->SetMass2(p_ms->Mass2(leg->Flav()));
    singlet->push_back(parton);
    parton->SetSing(singlet);
    
    //msg_Out()<<"Add parton "<<parton->Id()<<",conn="<<leg->Connected()<<") "
    //	     <<parton->GetFlavour()<<" "<<parton->Momentum()
    //	     <<" ["<<parton->GetFlow(1)<<", "<<parton->GetFlow(2)<<"]: "
    //	     <<"start = "<<parton->KtStart()<<" --> "
    //	     <<parton->KtPrev()<<"\n";
  }
  for (std::map<Parton*,Cluster_Leg*>::iterator pit=lmap.begin();
       pit!=lmap.end();pit++) {
    parton = pit->first;
    leg    = pit->second;
    EstablishRelations(parton,leg,pmap);
  }
  p_shower->SetMS(p_ms);

  pmap.clear();
  lmap.clear();
  //msg_Tracking()<<(*singlet)<<"\n";
  m_allsinglets.push_back(singlet);
  return true;
}

bool CS_Shower::PrepareStandardShower(Cluster_Amplitude *const ampl)
439 440 441 442 443 444
{
  CleanUp();
  msg_Debugging()<<METHOD<<"(): {\n";
  msg_Indent();
  p_rampl=ampl;
  p_ms=ampl->MS();
445 446
  KT2X_Map kt2xmap;
  GetKT2Min(ampl,kt2xmap);
447 448 449 450 451 452 453 454 455 456
  p_next->clear();
  All_Singlets allsinglets;
  std::map<size_t,Parton*> kmap;
  std::map<Parton*,Cluster_Leg*> almap;
  std::map<Cluster_Leg*,Parton*> apmap;
  for (Cluster_Amplitude *campl(ampl);campl;campl=campl->Next()) {
    msg_Debugging()<<*campl<<"\n";
    Parton *split(NULL);
    std::map<Parton*,Cluster_Leg*> lmap;
    std::map<Cluster_Leg*,Parton*> pmap;
457
    Singlet *sing(TranslateAmplitude(campl,pmap,lmap,kt2xmap));
458 459 460 461 462 463 464 465 466 467 468 469 470 471 472 473 474 475 476 477 478 479 480 481 482 483 484 485 486 487 488 489 490 491 492 493 494 495 496 497 498 499 500 501
    allsinglets.push_back(sing);
    for (size_t i(0);i<campl->Legs().size();++i) {
      Cluster_Leg *cl(campl->Leg(i));
      if (pmap.find(cl)==pmap.end()) continue;
      Parton *k(pmap[cl]);
      k->SetId(cl->Id());
      almap[apmap[cl]=k]=cl;
      std::map<size_t,Parton*>::iterator cit(kmap.find(cl->Id()));
      if (cit!=kmap.end()) {
	if (k->GetNext()!=NULL) 
	  THROW(fatal_error,"Invalid tree structure");
	k->SetNext(cit->second);
	cit->second->SetPrev(k);
	k->SetStat(1);
	kmap[cl->Id()]=k;
	continue;
      }
      for (std::map<size_t,Parton*>::iterator 
	     kit(kmap.begin());kit!=kmap.end();)
	if ((kit->first&cl->Id())!=kit->first) {
	  ++kit;
	}
	else {
	  if (k->GetSing()->GetLeft()==NULL) k->GetSing()->SetLeft(kit->second);
	  else if (k->GetSing()->GetRight()==NULL) {
	    if (kit->second->GetType()==pst::IS &&
		k->GetSing()->GetLeft()->GetType()==pst::FS) {
	      k->GetSing()->SetRight(k->GetSing()->GetLeft());
	      k->GetSing()->SetLeft(kit->second);
	    }
	    else {
	      k->GetSing()->SetRight(kit->second);
	    }
	  }
	  else THROW(fatal_error,"Invalid tree structure");
	  kit->second->SetPrev(k);
	  k->SetStat(1);
	  kmap.erase(kit);
	  kit=kmap.begin();
	  split=k;
	}
      kmap[cl->Id()]=k;
    }
    if (sing->GetSpec()) {
502 503
      split->SetOldMomentum(split->Momentum());
      Vec4D fixspec(sing->GetSpec()->Momentum());
Stefan Hoeche's avatar
Stefan Hoeche committed
504
      fixspec[0]=fixspec[0]<0.0?-fixspec.PSpat():fixspec.PSpat();
505
      split->SetFixSpec(fixspec);
506 507 508 509 510 511 512
      sing->SetSpec(sing->GetSpec()->GetNext());
      if (split==NULL) THROW(fatal_error,"Invalid tree structure");
      sing->SetSplit(split);
      Parton *l(sing->GetLeft()), *r(sing->GetRight()), *s(sing->GetSpec());
      almap[l]->SetMom(almap[l]->Id()&3?-l->Momentum():l->Momentum());
      almap[r]->SetMom(almap[r]->Id()&3?-r->Momentum():r->Momentum());
      almap[s]->SetMom(almap[s]->Id()&3?-s->Momentum():s->Momentum());
Stefan Hoeche's avatar
Stefan Hoeche committed
513
      split->SetKin(campl->Kin());
514
      CS_Parameters cp(p_cluster->KT2
515 516
		       (campl->Prev(),almap[l],almap[r],almap[s],
			split->GetType()==pst::FS?split->GetFlavour():
517
			split->GetFlavour().Bar(),p_ms,split->Kin()));
518 519 520 521 522
      l->SetTest(cp.m_kt2,cp.m_z,cp.m_y,cp.m_phi);
      l->SetStart(cp.m_kt2);
      r->SetStart(cp.m_kt2);
      msg_Debugging()<<"Set reco params: kt = "<<sqrt(cp.m_kt2)<<", z = "
		     <<cp.m_z<<", y = "<<cp.m_y<<", phi = "<<cp.m_phi
523
		     <<", mode = "<<cp.m_mode<<", scheme = "<<split->Kin()<<"\n";
524
      sing->SetAll(p_next);
525
      if (m_recocheck&1) {
526 527
      std::cout.precision(12);
      Vec4D oldl(l->Momentum()), oldr(r->Momentum()), olds(s->Momentum());
528 529
      Vec4D oldfl(l->FixSpec()), oldfr(r->FixSpec()), oldfs(s->FixSpec());
      sing->BoostBackAllFS(l,r,s,split,split->GetFlavour(),cp.m_mode|4);
530
      p_shower->ReconstructDaughters(sing,1);
531 532 533 534
      almap[l]->SetMom(almap[l]->Id()&3?-l->Momentum():l->Momentum());
      almap[r]->SetMom(almap[r]->Id()&3?-r->Momentum():r->Momentum());
      almap[s]->SetMom(almap[s]->Id()&3?-s->Momentum():s->Momentum());
      CS_Parameters ncp(p_cluster->KT2
535 536
			(campl->Prev(),almap[l],almap[r],almap[s],
			 split->GetType()==pst::FS?split->GetFlavour():
537
			 split->GetFlavour().Bar(),p_ms,split->Kin()));
538
      msg_Debugging()<<"New reco params: kt = "<<sqrt(ncp.m_kt2)<<", z = "
539 540
		     <<ncp.m_z<<", y = "<<ncp.m_y<<", phi = "<<ncp.m_phi
		     <<", kin = "<<ncp.m_kin<<"\n";
541
      msg_Debugging()<<"            vs.: kt = "<<sqrt(cp.m_kt2)<<", z = "
542 543
		     <<cp.m_z<<", y = "<<cp.m_y<<", phi = "<<cp.m_phi
		     <<", kin = "<<cp.m_kin<<"\n";
544 545 546 547 548 549 550 551
      if (!IsEqual(ncp.m_kt2,cp.m_kt2,1.0e-6) || 
	  !IsEqual(ncp.m_z,cp.m_z,1.0e-6) || 
	  !IsEqual(ncp.m_y,cp.m_y,1.0e-6) || 
	  !IsEqual(ncp.m_phi,cp.m_phi,1.0e-6) ||
	  !IsEqual(oldl,l->Momentum(),1.0e-6) || 
	  !IsEqual(oldr,r->Momentum(),1.0e-6) || 
	  !IsEqual(olds,s->Momentum(),1.0e-6)) {
	msg_Error()<<"\nFaulty reco params: kt = "<<sqrt(ncp.m_kt2)<<", z = "
552
		   <<ncp.m_z<<", y = "<<ncp.m_y<<", phi = "<<ncp.m_phi
553
		   <<", mode = "<<ncp.m_mode<<", scheme = "<<ncp.m_kin<<"\n";
554
	msg_Error()<<"               vs.: kt = "<<sqrt(cp.m_kt2)<<", z = "
555
		   <<cp.m_z<<", y = "<<cp.m_y<<", phi = "<<cp.m_phi
556
		   <<", mode = "<<cp.m_mode<<", scheme = "<<cp.m_kin<<"\n";
557 558 559
	msg_Error()<<"  "<<oldl<<" "<<oldr<<" "<<olds<<"\n";
	msg_Error()<<"  "<<l->Momentum()<<" "<<r->Momentum()
		   <<" "<<s->Momentum()<<"\n";
560
	if (m_recocheck&2) abort();
561
      }
562 563 564
      l->SetMomentum(oldl);
      r->SetMomentum(oldr);
      s->SetMomentum(olds);
565 566 567 568 569 570
      l->SetOldMomentum(oldl);
      r->SetOldMomentum(oldr);
      s->SetOldMomentum(olds);
      l->SetFixSpec(oldfl);
      r->SetFixSpec(oldfr);
      s->SetFixSpec(oldfs);
571
      }
572
      sing->BoostBackAllFS(l,r,s,split,split->GetFlavour(),cp.m_mode|4);
573
    }
574 575
    double kt2prev(campl->Next()?campl->KT2():kt2xmap[1].second);
    double kt2next(campl->Prev()?campl->Prev()->KT2():0.0);
576 577 578 579
    for (size_t i(0);i<campl->Legs().size();++i) {
      std::map<Cluster_Leg*,Parton*>::const_iterator 
	pit(apmap.find(campl->Leg(i)));
      if (pit!=apmap.end()) {
580
	pit->second->SetKtPrev(Min(kt2prev,m_kt2prev));
581 582 583 584 585 586 587 588 589 590 591 592 593
	pit->second->SetKtNext(kt2next);
      }
    }
    if (ampl->NIn()==1 && ampl->Leg(0)->Flav().IsHadron()) break;
    p_next->push_back(sing);
  }
  p_next->clear();
  for (All_Singlets::reverse_iterator
	 asit(allsinglets.rbegin());asit!=allsinglets.rend();++asit) {
    m_allsinglets.push_back(*asit);
    p_next->push_back(*asit);
  }
  msg_Debugging()<<"\nSinglet lists:\n\n";
594 595
  Cluster_Amplitude *campl(ampl);
  while (campl->Next()) campl=campl->Next();
596 597 598 599 600 601 602
  for (All_Singlets::const_iterator 
	 sit(m_allsinglets.begin());sit!=m_allsinglets.end();++sit) {
      for (Singlet::const_iterator 
	     pit((*sit)->begin());pit!=(*sit)->end();++pit) {
	if ((*pit)->GetPrev()) {
	  if ((*pit)->GetPrev()->GetNext()==*pit) 
	    (*pit)->SetStart((*pit)->GetPrev()->KtStart());
603
	  (*pit)->SetKtPrev(Min((*pit)->GetPrev()->KtNext(),m_kt2prev));
604 605 606
	}
      }
      (*sit)->SetJF(ampl->JF<PHASIC::Jet_Finder>());
607
      (*sit)->SetDecays(campl->Decays());
608 609 610 611 612 613 614 615 616 617 618 619
      (*sit)->SetAll(p_next);
      msg_Debugging()<<**sit;
    msg_Debugging()<<"\n";
  }
  msg_Debugging()<<"}\n";
  p_shower->SetMS(p_ms);
  return true;
}

Singlet *CS_Shower::TranslateAmplitude
(Cluster_Amplitude *const ampl,
 std::map<Cluster_Leg*,Parton*> &pmap,std::map<Parton*,Cluster_Leg*> &lmap,
620
 const KT2X_Map &kt2xmap)
621
{
622
  PHASIC::Jet_Finder *jf(ampl->JF<PHASIC::Jet_Finder>());
623
  double ktveto2(jf?jf->Ycut()*sqr(rpa->gen.Ecms()):4.0*ampl->MuR2());
624 625
  Singlet *singlet(new Singlet());
  singlet->SetMS(p_ms);
626
  singlet->SetNLO(ampl->NLO());
627 628 629 630
  for (size_t i(0);i<ampl->Legs().size();++i) {
    Cluster_Leg *cl(ampl->Leg(i));
    if (cl->Flav().IsHadron() && cl->Id()&((1<<ampl->NIn())-1)) continue;
    bool is(cl->Id()&((1<<ampl->NIn())-1));
631
    Particle p(1,is?cl->Flav().Bar():cl->Flav(),is?-cl->Mom():cl->Mom());
632
    if (is) {
633 634
      p.SetFlow(2,cl->Col().m_i);
      p.SetFlow(1,cl->Col().m_j);
635 636
    }
    else {
637 638
      p.SetFlow(1,cl->Col().m_i);
      p.SetFlow(2,cl->Col().m_j);
639
    }
640
    Parton *parton(new Parton(&p,is?pst::IS:pst::FS));
641 642 643
    parton->SetMass2(p_ms->Mass2(p.Flav()));
    if (!is && parton->Mass2()>10.0 && !p.Flav().Strong())
      parton->SetMass2(parton->Momentum().Abs2());
644 645 646
    pmap[cl]=parton;
    lmap[parton]=cl;
    parton->SetRFlow();
647
    parton->SetKin(p_shower->KinScheme());
648
    if (ampl->NLO()==1) parton->SetTMin(std::numeric_limits<double>::max());
Stefan Hoeche's avatar
Stefan Hoeche committed
649
    if (is) parton->SetBeam(i);
650
    KT2X_Map::const_iterator xit(kt2xmap.find(cl->Id()));
651
    parton->SetStart(xit->second.second);
652
    parton->SetKtMax(xit->second.first);
653 654 655 656 657 658 659
    parton->SetVeto(ktveto2);
    singlet->push_back(parton);
    parton->SetSing(singlet);
  }
  for (Singlet::const_iterator sit(singlet->begin());
       sit!=singlet->end();++sit) {
    int flow[2]={(*sit)->GetFlow(1),(*sit)->GetFlow(2)};
660 661 662
    if (!(*sit)->GetFlavour().Strong()) continue;
    if ((flow[0]==0 || (*sit)->GetLeft()!=NULL) &&
	(flow[1]==0 || (*sit)->GetRight()!=NULL)) continue;
663 664 665 666 667 668
    if (flow[0]) {
      for (Singlet::const_iterator tit(singlet->begin());
	   tit!=singlet->end();++tit)
	if (tit!=sit && (*tit)->GetFlow(2)==flow[0]) {
	  (*sit)->SetLeft(*tit);
	  (*tit)->SetRight(*sit);
669 670 671
	  break;
	}
    }
672 673 674 675 676 677 678
    if (flow[1]) {
      for (Singlet::const_iterator tit(singlet->begin());
	   tit!=singlet->end();++tit)
	if (tit!=sit && (*tit)->GetFlow(1)==flow[1]) {
	  (*sit)->SetRight(*tit);
	  (*tit)->SetLeft(*sit);
	  break;
679 680
	}
    }
681
    if ((flow[0] && (*sit)->GetLeft()==NULL) ||
682
	(flow[1] && (*sit)->GetRight()==NULL))
683
      THROW(fatal_error,"Missing colour partner");
684 685 686 687 688 689 690 691 692 693 694 695
  }
  for (size_t i(0);i<ampl->Legs().size();++i)
    if (ampl->Leg(i)->K()) {
      Singlet *sing(pmap[ampl->Leg(i)]->GetSing());
      sing->SetSpec(pmap[ampl->IdLeg(ampl->Leg(i)->K())]);
      break;
    }
  return singlet;
}

double CS_Shower::HardScale(const Cluster_Amplitude *const ampl)
{
696 697
  if (ampl->Next()) {
    Cluster_Amplitude *next(ampl->Next());
698
    if (next->OrderQCD()<ampl->OrderQCD()) return ampl->KT2();
699
    return HardScale(next);
700
  }
701
  return ampl->KT2();
702 703
}

704 705 706 707 708 709 710 711 712 713
double CS_Shower::CplFac(const ATOOLS::Flavour &fli,const ATOOLS::Flavour &flj,
                         const ATOOLS::Flavour &flk,const int type,
			 const int cpl,const double &mu2) const
{
  cstp::code stp((type&1)?
		 (type&2)?cstp::II:cstp::IF:
		 (type&2)?cstp::FI:cstp::FF);
  return p_shower->GetSudakov()->CplFac(fli, flj, flk,stp,cpl,mu2);
}

Stefan Hoeche's avatar
Stefan Hoeche committed
714 715
bool CS_Shower::JetVeto(ATOOLS::Cluster_Amplitude *const ampl,
			const int mode)
716 717 718 719
{
  DEBUG_FUNC("");
  msg_Debugging()<<*ampl<<"\n";
  PHASIC::Jet_Finder *jf(ampl->JF<PHASIC::Jet_Finder>());
720
  double q2cut(jf->Ycut()*sqr(rpa->gen.Ecms()));
721
  bool his(false);
722 723 724 725 726 727 728
  size_t noem(0), nospec(0);
  for (size_t i(0);i<ampl->Decays().size();++i) {
    noem|=ampl->Decays()[i]->m_id;
    if (!ampl->Decays()[i]->m_fl.Strong())
      nospec|=ampl->Decays()[i]->m_id;
  }
  msg_Debugging()<<"noem = "<<ID(noem)<<", nospec = "<<ID(nospec)<<"\n";
Stefan Hoeche's avatar
Stefan Hoeche committed
729 730 731
  double q2min(std::numeric_limits<double>::max());
  size_t imin(0), jmin(0), kmin(0);
  Flavour mofl;
732 733
  for (size_t i(0);i<ampl->Legs().size();++i) {
    Cluster_Leg *li(ampl->Leg(i));
734
    if (li->Id()&noem) continue;
735 736 737 738
    Flavour fi(i<ampl->NIn()?li->Flav().Bar():li->Flav());
    if (i<ampl->NIn() && fi.Resummed()) his=true;
    for (size_t j(Max(i+1,ampl->NIn()));j<ampl->Legs().size();++j) {
      Cluster_Leg *lj(ampl->Leg(j));
739
      if (lj->Id()&noem) continue;
740 741 742 743
      Flavour fj(j<ampl->NIn()?lj->Flav().Bar():lj->Flav());
      for (size_t k(0);k<ampl->Legs().size();++k) {
	if (k==i || k==j) continue;
	Cluster_Leg *lk(ampl->Leg(k));
744
	if (lk->Id()&nospec) continue;
745 746 747 748 749 750 751 752 753
	Flavour fk(k<ampl->NIn()?lk->Flav().Bar():lk->Flav());
	cstp::code et((i<ampl->NIn()||j<ampl->NIn())?
		      (k<ampl->NIn()?cstp::II:cstp::IF):
		      (k<ampl->NIn()?cstp::FI:cstp::FF));
	if (p_shower->GetSudakov()->HasKernel(fi,fj,fk,et)) {
	  double q2ijk(PDF::Qij2(li->Mom(),lj->Mom(),lk->Mom(),
				 li->Flav(),lj->Flav(),jf->DR()));
 	  msg_Debugging()<<"Q_{"<<ID(li->Id())<<ID(lj->Id())
			 <<","<<ID(lk->Id())<<"} = "<<sqrt(q2ijk)<<"\n";
Stefan Hoeche's avatar
Stefan Hoeche committed
754 755 756 757 758 759 760 761 762 763 764 765 766 767
	  if (mode==0) {
	    if (q2ijk<q2cut) return false;
	  }
	  else {
	    if (q2ijk<q2min) {
	      q2min=q2ijk;
	      mofl=Flavour(kf_gluon);
	      if (li->Flav().IsGluon()) mofl=lj->Flav();
	      if (lj->Flav().IsGluon()) mofl=li->Flav();
	      imin=i;
	      jmin=j;
	      kmin=k;
	    }
	  }
768 769 770 771 772 773 774 775
	}
	else {
	  msg_Debugging()<<"No kernel for "<<fi<<" "<<fj
			 <<" <-> "<<fk<<" ("<<et<<")\n";
	}
      }
    }
  }
776
  if (mode!=0 && imin!=jmin) {
Stefan Hoeche's avatar
Stefan Hoeche committed
777 778 779 780 781 782 783 784 785 786 787 788 789 790 791 792 793 794 795 796
    Vec4D_Vector p=p_cluster->Combine(*ampl,imin,jmin,kmin,mofl,ampl->MS(),1);
    Cluster_Amplitude *bampl(Cluster_Amplitude::New());
    bampl->SetNIn(ampl->NIn());
    bampl->SetJF(ampl->JF<void>());
    for (int i(0), j(0);i<ampl->Legs().size();++i) {
      if (i==jmin) continue;
      if (i==imin) {
	bampl->CreateLeg(p[j],mofl,ampl->Leg(i)->Col());
	bampl->Legs().back()->SetId(ampl->Leg(imin)->Id()|ampl->Leg(jmin)->Id());
	bampl->Legs().back()->SetK(ampl->Leg(kmin)->Id());	
      }
      else {
	bampl->CreateLeg(p[j],ampl->Leg(i)->Flav(),ampl->Leg(i)->Col());
      }
      ++j;
    }
    bool res=JetVeto(bampl,0);
    bampl->Delete();
    return res;
  }
797 798 799 800
  if (his) {
    for (size_t i(ampl->NIn());i<ampl->Legs().size();++i)
      if (ampl->Leg(i)->Flav().Resummed())
	if (ampl->Leg(i)->Mom().PPerp2()<
801
	    p_shower->GetSudakov()->ISPT2Min()) {
802 803 804 805 806
 	  msg_Debugging()<<"p_T_{"<<ID(ampl->Leg(i)->Id())<<"} = "
			 <<ampl->Leg(i)->Mom().PPerp()<<"\n";
	  return false;
	}
  }
807 808
  msg_Debugging()<<"--- Jet veto ---\n";
  return true;
809 810
}

811 812 813 814 815 816 817 818 819 820 821 822 823 824 825
namespace PDF {

  DECLARE_GETTER(CSS_Getter,"CSS",Shower_Base,Shower_Key);

  Shower_Base *CSS_Getter::operator()(const Shower_Key &key) const
  {
    return new CS_Shower(key.p_isr,key.p_model,key.p_read);
  }

  void CSS_Getter::PrintInfo(std::ostream &str,const size_t width) const
  { 
    str<<"The CSS shower"; 
  }

}