source: fact/tools/rootmacros/poissonanalysis.C@ 18481

Last change on this file since 18481 was 12204, checked in by neise, 13 years ago
removed CRLF
File size: 8.6 KB
Line 
1#include <TROOT.h>
2#include <TCanvas.h>
3#include <TProfile.h>
4#include <TH1.h>
5#include <TH2.h>
6#include <TF1.h>
7#include <TGraph.h>
8
9#include <stdint.h>
10#include <cstdio>
11
12#include "fits.h"
13#include "FOpenDataFile.c"
14#include "FOpenCalibFile.c"
15#include "FPedestal.c"
16#include "FIntFixedPosAllPx.c"
17
18#define n_peaks 2
19
20int poissonanalysis(const char *name, const char *drsname, UInt_t integration_delay, UInt_t integration_size, UInt_t spect_max )
21{
22//******************************************************************************
23//Read a datafile and plot the DRS-calibrated data
24//ATTENTION: only works for ROI=1024
25// (array indices of the calibration wrong otherwise)
26//******************************************************************************
27
28 gROOT->SetStyle("Plain");
29// TFile f("spectra.root","RECREATE");
30//-------------------------------------------
31//Define the data variables
32//-------------------------------------------
33 vector<int16_t> data;
34 vector<int16_t> data_offset;
35 unsigned int data_num;
36 size_t data_n;
37 UInt_t data_px;
38 UInt_t data_roi;
39
40//-------------------------------------------
41//Get the DRS calibration
42//-------------------------------------------
43 size_t drs_n;
44 vector<float> drs_basemean;
45 vector<float> drs_gainmean;
46 vector<float> drs_triggeroffsetmean;
47 FOpenCalibFile(drsname, drs_basemean, drs_gainmean, drs_triggeroffsetmean, drs_n);
48
49//-------------------------------------------
50//Check the sizes of the data columns
51//-------------------------------------------
52// if(drs_n!=data_n)
53// {
54// cout << "Data and DRS file incompatible (Px*ROI disagree)" << endl;
55// return 1;
56// }
57
58//-------------------------------------------
59//Open the file
60//-------------------------------------------
61 fits datafile(name);
62 if (!datafile)
63 {
64 cout << "Couldn't properly open the datafile." << endl;
65 return 1;
66 }
67 FOpenDataFile(datafile, data, data_offset, data_num, data_n, data_roi, data_px);
68
69//-------------------------------------------
70//Create the canvas, plot and title
71//-------------------------------------------
72 char title_spect[500];
73 char name_spect[50];
74// std::sprintf(title_spect,"Spectrum of Px %i (data: %s, DRS: %s)",pixelnr, name,drsname);
75 TCanvas *canv_spect = new TCanvas( "canv_spect", "Spectrum", 950, 10, 700, 500 );
76 TH1F* spectrum[data_px];
77 for(int i=0; i<data_px; i++) {
78 std::sprintf(title_spect,"Spectrum of Px %i, %i-%i (data: %s, DRS: %s)",i,integration_delay,integration_delay+integration_size,name,drsname);
79 std::sprintf(name_spect,"spect%i",i);
80// spectrum[i] = new TH1F(name_spect,title_spect,400,-150,650);
81 spectrum[i] = new TH1F(name_spect,title_spect,400,-150,spect_max);
82 spectrum[i]->GetXaxis()->SetTitle("Pulse size (a.u.)");
83 spectrum[i]->GetYaxis()->SetTitle("Counts");
84 }
85
86// TCanvas *canv_gain= new TCanvas( "canv_gain", "Gain in arbitrary units", 50, 500, 1800, 800 );
87//-------------------------------------------
88//Define the fit
89//-------------------------------------------
90// int fit_start = 20;
91// int fit_end = 200;
92// TF1 *f_peaks = new TF1("total","[2]*TMath::Exp(-(x-[0])*(x-[0])/(2*[3]*[3])) + [4]*TMath::Exp(-(x-[0]-[1])*(x-[0]-[1])/(2*[5]*[5])) + [6]*TMath::Exp((-x-[7])/[8])",fit_start,fit_end);
93// Double_t par_init[5+2*n_peaks] = {70,90,100,15,10,15,10,30,30};
94// Double_t parread_tmp[5+2*n_peaks];
95// Double_t parread[data_px][5+2*n_peaks];
96//// TF1 *f_peaks = new TF1("total","[2]*TMath::Exp(-(x-[0])*(x-[0])/(2*[3]*[3])) + [4]*TMath::Exp(-(x-[0]-[1])*(x-[0]-[1])/(2*[5]*[5])) + [6]*TMath::Exp(-(x-[0]-2*[1])*(x-[0]-2*[1])/(2*[7]*[7])) + [8]*TMath::Exp(-(x-[0]-3*[1])*(x-[0]-3*[1])/(2*[9]*[9])) + [10]*TMath::Exp(-(x-[12])*(x-[12])/(2*[11]*[11]))",fit_start,fit_end);
97//// Double_t par_init[3+2*n_peaks] = {70,90,500000,15,50000,15,5000,15,500,15,800000,20,0};
98// f_peaks->SetParameters(par_init);
99// f_peaks->SetLineColor(2);
100//
101// TF1 *f_pedestal = new TF1("pedestal","[0]*TMath::Exp((-x-[1])/[2])",fit_start,fit_end);
102// f_pedestal->SetLineColor(2);
103// f_pedestal->SetLineWidth(1);
104//
105// TF1* f_peak[n_peaks];
106// for(int i=0; i<n_peaks; i++) {
107// f_peak[i] = new TF1("peak","gaus",fit_start,fit_end);
108// f_peak[i]->SetLineColor(2);
109// f_peak[i]->SetLineWidth(1);
110// }
111
112//-------------------------------------------
113//Find the baseline
114//-------------------------------------------
115// float pedestal_mean, pedestal_rms;
116//// FPedestal(title_base, datafile, data, data_offset, drs_basemean, drs_gainmean, drs_triggeroffsetmean, data_roi, pixelnr, pedestal_mean, pedestal_rms);
117// //Return values
118// //Value of maximal probability: -2.225 +- 4.3144
119// pedestal_mean = -2.225;
120// pedestal_rms = 4.3144;
121
122//-------------------------------------------
123//Oscilloscope
124//-------------------------------------------
125// float threshold = pedestal_mean+pedestal_rms;
126// float threshold = 2.5;
127
128//***********************
129// UInt_t integration_size = 5;
130// UInt_t integration_delay = 238;
131
132 FIntFixedPosAllPx(datafile, data, data_offset, data_num, data_px, drs_basemean, drs_gainmean, drs_triggeroffsetmean, data_roi, integration_size, integration_delay, spectrum);
133
134 canv_spect->cd();
135 spectrum[0]->Draw();
136 canv_spect->Modified();
137 canv_spect->Update();
138
139// Double_t parread_dist[data_px];
140// Double_t parread_x[data_px];
141// for(int i=0; i<data_px; i++) {
142// f_peaks->SetParameters(par_init);
143// if(!(spectrum[i]->Fit(f_peaks,"NR"))) {
144//// std::cout << "Fit successful." << std::endl;
145// f_peaks->GetParameters(&parread_tmp[0]);
146// for(int j=0; j<5+2*n_peaks; j++) {
147// parread[i][j] = parread_tmp[j];
148// }
149// parread_dist[i] = parread_tmp[1];
150// parread_x[i] = i;
151// }
152// else {
153// for(int j=0; j<5+2*n_peaks; j++) {
154// parread[i][j] = 0;
155// }
156// parread_dist[i] = 0;
157// parread_x[i] = i;
158// }
159// }
160//
161// canv_gain->cd();
162// TGraph* plot_gain = new TGraph(data_px,parread_x,parread_dist);
163// plot_gain->GetYaxis()->SetRangeUser(0,150);
164// plot_gain->Draw("A*");
165// canv_gain->Modified();
166// canv_gain->Update();
167//
168// canv_spect->cd();
169// canv_spect->SetLogy();
170
171//-------------------------------------------
172//Draw the data where the gain is out of limits
173//-------------------------------------------
174// char temp = 'x';
175// std::cout << "Plot the spectra out of limits: Enter for next, 'a' to abort." << std::endl;
176// for(int i=0; i<data_px; i++) {
177// if((parread_dist[i]<80)||(parread_dist[i]>105)) {
178// std::cout << "Gain out of limits: Px " << i << " at " << parread_dist[i] << std::endl;
179// spectrum[i]->Draw();
180//
181// f_peaks->SetParameters(&parread[i][0]);
182// f_peak[0]->SetParameters(parread[i][2],parread[i][0],parread[i][3]);
183// f_peak[1]->SetParameters(parread[i][4],parread[i][0]+parread[i][1],parread[i][5]);
184//// f_peak[2]->SetParameters(parread[i][6],parread[i][0]+2*parread[i][1],parread[i][7]);
185//// f_peak[3]->SetParameters(parread[i][8],parread[i][0]+3*parread[i][1],parread[i][9]);
186//// f_peak[4]->SetParameters(parread[i][10],parread[i][12],parread[i][11]);
187// f_pedestal->SetParameters(parread[i][6],parread[i][7],parread[i][8]);
188//
189// f_peaks->Draw("same");
190// for(int i=0; i<n_peaks; i++) {
191// f_peak[i]->Draw("same");
192// }
193// f_pedestal->Draw("same");
194//
195// canv_spect->Modified();
196// canv_spect->Update();
197// temp='x';
198// while ((temp!='\n') && (temp!='a'))
199// cin.get(temp);
200// if(temp=='a') break;
201// }
202// }
203//-------------------------------------------
204//Draw the data
205//-------------------------------------------
206 canv_spect->cd();
207 Int_t plotspect_nr;
208 std::cout << "Enter the number of the spectrum to plot (-1 to quit): ";
209 std::cin >> plotspect_nr;
210 while(plotspect_nr>=0) {
211 plotspect_nr = plotspect_nr % data_px;
212 spectrum[plotspect_nr]->Draw();
213
214// f_peaks->SetParameters(&parread[plotspect_nr][0]);
215// f_peak[0]->SetParameters(parread[plotspect_nr][2],parread[plotspect_nr][0],parread[plotspect_nr][3]);
216// f_peak[1]->SetParameters(parread[plotspect_nr][4],parread[plotspect_nr][0]+parread[plotspect_nr][1],parread[plotspect_nr][5]);
217//// f_peak[2]->SetParameters(parread[plotspect_nr][6],parread[plotspect_nr][0]+2*parread[plotspect_nr][1],parread[plotspect_nr][7]);
218//// f_peak[3]->SetParameters(parread[plotspect_nr][8],parread[plotspect_nr][0]+3*parread[plotspect_nr][1],parread[plotspect_nr][9]);
219//// f_peak[4]->SetParameters(parread[plotspect_nr][10],parread[plotspect_nr][12],parread[plotspect_nr][11]);
220// f_pedestal->SetParameters(parread[plotspect_nr][6],parread[plotspect_nr][7],parread[plotspect_nr][8]);
221//
222// f_peaks->Draw("same");
223// for(int i=0; i<n_peaks; i++) {
224// f_peak[i]->Draw("same");
225// }
226// f_pedestal->Draw("same");
227
228 canv_spect->Modified();
229 canv_spect->Update();
230
231 std::cout << "Enter the number of the spectrum to plot (-1 to quit): ";
232 std::cin >> plotspect_nr;
233 }
234
235// f->Write();
236 return 0;
237}
Note: See TracBrowser for help on using the repository browser.