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

void cleanComp()
{    
  // general settings
  gROOT->Reset();
  gStyle->SetCanvasColor(0);
  gStyle->SetCanvasBorderMode(0);
  gStyle->SetPadBorderMode(0);
  gStyle->SetFrameBorderMode(0);
  gStyle->SetOptTitle(1);
  //  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;
  ton= new TChain("Events");
  
  // data file(s)
  ton->Add("/local_disk/jrico/mc/cleanComp0-3.root");
  
  // define aliases
  ton->SetAlias("run","MRawRunHeader.fRunNumber");
  ton->SetAlias("event","MRawEvtHeader.fDAQEvtNumber");
  ton->SetAlias("energy","MMcEvt.fEnergy");

  // alpha0 is for the file without noise
  // alphaNOCL is for the file with noise and no cleaning
  ton->SetAlias("alpha0","MHillasSrc0.fAlpha");
  ton->SetAlias("alphaNOCL","MHillasSrcnoclean.fAlpha");
  ton->SetAlias("alpha3015","MHillasSrc3015.fAlpha");
  ton->SetAlias("alpha4025","MHillasSrc4025.fAlpha");
  ton->SetAlias("alpha5035","MHillasSrc5035.fAlpha");
  ton->SetAlias("alpha6045","MHillasSrc6045.fAlpha");
  
  //  ton->Scan("run:event","alpha0<-30 && energy>2500 && energy<5000");
  
  //  plotAlphaWidth(ton);
   plotAlphaWidthVsEnergy(ton);
  // plotEfficiencies(ton);
}

void plotAlphaWidth(TChain* ton)
{
  // define and fill histos
  const Int_t nbins=90;
  TH1F* halpha3015 = new TH1F("Alpha3015" ,"Alpha3015",nbins,-180,180);
  TH1F* halpha4025 = new TH1F("Alpha4025" ,"Alpha4025",nbins,-180,180);
  TH1F* halpha5035 = new TH1F("Alpha5035" ,"Alpha5035",nbins,-180,180);
  TH1F* halpha6045 = new TH1F("Alpha6045" ,"Alpha6045",nbins,-180,180);
  
  halpha3015->SetXTitle("alpha-alpha0 (deg)");
  halpha3015->SetYTitle("Entries");
  
  halpha3015->SetLineColor(1);
  halpha4025->SetLineColor(2);
  halpha5035->SetLineColor(3);
  halpha6045->SetLineColor(4);
  
  ton->Draw("alpha3015-alpha0>>Alpha3015","alpha3015!=0");
  ton->Draw("alpha4025-alpha0>>Alpha4025","alpha4025!=0");
  ton->Draw("alpha5035-alpha0>>Alpha5035","alpha5035!=0");
  ton->Draw("alpha6045-alpha0>>Alpha6045","alpha6045!=0");
  
  // plot
  TCanvas *c1 = new TCanvas("c1","c1",800,500);
  c1->cd(1);
  
  halpha3015->Draw("HIST");
  halpha4025->Draw("HISTSAMES");
  halpha5035->Draw("HISTSAMES");
  halpha6045->Draw("HISTSAMES");
  
  gPad->Update();
  TPaveStats* pvstat4025  = (TPaveStats*) halpha4025->GetListOfFunctions()->FindObject("stats");
  TPaveStats* pvstat5035  = (TPaveStats*) halpha5035->GetListOfFunctions()->FindObject("stats");
  TPaveStats* pvstat6045  = (TPaveStats*) halpha6045->GetListOfFunctions()->FindObject("stats");
  
  Float_t shifty = pvstat4025->GetY2NDC()-pvstat4025->GetY1NDC();
  
  pvstat4025->SetY1NDC(pvstat4025->GetY1NDC()-shifty);
  pvstat4025->SetY2NDC(pvstat4025->GetY2NDC()-shifty);
  pvstat5035->SetY1NDC(pvstat5035->GetY1NDC()-2*shifty);
  pvstat5035->SetY2NDC(pvstat5035->GetY2NDC()-2*shifty);
  pvstat6045->SetY1NDC(pvstat6045->GetY1NDC()-3*shifty);
  pvstat6045->SetY2NDC(pvstat6045->GetY2NDC()-3*shifty);
  
  gPad->Modified();
  gPad->Update();
  
  c1->Print("cleanComp.ps");
}

