| 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 <TBox.h>
 | 
|---|
| 10 | #include <TMath.h>
 | 
|---|
| 11 | #include <TFile.h>
 | 
|---|
| 12 | #include <TStyle.h>
 | 
|---|
| 13 | 
 | 
|---|
| 14 | #include <stdio.h>
 | 
|---|
| 15 | #include <stdint.h>
 | 
|---|
| 16 | #include <cstdio>
 | 
|---|
| 17 | 
 | 
|---|
| 18 | #define NPIX  1440
 | 
|---|
| 19 | #define NCELL 1024
 | 
|---|
| 20 | #define FAD_MAX_SAMPLES 1024
 | 
|---|
| 21 | 
 | 
|---|
| 22 | #define HAVE_ZLIB
 | 
|---|
| 23 | #include "fits.h"
 | 
|---|
| 24 | #include "FOpenCalibFile.c"
 | 
|---|
| 25 | 
 | 
|---|
| 26 | #include "discriminator.h"
 | 
|---|
| 27 | #include "discriminator.C"
 | 
|---|
| 28 | #include "zerosearch.h"
 | 
|---|
| 29 | #include "zerosearch.C"
 | 
|---|
| 30 | #include "factfir.C"
 | 
|---|
| 31 | 
 | 
|---|
| 32 | 
 | 
|---|
| 33 | vector<int16_t> AllPixelDataVector;
 | 
|---|
| 34 | vector<int16_t> StartCellVector;
 | 
|---|
| 35 | 
 | 
|---|
| 36 | unsigned int CurrentEventID;
 | 
|---|
| 37 | 
 | 
|---|
| 38 | bool breakout=false;
 | 
|---|
| 39 | 
 | 
|---|
| 40 | size_t ROIxNOP;
 | 
|---|
| 41 | UInt_t NumberOfPixels;
 | 
|---|
| 42 | UInt_t RegionOfInterest;
 | 
|---|
| 43 | int NEvents;
 | 
|---|
| 44 | 
 | 
|---|
| 45 | size_t drs_n;
 | 
|---|
| 46 | vector<float> drs_basemean;
 | 
|---|
| 47 | vector<float> drs_gainmean;
 | 
|---|
| 48 | vector<float> drs_triggeroffsetmean;
 | 
|---|
| 49 | 
 | 
|---|
| 50 | int FOpenDataFile( fits & );
 | 
|---|
| 51 | 
 | 
|---|
| 52 | 
 | 
|---|
| 53 | vector<float> Ameas(FAD_MAX_SAMPLES);  // copy of the data (measured amplitude
 | 
|---|
| 54 | vector<float> N1mean(FAD_MAX_SAMPLES); // mean of the +1 -1 ch neighbors
 | 
|---|
| 55 | vector<float> Vcorr(FAD_MAX_SAMPLES);  // corrected Values
 | 
|---|
| 56 | vector<float> Vdiff(FAD_MAX_SAMPLES);  // numerical derivative
 | 
|---|
| 57 | 
 | 
|---|
| 58 | vector<float> Vslide(FAD_MAX_SAMPLES);  // sliding average result
 | 
|---|
| 59 | vector<float> Vcfd(FAD_MAX_SAMPLES);    // CDF result
 | 
|---|
| 60 | vector<float> Vcfd2(FAD_MAX_SAMPLES);    // CDF result + 2nd sliding average
 | 
|---|
| 61 | 
 | 
|---|
| 62 | // not needed?
 | 
|---|
| 63 | //vector<float> N2mean(FAD_MAX_SAMPLES); // mean of the +2 -2 ch neighbors
 | 
|---|
| 64 | 
 | 
|---|
| 65 | 
 | 
|---|
| 66 | float getValue( int, int );
 | 
|---|
| 67 | float correctDrsOffset( int slice, int pixel );
 | 
|---|
| 68 | void computeN1mean( int );
 | 
|---|
| 69 | void removeSpikes( int );
 | 
|---|
| 70 | 
 | 
|---|
| 71 | // histograms
 | 
|---|
| 72 | const int Ntypes = 7;
 | 
|---|
| 73 | const unsigned int  // arranged by Dominik
 | 
|---|
| 74 |         tAmeas  = 0, 
 | 
|---|
| 75 |         tN1mean = 1, 
 | 
|---|
| 76 |         tVcorr  = 2,
 | 
|---|
| 77 |         tVtest  = 3,
 | 
|---|
| 78 |         tVslide = 4,
 | 
|---|
| 79 |         tVcfd   = 5,
 | 
|---|
| 80 |         tVcfd2  = 6;
 | 
|---|
| 81 | 
 | 
|---|
| 82 | TH1F* h;
 | 
|---|
| 83 | TH1F *oldh;
 | 
|---|
| 84 | TH2F* hStartCell;   // id of the DRS physical pipeline cell where readout starts
 | 
|---|
| 85 |                     // x = pixel id, y = DRS cell id
 | 
|---|
| 86 | TH2F hPixelCellData(
 | 
|---|
| 87 |                 "PixelPedestal", "PixelPedestal", NCELL, 0., NCELL, 200, -50., 150.);
 | 
|---|
| 88 | TH1F *hBaseline[ NPIX ]; // histograms for baseline extraction
 | 
|---|
| 89 | TH1F *hMeanBsl, *hpltMeanBsl;
 | 
|---|
| 90 | TH1F *hRmsBsl, *hpltRmsBsl; 
 | 
|---|
| 91 | TH2F * hAmplSpek_cfd;
 | 
|---|
| 92 | TH2F * hAmplSpek_discri;
 | 
|---|
| 93 | TObjArray hList;
 | 
|---|
| 94 | TObjArray hListBaseline;
 | 
|---|
| 95 | 
 | 
|---|
| 96 | void BookHistos( );
 | 
|---|
| 97 | void SaveHistograms( char * );
 | 
|---|
| 98 | 
 | 
|---|
| 99 | // Create a canvas
 | 
|---|
| 100 | TCanvas* CW;
 | 
|---|
| 101 | TCanvas* cFilter;
 | 
|---|
| 102 | 
 | 
|---|
| 103 | 
 | 
