Rev 172 | Go to most recent revision | Show entire file | Regard whitespace | Details | Blame | Last modification | View Log | RSS feed
Rev 172 | Rev 173 | ||
---|---|---|---|
Line 57... | Line 57... | ||
57 | resol2d->widgetNE[1]->SetNumber(40); |
57 | resol2d->widgetNE[1]->SetNumber(40); |
58 | } |
58 | } |
59 | else if(analTab->GetCurrent() == 1) // Relative PDE |
59 | else if(analTab->GetCurrent() == 1) // Relative PDE |
60 | { |
60 | { |
61 | relPde->widgetChBox[0]->SetState(kButtonDown); |
61 | relPde->widgetChBox[0]->SetState(kButtonDown); |
- | 62 | relPde->widgetChBox[1]->SetState(kButtonUp); |
|
62 | midPeak->widgetChBox[0]->SetState(kButtonUp); |
63 | midPeak->widgetChBox[0]->SetState(kButtonUp); |
63 | zeroAngle->widgetNE[0]->SetNumber(0.00); |
64 | zeroAngle->widgetNE[0]->SetNumber(0.00); |
64 | } |
65 | } |
65 | else if(analTab->GetCurrent() == 2) // Breakdown voltage |
66 | else if(analTab->GetCurrent() == 2) // Breakdown voltage |
66 | { |
67 | { |
Line 93... | Line 94... | ||
93 | analtype = 2; |
94 | analtype = 2; |
94 | else if( analtab == 2 ) |
95 | else if( analtab == 2 ) |
95 | analtype = 3; |
96 | analtype = 3; |
96 | else if( analtab == 3 ) |
97 | else if( analtab == 3 ) |
97 | analtype = 4; |
98 | analtype = 4; |
- | 99 | ||
- | 100 | // Save analysis settings any time we run a new analysis |
|
- | 101 | SaveAnalSettings(); |
|
98 | 102 | ||
99 | // Only integrate spectrum or make relative PDE |
103 | // Only integrate spectrum or make relative PDE |
100 | if(type == 0) |
104 | if(type == 0) |
101 | { |
105 | { |
102 | files = new TList(); |
106 | files = new TList(); |
Line 368... | Line 372... | ||
368 | gCanvas = tempAnalysisCanvas->GetCanvas(); |
372 | gCanvas = tempAnalysisCanvas->GetCanvas(); |
369 | 373 | ||
370 | double range[4]; |
374 | double range[4]; |
371 | TGraph2D *gScan2D; |
375 | TGraph2D *gScan2D; |
372 | gScan2D = new TGraph2D(); |
376 | gScan2D = new TGraph2D(); |
373 | gScan2D->SetName("edgescan"); |
- | |
374 | range[0] = TMath::MinElement(nrfiles, surfxy); |
377 | range[0] = TMath::MinElement(nrfiles, surfxy); |
375 | range[1] = TMath::MaxElement(nrfiles, surfxy); |
378 | range[1] = TMath::MaxElement(nrfiles, surfxy); |
376 | range[2] = TMath::MinElement(nrfiles, surfz); |
379 | range[2] = TMath::MinElement(nrfiles, surfz); |
377 | range[3] = TMath::MaxElement(nrfiles, surfz); |
380 | range[3] = TMath::MaxElement(nrfiles, surfz); |
378 | 381 | ||
Line 531... | Line 534... | ||
531 | remove_from_last((char*)files->At(filenr-1)->GetTitle(), '_', ctemp); |
534 | remove_from_last((char*)files->At(filenr-1)->GetTitle(), '_', ctemp); |
532 | sprintf(exportname, "%s_edge.pdf", ctemp); |
535 | sprintf(exportname, "%s_edge.pdf", ctemp); |
533 | 536 | ||
534 | if(edit == 0) |
537 | if(edit == 0) |
535 | gCanvas = new TCanvas("canv1d","canv1d",1200,900); |
538 | gCanvas = new TCanvas("canv1d","canv1d",1200,900); |
536 | else |
539 | else |
537 | gCanvas = tempAnalysisCanvas->GetCanvas(); |
540 | gCanvas = tempAnalysisCanvas->GetCanvas(); |
538 | 541 | ||
539 | // Fit the PDF with a gaussian |
542 | // Fit the PDF with a gaussian |
540 | gScan[0]->Fit("gaus","Q"); |
543 | gScan[0]->Fit("gaus","Q"); |
541 | gScan[0]->GetFunction("gaus")->SetNpx(400); |
544 | gScan[0]->GetFunction("gaus")->SetNpx(400); |
542 | 545 | ||
Line 548... | Line 551... | ||
548 | 551 | ||
549 | gCanvas->Modified(); |
552 | gCanvas->Modified(); |
550 | gCanvas->Update(); |
553 | gCanvas->Update(); |
551 | 554 | ||
552 | TPaveStats *stats = (TPaveStats*)gScan[0]->FindObject("stats"); |
555 | TPaveStats *stats = (TPaveStats*)gScan[0]->FindObject("stats"); |
553 | if(!cleanPlots) |
556 | if(!cleanPlots) |
554 | { |
557 | { |
555 | stats->SetX1NDC(0.86); stats->SetX2NDC(1.0); |
558 | stats->SetX1NDC(0.86); stats->SetX2NDC(1.0); |
556 | stats->SetY1NDC(0.87); stats->SetY2NDC(1.0); |
559 | stats->SetY1NDC(0.87); stats->SetY2NDC(1.0); |
557 | } |
560 | } |
558 | else |
561 | else |
559 | { |
562 | { |
Line 587... | Line 590... | ||
587 | { |
590 | { |
588 | if(posUnitsPlot->widgetCB->GetSelected() == 0) |
591 | if(posUnitsPlot->widgetCB->GetSelected() == 0) |
589 | gScan[1]->SetTitle("SiPM edge detection;X [table units];Normalized ADC integral"); |
592 | gScan[1]->SetTitle("SiPM edge detection;X [table units];Normalized ADC integral"); |
590 | else if(posUnitsPlot->widgetCB->GetSelected() == 1) |
593 | else if(posUnitsPlot->widgetCB->GetSelected() == 1) |
591 | gScan[1]->SetTitle("SiPM edge detection;X [#mum];Normalized ADC integral"); |
594 | gScan[1]->SetTitle("SiPM edge detection;X [#mum];Normalized ADC integral"); |
592 | } |
595 | } |
593 | else if(axis == 2) |
596 | else if(axis == 2) |
594 | { |
597 | { |
595 | if(posUnitsPlot->widgetCB->GetSelected() == 0) |
598 | if(posUnitsPlot->widgetCB->GetSelected() == 0) |
596 | gScan[1]->SetTitle("SiPM edge detection;Y [table units];Normalized ADC integral"); |
599 | gScan[1]->SetTitle("SiPM edge detection;Y [table units];Normalized ADC integral"); |
597 | else if(posUnitsPlot->widgetCB->GetSelected() == 1) |
600 | else if(posUnitsPlot->widgetCB->GetSelected() == 1) |
598 | gScan[1]->SetTitle("SiPM edge detection;Y [#mum];Normalized ADC integral"); |
601 | gScan[1]->SetTitle("SiPM edge detection;Y [#mum];Normalized ADC integral"); |
599 | } |
602 | } |
600 | } |
603 | } |
601 | else |
604 | else |
602 | { |
605 | { |
603 | if(axis == 1) |
606 | if(axis == 1) |
604 | { |
607 | { |
605 | if(posUnitsPlot->widgetCB->GetSelected() == 0) |
608 | if(posUnitsPlot->widgetCB->GetSelected() == 0) |
606 | gScan[1]->SetTitle(";X [table units];Normalized ADC integral"); |
609 | gScan[1]->SetTitle(";X [table units];Normalized ADC integral"); |
607 | else if(posUnitsPlot->widgetCB->GetSelected() == 1) |
610 | else if(posUnitsPlot->widgetCB->GetSelected() == 1) |
608 | gScan[1]->SetTitle(";X [#mum];Normalized ADC integral"); |
611 | gScan[1]->SetTitle(";X [#mum];Normalized ADC integral"); |
609 | } |
612 | } |
610 | else if(axis == 2) |
613 | else if(axis == 2) |
611 | { |
614 | { |
612 | if(posUnitsPlot->widgetCB->GetSelected() == 0) |
615 | if(posUnitsPlot->widgetCB->GetSelected() == 0) |
613 | gScan[1]->SetTitle(";Y [table units];Normalized ADC integral"); |
616 | gScan[1]->SetTitle(";Y [table units];Normalized ADC integral"); |
614 | else if(posUnitsPlot->widgetCB->GetSelected() == 1) |
617 | else if(posUnitsPlot->widgetCB->GetSelected() == 1) |
615 | gScan[1]->SetTitle(";Y [#mum];Normalized ADC integral"); |
618 | gScan[1]->SetTitle(";Y [#mum];Normalized ADC integral"); |
616 | } |
619 | } |
617 | } |
620 | } |
618 | gScan[1]->SetLineColor(kBlue); |
621 | gScan[1]->SetLineColor(kBlue); |
619 | gScan[0]->SetLineWidth(2); |
622 | gScan[0]->SetLineWidth(2); |
620 | gScan[1]->SetLineWidth(2); |
623 | gScan[1]->SetLineWidth(2); |
621 | 624 | ||
622 | gCanvas->Modified(); |
625 | gCanvas->Modified(); |
623 | gCanvas->Update(); |
626 | gCanvas->Update(); |
624 | 627 | ||
625 | gCanvas->SaveAs(exportname); |
628 | gCanvas->SaveAs(exportname); |
626 | 629 | ||
627 | // If 2D edge plot, delete the 1D edge plots as we go |
630 | // If 2D edge plot, delete the 1D edge plots as we go |
628 | if(edge2d) |
631 | if(edge2d) |
629 | { |
632 | { |
630 | delete gScan[0]; |
633 | delete gScan[0]; |
631 | delete gScan[1]; |
634 | delete gScan[1]; |
632 | if(edit == 0) |
635 | if(edit == 0) |
633 | delete gCanvas; |
636 | delete gCanvas; |
634 | } |
637 | } |
635 | else |
638 | else |
636 | { |
639 | { |
637 | if(edit == 0) |
640 | if(edit == 0) |
638 | { |
641 | { |
639 | delete gScan[0]; |
642 | delete gScan[0]; |
640 | delete gScan[1]; |
643 | delete gScan[1]; |
641 | delete gCanvas; |
644 | delete gCanvas; |
642 | } |
645 | } |
643 | } |
646 | } |
644 | } |
647 | } |
645 | 648 | ||
Line 677... | Line 680... | ||
677 | double meanparam, paramsigma; |
680 | double meanparam, paramsigma; |
678 | int sortindex[20]; |
681 | int sortindex[20]; |
679 | int adcpedestal[2]; |
682 | int adcpedestal[2]; |
680 | int zeromu = 0; |
683 | int zeromu = 0; |
681 | int darkhist = -1; |
684 | int darkhist = -1; |
- | 685 | int nopeaks = -1; |
|
682 | 686 | ||
683 | double pointest[12]; |
687 | double pointest[12]; |
684 | bool exclude = false; |
688 | bool exclude = false; |
685 | 689 | ||
686 | // Zero the parameter values |
690 | // Zero the parameter values |
687 | for(int i = 0; i < 20; i++) {meansel[i] = 0; sigmasel[i] = 0; } |
691 | for(int i = 0; i < 20; i++) {meansel[i] = 0; sigmasel[i] = 0; } |
688 | 692 | ||
689 | float progVal = 0; |
693 | float progVal = 0; |
690 | analysisProgress->widgetPB->SetPosition(progVal); |
694 | analysisProgress->widgetPB->SetPosition(progVal); |
691 | gVirtualX->Update(1); |
695 | gVirtualX->Update(1); |
692 | 696 | ||
693 | // |
697 | // Check if the checkbox for no peaks is selected - TODO: Still need the situation when we do not have a peaked ADC spectrum |
694 | if( |
698 | if(relPde->widgetChBox[1]->IsDown()) |
695 | { |
699 | { |
696 |
|
700 | printf("PhotonMu(): ADC spectrum has no peak structure.\n"); |
697 |
|
701 | nopeaks = 1; |
698 | if( (nrfiles == 1) || (!multiSelect->widgetChBox[0]->IsDown()) ) |
702 | |
699 |
|
703 | // Error if there is no darkhist |
700 |
|
704 | if(strcmp("", darkRun->widgetTE->GetText()) == 0) |
701 | 705 | { |
|
702 |
|
706 | printf("PhotonMu(): Error! The no peak structure option needs a dark histogram.\n"); |
703 |
|
707 | delete[] integralCount; |
704 |
|
708 | delete[] integralPedestal; |
705 |
|
709 | delete[] angle; |
706 | 710 | delete[] pdeval; |
|
707 |
|
711 | delete[] muval; |
708 |
|
712 | return; |
709 |
|
713 | } |
710 | double par[300]; |
714 | |
711 |
|
715 | } |
712 | // Find spectrum background |
716 | |
713 |
|
717 | // Start if we select at least one file |
714 |
|
718 | if(nrfiles > 0) |
715 |
|
719 | { |
716 |
|
720 | // Find the pedestal peak for the dark histogram, and use it for all if there are no peaks |
717 |
|
721 | if(nopeaks != -1) |
718 |
|
722 | { |
719 |
|
723 | // Replot the spectrum on analysisCanvas and do not close the input file |
720 |
|
724 | DisplayHistogram( (char*)(darkRun->widgetTE->GetText()), 0, 1); |
721 |
|
725 | analysisCanvas->GetCanvas()->Modified(); |
722 | 726 | analysisCanvas->GetCanvas()->Update(); |
|
723 |
|
727 | |
724 |
|
728 | // Get the spectrum |
725 |
|
729 | histtemp = (TH1F*)analysisCanvas->GetCanvas()->GetPrimitive(histname); |
726 |
|
730 | npeaks = 15; |
727 |
|
731 | double par[300]; |
728 |
|
732 | spec = new TSpectrum(npeaks); |
729 |
|
733 | // Find spectrum background |
730 |
|
734 | histback = spec->Background(histtemp, (int)fitInter->widgetNE[0]->GetNumber(), "same"); |
731 |
|
735 | // Clone histogram and subtract background from it if we select that option |
732 |
|
736 | h2 = (TH1F*)histtemp->Clone("h2"); |
733 |
|
737 | if(fitChecks->widgetChBox[0]->IsDown()) |
734 | 738 | h2->Add(histback, -1); |
|
735 |
|
739 | // Search for the peaks |
736 |
|
740 | int found = spec->Search(h2, fitSigma->widgetNE[0]->GetNumber(), "goff", fitTresh->widgetNE[0]->GetNumber() ); |
737 |
|
741 | printf("PhotonMu(): Found %d candidates to fit.\n",found); |
738 |
|
742 | npeaks = found; |
739 |
|
743 | |
740 |
|
744 | // Set initial peak parameters |
741 |
|
745 | xpeaks = spec->GetPositionX(); |
742 |
|
746 | for(j = 0; j < found; j++) |
743 |
|
747 | { |
744 |
|
748 | float xp = xpeaks[j]; |
745 | 749 | int bin = h2->GetXaxis()->FindBin(xp); |
|
746 |
|
750 | float yp = h2->GetBinContent(bin); |
747 |
|
751 | par[3*j] = yp; |
748 |
|
752 | par[3*j+1] = xp; |
749 |
|
753 | par[3*j+2] = (double)fitSigma->widgetNE[0]->GetNumber(); |
750 |
|
754 | } |
751 |
|
755 | |
752 |
|
756 | // Fit the histogram |
753 |
|
757 | fit = new TF1("fit", FindPeaks, adcRange->widgetNE[0]->GetNumber(), adcRange->widgetNE[1]->GetNumber(), 3*npeaks); |
754 |
|
758 | TVirtualFitter::Fitter(histtemp, 3*npeaks); |
755 |
|
759 | fit->SetParameters(par); |
- | 760 | fit->SetNpx(300); |
|
- | 761 | h2->Fit("fit","Q"); |
|
- | 762 | // Get the fitted parameters |
|
- | 763 | fittingfunc = h2->GetFunction("fit"); |
|
- | 764 | fparam = fittingfunc->GetParameters(); |
|
- | 765 | fparamerr = fittingfunc->GetParErrors(); |
|
- | 766 | ||
- | 767 | // Gather the parameters (mean peak value for now) |
|
- | 768 | int j = 1; |
|
- | 769 | int nrfit = 0; |
|
- | 770 | while(1) |
|
- | 771 | { |
|
- | 772 | if( (fparam[j] < 1.E-30) || (nrfit > 8) ) |
|
- | 773 | break; |
|
- | 774 | else |
|
- | 775 | { |
|
- | 776 | // Check if pedestal is above the lower limit and sigma is smaller than the mean |
|
- | 777 | if( (fparam[j] > pedesLow->widgetNE[0]->GetNumber()) && ((double)fparamerr[j]/fparam[j] < accError->widgetNE[0]->GetNumber()) ) |
|
- | 778 | { |
|
- | 779 | // With the additional ADC offset, we can shift the mean values slightly, so they are not close to the X.5, but to the X.0 values |
|
- | 780 | meansel[nrfit] = fparam[j]+(adcOffset->widgetNE[0]->GetNumber()); |
|
- | 781 | sigmasel[nrfit] = fparam[j+1]; |
|
- | 782 | nrfit++; |
|
- | 783 | } |
|
- | 784 | } |
|
- | 785 | ||
- | 786 | j+=3; |
|
- | 787 | } |
|
- | 788 | TMath::Sort(nrfit, meansel, sortindex, kFALSE); |
|
- | 789 | ||
- | 790 | meanparam = meansel[sortindex[0]]; |
|
- | 791 | paramsigma = sigmasel[sortindex[0]]; |
|
- | 792 | ||
- | 793 | for(j = 0; j < nrfit; j++) |
|
- | 794 | if(DBGSIG) |
|
- | 795 | printf("PhotonMu(): %d: peak mean = %lf\n", j, meansel[sortindex[j]]); |
|
- | 796 | ||
- | 797 | j = 0; |
|
- | 798 | adcpedestal[0] = 0; |
|
- | 799 | adcpedestal[1] = -1; |
|
- | 800 | ||
- | 801 | while(1) |
|
- | 802 | { |
|
- | 803 | int bin = histtemp->GetXaxis()->FindBin((int)(j+meanparam+paramsigma)); |
|
- | 804 | int yp = histtemp->GetBinContent(bin); |
|
- | 805 | ||
- | 806 | // Check where we get to first minimum after pedestal peak or where we get to the half maximum of the pedestal peak (in case there is only a pedestal peak) |
|
- | 807 | if(adcpedestal[1] == -1) |
|
- | 808 | { |
|
- | 809 | adcpedestal[0] = j+meanparam+paramsigma; |
|
- | 810 | adcpedestal[1] = yp; |
|
- | 811 | } |
|
- | 812 | else |
|
- | 813 | { |
|
- | 814 | if( (npeaks > 1) && (adcpedestal[1] >= yp) ) |
|
- | 815 | { |
|
- | 816 | adcpedestal[0] = j+meanparam+paramsigma; |
|
- | 817 | adcpedestal[1] = yp; |
|
- | 818 | } |
|
- | 819 | else if( (npeaks == 1) && (adcpedestal[0] < meanparam+5*paramsigma) ) // TODO -> Determining the pedestal when only one peak |
|
- | 820 | { |
|
- | 821 | adcpedestal[0] = j+meanparam+paramsigma; |
|
- | 822 | adcpedestal[1] = yp; |
|
- | 823 | } |
|
- | 824 | else |
|
- | 825 | break; |
|
- | 826 | } |
|
- | 827 | ||
- | 828 | j++; |
|
- | 829 | if(j > 50) break; |
|
- | 830 | } |
|
- | 831 | ||
- | 832 | if(midPeak->widgetChBox[0]->IsDown()) |
|
- | 833 | { |
|
- | 834 | if( (meanparam - (int)meanparam >= 0.) && (meanparam - (int)meanparam < 0.5) ) |
|
- | 835 | m = TMath::Floor(meanparam); |
|
- | 836 | else if( (meanparam - (int)meanparam >= 0.5) && (meanparam - (int)meanparam < 1.) ) |
|
- | 837 | m = TMath::Ceil(meanparam); |
|
- | 838 | int bin = histtemp->GetXaxis()->FindBin(m); |
|
- | 839 | adcpedestal[0] = m; |
|
- | 840 | printf("midpeak x = %d, ", adcpedestal[0]); |
|
- | 841 | adcpedestal[1] = histtemp->GetBinContent(bin); |
|
- | 842 | } |
|
- | 843 | ||
- | 844 | // Option to show the fit |
|
- | 845 | fittingfunc->Draw("L SAME"); |
|
- | 846 | analysisCanvas->GetCanvas()->Modified(); |
|
- | 847 | analysisCanvas->GetCanvas()->Update(); |
|
- | 848 | ||
- | 849 | printf("Pedestal ends = %d and nr. of counts = %d\n", adcpedestal[0], adcpedestal[1]); |
|
- | 850 | ||
- | 851 | // Delete the opened histogram and spectrum |
|
- | 852 | delete spec; |
|
- | 853 | delete inroot; |
|
- | 854 | ||
- | 855 | // return; |
|
- | 856 | } |
|
- | 857 | ||
- | 858 | printf("PhotonMu(): Continuing with the rest of the spectra.\n"); |
|
- | 859 | ||
- | 860 | // Check all histograms for pedestal peak values |
|
- | 861 | for(int i = 0; i < (int)nrfiles; i++) |
|
- | 862 | { |
|
- | 863 | if( (nrfiles == 1) || (!multiSelect->widgetChBox[0]->IsDown()) ) |
|
- | 864 | { |
|
- | 865 | printf("PhotonMu(): Only one file selected. Not running analysis, just showing the fit.\n"); |
|
- | 866 | ||
- | 867 | // Replot the spectrum on analysisCanvas and do not close the input file |
|
- | 868 | DisplayHistogram( (char*)(files->At(i)->GetTitle()), 0, 1); |
|
- | 869 | analysisCanvas->GetCanvas()->Modified(); |
|
- | 870 | analysisCanvas->GetCanvas()->Update(); |
|
- | 871 | ||
- | 872 | // Get the spectrum |
|
- | 873 | histtemp = (TH1F*)analysisCanvas->GetCanvas()->GetPrimitive(histname); |
|
- | 874 | npeaks = 15; |
|
- | 875 | double par[300]; |
|
- | 876 | spec = new TSpectrum(npeaks); |
|
- | 877 | // Find spectrum background |
|
- | 878 | histback = spec->Background(histtemp, (int)fitInter->widgetNE[0]->GetNumber(), "same"); |
|
- | 879 | // Clone histogram and subtract background from it if we select that option |
|
- | 880 | h2 = (TH1F*)histtemp->Clone("h2"); |
|
- | 881 | if(fitChecks->widgetChBox[0]->IsDown()) |
|
- | 882 | h2->Add(histback, -1); |
|
- | 883 | // Search for the peaks |
|
- | 884 | int found = spec->Search(h2, fitSigma->widgetNE[0]->GetNumber(), "goff", fitTresh->widgetNE[0]->GetNumber() ); |
|
- | 885 | printf("PhotonMu(): Found %d candidates to fit.\n",found); |
|
- | 886 | npeaks = found; |
|
- | 887 | ||
- | 888 | // Set initial peak parameters |
|
- | 889 | xpeaks = spec->GetPositionX(); |
|
- | 890 | for(j = 0; j < found; j++) |
|
- | 891 | { |
|
- | 892 | float xp = xpeaks[j]; |
|
- | 893 | int bin = h2->GetXaxis()->FindBin(xp); |
|
- | 894 | float yp = h2->GetBinContent(bin); |
|
- | 895 | par[3*j] = yp; |
|
- | 896 | par[3*j+1] = xp; |
|
- | 897 | par[3*j+2] = (double)fitSigma->widgetNE[0]->GetNumber(); |
|
- | 898 | } |
|
- | 899 | ||
- | 900 | // Fit the histogram |
|
- | 901 | fit = new TF1("fit", FindPeaks, adcRange->widgetNE[0]->GetNumber(), adcRange->widgetNE[1]->GetNumber(), 3*npeaks); |
|
- | 902 | TVirtualFitter::Fitter(histtemp, 3*npeaks); |
|
- | 903 | fit->SetParameters(par); |
|
- | 904 | fit->SetNpx(300); |
|
- | 905 | h2->Fit("fit","Q"); |
|
- | 906 | // Get the fitted parameters |
|
- | 907 | fittingfunc = h2->GetFunction("fit"); |
|
- | 908 | fparam = fittingfunc->GetParameters(); |
|
- | 909 | fparamerr = fittingfunc->GetParErrors(); |
|
- | 910 | ||
- | 911 | // Gather the parameters (mean peak value for now) |
|
- | 912 | int j = 1; |
|
- | 913 | int nrfit = 0; |
|
- | 914 | while(1) |
|
- | 915 | { |
|
- | 916 | if( (fparam[j] < 1.E-30) || (nrfit > 8) ) |
|
- | 917 | break; |
|
- | 918 | else |
|
- | 919 | { |
|
- | 920 | // Check if pedestal is above the lower limit and sigma is smaller than the mean |
|
756 | if( (fparam[j] > pedesLow->widgetNE[0]->GetNumber()) && ((double)fparamerr[j]/fparam[j] < accError->widgetNE[0]->GetNumber()) ) |
921 | if( (fparam[j] > pedesLow->widgetNE[0]->GetNumber()) && ((double)fparamerr[j]/fparam[j] < accError->widgetNE[0]->GetNumber()) ) |
757 | { |
922 | { |
758 | // With the additional ADC offset, we can shift the mean values slightly, so they are not close to the X.5, but to the X.0 values |
923 | // With the additional ADC offset, we can shift the mean values slightly, so they are not close to the X.5, but to the X.0 values |
759 | meansel[nrfit] = fparam[j]+(adcOffset->widgetNE[0]->GetNumber()); |
924 | meansel[nrfit] = fparam[j]+(adcOffset->widgetNE[0]->GetNumber()); |
760 | sigmasel[nrfit] = fparam[j+1]; |
925 | sigmasel[nrfit] = fparam[j+1]; |
761 | nrfit++; |
926 | nrfit++; |
762 | } |
927 | } |
763 | } |
928 | } |
764 | 929 | ||
765 | j+=3; |
930 | j+=3; |
766 | } |
931 | } |
767 | TMath::Sort(nrfit, meansel, sortindex, kFALSE); |
932 | TMath::Sort(nrfit, meansel, sortindex, kFALSE); |
768 | 933 | ||
769 | fittingfunc->Draw("SAME"); |
934 | fittingfunc->Draw("SAME"); |
770 | analysisCanvas->GetCanvas()->Modified(); |
935 | analysisCanvas->GetCanvas()->Modified(); |
771 | analysisCanvas->GetCanvas()->Update(); |
936 | analysisCanvas->GetCanvas()->Update(); |
772 | 937 | ||
773 | meanparam = meansel[sortindex[0]]; |
938 | meanparam = meansel[sortindex[0]]; |
774 | paramsigma = sigmasel[sortindex[0]]; |
939 | paramsigma = sigmasel[sortindex[0]]; |
775 | 940 | ||
776 | for(j = 0; j < nrfit; j++) |
941 | for(j = 0; j < nrfit; j++) |
777 | printf("PhotonMu(): %d: peak mean = %lf\n", j, meansel[sortindex[j]]); |
942 | printf("PhotonMu(): %d: peak mean = %lf\n", j, meansel[sortindex[j]]); |
778 | 943 | ||
779 | 944 | ||
780 | 945 | ||
781 | return; |
946 | return; |
782 | } |
947 | } |
783 | if(files->At(i)) |
948 | if(files->At(i)) |
784 | { |
949 | { |
785 | if(strcmp(files->At(i)->GetTitle(),darkRun->widgetTE->GetText()) == 0) |
950 | if(strcmp(files->At(i)->GetTitle(),darkRun->widgetTE->GetText()) == 0) |
786 | { |
951 | { |
787 | printf("PhotonMu(): %s is the dark histogram file.\n", files->At(i)->GetTitle()); |
952 | printf("PhotonMu(): %s is the dark histogram file.\n", files->At(i)->GetTitle()); |
788 | darkhist = i; |
953 | darkhist = i; |
789 | } |
954 | } |
790 | 955 | ||
- | 956 | if(nopeaks == -1) |
|
- | 957 | { |
|
791 | // Replot the spectrum on analysisCanvas and do not close the input file |
958 | // Replot the spectrum on analysisCanvas and do not close the input file |
792 | DisplayHistogram( (char*)(files->At(i)->GetTitle()), 0, 1); |
959 | DisplayHistogram( (char*)(files->At(i)->GetTitle()), 0, 1); |
793 | analysisCanvas->GetCanvas()->Modified(); |
960 | analysisCanvas->GetCanvas()->Modified(); |
794 | analysisCanvas->GetCanvas()->Update(); |
961 | analysisCanvas->GetCanvas()->Update(); |
795 | 962 | ||
Line 836... | Line 1003... | ||
836 | int j = 1; |
1003 | int j = 1; |
837 | int nrfit = 0; |
1004 | int nrfit = 0; |
838 | while(1) |
1005 | while(1) |
839 | { |
1006 | { |
840 | if( (fparam[j] < 1.E-30) || (nrfit > 8) ) |
1007 | if( (fparam[j] < 1.E-30) || (nrfit > 8) ) |
841 | break; |
1008 | break; |
842 | else |
1009 | else |
843 | { |
1010 | { |
844 | // Check if pedestal is above the lower limit and sigma is smaller than the mean |
1011 | // Check if pedestal is above the lower limit and sigma is smaller than the mean |
845 | if( (fparam[j] > pedesLow->widgetNE[0]->GetNumber()) && ((double)fparamerr[j]/fparam[j] < accError->widgetNE[0]->GetNumber()) ) |
1012 | if( (fparam[j] > pedesLow->widgetNE[0]->GetNumber()) && ((double)fparamerr[j]/fparam[j] < accError->widgetNE[0]->GetNumber()) ) |
846 | { |
1013 | { |
847 | // With the additional ADC offset, we can shift the mean values slightly, so they are not close to the X.5, but to the X.0 values |
1014 | // With the additional ADC offset, we can shift the mean values slightly, so they are not close to the X.5, but to the X.0 values |
848 | meansel[nrfit] = fparam[j]+(adcOffset->widgetNE[0]->GetNumber()); |
1015 | meansel[nrfit] = fparam[j]+(adcOffset->widgetNE[0]->GetNumber()); |
849 | sigmasel[nrfit] = fparam[j+1]; |
1016 | sigmasel[nrfit] = fparam[j+1]; |
850 | nrfit++; |
1017 | nrfit++; |
851 | } |
1018 | } |
852 | } |
1019 | } |
853 | 1020 | ||
854 | j+=3; |
1021 | j+=3; |
855 | } |
1022 | } |
856 | TMath::Sort(nrfit, meansel, sortindex, kFALSE); |
1023 | TMath::Sort(nrfit, meansel, sortindex, kFALSE); |
857 | 1024 | ||
858 | meanparam = meansel[sortindex[0]]; |
1025 | meanparam = meansel[sortindex[0]]; |
859 | paramsigma = sigmasel[sortindex[0]]; |
1026 | paramsigma = sigmasel[sortindex[0]]; |
860 | 1027 | ||
861 | for(j = 0; j < nrfit; j++) |
1028 | for(j = 0; j < nrfit; j++) |
862 | if(DBGSIG) |
1029 | if(DBGSIG) |
863 | printf("PhotonMu(): %d: peak mean = %lf\n", j, meansel[sortindex[j]]); |
1030 | printf("PhotonMu(): %d: peak mean = %lf\n", j, meansel[sortindex[j]]); |
864 | 1031 | ||
865 | j = 0; |
1032 | j = 0; |
866 | adcpedestal[0] = 0; |
1033 | adcpedestal[0] = 0; |
867 | adcpedestal[1] = -1; |
1034 | adcpedestal[1] = -1; |
868 | 1035 | ||
869 | while(1) |
1036 | while(1) |
Line 895... | Line 1062... | ||
895 | 1062 | ||
896 | j++; |
1063 | j++; |
897 | if(j > 50) break; |
1064 | if(j > 50) break; |
898 | } |
1065 | } |
899 | 1066 | ||
900 | if |
1067 | if(npeaks > 1) |
901 | { |
1068 | { |
902 | int bin = histtemp->GetXaxis()->FindBin((int)(meanparam+meansel[sortindex[1]])/2); |
1069 | int bin = histtemp->GetXaxis()->FindBin((int)(meanparam+meansel[sortindex[1]])/2); |
903 | adcpedestal[0] = (meanparam+meansel[sortindex[1]])/2; |
1070 | adcpedestal[0] = (meanparam+meansel[sortindex[1]])/2; |
904 | printf("PhotonMu(): multipeak x = %d, ", adcpedestal[0]); |
1071 | printf("PhotonMu(): multipeak x = %d, ", adcpedestal[0]); |
905 | adcpedestal[1] = histtemp->GetBinContent(bin); |
1072 | adcpedestal[1] = histtemp->GetBinContent(bin); |
Line 925... | Line 1092... | ||
925 | printf("Pedestal ends = %d and nr. of counts = %d\n", adcpedestal[0], adcpedestal[1]); |
1092 | printf("Pedestal ends = %d and nr. of counts = %d\n", adcpedestal[0], adcpedestal[1]); |
926 | 1093 | ||
927 | // Delete the opened histogram and spectrum |
1094 | // Delete the opened histogram and spectrum |
928 | delete spec; |
1095 | delete spec; |
929 | delete inroot; |
1096 | delete inroot; |
- | 1097 | } |
|
930 | 1098 | ||
931 | // Open the input file and read header, ADC and TDC values |
1099 | // Open the input file and read header, ADC and TDC values |
932 | sprintf(ctemp, "%s", files->At(i)->GetTitle()); |
1100 | sprintf(ctemp, "%s", files->At(i)->GetTitle()); |
933 | inroot = new TFile(ctemp, "READ"); |
1101 | inroot = new TFile(ctemp, "READ"); |
934 | 1102 | ||
Line 992... | Line 1160... | ||
992 | integralPedestal[i] += (double)m2; |
1160 | integralPedestal[i] += (double)m2; |
993 | printf("PhotonMu(): %lf: Pedestal integral (%d evts) = %lf\n", angle[i], k2, integralPedestal[i]); |
1161 | printf("PhotonMu(): %lf: Pedestal integral (%d evts) = %lf\n", angle[i], k2, integralPedestal[i]); |
994 | if( (angle[i] == zeroAngle->widgetNE[0]->GetNumber()) && (darkhist != i) ) |
1162 | if( (angle[i] == zeroAngle->widgetNE[0]->GetNumber()) && (darkhist != i) ) |
995 | zeromu = i; |
1163 | zeromu = i; |
996 | 1164 | ||
997 | // Checking for errors when fitting a histogram |
- | |
998 | if(k2 == 0) |
- | |
999 | { |
- | |
1000 | printf("PhotonMu(): No pedestal entries found. Check the fitting results.\n"); |
- | |
1001 | muval[i] = -1; |
- | |
1002 | } |
- | |
1003 | else |
- | |
1004 | muval[i] = -TMath::Log((double)k2/(double)k); |
1165 | muval[i] = -TMath::Log((double)k2/(double)k); |
1005 | printf("PhotonMu(): %lf: muval = %lf\n", angle[i], muval[i]); |
1166 | printf("PhotonMu(): %lf: muval = %lf\n", angle[i], muval[i]); |
1006 | 1167 | ||
1007 | inroot->Close(); |
1168 | inroot->Close(); |
1008 | delete inroot; |
1169 | delete inroot; |
1009 | } |
1170 | } |
1010 | 1171 | ||
1011 | // Update the progress bar |
1172 | // Update the progress bar |
1012 | progVal = (float)(90.00/nrfiles)*i; |
1173 | progVal = (float)(90.00/nrfiles)*i; |
1013 | analysisProgress->widgetPB->SetPosition(progVal); |
1174 | analysisProgress->widgetPB->SetPosition(progVal); |
1014 | gVirtualX->Update(1); |
1175 | gVirtualX->Update(1); |
1015 | } |
1176 | } |
1016 | 1177 | ||
1017 | printf("PhotonMu(): %d files were selected.\n", nrfiles); |
1178 | printf("PhotonMu(): %d files were selected.\n", nrfiles); |
1018 | 1179 | ||
1019 | printf("PhotonMu(): angle\tmu\trelative PDE\n"); |
1180 | printf("PhotonMu(): angle\tmu\trelative PDE\n"); |
1020 | m = 0; |
1181 | m = 0; |
1021 | 1182 | ||
1022 | // Set the 0 degree muval, reuse meansel[1] |
1183 | // Set the 0 degree muval, reuse meansel[1] |
1023 | meansel[1] = muval[zeromu]; |
1184 | meansel[1] = muval[zeromu]; |
1024 | printf("Zero value (id=%d, angle=%lf) = %lf\n", zeromu, angle[zeromu], meansel[1]); |
1185 | printf("Zero value (id=%d, angle=%lf) = %lf\n", zeromu, angle[zeromu], meansel[1]); |
1025 | 1186 | ||
1026 | // TODO - point estimation still not working correctly! |
1187 | // TODO - point estimation still not working correctly! |
Line 1033... | Line 1194... | ||
1033 | exclude = false; |
1194 | exclude = false; |
1034 | 1195 | ||
1035 | // Get next point values (if zero value -> need to add the dark hist value again) |
1196 | // Get next point values (if zero value -> need to add the dark hist value again) |
1036 | pointest[10] = angle[i]; |
1197 | pointest[10] = angle[i]; |
1037 | pointest[11] = muval[i]; |
1198 | pointest[11] = muval[i]; |
1038 | 1199 | ||
1039 | // Check if next point has larger error than acceptable (if yes, set exclude signal to true), reuse meansel[0] |
1200 | // Check if next point has larger error than acceptable (if yes, set exclude signal to true), reuse meansel[0] |
1040 | meansel[0] = PointEstimate(5, pointest); // PointEstimate only works with very small step size |
1201 | meansel[0] = PointEstimate(5, pointest); // PointEstimate only works with very small step size |
1041 | if(meansel[0] > accError->widgetNE[0]->GetNumber()) |
1202 | if(meansel[0] > accError->widgetNE[0]->GetNumber()) |
1042 | { |
1203 | { |
1043 | printf("PhotonMu(): Point (%lf, %lf) excluded due to error: %lf\n", pointest[10], pointest[11], meansel[0]); |
1204 | printf("PhotonMu(): Point (%lf, %lf) excluded due to error: %lf\n", pointest[10], pointest[11], meansel[0]); |
1044 | exclude = true; |
1205 | exclude = true; |
1045 | } |
1206 | } |
1046 | 1207 | ||
1047 | // Value with 0 angle and dark histogram are always needed, so should not be excluded |
1208 | // Value with 0 angle and dark histogram are always needed, so should not be excluded |
1048 | if(i == darkhist) |
1209 | if(i == darkhist) |
1049 | exclude = false; |
1210 | exclude = false; |
1050 | - | ||
1051 | // Wrong fit |
- | |
1052 | if(muval[i] == -1) |
- | |
1053 | exclude = true; |
- | |
1054 | 1211 | ||
1055 | // If nothing excluded, pass the points in pointest variable like in a FIFO |
1212 | // If nothing excluded, pass the points in pointest variable like in a FIFO |
1056 | if(!exclude) |
1213 | if(!exclude) |
1057 | { |
1214 | { |
1058 | for(int j = 0; j < 10; j++) |
1215 | for(int j = 0; j < 10; j++) |
Line 1310... | Line 1467... | ||
1310 | // Check if pedestal is above the lower limit and sigma is smaller than the mean |
1467 | // Check if pedestal is above the lower limit and sigma is smaller than the mean |
1311 | if( (fparam[j] > pedesLow->widgetNE[0]->GetNumber()) && ((double)fparamerr[j]/fparam[j] < accError->widgetNE[0]->GetNumber()) ) |
1468 | if( (fparam[j] > pedesLow->widgetNE[0]->GetNumber()) && ((double)fparamerr[j]/fparam[j] < accError->widgetNE[0]->GetNumber()) ) |
1312 | { |
1469 | { |
1313 | // With the additional ADC offset, we can shift the mean values slightly, so they are not close to the X.5, but to the X.0 values |
1470 | // With the additional ADC offset, we can shift the mean values slightly, so they are not close to the X.5, but to the X.0 values |
1314 | meansel[nrfit] = fparam[j]+(adcOffset->widgetNE[0]->GetNumber()); |
1471 | meansel[nrfit] = fparam[j]+(adcOffset->widgetNE[0]->GetNumber()); |
- | 1472 | meanselerr[nrfit] = fparamerr[j]; |
|
1315 | sigmasel[nrfit] = fparam[j+1]; |
1473 | sigmasel[nrfit] = fparam[j+1]; |
1316 | nrfit++; |
1474 | nrfit++; |
1317 | } |
1475 | } |
1318 | } |
1476 | } |
1319 | 1477 | ||
Line 1485... | Line 1643... | ||
1485 | else |
1643 | else |
1486 | { |
1644 | { |
1487 | printf("BreakdownVolt(): The current separation measurements does not fall within acceptable errors.\n"); |
1645 | printf("BreakdownVolt(): The current separation measurements does not fall within acceptable errors.\n"); |
1488 | exclude = false; |
1646 | exclude = false; |
1489 | } |
1647 | } |
- | 1648 | ||
- | 1649 | printf("BreakdownVolt(): Calculated separation between peaks %d and %d: sep = %lf, seperr = %lf\n", (int)peakSepCalc->widgetNE[0]->GetNumber()-1, (int)peakSepCalc->widgetNE[0]->GetNumber(), sep[(int)peakSepCalc->widgetNE[0]->GetNumber()-1][p], seperr[(int)peakSepCalc->widgetNE[0]->GetNumber()-1][p]); |
|
1490 | 1650 | ||
1491 | // Write out parameters to a file |
1651 | // Write out parameters to a file |
1492 | fp = fopen(paramname, "a"); |
1652 | fp = fopen(paramname, "a"); |
1493 | if(exclude) |
1653 | if(exclude) |
1494 | { |
1654 | { |
Line 1538... | Line 1698... | ||
1538 | k = peakSepCalc->widgetNE[0]->GetNumber(); |
1698 | k = peakSepCalc->widgetNE[0]->GetNumber(); |
1539 | if(k < 4) |
1699 | if(k < 4) |
1540 | bdplot = new TGraphErrors(p, volt, sep[k-1], volterr, seperr[k-1]); |
1700 | bdplot = new TGraphErrors(p, volt, sep[k-1], volterr, seperr[k-1]); |
1541 | else |
1701 | else |
1542 | { |
1702 | { |
1543 | printf(" |
1703 | printf("BreakdownVolt(): Unsupported peak separation selected (%d).\n", k); |
1544 | return; |
1704 | return; |
1545 | } |
1705 | } |
1546 | 1706 | ||
1547 | bdplot->SetMarkerStyle(20); |
1707 | bdplot->SetMarkerStyle(20); |
1548 | bdplot->SetMarkerSize(0.4); |
1708 | bdplot->SetMarkerSize(0.4); |
Line 1575... | Line 1735... | ||
1575 | if(!cleanPlots) |
1735 | if(!cleanPlots) |
1576 | { |
1736 | { |
1577 | sprintf(ctemp, "#splitline{#Delta_{p}(U) = (%.2lf #pm %.2lf)#timesU + (%.2lf #pm %.3lf)}{U_{0} = %.2lf #pm %.3lf}", meansel[0], meanselerr[0], meansel[1], meanselerr[1], meansel[2], meansel[2]*(TMath::Abs(meanselerr[0]/meansel[0]) + TMath::Abs(meanselerr[1]/meansel[1])) ); |
1737 | sprintf(ctemp, "#splitline{#Delta_{p}(U) = (%.2lf #pm %.2lf)#timesU + (%.2lf #pm %.3lf)}{U_{0} = %.2lf #pm %.3lf}", meansel[0], meanselerr[0], meansel[1], meanselerr[1], meansel[2], meansel[2]*(TMath::Abs(meanselerr[0]/meansel[0]) + TMath::Abs(meanselerr[1]/meansel[1])) ); |
1578 | latex = new TLatex(); |
1738 | latex = new TLatex(); |
1579 | latex->SetTextSize(0.039); |
1739 | latex->SetTextSize(0.039); |
1580 | latex->DrawLatex(volt[0], 0.97*sep[0][sortindex[ |
1740 | latex->DrawLatex(volt[0], 0.97*sep[0][sortindex[1]], ctemp); |
1581 | printf("#Delta_{p}(U) = (%.6lf #pm %.8lf)#timesU + (%.6lf #pm %.8lf)}{U_{0} = %.6lf #pm %.8lf\n", meansel[0], meanselerr[0], meansel[1], meanselerr[1], meansel[2], meansel[2]*(TMath::Abs(meanselerr[0]/meansel[0]) + TMath::Abs(meanselerr[1]/meansel[1])) ); |
1741 | printf("#Delta_{p}(U) = (%.6lf #pm %.8lf)#timesU + (%.6lf #pm %.8lf)}{U_{0} = %.6lf #pm %.8lf\n", meansel[0], meanselerr[0], meansel[1], meanselerr[1], meansel[2], meansel[2]*(TMath::Abs(meanselerr[0]/meansel[0]) + TMath::Abs(meanselerr[1]/meansel[1])) ); |
1582 | } |
1742 | } |
1583 | else |
1743 | else |
1584 | printf("#Delta_{p}(U) = (%.2lf #pm %.2lf)#timesU + (%.2lf #pm %.3lf)}{U_{0} = %.2lf #pm %.3lf\n", meansel[0], meanselerr[0], meansel[1], meanselerr[1], meansel[2], meansel[2]*(TMath::Abs(meanselerr[0]/meansel[0]) + TMath::Abs(meanselerr[1]/meansel[1])) ); |
1744 | printf("#Delta_{p}(U) = (%.2lf #pm %.2lf)#timesU + (%.2lf #pm %.3lf)}{U_{0} = %.2lf #pm %.3lf\n", meansel[0], meanselerr[0], meansel[1], meanselerr[1], meansel[2], meansel[2]*(TMath::Abs(meanselerr[0]/meansel[0]) + TMath::Abs(meanselerr[1]/meansel[1])) ); |
1585 | 1745 | ||
Line 1750... | Line 1910... | ||
1750 | gCanvas = tempAnalysisCanvas->GetCanvas(); |
1910 | gCanvas = tempAnalysisCanvas->GetCanvas(); |
1751 | 1911 | ||
1752 | double range[4]; |
1912 | double range[4]; |
1753 | TGraph2D *gScan2D; |
1913 | TGraph2D *gScan2D; |
1754 | gScan2D = new TGraph2D(); |
1914 | gScan2D = new TGraph2D(); |
1755 | gScan2D->SetName("surfscan"); |
- | |
1756 | range[0] = TMath::MinElement(nrfiles, surfx); |
1915 | range[0] = TMath::MinElement(nrfiles, surfx); |
1757 | range[1] = TMath::MaxElement(nrfiles, surfx); |
1916 | range[1] = TMath::MaxElement(nrfiles, surfx); |
1758 | range[2] = TMath::MinElement(nrfiles, surfy); |
1917 | range[2] = TMath::MinElement(nrfiles, surfy); |
1759 | range[3] = TMath::MaxElement(nrfiles, surfy); |
1918 | range[3] = TMath::MaxElement(nrfiles, surfy); |
1760 | 1919 | ||
Line 1853... | Line 2012... | ||
1853 | } |
2012 | } |
1854 | else if(edit == 1) |
2013 | else if(edit == 1) |
1855 | { |
2014 | { |
1856 | gCanvas->Modified(); |
2015 | gCanvas->Modified(); |
1857 | gCanvas->Update(); |
2016 | gCanvas->Update(); |
1858 | - | ||
1859 | UpdateIntegrateSurface(-1); |
- | |
1860 | } |
2017 | } |
1861 | } |
2018 | } |
1862 | } |
2019 | } |