MiniballSort
Loading...
Searching...
No Matches
MiniballAngleFitter.cc
Go to the documentation of this file.
1#include "MiniballAngleFitter.hh"
2
3MiniballAngleFunction::MiniballAngleFunction( std::shared_ptr<MiniballSettings> _myset, std::shared_ptr<MiniballReaction> _myreact ) {
4
5 // Settings and reaction files
6 myset = _myset;
7 myreact = _myreact;
8
9 // Now initialise everything
10 Initialise();
11
12};
13
14void MiniballAngleFunction::Initialise() {
15
16 // Loop over clusters
17 present.resize( myset->GetNumberOfMiniballClusters() );
18 phiconst.resize( myset->GetNumberOfMiniballClusters() );
19 energy.resize( myset->GetNumberOfMiniballClusters() );
20 err.resize( myset->GetNumberOfMiniballClusters() );
21 phic.resize( myset->GetNumberOfMiniballClusters() );
22 phie.resize( myset->GetNumberOfMiniballClusters() );
23 for( unsigned int clu = 0; clu < myset->GetNumberOfMiniballClusters(); ++clu ) {
24
25 cluster.push_back(false);
26 present[clu].resize( myset->GetNumberOfMiniballCrystals() );
27 phiconst[clu].resize( myset->GetNumberOfMiniballCrystals() );
28 energy[clu].resize( myset->GetNumberOfMiniballCrystals() );
29 err[clu].resize( myset->GetNumberOfMiniballCrystals() );
30 phic[clu].resize( myset->GetNumberOfMiniballCrystals() );
31 phie[clu].resize( myset->GetNumberOfMiniballCrystals() );
32
33 // Loop over crystals
34 for( unsigned int cry = 0; cry < myset->GetNumberOfMiniballCrystals(); ++cry ) {
35
36 present[clu][cry].resize( myset->GetNumberOfMiniballSegments() );
37 phiconst[clu][cry].resize( myset->GetNumberOfMiniballSegments() );
38 energy[clu][cry].resize( myset->GetNumberOfMiniballSegments() );
39 err[clu][cry].resize( myset->GetNumberOfMiniballSegments() );
40 phic[clu][cry].resize( myset->GetNumberOfMiniballSegments() );
41 phie[clu][cry].resize( myset->GetNumberOfMiniballSegments() );
42
43 // Loop over segments
44 for( unsigned int seg = 0; seg < myset->GetNumberOfMiniballSegments(); ++seg ) {
45
46 present[clu][cry].push_back(false);
47 phiconst[clu][cry].push_back(false);
48 energy[clu][cry].push_back(0.0);
49 err[clu][cry].push_back(9.9);
50 phic[clu][cry].push_back(0.0);
51 phie[clu][cry].push_back(9.9);
52
53 } // seg
54
55 } // cry
56
57 } // clu
58
59 // Set the user z
60 user_z = myreact->GetOffsetZ();
61
62
63}
64
65bool MiniballAngleFunction::FitPeak( TH1D *h, double &en, double &er ){
66
67 // Work out some limits
68 double low_lim = en*0.92 - 5.;
69 double upp_lim = en*1.08 + 5.;
70 if( eref < 500. && upp_lim > 500. ) upp_lim = 500.;
71
72 // Draw the thing
73 auto c1 = std::make_unique<TCanvas>();
74 h->GetXaxis()->SetRangeUser( low_lim - 20., upp_lim + 20. );
75 h->Draw();
76
77 // Make a TF1
78 auto peakfit = std::make_unique<TF1>( "peakfit", "gaus(0)+pol1(3)", low_lim, upp_lim );
79
80 // Set initial parameters
81 double en_est = h->GetBinCenter( h->GetMaximumBin() );
82 double sig_est = 1.7; // 1.7 keV guess for sigma = 4.0 keV for FWHM
83 double bg_est = h->GetBinContent( h->FindBin( en_est - 5.0 * sig_est ) );
84 double amp_est = h->GetBinContent( h->GetMaximumBin() );
85 amp_est -= bg_est;
86 if( amp_est < 1.0 ) amp_est = h->GetBinContent( h->GetMaximumBin() ) + 1.0;
87 peakfit->SetParameter( 0, amp_est ); // amplitude
88 peakfit->SetParameter( 1, en_est ); // centroid
89 peakfit->SetParameter( 2, sig_est ); // sigma width
90 peakfit->SetParameter( 3, bg_est ); // bg const
91 peakfit->SetParameter( 4, 1e-9 ); // bg gradient
92
93 // Parameter limits
94 double integral = h->Integral( h->FindBin(low_lim), h->FindBin(upp_lim) );
95 peakfit->SetParLimits( 0, 1.0, integral ); // amplitude limit
96 peakfit->SetParLimits( 1, low_lim, upp_lim ); // centroid limit
97 peakfit->SetParLimits( 2, sig_est*0.25, sig_est*4.0 ); // sigma limit
98
99 // Make fit and check for success
100 gErrorIgnoreLevel = kError;
101 auto res = h->Fit( peakfit.get(), "LRSQ" );
102 c1->Print("peak_fits.pdf");
103 gErrorIgnoreLevel = kInfo;
104 if( res == 0 && integral > 150. ) {
105
106 // Set parameters back
107 en = peakfit->GetParameter(1);
108 er = peakfit->GetParError(1);
109 return true;
110
111 }
112
113 else return false;
114
115}
116
117void MiniballAngleFunction::FitSegmentEnergies( TFile *infile ){
118
119 // Names of the spectra in the events file
120 std::string hname, folder = "/miniball/cluster_", base = "mb_en_core_seg_";
121
122 // Histogram objects
123 TH1D *h1 = nullptr;
124 TH2D *h2 = nullptr;
125
126 // Open the pdf file for peak fits
127 auto c1 = std::make_unique<TCanvas>();
128 gErrorIgnoreLevel = kWarning;
129 c1->Print("peak_fits.pdf(");
130
131 // Loop over all clusters
132 for( unsigned int clu = 0; clu < myset->GetNumberOfMiniballClusters(); ++clu ) {
133
134 // Loop over crystals
135 for( unsigned int cry = 0; cry < myset->GetNumberOfMiniballCrystals(); ++cry ) {
136
137 // Loop over segments
138 for( unsigned int seg = 1; seg < myset->GetNumberOfMiniballSegments(); ++seg ) {
139
140 // Build up the 2D histogram name
141 hname = folder + std::to_string(clu) + "/";
142 hname += base + std::to_string(clu) + "_";
143 hname += std::to_string(cry) + "_ebis_on";
144
145 // Get 2D histogram from file
146 h2 = static_cast<TH2D*>( infile->Get( hname.data() ) );
147
148 // Build up the 1D histogram name
149 hname = std::to_string(clu) + "_";
150 hname += std::to_string(cry) + "_";
151 hname += std::to_string(seg) + "_ebis_on";
152
153 // Project a 1D histogram from this
154 h1 = static_cast<TH1D*>( h2->ProjectionY( hname.data(), seg+1, seg+1 ) );
155
156 // Check if we have some counts
157 if( h1->Integral() > 0 ) {
158 cluster[clu] = true;
159 present[clu][cry][seg] = true;
160 }
161 else {
162 present[clu][cry][seg] = false;
163 continue;
164 }
165
166 // Predict the centroid from intial guesses
167 double en_init = myreact->DopplerShift( eref, myreact->GetBeta(),
168 TMath::Cos( myreact->GetGammaTheta( clu, cry, seg ) ) );
169
170 // Fit peak and update the energy and error
171 energy[clu][cry][seg] = en_init;
172 if( !FitPeak( h1, energy[clu][cry][seg], err[clu][cry][seg] ) )
173 present[clu][cry][seg] = false;
174
175 // Don't include data with big errors, probably a fitting issue
176 if( err[clu][cry][seg] > 0.6 )
177 present[clu][cry][seg] = false;
178
179 } // seg
180
181 } // cry
182
183 } // clu
184
185 // Close the pdf file for peak fits
186 c1->Print("peak_fits.pdf)");
187
188 // Clean up
189 if (h1) delete h1;
190 if (h2) delete h2;
191
192 gErrorIgnoreLevel = kInfo;
193
194}
195
196void MiniballAngleFunction::LoadExpEnergies( std::string energy_file ){
197
198 // Open file
199 std::ifstream energyfile( energy_file );
200
201 // Loop over the whole file line-by-line
202 std::string data_line;
203 while( std::getline( energyfile, data_line ) ) {
204
205 // Make a string stream to get the data
206 std::stringstream data_stream( data_line );
207
208 // And the package it up in to a vector
209 int parsed_int;
210 double parsed_double;
211 std::vector<int> detid;
212 std::vector<double> data;
213
214 // Just the detector id first
215 while( data_stream >> parsed_int ) {
216
217 detid.push_back( parsed_int );
218 if( detid.size() == 3 ) break;
219
220 }
221
222 // If we didn't get a full ID, skip it
223 if( detid.size() != 3 ) continue;
224
225 // Otherwise we can split the id
226 int clu = detid[0], cry = detid[1], seg = detid[2];
227
228 // Check the cluster number is sensible
229 if( clu >= (int)myset->GetNumberOfMiniballClusters() || clu < 0 ) {
230
231 std::cerr << "Bad cluster number = " << clu << std::endl;
232 continue;
233
234 }
235
236 // Check the crystal number is sensible
237 if( cry >= (int)myset->GetNumberOfMiniballCrystals() || cry < 0 ) {
238
239 std::cerr << "Bad crystal number = " << cry << std::endl;
240 continue;
241
242 }
243
244 // Check the segment number is sensible
245 if( seg >= (int)myset->GetNumberOfMiniballSegments() || seg < 0 ) {
246
247 std::cerr << "Bad segment number = " << seg << std::endl;
248 continue;
249
250 }
251
252 // Then let's take the data for that segment
253 while( data_stream >> parsed_double )
254 data.push_back( parsed_double );
255
256 // We should have 2 or 4 items
257 if( data.size() != 2 && data.size() != 4 ) continue;
258
259 // Then we at least have the energy
260 present[clu][cry][seg] = true;
261 energy[clu][cry][seg] = data[0];
262 err[clu][cry][seg] = data[1];
263 cluster[clu] = true;
264
265 // If we have a phi constraint, we have 4 data items
266 if( data.size() == 4 ) {
267
268 phiconst[clu][cry][seg] = true;
269 phic[clu][cry][seg] = data[2];
270 phie[clu][cry][seg] = data[3];
271
272 }
273
274 }
275
276 energyfile.close();
277
278};
279
280void MiniballAngleFunction::SaveExpEnergies( std::string energy_file ){
281
282 // Open file
283 std::ofstream energyfile( energy_file, std::ios::out );
284
285 // Loop over all clusters
286 for( unsigned int clu = 0; clu < myset->GetNumberOfMiniballClusters(); ++clu ) {
287
288 // Loop over crystals
289 for( unsigned int cry = 0; cry < myset->GetNumberOfMiniballCrystals(); ++cry ) {
290
291 // Loop over segments
292 for( unsigned int seg = 1; seg < myset->GetNumberOfMiniballSegments(); ++seg ) {
293
294 if( present[clu][cry][seg] ){
295
296 energyfile << clu << "\t" << cry << "\t" << seg << "\t";
297 energyfile << energy[clu][cry][seg] << "\t" << err[clu][cry][seg] << std::endl;
298
299 }
300
301 } // seg
302
303 } // cry
304
305 } // clu
306
307 energyfile.close();
308
309}
310
311double MiniballAngleFunction::operator() ( const double *p ) {
312
313 // change angles
314 // p[0] = beta
315 // p[1] = theta1, p[2] = phi1, p[3] = alpha1, p[4] = r1,
316 // p[5] = theta2, p[6] = phi2, p[7] = alpha2, p[8] = r2, etc.
317
318 // Loop over clusters
319 for( unsigned int clu = 0; clu < myset->GetNumberOfMiniballClusters(); ++clu ) {
320
321 int indx = 1 + 4*clu;
322 myreact->SetupCluster( clu, p[indx], p[indx+1], p[indx+2], p[indx+3], user_z );
323
324 }
325
326 double chisq = 0;
327
328 // Loop over clusters
329 for( unsigned int clu = 0; clu < myset->GetNumberOfMiniballClusters(); ++clu ) {
330
331 // Loop over crystals
332 for( unsigned int cry = 0; cry < myset->GetNumberOfMiniballCrystals(); ++cry ) {
333
334 // Loop over segments
335 for( unsigned int seg = 0; seg < myset->GetNumberOfMiniballSegments(); ++seg ) {
336
337 // Only include segments with data
338 if( !present[clu][cry][seg] ) continue;
339
340 // Get the theta and phi of the reaction
341 double theta = myreact->GetGammaTheta( clu, cry, seg );
342 double phi = myreact->GetGammaPhi( clu, cry, seg );
343
344 // Phi should be compared to user input, which will be 0Ëš-360Ëš
345 phi *= TMath::RadToDeg();
346 if( phi < 0 ) phi += 360.;
347
348 // Zeroth parameter is beta
349 double beta = p[0];
350
351 // Doppler corrected energy
352 double edop = myreact->DopplerShift( eref, beta, TMath::Cos(theta) );
353
354 // Compare to fitted energy, increase chisq value
355 chisq += TMath::Power( ( energy[clu][cry][seg] - edop ) / err[clu][cry][seg], 2.0 );
356
357 // Check if we have a phi constraint for this segment
358 if( !phiconst[clu][cry][seg] ) continue;
359
360 // If so, then increase chisq value
361 chisq += TMath::Power( ( phic[clu][cry][seg] - phi ) / phie[clu][cry][seg], 2.0 );
362
363 }
364
365 }
366
367 }
368
369 return chisq;
370
371}
372
373MiniballAngleFitter::MiniballAngleFitter() {
374
375 // Just use dummy files
376 MiniballAngleFitter( "default", "default" );
377
378}
379
380MiniballAngleFitter::MiniballAngleFitter( std::string settings_file, std::string reaction_file ) {
381
382 // Settings and reaction files
383 myset = std::make_shared<MiniballSettings>( settings_file );
384 myreact = std::make_shared<MiniballReaction>( reaction_file, myset );
385
386 // Now call main function
387 MiniballAngleFitter( myset, myreact );
388
389}
390
391MiniballAngleFitter::MiniballAngleFitter( std::shared_ptr<MiniballSettings> _myset, std::shared_ptr<MiniballReaction> _myreact ) {
392
393 // Settings and reaction files
394 myset = _myset ;
395 myreact = _myreact;
396
397 // Now initialise
398 Initialise();
399
400}
401
402// Input ROOT file
403bool MiniballAngleFitter::SetInputROOTFile( std::string fname ){
404
405 // Open input file
406 input_root_file = new TFile( fname.data(), "read" );
407 if( input_root_file->IsZombie() ) {
408
409 std::cout << "Cannot open " << fname << std::endl;
410 return false;
411
412 }
413
414 // Set the flag for fitting the peaks from the ROOT file
415 flag_fit_peaks = true;
416
417 // Of course, when it works, we should return true here
418 return true;
419
420}
421
422// Input data file
423bool MiniballAngleFitter::SetInputEnergiesFile( std::string fname ){
424
425 // Open input file.
426 input_data_filename = fname;
427 std::ifstream ftest;
428 ftest.open( input_data_filename );
429 if( !ftest.is_open() ) {
430
431 std::cout << "Cannot open " << input_data_filename << std::endl;
432 return false;
433
434 }
435 ftest.close();
436
437 // Set the flag for taking the data from a list file
438 flag_fit_peaks = false;
439
440 return true;
441
442}
443
444void MiniballAngleFitter::Initialise() {
445
446 // Setup the fit function
447 ff = MiniballAngleFunction( myset, myreact );
448
449 // Resize the vectors for parameters and limits
450 npars = 1 + 4 * myset->GetNumberOfMiniballClusters();
451 pars.resize( npars );
452 names.resize( npars );
453 LL.resize( npars );
454 UL.resize( npars );
455
456 // set the initial velocity from the reaction file
457 pars[0] = myreact->GetBeta();
458 names[0] = "beta";
459 LL[0] = myreact->GetBeta() * 0.5;
460 UL[0] = myreact->GetBeta() * 1.2;
461
462 // Set the initial guesses for angles of each cluster from the reaction file
463 for( unsigned int clu = 0; clu < myset->GetNumberOfMiniballClusters(); ++clu ) {
464
465 // work out the parameter numbers
466 int indx = 1 + 4 * clu;
467
468 // get starting guesses from the reaction file
469 pars[indx+0] = myreact->GetMiniballTheta(clu) * TMath::RadToDeg();
470 pars[indx+1] = myreact->GetMiniballPhi(clu) * TMath::RadToDeg();
471 pars[indx+2] = myreact->GetMiniballAlpha(clu) * TMath::RadToDeg();
472 pars[indx+3] = myreact->GetMiniballR(clu);
473
474 // Names of the parameters
475 names[indx+0] = "theta_" + std::to_string(clu);
476 names[indx+1] = "phi_" + std::to_string(clu);
477 names[indx+2] = "alpha_" + std::to_string(clu);
478 names[indx+3] = "r_" + std::to_string(clu);
479
480 // Lower limits
481 // these may cause issue if we need to cross over the 0/360 boundary
482 LL[indx+0] = pars[indx+0] - 60.0;
483 LL[indx+1] = pars[indx+1] - 60.0;
484 LL[indx+2] = 0.0;
485 LL[indx+3] = 50.0;
486
487 // Upper limits
488 UL[indx+0] = pars[indx+0] + 60.0;
489 UL[indx+1] = pars[indx+1] + 60.0;
490 UL[indx+2] = 360.0;
491 UL[indx+3] = 180.0;
492
493 // Couple of consistency checks
494 if( LL[indx+0] < 0.0 ) LL[indx+0] = 0.0;
495 if( LL[indx+1] < 0.0 ) LL[indx+1] = 0.0;
496 if( UL[indx+0] > 360.0 ) UL[indx+0] = 360.0;
497 if( UL[indx+1] > 360.0 ) UL[indx+1] = 360.0;
498
499 }
500
501 // Some colours and marker styles
502 mystyles.push_back( kFullCircle );
503 mystyles.push_back( kFullSquare );
504 mystyles.push_back( kFullTriangleUp );
505 mystyles.push_back( kFullTriangleDown );
506 mystyles.push_back( kFullCross );
507 mystyles.push_back( kFullStar );
508 mystyles.push_back( kFullDiamond );
509 mystyles.push_back( kFullCrossX );
510 mystyles.push_back( kFullFourTrianglesX );
511 mystyles.push_back( kFullThreeTriangles );
512 mystyles.push_back( kFullDoubleDiamond );
513 mystyles.push_back( kFourSquaresX );
514 mystyles.push_back( kFourSquaresPlus );
515 mycolors.push_back( kRed+1 );
516 mycolors.push_back( kBlue+1 );
517 mycolors.push_back( kGreen+1 );
518 mycolors.push_back( kMagenta+1 );
519 mycolors.push_back( kYellow+1 );
520 mycolors.push_back( kCyan+1 );
521
522}
523
524void MiniballAngleFitter::DoFit() {
525
526 // First we need to fit all of the segment energies
527 // or read them in from a file
528 if( flag_fit_peaks ) ff.FitSegmentEnergies( input_root_file );
529 else ff.LoadExpEnergies( input_data_filename );
530
531 // Create minimizer
532 ROOT::Math::Minimizer *min =
533 ROOT::Math::Factory::CreateMinimizer("Minuit2", "Combined");
534
535 // Set print level
536 //min->SetPrintLevel(1);
537
538 // Create function for the fitting
539 ROOT::Math::Functor f_init( ff, npars );
540
541 // Some fit controls
542 min->SetErrorDef(1.);
543 min->SetMaxFunctionCalls(1e7);
544 min->SetMaxIterations(1e8);
545 min->SetPrecision(1e-12);
546 min->SetTolerance(1e-9);
547 min->SetStrategy(2); // 0: low, 1: medium, 2: high
548 min->SetFunction(f_init);
549
550 // Set limits in fit
551 for( unsigned int i = 0; i < npars; ++i ) {
552 min->SetLimitedVariable( i, names.at(i), pars.at(i), 0.0001, LL.at(i), UL.at(i) );
553 min->SetVariableStepSize( i, 1.0 );
554 }
555
556 // Set the variable step size for beta
557 min->SetVariableStepSize( 0, 0.001 );
558
559 // If a cluster is missing, fix the parameters
560 for( unsigned int clu = 0; clu < myset->GetNumberOfMiniballClusters(); ++clu ) {
561
562 int indx = 1 + 4*clu;
563 if( !ff.IsPresent(clu) ){
564
565 min->FixVariable( indx+0 );
566 min->FixVariable( indx+1 );
567 min->FixVariable( indx+2 );
568 min->FixVariable( indx+3 );
569
570 }
571
572 // Uncomment below to fix phi
573 //min->FixVariable( indx+1 );
574
575 } // clu
576
577 // Call the minimisation procedure
578 min->Minimize();
579
580 // Print the results to the terminal
581 min->PrintResults();
582
583 // Copy results into reaction object
584 double beta = min->X()[0];
585 for( unsigned int clu = 0; clu < myset->GetNumberOfMiniballClusters(); ++clu ) {
586
587 // work out the parameter numbers
588 int indx = 1 + 4 * clu;
589
590 // Setup the cluster again
591 myreact->SetupCluster( clu, min->X()[indx], min->X()[indx+1], min->X()[indx+2], min->X()[indx+3], 0.0 );
592
593 }
594
595 // print fit + residuals to pdf
596 TGraphErrors *engraph = new TGraphErrors();
597 engraph->SetName("engraph");
598 engraph->SetTitle("Experimental energies; Channel index; Energy [keV]");
599
600 TGraph *calcgraph = new TGraph();
601 calcgraph->SetName("calcgraph");
602 calcgraph->SetTitle("Calculated energies; Channel index; Energy [keV]");
603
604 TGraphErrors *resgraph = new TGraphErrors();
605 resgraph->SetName("resgraph");
606 resgraph->SetTitle("Residuals; Channel index; Energy [keV]");
607
608 TGraphErrors *corrgraph = new TGraphErrors();
609 corrgraph->SetName("corrgraph");
610 corrgraph->SetTitle("Doppler-corrected energies; Channel index; Energy [keV]");
611
612 TGraphErrors *phigraph = new TGraphErrors();
613 phigraph->SetName("phigraph");
614 phigraph->SetTitle("Phi residuals; Channel index; Energy [keV]");
615
616
617 // Loop over clusters
618 for( unsigned int clu = 0; clu < myset->GetNumberOfMiniballClusters(); ++clu ) {
619
620 // Loop over crystals
621 for( unsigned int cry = 0; cry < myset->GetNumberOfMiniballCrystals(); ++cry ) {
622
623 // Loop over segments
624 for( unsigned int seg = 0; seg < myset->GetNumberOfMiniballSegments(); ++seg ) {
625
626 // Check if it's present
627 if( !ff.IsPresent( clu, cry, seg ) ) continue;
628
629 // Get Doppler shifted energy, theta and phi
630 double theta = myreact->GetGammaTheta( clu, cry, seg );
631 double phi = myreact->GetGammaPhi( clu, cry, seg );
632 phi *= TMath::RadToDeg();
633 if( phi < 0.0 ) phi += 360.0;
634 double edop = myreact->DopplerShift( ff.GetReferenceEnergy(),
635 beta, TMath::Cos(theta) );
636 double corr = ff.GetReferenceEnergy() / edop;
637
638 // channel index
639 double indx = seg;
640 indx += cry * myset->GetNumberOfMiniballSegments();
641 indx += clu * myset->GetNumberOfMiniballCrystals() * myset->GetNumberOfMiniballSegments();
642
643 engraph->SetPoint(engraph->GetN(), indx, ff.GetExpEnergy( clu, cry, seg ) );
644 engraph->SetPointError( engraph->GetN()-1, 0, ff.GetExpError( clu, cry, seg ) );
645
646 corrgraph->SetPoint(corrgraph->GetN(), indx, ff.GetExpEnergy( clu, cry, seg ) * corr );
647 corrgraph->SetPointError( engraph->GetN()-1, 0, ff.GetExpError( clu, cry, seg ) );
648
649 calcgraph->SetPoint(calcgraph->GetN(), indx, edop );
650
651 resgraph->SetPoint(resgraph->GetN(), indx, edop - ff.GetExpEnergy( clu, cry, seg ) );
652 resgraph->SetPointError( resgraph->GetN()-1, 0, ff.GetExpError( clu, cry, seg ) );
653
654 // Get experimental phi
655 if( !ff.HasPhiConstraint( clu, cry, seg ) ) continue;
656 phigraph->SetPoint(phigraph->GetN(), indx, phi - ff.GetExpPhi( clu, cry, seg ) );
657 phigraph->SetPointError( phigraph->GetN()-1, 0, ff.GetExpPhiError( clu, cry, seg ) );
658
659 } // seg
660
661 } // cry
662
663 } // clu
664
665 // Draw the absolute experimental/calculated energies
666 auto c1 = std::make_unique<TCanvas>();
667 engraph->SetMarkerStyle(kFullCircle);
668 engraph->SetMarkerSize(0.5);
669 engraph->SetMarkerColor(kRed);
670 engraph->SetLineColor(kRed);
671 engraph->Draw("AP");
672 //engraph->Write();
673 calcgraph->SetMarkerStyle(kFullCircle);
674 calcgraph->SetMarkerSize(0.5);
675 calcgraph->Draw("P");
676 //calcgraph->Write();
677
678 // Add a legend
679 auto leg = std::make_unique<TLegend>(0.1,0.75,0.3,0.9);
680 leg->AddEntry(engraph, "Experimental energies");
681 leg->AddEntry(calcgraph, "Calculated energies");
682 leg->Draw();
683
684 // Save first plot as a PDF
685 gErrorIgnoreLevel = kWarning;
686 c1->Print("position_cal.pdf(","pdf");
687
688 // Draw the residuals
689 resgraph->SetMarkerStyle(kFullCircle);
690 resgraph->SetMarkerSize(0.5);
691 resgraph->Draw("AP");
692 //resgraph->Write();
693 auto func = std::make_unique<TF1>("func", "[0]", 0, 168);
694 func->SetParameter( 0, 0.0 );
695 func->Draw("same");
696
697 // Save second plot as a PDF
698 c1->Print("position_cal.pdf","pdf");
699
700 // Draw the corrected energies compared to reference
701 corrgraph->SetMarkerStyle(kFullCircle);
702 corrgraph->SetMarkerSize(0.5);
703 corrgraph->Draw("AP");
704 //corrgraph->Write();
705 func->SetParameter( 0, ff.GetReferenceEnergy() );
706 func->Draw("same");
707
708 // Save third plot as a PDF
709 c1->Print("position_cal.pdf","pdf");
710
711 // Draw the phi residuals
712 phigraph->SetMarkerStyle(kFullCircle);
713 phigraph->SetMarkerSize(0.5);
714 phigraph->Draw("AP");
715 func->SetParameter( 0, 0.0 );
716 func->Draw("same");
717
718 // Save fourth plot as a PDF
719 c1->Print("position_cal.pdf","pdf");
720
721 // A multi-graph showing all the positions in theta-phi, xy, xz, etc
722 auto tp_mg = std::make_unique<TMultiGraph>();
723 auto xy_f_mg = std::make_unique<TMultiGraph>();
724 auto xy_b_mg = std::make_unique<TMultiGraph>();
725 auto xz_r_mg = std::make_unique<TMultiGraph>();
726 auto xz_l_mg = std::make_unique<TMultiGraph>();
727 tp_mg->SetName("theta_phi_map");
728 xy_f_mg->SetName("xy_f_mg");
729 xy_b_mg->SetName("xy_b_mg");
730 xz_r_mg->SetName("xz_r_mg");
731 xz_l_mg->SetName("xz_l_mg");
732 std::vector< std::vector<TGraph*> > theta_phi;
733 std::vector< std::vector<TGraph*> > xy_f;
734 std::vector< std::vector<TGraph*> > xy_b;
735 std::vector< std::vector<TGraph*> > xz_l;
736 std::vector< std::vector<TGraph*> > xz_r;
737 theta_phi.resize( myset->GetNumberOfMiniballClusters() );
738 xy_f.resize( myset->GetNumberOfMiniballClusters() );
739 xy_b.resize( myset->GetNumberOfMiniballClusters() );
740 xz_l.resize( myset->GetNumberOfMiniballClusters() );
741 xz_r.resize( myset->GetNumberOfMiniballClusters() );
742
743 // Add a legend and make space for it on the pad
744 std::string leg_lab;
745 auto leg2 = std::make_unique<TLegend>( 0.79, 0.15, 0.95, 0.90 );
746 auto p1 = c1->cd();
747 p1->SetMargin( 0.15, 0.21, 0.15, 0.10 );
748
749 // Loop over clusters
750 for( unsigned int clu = 0; clu < myset->GetNumberOfMiniballClusters(); ++clu ) {
751
752 theta_phi[clu].resize( myset->GetNumberOfMiniballCrystals() );
753 xy_f[clu].resize( myset->GetNumberOfMiniballCrystals() );
754 xy_b[clu].resize( myset->GetNumberOfMiniballCrystals() );
755 xz_l[clu].resize( myset->GetNumberOfMiniballCrystals() );
756 xz_r[clu].resize( myset->GetNumberOfMiniballCrystals() );
757
758 // Loop over crystals
759 for( unsigned int cry = 0; cry < myset->GetNumberOfMiniballCrystals(); ++cry ) {
760
761 // Make the graphs to prevent any issues
762 theta_phi[clu][cry] = new TGraph();
763 xy_f[clu][cry] = new TGraph();
764 xy_b[clu][cry] = new TGraph();
765 xz_l[clu][cry] = new TGraph();
766 xz_r[clu][cry] = new TGraph();
767
768 leg_lab = "MB" + std::to_string(clu) + static_cast<char>(cry+65);
769 leg2->AddEntry( theta_phi[clu][cry], leg_lab.data() );
770
771 // But then skip if there's no data
772 if( !ff.IsPresent( clu ) ) continue;
773
774 // Loop over segments
775 for( unsigned int seg = 0; seg < myset->GetNumberOfMiniballSegments(); ++seg ) {
776
777 // Get position
778 double theta = myreact->GetGammaTheta( clu, cry, seg ) * TMath::RadToDeg();
779 double phi = myreact->GetGammaPhi( clu, cry, seg ) * TMath::RadToDeg();
780 double x = myreact->GetGammaX( clu, cry, seg );
781 double y = myreact->GetGammaY( clu, cry, seg );
782 double z = myreact->GetGammaZ( clu, cry, seg );
783
784 // Plot positions
785 theta_phi[clu][cry]->SetPoint(theta_phi[clu][cry]->GetN(), theta, phi );
786 if( z > 0 ) xy_f[clu][cry]->SetPoint(xy_f[clu][cry]->GetN(), y, x );
787 else xy_b[clu][cry]->SetPoint(xy_b[clu][cry]->GetN(), y, x );
788 if( y > 0 ) xz_r[clu][cry]->SetPoint(xz_r[clu][cry]->GetN(), z, x );
789 else xz_l[clu][cry]->SetPoint(xz_l[clu][cry]->GetN(), z, x );
790
791 } // seg
792
793 // Set colours etc
794 theta_phi[clu][cry]->SetMarkerSize(1.0);
795 theta_phi[clu][cry]->SetMarkerStyle(mystyles[clu]);
796 theta_phi[clu][cry]->SetMarkerColor(mycolors[cry]);
797 xy_f[clu][cry]->SetMarkerSize(1.0);
798 xy_f[clu][cry]->SetMarkerStyle(mystyles[clu]);
799 xy_f[clu][cry]->SetMarkerColor(mycolors[cry]);
800 xy_b[clu][cry]->SetMarkerSize(1.0);
801 xy_b[clu][cry]->SetMarkerStyle(mystyles[clu]);
802 xy_b[clu][cry]->SetMarkerColor(mycolors[cry]);
803 xz_r[clu][cry]->SetMarkerSize(1.0);
804 xz_r[clu][cry]->SetMarkerStyle(mystyles[clu]);
805 xz_r[clu][cry]->SetMarkerColor(mycolors[cry]);
806 xz_l[clu][cry]->SetMarkerSize(1.0);
807 xz_l[clu][cry]->SetMarkerStyle(mystyles[clu]);
808 xz_l[clu][cry]->SetMarkerColor(mycolors[cry]);
809
810 // Add to multi-graph
811 if( theta_phi[clu][cry]->GetN() > 0 ) tp_mg->Add( theta_phi[clu][cry] );
812 if( xy_f[clu][cry]->GetN() > 0 ) xy_f_mg->Add( xy_f[clu][cry] );
813 if( xy_b[clu][cry]->GetN() > 0 ) xy_b_mg->Add( xy_b[clu][cry] );
814 if( xz_r[clu][cry]->GetN() > 0 ) xz_r_mg->Add( xz_r[clu][cry] );
815 if( xz_l[clu][cry]->GetN() > 0 ) xz_l_mg->Add( xz_l[clu][cry] );
816
817 } // cry
818
819 // Skip if there's no data
820 if( !ff.IsPresent( clu ) ) continue;
821
822 } // clu
823
824 // Draw the multigraph for theta-phi
825 tp_mg->GetXaxis()->SetLimits(0,180);
826 tp_mg->GetYaxis()->SetLimits(-180,180);
827 tp_mg->GetXaxis()->SetTitle("Reaction Theta [deg]");
828 tp_mg->GetYaxis()->SetTitle("Reaction Phi [deg]");
829 tp_mg->Draw("AP");
830 leg2->Draw();
831 //tp_mg->Write();
832 c1->Print("position_cal.pdf","pdf");
833
834 // Draw the multigraph for xy-forward
835 xy_f_mg->GetYaxis()->SetTitle("x [mm]");
836 xy_f_mg->GetXaxis()->SetTitle("y [mm]");
837 xy_f_mg->Draw("AP");
838 leg2->Draw();
839 //xy_f_mg->Write();
840 c1->Print("position_cal.pdf","pdf");
841
842 // Draw the multigraph for xy-backwards
843 xy_b_mg->GetYaxis()->SetTitle("x [mm]");
844 xy_b_mg->GetXaxis()->SetTitle("y [mm]");
845 xy_b_mg->Draw("AP");
846 leg2->Draw();
847 //xy_b_mg->Write();
848 c1->Print("position_cal.pdf","pdf");
849
850 // Draw the multigraph for xz-right
851 xz_r_mg->GetYaxis()->SetTitle("x [mm]");
852 xz_r_mg->GetXaxis()->SetTitle("z [mm]");
853 xz_r_mg->Draw("AP");
854 leg2->Draw();
855 //xz_r_mg->Write();
856 c1->Print("position_cal.pdf","pdf");
857
858 // Draw the multigraph for xz-left
859 xz_l_mg->GetYaxis()->SetTitle("x [mm]");
860 xz_l_mg->GetXaxis()->SetTitle("z [mm]");
861 xz_l_mg->Draw("AP");
862 leg2->Draw();
863 //xz_l_mg->Write();
864 c1->Print("position_cal.pdf)","pdf");
865 gErrorIgnoreLevel = kInfo;
866
867 // Print final results to terminal
868 std::cout << "fitted beta = " << beta << std::endl;
869 printf(" theta phi alpha R\n");
870 for( unsigned int clu = 0; clu < myset->GetNumberOfMiniballClusters(); ++clu ) {
871
872 if( !ff.IsPresent(clu) ) continue;
873 int indx = 1 + 4*clu;
874 printf("MB%i %6.2f %6.2f %6.2f %6.2f\n", clu, min->X()[indx], min->X()[indx+1], min->X()[indx+2], min->X()[indx+3]);
875
876 }
877 std::cout << std::endl << std::endl;
878
879 // Print the angles in the reaction file format
880 myreact->PrintReaction( std::cout, "a" );
881
882 return;
883
884}
885
886void MiniballAngleFitter::SaveReactionFile( std::string fname ){
887
888 // Output
889 std::ofstream react_file;
890 react_file.open( fname, std::ios::out );
891 if( react_file.is_open() ) {
892
893 myreact->PrintReaction( react_file, "ae" );
894 react_file.close();
895
896 }
897
898 else std::cerr << "Couldn't open " << fname << std::endl;
899
900}
MiniballAngleFitter::UL
std::vector< double > UL
Definition MiniballAngleFitter.hh:179
MiniballAngleFitter::mystyles
std::vector< int > mystyles
Definition MiniballAngleFitter.hh:169
MiniballAngleFitter::LL
std::vector< double > LL
Definition MiniballAngleFitter.hh:178
MiniballAngleFitter::pars
std::vector< double > pars
Definition MiniballAngleFitter.hh:176
MiniballAngleFitter::names
std::vector< std::string > names
Definition MiniballAngleFitter.hh:177
MiniballAngleFitter::mycolors
std::vector< int > mycolors
Definition MiniballAngleFitter.hh:169
MiniballAngleFitter::myset
std::shared_ptr< MiniballSettings > myset
Definition MiniballAngleFitter.hh:153
MiniballAngleFitter::SetInputEnergiesFile
bool SetInputEnergiesFile(std::string fname)
Definition MiniballAngleFitter.cc:423
MiniballAngleFitter::SetInputROOTFile
bool SetInputROOTFile(std::string fname)
Definition MiniballAngleFitter.cc:403
MiniballAngleFitter::myreact
std::shared_ptr< MiniballReaction > myreact
Definition MiniballAngleFitter.hh:154
MiniballAngleFitter::SaveReactionFile
void SaveReactionFile(std::string fname)
Definition MiniballAngleFitter.cc:886
MiniballAngleFitter::ff
MiniballAngleFunction ff
Definition MiniballAngleFitter.hh:163
MiniballAngleFunction::SaveExpEnergies
void SaveExpEnergies(std::string energy_file)
Definition MiniballAngleFitter.cc:280
MiniballAngleFunction::LoadExpEnergies
void LoadExpEnergies(std::string energy_file)
Definition MiniballAngleFitter.cc:196
MiniballAngleFunction::cluster
std::vector< bool > cluster
Definition MiniballAngleFitter.hh:101
MiniballAngleFunction::IsPresent
bool IsPresent(unsigned int clu)
Definition MiniballAngleFitter.hh:59
MiniballAngleFunction::energy
std::vector< std::vector< std::vector< double > > > energy
Definition MiniballAngleFitter.hh:97
MiniballAngleFunction::myset
std::shared_ptr< MiniballSettings > myset
Definition MiniballAngleFitter.hh:104
MiniballAngleFunction::HasPhiConstraint
bool HasPhiConstraint(unsigned int clu, unsigned int cry, unsigned int seg)
Definition MiniballAngleFitter.hh:65
MiniballAngleFunction::phic
std::vector< std::vector< std::vector< double > > > phic
Definition MiniballAngleFitter.hh:99
MiniballAngleFunction::phiconst
std::vector< std::vector< std::vector< bool > > > phiconst
Definition MiniballAngleFitter.hh:96
MiniballAngleFunction::GetExpError
double GetExpError(unsigned int clu, unsigned int cry, unsigned int seg)
Definition MiniballAngleFitter.hh:77
MiniballAngleFunction::present
std::vector< std::vector< std::vector< bool > > > present
Definition MiniballAngleFitter.hh:95
MiniballAngleFunction::FitPeak
bool FitPeak(TH1D *h, double &en, double &er)
Definition MiniballAngleFitter.cc:65
MiniballAngleFunction::FitSegmentEnergies
void FitSegmentEnergies(TFile *infile)
Definition MiniballAngleFitter.cc:117
MiniballAngleFunction::GetExpPhi
double GetExpPhi(unsigned int clu, unsigned int cry, unsigned int seg)
Definition MiniballAngleFitter.hh:82
MiniballAngleFunction::myreact
std::shared_ptr< MiniballReaction > myreact
Definition MiniballAngleFitter.hh:105
MiniballAngleFunction::phie
std::vector< std::vector< std::vector< double > > > phie
Definition MiniballAngleFitter.hh:100
MiniballAngleFunction::GetExpEnergy
double GetExpEnergy(unsigned int clu, unsigned int cry, unsigned int seg)
Definition MiniballAngleFitter.hh:74
MiniballAngleFunction::err
std::vector< std::vector< std::vector< double > > > err
Definition MiniballAngleFitter.hh:98
MiniballAngleFunction::operator()
double operator()(const double *p)
Definition MiniballAngleFitter.cc:311
MiniballAngleFunction::GetExpPhiError
double GetExpPhiError(unsigned int clu, unsigned int cry, unsigned int seg)
Definition MiniballAngleFitter.hh:85
Generated by doxygen 1.11.0