|---|
| 104 | int searchSinglesPeaks = 0;
 | 
|---|
| 105 | 
 | 
|---|
| 106 | int fpeak_discri( 
 | 
|---|
| 107 |         char *datafilename              = "../../20111011_055.fits.gz",
 | 
|---|
| 108 |         const char *drsfilename = "../../20111011_054.drs.fits.gz",
 | 
|---|
| 109 |         int PixelID                             = -1,
 | 
|---|
| 110 |         int firstevent                  = 0, 
 | 
|---|
| 111 |         int nevents                     = -1,
 | 
|---|
| 112 |         bool spikeDebug = false,
 | 
|---|
| 113 |         int verbosityLevel = 1 // different verbosity levels can be implemented here
 | 
|---|
| 114 |         )
 | 
|---|
| 115 | {
 | 
|---|
| 116 | 
 | 
|---|
| 117 | gStyle->SetPalette(1,0);
 | 
|---|
| 118 | gROOT->SetStyle("Plain");
 | 
|---|
| 119 | // read FACT raw data 
 | 
|---|
| 120 | //      * remove spikes
 | 
|---|
| 121 | //      * calculate baseline
 | 
|---|
| 122 | //      * find peaks (CFD and normal discriminator)
 | 
|---|
| 123 | //      * compute pulse height and pulse integral spektrum of the peaks
 | 
|---|
| 124 | 
 | 
|---|
| 125 | // sliding window filter settings
 | 
|---|
| 126 |         int k_slide = 16;
 | 
|---|
| 127 |         vector<double> a_slide(k_slide, 1);
 | 
|---|
| 128 |         double b_slide = k_slide;
 | 
|---|
| 129 | 
 | 
|---|
| 130 | // CFD filter settings
 | 
|---|
| 131 |         int k_cfd = 10;
 | 
|---|
| 132 |         vector<double> a_cfd(k_cfd, 0);
 | 
|---|
| 133 |         double b_cfd = 1.;
 | 
|---|
| 134 |         a_cfd[0]=-0.75;
 | 
|---|
| 135 |         a_cfd[k_cfd-1]=1.;
 | 
|---|
| 136 | 
 | 
|---|
| 137 | // 2nd slinding window filter
 | 
|---|
| 138 |         int ks2 = 16;
 | 
|---|
| 139 |         vector<double> as2(ks2, 1);
 | 
|---|
| 140 |         double bs2 = ks2;
 | 
|---|
| 141 |         
 | 
|---|
| 142 | // Open the data file
 | 
|---|
| 143 |         fits *datafile = new fits( datafilename );
 | 
|---|
| 144 |         if (!datafile) {
 | 
|---|
| 145 |             printf( "Could not open the file: %s\n", datafilename );
 | 
|---|
| 146 |                   return 1;
 | 
|---|
| 147 |         }
 | 
|---|
| 148 |         
 | 
|---|
| 149 | // access data
 | 
|---|
| 150 |         NEvents = FOpenDataFile( *datafile );
 | 
|---|
| 151 |         printf("number of events in file: %d\n", NEvents);
 | 
|---|
| 152 |   if ( nevents == -1 || nevents > NEvents ) nevents = NEvents; // -1 means all!
 | 
|---|
| 153 | 
 | 
|---|
| 154 | //Get the DRS calibration
 | 
|---|
| 155 |         FOpenCalibFile( drsfilename,
 | 
|---|
| 156 |                                         drs_basemean,
 | 
|---|
| 157 |                                         drs_gainmean,
 | 
|---|
| 158 |                                         drs_triggeroffsetmean,
 | 
|---|
| 159 |                                         drs_n);
 | 
|---|
| 160 | 
 | 
|---|
| 161 | //Check the sizes of the data columns
 | 
|---|
| 162 |         if (drs_n != 1474560) {
 | 
|---|
| 163 |                 cerr << "error: DRS calib file has wrong ...erm...size ... drs_n is: "
 | 
|---|
| 164 |                         << drs_n << endl;
 | 
|---|
| 165 |                 cerr << " Aborting." << endl;
 | 
|---|
| 166 |                 return 1;
 | 
|---|
| 167 |         }
 | 
|---|
| 168 | 
 | 
|---|
| 169 |         if(ROIxNOP != 1474560)
 | 
|---|
| 170 |         {
 | 
|---|
| 171 |                 cout << "warning: data_n should better be 1440x1024=1474560, but is "
 | 
|---|
| 172 |                         << ROIxNOP << endl;
 | 
|---|
| 173 |                 cout << "this script is not guaranteed to run under these "
 | 
|---|
| 174 |                         <<" circumstances....any way ... it is never guaranteed." << endl;
 | 
|---|
| 175 |         }
 | 
|---|
| 176 | // Book the histograms
 | 
|---|
| 177 | 
 | 
|---|
| 178 |         BookHistos( );
 | 
|---|
| 179 |         TCanvas * cSpektrum = new TCanvas();
 | 
|---|
| 180 |         cSpektrum->Divide(1,2);
 | 
|---|
| 181 |         cSpektrum->cd(1);
 | 
|---|
| 182 |         hAmplSpek_discri->Draw("COLZ");
 | 
|---|
| 183 |         cSpektrum->cd(2);
 | 
|---|
| 184 |         hAmplSpek_cfd->Draw("COLZ");
 | 
|---|
| 185 | //      TCanvas * cStartCell = new TCanvas();
 | 
|---|
| 186 | //      cStartCell->cd();
 | 
|---|
| 187 | //      hStartCell->Draw();
 | 
|---|
| 188 | //      hPixelCellData.Draw();
 | 
|---|
| 189 | 
 | 
|---|
| 190 |         float calibratedVoltage;
 | 
|---|
| 191 |         for ( int ev = firstevent; ev < firstevent + nevents; ev++) {
 | 
|---|
| 192 |                 datafile->GetRow( ev );
 | 
|---|
| 193 | 
 | 
|---|
| 194 | 
 | 
|---|
| 195 |                 for ( int pix = 0; pix < 1440; pix++ ){
 | 
|---|
| 196 |                         // This means: GetEvent
 | 
|---|
| 197 | 
 | 
|---|
| 198 |                         // this is a stupid hack ... these is more code at the
 | 
|---|
| 199 |                         // end of this loop to complete this hack ...
 | 
|---|
| 200 |                         if (PixelID != -1) {
 | 
|---|
| 201 |                                 pix = PixelID;
 | 
|---|
| 202 |                                 if (verbosityLevel > 0){
 | 
|---|
| 203 |                                         cout << "Processing Event number: " << CurrentEventID << "\t"
 | 
|---|
| 204 |                                                 << "Pixel number: "<< pix << endl;
 | 
|---|
| 205 | 
 | 
|---|
| 206 |                                 }
 | 
|---|
| 207 |                         }
 | 
|---|
| 208 | 
 | 
|---|
| 209 |                         if (verbosityLevel > 0){
 | 
|---|
| 210 |                                 if (pix % 20 ==0){
 | 
|---|
| 211 |                                         cout << "Processing Event number: " << CurrentEventID << "\t"
 | 
|---|
| 212 |                                                 << "Pixel number: "<< pix << endl;
 | 
|---|
| 213 |                                 }
 | 
|---|
| 214 |                         }
 | 
|---|
| 215 | 
 | 
|---|
| 216 |                         //hStartCell->Fill( pix, StartCellVector[pix] );
 | 
|---|
| 217 | 
 | 
|---|
| 218 |                         for ( unsigned int sl = 0; sl < RegionOfInterest; sl++){
 | 
|---|
| 219 |                                 //calibratedVoltage = getValue( sl, pix);
 | 
|---|
| 220 |                                 calibratedVoltage = correctDrsOffset( sl, pix);
 | 
|---|
| 221 |                                 if (verbosityLevel > 10){
 | 
|---|
| 222 |                                         printf("calibratedVoltage = %f\n", calibratedVoltage);
 | 
|---|
| 223 |                                 }
 | 
|---|
| 224 |                                 Ameas[ sl ] = calibratedVoltage;
 | 
|---|
| 225 |                                 if (spikeDebug){
 | 
|---|
| 226 |                                         oldh[tAmeas].SetBinContent(sl, calibratedVoltage);
 | 
|---|
| 227 |                                 }
 | 
|---|
| 228 |                         }
 | 
|---|
| 229 | 
 | 
|---|
| 230 |                         computeN1mean( RegionOfInterest );
 | 
|---|
| 231 |                         // operates on Ameas[] and writes to N1mean[]
 | 
|---|
| 232 |                         removeSpikes( RegionOfInterest );
 | 
|---|
| 233 |                         // operates on Ameas[] and N1mean[], then writes to Vcorr[]
 | 
|---|
| 234 | 
 | 
|---|
| 235 |                         if (spikeDebug){
 | 
|---|
| 236 |                                 for ( unsigned int sl = 0; sl < RegionOfInterest; sl++){
 | 
|---|
| 237 |                                         hPixelCellData.Fill(sl, Vcorr[ sl ]);
 | 
|---|
| 238 |                                 }
 | 
|---|
| 239 |                         }
 | 
|---|
| 240 | 
 | 
|---|
| 241 |                         // filter Vcorr with sliding average using FIR filter function
 | 
|---|
| 242 |                         factfir(b_slide , a_slide, k_slide, Vcorr, Vslide);
 | 
|---|
| 243 |                         // filter Vslide with CFD using FIR filter function
 | 
|---|
| 244 |                         factfir(b_cfd , a_cfd, k_cfd, Vslide, Vcfd);
 | 
|---|
| 245 |                         // filter Vcfd with sliding average using FIR filter function
 | 
|---|
| 246 |                         factfir(bs2 , as2, ks2, Vcfd, Vcfd2);
 | 
|---|
| 247 | 
 | 
|---|
| 248 |                         float myTHR = 3.5;
 | 
|---|
| 249 |                         vector<Region> *Regions =  discriminator( Vslide, myTHR,true , 120);
 | 
|---|
| 250 |                         for (unsigned int p=0; p<Regions->size(); p++ ){
 | 
|---|
| 251 |                                         hAmplSpek_discri->Fill(pix , Regions->at(p).maxVal);
 | 
|---|
| 252 |                         }
 | 
|---|
| 253 | 
 | 
|---|
| 254 |                         // peaks in Ameas[] are found by searching for zero crossings
 | 
|---|
| 255 |                         // in Vcfd2
 | 
|---|
| 256 |                         // first Argument 1 means ... *rising* edge
 | 
|---|
| 257 |                         // second Argument 1 means ... search with stepsize 1 ... 10 is okay as well
 | 
|---|
| 258 |                         vector<Region> * zXings = zerosearch( Vcfd2 , 1 , 1);
 | 
|---|
| 259 | 
 | 
|---|
| 260 |                         // zXings means "zero cross ings"
 | 
|---|
| 261 |                         ShiftRegionBy(*zXings, -ks2/2);
 | 
|---|
| 262 |                         EnlargeRegion(*zXings, 10, 10);
 | 
|---|
| 263 |                         findAbsMaxInRegions(*zXings, Vslide);
 | 
|---|
| 264 | 
 | 
|---|
| 265 |                         if (zXings->size() != 0 ){
 | 
|---|
| 266 |                                 for (unsigned int i=0; i<zXings->size(); i++){
 | 
|---|
| 267 |                                         if (verbosityLevel > 1){
 | 
|---|
| 268 |                                                 cout << zXings->at(i).maxPos << ":\t"
 | 
|---|
| 269 |                                                         << zXings->at(i).maxVal <<endl;
 | 
|---|
| 270 |                                         }
 | 
|---|
| 271 |                                         hAmplSpek_cfd->Fill(pix, zXings->at(i).maxVal);
 | 
|---|
| 272 |                                 }
 | 
|---|
| 273 |                         }
 | 
|---|
| 274 | 
 | 
|---|
| 275 |                         if ( spikeDebug ){
 | 
|---|
| 276 |                                 for ( unsigned int sl = 0; sl < RegionOfInterest; sl++){
 | 
|---|
| 277 |                                         oldh[tVslide].SetBinContent( sl, Vslide[sl] );
 | 
|---|
| 278 |                                         oldh[tVcfd].SetBinContent( sl, Vcfd[sl] );
 | 
|---|
| 279 |                                         oldh[tVcfd2].SetBinContent( sl, Vcfd2[sl] );
 | 
|---|
| 280 |                                 }
 | 
|---|
| 281 | 
 | 
|---|
| 282 |                                 CW->cd( 1);
 | 
|---|
| 283 |                                 gPad->SetGrid();
 | 
|---|
| 284 |                                 oldh[tAmeas].Draw();
 | 
|---|
| 285 | 
 | 
|---|
| 286 | //                              CW->cd( tN1mean + 1);
 | 
|---|
| 287 | //                              gPad->SetGrid();
 | 
|---|
| 288 | //                              oldh[tN1mean].Draw();
 | 
|---|
| 289 | 
 | 
|---|
| 290 |                                 CW->cd( 2);
 | 
|---|
| 291 |                                 gPad->SetGrid();
 | 
|---|
| 292 |                                 oldh[tVcorr].Draw();
 | 
|---|
| 293 | 
 | 
|---|
| 294 |                                 //cFilter->cd( Ntypes - tVslide );
 | 
|---|
| 295 |                                 cFilter->cd(1);
 | 
|---|
| 296 |                                 gPad->SetGrid();
 | 
|---|
| 297 |                                 oldh[tVslide].Draw();
 | 
|---|
| 298 |                                 TLine thrLine(0, myTHR, 1024, myTHR);
 | 
|---|
| 299 |                                 thrLine.SetLineColor(kRed);
 | 
|---|
| 300 |                                 thrLine.SetLineWidth(2);
 | 
|---|
| 301 |                                 thrLine.Draw();
 | 
|---|
| 302 |                                 TLine * OneLine;
 | 
|---|
| 303 |                                 vector<TLine*> MyLines;
 | 
|---|
| 304 |                                 for (unsigned int p=0; p<Regions->size(); p++ ){
 | 
|---|
| 305 |                                         OneLine = new TLine(
 | 
|---|
| 306 |                                                 Regions->at(p).begin ,
 | 
|---|
| 307 |                                                 Regions->at(p).maxVal,
 | 
|---|
| 308 |                                                 Regions->at(p).end,
 | 
|---|
| 309 |                                                 Regions->at(p).maxVal);
 | 
|---|
| 310 |                                         OneLine->SetLineColor(kRed);
 | 
|---|
| 311 |                                         OneLine->SetLineWidth(2);
 | 
|---|
| 312 |                                         MyLines.push_back(OneLine);
 | 
|---|
| 313 |                                         OneLine->Draw();
 | 
|---|
| 314 |                                 }
 | 
|---|
| 315 |                                 TBox *OneBox;
 | 
|---|
| 316 |                                 vector<TBox*> MyBoxes;
 | 
|---|
| 317 |                                 for (unsigned int i=0; i<zXings->size(); i++){
 | 
|---|
| 318 |                                         OneBox = new TBox(
 | 
|---|
| 319 |                                                 zXings->at(i).maxPos -10 ,
 | 
|---|
| 320 |                                                 zXings->at(i).maxVal -0.5,
 | 
|---|
| 321 |                                                 zXings->at(i).maxPos +10 ,
 | 
|---|
| 322 |                                                 zXings->at(i).maxVal +0.5);
 | 
|---|
| 323 |                                         OneBox->SetLineColor(kBlue);
 | 
|---|
| 324 |                                         OneBox->SetLineWidth(1);
 | 
|---|
| 325 |                                         OneBox->SetFillStyle(0);
 | 
|---|
| 326 |                                         OneBox->SetFillColor(kRed);
 | 
|---|
| 327 |                                         MyBoxes.push_back(OneBox);
 | 
|---|
| 328 |                                         OneBox->Draw();
 | 
|---|
| 329 |                                 }
 | 
|---|
| 330 | 
 | 
|---|
| 331 | //                              cFilter->cd( Ntypes - tVcfd );
 | 
|---|
| 332 | //                              cFilter->cd(2);
 | 
|---|
| 333 | //                              gPad->SetGrid();
 | 
|---|
| 334 | //                              oldh[tVcfd].Draw();
 | 
|---|
| 335 | //                              TLine zeroline(0, 0, 1024, 0);
 | 
|---|
| 336 | //                              zeroline.SetLineColor(kBlue);
 | 
|---|
| 337 | //                              zeroline.Draw();
 | 
|---|
| 338 | 
 | 
|---|
| 339 |                                 //cFilter->cd( Ntypes - tVcfd2 );
 | 
|---|
| 340 |                                 cFilter->cd(2);
 | 
|---|
| 341 |                                 gPad->SetGrid();
 | 
|---|
| 342 |                                 oldh[tVcfd2].Draw();
 | 
|---|
| 343 |                                 TLine zeroline(0, 0, 1024, 0);
 | 
|---|
| 344 |                                 zeroline.SetLineColor(kBlue);
 | 
|---|
| 345 |                                 zeroline.Draw();
 | 
|---|
| 346 | 
 | 
|---|
| 347 | 
 | 
|---|
| 348 | 
 | 
|---|
| 349 | 
 | 
|---|
| 350 | 
 | 
|---|
| 351 |                                 CW->Update();
 | 
|---|
| 352 |                                 cFilter->Update();
 | 
|---|
| 353 | 
 | 
|---|
| 354 |                                 //Process gui events asynchronously during input
 | 
|---|
| 355 |                                 TTimer timer("gSystem->ProcessEvents();", 50, kFALSE);
 | 
|---|
| 356 |                                 timer.TurnOn();
 | 
|---|
| 357 |                                 TString input = Getline("Type 'q' to exit, <return> to go on: ");
 | 
|---|
| 358 |                                 timer.TurnOff();
 | 
|---|
| 359 |                                 if (input=="q\n") {
 | 
|---|
| 360 |                                         breakout=true;
 | 
|---|
| 361 |                                         delete OneLine;
 | 
|---|
| 362 | //                                      for (unsigned int i=0; i<MyLines.size(); i++){
 | 
|---|
| 363 | //                                              delete MyLines[i];
 | 
|---|
| 364 | //                                      }
 | 
|---|
| 365 | //                                      for (unsigned int i=0; i<MyBoxes.size(); i++){
 | 
|---|
| 366 | //                                              delete MyBoxes[i];
 | 
|---|
| 367 | //                                      }
 | 
|---|
| 368 |                                 }
 | 
|---|
| 369 |                         }// end of if(spikeDebug)
 | 
|---|
| 370 |                         delete Regions;
 | 
|---|
| 371 |                         delete zXings;
 | 
|---|
| 372 | 
 | 
|---|
| 373 |                         if (PixelID != -1){
 | 
|---|
| 374 |                                 pix = 2000;
 | 
|---|
| 375 |                         }
 | 
|---|
| 376 | 
 | 
|---|
| 377 |                         if (breakout)
 | 
|---|
| 378 |                                 break;
 | 
|---|
| 379 | 
 | 
|---|
| 380 |                 } // end of loop over pixels
 | 
|---|
| 381 | 
 | 
|---|
| 382 |         if (ev % 10 ==0){
 | 
|---|
| 383 |                 cSpektrum->Modified();
 | 
|---|
| 384 |                 cSpektrum->Update();
 | 
|---|
| 385 | //              cStartCell->Modified();
 | 
|---|
| 386 | //              cStartCell->Update();
 | 
|---|
| 387 |         }
 | 
|---|
| 388 |                 if (breakout)
 | 
|---|
| 389 |                         break;
 | 
|---|
| 390 |         }       // end of loop over pixels
 | 
|---|
| 391 | 
 | 
|---|
| 392 | 
 | 
|---|
| 393 | //      delete cStartCell;
 | 
|---|
| 394 |         return( 0 );
 | 
|---|
| 395 | }
 | 
