NcCollider.cxx

Go to the documentation of this file.

 
#include "NcCollider.h"
#include "Riostream.h"
 
ClassImp(NcCollider); // Class implementation to enable ROOT I/O
 
NcCollider::NcCollider() : TPythia6()
{
 
 fVertexmode=0;    // No vertex structure creation
 fResolution=1e-7; // Standard resolution is 0.1 micron
 fRunnum=0;
 fEventnum=0;
 fPrintfreq=1;
 fUserctrl=0;  // Automatic optimisation of some MC parameters 
 fElastic=0;   // No elastic and diffractive processes
 fMultiple=1;  // Include multiple interactions
 fEcmsmin=2.7; // Minimal CMS energy (in GeV) for events to get generated
 
 fEvent=0;
 
 fSpecpmin=0;
 
 fFrame="none";
 fWin=-1;
 fWxsec=0;
 
 fNucl=0;
 fZproj=0;
 fAproj=0;
 fZtarg=0;
 fAtarg=0;
 fFracpp=0;
 fFracnp=0;
 fFracpn=0;
 fFracnn=0;
 
 fOutFile=0;
 fOutTree=0;
 fMktree=0;
 fJob=0;
 fEvtuser.AddNamedSlot("BeamP");
 fEvtuser.AddNamedSlot("BeamTheta");
 fEvtuser.AddNamedSlot("BeamPhi");
 fEvtuser.AddNamedSlot("TargetP");
 fEvtuser.AddNamedSlot("TargetTheta");
 fEvtuser.AddNamedSlot("TargetPhi");
 
 fSelections=0;
 fSelect=0;
 
 fJetPpmin=0;
 fJetPpmax=0;
 fJetGpmin=0;
 fJetGpmax=0;
 fJetPspectrum=0;
 fJetPscale=0;
 fJetGspectrum=0;
 fJetGscale=0;
 
 TString s=GetName();
 s+=" (NcCollider)";
 SetName(s.Data());
 SetTitle("");
}

NcCollider::~NcCollider()
{
 
 if (fEvent)
 {
  delete fEvent;
  fEvent=0;
 }
 
 if (fSelections)
 {
  delete fSelections;
  fSelections=0;
 }
 
 // Deletion of the NcEvent output file also deletes the corresponding Tree
 if (fOutFile)
 {
  delete fOutFile;
  fOutFile=0;
  fOutTree=0;
 }
 
 if (fJob)
 {
  delete fJob;
  fJob=0;
 }
 
 if (fMktree)
 {
  delete fMktree;
  fMktree=0;
 }
 
 if (fJetPspectrum)
 {
  delete fJetPspectrum;
  fJetPspectrum=0;
 }
 
 if (fJetGspectrum)
 {
  delete fJetGspectrum;
  fJetGspectrum=0;
 }
}

NcTreeMaker* NcCollider::SetOutputFile(TString name,Int_t mode)
{
 
 // Expand the path name of the provided output filename
 name=gSystem->ExpandPathName(name.Data());
 
 // Flush and delete the current existing output file (if any) for the NcEvent data structures.
 // This will also delete the existing output tree connected to this file
 if (fOutFile)
 {
  if (fOutFile->IsOpen())
  {
   fOutFile->Write();
   delete fOutFile;
   fOutFile=0;
   fOutTree=0;
  }
 }
 
 // Close and delete the current existing output file (if any) for the plain ROOT tree data structures
 if (fMktree)
 {
  fMktree->CloseTree();
  delete fMktree;
  fMktree=0;
 }
 
 // Delete the NcJob environment
 if (fJob)
 {
  delete fJob;
  fJob=0;
 }
 
 TString fname;
 
 // Create the output file for the NcTreeMaker data structures
 if (mode>0)
 {
  fname=name;
  // Generate the corresponding file extension
  if (mode==2) fname+=".root";
 
  fMktree=new NcTreeMaker();
  fMktree->SetOutputFile(fname,"NcCollider event/track data in plain ROOT tree format");
 
  fJob=new NcJob("NcJob","NcCollider job (task) environment");
  fJob->Add(fMktree);
 
  printf(" *%-s::SetOutputFile* Plain ROOT tree event/track data will be written to output file: %-s \n",ClassName(),fname.Data());
 }
 
 // Create the output file for the NcEvent data structures
 if (mode==0 || mode==2)
 {
  fname=name;
  // Generate the corresponding file extension
  if (mode==2) fname+=".ncfspack";
 
  // Create a new NcEvent structure 
  if (fEvent)
  {
   delete fEvent;
   fEvent=0;
  }
  fEvent=new NcEvent();
  fEvent->SetOwner();
  fEvent->SetName(GetName());
  fEvent->SetTitle(GetTitle());
 
  // Create a new output file
  fOutFile=new TFile(fname,"RECREATE","NcCollider NcEvent data");
 
  // Create a new output Tree in the output file
  fOutTree=new TTree("T","NcCollider NcEvent data");
  Int_t bsize=32000;
  Int_t split=0;
  fOutTree->Branch("Events","NcEvent",&fEvent,bsize,split);
 
  printf(" *%-s::SetOutputFile* NcEvent data structures will be written to output file: %-s \n",ClassName(),fname.Data());
 }
 
 gROOT->cd(); // Make sure to work in the memory
 
 cout << endl;
 cout << endl;
 
 return fMktree;
}

void NcCollider::SetVertexMode(Int_t mode)
{
 
 if (mode<0 || mode >3)
 {
  cout << " *NcCollider::SetVertexMode* Invalid argument mode : " << mode << endl;
  fVertexmode=0;
 }
 else
 {
  fVertexmode=mode;
 }
}

Int_t NcCollider::GetVertexMode() const
{
 
 return fVertexmode;
}

void NcCollider::SetResolution(Double_t res)
{
 
 fResolution=fabs(res);
}

Double_t NcCollider::GetResolution() const
{
 
 return fResolution;
}

void NcCollider::SetRunNumber(Int_t run)
{
 
 fRunnum=run;
}

Int_t NcCollider::GetRunNumber() const
{
 
 return fRunnum;
}

void NcCollider::SetPrintFreq(Int_t n)
{
 
 fPrintfreq=n;
}

Int_t NcCollider::GetPrintFreq() const
{
 
 return fPrintfreq;
}

void NcCollider::SetUserControl(Int_t flag)
{
 
 fUserctrl=flag;
}

Int_t NcCollider::GetUserControl() const
{
 
 return fUserctrl;
}

void NcCollider::SetElastic(Int_t flag)
{
 
 fElastic=flag;
}

Int_t NcCollider::GetElastic() const
{
 
 return fElastic;
}

void NcCollider::SetMultiple(Int_t flag)
{
 
 fMultiple=flag;
}

Int_t NcCollider::GetMultiple() const
{
 
 return fMultiple;
}

void NcCollider::SetEcmsMin(Double_t ecms)
{
 
 fEcmsmin=ecms;
}

Double_t NcCollider::GetEcmsMin() const
{
 
 return fEcmsmin;
}

void NcCollider::SetRandomSeed(Int_t iseed)
{
 
 if (iseed>900000000) return;
 
 if (iseed<0)
 {
  NcTimestamp ts;
  Int_t mjd,sec,ns;
  ts.GetMJD(mjd,sec,ns);
  iseed=10000*sec+(mjd%10000);
 }
 SetMRPY(1,iseed);
 SetMRPY(2,0);
}

Int_t NcCollider::GetRandomSeed()
{
 
 Int_t iseed=GetMRPY(1);
 return iseed;
}

Float_t NcCollider::GetWin() const
{
 
 return fWin;
}

