455 const float gainSeed = 55.0;
456 const Double_t tSpectrumSigma = 1.75;
457 const Double_t tSpectrumThreshold = 0.18;
458 const size_t nPeakMax = 5;
461 h->GetXaxis()->SetRangeUser(350,700);
464 TCanvas* cspec =
new TCanvas();
467 TSpectrum *
s =
new TSpectrum();
469 int nPeak = s->Search(
h,tSpectrumSigma,
"",tSpectrumThreshold);
475 float* peakXsqr =
new float[5];
476 float* peakNdf =
new float[5];
477 float* peakSigma =
new float[5];
478 float* peakSigmaErr =
new float[5];
479 float* peakNorm =
new float[5];
480 float* peakNormErr =
new float[5];
481 float* peakMean =
new float[5];
482 float* peakMeanErr =
new float[5];
486 for(
size_t j=0; j<5; j++){
488 peakXsqr[j] = FLT_MAX;
489 peakNdf[j] = FLT_MAX;
490 peakSigma[j] = FLT_MAX;
491 peakSigmaErr[j] = FLT_MAX;
492 peakNorm[j] = FLT_MAX;
493 peakNormErr[j] = FLT_MAX;
494 peakMean[j] = FLT_MAX;
495 peakMeanErr[j] = FLT_MAX;
505 Double_t *peaks = s->GetPositionX();
514 for(
int p=0;
p<nPeak-1;
p++)
515 if(peaks[
p]>peaks[
p+1])
522 }
while( nchanges != 0 );
532 float gy[11], gey[11];
533 int peak_offset = round((peaks[0]-peds[chan])/gainSeed);
535 for (
int j=0; j<nPeak&&j<10; j++)
537 if(peaks[j]<0)
std::cout <<
"bad peak value: " << peaks[j] << std::endl;
538 x[j] = j+peak_offset;
549 gey[0] = pwidths[chan];
556 for (
int g=0 ;
g<nPeak&&
g<(int)nPeakMax;
g++)
559 TF1 *gfit =
new TF1(
"gfit",
"gaus",y[
g]-20,y[
g]+20);
560 gfit->SetParameter(0,
h->GetBinContent(
h->FindBin(y[
g])));
561 gfit->SetParameter(1,y[g]);
562 gfit->SetParameter(2,12);
563 gfit->SetParLimits(0,0,20000);
564 gfit->SetParLimits(1,y[g]-15,y[g]+15);
565 gfit->SetParLimits(2,8,40);
574 if( y[g]<
h->GetBinLowEdge(1)+15) nplow++;
575 if( y[g]>
h->GetBinLowEdge(
h->GetNbinsX())-15) nphigh++;
578 if( y[g]>
h->GetBinLowEdge(1)+15 && y[
g]<
h->GetBinLowEdge(
h->GetNbinsX())-15&&
abs(y[g] - gfit->GetParameter(1)) < 15)
581 if(g!=0 && g!=nPeak-1
582 && (y[g+1]-y[g]<gainSeed*0.7 || y[g]-y[g-1]<gainSeed*0.7))
585 for (
int j=g; j<nPeak; j++) x[j] = x[j]-1;
590 if(g!=0&&g!=nPeak-1&&y[g+1]-y[g]<gainSeed*1.2&&y[g]-y[g-1]>gainSeed*1.5)
592 for (
int j=g; j<nPeak; j++) x[j] = x[j]+1;
596 if(g==nPeak-1 && (y[g]-y[g-1])/gainSeed >1.5) {
597 x[
g]+=(int)((y[g]-y[g-1])/gainSeed);
601 peakXsqr[gg] = gfit->GetChisquare();
602 peakNdf[gg] = gfit->GetNDF();
603 peakMean[gg] = gfit->GetParameter(1);
604 peakMeanErr[gg] = gfit->GetParError(1);
605 peakNorm[gg] = gfit->GetParameter(0);
606 peakNormErr[gg] = gfit->GetParError(0);
607 peakSigma[gg] = gfit->GetParameter(2);
608 peakSigmaErr[gg] = gfit->GetParError(2);
611 gy[gg] = gfit->GetParameter(1);
612 gey[gg] = sqrt(gx[gg]+gey[0]*gey[0]);
620 if (nplow>0)
for (
int i=1; i<gg; i++) gx[i] = gx[i]-(nplow-1);
623 TGraphErrors* gr_mean =
new TGraphErrors(gg,gx,gy,0,gey);
626 TF1 *fit =
new TF1(
"fit",
"[0] + [1]*x",gx[0]-0.25,gx[gg-1]+0.25);
629 fit->SetParName(1,
"Gain");
630 fit->SetParName(0,
"Pedestal");
631 fit->SetParameter(1,gainSeed);
632 fit->SetParameter(0,peds[chan]);
633 fit->SetParLimits(1,gainSeed-20,gainSeed+20);
634 fit->SetParLimits(0,peds[chan]*0.8,peds[chan]*1.2);
638 gr_mean->Fit(fit,
"QR");
642 if (fit->GetChisquare()/fit->GetNDF()>5.0)
645 chisqr=fit->GetChisquare();
646 for(
int i=1; i<gg; i++) gx[i]+=1;
647 gr_mean =
new TGraphErrors(gg,gx,gy,0,gey);
648 fit->SetRange(gx[0]-0.25,gx[gg-1]+0.25);
649 gr_mean->Fit(fit,
"QR");
650 if (fit->GetChisquare()<chisqr) chisqr=fit->GetChisquare();
653 for(
int i=1; i<gg; i++) gx[i]-=2;
654 gr_mean =
new TGraphErrors(gg,gx,gy,0,gey);
655 fit->SetRange(gx[0]-0.25,gx[gg-1]+0.25);
656 gr_mean->Fit(fit,
"QR");
657 if (fit->GetChisquare()<chisqr) chisqr=fit->GetChisquare();
660 for(
int i=1; i<gg; i++) gx[i]+=1;
661 gr_mean =
new TGraphErrors(gg,gx,gy,0,gey);
662 fit->SetRange(gx[0]-0.25,gx[gg-1]+0.25);
663 gr_mean->Fit(fit,
"QR");
672 gr_mean->SetTitle(grtitle.c_str());
673 gr_mean->GetXaxis()->SetTitle(
"Peak #");
674 gr_mean->GetYaxis()->SetTitle(
"ADC value");
675 gr_mean->SetMarkerColor(4);
676 gr_mean->SetMarkerStyle(20);
677 gr_mean->SetFillColor(0);
679 gr_mean->GetXaxis()->SetRangeUser(gx[0]-0.5,gx[gg-1]+0.5);
684 gStyle->SetOptStat(0100);
685 gStyle->SetOptFit(1111);
689 TCanvas *c2 =
new TCanvas();
693 gStyle->SetStatX(0.5);
694 gStyle->SetStatY(0.9);
695 gStyle->SetStatH(0.15);
696 gStyle->SetStatW(0.2);
699 gr_mean->Draw(
"ALP");
705 statsarr[0][0] = fit->GetParameter(1);
706 statsarr[1][0] = fit->GetParError(1);
707 statsarr[2][0] = fit->GetParameter(0);
708 statsarr[3][0] = fit->GetParError(0);
709 statsarr[4][0] = fit->GetChisquare();
710 statsarr[5][0] = fit->GetNDF();
711 statsarr[6][0] = nPeak;
712 for(
size_t i=0; i<5; i++) {
713 statsarr[7][i] = peakXsqr[i];
714 statsarr[8][i] = peakMean[i];
715 statsarr[9][i] = peakMeanErr[i];
716 statsarr[10][i] = peakNorm[i];
717 statsarr[11][i] = peakNormErr[i];
718 statsarr[12][i] = peakSigma[i];
719 statsarr[13][i] = peakSigmaErr[i];
720 statsarr[14][i] = peakNdf[i];
725 delete[] peakMeanErr;
727 delete[] peakNormErr;
729 delete[] peakSigmaErr;
736 TString suff =
".png";
738 TString fname2 =
fname;
743 if(chan<10) fname+=
"0";
744 if(chan<10) fname2+=
"0";
746 fname2+=
"_peak_fit_spec";
750 while (!gSystem->AccessPathName(fname))
752 fname.Remove(fname.Length()-suff.Length(),fname.Length());
753 fname.Remove(fname.Length()-1,fname.Length());
758 while (!gSystem->AccessPathName(fname2))
760 fname2.Remove(fname2.Length()-suff.Length(),fname2.Length());
761 fname2.Remove(fname2.Length()-1,fname2.Length());
767 TImage *img = TImage::Create();
769 img->WriteImage(fname);
773 img->WriteImage(fname2);
process_name opflash particleana ie x
process_name opflash particleana ie ie y
float * GetPedSigma() const
std::string to_string(WindowPattern const &pattern)
then echo File list $list not found else cat $list while read file do echo $file sed s
BEGIN_PROLOG could also be cout