Subversion Repositories f9daq

Rev

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
   // Start if we select at least one file
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(nrfiles > 0)
698
   if(relPde->widgetChBox[1]->IsDown())
695
   {
699
   {
696
      for(int i = 0; i < (int)nrfiles; i++)
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
            printf("PhotonMu(): Only one file selected. Not running analysis, just showing the fit.\n");
704
      if(strcmp("", darkRun->widgetTE->GetText()) == 0)
701
 
705
      {
702
            // Replot the spectrum on analysisCanvas and do not close the input file
706
         printf("PhotonMu(): Error! The no peak structure option needs a dark histogram.\n");
703
            DisplayHistogram( (char*)(files->At(i)->GetTitle()), 0, 1);
707
         delete[] integralCount;
704
            analysisCanvas->GetCanvas()->Modified();
708
         delete[] integralPedestal;
705
            analysisCanvas->GetCanvas()->Update();
709
         delete[] angle;
706
       
710
         delete[] pdeval;
707
            // Get the spectrum
711
         delete[] muval;
708
            histtemp = (TH1F*)analysisCanvas->GetCanvas()->GetPrimitive(histname);
712
         return;
709
            npeaks = 15;
713
      }
710
            double par[300];
714
 
711
            spec = new TSpectrum(npeaks);
715
   }
712
            // Find spectrum background
716
 
713
            histback = spec->Background(histtemp, (int)fitInter->widgetNE[0]->GetNumber(), "same");
717
   // Start if we select at least one file
714
            // Clone histogram and subtract background from it if we select that option
718
   if(nrfiles > 0)
715
            h2 = (TH1F*)histtemp->Clone("h2");
719
   {
716
            if(fitChecks->widgetChBox[0]->IsDown())
720
      // Find the pedestal peak for the dark histogram, and use it for all if there are no peaks
717
               h2->Add(histback, -1);
721
      if(nopeaks != -1)
718
            // Search for the peaks
722
      {
719
            int found = spec->Search(h2, fitSigma->widgetNE[0]->GetNumber(), "goff", fitTresh->widgetNE[0]->GetNumber() );
723
         // Replot the spectrum on analysisCanvas and do not close the input file
720
            printf("PhotonMu(): Found %d candidates to fit.\n",found);
724
         DisplayHistogram( (char*)(darkRun->widgetTE->GetText()), 0, 1);
721
            npeaks = found;
725
         analysisCanvas->GetCanvas()->Modified();
722
   
726
         analysisCanvas->GetCanvas()->Update();
723
            // Set initial peak parameters
727
       
724
            xpeaks = spec->GetPositionX();
728
         // Get the spectrum
725
            for(j = 0; j < found; j++)
729
         histtemp = (TH1F*)analysisCanvas->GetCanvas()->GetPrimitive(histname);
726
            {
730
         npeaks = 15;
727
               float xp = xpeaks[j];
731
         double par[300];
728
               int bin = h2->GetXaxis()->FindBin(xp);
732
         spec = new TSpectrum(npeaks);
729
               float yp = h2->GetBinContent(bin);
733
         // Find spectrum background
730
               par[3*j] = yp;
734
         histback = spec->Background(histtemp, (int)fitInter->widgetNE[0]->GetNumber(), "same");
731
               par[3*j+1] = xp;
735
         // Clone histogram and subtract background from it if we select that option
732
               par[3*j+2] = (double)fitSigma->widgetNE[0]->GetNumber();
736
         h2 = (TH1F*)histtemp->Clone("h2");
733
            }
737
         if(fitChecks->widgetChBox[0]->IsDown())
734
         
738
            h2->Add(histback, -1);
735
            // Fit the histogram
739
         // Search for the peaks
736
            fit = new TF1("fit", FindPeaks, adcRange->widgetNE[0]->GetNumber(), adcRange->widgetNE[1]->GetNumber(), 3*npeaks);
740
         int found = spec->Search(h2, fitSigma->widgetNE[0]->GetNumber(), "goff", fitTresh->widgetNE[0]->GetNumber() );
737
            TVirtualFitter::Fitter(histtemp, 3*npeaks);
741
         printf("PhotonMu(): Found %d candidates to fit.\n",found);
738
            fit->SetParameters(par);
742
         npeaks = found;
739
            fit->SetNpx(300);
743
   
740
            h2->Fit("fit","Q");
744
         // Set initial peak parameters
741
            // Get the fitted parameters
745
         xpeaks = spec->GetPositionX();
742
            fittingfunc = h2->GetFunction("fit");
746
         for(j = 0; j < found; j++)
743
            fparam = fittingfunc->GetParameters();
747
         {
744
            fparamerr = fittingfunc->GetParErrors();
748
            float xp = xpeaks[j];
745
   
749
            int bin = h2->GetXaxis()->FindBin(xp);
746
            // Gather the parameters (mean peak value for now)
750
            float yp = h2->GetBinContent(bin);
747
            int j = 1;
751
            par[3*j] = yp;
748
            int nrfit = 0;
752
            par[3*j+1] = xp;
749
            while(1)
753
            par[3*j+2] = (double)fitSigma->widgetNE[0]->GetNumber();
750
            {
754
         }
751
               if( (fparam[j] < 1.E-30) || (nrfit > 8) )
755
         
752
                  break;
756
         // Fit the histogram
753
               else
757
         fit = new TF1("fit", FindPeaks, adcRange->widgetNE[0]->GetNumber(), adcRange->widgetNE[1]->GetNumber(), 3*npeaks);
754
               {
758
         TVirtualFitter::Fitter(histtemp, 3*npeaks);
755
                  // Check if pedestal is above the lower limit and sigma is smaller than the mean
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( (npeaks > 1) && (nrfit > 1) )
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("BreakdownVold(): Unsupported peak separation selected (%d).\n", k);
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[p-1]], ctemp);
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
}