#include "SHRiMPS/Main/Cluster_Algorithm.H" #include "SHRiMPS/Event_Generation/Inelastic_Event_Generator.H" #include "SHRiMPS/Cross_Sections/Sigma_Inelastic.H" #include "ATOOLS/Phys/Particle.H" #include "ATOOLS/Org/Run_Parameter.H" #include "ATOOLS/Org/Message.H" #include "ATOOLS/Math/Random.H" #include "ATOOLS/Math/Histogram.H" using namespace SHRIMPS; using namespace ATOOLS; using namespace std; Inelastic_Event_Generator::Inelastic_Event_Generator() { Initialise(); } Inelastic_Event_Generator::~Inelastic_Event_Generator() { Reset(); while (!m_Bgrids.empty()) { delete m_Bgrids.begin()->second; m_Bgrids.erase(m_Bgrids.begin()); } m_Bgrids.clear(); } void Inelastic_Event_Generator::Initialise() { m_sigma = 0.; Sigma_Inelastic sigma; list * eikonals(MBpars.GetEikonals()); for (list::iterator eikonal=eikonals->begin(); eikonal!=eikonals->end(); eikonal++) { m_Bgrids[(*eikonal)] = sigma.FillBGrid((*eikonal)); m_sigma += m_xsecs[(*eikonal)] = m_Bgrids[(*eikonal)]->back()*rpa->Picobarn(); } msg_Info()< m_Nladders) m_done = true; return 1; } bool Inelastic_Event_Generator::InitInelasticEvent(Blob_List * blobs) { if (m_init) return true; Blob * blob(blobs->FindFirst(btp::Soft_Collision)); if (!(blob && blob->Status()==blob_status::needs_minBias) || !SetUpEvent()) return false; m_init = true; m_laddergenerator.InitCollision(p_eikonal,m_B,m_Nladders); return true; } bool Inelastic_Event_Generator::SetUpEvent() { p_eikonal = 0; m_B = -1; m_Nladders = m_Nprim = m_Ngen = 0; size_t trials=0; do { SelectEikonal(); SelectB(); m_Nladders = ran->Poissonian(2.*(*p_eikonal)(m_B)); } while (trials++<1000 && m_Nladders<1); return (trials<1000); } bool Inelastic_Event_Generator::SelectEikonal() { p_eikonal = 0; while (p_eikonal==NULL) { double disc = ran->Get()*m_sigma; for (std::map::iterator eikiter=m_xsecs.begin(); eikiter!=m_xsecs.end();eikiter++) { disc-=eikiter->second; if (disc<=1.e-12) { p_eikonal = eikiter->first; break; } } } return (p_eikonal!=0); } bool Inelastic_Event_Generator::SelectB() { if (p_eikonal==0) { msg_Error()<<"Error in "< * grid = m_Bgrids[p_eikonal]; double deltaB(p_eikonal->DeltaB()); m_B = -1.; do { double random = ran->Get()*(*grid)[grid->size()-1]; size_t bin(0); while (binsize()-1 && (random-(*grid)[bin]>=0)) bin++; if (bin>=grid->size()) continue; double inthigh((*grid)[bin]), intlow((*grid)[bin-1]); double Bhigh(bin*deltaB), Blow((bin-1)*deltaB); m_B = (Blow*(random-intlow)+Bhigh*(inthigh-random))/(inthigh-intlow); } while (m_B<0.); return (m_B>=0.); } int Inelastic_Event_Generator::AddScatter(Blob_List * blobs) { Blob * blob(blobs->FindFirst(btp::Soft_Collision)); if (!(blob && blob->Status()==blob_status::needs_minBias)) return -1; if (m_Ngen>0) blob = CreateBlob(); SetBlobType(blob); if (!m_laddergenerator.MakePrimaryLadder(blob,m_Ngen==0)) { delete blob; return 0; } if (m_Ngen!=0) blobs->push_back(blob); p_cluster->SetTMax(4.*m_laddergenerator.TMax()); p_cluster->SetMinKT2(64.); //msg_Out()<SetExternalWeight(m_sigma); return 1; } Blob * Inelastic_Event_Generator::CreateBlob() { Blob * blob = new Blob(); blob->SetId(); return blob; } void Inelastic_Event_Generator::SetBlobType(Blob * blob) { blob->SetType(btp::Hard_Collision); blob->SetTypeSpec("MinBias"); blob->SetStatus(blob_status::needs_showers|blob_status::needs_beams); } void Inelastic_Event_Generator::Test(const std::string & dirname) { TestSelectB(dirname); TestNumberOfLadders(dirname); m_laddergenerator.Test(dirname); } void Inelastic_Event_Generator::TestSelectB(const std::string & dirname) { Histogram histo_in(0,0.,MBpars.GetEikonalParameters().bmax,100); Histogram histo_out(0,0.,MBpars.GetEikonalParameters().bmax,100); SelectEikonal(); Sigma_Inelastic sigma; sigma.SetEikonal(p_eikonal); double deltaB(p_eikonal->DeltaB()), B(deltaB/2.), val; do { val = 2.*M_PI*B*sigma.GetValue(B); histo_in.Insert(B,val); B += deltaB; } while (BPoissonian(arg))+0.5); ostringstream converter; converter << i; string istr(converter.str()); string name(string("NLadders_B_")+istr+string("_GeV.dat")); histo.Finalize(); histo.Output(dirname+"/"+name); } } // bool Inelastic_Event_Generator::DressShowerBlob(Blob * blob) { // msg_Error()<Id()<<", "<Type()<<"].\n"; // Abort(); // } // void Inelastic_Event_Generator::FixEikonalAndImpact() { // p_eikonal = m_sigma.SelectEikonal(); // m_B = m_sigma.SelectB(); // } // int Inelastic_Event_Generator:: // AddScatter(Blob_List * blobs,const double & xsec) // { // msg_Tracking()<GetCompensatorBlob())) { // m_connectblobs++; // return -1; // } // if (m_Nprim>0 && m_Nprim<=m_Nladders) { // m_rescatterhandler.UpdateCollision(blobs); // if (m_rescatterhandler.SelectRescatter(part1,part2)) { // m_Nsec++; // //m_laddergenerator.SetLadderGeneration(1+m_Nsec); // p_ladder = m_laddergenerator(part1,part2,true); // if (!p_ladder) return -1; // p_beams->SetInitials(part1,part2); // if (!p_beams->UpdateColours(p_ladder,false)) { // m_updatecols++; // return -1; // } // if (!CreateBlob(blobs,1.)) { // m_laddercols++; // return -1; // } // m_rescatterhandler.Map(part1,p_ladder->GetIn1()->GetParticle()); // m_rescatterhandler.Map(part2,p_ladder->GetIn2()->GetParticle()); // const double b1(m_laddergenerator.B1()), b2(m_laddergenerator.B2()); // if (m_analyse) { // m_histograms[string("B1_all")]->Insert(b1); // m_histograms[string("B2_all")]->Insert(b2); // } // m_Ngen++; // return 1; // } // } // if (m_NprimNextIS(part1,part2)) return -1; // m_Nsec = 0; // //m_laddergenerator.SetLadderGeneration(1); // p_ladder = m_laddergenerator(part1,part2,false,m_Nprim==0,m_weighted); // if (!p_ladder) return -1; // if (!p_beams->UpdateColours(p_ladder,m_Nprim+1==m_Nladders)) { // m_updatecols++; // return -1; // } // m_rescatterhandler.ResetRescatter(m_Nprim==0); // if (!CreateBlob(blobs,m_weighted?xsec*p_ladder->Weight():1.)) { // m_laddercols++; // return -1; // } // const double b1(m_laddergenerator.B1()), b2(m_laddergenerator.B2()); // if (m_analyse) { // m_histograms[string("B1_prim")]->Insert(b1); // m_histograms[string("B1_all")]->Insert(b1); // m_histograms[string("B2_prim")]->Insert(b2); // m_histograms[string("B2_all")]->Insert(b2); // } // m_Nprim++; m_Ngen++; // return 1; // } // else { // msg_Tracking() // <<"##################################################################\n" // <<" Out of event with "<Insert(m_Nprim); // m_histograms[string("N_ladder_true")]->Insert(m_Ngen); // m_histograms[string("N_ladder_sec")]->Insert(m_Ngen-m_Nprim); // m_histograms[string("N_ladder1_B")]->Insert(m_B,m_Nprim); // m_histograms[string("N_ladder_all_B")]->Insert(m_B,m_Ngen); // } // return 0; // } // msg_Tracking()<<" @@@ undefined, Nprim = "<Position()*rpa->hBar()*rpa->c()); // Blob * blob(blobs->FindFirst(btp::Soft_Collision)); // //msg_Out()<size()<<", blob = "<size()<=2) msg_Out()<<(*blobs)<<"\n"; // bool add(true); // if (blob && blob->Status()==blob_status::needs_minBias) { // if (blob->NInP()>0) { // msg_Error()<<"Error in "<DeleteInParticles(); // } // if (blob->NOutP()>0) { // msg_Error()<<"Error in "<DeleteOutParticles(); // } // blob->UnsetStatus(blob_status::needs_minBias); // add = false; // } // else { // blob = new Blob(); // blob->SetId(); // blobs->push_back(blob); // } // blob->SetType(btp::Hard_Collision); // if (m_isUE) blob->SetTypeSpec("UnderlyingEvent"); // else blob->SetTypeSpec("MinBias"); // blob->SetStatus(blob_status::needs_showers); // blob->SetPosition(pos); // Blob_Data_Base *winfo((*blob)["Weight"]); // if (!winfo) blob->AddData("Weight",new ATOOLS::Blob_Data(1.)); // Blob_Data_Base *wninfo((*blob)["Weight_Norm"]); // if (!wninfo) blob->AddData("Weight_Norm",new ATOOLS::Blob_Data(1.)); // Blob_Data_Base *tinfo((*blob)["Trials"]); // if (!tinfo) blob->AddData("Trials",new ATOOLS::Blob_Data(1.)); // Particle * part; // for (LadderMap::iterator liter=p_ladder->GetEmissionsBegin(); // liter!=p_ladder->GetEmissionsEnd();liter++) { // part = liter->second.GetParticle(); // blob->AddToOutParticles(part); // } // double shat((p_ladder->GetIn1()->GetParticle()->Momentum() // +p_ladder->GetIn2()->GetParticle()->Momentum()).Abs2()); // m_rescatterhandler.FillInitialStateIntoBlob(blob,p_ladder); // blob->SetCMS(); // if (blob->CheckMomentumConservation().Abs2()/shat>1.e-6 || // blob->CheckMomentumConservation()[0]/sqrt(shat)>1.e-3 || // blob->CheckMomentumConservation()[3]/sqrt(shat)>1.e-3) { // msg_Error()<<"Problem in "<Id()<<") seems fishy: " // <CheckMomentumConservation()<<".\n" // <<(*blob)<<"\n"<<(*p_ladder)<<"\n"; // } // if (!blob->CheckColour()) { // msg_Error()<<"Problem in "<Id()<<") seems fishy: " // <<"Bad colour configuration.\n" // <<(*blob)<<"\n"<<(*p_ladder)<<"\n"; // return false; // } // //msg_Out()< Hand over to shower now.\n" // // <<"===============================================\n"; // return true; // } // double Inelastic_Event_Generator::Smin() const { // double smin(m_luminosity.Smin()*m_Nladders); // if (!p_ladder) return smin; // // smin *= m_kt2fac; // if (p_ladder->IsHardDiffractive() && p_ladder->Size()==2) smin *= m_difffac; // return smin; // } // bool Inelastic_Event_Generator::IsLastRescatter() const { // if (!p_ladder) return false; // return p_ladder->IsRescatter(); // } ///////////////////////////////////////////////////////////////////////// //void Inelastic_Event_Generator:: //TestNumberOfLadders(Omega_ik * eikonal,const double & B){ // int nval(10000); // double mcmean(0.),anamean(0.),a,c,value; // double Y(p_sigma->Y()); // double kappa(eikonal->Kappa_i()); // double beta0(eikonal->FF1()->Beta0()); // double Lambda2(eikonal->Lambda2()); // double Delta(eikonal->Delta()); // a = Lambda2/(8.*(1.+kappa)); // c = sqr(beta0)*Lambda2*(1.+kappa)*exp(2.*Delta*Y)/(8.*M_PI); // anamean = c*exp(-a*sqr(B)); // for(int i=0; iPoissonian((*eikonal)(B))); // mcmean += value/nval; // } // msg_Tracking()<<"In "<< METHOD <<" mean number of ladders: "<