| 1 | #include <TROOT.h>
|
|---|
| 2 | #include <TCanvas.h>
|
|---|
| 3 | #include <TProfile.h>
|
|---|
| 4 | #include <TTimer.h>
|
|---|
| 5 | #include <TH1F.h>
|
|---|
| 6 | #include <TH2F.h>
|
|---|
| 7 | #include <Getline.h>
|
|---|
| 8 | #include <TLine.h>
|
|---|
| 9 | #include <TMath.h>
|
|---|
| 10 |
|
|---|
| 11 | #include <stdint.h>
|
|---|
| 12 | #include <cstdio>
|
|---|
| 13 |
|
|---|
| 14 | #define HAVE_ZLIB
|
|---|
| 15 | #include "fits.h"
|
|---|
| 16 |
|
|---|
| 17 | #include "openFits.h"
|
|---|
| 18 | #include "openFits.c"
|
|---|
| 19 |
|
|---|
| 20 | #include "DrsCalibration.h"
|
|---|
| 21 | #include "DrsCalibration.C"
|
|---|
| 22 |
|
|---|
| 23 | #include "SpikeRemoval.h"
|
|---|
| 24 | #include "SpikeRemoval.C"
|
|---|
| 25 |
|
|---|
| 26 | #define NPIX 1440
|
|---|
| 27 | #define NCELL 1024
|
|---|
| 28 |
|
|---|
| 29 | // data access and treatment
|
|---|
| 30 | #define FAD_MAX_SAMPLES 1024
|
|---|
| 31 |
|
|---|
| 32 | int NEvents;
|
|---|
| 33 | vector<int16_t> Data; // vector, which will be filled with raw data
|
|---|
| 34 | vector<int16_t> StartCells; // vector, which will be filled with DRS start positions
|
|---|
| 35 | unsigned int EventID; // index of the current event
|
|---|
| 36 | UInt_t RegionOfInterest; // Width of the Region, read out of the DRS
|
|---|
| 37 | UInt_t NumberOfPixels; // Total number of pixel, read out of the camera
|
|---|
| 38 | size_t PXLxROI; // Size of column "Data" = #Pixel x ROI
|
|---|
| 39 |
|
|---|
| 40 | int NBSLeve = 1000;
|
|---|
| 41 |
|
|---|
| 42 | size_t TriggerOffsetROI, RC;
|
|---|
| 43 | vector<float> Offset, Gain, TriggerOffset;
|
|---|
| 44 |
|
|---|
| 45 |
|
|---|
| 46 | vector<float> Ameas(FAD_MAX_SAMPLES); // copy of the data (measured amplitude
|
|---|
| 47 | vector<float> N1mean(FAD_MAX_SAMPLES); // mean of the +1 -1 ch neighbors
|
|---|
| 48 | vector<float> N2mean(FAD_MAX_SAMPLES); // mean of the +2 -2 ch neighbors
|
|---|
| 49 | vector<float> SpikeEst(FAD_MAX_SAMPLES); // corrected Values
|
|---|
| 50 | vector<float> Vcorr(FAD_MAX_SAMPLES); // corrected Values
|
|---|
| 51 | vector<float> Vdiff(FAD_MAX_SAMPLES); // numerical derivative
|
|---|
| 52 |
|
|---|
| 53 | vector<float> Vslide(FAD_MAX_SAMPLES); // sliding average result
|
|---|
| 54 | vector<float> Vcfd(FAD_MAX_SAMPLES); // CDF result
|
|---|
| 55 | vector<float> Vcfd2(FAD_MAX_SAMPLES); // CDF result + 2nd sliding average
|
|---|
| 56 |
|
|---|
| 57 | #include "factfir.C"
|
|---|
| 58 |
|
|---|
| 59 | void computeSpikeEstimator( int );
|
|---|
| 60 |
|
|---|
| 61 | // histograms
|
|---|
| 62 | const int Ntypes = 7;
|
|---|
| 63 | const unsigned int // arranged by Dominik
|
|---|
| 64 | tAmeas = 0,
|
|---|
| 65 | tSpikeEst = 1,
|
|---|
| 66 | tVcorr = 2,
|
|---|
| 67 | tVtest = 3,
|
|---|
| 68 | tVslide = 4,
|
|---|
| 69 | tVcfd = 5,
|
|---|
| 70 | tVcfd2 = 6;
|
|---|
| 71 |
|
|---|
| 72 | TH1F* h;
|
|---|
| 73 | //TH2F* hStartCell; // id of the DRS physical pipeline cell where readout starts
|
|---|
| 74 | // x = pixel id, y = DRS cell id
|
|---|
| 75 | TH2F hPixelCellData("PixelPedestal", "PixelPedestal", NCELL, 0., NCELL, 200, -50., 150.);
|
|---|
| 76 |
|
|---|
| 77 | void BookHistos( int );
|
|---|
| 78 |
|
|---|
| 79 | // Create a canvas
|
|---|
| 80 | TCanvas* CW;
|
|---|
| 81 | TCanvas* cFilter;
|
|---|
| 82 |
|
|---|
| 83 | int spikeDebug = 1;
|
|---|
| 84 |
|
|---|
| 85 |
|
|---|
| 86 | int fana2(
|
|---|
| 87 | char *datafilename = "../raw/20110916_025.fits",
|
|---|
| 88 | const char *drsfilename = "../raw/20110916_024.drs.fits",
|
|---|
| 89 | int pixelnr = 0,
|
|---|
| 90 | int firstevent = 0,
|
|---|
| 91 | int nevents = -1 ){
|
|---|
| 92 | // read and analyze FACT raw data
|
|---|
| 93 |
|
|---|
| 94 | // sliding window filter settings
|
|---|
| 95 | int k_slide = 16;
|
|---|
| 96 | vector<double> a_slide(k_slide, 1);
|
|---|
| 97 | double b_slide = k_slide;
|
|---|
| 98 |
|
|---|
| 99 | // CFD filter settings
|
|---|
| 100 | int k_cfd = 10;
|
|---|
| 101 | vector<double> a_cfd(k_cfd, 0);
|
|---|
| 102 | double b_cfd = 1.;
|
|---|
| 103 | a_cfd[0]=0.75;
|
|---|
| 104 | a_cfd[k_cfd-1]=-1.;
|
|---|
| 105 |
|
|---|
| 106 | // 2nd slinding window filter
|
|---|
| 107 | //int ks2 = 16;
|
|---|
| 108 | //vector<double> as2(ks2, 1);
|
|---|
| 109 | //double bs2 = ks2;
|
|---|
| 110 | gROOT->SetStyle("Plain");
|
|---|
| 111 |
|
|---|
| 112 | //-------------------------------------------
|
|---|
| 113 | // Open the file
|
|---|
| 114 | //-------------------------------------------
|
|---|
| 115 | fits * datafile = NULL;
|
|---|
| 116 | NEvents = openDataFits(
|
|---|
| 117 | datafilename,
|
|---|
| 118 | &datafile,
|
|---|
| 119 | Data,
|
|---|
| 120 | StartCells,
|
|---|
| 121 | EventID,
|
|---|
| 122 | RegionOfInterest,
|
|---|
| 123 | NumberOfPixels,
|
|---|
| 124 | PXLxROI );
|
|---|
| 125 | -
|
|---|
| 126 | printf( "number of events in file: %d\n", NEvents );
|
|---|
| 127 | if (NEvents == 0){
|
|---|
| 128 | cout << "return code of openDataFits:" << datafilename<< endl;
|
|---|
| 129 | cout << "is zero -> aborting." << endl;
|
|---|
| 130 | return 1;
|
|---|
| 131 | }
|
|---|
| 132 |
|
|---|
| 133 | // compare the number of events in the data file with the nevents the
|
|---|
| 134 | // the user would like to read. nevents = -1 means: read all
|
|---|
| 135 | if ( nevents == -1 || nevents > NEvents ) nevents = NEvents;
|
|---|
| 136 |
|
|---|
| 137 |
|
|---|
| 138 | //-------------------------------------------
|
|---|
| 139 | //Get the DRS calibration
|
|---|
| 140 | //-------------------------------------------
|
|---|
| 141 |
|
|---|
| 142 | RC = openCalibFits( drsfilename, Offset, Gain, TriggerOffset, TriggerOffsetROI, 3);
|
|---|
| 143 | if (RC == 0){
|
|---|
| 144 | cout << "return code of openCalibFits:" << drsfilename << endl;
|
|---|
| 145 | cout << "is zero -> aborting." << endl;
|
|---|
| 146 | return 1;
|
|---|
| 147 | }
|
|---|
| 148 |
|
|---|
| 149 | // Book the histograms
|
|---|
| 150 | BookHistos( RegionOfInterest );
|
|---|
| 151 |
|
|---|
| 152 | // Loop over events
|
|---|
| 153 | cout << "--------------------- Data --------------------" << endl;
|
|---|
| 154 |
|
|---|
| 155 | // float value;
|
|---|
| 156 |
|
|---|
| 157 | // TH1F * sp = new TH1F("spektrum", "test of Stepktrum", 256, -0.5, 63.5);
|
|---|
| 158 |
|
|---|
| 159 | size_t calib_RC = 1;
|
|---|
| 160 | for ( int ev = firstevent; ev < firstevent + nevents; ev++) {
|
|---|
| 161 |
|
|---|
| 162 | datafile->GetRow( ev );
|
|---|
| 163 | if (ev % 50 ==0){
|
|---|
| 164 | cout << "Event number: " << EventID << endl;
|
|---|
| 165 | }
|
|---|
| 166 |
|
|---|
| 167 | // get the data of this pixel from the Data vector
|
|---|
| 168 | // apply the Drs Calibration and cut off 12 slices at the beginning
|
|---|
| 169 | // and at the end.
|
|---|
| 170 | calib_RC = applyDrsCalibration( Ameas, pixelnr,12,12, Offset, Gain, TriggerOffset,
|
|---|
| 171 | RegionOfInterest, Data, StartCells);
|
|---|
| 172 | if (calib_RC == 0){
|
|---|
| 173 | break;
|
|---|
| 174 | }
|
|---|
| 175 |
|
|---|
| 176 | computeSpikeEstimator( RegionOfInterest );
|
|---|
| 177 | removeSpikes( Ameas, Vcorr );
|
|---|
| 178 | sliding_avg( Vcorr, Vslide, 8 );
|
|---|
| 179 |
|
|---|
| 180 | // for ( unsigned int sl = 0; sl < RegionOfInterest; sl++){
|
|---|
| 181 | // hPixelCellData.Fill(sl, Vcorr[ sl ]);
|
|---|
| 182 | //}
|
|---|
| 183 |
|
|---|
| 184 | // filter Vcorr with sliding average using FIR filter function
|
|---|
| 185 | factfir(b_slide , a_slide, k_slide, Vcorr, Vslide);
|
|---|
| 186 |
|
|---|
| 187 | // filter Vslide with CFD using FIR filter function
|
|---|
| 188 | factfir(b_cfd , a_cfd, k_cfd, Vslide, Vcfd);
|
|---|
| 189 | // filter Vcfd with sliding average using FIR filter function
|
|---|
| 190 | factfir(b_slide , a_slide, k_slide, Vcfd, Vcfd2);
|
|---|
| 191 |
|
|---|
| 192 | if ( spikeDebug ){
|
|---|
| 193 | for ( unsigned int sl = 0; sl < RegionOfInterest; sl++ ){
|
|---|
| 194 |
|
|---|
| 195 | h[tAmeas].SetBinContent( sl, Ameas[ sl ] );
|
|---|
| 196 | h[tVcorr].SetBinContent( sl, Vcorr[ sl ] );
|
|---|
| 197 |
|
|---|
| 198 | h[tVslide].SetBinContent( sl, Vslide[ sl ] );
|
|---|
| 199 | h[tVcfd].SetBinContent( sl, Vcfd[ sl ] );
|
|---|
| 200 | h[tVcfd2].SetBinContent( sl, Vcfd2[ sl ] );
|
|---|
| 201 | }
|
|---|
| 202 | }
|
|---|
| 203 |
|
|---|
| 204 | if ( spikeDebug ){
|
|---|
| 205 | CW->cd( tAmeas + 1);
|
|---|
| 206 | gPad->SetGrid();
|
|---|
| 207 | h[tAmeas].Draw();
|
|---|
| 208 |
|
|---|
| 209 | CW->cd( tSpikeEst + 1);
|
|---|
| 210 | gPad->SetGrid();
|
|---|
| 211 | h[tSpikeEst].Draw();
|
|---|
| 212 |
|
|---|
| 213 | CW->cd( tVcorr + 1);
|
|---|
| 214 | gPad->SetGrid();
|
|---|
| 215 | h[tVcorr].Draw();
|
|---|
| 216 |
|
|---|
| 217 | cFilter->cd( Ntypes - tVslide );
|
|---|
| 218 | cFilter->cd(1);
|
|---|
| 219 | gPad->SetGrid();
|
|---|
| 220 | h[tVslide].Draw();
|
|---|
| 221 | cFilter->cd( Ntypes - tVcfd );
|
|---|
| 222 | cFilter->cd(2);
|
|---|
| 223 | gPad->SetGrid();
|
|---|
| 224 | h[tVcfd].Draw();
|
|---|
| 225 | TLine zeroline(0, 0, 1024, 0);
|
|---|
| 226 | zeroline.SetLineColor(kBlue);
|
|---|
| 227 | zeroline.Draw();
|
|---|
| 228 |
|
|---|
| 229 | cFilter->cd( Ntypes - tVcfd2 );
|
|---|
| 230 | cFilter->cd(3);
|
|---|
| 231 | gPad->SetGrid();
|
|---|
| 232 | h[tVcfd2].Draw();
|
|---|
| 233 |
|
|---|
| 234 | zeroline.Draw();
|
|---|
| 235 |
|
|---|
| 236 | CW->Update();
|
|---|
| 237 | cFilter->Update();
|
|---|
| 238 |
|
|---|
| 239 | //Process gui events asynchronously during input
|
|---|
| 240 | TTimer timer("gSystem->ProcessEvents();", 50, kFALSE);
|
|---|
| 241 | timer.TurnOn();
|
|---|
| 242 | TString input = Getline("Type 'q' to exit, <return> to go on: ");
|
|---|
| 243 | timer.TurnOff();
|
|---|
| 244 | if (input=="q\n") break;
|
|---|
| 245 | }
|
|---|
| 246 |
|
|---|
| 247 |
|
|---|
| 248 | }
|
|---|
| 249 |
|
|---|
| 250 | return( 0 );
|
|---|
| 251 | }
|
|---|
| 252 |
|
|---|
| 253 | void computeSpikeEstimator( int Samples ){
|
|---|
| 254 |
|
|---|
| 255 | // compute the mean of the left and right neighbors of a channel
|
|---|
| 256 |
|
|---|
| 257 | for( int i = 1; i < Samples-1; i++){
|
|---|
| 258 |
|
|---|
| 259 | N1mean[ i ] = ( Ameas[i-1] + Ameas[i+1] ) / 2.;
|
|---|
| 260 | SpikeEst[ i ] = Ameas[ i ] - N1mean[ i ];
|
|---|
| 261 | h[tSpikeEst].SetBinContent( i, SpikeEst[ i ] );
|
|---|
| 262 | }
|
|---|
| 263 | } // end of computeN1mean computation
|
|---|
| 264 |
|
|---|
| 265 |
|
|---|
| 266 | void BookHistos( int Samples ){
|
|---|
| 267 | // booking and parameter settings for all histos
|
|---|
| 268 |
|
|---|
| 269 | h = new TH1F[ Ntypes ];
|
|---|
| 270 |
|
|---|
| 271 | for ( int type = 0; type < Ntypes; type++){
|
|---|
| 272 |
|
|---|
| 273 | h[ type ].SetBins(Samples, 0, Samples);
|
|---|
| 274 | h[ type ].SetLineColor(1);
|
|---|
| 275 | h[ type ].SetLineWidth(2);
|
|---|
| 276 |
|
|---|
| 277 | // set X axis paras
|
|---|
| 278 | h[ type ].GetXaxis()->SetLabelSize(0.1);
|
|---|
| 279 | h[ type ].GetXaxis()->SetTitleSize(0.1);
|
|---|
| 280 | h[ type ].GetXaxis()->SetTitleOffset(1.2);
|
|---|
| 281 | h[ type ].GetXaxis()->SetTitle(Form("Time slice (%.1f ns/slice)", 1./2.));
|
|---|
| 282 |
|
|---|
| 283 | // set Y axis paras
|
|---|
| 284 | h[ type ].GetYaxis()->SetLabelSize(0.1);
|
|---|
| 285 | h[ type ].GetYaxis()->SetTitleSize(0.1);
|
|---|
| 286 | h[ type ].GetYaxis()->SetTitleOffset(0.3);
|
|---|
| 287 | h[ type ].GetYaxis()->SetTitle("Amplitude (a.u.)");
|
|---|
| 288 | }
|
|---|
| 289 | CW = new TCanvas("CW","DRS Waveform",10,10,800,600);
|
|---|
| 290 | CW->Divide(1, 3);
|
|---|
| 291 | cFilter = new TCanvas("cFilter","filtered DRS Waveforms",10,10,800,600);
|
|---|
| 292 | cFilter->Divide(1, 3);
|
|---|
| 293 |
|
|---|
| 294 | // hStartCell = new TH2F("StartCell", "StartCell", 1440, 0., 1440., 1024, 0., 1024);
|
|---|
| 295 |
|
|---|
| 296 | }
|
|---|