Int_t NcCollider::Init(TString frame,TString beam,TString target,Float_t win,Nc3Vector* pbeam,Nc3Vector* ptarget,Int_t wxsec,Double_t fact)
{
 
 Int_t ier=0;
 if (frame!="cms" && frame!="fixt" && frame!="free") ier=1;
 if (frame!="free" && win<=0) ier=1;
 if (frame=="free" && (!pbeam || !ptarget)) ier=1;
 if (frame=="free" && (beam.Contains("gamma/") || target.Contains("gamma/"))) ier=1;
 if (ier)
 {
  cout << " *NcCollider::Init* Standard Pythia initialisation ==> Inconsistent input." << endl;
  return ier;
 }
 
 if (!fUserctrl) // Optimisation of some MC parameters
 {
  if (beam=="gamma" || target=="gamma")
  {
   SetMSTP(14,10);      // Real photons for photon beams or targets
  }
  SetPARP(2,0);         // No minimum CMS energy required at initialisation
  SetMSTP(33,2);        // Activate K factor. Separate for ordinary and color annih. graphs
  SetPMAS(25,1,125.09); // Setting the Higgs mass
 }
 
 fWxsec=0;
 if (frame=="free")
 {
  SetMSTP(171,1);  // Enable variation of beam and/or target momenta
  SetMSTP(82,1);   // Select abrupt pt_min cut-off (see Pythia parameter PARP(81)) for multiple interactions 
  if (wxsec)
  {
   SetMSTP(172,2); // Weight event production by cross section
   fWxsec=1;
  }
  else
  {
   SetMSTP(172,1); // Just always generate an event at the requested energy
  }
 }
 else
 {
  SetMSTP(171,0); // Disable variation of beam and/or target momenta
 }
 
 if (fElastic) SetMSEL(2); // Include low-Pt, elastic and diffractive events
 
 if (!fMultiple) // Disable multiple interactions
 {
  SetMSTP(81,0);
  SetMSTP(82,1);
 }
 
 fEventnum=0;
 fNucl=0;
 fFrame=frame;
 fWin=-1;
 if (win>0) fWin=win;
 Double_t s;
 Double_t ecms;
 fBeam.SetNameTitle(beam.Data(),"Beam");
 fTarget.SetNameTitle(target.Data(),"Target");
 
 // Set initial beam and target specifications for user defined system via "pbeam" and "ptarget"
 if (fFrame=="free")
 {
  Double_t bmass=0;
  Double_t tmass=0;
 
  if (beam=="e-" || beam=="e+") bmass=0.510998928e-3;
  if (beam=="mu-" || beam=="mu+") bmass=105.6583715e-3;
  if (beam=="tau-" || beam=="tau+") bmass=1.77686;
  if (beam=="pi+" || beam=="pi-") bmass=139.57018e-3;
  if (beam=="pi0") bmass=134.9766e-3;
  if (beam=="K+" || beam=="K-") bmass=493.677e-3;
  if (beam=="KS0" || beam=="KL0") bmass=497.611e-3;
  if (beam=="p") bmass=938.272046e-3;
  if (beam=="n") bmass=939.565379e-3;
  if (beam=="Lambda0") bmass=1.115683;
  if (beam=="Sigma+") bmass=1.18937;
  if (beam=="Sigma0") bmass=1.192642;
  if (beam=="Sigma-") bmass=1.197449;
  if (beam=="Xi-") bmass=1.32171;
  if (beam=="Xi0") bmass=1.31486;
  if (beam=="Omega-") bmass=1.67245;
 
  if (target=="e-" || target=="e+") tmass=0.510998928e-3;
  if (target=="mu-" || target=="mu+") tmass=105.6583715e-3;
  if (target=="tau-" || target=="tau+") tmass=1.77686;
  if (target=="pi+" || target=="pi-") tmass=139.57018e-3;
  if (target=="pi0") tmass=134.9766e-3;
  if (target=="K+" || target=="K-") tmass=493.677e-3;
  if (target=="KS0" || target=="KL0") tmass=497.611e-3;
  if (target=="p") tmass=938.272046e-3;
  if (target=="n") tmass=939.565379e-3;
  if (target=="Lambda0") tmass=1.115683;
  if (target=="Sigma+") tmass=1.18937;
  if (target=="Sigma0") tmass=1.192642;
  if (target=="Sigma-") tmass=1.197449;
  if (target=="Xi-") tmass=1.32171;
  if (target=="Xi0") tmass=1.31486;
  if (target=="Omega-") tmass=1.67245;
 
  fBeam.SetMass(bmass);
  fTarget.SetMass(tmass);
 
  Nc3Vector vinit;
  if (fact>0 && fWin<0) // Modify beam 3-momentum for initialisation only
  {
   vinit.Load(*pbeam);
   vinit*=fact;
   SetMomentum(vinit,1);
   SetMomentum(*ptarget,2);
   cout << endl;
   cout << " ************************************* NcCollider::Init ************************************" << endl;
   cout << " *** Beam momentum artificially increased for intialisation only. Increase factor : " << fact << endl;
   cout << " *******************************************************************************************" << endl;
   cout << endl;
  }
  else if (fact<0 && fWin<0) // Modify target 3-momentum for initialisation only
  {
   vinit.Load(*ptarget);
   vinit*=fabs(fact);
   SetMomentum(*pbeam,1);
   SetMomentum(vinit,2);
   cout << endl;
   cout << " ************************************** NcCollider::Init *************************************" << endl;
   cout << " *** Target momentum artificially increased for intialisation only. Increase factor : " << fabs(fact) << endl;
   cout << " *********************************************************************************************" << endl;
   cout << endl;
  }
  else // Use the provided beam and target 3-momenta for initialisation
  {
   SetMomentum(*pbeam,1);
   SetMomentum(*ptarget,2);
  }
  frame="3mom"; // Use the Pythia convention for the frame name
  if (win>=0)   // Forced event generation in the CMS
  {
   frame="cms";
   if (!win) // Set CM energy according to "pbeam", "ptarget" and "fact"
   {
    s=fBeam.GetInvariant()+fTarget.GetInvariant()+2.*fBeam.Dot(fTarget);
    ecms=sqrt(s);
    fWin=ecms;
   }
  }
 }
 
 // Prevent title overwriting by Initialize()
 TString title=GetTitle();
 
 Initialize(frame.Data(),beam.Data(),target.Data(),fWin);
 
 SetTitle(title.Data());
 
 // Use the Pythia beam and target specifications for consistency
 fBeam.SetMass(GetP(1,5));
 fTarget.SetMass(GetP(2,5));
 
 if (fFrame=="cms")
 {
  fBeam.Set3Vector(GetP(1,1),GetP(1,2),GetP(1,3),"car");
  fTarget.Set3Vector(GetP(2,1),GetP(2,2),GetP(2,3),"car");
 }
 if (fFrame=="fixt")
 {
  fBeam.Set3Vector(0,0,win,"car");
  fTarget.Set3Vector(0,0,0,"car");
 }
 if (fFrame=="free")
 {
  SetMomentum(*pbeam,1);
  SetMomentum(*ptarget,2);
 }
 
 TString sweight="Yes";
 if (!wxsec) sweight="No";
 
 s=fBeam.GetInvariant()+fTarget.GetInvariant()+2.*fBeam.Dot(fTarget);
 ecms=sqrt(s);
 
 cout << endl;
 cout << endl;
 cout << " ********************************************************" << endl;
 cout << " *** NcCollider::Init  Standard Pythia initialisation ***" << endl;
 cout << " ********************************************************" << endl;
 cout << " *** Beam particle: " << beam.Data() << " Target particle: " << target.Data() << " Frame: " << fFrame.Data() << endl;
 if (fFrame=="cms") cout << " *** Total CMS energy: " << win << " GeV" << endl;
 if (fFrame=="fixt") cout << " *** Beam particle momentum: " << win << " GeV/c" << endl;
 if (fFrame=="free") cout << " *** Event weighting by cross section: " << sweight.Data() << endl;
 cout << " *** Beam   particle 3-momentum (GeV/c): px=" << fBeam.GetX(1,"car") << " py=" << fBeam.GetX(2,"car") << " pz=" << fBeam.GetX(3,"car") << endl; 
 cout << " *** Target particle 3-momentum (GeV/c): px=" << fTarget.GetX(1,"car") << " py=" << fTarget.GetX(2,"car") << " pz=" << fTarget.GetX(3,"car") << endl; 
 if (fFrame!="cms") cout << " *** Total CMS energy: " << ecms << " GeV" << endl;
 if (fFrame=="free" && fWin>0) cout << " *** Forced CMS processing. Cross sections initialised for a CMS energy of " << fWin << " GeV" << endl;
 if (fOutFile) cout << " *** NcEvent data structures will be written to output file: " << fOutFile->GetName() << endl;
 cout << endl;
 cout << endl;
 
 return ier;
}

