« back to all changes in this revision
Viewing changes to tools/ParticleIdentification.cc
- browse files at revision 10
- compare with another revision
- download diff
- download tarball
- view history from revision 10
- Committer: Francois Drielsma
- Date: 2017-11-30 19:38:46 UTC
- Revision ID: francois.drielsma@gmail.com-20171130193846-i0mwvz79m5cee3sn
Extensive upgrade
- files added:
- src/AlignStations.cc
- files modified:
- README.txt
added
removed
tools/ParticleIdentification.cc
21
21
_method(""), _tof_ids(0), _tof_params(0), _tof_dz(0.) {
22
22
}
23
23
24
ParticleIdentification::ParticleIdentification(const std::vector<std::string>& data_files,
25
const std::vector<size_t>& ids,
26
const size_t nentries) :
27
_method(""), _tof_ids(ids), _tof_params(0), _tof_dz(0.) {
24
ParticleIdentification::ParticleIdentification(const std::vector<double>& tofs,
25
const std::vector<size_t>& tof_ids,
26
double tof_dz,
27
InfoBox* info) :
28
_method("tof"), _tof_ids(tof_ids), _tof_params(0), _tof_dz(tof_dz),
29
_muon_pdf(), _pion_pdf() {
28
30
29
31
30
32
try {
31
InitializeTOF(data_files, nentries);
33
InitializeTOF(tofs, info);
32
34
} catch ( Exceptions::Exception& e ) {
33
35
throw(Exceptions::Exception(Exceptions::nonRecoverable,
34
36
"Could not initialize the TOF method"+std::string(e.what()),
48
50
_tof_ids = pid._tof_ids;
49
51
_tof_params = pid._tof_params;
50
52
_tof_dz = pid._tof_dz;
53
_muon_pdf = pid._muon_pdf;
54
_pion_pdf = pid._pion_pdf;
51
55
52
56
return *this;
53
57
}
54
58
55
ParticleIdentification::~ParticleIdentification () {}
56
57
void ParticleIdentification::InitializeTOF(std::vector<std::string> data_files, size_t nentries) {
58
59
// Set the method
60
Pitch::print(Pitch::info, "Fitting the time-of-flight TOF"+std::to_string(_tof_ids[0])
61
+"->"+std::to_string(_tof_ids[1])+" for PID ("
62
+std::to_string(nentries)+" requested)");
63
_method = "tof";
64
65
// Counter for the number of particles found so far
66
ProgressBar pbar;
67
size_t n(0);
68
69
// Container for all the time of flights
70
std::vector<double> tof;
71
_tof_dz = 0.;
72
73
for (size_t iFile = 0; iFile < data_files.size(); iFile++) {
74
// Set up the ROOT file and data pointer
75
TFile data_file(data_files[iFile].c_str()); // Load the MAUS output file
76
if ( !data_file.IsOpen() )
77
throw(Exceptions::Exception(Exceptions::nonRecoverable,
78
"Data file not found: "+data_files[iFile],
79
"ParticleIdentification::InitializeTOF"));
80
TTree *T = (TTree*)data_file.Get("Spill"); // Pull out the TTree
81
MAUS::Data *data_ptr = new MAUS::Data(); // A variable to store the Data from each spill
82
T->SetBranchAddress("data", &data_ptr); // Set the address of data_ptr
83
84
// Loop over the spills, get the recon events
85
for (size_t i = 0; i < (size_t)T->GetEntries(); ++i) {
86
87
// Update the spill pointed to by data_ptr
88
T->GetEntry(i);
89
MAUS::Spill* spill = data_ptr->GetSpill();
90
if (spill == NULL || !(spill->GetDaqEventType() == "physics_event"))
91
continue;
92
std::vector<MAUS::ReconEvent*>* revts = spill->GetReconEvents();
93
94
// Loop over recon events in spill
95
for ( size_t ev = 0; ev < revts->size(); ++ev ) {
96
if ( !revts->at(ev) )
97
continue;
98
99
// Etract TOF space points
100
MAUS::TOFEvent *tofevent = revts->at(ev)->GetTOFEvent();
101
MAUS::TOFEventSpacePoint tofsp = tofevent->GetTOFEventSpacePoint();
102
std::vector<std::vector<MAUS::TOFSpacePoint>> tofsps = {tofsp.GetTOF0SpacePointArray(),
103
tofsp.GetTOF1SpacePointArray(),
104
tofsp.GetTOF2SpacePointArray()};
105
106
// Get the times and positions in each TOF, return if not a single SP
107
std::vector<double> t(2);
108
bool found(true);
109
size_t i;
110
for (i = 0 ; i < _tof_ids.size(); i++) {
111
if ( tofsps[_tof_ids[i]].size() == 1 ) {
112
t[i] = tofsps[_tof_ids[i]][0].GetTime();
113
} else {
114
found = false;
115
}
116
}
117
if ( !found )
118
continue;
119
120
// Skip absurd values of the time-of-flight (misreconstruction)
121
if ( t[1] - t[0] < 0 || t[1] - t[0] > 50 )
122
continue;
123
124
// If the distance between TOFs is not yet known, extract it
125
if ( !_tof_dz )
126
_tof_dz = tofsps[_tof_ids[1]][0].GetGlobalPosZ()
127
- tofsps[_tof_ids[0]][0].GetGlobalPosZ();
128
129
// Display the progress in %
130
pbar.GetProgress(n, nentries);
131
132
// Fill the TOF histogram and increment the count
133
tof.push_back(t[1] - t[0]);
134
135
if ( ++n > nentries-1 )
136
break;
137
}
138
if ( n > nentries-1 )
139
break;
140
}
141
data_file.Close();
142
if ( n > nentries-1 )
143
break;
144
}
145
if ( n < nentries )
146
Pitch::print(Pitch::warning, "Requested number of particles not met ("+
147
std::to_string(n)+"/"+std::to_string(nentries)+"), will attempt to fit");
59
ParticleIdentification::~ParticleIdentification () {
60
61
}
62
63
void ParticleIdentification::InitializeTOF(const std::vector<double>& tofs,
64
InfoBox* info) {
148
65
149
66
// Find the time taken by a particle travelling at the speed of light
150
67
double c = 299792458; // [m/s]
151
double ct = 1e9*_tof_dz*1e-3/c; // [ns]
68
double ct = 1e9*_tof_dz/c; // [ns]
152
69
153
70
// Define the limits of the TOF histogram from the recorded data
154
double mean = TMath::Mean(tof.begin(), tof.end());
155
double sig = TMath::RMS(tof.begin(), tof.end());
71
double mean = TMath::Mean(tofs.begin(), tofs.end());
72
double sig = TMath::RMS(tofs.begin(), tofs.end());
156
73
double min = std::max(mean-3*sig, ct+.5); // A bit slower than light or far enough from mus
157
74
double max = mean+3*sig; // Far enough from the mean to encompass everything
158
size_t nbins = (max-min)/.12; // 120 ns bins (2*resolution)
75
size_t nbins = (max-min)/.06; // 60 ps bins (resolution)
159
76
160
// Define a histogram that will contain the time of flight profile
161
// Bin width ~ 250 ps (4*res)
77
// Define a histogram that contains the time of flight profile
162
78
TH1F* htof = new TH1F(TString::Format("tof%d%d", (int)_tof_ids[0], (int)_tof_ids[1]),
163
TString::Format("Time-of-flight;t_{%d%d} [ns];PDF",
79
TString::Format(";t_{%d%d} [ns]",
164
80
(int)_tof_ids[0], (int)_tof_ids[1]), nbins, min, max);
165
htof->FillN(tof.size(), &(tof[0]), NULL);
81
htof->FillN(tofs.size(), &(tofs[0]), NULL);
166
82
htof->SetLineWidth(2);
167
83
84
TCanvas *canv0 = new TCanvas("c", "c", 1200, 800);
85
htof->Draw("E");
86
if ( info )
87
info->Draw();
88
canv0->SaveAs(TString::Format("tof%d%d.pdf", (int)_tof_ids[0], (int)_tof_ids[1]));
89
delete canv0;
90
91
htof->GetYaxis()->SetTitle("PDF [ns^{-1}]");
92
168
93
// Intitialize the TSpectrum object
169
94
const size_t npeaks = 2; // muons, pions
170
95
TSpectrum *spc = new TSpectrum(npeaks);
171
96
172
97
// Identify and substract the background
173
TH1F* background = (TH1F*)spc->Background(htof, 4, "nosmoothing");
98
TH1F* background = (TH1F*)spc->Background(htof, 7, "nosmoothing");
174
99
TH1F* substracted = (TH1F*)htof->Clone();
175
100
substracted->Add(background, -1);
176
101
177
102
// Normalise the histograms to 1 (make it a PDF)
178
103
htof->Sumw2();
179
htof->Scale(1./.12/substracted->GetEntries());
104
htof->Scale(1./htof->GetBinWidth(1)/htof->GetEntries());
180
105
181
106
background->Sumw2();
182
background->Scale(1./.12/substracted->GetEntries());
107
background->Scale(1./htof->GetBinWidth(1)/htof->GetEntries());
183
108
background->SetLineWidth(2);
184
109
185
110
substracted->Sumw2();
186
substracted->Scale(1./.12/substracted->GetEntries());
111
substracted->Scale(1./htof->GetBinWidth(1)/htof->GetEntries());
187
112
substracted->SetLineColor(kGreen+1);
188
113
substracted->SetLineWidth(2);
189
114
190
115
THStack* stack = new THStack("stack", TString::Format("%s;%s;%s",
191
116
htof->GetTitle(), htof->GetXaxis()->GetTitle(), htof->GetYaxis()->GetTitle()));
192
stack->Add(htof);
193
stack->Add(background);
194
stack->Add(substracted);
117
TLegend* leg = new TLegend(.7, .7, .89, .89);
118
leg->SetLineColorAlpha(0, 0);
119
leg->SetFillColor(0);
195
120
196
121
// Search the cleaned up sample
197
122
size_t nfound = spc->Search(substracted);
206
131
207
132
// Define a profile function to fit to the substracted distribution
208
133
TF1 *fpeaks = SimpleTofProfile(posx[0], posx[1], min, max);
209
// TF1 *fpeaks = TofProfile(posx[0], posx[1], min, max);
210
134
fpeaks->SetLineColor(substracted->GetLineColor());
211
135
212
136
// Fit peaks with the profile function and extract the relevant parameters
213
137
substracted->Fit("fpeaks", "RQ");
214
std::vector<double> pars;
138
std::vector<double> pars, errs;
215
139
for(size_t p = 0; p < npeaks; p++) {
140
pars.push_back(fpeaks->GetParameter(3*p));
216
141
pars.push_back(fpeaks->GetParameter(3*p+1));
217
142
pars.push_back(fpeaks->GetParameter(3*p+2));
143
errs.push_back(fpeaks->GetParError(3*p));
144
errs.push_back(fpeaks->GetParError(3*p+1));
145
errs.push_back(fpeaks->GetParError(3*p+2));
218
146
}
147
double lim = (pars[1]*pars[5]+pars[4]*pars[2])/(pars[2]+pars[5]);
148
149
// Build a muon and a pion PDF
150
_pion_pdf = TH1F("pion_pdf", "Pion PDF", nbins, min, max);
151
_pion_pdf.SetLineWidth(2);
152
_pion_pdf.SetLineColor(kRed+2);
153
_pion_pdf.SetLineStyle(2);
154
for (size_t i = 0; i < nbins; i++)
155
if ( substracted->GetBinCenter(i) > lim-.06 ) {
156
_pion_pdf.SetBinContent(i+1, substracted->GetBinContent(i+1));
157
_pion_pdf.SetBinError(i+1, substracted->GetBinError(i+1));
158
}
159
_pion_pdf.Sumw2();
160
161
_muon_pdf = TH1F(*(TH1F*)htof->Clone());
162
_muon_pdf.SetName("muon_pdf");
163
_muon_pdf.SetTitle("Muon PDF");
164
_muon_pdf.Add(&_pion_pdf, -1);
165
_muon_pdf.SetLineColor(kGreen+2);
166
167
stack->Add(&_muon_pdf);
168
stack->Add(&_pion_pdf);
169
170
leg->AddEntry(&_muon_pdf, _muon_pdf.GetTitle(), "lf");
171
leg->AddEntry(&_muon_pdf, TString::Format("MPV: %0.2f ns", pars[1]), "");
172
leg->AddEntry(&_pion_pdf, _pion_pdf.GetTitle(), "lf");
173
leg->AddEntry(&_pion_pdf, TString::Format("MPV: %0.2f ns", pars[4]), "");
219
174
220
175
// The first PID parameter is the limit between the positron and muon peaks
221
_tof_params.push_back(std::max(pars[0]-5*pars[1], ct+.5));
176
_tof_params.push_back(std::max(pars[1]-5*pars[2], ct+.5));
222
177
223
178
// The second PID parameter is the limit between the muon and pion peaks
224
_tof_params.push_back((pars[0]*pars[3]+pars[2]*pars[1])/(pars[1]+pars[3]));
179
_tof_params.push_back((pars[1]*pars[5]+pars[4]*pars[2])/(pars[2]+pars[5]));
225
180
}
226
181
227
182
// Plot the histogram
229
184
gStyle->SetStatY(.9);
230
185
gStyle->SetOptFit(1);
231
186
TCanvas *canv = new TCanvas("c", "c", 1200, 800);
232
stack->Draw("NOSTACK");
233
canv->SaveAs(TString::Format("tof%d%d.pdf", (int)_tof_ids[0], (int)_tof_ids[1]));
187
stack->Draw("HIST LF");
188
leg->Draw("SAME");
189
if ( info )
190
info->Draw();
191
canv->SaveAs(TString::Format("tof%d%d_pdf.pdf", (int)_tof_ids[0], (int)_tof_ids[1]));
234
192
delete canv;
235
193
delete htof;
236
194
241
199
"ParticleIdentification::InitializeTOF"));
242
200
}
243
201
244
int ParticleIdentification::GetID(MAUS::ReconEvent* event) const {
245
246
if ( _method == "tof" ) {
247
return GetTofID(event);
248
}
249
250
throw(Exceptions::Exception(Exceptions::recoverable,
251
"Particle identification method not initialized: "+_method,
252
"ParticleIdentification::GetID"));
253
return 0;
254
}
255
256
int ParticleIdentification::GetTofID(MAUS::ReconEvent* event) const {
257
258
// Extract TOF space points
259
MAUS::TOFEvent *tofevent = event->GetTOFEvent();
260
MAUS::TOFEventSpacePoint tofsp = tofevent->GetTOFEventSpacePoint();
261
std::vector<std::vector<MAUS::TOFSpacePoint>> tofsps = {tofsp.GetTOF0SpacePointArray(),
262
tofsp.GetTOF1SpacePointArray(),
263
tofsp.GetTOF2SpacePointArray()};
264
265
// Get the times and positions in each TOF, return if not a single SP
266
std::vector<double> t(2);
267
size_t i;
268
for (i = 0 ; i < _tof_ids.size(); i++) {
269
if ( tofsps[_tof_ids[i]].size() == 1 ) {
270
t[i] = tofsps[_tof_ids[i]][0].GetTime();
271
} else {
272
return 0;
273
}
274
}
275
276
return GetTofID(t[_tof_ids[1]]-t[_tof_ids[0]]);
277
}
278
279
int ParticleIdentification::GetTofID(double tof) const {
280
281
// Electron criterion, either 5 sigma below the muon peak or close to the speed of light
202
int ParticleIdentification::GetTofID(double tof) {
203
204
// If the muon pdf is 0, it is a positron
282
205
if ( tof < _tof_params[0] )
283
206
return 11;
284
207
285
// Perform pid based on the time-of-flight (muon-pion separation)
286
if ( tof < _tof_params[1] )
208
// If the pion pdf is 0, it is a muon
209
double p_pi = _pion_pdf.Interpolate(tof);
210
if ( !p_pi )
287
211
return 13;
288
212
289
return 211;
213
// Perform pid based on the log-likelihood otherwise
214
double p_mu = _muon_pdf.Interpolate(tof);
215
double llmu = log(p_mu)-log(p_pi);
216
double pmu = exp(llmu)/(exp(llmu)+exp(-llmu));
217
if ( pmu < .1 )
218
return 211;
219
220
return 13;
221
}
222
223
double ParticleIdentification::GetMass(int pid) const {
224
225
if ( fabs(pid) == 11 )
226
return 0.511;
227
228
if ( fabs(pid) == 13 )
229
return 105.66;
230
231
if ( fabs(pid) == 211 )
232
return 139.57;
233
234
throw(Exceptions::Exception(Exceptions::nonRecoverable,
235
"Geant4 particle ID not recognized: "+std::to_string(pid),
236
"ParticleIdentification::GetMass"));
290
237
}
291
238
292
239
TF1* ParticleIdentification::SimpleTofProfile(double mumu, double mupi, double min, double max) const {
293
240
294
241
// Define a two peaks gaussian
295
242
TF1 *fpeaks = new TF1("fpeaks",
296
"[0]*TMath::Gaus(x, [1], [2])/sqrt(2*TMath::Pi()*[2]*[2])+"
297
"[3]*TMath::Gaus(x, [4], [5])/sqrt(2*TMath::Pi()*[5]*[5])",
243
"[0]*TMath::Gaus(x, [1], [2], 1)+[3]*TMath::Gaus(x, [4], [5], 1)",
298
244
min, max);
299
245
fpeaks->SetParNames("A_{#mu}", "#mu_{#mu}", "#sigma_{#mu}",
300
246
"A_{#pi}", "#mu_{#pi}", "#sigma_{#pi}");
Loggerhead is a web-based interface for Breezy
Version: 2.0.1