//////////////////////////////////////////////////////////////////////////
//                                                                      //
// ATLFast MCMaker class.                                               //
//                                                                      //
// This class generates MC events or read events from a file.           //
//////////////////////////////////////////////////////////////////////////

#include <TFile.h>
#include <TTree.h>
#include <TPythia.h>
#include <TMCParticle.h>

#include "ATLFVirtualDisplay.h"
#include "ATLFMCMaker.h"
#include "ATLFast.h"

const Int_t kMAXPART = 10000;
const Double_t kPI = TMath::Pi();

#ifdef GENZ
# define genzinit genzinit_
# define genzread genzread_
# define genzend  genzend_
# define type_of_call
extern "C" void type_of_call genzinit();
extern "C" Int_t type_of_call genzread();
extern "C" void type_of_call genzend();
#endif


ClassImp(ATLFMCMaker)

//_____________________________________________________________________________
 ATLFMCMaker::ATLFMCMaker()
{
   m_Nparticles = 0;
   m_Tree   = 0;
   m_File   = 0;
   m_PT     = new Float_t[kMAXPART];
   m_PT2    = new Float_t[kMAXPART];
   m_Eta    = new Float_t[kMAXPART];
   m_Phi    = new Float_t[kMAXPART];
   m_Index  = new Int_t[kMAXPART];
   m_IsClonable = kFALSE;
}

//_____________________________________________________________________________
 ATLFMCMaker::ATLFMCMaker(const char *name, const char *title)
            :ATLFMaker(name,title)
{

   m_BranchName = "Particles";
   m_Nparticles = 0;
   m_Tree       = 0;
   m_File       = 0;
   m_IsClonable = kFALSE;

   SetHadronisationMode();
   SetProductionMode();
   SetRndNumber();

     // create arrays to store some intermediate particles parameters
   m_PT     = new Float_t[kMAXPART];
   m_PT2    = new Float_t[kMAXPART];
   m_Eta    = new Float_t[kMAXPART];
   m_Phi    = new Float_t[kMAXPART];
   m_Index  = new Int_t[kMAXPART];
}

//_____________________________________________________________________________
 ATLFMCMaker::~ATLFMCMaker()
{
   delete [] m_PT;
   delete [] m_PT2;
   delete [] m_Eta;
   delete [] m_Phi;
   delete [] m_Index;
}

//_____________________________________________________________________________
 void ATLFMCMaker::CleanParticles()
{
//    Clean particles array before writing to file
//   if m_Save = 1 (default), save only hard scattered particles
//   if m_Save = 2 save all particles (used by event display)

   if (m_Save != 1) return;
   TClonesArray *particles = gATLFast->MCMaker()->Fruits();
   particles->SetLast(m_Nstop);
}

//_____________________________________________________________________________
 void ATLFMCMaker::Clear(Option_t *option)
{
//    Reset Cell Maker

   m_Nparticles = 0;
   ATLFMaker::Clear(option);
}

//_____________________________________________________________________________
 void ATLFMCMaker::Draw(Option_t *)
{
//    Add itself to the list of graphics primitives in pad

   AppendPad();
}