Int_t NcCollider::Init(TString frame,Int_t zp,Int_t ap,Int_t zt,Int_t at,Float_t win,Nc3Vector* pbeam,Nc3Vector* ptarget,Int_t wxsec)
{
 
 Int_t ier=0;
 if (frame!="cms" && frame!="fixt" && frame!="free") ier=1;
 if (frame!="free" && win<=0) ier=1;
 if (frame=="free" && (!pbeam || !ptarget)) ier=1;
 if (ier)
 {
  cout << " *NcCollider::Init* Nucleus-Nucleus generator initialisation ==> Inconsistent input." << endl;
  return ier;
 }
 
 if (!fUserctrl) // Optimisation of some MC parameters
 {
  SetPARP(2,0);         // No minimum CMS energy required at initialisation
  SetMSTP(33,2);        // Activate K factor. Separate for ordinary and color annih. graphs
  SetPMAS(25,1,125.09); // Setting the Higgs mass
 }
 
 fWxsec=0;
 if (frame=="free")
 {
  SetMSTP(171,1);  // Enable variation of beam and/or target momenta
  SetMSTP(82,1);   // Select abrupt pt_min cut-off (see Pythia parameter PARP(81)) for multiple interactions 
  if (wxsec)
  {
   SetMSTP(172,2); // Weight event production by cross section
   fWxsec=1;
  }
  else
  {
   SetMSTP(172,1); // Just always generate an event at the requested energy
  }
 }
 else
 {
  SetMSTP(171,0); // Disable variation of beam and/or target momenta
 }
 
 if (fElastic) SetMSEL(2); // Include low-Pt, elastic and diffractive events
 
 if (!fMultiple) // Disable multiple interactions
 {
  SetMSTP(81,0);
  SetMSTP(82,1);
 }
 
 fEventnum=0;
 fNucl=1;
 fFrame=frame;
 fWin=-1;
 if (win>0) fWin=win;
 fZproj=0;
 fAproj=0;
 fZtarg=0;
 fAtarg=0;
 fFracpp=0;
 fFracnp=0;
 fFracpn=0;
 fFracnn=0;
 
 if (ap<1 || at<1 || zp>ap || zt>at)
 {
  cout << endl;
  cout << " *NcCollider::Init* Invalid input value(s). Zproj = " << zp
       << " Aproj = " << ap << " Ztarg = " << zt << " Atarg = " << at << endl;
  return 1;
 }
 
 fZproj=zp;
 fAproj=ap;
 fZtarg=zt;
 fAtarg=at;
 
 TString beam="(Z=";
 beam+=zp;
 beam+=",A=";
 beam+=ap;
 beam+=")";
 fBeam.SetNameTitle(beam.Data(),"Beam");
 
 TString target="(Z=";
 target+=zt;
 target+=",A=";
 target+=at;
 target+=")";
 fTarget.SetNameTitle(target.Data(),"Target");
 
 Float_t mnuc=(0.9382720+0.9395654)/2.;
 fBeam.SetMass(mnuc);
 fTarget.SetMass(mnuc);
 
 if (fFrame=="cms")
 {
  Float_t pcms=sqrt((win*win/4.)-mnuc*mnuc);
  fBeam.Set3Vector(0,0,pcms,"car");
  fTarget.Set3Vector(0,0,-pcms,"car");
 }
 if (fFrame=="fixt")
 {
  fBeam.Set3Vector(0,0,win,"car");
  fTarget.Set3Vector(0,0,0,"car");
 }
 if (fFrame=="free")
 {
  SetMomentum(*pbeam,1);
  SetMomentum(*ptarget,2);
 }
 
 TString sweight="Yes";
 if (!wxsec) sweight="No";
 
 Double_t s=fBeam.GetInvariant()+fTarget.GetInvariant()+2.*fBeam.Dot(fTarget);
 Double_t ecms=sqrt(s);
 
 // Set CM energy according to "pbeam" and "ptarget" for forced CMS event generation
 if (fFrame=="free" && !win) fWin=ecms;
 
 cout << endl;
 cout << endl;
 cout << " ******************************************************************" << endl;
 cout << " *** NcCollider::Init  Nucleus-Nucleus generator initialisation ***" << endl;
 cout << " ******************************************************************" << endl;
 cout << " *** Beam nucleus: " << fBeam.GetName() << " Target nucleus: " << fTarget.GetName() << " Frame: " << fFrame.Data() << endl;
 if (fFrame=="cms") cout << " *** Total CMS energy per nucleon-nucleon collision: " << win << " GeV" << endl;
 if (fFrame=="fixt") cout << " *** Beam momentum: " << win << " GeV/c per nucleon" << endl;
 if (fFrame=="free") cout << " *** Event weighting by cross section: " << sweight.Data() << endl;
 cout << " *** Beam   3-momentum in GeV/c per nucleon: px=" << fBeam.GetX(1,"car") << " py=" << fBeam.GetX(2,"car") << " pz=" << fBeam.GetX(3,"car") << endl; 
 cout << " *** Target 3-momentum in GeV/c per nucleon: px=" << fTarget.GetX(1,"car") << " py=" << fTarget.GetX(2,"car") << " pz=" << fTarget.GetX(3,"car") << endl; 
 if (fFrame!="cms") cout << " *** Total CMS energy per nucleon-nucleon collision: " << ecms << " GeV" << endl;
 if (fFrame=="free" && fWin>0) cout << " *** Forced CMS processing. Cross sections initialised for a nucleon-nucleon CMS energy of " << fWin << " GeV" << endl;
 if (fOutFile) cout << " *** Event data will be written to output file: " << fOutFile->GetName() << endl;
 cout << endl;
 cout << endl;
 
 return ier;
}

void NcCollider::GetFractions(Float_t zp,Float_t ap,Float_t zt,Float_t at)
{
 
 if (zp<0) zp=0;
 if (zt<0) zt=0;
 
 fFracpp=0;
 fFracnp=0;
 fFracpn=0;
 fFracnn=0;
 
 if (ap>0 && at>0)
 {
  fFracpp=(zp/ap)*(zt/at);
  fFracnp=(1.-zp/ap)*(zt/at);
  fFracpn=(zp/ap)*(1.-zt/at);
  fFracnn=(1.-zp/ap)*(1.-zt/at);
 }
}

void NcCollider::SetMomentum(Nc3Vector& p,Int_t mode)
{
 
 if (mode!=1 && mode!=2)
 {
  cout << " *NcCollider::SetMomentum* Invalid input value. mode = " << mode << endl;
  return;  
 }
 
 if (fFrame!="free")
 {
  cout << " *NcCollider::SetMomentum* Not valid for frame = " << fFrame.Data() << endl;
  return;  
 }
 
 if (mode==1) fBeam.Set3Momentum(p);
 if (mode==2) fTarget.Set3Momentum(p);
 
 
 if (fWin<0) // Update the internal Pythia beam c.q. target momentum
 {
  SetP(mode,1,p.GetX(1,"car"));
  SetP(mode,2,p.GetX(2,"car"));
  SetP(mode,3,p.GetX(3,"car"));
 }
 else // Update the corresponding Ecms scaling factor in case of forced CMS processing
 {
  Double_t s=fBeam.GetInvariant()+fTarget.GetInvariant()+2.*fBeam.Dot(fTarget);
  Double_t ecms=sqrt(s);
  SetPARP(171,ecms/fWin);
 }
}