|---|
| 396 | 
 | 
|---|
| 397 | void removeSpikes(int Samples){
 | 
|---|
| 398 | 
 | 
|---|
| 399 |     const float fract = 0.8;
 | 
|---|
| 400 |     float x, xp, xpp, x3p;
 | 
|---|
| 401 | 
 | 
|---|
| 402 |     // assume that there are no spikes
 | 
|---|
| 403 |     for ( int i = 0; i <  Samples; i++) Vcorr[i] = Ameas[i];
 | 
|---|
| 404 | 
 | 
|---|
| 405 | // find the spike and replace it by mean value of neighbors
 | 
|---|
| 406 |     for ( int i = 0; i < Samples; i++) {
 | 
|---|
| 407 | 
 | 
|---|
| 408 |     // printf("Vcorr[%d] = %f, Ameas[%d] = %f\n", i, Vcorr[ i ], i, Ameas[ i ] );
 | 
|---|
| 409 | 
 | 
|---|
| 410 |     x = Ameas[i] - N1mean[i];
 | 
|---|
| 411 | 
 | 
|---|
| 412 |         if ( x < -5. ){ // a spike candidate
 | 
|---|
| 413 |             // check consistency with a single channel spike
 | 
|---|
| 414 |             xp = Ameas[i+1] - N1mean[i+1];
 | 
|---|
| 415 |             xpp = Ameas[i+2] - N1mean[i+2];
 | 
|---|
| 416 |             x3p = Ameas[i+3] - N1mean[i+3];
 | 
|---|
| 417 | 
 | 
|---|
| 418 |             // printf("candidates x[%d] = %f; xp = %f; xpp = %f, x3p = %f\n", i, x, xp, xpp, x3p);
 | 
|---|
| 419 | 
 | 
|---|
| 420 |             if ( Ameas[i+2] - ( Ameas[i] + Ameas[i+3] )/2. > 10. ){
 | 
|---|
| 421 |                 // printf("double spike candidate\n");
 | 
|---|
| 422 |                 Vcorr[i+1] = ( Ameas[i] + Ameas[i+3] )/2.;
 | 
|---|
| 423 |                 Vcorr[i+2] = ( Ameas[i] + Ameas[i+3] )/2.;
 | 
|---|
| 424 |                 // printf("Vcorr[%d] = %f %f %f %f\n", i, Vcorr[i], Vcorr[i+1], Vcorr[i+2], Vcorr[ i+3 ]);
 | 
|---|
| 425 |                 // printf("Ameas[%d] = %f %f %f %f\n", i, Ameas[ i ], Ameas[ i+1 ], Ameas[ i+2 ], Ameas[i+3]);
 | 
|---|
| 426 |                 i = i + 3; 
 | 
|---|
| 427 |             }
 | 
|---|
| 428 |             else{
 | 
|---|
| 429 |     
 | 
|---|
| 430 |                 if ( ( xp > -2.*x*fract ) && ( xpp < -10. ) ){
 | 
|---|
| 431 |                     Vcorr[i+1] = N1mean[i+1];
 | 
|---|
| 432 |                     // printf("Vcorr[%d] = %f %f %f\n", i, Vcorr[i], Vcorr[i+1], Vcorr[i+2]);
 | 
|---|
| 433 |                     // N1mean[i+1] = (Ameas[i] - Ameas[i+2] / 2.);
 | 
|---|
| 434 |                     N1mean[i+2] = (Ameas[i+1] - Ameas[i+3] / 2.);
 | 
|---|
| 435 |                     i = i + 2;//do not care about the next sample it was the spike
 | 
|---|
| 436 |                 }
 | 
|---|
| 437 |                 // treatment for the end of the pipeline must be added !!!
 | 
|---|
| 438 |             }
 | 
|---|
| 439 |         }
 | 
|---|
| 440 |         else{
 | 
|---|
| 441 |              // do nothing
 | 
|---|
| 442 |         }
 | 
|---|
| 443 |     } // end of spike search and correction
 | 
|---|
| 444 |     for ( int i = 0; i < Samples; i++ ) oldh[ tVcorr ].SetBinContent( i, Vcorr[i] );
 | 
|---|
| 445 | }
 | 
