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


void Einfangzeiten(double xmin = 300., double xmax = 20000., 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.);

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

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

    TH1D *hEinfang = new TH1D(*ho);
    hEinfang->Add(hu, 1.);
    hEinfang->SetTitle("Einfangzeit");
    hEinfang->SetName("hE");

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

    TF1 *f1 = new TF1("f1", fitFunc, 0., 20000., 5);    // (Name, Titel, xmin, xmax, Anzahl der Parameter)
    f1->SetParameters(2000., 1.3, 5000., 880.,30.);       // Startwerte fuer Parameter festlegen
    f1->SetParName(0, "#tau_{0}");                      // Parameter benennen
    f1->SetParName(1, "f");
    f1->SetParName(2, "# #mu^{+}");
    f1->SetParName(3, "#tau_{c}");
    f1->SetParName(4, "BG");
    f1->SetNpx(200);                                    // Anzahl der Punkte auf denen geschrieben wird

    // Hier können die Bereiche, in denen die Parameter beim variiert werden
    // angegeben werden. 
    f1->SetParLimits(0, 0., 10000.);  // ---> Bereich f\"ur tau_0: 0ns ... 10000ns 
    f1->SetParLimits(1, 0., 10000.);   // ---> Bereich f\"ur f: 0ns ... 10000ns                  
    f1->SetParLimits(3, 0., 10000.); // ---> Bereich f\"ur tau_c: 0ns ... 10000ns

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

    c0->Print("Einfangzeiten.ps(");

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

    c1->Print("ZeitspektrumGesamt.ps");

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

    c2->Print("ZeitspektrumGesamt.ps");

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

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

    c3->Print("Einfangzeiten.ps");

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

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

    c4->Print("Einfangzeiten.ps");

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

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

    c5->Print("Einfangzeiten.ps)");
}

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

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

    return value;
}