fillSpaceChargeMaps/macros/simulateTimestamps.C

0001 #include <TFile.h>
0002 #include <TTree.h>
0003
0004 #include <gsl/gsl_randist.h>
0005 #include <gsl/gsl_rng.h>
0006
0007 #include <Rtypes.h>
0008
0009 #include <format>
0010 #include <fstream>
0011 #include <iostream>
0012
0013 R__LOAD_LIBRARY(libgslcblas.so)
0014 R__LOAD_LIBRARY(libgsl.so)
0015
0016 //_______________________________________________________________
0017 void simulateTimestamps()
0018 {
0019   // random generator
0020   auto *rng = gsl_rng_alloc(gsl_rng_mt19937);
0021   // gsl_rng_set( rng, time(0) );//time(0) is the seed
0022   gsl_rng_set(rng, 1);
0023
0024   // time interval (s) between two bunch crossing
0025   /* value copied from generators/phhepmc/Fun4AllHepMCPileupInputManager.cc */
0026   static constexpr double deltat_crossing = 106e-9;
0027
0028   // triggered rate (Hz)
0029   // static constexpr double trigger_rate = 1e5; // Hijing 100kHz
0030   static constexpr double trigger_rate = 5e4;  // Hijing 50kHz
0031   // static constexpr double trigger_rate = 3e6; // pp collisions
0032
0033   // mean number of collision per crossing
0034   static constexpr double mu = trigger_rate * deltat_crossing;
0035
0036   // print configuration
0037   std::cout << "SimulateCollisions - deltat_crossing: " << deltat_crossing << std::endl;
0038   std::cout << "SimulateCollisions - trigger_rate: " << trigger_rate << std::endl;
0039   std::cout << "SimulateCollisions - mu: " << mu << std::endl;
0040
0041   // number of requested triggers (should correspond approximately to 1/10 second of data taking)
0042   // static constexpr uint ntrigtot = 1e4;
0043   // static constexpr uint ntrigtot = 6e6;
0044   static constexpr uint ntrigtot = 1e6;
0045   // write all timestamps to a file
0046   std::ofstream out("./data/timestamps_50kHz_1M.txt");
0047   out << "// bunchcrossin id; time (ns)" << std::endl;
0048   out << "// assuming 106ns between bunches" << std::endl;
0049
0050   // running collision time
0051   int64_t bunchcrossing = 0;
0052   float time = 0;
0053   float timeD = 0;
0054   TTree *timestamps = new TTree("timestamps", "beamcrossings timestamps");
0055   timestamps->Branch("bc", &bunchcrossing, "bc/I");
0056   timestamps->Branch("t", &time, "t/F");
0057   timestamps->Branch("dt", &timeD, "dt/F");
0058
0059   // keep track of the last collision time
0060   double previous_trigger_time = 0;
0061
0062   // generate triggers
0063   for (int itrig = 0; itrig < ntrigtot;)
0064   {
0065     ++bunchcrossing;
0066     time += deltat_crossing;
0067
0068     auto ntrig = gsl_ran_poisson(rng, mu);
0069     for (uint i = 0; i < ntrig; ++i)
0070     {
0071       out << " " << static_cast<int>(bunchcrossing) << " " << static_cast<int>(time * 1e9) << std::endl;
0072       timeD = (time - previous_trigger_time) * 1e9;
0073       timestamps->Fill();
0074       previous_trigger_time = time;
0075       ++itrig;
0076       if (itrig % 100000 == 0)
0077       {
0078         std::cout << "SimulateCollisions - itrig: " << itrig << std::endl;
0079       }
0080     }
0081   }
0082
0083   // printout last trigger time
0084   std::cout << "SimulateCollisions - last trigger time: " << time << std::endl;
0085
0086   out.close();
0087   // out.close();
0088   gsl_rng_free(rng);
0089
0090   TFile *outf = new TFile("./data/timestamps_50kHz_1M_t.root", "recreate");
0091   timestamps->Write();
0092
0093   outf->Write();
0094   outf->Close();
0095 }