source: fact/tools/rootmacros/gainfit2.C@ 12677

#include <TFile.h>
#include <TH2F.h>
#include <TH1D.h>
#include <TROOT.h>
#include <TF1.h>
#include <TTimer.h>
#include <TString.h>
#include <Getline.h>
#include <TCanvas.h>
#include <TGraph.h>
#include <TMath.h>
#include <TMultiGraph.h>
#include <TH1.h>
#include <TProfile.h>
#include <TStyle.h>

#include <iostream>
#include <fstream>
using namespace std;

//-----------------------------------------------------------------------------
// Decleration of variables
//-----------------------------------------------------------------------------
UInt_t NumberOfPixels;

  typedef struct{
  UInt_t PID;
  vector<float> fitParameters[9];
  float amplSD;
  float amplDT;
  float amplMean;
  } proj;
vector<proj> proj_of_pixel;

//-----------------------------------------------------------------------------
// Functions
//-----------------------------------------------------------------------------


//-----------------------------------------------------------------------------
// Main function
//-----------------------------------------------------------------------------


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
        bool showHistos = false,
        bool debug = false)
{

//-----------------------------------------------------------------------------
// Histograms, Canvases and Fit Functions
//-----------------------------------------------------------------------------

    TH1D *hPixelProjection = NULL;
    TFile * file = TFile::Open( InputRootFileName );
    TH2F *h = (TH2F*)file->Get("hAmplSpek_cfd");
    

    Double_t par[9];
    TF1 *g1    = new TF1("g1","gaus",5,15);
    TF1 *g2    = new TF1("g2","gaus",15,25);
    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);
    cGainFit->Clear();
    cGainFit->Divide(2, 2);

//-----------------------------------------------------------------------------
// Loop over all Pixels and fit Distribution of singles, doubles and tripples
//-----------------------------------------------------------------------------

    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");
    proj_of_pixel.resize(h->GetNbinsX());
    NumberOfPixels = h->GetNbinsX();
    cout << "Number of Pixels: " << NumberOfPixels << endl;

    // Begin of Loop
    for (UInt_t pixel = 0; pixel < NumberOfPixels; pixel++){

        cout << "-----------------------------------------------------------" << endl;
        cout << " PID: " << pixel+1 << endl;
        cout << "-----------------------------------------------------------" << endl;
        
        hPixelProjection = h->ProjectionY( "# " , pixel+1, pixel+1);     //Projectipon of each Pixel out of Ampl.Spectrum
        hPixelProjection->SetTitle("Amplitude Spectrum of one Pixels");
        hPixelProjection->SetAxisRange(0, 40, "X");

//-----------------------------------------------------------------------------
// Single Gaussian Fits and Tripple-Gaussian Fit Funciton
//-----------------------------------------------------------------------------

        total->SetLineColor(2);
        cGainFit->cd(2);
        gPad->SetLogy(1);
        hPixelProjection->SetYTitle("Counts");
        hPixelProjection->Fit(g1,"R");
        hPixelProjection->Fit(g2,"R+");
        hPixelProjection->Fit(g3,"R+");
        g1->GetParameters(&par[0]);
        g2->GetParameters(&par[3]);
        g3->GetParameters(&par[6]);
        total->SetParameters(par);
        if( showHistos ){
        hPixelProjection->Draw();
        hPixelProjection->Fit(total,"R+");
        }

        //TODO: Logarithmic Scale

//-----------------------------------------------------------------------------
// Save Parameters into vetor
//-----------------------------------------------------------------------------

//        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);

//-----------------------------------------------------------------------------
// 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) ;
        hAmplDistribution2->Fit(g4);

    if( showHistos ){
        cGainFit->cd(3);
        hAmplDistribution->Draw();
        cGainFit->cd(4);
        hAmplDistribution2->Draw();
        cGainFit->Modified();
        cGainFit->Update();
    }
        if(debug){

            TTimer timer("gSystem->ProcessEvents();", 50, kFALSE);
            timer.TurnOn();
            TString input = Getline("Type 'q' to exit, <return> to go on: ");
            timer.TurnOff();
            if (input=="q\n") {
              break;
            }
        }
    }
//-----------------------------------------------------------------------------
// End of Loop
//-----------------------------------------------------------------------------

    cGainFit->cd(2);
    hPixelProjection->Draw();
    hPixelProjection->Fit(total,"R+");

    cGainFit->cd(3);
    hAmplDistribution->Draw();

    cGainFit->cd(4);
    hAmplDistribution2->Draw();

    cGainFit->Modified();
    cGainFit->Update();

//-----------------------------------------------------------------------------
// Write into logfile
//-----------------------------------------------------------------------------
    string filename =  InputRootFileName;
    size_t lastslash = filename.find_last_of('/');
    filename = filename.substr(lastslash+1);
    size_t firstdot = filename.find_first_of('.');
    filename = filename.substr(0, firstdot);

    string::iterator it;
    for (it=filename.begin(); it < filename.end(); ){

            if (!isdigit(*it)){
                    it=filename.erase(it);
            }
            else
                    ++it;
    }
    cout << filename << endl;

    ofstream out;
    out.open( OutTextFileName );
    out << "pixelID/I:gainSD/F:gainDT:filename/l" << endl;
    out << "# ---- the 1st line helps to use TTree::ReadFile, for reading again this txt file with root. " << endl;
    out << "#" << InputRootFileName << endl;
    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" ;
            out << filename << endl;

    }
    out.close();

return( 0 );
}