Changeset 13422
- Timestamp:
- 04/24/12 08:22:49 (13 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
fact/tools/rootmacros/FPulseTemplate.C
r13378 r13422 35 35 #include <stdint.h> 36 36 #include <cstdio> 37 #include <iostream> 38 #include <fstream> 39 using namespace std; 37 40 38 41 #define NPIX 1440 … … 94 97 95 98 typedef struct{ 96 floatmaxAmpl;97 float countOfMax;99 double maxAmpl; 100 int countOfMax; 98 101 } OverlayMaximum; 99 102 … … 146 149 TH1D* hProjPeak = NULL; 147 150 TH1F* hTesthisto = NULL; 151 TH2F* hTesthisto2 = NULL; 148 152 149 153 //Pixelwise Histograms 150 154 TH2F* hPixelOverlay[MAX_PULS_ORDER]= {NULL};//histogrammm for overlay of detected Peaks 155 TH2F* hPixelEdgeOverlay[MAX_PULS_ORDER] = {NULL}; 151 156 TProfile* hPixelProfile[MAX_PULS_ORDER] = {NULL};//histogrammm for Profile of detected Peaks 157 TProfile* hPixelProfile2[MAX_PULS_ORDER] = {NULL};//histogrammm for Profile of detected Peaks 152 158 TH1F* hPixelMax[MAX_PULS_ORDER] = {NULL}; 153 TH2F* hPixelEdgeOverlay[MAX_PULS_ORDER] = {NULL}; 159 TH1F* hPixelMedian[MAX_PULS_ORDER] = {NULL}; 160 TH1F* hPixelMean[MAX_PULS_ORDER] = {NULL}; 154 161 155 162 //All Pixel Histograms 156 163 TH2F* hAllPixelOverlay[MAX_PULS_ORDER] = {NULL}; 157 164 TProfile* hAllPixelProfile[MAX_PULS_ORDER] = {NULL}; 165 TProfile* hAllPixelProfile2[MAX_PULS_ORDER] = {NULL}; 158 166 TH1F* hAllPixelMax[MAX_PULS_ORDER] = {NULL}; 159 TH1F* hAllPixelMaxGaus[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}; 160 170 161 171 //Histogram Parameters … … 183 193 void DeletePixelHistos(int ); 184 194 void SaveHistograms( const char*, const char*, TObjArray*, int ); 185 void FillHistograms(vector<Region>*, int, int );195 void FillHistograms(vector<Region>*, int, int, int); 186 196 void PlotPulseShapeOfMaxPropability(unsigned int, int, int); 197 187 198 void DrawPulseHistograms(int, int); 188 199 void FitMaxPropabilityPuls( TH1F* , int ); 189 200 void DrawMaximumHistograms( int, int ); 190 201 void DrawTestHistograms( int); 202 Double_t MedianOfH1( TH1 *h1 ); 203 void PlotMedianEachSliceOfPulse(unsigned int, int, int); 204 191 205 192 206 void UpdateCanvases( int, int); … … 198 212 TString CreateSubDirName(int); //creates a String containing the path to the subdirectory in root file 199 213 TString CreateSubDirName(const char* ); 214 void WritePixelTemplateToCsv( TString, const char*, const char*, int, int); 215 void WriteAllPixelTemplateToCsv( TString, const char*, const char*, int); 200 216 //void SetHistogrammSettings( const char*, int, int , int); 201 217 //void CreatePulseCanvas( TCanvas* , unsigned int, const char* , const char*, const char*, int); … … 212 228 const char* drsfilename = "../data/2011/11/09/20111109_003.drs.fits.gz", 213 229 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 216 233 bool spikeDebug = false, 217 bool debugPixel = true,234 bool debugPixel = false, 218 235 bool fitdata = false, 219 int verbosityLevel = 1,// different verbosity levels can be implemented here220 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 histograms236 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 225 242 float GainMean = 8, 226 float BSLMean = -1, //4 Histogramms will be drawn, decide how far pulsehights differ from eachother243 float BSLMean = -1, //4 Histogramms will be drawn, decide how far pulsehights differ from eachother 227 244 int avg1 = 8, 228 245 int avg2 = 8, … … 360 377 for ( int pixel = firstpixel; pixel < firstpixel + npixel; pixel++ ) 361 378 { 362 if (verbosityLevel > 0)379 if (verbosityLevel >= 0) 363 380 { 364 381 cout << "------------------------------------------------" << endl … … 440 457 removeRegionWithMaxOnEdge( *pZXings, 2); 441 458 removeRegionOnFallingEdge( *pZXings, 100); 442 findTimeOfHalfMaxLeft(*pZXings, V slide, gBSLMean, 5, 10, verbosityLevel );459 findTimeOfHalfMaxLeft(*pZXings, Vcorr, gBSLMean, 5, 10, verbosityLevel ); 443 460 if (verbosityLevel > 2) cout << "...done" << endl; 444 461 … … 446 463 // Fill Overlay Histos 447 464 //----------------------------------------------------------------------------- 448 FillHistograms(pZXings, AmplWindowWidth, verbosityLevel );465 FillHistograms(pZXings, AmplWindowWidth, ev, verbosityLevel ); 449 466 450 467 //----------------------------------------------------------------------------- … … 557 574 //------------------------------------- 558 575 576 577 578 559 579 //------------------------------------- 560 580 // Histogramms of Maximas in Overlay Spectra … … 567 587 ); 568 588 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 ); 569 602 570 603 if (verbosityLevel > 2) cout << "...done" << endl; … … 599 632 { 600 633 //Process gui events asynchronously during input 634 cout << endl; 601 635 TTimer timer("gSystem->ProcessEvents();", 50, kFALSE); 602 636 timer.TurnOn(); … … 647 681 } 648 682 683 WriteAllPixelTemplateToCsv( 684 OutPutPath, 685 "PulseTemplate_PointSet", 686 "Maximum", 687 verbosityLevel 688 ); 689 649 690 DeletePixelCanvases( verbosityLevel ); 650 691 return( 0 ); … … 662 703 663 704 664 void BookPixelHistos( int verbosityLevel) 705 void 706 BookPixelHistos( int verbosityLevel) 665 707 { 666 708 if (verbosityLevel > 2) cout << endl << "...book pixel histograms" << endl; … … 686 728 512, 687 729 -55.5, 688 300.5730 200.5 689 731 ); 690 732 hPixelOverlay[order]->SetAxisRange( … … 722 764 hList->Add( hPixelMax[order] ); 723 765 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 //------------------------------------------------------------------------ 725 819 // SetHistogramAttributes( 726 820 // "PeakMaxGaus", … … 740 834 512, 741 835 -55.5, 742 300.5836 200.5 743 837 ); 744 838 … … 771 865 300.5, //yup 772 866 "s"); //option 773 hPixel Max[order]->SetAxisRange(867 hPixelProfile[order]->SetAxisRange( 774 868 gBSLMean - 5, 775 869 (gGainMean*(order+1)) + 10, … … 779 873 //hPixelProfile->SetBit(TH2F::kCanRebin); 780 874 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 781 909 } 782 910 if (verbosityLevel > 2) cout << "...done" << endl; … … 795 923 if (verbosityLevel > 3) cout << endl << "...deleting hPixelMax" << order; 796 924 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]; 797 929 if (verbosityLevel > 3) cout << endl << "...deleting hPixelEdgeOverlay" << order; 798 930 delete hPixelEdgeOverlay[order]; 799 931 if (verbosityLevel > 3) cout << endl << "...deleting hPixelProfile" << order; 800 932 delete hPixelProfile[order]; 933 if (verbosityLevel > 3) cout << endl << "...deleting hPixelProfile2" << order; 934 delete hPixelProfile2[order]; 801 935 } 802 936 if (verbosityLevel > 3) cout << endl << "...deleting hList"; … … 814 948 "hTesthisto", 815 949 "Deviation of rising edge and maximum", 816 200,817 - 300,818 300950 600, 951 -10.1, 952 10.1 819 953 ); 820 954 hTesthisto->GetXaxis()->SetTitle( "Timeslices [a.u.]" ); 821 955 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 ); 822 970 823 971 if (verbosityLevel > 2) cout << "...done" << endl; … … 910 1058 hAllPixelList->Add( hAllPixelMax[order] ); 911 1059 912 1060 //------------------------------------------------------------------------ 913 1061 // 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", 915 1115 // order, 916 1116 // gMaxGausAttrib, 917 1117 // MaxAmplOfFirstPulse 918 1118 // ); 919 histoname = "hAllPixel MaxGaus";1119 histoname = "hAllPixelEdgeOverlay"; 920 1120 histoname += order; 921 1121 if (verbosityLevel > 3) cout << "...booking " << histoname << endl; 922 hAllPixel MaxGaus[order] = new TH1F(1122 hAllPixelEdgeOverlay[order] = new TH2F( 923 1123 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 ", 925 1125 gPixelOverlayXaxisLeft + gPixelOverlayXaxisRight , 926 1126 (-1*gPixelOverlayXaxisLeft)-0.5, 927 gPixelOverlayXaxisRight-0.5 1127 gPixelOverlayXaxisRight-0.5 , 1128 512, 1129 -55.5, 1130 200.5 928 1131 ); 929 hAllPixelMaxGaus[order]->SetAxisRange( 1132 1133 hAllPixelEdgeOverlay[order]->SetAxisRange( 930 1134 gBSLMean - 5, 931 1135 (gGainMean*(order+1)) + 10, 932 1136 "Y"); 933 hAllPixel MaxGaus[order]->GetXaxis()->SetTitle( "Timeslices [a.u.]" );934 hAllPixel MaxGaus[order]->GetYaxis()->SetTitle( "Amplitude [mV]" );935 hAllPixelList->Add( hAllPixel MaxGaus[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 //------------------------------------------------------------------------ 938 1142 // SetHistogramAttributes( 939 1143 // "AllPixelPeakProfile", … … 964 1168 hAllPixelList->Add( hAllPixelProfile[order] ); 965 1169 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 966 1200 967 1201 … … 1077 1311 vector<Region>* pZXings, 1078 1312 int AmplWindowWidth, 1313 int eventNumber, 1079 1314 int verbosityLevel 1080 1315 ) … … 1103 1338 1104 1339 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 ) ; 1108 1347 1109 1348 if (verbosityLevel > 2) cout << endl << "\t...choosing Histogram" ; … … 1132 1371 if (verbosityLevel > 2) cout << "...filling" ; 1133 1372 for ( int pos = Left; pos < Right; pos++){ 1373 // if (); 1134 1374 hPixelOverlay[order_of_pulse]->Fill( pos - (reg->maxPos), Ameas[pos]) ; 1135 1375 hPixelProfile[order_of_pulse]->Fill( pos - (reg->maxPos), Ameas[pos]) ; … … 1140 1380 // cout << endl << "###filling edge histo###" << endl; 1141 1381 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]) ; 1142 1385 // cout << endl << "######" << endl; 1143 1386 } … … 1161 1404 cgpTestHistos->cd(1); 1162 1405 hTesthisto->Draw(); 1406 cgpTestHistos->cd(2); 1407 hTesthisto2->Draw("BOX"); 1163 1408 } 1164 1409 // end of DrawTestHistograms … … 1175 1420 cgpPixelPulses[pulse_order]->cd(1); 1176 1421 hPixelOverlay[pulse_order]->Draw("COLZ"); 1422 cgpPixelPulses[pulse_order]->cd(3); 1423 hPixelProfile[pulse_order]->Draw("COLZ"); 1424 hPixelProfile2[pulse_order]->Draw("SAME"); 1425 1177 1426 cgpPixelPulses[pulse_order]->cd(2); 1178 hPixelProfile[pulse_order]->Draw("COLZ");1179 cgpPixelPulses[pulse_order]->cd(3);1180 1427 hPixelEdgeOverlay[pulse_order]->Draw("COLZ"); 1181 1428 // cgpPixelPulses[pulse_order]->cd(4); … … 1185 1432 hAllPixelOverlay[pulse_order]->Draw("COLZ"); 1186 1433 cgpAllPixelPulses[pulse_order]->cd(2); 1434 hAllPixelEdgeOverlay[pulse_order]->Draw("COLZ"); 1435 1436 cgpAllPixelPulses[pulse_order]->cd(3); 1187 1437 hAllPixelProfile[pulse_order]->Draw("COLZ"); 1438 hAllPixelProfile2[pulse_order]->Draw("SAME"); 1188 1439 1189 1440 } … … 1202 1453 { 1203 1454 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"); 1205 1458 1206 1459 cgpAllPixelPulses[pulse_order]->cd(4); 1207 1460 hAllPixelMax[pulse_order]->Draw(); 1461 hAllPixelMedian[pulse_order]->Draw("SAME"); 1462 hAllPixelMean[pulse_order]->Draw("SAME"); 1208 1463 // cgpAllPixelPulses[pulse_order]->cd(3); 1209 1464 // hAllPixelMaxGaus[pulse_order]->Draw("COLZ"); … … 1236 1491 1237 1492 void 1493 PlotMedianEachSliceOfPulse( 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 1535 Double_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 1555 void 1238 1556 PlotPulseShapeOfMaxPropability( 1239 1557 unsigned int max_pulse_order, … … 1247 1565 for (unsigned int pulse_order = 0 ; pulse_order < max_pulse_order ; pulse_order ++) 1248 1566 { 1249 if (verbosityLevel > 2) cout << "\t...calculati ngof "1567 if (verbosityLevel > 2) cout << "\t...calculation of " 1250 1568 << "pulse order " << pulse_order; 1251 1569 // 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++) 1258 1575 { 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; 1262 1577 //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; 1267 1587 } 1268 1588 … … 1282 1602 void 1283 1603 FitMaxPropabilityPuls( 1284 TH1F* hMaximumTemp,1285 int verbosityLevel1604 TH1F* hMaximumTemp, 1605 int verbosityLevel 1286 1606 ) 1287 1607 { … … 1415 1735 // histoAttrib[pulse_order].yTitle += "Amplitude [mV]"; 1416 1736 //} 1737 1738 void 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 1785 void 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.