|---|
| 446 | /*
 | 
|---|
| 447 | void computeN1mean( int Samples ){
 | 
|---|
| 448 | cout << "In compute N1mean" << endl;
 | 
|---|
| 449 | // compute the mean of the left and right neighbors of a channel
 | 
|---|
| 450 | 
 | 
|---|
| 451 |     for( int i = 0; i < Samples; i++){
 | 
|---|
| 452 |         if (i == 0){ // use right sample es mean
 | 
|---|
| 453 |             N1mean[i] = Ameas[i+1];
 | 
|---|
| 454 |         }
 | 
|---|
| 455 |         else if ( i == Samples-1 ){ //use left sample as mean
 | 
|---|
| 456 |             N1mean[i] = Ameas[i-1];
 | 
|---|
| 457 |         }
 | 
|---|
| 458 |         else{
 | 
|---|
| 459 |             N1mean[i] = ( Ameas[i-1] + Ameas[i+1] ) / 2.;
 | 
|---|
| 460 |         }
 | 
|---|
| 461 |         h[tN1mean].SetBinContent(i, Ameas[i] - N1mean[i]);
 | 
|---|
| 462 |     }
 | 
|---|
| 463 | } // end of computeN1mean computation
 | 
|---|
| 464 | */
 | 
|---|
| 465 | 
 | 
