« back to all changes in this revision
Viewing changes to align_trackers.cpp
- browse files at revision 1
- compare with another revision
- compare with revision 1
- download diff
- download tarball
- view history from revision 1
- Committer: francois.drielsma at gmail
- Date: 2016-07-09 16:52:28 UTC
- Revision ID: francois.drielsma@gmail.com-20160709165228-box2rjzdtb7katwo
Alignment code 0.1
- files added:
- build
added
removed
align_trackers.cpp
1
// Cpp includes
2
#include <ctime>
3
#include <cmath>
4
#include <iostream>
5
#include <fstream>
6
#include <algorithm>
7
#include <regex>
8
9
// Root includes
10
#include "TROOT.h"
11
#include "TSystem.h"
12
#include "TStyle.h"
13
#include "TObject.h"
14
#include "TString.h"
15
#include "TFile.h"
16
#include "TTree.h"
17
#include "TCanvas.h"
18
#include "TLegend.h"
19
#include "TLatex.h"
20
#include "TH1D.h"
21
#include "TH2D.h"
22
#include "TF1.h"
23
#include "TF2.h"
24
#include "TVector3.h"
25
#include "TMath.h"
26
#include "TSpectrum.h"
27
28
// MAUS includes
29
#include "DataStructure/Spill.hh"
30
#include "DataStructure/Data.hh"
31
#include "DataStructure/ReconEvent.hh"
32
#include "DataStructure/TOFEvent.hh"
33
#include "DataStructure/SciFiEvent.hh"
34
35
// Returns the global coordinates of the detectors (fetched in the geometry)
36
class DetectorGlobals {
37
public:
38
39
// Default constructor, initilizes everything to 0
40
DetectorGlobals() : _globals(0., 0., 0.), _angles (0., 0., 0.) {}
41
42
// Proper constructor
43
DetectorGlobals(std::string filename, std::string name) {
44
this->Load(filename, name);
45
}
46
47
// Copy constructor
48
DetectorGlobals(const DetectorGlobals& detector_globals) {
49
*this = detector_globals;
50
}
51
52
// Equality operator
53
DetectorGlobals operator=(const DetectorGlobals& detector_globals) {
54
if ( this == &detector_globals )
55
return *this;
56
57
this->SetDetectorGlobals(detector_globals._globals);
58
this->SetDetectorRotations(detector_globals._angles);
59
60
return *this;
61
}
62
63
// Destructor
64
~DetectorGlobals () {}
65
66
// Loads the global constants from the parent geometry file
67
bool Load(std::string filename, std::string name) {
68
std::ifstream infile(filename.c_str());
69
if ( !infile.is_open() ) {
70
std::cerr << "The geometry file could not be found: " << filename << std::endl;
71
return false;
72
}
73
74
bool det_found(false);
75
while ( !infile.eof() ) {
76
std::string word;
77
infile >> word;
78
if ( word.find(name) != std::string::npos ) {
79
det_found = true;
80
while ( word != "}" ) {
81
infile >> word;
82
if ( word == "Position" ) {
83
infile >> word;
84
_globals.SetX(std::atof(word.c_str()));
85
infile >> word;
86
_globals.SetY(std::atof(word.c_str()));
87
infile >> word;
88
_globals.SetZ(std::atof(word.c_str()));
89
} else if ( word == "Rotation" ) {
90
infile >> word;
91
_angles.SetX(std::atof(word.c_str()));
92
infile >> word;
93
_angles.SetY(std::atof(word.c_str()));
94
infile >> word;
95
_angles.SetZ(std::atof(word.c_str()));
96
}
97
}
98
99
break;
100
}
101
}
102
103
infile.close();
104
if ( !det_found ) {
105
std::cerr << "Detector not found in geometry file: " << name << std::endl;
106
return false;
107
}
108
109
return true;
110
}
111
112
// Setters and getters
113
void SetDetectorGlobals(TVector3 globals) { _globals = globals; }
114
TVector3 GetDetectorGlobals() { return _globals; }
115
void SetDetectorRotations(TVector3 angles) { _angles = angles; }
116
TVector3 GetDetectorRotations() { return _angles; }
117
double x() { return _globals.x(); }
118
double y() { return _globals.y(); }
119
double z() { return _globals.z(); }
120
double alpha() { return _angles.x(); }
121
double beta() { return _angles.y(); }
122
double gamma() { return _angles.z(); }
123
124
private:
125
TVector3 _globals;
126
TVector3 _angles;
127
};
128
129
// Computes the geometrical distance between two detector space points
130
double distance (const TVector3& a, const TVector3& b) {
131
132
return sqrt( pow(b.x()-a.x(), 2) +
133
pow(b.y()-a.y(), 2) +
134
pow(b.z()-a.z(), 2) );
135
}
136
137
138
139
140
// Main alignment code
141
int main(int argc, char ** argv) {
142
143
// Return if there is no input file
144
if ( argc < 2 ) {
145
std::cerr << "No data card specified, please execute in the following fashion:" << std::endl;
146
std::cerr << "./align_trackers datacard.txt data0.root [... dataN.root]" << std::endl;
147
return 2;
148
}
149
150
if ( argc == 2 ) {
151
std::cerr << "No data file specified, please execute in the following fashion:" << std::endl;
152
std::cerr << "./align_trackers datacard.txt data0.root [... dataN.root]" << std::endl;
153
return 2;
154
}
155
156
// Get the run name from the file name (remove superfluous with regex)
157
std::string run_name = std::regex_replace(std::string(argv[2]),
158
std::regex("(.*/)|(.root)"), "");
159
160
// Fitting algorithm parameters
161
std::string geo_file;
162
size_t iterations(0);
163
bool do_pid(false);
164
bool do_roll(false);
165
166
std::ifstream datacard(argv[1]);
167
if ( !datacard.is_open() ) {
168
std::cerr << "Datacards not found, please verify path: " << argv[1] << std::endl;
169
return 2;
170
}
171
while ( !datacard.eof() ) {
172
std::string variable, value;
173
datacard >> variable >> value;
174
if ( variable == "geometry_filename" )
175
geo_file = value;
176
else if ( variable == "number_of_iterations" )
177
iterations = atoi(value.c_str());
178
else if ( variable == "do_pid" )
179
do_pid = (bool)atoi(value.c_str());
180
else if ( variable == "do_roll" )
181
do_roll = (bool)atoi(value.c_str());
182
else if ( variable.size() )
183
std::cerr << "Unrecognised argument: " << variable << std::endl;
184
}
185
186
187
// Position of the detectors of interest in the beam line
188
DetectorGlobals tof1_globals(geo_file, "TOF1");
189
DetectorGlobals tof2_globals(geo_file, "TOF2");
190
DetectorGlobals tku_globals(geo_file, "Tracker0");
191
DetectorGlobals tkd_globals(geo_file, "Tracker1");
192
193
// Trackers alignment constants
194
double alphau(0.), betau(0.), gammau(0.), xu(0.), yu(0.);
195
double alphad(0.), betad(0.), gammad(0.), xd(0.), yd(0.);
196
197
// METHOD THAT DOES NOT INCLUDE THE ROLL OF THE TRACKERS
198
// Residuals distributions to fit at the centres (not roll sensitive)
199
TH2D* htku_resxy = new TH2D("tku_resxy",
200
"Residuals TKU/TOF12;"
201
"x_{u}-#xi_{12}^{u} [mm];y_{u}-#upsilon_{12}^{u} [mm]",
202
50, -250, 250, 50, -250, 250);
203
TH2D* htku_resxpyp = new TH2D("tku_resxpyp",
204
"Angular residuals TKU/TOF12;"
205
"x'_{u}-#xi'_{12} [mrad];y'_{u}-#upsilon'_{12} [mrad]",
206
50, -100, 100, 50, -100, 100);
207
208
TH2D* htkd_resxy = new TH2D("tkd_resxy",
209
"Residuals TKD/TOF12;"
210
"x_{d}-#xi_{12}^{d} [mm];y_{d}-#upsilon_{12}^{d} [mm]",
211
50, -250, 250, 50, -250, 250);
212
TH2D* htkd_resxpyp = new TH2D("tkd_resxpyp",
213
"Angular residuals TKD/TOF12;"
214
"x'_{d}-#xi'_{12} [mrad];y'_{d}-#upsilon'_{12} [mrad]",
215
50, -100, 100, 50, -100, 100);
216
217
std::vector<TH2D*> res_hists;
218
res_hists.push_back(htku_resxy);
219
res_hists.push_back(htku_resxpyp);
220
res_hists.push_back(htkd_resxy);
221
res_hists.push_back(htkd_resxpyp);
222
223
// METHOD THAT INCLUDES THE ROLL OF THE TRACKERS
224
// Residuals distributions as a function the position (roll sensitive)
225
TH2D* htku_resy_x = new TH2D("tku_resy_x",
226
"res_{y} vs x in TKU; x_{u} [mm]; y_{u}-#upsilon_{12}^{u} [mm]",
227
40, -200, 200, 50, -250, 250);
228
TH2D* htku_resx_y = new TH2D("tku_resx_y",
229
"res_{x} vs y in TKU; y_{u} [mm]; x_{u}-#xi_{12}^{u} [mm]",
230
40, -200, 200, 50, -250, 250);
231
TH2D* htku_resyp_xp = new TH2D("tku_resyp_xp",
232
"res_{y'} vs x' in TKU; x'_{u} [mrad]; y'_{u}-#upsilon'_{12} [mrad]",
233
25, -50, 50, 50, -100, 100);
234
TH2D* htku_resxp_yp = new TH2D("tku_resxp_yp",
235
"res_{x'} vs y' in TKU; y'_{u} [mrad]; x'_{u}-#xi'_{12} [mrad]",
236
25, -50, 50, 50, -100, 100);
237
238
TH2D* htkd_resy_x = new TH2D("tkd_resy_x",
239
"res_{y} vs x in TKD; x_{d} [mm]; y_{d}-#upsilon_{12}^{d} [mm]",
240
40, -200, 200, 50, -250, 250);
241
TH2D* htkd_resx_y = new TH2D("tkd_resx_y",
242
"res_{x} vs y in TKD; y_{d} [mm]; x_{d}-#xi_{12}^{d} [mm]",
243
40, -200, 200, 50, -250, 250);
244
TH2D* htkd_resyp_xp = new TH2D("tkd_resyp_xp",
245
"res_{y'} vs x' in TKD; x'_{d} [mrad]; y'_{d}-#upsilon'_{12} [mrad]",
246
25, -50, 50, 50, -100, 100);
247
TH2D* htkd_resxp_yp = new TH2D("tkd_resxp_yp",
248
"res_{x'} vs y' in TKD; y'_{d} [mrad]; x'_{d}-#xi'_{12} [mrad]",
249
25, -50, 50, 50, -100, 100);
250
251
std::vector<TH2D*> res_hists_roll;
252
res_hists_roll.push_back(htku_resx_y);
253
res_hists_roll.push_back(htku_resy_x);
254
res_hists_roll.push_back(htku_resxp_yp);
255
res_hists_roll.push_back(htku_resyp_xp);
256
res_hists_roll.push_back(htkd_resx_y);
257
res_hists_roll.push_back(htkd_resy_x);
258
res_hists_roll.push_back(htkd_resxp_yp);
259
res_hists_roll.push_back(htkd_resyp_xp);
260
261
// Iterative histograms
262
TH1D* htku_optx = new TH1D("tku_optx", "TKU optimal x;;x_{TKU}^{*} [mm]",
263
iterations, .5, iterations+.5);
264
TH1D* htku_opty = new TH1D("tku_opty", "TKU optimal y;;y_{TKU}^{*} [mm]",
265
iterations, .5, iterations+.5);
266
TH1D* htku_optalpha = new TH1D("tku_optalpha", "TKU optimal #alpha;;#alpha_{TKU}^{*} [mrad]",
267
iterations, .5, iterations+.5);
268
TH1D* htku_optbeta = new TH1D("tku_optbeta", "TKU optimal #beta;;#beta_{TKU}^{*} [mrad]",
269
iterations, .5, iterations+.5);
270
TH1D* htku_optgamma = new TH1D("tku_optgamma", "TKU optimal #gamma;;gamma_{TKU}^{*} [mrad]",
271
iterations, .5, iterations+.5);
272
273
TH1D* htkd_optx = new TH1D("tkd_optx", "TKD optimal x;;x_{TKD}^{*} [mm]",
274
iterations, .5, iterations+.5);
275
TH1D* htkd_opty = new TH1D("tkd_opty", "TKD optimal y;;y_{TKD}^{*} [mm]",
276
iterations, .5, iterations+.5);
277
TH1D* htkd_optalpha = new TH1D("tkd_optalpha", "TKD optimal #alpha;;#alpha_{TKD}^{*} [mrad]",
278
iterations, .5, iterations+.5);
279
TH1D* htkd_optbeta = new TH1D("tkd_optbeta", "TKD optimal #beta;;#beta_{TKD}^{*} [mrad]",
280
iterations, .5, iterations+.5);
281
TH1D* htkd_optgamma = new TH1D("tkd_optgamma", "TKD optimal #gamma;;#gamma_{TKD}^{*} [mrad]",
282
iterations, .5, iterations+.5);
283
284
std::vector<TH1D*> ite_hists;
285
ite_hists.push_back(htku_optx);
286
ite_hists.push_back(htku_opty);
287
ite_hists.push_back(htku_optalpha);
288
ite_hists.push_back(htku_optbeta);
289
ite_hists.push_back(htku_optgamma);
290
ite_hists.push_back(htkd_optx);
291
ite_hists.push_back(htkd_opty);
292
ite_hists.push_back(htkd_optalpha);
293
ite_hists.push_back(htkd_optbeta);
294
ite_hists.push_back(htkd_optgamma);
295
296
for (size_t i = 0; i < ite_hists.size(); i++)
297
ite_hists[i]->GetXaxis()->SetTitle("Iteration ID");
298
299
// PARTICLE IDENTIFICATION
300
double mu[2] = {0., 0.}; // Peak means
301
double sig[2] = {0., 0.}; // Peak spreads
302
303
// Here the muon and pion peaks are identified and fitted (for pion beams only)
304
if ( do_pid ) {
305
// Number of particles to look for before trying to fit the distribution
306
int n(0);
307
const int nentries(10000);
308
309
std::cerr << "Fitting the time-of-flight TOF1->2 for PID (";
310
std::cerr << nentries << " entries requested)" << std::endl;
311
312
// Find the time taken by a particle travelling at the speed of light
313
double c = 299792458; // [m/s]
314
double s12 = tof2_globals.z()-tof1_globals.z(); // [mm]
315
double ct12 = 1e9*s12*1e-3/c; // [ns]
316
317
// Define a histogram that will contain the time of flight profile
318
// Bin width ~ 250 ps (5*res)
319
// Min TOF = ct12+.5 ns (a bit slower than light)
320
// Max TOF = 50 ns (70 MeV/c muons)
321
TH1F* htof12 = new TH1F("tof12", "Time-of-flight", 100, ct12+0.5, 50);
322
htof12->GetXaxis()->SetTitle("TOF_{12} [ns]");
323
324
for (size_t iFile = 2; iFile < (size_t)argc; iFile++) {
325
// Set up the ROOT file and data pointer
326
TFile f1(argv[iFile]); // Load the MAUS output file
327
TTree *T = (TTree*)f1.Get("Spill"); // Pull out the TTree
328
MAUS::Data *data_ptr = new MAUS::Data(); // A variable to store the Data from each spill
329
T->SetBranchAddress("data", &data_ptr); // Set the address of data_ptr
330
331
// Loop over the spills, get the recon events
332
for (size_t i = 0; i < (size_t)T->GetEntries(); ++i) {
333
334
// Update the spill pointed to by data_ptr
335
T->GetEntry(i);
336
MAUS::Spill* spill = data_ptr->GetSpill();
337
if (spill == NULL || !(spill->GetDaqEventType() == "physics_event"))
338
continue;
339
std::vector<MAUS::ReconEvent*>* revts = spill->GetReconEvents();
340
341
// Loop over recon events in spill
342
for ( size_t ev = 0; ev < revts->size(); ++ev ) {
343
if ( !revts->at(ev) )
344
continue;
345
346
MAUS::TOFEvent *tofevent = revts->at(ev)->GetTOFEvent();
347
MAUS::TOFEventSpacePoint tofsp = tofevent->GetTOFEventSpacePoint();
348
std::vector<MAUS::TOFSpacePoint> tof1sps = tofsp.GetTOF1SpacePointArray();
349
std::vector<MAUS::TOFSpacePoint> tof2sps = tofsp.GetTOF2SpacePointArray();
350
351
double t1(0);
352
if ( tof1sps.size() == 1 ) {
353
t1 = tof1sps[0].GetTime();
354
} else {
355
continue;
356
}
357
358
double t2(0);
359
if ( tof2sps.size() == 1 ) {
360
t2 = tof2sps[0].GetTime();
361
} else {
362
continue;
363
}
364
365
// Display the progress in %
366
size_t j = (50*n)/(nentries-1);
367
if ( j > (size_t)(50*(n-1))/(nentries-1) || !i ) {
368
std::cerr << "[";
369
for (size_t k = 0; k < 51; k++) {
370
if ( k < j ) std::cerr << "=";
371
if ( k == j ) std::cerr << ">";
372
if ( k > j ) std::cerr << " ";
373
}
374
std::cerr << "] " << j*2;
375
if ( n != (size_t)(nentries-1) ) {
376
std::cerr << " %\r";
377
std::cerr.flush();
378
} else {
379
std::cerr << std::endl;
380
}
381
}
382
383
// Fill the TOF histogram and increment the count
384
htof12->Fill(t2 - t1);
385
n++;
386
}
387
if ( n > nentries )
388
break;
389
}
390
if ( n > nentries )
391
break;
392
}
393
if ( n < nentries )
394
std::cerr << "\n Not enough entries in the files to fit the TOF distributions" << std::endl;
395
396
// Find the peaks using the root TSpectrum code
397
const int npeaks = 2; // muons, pions
398
TSpectrum *spc = new TSpectrum(npeaks);
399
int nfound = spc->Search(htof12);
400
401
if ( npeaks != nfound ) {
402
std::cerr << "Wrong number of peaks found (" << nfound << "), aborting..." << std::endl;
403
} else {
404
float *xpeaks = spc->GetPositionX(); // array of npeaks doubles
405
std::vector<double> posx;
406
for (size_t i = 0; i < npeaks; i++)
407
posx.push_back((double)xpeaks[i]);
408
std::sort(posx.begin(), posx.end());
409
410
// Define a two peaks gaussian to fit to the TOF distribution
411
TF1 *fpeaks = new TF1("fpeaks",
412
"[0]*exp(-0.5*((x-[1])/[2])*((x-[1])/[2]))+"
413
"[3]*exp(-0.5*((x-[4])/[5])*((x-[4])/[5]))", 25, 50);
414
fpeaks->SetParNames("A_{#mu}","#mu_{#mu}","#sigma_{#mu}",
415
"A_{#pi}","#mu_{#pi}","#sigma_{#pi}");
416
417
// Set first approximation of parameters
418
for(size_t p = 0; p < npeaks; p++) {
419
int bin = htof12->GetXaxis()->FindBin(posx[p]);
420
double yp = htof12->GetBinContent(bin);
421
422
fpeaks->FixParameter(3*p, yp);
423
fpeaks->FixParameter(3*p+1, posx[p]);
424
fpeaks->SetParameter(3*p+2, .1);
425
}
426
427
// Fit peaks with the two peak gaussian and extract the parameters
428
htof12->Fit("fpeaks", "RQ");
429
for(size_t p = 0; p < npeaks; p++) {
430
mu[p] = fpeaks->GetParameter(3*p+1);
431
sig[p] = fpeaks->GetParameter(3*p+2);
432
}
433
434
// Now plot the histogram
435
gStyle->SetStatX(.9);
436
gStyle->SetStatY(.9);
437
gStyle->SetOptFit(1);
438
TCanvas *c = new TCanvas("tof12", "tof12", 1200, 800);
439
htof12->Draw();
440
c->SaveAs("tof12.png");
441
delete c;
442
}
443
}
444
445
// TRACKER ALIGNMENT
446
// Here the offsets, pitches and yaws (and roll) of the two trackers are optimized
447
for (size_t it = 0; it < iterations; it++) {
448
std::cerr << "Iteration " << it+1 << " of the tracker alignment algorithm" << std::endl;
449
450
for (size_t iFile = 2; iFile < (size_t)argc; iFile++) {
451
// Set up the ROOT file and data pointer
452
TFile f1(argv[iFile]); // Load the MAUS output file
453
TTree *T = (TTree*)f1.Get("Spill"); // Pull out the TTree
454
MAUS::Data *data_ptr = new MAUS::Data(); // A variable to store the Data from each spill
455
T->SetBranchAddress("data", &data_ptr); // Set the address of data_ptr
456
457
// Number of spills to process
458
std::cerr << "Processing \"" << argv[iFile] << "\"" << std::endl;
459
460
// Loop over the spills, get the recon events
461
for (size_t i = 0; i < (size_t)T->GetEntries(); ++i) {
462
463
// Display the progress in %
464
size_t j = (50*i)/(T->GetEntries()-1);
465
if ( j > (size_t)(50*(i-1))/(T->GetEntries()-1) || !i ) {
466
std::cerr << "[";
467
for (size_t k = 0; k < 51; k++) {
468
if ( k < j ) std::cerr << "=";
469
if ( k == j ) std::cerr << ">";
470
if ( k > j ) std::cerr << " ";
471
}
472
std::cerr << "] " << j*2;
473
if ( i != (size_t)(T->GetEntries()-1) ) {
474
std::cerr << " %\r";
475
std::cerr.flush();
476
} else {
477
std::cerr << std::endl;
478
}
479
}
480
481
// Update the spill pointed to by data_ptr
482
T->GetEntry(i);
483
MAUS::Spill* spill = data_ptr->GetSpill();
484
if (spill == NULL || !(spill->GetDaqEventType() == "physics_event"))
485
continue;
486
std::vector<MAUS::ReconEvent*>* revts = spill->GetReconEvents();
487
488
// Loop over recon events in spill
489
for ( size_t ev = 0; ev < revts->size(); ++ev ) {
490
if ( !revts->at(ev) )
491
continue;
492
493
//////////////////////////////
494
////// TOF1 COORDINATES //////
495
//////////////////////////////
496
TVector3 threeve1;
497
double t1(0);
498
MAUS::TOFEvent *tofevent = revts->at(ev)->GetTOFEvent();
499
MAUS::TOFEventSpacePoint tofsp = tofevent->GetTOFEventSpacePoint();
500
std::vector<MAUS::TOFSpacePoint> tof1sps = tofsp.GetTOF1SpacePointArray();
501
502
if ( tof1sps.size() == 1 ) {
503
threeve1.SetX(tof1sps[0].GetGlobalPosX());
504
threeve1.SetY(tof1sps[0].GetGlobalPosY());
505
threeve1.SetZ(tof1sps[0].GetGlobalPosZ());
506
t1 = tof1sps[0].GetTime();
507
} else {
508
continue;
509
}
510
511
512
//////////////////////////////
513
////// TOF2 COORDINATES //////
514
//////////////////////////////
515
TVector3 threeve2;
516
double t2(0);
517
std::vector<MAUS::TOFSpacePoint> tof2sps = tofsp.GetTOF2SpacePointArray();
518
519
if ( tof2sps.size() == 1 ) {
520
threeve2.SetX(tof2sps[0].GetGlobalPosX());
521
threeve2.SetY(tof2sps[0].GetGlobalPosY());
522
threeve2.SetZ(tof2sps[0].GetGlobalPosZ());
523
t2 = tof2sps[0].GetTime();
524
} else {
525
continue;
526
}
527
528
529
//////////////////////////////
530
//////// TOFs ANGLES /////////
531
//////////////////////////////
532
// Reconstruct the angles between the TOFs
533
double ax12 = (threeve2.x()-threeve1.x())/(threeve2.z()-threeve1.z());
534
double ay12 = (threeve2.y()-threeve1.y())/(threeve2.z()-threeve1.z());
535
536
537
//////////////////////////////
538
////////// TOF PID ///////////
539
//////////////////////////////
540
// Evaluate the time between TOF1 and TOF2 and use to reject electrons/positrons
541
// as they have a radically different scattering pattern from mu and pi
542
double c = 299792458; // [m/s]
543
double t12 = t2 - t1; // [ns]
544
double s12 = distance(threeve1, threeve2); // [mm]
545
double ct12 = 1e9*s12*1e-3/c; // [ns]
546
if ( t12 < ct12 + .5 ) // 5 sigma electron rejection
547
continue;
548
549
// Develop a momentum based fiducial area downstream
550
// Calculate the momentum and Lorentz beta of the particle assuming it's a muon
551
// Perform pid based on the time-of-flight if requested
552
double m(0.);
553
double mmu = 105.66; // [MeV/c^2]
554
double mpi = 139.57; // [MeV/c^2]
555
if ( do_pid ) {
556
if ( fabs(t12-mu[0])/sig[0] < fabs(t12-mu[1])/sig[1] ) {
557
m = mmu;
558
} else {
559
m = mpi;
560
}
561
} else {
562
m = mmu;
563
}
564
565
double p12 = m/sqrt(pow(t12*1e-9*c/(s12*1e-3), 2)-1); // [MeV/c]
566
double beta12 = 1./sqrt(pow(m/p12, 2)+1);
567
568
// Calculate the width of the scattering distribution in angle
569
// The form factor, g(x/X0), has been studied for run 07418 and yielded 0.212
570
double th0 = 2.88/(p12*beta12); // [rad]
571
572
// Define the cut based on the probability to be contained in TKD
573
// The effective distance, Delta z, has been studied for run 07418 and yielded 2935.37 mm
574
double dz_eff = 2935.37; // [mm]
575
double rad_cut = 150.-2*dz_eff*th0; // 2 sigma containment
576
577
578
//////////////////////////////
579
////// TKD&U COORDINATES /////
580
//////////////////////////////
581
TVector3 threeved, threeveu;
582
double mx[2] = {-999, -999}; // Gradients of the linear fits
583
double my[2] = {-999, -999};
584
double dz_trackers = tkd_globals.z()-tku_globals.z(); // [mm]
585
double dz_tku_tof1 = tku_globals.z()-threeve1.z();
586
double dz_tkd_tof2 = tkd_globals.z()-threeve2.z();
587
double dz_ref = -549.9513; // [mm]
588
bool fiducial(false);
589
590
// PATTERN REGOGNITION METHOD //
591
MAUS::SciFiEvent* sfevt = revts->at(ev)->GetSciFiEvent(); // Pull out scifi event
592
MAUS::SciFiStraightPRTrackPArray sfst = sfevt->straightprtracks();
593
MAUS::SciFiTrackPArray sft = sfevt->scifitracks();
594
595
if ( sfst.size() != 2 )
596
continue;
597
for ( size_t st = 0; st < sfst.size(); ++st ) { // Loop over straight tracks
598
if ( sfst[st]->get_tracker() == 0 ) { // Upstream
599
600
// Fetch the local fitting parameters (given at the ref. plane)
601
threeveu.SetX(-sfst[st]->get_x0()); // x flipped
602
threeveu.SetY( sfst[st]->get_y0()); // y conserved
603
threeveu.SetZ(-dz_ref); // z flipped
604
mx[0] = sfst[st]->get_mx(); // x flipped, z flipped
605
my[0] = -sfst[st]->get_my(); // y conserved, z flipped
606
607
// Get the global position in the centre of TKU
608
TVector3 vtku(threeveu.x()+mx[0]*dz_ref, threeveu.y()+my[0]*dz_ref, 0);
609
vtku.RotateX(alphau*1e-3);
610
vtku.RotateY(betau*1e-3);
611
vtku.RotateZ(gammau*1e-3);
612
vtku += TVector3(xu, yu, 0);
613
614
// Get the global gradients of the track in TKU
615
TVector3 gtku(mx[0], my[0], 1);
616
gtku.RotateX(alphau*1e-3);
617
gtku.RotateY(betau*1e-3);
618
gtku.RotateZ(gammau*1e-3);
619
620
// Estimate where we expect the particle to hit the centre of TKD
621
double xp = vtku.x() + dz_trackers*gtku.x()/gtku.z();
622
double yp = vtku.y() + dz_trackers*gtku.y()/gtku.z();
623
624
// For the alignment, only keep the particles that do not scatter out
625
if ( sqrt(pow(xp, 2) + pow(yp, 2)) < rad_cut ) {
626
fiducial = true;
627
}
628
} else if ( sfst[st]->get_tracker() == 1 ) { // Downstream
629
630
// Fetch the local fitting parameters (given at the ref. plane)
631
threeved.SetX(sfst[st]->get_x0());
632
threeved.SetY(sfst[st]->get_y0());
633
threeved.SetZ(dz_ref);
634
mx[1] = sfst[st]->get_mx();
635
my[1] = sfst[st]->get_my();
636
}
637
}
638
639
640
//////////////////////////////
641
////// FILL HISTOGRAMS ///////
642
//////////////////////////////
643
644
if ( fiducial ) {
645
if ( !do_roll ) {
646
htku_resxpyp->Fill(1000*(mx[0]-ax12), 1000*(my[0]-ay12));
647
htku_resxy->Fill(threeveu.x()+dz_ref*mx[0] - (threeve1.x()+ax12*dz_tku_tof1),
648
threeveu.y()+dz_ref*my[0] - (threeve1.y()+ay12*dz_tku_tof1));
649
650
htkd_resxpyp->Fill(1000*(mx[1]-ax12), 1000*(my[1]-ay12));
651
htkd_resxy->Fill(threeved.x()-dz_ref*mx[1] - (threeve2.x()+ax12*dz_tkd_tof2),
652
threeved.y()-dz_ref*my[1] - (threeve2.y()+ay12*dz_tkd_tof2));
653
} else {
654
htku_resx_y->Fill(threeveu.y()+dz_ref*my[0],
655
threeveu.x()+dz_ref*mx[0] - (threeve1.x()+ax12*dz_tku_tof1));
656
htku_resy_x->Fill(threeveu.x()+dz_ref*mx[0],
657
threeveu.y()+dz_ref*my[0] - (threeve1.y()+ay12*dz_tku_tof1));
658
htku_resxp_yp->Fill(1000*my[0], 1000*(mx[0]-ax12));
659
htku_resyp_xp->Fill(1000*mx[0], 1000*(my[0]-ay12));
660
661
htkd_resx_y->Fill(threeved.y()-dz_ref*my[1],
662
threeved.x()-dz_ref*mx[1] - (threeve2.x()+ax12*dz_tkd_tof2));
663
htkd_resy_x->Fill(threeved.x()-dz_ref*mx[1],
664
threeved.y()-dz_ref*my[1] - (threeve2.y()+ay12*dz_tkd_tof2));
665
htkd_resxp_yp->Fill(1000*my[1], 1000*(mx[1]-ax12));
666
htkd_resyp_xp->Fill(1000*mx[1], 1000*(my[1]-ay12));
667
}
668
}
669
} // End of the list of recEvents
670
} // End of the list of TTree entries
671
f1.Close();
672
} // End of the list of files
673
674
if ( !do_roll ) {
675
// Fit the distributions and get the constants from the peak
676
TF2 *gausxy = new TF2("gausxy", "[0]*TMath::Gaus(x, [1], [2])*TMath::Gaus(y, [3], [4])",
677
-250, 250, -250, 250);
678
gausxy->SetParNames("A", "#mu_{x}", "#sigma_{x}", "#mu_{y}", "#sigma_{y}");
679
680
// Upstream tracker paramters, fit the peak
681
gausxy->SetRange(-40, 40, -40, 40);
682
gausxy->SetParameters(htku_resxy->GetMaximum(), 0, 10, 0, 10);
683
htku_resxy->Fit("gausxy", "RQ0");
684
xu = -gausxy->GetParameter(1);
685
yu = -gausxy->GetParameter(3);
686
htku_optx->SetBinContent(it+1, xu);
687
htku_opty->SetBinContent(it+1, yu);
688
htku_optx->SetBinError(it+1, gausxy->GetParError(1));
689
htku_opty->SetBinError(it+1, gausxy->GetParError(3));
690
691
gausxy->SetRange(-20, 20, -20, 20);
692
gausxy->SetParameters(htku_resxpyp->GetMaximum(), 0, 10, 0, 10);
693
htku_resxpyp->Fit("gausxy", "RQ0");
694
alphau = gausxy->GetParameter(3);
695
betau = -gausxy->GetParameter(1);
696
htku_optalpha->SetBinContent(it+1, alphau);
697
htku_optbeta->SetBinContent(it+1, betau);
698
htku_optalpha->SetBinError(it+1, gausxy->GetParError(3));
699
htku_optbeta->SetBinError(it+1, gausxy->GetParError(1));
700
701
// Downstream tracker parameters, fit the peak
702
gausxy->SetRange(-40, 40, -40, 40);
703
gausxy->SetParameters(htkd_resxy->GetMaximum(), 0, 10, 0, 10);
704
htkd_resxy->Fit("gausxy", "RQ0");
705
xd = -gausxy->GetParameter(1);
706
yd = -gausxy->GetParameter(3);
707
htkd_optx->SetBinContent(it+1, xd);
708
htkd_opty->SetBinContent(it+1, yd);
709
htkd_optx->SetBinError(it+1, gausxy->GetParError(1));
710
htkd_opty->SetBinError(it+1, gausxy->GetParError(3));
711
712
gausxy->SetRange(-20, 20, -20, 20);
713
gausxy->SetParameters(htkd_resxpyp->GetMaximum(), 0, 10, 0, 10);
714
htkd_resxpyp->Fit("gausxy", "RQ0");
715
alphad = gausxy->GetParameter(3);
716
betad = -gausxy->GetParameter(1);
717
htkd_optalpha->SetBinContent(it+1, alphad);
718
htkd_optbeta->SetBinContent(it+1, betad);
719
htkd_optalpha->SetBinError(it+1, gausxy->GetParError(3));
720
htkd_optbeta->SetBinError(it+1, gausxy->GetParError(1));
721
722
} else {
723
724
// Vectors to hold the different measurements of the roll, gamma
725
std::vector<double> gamu, gamue;
726
std::vector<double> gamd, gamde;
727
728
// Get the fitting parameters for the upstream tracker
729
// res_x vs y in TKU, measures -x_U and gamma_U
730
TH1D *htku_resx_y_1D = new TH1D("", "", htku_resx_y->GetNbinsX(),
731
htku_resx_y->GetXaxis()->GetXmin(),
732
htku_resx_y->GetXaxis()->GetXmax());
733
for (int k = 0; k < htku_resx_y->GetNbinsX(); k++) {
734
TH1D *hbinsy = htku_resx_y->ProjectionY(Form("bin%d",k+1), k+1, k+1);
735
double mean = hbinsy->GetMean();
736
double meanerr = hbinsy->GetMeanError();
737
738
htku_resx_y_1D->SetBinContent(k+1, mean);
739
htku_resx_y_1D->SetBinError(k+1, meanerr);
740
}
741
742
htku_resx_y_1D->Fit("pol1", "Q", "",
743
htku_resx_y->GetMean(1)-2*htku_resx_y->GetRMS(1),
744
htku_resx_y->GetMean(1)+2*htku_resx_y->GetRMS(1));
745
xu = -htku_resx_y_1D->GetFunction("pol1")->GetParameter(0);
746
htku_optx->SetBinContent(it+1, xu);
747
htku_optx->SetBinError(it+1, htku_resx_y_1D->GetFunction("pol1")->GetParError(0));
748
gamu.push_back(htku_resx_y_1D->GetFunction("pol1")->GetParameter(1));
749
gamue.push_back(htku_resx_y_1D->GetFunction("pol1")->GetParError(1));
750
delete htku_resx_y_1D;
751
752
// res_y vs x in TKU, measures -y_U and -gamma_U
753
TH1D *htku_resy_x_1D = new TH1D("", "", htku_resy_x->GetNbinsX(),
754
htku_resy_x->GetXaxis()->GetXmin(),
755
htku_resy_x->GetXaxis()->GetXmax());
756
for (int k = 0; k <htku_resy_x->GetNbinsX(); k++) {
757
TH1D *hbinsy = htku_resy_x->ProjectionY(Form("bin%d",k+1), k+1, k+1);
758
double mean = hbinsy->GetMean();
759
double meanerr = hbinsy->GetMeanError();
760
761
htku_resy_x_1D->SetBinContent(k+1, mean);
762
htku_resy_x_1D->SetBinError(k+1, meanerr);
763
}
764
765
htku_resy_x_1D->Fit("pol1", "Q", "",
766
htku_resy_x->GetMean(1)-2*htku_resy_x->GetRMS(1),
767
htku_resy_x->GetMean(1)+2*htku_resy_x->GetRMS(1));
768
yu = -htku_resy_x_1D->GetFunction("pol1")->GetParameter(0);
769
htku_opty->SetBinContent(it+1, yu);
770
htku_opty->SetBinError(it+1, htku_resy_x_1D->GetFunction("pol1")->GetParError(0));
771
gamu.push_back(-htku_resy_x_1D->GetFunction("pol1")->GetParameter(1));
772
gamue.push_back(htku_resy_x_1D->GetFunction("pol1")->GetParError(1));
773
delete htku_resy_x_1D;
774
775
// res_xp vs yp in TKU, measures -beta_U and gamma_U
776
TH1D *htku_resxp_yp_1D = new TH1D("", "", htku_resxp_yp->GetNbinsX(),
777
htku_resxp_yp->GetXaxis()->GetXmin(),
778
htku_resxp_yp->GetXaxis()->GetXmax());
779
for (int k = 0; k <htku_resxp_yp->GetNbinsX(); k++) {
780
TH1D *hbinsy = htku_resxp_yp->ProjectionY(Form("bin%d",k+1), k+1, k+1);
781
double mean = hbinsy->GetMean();
782
double meanerr = hbinsy->GetMeanError();
783
784
htku_resxp_yp_1D->SetBinContent(k+1, mean);
785
htku_resxp_yp_1D->SetBinError(k+1, meanerr);
786
}
787
788
htku_resxp_yp_1D->Fit("pol1", "Q", "",
789
htku_resxp_yp->GetMean(1)-2*htku_resxp_yp->GetRMS(1),
790
htku_resxp_yp->GetMean(1)+2*htku_resxp_yp->GetRMS(1));
791
betau = -htku_resxp_yp_1D->GetFunction("pol1")->GetParameter(0);
792
htku_optbeta->SetBinContent(it+1, betau);
793
htku_optbeta->SetBinError(it+1, htku_resxp_yp_1D->GetFunction("pol1")->GetParError(0));
794
gamu.push_back(htku_resxp_yp_1D->GetFunction("pol1")->GetParameter(1));
795
gamue.push_back(htku_resxp_yp_1D->GetFunction("pol1")->GetParError(1));
796
delete htku_resxp_yp_1D;
797
798
// res_yp vs xp in TKU, measures alpha_U and -gamma_U
799
TH1D *htku_resyp_xp_1D = new TH1D("", "", htku_resyp_xp->GetNbinsX(),
800
htku_resyp_xp->GetXaxis()->GetXmin(),
801
htku_resyp_xp->GetXaxis()->GetXmax());
802
for (int k = 0; k <htku_resyp_xp->GetNbinsX(); k++) {
803
TH1D *hbinsy = htku_resyp_xp->ProjectionY(Form("bin%d",k+1), k+1, k+1);
804
double mean = hbinsy->GetMean();
805
double meanerr = hbinsy->GetMeanError();
806
807
htku_resyp_xp_1D->SetBinContent(k+1, mean);
808
htku_resyp_xp_1D->SetBinError(k+1, meanerr);
809
}
810
811
htku_resyp_xp_1D->Fit("pol1", "Q", "",
812
htku_resyp_xp->GetMean(1)-2*htku_resyp_xp->GetRMS(1),
813
htku_resyp_xp->GetMean(1)+2*htku_resyp_xp->GetRMS(1));
814
alphau = htku_resyp_xp_1D->GetFunction("pol1")->GetParameter(0);
815
htku_optalpha->SetBinContent(it+1, alphau);
816
htku_optalpha->SetBinError(it+1, htku_resyp_xp_1D->GetFunction("pol1")->GetParError(0));
817
gamu.push_back(-htku_resyp_xp_1D->GetFunction("pol1")->GetParameter(1));
818
gamue.push_back(htku_resyp_xp_1D->GetFunction("pol1")->GetParError(1));
819
delete htku_resyp_xp_1D;
820
821
// res_x vs y in TKD, measures -x_D and gamma_D
822
TH1D *htkd_resx_y_1D = new TH1D("", "", htkd_resx_y->GetNbinsX(),
823
htkd_resx_y->GetXaxis()->GetXmin(),
824
htkd_resx_y->GetXaxis()->GetXmax());
825
for (int k = 0; k <htkd_resx_y->GetNbinsX(); k++) {
826
TH1D *hbinsy = htkd_resx_y->ProjectionY(Form("bin%d",k+1), k+1, k+1);
827
double mean = hbinsy->GetMean();
828
double meanerr = hbinsy->GetMeanError();
829
830
htkd_resx_y_1D->SetBinContent(k+1, mean);
831
htkd_resx_y_1D->SetBinError(k+1, meanerr);
832
}
833
834
htkd_resx_y_1D->Fit("pol1", "Q", "",
835
htkd_resx_y->GetMean(1)-2*htkd_resx_y->GetRMS(1),
836
htkd_resx_y->GetMean(1)+2*htkd_resx_y->GetRMS(1));
837
xd = -htkd_resx_y_1D->GetFunction("pol1")->GetParameter(0);
838
htkd_optx->SetBinContent(it+1, xd);
839
htkd_optx->SetBinError(it+1, htkd_resx_y_1D->GetFunction("pol1")->GetParError(0));
840
gamd.push_back(htkd_resx_y_1D->GetFunction("pol1")->GetParameter(1));
841
gamde.push_back(htkd_resx_y_1D->GetFunction("pol1")->GetParError(1));
842
delete htkd_resx_y_1D;
843
844
// res_y vs x in TKD, measures -x_D and gamma_D
845
TH1D *htkd_resy_x_1D = new TH1D("", "", htkd_resy_x->GetNbinsX(),
846
htkd_resy_x->GetXaxis()->GetXmin(),
847
htkd_resy_x->GetXaxis()->GetXmax());
848
for (int k = 0; k < htkd_resy_x->GetNbinsX(); k++) {
849
TH1D *hbinsy = htkd_resy_x->ProjectionY(Form("bin%d",k+1), k+1, k+1);
850
double mean = hbinsy->GetMean();
851
double meanerr = hbinsy->GetMeanError();
852
853
htkd_resy_x_1D->SetBinContent(k+1, mean);
854
htkd_resy_x_1D->SetBinError(k+1, meanerr);
855
}
856
857
htkd_resy_x_1D->Fit("pol1", "Q", "",
858
htkd_resy_x->GetMean(1)-2*htkd_resy_x->GetRMS(1),
859
htkd_resy_x->GetMean(1)+2*htkd_resy_x->GetRMS(1));
860
yd = -htkd_resy_x_1D->GetFunction("pol1")->GetParameter(0);
861
htkd_opty->SetBinContent(it+1, yd);
862
htkd_opty->SetBinError(it+1, htkd_resy_x_1D->GetFunction("pol1")->GetParError(0));
863
gamd.push_back(-htkd_resy_x_1D->GetFunction("pol1")->GetParameter(1));
864
gamde.push_back(htkd_resy_x_1D->GetFunction("pol1")->GetParError(1));
865
delete htkd_resy_x_1D;
866
867
// res_xp vs yp in TKD, measures -x_D and gamma_D
868
TH1D *htkd_resxp_yp_1D = new TH1D("", "", htkd_resxp_yp->GetNbinsX(),
869
htkd_resxp_yp->GetXaxis()->GetXmin(),
870
htkd_resxp_yp->GetXaxis()->GetXmax());
871
for (int k = 0; k < htkd_resxp_yp->GetNbinsX(); k++) {
872
TH1D *hbinsy = htkd_resxp_yp->ProjectionY(Form("bin%d",k+1), k+1, k+1);
873
double mean = hbinsy->GetMean();
874
double meanerr = hbinsy->GetMeanError();
875
876
htkd_resxp_yp_1D->SetBinContent(k+1, mean);
877
htkd_resxp_yp_1D->SetBinError(k+1, meanerr);
878
}
879
htkd_resxp_yp_1D->Fit("pol1", "Q", "",
880
htkd_resxp_yp->GetMean(1)-2*htkd_resxp_yp->GetRMS(1),
881
htkd_resxp_yp->GetMean(1)+2*htkd_resxp_yp->GetRMS(1));
882
betad = -htkd_resxp_yp_1D->GetFunction("pol1")->GetParameter(0);
883
htkd_optbeta->SetBinContent(it+1, betad);
884
htkd_optbeta->SetBinError(it+1, htkd_resxp_yp_1D->GetFunction("pol1")->GetParError(0));
885
gamd.push_back(htkd_resxp_yp_1D->GetFunction("pol1")->GetParameter(1));
886
gamde.push_back(htkd_resxp_yp_1D->GetFunction("pol1")->GetParError(1));
887
delete htkd_resxp_yp_1D;
888
889
// res_yp vs xp in TKD, measures -x_D and gamma_D
890
TH1D *htkd_resyp_xp_1D = new TH1D("", "", htkd_resyp_xp->GetNbinsX(),
891
htkd_resyp_xp->GetXaxis()->GetXmin(),
892
htkd_resyp_xp->GetXaxis()->GetXmax());
893
for (int k = 0; k < htkd_resyp_xp->GetNbinsX(); k++) {
894
TH1D *hbinsy = htkd_resyp_xp->ProjectionY(Form("bin%d",k+1), k+1, k+1);
895
double mean = hbinsy->GetMean();
896
double meanerr = hbinsy->GetMeanError();
897
898
htkd_resyp_xp_1D->SetBinContent(k+1, mean);
899
htkd_resyp_xp_1D->SetBinError(k+1, meanerr);
900
}
901
902
htkd_resyp_xp_1D->Fit("pol1", "Q", "",
903
htkd_resyp_xp->GetMean(1)-2*htkd_resyp_xp->GetRMS(1),
904
htkd_resyp_xp->GetMean(1)+2*htkd_resyp_xp->GetRMS(1));
905
alphad = htkd_resyp_xp_1D->GetFunction("pol1")->GetParameter(0);
906
htkd_optalpha->SetBinContent(it+1, alphad);
907
htkd_optalpha->SetBinError(it+1, htkd_resyp_xp_1D->GetFunction("pol1")->GetParError(0));
908
gamd.push_back(-htkd_resyp_xp_1D->GetFunction("pol1")->GetParameter(1));
909
gamde.push_back(htkd_resyp_xp_1D->GetFunction("pol1")->GetParError(1));
910
delete htkd_resyp_xp_1D;
911
912
// Do a weighted average to find the best fit to the roll
913
double gamu_mean(0.), gamd_mean(0.), gamu_error(0.), gamd_error(0.);
914
for (size_t sample = 0; sample < gamu.size(); sample++) {
915
gamu_mean += gamu[sample]*pow(gamue[sample], -2);
916
gamd_mean += gamd[sample]*pow(gamde[sample], -2);
917
918
gamu_error += pow(gamue[sample], -2);
919
gamd_error += pow(gamde[sample], -2);
920
}
921
922
gamu_mean /= gamu_error;
923
gamd_mean /= gamd_error;
924
925
gammau = 1000*gamu_mean;
926
htku_optgamma->SetBinContent(it+1, gammau);
927
htku_optgamma->SetBinError(it+1, 1000./sqrt(gamu_error));
928
929
gammad = 1000*gamd_mean;
930
htkd_optgamma->SetBinContent(it+1, gammad);
931
htkd_optgamma->SetBinError(it+1, 1000./sqrt(gamd_error));
932
}
933
934
// Scale the histograms to 0 for the next round
935
if ( it < iterations-1 ) {
936
if ( !do_roll ) {
937
htku_resxy->Reset();
938
htkd_resxy->Reset();
939
htku_resxpyp->Reset();
940
htkd_resxpyp->Reset();
941
} else {
942
htku_resx_y->Reset();
943
htku_resy_x->Reset();
944
htku_resxp_yp->Reset();
945
htku_resyp_xp->Reset();
946
947
htkd_resx_y->Reset();
948
htkd_resy_x->Reset();
949
htkd_resxp_yp->Reset();
950
htkd_resyp_xp->Reset();
951
}
952
}
953
} // End of the fitting iterations
954
955
// Set the Stat box to be along the top right corner of the TPad in the TCanvas
956
gStyle->SetStatX(.9);
957
gStyle->SetStatY(.9);
958
959
// Labels of the fits for the plots that include the roll
960
std::vector<std::string> yint_labels;
961
yint_labels.push_back("-x_{U}");
962
yint_labels.push_back("-y_{U}");
963
yint_labels.push_back("-#beta_{U}");
964
yint_labels.push_back("#alpha_{U}");
965
yint_labels.push_back("-x_{D}");
966
yint_labels.push_back("-y_{D}");
967
yint_labels.push_back("-#beta_{D}");
968
yint_labels.push_back("#alpha_{D}");
969
std::vector<std::string> grad_labels;
970
grad_labels.push_back("#gamma_{U}");
971
grad_labels.push_back("-#gamma_{U}");
972
grad_labels.push_back("#gamma_{U}");
973
grad_labels.push_back("-#gamma_{U}");
974
grad_labels.push_back("#gamma_{D}");
975
grad_labels.push_back("-#gamma_{D}");
976
grad_labels.push_back("#gamma_{D}");
977
grad_labels.push_back("-#gamma_{D}");
978
979
// Plot the graphs and their fit, store them to a root file
980
std::string output_name = run_name+"_alignment.root";
981
TFile* output = new TFile(output_name.c_str(), "RECREATE");
982
for (size_t id = 0; id < res_hists.size(); id++) {
983
if ( !do_roll ) {
984
gStyle->SetOptFit(1);
985
TF2 *gausxy = new TF2("gausxy", "[0]*TMath::Gaus(x, [1], [2])*TMath::Gaus(y, [3], [4])",
986
-500, 500, -500, 500);
987
gausxy->SetParNames("A", "#mu_{x}", "#sigma_{x}", "#mu_{y}", "#sigma_{y}");
988
TCanvas *cd = new TCanvas("cd", "cd", 1200, 900);
989
res_hists[id]->Draw("COLZ");
990
res_hists[id]->Fit("gausxy");
991
cd->SaveAs(TString(res_hists[id]->GetName())+".png");
992
delete cd;
993
994
res_hists[id]->Write();
995
} else {
996
double minx = res_hists_roll[id]->GetXaxis()->GetXmin();
997
double maxx = res_hists_roll[id]->GetXaxis()->GetXmax();
998
999
TString name_1D = TString(res_hists_roll[id]->GetName())+"_1D";
1000
TH1D *res_hists_1D = new TH1D(name_1D.Data(), res_hists_roll[id]->GetTitle(),
1001
res_hists_roll[id]->GetNbinsX(), minx, maxx);
1002
res_hists_1D->GetXaxis()->SetTitle(res_hists_roll[id]->GetXaxis()->GetTitle());
1003
res_hists_1D->GetYaxis()->SetTitle(res_hists_roll[id]->GetYaxis()->GetTitle());
1004
1005
for (int k = 0; k < res_hists[id]->GetNbinsX(); k++) {
1006
TH1D *hbinsy = res_hists[id]->ProjectionY(Form("bin%d",k+1), k+1, k+1);
1007
1008
double mean = hbinsy->GetMean();
1009
double meanerr = hbinsy->GetMeanError();
1010
1011
res_hists_1D->SetBinContent(k+1, mean);
1012
res_hists_1D->SetBinError(k+1, meanerr);
1013
}
1014
1015
// Draw the distributions
1016
gStyle->SetOptStat(10);
1017
gStyle->SetOptFit(1);
1018
TCanvas *cd = new TCanvas("cd", "cd", 1200, 900);
1019
res_hists_roll[id]->Draw("COLZ");
1020
cd->SaveAs(TString(res_hists_roll[id]->GetName())+".png");
1021
delete cd;
1022
1023
// Draw the mean RMS and the best fit
1024
cd = new TCanvas("cd", "cd", 1200, 900);
1025
res_hists_1D->Draw("");
1026
res_hists_1D->Fit("pol1", "Q", "",
1027
res_hists_roll[id]->GetMean(1)-2*res_hists_roll[id]->GetRMS(1),
1028
res_hists_roll[id]->GetMean(1)+2*res_hists_roll[id]->GetRMS(1));
1029
res_hists_1D->GetFunction("pol1")->SetParNames(yint_labels[id].c_str(),
1030
grad_labels[id].c_str());
1031
cd->SaveAs(name_1D+".png");
1032
delete cd;
1033
1034
res_hists_roll[id]->Write();
1035
res_hists_1D->Write();
1036
delete res_hists_1D;
1037
}
1038
}
1039
1040
// Sensitivity of the fit to the initial parameters
1041
gStyle->SetOptStat(0);
1042
for (size_t id = 0; id < ite_hists.size(); id++) {
1043
double mean = ite_hists[id]->GetBinContent(iterations);
1044
double err = ite_hists[id]->GetBinError(iterations);
1045
ite_hists[id]->SetTitle(TString(ite_hists[id]->GetTitle())+
1046
TString::Format(": %1.3f +/- % 1.3f", mean, err));
1047
TCanvas *c = new TCanvas("c", "c", 1200, 900);
1048
ite_hists[id]->Draw("E");
1049
c->SaveAs(TString(ite_hists[id]->GetName())+".png");
1050
delete c;
1051
1052
ite_hists[id]->Write();
1053
}
1054
1055
output->Close();
1056
}
Loggerhead is a web-based interface for Breezy
Version: 2.0.1