Changeset 12690 for fact/tools/rootmacros/gainfit2.C

Timestamp:

12/02/11 18:43:49 (13 years ago)

Author:

Jens Buss

Message:

added plot Amplification vs. Channel

File:

• fact/tools/rootmacros/gainfit2.C (modified) (11 diffs)

fact/tools/rootmacros/gainfit2.C

@@ -31 +31 @@
   float amplMean;
   } proj;
-vector<proj> proj_of_pixel;
+vector<proj> 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 +52 @@
 //-----------------------------------------------------------------------------
 
+void BookHistos();
+void SaveHistograms(const char * );
+//Double_t spek_fit(Double_t *x, Double_t *par)
+
 
 //-----------------------------------------------------------------------------
@@ -44 +63 @@
 
 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 +75 @@
 //-----------------------------------------------------------------------------
 
+//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 +93 @@
     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 +110 @@
     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; i<NumberOfPixels; i++){
+       pixelvec[i]=i;
+       simple_gain[i]=0;
+       gain[i]=0;
+    }
+
+    TGraph *gr1 = new TGraph(NumberOfPixels,pixelvec,simple_gain);
+    TGraph *gr2 = new TGraph(NumberOfPixels,pixelvec,gain);
+
+    gr1->SetMarkerStyle(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 +167 @@
         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 +176 @@
 
         total->SetLineColor(2);
-        cGainFit->cd(2);
+        cGainFit->cd(3);
         gPad->SetLogy(1);
         hPixelProjection->SetYTitle("Counts");
@@ -126 +197 @@
 //-----------------------------------------------------------------------------
 
-//        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 +266 @@
 //-----------------------------------------------------------------------------
 
-    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 +316 @@
     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; pixel<NumberOfPixels; pixel++){
+            out << channel[ pixel ].PID << ";\t"    ;
+            out << channel[ pixel ].amplSD << ";\t" ;
+            out << channel[ pixel ].amplDT << ";\t" ;
             out << filename << endl;
 
     }
     out.close();
-
+SaveHistograms( RootOutFileName );
 return( 0 );
 }
+
+//-----------------------------------------------------------------------------
+// End of Main
+//-----------------------------------------------------------------------------
+
+
+
+
+//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+//+                                                                           +
+//+                     Function definitions                                  +
+//+                                                                           +
+//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+
+
+//-----------------------------------------------------------------------------
+// Function: Book Histograms
+//-----------------------------------------------------------------------------
+
+void BookHistos(){
+
+  //Amplification Distribution - Single2Double
+  hAmplDistribution = new TH1F("hAmplDistribution","Distribution of G-APD Amplification (1+2)", 500, 0 , 20 );
+  hAmplDistribution->SetXTitle("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;
+        }
+*/
+