/*******************************************************************************/
void plotAlphaWidthVsEnergy(TChain* ton)
{
  const Int_t nbins = 8;
  Float_t bound[nbins+1] ={0,30,50,80,150,300,600,1200,3000}; // boundary of considered size bins
  
  // define the arrays that will feed the graphs
  // value of rms of the alpha-alpha0 distribution, and errors
  Float_t alpha0[nbins];
  Float_t alpha3015[nbins];
  Float_t alpha4025[nbins];
  Float_t alpha5035[nbins];
  Float_t alpha6045[nbins];

  Float_t ealpha0[nbins];
  Float_t ealpha3015[nbins];
  Float_t ealpha4025[nbins];
  Float_t ealpha5035[nbins];
  Float_t ealpha6045[nbins];

  // x axis and error
  Float_t x[nbins];
  Float_t ex[nbins];

  TH1F* hp[nbins];

  // define the function that will be used to fit the histograms
  TF1 f2("f2","[0]+[1]/[3]*exp(-(x-[2])*(x-[2])/(2*[3]*[3]))",-180.,180.);
  f2.SetParameters(0,100,0,02);
  f2.FixParameter(0,0);
  f2.SetLineWidth(1);

  // fill the arrays
  for(Int_t i=0;i<nbins;i++)
    {
      //      cout << "Bin " << i << endl;
      TH1F dummy0("dummy0" ,"dummy0",100,-90,90);
      TH1F dummy3015("dummy3015" ,"dummy3015",100,-180,180);
      TH1F dummy4025("dummy4025" ,"dummy4025",100,-180,180);
      TH1F dummy5035("dummy5035" ,"dummy5035",100,-180,180);
      TH1F dummy6045("dummy6045" ,"dummy6045",100,-180,180);

      TH1F xdummy("xdummy" ,"xdummy",100,bound[i],bound[i+1]);

      // define the cuts for each of the cleaning levels
      char cut[100];
      sprintf(cut,"energy>%f && energy<%f && alpha0!=0",bound[i],bound[i+1]);
      char cut3015[200];
      sprintf(cut3015,"%s && alpha3015!=0",cut);
      char cut4025[200];
      sprintf(cut4025,"%s && alpha4025!=0",cut);
      char cut5035[200];
      sprintf(cut5035,"%s && alpha5035!=0",cut);
      char cut6045[200];
      sprintf(cut6045,"%s && alpha6045!=0",cut);

      TF1 f1("f1","[0]+[1]/[3]*exp(-(x-[2])*(x-[2])/(2*[3]*[3]))",-180.,180.);
      f1.SetParameters(0,100,0,02);
      f1.FixParameter(0,0);


      // fill the histograms to extract later the information from rms, entries, means, etc
      ton->Draw("alpha0>>dummy0",cut);
      ton->Draw("alpha3015-alpha0>>dummy3015",cut3015);
      ton->Draw("alpha4025-alpha0>>dummy4025",cut4025);
      ton->Draw("alpha5035-alpha0>>dummy5035",cut5035);
      ton->Draw("alpha6045-alpha0>>dummy6045",cut6045);

      ton->Draw("energy>>xdummy",cut);

      // fit to a gaussian and get the sigma to be plotted later
      dummy0.Fit("f1");
      alpha0[i]   =f1.GetParameter(3);
      dummy3015.Fit("f1");
      alpha3015[i]=f1.GetParameter(3);
      dummy4025.Fit("f1");
      alpha4025[i]=f1.GetParameter(3);
      dummy5035.Fit("f1");
      alpha5035[i]=f1.GetParameter(3);
      dummy6045.Fit("f1");
      alpha6045[i]=f1.GetParameter(3);

      Float_t nevt0    = dummy0.GetEntries();
      Float_t nevt3015 = dummy3015.GetEntries();
      Float_t nevt4025 = dummy4025.GetEntries();
      Float_t nevt5035 = dummy5035.GetEntries();
      Float_t nevt6045 = dummy6045.GetEntries();

      // since gaussian is not so good aproximation, compute error like follows
      ealpha0[i]   = alpha0[i]/TMath::Sqrt(2*nevt0);
      ealpha3015[i]= alpha3015[i]/TMath::Sqrt(2*nevt3015);
      ealpha4025[i]= alpha4025[i]/TMath::Sqrt(2*nevt4025);
      ealpha5035[i]= alpha5035[i]/TMath::Sqrt(2*nevt5035);
      ealpha6045[i]= alpha6045[i]/TMath::Sqrt(2*nevt6045);

      // copy some histos to be plotted later
      Char_t nom[10];
      sprintf(nom,"hp[%02d]",i);
      hp[i] = new TH1F(nom,nom,100,-180,180);
      hp[i]->GetXaxis()->SetTitle("alpha (deg)");
      hp[i]->Add(&dummy6045);
      hp[i]->Fit("f2");

      //      cout << xdummy.GetMean() << " " << xdummy.GetRMS()<< " " << xdummy.GetEntries() << endl;
      x[i]  = xdummy.GetMean();            
      //      Float_t rms = xdummy.GetRMS();
      ex[i] = 0;//rms*rms/xdummy.GetEntries();
    }

  //#if 0
  TCanvas* cdum = new TCanvas("cdum","cdum",600,900);
  cdum->Divide(2,nbins/2+1);
  for(Int_t i=0;i<nbins;i++)
    {
      cdum->cd(i+1);
      hp[i]->Draw();
      gPad->Modified();
      gPad->Update();
    }
  //#endif
  
  // create the graphs
  TGraphErrors* gr0    = new TGraphErrors(nbins,x,alpha0,ex,ealpha0);
  TGraphErrors* gr3015 = new TGraphErrors(nbins,x,alpha3015,ex,ealpha3015);
  TGraphErrors* gr4025 = new TGraphErrors(nbins,x,alpha4025,ex,ealpha4025);
  TGraphErrors* gr5035 = new TGraphErrors(nbins,x,alpha5035,ex,ealpha5035);
  TGraphErrors* gr6045 = new TGraphErrors(nbins,x,alpha6045,ex,ealpha6045);  

  // marker/line colors
  gr0->SetMarkerStyle(20);
  gr3015->SetMarkerStyle(24);
  gr4025->SetMarkerStyle(21);
  gr5035->SetMarkerStyle(25);
  gr6045->SetMarkerStyle(22);

  gr0->SetMarkerColor(1);
  gr3015->SetMarkerColor(2);
  gr4025->SetMarkerColor(3);
  gr5035->SetMarkerColor(4);
  gr6045->SetMarkerColor(6);

  gr0->SetLineColor(1);
  gr3015->SetLineColor(2);
  gr4025->SetLineColor(3);
  gr5035->SetLineColor(4);
  gr6045->SetLineColor(6);

  //////////////////////
  // plot alpha width //
  //////////////////////
  TCanvas *c1 = new TCanvas("c1","c1",800,500);
  c1->cd(1);
  gPad->SetGridx();  
  gPad->SetGridy();  
  gPad->SetLogx();

  TGraphErrors* grfirst = gr4025;
  grfirst->Draw("APL");
  gr0->Draw("PL");
  gr3015->Draw("PL");
  gr4025->Draw("PL");
  gr5035->Draw("PL");
  gr6045->Draw("PL");

  grfirst->GetHistogram()->SetTitle("Cleaning comparison (alpha)");
  grfirst->GetHistogram()->GetXaxis()->SetTitle("E_{MC} [GeV]");
  grfirst->GetHistogram()->GetYaxis()->SetTitle("#sigma (alpha-alpha_{0}) [deg]");
  grfirst->GetHistogram()->SetMaximum(80);

  TPaveLabel* label = new TPaveLabel(0.25,0.75,0.6,0.85,"MC Gammas - zbins 0-3","NDC");
  label->SetBorderSize(0);
  label->SetFillColor(0);
  label->Draw();

  // legend
  TLegend* leg = new TLegend(.6,.5,.8,.8);
  leg->SetHeader("");
  leg->SetFillColor(0);
  leg->SetLineColor(0);
  leg->SetBorderSize(0);
  leg->AddEntry(gr0,"No noise","p");
  leg->AddEntry(gr3015,"(3.0,1.5)","p");
  leg->AddEntry(gr4025,"(4.0,2.5)","p");
  leg->AddEntry(gr5035,"(5.0,3.5)","p");
  leg->AddEntry(gr6045,"(6.0,4.5)","p");
  leg->Draw();

  gPad->Modified();
  gPad->Update();
  // print to file
  c1->Print("cleanAlphaWidthVSEnergy.ps");
}

