WebSVN – f9daq – Diff – /lab/sipmscan/trunk/src/analysis.cpp

       resol2d->widgetNE[1]->SetNumber(40);
+   }
    else if(analTab->GetCurrent() == 1)  // Relative PDE
+   {
       relPde->widgetChBox[0]->SetState(kButtonDown);
-      relPde->widgetChBox[1]->SetState(kButtonUp);
       midPeak->widgetChBox[0]->SetState(kButtonUp);
       zeroAngle->widgetNE[0]->SetNumber(0.00);
+   }
    else if(analTab->GetCurrent() == 2)  // Breakdown voltage
+   {
    float progVal = 0;
    analysisProgress->widgetPB->SetPosition(progVal);
    gVirtualX->Update(1);
    // Check if the checkbox for no peaks is selected - TODO: Still need the situation when we do not have a peaked ADC spectrum
-   if(relPde->widgetChBox[1]->IsDown())
+/*   if(relPde->widgetChBox[1]->IsDown())
+   {
       printf("PhotonMu(): ADC spectrum has no peak structure.\n");
       nopeaks = 1;
       // Error if there is no darkhist
          delete[] pdeval;
          delete[] muval;
          return;
+      }
-   }
+   }*/
    // Start if we select at least one file
    if(nrfiles > 0)
    {
       // Find the pedestal peak for the dark histogram, and use it for all if there are no peaks
       if(nopeaks != -1)
       {
          // Replot the spectrum on analysisCanvas and do not close the input file
          DisplayHistogram( (char*)(darkRun->widgetTE->GetText()), 0, 1);
          analysisCanvas->GetCanvas()->Modified();
          analysisCanvas->GetCanvas()->Update();
          // Get the spectrum
          histtemp = (TH1F*)analysisCanvas->GetCanvas()->GetPrimitive(histname);
          npeaks = 15;
          double par[300];
          spec = new TSpectrum(npeaks);
          // Find spectrum background
          histback = spec->Background(histtemp, (int)fitInter->widgetNE[0]->GetNumber(), "same");
          // Clone histogram and subtract background from it if we select that option
          h2 = (TH1F*)histtemp->Clone("h2");
          if(fitChecks->widgetChBox[0]->IsDown())
             h2->Add(histback, -1);
          // Search for the peaks
          int found = spec->Search(h2, fitSigma->widgetNE[0]->GetNumber(), "goff", fitTresh->widgetNE[0]->GetNumber() );
          printf("PhotonMu(): Found %d candidates to fit.\n",found);
          npeaks = found;
          // Set initial peak parameters
          xpeaks = spec->GetPositionX();
          for(j = 0; j < found; j++)
          {
             float xp = xpeaks[j];
             int bin = h2->GetXaxis()->FindBin(xp);
             float yp = h2->GetBinContent(bin);
             par[3*j] = yp;
             par[3*j+1] = xp;
             par[3*j+2] = (double)fitSigma->widgetNE[0]->GetNumber();
          }
          // Fit the histogram
          fit = new TF1("fit", FindPeaks, adcRange->widgetNE[0]->GetNumber(), adcRange->widgetNE[1]->GetNumber(), 3*npeaks);
          TVirtualFitter::Fitter(histtemp, 3*npeaks);
          fit->SetParameters(par);
          fit->SetNpx(300);
          h2->Fit("fit","Q");
          // Get the fitted parameters
          fittingfunc = h2->GetFunction("fit");
          fparam = fittingfunc->GetParameters();
          fparamerr = fittingfunc->GetParErrors();
          // Gather the parameters (mean peak value for now)
          int j = 1;
          int nrfit = 0;
          while(1)
          {
             if( (fparam[j] < 1.E-30) || (nrfit > 8) )
                break;
             else
             {
                // Check if pedestal is above the lower limit and sigma is smaller than the mean
                if( (fparam[j] > pedesLow->widgetNE[0]->GetNumber()) && ((double)fparamerr[j]/fparam[j] < accError->widgetNE[0]->GetNumber()) )
                {
                   // 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
                   meansel[nrfit] = fparam[j]+(adcOffset->widgetNE[0]->GetNumber());
                   sigmasel[nrfit] = fparam[j+1];
                   nrfit++;
+               }
+            }
+      {
          if( (nrfiles == 1) || (!multiSelect->widgetChBox[0]->IsDown()) )
+         {
             printf("PhotonMu(): Only one file selected. Not running analysis, just showing the fit.\n");
             // Replot the spectrum on analysisCanvas and do not close the input file
             DisplayHistogram( (char*)(files->At(i)->GetTitle()), 0, 1);
             analysisCanvas->GetCanvas()->Modified();
             analysisCanvas->GetCanvas()->Update();
             // Get the spectrum
             histtemp = (TH1F*)analysisCanvas->GetCanvas()->GetPrimitive(histname);
             npeaks = 15;
             double par[300];
             spec = new TSpectrum(npeaks);
             // Find spectrum background
             histback = spec->Background(histtemp, (int)fitInter->widgetNE[0]->GetNumber(), "same");
             // Clone histogram and subtract background from it if we select that option
             h2 = (TH1F*)histtemp->Clone("h2");
             if(fitChecks->widgetChBox[0]->IsDown())
                h2->Add(histback, -1);
             // Search for the peaks
             int found = spec->Search(h2, fitSigma->widgetNE[0]->GetNumber(), "goff", fitTresh->widgetNE[0]->GetNumber() );
             printf("PhotonMu(): Found %d candidates to fit.\n",found);
             npeaks = found;
             // Set initial peak parameters
             xpeaks = spec->GetPositionX();
             for(j = 0; j < found; j++)
             {
                float xp = xpeaks[j];
                int bin = h2->GetXaxis()->FindBin(xp);
                float yp = h2->GetBinContent(bin);
                par[3*j] = yp;
                par[3*j+1] = xp;
                par[3*j+2] = (double)fitSigma->widgetNE[0]->GetNumber();
             }
             // Fit the histogram
             fit = new TF1("fit", FindPeaks, adcRange->widgetNE[0]->GetNumber(), adcRange->widgetNE[1]->GetNumber(), 3*npeaks);
             TVirtualFitter::Fitter(histtemp, 3*npeaks);
             fit->SetParameters(par);
             fit->SetNpx(300);
             h2->Fit("fit","Q");
             // Get the fitted parameters
             fittingfunc = h2->GetFunction("fit");
             fparam = fittingfunc->GetParameters();
             fparamerr = fittingfunc->GetParErrors();
             // Gather the parameters (mean peak value for now)
             int j = 1;
             int nrfit = 0;
             while(1)
+            {
             exclude = false;
             // Get next point values (if zero value -> need to add the dark hist value again)
             pointest[10] = angle[i];
             pointest[11] = muval[i];
+//          printf("Last point = [%lf,%lf]\n", pointest[10], pointest[11]);
             // Check if next point has larger error than acceptable (if yes, set exclude signal to true), reuse meansel[0]
             meansel[0] = PointEstimate(5, pointest);    // PointEstimate only works with very small step size
             if(meansel[0] > accError->widgetNE[0]->GetNumber())
+            {
                printf("PhotonMu(): Point (%lf, %lf) excluded due to error: %lf\n", pointest[10], pointest[11], meansel[0]);
+               exclude = true;
+            }
+            if(isinf(meansel[0]))
+            {
+               printf("PhotonMu(): Point (%lf, %lf) excluded due to being infinite: %lf\n", pointest[10], pointest[11], meansel[0]);
                exclude = true;
+            }
             // Value with 0 angle and dark histogram are always needed, so should not be excluded
             if(i == darkhist)
                exclude = false;
             // If nothing excluded, pass the points in pointest variable like in a FIFO
             if(!exclude)
             {
                for(int j = 0; j < 10; j++)
+               {
                   if(DBGSIG) printf("PhotonMu(): j = %d: Old X value = %lf\n", j, pointest[j]);
                   pointest[j] = pointest[j+2];
                }
             }
             else
+            {
                for(int j = 0; j < 10; j++)
                   if(DBGSIG) printf("PhotonMu(): No j = %d: Old X value = %lf\n", j, pointest[j]);
+            }
+         {
             // First 5 points act as estimator points for next one
             pointest[2*m] = angle[i];
             pointest[2*m+1] = muval[i];
+         }
+         printf("Exclude signal = %d\n", (int)exclude);
          // Run only if we have a dark run histogram and middle pedestal peak estimation
          if( (darkhist != -1) && midPeak->widgetChBox[0]->IsDown() )
+         {
             if(DBGSIG) printf("PhotonMu(): m = %d, i = %d: muval = %lf, ", m, i, muval[i]);

Subversion Repositories f9daq

(root)/lab/sipmscan/trunk/src/analysis.cpp – Rev 173 → 181