Int_t NcCollider::MakeEvent(Int_t npt,Int_t mlist,Int_t medit)
{
 
 fEventnum++; 
 
 Int_t specmode=1;
 if (npt<0)
 {
  specmode=0;
  npt=abs(npt);
 }
 
 // Counters for the various (proj,targ) combinations : p+p, n+p, p+n and n+n
 Int_t ncols[4]={0,0,0,0};
 
 Int_t zp=0;
 Int_t ap=0;
 Int_t zt=0;
 Int_t at=0;
 
 Int_t ncol=1;
 
 if (fNucl)
 {
  if (npt<2 || npt>(fAproj+fAtarg))
  {
   cout << " *NcCollider::MakeEvent* Invalid input value. npt = " << npt
        << " Aproj = " << fAproj << " Atarg = " << fAtarg << endl;
   return -1;
  }
 
  // Determine the number of nucleon-nucleon collisions
  ncol=npt/2;
  if (npt%2 && fRan.Uniform()>0.5) ncol+=1;
 
  // Determine the number of the various types of N+N interactions
  zp=fZproj;
  ap=fAproj;
  zt=fZtarg;
  at=fAtarg;
  Int_t maxa=2; // Indicator whether proj (1) or target (2) has maximal A left
  if (ap>at) maxa=1;
  Float_t* rans=new Float_t[ncol];
  fRan.Uniform(rans,ncol);
  Float_t rndm=0;
  for (Int_t i=0; i<ncol; i++)
  {
   GetFractions(zp,ap,zt,at);
   rndm=rans[i];
   if (rndm<=fFracpp) // p+p interaction
   {
    ncols[0]++;
    if (maxa==2)
    {
     at--;
     zt--;
    } 
    else
    {
     ap--;
     zp--;
    }
   }
   if (rndm>fFracpp && rndm<=(fFracpp+fFracnp)) // n+p interaction
   {
    ncols[1]++;
    if (maxa==2)
    {
     at--;
     zt--;
    } 
    else
    {
     ap--;
    }
   }
   if (rndm>(fFracpp+fFracnp) && rndm<=(fFracpp+fFracnp+fFracpn)) // p+n interaction
   {
    ncols[2]++;
    if (maxa==2)
    {
     at--;
    } 
    else
    {
     ap--;
     zp--;
    }
   }
   if (rndm>(fFracpp+fFracnp+fFracpn)) // n+n interaction
   {
    ncols[3]++; 
    if (maxa==2)
    {
     at--;
    } 
    else
    {
     ap--;
    }
   }
  }
  delete [] rans;
 }
 
 TString sweight="Yes";
 if (!fWxsec) sweight="No";
 
 Double_t s=fBeam.GetInvariant()+fTarget.GetInvariant()+2.*fBeam.Dot(fTarget);
 Double_t ecms=sqrt(s);
 
 if (fPrintfreq)
 {
  if (!(fEventnum%fPrintfreq))
  {
   cout << " *NcCollider::MakeEvent* Run : " << fRunnum << " Event : " << fEventnum << endl;
   if (fFrame=="free") cout << "  Event weighting by cross section: " << sweight.Data() << endl;
 
   if (fNucl)
   {
    cout << "  Beam nucleus: " << fBeam.GetName() << " Target nucleus: " << fTarget.GetName() << " Frame: " << fFrame.Data() << endl;
    cout << "  Beam   3-momentum in GeV/c per nucleon: px=" << fBeam.GetX(1,"car") << " py=" << fBeam.GetX(2,"car") << " pz=" << fBeam.GetX(3,"car") << endl; 
    cout << "  Target 3-momentum in GeV/c per nucleon: px=" << fTarget.GetX(1,"car") << " py=" << fTarget.GetX(2,"car") << " pz=" << fTarget.GetX(3,"car") << endl; 
    cout << "  Total CMS energy per nucleon-nucleon collision: " << ecms << " GeV" << endl;
    cout << "  # participants and collisions: npart=" << npt << " ncol=" << ncol 
         << " ncolpp=" << ncols[0] << " ncolnp=" << ncols[1] << " ncolpn=" << ncols[2] << " ncolnn=" << ncols[3] << endl;
   }
   else
   {
    cout << "  Beam particle: " << fBeam.GetName() << " Target particle: " << fTarget.GetName() << " Frame: " << fFrame.Data() << endl;
    cout << "  Beam   particle 3-momentum (GeV/c): px=" << fBeam.GetX(1,"car") << " py=" << fBeam.GetX(2,"car") << " pz=" << fBeam.GetX(3,"car") << endl; 
    cout << "  Target particle 3-momentum (GeV/c): px=" << fTarget.GetX(1,"car") << " py=" << fTarget.GetX(2,"car") << " pz=" << fTarget.GetX(3,"car") << endl; 
    cout << "  Total CMS energy: " << ecms << " GeV" << endl;
   }
 
   if (ecms<fEcmsmin) printf("  *** No event generated. Ecms is below the minimal requirement of : %-g GeV. \n \n",fEcmsmin);
  }
 }
 
 if (ecms<fEcmsmin) return 0;
 
 if (!fEvent)
 {
  fEvent=new NcEvent();
  fEvent->SetOwner();
  fEvent->SetName(GetName());
  fEvent->SetTitle(GetTitle());
 }
 else
 {
  fEvent->Reset();
 }
 fEvent->SetRunNumber(fRunnum);
 fEvent->SetEventNumber(fEventnum);
 
 // Set event title text if not provided by the user
 TString title=GetTitle();
 if (title=="")
 {
  title=fBeam.GetName();
  title+=" on ";
  title+=fTarget.GetName();
  title+=" collision";
  fEvent->SetTitle(title.Data());
 }
 
 NcTrack t;
 Nc3Vector p;
 NcPosition r,rx;
 Float_t v[3];
 NcVertex vert;
 Nc3Vector pproj,ptarg;
 Nc3Vector v3;
 Nc4Vector v4;
 Double_t ct=0;
 
 // The Lorentz boost for the case of forced CMS processing
 Nc4Vector ptot;
 ptot.SetInvariant(s);
 p=fBeam.Get3Momentum();
 v3=fTarget.Get3Momentum();
 p+=v3;
 ptot.Set3Vector(p);
 fLorbo.Set4Momentum(ptot);
 
 if (fVertexmode)
 {
  // Make sure the primary vertex gets correct location and Id=1
  v[0]=0;
  v[1]=0;
  v[2]=0;
  r.SetPosition(v,"car");
  v[0]=fResolution;
  v[1]=fResolution;
  v[2]=fResolution;
  r.SetPositionErrors(v,"car");
 
  vert.SetId(1);
  vert.SetTrackCopy(0);
  vert.SetVertexCopy(0);
  vert.SetPosition(r);
  fEvent->AddVertex(vert,0);
 }
 
 Int_t kf=0;
 Float_t charge=0,mass=0;
 TString name;
 
 Int_t ntypes=4;
 
 // Singular settings for a normal Pythia elementary particle interation 
 if (!fNucl)
 {
  ntypes=1;
  ncols[0]=1;
 }
 
 // Generate all the various collisions
 fSelect=0;      // Flag to indicate whether the total event is selected or not 
 Int_t select=0; // Flag to indicate whether the sub-event is selected or not 
 Int_t first=1;  // Flag to indicate the first collision process
 TString frame=fFrame;
 if (fFrame=="free")
 {
  frame="3mom";            // Use the Pythia convention for the frame name
  if (fWin>0) frame="cms"; // Forced event generation in the CMS
 }
 pproj=fBeam.Get3Momentum();
 ptarg=fTarget.Get3Momentum();
 Int_t npart=0,ntk=0;
 Double_t dist=0;
 for (Int_t itype=0; itype<ntypes; itype++)
 {
  if (fNucl)
  {
   if (fFrame=="free")
   {
    SetMomentum(pproj,1);
    SetMomentum(ptarg,2);
   }
 
   if (itype==0 && ncols[itype]) Initialize(frame.Data(),"p","p",fWin);
   if (itype==1 && ncols[itype]) Initialize(frame.Data(),"n","p",fWin);
   if (itype==2 && ncols[itype]) Initialize(frame.Data(),"p","n",fWin);
   if (itype==3 && ncols[itype]) Initialize(frame.Data(),"n","n",fWin);
 
   if (fFrame=="free")
   {
    SetMomentum(pproj,1);
    SetMomentum(ptarg,2);
   }
  }
 
  for (Int_t jcol=0; jcol<ncols[itype]; jcol++)
  {
 
   GenerateEvent();
 
   select=IsSelected();
   if (select) fSelect=1;
 
   if (first) // Store generator parameter information in the event structure
   {
    // Enter generator parameters as a device in the event
 
    NcDevice params;
    params.SetNameTitle("NcCollider","NcCollider generator parameters");
    params.SetSlotName("Medit",1);
    params.SetSlotName("Vertexmode",2);
    params.SetSlotName("Resolution",3);
    params.SetSlotName("Userctrl",4);
    params.SetSlotName("Elastic",5);
    params.SetSlotName("Multiple",6);
    params.SetSlotName("Wxsec",7);
    params.SetSlotName("Ecms",8);
 
    params.SetSignal(medit,1);
    params.SetSignal(fVertexmode,2);
    params.SetSignal(fResolution,3);
    params.SetSignal(fUserctrl,4);
    params.SetSignal(fElastic,5);
    params.SetSignal(fMultiple,6);
    params.SetSignal(fWxsec,7);
    params.SetSignal(ecms,8);
 
    // Store projectile and target information in the event structure
    if (fNucl)
    {
     fEvent->SetProjectile(fAproj,fZproj,pproj);
     fEvent->SetTarget(fAtarg,fZtarg,ptarg);
 
     params.AddNamedSlot("specmode");
     params.AddNamedSlot("Specpmin");
     params.AddNamedSlot("npart");
     params.AddNamedSlot("ncolpp");
     params.AddNamedSlot("ncolnp");
     params.AddNamedSlot("ncolpn");
     params.AddNamedSlot("ncolnn");
 
     params.SetSignal(specmode,"specmode");
     params.SetSignal(fSpecpmin,"Specpmin");
     params.SetSignal(npt,"npart");
     params.SetSignal(ncols[0],"ncolpp");
     params.SetSignal(ncols[1],"ncolnp");
     params.SetSignal(ncols[2],"ncolpn");
     params.SetSignal(ncols[3],"ncolnn");
    }
    else
    {
     kf=GetK(1,2);
     fEvent->SetProjectile(0,0,pproj,kf);
     kf=GetK(2,2);
     fEvent->SetTarget(0,0,ptarg,kf);
    }
 
    fEvent->AddDevice(params);
 
    first=0;
   }
 
   if (medit >= 0) Pyedit(medit); // Define which particles are to be kept
 
   if (mlist>=0 && select)
   {
    Pylist(mlist);
    cout << endl;
   }
 
   npart=GetN();
   for (Int_t jpart=1; jpart<=npart; jpart++)
   {
    kf=GetK(jpart,2);
    charge=Pychge(kf)/3.;
    mass=GetP(jpart,5);
    name=GetPyname(kf);
 
    // 3-momentum in GeV/c
    v[0]=GetP(jpart,1);
    v[1]=GetP(jpart,2);
    v[2]=GetP(jpart,3);
    p.SetVector(v,"car");
 
    // Production location in meter.
    v[0]=GetV(jpart,1)/1000.;
    v[1]=GetV(jpart,2)/1000.;
    v[2]=GetV(jpart,3)/1000.;
    r.SetPosition(v,"car");
 
    ct=GetV(jpart,4)/1000.;
 
    // Boost the track momentum and vertex location into the user defined frame in case of forced CMS processing
    if (fFrame=="free" && fWin>0)
    {
     // Momentum boost
     v4.SetInvariant(mass*mass);
     v4.Set3Vector(p);
     v4=fLorbo.Inverse(v4);
     p=v4.Get3Vector();
 
     // Vertex location boost
     v4.SetInvariant(pow(ct,2)-r.Dot(r));
     v4.Set3Vector(r);
     v4=fLorbo.Inverse(v4);
     v3=v4.Get3Vector();
     r.SetPosition(v3);
    }    
 
    ntk++;
 
    t.Reset();
    t.SetId(ntk);
    t.SetParticleCode(kf);
    t.SetName(name.Data());
    t.SetCharge(charge);
    t.SetMass(mass);
    t.Set3Momentum(p);
    t.SetBeginPoint(r);
 
    fEvent->AddTrack(t);
 
    // Build the vertex structures if requested
    if (fVertexmode)
    {
     // Check if track belongs within the resolution to an existing vertex
     Int_t add=0;  
     for (Int_t jv=1; jv<=fEvent->GetNvertices(); jv++)
     {
      NcVertex* vx=fEvent->GetVertex(jv);
      if (vx)
      {
       rx=vx->GetPosition();
       dist=rx.GetDistance(r);
       if (dist < fResolution)
       {
        NcTrack* tx=fEvent->GetIdTrack(ntk);
        if (tx)
        { 
         vx->AddTrack(tx);
         add=1;
        }
       }
      }
      if (add) break; // No need to look further for vertex candidates
     }
 
     // If track was not close enough to an existing vertex
     // a new secondary vertex is created      
     if (!add && fVertexmode>1)
     {
      NcTrack* tx=fEvent->GetIdTrack(ntk);
      if (tx)
      {
       v[0]=fResolution;
       v[1]=fResolution;
       v[2]=fResolution;
       r.SetPositionErrors(v,"car");
       vert.Reset();
       vert.SetTrackCopy(0);
       vert.SetVertexCopy(0);
       vert.SetId((fEvent->GetNvertices())+1);
       vert.SetPosition(r);
       vert.AddTrack(tx);
       fEvent->AddVertex(vert,0);
      } 
     }
    }
   } // End of loop over the produced particles for each collision
  } // End of loop over number of collisions for each type
 } // End of loop over collision types
 
 // Link sec. vertices to primary if requested
 // Note that also the connecting tracks are automatically created
 if (fVertexmode>2)
 {
  NcVertex* vp=fEvent->GetIdVertex(1); // Primary vertex
  if (vp)
  {
   for (Int_t i=2; i<=fEvent->GetNvertices(); i++)
   {
    NcVertex* vx=fEvent->GetVertex(i);
    if (vx)
    {
     if (vx->GetId() != 1) vp->AddVertex(vx);
    }
   }
  }
 }
 
 // Include the spectator tracks in the event structure.
 if (fNucl && specmode)
 {
  v[0]=0;
  v[1]=0;
  v[2]=0;
  r.SetPosition(v,"car");
 
  zp=fZproj-(ncols[0]+ncols[2]);
  if (zp<0) zp=0;
  ap=fAproj-(ncols[0]+ncols[1]+ncols[2]+ncols[3]);
  if (ap<0) ap=0;
  zt=fZtarg-(ncols[0]+ncols[1]);
  if (zt<0) zt=0;
  at=fAtarg-(ncols[0]+ncols[1]+ncols[2]+ncols[3]);
  if (at<0) at=0;
 
  Int_t nspec=0;
 
  if (pproj.GetNorm() > fSpecpmin)
  {
   kf=2212; // Projectile spectator protons
   charge=Pychge(kf)/3.;
   mass=GetPMAS(Pycomp(kf),1);
   name=GetPyname(kf);
   for (Int_t iprojp=1; iprojp<=zp; iprojp++)
   {
    nspec++;
    t.Reset();
    t.SetId(-nspec);
    t.SetParticleCode(kf);
    t.SetName(name.Data());
    t.SetTitle("Projectile spectator proton");
    t.SetCharge(charge);
    t.SetMass(mass);
    t.Set3Momentum(pproj);
    t.SetBeginPoint(r);
 
    fEvent->AddTrack(t);
   }
 
   kf=2112; // Projectile spectator neutrons
   charge=Pychge(kf)/3.;
   mass=GetPMAS(Pycomp(kf),1);
   name=GetPyname(kf);
   for (Int_t iprojn=1; iprojn<=(ap-zp); iprojn++)
   {
    nspec++;
    t.Reset();
    t.SetId(-nspec);
    t.SetParticleCode(kf);
    t.SetName(name.Data());
    t.SetTitle("Projectile spectator neutron");
    t.SetCharge(charge);
    t.SetMass(mass);
    t.Set3Momentum(pproj);
    t.SetBeginPoint(r);
 
    fEvent->AddTrack(t);
   }
  }
 
  if (ptarg.GetNorm() > fSpecpmin)
  {
   kf=2212; // Target spectator protons
   charge=Pychge(kf)/3.;
   mass=GetPMAS(Pycomp(kf),1);
   name=GetPyname(kf);
   for (Int_t itargp=1; itargp<=zt; itargp++)
   {
    nspec++;
    t.Reset();
    t.SetId(-nspec);
    t.SetParticleCode(kf);
    t.SetName(name.Data());
    t.SetTitle("Target spectator proton");
    t.SetCharge(charge);
    t.SetMass(mass);
    t.Set3Momentum(ptarg);
    t.SetBeginPoint(r);
 
    fEvent->AddTrack(t);
   }
 
   kf=2112; // Target spectator neutrons
   charge=Pychge(kf)/3.;
   mass=GetPMAS(Pycomp(kf),1);
   name=GetPyname(kf);
   for (Int_t itargn=1; itargn<=(at-zt); itargn++)
   {
    nspec++;
    t.Reset();
    t.SetId(-nspec);
    t.SetParticleCode(kf);
    t.SetName(name.Data());
    t.SetTitle("Target spectator neutron");
    t.SetCharge(charge);
    t.SetMass(mass);
    t.Set3Momentum(ptarg);
    t.SetBeginPoint(r);
 
    fEvent->AddTrack(t);
   }
  }
 
  // Link the spectator tracks to the primary vertex.
  if (fVertexmode)
  {
   NcVertex* vp=fEvent->GetIdVertex(1);
   if (vp)
   {
    for (Int_t ispec=1; ispec<=nspec; ispec++)
    {
     NcTrack* tx=fEvent->GetIdTrack(-ispec);
     if (tx) vp->AddTrack(tx);
    }
   }
  }
 }
 
 if (fPrintfreq)
 {
  if (!(fEventnum%fPrintfreq) && (mlist || fEvent))
  {
   if (fEvent) printf("  Number of tracks in the event structure : %-i \n",fEvent->GetNtracks());
   printf("\n"); // Create empty output line after the event
  }
 }
 
 // Record the actual beam and target momenta as user data in the event structure
 if (fSelect)
 {
  fEvtuser.SetSignal(fBeam.GetX(1,"sph"),"BeamP");
  fEvtuser.SetSignal(fBeam.GetX(2,"sph","rad"),"BeamTheta");
  fEvtuser.SetSignal(fBeam.GetX(3,"sph","rad"),"BeamPhi");
  fEvtuser.SetSignal(fTarget.GetX(1,"sph"),"TargetP");
  fEvtuser.SetSignal(fTarget.GetX(2,"sph","rad"),"TargetTheta");
  fEvtuser.SetSignal(fTarget.GetX(3,"sph","rad"),"TargetPhi");
  fEvent->SetUserData(fEvtuser);
 }
 
 if (fOutTree && fSelect) fOutTree->Fill();
 
 if (fJob && fSelect) fJob->ProcessObject(fEvent);
 
 return fSelect;
}

