Changeset 13422 for fact/tools


Ignore:
Timestamp:
04/24/12 08:22:49 (13 years ago)
Author:
Jens Buss
Message:
now with median, mean and max of slice
File:
1 edited

Legend:

Unmodified
Added
Removed

  • fact/tools/rootmacros/FPulseTemplate.C

    r13378 r13422  
    3535#include <stdint.h>
    3636#include <cstdio>
     37#include <iostream>
     38#include <fstream>
     39using namespace std;
    3740
    3841#define NPIX  1440
    … …  
    9497
    9598typedef struct{
    96   float maxAmpl;
    97   float countOfMax;
     99  double maxAmpl;
     100  int countOfMax;
    98101  } OverlayMaximum;
    99102
    … …  
    146149TH1D*       hProjPeak = NULL;
    147150TH1F*       hTesthisto = NULL;
     151TH2F*       hTesthisto2 = NULL;
    148152
    149153//Pixelwise Histograms
    150154TH2F*       hPixelOverlay[MAX_PULS_ORDER]= {NULL};//histogrammm for overlay of detected Peaks
     155TH2F*       hPixelEdgeOverlay[MAX_PULS_ORDER] = {NULL};
    151156TProfile*   hPixelProfile[MAX_PULS_ORDER] = {NULL};//histogrammm for Profile of detected Peaks
     157TProfile*   hPixelProfile2[MAX_PULS_ORDER] = {NULL};//histogrammm for Profile of detected Peaks
    152158TH1F*       hPixelMax[MAX_PULS_ORDER] = {NULL};
    153 TH2F*       hPixelEdgeOverlay[MAX_PULS_ORDER] = {NULL};
     159TH1F*       hPixelMedian[MAX_PULS_ORDER] = {NULL};
     160TH1F*       hPixelMean[MAX_PULS_ORDER] = {NULL};
    154161
    155162//All Pixel Histograms
    156163TH2F*       hAllPixelOverlay[MAX_PULS_ORDER] = {NULL};
    157164TProfile*   hAllPixelProfile[MAX_PULS_ORDER] = {NULL};
     165TProfile*   hAllPixelProfile2[MAX_PULS_ORDER] = {NULL};
    158166TH1F*       hAllPixelMax[MAX_PULS_ORDER]     = {NULL};
    159 TH1F*       hAllPixelMaxGaus[MAX_PULS_ORDER] = {NULL};
     167TH1F*       hAllPixelMedian[MAX_PULS_ORDER] = {NULL};
     168TH1F*       hAllPixelMean[MAX_PULS_ORDER] = {NULL};
     169TH2F*       hAllPixelEdgeOverlay[MAX_PULS_ORDER] = {NULL};
    160170
    161171//Histogram Parameters
    … …  
    183193void DeletePixelHistos(int );
    184194void SaveHistograms( const char*, const char*, TObjArray*, int );
    185 void FillHistograms(vector<Region>*, int, int);
     195void FillHistograms(vector<Region>*, int, int, int);
    186196void PlotPulseShapeOfMaxPropability(unsigned int, int, int);
     197
    187198void DrawPulseHistograms(int, int);
    188199void FitMaxPropabilityPuls( TH1F* , int );
    189200void DrawMaximumHistograms( int, int );
    190201void DrawTestHistograms( int);
     202Double_t MedianOfH1( TH1 *h1 );
     203void PlotMedianEachSliceOfPulse(unsigned int, int, int);
     204
    191205
    192206void UpdateCanvases( int, int);
    … …  
    198212TString CreateSubDirName(int); //creates a String containing the path to the subdirectory in root file
    199213TString CreateSubDirName(const char* );
     214void WritePixelTemplateToCsv( TString, const char*, const char*, int, int);
     215void WriteAllPixelTemplateToCsv( TString, const char*, const char*, int);
    200216//void SetHistogrammSettings( const char*, int, int , int);
    201217//void CreatePulseCanvas( TCanvas* , unsigned int, const char* , const char*, const char*, int);
    … …  
    212228  const char*   drsfilename         = "../data/2011/11/09/20111109_003.drs.fits.gz",
    213229  const char*   OutRootFileName     = "../analysis/FPulseTemplate/20111109_006/20111109_006.pulses.root",
    214   bool          ProduceGraphic      = true,
    215   int           refresh_rate        = 500,   //refresh rate for canvases
     230  const char*   OutPutPath          = "../analysis/FPulseTemplate/20111109_006/Templates/",
     231  bool          ProduceGraphic      = false,
     232  int           refresh_rate        = 500,      //refresh rate for canvases
    216233  bool          spikeDebug          = false,
    217   bool          debugPixel          = true,
     234  bool          debugPixel          = false,
    218235  bool          fitdata             = false,
    219   int           verbosityLevel      = 1, // different verbosity levels can be implemented here
    220   int           firstevent          = 272,
    221   int           nevents             = 500,
    222   int           firstpixel          = 1100,
    223   int           npixel              = 1,
    224   int           AmplWindowWidth     = 14,  //Width of Window for selection of pulses to histograms
     236  int           verbosityLevel      = 0,        // different verbosity levels can be implemented here
     237  int           firstevent          = 0,
     238  int           nevents             = -1,
     239  int           firstpixel          = 0,
     240  int           npixel              = -1,
     241  int           AmplWindowWidth     = 14,       //Width of Window for selection of pulses to histograms
    225242  float         GainMean           = 8,
    226   float         BSLMean            = -1, //4 Histogramms will be drawn, decide how far pulsehights differ from eachother
     243  float         BSLMean            = -1,        //4 Histogramms will be drawn, decide how far pulsehights differ from eachother
    227244  int           avg1                = 8,
    228245  int           avg2                = 8,
    … …  
    360377    for ( int pixel = firstpixel; pixel < firstpixel + npixel; pixel++ )
    361378    {
    362         if (verbosityLevel > 0)
     379        if (verbosityLevel >= 0)
    363380        {
    364381            cout << "------------------------------------------------" << endl
    … …  
    440457            removeRegionWithMaxOnEdge( *pZXings, 2);
    441458            removeRegionOnFallingEdge( *pZXings, 100);
    442             findTimeOfHalfMaxLeft(*pZXings, Vslide, gBSLMean, 5, 10, verbosityLevel );
     459            findTimeOfHalfMaxLeft(*pZXings, Vcorr, gBSLMean, 5, 10, verbosityLevel );
    443460            if (verbosityLevel > 2) cout << "...done" << endl;
    444461
    … …  
    446463// Fill Overlay Histos
    447464//-----------------------------------------------------------------------------
    448         FillHistograms(pZXings, AmplWindowWidth, verbosityLevel );
     465        FillHistograms(pZXings, AmplWindowWidth, ev, verbosityLevel );
    449466
    450467//-----------------------------------------------------------------------------
    … …  
    557574//-------------------------------------
    558575
     576
     577
     578
    559579//-------------------------------------
    560580// Histogramms of Maximas in Overlay Spectra
    … …  
    567587                 );
    568588
     589         PlotMedianEachSliceOfPulse(
     590                     MAX_PULS_ORDER,
     591                     fitdata,
     592                     verbosityLevel
     593                 );
     594
     595         WritePixelTemplateToCsv(
     596                     OutPutPath,
     597                     "PulseTemplate_PointSet",
     598                     "Maximum",
     599                     pixel,
     600                     verbosityLevel
     601                     );
    569602
    570603       if (verbosityLevel > 2) cout << "...done" << endl;
    … …  
    599632       {
    600633           //Process gui events asynchronously during input
     634           cout << endl;
    601635           TTimer timer("gSystem->ProcessEvents();", 50, kFALSE);
    602636           timer.TurnOn();
    … …  
    647681    }
    648682
     683    WriteAllPixelTemplateToCsv(
     684                OutPutPath,
     685                "PulseTemplate_PointSet",
     686                "Maximum",
     687                verbosityLevel
     688                );
     689
    649690    DeletePixelCanvases( verbosityLevel );
    650691    return( 0 );
    … …  
    662703
    663704
    664 void BookPixelHistos( int verbosityLevel)
     705void
     706BookPixelHistos( int verbosityLevel)
    665707{
    666708    if (verbosityLevel > 2) cout << endl << "...book pixel histograms" << endl;
    … …  
    686728                    512,
    687729                    -55.5,
    688                     300.5
     730                    200.5
    689731                    );
    690732        hPixelOverlay[order]->SetAxisRange(
    … …  
    722764        hList->Add( hPixelMax[order] );
    723765
    724     //------------------------------------------------------------------------
     766//------------------------------------------------------------------------
     767    //        SetHistogramAttributes(
     768    //                    "PeakMaximum",
     769    //                    order,
     770    //                    gMaximumAttrib,
     771    //                    MaxAmplOfFirstPulse
     772    //                    );
     773            histoname = "hPixelMedian";
     774            histoname += order;
     775            if (verbosityLevel > 3) cout << "...booking " << histoname << endl;
     776            hPixelMedian[order] = new TH1F(
     777                        histoname,
     778                        "Median of each slice in overlay plot",
     779                        gPixelOverlayXaxisLeft + gPixelOverlayXaxisRight ,
     780                        (-1*gPixelOverlayXaxisLeft)-0.5,
     781                        gPixelOverlayXaxisRight-0.5
     782                        );
     783            hPixelMedian[order]->SetAxisRange(
     784                        gBSLMean - 5,
     785                        (gGainMean*(order+1)) + 10,
     786                        "Y");
     787            hPixelMedian[order]->SetLineColor(kRed);
     788            hPixelMedian[order]->GetXaxis()->SetTitle( "Timeslices [a.u.]" );
     789            hPixelMedian[order]->GetYaxis()->SetTitle( "Amplitude [mV]" );
     790            hList->Add( hPixelMedian[order] );
     791
     792//------------------------------------------------------------------------
     793    //        SetHistogramAttributes(
     794    //                    "PeakMaximum",
     795    //                    order,
     796    //                    gMaximumAttrib,
     797    //                    MaxAmplOfFirstPulse
     798    //                    );
     799            histoname = "hPixelMean";
     800            histoname += order;
     801            if (verbosityLevel > 3) cout << "...booking " << histoname << endl;
     802            hPixelMean[order] = new TH1F(
     803                        histoname,
     804                        "Mean of each slice in overlay plot",
     805                        gPixelOverlayXaxisLeft + gPixelOverlayXaxisRight ,
     806                        (-1*gPixelOverlayXaxisLeft)-0.5,
     807                        gPixelOverlayXaxisRight-0.5
     808                        );
     809            hPixelMean[order]->SetAxisRange(
     810                        gBSLMean - 5,
     811                        (gGainMean*(order+1)) + 10,
     812                        "Y");
     813            hPixelMean[order]->SetLineColor(kBlue);
     814            hPixelMean[order]->GetXaxis()->SetTitle( "Timeslices [a.u.]" );
     815            hPixelMean[order]->GetYaxis()->SetTitle( "Amplitude [mV]" );
     816            hList->Add( hPixelMean[order] );
     817
     818//------------------------------------------------------------------------
    725819//        SetHistogramAttributes(
    726820//                    "PeakMaxGaus",
    … …  
    740834                    512,
    741835                    -55.5,
    742                     300.5
     836                    200.5
    743837                    );
    744838
    … …  
    771865                    300.5,                                          //yup
    772866                    "s");                                           //option
    773         hPixelMax[order]->SetAxisRange(
     867        hPixelProfile[order]->SetAxisRange(
    774868                    gBSLMean - 5,
    775869                    (gGainMean*(order+1)) + 10,
    … …  
    779873        //hPixelProfile->SetBit(TH2F::kCanRebin);
    780874        hList->Add( hPixelProfile[order] );
     875
     876
     877    //------------------------------------------------------------------------
     878    //        SetHistogramAttributes(
     879    //                    "PeakProfile",
     880    //                    order,
     881    //                    gProfileAttrib,
     882    //                    MaxAmplOfFirstPulse
     883    //                    );
     884            histoname = "hPixelProfile2";
     885            histoname += order;
     886            if (verbosityLevel > 3) cout << "...booking " << histoname << endl;
     887            hPixelProfile2[order] = new TProfile(
     888                        histoname,
     889                        "Mean value of each slice in overlay plot (Tprofile)",
     890                        gPixelOverlayXaxisLeft + gPixelOverlayXaxisRight ,//nbinsx
     891                        (-1*gPixelOverlayXaxisLeft)-0.5,                 //xlow
     892                        gPixelOverlayXaxisRight-0.5 ,                    //xup
     893        //                512,                                            //nbinsy
     894                        -55.5,                                          //ylow
     895                        300.5,                                          //yup
     896                        "s");                                           //option
     897            hPixelProfile2[order]->SetLineColor(kRed);
     898            hPixelProfile2[order]->SetAxisRange(
     899                        gBSLMean - 5,
     900                        (gGainMean*(order+1)) + 10,
     901                        "Y");
     902            hPixelProfile2[order]->GetXaxis()->SetTitle( "Timeslices [a.u.]" );
     903            hPixelProfile2[order]->GetYaxis()->SetTitle( "Amplitude [mV]" );
     904            //hPixelProfile->SetBit(TH2F::kCanRebin);
     905
     906
     907
     908
    781909    }
    782910    if (verbosityLevel > 2) cout << "...done" << endl;
    … …  
    795923        if (verbosityLevel > 3) cout << endl << "...deleting hPixelMax" << order;
    796924        delete hPixelMax[order];
     925        if (verbosityLevel > 3) cout << endl << "...deleting hPixelMedian" << order;
     926        delete hPixelMedian[order];
     927        if (verbosityLevel > 3) cout << endl << "...deleting hPixelMean" << order;
     928        delete hPixelMean[order];
    797929        if (verbosityLevel > 3) cout << endl << "...deleting hPixelEdgeOverlay" << order;
    798930        delete hPixelEdgeOverlay[order];
    799931        if (verbosityLevel > 3) cout << endl << "...deleting hPixelProfile" << order;
    800932        delete hPixelProfile[order];
     933        if (verbosityLevel > 3) cout << endl << "...deleting hPixelProfile2" << order;
     934        delete hPixelProfile2[order];
    801935    }
    802936    if (verbosityLevel > 3) cout << endl << "...deleting hList";
    … …  
    814948                "hTesthisto",
    815949                "Deviation of rising edge and maximum",
    816                 200,
    817                 -300,
    818                 300
     950                600,
     951                -10.1,
     952                10.1
    819953                );
    820954    hTesthisto->GetXaxis()->SetTitle( "Timeslices [a.u.]" );
    821955    hTesthisto->GetYaxis()->SetTitle( "counts" );
     956
     957    hTesthisto2 = new TH2F (
     958                "hTesthisto2",
     959                "Deviation of rising edge and maximum by event #",
     960                gNEvents,
     961                250,
     962                800,
     963                600,
     964                -10.1,
     965                10.1
     966                );
     967//    hTesthisto2->GetXaxis()->SetTitle( "Timeslices [a.u.]" );
     968//    hTesthisto2->GetYaxis()->SetTitle( "counts" );
     969//    hTesthisto2->SetMarkerStyle( 2 );
    822970
    823971    if (verbosityLevel > 2) cout << "...done" << endl;
    … …  
    9101058        hAllPixelList->Add( hAllPixelMax[order] );
    9111059
    912     //------------------------------------------------------------------------
     1060//------------------------------------------------------------------------
    9131061//        SetHistogramAttributes(
    914 //                    "AllPixelPeakMaxGaus",
     1062//                    "PeakMaximum",
     1063//                    order,
     1064//                    gMaximumAttrib,
     1065//                    MaxAmplOfFirstPulse
     1066//                    );
     1067        histoname = "hAllPixelMedian";
     1068        histoname += order;
     1069        if (verbosityLevel > 3) cout << "...booking " << histoname << endl;
     1070        hAllPixelMedian[order] = new TH1F(
     1071                    histoname,
     1072                    "Median of each slice in overlay plot for All Pixels",
     1073                    gPixelOverlayXaxisLeft + gPixelOverlayXaxisRight ,
     1074                    (-1*gPixelOverlayXaxisLeft)-0.5,
     1075                    gPixelOverlayXaxisRight-0.5
     1076                    );
     1077        hAllPixelMedian[order]->SetAxisRange(
     1078                    gBSLMean - 5,
     1079                    (gGainMean*(order+1)) + 10,
     1080                    "Y");
     1081        hAllPixelMedian[order]->SetLineColor(kRed);
     1082        hAllPixelMedian[order]->GetXaxis()->SetTitle( "Timeslices [a.u.]" );
     1083        hAllPixelMedian[order]->GetYaxis()->SetTitle( "Amplitude [mV]" );
     1084        hAllPixelList->Add( hAllPixelMedian[order] );
     1085        if (verbosityLevel > 3) cout << "...BREAKPOINT in " << histoname << endl;
     1086//------------------------------------------------------------------------
     1087//        SetHistogramAttributes(
     1088//                    "PeakMaximum",
     1089//                    order,
     1090//                    gMaximumAttrib,
     1091//                    MaxAmplOfFirstPulse
     1092//                    );
     1093        histoname = "hAllPixelMean";
     1094        histoname += order;
     1095        if (verbosityLevel > 3) cout << "...booking " << histoname << endl;
     1096        hAllPixelMean[order] = new TH1F(
     1097                    histoname,
     1098                    "Mean of each slice in overlay plot for All Pixels",
     1099                    gPixelOverlayXaxisLeft + gPixelOverlayXaxisRight ,
     1100                    (-1*gPixelOverlayXaxisLeft)-0.5,
     1101                    gPixelOverlayXaxisRight-0.5
     1102                    );
     1103        hAllPixelMean[order]->SetAxisRange(
     1104                    gBSLMean - 5,
     1105                    (gGainMean*(order+1)) + 10,
     1106                    "Y");
     1107        hAllPixelMean[order]->SetLineColor(kBlue);
     1108        hAllPixelMean[order]->GetXaxis()->SetTitle( "Timeslices [a.u.]" );
     1109        hAllPixelMean[order]->GetYaxis()->SetTitle( "Amplitude [mV]" );
     1110        hAllPixelList->Add( hAllPixelMean[order] );
     1111
     1112//------------------------------------------------------------------------
     1113//        SetHistogramAttributes(
     1114//                    "PeakMaxGaus",
    9151115//                    order,
    9161116//                    gMaxGausAttrib,
    9171117//                    MaxAmplOfFirstPulse
    9181118//                    );
    919         histoname = "hAllPixelMaxGaus";
     1119        histoname = "hAllPixelEdgeOverlay";
    9201120        histoname += order;
    9211121        if (verbosityLevel > 3) cout << "...booking " << histoname << endl;
    922         hAllPixelMaxGaus[order] = new TH1F(
     1122        hAllPixelEdgeOverlay[order] = new TH2F(
    9231123                    histoname,
    924                     "Mean value of each slice in overlay plot for All Pixels",
     1124                    "Overlay at rising edge of detected pulses of one pulse order for All Pixels ",
    9251125                    gPixelOverlayXaxisLeft + gPixelOverlayXaxisRight ,
    9261126                    (-1*gPixelOverlayXaxisLeft)-0.5,
    927                     gPixelOverlayXaxisRight-0.5
     1127                    gPixelOverlayXaxisRight-0.5 ,
     1128                    512,
     1129                    -55.5,
     1130                    200.5
    9281131                    );
    929         hAllPixelMaxGaus[order]->SetAxisRange(
     1132
     1133        hAllPixelEdgeOverlay[order]->SetAxisRange(
    9301134                    gBSLMean - 5,
    9311135                    (gGainMean*(order+1)) + 10,
    9321136                    "Y");
    933         hAllPixelMaxGaus[order]->GetXaxis()->SetTitle( "Timeslices [a.u.]" );
    934         hAllPixelMaxGaus[order]->GetYaxis()->SetTitle( "Amplitude [mV]" );
    935         hAllPixelList->Add( hAllPixelMaxGaus[order] );
    936 
    937     //------------------------------------------------------------------------
     1137        hAllPixelEdgeOverlay[order]->GetXaxis()->SetTitle( "Timeslices [a.u.]" );
     1138        hAllPixelEdgeOverlay[order]->GetYaxis()->SetTitle( "Amplitude [mV]" );
     1139        hAllPixelList->Add( hAllPixelEdgeOverlay[order] );
     1140
     1141//------------------------------------------------------------------------
    9381142//        SetHistogramAttributes(
    9391143//                    "AllPixelPeakProfile",
    … …  
    9641168        hAllPixelList->Add( hAllPixelProfile[order] );
    9651169
     1170//------------------------------------------------------------------------
     1171//        SetHistogramAttributes(
     1172//                    "AllPixelPeakProfile",
     1173//                    order,
     1174//                    gProfileAttrib,
     1175//                    MaxAmplOfFirstPulse
     1176//                    );
     1177        histoname = "hAllPixelProfile2";
     1178        histoname += order;
     1179        if (verbosityLevel > 3) cout << "...booking " << histoname << endl;
     1180        hAllPixelProfile2[order] = new TProfile(
     1181                    histoname,
     1182                    "Mean value of each slice in overlay plot for All Pixels (Tprofile)",
     1183                    gPixelOverlayXaxisLeft + gPixelOverlayXaxisRight ,//nbinsx
     1184                    (-1*gPixelOverlayXaxisLeft)-0.5,                 //xlow
     1185                    gPixelOverlayXaxisRight-0.5 ,                    //xup
     1186    //                512,                                            //nbinsy
     1187                    -55.5,                                          //ylow
     1188                    300.5,                                          //yup
     1189                    "s");                                           //option
     1190        hAllPixelProfile2[order]->SetAxisRange(
     1191                    gBSLMean - 5,
     1192                    (gGainMean*(order+1)) + 10,
     1193                    "Y");
     1194        hAllPixelProfile2[order]->SetLineColor(kRed);
     1195        hAllPixelProfile2[order]->GetXaxis()->SetTitle( "Timeslices [a.u.]" );
     1196        hAllPixelProfile2[order]->GetYaxis()->SetTitle( "Amplitude [mV]" );
     1197        //hPixelProfile->SetBit(TH2F::kCanRebin);
     1198        hAllPixelList->Add( hAllPixelProfile2[order] );
     1199
    9661200
    9671201
    … …  
    10771311        vector<Region>* pZXings,
    10781312        int AmplWindowWidth,
     1313        int eventNumber,
    10791314        int verbosityLevel
    10801315        )
    … …  
    11031338
    11041339        if (Left < 0) Left =0;
    1105         if (Right > (int)Vcorr.size() ) Right=Vcorr.size() ;
    1106 
    1107         hTesthisto->Fill( reg->maxPos - reg->halfRisingEdgePos ) ;
     1340        if (Right > (int)Ameas.size() ) Right=Ameas.size() ;
     1341        if (EdgeLeft < 0) EdgeLeft =0;
     1342        if (EdgeRight > (int)Ameas.size() ) EdgeRight=Ameas.size() ;
     1343
     1344
     1345        hTesthisto->Fill( reg->slopeOfRisingEdge ) ;
     1346        hTesthisto2->Fill( eventNumber, reg->slopeOfRisingEdge ) ;
    11081347
    11091348        if (verbosityLevel > 2) cout << endl << "\t...choosing Histogram" ;
    … …  
    11321371            if (verbosityLevel > 2) cout << "...filling" ;
    11331372            for ( int pos = Left; pos < Right; pos++){
     1373//                if ();
    11341374                hPixelOverlay[order_of_pulse]->Fill( pos - (reg->maxPos), Ameas[pos]) ;
    11351375                hPixelProfile[order_of_pulse]->Fill( pos - (reg->maxPos), Ameas[pos]) ;
    … …  
    11401380//                cout << endl << "###filling edge histo###" << endl;
    11411381                hPixelEdgeOverlay[order_of_pulse]->Fill( pos - (reg->halfRisingEdgePos), Ameas[pos]) ;
     1382                hPixelProfile2[order_of_pulse]->Fill( pos - (reg->halfRisingEdgePos), Ameas[pos]) ;
     1383                hAllPixelEdgeOverlay[order_of_pulse]->Fill( pos - (reg->halfRisingEdgePos), Ameas[pos]) ;
     1384                hAllPixelProfile2[order_of_pulse]->Fill( pos - (reg->halfRisingEdgePos), Ameas[pos]) ;
    11421385//                cout << endl << "######" << endl;
    11431386            }
    … …  
    11611404    cgpTestHistos->cd(1);
    11621405    hTesthisto->Draw();
     1406    cgpTestHistos->cd(2);
     1407    hTesthisto2->Draw("BOX");
    11631408}
    11641409// end of DrawTestHistograms
    … …  
    11751420        cgpPixelPulses[pulse_order]->cd(1);
    11761421        hPixelOverlay[pulse_order]->Draw("COLZ");
     1422        cgpPixelPulses[pulse_order]->cd(3);
     1423        hPixelProfile[pulse_order]->Draw("COLZ");
     1424        hPixelProfile2[pulse_order]->Draw("SAME");
     1425
    11771426        cgpPixelPulses[pulse_order]->cd(2);
    1178         hPixelProfile[pulse_order]->Draw("COLZ");
    1179         cgpPixelPulses[pulse_order]->cd(3);
    11801427        hPixelEdgeOverlay[pulse_order]->Draw("COLZ");
    11811428//        cgpPixelPulses[pulse_order]->cd(4);
    … …  
    11851432        hAllPixelOverlay[pulse_order]->Draw("COLZ");
    11861433        cgpAllPixelPulses[pulse_order]->cd(2);
     1434        hAllPixelEdgeOverlay[pulse_order]->Draw("COLZ");
     1435
     1436        cgpAllPixelPulses[pulse_order]->cd(3);
    11871437        hAllPixelProfile[pulse_order]->Draw("COLZ");
     1438        hAllPixelProfile2[pulse_order]->Draw("SAME");
    11881439
    11891440    }
    … …  
    12021453    {
    12031454        cgpPixelPulses[pulse_order]->cd(4);
    1204         hPixelMax[pulse_order]->Draw("LF2");
     1455        hPixelMax[pulse_order]->Draw();
     1456        hPixelMedian[pulse_order]->Draw("SAME");
     1457        hPixelMean[pulse_order]->Draw("SAME");
    12051458
    12061459        cgpAllPixelPulses[pulse_order]->cd(4);
    12071460        hAllPixelMax[pulse_order]->Draw();
     1461        hAllPixelMedian[pulse_order]->Draw("SAME");
     1462        hAllPixelMean[pulse_order]->Draw("SAME");
    12081463//        cgpAllPixelPulses[pulse_order]->cd(3);
    12091464     //   hAllPixelMaxGaus[pulse_order]->Draw("COLZ");
    … …  
    12361491
    12371492void
     1493PlotMedianEachSliceOfPulse(
     1494        unsigned int    max_pulse_order,
     1495        int             fitdata,
     1496        int             verbosityLevel
     1497        )
     1498{
     1499    if (verbosityLevel > 2) cout << endl
     1500                                 << "...calculating pulse shape of slice's Median"
     1501                                 << endl;
     1502    for (unsigned int pulse_order = 0 ; pulse_order < max_pulse_order ; pulse_order ++)
     1503    {
     1504        if (verbosityLevel > 2) cout << "\t...calculation of "
     1505                                     << "pulse order " << pulse_order;
     1506
     1507        Int_t nbins = hPixelOverlay[pulse_order]->GetXaxis()->GetNbins();
     1508
     1509        for (Int_t TimeSlice=1;TimeSlice<=nbins;TimeSlice++) {
     1510           TH1 *h1 = hPixelOverlay[pulse_order]->ProjectionY("",TimeSlice,TimeSlice);
     1511           float median = MedianOfH1(h1);
     1512           float mean = h1->GetMean();
     1513           if (verbosityLevel > 2) printf("Median of Slice %d, Median=%g\n",TimeSlice,median);
     1514           delete h1;
     1515           hPixelMedian[pulse_order]->SetBinContent(TimeSlice, median );
     1516           hPixelMean[pulse_order]->SetBinContent(TimeSlice, mean );
     1517           hAllPixelMedian[pulse_order]->SetBinContent(TimeSlice, median );
     1518           hAllPixelMean[pulse_order]->SetBinContent(TimeSlice, mean );
     1519        }
     1520
     1521        if (verbosityLevel > 2) cout << "\t...done" << endl;
     1522
     1523        if (fitdata)
     1524        {
     1525           FitMaxPropabilityPuls(
     1526                    hPixelMedian[pulse_order],
     1527                    verbosityLevel
     1528                    );
     1529        }
     1530    }
     1531}
     1532// end of PlotMedianEachSliceOfPulse
     1533//----------------------------------------------------------------------------
     1534
     1535Double_t MedianOfH1 (
     1536        TH1*            h1
     1537        )
     1538{
     1539   //compute the median for 1-d histogram h1
     1540   Int_t nbins = h1->GetXaxis()->GetNbins();
     1541   Double_t *x = new Double_t[nbins];
     1542   Double_t *y = new Double_t[nbins];
     1543   for (Int_t i=0;i<nbins;i++) {
     1544      x[i] = h1->GetXaxis()->GetBinCenter(i+1);
     1545      y[i] = h1->GetBinContent(i+1);
     1546   }
     1547   Double_t median = TMath::Median(nbins,x,y);
     1548   delete [] x;
     1549   delete [] y;
     1550   return median;
     1551}
     1552// end of PlotMedianEachSliceOfPulse
     1553//----------------------------------------------------------------------------
     1554
     1555void
    12381556PlotPulseShapeOfMaxPropability(
    12391557        unsigned int    max_pulse_order,
    … …  
    12471565    for (unsigned int pulse_order = 0 ; pulse_order < max_pulse_order ; pulse_order ++)
    12481566    {
    1249         if (verbosityLevel > 2) cout << "\t...calculating of "
     1567        if (verbosityLevel > 2) cout << "\t...calculation of "
    12501568                                     << "pulse order " << pulse_order;
    12511569        //  vector max_value_of to number of timeslices in Overlay Spectra
    1252         vector<OverlayMaximum> max_value_of;
    1253         max_value_of.resize(hPixelOverlay[pulse_order]->GetNbinsX());
    1254 
    1255         for (int TimeSlice = 0;
    1256              TimeSlice <= hPixelOverlay[pulse_order]->GetNbinsX();
    1257              TimeSlice++ )
     1570        float max_value_of_slice;
     1571
     1572        Int_t nbins = hPixelOverlay[pulse_order]->GetXaxis()->GetNbins();
     1573        if (verbosityLevel > 2) cout << "...generating projections" << endl;
     1574        for (Int_t TimeSlice=1;TimeSlice<=nbins;TimeSlice++)
    12581575        {
    1259             histotitle = "hproj_py";    //generate title of histogram of which MaxVal
    1260             histotitle += pulse_order ; //will be calculated
    1261 
     1576            if (verbosityLevel > 2) cout << "...Timeslice: " << TimeSlice;
    12621577            //put maximumvalue of every Y-projection of every timeslice into vector
    1263             hProjPeak = hPixelOverlay[pulse_order]->ProjectionY(histotitle, TimeSlice, TimeSlice, "");
    1264             max_value_of[ TimeSlice ].maxAmpl = (hProjPeak->GetMaximumBin() * 0.5) - 3.5;
    1265             max_value_of[ TimeSlice ].countOfMax = hProjPeak->GetBinContent( hProjPeak->GetMaximumBin() );
    1266             hPixelMax[pulse_order]->SetBinContent(TimeSlice, max_value_of[ TimeSlice ].maxAmpl );
     1578            TH1D *h1 =  hPixelOverlay[pulse_order]->ProjectionY( "", TimeSlice,TimeSlice);
     1579
     1580            max_value_of_slice = h1->GetBinCenter( h1->GetMaximumBin() );
     1581
     1582            if (verbosityLevel > 2) cout << " with max value "
     1583                                         << max_value_of_slice << endl;
     1584            hPixelMax[pulse_order]->SetBinContent(TimeSlice, max_value_of_slice );
     1585            hAllPixelMax[pulse_order]->SetBinContent(TimeSlice, max_value_of_slice );
     1586            delete h1;
    12671587        }
    12681588
    … …  
    12821602void
    12831603FitMaxPropabilityPuls(
    1284         TH1F*   hMaximumTemp,
    1285         int     verbosityLevel
     1604        TH1F*       hMaximumTemp,
     1605        int         verbosityLevel
    12861606        )
    12871607    {
    … …  
    14151735//    histoAttrib[pulse_order].yTitle += "Amplitude [mV]";
    14161736//}
     1737
     1738void WritePixelTemplateToCsv(
     1739        TString         path,
     1740        const char*     csv_file_name,
     1741        const char*     overlay_method,
     1742        int             pixel,
     1743        int             verbosityLevel
     1744        )
     1745{
     1746    path += csv_file_name;
     1747    path += "_";
     1748    path += pixel;
     1749    path += ".csv";
     1750
     1751    Int_t nbins = hPixelMax[0]->GetXaxis()->GetNbins();
     1752    if (verbosityLevel > 0)
     1753    {
     1754        cout << "writing point-set to csv file: " ;
     1755        cout << path << endl;
     1756        cout << "...opening file" << endl;
     1757    }
     1758    if (verbosityLevel > 2) cout << "...number of bins " << nbins << endl;
     1759    ofstream out;
     1760    out.open( path );
     1761    out << "### point-set of a single photon pulse template" << endl;
     1762    out << "### template determined with pulse overlay at: "
     1763        << overlay_method << endl;
     1764    out << "### Slice's Amplitude determined by calculating the " << endl
     1765        << "### value of maximum propability of slice -> Amplitude1 " << endl
     1766        << "### mean of slice -> Amplitude2 " << endl
     1767        << "### median of slice -> Amplitude3 " << endl
     1768        << "### for each slice" << endl;
     1769    out << "### Pixel number (CHid): " << pixel << endl
     1770        << endl;
     1771
     1772    out << "time [slices],Amplitude1 [mV],Amplitude2 [mV],Amplitude3 [mV]" << endl;
     1773
     1774    for (int TimeSlice=1;TimeSlice<=nbins;TimeSlice++)
     1775    {
     1776        out << TimeSlice << "," ;
     1777        out << hPixelMax[0]->GetBinContent(TimeSlice) << ",";
     1778        out << hPixelMean[0]->GetBinContent(TimeSlice) << ",";
     1779        out << hPixelMedian[0]->GetBinContent(TimeSlice) << endl;
     1780    }
     1781    out.close();
     1782    if (verbosityLevel > 0) cout << "...file closed" << endl;
     1783}
     1784
     1785void WriteAllPixelTemplateToCsv(
     1786        TString         path,
     1787        const char*     csv_file_name,
     1788        const char*     overlay_method,
     1789        int             verbosityLevel
     1790        )
     1791{
     1792    path += csv_file_name;
     1793    path += "_AllPixel.csv";
     1794
     1795    Int_t nbins = hAllPixelMax[0]->GetXaxis()->GetNbins();
     1796    if (verbosityLevel > 0)
     1797    {
     1798        cout << "writing point-set to csv file: " ;
     1799        cout << path << endl;
     1800        cout << "...opening file" << endl;
     1801    }
     1802    if (verbosityLevel > 2) cout << "...number of bins " << nbins << endl;
     1803    ofstream out;
     1804    out.open( path );
     1805    out << "### point-set of a single photon pulse template" << endl;
     1806    out << "### template determined with pulse overlay at: "
     1807        << overlay_method << "of all Pixels";
     1808    out << "### Slice's Amplitude determined by calculating the " << endl
     1809        << "### value of maximum propability of slice -> Amplitude1 " << endl
     1810        << "### mean of slice -> Amplitude2 " << endl
     1811        << "### median of slice -> Amplitude3 " << endl
     1812        << "### for each slice" << endl
     1813        << endl;
     1814
     1815    out << "time [slices],Amplitude1 [mV],Amplitude2 [mV],Amplitude3 [mV]" << endl;
     1816
     1817    for (int TimeSlice=1;TimeSlice<=nbins;TimeSlice++)
     1818    {
     1819        out << TimeSlice << "," ;
     1820        out << hAllPixelMax[0]->GetBinContent(TimeSlice) << ",";
     1821        out << hAllPixelMean[0]->GetBinContent(TimeSlice) << ",";
     1822        out << hAllPixelMedian[0]->GetBinContent(TimeSlice) << endl;
     1823    }
     1824    out.close();
     1825    if (verbosityLevel > 0) cout << "...file closed" << endl;
     1826}
Note: See TracChangeset for help on using the changeset viewer.