ProtoDUNEBeamInstrument.cxx
Go to the documentation of this file.
1 //// Created and Modified by Pablo and Leigh H. Howard,
2 ////Smear important variables of beam monitors, mimic Cherenkov monitors response
3 ////and store some histograms through art utilities
4 //// pablo.fer@cern.ch
5 //// July 2018
6 /////////////////////////////////////////////////////////////
7 
8 #include <iostream>
10 
11 namespace sim{
12 
13  //--------------------------constructors-----------------------------------------------
15  fInstrumentName = "default";
16  fX = 0.0;
17  fY = 0.0;
18  fZ = 0.0;
19  fT = 0.0;
20  fPx = 0.0;
21  fPy = 0.0;
22  fPz = 0.0;
23  fPDGid = 0;
24  fEventID = 0;
25  fTrackID = 0;
26  fSmearedVar1 = 0.0;
27  fSmearedVar2 = 0.0;
28  fResolution = 0.0;
29  }
30 
31  //-------------------------------default destructor------------------------------------------
33 
34  }
35 
37  Float_t x,
38  Float_t y,
39  Float_t z,
40  Float_t t,
41  Float_t Px,
42  Float_t Py,
43  Float_t Pz,
44  Int_t PDGid,
45  Int_t EventID,
46  Int_t TrackID,
47  Float_t Resolution
48  ){
50  fX = x; fY = y; fZ = z; fT = t;
51  fPx = Px; fPy = Py; fPz = Pz;
52  fPDGid = PDGid;
53  fEventID = EventID;
54  fTrackID = TrackID;
55  fResolution = Resolution;
56  srand (static_cast <unsigned> (time(NULL)));
57  if (name == "TOF1" || name == "TRIG2"){
58  Float_t delta_t = 20.;
59 // if (name == "TOF1") srand (static_cast <unsigned> (time(0)*9));
60  Float_t t_test = t - delta_t/2. + static_cast <float> (rand()) /( static_cast <float> (RAND_MAX/(delta_t)));
61  Float_t p_test = static_cast <float> (rand()) /( static_cast <float> (RAND_MAX));
62  Float_t p_gauss = ProtoDUNEBeamInstrument::UnitGauss(t,t_test,fResolution);
63  while (p_test > p_gauss){
64  p_test = static_cast <float> (rand()) /( static_cast <float> (RAND_MAX));
65  t_test = t -delta_t/2. + static_cast <float> (rand()) /( static_cast <float> (RAND_MAX/(delta_t)));
67  }
68  fSmearedVar1 = t_test;
69  }
70  if (name == "BPROFEXT" || name == "BPROF4"){
71  Float_t delta_x = 40.;
72 // srand (static_cast <unsigned> (time(0)*99));
73  Float_t x_test = x -delta_x/2. + static_cast <float> (rand()) /( static_cast <float> (RAND_MAX/(delta_x)));
74  Float_t p_test = static_cast <float> (rand()) /( static_cast <float> (RAND_MAX));
75  Float_t p_gauss = ProtoDUNEBeamInstrument::UnitGauss(x,x_test,fResolution);
76  while (p_test > p_gauss){
77  p_test = static_cast <float> (rand()) /( static_cast <float> (RAND_MAX));
78  x_test = x -delta_x/2. + static_cast <float> (rand()) /( static_cast <float> (RAND_MAX/(delta_x)));
80  }
81  fSmearedVar1 = x_test;
82 
83  Float_t delta_y = 40.;
84 // srand (static_cast <unsigned> (time(0)*7));
85  Float_t y_test = y -delta_y/2. + static_cast <float> (rand()) /( static_cast <float> (RAND_MAX/(delta_y)));
86  p_test = static_cast <float> (rand()) /( static_cast <float> (RAND_MAX));
88  while (p_test > p_gauss){
89  p_test = static_cast <float> (rand()) /( static_cast <float> (RAND_MAX));
90  y_test = y -delta_y + static_cast <float> (rand()) /( static_cast <float> (RAND_MAX/(delta_y)));
92  }
94 }
95  if (name == "CHERENKOV1"){
96 // srand (static_cast <unsigned> (time(0)*77));
97  Float_t p_test = static_cast <float> (rand()) /( static_cast <float> (RAND_MAX));
98  fSmearedVar1 = 0;
99  if (p_test <= fResolution){
100  Float_t Ptot = pow(pow(Px,2)+pow(Py,2)+pow(Pz,2),0.5)/1000.;
101  if (Ptot <= 2.0){
102  if (abs(PDGid) == 2212 || abs(PDGid) == 310) fSmearedVar1 = 0;
103  if (abs(PDGid) == 11) fSmearedVar1 = 1;
104  if (abs(PDGid) == 211) fSmearedVar1 = 1;
105 }
106  if (Ptot <= 3.0 && Ptot >2.0){
107  if (abs(PDGid) == 2212 || abs(PDGid) == 310) fSmearedVar1 = 0;
108  if (abs(PDGid) == 11) fSmearedVar1 = 1;
109  if (abs(PDGid) == 211) fSmearedVar1 = 1;
110 }
111  if (Ptot <= 5.0 && Ptot > 3.0){
112  if (abs(PDGid) == 2212 || abs(PDGid) == 310) fSmearedVar1 = 0;
113  if (abs(PDGid) == 11) fSmearedVar1 = 1;
114  if (abs(PDGid) == 211) fSmearedVar1 = 1;
115 }
116  if (Ptot > 5.0){
117  if (abs(PDGid) == 2212 || abs(PDGid) == 310) fSmearedVar1 = 0;
118  if (abs(PDGid) == 11) fSmearedVar1 = 0;
119  if (abs(PDGid) == 211) fSmearedVar1 = 1;
120 }
121 }
122 std::cout << "Particle ID: " << fPDGid << std::endl;
123 std::cout << "Cherenkov 1: " << fSmearedVar1 << std::endl;
124 }
125  if (name == "CHERENKOV2"){
126 // srand (static_cast <unsigned> (time(0)*97));
127  Float_t p_test = static_cast <float> (rand()) /( static_cast <float> (RAND_MAX));
128  fSmearedVar1 = 0;
129  if (p_test <= fResolution){
130  Float_t Ptot = pow(pow(Px,2)+pow(Py,2)+pow(Pz,2),0.5)/1000.;
131  if (Ptot <= 2.0){
132  if (abs(PDGid) == 2212 || abs(PDGid) == 310) fSmearedVar1 = 0;
133  if (abs(PDGid) == 11) fSmearedVar1 = 0;
134  if (abs(PDGid) == 211) fSmearedVar1 = 1;
135 }
136  if (Ptot <= 3.0 && Ptot >2.0){
137  if (abs(PDGid) == 2212 || abs(PDGid) == 310) fSmearedVar1 = 0;
138  if (abs(PDGid) == 11) fSmearedVar1 = 0;
139  if (abs(PDGid) == 211) fSmearedVar1 = 1;
140 }
141  if (Ptot <= 5.0 && Ptot > 3.0){
142  if (abs(PDGid) == 2212 || abs(PDGid) == 310) fSmearedVar1 = 0;
143  if (abs(PDGid) == 11) fSmearedVar1 = 0;
144  if (abs(PDGid) == 211) fSmearedVar1 = 1;
145 }
146  if (Ptot > 5.0){
147  if (abs(PDGid) == 2212 || abs(PDGid) == 310) fSmearedVar1 = 0;
148  if (abs(PDGid) == 11) fSmearedVar1 = 0;
149  if (abs(PDGid) == 211) fSmearedVar1 = 1;
150 }
151 std::cout << "Cherenkov 2: " << fSmearedVar1 << std::endl;
152 }
153 }
154 }
155 
156  /// Vector-based constructor
158  std::vector<Float_t> position,
159  Float_t t,
160  std::vector<Float_t> momentum,
161  Int_t PDGid,
162  Int_t EventID,
163  Int_t TrackID,
164  Float_t Resolution
165  ){
167  fX = position[0]; fY = position[1]; fZ = position[2]; fT = t;
168  fPx = momentum[0]; fPy = momentum[1]; fPz = momentum[2];
169  fPDGid = PDGid;
170  fEventID = EventID;
171  fTrackID = TrackID;
172  fResolution = Resolution;
173 
174  }
175 
176  // TVector3-based constructor
178  TVector3 position,
179  Float_t t,
180  TVector3 momentum,
181  Int_t PDGid,
182  Int_t EventID,
183  Int_t TrackID,
184  Float_t Resolution
185  ){
187  fX = position.X(); fY = position.Y(); fZ = position.Z(); fT = t;
188  fPx = momentum.X(); fPy = momentum.Y(); fPz = momentum.Z();
189  fPDGid = PDGid;
190  fEventID = EventID;
191  fTrackID = TrackID;
192  fResolution = Resolution;
193 
194  }
195 
197  fX = rhs.fX; fY = rhs.fY; fZ = rhs.fZ; fT = rhs.fT;
198  fPx = rhs.fPx; fPy = rhs.fPy; fPz = rhs.fPz;
199  fPDGid = rhs.fPDGid;
200  fEventID = rhs.fEventID;
201  fTrackID = rhs.fTrackID;
205  fResolution = rhs.fResolution;
206  }
207 
208 }
209 
static QCString name
Definition: declinfo.cpp:673
std::string string
Definition: nybbler.cc:12
constexpr T pow(T x)
Definition: pow.h:72
T abs(T value)
Code to link reconstructed objects back to the MC truth information.
list x
Definition: train.py:276
def momentum(x1, x2, x3, scale=1.)
QTextStream & endl(QTextStream &s)
Float_t UnitGauss(Float_t mean, Float_t value, Float_t sigma)