#include "TString.h"
#include "TMultiGraph.h"
#include "TGraphErrors.h"
#include "TStyle.h"
#include "TCanvas.h"
#include "TLegend.h"

#include "MParList.h"
#include "MTaskList.h"
#include "MReadMarsFile.h"
#include "MGeomApply.h"
#include "MEvtLoop.h"
#include "MPedCalcPedRun.h"
#include "MPedestalCam.h"
#include "MPedestalPix.h"

#include <iostream>
#include <fstream>
#include <stdlib.h>

using namespace std;


const UInt_t nmaxslice = 30;
const UInt_t step = 2;

void pedvsslices2(TString pedname);

int main()
//void pedvsslices()
{
  pedvsslices2("~markus/Magic/Soft/Mars/20040303_20123_P_NewCalBoxTestLidOpen_E.root");
}

void pedvsslices2(TString pedname)
{
  const Int_t nbad = 22;
  Int_t bad[nbad]={0,10,47,54,124,125,159,162,222,300,388,395,464,558,559,543,89,510,576,459,550,561};

  static const UInt_t rnmax   = (nmaxslice/step)*step;
  static const UInt_t nslices = rnmax/step;
  static const UInt_t nloops  = rnmax*(rnmax-step+2)/(2*step);

  cout << "nslices = " << nslices << endl;

  // general settings
  gStyle->SetCanvasColor(0);
  gStyle->SetStatColor(0);
  gStyle->SetStatBorderSize(1);
  gStyle->SetMarkerSize(0.6);
  gStyle->SetPadBorderMode(0);

  // Define the graphs and their atributes
  TGraphErrors grmeanIn[nslices];
  TGraphErrors grrmsIn[nslices];
  TGraphErrors grmeanOut[nslices];
  TGraphErrors grrmsOut[nslices];

  // Null graph to adjust plot limits
  Float_t nullx[rnmax+1];
  Float_t nullyM[rnmax+1];
  Float_t nullyR[rnmax+1];
  for(UInt_t i=0;i<=rnmax;i++)
    {
      nullx[i]=(Float_t)i;
      nullyM[i]=9+2./rnmax*i;
      nullyR[i]=3.5+5./rnmax*i;
    }
  TGraph grnullM(rnmax+1,nullx,nullyM);
  TGraph grnullR(rnmax+1,nullx,nullyR);
  grnullM.SetMarkerSize(0);
  grnullR.SetMarkerSize(0);
  grnullM.SetMarkerColor(0);
  grnullR.SetMarkerColor(0);
  grnullM.GetXaxis()->SetTitle("First used time slice");
  grnullR.GetXaxis()->SetTitle("First used time slice");
  grnullM.GetYaxis()->SetTitle("Average Pedestal mean (ADC cnts)");
  grnullR.GetYaxis()->SetTitle("Average Pedestal rms (ADC cnts)");

  for(UInt_t i=0;i<nslices;i++)
    {
      grmeanIn[i].SetMarkerColor(i%9+1);
      grrmsIn[i].SetMarkerColor(i%9+1);
      grmeanOut[i].SetMarkerColor(i%9+1);
      grrmsOut[i].SetMarkerColor(i%9+1);

      grmeanIn[i].SetMarkerStyle(20+i%19+1);
      grrmsIn[i].SetMarkerStyle(20+i%19+1);
      grmeanOut[i].SetMarkerStyle(20+i%19+1);
      grrmsOut[i].SetMarkerStyle(20+i%19+1);
    }


  // Fill the graphs
  UInt_t islice=0;
  UInt_t fslice=step-1;    
  for(UInt_t i=0;i<nloops;i++)
    {
      MParList plist;
      MTaskList tlist;
      MPedestalCam pedcam;
  
      plist.AddToList(&tlist);
      plist.AddToList(&pedcam);
      
      MReadMarsFile read("Events", pedname);
      read.DisableAutoScheme();
      MGeomApply     geomapl;
      MPedCalcPedRun pedcalc;
      cout << "Checking slices from " << islice << " to " << fslice << endl;
      pedcalc.SetRange(islice,fslice,0,0);
      pedcalc.SetWindowSize(fslice-islice+1);

      tlist.AddToList(&read);
      tlist.AddToList(&geomapl);
      tlist.AddToList(&pedcalc);
      
      MEvtLoop evtloop;
      evtloop.SetParList(&plist);
      if (!evtloop.Eventloop())
	return;

      // compute the pedestal mean and rms averaged in the camera
      Float_t meansumIn =0;
      Float_t meansum2In=0;
      Float_t rmssumIn  =0;
      Float_t rmssum2In =0;
      Float_t meansumOut =0;
      Float_t meansum2Out=0;
      Float_t rmssumOut  =0;
      Float_t rmssum2Out =0;

      const Int_t ninn = 397;
      Int_t rninn=0;
      Int_t rnout=0;
      for(Int_t j=0;j<pedcam.GetSize();j++)
	{
	  // skip bad pixels from list above
	  Bool_t isBad=kFALSE;
	  for(Int_t k=0;k<nbad;k++)
	    if(bad[k]==j)
	      {
		isBad=kTRUE;
		break;
	      }
	  if(isBad)
	    continue;

	  // keep the sum of values and of squares
	  Float_t ped = pedcam[j].GetPedestal();
	  Float_t rms = pedcam[j].GetPedestalRms();
	  if(j<ninn)
	    {
	      meansumIn+= ped;
	      meansum2In+=ped*ped;
	      rmssumIn+=rms;
	      rmssum2In+=rms*rms;
	      rninn++;
	    }
	  else
	    {
	      meansumOut+= ped;
	      meansum2Out+=ped*ped;
	      rmssumOut+=rms;
	      rmssum2Out+=rms*rms;
	      rnout++;
	    }	    
	}
      
      Float_t avmeanIn  = meansumIn/rninn;
      Float_t sigmeanIn = (meansum2In-meansumIn*meansumIn/rninn)/(rninn-1);
      Float_t avrmsIn   = rmssumIn/rninn;
      Float_t sigrmsIn  = (rmssum2In-rmssumIn*rmssumIn/rninn)/(rninn-1);
       
      Float_t avmeanOut  = meansumOut/rnout;
      Float_t sigmeanOut = (meansum2Out-meansumOut*meansumOut/rnout)/(rnout-1);
      Float_t avrmsOut   = rmssumOut/rnout;
      Float_t sigrmsOut  = (rmssum2Out-rmssumOut*rmssumOut/rnout)/(rnout-1);
      
      // save the values into the graphs
      grmeanIn[(fslice-islice+1)/step-1].SetPoint(islice,islice,avmeanIn);
      grmeanIn[(fslice-islice+1)/step-1].SetPointError(islice,0,sigmeanIn/rninn);
      grmeanOut[(fslice-islice+1)/step-1].SetPoint(islice,islice,avmeanOut);
      grmeanOut[(fslice-islice+1)/step-1].SetPointError(islice,0,sigmeanOut/rnout);

      grrmsIn[(fslice-islice+1)/step-1].SetPoint(islice,islice,avrmsIn);
      grrmsIn[(fslice-islice+1)/step-1].SetPointError(islice,0,sigrmsIn/rninn);
      grrmsOut[(fslice-islice+1)/step-1].SetPoint(islice,islice,avrmsOut);
      grrmsOut[(fslice-islice+1)/step-1].SetPointError(islice,0,sigrmsOut/rnout);

      // compute next initial and final slice values
      islice++;
      fslice++;
      if(fslice>=nmaxslice)
	{
	  fslice=fslice-islice+step;
	  islice=0;
	}
    }

  TCanvas* myC = new TCanvas("myC","pedestal studies",900,1250);
  myC->Divide(2,2);

  TLegend* leg1 = new TLegend(.2,.2,.4,.5);
  TLegend* leg2 = new TLegend(.44,.2,.88,.5);
  leg1->SetEntrySeparation(0.1);
  leg1->SetHeader("");
  leg2->SetEntrySeparation(0.1);
  leg2->SetHeader("");

  myC->cd(1);
  grnullM.SetTitle("Average pedestal mean for inner pixels");
  grnullM.DrawClone("AP");
  for(UInt_t i=0;i<nslices;i++)
    {      
      grmeanIn[i].DrawClone("P");
      char leglab[100];
      sprintf(leglab,"%d slices",(i+1)*step);
      if(i<nslices/2)
	leg1->AddEntry(&grmeanIn[i],leglab,"p");
      else
	leg2->AddEntry(&grmeanIn[i],leglab,"p");
    }
  leg1->SetFillColor(0);
  leg1->SetLineColor(0);
  leg1->SetBorderSize(0);
  leg1->Draw();
  leg2->SetFillColor(0);
  leg2->SetLineColor(0);
  leg2->SetBorderSize(0);
  leg2->Draw();

  myC->cd(2);
  grnullR.SetTitle("Average pedestal rms for inner pixels");
  grnullR.DrawClone("AP");
  for(UInt_t i=0;i<nslices;i++)
    grrmsIn[i].DrawClone("P");
  
  myC->cd(3);
  grnullM.SetTitle("Average pedestal mean for outer pixels");
  grnullM.DrawClone("AP");
  for(UInt_t i=0;i<nslices;i++)
    grmeanOut[i].DrawClone("P");
  
  myC->cd(4);
  grnullR.SetTitle("Average pedestal rms for outer pixels");
  grnullR.DrawClone("AP");
  for(UInt_t i=0;i<nslices;i++)
    grrmsOut[i].DrawClone("P");

  myC->Print("pedestaSliceDependence.ps");

}