Calibration_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 // Class: Calibration
3 // Module Type: analyzer
4 // File: Calibration_module.cc
5 //
6 // Generated at Wed Nov 11 14:47:12 2015 by Andrey Sheshukov using artmod
7 // from cetpkgsupport v1_08_07.
8 ////////////////////////////////////////////////////////////////////////
9 
10 #include <iostream>
20 #include <TH2I.h>
21 #include <TProfile.h>
22 #include <TMath.h>
23 #include "fhiclcpp/ParameterSet.h"
31 
32 
33 // This analyzer's goal is to collect data from events into calibration histograms.
34 // It will collect EventsForCalibration before analizing this data and store for further usage.
35 
36 
37 namespace novaddt {
38  class Calibration;
39 }
40 
41 
43 public:
44  explicit Calibration(fhicl::ParameterSet const & p);
45  // The destructor generated by the compiler is fine for classes
46  // without bare pointers or other resource use.
47 
48  // Plugins should not be copied or assigned.
49  Calibration(Calibration const &) = delete;
50  Calibration(Calibration &&) = delete;
51  Calibration & operator = (Calibration const &) = delete;
52  Calibration & operator = (Calibration &&) = delete;
53 
54  // Required functions.
55  void analyze(art::Event const & e) override;
56 
57  // Selected optional functions.
58  void beginJob() override;
59  void endJob() override;
60 
61 private:
62  void accumulate(art::Event const& e);
63  void recalibrate();
66 
67 private:
68  // Declare member data here.
69  unsigned _EventCounter = 0;
70  unsigned _EventsToCollect = 50; //number of events to accumulate for calibration
71  std::vector<std::string> _hits_id; //[label,instance] of hits
72  std::vector<std::string> _tracks_id; //[label, instance] of tracks
73  TProfile* _histAtt; //attenuation study histogram
74 };
75 
76 
77 
78 double FindW(const novaddt::Track3D& t,const novaddt::DAQHit& h){
79  TVector3 dir=(t.End()-t.Start());
80  dir*=1./dir.Z();
81  TVector3 pos=t.Start()+dir*(h.Plane().val-t.Start().Z());
82 /* std::cout<<"Z: "<<t.Start().Z()<<" - "<<t.End().Z()
83  <<"pos["<<h.Plane().val<<"] = "
84  <<pos.X()<<", "<<pos.Y()<<", "<<pos.Z()
85  <<std::endl;*/
86  if(h.View().val==1)return pos.Y();
87  if(h.View().val==2)return pos.X();
88  return -999;
89 }
90 double FindE(const novaddt::Track3D& t,const novaddt::DAQHit& h){
91  TVector3 dir=(t.End()-t.Start()).Unit();
92  //geometric factor = sin(angle between dir and cell orientation)
93  double Factor=1;
94  if(h.View().val==1)Factor=TMath::Sqrt(1-dir.X()*dir.X());
95  if(h.View().val==2)Factor=TMath::Sqrt(1-dir.Y()*dir.Y());
96  return h.ADC().val*Factor;
97 }
98 
100  EDAnalyzer(p),
101  _hits_id(p.get<std::vector<std::string>>("hits_id")),
102  _tracks_id(p.get<std::vector<std::string>>("tracks_id"))
103 {
104 }
105 
107 {
108  ++_EventCounter;
109  // Hitlist from event (input)
111  e.getByLabel(_hits_id[0],_hits_id[1], hitsHandle);
112 
113  //Create a set of all hits in tracks
114  novaddt::HitSet allTrackHits=fetchAllTrackHits(e);
115  novaddt::HitSet freeHits(*hitsHandle);
116  freeHits.subtract(allTrackHits);
117 
118  //get access to the shared hit map
120  //increment event counter
121  std::cout<<"Counter="<<cms->getCounter()<<std::endl;
122  cms->addCounter();
123  //obtain sharable lock, to indicate we're currently R\W data in map
124  auto lock = cms->GetLock();
125  auto &NoiseMap = cms->GetWMap();
126  //fill the hits map
127  for(const auto &hit: freeHits)
128  NoiseMap[hit]+=1;
129 
130  //now try reading the previous map
131  if(cms->getCounter()==1){
132  MakeNoiseHisto(cms->GetRMap());
133  }
134  //do the attenuation calibration:
135  //read tracks from event
137  e.getByLabel(_tracks_id[0],_tracks_id[1], tracksHandle);
138  //get associated hits
139  art::FindOneP<novaddt::HitList> fohl(tracksHandle, e, _tracks_id[0]);
140  for(size_t nt = 0; nt < tracksHandle->size(); ++nt){
141  auto track=tracksHandle->at(nt);
142  //skip non-3D tracks
143  if(track.Is3D() == false)continue;
144  //skip too vertical tracks
145  if(std::abs(track.End().Z()-track.Start().Z())<30)continue;
146  //track direction cosine
147  //double Cz = (track.End()-track.Start()).CosTheta();
148  //get all hits for this track
149  art::Ptr<novaddt::HitList> trackHits = fohl.at(nt);
150 
151  for(const auto& hit: *trackHits){
152  //get the distance from hit to the center of detector
153  double W=FindW(track,hit);
154  double E=FindE(track,hit);
155  /* std::cout<<"hit (ADC="<<hit.ADC().val
156  <<", cell="<<hit.Cell().val
157  <<", plane="<<hit.Plane().val
158  <<",view="<<hit.View().val<<")"
159  <<" W = "<<W
160  <<" E = "<<E
161  <<std::endl;*/
162  _histAtt->Fill(W,E);
163  }
164 
165  }
166 
167 }
168 
170  //Create a set of all hits, that belong to tracks
171  novaddt::HitSet allTrackHits;
172  //Fill this set from tracks
174  event.getByLabel(_tracks_id[0],_tracks_id[1],trackHitLists);
175  for(const auto& trackHits: *trackHitLists)
176  allTrackHits.add(trackHits);
177  return allTrackHits;
178 }
179 
180 
183  auto hist=tfs->make<TH2I>(Form("h_%d",_EventCounter),
184  Form("Noise, event #%d",_EventCounter),
185  m.Nplanes(),0,m.Nplanes(),
186  m.Ncells(),0,m.Ncells());
187  uint64_t sum=0;
188  for(size_t nc=0;nc<m.Ncells();++nc)
189  for(size_t np=0;np<m.Nplanes();++np){
190  auto val=m(novaddt::Plane(np),novaddt::Cell(nc));
191  sum+=val;
192  hist->SetBinContent(np,nc,val);
193  }
194  std::cout<<"Made hist with "<<sum<<" hits"<<std::endl;
195  return hist;
196 }
197 
199 {
200  std::cout<<"---- Calibration analyzer options -----"<<std::endl;
201  std::cout<<"\tEventsToCollect = "<<_EventsToCollect <<std::endl;
202  std::cout<<"\tInput hits label = ["<<_hits_id[0] <<" : "<<_hits_id[1] <<"]"<<std::endl;
203  std::cout<<"\tInput tracks label = ["<<_tracks_id[0]<<" : "<<_tracks_id[1]<<"]"<<std::endl;
205  _histAtt=tfs->make<TProfile>("hAtt","Attenuation;#cell;ADC_{corr}",
206  200,0,384);
207 }
208 
210 {
211  // Implementation of optional member function here.
212 }
213 
std::vector< std::string > _tracks_id
TH2I * MakeNoiseHisto(const novaddt::calib::FDHitMap &m)
novaddt::Plane const & Plane() const
Definition: DAQHit.h:70
const char * p
Definition: xmltok.h:285
Definition: event.h:19
DEFINE_ART_MODULE(TestTMapFile)
float abs(float number)
Definition: d0nt_math.hpp:39
TVector3 const & End() const
Definition: Track3D.h:46
novaddt::HitSet fetchAllTrackHits(const art::Event &event)
void add(const AnyCollection &v)
Definition: HitSet.h:24
void analyze(art::Event const &e) override
Calibration & operator=(Calibration const &)=delete
static size_t Ncells()
Definition: HitMap.h:35
static size_t Nplanes()
Definition: HitMap.h:34
Float_t E
Definition: plot.C:20
novaddt::ADC const & ADC() const
Definition: DAQHit.h:73
novaddt::View const & View() const
Definition: DAQHit.h:72
TVector3 const & Start() const
Definition: Track3D.h:45
value_type val
Definition: BaseProducts.h:84
EDAnalyzer(Table< Config > const &config)
Definition: EDAnalyzer.h:100
TVector3 Unit() const
OStream cout
Definition: OStream.cxx:6
T * make(ARGS...args) const
TDirectory * dir
Definition: macro.C:5
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
Var Sqrt(const Var &v)
Use to take sqrt of a var.
Definition: Var.cxx:326
Definition: structs.h:12
value_type val
Definition: BaseProducts.h:137
enum BeamMode nc
value_type val
Definition: BaseProducts.h:65
void accumulate(art::Event const &e)
Double_t sum
Definition: plot.C:31
Float_t e
Definition: plot.C:35
TNtuple * nt
Definition: drawXsec.C:2
#define W(x)
Calibration(fhicl::ParameterSet const &p)
double FindW(const novaddt::Track3D &t, const novaddt::DAQHit &h)
void subtract(const AnyCollection &v)
Definition: HitSet.h:30
Definition: fwd.h:28
double FindE(const novaddt::Track3D &t, const novaddt::DAQHit &h)
std::vector< std::string > _hits_id
enum BeamMode string