EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
SmearBeAST_0_1.cxx
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file SmearBeAST_0_1.cxx
1 // NOT intended tobe representative of the most recent
2 // BeAST developments (coming soon)
3 // Rather an example, similar to what was used for the design study,
4 // using more complicated parameterizations than the Handbook
5 
6 #include "eicsmear/erhic/VirtualParticle.h"
7 #include "eicsmear/smear/Acceptance.h"
8 #include "eicsmear/smear/Device.h"
9 #include "eicsmear/smear/Detector.h"
10 #include "eicsmear/smear/Smearer.h"
11 #include "eicsmear/smear/ParticleMCS.h"
12 #include "eicsmear/smear/PerfectID.h"
13 #include <eicsmear/smear/Smear.h>
14 #include <eicsmear/erhic/ParticleMC.h>
15 #include "Math/Vector4D.h"
16 
17 // declare static --> local to this file, won't clash with others
18 static double ThetaFromEta( const double eta );
19 
20 Smear::Detector BuildBeAST_0_1() {
21 
22  gSystem->Load("libeicsmear");
23 
24  // Create the detector object to hold all devices
25  Smear::Detector det;
26  det.SetEventKinematicsCalculator("NM JB DA"); // The detector will calculate event kinematics from smeared values
27 
28  // For proper kinematics calculations,
29  // phi and theta need to be smeared (at least with "0") for
30  // all particles that see a kE or kP smearer
31 
32  // Phi and theta smearing for all particles
33  // ----------------------------------------
34  // eta = -3.5 -- 3.5
35  // charged particles in the tracker
36  // includes electrons
37  Smear::Acceptance::Zone TrackZone(ThetaFromEta ( 3.5 ),ThetaFromEta ( -3.5 ));
38  Smear::Device trackTheta(Smear::kTheta, "((1.0/(1.0*P))*(0.752935 + 0.280370*pow((-log(tan(theta/2.0))), 2) - 0.0359713*pow((-log(tan(theta/2.0))), 4) + 0.00200623*pow((-log(tan(theta/2.0))), 6)) + 0.0282315 - 0.00998623*pow((-log(tan(theta/2.0))), 2) + 0.00117487*pow((-log(tan(theta/2.0))), 4) - 0.0000443918*pow((-log(tan(theta/2.0))), 6))*0.001");
39  trackTheta.Accept.AddZone(TrackZone);
40  trackTheta.Accept.SetCharge(Smear::kCharged);
41  det.AddDevice(trackTheta);
42  Smear::Device trackPhi(Smear::kPhi, "((1.0/(1.0*P))*(0.743977 + 0.753393*pow((-log(tan(theta/2.0))), 2) + 0.0634184*pow((-log(tan(theta/2.0))), 4) + 0.0128001*pow((-log(tan(theta/2.0))), 6)) + 0.0308753 + 0.0480770*pow((-log(tan(theta/2.0))), 2) - 0.0129859*pow((-log(tan(theta/2.0))), 4) + 0.00109374*pow((-log(tan(theta/2.0))), 6))*0.001");
43  trackPhi.Accept.AddZone(TrackZone);
44  trackPhi.Accept.SetCharge(Smear::kCharged);
45  det.AddDevice(trackPhi);
46 
47  // emcal stretches to -4.5 < eta < 4.5
48  // includes gammas
49  // Assuming 1 mrad
50  Smear::Acceptance::Zone AngleZoneEmcalB(ThetaFromEta ( -3.5 ),ThetaFromEta ( -4.5 ));
51  Smear::Device SmearThetaEmcalB(Smear::kTheta, "0.001");
52  SmearThetaEmcalB.Accept.AddZone(AngleZoneEmcalB);
53  SmearThetaEmcalB.Accept.SetGenre(Smear::kElectromagnetic);
54  det.AddDevice(SmearThetaEmcalB);
55 
56  Smear::Device SmearPhiEmcalB(Smear::kPhi, "0.001");
57  SmearPhiEmcalB.Accept.AddZone(AngleZoneEmcalB);
58  SmearPhiEmcalB.Accept.SetGenre(Smear::kElectromagnetic);
59  det.AddDevice(SmearPhiEmcalB);
60 
61  Smear::Acceptance::Zone AngleZoneEmcalF(ThetaFromEta ( 4.5 ),ThetaFromEta ( 3.5 ));
62  Smear::Device SmearThetaEmcalF(Smear::kTheta, "0.001");
63  SmearThetaEmcalF.Accept.AddZone(AngleZoneEmcalF);
64  SmearThetaEmcalF.Accept.SetGenre(Smear::kElectromagnetic);
65  det.AddDevice(SmearThetaEmcalF);
66 
67  Smear::Device SmearPhiEmcalF(Smear::kPhi, "0.001");
68  SmearPhiEmcalF.Accept.AddZone(AngleZoneEmcalF);
69  SmearPhiEmcalF.Accept.SetGenre(Smear::kElectromagnetic);
70  det.AddDevice(SmearPhiEmcalF);
71 
72  // Haven't covered gammas in the barrel yet
73  Smear::Acceptance::Zone AngleZoneGamma(ThetaFromEta ( 3.5 ),ThetaFromEta ( -3.5 ));
74  Smear::Device SmearThetaGamma(Smear::kTheta, "0.001");
75  SmearThetaGamma.Accept.AddZone(AngleZoneGamma);
76  SmearThetaGamma.Accept.AddParticle(22);
77  det.AddDevice(SmearThetaGamma);
78 
79  Smear::Device SmearPhiGamma(Smear::kPhi, "0.001");
80  SmearPhiGamma.Accept.AddZone(AngleZoneGamma);
81  SmearPhiGamma.Accept.AddParticle(22);
82  det.AddDevice(SmearPhiGamma);
83 
84  // Wha'ts left is neutral hadrons.
85  // HCal in this example covers +/-4.5
86  // Using the same, very optimistic, 1 mrad
87  Smear::Acceptance::Zone AngleZoneHCal(ThetaFromEta ( 4.5 ),ThetaFromEta ( -4.5 ));
88  Smear::Device SmearThetaHCal(Smear::kTheta, "0.001");
89  SmearThetaHCal.Accept.AddZone(AngleZoneHCal);
90  SmearThetaHCal.Accept.SetGenre(Smear::kHadronic);
91  SmearThetaHCal.Accept.SetCharge(Smear::kNeutral);
92  det.AddDevice(SmearThetaHCal);
93 
94  Smear::Device SmearPhiHCal(Smear::kPhi, "0.001");
95  SmearPhiHCal.Accept.AddZone(AngleZoneHCal);
96  SmearPhiHCal.Accept.SetGenre(Smear::kHadronic);
97  SmearPhiHCal.Accept.SetCharge(Smear::kNeutral);
98  det.AddDevice(SmearPhiHCal);
99 
100  // Tracking
101  // --------
102  // Note: Smear::kCharged checks pdg charge, so includes muons (good)
103 
104  // The momentum parametrization (a*p + b) gives sigma_P/P in percent.
105  // So Multiply through by P and divide by 100 to get absolute sigma_P
106  // Theta and Phi parametrizations give absolute sigma in miliradians
107 
108  // Track Momentum
109  // eta = -3.5 -- 3.5
110  // sigma_p/p ~ 0.1% p+2.0%
111  Smear::Device momentum(Smear::kP, "(P*P*(0.0182031 + 0.00921047*pow((-log(tan(theta/2.0))), 2) - 0.00291243*pow((-log(tan(theta/2.0))), 4) + 0.000264353*pow((-log(tan(theta/2.0))), 6)) + P*(0.209681 + 0.275144*pow((-log(tan(theta/2.0))), 2) - 0.0436536*pow((-log(tan(theta/2.0))), 4) + 0.00367412*pow((-log(tan(theta/2.0))), 6)))*0.01");
112  momentum.Accept.SetCharge(Smear::kCharged);
113  momentum.Accept.AddZone(TrackZone);
114  det.AddDevice(momentum);
115 
116  // EMCal
117  // -----
118  // Calorimeter resolution usually given as sigma_E/E = const% + stocastic%/Sqrt{E}
119  // EIC Smear needs absolute sigma: sigma_E = Sqrt{const*const*E*E + stoc*stoc*E}
120  // eta = -4.5 -- -2
121  Smear::Acceptance::Zone emBck(ThetaFromEta ( -2 ),ThetaFromEta ( -4.5 ));
122  Smear::Device emcalBck(Smear::kE, "sqrt(0.01*0.01*E*E + 0.015*0.015*E)");
123  emcalBck.Accept.AddZone(emBck);
124  emcalBck.Accept.SetGenre(Smear::kElectromagnetic);
125  det.AddDevice(emcalBck);
126 
127  // eta = -2 -- -1
128  Smear::Acceptance::Zone emMidBck(ThetaFromEta ( -1 ),ThetaFromEta ( -2 ));
129  Smear::Device emcalMidBck(Smear::kE, "sqrt(0.01*0.01*E*E + 0.07*0.07*E)");
130  emcalMidBck.Accept.AddZone(emMidBck);
131  emcalMidBck.Accept.SetGenre(Smear::kElectromagnetic);
132  det.AddDevice(emcalMidBck);
133 
134  // eta = -1 -- 1
135  Smear::Acceptance::Zone emMid(ThetaFromEta ( 1 ),ThetaFromEta ( -1 ));
136  Smear::Device emcalMid(Smear::kE, "sqrt(0.01*0.01*E*E + 0.10*0.10*E)");
137  emcalMid.Accept.AddZone(emMid);
138  emcalMid.Accept.SetGenre(Smear::kElectromagnetic);
139  det.AddDevice(emcalMid);
140 
141  // eta = 1 -- 4.5
142  Smear::Acceptance::Zone emFwd(ThetaFromEta ( 4.5 ),ThetaFromEta ( 1 ));
143  Smear::Device emcalFwd(Smear::kE, "sqrt(0.01*0.01*E*E + 0.07*0.07*E)");
144  emcalFwd.Accept.AddZone(emFwd);
145  emcalFwd.Accept.SetGenre(Smear::kElectromagnetic);
146  det.AddDevice(emcalFwd);
147 
148  // HCal
149  // -----
150  // Create the Forward/Backward Hadron Calorimeter
151  // eta = -4.5 -- -1
152  Smear::Acceptance::Zone hBck(ThetaFromEta ( -1 ),ThetaFromEta ( -4.5 ));
153  Smear::Device hcalBck(Smear::kE, "sqrt(0.015*0.015*E*E + 0.50*0.50*E)");
154  hcalBck.Accept.AddZone(hBck);
155  hcalBck.Accept.SetGenre(Smear::kHadronic);
156  det.AddDevice(hcalBck);
157 
158  // eta = 1 -- 4.5
159  Smear::Acceptance::Zone hFwd(ThetaFromEta ( 4.5 ),ThetaFromEta ( 1 ));
160  Smear::Device hcalFwd(Smear::kE, "sqrt(0.015*0.015*E*E + 0.50*0.50*E)");
161  hcalFwd.Accept.AddZone(hFwd);
162  hcalFwd.Accept.SetGenre(Smear::kHadronic);
163  det.AddDevice(hcalFwd);
164 
165  // Create the Hypothetical Mid Rap Calorimeter
166  // eta = -1 -- 1
167  Smear::Acceptance::Zone hMid(ThetaFromEta ( 1 ),ThetaFromEta ( -1 ));
168  //Smear::Device hcalMid(Smear::kE, "sqrt(0.07*0.07*E*E + 0.85*0.85*E)"); // ~CMS
169  Smear::Device hcalMid(Smear::kE, "sqrt(0.02*0.02*E*E + 0.35*0.35*E)"); // ~Zeus
170  hcalMid.Accept.AddZone(hMid);
171  hcalMid.Accept.SetGenre(Smear::kHadronic);
172  det.AddDevice(hcalMid);
173 
174  // // Create PID based on Hermes RICH.
175  // Smear::ParticleID rich("PIDMatrix.dat");
176  // det.AddDevice(rich);
177 
178  return det;
179 }
180 
181 // -------------------------------------------------------------------
182 double ThetaFromEta( const double eta ) {
183  if ( !std::isnan(eta) && !std::isinf(eta) ) {
184  return 2.0 * atan( exp( -eta ));
185  }
186  throw std::runtime_error("ThetaFromEta called with NaN or Inf");
187  return -1;
188 }