| 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 |
|
|---|
| 20 | int 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 | }
|
|---|