//_____________________________________________________________________________
 void ATLFMCMaker::Init()
{
//.............................................
//  Initialize the Event generator
//  The current implementation supports two operational modes
//     A- Generation on the fly using Pythia with the following 10 modes
//        1- H-->gam gam
//        2- Z-->ee
//        3- Z-->mumu
//        4- H-->ZZ-->4e
//        5- H-->ZZ-->4mu
//        6- H-->ZZ-->2mu2e
//        7- WH, H-->bb, W-->lnu
//        8- WH, H-->gg, W-->lnu
//        9- WH, H-->jj, W-->lnu
//       10- bbA, A--> tau tau
//
//    B- Read a Root Tree file previously generated with one of the modes above
//      101- Read file pythia1.root with H-->gam gam events
//      102- Read file pythia2.root with Z-->ee events
//      .....
//      110  Read file pythia10.root with bbA, A--> tau tau events
//
//    C- Read a GENZ file
//      201
//
//   The operational mode is selected via ATLFMCMaker->SetProductionMode.
//.............................................

   TPythia *pythia;

#ifdef GENZ
//.............................................
//... Reading a GENZ file
//.............................................
   if (m_ProductionMode > 200) {
      genzinit();
      pythia = new TPythia();
      pythia->SetTitle("GENZ events");
      SetTitle("GENZ events");
      m_Generator = pythia;
      m_Fruits    = pythia->GetListOfPrimaries();
      return;
   }
#endif

//.............................................
//... Reading a Pythia file
//.............................................
   Text_t fname[20];
   if (m_ProductionMode > 100) {
      m_Fruits     = new TClonesArray("TMCParticle",kMAXPART, kFALSE);
      TDirectory *dirsav  = gDirectory;
      sprintf(fname,"pythia%i.root",m_ProductionMode%100);
        // Open the Root/Pythia file:
      m_File = new TFile(fname);

        // Get the Tree called "t" from the file
      m_Tree = (TTree*)m_File->Get("t");

     // Data will be read in in the standard ClonesArray
      m_Tree->SetBranchAddress("particles",m_Fruits);
      if (dirsav) dirsav->cd();
      return;
   }

//.............................................
//...  On the fly generation with Pythia
//.............................................

   pythia=new TPythia();
   m_Generator = pythia;
   m_Fruits    = pythia->GetListOfPrimaries();
   pythia->SetMSTP(125,1);   
   pythia->SetMSTP(51,9);
   pythia->SetMRLU(1,m_RndNumber);

   Float_t HiggsMass = 130.0;
   char ptitle[80];
   char *mode = "";
   Int_t i;

   switch(m_ProductionMode){
   
      case 1:
            mode = "H-->gam gam";
            pythia->SetMSEL(16);
            pythia->SetMDCY(25,1,0);
            for (i=192;i<=208;i++) pythia->SetMDME(i,1,0);
            pythia->SetMDME(205,1,1);
            HiggsMass = 100.0;
            break;
				
      case 2:
            mode = "Z-->ee";
            pythia->SetMSEL(0);  
            pythia->SetCKIN(1,80.0);
            pythia->SetMSUB(1,1);
            pythia->SetMDCY(23,1,0);
            for (i=156;i<=171;i++) pythia->SetMDME(i,1,0);
            pythia->SetMDME(164,1,1);
            break;
				
      case 3:
            mode = "Z-->mumu";
            pythia->SetMSEL(0);  
            pythia->SetCKIN(1,80.0);
            pythia->SetMSUB(1,1);
            pythia->SetMDCY(23,1,0);
            for (i=156;i<=171;i++) pythia->SetMDME(i,1,0);
            pythia->SetMDME(166,1,1);
            break;
				
      case 4: 
            mode = "H-->ZZ-->4e";
            pythia->SetMSEL(16);
            pythia->SetMDCY(25,1,0);
            for (i=192;i<=208;i++) pythia->SetMDME(i,1,0);
            pythia->SetMDME(207,1,1);
            pythia->SetMDCY(23,1,0);
            for (i=156;i<=171;i++) pythia->SetMDME(i,1,0);
            pythia->SetMDME(164,1,1);
            break;
				
      case 5:
            mode = "H-->ZZ-->4mu";
            pythia->SetMSEL(16);
            pythia->SetMDCY(23,1,0);
            pythia->SetMDCY(25,1,0);
            for (i=192;i<=208;i++) pythia->SetMDME(i,1,0);
            pythia->SetMDME(207,1,1);
            for (i=156;i<=171;i++) pythia->SetMDME(i,1,0);
            pythia->SetMDME(166,1,1);
            break;
				
      case 6: 
            mode = "H-->ZZ-->2mu2e";
            pythia->SetMSEL(16);
            pythia->SetMDCY(23,1,0);
            pythia->SetMDCY(25,1,0);
            for (i=192;i<=208;i++) pythia->SetMDME(i,1,0);
            pythia->SetMDME(207,1,1);
            for (i=156;i<=171;i++) pythia->SetMDME(i,1,0);
            pythia->SetMDME(164,1,5);
            pythia->SetMDME(166,1,4);
            break;
				
      case 7:
            mode = "WH, H-->bb, W-->lnu";
            pythia->SetMSEL(0);  
            pythia->SetMSUB(26,1);
            pythia->SetMDCY(24,1,0);
            for (i=172;i<=191;i++) pythia->SetMDME(i,1,0);
            pythia->SetMDME(188,1,1);
            pythia->SetMDME(189,1,1);
            pythia->SetMDCY(25,1,0);
            for (i=192;i<=208;i++) pythia->SetMDME(i,1,0);
            pythia->SetMDME(196,1,1);
            HiggsMass = 100;            
            break;
				
      case 8:
            mode = "WH, H-->gg, W-->lnu";
            pythia->SetMSEL(0);  
            pythia->SetMSUB(26,1);
            pythia->SetMDCY(24,1,0);
            for (i=172;i<=191;i++) pythia->SetMDME(i,1,0);
            pythia->SetMDME(188,1,1);
            pythia->SetMDME(189,1,1);
            pythia->SetMDCY(25,1,0);
            for (i=192;i<=208;i++) pythia->SetMDME(i,1,0);
            pythia->SetMDME(204,1,1);
            HiggsMass = 100;
            break;
				
      case 9:
            mode = "WH, H-->jj, W-->lnu";
            pythia->SetMSEL(0);  
            pythia->SetMSUB(26,1);
            pythia->SetMDCY(24,1,0);
            for (i=172;i<=191;i++) pythia->SetMDME(i,1,0);
            pythia->SetMDME(188,1,1);
            pythia->SetMDME(189,1,1);
            pythia->SetMDCY(25,1,0);
            for (i=192;i<=208;i++) pythia->SetMDME(i,1,0);
            pythia->SetMDME(192,1,1);
            pythia->SetMDME(193,1,1);
            HiggsMass = 100;
            break;
				
      case 10:
            mode = "bbA, A--> tau tau";
            pythia->SetMSEL(0);  
            pythia->SetMSUB(186,1);
            pythia->SetMSUB(187,1);
            pythia->SetKFPR(186,2,5);
            pythia->SetKFPR(187,2,5);
            pythia->SetPMAS(36,1,150.0);
            pythia->SetMDCY(36,1,0);
            for (i=274;i<=291;i++) pythia->SetMDME(i,1,0);
            pythia->SetMDME(284,1,1);
            for (i=83;i<=123;i++) pythia->SetMDME(i,1,2);
            pythia->SetMDME(83,1,3);
            pythia->SetMDME(84,1,3);
            HiggsMass = 150;
            break;
      case 11:
            mode = "WH, H-->bb, W-->lnu";
            pythia->SetMSEL(0);  
            pythia->SetMSUB(26,1);
            pythia->SetMDCY(24,1,0);
            for (i=172;i<=191;i++) pythia->SetMDME(i,1,0);
            pythia->SetMDME(188,1,1);
            pythia->SetMDME(189,1,1);
            pythia->SetMDCY(25,1,0);
            for (i=192;i<=208;i++) pythia->SetMDME(i,1,0);
            pythia->SetMDME(195,1,1);
            HiggsMass = 100;            
            break;

      default:
            printf("NO predefined ProductionMode selected!n");
            break;
   }
				
   switch(m_HadronisationMode){
   
      case 0:
            pythia->SetMSTP(61,0);
            pythia->SetMSTP(71,0);
            pythia->SetMSTP(81,0);
            pythia->SetMSTP(111,0);
            break;
      case 1:
            pythia->SetMSTP(71,0);
            pythia->SetMSTP(111,0);
            break;
      case 2:
            pythia->SetMSTP(111,0);
            break;
      case 3:
            break;
      case 4:
            pythia->SetMSTP(61,0);
            pythia->SetMSTP(81,0);
            pythia->SetMSTP(111,0);
            break;
      case 5:
            pythia->SetMSTP(61,0);
            pythia->SetMSTP(81,0);
            break;
                                
      default:
            printf("NO valid HadronisationMode selected!n");
            break;
   }
   pythia->SetMDCY(6,1,1);
   pythia->SetMSTP(48,1);

   pythia->SetPMAS(6,1,175.0); // Mass of top
   pythia->SetPMAS(25,1,HiggsMass);  // Mass of Higgs

   pythia->SetMSTJ(41,1); //only QCD showering

   pythia->SetMSTJ(11,3);  //Peterson fragmentation function for b and c:
   pythia->SetPARJ(54,-0.07);
   pythia->SetPARJ(55,-0.006);

// default (KEYHAD=3):
// no additional settings or switches required...


//Initializing Pythia:
   pythia->Initialize("CMS","p","p",14000);
   sprintf(ptitle,"%s: %s",pythia->GetTitle(),mode);
   pythia->SetTitle(ptitle);
   SetTitle(mode);
}

