Project

General

Profile

HitDumper_module.cc

Michelle Stancari, 04/28/2015 03:36 PM

 
1
////////////////////////////////////////////////////////////////////////
2
// Class:       Hitdumper
3
// Module Type: analyzer
4
// File:        Hitdumper_module.cc
5
//
6
// Generated at Sun Mar 24 09:05:02 2013 by Tingjun Yang using artmod
7
// from art v1_02_06.
8
////////////////////////////////////////////////////////////////////////
9
#ifndef Hitdumper_Module
10
#define Hitdumper_Module
11

    
12
// Framework includes
13
#include "art/Framework/Core/EDAnalyzer.h"
14
#include "art/Framework/Core/ModuleMacros.h" 
15
#include "art/Framework/Principal/Event.h" 
16
#include "fhiclcpp/ParameterSet.h" 
17
#include "art/Framework/Principal/Handle.h" 
18
#include "art/Persistency/Common/Ptr.h" 
19
#include "art/Persistency/Common/PtrVector.h" 
20
#include "art/Framework/Services/Registry/ServiceHandle.h" 
21
#include "art/Framework/Services/Optional/TFileService.h" 
22
#include "art/Framework/Services/Optional/TFileDirectory.h" 
23
#include "messagefacility/MessageLogger/MessageLogger.h" 
24

    
25
// LArSoft includes
26
#include "Geometry/Geometry.h"
27
#include "Geometry/PlaneGeo.h"
28
#include "Geometry/WireGeo.h"
29
#include "RecoBase/Hit.h"
30
#include "Utilities/LArProperties.h"
31
#include "Utilities/DetectorProperties.h"
32
#include "Utilities/AssociationUtil.h"
33
#include "Utilities/TimeService.h"
34
#include "Simulation/AuxDetSimChannel.h"
35
#include "RawData/ExternalTrigger.h"
36

    
37
// ROOT includes
38
#include "TTree.h"
39
#include "TTimeStamp.h"
40

    
41
const int kMaxHits       = 10000; //maximum number of hits
42
const int kMaxAuxDets = 100;
43
//const unsigned short kMaxTkIDs = 100;
44
namespace MyHitdumper {
45

    
46
class Hitdumper : public art::EDAnalyzer {
47
public:
48
  explicit Hitdumper(fhicl::ParameterSet const & p);
49
  virtual ~Hitdumper();
50

    
51
   // This method is called once, at the start of the job. In this
52
    // example, it will define the histograms and n-tuples we'll write.
53
    void beginJob();
54

    
55
    // This method is called once, at the start of each run. It's a
56
    // good place to read databases or files that may have
57
    // run-dependent information.
58
  //    void beginRun(const art::Run& run);
59

    
60
    // This method reads in any parameters from the .fcl files. This
61
    // method is called 'reconfigure' because it might be called in the
62
    // middle of a job; e.g., if the user changes parameter values in an
63
    // interactive event display.
64
    void reconfigure(fhicl::ParameterSet const& pset);
65

    
66
    // The analysis routine, called once per event. 
67
    void analyze (const art::Event& evt); 
68

    
69
private:
70

    
71
  void ResetVars();
72
  
73
  TTree* fTree;
74
  //run information
75
  int run;
76
  int subrun;
77
  int event;
78
  double evttime;
79
  double efield[3];
80
  int t0;
81

    
82
  int    nhits;
83
  int    hit_cryostat[kMaxHits];
84
  int    hit_tpc[kMaxHits];
85
  int    hit_plane[kMaxHits];
86
  int    hit_wire[kMaxHits];
87
  int    hit_channel[kMaxHits];
88
  double hit_peakT[kMaxHits];
89
  double hit_charge[kMaxHits];
90
  double hit_ph[kMaxHits];
91

    
92
  int counter_flag[kMaxAuxDets];
93

    
94
  std::string fHitsModuleLabel;
95
  std::string fLArG4ModuleLabel;
96
  std::string fCounterModuleLabel;
97
  double fCombinedTimeDelay;
98

    
99
};
100

    
101

    
102
Hitdumper::Hitdumper(fhicl::ParameterSet const& pset)
103
  : EDAnalyzer(pset)
104
{
105
   // Read in the parameters from the .fcl file.
106
    this->reconfigure(pset);
107
}
108
  
109
void Hitdumper::reconfigure(fhicl::ParameterSet const& p)
110
{   // Read parameters from the .fcl file. The names in the arguments
111
  // to p.get<TYPE> must match names in the .fcl file.
112

    
113
  fHitsModuleLabel = p.get< std::string >("HitsModuleLabel");
114
  fCounterModuleLabel = p.get< std::string >("CounterModuleLabel");
115
  //    fLArG4ModuleLabel(p.get< std::string >("LArGeantModuleLabel", "largeant"));
116
  fLArG4ModuleLabel = "largeant";
117
  fCombinedTimeDelay = p.get< double >("CombinedTimeDelay");
118

    
119
  
120

    
121
  return;
122
}
123
  
124

    
125

    
126

    
127
Hitdumper::~Hitdumper()
128
{
129
  // Clean up dynamic memory and other resources here.
130
}
131

    
132
void Hitdumper::analyze(const art::Event& evt)
133
{
134
  // Implementation of required member function here.
135
  ResetVars();
136

    
137
  art::ServiceHandle<geo::Geometry> geom;
138
  art::ServiceHandle<util::LArProperties> larprop;
139
  art::ServiceHandle<util::DetectorProperties> detprop;
140
  art::ServiceHandle<util::TimeService> timeserv;
141

    
142

    
143
  run = evt.run();
144
  subrun = evt.subRun();
145
  event = evt.id().event();
146

    
147
  /* old "hacker" way to check counter information */
148
  /*
149
  std::vector<const sim::AuxDetSimChannel*> fAuxDetSimChannels;
150
  evt.getView(fLArG4ModuleLabel, fAuxDetSimChannels);
151
  // int nauxdets = fAuxDetSimChannels.size();
152
  // std::cout << "Nauxdets=" << nauxdets << std::endl;
153
  for (size_t i = 0; i<fAuxDetSimChannels.size(); ++i){
154
    const sim::AuxDetSimChannel* c =  fAuxDetSimChannels[i];
155
    int auxdetid = c->AuxDetID();
156
    const std::vector<sim::AuxDetIDE>& setOfIDEs = c->AuxDetIDEs();
157
    if (setOfIDEs.size()>0) {
158
      float totE = 0.0;
159
      for (size_t j = 0; j<setOfIDEs.size(); ++j) {
160
        totE+=setOfIDEs[j].energyDeposited;
161
        float htime = setOfIDEs[j].entryT;
162
        std::cout << "counter " << auxdetid << " time  " << htime 
163
                  << " energy " << setOfIDEs[j].energyDeposited << std::endl;
164
      }
165
      if (totE>0.001) counter_flag[auxdetid]=1;
166
    }
167
  }
168
  */
169

    
170
  // New way to get counter hits
171
  art::Handle< std::vector< raw::ExternalTrigger> > externalTriggerListHandle;
172
  evt.getByLabel(fCounterModuleLabel, externalTriggerListHandle);
173
  std::vector< art::Ptr< raw::ExternalTrigger> > trigs;
174
  art::fill_ptr_vector(trigs,externalTriggerListHandle);
175

    
176
  t0 = detprop->TriggerOffset();  // units of TPC ticks
177
  // int nhits = externalTriggerListHandle->size();
178
  // std::cout << "number of counter hits is " << nhits << std::endl;
179

    
180
  for (auto const& trig : trigs) {
181
    int auxdetid = trig->GetTrigID();
182
    int tick = trig->GetTrigTime();
183
    std::cout << "Counter hit " << auxdetid << " " << tick << std::endl;
184
    // convert Penn Board ticks to TPC ticks
185
    // fCombinedTimeDelay=160 ns is the cable delay offset used in counter simulation
186
    // detprop->SamplingRate();    //time sample in ns
187
    // fTimeService->TriggerOffsetTPC()*1.e3; // ns
188

    
189
    float realtime = (tick+0.5)*32.0 - fCombinedTimeDelay;
190
    float TPCticks = realtime/detprop->SamplingRate()-t0;
191
    // This timing cut only works if there was no timing offset in the 
192
    //   simulation - won't work on "untriggered" or CRY MC samples
193
    // i.e. physics.producers.generator.T0:   [0.] # in ns     
194
    if (TPCticks<=1 && TPCticks>=-1) 
195
      { 
196
        std::cout << "good counter hit " << auxdetid << std::endl;
197
        counter_flag[auxdetid]=1;
198
      }
199
  }
200

    
201
  // run = evt.run();
202
  // subrun = evt.subRun();
203
  // event = evt.id().event();
204

    
205
  art::Timestamp ts = evt.time();
206
  TTimeStamp tts(ts.timeHigh(), ts.timeLow());
207
  evttime = tts.AsDouble();
208

    
209

    
210
  efield[0] = larprop->Efield(0);
211
  efield[1] = larprop->Efield(1);
212
  efield[2] = larprop->Efield(2);
213

    
214
  // moved earlier
215
  //  t0 = detprop->TriggerOffset();
216

    
217
  art::Handle< std::vector<recob::Hit> > hitListHandle;
218
  std::vector<art::Ptr<recob::Hit> > hitlist;
219
  if (evt.getByLabel(fHitsModuleLabel,hitListHandle))
220
    art::fill_ptr_vector(hitlist, hitListHandle);
221
  nhits = hitlist.size();
222

    
223
  for (size_t i = 0; i<hitlist.size(); ++i){
224
    geo::WireID wireid = hitlist[i]->WireID();
225
    hit_cryostat[i]   = wireid.Cryostat;
226
    hit_tpc[i]   = wireid.TPC;
227
    hit_plane[i]   = wireid.Plane;
228
    hit_wire[i]    = wireid.Wire;
229
    hit_channel[i] = hitlist[i]->Channel();
230
    // peak time needs plane dependent offset correction applied.
231
    hit_peakT[i]   = hitlist[i]->PeakTime();
232
    hit_charge[i]  = hitlist[i]->Integral();
233
    hit_ph[i]      = hitlist[i]->PeakAmplitude();
234
  }
235
  fTree->Fill();
236

    
237
}
238

    
239
void Hitdumper::beginJob()
240
{
241
  // Implementation of optional member function here.
242
  art::ServiceHandle<art::TFileService> tfs;
243
  fTree = tfs->make<TTree>("hitdumper","analysis tree");
244
  fTree->Branch("run",&run,"run/I");
245
  fTree->Branch("subrun",&subrun,"subrun/I");
246
  fTree->Branch("event",&event,"event/I");
247
  fTree->Branch("evttime",&evttime,"evttime/D");
248
  fTree->Branch("efield",efield,"efield[3]/D");
249
  fTree->Branch("t0",&t0,"t0/I");
250
  fTree->Branch("nhits",&nhits,"nhits/I");
251
  fTree->Branch("hit_cryostat",hit_cryostat,"hit_cryostat[nhits]/I");
252
  fTree->Branch("hit_tpc",hit_tpc,"hit_tpc[nhits]/I");
253
  fTree->Branch("hit_plane",hit_plane,"hit_plane[nhits]/I");
254
  fTree->Branch("hit_wire",hit_wire,"hit_wire[nhits]/I");
255
  fTree->Branch("hit_channel",hit_channel,"hit_channel[nhits]/I");
256
  fTree->Branch("hit_peakT",hit_peakT,"hit_peakT[nhits]/D");
257
  fTree->Branch("hit_charge",hit_charge,"hit_charge[nhits]/D");
258
  fTree->Branch("hit_ph",hit_ph,"hit_ph[nhits]/D");
259
  fTree->Branch("counter_flag",counter_flag,"counter_flag[100]/I");
260
}
261

    
262
//void Hitdumper::reconfigure(fhicl::ParameterSet const & p)
263
//{
264
//  // Implementation of optional member function here.
265
//}
266
void Hitdumper::ResetVars(){
267

    
268
  run = -99999;
269
  subrun = -99999;
270
  event = -99999;
271
  evttime = -99999;
272
  for (int i = 0; i<3; ++i){
273
    efield[i] = -99999;
274
  }
275
  t0 = -99999;
276
  nhits = -99999;
277
  for (int i = 0; i<kMaxHits; ++i){
278
    hit_cryostat[i] = -99999;
279
    hit_tpc[i] = -99999;
280
    hit_plane[i] = -99999;
281
    hit_wire[i] = -99999;
282
    hit_channel[i] = -99999;
283
    hit_peakT[i] = -99999;
284
    hit_charge[i] = -99999;
285
    hit_ph[i] = -99999;
286
  }
287
  for (int i=0;i<kMaxAuxDets;++i) counter_flag[i]=0;
288
  
289
}
290

    
291
DEFINE_ART_MODULE(Hitdumper)
292

    
293
  }  // namespace Hitdumper
294

    
295
#endif // Hitdumper_Module