|---|
| 466 | void computeN1mean( int Samples ){
 | 
|---|
| 467 | // compute the mean of the left and right neighbors of a channel
 | 
|---|
| 468 | 
 | 
|---|
| 469 |     for( int i = 2; i < Samples - 2; i++){
 | 
|---|
| 470 | /*        if (i == 0){ // use right sample es mean
 | 
|---|
| 471 |             N1mean[i] = Ameas[i+1];
 | 
|---|
| 472 |         }
 | 
|---|
| 473 |         else if ( i == Samples-1 ){ //use left sample as mean
 | 
|---|
| 474 |             N1mean[i] = Ameas[i-1];
 | 
|---|
| 475 |         }
 | 
|---|
| 476 |         else{
 | 
|---|
| 477 |             N1mean[i] = ( Ameas[i-1] + Ameas[i+1] ) / 2.;
 | 
|---|
| 478 |         }
 | 
|---|
| 479 | */
 | 
|---|
| 480 |         N1mean[i] = ( Ameas[i-1] + Ameas[i+1] ) / 2.;
 | 
|---|
| 481 |     }
 | 
|---|
| 482 | } // end of computeN1mean computation
 | 
|---|
| 483 | 
 | 
|---|
| 484 | float getValue( int slice, int pixel ){
 | 
|---|
| 485 |         const float dconv = 2000/4096.0;
 | 
|---|
| 486 | 
 | 
|---|
| 487 |         float vraw, vcal;
 | 
|---|
| 488 | 
 | 
|---|
| 489 |         unsigned int pixel_pt;
 | 
|---|
| 490 |         unsigned int slice_pt;
 | 
|---|
| 491 |         unsigned int cal_pt;
 | 
|---|
| 492 |         unsigned int drs_cal_offset;
 | 
|---|
| 493 | 
 | 
|---|
| 494 |         // printf("pixel = %d, slice = %d\n", slice, pixel);
 | 
|---|
| 495 | 
 | 
|---|
| 496 |         pixel_pt = pixel * RegionOfInterest;
 | 
|---|
| 497 |         slice_pt = pixel_pt + slice;
 | 
|---|
| 498 |         drs_cal_offset = ( slice + StartCellVector[ pixel ] )%RegionOfInterest;
 | 
|---|
| 499 |         cal_pt    = pixel_pt + drs_cal_offset;
 | 
|---|
| 500 | 
 | 
|---|
| 501 |         vraw = AllPixelDataVector[ slice_pt ] * dconv;
 | 
|---|
| 502 |         vcal = ( vraw - drs_basemean[ cal_pt ] - drs_triggeroffsetmean[ slice_pt ] ) / drs_gainmean[ cal_pt ]*1907.35;
 | 
|---|
| 503 | 
 | 
|---|
| 504 |         return( vcal );
 | 
|---|
| 505 | }
 | 
