Index: fact/tools/rootmacros/gainfit2.C =================================================================== --- fact/tools/rootmacros/gainfit2.C (revision 12690) +++ fact/tools/rootmacros/gainfit2.C (revision 12712) @@ -27,10 +27,18 @@ UInt_t PID; vector fitParameters[9]; - float amplSD; - float amplDT; - float amplMean; + Double_t amplSD; + Double_t amplDT; + Double_t amplMean; + Double_t ampl_fspek_fit; + bool crazy; } proj; vector channel; +typedef struct{ + Double_t constant; + Double_t mean; + Double_t sigma; +} dist_t; +dist_t distribution; //----------------------------------------------------------------------------- @@ -41,19 +49,19 @@ TH1F *hAmplDistribution = NULL; TH1F *hAmplDistribution2 = NULL; -TH1F *hAmplVsChannel = NULL; - -//Lists of Histograms taht are to be saved +TH1F *hPixelGainFit = NULL; + +//Lists of Histograms that have to be saved TObjArray hList; TObjArray hListPeakSpektra; - //----------------------------------------------------------------------------- // Functions //----------------------------------------------------------------------------- +Double_t spek_fit(Double_t *, Double_t *); +Double_t spek_fit2(Double_t *, Double_t *); void BookHistos(); void SaveHistograms(const char * ); -//Double_t spek_fit(Double_t *x, Double_t *par) - +bool SearchAnormal(UInt_t, Double_t ); //----------------------------------------------------------------------------- @@ -79,13 +87,18 @@ TFile * file = TFile::Open( InputRootFileName ); - - TH1D *hPixelProjection = NULL; //Baseline TH2F * hbaseline = (TH2F*)baseline->Get("histo_mean"); //Amplitude spectrum - TH2F *h = (TH2F*)file->Get("hAmplSpek_cfd"); - + TH2F *hAmplSpek = (TH2F*)file->Get("hAmplSpek_cfd"); + NumberOfPixels = hAmplSpek->GetNbinsX(); + channel.resize(hAmplSpek->GetNbinsX()); + +//Amplitude Spectrum of a Pixel + TH1D* hPixelAmplSpek[1440]; ///++++ TODO: Warning hardcoded pixelnumber, change that +++ + +//Book Histograms BookHistos(); +//Define Fit Functions Double_t par[9]; TF1 *g1 = new TF1("g1","gaus",5,15); @@ -94,22 +107,24 @@ TF1 *g4 = new TF1("g4","gaus"); TF1 *total = new TF1("total","gaus(0)+gaus(3)+gaus(6)",7,32); - - + TF1 *fspek_fit = new TF1("fspek_fit", spek_fit, 5, 60, 5); + TF1 *fspek_fit2 = new TF1("fspek_fit2", spek_fit2, 5, 60, 5); + +//Create Canvases TCanvas *cGainFit = new TCanvas("cGainFit", "Fit of Amplitude Spektra", 0,0, 1200,400); cGainFit->Clear(); cGainFit->Divide(3, 2); - TCanvas *cGraphs = new TCanvas("cGraphs", "Horst", 0,401, 1200,400); + TCanvas *cGraphs = new TCanvas("cGraphs", "gain vs. channel", 0,401, 1200,400); cGraphs->Clear(); cGraphs->Divide(1, 2); +//Draw imported Histograms cGainFit->cd(1); gStyle->SetPalette(1,0); gROOT->SetStyle("Plain"); - h->SetTitle("Amplitude Spectrum of all Pixels"); - h->SetAxisRange(0, 100, "Y"); - h->Draw("COLZ"); - channel.resize(h->GetNbinsX()); - NumberOfPixels = h->GetNbinsX(); + hAmplSpek->SetTitle("Amplitude Spectrum of all Pixels"); + hAmplSpek->SetAxisRange(0, 100, "Y"); + hAmplSpek->Draw("COLZ"); + cGainFit->cd(2); @@ -146,7 +161,4 @@ mg->Add(gr2); - - - //----------------------------------------------------------------------------- // Loop over all Pixels and fit Distribution of singles, doubles and tripples @@ -158,16 +170,23 @@ //Title - TString pixelnumber; - pixelnumber += "hPixelProj"; - pixelnumber += "_"; - pixelnumber += pixel + 1 ; + TString histotitle; + histotitle += "hPixelPro_"; + histotitle += pixel + 1 ; cout << "-----------------------------------------------------------" << endl; cout << " PID: " << pixel+1 << endl; cout << "-----------------------------------------------------------" << endl; - - hPixelProjection = h->ProjectionY( pixelnumber , pixel+1, pixel+1); //Projectipon of each Pixel out of Ampl.Spectrum - hPixelProjection->SetTitle("Amplitude Spectrum of one Pixels"); - hPixelProjection->SetAxisRange(0, 40, "X"); + //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"); + + histotitle = "Amplitude Spectrum of Pixel #"; + histotitle += pixel + 1 ; + + hPixelAmplSpek[ pixel ]->SetTitle(histotitle); + hPixelAmplSpek[ pixel ]->SetYTitle("Counts [a.u.]"); + hListPeakSpektra.Add(hPixelAmplSpek[ pixel ]); + //----------------------------------------------------------------------------- @@ -176,20 +195,30 @@ total->SetLineColor(2); + fspek_fit->SetLineColor(3); + fspek_fit2->SetLineColor(kBlue); + fspek_fit2->SetLineStyle(2); + cGainFit->cd(3); gPad->SetLogy(1); - hPixelProjection->SetYTitle("Counts"); - hPixelProjection->Fit(g1,"R"); - hPixelProjection->Fit(g2,"R+"); - hPixelProjection->Fit(g3,"R+"); + + hPixelAmplSpek[ pixel ]->Fit(g1,"R"); + hPixelAmplSpek[ pixel ]->Fit(g2,"R+"); + hPixelAmplSpek[ pixel ]->Fit(g3,"R+"); g1->GetParameters(&par[0]); g2->GetParameters(&par[3]); g3->GetParameters(&par[6]); + + fspek_fit->SetParameters(par); + fspek_fit2->SetParameters(par); total->SetParameters(par); + hPixelAmplSpek[ pixel ]->Fit(total,"R+"); + hPixelAmplSpek[ pixel ]->Fit(fspek_fit,"R+"); + hPixelAmplSpek[ pixel ]->Fit(fspek_fit2,"R+"); + if( showHistos ){ - hPixelProjection->Draw(); - hPixelProjection->Fit(total,"R+"); + hPixelAmplSpek[ pixel ]->Draw(); } - //TODO: Logarithmic Scale + //----------------------------------------------------------------------------- @@ -199,5 +228,8 @@ // 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); @@ -211,35 +243,43 @@ gPad->SetLogy(1); g4->SetLineColor(2); - channel[ pixel ].amplSD = total->GetParameter(4) - total->GetParameter(1); hAmplDistribution->Fill( channel[ pixel ].amplSD ) ; hAmplDistribution->Fit(g4); + distribution.sigma = g4->GetParameter(2); + distribution.mean = g4->GetParameter(1); + distribution.constant = g4->GetParameter(0); cGainFit->cd(5); gPad->SetLogy(1); g4->SetLineColor(2); - channel[ pixel ].amplDT = total->GetParameter(7) - total->GetParameter(4); - hAmplDistribution2->Fill(channel[ pixel ].amplDT) ; + hAmplDistribution2->Fill(channel[ pixel ].amplMean) ; hAmplDistribution2->Fit(g4); cGainFit->cd(6); - hAmplVsChannel->SetBinContent(pixel, channel[ pixel ].amplSD); + hPixelGainFit->SetBinContent(pixel, channel[ pixel ].amplSD); + + if( showHistos ){ + if (pixel%40==0){ + cGainFit->cd(4); + hAmplDistribution->Draw(); + cGainFit->cd(5); + hAmplDistribution2->Draw(); + } + + cGainFit->cd(6); + hPixelGainFit->Draw(); + } + +//----------------------------------------------------------------------------- +// Draw Graphs of Distributions +//----------------------------------------------------------------------------- gr1->SetPoint(pixel, channel[pixel].PID, channel[pixel].amplSD); - gr2->SetPoint(pixel, channel[pixel].PID, channel[pixel].amplDT); + gr2->SetPoint(pixel, channel[pixel].PID, channel[pixel].amplMean); if( showHistos ){ - if (pixel%40==0){ - cGainFit->cd(4); - hAmplDistribution->Draw(); - cGainFit->cd(5); - hAmplDistribution2->Draw(); - } - - cGainFit->cd(6); - hAmplVsChannel->Draw(); if (pixel%60==0){ cGraphs->cd(2); - mg->Draw("APL"); + gr2->Draw("APL"); cGraphs->cd(1); gr1->Draw("APL"); @@ -251,4 +291,8 @@ cGraphs->Update(); } + + //----------------------------------------------------------------------------- + // Debug Modus + //----------------------------------------------------------------------------- if(debug){ @@ -266,7 +310,7 @@ //----------------------------------------------------------------------------- - cGainFit->cd(3); - hPixelProjection->Draw(); - hPixelProjection->Fit(total,"R+"); + +// ( cGainFit->cd(3); +// hPixelAmplSpek[ NumberOfPixels ]->Draw(); cGainFit->cd(4); @@ -277,5 +321,5 @@ cGainFit->cd(6); - hAmplVsChannel->Draw(); + hPixelGainFit->Draw(); cGainFit->Modified(); @@ -283,5 +327,5 @@ cGraphs->cd(2); - mg->Draw("APL"); + gr2->Draw("APL"); cGraphs->cd(1); gr1->Draw("APL"); @@ -320,4 +364,9 @@ out << channel[ pixel ].amplSD << ";\t" ; out << channel[ pixel ].amplDT << ";\t" ; + out << channel[ pixel ].ampl_fspek_fit << "\t" ; + + channel[ pixel ].crazy = SearchAnormal(pixel, 4); //check if pixel had a strange gain + out << channel[ pixel ].crazy << "\t" ; + out << filename << endl; @@ -364,9 +413,10 @@ //Amplification vs Channels - hAmplVsChannel = new TH1F("hAmplVsChannel","Amplifications (SD) of each Channel", 1440, -0.5 , 1439.5 ); - hAmplVsChannel->SetYTitle("Amplification [mV]"); - hAmplVsChannel->SetXTitle("Channels"); - hAmplVsChannel->SetLineColor(2); - hList.Add( hAmplVsChannel ); + 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 ); + } @@ -389,4 +439,68 @@ // MulitOrder MultiGauß Function //----------------------------------------------------------------------------- + +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 = 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]); + 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 + Double_t check = channel[ pixel ].ampl_fspek_fit; + check = TMath::Sqrt(check*check); + if (check < (distance * distribution.sigma + distribution.mean)){ + return false; + } + else{ + cout << "Not normal behaviotruer in Pixel: " << pixel +1 << endl; + return true; + } +} + +/* + + + + + + + */ + /* @@ -410,3 +524,2 @@ } */ -