Skip to content

Commit d074a9d

Browse files
gianipezgianipez
authored
Merge pull request #30 from Mu2e/pgirotti/Gr4Analysis
Adding modules and macros to analyze GR4 calo data
2 parents da9a1a7 + 05e498f commit d074a9d

File tree

4 files changed

+334
-85
lines changed

4 files changed

+334
-85
lines changed

‎otsdaq-mu2e-calorimeter/ArtModules/CaloiercAnalyzer_module.cc‎

Lines changed: 176 additions & 82 deletions
Original file line numberDiff line numberDiff line change
@@ -10,36 +10,35 @@
1010
#include "canvas/Utilities/InputTag.h"
1111

1212
#include "Offline/RecoDataProducts/inc/CaloDigi.hh"
13-
14-
#include <iomanip>
15-
#include <sstream>
16-
#include <vector>
17-
18-
#include <TDirectory.h>
19-
#include <TFile.h>
20-
#include <TGraph.h>
21-
#include <TGraphErrors.h>
22-
#include <TH1F.h>
23-
#include <TH2F.h>
24-
#include <TMath.h>
25-
#include <iostream>
26-
#include "TH1.h"
27-
#include "TH2.h"
2813
#include "art_root_io/TFileService.h"
2914

30-
#include <TCanvas.h>
31-
#include <TF1.h>
32-
#include <TLine.h>
33-
#include <TProfile.h>
34-
#include <TSpline.h>
35-
#include <TTree.h>
36-
#include <TVirtualFitter.h>
3715
#include <algorithm>
3816
#include <cmath>
3917
#include <fstream>
40-
#include <iomanip> // Libreria per manipolare il formato di output
18+
#include <iomanip>
4119
#include <iostream>
4220
#include <map>
21+
#include <sstream>
22+
#include <vector>
23+
24+
#include "TCanvas.h"
25+
#include "TDirectory.h"
26+
#include "TF1.h"
27+
#include "TFile.h"
28+
#include "TFitResult.h"
29+
#include "TGraph.h"
30+
#include "TGraphErrors.h"
31+
#include "TH1.h"
32+
#include "TH1F.h"
33+
#include "TH2.h"
34+
#include "TH2F.h"
35+
#include "TKey.h"
36+
#include "TLine.h"
37+
#include "TMath.h"
38+
#include "TProfile.h"
39+
#include "TSpline.h"
40+
#include "TTree.h"
41+
#include "TVirtualFitter.h"
4342

4443
namespace mu2e {
4544
class CaloiercAnalyzer : public art::EDAnalyzer {
@@ -63,35 +62,45 @@ class CaloiercAnalyzer : public art::EDAnalyzer {
6362
false};
6463
fhicl::Atom<int> skipAfterN{
6564
fhicl::Name("skipAfterN"), fhicl::Comment("Don't fit after N hits"), -1};
65+
fhicl::Atom<bool> produceTree{fhicl::Name("produceTree"),
66+
fhicl::Comment("Produce tree with fit results"),
67+
true};
6668
};
6769

6870
explicit CaloiercAnalyzer(const art::EDAnalyzer::Table<Config>& config);
6971
void analyze(art::Event const& event) override;
7072
void endJob() override;
7173

7274
private:
75+
int load_templates();
76+
TFitResultPtr fit_waveform(TGraph* graph, int sipm_id);
77+
7378
std::string caloDigiModuleLabel_;
7479
std::string caloDigiInstanceLabel_;
7580
int verbosity_;
7681
std::string splineFilename_;
7782
bool uset0_;
7883
int skipAfterN_;
84+
bool produceTree_;
7985

80-
std::map<int, std::vector<int>> digiMap;
81-
std::map<int, std::vector<float>> timeMap;
82-
std::map<int, std::vector<float>> chi2Map;
83-
std::map<int, std::vector<float>> chi2rMap;
86+
std::map<int, std::vector<int>> digiMap;
87+
std::map<int, std::vector<float>> timeMap;
88+
std::map<int, std::vector<float>> chi2Map;
89+
std::map<int, std::vector<float>> chi2rMap;
90+
std::map<int, std::map<int, TSpline3*>> templateMap;
8491

8592
int total_events;
8693
int last_event;
8794
int total_hits;
8895
int total_unsorted;
8996
int total_badfits;
9097

98+
// Fit parameters
99+
const static int FITSCALE = 3850.;
100+
const static int FITOFF = 0.;
101+
const static int FITPED = 0.;
102+
91103
art::ServiceHandle<art::TFileService> tfs;
92-
art::TFileDirectory* sameChanDir;
93-
art::TFileDirectory* sameBoardDir;
94-
art::TFileDirectory* diffBoardDir;
95104

96105
std::map<int, TH1F*> map_h1_dt_singlechan;
97106
std::map<int, TGraph*> map_g_dtevt_singlechan;
@@ -102,8 +111,11 @@ class CaloiercAnalyzer : public art::EDAnalyzer {
102111
std::map<int, TH1F*> map_h1_dt_diffboard;
103112
std::map<int, TGraph*> map_g_dtevt_diffboard;
104113

105-
TSpline3* spline;
106-
TF1* f_spline;
114+
TTree* tree;
115+
long int t_ewt;
116+
std::vector<int>* t_SiPMID = 0;
117+
std::vector<std::vector<float>>* t_times = 0;
118+
std::vector<std::vector<float>>* t_chi2ndf = 0;
107119
};
108120
} // namespace mu2e
109121