|---|
| 506 | float correctDrsOffset( int slice, int pixel ){
 | 
|---|
| 507 | 
 | 
|---|
| 508 |         const float dconv = 2000/4096.0;
 | 
|---|
| 509 | 
 | 
|---|
| 510 |         // here 1024 is not the RegionOfInterest, but really the lenth of the pipeline
 | 
|---|
| 511 |         unsigned int physical_slice = ( slice + StartCellVector[ pixel ] ) % 1024;
 | 
|---|
| 512 | 
 | 
|---|
| 513 |         unsigned int slice_pt;
 | 
|---|
| 514 |         unsigned int physical_slice_pt;
 | 
|---|
| 515 |         slice_pt                        = pixel * RegionOfInterest + slice;
 | 
|---|
| 516 |         physical_slice_pt       = pixel * RegionOfInterest + physical_slice;
 | 
|---|
| 517 | 
 | 
|---|
| 518 |         float vcal = AllPixelDataVector[ slice_pt ] *
 | 
|---|
| 519 |                 dconv - drs_basemean[ physical_slice_pt ];
 | 
|---|
| 520 |         return( vcal );
 | 
|---|
| 521 | }
 | 
|---|
| 522 | 
 | 
|---|
| 523 | void BookHistos( ){
 | 
|---|
| 524 | // booking and parameter settings for all histos
 | 
|---|
| 525 | 
 | 
|---|
| 526 |         // histograms for baseline extraction
 | 
|---|
| 527 |         char hName[500];
 | 
|---|
| 528 |         char hTitle[500];
 | 
|---|
| 529 | 
 | 
|---|
| 530 |         TH1F *h;
 | 
|---|
| 531 | 
 | 
|---|
| 532 |         for( int i = 0; i < NPIX; i++ ) {
 | 
|---|
| 533 |                 sprintf(&hTitle[0],"all events all slices of pixel %d", i);
 | 
|---|
| 534 |                 sprintf(&hName[0],"base%d", i);
 | 
|---|
| 535 | 
 | 
|---|
| 536 |                 h = new TH1F( hName, hTitle, 400, -99.5 ,100.5 );
 | 
|---|
| 537 | 
 | 
|---|
| 538 |                 h->GetXaxis()->SetTitle( "Sample value (mV)" );
 | 
|---|
| 539 |                 h->GetYaxis()->SetTitle( "Entries / 0.5 mV" );
 | 
|---|
| 540 |                 hListBaseline.Add( h );
 | 
|---|
| 541 |                 hBaseline[i] = h;
 | 
|---|
| 542 |         }
 | 
|---|
| 543 | 
 | 
|---|
| 544 |     hMeanBsl = new TH1F("histo_mean","Value of maximal probability",400,-99.5,100.5);
 | 
|---|
| 545 |     hMeanBsl->GetXaxis()->SetTitle( "max value (mV)" );
 | 
|---|
| 546 |     hMeanBsl->GetYaxis()->SetTitle( "Entries / 0.5 mV" );
 | 
|---|
| 547 |     hList.Add( hMeanBsl );
 | 
|---|
| 548 | 
 | 
|---|
| 549 |     hpltMeanBsl = new TH1F("hplt_mean","Value of maximal probability",1440,-0.5,1439.5);
 | 
|---|
| 550 |     hpltMeanBsl->GetXaxis()->SetTitle( "pixel" );
 | 
|---|
| 551 |     hpltMeanBsl->GetYaxis()->SetTitle( "max value in mV" );
 | 
|---|
| 552 |     hList.Add( hpltMeanBsl );
 | 
|---|
| 553 | 
 | 
|---|
| 554 |     hRmsBsl = new TH1F("histo_rms","RMS in mV",2000,-99.5,100.5);
 | 
|---|
| 555 |     hRmsBsl->GetXaxis()->SetTitle( "RMS (mV)" );
 | 
|---|
| 556 |     hRmsBsl->GetYaxis()->SetTitle( "Entries / 0.5 mV" );
 | 
|---|
| 557 |     hList.Add( hRmsBsl );
 | 
|---|
| 558 | 
 | 
|---|
| 559 |     hpltRmsBsl = new TH1F("hplt_rms","Value of maximal probability",1440,-0.5,1439.5);
 | 
|---|
| 560 |     hpltRmsBsl->GetXaxis()->SetTitle( "pixel" );
 | 
|---|
| 561 |     hpltRmsBsl->GetYaxis()->SetTitle( "RMS in mV" );
 | 
|---|
| 562 |     hList.Add( hpltRmsBsl );
 | 
|---|
| 563 | 
 | 
|---|
| 564 |         hAmplSpek_cfd = new TH2F("hAmplSpek_cfd","amplitude spektrum - CFD",1440,-0.5,1439.5, 256, -27.5, 100.5);
 | 
|---|
| 565 |         hAmplSpek_cfd->GetXaxis()->SetTitle( "pixel" );
 | 
|---|
| 566 |         hAmplSpek_cfd->GetYaxis()->SetTitle( "amplitude in mV" );
 | 
|---|
| 567 |         hList.Add( hAmplSpek_cfd );
 | 
|---|
| 568 | 
 | 
|---|
| 569 |         hAmplSpek_discri = new TH2F("hAmplSpek_discri","amplitude spektrum - std discriminator",1440,-0.5,1439.5, 256, -27.5, 100.5);
 | 
|---|
| 570 |         hAmplSpek_discri->GetXaxis()->SetTitle( "pixel" );
 | 
|---|
| 571 |         hAmplSpek_discri->GetXaxis()->SetTitle( "amplitude in mV" );
 | 
|---|
| 572 |         hList.Add( hAmplSpek_discri );
 | 
|---|
| 573 | 
 | 
|---|
| 574 |         hStartCell = new TH2F("StartCell", "StartCell", 1440, 0., 1440., 1024, 0., 1024);
 | 
|---|
| 575 |         hStartCell->GetXaxis()->SetTitle( "pixel" );
 | 
|---|
| 576 |         hStartCell->GetXaxis()->SetTitle( "slice" );
 | 
|---|
| 577 |         hList.Add( hStartCell );
 | 
|---|
| 578 | 
 | 
|---|
| 579 |     oldh = new TH1F[ Ntypes ];
 | 
|---|
| 580 | 
 | 
|---|
| 581 |     for ( int type = 0; type < Ntypes; type++){
 | 
|---|
| 582 | 
 | 
|---|
| 583 |         oldh[ type ].SetBins(1024, 0, 1024);
 | 
|---|
| 584 |         oldh[ type ].SetLineColor(1);
 | 
|---|
| 585 |         oldh[ type ].SetLineWidth(2);
 | 
|---|
| 586 | 
 | 
|---|
| 587 |         // set X axis paras
 | 
|---|
| 588 |         oldh[ type ].GetXaxis()->SetLabelSize(0.1);
 | 
|---|
| 589 |         oldh[ type ].GetXaxis()->SetTitleSize(0.1);
 | 
|---|
| 590 |         oldh[ type ].GetXaxis()->SetTitleOffset(1.2);
 | 
|---|
| 591 |         oldh[ type ].GetXaxis()->SetTitle(Form("Time slice (%.1f ns/slice)", 1./2.));
 | 
|---|
| 592 | 
 | 
|---|
| 593 |         // set Y axis paras
 | 
|---|
| 594 |         oldh[ type ].GetYaxis()->SetLabelSize(0.1);
 | 
|---|
| 595 |         oldh[ type ].GetYaxis()->SetTitleSize(0.1);
 | 
|---|
| 596 |         oldh[ type ].GetYaxis()->SetTitleOffset(0.3);
 | 
|---|
| 597 |         oldh[ type ].GetYaxis()->SetTitle("Amplitude (a.u.)");
 | 
|---|
| 598 |     }
 | 
|---|
| 599 | 
 | 
|---|
| 600 |     CW = new TCanvas("CW","DRS Waveform",10,10,800,600);
 | 
|---|
| 601 |         CW->Divide(1, 2);
 | 
|---|
| 602 |     cFilter = new TCanvas("cFilter","filtered DRS Waveforms",10,10,800,600);
 | 
|---|
| 603 |         cFilter->Divide(1, 2);
 | 
|---|
| 604 | }
 | 
