Index: fact/tools/rootmacros/gainfit2.C =================================================================== --- fact/tools/rootmacros/gainfit2.C (revision 12683) +++ fact/tools/rootmacros/gainfit2.C (revision 12690) @@ -31,5 +31,20 @@ float amplMean; } proj; -vector proj_of_pixel; +vector channel; + + +//----------------------------------------------------------------------------- +// Decleration of Histograms +//----------------------------------------------------------------------------- + +//Amplification Distribution +TH1F *hAmplDistribution = NULL; +TH1F *hAmplDistribution2 = NULL; +TH1F *hAmplVsChannel = NULL; + +//Lists of Histograms taht are to be saved +TObjArray hList; +TObjArray hListPeakSpektra; + //----------------------------------------------------------------------------- @@ -37,4 +52,8 @@ //----------------------------------------------------------------------------- +void BookHistos(); +void SaveHistograms(const char * ); +//Double_t spek_fit(Double_t *x, Double_t *par) + //----------------------------------------------------------------------------- @@ -44,7 +63,8 @@ int gainfit2( - const char * InputRootFileName = "../analysis/fpeak/20111109_006_fpeak.root", - const char * OutTextFileName = "../analysis/fpeak/20111109_006_gainfit.txt", - //TODO: Output into txt-file + 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", bool showHistos = false, bool debug = false) @@ -55,8 +75,16 @@ //----------------------------------------------------------------------------- +//Open Rootfiles Files + TFile * baseline = TFile::Open( InputBaselineFileName ); + TFile * file = TFile::Open( InputRootFileName ); + + TH1D *hPixelProjection = NULL; - TFile * file = TFile::Open( InputRootFileName ); +//Baseline + TH2F * hbaseline = (TH2F*)baseline->Get("histo_mean"); +//Amplitude spectrum TH2F *h = (TH2F*)file->Get("hAmplSpek_cfd"); + BookHistos(); Double_t par[9]; @@ -65,17 +93,14 @@ 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)",5,35); - TH1F *hAmplDistribution = new TH1F("hAmplDistribution","Distribution of G-APD Amplification (1+2)", 500, 0 , 20 ); - TH1F *hAmplDistribution2 = new TH1F("hAmplDistribution2","Distribution of G-APD Amplification(2+3)", 500, 0 , 20 ); - TH1F *hAmplDistribution3 = new TH1F("hAmplDistribution3","Distribution of G-APD Amplification(2+3)", 500, 0 , 20 ); - - TCanvas *cGainFit = new TCanvas("cGainFit", "Fit of Amplitude Spektra", 0,0, 800,800); - //TCanvas *cGainDist = new TCanvas("cGainDist", "cGainDist", 0,501, 1000,500); + TF1 *total = new TF1("total","gaus(0)+gaus(3)+gaus(6)",7,32); + + + TCanvas *cGainFit = new TCanvas("cGainFit", "Fit of Amplitude Spektra", 0,0, 1200,400); cGainFit->Clear(); - cGainFit->Divide(2, 2); - -//----------------------------------------------------------------------------- -// Loop over all Pixels and fit Distribution of singles, doubles and tripples -//----------------------------------------------------------------------------- + cGainFit->Divide(3, 2); + + TCanvas *cGraphs = new TCanvas("cGraphs", "Horst", 0,401, 1200,400); + cGraphs->Clear(); + cGraphs->Divide(1, 2); cGainFit->cd(1); @@ -85,10 +110,56 @@ h->SetAxisRange(0, 100, "Y"); h->Draw("COLZ"); - proj_of_pixel.resize(h->GetNbinsX()); + channel.resize(h->GetNbinsX()); NumberOfPixels = h->GetNbinsX(); + + cGainFit->cd(2); + hbaseline->SetTitle("Baselinedistribution of all Pixels"); + hbaseline->SetAxisRange(-5, 5, "X"); + //hbaseline->SetAxisRange(0, 100, "Y"); + hbaseline->Fit("gaus"); + hbaseline->Draw(); + + //----------------------------------------------------------------------------- + // Graphs + //----------------------------------------------------------------------------- + + Double_t simple_gain[NumberOfPixels]; + Double_t gain[NumberOfPixels]; + Double_t pixelvec[NumberOfPixels]; + for (UInt_t i=0; iSetMarkerStyle(5); + gr2->SetMarkerStyle(2); + gr1->SetMarkerColor(kRed); + gr2->SetMarkerColor(kBlack); + gr1->SetMarkerSize(0.5); + gr2->SetMarkerSize(0.5); + TMultiGraph *mg = new TMultiGraph(); + mg->Add(gr1); + mg->Add(gr2); + + + + +//----------------------------------------------------------------------------- +// Loop over all Pixels and fit Distribution of singles, doubles and tripples +//----------------------------------------------------------------------------- + cout << "Number of Pixels: " << NumberOfPixels << endl; - // Begin of Loop for (UInt_t pixel = 0; pixel < NumberOfPixels; pixel++){ + + //Title + TString pixelnumber; + pixelnumber += "hPixelProj"; + pixelnumber += "_"; + pixelnumber += pixel + 1 ; cout << "-----------------------------------------------------------" << endl; @@ -96,5 +167,5 @@ cout << "-----------------------------------------------------------" << endl; - hPixelProjection = h->ProjectionY( "# " , pixel+1, pixel+1); //Projectipon of each Pixel out of Ampl.Spectrum + 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"); @@ -105,5 +176,5 @@ total->SetLineColor(2); - cGainFit->cd(2); + cGainFit->cd(3); gPad->SetLogy(1); hPixelProjection->SetYTitle("Counts"); @@ -126,42 +197,57 @@ //----------------------------------------------------------------------------- -// total->GetParameters(proj_of_pixel[ pixel ].fitParameters); - proj_of_pixel[ pixel ].PID = pixel+1; - - -// proj_of_pixel[ pixel ].tripple_peakSigm = total->GetParameter(8); -// proj_of_pixel[ pixel ].ampl = total->GetParameter(1) - total->GetParameter(4); +// total->GetParameters(channel[ pixel ].fitParameters); + channel[ pixel ].PID = pixel+1; + + +// channel[ pixel ].tripple_peakSigm = total->GetParameter(8); +// channel[ pixel ].ampl = total->GetParameter(1) - total->GetParameter(4); //----------------------------------------------------------------------------- // Draw Histograms //----------------------------------------------------------------------------- - - cGainFit->cd(3); - gPad->SetLogy(1); - g4->SetLineColor(2); - hAmplDistribution->SetXTitle("Amplification [mV]"); - hAmplDistribution->SetYTitle("Counts"); - hAmplDistribution->SetAxisRange(7.5, 12.5, "X"); - proj_of_pixel[ pixel ].amplSD = total->GetParameter(4) - total->GetParameter(1); - hAmplDistribution->Fill( proj_of_pixel[ pixel ].amplSD ) ; - hAmplDistribution->Fit(g4); cGainFit->cd(4); gPad->SetLogy(1); g4->SetLineColor(2); - hAmplDistribution2->SetXTitle("Amplification [mV]"); - hAmplDistribution2->SetYTitle("Counts"); - hAmplDistribution2->SetAxisRange(7.5, 12.5, "X"); - proj_of_pixel[ pixel ].amplDT = total->GetParameter(7) - total->GetParameter(4); - hAmplDistribution2->Fill(proj_of_pixel[ pixel ].amplDT) ; + channel[ pixel ].amplSD = total->GetParameter(4) - total->GetParameter(1); + hAmplDistribution->Fill( channel[ pixel ].amplSD ) ; + hAmplDistribution->Fit(g4); + + cGainFit->cd(5); + gPad->SetLogy(1); + g4->SetLineColor(2); + channel[ pixel ].amplDT = total->GetParameter(7) - total->GetParameter(4); + hAmplDistribution2->Fill(channel[ pixel ].amplDT) ; hAmplDistribution2->Fit(g4); + cGainFit->cd(6); + hAmplVsChannel->SetBinContent(pixel, channel[ pixel ].amplSD); + + gr1->SetPoint(pixel, channel[pixel].PID, channel[pixel].amplSD); + gr2->SetPoint(pixel, channel[pixel].PID, channel[pixel].amplDT); + + if( showHistos ){ - cGainFit->cd(3); + if (pixel%40==0){ + cGainFit->cd(4); hAmplDistribution->Draw(); - cGainFit->cd(4); + cGainFit->cd(5); hAmplDistribution2->Draw(); + } + + cGainFit->cd(6); + hAmplVsChannel->Draw(); + if (pixel%60==0){ + cGraphs->cd(2); + mg->Draw("APL"); + cGraphs->cd(1); + gr1->Draw("APL"); + } + cGainFit->Modified(); cGainFit->Update(); + cGraphs->Modified(); + cGraphs->Update(); } if(debug){ @@ -180,16 +266,26 @@ //----------------------------------------------------------------------------- - cGainFit->cd(2); + cGainFit->cd(3); hPixelProjection->Draw(); hPixelProjection->Fit(total,"R+"); - cGainFit->cd(3); + cGainFit->cd(4); hAmplDistribution->Draw(); - cGainFit->cd(4); + cGainFit->cd(5); hAmplDistribution2->Draw(); + + cGainFit->cd(6); + hAmplVsChannel->Draw(); cGainFit->Modified(); cGainFit->Update(); + + cGraphs->cd(2); + mg->Draw("APL"); + cGraphs->cd(1); + gr1->Draw("APL"); + cGraphs->Modified(); + cGraphs->Update(); //----------------------------------------------------------------------------- @@ -220,13 +316,97 @@ out << "# PixelID : gain between Single and Double Peak : gain between Double and Tripple Peak " << endl; out << "# gain is derived from the parameters of a three gaussian fit function " << endl; - for (int pixel=0; pixel<1440; pixel++){ - out << proj_of_pixel[ pixel ].PID << ";\t" ; - out << proj_of_pixel[ pixel ].amplSD << ";\t" ; - out << proj_of_pixel[ pixel ].amplDT << ";\t" ; + for (UInt_t pixel=0; pixelSetXTitle("Amplification [mV]"); + hAmplDistribution->SetYTitle("Counts"); + hAmplDistribution->SetAxisRange(7.5, 12.5, "X"); + hList.Add( hAmplDistribution ); + + //Amplification Distribution2 - Double2Tripple + hAmplDistribution2 = new TH1F("hAmplDistribution2","Distribution of G-APD Amplification(2+3)", 500, 0 , 20 ); + hAmplDistribution2->SetXTitle("Amplification [mV]"); + hAmplDistribution2->SetYTitle("Counts"); + hAmplDistribution2->SetAxisRange(7.5, 12.5, "X"); + hList.Add( hAmplDistribution2 ); + + //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 ); +} + +//----------------------------------------------------------------------------- +// Save Histograms +//----------------------------------------------------------------------------- + +void SaveHistograms( const char * loc_fname){ + // create the histogram file (replace if already existing) + TFile tf( loc_fname, "RECREATE"); + + hList.Write(); // write the major histograms into the top level directory + tf.mkdir("PeakSpektra"); tf.cd("PeakSpektra"); // go to new subdirectory + hListPeakSpektra.Write(); // write histos into subdirectory + + tf.Close(); // close the file +} + +//----------------------------------------------------------------------------- +// MulitOrder MultiGauß Function +//----------------------------------------------------------------------------- + +/* +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; + } +*/ +