EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
SmearJLEIC_0_1.cxx
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file SmearJLEIC_0_1.cxx
1 // Adapted from https://gitlab.com/eic/escalate/ejana/-/blob/master/src/plugins/reco/eic_smear/SS_DetectorJleic_v1_0_2.cc
2 // Work in progress - no guarantees for correctness!
3 // Characteristics:
4 // -- Calorimetry for CHARGED and PHOTONS: sigma (E) = a *sqrt(E) (+) b*E (+) c
5 // ---- in (Eta() > -3.5 && Eta() < -1.1): a= 0.02, b=0.001, c = 0.005
6 // ---- in (Eta() > -1.1 && Eta() < 3.5): a= 0.1, b=0.01, c = 0.005
7 // -- Calorimetry for NEUTRALS:
8 // 100% / sqrt(p.E()) with NO acceptance cuts!
9 // technically, special case with the same 100% /sqrt(p.E()) for neutrons in Theta() < 0.01
10 // -- Tracking for all CHARGED sigma_PT
11 // ---- in (Eta() > -3.5 && Eta() < 3.5): sigma = 0.01 Pt()^2 (+) 0.005 Pt()
12 // with angular resolution = 0.001 in phi and theta
13 // ---- in Theta() < 0.01 : sigma = 0.02
14 // ---- in 0.01 < Theta()< 0.05 : sigma = 0.05
15 
16 
17 #include "eicsmear/erhic/VirtualParticle.h"
18 #include "eicsmear/smear/Acceptance.h"
19 #include "eicsmear/smear/Device.h"
20 #include "eicsmear/smear/Detector.h"
21 #include "eicsmear/smear/Smearer.h"
22 #include "eicsmear/smear/ParticleMCS.h"
23 #include "eicsmear/smear/PerfectID.h"
24 #include <eicsmear/smear/Smear.h>
25 #include <eicsmear/erhic/ParticleMC.h>
26 #include "Math/Vector4D.h"
27 
28 using std::cout;
29 using std::endl;
30 
31 // declare static --> local to this file, won't clash with others
32 static double ThetaFromEta( const double eta );
33 
34 Smear::Detector BuildJLEIC_0_1() {
35 
36  gSystem->Load("libeicsmear");
37 
38  // Create the detector object to hold all devices
39  Smear::Detector det;
40  det.SetEventKinematicsCalculator("NM JB DA"); // The detector will calculate event kinematics from smeared values
41 
42  // Phi and theta smearing for all particles
43  // ----------------------------------------
44  // Eic-smear rule: if it is smeared, it needs phi and theta information
45  // Here: Finite resolution only given for charged particles in +/-3.5
46  // But neutral hadrons are smeared everywhere, so endow them with perfect resolution
47 
48  // Neutral hadrons - assume perfect resolution
49  Smear::Acceptance::Zone AngleZoneNeutralHadrons (ThetaFromEta(15.),ThetaFromEta(-15.));
50  Smear::Device SmearThetaNeutralHadrons(Smear::kTheta, "0");
51  SmearThetaNeutralHadrons.Accept.AddZone(AngleZoneNeutralHadrons);
52  SmearThetaNeutralHadrons.Accept.SetGenre(Smear::kHadronic);
53  SmearThetaNeutralHadrons.Accept.SetCharge(Smear::kNeutral);
54  det.AddDevice(SmearThetaNeutralHadrons);
55 
56  Smear::Device SmearPhiNeutralHadrons(Smear::kPhi, "0");
57  SmearPhiNeutralHadrons.Accept.AddZone(AngleZoneNeutralHadrons);
58  SmearPhiNeutralHadrons.Accept.SetGenre(Smear::kHadronic);
59  SmearPhiNeutralHadrons.Accept.SetCharge(Smear::kNeutral);
60  det.AddDevice(SmearPhiNeutralHadrons);
61 
62  // Charged particles, barrel+endcap
63  Smear::Acceptance::Zone AngleZoneChargedBarrel(ThetaFromEta ( 3.5 ),ThetaFromEta ( -3.5 ));
64  Smear::Device SmearThetaChargedBarrel(Smear::kTheta, "0.001");
65  SmearThetaChargedBarrel.Accept.AddZone(AngleZoneChargedBarrel);
66  SmearThetaChargedBarrel.Accept.SetCharge(Smear::kCharged);
67  det.AddDevice(SmearThetaChargedBarrel);
68 
69  Smear::Device SmearPhiChargedBarrel(Smear::kPhi, "0.001");
70  SmearPhiChargedBarrel.Accept.AddZone(AngleZoneChargedBarrel);
71  SmearPhiChargedBarrel.Accept.SetCharge(Smear::kCharged);
72  det.AddDevice(SmearPhiChargedBarrel);
73 
74  // barrel + endcap calorimetry also accepts photons
75  Smear::Device SmearThetaCalPhotons(Smear::kTheta, "0.001");
76  SmearThetaCalPhotons.Accept.AddZone(AngleZoneChargedBarrel);
77  SmearThetaCalPhotons.Accept.AddParticle(22);
78  det.AddDevice(SmearThetaCalPhotons);
79 
80  Smear::Device SmearPhiCalPhotons(Smear::kPhi, "0.001");
81  SmearPhiCalPhotons.Accept.AddZone(AngleZoneChargedBarrel);
82  SmearPhiCalPhotons.Accept.AddParticle(22);
83  det.AddDevice(SmearPhiCalPhotons);
84 
85 
86  // Charged particles, far forward - assume perfect resolution
87  // keep theta finite to avoid division by zero
88  Smear::Acceptance::Zone AngleZoneChargedFarForward( 1e-12, 0.05);
89  Smear::Device SmearThetaChargedFarForward(Smear::kTheta, "0");
90  SmearThetaChargedFarForward.Accept.AddZone(AngleZoneChargedFarForward);
91  SmearThetaChargedFarForward.Accept.SetCharge(Smear::kCharged);
92  det.AddDevice(SmearThetaChargedFarForward);
93 
94  Smear::Device SmearPhiChargedFarForward(Smear::kPhi, "0");
95  SmearPhiChargedFarForward.Accept.AddZone(AngleZoneChargedFarForward);
96  SmearPhiChargedFarForward.Accept.SetCharge(Smear::kCharged);
97  det.AddDevice(SmearPhiChargedFarForward);
98 
99  // Calorimetry, Neutral hadrons
100  // ----------------------------
101  // Not separating out the ZDC since it logically does nothing on top of this.
102  Smear::Device CalNeutralHadrons(Smear::kE, "1*sqrt(E)" );
103  CalNeutralHadrons.Accept.AddZone(AngleZoneNeutralHadrons);
104  CalNeutralHadrons.Accept.SetGenre(Smear::kHadronic);
105  CalNeutralHadrons.Accept.SetCharge(Smear::kNeutral);
106  det.AddDevice( CalNeutralHadrons );
107 
108  // Calorimetry, Charged and photons
109  // --------------------------------
110  // There's no elegant way to accept charged hadrons, leptons, and photons
111  // so we make two devices
112  // eta = -3.5 -- -1.1
113  // sigma^2 = (0.02*sqrt(E))^2 + (0.001*E)^2 + 0.005^2;
114  Smear::Acceptance::Zone CalBackZone(ThetaFromEta ( -1.1 ), ThetaFromEta ( -3.5 ));
115  Smear::Device CalBack(Smear::kE, "sqrt( pow(0.02*sqrt(E),2) + pow(0.001*E,2) + pow(0.005,2) )" );
116  CalBack.Accept.AddZone(CalBackZone);
117  CalBack.Accept.SetCharge(Smear::kCharged);
118  CalBack.Accept.SetGenre(Smear::kHadronic);
119  det.AddDevice(CalBack);
120 
121  // eta = -1.1 -- 3.5
122  // sigma^2 = (0.1*sqrt(E))^2 + (0.01*E)^2 + 0.005^2;
123  Smear::Acceptance::Zone CalForwardZone(ThetaFromEta ( 3.5 ), ThetaFromEta ( -1.1 ));
124  Smear::Device CalForward(Smear::kE, "sqrt( pow(0.1*sqrt(E),2) + pow(0.01*E,2) + pow(0.005,2) )" );
125  CalForward.Accept.AddZone(CalForwardZone);
126  CalForward.Accept.SetCharge(Smear::kCharged);
127  CalForward.Accept.SetGenre(Smear::kHadronic);
128  det.AddDevice(CalForward);
129 
130  // eta = -3.5 -- -1.1
131  // sigma^2 = (0.02*sqrt(E))^2 + (0.001*E)^2 + 0.005^2;
132  Smear::Device PhotonCalBack(Smear::kE, "sqrt( pow(0.02*sqrt(E),2) + pow(0.001*E,2) + pow(0.005,2) )" );
133  PhotonCalBack.Accept.AddZone(CalBackZone);
134  PhotonCalBack.Accept.SetGenre(Smear::kElectromagnetic);
135  det.AddDevice(PhotonCalBack);
136 
137  // eta = -1.1 -- 3.5
138  // sigma^2 = (0.1*sqrt(E))^2 + (0.01*E)^2 + 0.005^2;
139  Smear::Device PhotonCalForward(Smear::kE, "sqrt( pow(0.1*sqrt(E),2) + pow(0.01*E,2) + pow(0.005,2) )" );
140  PhotonCalForward.Accept.AddZone(CalForwardZone);
141  PhotonCalForward.Accept.SetGenre(Smear::kElectromagnetic);
142  det.AddDevice(PhotonCalForward);
143 
144  // Tracking
145  // --------
146  // Note: Original implementation has vertex smearing,
147  // which (at least for now) is not supported here
148 
149  // eta = -3.5 -- 3.5
150  Smear::Acceptance::Zone TrackBarrelZone(ThetaFromEta ( 3.5 ),ThetaFromEta ( -3.5 ));
151  Smear::Device TrackBarrel(Smear::kPt, "sqrt ( pow( 0.01 * pow(pT,2), 2) + pow(0.005 * pT,2))");
152  TrackBarrel.Accept.AddZone(TrackBarrelZone);
153  TrackBarrel.Accept.SetCharge(Smear::kCharged);
154  det.AddDevice(TrackBarrel);
155 
156 
157  // Roman Pots, Theta < 0.01
158  Smear::Acceptance::Zone TrackRPZone(1e-7,0.01);
159  Smear::Device TrackRP(Smear::kPt, "0.02");
160  TrackRP.Accept.AddZone(TrackRPZone);
161  TrackRP.Accept.SetCharge(Smear::kCharged);
162  det.AddDevice(TrackRP);
163 
164  // D1, 0.01 < Theta < 0.05
165  Smear::Acceptance::Zone TrackD1Zone(0.01,0.05);
166  Smear::Device TrackD1(Smear::kPt, "0.01");
167  TrackD1.Accept.AddZone(TrackD1Zone);
168  TrackD1.Accept.SetCharge(Smear::kCharged);
169  det.AddDevice(TrackD1);
170 
171  return det;
172 }
173 
174 // -------------------------------------------------------------------
175 double ThetaFromEta( const double eta ) {
176  if ( !std::isnan(eta) && !std::isinf(eta) ) {
177  return 2.0 * atan( exp( -eta ));
178  }
179  throw std::runtime_error("ThetaFromEta called with NaN or Inf");
180  return -1;
181 }
182