|---|
| 605 | 
 | 
|---|
| 606 | void SaveHistograms( char * loc_fname ){
 | 
|---|
| 607 | 
 | 
|---|
| 608 |         TString fName; // name of the histogram file
 | 
|---|
| 609 | 
 | 
|---|
| 610 |         // create the filename for the histogram file
 | 
|---|
| 611 |         fName = loc_fname; // use the name of the tree file
 | 
|---|
| 612 |         // TODO ... next statement doesn't work for ".fits.gz"
 | 
|---|
| 613 |         fName.Remove(fName.Length() - 5); // remove the extension .fits
 | 
|---|
| 614 |         fName += "_discri.root"; // add the new extension
 | 
|---|
| 615 | 
 | 
|---|
| 616 |         // create the histogram file (replace if already existing)
 | 
|---|
| 617 |         TFile tf( fName, "RECREATE");
 | 
|---|
| 618 | 
 | 
|---|
| 619 |         hList.Write(); // write the major histograms into the top level directory
 | 
|---|
| 620 |         tf.mkdir("BaselineHisto");
 | 
|---|
| 621 |         tf.cd("BaselineHisto"); // go to new subdirectory
 | 
|---|
| 622 |         hListBaseline.Write(); // write histos into subdirectory
 | 
|---|
| 623 | 
 | 
|---|
| 624 |         tf.Close(); // close the file
 | 
|---|
| 625 | } // end of SaveHistograms( char * loc_fname )
 | 