NcEvent* NcCollider::GetEvent(Int_t select) const
{
 
 if (!select || fSelect)
 {
  return fEvent;
 }
 else
 {
  return 0;
 }
}

void NcCollider::EndRun()
{
 
 if (!fOutFile && !fMktree) return;
 
 // Deleting the NcEvent output file also deletes the corresponding Tree
 if (fOutFile)
 {
  if (fOutFile->IsOpen())
  {
   fOutFile->Write();
   fOutFile->Close();
   delete fOutFile;
   fOutFile=0;
   fOutTree=0;
  }
 }
 
 if (fMktree)
 {
  fMktree->CloseTree();
  delete fMktree;
  fMktree=0;
 }
 
 printf(" *%-s::EndRun* Output file(s) correctly written and closed. \n",ClassName());
}

void NcCollider::SetStable(Int_t id,Int_t mode,Int_t cls)
{
 
 if (mode==0 || mode==1)
 {
  Int_t kc=0;
  Int_t decay=1-mode;
  
  // Introduce additional c.q. missing decay channels
  Int_t idc=7000;   // Start entry point in the decay table for the new decay processes
  SetPARJ(64,5e-4); // Reduce the allowed minimal mass difference in decays
 
  // Specification of neutron decay data
  Double_t ctau=2.6391e14; // Average ctau in mm
  kc=Pycomp(2112);         // kc code for (anti)neutrons
  SetPMAS(kc,4,ctau);      // Set average lifetime via ctau in mm
  SetMDCY(kc,2,idc);       // Set the idc entry point for this particle decay
  SetMDCY(kc,3,1);         // Number of decay modes
  SetMDME(idc,1,1);        // Activate this decay channel
  SetMDME(idc,2,0);        // Use a normal decay matrix element for this decay channel
  SetBRAT(idc,1);          // Set branching ratio for this decay channel
  SetKFDP(idc,1,2212);     // The "p" decay product
  SetKFDP(idc,2,11);       // The "e-" decay product
  SetKFDP(idc,3,-12);      // The "nue_bar" decay product
 
  if (!id) // Settings for a class of particles
  {
   if (cls>0 && cls<6)
   {
    for (Int_t i=1; i<1000; i++)
    {
     kc=Pycomp(i);
     if (kc>0)
     {
      if (cls==1 && i<10) SetMDCY(kc,1,decay);
      if (cls==2 && i>10 && i<20) SetMDCY(kc,1,decay);
      if (cls==3 && i>20 && i<26) SetMDCY(kc,1,decay);
      if (cls==4 && i>100 && i<1000) SetMDCY(kc,1,decay);
      if (cls==5 && i>1000) SetMDCY(kc,1,decay);
     }
    }
    // Special entries for Psi' and Ypsilon'
    if (cls==4)
    {
     kc=Pycomp(100443); 
     if (kc>0) SetMDCY(kc,1,decay);
     kc=Pycomp(100553); 
     if (kc>0) SetMDCY(kc,1,decay);
    }
   }
   else
   {
    cout << " *NcCollider::SetStable* Invalid parameter. cls = " << cls << endl;
   }
  }
  else // Settings for individual particles
  {
   kc=Pycomp(id);
   if (kc>0)
   {
    SetMDCY(kc,1,decay);
   }
   else 
   {
    cout << " *NcCollider::SetStable* Unknown particle code. id = " << id << endl;
   }
  }
 } 
 else
 {
  cout << " *NcCollider::SetStable* Invalid parameter. mode = " << mode << endl;
 }
}

