ADCShapeFitAna_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 // \file ADCShapeFitAna_module.cc
3 // \brief TODO
4 // \author Christopher Backhouse - bckhouse@caltech.edu
5 ////////////////////////////////////////////////////////////////////////
6 
7 // Framework includes
10 #include "art_root_io/TFileService.h"
13 
15 #include "Geometry/Geometry.h"
16 #include "RawData/RawDigit.h"
17 
18 #include "TCanvas.h"
19 #include "TGraph.h"
20 #include "TH2.h"
21 
22 namespace calib
23 {
24  /// \brief TODO
26  {
27  public:
28  explicit ADCShapeFitAna(fhicl::ParameterSet const& pset);
30 
31  virtual void analyze(art::Event const& evt);
32 
33  void reconfigure(fhicl::ParameterSet const& pset);
34 
35  virtual void beginJob();
36  protected:
38  TH2* fAvg;
39 
40  int fNCanvs;
41  };
42 
43  //......................................................................
45  : EDAnalyzer(pset)
46  , fShapeFit(0)
47  , fAvg(0)
48  , fNCanvs(0)
49  {
50  reconfigure(pset);
51  }
52 
53  //......................................................................
55  {
56  }
57 
58  //......................................................................
59 
61  {
62  }
63 
64  //......................................................................
66  {
67  }
68 
69  //......................................................................
71  {
73 
74  const bool isND = geom->DetId() == novadaq::cnv::kNEARDET;
75 
76  if(!fShapeFit){
78  cet::search_path sp("FW_SEARCH_PATH");
79 
80  if(isND){
81  sp.find_file("Calibrator/adc_shape_fit_table_nd_140_4500.root", fname);
82  }
83  else{
84  sp.find_file("Calibrator/adc_shape_fit_table_fd_380_7000.root", fname);
85  }
86  assert(!fname.empty());
87 
88  // TODO add proper configuration of the mode parameter here, and the table files above
89  fShapeFit = new ADCShapeFitTable(fname, isND, !evt.isRealData(),0);
90  }
91 
92  const double sampleTime = isND ? 125 : 500; // ns
93  double riseTime = (isND ? 140 : 380); // ns
94  const double fallTime = (isND ? 4500 : 7000); // ns
95 
96  // ND data and MC are different, the ND data seems to have a ~140ns rise
97  // time rather than 80ns. This solution isn't very scaleable, but it'll do
98  // for now.
99  if(isND && evt.isRealData()) riseTime = 140;
100 
101  static const double peak =
102  pow(riseTime/(riseTime+fallTime), riseTime/fallTime)-
103  pow(riseTime/(riseTime+fallTime), 1+riseTime/fallTime);
104 
105 
106  if(!fAvg){
108  fAvg = tfs->make<TH2F>("avg", ";TNS", 64*9, -sampleTime, 8*sampleTime, 140, -.2, 1.2);
109 
110  TGraph* nominal = tfs->make<TGraph>();
111 
112  nominal->SetPoint(0, -sampleTime, 0);
113  for(double t = 0; t < 8.1*sampleTime; t += 10){
114  nominal->SetPoint(nominal->GetN(), t, exp(-t/fallTime)*(1-exp(-t/riseTime)));
115  }
116  nominal->SetLineWidth(2);
117  nominal->SetLineColor(kRed);
118 
119  nominal->Write("nominal");
120  }
121 
123  evt.getByLabel("daq", digcol);
124 
125  for(unsigned int digIdx = 0; digIdx < digcol->size(); ++digIdx){
126  const art::Ptr<rawdata::RawDigit> dig(digcol, digIdx);
127 
128  if(dig->NADC() == 1) continue;
129 
130  double adcpeak, base;
131  bool goodTime;
132  const double tns = fShapeFit->TNS(0,
133  dig->ADC(1) - dig->ADC(0),
134  dig->ADC(2) - dig->ADC(0),
135  dig->ADC(3) - dig->ADC(0),
136  adcpeak, goodTime, base);
137 
138  base += dig->ADC(0);
139 
140  if(!goodTime) continue;
141 
142  for(unsigned int i = 0; i < dig->NADC(); ++i){
143  fAvg->Fill(sampleTime*i-tns, (dig->ADC(i)-base)/adcpeak*peak);
144  }
145 
146  if(fNCanvs < 100){
147  ++fNCanvs;
148 
150  TCanvas* c = tfs->make<TCanvas>(TString::Format("trace%d", fNCanvs),
151  TString::Format("trace%d", fNCanvs));
152  c->cd();
153 
154  TH2* axes = new TH2F("", ";TNS;ADC",
155  100, -sampleTime, sampleTime*dig->NADC(),
156  100, dig->ADC(0)-.2*adcpeak, dig->ADC(3)+.2*adcpeak);
157  axes->Draw();
158 
159  TGraph* g = new TGraph;
160  for(unsigned int i = 0; i < dig->NADC(); ++i){
161  g->SetPoint(i, i*sampleTime, dig->ADC(i));
162  }
163 
164  g->SetMarkerStyle(kFullCircle);
165  g->Draw("p same");
166 
167  TGraph* fit = new TGraph;
168 
169  fit->SetPoint(0, -sampleTime, base);
170 
171  for(double t = tns; t < dig->NADC()*sampleTime; t += 10){
172  fit->SetPoint(fit->GetN(), t, base+adcpeak/peak*exp(-(t-tns)/fallTime)*(1-exp(-(t-tns)/riseTime)));
173  }
174 
175  fit->SetLineWidth(2);
176  fit->SetLineColor(kRed);
177  fit->Draw("l same");
178 
179  c->Write();
180  }
181  } // end for digIdx
182  }
183 
185 
186 } // end namespace
187 ////////////////////////////////////////////////////////////////////////
188 
enum BeamMode kRed
constexpr T pow(T x)
Definition: pow.h:72
void reconfigure(fhicl::ParameterSet const &pset)
virtual void analyze(art::Event const &evt)
EDAnalyzer(fhicl::ParameterSet const &pset)
Definition: EDAnalyzer.h:27
DEFINE_ART_MODULE(TestTMapFile)
ADCShapeFitAna(fhicl::ParameterSet const &pset)
std::string find_file(std::string const &filename) const
const XML_Char int const XML_Char int const XML_Char * base
Definition: expat.h:331
bool isRealData() const
CDPStorage service.
bool getByLabel(std::string const &label, std::string const &instance, Handle< PROD > &result) const
Definition: DataViewImpl.h:446
novadaq::cnv::DetId DetId() const
Prefer ds::DetectorService::DetId() instead.
Definition: GeometryBase.h:243
double TNS(double tns0, int16_t adc1, int16_t adc2, int16_t adc3, double &adcpeak, bool &goodTime, double &base) const
int evt
Near Detector in the NuMI cavern.
calib::ADCShapeFitTable * fShapeFit
int16_t ADC(uint32_t i) const
Definition: RawDigit.cxx:58
uint32_t NADC() const
Definition: RawDigit.h:81
void geom(int which=0)
Definition: geom.C:163
assert(nhit_max >=nhit_nbins)
void Format(TGraph *gr, int lcol, int lsty, int lwid, int mcol, int msty, double msiz)
Definition: Style.cxx:154
Definition: fwd.h:29
Encapsulate the geometry of one entire detector (near, far, ndos)
enum BeamMode string