Index: fact/tools/rootmacros/gainfit2.C =================================================================== --- fact/tools/rootmacros/gainfit2.C (revision 12712) +++ fact/tools/rootmacros/gainfit2.C (revision 12713) @@ -1,2 +1,22 @@ +//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +//+ Root Macro: GainFit2 for FACT + +//+---------------------------------------------------------------------------+ +//+ by: jens.buss@tu-dortmund.de + +//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +// +// - analyse gain of pixels from Dark-Count-Photon-Spektrum (pedestal on) +// - read a root-file produced by "fpeak.C" and get Amplitudespektrum +// of all pixel +// - read a root-file produced by "fbsl.C" and get baseline +// - calculate pixelgain from multi-gaussian-fit +// over each pixels amplitudespektrum +// - draw distribution of gains +// - draw gain vs. pixel +// - find Pixels with deviating gainvalue (e.g. dead , crazy Pixels) +// and write their number to log-file (e.g. date_run_gainfit.txt) +// - save histograms into root-file (e.g. date_run_gainfit.root) +// +//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + #include #include @@ -29,6 +49,8 @@ Double_t amplSD; Double_t amplDT; - Double_t amplMean; + Double_t cross; Double_t ampl_fspek_fit; + Double_t ampl_fspek_fit2; + Double_t bsl; bool crazy; } proj; @@ -42,4 +64,6 @@ dist_t distribution; +//Verbosity for Histograms +TString verbos[3]={"Q","", ""}; //----------------------------------------------------------------------------- // Decleration of Histograms @@ -49,9 +73,10 @@ TH1F *hAmplDistribution = NULL; TH1F *hAmplDistribution2 = NULL; -TH1F *hPixelGainFit = NULL; +TH1F *hGainFitBsl = NULL; //Lists of Histograms that have to be saved TObjArray hList; TObjArray hListPeakSpektra; +TObjArray hListCrazy; //----------------------------------------------------------------------------- @@ -63,5 +88,5 @@ void BookHistos(); void SaveHistograms(const char * ); -bool SearchAnormal(UInt_t, Double_t ); +bool SearchAnormal(UInt_t, Double_t, int verbosityLevel ); //----------------------------------------------------------------------------- @@ -73,8 +98,9 @@ const char * InputRootFileName = "../analysis/fpeak/20111109_006/20111109_006_fpeak.root", const char * InputBaselineFileName = "../analysis/fpeak/20111109_006/20111109_006_fbsl.root", - const char * OutTextFileName = "../analysis/fpeak/20111109_006_gainfit_test.txt", - const char * RootOutFileName = "../analysis/fpeak/20111109_006_gainfit_test.root", + const char * OutTextFileName = "../analysis/fpeak/20111109_006/20111109_006_gainfit.txt", + const char * RootOutFileName = "../analysis/fpeak/20111109_006/20111109_006_gainfit.root", bool showHistos = false, - bool debug = false) + bool debug = false, + int verbosityLevel = 3) { @@ -106,7 +132,7 @@ TF1 *g3 = new TF1("g3","gaus",25,35); TF1 *g4 = new TF1("g4","gaus"); - TF1 *total = new TF1("total","gaus(0)+gaus(3)+gaus(6)",7,32); + //TF1 *total = new TF1("total","gaus(0)+gaus(3)+gaus(6)",7,60); TF1 *fspek_fit = new TF1("fspek_fit", spek_fit, 5, 60, 5); - TF1 *fspek_fit2 = new TF1("fspek_fit2", spek_fit2, 5, 60, 5); + TF1 *fspek_fit2 = new TF1("fspek_fit2", spek_fit2, 5, 60, 4); //Create Canvases @@ -128,9 +154,9 @@ - cGainFit->cd(2); + cGainFit->cd(3); hbaseline->SetTitle("Baselinedistribution of all Pixels"); - hbaseline->SetAxisRange(-5, 5, "X"); - //hbaseline->SetAxisRange(0, 100, "Y"); - hbaseline->Fit("gaus"); + Double_t Xrange = hbaseline->GetMean(1); + hbaseline->SetAxisRange(Xrange-1.5, Xrange + 1.5, "X"); + hbaseline->SetAxisRange(0, (hbaseline->GetBinContent(hbaseline->GetMaximumBin())+300), "Y"); hbaseline->Draw(); @@ -154,7 +180,9 @@ gr2->SetMarkerStyle(2); gr1->SetMarkerColor(kRed); - gr2->SetMarkerColor(kBlack); + gr2->SetMarkerColor(kBlue); gr1->SetMarkerSize(0.5); gr2->SetMarkerSize(0.5); + gr1->SetTitle("G-APD Amplification vs. Channel (with Baseline Parameter)"); + gr2->SetTitle("G-APD Amplification vs. Channel (without Baseline Parameter)"); TMultiGraph *mg = new TMultiGraph(); mg->Add(gr1); @@ -165,23 +193,27 @@ //----------------------------------------------------------------------------- - cout << "Number of Pixels: " << NumberOfPixels << endl; + if (verbosityLevel > 1) cout << "Number of Pixels: " << NumberOfPixels << endl; // Begin of Loop for (UInt_t pixel = 0; pixel < NumberOfPixels; pixel++){ - //Title + channel[ pixel ].PID = pixel+0; + + //Title TString histotitle; histotitle += "hPixelPro_"; - histotitle += pixel + 1 ; - + histotitle += pixel+0 ; + + if (verbosityLevel > 1){ cout << "-----------------------------------------------------------" << endl; - cout << " PID: " << pixel+1 << endl; + cout << " PID: " << pixel+0 << endl; cout << "-----------------------------------------------------------" << endl; + } //Amplitude Spectrum of a Pixel hPixelAmplSpek[ pixel ] = new TH1D; - hPixelAmplSpek[ pixel ] = hAmplSpek->ProjectionY( histotitle , pixel+1, pixel+1); //Projectipon of each Pixel out of Ampl.Spectrum - hPixelAmplSpek[ pixel ]->SetAxisRange(0, 40, "X"); + hPixelAmplSpek[ pixel ] = hAmplSpek->ProjectionY( histotitle , pixel+0, pixel+0); //Projectipon of each Pixel out of Ampl.Spectrum + hPixelAmplSpek[ pixel ]->SetAxisRange(0, 60, "X"); histotitle = "Amplitude Spectrum of Pixel #"; - histotitle += pixel + 1 ; + histotitle += pixel+0 ; hPixelAmplSpek[ pixel ]->SetTitle(histotitle); @@ -194,26 +226,39 @@ //----------------------------------------------------------------------------- - total->SetLineColor(2); - fspek_fit->SetLineColor(3); + //total->SetLineColor(2); + fspek_fit->SetLineColor(kRed); fspek_fit2->SetLineColor(kBlue); fspek_fit2->SetLineStyle(2); - cGainFit->cd(3); + cGainFit->cd(4); gPad->SetLogy(1); - hPixelAmplSpek[ pixel ]->Fit(g1,"R"); - hPixelAmplSpek[ pixel ]->Fit(g2,"R+"); - hPixelAmplSpek[ pixel ]->Fit(g3,"R+"); + hPixelAmplSpek[ pixel ]->Fit(g1,verbos[verbosityLevel-1]+"R"); + hPixelAmplSpek[ pixel ]->Fit(g2,verbos[verbosityLevel-1]+"0R+"); + hPixelAmplSpek[ pixel ]->Fit(g3,verbos[verbosityLevel-1]+"0R+"); g1->GetParameters(&par[0]); g2->GetParameters(&par[3]); g3->GetParameters(&par[6]); + //total->SetParameters(par); + + //hPixelAmplSpek[ pixel ]->Fit(total,verbos[verbosityLevel-1]+"R+"); + //channel[ pixel ].amplSD = total->GetParameter(4) - total->GetParameter(1); + //channel[ pixel ].amplDT = total->GetParameter(7) - total->GetParameter(4); + //channel[ pixel ].amplMean = fspek_fit->GetParameter(1); + + par[4] = -1; fspek_fit->SetParameters(par); + cout << "spek_fit" << endl; + hPixelAmplSpek[ pixel ]->Fit(fspek_fit,verbos[verbosityLevel-1]+"R+"); + channel[ pixel ].ampl_fspek_fit = fspek_fit->GetParameter(1); + channel[ pixel ].bsl = fspek_fit->GetParameter(4); + cout << "Baseline:" << channel[ pixel ].bsl << endl; + fspek_fit2->SetParameters(par); - total->SetParameters(par); - hPixelAmplSpek[ pixel ]->Fit(total,"R+"); - hPixelAmplSpek[ pixel ]->Fit(fspek_fit,"R+"); - hPixelAmplSpek[ pixel ]->Fit(fspek_fit2,"R+"); - + cout << "fspek_fit2" << endl; + hPixelAmplSpek[ pixel ]->Fit(fspek_fit2,verbos[verbosityLevel-1]+"R+"); + channel[ pixel ].ampl_fspek_fit2 = fspek_fit2->GetParameter(1); + if( showHistos ){ hPixelAmplSpek[ pixel ]->Draw(); @@ -222,17 +267,6 @@ -//----------------------------------------------------------------------------- -// Save Parameters into vetor -//----------------------------------------------------------------------------- - -// total->GetParameters(channel[ pixel ].fitParameters); - channel[ pixel ].PID = pixel+1; - channel[ pixel ].amplSD = total->GetParameter(4) - total->GetParameter(1); - channel[ pixel ].amplDT = total->GetParameter(7) - total->GetParameter(4); - channel[ pixel ].amplMean = fspek_fit->GetParameter(1); - channel[ pixel ].ampl_fspek_fit = fspek_fit->GetParameter(1); - -// channel[ pixel ].tripple_peakSigm = total->GetParameter(8); -// channel[ pixel ].ampl = total->GetParameter(1) - total->GetParameter(4); + + //----------------------------------------------------------------------------- @@ -240,9 +274,9 @@ //----------------------------------------------------------------------------- - cGainFit->cd(4); + cGainFit->cd(2); gPad->SetLogy(1); g4->SetLineColor(2); - hAmplDistribution->Fill( channel[ pixel ].amplSD ) ; - hAmplDistribution->Fit(g4); + hAmplDistribution->Fill( channel[ pixel ].ampl_fspek_fit ) ; + hAmplDistribution->Fit(g4, verbos[verbosityLevel-1]); distribution.sigma = g4->GetParameter(2); distribution.mean = g4->GetParameter(1); @@ -251,14 +285,14 @@ cGainFit->cd(5); gPad->SetLogy(1); - g4->SetLineColor(2); - hAmplDistribution2->Fill(channel[ pixel ].amplMean) ; - hAmplDistribution2->Fit(g4); - - cGainFit->cd(6); - hPixelGainFit->SetBinContent(pixel, channel[ pixel ].amplSD); - + g4->SetLineColor(kBlue); + hAmplDistribution2->Fill(channel[ pixel ].ampl_fspek_fit2) ; + hAmplDistribution2->Fit(g4, verbos[verbosityLevel-1]); + + cGainFit->cd(3); + hGainFitBsl->Fill(channel[ pixel ].bsl); + hGainFitBsl->Draw("same"); if( showHistos ){ if (pixel%40==0){ - cGainFit->cd(4); + cGainFit->cd(2); hAmplDistribution->Draw(); cGainFit->cd(5); @@ -267,5 +301,5 @@ cGainFit->cd(6); - hPixelGainFit->Draw(); + hGainFitBsl->Draw(); } @@ -274,13 +308,13 @@ //----------------------------------------------------------------------------- - gr1->SetPoint(pixel, channel[pixel].PID, channel[pixel].amplSD); - gr2->SetPoint(pixel, channel[pixel].PID, channel[pixel].amplMean); + gr1->SetPoint(pixel, channel[pixel].PID, channel[pixel].ampl_fspek_fit); + gr2->SetPoint(pixel, channel[pixel].PID, channel[pixel].ampl_fspek_fit2); if( showHistos ){ if (pixel%60==0){ + cGraphs->cd(1); + gr2->Draw("APL"); cGraphs->cd(2); - gr2->Draw("APL"); - cGraphs->cd(1); gr1->Draw("APL"); } @@ -311,8 +345,8 @@ -// ( cGainFit->cd(3); +// ( cGainFit->cd(4); // hPixelAmplSpek[ NumberOfPixels ]->Draw(); - cGainFit->cd(4); + cGainFit->cd(2); hAmplDistribution->Draw(); @@ -321,12 +355,13 @@ cGainFit->cd(6); - hPixelGainFit->Draw(); + hGainFitBsl->Draw(); + hbaseline->Draw("same"); cGainFit->Modified(); cGainFit->Update(); + cGraphs->cd(1); + gr2->Draw("APL"); cGraphs->cd(2); - gr2->Draw("APL"); - cGraphs->cd(1); gr1->Draw("APL"); cGraphs->Modified(); @@ -351,5 +386,5 @@ ++it; } - cout << filename << endl; +if (verbosityLevel > 0) cout << filename << endl; ofstream out; @@ -362,9 +397,11 @@ for (UInt_t pixel=0; pixelSetXTitle("Amplification [mV]"); hAmplDistribution->SetYTitle("Counts"); @@ -406,5 +443,5 @@ //Amplification Distribution2 - Double2Tripple - hAmplDistribution2 = new TH1F("hAmplDistribution2","Distribution of G-APD Amplification(2+3)", 500, 0 , 20 ); + hAmplDistribution2 = new TH1F("hAmplDistribution2","Distribution of G-APD Amplification(without Baseline Parameter)", 500, 0 , 20 ); hAmplDistribution2->SetXTitle("Amplification [mV]"); hAmplDistribution2->SetYTitle("Counts"); @@ -412,10 +449,10 @@ hList.Add( hAmplDistribution2 ); - //Amplification vs Channels - hPixelGainFit = new TH1F("hPixelGainFit","Amplifications (SD) of each Channel", 1440, -0.5 , 1439.5 ); - hPixelGainFit->SetYTitle("Amplification [mV]"); - hPixelGainFit->SetXTitle("Channels"); - hPixelGainFit->SetLineColor(2); - hList.Add( hPixelGainFit ); + //cross vs Channels + hGainFitBsl = new TH1F("hGainFitBsl","BaselineDistribution from plot", 400,-99.5,100.5); + hGainFitBsl->SetYTitle("Counts"); + hGainFitBsl->SetXTitle("Voltage [mV]"); + hGainFitBsl->SetLineColor(2); + hList.Add( hGainFitBsl ); } @@ -432,13 +469,36 @@ tf.mkdir("PeakSpektra"); tf.cd("PeakSpektra"); // go to new subdirectory hListPeakSpektra.Write(); // write histos into subdirectory - + tf.mkdir("CrazyPixels"); tf.cd("CrazyPixels"); // go to new subdirectory + hListCrazy.Write(); // write histos into subdirector tf.Close(); // close the file } //----------------------------------------------------------------------------- -// MulitOrder MultiGauß Function +// Fit-Function: Sum of several Gauß-Function with falling Maximum Amplitude //----------------------------------------------------------------------------- Double_t spek_fit(Double_t *x, Double_t *par){ + Double_t arg = 0; + Double_t arg2 = 0; + Double_t arg3 = 0; + Double_t arg4 = 0; + Double_t cross = 0; + if (par[0] != 0) cross = par[3]/par[0]; + if (par[2] != 0) arg = (x[0] - par[1]-par[4])/par[2]; + if (par[2] != 0) arg2 = (x[0] - 2*par[1]-par[4])/par[2]; + if (par[2] != 0) arg3 = (x[0] - 3*par[1]-par[4])/par[2]; + if (par[2] != 0) arg4 = (x[0] - 4*par[1]-par[4])/par[2]; + Double_t peak1 = par[0]*TMath::Exp(-0.5 * arg * arg); + Double_t peak2 = cross*par[0]*TMath::Exp(-0.5 * arg2 * arg2 /2); + Double_t peak3 = cross*cross*par[0]*TMath::Exp(-0.5 * arg3 * arg3 /3); + Double_t peak4 = cross*cross*cross*par[0]*TMath::Exp(-0.5 * arg4 * arg4 /4); + return peak1 + peak2 + peak3 + peak4; +} + +//----------------------------------------------------------------------------- +// Fit-Function: Sum of several Gauß-Function with overlying potential decay +//----------------------------------------------------------------------------- + +Double_t spek_fit2(Double_t *x, Double_t *par){ Double_t arg = 0; @@ -447,40 +507,21 @@ Double_t arg4 = 0; Double_t cross = 1; - if (par[0] != 0) cross = par[4]/par[0]; + if (par[0] != 0) cross = par[3]/par[0]; if (par[2] != 0) arg = (x[0] - par[1])/par[2]; if (par[2] != 0) arg2 = (x[0] - 2*par[1])/par[2]; if (par[2] != 0) arg3 = (x[0] - 3*par[1])/par[2]; if (par[2] != 0) arg4 = (x[0] - 4*par[1])/par[2]; - cross = par[4]/par[0]; Double_t peak1 = par[0]*TMath::Exp(-0.5 * arg * arg); Double_t peak2 = cross*par[0]*TMath::Exp(-0.5 * arg2 * arg2 /2); Double_t peak3 = cross*cross*par[0]*TMath::Exp(-0.5 * arg3 * arg3 /3); Double_t peak4 = cross*cross*cross*par[0]*TMath::Exp(-0.5 * arg4 * arg4 /4); - //Double_t decayoverlay = TMath::Exp(-x[0]); return peak1 + peak2 + peak3 + peak4; } -Double_t spek_fit2(Double_t *x, Double_t *par){ - - Double_t arg = 0; - Double_t arg2 = 0; - Double_t arg3 = 0; - Double_t arg4 = 0; - Double_t cross = 1; - if (par[0] != 0) cross = par[4]/par[0]; - if (par[2] != 0) arg = (x[0] - par[1])/par[2]; - if (par[2] != 0) arg2 = (x[0] - 2*par[1])/par[2]; - if (par[2] != 0) arg3 = (x[0] - 3*par[1])/par[2]; - if (par[2] != 0) arg4 = (x[0] - 4*par[1])/par[2]; - cross = par[4]/par[0]; - Double_t peak1 = par[0]*TMath::Exp(-0.5 * arg * arg); - Double_t peak2 = cross*par[0]*TMath::Exp(-0.5 * arg2 * arg2 /2); - Double_t peak3 = cross*cross*par[0]*TMath::Exp(-0.5 * arg3 * arg3 /3); - Double_t peak4 = cross*cross*cross*par[0]*TMath::Exp(-0.5 * arg4 * arg4 /4); - Double_t decayoverlay = TMath::Exp(-x[0]*TMath::Log(1.0001)); - return peak1 + peak2 + peak3 + peak4 + decayoverlay; -} - -bool SearchAnormal(UInt_t pixel, Double_t distance){ //// MACHT ALLES KEINEN SINN WENN MAN KEINE BETRÄGE BETRACHTET +//----------------------------------------------------------------------------- +// Function: Search for Pixel that do not have a compareable gain to the Other +//----------------------------------------------------------------------------- + +bool SearchAnormal(UInt_t pixel, Double_t distance, int verbosityLevel){ Double_t check = channel[ pixel ].ampl_fspek_fit; check = TMath::Sqrt(check*check); @@ -489,37 +530,7 @@ } else{ - cout << "Not normal behaviotruer in Pixel: " << pixel +1 << endl; + if (verbosityLevel > 2) cout << "Not normal behaviotruer in Pixel: " << pixel +1 << endl; return true; } } -/* - - - - - - - */ - - -/* -Double_t spek_fit(Double_t *x, Double_t *par, Uint_t *maxorder) - { - Double_t exp = 0; - - //Double_t arg2 = 0; - //if (par[2] != 0) arg = (v[0] - par[1])/par[2]; - //if (par[5] != 0) arg2 = (v[0] - par[4])/par[5]; - - //Double_t fitval = par[0]*TMath::Exp(-0.5*arg*arg) + par[3]*TMath::Exp(-0.5*arg2*arg2); - - for(int order = 0; order < maxorder; order++){ - peak[order]=par[0]*par[4]/(order+1)*((x[0]-(order+1)*par[1])*(x[0]-(order+1)*par[1]))/(par[2]*par[2]); - exp = exp + peak[order]; - } - - Double_t fitval = TMath::Exp(-0.5*exp); - return fitval; - } -*/