void NcCollider::SelectEvent(Int_t id)
{
 
 if (!id)
 {
  if (fSelections)
  {
   delete fSelections;
   fSelections=0;
  }
 }
 else
 {
  Int_t kc=Pycomp(id);
  if (!fSelections)
  {
   fSelections=new TArrayI(1);
   fSelections->AddAt(kc,0);
  }
  else
  {
   Int_t exist=0;
   Int_t size=fSelections->GetSize();
   for (Int_t i=0; i<size; i++)
   {
    if (kc==fSelections->At(i))
    {
     exist=1;
     break;
    }
   }
  
   if (!exist)
   {
    fSelections->Set(size+1);
    fSelections->AddAt(kc,size);
   }
  }
 }
}

Int_t NcCollider::GetSelectionFlag() const
{
 
 return fSelect;
}

Int_t NcCollider::IsSelected()
{
 
 if (GetMSTI(61)) return 0;
 
 if (!fSelections) return 1; 
 
 Int_t nsel=fSelections->GetSize();
 Int_t npart=GetN();
 Int_t kf,kc;
 
 Int_t select=0;
 for (Int_t jpart=1; jpart<=npart; jpart++)
 {
  kf=GetK(jpart,2);
  kc=Pycomp(kf);
  for (Int_t i=0; i<nsel; i++)
  {
   if (kc==fSelections->At(i))
   {
    select=1;
    break;
   }
  }
  if (select) break;
 }
 return select;
}

