NuMISpillTimeStructure.cxx
Go to the documentation of this file.
1 // #include "EventGenerator/GENIE/NuMISpillTimeStructure.h"
3 #include <cassert>
4 #include <iostream>
5 #include "TMath.h"
6 #include "TRandom.h"
7 using namespace evgen;
8 
10  //
11  // These numbers specify the spill structure for NuMI
12  //
13  // * There are 84 RF buckets in a single batch
14  // * There are 6 filled buckets in a single spill
15  // * Within a batch, the first 81 buckets are filled
16  //
17  // * 53.103 MHz is the RF frequency of Main Injector at 120 GeV
18  //
19  // * The RF sine wave localizes the protons to the center of the
20  // bucket with a ~gaussian time profile of width ~750 ps (e-mail
21  // from Phil Adamson (May 28, 2015)
22  //
23  // * There is some variation bucket-to-bucket in the number of
24  // protons. This is estimated to be about +/- 15% with a gasussian
25  // spread based on a quick look at some of the data from the
26  // resistive wall monitors on the fast NuMI DAQ
27  //
28  // messier@indiana.edu (June 1, 2015)
29  //
30  fNbatch (6),
31  fNbucketPerBatch(84),
32  fNfilledBucket (81),
33  fTbucket (1e9/53.103E6),
34  fBucketSigmaT (0.750),
35  fBucketSigmaPOT (0.15),
36  fNbinPerBucket (50),
37  fNbin (2*fNbinPerBucket+fNbucketPerBatch*fNbatch*fNbinPerBucket),
38  fPOTbyBatch (fNbatch),
39  fSpillPOT (fNbin),
40  fSpillPOTI (fNbin)
41 {
42  fPOTbyBatch[0] = 1.9;
43  fPOTbyBatch[1] = 1.9;
44  fPOTbyBatch[2] = 1.0;
45  fPOTbyBatch[3] = 1.0;
46  fPOTbyBatch[4] = 1.0;
47  fPOTbyBatch[5] = 1.0;
48 
49  for (unsigned int i=0; i<fSpillPOT. size(); ++i) fSpillPOT[i] = 0;
50  for (unsigned int i=0; i<fSpillPOTI.size(); ++i) fSpillPOTI[i] = 0;
51 }
52 
53 //.....................................................................
54 
55 double NuMISpillTimeStructure::BinToT(unsigned int i) const
56 {
57  return fTbucket*((float)i-(float)fNbinPerBucket+0.5)/(float)fNbinPerBucket;
58 }
59 
60 //.....................................................................
61 
63 {
64  unsigned int i, j, k;
65  unsigned int it, it0;
66  double t, t0;
67  double pot, fpot;
68 
69  //
70  // Figure out the relative number of protons in each time bin
71  //
72  it = fNbinPerBucket;
73  for (i=0; i<fNbatch; ++i) {
74  pot = fPOTbyBatch[i];
75  for (j=0; j<fNbucketPerBatch; ++j) {
76  it0 = it+fNbinPerBucket/2;
77  t0 = this->BinToT(it0);
78  fpot = gRandom->Gaus(1.0, fBucketSigmaPOT);
79  for (k=0; k<fNbinPerBucket; ++k) {
80  if (j>fNfilledBucket) {
81  fSpillPOT[it] = 0.0;
82  }
83  else {
84  t = this->BinToT(it);
85  fSpillPOT[it] = pot*fpot*TMath::Gaus(t-t0,0,fBucketSigmaT);
86  }
87  ++it;
88  assert(it<=fNbin);
89  }
90  }
91  }
92 
93  //
94  // Now compute the integral and normalize it to 1 so we can sample
95  // the timing distribution
96  //
97  fSpillPOTI[0] = fSpillPOT[i];
98  for (i=1; i<fSpillPOTI.size(); ++i) {
99  fSpillPOTI[i] = fSpillPOTI[i-1] + fSpillPOT[i];
100  }
101  for (i=0; i<fSpillPOTI.size(); ++i) {
102  fSpillPOTI[i] /= fSpillPOTI[fSpillPOTI.size()-1];
103  }
104 }
105 
106 //......................................................................
107 
109 {
110  unsigned int i;
111  double r;
112  r = gRandom->Uniform();
113  for (i=0; i<fSpillPOTI.size(); ++i) {
114  if (fSpillPOTI[i]>r) return this->BinToT(i);
115  }
116  return 0.0;
117 }
118 
119 ////////////////////////////////////////////////////////////////////////
unsigned int fNfilledBucket
Number of buckets that are filled.
std::vector< double > fSpillPOT
POT/bin.
set< int >::iterator it
unsigned int fNbucketPerBatch
Number of RF buckets per batch.
std::vector< double > fSpillPOTI
Integrated POT.
unsigned int fNbatch
Number of proton batches.
std::vector< double > fPOTbyBatch
Relative number of POT in a batch (arb.)
unsigned int fNbinPerBucket
Number of time bins to use per bucket.
unsigned int fNbin
Total number of bins in time profile.
#define pot
double Gaus(TH1D *h, double &err, bool isTruth)
Definition: absCal.cxx:489
const double j
Definition: BetheBloch.cxx:29
double fBucketSigmaT
Width of protons inside an RF time bucket (ns)
double fBucketSigmaPOT
bucket-to-bucket variation in POT (fractional)
assert(nhit_max >=nhit_nbins)
TRandom3 r(0)
double BinToT(unsigned int i) const
double fTbucket
How many ns are there per bucket?
Module to generate only pions from cosmic rays.