// Usage:
// .L ../Lebensdauer.C
// Lebensdauer(xmin, xmax, "../fp13.root")


void Lebensdauer(double xmin = 300., double xmax = 10000., const char *filename = "fp13.root") {

//---------------------------------------------------------------------//
// Laden der Histogramme                                               //
//---------------------------------------------------------------------//

    gStyle->SetOptFit(111);
    gStyle->SetOptStat("ni");
    gStyle->SetStatFormat("g");
    gStyle->SetMarkerStyle(20);
    gStyle->SetMarkerSize(0.5);
    gStyle->SetHistLineWidth(1);
    gROOT->ForceStyle();

    TFile *f = new TFile(filename);

    TH1D *a1 = (TH1D*)f->Get("a1");                         
    TH1D *a2 = (TH1D*)f->Get("a2");                              
    TH1D *a3 = (TH1D*)f->Get("a3");                              
    TH1D *a4 = (TH1D*)f->Get("a4");                             
    TH1D *a5 = (TH1D*)f->Get("a5");
    TH1D *a6 = (TH1D*)f->Get("a6");

    TH1D *b1 = (TH1D*)f->Get("b1");                         
    TH1D *b2 = (TH1D*)f->Get("b2");                              
    TH1D *b3 = (TH1D*)f->Get("b3");                              
    TH1D *b4 = (TH1D*)f->Get("b4");                             
    TH1D *b5 = (TH1D*)f->Get("b5");  
    TH1D *b6 = (TH1D*)f->Get("b6");  

    TH1D *x1 = (TH1D*)f->Get("x1");                              
    TH1D *x2 = (TH1D*)f->Get("x2");                              
    TH1D *x3 = (TH1D*)f->Get("x3");                              
    TH1D *x4 = (TH1D*)f->Get("x4");                              
    TH1D *x5 = (TH1D*)f->Get("x5");                              
    TH1D *x6 = (TH1D*)f->Get("x6");  

//---------------------------------------------------------------------//
// Bereinigung der Histogramme von Nachpulse                           //
//---------------------------------------------------------------------//
    
    a1->Sumw2();              // speichert die Fehler als Wurzel des Eintrages
    a2->Sumw2();              // um diese spaeter richtig zu berechnen
    a3->Sumw2();
    a4->Sumw2();
    a6->Sumw2();
    b1->Sumw2();
    b2->Sumw2();
    b3->Sumw2();
    b4->Sumw2();
    b5->Sumw2();
    x1->Sumw2();
    x2->Sumw2();
    x3->Sumw2();
    x4->Sumw2();
    x5->Sumw2();
    x6->Sumw2();

    a1->Add(x1,0.);
    a2->Add(x2,0.);
    a3->Add(x3,0.);
    a4->Add(x4,0.);
    a5->Add(x5,0.);
    a6->Add(x6,0.);

    b1->Add(x2,0.);
    b2->Add(x3,0.);
    b3->Add(x4,0.);
    b4->Add(x5,0.);
    b5->Add(x6,0.);

//---------------------------------------------------------------------//  
// Laden und setzen der Fitfunktion                                    //
//---------------------------------------------------------------------//

    TF1 *f1 = new TF1("f1", fitFunc, 0., 20000., 3);       // ( Name, Titel, xmin, xmax, Anzahl der Parameter)
    f1->SetParameters(2000., 5000., 10.);                   // Startwerte fuer Parameter festlegen
    f1->SetParName(0, "#tau_{0}");                         // Parameter benennen
    f1->SetParName(1, "# #mu");
    f1->SetParName(2, "BG");
    f1->SetNpx(200);                                       // Anzahl der Punkte auf denen geschrieben wird
    f1->SetParLimits(1, 0., 100000.);
    f1->SetParLimits(2, 0.2, 100000.);

//---------------------------------------------------------------------//
// Zusammenfuegen der Histogramme                                      //
//---------------------------------------------------------------------//

    TH1D *ho = new TH1D(*a1);
    ho->Add(a2, 1.);
    ho->Add(a3, 1.);
    ho->Add(a4, 1.);
    ho->Add(a5, 1.);
    ho->Add(a6, 1.);
    ho->SetTitle("Zerfall nach oben");
    ho->SetName("ho");

    printf("\n\n   ---> Fit Lebensdauer: Zerfaelle nach oben \n\n\n ");
    ho->Fit(f1, "", "e", xmin, xmax);

    TH1D *hu = new TH1D(*b1);
    hu->Add(b2, 1.);
    hu->Add(b3, 1.);
    hu->Add(b4, 1.);
    hu->Add(b5, 1.);
    hu->Add(b6, 1.);
    hu->SetTitle("Zerfall nach unten");
    hu->SetName("hu");

    printf("\n\n   ---> Fit Lebensdauer: Zerfaelle nach oben \n\n\n ");
    hu->Fit(f1, "", "e", xmin, xmax);

    TH1D *hLebensdauer = new TH1D(*ho);
    hLebensdauer->Add(hu, 1.);
    hLebensdauer->SetTitle("Lebensdauer");
    hLebensdauer->SetName("hL");

//---------------------------------------------------------------------//
// Histogramme zeichnen                                                //
//---------------------------------------------------------------------//

    TCanvas *c0 = new TCanvas("c0");
    c0->Clear();
    c0->Divide(3,4);
    gPad->SetLogy();

    c0->cd(1);
    gPad->SetLogy();
    a1->Draw("e");

    c0->cd(2);
    gPad->SetLogy();
    a2->Draw("e");

    c0->cd(3);
    gPad->SetLogy();
    a3->Draw("e");

    c0->cd(4);
    gPad->SetLogy();
    a4->Draw("e");

    c0->cd(5);
    gPad->SetLogy();
    a5->Draw("e");

    c0->cd(6);
    gPad->SetLogy();
    a6->Draw("e");

    c0->cd(7);
    gPad->SetLogy();
    b1->Draw("e");

    c0->cd(8);
    gPad->SetLogy();
    b2->Draw("e");

    c0->cd(9);
    gPad->SetLogy();
    b3->Draw("e");

    c0->cd(10);
    gPad->SetLogy();
    b4->Draw("e");

    c0->cd(11);
    gPad->SetLogy();
    b5->Draw("e");

    c0->cd(12);
    gPad->SetLogy();
    b6->Draw("e");



    TCanvas *c1 = new TCanvas("c1");
    c1->Clear();
    c1->Divide(2,3);
    gPad->SetLogy();

    c1->cd(1);
    gPad->SetLogy();
    x1->Draw("e");

    c1->cd(2);
    gPad->SetLogy();
    x2->Draw("e");

    c1->cd(3);
    gPad->SetLogy();
    x3->Draw("e");

    c1->cd(4);
    gPad->SetLogy();
    x4->Draw("e");

    c1->cd(5);
    gPad->SetLogy();
    x5->Draw("e");

    c1->cd(6);
    gPad->SetLogy();
    x6->Draw("e");


    TCanvas *c2 = new TCanvas("c2");
    c2->Clear();
    c2->Divide(1,2);
    gPad->SetLogy();

    c2->cd(1);
    gPad->SetLogy();
    ho->Draw("e");

    c2->cd(2);
    gPad->SetLogy();
    hu->Draw("e");



    TCanvas *c3 = new TCanvas("c3");
    c3->Clear();
    gPad->SetLogy();

    printf("\n\n   ---> Fit Lebensdauer: Alle Zerfaelle \n\n\n ");
    hLebensdauer->Fit(f1, "", "e", xmin, xmax);
}

//---------------------------------------------------------------------//
// definieren der Fitfunktion                                          //
//---------------------------------------------------------------------//

double fitFunc(double *x, double *par) {
  double value =  par[1]*(TMath::Exp(-x[0]/par[0])) + par[2];
  //double value =  par[1]*(TMath::Exp(-x[0]/par[0]));
    return value;
}