//_____________________________________________________________________________
 void ATLFMCMaker::Finish()
{
// Function called when all makers for one event have been called

#ifdef GENZ
   if (m_ProductionMode > 200) {
      genzend_();
   }
#endif
}

//_____________________________________________________________________________
 Int_t ATLFMCMaker::Make()
{
//.............................................
//  This function generates one MonteCarlo event
//  Read from a file if production mode > 100
//  else generation on the fly


   TPythia *pythia;
   if (m_ProductionMode > 200) {
#ifdef GENZ
//.............................................
//... Reading a GENZ file
//.............................................
      Int_t ret = genzread();
      if (ret < 0) return ret;
      pythia = (TPythia*)m_Generator;
      pythia->GetPrimaries();
#endif
   } else if (m_ProductionMode > 100) {
      Int_t event = gATLFast->Event();
      m_Tree->GetEvent(event);

   } else {
      pythia = (TPythia*)m_Generator;
      pythia->GenerateEvent();
   }

   TClonesArray *particles = (TClonesArray*)m_Fruits;
   m_Nparticles = particles->GetEntriesFast();
   TMCParticle *part;
   if (m_Nparticles >= kMAXPART) {
      Error("Make","Too many particles generated=%d, buffer size=%d",m_Nparticles,kMAXPART);
   }

//  Find last original parton number
   Int_t i;
   for (i=0;i<m_Nparticles;i++) {
      part = (TMCParticle*)particles->UncheckedAt(i);     
      if (part->GetKS() != 21) {
         m_Nstop  = i;
         m_Nstart = i+1;
         break;
      }
   }

//  Fill additional arrays of particles parameters
   m_Nstables = 0;
   Bool_t compute;
   for (i=0;i<m_Nparticles;i++) {
      part = (TMCParticle*)particles->UncheckedAt(i);     
      if (i < m_Nstop) compute = kTRUE;
      else             compute = kFALSE;
      if (part->GetKS() >=0 && part->GetKS() < 11) {
         m_Index[m_Nstables] = i;
         m_Nstables++;
         compute = kTRUE;
      }
      if (compute) {
         m_PT2[i] = part->GetPx()*part->GetPx() + part->GetPy()*part->GetPy();
         m_PT[i]  = TMath::Sqrt(m_PT2[i]);
         m_Eta[i] = Rapidity(m_PT[i], part->GetPz());
         m_Phi[i] = Angle(part->GetPx(), part->GetPy());
      }
   }
   return 0;
}

