source: trunk/MagicSoft/Mars/mtemp/mberlin/macros/AlphaPlot.C@ 5138

#include "TString.h"
#include "TChain.h"
#include "TCanvas.h"
#include "TH1F.h"
#include "TF1.h"
#include "TLegend.h"
#include "TArrayF.h"
#include "TGraphErrors.h"
#include "MRawRunHeader.h"
#include "TPad.h"
#include "TPaveText.h"
#include "iostream.h"
#include "TObjArray.h"
#include "TFile.h"

//--------------------------------------------------------------------------------------
// MAIN INPUT/OUTPUT

#include "IOMkn421.h"
//--------------------------------------------------------------------------------------

Double_t signal(Double_t *x, Double_t *par)
{
    return par[0]/(TMath::Sqrt(2.*TMath::Pi())*par[1])*TMath::Exp(-(x[0]-par[2])*(x[0]-par[2])/2./par[1]/par[1]);
}

Double_t background(Double_t *x, Double_t *par)
{
    return par[0]+par[1]*x[0]+par[2]*x[0]*x[0]+par[3]*x[0]*x[0]*x[0];
}

Double_t fitFunction(Double_t *x, Double_t *par)
{
    return signal(x,par) + background(x,&par[3]);
}

void AlphaPlot()
{
    DistLo*=deg2mm;
    DistUp*=deg2mm;

    Double_t acut = 10; // alpha region of excess
    Int_t nbins   = 36;
    Int_t nbin_flat= nbins/3;

    // input file chain
    TChain *fEvents=new TChain("Events");
    fEvents->Add(fileOn1.Data());
    fEvents->Add(fileOn2.Data());

    TChain *fEventsOff=new TChain("Events");
    fEventsOff->Add(fileOff1.Data());
    fEventsOff->Add(fileOff2.Data());


    cout<<"Non="<<fEvents->GetEntries()<<"  Noff="<<fEventsOff->GetEntries()<<endl;
    //************************************************************************
    // alpha plot + fit

    // make alpha plot (after cuts)
    TCanvas *c1 = new TCanvas("c1","Alpha distribution after g/h separation",200,10,800,500);
    c1->SetFillColor(0);
    c1->SetBorderMode(0);
    c1->cd();

    TH1F *halpha=new TH1F("halpha","",nbins,0.,90.);
    halpha->GetXaxis()->SetTitle("Alpha [deg]");
    halpha->GetYaxis()->SetTitle("Counts");

    TH1F *halphaOff=new TH1F("halphaOff","",nbins,0.,90.);
    halpha->SetStats(kFALSE);
    halphaOff->SetStats(kFALSE);

    TString strcut=Form("MHillas.fSize>%f",SizeLo);

    strcut+=Form("&& MHillasSrc.fDist>%f",DistLo);
    strcut+=Form("&& MHillasSrc.fDist<%f",DistUp);

    strcut+=Form("&& MHillas.fWidth>%f",WidthLo);
    strcut+=Form("&& MHillas.fWidth<%f",WidthUp);

    strcut+=Form("&& MHillas.fLength>%f",LengthLo);
    strcut+=Form("&& MHillas.fLength<%f",LengthUp);

    // fill histograms
    fEvents->Project("halpha","TMath::Abs(MHillasSrc.fAlpha)",strcut.Data());
    halpha->SetLineWidth(2);
    halpha->DrawCopy("e");

    fEventsOff->Project("halphaOff","TMath::Abs(MHillasSrc.fAlpha)",strcut.Data());
    halphaOff->SetLineColor(2);
    halphaOff->SetLineWidth(2);


    // scale Off histogram to On data using bins nbin_flat(first bin after peak) - nbins
    Double_t nsample_on=0;
    Double_t nsample_off=0;
    for(Int_t i=nbin_flat;i<nbins;i++)
        nsample_on+=halpha->GetBinContent(i);
    for(Int_t i=nbin_flat;i<nbins;i++)
        nsample_off+=halphaOff->GetBinContent(i);

    Double_t scal=nsample_on/nsample_off;
    halphaOff->Sumw2();
    halphaOff->Scale(scal);
    halphaOff->Draw("e same");

    gPad->SetGridx();
    gPad->SetGridy();

    //--------------------------------------------------------------------------
    // fit alpha hist
    TF1 *backFcn = new TF1("backFcn",background,20.,89.5,4);
    backFcn->SetLineWidth(2);
    backFcn->SetLineColor(kRed);

    TF1 *signalFcn = new TF1("signalFcn",signal,0.,89.5,3);
    signalFcn->SetLineWidth(2);
    signalFcn->SetLineColor(kBlue);
    signalFcn->SetNpx(500);

    TF1 *fitFcn = new TF1("fitFcn",fitFunction,0.,89.5,7);
    fitFcn->SetNpx(500);
    fitFcn->SetLineWidth(2);
    fitFcn->SetLineColor(kMagenta);

    // start-parameters
    Double_t p3start=halpha->GetBinContent(halpha->GetNbinsX()-1);
    Double_t p4start=0.;
    Double_t p5start=0.;
    Double_t p6start=0.;

    Double_t p0start=(halpha->GetEntries()-p3start*nbins)*halpha->GetBinWidth(1);
    Double_t p1start=5.;
    Double_t p2start=0.;

    cout<<"Start values for fit:"<<endl;

    cout<<" Gaussian:"<<endl;
    cout<<"  p0 = "<<p0start<<endl;
    cout<<"  p1 = "<<p1start<<endl;
    cout<<"  p2 = "<<p2start<<endl;

    cout<<" Pol3:"<<endl;
    cout<<"  p0 = "<<p3start<<endl;
    cout<<"  p1 = "<<p4start<<endl;
    cout<<"  p2 = "<<p5start<<endl;
    cout<<"  p3 = "<<p6start<<endl<<endl<<endl;

    // background fit
    Double_t parb[4];
    backFcn->SetParameters(p3start,p4start,p5start,p6start);
    backFcn->FixParameter(1,0);   // deriv. at zero equal zero
    //backFcn->FixParameter(2,0); // no 2nd order term
    backFcn->FixParameter(3,0);   // no 3rd order term
    halpha->Fit ("backFcn","RN");
    backFcn->GetParameters(parb);

    // total fit
    Double_t par[7];
    fitFcn->SetParameters(p0start,p1start,p2start,
                          parb[0],parb[1],parb[2],parb[3]);

    fitFcn->FixParameter(2,0.);
    fitFcn->FixParameter(3,parb[0]);
    fitFcn->FixParameter(4,parb[1]);
    fitFcn->FixParameter(5,parb[2]);
    fitFcn->FixParameter(6,parb[3]);
    halpha->Fit ("fitFcn","RN");

    // draw fit results + legend
    TLegend *legend=new TLegend(0.4,0.5,0.93,0.65);
    legend->SetTextFont(40);
    legend->SetTextSize(0.04);
    legend->SetFillColor(19);
    legend->AddEntry(halpha,"Data","lp");
    TString strpol("");
    strpol=Form("Polynom");
    legend->AddEntry(backFcn,strpol.Data(),"l");
    legend->AddEntry(fitFcn,"Global Fit","l");
    //legend->Draw();

    fitFcn->GetParameters(par);
    fitFcn->SetRange(-90.,90.);
    fitFcn->SetLineWidth(2);
    fitFcn->Draw("same");

    signalFcn->SetParameters(par);
    //signalFcn->SetLineWidth(2);
    //signalFcn->Draw("same");

    backFcn->SetParameters(&par[3]);
    backFcn->SetRange(-90,90);
    backFcn->SetLineWidth(2);
    backFcn->Draw("same");

    //-----------------------------------------------------------------------
    // calculate significance

    // analytical method
    Double_t nSig=signalFcn->Integral(0.,acut)/halpha->GetBinWidth(1);
    Double_t nOffScaled=backFcn->Integral(0.,acut)/halpha->GetBinWidth(1);

    Double_t nOn=nSig+nOffScaled;

    // significance after Li/Ma
    Double_t theta=1.; // should be = Ton/Toff ???
    Double_t nOff=nOffScaled/theta;
    Double_t siglima = sqrt(2*(nOn*TMath::Log((1+theta)*nOn/(theta*(nOn+nOff)))
                               +nOff*TMath::Log((1+theta)*(nOff/(nOn+nOff)))));

    // "standard" significance
    Double_t sigstd  = nSig/sqrt(nOn+nOffScaled);

    TPaveText *pt = new TPaveText(0.25,0.75,0.56,0.98,"brNDC");
    pt->SetTextFont(60);
    pt->SetTextSize(0.03);
    pt->SetTextAlign(12);

    TString strcut1(Form("cuts:"));
    TString strcut3(Form(" |alpha| <  %4.2f ",acut));
    pt->AddText(strcut3.Data());
    TString str("");
    pt->AddText(str);

    TString strsig(Form("NExcess = %d",Int_t(nSig)));
    pt->AddText(strsig.Data());
    TString strbak(Form("NOff    = %d",Int_t(nOffScaled)));
    pt->AddText(strbak.Data());
    pt->AddText(str);

    TString strsg(Form("Significance (Li/Ma) = %4.2f",siglima));
    pt->AddText(strsg.Data());
    strsg=Form("Significance (std) = %4.2f",sigstd);
    pt->AddText(strsg.Data());
    pt->Draw();

    cout<<endl<<"Results for events with |alpha| < "<<acut<<" degrees :"<<endl;
    cout<<" noff         ="<<nOffScaled<<endl;
    cout<<" nexcess      ="<<nSig<<endl;
    cout<<" sig. (Li/Ma) ="<<siglima<<endl;
    cout<<" sig.-std     ="<<sigstd<<endl<<endl<<endl;

    gPad->SetGridx();
    gPad->SetGridy();

    c1->Update();

    return;
}