/* ======================================================================== *\!
!   Author(s): Oscar Blanch & Javier Rico
\* ======================================================================== */

void alpha(
    Float_t cutsize=2000,
    Float_t upcutsize=9999999,
    Float_t lengthlowcut=0.,
    Float_t lengthcut=0.26,
    Float_t widthlowcut=0.,
    Float_t widthcut=0.12
)
{    
    // constants
    const Int_t nbins = 9;
    TString psname("alphaMrk.gif");
    // general settings
    gROOT->Reset();
    gStyle->SetCanvasColor(0);
    gStyle->SetCanvasBorderMode(0);
    gStyle->SetPadBorderMode(0);
    gStyle->SetFrameBorderMode(0);
    gStyle->SetOptTitle(0);
    gStyle->SetTitleOffset(1.7,"y");
    gStyle->SetPadLeftMargin(0.15);
    gStyle->SetOptStat(111110);
    gStyle->SetStatColor(0);
    gStyle->SetStatBorderSize(1);
    gStyle->SetStatW(0.2);
    gStyle->SetStatH(0.1);
    gStyle->SetStatX(0.9);
    gStyle->SetStatY(0.9);
    
    // define files, chain, trees, canvas...
    TChain* ton= new TChain("Parameters");
    
    // ON data
    ton->Add("/mnt/users/jrico/magic/mars/mars/mtemp/mifae/hillas/crab20040215OnRotateCalA-D.root");

    TChain* toff= new TChain("Parameters");

    // OFF data
    toff->Add("/mnt/users/jrico/magic/mars/mars/mtemp/mifae/hillas/crab20040215OffRotateCalA-H.root");
    
    TCanvas *c1 = new TCanvas("c1","c1",800,500);
    c1->cd(1);

    // define aliases
    ton->SetAlias("length","MHillas.fLength*0.6/189");
    ton->SetAlias("width","MHillas.fWidth*0.6/189");
    ton->SetAlias("dist","MHillasSrc.fDist*0.6/189");
    ton->SetAlias("conc","MNewImagePar.fConc");
    ton->SetAlias("size","MHillas.fSize");
    ton->SetAlias("event","MRawEvtHeader.fDAQEvtNumber");
    ton->SetAlias("alpha","abs(MHillasSrc.fAlpha)");
    ton->SetAlias("leakage","MNewImagePar.fInnerLeakage1");
    
    toff->SetAlias("length","MHillas.fLength*0.6/189");
    toff->SetAlias("width","MHillas.fWidth*0.6/189");
    toff->SetAlias("dist","MHillasSrc.fDist*0.6/189");
    toff->SetAlias("conc","MNewImagePar.fConc");
    toff->SetAlias("size","MHillas.fSize");
    toff->SetAlias("event","MRawEvtHeader.fDAQEvtNumber");
    toff->SetAlias("alpha","abs(MHillasSrc.fAlpha)");
    toff->SetAlias("leakage","MNewImagePar.fInnerLeakage1");

    // define  cut(s)
//    const char* eventcut="event%4==0 &&";
    const char* eventcut="";
    const char varcut[512];
    //    sprintf(varcut,"size>%f && size<%f && length>%f  &&length<%f && width>%f &&  width<%f",cutsize,upcutsize,lengthlowcut,lengthcut,widthlowcut,widthcut);
    sprintf(varcut,"size>%f && size<%f  && dist>0.2 && dist<1.1 && leakage==0 && length>%f  &&length<%f && width>%f &&  width<%f",cutsize,upcutsize,lengthlowcut,lengthcut,widthlowcut,widthcut);
//    sprintf(varcut,"size>%f && dist>0.2 && dist<0.8 && length<0.21 && width<0.125",cutsize,upcutsize,lengthcut,widthcut);
//    sprintf(varcut,"size>%f && size<%f  && dist>0.2 && dist<0.8 && length<0.21 && width<0.125",cutsize,upcutsize);
//    sprintf(varcut,"size>%f && size<%f  && dist>0.2 && dist<0.8 && length<0.21",cutsize,upcutsize);
//    sprintf(varcut,"size>%f && size<%f  && dist>0.2 && dist<0.8",cutsize,upcutsize);
//    sprintf(varcut,"size>%f && size<%f",cutsize,upcutsize)
//    sprintf(varcut,"");
    Char_t cut[256];
    Char_t cutoff[256];
    strcat(cut,eventcut);
    strcat(cut,varcut);
    strcat(cutoff,varcut);

    // fill on/off histos
    TH1F* onhisto       = new TH1F("OnHist" ,"Alpha Plot",nbins,0,90);
    TH1F* offhisto      = new TH1F("OffHist","Alpha Plot",nbins,0,90);    
    ton->Draw("alpha>>OnHist",cut);
    toff->Draw("alpha>>OffHist",cutoff);
    cout << "Applied cut: " << cut << endl;
    
    // line/normalization
    const Int_t inibin = 20./90.*nbins+1;
    Float_t level=0;
    Float_t leveloff=0;
    Float_t levelofferror=0;
    for(Int_t ibin = inibin; ibin<=nbins;ibin++)
    {
	level+=onhisto->GetBinContent(ibin);
	leveloff+=offhisto->GetBinContent(ibin);
    }
    level/=(nbins-inibin+1);
    leveloff/=(nbins-inibin+1);
    
    // normalize on/off
    offhisto->Sumw2(); // needed to compute correct errors after normalization
    const Float_t norm = level/leveloff;
    cout << "Normalizing by factor " << norm <<endl;
    offhisto->Scale(norm);    

    // significance:
    Float_t sig=0,bg=0,esig=0,ebg=0;
    Float_t significance=0;
    const Int_t signbins = inibin-1;
    for(Int_t ibin = 1; ibin<=signbins;ibin++)
    {
//	Float_t sigma = (onhisto->GetBinContent(ibin)-level)/onhisto->GetBinError(ibin);
//	significance+=sigma*sigma;
	sig  += onhisto->GetBinContent(ibin);
	esig += onhisto->GetBinError(ibin)*onhisto->GetBinError(ibin);
	bg   += offhisto->GetBinContent(ibin);
	ebg  += offhisto->GetBinError(ibin)*offhisto->GetBinError(ibin);
    }
    Float_t error= TMath::Sqrt(esig+ebg);
    significance = (sig-bg)/error;
    
    cout << "Excess: " << sig-bg << endl;
    cout << "N bkg: "  << bg << endl; 
    cout << "Significance: " << significance << endl;

    // plot
    onhisto->SetXTitle("alpha (deg)");
    onhisto->SetYTitle("Entries");
    onhisto->SetMarkerColor(4);
    onhisto->SetMarkerStyle(20);
    onhisto->SetMarkerSize(.7);
    onhisto->SetLineColor(4);
    onhisto->SetMaximum(onhisto->GetBinContent(nbins)*1.5);

    offhisto->SetFillColor(2);
    offhisto->SetMaximum(offhisto->GetBinContent(nbins)*1.5);
    offhisto->SetXTitle("alpha (deg)");
    offhisto->SetYTitle("Entries");

    offhisto->Draw("HIST");
    offhisto->Draw("ESAME");
    onhisto->Draw("ESAMES");
    //    onhisto->Draw("E");

    // move stat box to make them all visible
    gPad->Update();
    TPaveStats* offpavstat = (TPaveStats*) offhisto->GetListOfFunctions()->FindObject("stats");
    if(offpavstat)
    {
	Float_t shiftx = offpavstat->GetX2NDC()-offpavstat->GetX1NDC();
	offpavstat->SetX1NDC(offpavstat->GetX1NDC()-shiftx);
	offpavstat->SetX2NDC(offpavstat->GetX2NDC()-shiftx);
    }

    // draw line
    TLine* line = new TLine(0.,level,90.,level);
    line->SetLineColor(2);
    line->SetLineStyle(2);
//    line->Draw();

    gPad->Modified();
    gPad->Update();

    // produce ps
    //  char psname[50];
    //  sprintf(psname,"%s.ps",var);
    c1->Print(psname);
}