|---|
| 626 | 
 | 
|---|
| 627 | int FOpenDataFile(fits &datafile){
 | 
|---|
| 628 | 
 | 
|---|
| 629 | //      cout << "-------------------- Data Header -------------------" << endl;
 | 
|---|
| 630 | //      datafile.PrintKeys();
 | 
|---|
| 631 | //      cout << "------------------- Data Columns -------------------" << endl;
 | 
|---|
| 632 | //      datafile.PrintColumns();
 | 
|---|
| 633 | 
 | 
|---|
| 634 |         //Get the size of the data column
 | 
|---|
| 635 |         RegionOfInterest        = datafile.GetUInt("NROI");
 | 
|---|
| 636 |         NumberOfPixels          = datafile.GetUInt("NPIX");
 | 
|---|
| 637 |         //Size of column "Data" = #Pixel x ROI
 | 
|---|
| 638 |         ROIxNOP                         = datafile.GetN("Data");
 | 
|---|
| 639 | 
 | 
|---|
| 640 |         //Set the sizes of the data vectors
 | 
|---|
| 641 |         AllPixelDataVector.resize(ROIxNOP,0);
 | 
|---|
| 642 |         StartCellVector.resize(NumberOfPixels,0);
 | 
|---|
| 643 | 
 | 
|---|
| 644 |         //Link the data to variables
 | 
|---|
| 645 |         datafile.SetRefAddress("EventNum", CurrentEventID);
 | 
|---|
| 646 |         datafile.SetVecAddress("Data", AllPixelDataVector);
 | 
|---|
| 647 |         datafile.SetVecAddress("StartCellData", StartCellVector);
 | 
|---|
| 648 |         datafile.GetRow(0);
 | 
|---|
| 649 | 
 | 
|---|
| 650 |         return datafile.GetNumRows() ;
 | 
|---|
| 651 | }
 | 
|---|
| 652 | 
 | 
|---|
| 653 | 
 | 
|---|
| 654 | /////////////////////////////////////////////////////////////////////////////
 | 
|---|
| 655 | /// old BookHistos
 | 
|---|
| 656 | /*
 | 
|---|
| 657 | void BookHistos( int Samples ){
 | 
|---|
| 658 | // booking and parameter settings for all histos
 | 
|---|
| 659 | 
 | 
|---|
| 660 |     h = new TH1F[ Ntypes ];
 | 
|---|
| 661 | 
 | 
|---|
| 662 |     for ( int type = 0; type < Ntypes; type++){
 | 
|---|
| 663 | 
 | 
|---|
| 664 |         h[ type ].SetBins(Samples, 0, Samples);
 | 
|---|
| 665 |         h[ type ].SetLineColor(1);
 | 
|---|
| 666 |         h[ type ].SetLineWidth(2);
 | 
|---|
| 667 | 
 | 
|---|
| 668 |         // set X axis paras
 | 
|---|
| 669 |         h[ type ].GetXaxis()->SetLabelSize(0.1);
 | 
|---|
| 670 |         h[ type ].GetXaxis()->SetTitleSize(0.1);
 | 
|---|
| 671 |         h[ type ].GetXaxis()->SetTitleOffset(1.2);
 | 
|---|
| 672 |         h[ type ].GetXaxis()->SetTitle(Form("Time slice (%.1f ns/slice)", 1./2.));
 | 
|---|
| 673 | 
 | 
|---|
| 674 |         // set Y axis paras
 | 
|---|
| 675 |         h[ type ].GetYaxis()->SetLabelSize(0.1);
 | 
|---|
| 676 |         h[ type ].GetYaxis()->SetTitleSize(0.1);
 | 
|---|
| 677 |         h[ type ].GetYaxis()->SetTitleOffset(0.3);
 | 
|---|
| 678 |         h[ type ].GetYaxis()->SetTitle("Amplitude (a.u.)");
 | 
|---|
| 679 |     }
 | 
|---|
| 680 |     CW = new TCanvas("CW","DRS Waveform",10,10,800,600);
 | 
|---|
| 681 |     CW->Divide(1, 3);
 | 
|---|
| 682 |     cFilter = new TCanvas("cFilter","filtered DRS Waveforms",10,10,800,600);
 | 
|---|
| 683 |     cFilter->Divide(1, 3);
 | 
|---|
| 684 | 
 | 
|---|
| 685 |     hStartCell = new TH2F("StartCell", "StartCell", 1440, 0., 1440., 1024, 0., 1024);
 | 
|---|
| 686 | 
 | 
|---|
| 687 | }
 | 
|---|
| 688 | */
 | 
|---|
| 689 | 
 | 
|---|