//_____________________________________________________________________________
 void ATLFMCMaker::Paint(Option_t *option)
{
//    Paint generated particles

   gATLFast->Display()->PaintParticles(option);
}

//______________________________________________________________________________
 void ATLFMCMaker::Sizeof3D() const
{
   gATLFast->Display()->SizeParticles();
}

//_____________________________________________________________________________
 void ATLFMCMaker::PrintInfo()
{
   ATLFMaker::PrintInfo();
}

//_____________________________________________________________________________
 Float_t ATLFMCMaker::Angle(Float_t x, Float_t y)
{
//   Compute phi angle of particle
// ... this is a copy of function ULANGL
//  .. sign(a,b) = -abs(a) if b <  0
//               =  abs(a) if b >= 0

   Float_t angle = 0;
   Float_t r = TMath::Sqrt(x*x + y*y);
   if (r < 1e-20) return angle;
   if (TMath::Abs(x)/r < 0.8) {
      angle = TMath::Sign((Float_t)TMath::Abs(TMath::ACos(x/r)), y);
   } else {
      angle = TMath::ASin(y/r);
      if (x < 0 ) {
         if(angle >= 0) angle =  kPI - angle;
         else           angle = -kPI - angle;
      }
   }
   return angle;
}

//_____________________________________________________________________________
 Int_t ATLFMCMaker::Charge(Int_t kf)
{
//... this is copy of function LUCHGE 
//...Purpose: to give three times the charge for a particle/parton. 

  static Int_t kchg[500] = { -1,2,-1,2,-1,2,-1,2,0,0,-3,0,-3,0,-3,0,-3,0,
        0,0,0,0,0,3,0,0,0,0,0,0,0,0,0,3,0,0,3,0,-1,0,0,0,0,0,0,0,0,0,0,0,
       0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
       2,-1,2,-1,2,3,0,0,0,0,0,0,0,0,0,0,0,3,0,3,3,0,3,0,3,0,3,0,0,0,0,0,
       0,0,0,0,0,3,0,3,3,0,3,0,3,0,3,0,0,0,0,0,0,0,0,0,0,3,0,3,3,0,3,0,3,
       0,3,0,0,0,0,0,0,0,0,0,0,3,0,3,3,0,3,0,3,0,3,0,0,0,0,0,0,0,0,0,0,3,
       0,3,3,0,3,0,3,0,3,0,0,0,0,0,0,0,0,0,0,3,0,3,3,0,3,0,3,0,3,0,0,0,0,
       0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
       0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
       3,0,3,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,3,0,
       0,3,0,0,0,0,0,0,0,0,-3,0,0,0,0,0,0,0,0,3,0,-3,0,3,-3,0,0,0,3,6,0,
       3,0,0,0,0,0,-3,0,3,-3,0,-3,0,0,0,0,-3,0,3,6,-3,0,3,-3,0,-3,0,3,6,
       0,3,0,0,0,0,0,-3,0,3,-3,0,-3,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
       0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
       0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
       0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 };

//    extern integer kfcomp_(integer *);
  Int_t ipower;
  Int_t ret = 0;
  Int_t kfa = TMath::Abs(kf);
  Int_t kc  = Compress(kfa);

//...Initial values. Simple case of direct readout. 
  if (kc == 0) {
  } else if (kfa <= 100 || kc <= 80 || kc > 100) {
     ret = kchg[kc-1];

// ...Construction from quark content for heavy meson, diquark, baryon.
  } else if (kfa/1000 % 10 == 0) {
	ipower = kfa/100 % 10;
	ret = (kchg[kfa / 100 % 10 - 1] - kchg[kfa/10 % 10 - 1])*Int_t(TMath::Power(-1, ipower));
  } else if (kfa / 10 % 10 == 0) {
	ret = kchg[kfa/1000 % 10 - 1] + kchg[kfa/100 % 10 - 1];
  } else {
	ret = kchg[kfa/1000 % 10 - 1] + kchg[kfa/100 % 10 - 1] + kchg[kfa/10 % 10 - 1];
  }

// ...Add on correct sign.
  if (kf > 0) return ret;
  else        return -ret;
}