void NcCollider::SetSpectatorPmin(Float_t pmin)
{
 
 fSpecpmin=pmin;
}

Float_t NcCollider::GetSpectatorPmin() const
{
 
 return fSpecpmin;
}

TString NcCollider::GetPyname(Int_t kf)
{
 
 char name[16];
 TString sname;
 Pyname(kf,name);
 sname=name[0];
 for (Int_t i=1; i<16; i++)
 {
  if (name[i]==' ') break;
  sname=sname+name[i];
 }
 return sname;
}

void NcCollider::SetJetProtonSpectrum(Double_t pmin,Double_t pmax,TF1* fspec,TH1* hspec,Int_t mode)
{
 
 gROOT->cd(); // Make sure to work in memory
 
 fJetPpmin=0;
 fJetPpmax=0;
 
 if (fJetPspectrum)
 {
  delete fJetPspectrum;
  fJetPspectrum=0;
 }
 
 if (pmax<=pmin && pmin<=0)
 {
  printf(" *%-s::SetJetProtonSpectrum* Inconsistent input pmin=%-g pmax=%-g \n",ClassName(),pmin,pmax);
  return;
 }
 
 if (pmax>pmin && !fspec && !hspec)
 {
  printf(" *%-s::SetJetProtonSpectrum* Inconsistent input pmin=%-g pmax=%-g fspec=0 hspec=0 \n",ClassName(),pmin,pmax);
  return;
 }
 
 fJetPpmin=pmin;
 fJetPpmax=pmax;
 if (pmax<=pmin) fJetPpmax=pmin;
 fJetPscale=0;
 
 // Momentum distribution specified by a histogram
 if (hspec)
 {
  fJetPscale=mode;
  if (mode==1)
  {
   pmin=log10(pmin);
   pmax=log10(pmax);
  }
  if (mode==2)
  {
   pmin=log(pmin);
   pmax=log(pmax);
  }
 
  if (!fspec) // Histogram contains N vs. p distribution
  {
   fJetPspectrum=(TH1*)hspec->Clone();
  }
  else // Histogram contains dN/dp distribution scaled by "fspec"
  {
   TH1F hpbeam=fLab.GetCountsHistogram(*hspec,mode,"",fspec);
   fJetPspectrum=(TH1*)hpbeam.Clone();
   fJetPspectrum->SetName("JetProton");
  }
 
  Int_t nbins=fJetPspectrum->GetNbinsX();
  Int_t ibinlow=fJetPspectrum->FindFixBin(pmin);
  Int_t ibinup=fJetPspectrum->FindFixBin(pmax);
  // Only keep the histogram contents for the momentum range [pmin,pmax]
  for (Int_t i=1; i<=nbins; i++)
  {
   if (i<ibinlow || i>ibinup) fJetPspectrum->SetBinContent(i,0);
  }
 }
 
 // dN/dp spectrum specified by a function 
 if (fspec && !hspec && pmax>pmin)
 {
  fJetPscale=mode;
  TString title="Jet proton (beam) dN/dp=";
  title+=fspec->GetExpFormula("p");
  title.ReplaceAll("x","p");
  title+=" spectrum";
  if (mode==0) title+=";Momentum [GeV/c];Counts";
  if (mode==1)
  {
   title+=";^{10}Log(Momentum) [GeV/c];Counts";
   pmin=log10(pmin);
   pmax=log10(pmax);
  }
  if (mode==2)
  {
   title+=";Ln(Momentum) [GeV/c];Counts";
   pmin=log(pmin);
   pmax=log(pmax);
  }
  TH1F hpbeam=fLab.GetCountsHistogram(*fspec,1000,pmin,pmax,mode);
  fJetPspectrum=(TH1*)hpbeam.Clone();
  fJetPspectrum->SetTitle(title);
  fJetPspectrum->SetName("JetProton");
 }
}

void NcCollider::SetJetGammaSpectrum(Double_t pmin,Double_t pmax,TF1* fspec,TH1* hspec,Int_t mode)
{
 
 gROOT->cd(); // Make sure to work in memory
 
 fJetGpmin=0;
 fJetGpmax=0;
 
 if (fJetGspectrum)
 {
  delete fJetGspectrum;
  fJetGspectrum=0;
 }
 
 if (pmax<=pmin && pmin<=0)
 {
  printf(" *%-s::SetJetGammaSpectrum* Inconsistent input pmin=%-g pmax=%-g \n",ClassName(),pmin,pmax);
  return;
 }
 
 if (pmax>pmin && !fspec && !hspec)
 {
  printf(" *%-s::SetJetGammaSpectrum* Inconsistent input pmin=%-g pmax=%-g fspec=0 hspec=0 \n",ClassName(),pmin,pmax);
  return;
 }
 
 fJetGpmin=pmin;
 fJetGpmax=pmax;
 if (pmax<=pmin) fJetGpmax=pmin;
 fJetGscale=0;
 
 // Momentum distribution specified by a histogram
 if (hspec)
 {
  fJetGscale=mode;
  if (mode==1)
  {
   pmin=log10(pmin);
   pmax=log10(pmax);
  }
  if (mode==2)
  {
   pmin=log(pmin);
   pmax=log(pmax);
  }
 
  if (!fspec) // Histogram contains N vs. p distribution
  {
   fJetGspectrum=(TH1*)hspec->Clone();
  }
  else // Histogram contains dN/dp distribution scaled by "fspec"
  {
   TH1F hptarget=fLab.GetCountsHistogram(*hspec,mode,"",fspec);
   fJetGspectrum=(TH1*)hptarget.Clone();
   fJetGspectrum->SetName("JetGamma");
  }
 
  Int_t nbins=fJetGspectrum->GetNbinsX();
  Int_t ibinlow=fJetGspectrum->FindFixBin(pmin);
  Int_t ibinup=fJetGspectrum->FindFixBin(pmax);
  // Only keep the histogram contents for the momentum range [pmin,pmax]
  for (Int_t i=1; i<=nbins; i++)
  {
   if (i<ibinlow || i>ibinup) fJetGspectrum->SetBinContent(i,0);
  }
 }
 
 // dN/dp spectrum specified by a function 
 if (fspec && !hspec && pmax>pmin)
 {
  fJetGscale=mode;
  TString title="Jet gamma (target) dN/dp=";
  title+=fspec->GetExpFormula("p");
  title.ReplaceAll("x","p");
  title+=" spectrum";
  if (mode==0) title+=";Momentum [GeV/c];Counts";
  if (mode==1)
  {
   title+=";^{10}Log(Momentum) [GeV/c];Counts";
   pmin=log10(pmin);
   pmax=log10(pmax);
  }
  if (mode==2)
  {
   title+=";Ln(Momentum) [GeV/c];Counts";
   pmin=log(pmin);
   pmax=log(pmax);
  }
  TH1F hptarget=fLab.GetCountsHistogram(*fspec,1000,pmin,pmax,0);
  fJetGspectrum=(TH1*)hptarget.Clone();
  fJetGspectrum->SetTitle(title);
  fJetGspectrum->SetName("JetGamma");
 }
}

TH1* NcCollider::GetJetProtonSpectrum(Double_t* pmin,Double_t* pmax)
{
 
 if (pmin) *pmin=fJetPpmin;
 if (pmax) *pmax=fJetPpmax;
 
 return fJetPspectrum;
}

