Changeset 13641
- Timestamp:
- 05/10/12 16:44:09 (13 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
fact/tools/rootmacros/PulseTemplates/FCalcPulseTemplate.C
r13465 r13641 1 1 /********************** FPulseTemplate *********************** 2 * read FACT raw data 3 * remove spikes 4 * calculate baseline 5 * find peaks (CFD and normal discriminator) 2 * open a root-file with pulse overlay histograms 3 * calculate most propable Pulse-Shape 6 4 * compute pulse height and pulse integral spektrum of the peaks 7 + search for Peaks in data8 + put peaks into different histos depending und Amplitude9 5 + draw histos for single, double, tripple, ....above 100mV Photonpulses 10 6 + draw Parameterdevelopment depending on photon quantity … … 19 15 #include <TROOT.h> 20 16 #include <TCanvas.h> 21 #include <TProfile.h>22 17 #include <TTimer.h> 23 #include <TH1F.h>24 #include <TH2F.h>25 18 #include <Getline.h> 26 #include <TLine.h>27 #include <TBox.h>28 19 #include <TMath.h> 29 20 #include <TFile.h> 30 21 #include <TStyle.h> 31 22 #include <TString.h> 32 #include <TProfile.h>33 23 34 24 #include <stdio.h> … … 48 38 // rootmacros 49 39 //---------------------------------------------------------------------------- 50 51 #include "../fits.h" 52 53 #include "../openFits.h" 54 #include "../openFits.c" 55 56 #include "../discriminator.h" 57 #include "../discriminator.C" 58 #include "../zerosearch.h" 59 #include "../zerosearch.C" 60 #include "../factfir.C" 61 62 #include "../DrsCalibration.C" 63 #include "../DrsCalibration.h" 64 65 #include "../SpikeRemoval.h" 66 #include "../SpikeRemoval.C" 40 #include "rootfilehandler.h" 41 #include "pixel.h" 42 #include "templateextractors.h" 67 43 68 44 //---------------------------------------------------------------------------- 69 45 // Decleration of global variables 70 46 //---------------------------------------------------------------------------- 71 72 47 bool breakout=false; 73 48 74 int gNEvents; 75 vector<int16_t> AllPixelDataVector; 76 vector<int16_t> StartCellVector; 77 unsigned int CurrentEventID; 78 size_t PXLxROI; 79 UInt_t gRegionOfInterest; 80 UInt_t NumberOfPixels; 81 TString histotitle; 82 83 size_t TriggerOffsetROI, RC; 84 size_t drs_n; 85 vector<float> drs_basemean; 86 vector<float> drs_gainmean; 87 vector<float> drs_triggeroffsetmean; 88 89 vector<float> Ameas(FAD_MAX_SAMPLES); // copy of the data (measured amplitude 90 vector<float> Vcorr(FAD_MAX_SAMPLES); // corrected Values 91 vector<float> Vslide(FAD_MAX_SAMPLES); // sliding average result 92 vector<float> Vcfd(FAD_MAX_SAMPLES); // CDF result 93 vector<float> Vcfd2(FAD_MAX_SAMPLES); // CDF result + 2nd sliding average 94 95 float gGainMean = 9; 96 float gBSLMean = 0; 97 98 typedef struct{ 99 double maxAmpl; 100 int countOfMax; 101 } OverlayMaximum; 102 103 //typedef struct{ 104 // TString name; 105 // TString title; 106 // TString xTitle; 107 // TString yTitle; 108 // float yMax; 109 // float yMin; 110 //} histAttrib_t; 111 112 //histAttrib_t* gHistoAttrib[MAX_PULS_ORDER]; 113 //histAttrib_t* gMaximumAttrib[MAX_PULS_ORDER]; 114 //histAttrib_t* gMaxGausAttrib[MAX_PULS_ORDER]; 115 //histAttrib_t* gProfileAttrib[MAX_PULS_ORDER]; 116 117 118 // histograms 119 const int Number_Of_Debug_Histo_Types = 7; 120 121 const unsigned int 122 Ameas_ = 0, 123 N1mean_ = 1, 124 Vcorr_ = 2, 125 Vtest_ = 3, 126 Vslide_ = 4, 127 Vcfd_ = 5, 128 Vcfd2_ = 6; 129 130 //const char* gHistoTypes[8] = { 131 // "hPixelOverlay", 132 // "hPixelMax", 133 // "hPixelEdgeOverlay", 134 // "hPixelProfile", 135 // "hAllPixelOverlay", 136 // "hAllPixelMax", 137 // "hAllPixelMaxGaus", 138 // "hAllPixelProfile" 139 //}; 140 141 //---------------------------------------------------------------------------- 142 // Initialisation of histograms 49 //---------------------------------------------------------------------------- 50 // Initialisation of Root Objects 143 51 //---------------------------------------------------------------------------- 144 52 145 53 // Temporary Objects 146 TH1F* debugHistos = NULL;147 TH2F* hOverlayTemp = NULL;148 TH1F* hMaximumTemp = NULL;149 TH1D* hProjPeak = NULL;150 TH1F* hTesthisto = NULL;151 TH2F* hTesthisto2 = NULL;152 54 153 55 //Pixelwise Histograms 154 TH2F* hPixelOverlay[MAX_PULS_ORDER]= {NULL};//histogrammm for overlay of detected Peaks155 TH2F* hPixelEdgeOverlay[MAX_PULS_ORDER] = {NULL};156 TProfile* hPixelProfile[MAX_PULS_ORDER] = {NULL};//histogrammm for Profile of detected Peaks157 TProfile* hPixelProfile2[MAX_PULS_ORDER] = {NULL};//histogrammm for Profile of detected Peaks158 TH1F* hPixelMax[MAX_PULS_ORDER] = {NULL};159 TH1F* hPixelMedian[MAX_PULS_ORDER] = {NULL};160 TH1F* hPixelMean[MAX_PULS_ORDER] = {NULL};161 56 162 57 //All Pixel Histograms 163 TH2F* hAllPixelOverlay[MAX_PULS_ORDER] = {NULL};164 TProfile* hAllPixelProfile[MAX_PULS_ORDER] = {NULL};165 TProfile* hAllPixelProfile2[MAX_PULS_ORDER] = {NULL};166 TH1F* hAllPixelMax[MAX_PULS_ORDER] = {NULL};167 TH1F* hAllPixelMedian[MAX_PULS_ORDER] = {NULL};168 TH1F* hAllPixelMean[MAX_PULS_ORDER] = {NULL};169 TH2F* hAllPixelEdgeOverlay[MAX_PULS_ORDER] = {NULL};170 171 //Histogram Parameters172 Int_t gPixelOverlayXaxisLeft = 0;173 Int_t gPixelOverlayXaxisRight = 0;174 58 175 59 //Root-File Objects 176 TObjArray* hList = NULL; 177 TObjArray* hAllPixelList = NULL; 178 TObjArray* hRootList = NULL; 60 TObjArray* hAllPixelList = NULL; 61 TObjArray* hRootList = NULL; 179 62 180 63 //Canvases 181 TCanvas* cgpPixelPulses[MAX_PULS_ORDER] = {NULL};182 TCanvas* cgpAllPixelPulses[MAX_PULS_ORDER] = {NULL};183 TCanvas* cgpTestHistos= NULL;184 //TCanvas* gpcDevelopment= NULL;64 TCanvas** cgpPixelPulses = NULL; 65 TCanvas** cgpAllPixelPulses = NULL; 66 TCanvas** cgpDistributions = NULL; 67 //TCanvas* cgpTestHistos = NULL; 185 68 186 69 //---------------------------------------------------------------------------- 187 70 // Functions 188 71 //---------------------------------------------------------------------------- 189 190 void BookHistos(int ); 191 void BookPixelHistos(int ); 192 void BookTestHistos( int ); 193 194 void DeletePixelHistos(int ); 195 void SaveHistograms( const char*, const char*, TObjArray*, int ); 196 197 void FillHistograms(vector<Region>*, int, int, int); 198 void DrawPulseHistograms(int, int); 199 void DrawMaximumHistograms( int, int ); 200 void DrawTestHistograms( int); 201 202 void PlotPulseShapeOfMaxPropability(unsigned int, int, int); 203 void FitMaxPropabilityPuls( TH1F* , int ); 204 205 Double_t MedianOfH1( TH1 *h1 ); 206 void PlotMedianEachSliceOfPulse(unsigned int, int, int); 207 208 void UpdateCanvases( int, int); 209 void DeletePixelCanvases( int ); 210 211 void CreateRootFile(const char*, int ); 212 TFile *OpenRootFile(const char*, int ); 213 void CloseRootFile(TFile *tf); 214 215 TString CreateSubDirName(int); //creates a String containing the path to the subdirectory in root file 216 TString CreateSubDirName(const char* ); 217 218 void WritePixelTemplateToCsv( TString, const char*, const char*, int, int); 219 void WriteAllPixelTemplateToCsv( TString, const char*, const char*, int); 220 //void SetHistogrammSettings( const char*, int, int , int); 221 //void CreatePulseCanvas( TCanvas* , unsigned int, const char* , const char*, const char*, int); 222 //void DeletePulseCanvas( TCanvas* , unsigned int); 223 224 72 void DeletePixelCanvases( int, int ); 73 void 74 UpdateCanvases( 75 int verbosityLevel, 76 int max_pulse_order, 77 bool testmode 78 ); 225 79 //---------------------------------------------------------------------------- 226 80 //---------------------------------------------------------------------------- … … 229 83 //---------------------------------------------------------------------------- 230 84 int FCalcPulseTemplate( 231 char* datafilename = "../data/2011/11/09/20111109_006.fits.gz", 232 const char* drsfilename = "../data/2011/11/09/20111109_003.drs.fits.gz", 233 const char* OutRootFileName = "../analysis/FPulseTemplate/20111109_006/20111109_006.pulses.root", 234 const char* OutPutPath = "../analysis/FPulseTemplate/20111109_006/Templates/", 235 bool ProduceGraphic = false, 236 int refresh_rate = 500, //refresh rate for canvases 237 bool spikeDebug = false, 238 bool debugPixel = false, 239 bool fitdata = false, 240 int verbosityLevel = 0, // different verbosity levels can be implemented here 241 int firstevent = 0, 242 int nevents = -1, 243 int firstpixel = 0, 244 int npixel = -1, 245 int sampleSetSize = 36, 246 int AmplWindowWidth = 14, //Width of Window for selection of pulses to histograms 247 float GainMean = 8, 248 float BSLMean = -1, //4 Histogramms will be drawn, decide how far pulsehights differ from eachother 249 int avg1 = 8, 250 int avg2 = 8, 251 int OverlayWindowLeft = 70, 252 int OverlayWindowRight = 230 85 TString InRootFileName = "20111109_006.pulses.root", 86 TString InputPath = "../analysis/FPulseTemplate/20111109_006/", 87 TString OutputRootFileName = "20111109_006.pulses.root", 88 TString OutPutPath = "../analysis/FPulseTemplate/20111109_006/Templates/", 89 bool ProduceGraphic = false, 90 bool fitdata = false, 91 bool stats = false, 92 bool debugPixel = false, 93 int pixelSetSize = 40, 94 int maxPulseOrder = 3, 95 int refresh_rate = 500, //refresh rate for canvases 96 int verbosityLevel = 0, // different verbosity levels can be implemented here 97 int firstpixel = 0, 98 int npixel = -1 253 99 ) 254 100 { 255 gGainMean = GainMean; 256 gBSLMean = BSLMean; 257 //---------------------------------------------------------------------------- 258 // Save-Root-File Settings 259 //---------------------------------------------------------------------------- 260 CreateRootFile( OutRootFileName, verbosityLevel ); 261 101 102 //---------------------------------------------------------------------------- 103 // Open-Root-File Settings 104 //---------------------------------------------------------------------------- 105 TFile * inputRootFile = OpenRootFile( InputPath, InRootFileName, verbosityLevel ); 106 107 //---------------------------------------------------------------------------- 108 // global variable Settings 109 //---------------------------------------------------------------------------- 110 if ( npixel == -1 ) 111 { 112 npixel = NPIX - firstpixel; 113 } 114 115 116 if ( pixelSetSize == -1 ) 117 { 118 pixelSetSize = firstpixel +npixel; 119 } 120 // float GainMean = GainMean; // this has to be extracted from root files 121 // float BSLMean = BSLMean; // this has to be extracted from root files 262 122 263 123 //---------------------------------------------------------------------------- 264 124 // root global Settings 265 125 //---------------------------------------------------------------------------- 266 267 gPixelOverlayXaxisLeft = OverlayWindowLeft;268 gPixelOverlayXaxisRight = OverlayWindowRight;269 270 126 gStyle->SetPalette(1,0); 271 127 gROOT->SetStyle("Plain"); 272 // gPad->SetGrid(); 273 128 gPad->SetGrid(); 129 130 131 132 //---------------------------------------------------------------------------- 133 // root Canvas Settings 134 //---------------------------------------------------------------------------- 135 136 //Canvas Pad numbering 137 int pulsesCanvasFrameNrs[4] = { 138 1, // Top left 139 2, // Top right 140 3, // bottom left 141 4 // bootom right 142 }; 143 144 //Canvas Pad numbering 145 int distributionCanvasFrameNrs[4] = { 146 1, // Top left 147 2, // Top right 148 3, // bottom left 149 4 // bootom right 150 }; 151 152 if (ProduceGraphic) 153 { 154 155 //Canvases 156 cgpPixelPulses = new TCanvas*[maxPulseOrder]; 157 cgpDistributions = new TCanvas*[maxPulseOrder]; 158 159 //TCanvas* gpcDevelopment = NULL; 160 TString cName = ""; 161 TString cTitle = ""; 162 163 //naming of pulse canvases 274 164 for ( 275 int pulse_order = MAX_PULS_ORDER- 1;165 int pulse_order = maxPulseOrder - 1; 276 166 pulse_order >= 0 ; 277 167 pulse_order-- 278 168 ) 279 169 { 280 TString cName ="cgpPixelPulses"; 281 cName += pulse_order; 282 TString cTitle ="Pulses of Order: "; 283 cTitle += pulse_order; 284 cgpPixelPulses[pulse_order] = new TCanvas(cName,cTitle, 0,pulse_order*20,800,800); 170 cName ="cgpDistributions"; 171 cName += pulse_order; 172 cTitle ="Distributions of Pulses with Order of: "; 173 cTitle += pulse_order; 174 cgpDistributions[pulse_order] = new TCanvas(cName,cTitle, 720,pulse_order*20,720,720); 175 cgpDistributions[pulse_order]->Divide(2, 2); 176 cName ="cgpPixelPulses"; 177 cName += pulse_order; 178 cTitle ="Overlays of Pulses with Order of: "; 179 cTitle += pulse_order; 180 cgpPixelPulses[pulse_order] = new TCanvas(cName,cTitle, 0,pulse_order*20,720,720); 285 181 cgpPixelPulses[pulse_order]->Divide(2, 2); 286 cName = "cgpAllPixelPulses";287 cName += pulse_order;288 cTitle ="AllPixel average of puls shapes of Order: ";289 cTitle += pulse_order;290 cgpAllPixelPulses[pulse_order] = new TCanvas( cName, cTitle, 801, pulse_order*20, 800, 800 );291 cgpAllPixelPulses[pulse_order]->Divide(2, 2);292 182 } 183 293 184 294 185 // Create (pointer to) Canvases, which are used in every run, … … 297 188 // the pointers anyway ... 298 189 299 TCanvas *cFiltered = NULL; 300 if (spikeDebug) 301 { 302 cFiltered = new TCanvas("cFiltered","filtered DRS Waveforms", 1,310,400,300); 303 cFiltered->Divide(1, 3); 304 } 305 306 cgpTestHistos = new TCanvas( "cgpTestHistos", "Test Histograms", 801, 0, 800, 800 ); 307 cgpTestHistos->Divide(2,0); 190 // if (testmode) 191 // { 192 // //additional Test histograms 193 // cgpTestHistos = new TCanvas( "cgpTestHistos", "Test Histograms", 360, 420, 360, 360 ); 194 // cgpTestHistos->Divide(2,0); 195 // } 196 } 197 308 198 //----------------------------------------------------------------------------- 309 199 // Filter-Settings 310 200 //----------------------------------------------------------------------------- 311 // CFD filter settings 312 int k_cfd = 10; 313 vector<double> a_cfd(k_cfd, 0); 314 double b_cfd = 1.; 315 a_cfd[0]=-0.75; 316 a_cfd[k_cfd-1]=1.; 201 317 202 318 203 //----------------------------------------------------------------------------- … … 320 205 //----------------------------------------------------------------------------- 321 206 322 // Open the data file 323 324 fits * datafile; 325 // Opens the raw data file and 'binds' the variables given as 326 // Parameters to the data file. So they are filled with 327 // raw data as soon as datafile->GetRow(int) is called. 328 gNEvents = openDataFits( 329 datafilename, 330 &datafile, 331 AllPixelDataVector, 332 StartCellVector, 333 CurrentEventID, 334 gRegionOfInterest, 335 NumberOfPixels, 336 PXLxROI, 337 verbosityLevel 338 ); 339 340 if (gNEvents == 0) 341 { 342 cout << "return code of OpenDataFile:" << datafilename<< endl; 343 cout << "is zero -> aborting." << endl; 344 return 1; 345 } 346 347 if (verbosityLevel > 0) 348 { 349 cout << endl <<"number of events in file: "<< gNEvents << "\t"; 350 } 351 352 if ( nevents == -1 || nevents > gNEvents ) nevents = gNEvents; // -1 means all! 353 354 if (verbosityLevel > 0) 355 { 356 cout <<"of, which "<< nevents << " will be processed"<< endl; 357 cout <<"Total # of Pixel: "<< NumberOfPixels << "\t"; 358 } 359 360 if ( npixel == -1 || npixel > (int)NumberOfPixels ) npixel = (int)NumberOfPixels; // -1 means all! 361 362 if (verbosityLevel > 0) 363 { 364 cout <<"of, which "<< npixel << " will be processed"<< endl; 365 } 366 367 //Get the DRS calibration 368 RC = openCalibFits( 369 drsfilename, 370 drs_basemean, 371 drs_gainmean, 372 drs_triggeroffsetmean, 373 TriggerOffsetROI 374 ); 375 376 if (RC == 0) 377 { 378 cout << "return code of openCalibFits:" << drsfilename << endl; 379 cout << "is zero -> aborting." << endl; 380 return 1; 381 } 382 383 // Book the histograms 384 BookHistos(verbosityLevel ); 385 BookTestHistos( verbosityLevel ); 207 //---------------------------------------------------------------------------- 208 // Initialize Pixel 209 //---------------------------------------------------------------------------- 210 Pixel** pixel = new Pixel*[NPIX]; 211 212 for (int i = 0 ; i < NPIX; i++) 213 { 214 pixel[i] = NULL; 215 } 216 386 217 //------------------------------------- 387 218 // Loop over Pixel Sets 388 219 //------------------------------------- 389 for ( int firstPixelOfSample = firstpixel; 390 firstPixelOfSample < firstpixel + npixel; 391 firstPixelOfSample += sampleSetSize ) 392 { 393 220 for ( int firstPixelOfSet = firstpixel; 221 firstPixelOfSet < firstpixel + npixel; 222 firstPixelOfSet += pixelSetSize ) 223 { 394 224 if (verbosityLevel >= 0) 395 225 { 396 226 cout << "------------------------------------------------" << endl 397 << "...processing Sample fromPixel: "398 << firstPixelOfS ample;227 << "...processing Pixel: " 228 << firstPixelOfSet 399 229 << "to Pixel: " 400 << firstPixelOfS ample+sampleSetSize-1 << endl;230 << firstPixelOfSet+pixelSetSize-1 << endl; 401 231 } 402 BookPixelHistos(verbosityLevel ); 403 404 //-------------------------------------------------------------------- 405 // Loops over Every Event of Pixel 406 //-------------------------------------------------------------------- 407 for ( int ev = firstevent; ev < firstevent + nevents; ev++) 232 233 //-------------------------------------------------------------------- 234 // Loops over every Pixel of a Set of Pixels 235 //-------------------------------------------------------------------- 236 for ( int pixelID = firstPixelOfSet; 237 pixelID < firstPixelOfSet + pixelSetSize 238 && pixelID < firstpixel + npixel; 239 pixelID++ ) 408 240 { 409 // Get an Event --> consists of 1440 Pixel ...erm....data 410 datafile->GetRow( ev ); 411 412 if (verbosityLevel > 0) 241 242 if (verbosityLevel > 1) 413 243 { 414 244 cout << "-------------------------------------" << endl 415 << "...processing Event: " << CurrentEventID 416 << "/" << nevents << " of Pixel: " << pixel << endl; 245 << "...processing Set from Pixel " 246 << firstPixelOfSet 247 << " to Pixel " 248 << firstPixelOfSet+pixelSetSize-1 249 << " Pixel: " << pixelID 250 << "/" << firstpixel + npixel -1 << endl; 417 251 } 418 252 419 //-------------------------------------------------------------------- 420 // Loops over every Pixel of a Set of Pixels 421 //-------------------------------------------------------------------- 422 for ( int pixel = firstPixelOfSample; 423 pixel < sampleSetSize || pixel < firstpixel + npixel; 424 pixel++ ) 425 { 426 if (verbosityLevel > 0) 427 { 428 cout << "-------------------------------------" << endl 429 << "...processing Pixel: " << pixel 430 << "/" << nPixelSample << endl; 431 } 432 433 //------------------------------------- 434 // Apply Calibration 435 //------------------------------------- 436 if (verbosityLevel > 2) cout << "...applying DrsCalibration"; 437 applyDrsCalibration( 438 Ameas, 439 pixel, 440 12, 441 12, 442 drs_basemean, 443 drs_gainmean, 444 drs_triggeroffsetmean, 445 gRegionOfInterest, 446 AllPixelDataVector, 447 StartCellVector 253 //------------------------------------- 254 // Create individual Pixel 255 //------------------------------------- 256 pixel[pixelID] = new Pixel( 257 pixelID, 258 maxPulseOrder, ///One should get this from the file 259 verbosityLevel, 260 stats, 261 inputRootFile 262 ); 263 264 if (breakout) break; 265 266 //------------------------------------- 267 // Histogramms of Maximas in Overlay Spectra 268 //------------------------------------- 269 270 for ( int pulse_order = 0 ; 271 pulse_order < pixel[pixelID]->mMaxPulseOrder ; 272 pulse_order ++) 273 { 274 if (verbosityLevel > 2) 275 { 276 cout << "-------------------------------------" << endl 277 << "...processing Set from Pixel " 278 << firstPixelOfSet 279 << " to Pixel " 280 << firstPixelOfSet+pixelSetSize-1 281 << " Pixel: " << pixelID 282 << "/" << firstpixel + npixel -1 << endl 283 << " Pulse-Order: " << pulse_order; 284 } 285 286 // Calculate Max Prop. Value of each slice 287 //------------------------------------- 288 289 //from Maximum Overlay 290 ExtractPulseTemplate( 291 pixel[pixelID], 292 "Maximum", 293 pulse_order, 294 verbosityLevel 295 ); 296 297 //from Edge Overlay 298 ExtractPulseTemplate( 299 pixel[pixelID], 300 "Edge", 301 pulse_order, 302 verbosityLevel 303 ); 304 305 WritePixelTemplateToCsv( 306 pixel[pixelID], 307 OutPutPath, 308 "Maximum", 309 pulse_order, 310 verbosityLevel 311 ); 312 313 WritePixelTemplateToCsv( 314 pixel[pixelID], 315 OutPutPath, 316 "Edge", 317 pulse_order, 318 verbosityLevel 319 ); 320 321 pixel[pixelID]->DrawTemplateHistograms( 322 cgpPixelPulses, 323 pulsesCanvasFrameNrs 324 ); 325 326 pixel[pixelID]->DrawEdgeTemplateHistograms( 327 cgpPixelPulses, 328 pulsesCanvasFrameNrs 329 ); 330 331 332 333 } 334 // End of Loop over pulsorder of current Pixel 335 336 if (ProduceGraphic) 337 { 338 UpdateCanvases( 339 verbosityLevel, 340 MAX_PULS_ORDER, 341 false 448 342 ); 449 if (verbosityLevel > 2) cout << "...done " << endl; 450 451 //------------------------------------- 452 // Apply Filters 453 //------------------------------------- 454 // finds spikes in the raw data, and interpolates the value 455 // spikes are: 1 or 2 slice wide, positive non physical artifacts 456 if (verbosityLevel > 2) cout << "...removeing Spikes"; 457 removeSpikes (Ameas, Vcorr); 458 if (verbosityLevel > 2) cout << "...done " << endl; 459 460 // filter Vcorr with sliding average using FIR filter function 461 if (verbosityLevel > 2) cout << "...applying sliding average filter"; 462 sliding_avg(Vcorr, Vslide, avg1); 463 if (verbosityLevel > 2) cout << "...done " << endl; 464 465 // filter Vslide with CFD using FIR filter function 466 if (verbosityLevel > 2) cout << "...apllying factfir filter"; 467 factfir(b_cfd , a_cfd, k_cfd, Vslide, Vcfd); 468 if (verbosityLevel > 2) cout << "...done " << endl; 469 470 // filter Vcfd with sliding average using FIR filter function 471 if (verbosityLevel > 2) cout << "...applying 2nd sliding average filter"; 472 sliding_avg(Vcfd, Vcfd2, avg2); 473 if (verbosityLevel > 2) cout << "...done " << endl; 474 475 //------------------------------------- 476 // Search vor Zero crossings 477 //------------------------------------- 478 if (verbosityLevel > 2) cout << endl << "...searching zero crossings" ; 479 // peaks in Ameas[] are found by searching for zero crossings 480 // in Vcfd2 481 // first Argument 1 means ... *rising* edge 482 // second Argument 1 means ... search with stepsize 1 ... 10 is okay as well 483 vector<Region>* pZXings = zerosearch( Vcfd2 , 1 , 8); 484 // pZXings means "zero cross ings" 485 EnlargeRegion(*pZXings, 10, 10); 486 findAbsMaxInRegions(*pZXings, Vslide); 487 removeMaximaBelow( *pZXings, 3.0); 488 removeRegionWithMaxOnEdge( *pZXings, 2); 489 removeRegionOnFallingEdge( *pZXings, 100); 490 findTimeOfHalfMaxLeft(*pZXings, Vcorr, gBSLMean, 5, 10, verbosityLevel ); 491 if (verbosityLevel > 2) cout << "...done" << endl; 492 493 //----------------------------------------------------------------------------- 494 // Fill Overlay Histos 495 //----------------------------------------------------------------------------- 496 FillHistograms(pZXings, AmplWindowWidth, ev, verbosityLevel ); 497 498 //----------------------------------------------------------------------------- 499 // Spike Debug 500 //----------------------------------------------------------------------------- 501 if ( spikeDebug ) 502 { 503 // TODO do this correct. The vectors should be the rigt ones... this is just luck 504 debugHistos[Ameas_].GetXaxis()->Set(Ameas.size() , -0.5 , Ameas.size()-0.5); 505 debugHistos[Vcorr_].GetXaxis()->Set(Ameas.size() , -0.5 , Ameas.size()-0.5); 506 debugHistos[Vslide_].GetXaxis()->Set(Ameas.size() , -0.5 , Ameas.size()-0.5); 507 debugHistos[Vcfd_].GetXaxis()->Set(Ameas.size() , -0.5 , Ameas.size()-0.5); 508 debugHistos[Vcfd2_].GetXaxis()->Set(Ameas.size() , -0.5 , Ameas.size()-0.5); 509 510 for ( unsigned int sl = 0; sl < gRegionOfInterest; sl++) 343 } 344 345 if (debugPixel) 511 346 { 512 debugHistos[Ameas_].SetBinContent(sl, Ameas[sl]); 513 debugHistos[Vcorr_].SetBinContent(sl, Vcorr[sl]); 514 debugHistos[Vslide_].SetBinContent( sl, Vslide[sl] ); 515 debugHistos[Vcfd_].SetBinContent( sl, Vcfd[sl] ); 516 debugHistos[Vcfd2_].SetBinContent( sl, Vcfd2[sl] ); 347 //Process gui events asynchronously during input 348 cout << endl; 349 TTimer timer("gSystem->ProcessEvents();", 50, kFALSE); 350 timer.TurnOn(); 351 TString input = Getline("Type 'q' to exit, <return> to go on: "); 352 timer.TurnOff(); 353 if (input=="q\n") 354 { 355 break; 356 } 517 357 } 518 358 519 520 cFiltered->cd(1); 521 gPad->SetGrid(); 522 debugHistos[Ameas_].Draw(); 523 524 cFiltered->cd(2); 525 gPad->SetGrid(); 526 debugHistos[Vcorr_].Draw(); 527 528 cFiltered->cd(3); 529 gPad->SetGrid(); 530 debugHistos[Vslide_].Draw(); 531 532 TBox *OneBox; 533 vector<TBox*> MyBoxes; 534 for (unsigned int i=0; i<pZXings->size(); i++){ 535 OneBox = new TBox( 536 pZXings->at(i).maxPos -10 , 537 pZXings->at(i).maxVal -0.5, 538 pZXings->at(i).maxPos +10 , 539 pZXings->at(i).maxVal +0.5); 540 OneBox->SetLineColor(kBlue); 541 OneBox->SetLineWidth(1); 542 OneBox->SetFillStyle(0); 543 OneBox->SetFillColor(kRed); 544 MyBoxes.push_back(OneBox); 545 OneBox->Draw(); 546 } 547 548 // cFiltered->cd(3); 549 // gPad->SetGrid(); 550 // debugHistos[Vcfd2_].Draw(); 551 // TLine *zeroline = new TLine(0, 0, 1024, 0); 552 // zeroline->SetLineColor(kBlue); 553 // zeroline->Draw(); 554 555 cFiltered->Update(); 556 557 558 //Process gui events asynchronously during input 559 TTimer timer("gSystem->ProcessEvents();", 50, kFALSE); 560 timer.TurnOn(); 561 TString input = Getline("Type 'q' to exit, <return> to go on: "); 562 timer.TurnOff(); 563 if (input=="q\n") { 564 breakout=true; 565 } 566 567 //TODO!!!!!!!!! 568 // do some Garbage collection ... 569 // all the Objects on the heap should be deleted here. 570 571 }// end of if(spikeDebug) 572 573 delete pZXings; 574 if (breakout) break; 575 //------------------------------------- 576 // Draw 1. Set of Pixel Histograms 577 //------------------------------------- 578 if ((ev % refresh_rate) == 0) 359 if (verbosityLevel > 2) 579 360 { 580 DrawPulseHistograms( 581 verbosityLevel, 582 MAX_PULS_ORDER 583 ); 584 DrawTestHistograms( 585 verbosityLevel 586 ); 587 if (ProduceGraphic) 588 { 589 UpdateCanvases( 590 verbosityLevel, 591 MAX_PULS_ORDER 592 ); 593 } 594 } 595 596 if (breakout) break; 597 598 if (verbosityLevel > 2) 599 { 600 cout << endl << "-------------------------------------"<< endl; 601 } 602 603 604 605 606 }//End of Loop over Sample 607 //------------------------------------- 608 // end of Loops over Sample 609 //------------------------------------- 610 611 }//End of Loop over Events 612 //------------------------------------- 613 // end of Loops over Events 614 //------------------------------------- 615 616 617 618 619 //------------------------------------- 620 // Histogramms of Maximas in Overlay Spectra 621 //------------------------------------- 622 623 PlotPulseShapeOfMaxPropability( 624 MAX_PULS_ORDER, 625 fitdata, 626 verbosityLevel 627 ); 628 629 for (unsigned int pulse_order = 0 ; pulse_order < max_pulse_order ; pulse_order ++) 630 { 631 if (verbosityLevel > 2) cout << "\t...plotting most probable pulse of" 632 << "pulse order " << pulse_order; 633 PlotMedianEachSliceOfPulse( 634 hPixelOverlay[pulse_order], 635 hPixelMedian[pulse_order], 636 fitdata, 637 verbosityLevel 638 ); 639 640 PlotMeanEachSliceOfPulse( 641 hPixelOverlay[pulse_order], 642 hPixelMean[pulse_order], 643 fitdata, 644 verbosityLevel 645 ); 646 } 647 WritePixelTemplateToCsv( 648 OutPutPath, 649 "PulseTemplate_PointSet", 650 "Maximum", 651 pixel, 652 verbosityLevel 653 ); 654 655 if (verbosityLevel > 2) cout << "...done" << endl; 656 //here is what happends at the end of each loop over all pixels 657 //Save Histograms of Pixel into Output rootfile 658 DrawPulseHistograms( 659 verbosityLevel, 660 MAX_PULS_ORDER 661 ); 662 663 DrawMaximumHistograms( 664 verbosityLevel, 665 MAX_PULS_ORDER 666 ); 667 668 if (ProduceGraphic) 669 { 670 UpdateCanvases( 671 verbosityLevel, 672 MAX_PULS_ORDER 673 ); 674 } 675 676 // SaveHistograms( 677 // OutRootFileName, 678 // CreateSubDirName(pixel), 679 // hList, 680 // verbosityLevel 681 // ); 682 683 if (debugPixel) 684 { 685 //Process gui events asynchronously during input 686 cout << endl; 687 TTimer timer("gSystem->ProcessEvents();", 50, kFALSE); 688 timer.TurnOn(); 689 TString input = Getline("Type 'q' to exit, <return> to go on: "); 690 timer.TurnOff(); 691 if (input=="q\n") { 692 break; 693 } 694 } 695 696 DeletePixelHistos(verbosityLevel); 697 698 if (verbosityLevel > 2) 699 { 700 cout << endl << "...End of Pixel" 361 cout << endl << "...End of Pixel" 701 362 << endl << "------------------------------------------------" 702 363 << endl; 703 } 704 } 705 // End of Loop over all Pixels 364 } 365 } 366 // End of Loop over all Pixels of set 367 368 if (verbosityLevel > 1) 369 { 370 cout << endl << "...End of Loop over all Pixels of set" 371 << endl << "------------------------------------------------" 372 << endl; 373 } 374 } 375 // End of Loop over all Pixelsets 376 377 if (verbosityLevel > 0) 378 { 379 cout << endl << "...End of Loop over all Pixelsets" 380 << endl << "------------------------------------------------" 381 << endl; 382 } 706 383 707 384 //------------------------------------- 708 // Draw Histograms385 // Draw All Pixel Histograms 709 386 //------------------------------------- 710 387 711 SaveHistograms( //save histograms of all pixel into output root file712 OutRootFileName,713 CreateSubDirName("All"),714 hAllPixelList,715 verbosityLevel716 );388 // SaveHistograms( //save histograms of all pixel into output root file 389 // OutInRootFileName, 390 // CreateSubDirName("All"), 391 // hAllPixelList, 392 // verbosityLevel 393 // ); 717 394 718 395 // SaveHistograms( //save histograms of generell results into output root file 719 // Out RootFileName,396 // OutInRootFileName, 720 397 // "root", 721 398 // hRootList, … … 725 402 726 403 727 if (ProduceGraphic)728 {729 UpdateCanvases(730 verbosityLevel,731 MAX_PULS_ORDER732 );733 }734 735 WriteAllPixelTemplateToCsv(736 OutPutPath,737 "PulseTemplate_PointSet",738 "Maximum",739 verbosityLevel740 );741 742 DeletePixelCanvases( verbosityLevel );404 // if (ProduceGraphic) 405 // { 406 // UpdateCanvases( 407 // verbosityLevel, 408 // MAX_PULS_ORDER 409 // ); 410 // } 411 412 // WriteAllPixelTemplateToCsv( 413 // OutPutPath, 414 // "PulseTemplate_PointSet", 415 // "Maximum", 416 // verbosityLevel 417 // ); 418 419 DeletePixelCanvases( maxPulseOrder ,verbosityLevel ); 743 420 return( 0 ); 744 421 } … … 753 430 // Funktions 754 431 //----------------------------------------------------------------------------- 755 756 /*757 ich erzeuge sowas wie 36 pixel mit ca 20 Histogrammen pro pixel758 wie viel speicher braucht das?759 Ich brauche eine möglichkeit jedes dieser histogramme zu identifizieren760 am besten ein array oder eine liste oder einen abstracten datentyp mit den histogrammen als attribute.761 ne klasse wäre super762 763 */764 432 void 765 BookPixelHistos( 766 int sampleSize, 767 768 int verbosityLevel 769 ) 770 { 771 if (verbosityLevel > 2) cout << endl << "...book pixel histograms" << endl; 772 hList = new TObjArray; 773 774 TString histoname; 775 for (int order = 0; order < MAX_PULS_ORDER; order ++) 776 { 777 // SetHistogramAttributes( 778 // "PeakOverlay", 779 // order, 780 // gOverlayAttrib, 781 // MaxAmplOfFirstPulse 782 // ); 783 histoname = "hPixelOverlay"; 784 histoname += order; 785 if (verbosityLevel > 3) cout << "...booking " << histoname << endl; 786 hPixelOverlay[order] = new TH2F( 787 histoname, 788 "Overlay of detected pulses of one pulse order for one Pixel", 789 gPixelOverlayXaxisLeft + gPixelOverlayXaxisRight , 790 (-1*gPixelOverlayXaxisLeft)-0.5, 791 gPixelOverlayXaxisRight-0.5 , 792 512, 793 -55.5, 794 200.5 795 ); 796 hPixelOverlay[order]->SetAxisRange( 797 gBSLMean - 5, 798 (gGainMean*(order+1)) + 10, 799 "Y"); 800 hPixelOverlay[order]->GetXaxis()->SetTitle( "Timeslices [a.u.]" ); 801 hPixelOverlay[order]->GetYaxis()->SetTitle( "Amplitude [mV]" ); 802 //hPixelProfile->SetBit(TH2F::kCanRebin); 803 hList->Add( hPixelOverlay[order] ); 804 805 //------------------------------------------------------------------------ 806 // SetHistogramAttributes( 807 // "PeakMaximum", 808 // order, 809 // gMaximumAttrib, 810 // MaxAmplOfFirstPulse 811 // ); 812 histoname = "hPixelMax"; 813 histoname += order; 814 if (verbosityLevel > 3) cout << "...booking " << histoname << endl; 815 hPixelMax[order] = new TH1F( 816 histoname, 817 "Maximum value of each slice in overlay plot", 818 gPixelOverlayXaxisLeft + gPixelOverlayXaxisRight , 819 (-1*gPixelOverlayXaxisLeft)-0.5, 820 gPixelOverlayXaxisRight-0.5 821 ); 822 hPixelMax[order]->SetAxisRange( 823 gBSLMean - 5, 824 (gGainMean*(order+1)) + 10, 825 "Y"); 826 hPixelMax[order]->GetXaxis()->SetTitle( "Timeslices [a.u.]" ); 827 hPixelMax[order]->GetYaxis()->SetTitle( "Amplitude [mV]" ); 828 hList->Add( hPixelMax[order] ); 829 830 //------------------------------------------------------------------------ 831 // SetHistogramAttributes( 832 // "PeakMaximum", 833 // order, 834 // gMaximumAttrib, 835 // MaxAmplOfFirstPulse 836 // ); 837 histoname = "hPixelMedian"; 838 histoname += order; 839 if (verbosityLevel > 3) cout << "...booking " << histoname << endl; 840 hPixelMedian[order] = new TH1F( 841 histoname, 842 "Median of each slice in overlay plot", 843 gPixelOverlayXaxisLeft + gPixelOverlayXaxisRight , 844 (-1*gPixelOverlayXaxisLeft)-0.5, 845 gPixelOverlayXaxisRight-0.5 846 ); 847 hPixelMedian[order]->SetAxisRange( 848 gBSLMean - 5, 849 (gGainMean*(order+1)) + 10, 850 "Y"); 851 hPixelMedian[order]->SetLineColor(kRed); 852 hPixelMedian[order]->GetXaxis()->SetTitle( "Timeslices [a.u.]" ); 853 hPixelMedian[order]->GetYaxis()->SetTitle( "Amplitude [mV]" ); 854 hList->Add( hPixelMedian[order] ); 855 856 //------------------------------------------------------------------------ 857 // SetHistogramAttributes( 858 // "PeakMaximum", 859 // order, 860 // gMaximumAttrib, 861 // MaxAmplOfFirstPulse 862 // ); 863 histoname = "hPixelMean"; 864 histoname += order; 865 if (verbosityLevel > 3) cout << "...booking " << histoname << endl; 866 hPixelMean[order] = new TH1F( 867 histoname, 868 "Mean of each slice in overlay plot", 869 gPixelOverlayXaxisLeft + gPixelOverlayXaxisRight , 870 (-1*gPixelOverlayXaxisLeft)-0.5, 871 gPixelOverlayXaxisRight-0.5 872 ); 873 hPixelMean[order]->SetAxisRange( 874 gBSLMean - 5, 875 (gGainMean*(order+1)) + 10, 876 "Y"); 877 hPixelMean[order]->SetLineColor(kBlue); 878 hPixelMean[order]->GetXaxis()->SetTitle( "Timeslices [a.u.]" ); 879 hPixelMean[order]->GetYaxis()->SetTitle( "Amplitude [mV]" ); 880 hList->Add( hPixelMean[order] ); 881 882 //------------------------------------------------------------------------ 883 // SetHistogramAttributes( 884 // "PeakMaxGaus", 885 // order, 886 // gMaxGausAttrib, 887 // MaxAmplOfFirstPulse 888 // ); 889 histoname = "hPixelEdgeOverlay"; 890 histoname += order; 891 if (verbosityLevel > 3) cout << "...booking " << histoname << endl; 892 hPixelEdgeOverlay[order] = new TH2F( 893 histoname, 894 "Overlay at rising edge of detected pulses of one pulse order for one Pixel ", 895 gPixelOverlayXaxisLeft + gPixelOverlayXaxisRight , 896 (-1*gPixelOverlayXaxisLeft)-0.5, 897 gPixelOverlayXaxisRight-0.5 , 898 512, 899 -55.5, 900 200.5 901 ); 902 903 hPixelEdgeOverlay[order]->SetAxisRange( 904 gBSLMean - 5, 905 (gGainMean*(order+1)) + 10, 906 "Y"); 907 hPixelEdgeOverlay[order]->GetXaxis()->SetTitle( "Timeslices [a.u.]" ); 908 hPixelEdgeOverlay[order]->GetYaxis()->SetTitle( "Amplitude [mV]" ); 909 hList->Add( hPixelEdgeOverlay[order] ); 910 911 //------------------------------------------------------------------------ 912 // SetHistogramAttributes( 913 // "PeakProfile", 914 // order, 915 // gProfileAttrib, 916 // MaxAmplOfFirstPulse 917 // ); 918 histoname = "hPixelProfile"; 919 histoname += order; 920 if (verbosityLevel > 3) cout << "...booking " << histoname << endl; 921 hPixelProfile[order] = new TProfile( 922 histoname, 923 "Mean value of each slice in overlay plot (Tprofile)", 924 gPixelOverlayXaxisLeft + gPixelOverlayXaxisRight ,//nbinsx 925 (-1*gPixelOverlayXaxisLeft)-0.5, //xlow 926 gPixelOverlayXaxisRight-0.5 , //xup 927 // 512, //nbinsy 928 -55.5, //ylow 929 300.5, //yup 930 "s"); //option 931 hPixelProfile[order]->SetAxisRange( 932 gBSLMean - 5, 933 (gGainMean*(order+1)) + 10, 934 "Y"); 935 hPixelProfile[order]->GetXaxis()->SetTitle( "Timeslices [a.u.]" ); 936 hPixelProfile[order]->GetYaxis()->SetTitle( "Amplitude [mV]" ); 937 //hPixelProfile->SetBit(TH2F::kCanRebin); 938 hList->Add( hPixelProfile[order] ); 939 940 941 //------------------------------------------------------------------------ 942 // SetHistogramAttributes( 943 // "PeakProfile", 944 // order, 945 // gProfileAttrib, 946 // MaxAmplOfFirstPulse 947 // ); 948 histoname = "hPixelProfile2"; 949 histoname += order; 950 if (verbosityLevel > 3) cout << "...booking " << histoname << endl; 951 hPixelProfile2[order] = new TProfile( 952 histoname, 953 "Mean value of each slice in overlay plot (Tprofile)", 954 gPixelOverlayXaxisLeft + gPixelOverlayXaxisRight ,//nbinsx 955 (-1*gPixelOverlayXaxisLeft)-0.5, //xlow 956 gPixelOverlayXaxisRight-0.5 , //xup 957 // 512, //nbinsy 958 -55.5, //ylow 959 300.5, //yup 960 "s"); //option 961 hPixelProfile2[order]->SetLineColor(kRed); 962 hPixelProfile2[order]->SetAxisRange( 963 gBSLMean - 5, 964 (gGainMean*(order+1)) + 10, 965 "Y"); 966 hPixelProfile2[order]->GetXaxis()->SetTitle( "Timeslices [a.u.]" ); 967 hPixelProfile2[order]->GetYaxis()->SetTitle( "Amplitude [mV]" ); 968 //hPixelProfile->SetBit(TH2F::kCanRebin); 969 970 971 972 973 } 974 if (verbosityLevel > 2) cout << "...done" << endl; 975 } 976 //end of BookPixelHistos 977 //---------------------------------------------------------------------------- 978 979 980 void DeletePixelHistos( int verbosityLevel ) 981 { 982 if (verbosityLevel > 2) cout << endl << "...delete current pixel histograms" ; 983 for (int order = 0; order < MAX_PULS_ORDER; order ++) 984 { 985 if (verbosityLevel > 3) cout << endl << "...deleting hPixelOverlay" << order; 986 delete hPixelOverlay[order]; 987 if (verbosityLevel > 3) cout << endl << "...deleting hPixelMax" << order; 988 delete hPixelMax[order]; 989 if (verbosityLevel > 3) cout << endl << "...deleting hPixelMedian" << order; 990 delete hPixelMedian[order]; 991 if (verbosityLevel > 3) cout << endl << "...deleting hPixelMean" << order; 992 delete hPixelMean[order]; 993 if (verbosityLevel > 3) cout << endl << "...deleting hPixelEdgeOverlay" << order; 994 delete hPixelEdgeOverlay[order]; 995 if (verbosityLevel > 3) cout << endl << "...deleting hPixelProfile" << order; 996 delete hPixelProfile[order]; 997 if (verbosityLevel > 3) cout << endl << "...deleting hPixelProfile2" << order; 998 delete hPixelProfile2[order]; 999 } 1000 if (verbosityLevel > 3) cout << endl << "...deleting hList"; 1001 delete hList; 1002 if (verbosityLevel > 2) cout << endl << "...done" << endl; 1003 } 1004 //end of DeletePixelHistos 1005 //---------------------------------------------------------------------------- 1006 1007 void BookTestHistos( int verbosityLevel ) 1008 { 1009 if (verbosityLevel > 2) cout << endl << "...book pixel histograms" << endl; 1010 1011 hTesthisto = new TH1F ( 1012 "hTesthisto", 1013 "Deviation of rising edge and maximum", 1014 600, 1015 -10.1, 1016 10.1 1017 ); 1018 hTesthisto->GetXaxis()->SetTitle( "Timeslices [a.u.]" ); 1019 hTesthisto->GetYaxis()->SetTitle( "counts" ); 1020 1021 hTesthisto2 = new TH2F ( 1022 "hTesthisto2", 1023 "Deviation of rising edge and maximum by event #", 1024 gNEvents, 1025 250, 1026 800, 1027 600, 1028 -10.1, 1029 10.1 1030 ); 1031 // hTesthisto2->GetXaxis()->SetTitle( "Timeslices [a.u.]" ); 1032 // hTesthisto2->GetYaxis()->SetTitle( "counts" ); 1033 // hTesthisto2->SetMarkerStyle( 2 ); 1034 1035 if (verbosityLevel > 2) cout << "...done" << endl; 1036 } 1037 //end of BookTestHistos 1038 //---------------------------------------------------------------------------- 1039 1040 void BookHistos( int verbosityLevel ) 1041 { 1042 if (verbosityLevel > 2) cout << endl << "...book histograms" << endl; 1043 hAllPixelList = new TObjArray; 1044 hRootList = new TObjArray; 1045 debugHistos = new TH1F[ Number_Of_Debug_Histo_Types ]; 1046 for ( int type = 0; type < Number_Of_Debug_Histo_Types; type++){ 1047 debugHistos[ type ].SetBins(1024, 0, 1024); 1048 debugHistos[ type ].SetLineColor(1); 1049 debugHistos[ type ].SetLineWidth(2); 1050 debugHistos[ type ].SetStats(false); 1051 1052 // set X axis paras 1053 debugHistos[ type ].GetXaxis()->SetLabelSize(0.08); 1054 debugHistos[ type ].GetXaxis()->SetTitleSize(0.08); 1055 debugHistos[ type ].GetXaxis()->SetTitleOffset(1.0); 1056 debugHistos[ type ].GetXaxis()->SetTitle(Form("Time slice [%.1f ns/slice]", 1./2.)); 1057 1058 // set Y axis paras 1059 debugHistos[ type ].GetYaxis()->SetLabelSize(0.08); 1060 debugHistos[ type ].GetYaxis()->SetTitleSize(0.08); 1061 debugHistos[ type ].GetYaxis()->SetTitleOffset(0.3); 1062 debugHistos[ type ].GetYaxis()->SetTitle("Amplitude [mV]"); 1063 } 1064 1065 // All Pixel Histograms 1066 //------------------------------------------------------------------------ 1067 TString histoname; 1068 for (int order = 0; order < MAX_PULS_ORDER; order ++) 1069 { 1070 // SetHistogramAttributes( 1071 // "AllPixelPeakOverlay", 1072 // order, 1073 // gOverlayAttrib, 1074 // MaxAmplOfFirstPulse 1075 // ); 1076 histoname = "hAllPixelOverlay"; 1077 histoname += order; 1078 if (verbosityLevel > 3) cout << "...booking " << histoname << endl; 1079 hAllPixelOverlay[order] = new TH2F( 1080 histoname, 1081 "Overlay of detected Pulses of one Order for All Pixels", 1082 gPixelOverlayXaxisLeft + gPixelOverlayXaxisRight , 1083 (-1*gPixelOverlayXaxisLeft)-0.5, 1084 gPixelOverlayXaxisRight-0.5 , 1085 512, 1086 -55.5, 1087 300.5 1088 ); 1089 hAllPixelOverlay[order]->SetAxisRange( 1090 gBSLMean - 5, 1091 (gGainMean*(order+1)) + 10, 1092 "Y"); 1093 hAllPixelOverlay[order]->GetXaxis()->SetTitle( "Timeslices [a.u.]" ); 1094 hAllPixelOverlay[order]->GetYaxis()->SetTitle( "Amplitude [mV]" ); 1095 //hAllPixelOverlay->SetBit(TH2F::kCanRebin); 1096 hAllPixelList->Add( hAllPixelOverlay[order] ); 1097 1098 //------------------------------------------------------------------------ 1099 1100 // SetHistogramAttributes( 1101 // "AllPixelPeakMax", 1102 // order, 1103 // gMaximumAttrib, 1104 // MaxAmplOfFirstPulse 1105 // ); 1106 histoname = "hAllPixelMax"; 1107 histoname += order; 1108 if (verbosityLevel > 3) cout << "...booking " << histoname << endl; 1109 hAllPixelMax[order] = new TH1F( 1110 histoname, 1111 "Maximum value of each slice in overlay plot for All Pixels", 1112 gPixelOverlayXaxisLeft + gPixelOverlayXaxisRight , 1113 (-1*gPixelOverlayXaxisLeft)-0.5, 1114 gPixelOverlayXaxisRight-0.5 1115 ); 1116 hAllPixelMax[order]->SetAxisRange( 1117 gBSLMean - 5, 1118 (gGainMean*(order+1)) + 10, 1119 "Y"); 1120 hAllPixelMax[order]->GetXaxis()->SetTitle( "Timeslices [a.u.]" ); 1121 hAllPixelMax[order]->GetYaxis()->SetTitle( "Amplitude [mV]" ); 1122 hAllPixelList->Add( hAllPixelMax[order] ); 1123 1124 //------------------------------------------------------------------------ 1125 // SetHistogramAttributes( 1126 // "PeakMaximum", 1127 // order, 1128 // gMaximumAttrib, 1129 // MaxAmplOfFirstPulse 1130 // ); 1131 histoname = "hAllPixelMedian"; 1132 histoname += order; 1133 if (verbosityLevel > 3) cout << "...booking " << histoname << endl; 1134 hAllPixelMedian[order] = new TH1F( 1135 histoname, 1136 "Median of each slice in overlay plot for All Pixels", 1137 gPixelOverlayXaxisLeft + gPixelOverlayXaxisRight , 1138 (-1*gPixelOverlayXaxisLeft)-0.5, 1139 gPixelOverlayXaxisRight-0.5 1140 ); 1141 hAllPixelMedian[order]->SetAxisRange( 1142 gBSLMean - 5, 1143 (gGainMean*(order+1)) + 10, 1144 "Y"); 1145 hAllPixelMedian[order]->SetLineColor(kRed); 1146 hAllPixelMedian[order]->GetXaxis()->SetTitle( "Timeslices [a.u.]" ); 1147 hAllPixelMedian[order]->GetYaxis()->SetTitle( "Amplitude [mV]" ); 1148 hAllPixelList->Add( hAllPixelMedian[order] ); 1149 if (verbosityLevel > 3) cout << "...BREAKPOINT in " << histoname << endl; 1150 //------------------------------------------------------------------------ 1151 // SetHistogramAttributes( 1152 // "PeakMaximum", 1153 // order, 1154 // gMaximumAttrib, 1155 // MaxAmplOfFirstPulse 1156 // ); 1157 histoname = "hAllPixelMean"; 1158 histoname += order; 1159 if (verbosityLevel > 3) cout << "...booking " << histoname << endl; 1160 hAllPixelMean[order] = new TH1F( 1161 histoname, 1162 "Mean of each slice in overlay plot for All Pixels", 1163 gPixelOverlayXaxisLeft + gPixelOverlayXaxisRight , 1164 (-1*gPixelOverlayXaxisLeft)-0.5, 1165 gPixelOverlayXaxisRight-0.5 1166 ); 1167 hAllPixelMean[order]->SetAxisRange( 1168 gBSLMean - 5, 1169 (gGainMean*(order+1)) + 10, 1170 "Y"); 1171 hAllPixelMean[order]->SetLineColor(kBlue); 1172 hAllPixelMean[order]->GetXaxis()->SetTitle( "Timeslices [a.u.]" ); 1173 hAllPixelMean[order]->GetYaxis()->SetTitle( "Amplitude [mV]" ); 1174 hAllPixelList->Add( hAllPixelMean[order] ); 1175 1176 //------------------------------------------------------------------------ 1177 // SetHistogramAttributes( 1178 // "PeakMaxGaus", 1179 // order, 1180 // gMaxGausAttrib, 1181 // MaxAmplOfFirstPulse 1182 // ); 1183 histoname = "hAllPixelEdgeOverlay"; 1184 histoname += order; 1185 if (verbosityLevel > 3) cout << "...booking " << histoname << endl; 1186 hAllPixelEdgeOverlay[order] = new TH2F( 1187 histoname, 1188 "Overlay at rising edge of detected pulses of one pulse order for All Pixels ", 1189 gPixelOverlayXaxisLeft + gPixelOverlayXaxisRight , 1190 (-1*gPixelOverlayXaxisLeft)-0.5, 1191 gPixelOverlayXaxisRight-0.5 , 1192 512, 1193 -55.5, 1194 200.5 1195 ); 1196 1197 hAllPixelEdgeOverlay[order]->SetAxisRange( 1198 gBSLMean - 5, 1199 (gGainMean*(order+1)) + 10, 1200 "Y"); 1201 hAllPixelEdgeOverlay[order]->GetXaxis()->SetTitle( "Timeslices [a.u.]" ); 1202 hAllPixelEdgeOverlay[order]->GetYaxis()->SetTitle( "Amplitude [mV]" ); 1203 hAllPixelList->Add( hAllPixelEdgeOverlay[order] ); 1204 1205 //------------------------------------------------------------------------ 1206 // SetHistogramAttributes( 1207 // "AllPixelPeakProfile", 1208 // order, 1209 // gProfileAttrib, 1210 // MaxAmplOfFirstPulse 1211 // ); 1212 histoname = "hAllPixelProfile"; 1213 histoname += order; 1214 if (verbosityLevel > 3) cout << "...booking " << histoname << endl; 1215 hAllPixelProfile[order] = new TProfile( 1216 histoname, 1217 "Mean value of each slice in overlay plot for All Pixels (Tprofile)", 1218 gPixelOverlayXaxisLeft + gPixelOverlayXaxisRight ,//nbinsx 1219 (-1*gPixelOverlayXaxisLeft)-0.5, //xlow 1220 gPixelOverlayXaxisRight-0.5 , //xup 1221 // 512, //nbinsy 1222 -55.5, //ylow 1223 300.5, //yup 1224 "s"); //option 1225 hAllPixelProfile[order]->SetAxisRange( 1226 gBSLMean - 5, 1227 (gGainMean*(order+1)) + 10, 1228 "Y"); 1229 hAllPixelProfile[order]->GetXaxis()->SetTitle( "Timeslices [a.u.]" ); 1230 hAllPixelProfile[order]->GetYaxis()->SetTitle( "Amplitude [mV]" ); 1231 //hPixelProfile->SetBit(TH2F::kCanRebin); 1232 hAllPixelList->Add( hAllPixelProfile[order] ); 1233 1234 //------------------------------------------------------------------------ 1235 // SetHistogramAttributes( 1236 // "AllPixelPeakProfile", 1237 // order, 1238 // gProfileAttrib, 1239 // MaxAmplOfFirstPulse 1240 // ); 1241 histoname = "hAllPixelProfile2"; 1242 histoname += order; 1243 if (verbosityLevel > 3) cout << "...booking " << histoname << endl; 1244 hAllPixelProfile2[order] = new TProfile( 1245 histoname, 1246 "Mean value of each slice in overlay plot for All Pixels (Tprofile)", 1247 gPixelOverlayXaxisLeft + gPixelOverlayXaxisRight ,//nbinsx 1248 (-1*gPixelOverlayXaxisLeft)-0.5, //xlow 1249 gPixelOverlayXaxisRight-0.5 , //xup 1250 // 512, //nbinsy 1251 -55.5, //ylow 1252 300.5, //yup 1253 "s"); //option 1254 hAllPixelProfile2[order]->SetAxisRange( 1255 gBSLMean - 5, 1256 (gGainMean*(order+1)) + 10, 1257 "Y"); 1258 hAllPixelProfile2[order]->SetLineColor(kRed); 1259 hAllPixelProfile2[order]->GetXaxis()->SetTitle( "Timeslices [a.u.]" ); 1260 hAllPixelProfile2[order]->GetYaxis()->SetTitle( "Amplitude [mV]" ); 1261 //hPixelProfile->SetBit(TH2F::kCanRebin); 1262 hAllPixelList->Add( hAllPixelProfile2[order] ); 1263 1264 1265 1266 }//end of for over order 1267 if (verbosityLevel > 2) cout << "...done" << endl; 1268 } 1269 // end of BookHistos 1270 //---------------------------------------------------------------------------- 1271 1272 1273 void CreateRootFile(const char * loc_fname, int verbosityLevel) 1274 { 1275 TFile* tf = new TFile(loc_fname,"UPDATE"); 1276 if (tf->IsZombie()) 1277 { 1278 cout << "Error opening file" << endl; 1279 exit(-1); 1280 } 1281 else { 1282 if (verbosityLevel > 1) cout << "creating root-file successfull" << endl; 1283 tf->Close(); 1284 } 1285 } 1286 //end of CreateRootFile 1287 //---------------------------------------------------------------------------- 1288 1289 1290 TFile* OpenRootFile(const char * loc_fname, int verbosityLevel) 1291 { 1292 TFile* tf = new TFile(loc_fname,"UPDATE"); 1293 if (tf->IsZombie()) 1294 { 1295 cout << "Error opening file" << endl; 1296 exit(-1); 1297 } 1298 else { 1299 if (verbosityLevel > 1) cout << "...opening root-file:" 1300 << loc_fname 1301 << " successfull" << endl; 1302 } 1303 return tf; 1304 } 1305 //end of OpenRootFile 1306 //---------------------------------------------------------------------------- 1307 1308 1309 void CloseRootFile(TFile* tf) 1310 { 1311 if (tf->IsZombie()) 1312 { 1313 cout << "Error closing file" << endl; 1314 exit(-1); 1315 } else { 1316 tf->Close(); 1317 } 1318 } 1319 //end of CloseRootFile 1320 //---------------------------------------------------------------------------- 1321 1322 1323 void SaveHistograms( 1324 const char* loc_fname, 1325 const char* subdirectory, 1326 TObjArray* histList, 433 DeletePixelCanvases( 434 int maxPulseOrder, 1327 435 int verbosityLevel 1328 436 ) 1329 437 { 1330 if (verbosityLevel > 2) cout << endl 1331 << "...saving pixel histograms to subdirectory: " 1332 << subdirectory << endl ; 1333 1334 TFile* tf = OpenRootFile(loc_fname, verbosityLevel); 1335 if ( subdirectory != "root") 1336 { 1337 tf->cd(); 1338 tf->mkdir(subdirectory,subdirectory); 1339 tf->cd(subdirectory); 1340 } 1341 else 1342 { 1343 tf->cd(); 1344 } 1345 histList->Write(); // write the major histograms into the top level directory 1346 if (verbosityLevel > 3) tf->ls(); 1347 CloseRootFile( tf ); // close the file 1348 if (verbosityLevel > 2) cout << "...done" << endl; 438 if (verbosityLevel > 2) cout << endl << "...delete pixel Canvases" ; 439 for (int pulse_order = 0; pulse_order < maxPulseOrder; pulse_order++ ) 440 { 441 delete cgpPixelPulses[pulse_order]; 442 delete cgpDistributions[pulse_order]; 443 } 444 delete[] cgpPixelPulses; 445 delete[] cgpDistributions; 1349 446 } 1350 // end of SaveHistograms 1351 //---------------------------------------------------------------------------- 1352 1353 1354 TString CreateSubDirName(int pixel) 1355 { 1356 TString path = "Pixel_"; 1357 path += pixel; 1358 // cout << "path " << path << endl ; 1359 return path; 1360 } 1361 // end of CreateSubDirName 1362 //---------------------------------------------------------------------------- 1363 1364 TString CreateSubDirName(const char* title) 1365 { 1366 TString path = title; 1367 path += "_Pixel"; 1368 // cout << "path " << path << endl ; 1369 return path; 1370 } 1371 // end of CreateSubDirName 1372 //---------------------------------------------------------------------------- 1373 1374 void FillHistograms( 1375 vector<Region>* pZXings, 1376 int AmplWindowWidth, 1377 int eventNumber, 1378 int verbosityLevel 1379 ) 1380 { 1381 if (verbosityLevel > 2) cout << endl << "...filling pulse histograms" ; 1382 bool use_this_peak=false; 1383 int order_of_pulse=0; 1384 vector<Region>::iterator reg; 1385 1386 //Loop over all found zerocrossings in Region 1387 for (reg = pZXings->begin() ; reg < pZXings->end() ; reg++) 1388 { 1389 //skip those who are at the Rim of the ROI 1390 if (reg->maxPos < 12 || (unsigned int) reg->maxPos > gRegionOfInterest-12) 1391 { 1392 if (verbosityLevel > 2) cout << endl << "\t...out of range" << endl; 1393 continue; 1394 } 1395 //domstest->Fill(reg->maxPos); 1396 1397 // Define axis range of Histogram 1398 int Left = reg->maxPos - gPixelOverlayXaxisLeft; 1399 int Right = reg->maxPos + gPixelOverlayXaxisRight; 1400 int EdgeLeft = reg->halfRisingEdgePos - gPixelOverlayXaxisLeft; 1401 int EdgeRight = reg->halfRisingEdgePos + gPixelOverlayXaxisRight; 1402 1403 if (Left < 0) Left =0; 1404 if (Right > (int)Ameas.size() ) Right=Ameas.size() ; 1405 if (EdgeLeft < 0) EdgeLeft =0; 1406 if (EdgeRight > (int)Ameas.size() ) EdgeRight=Ameas.size() ; 1407 1408 1409 hTesthisto->Fill( reg->slopeOfRisingEdge ) ; 1410 hTesthisto2->Fill( eventNumber, reg->slopeOfRisingEdge ) ; 1411 1412 if (verbosityLevel > 2) cout << endl << "\t...choosing Histogram" ; 1413 1414 //determine order of pulse and which histogram shall be filled 1415 if (verbosityLevel > 2) cout << endl << "\t...choosing Histogram" ; 1416 for(int order = 0; order < MAX_PULS_ORDER; order++) 1417 { 1418 if (Ameas[reg->maxPos] >= ((gGainMean*(order+1)) - (AmplWindowWidth/2)) 1419 && Ameas[reg->maxPos] < ((gGainMean*(order+1)) + AmplWindowWidth/2)) 1420 { 1421 //hOverlayTemp = hPixelOverlay[order]; 1422 if (verbosityLevel > 2) cout << "...#" << order ; 1423 order_of_pulse = order; 1424 use_this_peak = true; 1425 break; 1426 } 1427 else if (order >= (MAX_PULS_ORDER - 1)){ 1428 use_this_peak = false; 1429 if (verbosityLevel > 2) cout << "...NONE" << endl ; 1430 } 1431 } 1432 1433 //Fill overlay und profile histograms 1434 if (use_this_peak){ 1435 if (verbosityLevel > 2) cout << "...filling" ; 1436 for ( int pos = Left; pos < Right; pos++){ 1437 // if (); 1438 hPixelOverlay[order_of_pulse]->Fill( pos - (reg->maxPos), Ameas[pos]) ; 1439 hPixelProfile[order_of_pulse]->Fill( pos - (reg->maxPos), Ameas[pos]) ; 1440 hAllPixelOverlay[order_of_pulse]->Fill( pos - (reg->maxPos), Ameas[pos]) ; 1441 hAllPixelProfile[order_of_pulse]->Fill( pos - (reg->maxPos), Ameas[pos]) ; 1442 } 1443 for ( int pos = EdgeLeft; pos < EdgeRight; pos++){ 1444 // cout << endl << "###filling edge histo###" << endl; 1445 hPixelEdgeOverlay[order_of_pulse]->Fill( pos - (reg->halfRisingEdgePos), Ameas[pos]) ; 1446 hPixelProfile2[order_of_pulse]->Fill( pos - (reg->halfRisingEdgePos), Ameas[pos]) ; 1447 hAllPixelEdgeOverlay[order_of_pulse]->Fill( pos - (reg->halfRisingEdgePos), Ameas[pos]) ; 1448 hAllPixelProfile2[order_of_pulse]->Fill( pos - (reg->halfRisingEdgePos), Ameas[pos]) ; 1449 // cout << endl << "######" << endl; 1450 } 1451 if (verbosityLevel > 2) cout << "...done" << endl; 1452 } 1453 //Breakout option 1454 if (breakout) 1455 break; 1456 } 1457 if (verbosityLevel > 2) cout << "...done" << endl; 1458 } 1459 // end of FillHistograms 1460 //---------------------------------------------------------------------------- 1461 1462 void DrawTestHistograms( 1463 int verbosityLevel 1464 ) 1465 { 1466 if (verbosityLevel > 2) cout << endl << "...drawing pulse histograms" ; 1467 1468 cgpTestHistos->cd(1); 1469 hTesthisto->Draw(); 1470 cgpTestHistos->cd(2); 1471 hTesthisto2->Draw("BOX"); 1472 } 1473 // end of DrawTestHistograms 1474 //---------------------------------------------------------------------------- 1475 1476 void DrawPulseHistograms( 447 448 void 449 UpdateCanvases( 1477 450 int verbosityLevel, 1478 int max_pulse_order 1479 ) 1480 { 1481 if (verbosityLevel > 2) cout << endl << "...drawing pulse histograms" ; 1482 for (int pulse_order = 0; pulse_order < max_pulse_order; pulse_order++) 1483 { 1484 cgpPixelPulses[pulse_order]->cd(1); 1485 hPixelOverlay[pulse_order]->Draw("COLZ"); 1486 cgpPixelPulses[pulse_order]->cd(3); 1487 hPixelProfile[pulse_order]->Draw("COLZ"); 1488 hPixelProfile2[pulse_order]->Draw("SAME"); 1489 1490 cgpPixelPulses[pulse_order]->cd(2); 1491 hPixelEdgeOverlay[pulse_order]->Draw("COLZ"); 1492 // cgpPixelPulses[pulse_order]->cd(4); 1493 // hTesthisto->Draw(); 1494 1495 cgpAllPixelPulses[pulse_order]->cd(1); 1496 hAllPixelOverlay[pulse_order]->Draw("COLZ"); 1497 cgpAllPixelPulses[pulse_order]->cd(2); 1498 hAllPixelEdgeOverlay[pulse_order]->Draw("COLZ"); 1499 1500 cgpAllPixelPulses[pulse_order]->cd(3); 1501 hAllPixelProfile[pulse_order]->Draw("COLZ"); 1502 hAllPixelProfile2[pulse_order]->Draw("SAME"); 1503 1504 } 1505 } 1506 // end of DrawPulseHistograms 1507 //---------------------------------------------------------------------------- 1508 1509 1510 void DrawMaximumHistograms( 1511 int verbosityLevel, 1512 int max_pulse_order 1513 ) 1514 { 1515 if (verbosityLevel > 2) cout << endl << "...drawing maximum histograms" ; 1516 for (int pulse_order = 0; pulse_order < max_pulse_order; pulse_order++) 1517 { 1518 cgpPixelPulses[pulse_order]->cd(4); 1519 hPixelMax[pulse_order]->Draw(); 1520 hPixelMedian[pulse_order]->Draw("SAME"); 1521 hPixelMean[pulse_order]->Draw("SAME"); 1522 1523 cgpAllPixelPulses[pulse_order]->cd(4); 1524 hAllPixelMax[pulse_order]->Draw(); 1525 hAllPixelMedian[pulse_order]->Draw("SAME"); 1526 hAllPixelMean[pulse_order]->Draw("SAME"); 1527 // cgpAllPixelPulses[pulse_order]->cd(3); 1528 // hAllPixelMaxGaus[pulse_order]->Draw("COLZ"); 1529 } 1530 } 1531 // end of DrawMaximumHistograms 1532 //---------------------------------------------------------------------------- 1533 1534 1535 void UpdateCanvases( 1536 int verbosityLevel, 1537 int max_pulse_order 451 int max_pulse_order, 452 bool testmode 1538 453 ) 1539 454 { … … 1543 458 cgpPixelPulses[pulse_order]->Modified(); 1544 459 cgpPixelPulses[pulse_order]->Update(); 1545 cgpAllPixelPulses[pulse_order]->Modified(); 1546 cgpAllPixelPulses[pulse_order]->Update(); 1547 cgpTestHistos->Modified(); 1548 cgpTestHistos->Update(); 460 cgpDistributions[pulse_order]->Modified(); 461 cgpDistributions[pulse_order]->Update(); 462 // if ( testmode ) 463 // { 464 // cgpTestHistos->Modified(); 465 // cgpTestHistos->Update(); 466 // } 1549 467 } 1550 468 } 1551 // end of UpdateCanvases 1552 //---------------------------------------------------------------------------- 1553 1554 1555 1556 void 1557 PlotMedianEachSliceOfPulse( 1558 unsigned int max_pulse_order, 1559 int fitdata, 1560 int verbosityLevel 1561 ) 469 //----------------------------------------------------------------------------- 470 // MAIN Funktion for C compilation 471 //----------------------------------------------------------------------------- 472 473 int main() 1562 474 { 1563 if (verbosityLevel > 2) cout << endl 1564 << "...calculating pulse shape of slice's Median" 1565 << endl; 1566 for (unsigned int pulse_order = 0 ; pulse_order < max_pulse_order ; pulse_order ++) 1567 { 1568 if (verbosityLevel > 2) cout << "\t...calculation of " 1569 << "pulse order " << pulse_order; 1570 1571 Int_t nbins = hPixelOverlay[pulse_order]->GetXaxis()->GetNbins(); 1572 1573 for (Int_t TimeSlice=1;TimeSlice<=nbins;TimeSlice++) { 1574 TH1 *h1 = hPixelOverlay[pulse_order]->ProjectionY("",TimeSlice,TimeSlice); 1575 float median = MedianOfH1(h1); 1576 float mean = h1->GetMean(); 1577 if (verbosityLevel > 2) printf("Median of Slice %d, Median=%g\n",TimeSlice,median); 1578 delete h1; 1579 hPixelMedian[pulse_order]->SetBinContent(TimeSlice, median ); 1580 hPixelMean[pulse_order]->SetBinContent(TimeSlice, mean ); 1581 hAllPixelMedian[pulse_order]->SetBinContent(TimeSlice, median ); 1582 hAllPixelMean[pulse_order]->SetBinContent(TimeSlice, mean ); 1583 } 1584 1585 if (verbosityLevel > 2) cout << "\t...done" << endl; 1586 1587 if (fitdata) 1588 { 1589 FitMaxPropabilityPuls( 1590 hPixelMedian[pulse_order], 1591 verbosityLevel 1592 ); 1593 } 1594 } 475 476 FCalcPulseTemplate(); 477 return 0; 478 1595 479 } 1596 // end of PlotMedianEachSliceOfPulse1597 //----------------------------------------------------------------------------1598 1599 Double_t MedianOfH1 (1600 TH1* h11601 )1602 {1603 //compute the median for 1-d histogram h11604 Int_t nbins = h1->GetXaxis()->GetNbins();1605 Double_t *x = new Double_t[nbins];1606 Double_t *y = new Double_t[nbins];1607 for (Int_t i=0;i<nbins;i++) {1608 x[i] = h1->GetXaxis()->GetBinCenter(i+1);1609 y[i] = h1->GetBinContent(i+1);1610 }1611 Double_t median = TMath::Median(nbins,x,y);1612 delete [] x;1613 delete [] y;1614 return median;1615 }1616 // end of PlotMedianEachSliceOfPulse1617 //----------------------------------------------------------------------------1618 1619 void1620 PlotPulseShapeOfMaxPropability(1621 unsigned int max_pulse_order,1622 int fitdata,1623 int verbosityLevel1624 )1625 {1626 if (verbosityLevel > 2) cout << endl1627 << "...calculating pulse shape of max. propability"1628 << endl;1629 for (unsigned int pulse_order = 0 ; pulse_order < max_pulse_order ; pulse_order ++)1630 {1631 if (verbosityLevel > 2) cout << "\t...calculation of "1632 << "pulse order " << pulse_order;1633 // vector max_value_of to number of timeslices in Overlay Spectra1634 float max_value_of_slice;1635 1636 Int_t nbins = hPixelOverlay[pulse_order]->GetXaxis()->GetNbins();1637 if (verbosityLevel > 2) cout << "...generating projections" << endl;1638 for (Int_t TimeSlice=1;TimeSlice<=nbins;TimeSlice++)1639 {1640 if (verbosityLevel > 2) cout << "...Timeslice: " << TimeSlice;1641 //put maximumvalue of every Y-projection of every timeslice into vector1642 TH1D *h1 = hPixelOverlay[pulse_order]->ProjectionY( "", TimeSlice,TimeSlice);1643 1644 max_value_of_slice = h1->GetBinCenter( h1->GetMaximumBin() );1645 1646 if (verbosityLevel > 2) cout << " with max value "1647 << max_value_of_slice << endl;1648 hPixelMax[pulse_order]->SetBinContent(TimeSlice, max_value_of_slice );1649 hAllPixelMax[pulse_order]->SetBinContent(TimeSlice, max_value_of_slice );1650 delete h1;1651 }1652 1653 if (verbosityLevel > 2) cout << "\t...done" << endl;1654 1655 if (fitdata)1656 {1657 FitMaxPropabilityPuls(1658 hPixelMax[pulse_order],1659 verbosityLevel1660 );1661 }1662 }1663 }1664 1665 1666 void1667 FitMaxPropabilityPuls(1668 TH1F* hMaximumTemp,1669 int verbosityLevel1670 )1671 {1672 if (verbosityLevel > 2) cout << "...fit Landau in histograms" ;1673 if (verbosityLevel > 2) cout << "\t...calculating Landaufit" ;1674 hMaximumTemp->Fit("landau", "", "", -50, 250);1675 if (verbosityLevel > 2) cout << "...done" << endl;1676 }1677 1678 1679 //void SetHistogrammSettings(1680 // const char* histoType,1681 // int max_puls_order,1682 // int MaxAmplOfFirstPulse,1683 // int verbosityLevel1684 // )1685 //{1686 // if (verbosityLevel > 2) cout << endl << "...set histogram Settings" << endl;1687 // for (int pulse_order = 1; pulse_order <= max_puls_order; pulse_order++)1688 // {1689 // qHistSetting[pulse_order].yMax = (MaxAmplOfFirstPulse * pulse_order);1690 // cout << "\t Max @ " << gHistSetting[pulse_order].yMax << endl;1691 // gHistSetting[pulse_order].name = "h";1692 // gHistSetting[pulse_order].name += histoType;1693 // gHistSetting[pulse_order].name += "PeakOverlay";1694 // gHistSetting[pulse_order].name += pulse_order;1695 // gHistSetting[pulse_order].title = "PeakOverlay with Amplitude around ";1696 // gHistSetting[pulse_order].title += gHistSetting[pulse_order].yMax;1697 // gHistSetting[pulse_order].title += " mV";1698 // gHistSetting[pulse_order].Mname = "h";1699 // gHistSetting[pulse_order].Mname += histoType;1700 // gHistSetting[pulse_order].Mname += "PeakMaximum";1701 // gHistSetting[pulse_order].Mname += pulse_order;1702 // gHistSetting[pulse_order].Mtitle = "Peak (approx) derived with Maximum of above Spektrum with Max of ";1703 // gHistSetting[pulse_order].Mtitle += gHistSetting[pulse_order].yMax;1704 // gHistSetting[pulse_order].Mtitle += " mV";1705 // }1706 // if (verbosityLevel > 2) cout << "...done" << endl;1707 //}1708 1709 1710 //void CreatePulseCanvas(1711 // TCanvas* array_of_canvases,1712 // unsigned int max_pulse_order,1713 // const char* canvas_name,1714 // const char* canvas_title,1715 // const char* pixel,1716 // int verbosityLevel1717 // )1718 //{1719 // if (verbosityLevel > 2) cout << "...building Canvases" ;1720 1721 // for (1722 // unsigned int pulse_order = 1;1723 // pulse_order <= max_pulse_order;1724 // pulse_order++1725 // )1726 // {1727 // TString name = canvas_name;1728 // name += " ";1729 // name += pulse_order;1730 // name += " ";1731 1732 // TString title = "Order ";1733 // title += pulse_order;1734 // title += " ";1735 // title += canvas_title;1736 // title += " of ";1737 // title += pixel;1738 // title += ". Pixel";1739 1740 // array_of_canvases[pulse_order] = new TCanvas(name, title, 1, 1, 1400, 700);1741 // array_of_canvases[pulse_order].Divide(2,2);1742 // }1743 // if (verbosityLevel > 2) cout << "...done" << endl;1744 //}1745 1746 1747 //void1748 //DeletePulseCanvas(1749 // TCanvas* array_of_canvases,1750 // unsigned int max_pulse_order1751 // )1752 //{1753 // for (1754 // unsigned int pulse_order = 1;1755 // pulse_order <= max_pulse_order;1756 // pulse_order++1757 // )1758 // {1759 // delete array_of_canvases[pulse_order];1760 // }1761 //}1762 1763 1764 1765 1766 void1767 DeletePixelCanvases( int verbosityLevel )1768 {1769 if (verbosityLevel > 2) cout << endl << "...delete pixel Canvas" ;1770 for (int pulse_order = 0; pulse_order < MAX_PULS_ORDER; pulse_order ++)1771 {1772 delete cgpPixelPulses[pulse_order];1773 }1774 }1775 1776 //----------------------------------------------------------------------------1777 //----------------------------------------------------------------------------1778 1779 //void SetHistogramAttributes(1780 // const char* histoType,1781 // int pulse_order,1782 // histAttrib_t* histoAttrib,1783 // int max_ampl_of_first_pulse1784 // )1785 //{1786 // histoAttrib[pulse_order].yMax = (max_ampl_of_first_pulse * pulse_order);1787 // histoAttrib[pulse_order].yMin = (histoAttrib[pulse_order].yMax/2);1788 // histoAttrib[pulse_order].name = "h";1789 // histoAttrib[pulse_order].name += histoType;1790 // histoAttrib[pulse_order].name = "_";1791 // histoAttrib[pulse_order].name += pulse_order;1792 1793 // histoAttrib[pulse_order].title += histoType;1794 // histoAttrib[pulse_order].title += " with Amplitude around ";1795 // histoAttrib[pulse_order].title += histoAttrib[pulse_order].yMax;1796 // histoAttrib[pulse_order].title += " mV";1797 1798 // histoAttrib[pulse_order].xTitle += "Timeslices [a.u.]";1799 // histoAttrib[pulse_order].yTitle += "Amplitude [mV]";1800 //}1801 1802 void WritePixelTemplateToCsv(1803 TString path,1804 const char* csv_file_name,1805 const char* overlay_method,1806 int pixel,1807 int verbosityLevel1808 )1809 {1810 path += csv_file_name;1811 path += "_";1812 path += pixel;1813 path += ".csv";1814 1815 Int_t nbins = hPixelMax[0]->GetXaxis()->GetNbins();1816 if (verbosityLevel > 0)1817 {1818 cout << "writing point-set to csv file: " ;1819 cout << path << endl;1820 cout << "...opening file" << endl;1821 }1822 if (verbosityLevel > 2) cout << "...number of bins " << nbins << endl;1823 ofstream out;1824 out.open( path );1825 out << "### point-set of a single photon pulse template" << endl;1826 out << "### template determined with pulse overlay at: "1827 << overlay_method << endl;1828 out << "### Slice's Amplitude determined by calculating the " << endl1829 << "### value of maximum propability of slice -> Amplitude1 " << endl1830 << "### mean of slice -> Amplitude2 " << endl1831 << "### median of slice -> Amplitude3 " << endl1832 << "### for each slice" << endl;1833 out << "### Pixel number (CHid): " << pixel << endl1834 << endl;1835 1836 out << "time [slices],Amplitude1 [mV],Amplitude2 [mV],Amplitude3 [mV]" << endl;1837 1838 for (int TimeSlice=1;TimeSlice<=nbins;TimeSlice++)1839 {1840 out << TimeSlice << "," ;1841 out << hPixelMax[0]->GetBinContent(TimeSlice) << ",";1842 out << hPixelMean[0]->GetBinContent(TimeSlice) << ",";1843 out << hPixelMedian[0]->GetBinContent(TimeSlice) << endl;1844 }1845 out.close();1846 if (verbosityLevel > 0) cout << "...file closed" << endl;1847 }1848 1849 void WriteAllPixelTemplateToCsv(1850 TString path,1851 const char* csv_file_name,1852 const char* overlay_method,1853 int verbosityLevel1854 )1855 {1856 path += csv_file_name;1857 path += "_AllPixel.csv";1858 1859 Int_t nbins = hAllPixelMax[0]->GetXaxis()->GetNbins();1860 if (verbosityLevel > 0)1861 {1862 cout << "writing point-set to csv file: " ;1863 cout << path << endl;1864 cout << "...opening file" << endl;1865 }1866 if (verbosityLevel > 2) cout << "...number of bins " << nbins << endl;1867 ofstream out;1868 out.open( path );1869 out << "### point-set of a single photon pulse template" << endl;1870 out << "### template determined with pulse overlay at: "1871 << overlay_method << "of all Pixels";1872 out << "### Slice's Amplitude determined by calculating the " << endl1873 << "### value of maximum propability of slice -> Amplitude1 " << endl1874 << "### mean of slice -> Amplitude2 " << endl1875 << "### median of slice -> Amplitude3 " << endl1876 << "### for each slice" << endl1877 << endl;1878 1879 out << "time [slices],Amplitude1 [mV],Amplitude2 [mV],Amplitude3 [mV]" << endl;1880 1881 for (int TimeSlice=1;TimeSlice<=nbins;TimeSlice++)1882 {1883 out << TimeSlice << "," ;1884 out << hAllPixelMax[0]->GetBinContent(TimeSlice) << ",";1885 out << hAllPixelMean[0]->GetBinContent(TimeSlice) << ",";1886 out << hAllPixelMedian[0]->GetBinContent(TimeSlice) << endl;1887 }1888 out.close();1889 if (verbosityLevel > 0) cout << "...file closed" << endl;1890 }
Note:
See TracChangeset
for help on using the changeset viewer.
