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

Timestamp:

12/08/11 00:28:06 (13 years ago)

Author:

Jens Buss

Message:

added differnt typee of fit function to see which fits best; detect pixel with not normal gain

File:

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

fact/tools/rootmacros/gainfit2.C

@@ -27 +27 @@
   UInt_t PID;
   vector<float> 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<proj> channel;
 
+typedef struct{
+  Double_t constant;
+  Double_t mean;
+  Double_t sigma;
+} dist_t;
+dist_t distribution;
 
 //-----------------------------------------------------------------------------
@@ -41 +49 @@
 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 +87 @@
     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 +107 @@
     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 +161 @@
     mg->Add(gr2);
 
-
-
-
 //-----------------------------------------------------------------------------
 // Loop over all Pixels and fit Distribution of singles, doubles and tripples
@@ -158 +170 @@
 
         //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 +195 @@
 
         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 +228 @@
 //        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 +243 @@
         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 +291 @@
         cGraphs->Update();
     }
+
+    //-----------------------------------------------------------------------------
+    // Debug Modus
+    //-----------------------------------------------------------------------------
         if(debug){
 
@@ -266 +310 @@
 //-----------------------------------------------------------------------------
 
-    cGainFit->cd(3);
-    hPixelProjection->Draw();
-    hPixelProjection->Fit(total,"R+");
+
+//   ( cGainFit->cd(3);
+//    hPixelAmplSpek[ NumberOfPixels ]->Draw();
 
     cGainFit->cd(4);
@@ -277 +321 @@
 
     cGainFit->cd(6);
-    hAmplVsChannel->Draw();
+    hPixelGainFit->Draw();
 
     cGainFit->Modified();
@@ -283 +327 @@
 
     cGraphs->cd(2);
-    mg->Draw("APL");
+    gr2->Draw("APL");
     cGraphs->cd(1);
     gr1->Draw("APL");
@@ -320 +364 @@
             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 +413 @@
 
   //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 +439 @@
 // 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 +524 @@
         }
 */
-