/////////////////////////////////////////////////////////////////////
//                                                                 //
// This program is used to calculate number of events involving    //
// B+ and B- mesons in the same event at LHC conditions            //
//                                                                 //
// To run at background:                                           //
//  $ ./ex3 >logfile.txt 2>logerr.txt &                            //
//                                                                 //
// July 2025                                                       //
/////////////////////////////////////////////////////////////////////

#include <iostream>
#include <ctime>
#include "Pythia8/Pythia.h"
#include "TThread.h"

using namespace std;
using namespace Pythia8;

//
// *** handle function for multi-threading ***
//
void *handle(void *ptr)
{
  int ith = (long) ptr;

  Pythia8::Pythia pythia;

  // Read command file
  pythia.readFile("ex3.cmnd");

  // We will use 4 threads each will start with different random seed
  pythia.readString("Random:setSeed = on");  // set random seed
  pythia.readString(Form("Random:seed = %d",ith));  // set seed value

  // inilize pythia
  pythia.init();

  int nevents = 5000, totalBevent = 0;

  cout << "Counting B+ and B- particles..." << endl;

  // Main event loop ------------------------------------
  for(int i=0; i<nevents; i++)
  {
    if (!pythia.next()) continue;

    int entries = pythia.event.size();
    int Bplus=0, Bminus=0;

    for(int j=0; j<entries; j++)
    {
       int    id = pythia.event[j].id();
       double px = pythia.event[j].px();
       double py = pythia.event[j].py();
       double pT = sqrt(px*px+py*py);
       if(pT<5.0) continue; // transverse momentum cut

       if(id ==  521) Bplus++;  // count B+ and
       if(id == -521) Bminus++; // count B-
    }

    if(Bplus >=1 && Bminus >=1) totalBevent++;
  }
  // ----------------------------------------------------

  cout << scientific;
  cout <<"Thread # " << ith
       <<": Total number of events with B+ B- " << totalBevent << " / " << nevents
       << " = " << double(totalBevent)/nevents << endl;

  return ptr;
}

//
// *** main program ***
//
int main()
{
  const int nt = 4; // number of threads
  TThread *th[nt];

  for(int i=0;i<nt;i++)
  {
     th[i] = new TThread(Form("th%d",i), handle, (void*)i);
     th[i]->Run();
  }

  for(int i=0; i<nt; i++) th[i]->Join();

  cout << "End of program\n";

  return 0;
}