/////////////////////////////////////////////////////////////////////////////
//                                                                         //
// analysis_read.C                                                         //
//                                                                         //
// Author(s): S.C. Commichau, 3/2004                                       //
//                                                                         //
/////////////////////////////////////////////////////////////////////////////


// This macro reads the output of analysis.C and plots some of the
// image parameters. It also produces an ALPHA plot and calculates
// the significance

void analysis_read(){

    gROOT->Reset();
    
    // Some parameters for the ALPHA-plot
    const Int_t inibin = 9;
    const Int_t nbins  = 18;

    // ON data chain
    TChain* CON  = new TChain("Parameters");

    // OFF data chain
    TChain* COFF = new TChain("Parameters");

    // Add the root file(s) to the chain and treat it like a tree
    CON->Add("~/Mars/on1.root");
    COFF->Add("~/Mars/off1.root");

    //CON->Add("~/no_back/hillas_on_all23.root");
    //COFF->Add("~/no_back/hillas_off_all23.root");

    TString title = "Mrk421 - 04/23/2004";

    // Set some global style options
    //gStyle->SetOptDate(10);
    //gStyle->SetDateX(.1);
    //gStyle->SetDateY(.905);
    gStyle->SetOptFit(0);
    gStyle->SetOptStat(110011);
    gStyle->SetFrameBorderMode(0);
    //gStyle->SetPalette(1);
    gStyle->SetPalette(1,0);
    gStyle->SetFrameBorderSize(0);
    gStyle->SetCanvasColor(0);
    gStyle->SetFrameFillColor(0);
    gStyle->SetTitleFillColor(0);
    gStyle->SetTitleBorderSize(0);
    gStyle->SetStatColor(0);
    gStyle->SetStatBorderSize(1);
    // Put the statistics box into the right corner!
    gStyle->SetStatX(.9);
    gStyle->SetStatY(.9);
    //gStyle->SetStatTextColor(2);

    // Set alias for each quantity to simplify the handling
    CON->SetAlias("length","MHillas.fLength*0.6/189.0");
    CON->SetAlias("width","MHillas.fWidth*0.6/189.0");
    CON->SetAlias("meanx","MHillas.fMeanX*0.6/189.0");
    CON->SetAlias("meany","MHillas.fMeanY*0.6/189.0");
    CON->SetAlias("size","MHillas.fSize");
    CON->SetAlias("dist","MHillasSrc.fDist*0.6/189.0");
    CON->SetAlias("alpha","MHillasSrc.fAlpha");
    CON->SetAlias("delta","abs(MHillas.fDelta)*180/3.1415");
    //CON->SetAlias("delta","MDCA.fDelta1*180/3.1415");
    //CON->SetAlias("dca","MDCA.fDCA*0.6/189.0");

    COFF->SetAlias("length","MHillas.fLength*0.6/189.0");
    COFF->SetAlias("width","MHillas.fWidth*0.6/189.0");
    COFF->SetAlias("meanx","MHillas.fMeanX*0.6/189.0");
    COFF->SetAlias("meany","MHillas.fMeanY*0.6/189.0");
    COFF->SetAlias("size","MHillas.fSize");
    COFF->SetAlias("dist","MHillasSrc.fDist*0.6/189.0");
    COFF->SetAlias("alpha","MHillasSrc.fAlpha");
    COFF->SetAlias("delta","abs(MHillas.fDelta)*180/3.1415");


    // Cuts...
    // Cut for ALPHA
    TString cut1 = "(dist>0.2) && (dist<.7) && (width>0.04) && (width<0.14) && (length>0.14) && (length<0.26) && (size>500)";
        
    TString cut2 = "";//size>1000 && size<2000";

    TString cut3 = "";

    TString cut4 = "";

    TString cut5 = "";
    

    // And now the draw area
    
    //************ Width vs Length ************

    TCanvas* CWL = new TCanvas("CWL","Canvas",0,0,300,300);
    CWL->cd();
    CWL->SetBorderMode(0); 
    gPad->SetBorderMode(0);
    gStyle->SetOptStat(0);
            
    TH2F* hWL = new TH2F("hWL","",200,0,0.7,200,0,1);
    
    // Leaf, Cuts, Draw option, use >>+ to append data to an ex. histogram
    CON->Draw("length:width>>hWL",cut2);

    hWL->SetTitle(title); 
    hWL->GetXaxis()->SetTitle("WIDTH [#circ]");
    hWL->GetYaxis()->SetTitle("LENGTH [#circ]");
    hWL->GetYaxis()->SetTitleOffset(1.2);
    hWL->SetMarkerColor(2);
     
    hWL->SetMaximum(150);
    hWL->Draw("COLZ");//CONT


    //************ Width and Length ************

    // Create Canvas and set Canvas some options
    TCanvas* Cwl = new TCanvas("Cwl","Canvas",325,0,300,300);
    Cwl->cd();
    Cwl->SetBorderMode(0);   
    gPad->SetBorderMode(0);
    gStyle->SetOptStat(110011);

    TH1F* hLength = new TH1F("hLength","",100,0,1.1);
        
    CON->Draw("length>>hLength",cut5,"");

    hLength->SetTitle(title); 
    hLength->GetXaxis()->SetTitle("LENGTH and WIDTH [#circ]");
    hLength->GetYaxis()->SetTitle("Entries");
    hLength->GetYaxis()->SetTitleOffset(1.6);
    hLength->SetFillColor(50);
    hLength->SetLineColor(50);
    hLength->SetLineWidth(2);
    hLength->SetFillStyle(3005);
    //hLength->SetMaximum(55000);
    hLength->SetName("LENGTH");

    TH1F* hWidth = new TH1F("hWidth","",100,0,1.1);

    CON->Draw("width>>hWidth",cut5,"");
 
    gStyle->SetStatTextColor(4);
    hWidth->SetFillColor(4);
    hWidth->SetLineColor(4);
    hWidth->SetFillStyle(3004);
    hWidth->SetLineWidth(2);
    hWidth->SetName("WIDTH");

    TLegend* leg = new TLegend(.59,.75,.8,.84);
    leg->AddEntry(hLength,"LENGTH","F");
    leg->AddEntry(hWidth,"WIDTH","F");
    leg->SetFillColor(0);
    leg->SetLineColor(1);
    leg->SetBorderSize(0);

    hLength->Draw();//DrawNormalized or DrawCopy are possible options
    hWidth->Draw("sames");
    //leg->Draw("same");

    // Beautify the statistic boxes!
    gPad->Update();
    TPaveStats* pavstat = (TPaveStats*) hLength->GetListOfFunctions()->FindObject("stats");
    if(pavstat)
    {
	Float_t shiftx = pavstat->GetX2NDC()-pavstat->GetX1NDC();
	pavstat->SetX1NDC(pavstat->GetX1NDC()-shiftx);
	pavstat->SetX2NDC(pavstat->GetX2NDC()-shiftx);
	pavstat->SetTextColor(50);
    }

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


    //************ Dist ************
    
    TCanvas* CDist = new TCanvas("CDist","Canvas",650,0,300,300);
    CDist->cd();
    CDist->SetBorderMode(0); 
    gStyle->SetOptStat(110011);
    //CDist->SetLogy();

    TH1F* hDist = new TH1F("hDist","",90,0,1.8);

    CON->Draw("dist>>hDist");

    hDist->GetXaxis()->SetTitle("DIST [#circ]");
    hDist->GetYaxis()->SetTitle("Entries");
    hDist->GetYaxis()->SetTitleOffset(1.6);
    //hDist->GetXaxis()->SetTitleOffset(1.3);
    hDist->SetLineColor(4);
    hDist->SetMarkerColor(4);
    hDist->SetFillColor(4);
    hDist->SetFillStyle(3004);
    hDist->SetLineWidth(2);
    hDist->SetTitle(title);
    hDist->SetName("DIST");
    
    hDist->Draw();


    //************ MeanY vs MeanX ************
    
    TCanvas* CXY = new TCanvas("CXY","Canvas",975,0,300,300);
    CXY->cd();
    CXY->SetBorderMode(0);  
    gStyle->SetOptStat(0);
    
    TH2F* hXY = new TH2F("hXY","",500,-1.5,1.5,500,-1.5,1.5);

    CON->Draw("meany:meanx>>hXY",cut3,"");

    TString titlexy = TString(title);

    hXY->SetTitle(titlexy); 
    hXY->GetXaxis()->SetTitle("MeanX [#circ]");
    hXY->GetYaxis()->SetTitle("MeanY [#circ]");
    hXY->GetYaxis()->SetTitleOffset(1.2);
    hXY->SetMaximum(15);
     
    hXY->Draw("COLZ");


    //************ DCA ************
    /*     
    TCanvas* CDCA = new TCanvas("CDCA","Canvas",10,400,300,300);
    CDCA->cd();
    CDCA->SetBorderMode(0); 
    gStyle->SetOptStat(110011);
    
    TH1F* hDCA = new TH1F("hDCA","",100,-1.5,1.5);

    CON->Draw("dca>>hDCA");

    hDCA->GetXaxis()->SetTitle("DCA [#circ]");
    hDCA->GetYaxis()->SetTitle("Entries");
    hDCA->GetYaxis()->SetTitleOffset(2);
    hDCA->SetLineColor(4);
    hDCA->SetFillColor(4);
    hDCA->SetFillStyle(3004);
    hDCA->SetTitle(title);
    hDCA->SetName("DCA");

    hDCA->Draw();
    */    


    //************ Delta ************
     
    TCanvas* CDelta = new TCanvas("CDelta","Canvas",0,335,300,300);
    CDelta->cd();
    CDelta->SetBorderMode(0); 
    gStyle->SetOptStat(110011);

    TH1F* hDelta = new TH1F("hDelta","",100,0,90);

    CON->Draw("delta>>hDelta",cut3);
    
    hDelta->GetXaxis()->SetTitle("#delta [#circ]");
    hDelta->GetYaxis()->SetTitle("Entries");
    hDelta->GetYaxis()->SetTitleOffset(1.4);
    hDelta->SetLineColor(4);
    hDelta->SetFillColor(4);
    hDelta->SetFillStyle(3004);
    hDelta->SetLineWidth(2);
    hDelta->SetTitle(title);
    hDelta->SetName("Delta");

    hDelta->Draw();


    //************ Size ************
     
    TCanvas* CSize = new TCanvas("CSize","Canvas",325,335,300,300);
    CSize->cd();
    CSize->SetBorderMode(0); 
    gStyle->SetOptStat(110011);

    TH1F* hSize = new TH1F("hSize","",100,0.1,5000);

    CON->Draw("size>>hSize");
    
    hSize->GetXaxis()->SetTitle("Size [Photons]");
    hSize->GetYaxis()->SetTitle("Entries");
    hSize->GetYaxis()->SetTitleOffset(1.6);
    hSize->GetXaxis()->SetTitleOffset(1.2);
    hSize->SetLineColor(4);
    hSize->SetFillColor(4);
    hSize->SetLineWidth(2);
    hSize->SetFillStyle(3004);
    hSize->SetTitle(title);
    hSize->SetName("SIZE");
    
    hSize->Draw();


    //************ Alpha ************
       
    TCanvas* CAlpha = new TCanvas("CAlpha","Canvas",650,335,300,300);
    CAlpha->cd();
    CAlpha->SetBorderMode(0); 
    CAlpha->SetGridx();
    CAlpha->SetGridy();
    gStyle->SetOptStat(11);
    
    TH1F *hAlpha = (TH1F*)new TH1F("hAlpha","Alpha ON",nbins,0,90);
    TH1F *hAlphaoff = (TH1F*)new TH1F("hAlphaoff","Alpha OFF",nbins,0,90);

    CON->Draw("abs(alpha)>>hAlpha",cut1);

    COFF->Draw("abs(alpha)>>hAlphaoff",cut1);


    //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 
    // Normalize ON and OFF data
    cout << "******************************************************" << endl;
    cout << " Bin error before scaling: " << hAlphaoff->GetBinError(2) << endl;

    // If TH1::Sumw2 has been called before filling, the sum of squares of
    // weights is also stored - I do not understand this!!
    hAlphaoff->Sumw2();
    
    Float_t level = 0, leveloff = 0;
    
    // Calculate the mean of entries beside the signal region
    for(Int_t ibin = inibin; ibin<=nbins; ibin++)
	level += hAlpha->GetBinContent(ibin);
    level /= (nbins-inibin+1);
    
    cout << " Mean beside signal: " << level << endl;

    for(Int_t ibin = inibin; ibin<=nbins; ibin++)
	leveloff += hAlphaoff->GetBinContent(ibin);
    leveloff /= (nbins-inibin+1);

    // Compare the mean values of both histograms
    if(leveloff>0)
	const Float_t norm = level/leveloff;//*nbins/hAlphaoff->GetEntries();
    cout << " Normalizing by factor " << norm <<endl;
    hAlphaoff->Scale(norm);
    
    cout << " Bin error after scaling: " << hAlphaoff->GetBinError(2) << endl;

    cout << " Cut on (ALPHA): " << cut1 << endl;
    //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 

   //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 
    // Significance calculation
    Float_t significance = 0;
    const Int_t signbins = 3;
    
    Int_t Non  = 0;
    Int_t Noff = 0;

    cout << " Significance Bins: " << signbins << endl;

    // Determine number of excess events and off events in the same region 
    for(Int_t ibin = 1; ibin<=signbins;ibin++)
    {
	Non +=hAlpha->GetBinContent(ibin);
	Noff+=hAlphaoff->GetBinContent(ibin);
    }

    cout << " Non = " << Non << " Noff = " << Noff << endl;

    significance = Significance(Non, Noff, 1.0, 0); //norm

    cout << " Significance: " << significance << endl;
    cout << "******************************************************" << endl;
    //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 

    // Set some draw options and draw the histograms
    hAlphaoff->GetXaxis()->SetTitle("ALPHA [#circ]");
    hAlphaoff->GetYaxis()->SetTitle("Entries");
    hAlphaoff->GetYaxis()->SetTitleOffset(1.4);
    hAlphaoff->GetXaxis()->SetTitleOffset(1.1);

    char legentry[128];
    sprintf(legentry,"  (Sign.: %2.3f #sigma,  Norm.: %2.3f,  |ALPHA| #leq %d#circ)", significance, norm, 90/nbins*signbins);

    hAlphaoff->SetTitle(title+legentry);
    hAlphaoff->SetName("OFF source");
    hAlphaoff->SetLineColor(4);
    hAlphaoff->SetLineWidth(2);
    hAlphaoff->SetMarkerColor(8);
    hAlphaoff->SetMarkerStyle(21);
    hAlphaoff->SetMarkerSize(0);
    hAlphaoff->SetFillColor(18);
    hAlphaoff->SetMaximum(hAlpha->GetMaximum()+200);

    hAlpha->SetLineColor(1);
    hAlpha->SetLineWidth(2);
    hAlpha->SetMarkerColor(50);
    hAlpha->SetMarkerStyle(20);
    hAlpha->SetMarkerSize(0.9);
    hAlpha->SetName("ON source");
    //hAlpha->SetMaximum(hAlpha->GetMaximum());

    hAlphaoff->Draw("eH");
    hAlpha->Draw("e1sames");

    // Beautify the statistic boxes!
    gPad->Update();
    TPaveStats* pavstat = (TPaveStats*) hAlphaoff->GetListOfFunctions()->FindObject("stats");
    if(pavstat)
    {
	Float_t shiftx = pavstat->GetX2NDC()-pavstat->GetX1NDC();
	pavstat->SetX1NDC(pavstat->GetX1NDC()-shiftx);
	pavstat->SetX2NDC(pavstat->GetX2NDC()-shiftx);
	pavstat->SetTextColor(4);
    }

    TPaveStats* pavstaton = (TPaveStats*) hAlpha->GetListOfFunctions()->FindObject("stats");
    if(pavstaton)
    {
	pavstaton->SetTextColor(50);
    }

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


    //************ Length/Size ************
    
    TCanvas* CLSiz = new TCanvas("CLSiz","Canvas",975,335,300,300);
    CLSiz->cd();
    CLSiz->SetBorderMode(0); 
    gStyle->SetOptStat(110011);

    TH1F* hLenSiz = new TH1F("hLenSiz","",100,0,0.005);

    CON->Draw("length/size>>hLenSiz",cut2,"");

    hLenSiz->GetXaxis()->SetTitle("LENGTH/SIZE [#circ/Photons]");
    hLenSiz->GetYaxis()->SetTitle("Entries");
    hLenSiz->GetYaxis()->SetTitleOffset(1.5);
    hLenSiz->GetXaxis()->SetTitleOffset(1.1);
    hLenSiz->SetLineColor(4);
    hLenSiz->SetFillColor(4);
    hLenSiz->SetLineWidth(2);
    hLenSiz->SetFillStyle(3004);
    hLenSiz->SetNdivisions(505);

    hLenSiz->SetTitle(title);
    hLenSiz->SetName("LENGTH/SIZE");

    hLenSiz->Draw();


    // Uncomment the following lines for MC only!

    //************ Energy ************
    /*
   
    TH1F *hEnergy = (TH1F*)gDirectory->Get("hEnergy");
    
    TCanvas* CEnergy = new TCanvas("CEnergy","Canvas",0,670,300,300);
    CEnergy->cd();
    CEnergy->SetBorderMode(0); 
    //CEnergy->SetLogy();
    gStyle->SetOptStat(110011);

    TH1F* hEnergy = new TH1F("hEnergy","",100,0,1000);

    CON->Draw("energy>>hEnergy","","");
     
    hEnergy->GetXaxis()->SetTitle("Energy [GeV]");
    hEnergy->GetYaxis()->SetTitle("Entries");
    hEnergy->GetYaxis()->SetTitleOffset(1.2);
    hEnergy->GetXaxis()->SetTitleOffset(1.1);
    hEnergy->SetLineColor(4);
    hEnergy->SetFillColor(4);
    hEnergy->SetLineWidth(2);
    hEnergy->SetFillStyle(3004);
    hEnergy->SetNdivisions(505);

    hEnergy->SetTitle(title);
    hEnergy->SetName("Energy");

    hEnergy->Draw();
    */


    //************ MeanX and MeanY ************

    // Create Canvas and set Canvas some options
    TCanvas* Cxy = new TCanvas("Cxy","Canvas",325,670,300,300);
    Cxy->cd();    
    Cxy->SetBorderMode(0);   
    gPad->SetBorderMode(0);
    gStyle->SetOptStat(111);

    TH1F* hMeanX = new TH1F("hMeanX","",100,-1.5,1.5);
    
    CON->Draw("meanx>>hMeanX",cut3,"");

    TH1F* hMeanY = new TH1F("hMeanY","",100,-1.5,1.5);

    CON->Draw("meany>>hMeanY",cut3,"");
    
    hMeanX->SetTitle(title); 
    hMeanX->GetXaxis()->SetTitle("MeanX and MeanY [#circ]");
    hMeanX->GetYaxis()->SetTitle("Entries");
    hMeanX->GetYaxis()->SetTitleOffset(1.6);
    hMeanX->SetFillColor(50);
    hMeanX->SetLineColor(50);
    hMeanX->SetLineWidth(2);
    hMeanX->SetFillStyle(3005);
    //hMeanX->SetMaximum(16000);
    //hMeanX->SetMaximum(6500);
    hMeanX->SetName("MeanX");
 
    gStyle->SetStatTextColor(4);
    hMeanY->SetFillColor(4);
    hMeanY->SetLineColor(4);
    hMeanY->SetFillStyle(3004);
    hMeanY->SetLineWidth(2);
    hMeanY->SetName("MeanY");

    hMeanX->Draw();//DrawNormalized or DrawCopy are possible options
    hMeanY->Draw("sames");//esames

    // Beautify the statistic boxes!
    gPad->Update();
    TPaveStats* pavstat = (TPaveStats*) hMeanX->GetListOfFunctions()->FindObject("stats");
    if(pavstat)
    {
	Float_t shiftx = pavstat->GetX2NDC()-pavstat->GetX1NDC();
	pavstat->SetX1NDC(pavstat->GetX1NDC()-shiftx);
	pavstat->SetX2NDC(pavstat->GetX2NDC()-shiftx);
	pavstat->SetTextColor(50);
    }

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

    //delete hWL;
    //delete hXY;
}


Double_t Significance(Float_t Non, Float_t Noff, Float_t alpha, Int_t type = 0)
{
    // Calculate significance after Li and Ma
    
    // Formula Nr. 17
    if(type == 0)
	return TMath::Sqrt(2*(Non*TMath::Log(((1+alpha)/alpha)*(Non/(Non+Noff))) + Noff*TMath::Log((1+alpha)*Noff/(Non+Noff))));

    // Formula Nr. 9
    if(type == 1)
	return (Non-alpha*Noff)/TMath::Sqrt(alpha*(Non+Noff));

}