TH1* NcCollider::GetJetGammaSpectrum(Double_t* pmin,Double_t* pmax)
{
 
 if (pmin) *pmin=fJetGpmin;
 if (pmax) *pmax=fJetGpmax;
 
 return fJetGspectrum;
}

void NcCollider::ProcessJet(Double_t np,Double_t gfrac,TString flux,Double_t dthmax,Int_t nlist,Int_t ntrymax,Int_t wxsec,Double_t finit,Int_t full)
{
 
 if (fJetPpmax<=0 || fJetGpmax<=0 || np<1 || gfrac<0 || ntrymax<1)
 {
  printf(" *%-s::ProcessJet* Inconsistent intialisation. \n",ClassName());
  printf(" Momentum range for (beam) protons : [%-g,%-g] GeV/c \n",fJetPpmin,fJetPpmax);
  printf(" Momentum range for (target) gammas : [%-g,%-g] GeV/c \n",fJetGpmin,fJetGpmax);
  printf(" Number of simulated (beam) protons : %-g \n",np);
  printf(" Fraction of (beam) protons used for p+gamma interactions : %-g \n",gfrac);
  printf(" Maximum number of phase-space trials per event : %-i \n",ntrymax);
  return;
 }
 
 printf(" *%-s::ProcessJet* Parameter settings for astrophysical Jet simulation \n",ClassName());  
 printf(" Multiple partonic interactions flag : %-i \n",GetMultiple());
 printf(" Low-Pt, Elastic and Diffractive scattering flag : %-i \n",GetElastic());
 printf(" Minimal CMS energy for event generation : %-g GeV \n",GetEcmsMin());
 printf(" Number of simulated (beam) protons : %-g \n",np);
 printf(" Fraction of (beam) protons used for p+gamma interactions : %-g \n",gfrac);
 printf(" Maximum number of phase-space trials per event : %-i \n",ntrymax);
 printf(" Final particle species that are recorded : %-s \n",flux.Data()); 
 if (!fJetPspectrum)
 {
  printf(" Proton (beam) momenta will be mono-energetic at %-g GeV/c \n",fJetPpmax); 
 }
 else
 {
  printf(" Momentum range for (beam) protons : [%-g,%-g] GeV/c \n",fJetPpmin,fJetPpmax);
 }
 if (!fJetGspectrum)
 {
  printf(" Gamma (target) momenta will be mono-energetic at %-g GeV/c \n",fJetGpmax); 
 }
 else
 {
  printf(" Momentum range for (target) gammas : [%-g,%-g] GeV/c \n",fJetGpmin,fJetGpmax);
 }
 
 if (fMktree)
 {
  fMktree->Select("event","jrun");
  fMktree->Select("event","jevt");
  fMktree->Select("event","user","BeamP");
  fMktree->Select("event","user","BeamTheta");
  fMktree->Select("event","user","BeamPhi");
  fMktree->Select("event","user","TargetP");
  fMktree->Select("event","user","TargetTheta");
  fMktree->Select("event","user","TargetPhi");
 
  fMktree->Select("track","p");
 
  if (flux.Contains("nu"))
  {
   fMktree->UseTracks("nu_mu");
   fMktree->UseTracks("nu_mubar");
   fMktree->UseTracks("nu_e");
   fMktree->UseTracks("nu_ebar");
   fMktree->UseTracks("nu_tau");
   fMktree->UseTracks("nu_taubar");
   fMktree->UseTracks("nu",-1,1);
  }
  if (flux.Contains("gamma")) fMktree->UseTracks("gamma");
  if (flux.Contains("neutron"))
  {
   fMktree->UseTracks("n0");
   fMktree->UseTracks("nbar0");
   fMktree->UseTracks("p+");
   fMktree->UseTracks("pbar-");
  }
 }

 // Generate both p+p (jrun>0) and p+gamma (jrun<0) processes. //
 
 Int_t jrun=0;
 Nc3Vector pbeam;
 Nc3Vector ptarget;
 Nc3Vector pfixed;
 pfixed.SetVector(0,0,0,"car");
 Int_t ier=0;     // Flag to identify initialisation error
 Int_t nevents=0; // Number of events to be generated for each process
 Int_t igen=0;    // Flag to denote successful generation of an event
 Int_t ievt=0;    // Counter for successfully generated events
 Int_t ntry=0;    // Counter for phase-space trials per event
 Float_t beamp=fJetPpmax;
 Float_t targetp=fJetGpmax;
 
 if (fJetPspectrum || fJetGspectrum) wxsec=1;
 
 // Initialisation of the processes
 for (Int_t k=0; k<2; k++)
 {
  pbeam.SetVector(0,0,fJetPpmax,"car");
  ptarget.SetVector(0,0,-fJetGpmax,"car");
 
  if (!k) // p+p process
  {   
   nevents=(1.-gfrac)*np;
   jrun=1;
   ier=Init("free","p","p",0,&pbeam,&pfixed,wxsec,finit);
  }
  else // p+gamma process
  {
   nevents=gfrac*np;
   jrun=-1;
   ier=Init("free","p","gamma",0,&pbeam,&ptarget,wxsec,finit);
  }
 
  if (ier) return;
 
  SetRunNumber(jrun);
 
  // Define several particles as (un)stable according to selected analysis mode
  SetStable(0,1,4);                                                          // Declare all mesons as stable
  if (flux.Contains("nu") || flux.Contains("gamma")) SetStable(0,0,4);       // Declare all mesons as unstable
  if (!(flux.Contains("gamma")))  SetStable(111,1);                          // Declare pi0 as stable 
  if (!(flux.Contains("nu")))  SetStable(211,1);                             // Declare pi+ and pi- as stable 
  if (flux.Contains("nu")) SetStable(13,0);                                  // Declare mu+ and mu- as unstable
  if (flux.Contains("nu") && !(flux.Contains("neutron"))) SetStable(2112,0); // Declare n and nbar as unstable
 
  // Generation of the events for this process
  ievt=0;
  ntry=0;
  while (ievt<nevents && ntry<ntrymax)
  {
   // Pick a proton momentum from the beam momentum distribution 
   if (fJetPspectrum) beamp=fJetPspectrum->GetRandom();
   // Convert to linear scale momentum value
   if (fJetPscale==1) beamp=pow(10,beamp);
   if (fJetPscale==2) beamp=exp(beamp);
   pbeam.SetVector(0,0,beamp,"car");
   SetMomentum(pbeam,1);
 
   // Pick a photon momentum from the target momentum distribution 
   if (fJetGspectrum) targetp=fJetGspectrum->GetRandom();
   // Convert to linear scale momentum value
   if (fJetGscale==1) targetp=pow(10,targetp);
   if (fJetGscale==2) targetp=exp(targetp);
   ptarget.SetVector(0,0,-targetp,"car");
   SetMomentum(ptarget,2);
 
   // Randomisation of the beam or target direction
   if (dthmax>0)
   {
    fLab.RandomPosition(pbeam,0,dthmax,0,360);
    SetMomentum(pbeam,1);
   }
   if (dthmax<0)
   {
    fLab.RandomPosition(ptarget,180+dthmax,180,0,360);
    SetMomentum(ptarget,2);
   }
 
   // Fixed target for p+p events
   if (!k) SetMomentum(pfixed,2);
 
   if (nlist && ievt<nlist) // Produce Pythia listing for the first "nlist" events of each sample
   {
    if (!full)
    {
     igen=MakeEvent(0,1);
    }
    else
    {
     igen=MakeEvent(0,1,-1);
    }
   }
   else // No Pythia event listing
   {
    igen=MakeEvent();
   }
 
   if (igen<0) break;   // Error
   if (!igen) // Event did not pass selection
   {
    ntry++;
    continue;
   }
 
   ievt++;
   ntry=0;
  } // End of loop over the events for this process
 
  // Printout Pythia statistics for this event sample
  Pystat(1);
 } // End of loop over the processes
 
 EndRun();
}