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


void Asymmetrie2(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("fp13_ohneBFeld.root");

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

      TFile *f = new TFile("fp13_mitBFeld.root");

    TH1D *am1 = (TH1D*)f->Get("a1");                         
    TH1D *am2 = (TH1D*)f->Get("a2");                              
    TH1D *am3 = (TH1D*)f->Get("a3");                              
    TH1D *am4 = (TH1D*)f->Get("a4");                             
    TH1D *am5 = (TH1D*)f->Get("a5");
    TH1D *am6 = (TH1D*)f->Get("a6");

    TH1D *bm1 = (TH1D*)f->Get("b1");                         
    TH1D *bm2 = (TH1D*)f->Get("b2");                              
    TH1D *bm3 = (TH1D*)f->Get("b3");                              
    TH1D *bm4 = (TH1D*)f->Get("b4");                             
    TH1D *bm5 = (TH1D*)f->Get("b5");  
    TH1D *bm6 = (TH1D*)f->Get("b6");  

    TH1D *xm1 = (TH1D*)f->Get("x1");                              
    TH1D *xm2 = (TH1D*)f->Get("x2");                              
    TH1D *xm3 = (TH1D*)f->Get("x3");                              
    TH1D *xm4 = (TH1D*)f->Get("x4");                              
    TH1D *xm5 = (TH1D*)f->Get("x5");                              
    TH1D *xm6 = (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();
    a5->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();

    am1->Sumw2();              // speichert die Fehler als Wurzel des Eintrages
    am2->Sumw2();              // um diese spaeter richtig zu berechnen
    am3->Sumw2();
    am4->Sumw2();
    am5->Sumw2();
    am6->Sumw2();
    bm1->Sumw2();
    bm2->Sumw2();
    bm3->Sumw2();
    bm4->Sumw2();
    bm5->Sumw2();
    xm1->Sumw2();
    xm2->Sumw2();
    xm3->Sumw2();
    xm4->Sumw2();
    xm5->Sumw2();
    xm6->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.);

    am1->Add(xm1,0.);
    am2->Add(xm2,0.);
    am3->Add(xm3,0.);
    am4->Add(xm4,0.);
    am5->Add(xm5,0.);
    am6->Add(xm6,0.);
    
    bm1->Add(xm2,.0);
    bm2->Add(xm3,0.);
    bm3->Add(xm4,0.);
    bm4->Add(xm5,0.);
    bm5->Add(xm6,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 ohne B-Feld");
    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 ohne B-Feld");
    hu->SetName("hu");

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

    TH1D *hmo = new TH1D(*am1);
    hmo->Add(am2, 1.);
    hmo->Add(am3, 1.);
    hmo->Add(am4, 1.);
    hmo->Add(am5, 1.);
    hmo->Add(am6, 1.);
    hmo->SetTitle("Zerfall nach oben mit B-Feld");
    hmo->SetName("hmo");

    TH1D *hmu = new TH1D(*bm1);
    hmu->Add(bm2, 1.);
    hmu->Add(bm3, 1.);
    hmu->Add(bm4, 1.);
    hmu->Add(bm5, 1.);
    hmu->Add(bm6, 1.);
    hmu->SetTitle("Zerfall nach unten mit B-Feld");
    hmu->SetName("hmu");

    TH1D *hEinfangm = new TH1D(*hmo);
    hEinfangm->Add(hmu, 1.);
    hEinfangm->SetTitle("Lebensdauer mit B-Feld");
    hEinfangm->SetName("hmE");


    TH1D *totsumu = new TH1D(*hu);
    totsumu->Add(hmu, 18999./20256.);
    totsumu->SetName("totsumu");
    totsumu->SetTitle("totsumu");

    TH1D *totdiffu = new TH1D(*hu);
    totdiffu->Add(hmu, -18999./20256.);
    totdiffu->SetName("totdiffu");
    totdiffu->SetTitle("totdiffu");

   
    totdiffu->Divide(totsumu);

    TH1D *totsumo = new TH1D(*ho);
    totsumo->Add(hmo, 26580./25379.);
    totsumo->SetName("totsumo");
    totsumo->SetTitle("totsumuo");

    TH1D *totdiffo = new TH1D(*ho);
    totdiffo->Add(hmo, -26580./25379.);
    totdiffo->SetName("totdiffo");
    totdiffo->SetTitle("totdiffo");
    
   
    totdiffo->Divide(totsumo);


    TH1D *hAsymo = new TH1D(*totdiffo);
    hAsymo->SetName("hAsymo");
    hAsymo->SetTitle("Asymmetrie oben");
    hAsymo->SetMinimum(-0.4);
    hAsymo->SetMaximum(0.4);

    TH1D *hAsymu = new TH1D(*totdiffu);
    hAsymu->SetName("hAsymu");
    hAsymu->SetTitle("Asymmetrie unten");
    hAsymu->SetMinimum(-0.4);
    hAsymu->SetMaximum(0.4);

    TH1D *hAsym = new TH1D(*hAsymu);
    hAsym->Add(hAsymo, -1.);
    hAsym->Scale(0.5);
    hAsymu->SetName("hAsym");
    hAsym->SetTitle("Asymmetrie gesamt");
    hAsym->SetMinimum(-0.4);
    hAsym->SetMaximum(0.4);


//---------------------------------------------------------------------//
// 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 *c2 = new TCanvas("c3");
    c3->Clear();
    c3->Divide(1,2);
    gPad->SetLogy();

    c3->cd(1);
    gPad->SetLogy();
    hmo->Draw("e");

    c3->cd(2);
    gPad->SetLogy();
    hmu->Draw("e");



//---------------------------------------------------------------------//  
// Laden und setzen der Fitfunktion Asymmetrie unten                   //
//---------------------------------------------------------------------//

    TF1 *f1 = new TF1("f1", fitFuncAsym, 0., 10000., 4);  // ( Name, Titel, xmin, xmax, Anzahl der Parameter) 
    f1->SetParameters(0.08, 4.e-3, 0.0, 0.0);           // Startwerte fuer Parameter festlegen
    f1->SetParName(0, "P*A");                             // Parameter benennen
    f1->SetParName(1, "#omega");
    f1->SetParName(2, "#varphi");
    f1->SetParName(3, "const");
    //    f1->SetParLimits(0, 0., 1.);

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

    hAsymu->Fit(f1, "", "e", xmin, xmax);
    TLine s;
    s.DrawLine(0., 0., 10000., 0.);

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


//---------------------------------------------------------------------//  
// Laden und setzen der Fitfunktion Asymmetrie oben                    //
//---------------------------------------------------------------------//

    TF1 *f1 = new TF1("f1", fitFuncAsym, 0., 10000., 4);  // ( Name, Titel, xmin, xmax, Anzahl der Parameter) 
    f1->SetParameters(0.08, 4e-3, 1.14, 0.0);           // Startwerte fuer Parameter festlegen
    f1->SetParName(0, "P*A");                             // Parameter benennen
    f1->SetParName(1, "#omega");
    f1->SetParName(2, "#varphi");
    f1->SetParName(3, "const");
//    f1->SetParLimits(0, 0., 1.);

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

    hAsymo->Fit(f1, "", "e", xmin, xmax);
    TLine s;
    s.DrawLine(0., 0., 10000., 0.);

    c5->Print("Asymmetrie.ps");

//---------------------------------------------------------------------//  
// Laden und setzen der Fitfunktion Asymmetrie gesamt                    //
//---------------------------------------------------------------------//

    TF1 *f1 = new TF1("f1", fitFuncAsym, 0., 10000., 4);  // ( Name, Titel, xmin, xmax, Anzahl der Parameter) 
    f1->SetParameters(0.08, 4e-3, 1.14, 0.0);           // Startwerte fuer Parameter festlegen
    f1->SetParName(0, "P*A");                             // Parameter benennen
    f1->SetParName(1, "#omega");
    f1->SetParName(2, "#varphi");
    f1->SetParName(3, "const");
//    f1->SetParLimits(0, 0., 1.);

    TCanvas *c6 = new TCanvas("c6");
    c6->Clear();

    hAsym->Fit(f1, "", "e", xmin, xmax);
    TLine s;
    s.DrawLine(0., 0., 10000., 0.);

    c6->Print("Asymmetrie.ps)");

//---------------------------------------------------------------------//
// Schreiben der Histogramme                                           //
//---------------------------------------------------------------------//
    TFile Open("Results.root","RECREATE");
    ho->Write();
    hu->Write();
    hmo->Write();
    hmu->Write();
    hEinfang->Write();
    hAsymo->Write();
    hAsymu->Write();
    hAsym->Write();
}
//---------------------------------------------------------------------//
// definieren der Fitfunktion                                          //
//---------------------------------------------------------------------//

double fitFuncAsym(double *x, double *par) {
    double value = 0.5 * par[0] * TMath::Cos(par[1] * x[0] + par[2]) + par[3];
    return value;
}