//_____________________________________________________________________________
 Int_t ATLFMCMaker::Compress(Int_t kf)
{
//... this is copy of function LUCOMP 
//...Purpose: to compress the standard KF codes for use in mass and decay 
//...arrays; also to check whether a given code actually is defined.
//     from BLOCK LUDATA
  static Int_t  kchg[500] = { 1,1,1,1,1,1,1,1,0,0,1,1,1,1,1,1,1,1,0,0,0,0,
        0,1,0,0,0,0,0,0,0,0,0,1,0,0,1,0,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
        0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,1,1,1,1,1,
        1,1,0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,1,1,1,0,0,0,0,0,0,0,0,0,0,1,
        1,1,1,1,1,1,1,1,1,0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,1,1,1,0,0,0,0,
        0,0,0,0,0,0,1,1,1,1,1,1,1,1,1,1,0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,
        1,1,1,0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,1,1,1,0,0,0,0,0,0,0,0,0,0,
        0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
        0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,0,0,0,
        0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,1,0,1,1,0,0,0,0,
        0,0,1,0,1,1,0,0,0,0,0,0,1,1,0,1,1,1,1,1,0,1,1,1,1,1,1,0,0,0,0,1,1,
        1,1,1,1,0,0,0,0,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,0,0,0,0,1,1,1,1,1,
        1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
        0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
        0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
        0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 };
  static Int_t kftab[25] = { 211,111,221,311,321,130,310,213,113,223,313,
	    323,2112,2212,210,2110,2210,110,220,330,440,30443,30553,0,0 };
  static Int_t kctab[25] = { 101,111,112,102,103,221,222,121,131,132,122,
	    123,332,333,281,282,283,284,285,286,287,231,235,0,0 };

  Int_t ret = 0;
  Int_t kfla, kflb, kflc, kflr, kfls, kfa, ikf;

  kfa = TMath::Abs(kf);
//...Simple cases: direct translation or table.
  if (kfa == 0 || kfa >= 100000) {
     return ret;
  } else if (kfa <= 100) {
     ret = kfa;
     if (kf < 0 && kchg[kfa - 1] == 0)  ret = 0;
     return ret;
  } else {
     for (ikf = 1; ikf <= 23; ++ikf) {
        if (kfa == kftab[ikf-1]) {
           ret = kctab[ikf-1];
           if (kf < 0 && kchg[ret-1] == 0)  ret = 0;
           return ret;
        }
     }
  }
// ...Subdivide KF code into constituent pieces.
  kfla = kfa / 1000%10;
  kflb = kfa / 100%10;
  kflc = kfa / 10%10;
  kfls = kfa%10;
  kflr = kfa / 10000%10;
// ...Mesons.
  if (kfa - kflr*10000 < 1000) {
     if (kflb == 0 || kflb == 9 || kflc == 0 || kflc == 9) {
     } else if (kflb < kflc) {
     } else if (kf < 0 && kflb == kflc) {
     } else if (kflb == kflc) {
        if (kflr == 0 && kfls == 1) {        ret = kflb + 110;
        } else if (kflr == 0 && kfls == 3) { ret = kflb + 130;
        } else if (kflr == 1 && kfls == 3) { ret = kflb + 150;
        } else if (kflr == 1 && kfls == 1) { ret = kflb + 170;
        } else if (kflr == 2 && kfls == 3) { ret = kflb + 190;
        } else if (kflr == 0 && kfls == 5) { ret = kflb + 210;
        }
     } else if (kflb <= 5) {
        if (kflr == 0 && kfls == 1) {        ret = (kflb-1)*(kflb-2)/2 + 100 + kflc;
        } else if (kflr == 0 && kfls == 3) { ret = (kflb-1)*(kflb-2)/2 + 120 + kflc;
        } else if (kflr == 1 && kfls == 3) { ret = (kflb-1)*(kflb-2)/2 + 140 + kflc;
        } else if (kflr == 1 && kfls == 1) { ret = (kflb-1)*(kflb-2)/2 + 160 + kflc;
        } else if (kflr == 2 && kfls == 3) { ret = (kflb-1)*(kflb-2)/2 + 180 + kflc;
        } else if (kflr == 0 && kfls == 5) { ret = (kflb-1)*(kflb-2)/2 + 200 + kflc;
        }
     } else if (kfls == 1 && kflr <= 1 || kfls == 3 && kflr <= 2 || kfls == 5 && kflr == 0) {
        ret = kflb + 80;
     }
// ...Diquarks.
  } else if ((kflr == 0 || kflr == 1) && kflc == 0) {
     if (kfls != 1 && kfls != 3) {
     } else if (kfla == 9 || kflb == 0 || kflb == 9) {
     } else if (kfla < kflb) {
     } else if (kfls == 1 && kfla == kflb) {
     } else { ret = 90;
     }
// ...Spin 1/2 baryons.
  } else if (kflr == 0 && kfls == 2) {
     if (kfla == 9 || kflb == 0 || kflb == 9 || kflc == 9) {
     } else if (kfla <= kflc || kfla < kflb) {
     } else if (kfla >= 6 || kflb >= 4 || kflc >= 4) {
         ret = kfla + 80;
     } else if (kflb < kflc) {
         ret = (kfla+1)*kfla*(kfla-1)/6 + 300 + kflc*(kflc-1)/2 + kflb;
     } else {
         ret = (kfla+1)*kfla*(kfla-1)/6 + 330 + kflb*(kflb-1)/2 + kflc;
     }
// ...Spin 3/2 baryons.
  } else if (kflr == 0 && kfls == 4) {
     if (kfla == 9 || kflb == 0 || kflb == 9 || kflc == 9) {
     } else if (kfla < kflb || kflb < kflc) {
     } else if (kfla >= 6 || kflb >= 4) {
         ret = kfla + 80;
     } else {
         ret = (kfla+1)*kfla*(kfla-1) / 6 + 360 + kflb*(kflb -1) / 2 + kflc;
     }
  }
    return ret;
}

//_____________________________________________________________________________
 Float_t ATLFMCMaker::Rapidity(Float_t pt, Float_t pz)
{
//   Compute rapidity

   Float_t etalog = TMath::Log((TMath::Sqrt(pt*pt + pz*pz) + TMath::Abs(pz))/pt);
   if (pz < 0 ) return -TMath::Abs(etalog);
   else         return  TMath::Abs(etalog);
}