Changeset 12713

Timestamp:

12/09/11 18:32:03 (13 years ago)

Author:

Jens Buss

Message:

 histos of pixels with to high gain in seperate folder, comments, fit functions improved  

File:

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

fact/tools/rootmacros/gainfit2.C

@@ +1 @@
+//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+//+   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 <TFile.h>
 #include <TH2F.h>
@@ -29 +49 @@
   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 +64 @@
 dist_t distribution;
 
+//Verbosity for Histograms
+TString verbos[3]={"Q","", ""};
 //-----------------------------------------------------------------------------
 // Decleration of Histograms
@@ -49 +73 @@
 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 +88 @@
 void BookHistos();
 void SaveHistograms(const char * );
-bool SearchAnormal(UInt_t, Double_t );
+bool SearchAnormal(UInt_t, Double_t, int verbosityLevel );
 
 //-----------------------------------------------------------------------------
@@ -73 +98 @@
         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 +132 @@
     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 +154 @@
 
 
-    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 +180 @@
     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 +193 @@
 //-----------------------------------------------------------------------------
 
-    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 +226 @@
 //-----------------------------------------------------------------------------
 
-        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 +267 @@
 
 
-//-----------------------------------------------------------------------------
-// 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 +274 @@
 //-----------------------------------------------------------------------------
 
-        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 +285 @@
         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 +301 @@
 
             cGainFit->cd(6);
-            hPixelGainFit->Draw();
+            hGainFitBsl->Draw();
           }
 
@@ -274 +308 @@
 //-----------------------------------------------------------------------------
 
-        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 +345 @@
 
 
-//   ( cGainFit->cd(3);
+//   ( cGainFit->cd(4);
 //    hPixelAmplSpek[ NumberOfPixels ]->Draw();
 
-    cGainFit->cd(4);
+    cGainFit->cd(2);
     hAmplDistribution->Draw();
 
@@ -321 +355 @@
 
     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 +386 @@
                     ++it;
     }
-    cout << filename << endl;
+if (verbosityLevel > 0)    cout << filename << endl;
 
     ofstream out;
@@ -362 +397 @@
     for (UInt_t pixel=0; pixel<NumberOfPixels; pixel++){
             out << channel[ pixel ].PID << ";\t"    ;
-            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 ].ampl_fspek_fit << ";\t" ;
+           // out << channel[ pixel ].amplDT << ";\t" ;
+            out << channel[ pixel ].ampl_fspek_fit2 << "\t" ;
+
+    //check if pixel gain is compareable to gain of the others and save result to logfile
+            channel[ pixel ].crazy = SearchAnormal(pixel, 4, verbosityLevel);
+            if (channel[ pixel ].crazy) hListCrazy.Add(hPixelAmplSpek[ pixel ]);
             out << channel[ pixel ].crazy << "\t" ;
 
@@ -399 +436 @@
 
   //Amplification Distribution - Single2Double
-  hAmplDistribution = new TH1F("hAmplDistribution","Distribution of G-APD Amplification (1+2)", 500, 0 , 20 );
+  hAmplDistribution = new TH1F("hAmplDistribution","Distribution of G-APD Amplification (with Baseline Parameter)", 500, 0 , 20 );
   hAmplDistribution->SetXTitle("Amplification [mV]");
   hAmplDistribution->SetYTitle("Counts");
@@ -406 +443 @@
 
   //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 +449 @@
   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 +469 @@
     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 +507 @@
   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 +530 @@
     }
   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;
-        }
-*/