Core/test/calcnattempts.C

0001 #include "TGraph.h"
0002 #include "TCanvas.h"
0003 #include "TMath.h"
0004
0005 //J.Bendavid
0006 //Script to optimize the number of hadronization steps for a given matching*filter efficiency
0007 //Optimization is exact under the assumption that every input event has the same underlying matching*filter efficiency
0008 //In practice this is not true for jet matching when mixing multiplicities, so the optimization becomes approximate.
0009 //This assumption is also not strictly true in case gen filters include pt/acceptance cuts which correlate the filter
0010 //efficiency with the parton kinematics from the hard event.
0011
0012 void calcnattempts() {
0013
0014   const double eff=0.273*0.00122; // (54924/201000)*(67/54924)
0015   const unsigned int ntmax=1000;
0016
0017   const double tlhe = 2.8;  //cpu time to generate an lhe event
0018   const double tps = 19.5e-3;  //cpu time to shower/hadronize an event
0019   const double ts = 70.6;   //cpu time to simulate an event
0020
0021   double tmin = (1./eff)*(tlhe+tps) + ts;
0022   int ntmin = 1;
0023   double ninmin = 0.;
0024   double noutmin = 0.;
0025
0026   printf("initial tmin = %5f\n",tmin);
0027
0028   TGraph *htime = new TGraph;
0029   for (int nt=1; nt<=ntmax; ++nt) {
0030     double pk0 = pow(1.-eff,nt);
0031     //double pk0 = TMath::Binomial(nt,0)*pow(eff,0)*pow(1.-eff,nt-0);
0032     double r1 = 0.;
0033     double r2 = 0.;
0034     bool trippednan = false;
0035     for (int k=0; k<=nt; ++k) {
0036       double pk = 0;
0037       if (!trippednan) {
0038     pk = TMath::Binomial(nt,k)*pow(eff,k)*pow(1.-eff,nt-k);
0039     if (std::isnan(pk)) trippednan = true;
0040     //assert(!trippednan);
0041       }
0042       if (trippednan) {
0043     //printf("binomialcoeff = %5f\n", TMath::Binomial(nt,k));
0044     double lambda = (double)nt*eff;
0045     pk = TMath::PoissonI(k,lambda);
0046       }
0047       r1 += pk*k*k;
0048       r2 += pk*k;
0049     }
0050     double r = r1/(r2*r2);
0051     double nin = r;
0052     double nout = r*(1.-pk0);
0053     double tcpu = r*(tlhe + nt*tps + (1.-pk0)*ts);
0054     //double tcpu = r*(tlhe + nt*tps + ts);
0055     if (tcpu<tmin) {
0056       tmin = tcpu;
0057       ntmin = nt;
0058       ninmin = nin;
0059       noutmin = nout;
0060     }
0061     printf("pk0 = %5f, r = %5f, r1 = %5f, r2 = %5f\n",pk0, r,r1,r2);
0062     printf("nt = %i, NE = 1, Nin = %5f, Nout = %5f, tcpu = %5f\n",nt,nin,nout,tcpu);
0063     htime->SetPoint(nt-1,double(nt),tcpu);
0064   }
0065   htime->Draw();
0066
0067   printf("Optimal nAttempts = %i, cpu time per unweighted event equivalent = %5f\n",ntmin, tmin);
0068   printf("NE = 1, Nin = %5f, Nout = %5f\n",ninmin,noutmin);
0069   printf("Nout/Nin = %5f\n",noutmin/ninmin);
0070   printf("NE/Nin (equivalent unweighted events per input event) = %5f\n",1./ninmin);
0071
0072 }