/********************** FPulseTemplate *********************** * read FACT raw data * remove spikes * calculate baseline * find peaks (CFD and normal discriminator) * compute pulse height and pulse integral spektrum of the peaks + search for Peaks in data + put peaks into different histos depending und Amplitude + draw histos for single, double, tripple, ....above 100mV Photonpulses + draw Parameterdevelopment depending on photon quantity + form a template (shape depending on number of photons) so it can be used for detection of other peaks *************************************************************/ //---------------------------------------------------------------------------- // root libraries //---------------------------------------------------------------------------- #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #define NPIX 1440 #define NCELL 1024 #define FAD_MAX_SAMPLES 1024 #define MAX_PULS_ORDER 1 #define HAVE_ZLIB //---------------------------------------------------------------------------- // rootmacros //---------------------------------------------------------------------------- #include "fits.h" #include "openFits.h" #include "openFits.c" #include "discriminator.h" #include "discriminator.C" #include "zerosearch.h" #include "zerosearch.C" #include "factfir.C" #include "DrsCalibration.C" #include "DrsCalibration.h" #include "SpikeRemoval.h" #include "SpikeRemoval.C" //---------------------------------------------------------------------------- // Decleration of global variables //---------------------------------------------------------------------------- bool breakout=false; int gNEvents; vector AllPixelDataVector; vector StartCellVector; unsigned int CurrentEventID; size_t PXLxROI; UInt_t gRegionOfInterest; UInt_t NumberOfPixels; TString histotitle; size_t TriggerOffsetROI, RC; size_t drs_n; vector drs_basemean; vector drs_gainmean; vector drs_triggeroffsetmean; vector Ameas(FAD_MAX_SAMPLES); // copy of the data (measured amplitude vector Vcorr(FAD_MAX_SAMPLES); // corrected Values vector Vslide(FAD_MAX_SAMPLES); // sliding average result vector Vcfd(FAD_MAX_SAMPLES); // CDF result vector Vcfd2(FAD_MAX_SAMPLES); // CDF result + 2nd sliding average float gGainMean = 9; float gBSLMean = 0; typedef struct{ float maxAmpl; float countOfMax; } OverlayMaximum; //typedef struct{ // TString name; // TString title; // TString xTitle; // TString yTitle; // float yMax; // float yMin; //} histAttrib_t; //histAttrib_t* gHistoAttrib[MAX_PULS_ORDER]; //histAttrib_t* gMaximumAttrib[MAX_PULS_ORDER]; //histAttrib_t* gMaxGausAttrib[MAX_PULS_ORDER]; //histAttrib_t* gProfileAttrib[MAX_PULS_ORDER]; // histograms const int Number_Of_Debug_Histo_Types = 7; const unsigned int Ameas_ = 0, N1mean_ = 1, Vcorr_ = 2, Vtest_ = 3, Vslide_ = 4, Vcfd_ = 5, Vcfd2_ = 6; //const char* gHistoTypes[8] = { // "hPixelOverlay", // "hPixelMax", // "hPixelEdgeOverlay", // "hPixelProfile", // "hAllPixelOverlay", // "hAllPixelMax", // "hAllPixelMaxGaus", // "hAllPixelProfile" //}; //---------------------------------------------------------------------------- // Initialisation of histograms //---------------------------------------------------------------------------- // Temporary Objects TH1F* debugHistos = NULL; TH2F* hOverlayTemp = NULL; TH1F* hMaximumTemp = NULL; TH1D* hProjPeak = NULL; TH1F* hTesthisto = NULL; //Pixelwise Histograms TH2F* hPixelOverlay[MAX_PULS_ORDER]= {NULL};//histogrammm for overlay of detected Peaks TProfile* hPixelProfile[MAX_PULS_ORDER] = {NULL};//histogrammm for Profile of detected Peaks TH1F* hPixelMax[MAX_PULS_ORDER] = {NULL}; TH2F* hPixelEdgeOverlay[MAX_PULS_ORDER] = {NULL}; //All Pixel Histograms TH2F* hAllPixelOverlay[MAX_PULS_ORDER] = {NULL}; TProfile* hAllPixelProfile[MAX_PULS_ORDER] = {NULL}; TH1F* hAllPixelMax[MAX_PULS_ORDER] = {NULL}; TH1F* hAllPixelMaxGaus[MAX_PULS_ORDER] = {NULL}; //Histogram Parameters Int_t gPixelOverlayXaxisLeft = 0; Int_t gPixelOverlayXaxisRight = 0; //Root-File Objects TObjArray* hList = NULL; TObjArray* hAllPixelList = NULL; TObjArray* hRootList = NULL; //Canvases TCanvas* cgpPixelPulses[MAX_PULS_ORDER] = {NULL}; TCanvas* cgpAllPixelPulses[MAX_PULS_ORDER] = {NULL}; TCanvas* cgpTestHistos = NULL; //TCanvas* gpcDevelopment = NULL; //---------------------------------------------------------------------------- // Functions //---------------------------------------------------------------------------- void BookHistos(int ); void BookPixelHistos(int ); void BookTestHistos( int ); void DeletePixelHistos(int ); void SaveHistograms( const char*, const char*, TObjArray*, int ); void FillHistograms(vector*, int, int); void PlotPulseShapeOfMaxPropability(unsigned int, int, int); void DrawPulseHistograms(int, int); void FitMaxPropabilityPuls( TH1F* , int ); void DrawMaximumHistograms( int, int ); void DrawTestHistograms( int); void UpdateCanvases( int, int); void DeletePixelCanvases( int ); void CreateRootFile(const char*, int ); TFile *OpenRootFile(const char*, int ); void CloseRootFile(TFile *tf); TString CreateSubDirName(int); //creates a String containing the path to the subdirectory in root file TString CreateSubDirName(const char* ); //void SetHistogrammSettings( const char*, int, int , int); //void CreatePulseCanvas( TCanvas* , unsigned int, const char* , const char*, const char*, int); //void DeletePulseCanvas( TCanvas* , unsigned int); //---------------------------------------------------------------------------- //---------------------------------------------------------------------------- // MAIN - Funtion //---------------------------------------------------------------------------- //---------------------------------------------------------------------------- int FPulseTemplate( char* datafilename = "../data/2011/11/09/20111109_006.fits.gz", const char* drsfilename = "../data/2011/11/09/20111109_003.drs.fits.gz", const char* OutRootFileName = "../analysis/FPulseTemplate/20111109_006/20111109_006.pulses.root", bool ProduceGraphic = true, int refresh_rate = 500, //refresh rate for canvases bool spikeDebug = false, bool debugPixel = true, bool fitdata = false, int verbosityLevel = 1, // different verbosity levels can be implemented here int firstevent = 272, int nevents = 500, int firstpixel = 1100, int npixel = 1, int AmplWindowWidth = 14, //Width of Window for selection of pulses to histograms float GainMean = 8, float BSLMean = -1, //4 Histogramms will be drawn, decide how far pulsehights differ from eachother int avg1 = 8, int avg2 = 8, int OverlayWindowLeft = 70, int OverlayWindowRight = 230 ) { gGainMean = GainMean; gBSLMean = BSLMean; //---------------------------------------------------------------------------- // Save-Root-File Settings //---------------------------------------------------------------------------- CreateRootFile( OutRootFileName, verbosityLevel ); //---------------------------------------------------------------------------- // root global Settings //---------------------------------------------------------------------------- gPixelOverlayXaxisLeft = OverlayWindowLeft; gPixelOverlayXaxisRight = OverlayWindowRight; gStyle->SetPalette(1,0); gROOT->SetStyle("Plain"); // gPad->SetGrid(); for ( int pulse_order = MAX_PULS_ORDER - 1; pulse_order >= 0 ; pulse_order-- ) { TString cName ="cgpPixelPulses"; cName += pulse_order; TString cTitle ="Pulses of Order: "; cTitle += pulse_order; cgpPixelPulses[pulse_order] = new TCanvas(cName,cTitle, 0,pulse_order*20,800,800); cgpPixelPulses[pulse_order]->Divide(2, 2); cName = "cgpAllPixelPulses"; cName += pulse_order; cTitle ="AllPixel average of puls shapes of Order: "; cTitle += pulse_order; cgpAllPixelPulses[pulse_order] = new TCanvas( cName, cTitle, 801, pulse_order*20, 800, 800 ); cgpAllPixelPulses[pulse_order]->Divide(2, 2); } // Create (pointer to) Canvases, which are used in every run, // also in 'non-debug' runs // Canvases only need if spike Debug, but I want to deklare // the pointers anyway ... TCanvas *cFiltered = NULL; if (spikeDebug) { cFiltered = new TCanvas("cFiltered","filtered DRS Waveforms", 1,310,400,300); cFiltered->Divide(1, 3); } cgpTestHistos = new TCanvas( "cgpTestHistos", "Test Histograms", 801, 0, 800, 800 ); cgpTestHistos->Divide(2,0); //----------------------------------------------------------------------------- // Filter-Settings //----------------------------------------------------------------------------- // CFD filter settings int k_cfd = 10; vector a_cfd(k_cfd, 0); double b_cfd = 1.; a_cfd[0]=-0.75; a_cfd[k_cfd-1]=1.; //----------------------------------------------------------------------------- // prepare datafile //----------------------------------------------------------------------------- // Open the data file fits * datafile; // Opens the raw data file and 'binds' the variables given as // Parameters to the data file. So they are filled with // raw data as soon as datafile->GetRow(int) is called. gNEvents = openDataFits( datafilename, &datafile, AllPixelDataVector, StartCellVector, CurrentEventID, gRegionOfInterest, NumberOfPixels, PXLxROI, verbosityLevel ); if (gNEvents == 0) { cout << "return code of OpenDataFile:" << datafilename<< endl; cout << "is zero -> aborting." << endl; return 1; } if (verbosityLevel > 0) { cout << endl <<"number of events in file: "<< gNEvents << "\t"; } if ( nevents == -1 || nevents > gNEvents ) nevents = gNEvents; // -1 means all! if (verbosityLevel > 0) { cout <<"of, which "<< nevents << " will be processed"<< endl; cout <<"Total # of Pixel: "<< NumberOfPixels << "\t"; } if ( npixel == -1 || npixel > (int)NumberOfPixels ) npixel = (int)NumberOfPixels; // -1 means all! if (verbosityLevel > 0) { cout <<"of, which "<< npixel << " will be processed"<< endl; } //Get the DRS calibration RC = openCalibFits( drsfilename, drs_basemean, drs_gainmean, drs_triggeroffsetmean, TriggerOffsetROI); if (RC == 0) { cout << "return code of openCalibFits:" << drsfilename << endl; cout << "is zero -> aborting." << endl; return 1; } // Book the histograms BookHistos(verbosityLevel ); BookTestHistos( verbosityLevel ); //------------------------------------- // Loop over every Pixel //------------------------------------- for ( int pixel = firstpixel; pixel < firstpixel + npixel; pixel++ ) { if (verbosityLevel > 0) { cout << "------------------------------------------------" << endl << "...processing Pixel: " << pixel << endl; } BookPixelHistos(verbosityLevel ); //-------------------------------------------------------------------- // Loops over Every Event of Pixel //-------------------------------------------------------------------- for ( int ev = firstevent; ev < firstevent + nevents; ev++) { // Get an Event --> consists of 1440 Pixel ...erm....data datafile->GetRow( ev ); if (verbosityLevel > 0) { cout << "-------------------------------------" << endl << "...processing Event: " << CurrentEventID << "/" << nevents << " of Pixel: " << pixel << endl; } //------------------------------------- // Apply Calibration //------------------------------------- if (verbosityLevel > 2) cout << "...applying DrsCalibration"; applyDrsCalibration( Ameas, pixel, 12, 12, drs_basemean, drs_gainmean, drs_triggeroffsetmean, gRegionOfInterest, AllPixelDataVector, StartCellVector ); if (verbosityLevel > 2) cout << "...done " << endl; //------------------------------------- // Apply Filters //------------------------------------- // finds spikes in the raw data, and interpolates the value // spikes are: 1 or 2 slice wide, positive non physical artifacts if (verbosityLevel > 2) cout << "...removeing Spikes"; removeSpikes (Ameas, Vcorr); if (verbosityLevel > 2) cout << "...done " << endl; // filter Vcorr with sliding average using FIR filter function if (verbosityLevel > 2) cout << "...applying sliding average filter"; sliding_avg(Vcorr, Vslide, avg1); if (verbosityLevel > 2) cout << "...done " << endl; // filter Vslide with CFD using FIR filter function if (verbosityLevel > 2) cout << "...apllying factfir filter"; factfir(b_cfd , a_cfd, k_cfd, Vslide, Vcfd); if (verbosityLevel > 2) cout << "...done " << endl; // filter Vcfd with sliding average using FIR filter function if (verbosityLevel > 2) cout << "...applying 2nd sliding average filter"; sliding_avg(Vcfd, Vcfd2, avg2); if (verbosityLevel > 2) cout << "...done " << endl; //------------------------------------- // Search vor Zero crossings //------------------------------------- if (verbosityLevel > 2) cout << endl << "...searching zero crossings" ; // peaks in Ameas[] are found by searching for zero crossings // in Vcfd2 // first Argument 1 means ... *rising* edge // second Argument 1 means ... search with stepsize 1 ... 10 is okay as well vector* pZXings = zerosearch( Vcfd2 , 1 , 8); // pZXings means "zero cross ings" EnlargeRegion(*pZXings, 10, 10); findAbsMaxInRegions(*pZXings, Vslide); removeMaximaBelow( *pZXings, 3.0); removeRegionWithMaxOnEdge( *pZXings, 2); removeRegionOnFallingEdge( *pZXings, 100); findTimeOfHalfMaxLeft(*pZXings, Vslide, gBSLMean, 5, 10, verbosityLevel ); if (verbosityLevel > 2) cout << "...done" << endl; //----------------------------------------------------------------------------- // Fill Overlay Histos //----------------------------------------------------------------------------- FillHistograms(pZXings, AmplWindowWidth, verbosityLevel ); //----------------------------------------------------------------------------- // Spike Debug //----------------------------------------------------------------------------- if ( spikeDebug ) { // TODO do this correct. The vectors should be the rigt ones... this is just luck debugHistos[Ameas_].GetXaxis()->Set(Ameas.size() , -0.5 , Ameas.size()-0.5); debugHistos[Vcorr_].GetXaxis()->Set(Ameas.size() , -0.5 , Ameas.size()-0.5); debugHistos[Vslide_].GetXaxis()->Set(Ameas.size() , -0.5 , Ameas.size()-0.5); debugHistos[Vcfd_].GetXaxis()->Set(Ameas.size() , -0.5 , Ameas.size()-0.5); debugHistos[Vcfd2_].GetXaxis()->Set(Ameas.size() , -0.5 , Ameas.size()-0.5); for ( unsigned int sl = 0; sl < gRegionOfInterest; sl++) { debugHistos[Ameas_].SetBinContent(sl, Ameas[sl]); debugHistos[Vcorr_].SetBinContent(sl, Vcorr[sl]); debugHistos[Vslide_].SetBinContent( sl, Vslide[sl] ); debugHistos[Vcfd_].SetBinContent( sl, Vcfd[sl] ); debugHistos[Vcfd2_].SetBinContent( sl, Vcfd2[sl] ); } cFiltered->cd(1); gPad->SetGrid(); debugHistos[Ameas_].Draw(); cFiltered->cd(2); gPad->SetGrid(); debugHistos[Vcorr_].Draw(); cFiltered->cd(3); gPad->SetGrid(); debugHistos[Vslide_].Draw(); TBox *OneBox; vector MyBoxes; for (unsigned int i=0; isize(); i++){ OneBox = new TBox( pZXings->at(i).maxPos -10 , pZXings->at(i).maxVal -0.5, pZXings->at(i).maxPos +10 , pZXings->at(i).maxVal +0.5); OneBox->SetLineColor(kBlue); OneBox->SetLineWidth(1); OneBox->SetFillStyle(0); OneBox->SetFillColor(kRed); MyBoxes.push_back(OneBox); OneBox->Draw(); } // cFiltered->cd(3); // gPad->SetGrid(); // debugHistos[Vcfd2_].Draw(); // TLine *zeroline = new TLine(0, 0, 1024, 0); // zeroline->SetLineColor(kBlue); // zeroline->Draw(); cFiltered->Update(); //Process gui events asynchronously during input TTimer timer("gSystem->ProcessEvents();", 50, kFALSE); timer.TurnOn(); TString input = Getline("Type 'q' to exit, to go on: "); timer.TurnOff(); if (input=="q\n") { breakout=true; } //TODO!!!!!!!!! // do some Garbage collection ... // all the Objects on the heap should be deleted here. }// end of if(spikeDebug) delete pZXings; if (breakout) break; //------------------------------------- // Draw 1. Set of Pixel Histograms //------------------------------------- if ((ev % refresh_rate) == 0) { DrawPulseHistograms( verbosityLevel, MAX_PULS_ORDER ); DrawTestHistograms( verbosityLevel ); if (ProduceGraphic) { UpdateCanvases( verbosityLevel, MAX_PULS_ORDER ); } } if (breakout) break; if (verbosityLevel > 2) { cout << endl << "-------------------------------------"<< endl; } }//End of Loop over Events //------------------------------------- // end of Loops over Events //------------------------------------- //------------------------------------- // Histogramms of Maximas in Overlay Spectra //------------------------------------- PlotPulseShapeOfMaxPropability( MAX_PULS_ORDER, fitdata, verbosityLevel ); if (verbosityLevel > 2) cout << "...done" << endl; //here is what happends at the end of each loop over all pixels //Save Histograms of Pixel into Output rootfile DrawPulseHistograms( verbosityLevel, MAX_PULS_ORDER ); DrawMaximumHistograms( verbosityLevel, MAX_PULS_ORDER ); if (ProduceGraphic) { UpdateCanvases( verbosityLevel, MAX_PULS_ORDER ); } // SaveHistograms( // OutRootFileName, // CreateSubDirName(pixel), // hList, // verbosityLevel // ); if (debugPixel) { //Process gui events asynchronously during input TTimer timer("gSystem->ProcessEvents();", 50, kFALSE); timer.TurnOn(); TString input = Getline("Type 'q' to exit, to go on: "); timer.TurnOff(); if (input=="q\n") { break; } } DeletePixelHistos(verbosityLevel); if (verbosityLevel > 2) { cout << endl << "...End of Pixel" << endl << "------------------------------------------------" << endl; } } // End of Loop over all Pixels //------------------------------------- // Draw Histograms //------------------------------------- SaveHistograms( //save histograms of all pixel into output root file OutRootFileName, CreateSubDirName("All"), hAllPixelList, verbosityLevel ); // SaveHistograms( //save histograms of generell results into output root file // OutRootFileName, // "root", // hRootList, // verbosityLevel // ); if (ProduceGraphic) { UpdateCanvases( verbosityLevel, MAX_PULS_ORDER ); } DeletePixelCanvases( verbosityLevel ); return( 0 ); } //---------------------------------------------------------------------------- // end of main function //----------------------------------------------------------------------------- //----------------------------------------------------------------------------- // Funktions //----------------------------------------------------------------------------- void BookPixelHistos( int verbosityLevel) { if (verbosityLevel > 2) cout << endl << "...book pixel histograms" << endl; hList = new TObjArray; TString histoname; for (int order = 0; order < MAX_PULS_ORDER; order ++) { // SetHistogramAttributes( // "PeakOverlay", // order, // gOverlayAttrib, // MaxAmplOfFirstPulse // ); histoname = "hPixelOverlay"; histoname += order; if (verbosityLevel > 3) cout << "...booking " << histoname << endl; hPixelOverlay[order] = new TH2F( histoname, "Overlay of detected pulses of one pulse order for one Pixel", gPixelOverlayXaxisLeft + gPixelOverlayXaxisRight , (-1*gPixelOverlayXaxisLeft)-0.5, gPixelOverlayXaxisRight-0.5 , 512, -55.5, 300.5 ); hPixelOverlay[order]->SetAxisRange( gBSLMean - 5, (gGainMean*(order+1)) + 10, "Y"); hPixelOverlay[order]->GetXaxis()->SetTitle( "Timeslices [a.u.]" ); hPixelOverlay[order]->GetYaxis()->SetTitle( "Amplitude [mV]" ); //hPixelProfile->SetBit(TH2F::kCanRebin); hList->Add( hPixelOverlay[order] ); //------------------------------------------------------------------------ // SetHistogramAttributes( // "PeakMaximum", // order, // gMaximumAttrib, // MaxAmplOfFirstPulse // ); histoname = "hPixelMax"; histoname += order; if (verbosityLevel > 3) cout << "...booking " << histoname << endl; hPixelMax[order] = new TH1F( histoname, "Maximum value of each slice in overlay plot", gPixelOverlayXaxisLeft + gPixelOverlayXaxisRight , (-1*gPixelOverlayXaxisLeft)-0.5, gPixelOverlayXaxisRight-0.5 ); hPixelMax[order]->SetAxisRange( gBSLMean - 5, (gGainMean*(order+1)) + 10, "Y"); hPixelMax[order]->GetXaxis()->SetTitle( "Timeslices [a.u.]" ); hPixelMax[order]->GetYaxis()->SetTitle( "Amplitude [mV]" ); hList->Add( hPixelMax[order] ); //------------------------------------------------------------------------ // SetHistogramAttributes( // "PeakMaxGaus", // order, // gMaxGausAttrib, // MaxAmplOfFirstPulse // ); histoname = "hPixelEdgeOverlay"; histoname += order; if (verbosityLevel > 3) cout << "...booking " << histoname << endl; hPixelEdgeOverlay[order] = new TH2F( histoname, "Overlay at rising edge of detected pulses of one pulse order for one Pixel ", gPixelOverlayXaxisLeft + gPixelOverlayXaxisRight , (-1*gPixelOverlayXaxisLeft)-0.5, gPixelOverlayXaxisRight-0.5 , 512, -55.5, 300.5 ); hPixelEdgeOverlay[order]->SetAxisRange( gBSLMean - 5, (gGainMean*(order+1)) + 10, "Y"); hPixelEdgeOverlay[order]->GetXaxis()->SetTitle( "Timeslices [a.u.]" ); hPixelEdgeOverlay[order]->GetYaxis()->SetTitle( "Amplitude [mV]" ); hList->Add( hPixelEdgeOverlay[order] ); //------------------------------------------------------------------------ // SetHistogramAttributes( // "PeakProfile", // order, // gProfileAttrib, // MaxAmplOfFirstPulse // ); histoname = "hPixelProfile"; histoname += order; if (verbosityLevel > 3) cout << "...booking " << histoname << endl; hPixelProfile[order] = new TProfile( histoname, "Mean value of each slice in overlay plot (Tprofile)", gPixelOverlayXaxisLeft + gPixelOverlayXaxisRight ,//nbinsx (-1*gPixelOverlayXaxisLeft)-0.5, //xlow gPixelOverlayXaxisRight-0.5 , //xup // 512, //nbinsy -55.5, //ylow 300.5, //yup "s"); //option hPixelMax[order]->SetAxisRange( gBSLMean - 5, (gGainMean*(order+1)) + 10, "Y"); hPixelProfile[order]->GetXaxis()->SetTitle( "Timeslices [a.u.]" ); hPixelProfile[order]->GetYaxis()->SetTitle( "Amplitude [mV]" ); //hPixelProfile->SetBit(TH2F::kCanRebin); hList->Add( hPixelProfile[order] ); } if (verbosityLevel > 2) cout << "...done" << endl; } //end of BookPixelHistos //---------------------------------------------------------------------------- void DeletePixelHistos( int verbosityLevel ) { if (verbosityLevel > 2) cout << endl << "...delete current pixel histograms" ; for (int order = 0; order < MAX_PULS_ORDER; order ++) { if (verbosityLevel > 3) cout << endl << "...deleting hPixelOverlay" << order; delete hPixelOverlay[order]; if (verbosityLevel > 3) cout << endl << "...deleting hPixelMax" << order; delete hPixelMax[order]; if (verbosityLevel > 3) cout << endl << "...deleting hPixelEdgeOverlay" << order; delete hPixelEdgeOverlay[order]; if (verbosityLevel > 3) cout << endl << "...deleting hPixelProfile" << order; delete hPixelProfile[order]; } if (verbosityLevel > 3) cout << endl << "...deleting hList"; delete hList; if (verbosityLevel > 2) cout << endl << "...done" << endl; } //end of DeletePixelHistos //---------------------------------------------------------------------------- void BookTestHistos( int verbosityLevel ) { if (verbosityLevel > 2) cout << endl << "...book pixel histograms" << endl; hTesthisto = new TH1F ( "hTesthisto", "Deviation of rising edge and maximum", 200, -300, 300 ); hTesthisto->GetXaxis()->SetTitle( "Timeslices [a.u.]" ); hTesthisto->GetYaxis()->SetTitle( "counts" ); if (verbosityLevel > 2) cout << "...done" << endl; } //end of BookTestHistos //---------------------------------------------------------------------------- void BookHistos( int verbosityLevel ) { if (verbosityLevel > 2) cout << endl << "...book histograms" << endl; hAllPixelList = new TObjArray; hRootList = new TObjArray; debugHistos = new TH1F[ Number_Of_Debug_Histo_Types ]; for ( int type = 0; type < Number_Of_Debug_Histo_Types; type++){ debugHistos[ type ].SetBins(1024, 0, 1024); debugHistos[ type ].SetLineColor(1); debugHistos[ type ].SetLineWidth(2); debugHistos[ type ].SetStats(false); // set X axis paras debugHistos[ type ].GetXaxis()->SetLabelSize(0.08); debugHistos[ type ].GetXaxis()->SetTitleSize(0.08); debugHistos[ type ].GetXaxis()->SetTitleOffset(1.0); debugHistos[ type ].GetXaxis()->SetTitle(Form("Time slice [%.1f ns/slice]", 1./2.)); // set Y axis paras debugHistos[ type ].GetYaxis()->SetLabelSize(0.08); debugHistos[ type ].GetYaxis()->SetTitleSize(0.08); debugHistos[ type ].GetYaxis()->SetTitleOffset(0.3); debugHistos[ type ].GetYaxis()->SetTitle("Amplitude [mV]"); } // All Pixel Histograms //------------------------------------------------------------------------ TString histoname; for (int order = 0; order < MAX_PULS_ORDER; order ++) { // SetHistogramAttributes( // "AllPixelPeakOverlay", // order, // gOverlayAttrib, // MaxAmplOfFirstPulse // ); histoname = "hAllPixelOverlay"; histoname += order; if (verbosityLevel > 3) cout << "...booking " << histoname << endl; hAllPixelOverlay[order] = new TH2F( histoname, "Overlay of detected Pulses of one Order for All Pixels", gPixelOverlayXaxisLeft + gPixelOverlayXaxisRight , (-1*gPixelOverlayXaxisLeft)-0.5, gPixelOverlayXaxisRight-0.5 , 512, -55.5, 300.5 ); hAllPixelOverlay[order]->SetAxisRange( gBSLMean - 5, (gGainMean*(order+1)) + 10, "Y"); hAllPixelOverlay[order]->GetXaxis()->SetTitle( "Timeslices [a.u.]" ); hAllPixelOverlay[order]->GetYaxis()->SetTitle( "Amplitude [mV]" ); //hAllPixelOverlay->SetBit(TH2F::kCanRebin); hAllPixelList->Add( hAllPixelOverlay[order] ); //------------------------------------------------------------------------ // SetHistogramAttributes( // "AllPixelPeakMax", // order, // gMaximumAttrib, // MaxAmplOfFirstPulse // ); histoname = "hAllPixelMax"; histoname += order; if (verbosityLevel > 3) cout << "...booking " << histoname << endl; hAllPixelMax[order] = new TH1F( histoname, "Maximum value of each slice in overlay plot for All Pixels", gPixelOverlayXaxisLeft + gPixelOverlayXaxisRight , (-1*gPixelOverlayXaxisLeft)-0.5, gPixelOverlayXaxisRight-0.5 ); hAllPixelMax[order]->SetAxisRange( gBSLMean - 5, (gGainMean*(order+1)) + 10, "Y"); hAllPixelMax[order]->GetXaxis()->SetTitle( "Timeslices [a.u.]" ); hAllPixelMax[order]->GetYaxis()->SetTitle( "Amplitude [mV]" ); hAllPixelList->Add( hAllPixelMax[order] ); //------------------------------------------------------------------------ // SetHistogramAttributes( // "AllPixelPeakMaxGaus", // order, // gMaxGausAttrib, // MaxAmplOfFirstPulse // ); histoname = "hAllPixelMaxGaus"; histoname += order; if (verbosityLevel > 3) cout << "...booking " << histoname << endl; hAllPixelMaxGaus[order] = new TH1F( histoname, "Mean value of each slice in overlay plot for All Pixels", gPixelOverlayXaxisLeft + gPixelOverlayXaxisRight , (-1*gPixelOverlayXaxisLeft)-0.5, gPixelOverlayXaxisRight-0.5 ); hAllPixelMaxGaus[order]->SetAxisRange( gBSLMean - 5, (gGainMean*(order+1)) + 10, "Y"); hAllPixelMaxGaus[order]->GetXaxis()->SetTitle( "Timeslices [a.u.]" ); hAllPixelMaxGaus[order]->GetYaxis()->SetTitle( "Amplitude [mV]" ); hAllPixelList->Add( hAllPixelMaxGaus[order] ); //------------------------------------------------------------------------ // SetHistogramAttributes( // "AllPixelPeakProfile", // order, // gProfileAttrib, // MaxAmplOfFirstPulse // ); histoname = "hAllPixelProfile"; histoname += order; if (verbosityLevel > 3) cout << "...booking " << histoname << endl; hAllPixelProfile[order] = new TProfile( histoname, "Mean value of each slice in overlay plot for All Pixels (Tprofile)", gPixelOverlayXaxisLeft + gPixelOverlayXaxisRight ,//nbinsx (-1*gPixelOverlayXaxisLeft)-0.5, //xlow gPixelOverlayXaxisRight-0.5 , //xup // 512, //nbinsy -55.5, //ylow 300.5, //yup "s"); //option hAllPixelProfile[order]->SetAxisRange( gBSLMean - 5, (gGainMean*(order+1)) + 10, "Y"); hAllPixelProfile[order]->GetXaxis()->SetTitle( "Timeslices [a.u.]" ); hAllPixelProfile[order]->GetYaxis()->SetTitle( "Amplitude [mV]" ); //hPixelProfile->SetBit(TH2F::kCanRebin); hAllPixelList->Add( hAllPixelProfile[order] ); }//end of for over order if (verbosityLevel > 2) cout << "...done" << endl; } // end of BookHistos //---------------------------------------------------------------------------- void CreateRootFile(const char * loc_fname, int verbosityLevel) { TFile* tf = new TFile(loc_fname,"UPDATE"); if (tf->IsZombie()) { cout << "Error opening file" << endl; exit(-1); } else { if (verbosityLevel > 1) cout << "creating root-file successfull" << endl; tf->Close(); } } //end of CreateRootFile //---------------------------------------------------------------------------- TFile* OpenRootFile(const char * loc_fname, int verbosityLevel) { TFile* tf = new TFile(loc_fname,"UPDATE"); if (tf->IsZombie()) { cout << "Error opening file" << endl; exit(-1); } else { if (verbosityLevel > 1) cout << "...opening root-file:" << loc_fname << " successfull" << endl; } return tf; } //end of OpenRootFile //---------------------------------------------------------------------------- void CloseRootFile(TFile* tf) { if (tf->IsZombie()) { cout << "Error closing file" << endl; exit(-1); } else { tf->Close(); } } //end of CloseRootFile //---------------------------------------------------------------------------- void SaveHistograms( const char* loc_fname, const char* subdirectory, TObjArray* histList, int verbosityLevel ) { if (verbosityLevel > 2) cout << endl << "...saving pixel histograms to subdirectory: " << subdirectory << endl ; TFile* tf = OpenRootFile(loc_fname, verbosityLevel); if ( subdirectory != "root") { tf->cd(); tf->mkdir(subdirectory,subdirectory); tf->cd(subdirectory); } else { tf->cd(); } histList->Write(); // write the major histograms into the top level directory if (verbosityLevel > 3) tf->ls(); CloseRootFile( tf ); // close the file if (verbosityLevel > 2) cout << "...done" << endl; } // end of SaveHistograms //---------------------------------------------------------------------------- TString CreateSubDirName(int pixel) { TString path = "Pixel_"; path += pixel; // cout << "path " << path << endl ; return path; } // end of CreateSubDirName //---------------------------------------------------------------------------- TString CreateSubDirName(const char* title) { TString path = title; path += "_Pixel"; // cout << "path " << path << endl ; return path; } // end of CreateSubDirName //---------------------------------------------------------------------------- void FillHistograms( vector* pZXings, int AmplWindowWidth, int verbosityLevel ) { if (verbosityLevel > 2) cout << endl << "...filling pulse histograms" ; bool use_this_peak=false; int order_of_pulse=0; vector::iterator reg; //Loop over all found zerocrossings in Region for (reg = pZXings->begin() ; reg < pZXings->end() ; reg++) { //skip those who are at the Rim of the ROI if (reg->maxPos < 12 || (unsigned int) reg->maxPos > gRegionOfInterest-12) { if (verbosityLevel > 2) cout << endl << "\t...out of range" << endl; continue; } //domstest->Fill(reg->maxPos); // Define axis range of Histogram int Left = reg->maxPos - gPixelOverlayXaxisLeft; int Right = reg->maxPos + gPixelOverlayXaxisRight; int EdgeLeft = reg->halfRisingEdgePos - gPixelOverlayXaxisLeft; int EdgeRight = reg->halfRisingEdgePos + gPixelOverlayXaxisRight; if (Left < 0) Left =0; if (Right > (int)Vcorr.size() ) Right=Vcorr.size() ; hTesthisto->Fill( reg->maxPos - reg->halfRisingEdgePos ) ; if (verbosityLevel > 2) cout << endl << "\t...choosing Histogram" ; //determine order of pulse and which histogram shall be filled if (verbosityLevel > 2) cout << endl << "\t...choosing Histogram" ; for(int order = 0; order < MAX_PULS_ORDER; order++) { if (Ameas[reg->maxPos] >= ((gGainMean*(order+1)) - (AmplWindowWidth/2)) && Ameas[reg->maxPos] < ((gGainMean*(order+1)) + AmplWindowWidth/2)) { //hOverlayTemp = hPixelOverlay[order]; if (verbosityLevel > 2) cout << "...#" << order ; order_of_pulse = order; use_this_peak = true; break; } else if (order >= (MAX_PULS_ORDER - 1)){ use_this_peak = false; if (verbosityLevel > 2) cout << "...NONE" << endl ; } } //Fill overlay und profile histograms if (use_this_peak){ if (verbosityLevel > 2) cout << "...filling" ; for ( int pos = Left; pos < Right; pos++){ hPixelOverlay[order_of_pulse]->Fill( pos - (reg->maxPos), Ameas[pos]) ; hPixelProfile[order_of_pulse]->Fill( pos - (reg->maxPos), Ameas[pos]) ; hAllPixelOverlay[order_of_pulse]->Fill( pos - (reg->maxPos), Ameas[pos]) ; hAllPixelProfile[order_of_pulse]->Fill( pos - (reg->maxPos), Ameas[pos]) ; } for ( int pos = EdgeLeft; pos < EdgeRight; pos++){ // cout << endl << "###filling edge histo###" << endl; hPixelEdgeOverlay[order_of_pulse]->Fill( pos - (reg->halfRisingEdgePos), Ameas[pos]) ; // cout << endl << "######" << endl; } if (verbosityLevel > 2) cout << "...done" << endl; } //Breakout option if (breakout) break; } if (verbosityLevel > 2) cout << "...done" << endl; } // end of FillHistograms //---------------------------------------------------------------------------- void DrawTestHistograms( int verbosityLevel ) { if (verbosityLevel > 2) cout << endl << "...drawing pulse histograms" ; cgpTestHistos->cd(1); hTesthisto->Draw(); } // end of DrawTestHistograms //---------------------------------------------------------------------------- void DrawPulseHistograms( int verbosityLevel, int max_pulse_order ) { if (verbosityLevel > 2) cout << endl << "...drawing pulse histograms" ; for (int pulse_order = 0; pulse_order < max_pulse_order; pulse_order++) { cgpPixelPulses[pulse_order]->cd(1); hPixelOverlay[pulse_order]->Draw("COLZ"); cgpPixelPulses[pulse_order]->cd(2); hPixelProfile[pulse_order]->Draw("COLZ"); cgpPixelPulses[pulse_order]->cd(3); hPixelEdgeOverlay[pulse_order]->Draw("COLZ"); // cgpPixelPulses[pulse_order]->cd(4); // hTesthisto->Draw(); cgpAllPixelPulses[pulse_order]->cd(1); hAllPixelOverlay[pulse_order]->Draw("COLZ"); cgpAllPixelPulses[pulse_order]->cd(2); hAllPixelProfile[pulse_order]->Draw("COLZ"); } } // end of DrawPulseHistograms //---------------------------------------------------------------------------- void DrawMaximumHistograms( int verbosityLevel, int max_pulse_order ) { if (verbosityLevel > 2) cout << endl << "...drawing maximum histograms" ; for (int pulse_order = 0; pulse_order < max_pulse_order; pulse_order++) { cgpPixelPulses[pulse_order]->cd(4); hPixelMax[pulse_order]->Draw("LF2"); cgpAllPixelPulses[pulse_order]->cd(4); hAllPixelMax[pulse_order]->Draw(); // cgpAllPixelPulses[pulse_order]->cd(3); // hAllPixelMaxGaus[pulse_order]->Draw("COLZ"); } } // end of DrawMaximumHistograms //---------------------------------------------------------------------------- void UpdateCanvases( int verbosityLevel, int max_pulse_order ) { if (verbosityLevel > 3) cout << endl << "...updating canvases" ; for (int pulse_order = 0; pulse_order < max_pulse_order; pulse_order++) { cgpPixelPulses[pulse_order]->Modified(); cgpPixelPulses[pulse_order]->Update(); cgpAllPixelPulses[pulse_order]->Modified(); cgpAllPixelPulses[pulse_order]->Update(); cgpTestHistos->Modified(); cgpTestHistos->Update(); } } // end of UpdateCanvases //---------------------------------------------------------------------------- void PlotPulseShapeOfMaxPropability( unsigned int max_pulse_order, int fitdata, int verbosityLevel ) { if (verbosityLevel > 2) cout << endl << "...calculating pulse shape of max. propability" << endl; for (unsigned int pulse_order = 0 ; pulse_order < max_pulse_order ; pulse_order ++) { if (verbosityLevel > 2) cout << "\t...calculating of " << "pulse order " << pulse_order; // vector max_value_of to number of timeslices in Overlay Spectra vector max_value_of; max_value_of.resize(hPixelOverlay[pulse_order]->GetNbinsX()); for (int TimeSlice = 0; TimeSlice <= hPixelOverlay[pulse_order]->GetNbinsX(); TimeSlice++ ) { histotitle = "hproj_py"; //generate title of histogram of which MaxVal histotitle += pulse_order ; //will be calculated //put maximumvalue of every Y-projection of every timeslice into vector hProjPeak = hPixelOverlay[pulse_order]->ProjectionY(histotitle, TimeSlice, TimeSlice, ""); max_value_of[ TimeSlice ].maxAmpl = (hProjPeak->GetMaximumBin() * 0.5) - 3.5; max_value_of[ TimeSlice ].countOfMax = hProjPeak->GetBinContent( hProjPeak->GetMaximumBin() ); hPixelMax[pulse_order]->SetBinContent(TimeSlice, max_value_of[ TimeSlice ].maxAmpl ); } if (verbosityLevel > 2) cout << "\t...done" << endl; if (fitdata) { FitMaxPropabilityPuls( hPixelMax[pulse_order], verbosityLevel ); } } } void FitMaxPropabilityPuls( TH1F* hMaximumTemp, int verbosityLevel ) { if (verbosityLevel > 2) cout << "...fit Landau in histograms" ; if (verbosityLevel > 2) cout << "\t...calculating Landaufit" ; hMaximumTemp->Fit("landau", "", "", -50, 250); if (verbosityLevel > 2) cout << "...done" << endl; } //void SetHistogrammSettings( // const char* histoType, // int max_puls_order, // int MaxAmplOfFirstPulse, // int verbosityLevel // ) //{ // if (verbosityLevel > 2) cout << endl << "...set histogram Settings" << endl; // for (int pulse_order = 1; pulse_order <= max_puls_order; pulse_order++) // { // qHistSetting[pulse_order].yMax = (MaxAmplOfFirstPulse * pulse_order); // cout << "\t Max @ " << gHistSetting[pulse_order].yMax << endl; // gHistSetting[pulse_order].name = "h"; // gHistSetting[pulse_order].name += histoType; // gHistSetting[pulse_order].name += "PeakOverlay"; // gHistSetting[pulse_order].name += pulse_order; // gHistSetting[pulse_order].title = "PeakOverlay with Amplitude around "; // gHistSetting[pulse_order].title += gHistSetting[pulse_order].yMax; // gHistSetting[pulse_order].title += " mV"; // gHistSetting[pulse_order].Mname = "h"; // gHistSetting[pulse_order].Mname += histoType; // gHistSetting[pulse_order].Mname += "PeakMaximum"; // gHistSetting[pulse_order].Mname += pulse_order; // gHistSetting[pulse_order].Mtitle = "Peak (approx) derived with Maximum of above Spektrum with Max of "; // gHistSetting[pulse_order].Mtitle += gHistSetting[pulse_order].yMax; // gHistSetting[pulse_order].Mtitle += " mV"; // } // if (verbosityLevel > 2) cout << "...done" << endl; //} //void CreatePulseCanvas( // TCanvas* array_of_canvases, // unsigned int max_pulse_order, // const char* canvas_name, // const char* canvas_title, // const char* pixel, // int verbosityLevel // ) //{ // if (verbosityLevel > 2) cout << "...building Canvases" ; // for ( // unsigned int pulse_order = 1; // pulse_order <= max_pulse_order; // pulse_order++ // ) // { // TString name = canvas_name; // name += " "; // name += pulse_order; // name += " "; // TString title = "Order "; // title += pulse_order; // title += " "; // title += canvas_title; // title += " of "; // title += pixel; // title += ". Pixel"; // array_of_canvases[pulse_order] = new TCanvas(name, title, 1, 1, 1400, 700); // array_of_canvases[pulse_order].Divide(2,2); // } // if (verbosityLevel > 2) cout << "...done" << endl; //} //void //DeletePulseCanvas( // TCanvas* array_of_canvases, // unsigned int max_pulse_order // ) //{ // for ( // unsigned int pulse_order = 1; // pulse_order <= max_pulse_order; // pulse_order++ // ) // { // delete array_of_canvases[pulse_order]; // } //} void DeletePixelCanvases( int verbosityLevel ) { if (verbosityLevel > 2) cout << endl << "...delete pixel Canvas" ; for (int pulse_order = 0; pulse_order < MAX_PULS_ORDER; pulse_order ++) { delete cgpPixelPulses[pulse_order]; } } //---------------------------------------------------------------------------- //---------------------------------------------------------------------------- //void SetHistogramAttributes( // const char* histoType, // int pulse_order, // histAttrib_t* histoAttrib, // int max_ampl_of_first_pulse // ) //{ // histoAttrib[pulse_order].yMax = (max_ampl_of_first_pulse * pulse_order); // histoAttrib[pulse_order].yMin = (histoAttrib[pulse_order].yMax/2); // histoAttrib[pulse_order].name = "h"; // histoAttrib[pulse_order].name += histoType; // histoAttrib[pulse_order].name = "_"; // histoAttrib[pulse_order].name += pulse_order; // histoAttrib[pulse_order].title += histoType; // histoAttrib[pulse_order].title += " with Amplitude around "; // histoAttrib[pulse_order].title += histoAttrib[pulse_order].yMax; // histoAttrib[pulse_order].title += " mV"; // histoAttrib[pulse_order].xTitle += "Timeslices [a.u.]"; // histoAttrib[pulse_order].yTitle += "Amplitude [mV]"; //}