@@ -114,44 +126,112 @@ mu2e::CaloiercAnalyzer::CaloiercAnalyzer(const art::EDAnalyzer::Table<Config>& c
114126
, verbosity_(config().verbosity())
115127
, splineFilename_(config().splineFilename())
116128
, uset0_(config().uset0())
117-
, skipAfterN_(config().skipAfterN()) {
118-
TFile* filetemp = TFile::Open(splineFilename_.c_str());
119-
if(!filetemp->IsOpen()) {
120-
std::cout << splineFilename_ << " NOT FOUND" << std::endl;
121-
}
122-
123-
char hname[200];
124-
sprintf(hname, "spline_0_1");
125-
if(filetemp->Get(hname))
126-
spline = (TSpline3*)filetemp->Get(hname);
127-
else {
128-
std::cout << "No template.." << std::endl;
129-
}
130-
131-
f_spline = new TF1(
132-
"f_spline",
133-
[this](double* x, double* par) {
134-
return par[0] * this->spline->Eval(x[0] - par[1]) + par[2];
135-
},
136-
0.,
137-
19.,
138-
3);
139-
f_spline->SetParNames("scale", "tpeak", "ped");
140-
f_spline->SetNpx(10000);
141-
f_spline->SetRange(spline->GetXmin(), spline->GetXmax());
129+
, skipAfterN_(config().skipAfterN())
130+
, produceTree_(config().produceTree()) {
131+
load_templates();
142132

143133
total_events = 0;
144134
total_hits = 0;
145135
total_unsorted = 0;
146136
total_badfits = 0;
147137
// TVirtualFitter::SetDefaultFitter("Minuit");
148138

149-
art::TFileDirectory t_sameChanDir = tfs->mkdir("sameChannel");
150-
art::TFileDirectory t_sameBoardDir = tfs->mkdir("sameBoard");
151-
art::TFileDirectory t_diffBoardDir = tfs->mkdir("diffBoard");
152-
sameChanDir = &t_sameChanDir;
153-
sameBoardDir = &t_sameBoardDir;
154-
diffBoardDir = &t_diffBoardDir;
139+
if(produceTree_) {
140+
tree = tfs->make<TTree>("tree", "Event tree with timing fit results");
141+
tree->Branch("ewt", &t_ewt, "t_ewt/L");
142+
tree->Branch("SiPMID", &t_SiPMID);
143+
tree->Branch("times", &t_times);
144+
tree->Branch("chi2ndf", &t_chi2ndf);
145+
}
146+
}
147+
148+
int mu2e::CaloiercAnalyzer::load_templates() {
149+
TFile template_file(splineFilename_.c_str(), "READ");
150+
if(!template_file.IsOpen()) {
151+
std::cout << splineFilename_ << " NOT FOUND" << std::endl;
152+
return 1;
153+
} else {
154+
std::cout << "LOADING SPLINES FROM " << splineFilename_ << std::endl;
155+
}
156+
157+
TIter nextkey(template_file.GetListOfKeys());
158+
while(auto key = (TKey*)nextkey()) {
159+
std::string templateName = key->GetName();
160+
int board, channel;
161+
if(sscanf(templateName.c_str(), "spline_%d_%d", &board, &channel) == 2) {
162+
// HARDCODED HACK TO MATCH FILE NAMES
163+
if(board == 2) {
164+
board = 1;
165+
} else if(board == 1) {
166+
board = 2;
167+
}
168+
//
169+
170+
templateMap[board][channel] =
171+
(TSpline3*)template_file.Get(templateName.c_str());
172+
std::cout << "Loaded template for Board " << board << " channel " << channel
173+
<< "\n";
174+
}
175+
}
176+
177+
return 0;
178+
}
179+
180+
TFitResultPtr mu2e::CaloiercAnalyzer::fit_waveform(TGraph* graph, int sipm_id) {
181+
int board = (sipm_id / 20) % 6;
182+
int channel = sipm_id % 20;
183+
TSpline3* thisSpline;
184+
if(templateMap.find(board) == templateMap.end()) {
185+
// No templates for this board! Using the first available...
186+
if(!templateMap.empty()) {
187+
auto board_pair = templateMap.begin();
188+
auto channel_pair = board_pair->second.begin();
189+
std::cout << "No available template for Board " << board << " channel "
190+
<< channel << "!\n";
191+
std::cout << "Using the one of Board " << board_pair->first << " channel "
192+
<< channel_pair->first << " instead!\n";
193+
thisSpline = channel_pair->second;
194+
} else {
195+
std::cout << "No templates? Skipping fit...\n";
196+
return -1;
197+
}
198+
} else if(templateMap[board].find(channel) == templateMap[board].end()) {
199+
// No template for this channel -- but board is known
200+
if(!templateMap[board].empty()) {
201+
auto channel_pair = templateMap[board].begin();
202+
std::cout << "No available template for Board " << board << " channel "
203+
<< channel << "!\n";
204+
std::cout << "Using the one of Board " << board << " channel "
205+
<< channel_pair->first << " instead!\n";
206+
thisSpline = channel_pair->second;
207+
} else {
208+
std::cout << "No templates for board " << board << "?? Skipping fit...\n";
209+
return -1;
210+
}
211+
} else { // all good!
212+
thisSpline = templateMap[board][channel];
213+
}
214+
215+
double xmin = thisSpline->GetXmin();
216+
double xmax = thisSpline->GetXmax();
217+
TF1* f_spline = new TF1(
218+
"f_spline",
219+
[&](double* x, double* par) {
220+
return par[0] * thisSpline->Eval(x[0] - par[1]) + par[2];
221+
},
222+
0.,
223+
19.,
224+
3);
225+
f_spline->SetParNames("scale", "tpeak offset", "ped offset", "roc", "chan");
226+
f_spline->SetNpx(10000);
227+
f_spline->SetRange(xmin, xmax);
228+
f_spline->SetParameter(0, FITSCALE);
229+
f_spline->SetParameter(1, FITOFF);
230+
f_spline->SetParameter(2, FITPED);
231+
TFitResultPtr fitResult = graph->Fit(f_spline, "QRNS");
232+
delete f_spline;
233+
234+
return fitResult;
155235
}
156236

157237
void mu2e::CaloiercAnalyzer::analyze(art::Event const& event) {
@@ -163,6 +243,11 @@ void mu2e::CaloiercAnalyzer::analyze(art::Event const& event) {
163243
const auto& caloDigis =
164244
*event.getValidHandle(consumes<mu2e::CaloDigiCollection>(caloDigiModuleLabel_));
165245

246+
if(caloDigis.size() == 0) {
247+
// std::cout << "No calodigis!\n";
248+
return;
249+
}
250+
166251
// Fill time-ordered map for each sipm
167252
digiMap.clear();
168253
for(uint ihit = 0; ihit < caloDigis.size(); ihit++) {
@@ -235,26 +320,21 @@ void mu2e::CaloiercAnalyzer::analyze(art::Event const& event) {
235320
gadc->SetPoint(gi, gi, this_waveform[gi]);
236321
// gadc->SetPointError(gi, 0., 1.);
237322
}
238-
// std::cout<<" ("<<this_waveform.size()<<" points)\n";
239-
f_spline->SetParameters(3850, 0, 0);
240-
// double xmin, xmax;
241-
// f_spline->GetRange(xmin, xmax);
242-
// std::cout<<"Now fitting in range ("<<xmin<<","<<xmax<<")\n";
243-
int fitStatus = gadc->Fit("f_spline", "QRN"); // FIT!
244-
if(fitStatus >= 0) {
245-
// std::cout<<"fitStatus: "<<fitStatus<<", ped:
246-
// "<<f_spline->GetParameter(2)<<", time:
247-
// "<<f_spline->GetParameter(1)<<", scale:
248-
// "<<f_spline->GetParameter(0)<<"\n"; std::cout<<"chi2:
249-
// "<<f_spline->GetChisquare()<<", ndf: "<<f_spline->GetNDF()<<",
250-
// chi2/ndf: "<<f_spline->GetChisquare()/f_spline->GetNDF()<<"\n";
323+
TFitResultPtr fitResult =
324+
fit_waveform(gadc, caloDigis[idx].SiPMID()); // FIT!
325+
if(fitResult >= 0) {
326+
// std::cout<<"fitResult: "<<fitResult<<", ped:
327+
// "<<fitResult->Parameter(2)<<", time:
328+
// "<<fitResult->Parameter(1)<<", scale:
329+
// "<<fitResult->Parameter(0)<<"\n"; std::cout<<"chi2:
330+
// "<<fitResult->Chi2()<<", ndf: "<<fitResult->Ndf()<<",
331+
// chi2/ndf: "<<fitResult->Chi2()/fitResult->Ndf()<<"\n";
251332
timeMap[pair.first].push_back(
252-
5. * (caloDigis[idx].t0() + f_spline->GetParameter(1)));
253-
chi2Map[pair.first].push_back(f_spline->GetChisquare());
254-
chi2rMap[pair.first].push_back(f_spline->GetChisquare() /
255-
f_spline->GetNDF());
333+
5. * (caloDigis[idx].t0() + fitResult->Parameter(1)));
334+
chi2Map[pair.first].push_back(fitResult->Chi2());
335+
chi2rMap[pair.first].push_back(fitResult->Chi2() / fitResult->Ndf());
256336
} else { // if bad fit, don't fill at all
257-
std::cout << "Bad fit status: " << fitStatus << " for sipmid "
337+
std::cout << "Bad fit status: " << fitResult << " for sipmid "
258338
<< pair.first << " hit " << idx << "\n";
259339
total_badfits++;
260340
}
@@ -285,6 +365,20 @@ void mu2e::CaloiercAnalyzer::analyze(art::Event const& event) {
285365
}
286366
}
287367

368+
//------- Fill the tree
369+
if(produceTree_) {
370+
t_ewt = this_eventNumber;
371+
t_SiPMID->clear();
372+
t_times->clear();
373+
t_chi2ndf->clear();
374+
for(auto pair : timeMap) {
375+
t_SiPMID->push_back(pair.first);
376+
t_times->push_back(pair.second);
377+
t_chi2ndf->push_back(chi2rMap[pair.first]);
378+
}
379+
tree->Fill();
380+
}
381+
288382
//--------Now we can fill hists
289383

290384
// Same channel
@@ -443,7 +537,7 @@ void mu2e::CaloiercAnalyzer::analyze(art::Event const& event) {
443537
}
444538

445539
void mu2e::CaloiercAnalyzer::endJob() {
446-
std::cout << "\n-- JOB SUMMARY --\n";
540+
std::cout << "\n-- CaloiercAnalyzer JOB SUMMARY --\n";
447541
std::cout << "Total events: " << total_events << " (last event: " << last_event
448542
<< ")\n";
449543
std::cout << "Total hits: " << total_hits << "\n";

0 commit comments

Comments
 (0)