LArSoft  v10_04_05
Liquid Argon Software toolkit - https://larsoft.org/
MagDriftAna_module.cc
Go to the documentation of this file.
1 //
3 // MagDriftAna class
4 //
5 // dmckee@phys.ksu.edu
6 //
8 
9 // C++ std library includes
10 #include <algorithm>
11 #include <string>
12 
13 // Framework includes
19 #include "art_root_io/TFileService.h"
21 #include "fhiclcpp/ParameterSet.h"
23 
24 // Root Includes
25 #include "TH2.h"
26 #include "TLine.h"
27 
28 // LArSoft includes
39 
41 namespace hit {
42 
44  class MagDriftAna : public art::EDAnalyzer {
45  public:
46  explicit MagDriftAna(fhicl::ParameterSet const& pset);
47 
48  private:
49  void analyze(const art::Event& evt) override;
50  void endJob() override;
51 
52  // intilize the histograms
53  //
54  // Can't be done in Begin job because I want to use LArProperties
55  // which used the database, so I test and run on each
56  // event. Wasteful and silly, but at least it *works*.
57  void ensureHists(art::Event const& evt, detinfo::DetectorClocksData const& clockData);
58 
60  std::string fLArG4ModuleLabel;
61 
62  // Flag for initialization done, because we set up histograms the
63  // first time through beginRun() so that we can use the
64  // database...
65  bool initDone{false};
66 
67  // Drift properties
68  double fDirCosY{0.};
69  double fDirCosZ{0.};
70 
71  TH1D* fChargeXpos{nullptr}; // << position of the MC Truth charge deposition
72  TH1D* fChargeYpos{nullptr};
73  TH1D* fChargeZpos{nullptr};
74  TH1D* fHitZpos{nullptr}; // << Z position of the recorded hit (from the
75  // z-sensitive wire)
76 
77  TH1D* fDriftDeltaZ{nullptr}; // << Difference in MC charge Z and recorded hit Z
78  TH1D* fDeltaZoverX{nullptr}; // << Delta Z as a function of drift distance
79  TH2D* fDeltaZvsX{nullptr};
80 
81  // Same as above, but only for long drift distances (greater than
82  // 4/5 of the detector)
83  TH1D* fDriftDeltaZAway{nullptr}; // << Difference in MC charge Z and recorded hit Z
84  TH1D* fDeltaZoverXAway{nullptr}; // << Delta Z as a function of drift distance
85 
86  }; // class MagdriftAna
87 
88  //-------------------------------------------------
90  : EDAnalyzer(pset)
91  , fFFTHitFinderModuleLabel{pset.get<std::string>("HitsModuleLabel")}
92  , fLArG4ModuleLabel{pset.get<std::string>("LArGeantModuleLabel")}
93  {}
94 
95  //-------------------------------------------------
97  {
98  if (initDone) return; // Bail if we've already done this.
99  initDone = true; // Insure that we bail later on
100 
102 
103  // Find magnetic field related corrections
104  auto const detProp =
106 
107  // art::ServiceHandle<mag::MagneticField const> MagField;
109  auto const* MagField = MagFieldHandle->provider();
110  double Efield = detProp.Efield();
111  double Temperature = detProp.Temperature();
112  double DriftVelocity = detProp.DriftVelocity(Efield, Temperature) / 1000.;
113 
114  // MagneticField::FieldAtPoint() returns (0, 0, 0) if there is no
115  // field at the requested point, so these direction cosines are
116  // 0 if there is no field
117  fDirCosY = -DriftVelocity * MagField->FieldAtPoint().z() / Efield;
118  fDirCosZ = +DriftVelocity * MagField->FieldAtPoint().y() / Efield;
119  MF_LOG_VERBATIM("MagDriftAna") << "Drift ratios: "
120  << "dY/dX = " << fDirCosY << ", "
121  << "dZ/dX = " << fDirCosZ;
122 
123  // geometry data.
125  // assumes all TPCs are the same
126  auto const& tpc_0 = geom->TPC();
127  double width = 2 * tpc_0.HalfWidth();
128  double halfHeight = tpc_0.HalfHeight();
129  double length = tpc_0.Length();
130 
131  double zScale = std::max(fDirCosZ / 2.0, 4e-4);
132 
133  // Assumes microboone dimensions. Ideally we'd fix this later...
134  fChargeXpos =
135  tfs->make<TH1D>("hChargeXpos", "MC X charge depositions; X (cm); Events", 101, 0.0, width);
136  fChargeYpos = tfs->make<TH1D>(
137  "hChargeYpos", "MC Y charge depositions; Y (cm); Events", 101, -halfHeight, halfHeight);
138  fChargeZpos =
139  tfs->make<TH1D>("hChargeZpos", "MC Z charge depositions; Z (cm); Events", 101, 0.0, length);
140  fHitZpos = tfs->make<TH1D>("hHitZpos", "Z charge collection; Z (cm); Events", 101, 0.0, length);
141 
142  fDriftDeltaZ = tfs->make<TH1D>("hDriftDeltaZ",
143  "Z drift of charge; delta Z (cm); Events",
144  101,
145  -5 * zScale * width,
146  5 * zScale * width);
147  fDeltaZoverX = tfs->make<TH1D>(
148  "hDeltaZoverX", "Z drift of charge; delta Z/X; Events", 51, -10 * zScale, 10 * zScale);
149  fDeltaZvsX = tfs->make<TH2D>("hDeltaZvsX",
150  "delta Z vs X; X (cm); delta Z (cm), Events",
151  51,
152  0.0,
153  width,
154  51,
155  -20 * zScale,
156  20 * zScale);
157 
158  // Some stats only when the Xdrift is large (more than 3/4)
159  fDriftDeltaZAway = tfs->make<TH1D>("hDriftDeltaZAway",
160  "Z drift of charge (long drift); delta Z (cm); Events",
161  101,
162  -5 * zScale * width,
163  5 * zScale * width);
164  fDeltaZoverXAway = tfs->make<TH1D>("hDeltaZoverXAway",
165  "Z drift of charge (long drift); delta Z/X; Events",
166  51,
167  -10 * zScale,
168  10 * zScale);
169  }
170 
171  //-------------------------------------------------
173  {
174  // Add a line on the deltaZ/X graph to denote the calculated value
175  // of the drift ration
176  TLine* l = new TLine(fDirCosZ, 0, fDirCosZ, 1.05 * fDeltaZoverX->GetMaximum());
177  l->SetLineColor(kRed);
178  l->SetLineStyle(kDotted);
179  fDeltaZoverX->GetListOfFunctions()->Add(l);
180 
181  // I know this looks like a memory leak, but each historgram needs
182  // it's own copy of the line to prevent double freeing by the
183  // framework...
184  l = new TLine(fDirCosZ, 0, fDirCosZ, 1.05 * fDeltaZoverX->GetMaximum());
185 
186  l->SetLineColor(kRed);
187  l->SetLineStyle(kDotted);
188  fDeltaZoverXAway->GetListOfFunctions()->Add(l);
189  }
190 
191  //-------------------------------------------------
193  {
194  if (evt.isRealData()) {
195  throw cet::exception("MagDriftAna: ") << "Not for use on Data yet...\n";
196  }
197 
198  auto const clockData = art::ServiceHandle<detinfo::DetectorClocksService const>()->DataFor(evt);
199 
200  ensureHists(evt, clockData);
201 
203  evt.getByLabel(fFFTHitFinderModuleLabel, hitHandle);
204 
205  auto const& wireReadoutGeom = art::ServiceHandle<geo::WireReadout>()->Get();
206 
207  // We're going to want to compare the reconstructed Z with the
208  // simulted Z. For that purpose we use the simultion backtracking.
209 
211 
212  std::vector<art::Ptr<recob::Hit>> hits;
213  art::fill_ptr_vector(hits, hitHandle);
214 
215  geo::WireID hitWireID;
216 
217  //++++++++++
218  // Loop over the hits (data) and fill histos
219  //++++++++++
220  for (auto itr : hits) {
221 
222  hitWireID = itr->WireID();
223  // By assumption the drift occurs only in the z-direction, so
224  // we can get all the info we need from the z-measug plane.
225  if (hitWireID.Plane != (wireReadoutGeom.Nplanes() - 1)) continue;
226 
227  // Charge collected at the wire
228  //
229  // Exactly once for each recob::Hit
230  auto const w0pos = wireReadoutGeom.Plane(hitWireID).Wire(0).GetCenter();
231  double HitZpos =
232  w0pos.Z() +
233  hitWireID.Wire * wireReadoutGeom.FirstPlane({hitWireID.TPC, hitWireID.Plane}).WirePitch();
234  double Charge = itr->Integral();
235  fHitZpos->Fill(HitZpos, Charge);
236 
237  // Charge deposition in the detector
238  std::vector<double> xyz = bt_serv->HitToXYZ(clockData, itr);
239  fChargeXpos->Fill(xyz[0], Charge);
240  fChargeYpos->Fill(xyz[1], Charge);
241  double ChargeZpos = xyz[2];
242  fChargeZpos->Fill(ChargeZpos, Charge);
243 
244  // Delta-Z
245  //
246  // Compares the collected z position from the wire to the
247  // simulated z position
248  double DeltaZ = HitZpos - ChargeZpos;
249  fDriftDeltaZ->Fill(DeltaZ, Charge);
250  // Delta Z correlation with X
251  fDeltaZoverX->Fill(DeltaZ / xyz[0], Charge);
252  fDeltaZvsX->Fill(xyz[0], DeltaZ, Charge);
253  // The X related histograms know the dimensions of the
254  // detector, so we use them to set the "away" limit
255  if (xyz[0] > (fChargeYpos->GetXaxis()->GetXmax() * 0.80)) {
256  fDriftDeltaZAway->Fill(DeltaZ, Charge);
257  fDeltaZoverXAway->Fill(DeltaZ / xyz[0], Charge);
258  }
259 
260  } // loop on Hits
261 
262  return;
263  } // end analyze method
264 
266 
267 } // end of hit namespace
void analyze(const art::Event &evt) override
Declaration of signal hit object.
std::string fLArG4ModuleLabel
std::string fFFTHitFinderModuleLabel
cout<< "Opened file "<< fin<< " ixs= "<< ixs<< endl;if(ixs==0) hhh=(TH1F *) fff-> Get("h1")
Definition: AddMC.C:8
WireID_t Wire
Index of the wire within its plane.
Definition: geo_types.h:430
bool isRealData() const
Definition: Event.cc:53
std::vector< double > HitToXYZ(detinfo::DetectorClocksData const &clockData, const recob::Hit &hit) const
EDAnalyzer(fhicl::ParameterSet const &pset)
Definition: EDAnalyzer.cc:6
Base class for creation of raw signals on wires.
void hits()
Definition: readHits.C:15
virtual const mag::MagneticField * provider() const =0
#define DEFINE_ART_MODULE(klass)
Definition: ModuleMacros.h:65
T get(std::string const &key) const
Definition: ParameterSet.h:314
void endJob() override
PlaneID_t Plane
Index of the plane within its TPC.
Definition: geo_types.h:373
MagDriftAna(fhicl::ParameterSet const &pset)
Detector simulation of raw signals on wires.
Encapsulate the geometry of a wire .
void ensureHists(art::Event const &evt, detinfo::DetectorClocksData const &clockData)
bool getByLabel(std::string const &label, std::string const &instance, Handle< PROD > &result) const
Encapsulate the construction of a single detector plane .
Contains all timing reference information for the detector.
#define MF_LOG_VERBATIM(category)
constexpr WireID()=default
Default constructor: an invalid TPC ID.
TCEvent evt
Definition: DataStructs.cxx:8
void fill_ptr_vector(std::vector< Ptr< T >> &ptrs, H const &h)
Definition: Ptr.h:306
TPCID_t TPC
Index of the TPC within its cryostat.
Definition: geo_types.h:315
TPCGeo const & TPC(TPCID const &tpcid=details::tpc_zero) const
Returns the specified TPC.
Definition: GeometryCore.h:448
Float_t e
Definition: plot.C:35
art framework interface to geometry description
double HalfWidth() const
Width is associated with x coordinate [cm].
Definition: TPCGeo.h:102
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
Encapsulate the construction of a single detector plane .