/*******************************************************************************/
void plotEfficiencies(TChain* ton)
{
  const Int_t nbins = 8;
  Float_t bound[nbins+1] ={0,30,50,80,150,300,600,1200,3000}; // boundary of considered size bins
  
  // define the arrays that will feed the graphs
  // gamma detection efficiency, and errors  
  Float_t effNOCL[nbins];
  Float_t eff3015[nbins];
  Float_t eff4025[nbins];
  Float_t eff5035[nbins];
  Float_t eff6045[nbins];

  Float_t eeffNOCL[nbins];
  Float_t eeff3015[nbins];
  Float_t eeff4025[nbins];
  Float_t eeff5035[nbins];
  Float_t eeff6045[nbins];

  // x axis and error
  Float_t x[nbins];
  Float_t ex[nbins];

  // fill the arrays
  for(Int_t i=0;i<nbins;i++)
    {
      //      cout << "Bin " << i << endl;
      TH1F dummyNOCL("dummyNOCL" ,"dummyNOCL",100,-90,90);
      TH1F dummy3015("dummy3015" ,"dummy3015",100,-90,90);
      TH1F dummy4025("dummy4025" ,"dummy4025",100,-90,90);
      TH1F dummy5035("dummy5035" ,"dummy5035",100,-90,90);
      TH1F dummy6045("dummy6045" ,"dummy6045",100,-90,90);

      TH1F xdummy("xdummy" ,"xdummy",100,bound[i],bound[i+1]);

      // define the cuts for each of the cleaning levels
      char cut[100];
      sprintf(cut,"energy>%f && energy<%f && alpha0!=0",bound[i],bound[i+1]);
      char cutNOCL[200];
      sprintf(cutNOCL,"%s && alphaNOCL!=0",cut);
      char cut3015[200];
      sprintf(cut3015,"%s && alpha3015!=0",cut);
      char cut4025[200];
      sprintf(cut4025,"%s && alpha4025!=0",cut);
      char cut5035[200];
      sprintf(cut5035,"%s && alpha5035!=0",cut);
      char cut6045[200];
      sprintf(cut6045,"%s && alpha6045!=0",cut);


      // fill the histograms to extract later the information from rms, entries, means, etc
      ton->Draw("alpha0>>dummyNOCL",cutNOCL);
      ton->Draw("alpha0>>dummy3015",cut3015);
      ton->Draw("alpha0>>dummy4025",cut4025);
      ton->Draw("alpha0>>dummy5035",cut5035);
      ton->Draw("alpha0>>dummy6045",cut6045);

      ton->Draw("energy>>xdummy",cut);

      Float_t nevt0    = dummyNOCL.GetEntries();
      Float_t nevt3015 = dummy3015.GetEntries();
      Float_t nevt4025 = dummy4025.GetEntries();
      Float_t nevt5035 = dummy5035.GetEntries();
      Float_t nevt6045 = dummy6045.GetEntries();

      // fill the array of efficiencies
      eff3015[i] = nevt3015/nevt0*100;
      eff4025[i] = nevt4025/nevt0*100;
      eff5035[i] = nevt5035/nevt0*100;
      eff6045[i] = nevt6045/nevt0*100;

      eeff3015[i] = TMath::Sqrt(nevt3015)/nevt0*100;
      eeff4025[i] = TMath::Sqrt(nevt4025)/nevt0*100;
      eeff5035[i] = TMath::Sqrt(nevt5035)/nevt0*100;
      eeff6045[i] = TMath::Sqrt(nevt6045)/nevt0*100;

      //      cout << xdummy.GetMean() << " " << xdummy.GetRMS()<< " " << xdummy.GetEntries() << endl;
      x[i]  = xdummy.GetMean();            
      //      Float_t rms = xdummy.GetRMS();
      ex[i] = 0;//rms*rms/xdummy.GetEntries();
    }

  // create the graphs
  TGraphErrors* efgr3015 = new TGraphErrors(nbins,x,eff3015,ex,eeff3015);
  TGraphErrors* efgr4025 = new TGraphErrors(nbins,x,eff4025,ex,eeff4025);
  TGraphErrors* efgr5035 = new TGraphErrors(nbins,x,eff5035,ex,eeff5035);
  TGraphErrors* efgr6045 = new TGraphErrors(nbins,x,eff6045,ex,eeff6045);  

  // marker/line colors
  efgr3015->SetMarkerStyle(24);
  efgr4025->SetMarkerStyle(21);
  efgr5035->SetMarkerStyle(25);
  efgr6045->SetMarkerStyle(22);

  efgr3015->SetMarkerColor(2);
  efgr4025->SetMarkerColor(3);
  efgr5035->SetMarkerColor(4);
  efgr6045->SetMarkerColor(6);

  efgr3015->SetLineColor(2);
  efgr4025->SetLineColor(3);
  efgr5035->SetLineColor(4);
  efgr6045->SetLineColor(6);

  ///////////////////////
  // plot efficiencies //
  ///////////////////////
  TCanvas *c1 = new TCanvas("c1","c1",800,500);
  c1->cd(1);
  gPad->SetGridx();  
  gPad->SetGridy();  
  gPad->SetLogx();

  TGraphErrors* efgrfirst = efgr6045;
  efgrfirst->Draw("APL");
  efgr3015->Draw("PL");
  efgr4025->Draw("PL");
  efgr5035->Draw("PL");
  efgr6045->Draw("PL");

  efgrfirst->GetHistogram()->SetTitle("Cleaning comparison (efficiency)");
  efgrfirst->GetHistogram()->GetXaxis()->SetTitle("E_{MC} [GeV]");
  efgrfirst->GetHistogram()->GetYaxis()->SetTitle("Detection efficiency (%)");

  // label
  TPaveLabel* eflabel = new TPaveLabel(0.5,0.49,0.85,0.59,"MC Gammas - zbins 0-3","NDC");
  eflabel->SetBorderSize(0);
  eflabel->SetFillColor(0);
  eflabel->Draw();
 
  // legend
  TLegend* efleg = new TLegend(.6,.2,.8,.5);
  efleg->SetHeader("");
  efleg->SetFillColor(0);
  efleg->SetLineColor(0);
  efleg->SetBorderSize(0);
  efleg->AddEntry(efgr3015,"(3.0,1.5)","p");
  efleg->AddEntry(efgr4025,"(4.0,2.5)","p");
  efleg->AddEntry(efgr5035,"(5.0,3.5)","p");
  efleg->AddEntry(efgr6045,"(6.0,4.5)","p");
  efleg->Draw();

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

  // print to file
  c1->Print("cleanEfficiency.ps");
}