/* ======================================================================== *\
!
! *
! * This file is part of MARS, the MAGIC Analysis and Reconstruction
! * Software. It is distributed to you in the hope that it can be a useful
! * and timesaving tool in analysing Data of imaging Cerenkov telescopes.
! * It is distributed WITHOUT ANY WARRANTY.
! *
! * Permission to use, copy, modify and distribute this software and its
! * documentation for any purpose is hereby granted without fee,
! * provided that the above copyright notice appear in all copies and
! * that both that copyright notice and this permission notice appear
! * in supporting documentation. It is provided "as is" without express
! * or implied warranty.
! *
!
!
!   Author(s): Markus Gaug, 04/2004 <mailto:markus@ifae.es>
!
!   Copyright: MAGIC Software Development, 2000-2004
!
!
\* ======================================================================== */
//////////////////////////////////////////////////////////////////////////////
//
// pedestalstudies_pisa.C
//
// macro to study the pedestal and pedestalRMS with the number of FADC 
// slices summed up. 
//
/////////////////////////////////////////////////////////////////////////////////
const TString pedfile = "/data0/Temp/20040422_23203_P_Mrk421_E.root";

void pedestalstudies_pisa(const TString pedname=pedfile, Int_t pixin=0, Int_t pixout=397,TString f2="",TString f3="",TString f4="",TString f5="")
{

  Int_t loops = 1;
  Int_t stepsize = 2;

  gStyle->SetOptStat(1111);
  gStyle->SetOptFit();
  
  TArrayF *hmeandiffinn = new TArrayF(loops);
  TArrayF *hrmsdiffinn  = new TArrayF(loops);
  TArrayF *hmeandiffout = new TArrayF(loops);
  TArrayF *hrmsdiffout  = new TArrayF(loops);
  TArrayF *hmeaninn  = new TArrayF(loops);
  TArrayF *hmeanout  = new TArrayF(loops);
  TArrayF *hrmsinn   = new TArrayF(loops);
  TArrayF *hrmsout   = new TArrayF(loops);
  TArrayF *hmuinn    = new TArrayF(loops);
  TArrayF *hmuout    = new TArrayF(loops);
  TArrayF *hsigmainn = new TArrayF(loops);
  TArrayF *hsigmaout = new TArrayF(loops);

  TArrayF *hmeandiffinnerr = new TArrayF(loops);
  TArrayF *hrmsdiffinnerr  = new TArrayF(loops);
  TArrayF *hmeandiffouterr = new TArrayF(loops);
  TArrayF *hrmsdiffouterr  = new TArrayF(loops);
  TArrayF *hmeaninnerr  = new TArrayF(loops);
  TArrayF *hmeanouterr  = new TArrayF(loops);
  TArrayF *hrmsinnerr   = new TArrayF(loops);
  TArrayF *hrmsouterr   = new TArrayF(loops);
  TArrayF *hmuinnerr    = new TArrayF(loops);
  TArrayF *hmuouterr    = new TArrayF(loops);
  TArrayF *hsigmainnerr = new TArrayF(loops);
  TArrayF *hsigmaouterr = new TArrayF(loops);


  MStatusDisplay *display = new MStatusDisplay;
  display->SetUpdateTime(500);
  display->Resize(850,700);
      
  //
  // Create a empty Parameter List and an empty Task List
  // The tasklist is identified in the eventloop by its name
  //
  MParList  plist;
  MTaskList tlist;
  plist.AddToList(&tlist);

  TString tabname="";

  for (Int_t samples=2; samples < stepsize*loops+1; samples=samples+stepsize)
    {

      plist.Reset();
      tlist.Reset();
      
      //
      // Now setup the tasks and tasklist for the pedestals:
      // ---------------------------------------------------
      //
      
      MReadMarsFile read("Events", pedname);
      read.DisableAutoScheme();
     
	Int_t neventi;
 	neventi=1000;

	Float_t samplingRate=500;	//Hz
	if (f2!=""){
	read.AddFile(f2);
	neventi += 1000;
	}

	if (f3!=""){
	read.AddFile(f3);
	neventi += 1000;
	}

	if (f4!=""){
	read.AddFile(f4);
	neventi += 1000;
	}

	if (f5!=""){
	read.AddFile(f5);
	neventi += 1000;
	}

	const Int_t NEvt = neventi;
	cout << "E" <<neventi<< " NEvt:" << NEvt << endl;
	Float_t timevt[NEvt];

      MGeomApply      geomapl;
      //
      // Set the extraction range higher:
      //		
      //MExtractFixedWindow sigcalc;
      //MExtractSignal3 sigcalc;
      	MExtractSlidingWindow sigcalc;
      sigcalc.SetRange(0,samples-1,0,1);

     
  
      MPedCalcPedRun pedcalc;
      pedcalc.SetRange(0,samples-1,0,0);
      pedcalc.SetWindowSize((Int_t)sigcalc.GetNumHiGainSamples());

      //MPedCalcFromData pedcalc;
      //pedcalc.SetfHiGainThreshold(100);

      //
      // Additionally to calculating the pedestals, 
      // you can fill histograms and look at them
      //
      MFillH fill("MHPedestalCam", "MExtractedSignalCam");
      tabname="PedCam ";
      tabname += samples;
      //fill.SetNameTab(Form("%s%2d","PedCam",samples));
      fill.SetNameTab(tabname);

      tlist.AddToList(&read);
      tlist.AddToList(&geomapl);
      tlist.AddToList(&sigcalc);
      tlist.AddToList(&pedcalc);
      tlist.AddToList(&fill);
      
      const MGeomCamMagic      geomcam;
      MPedestalCam       pedcam;
      MBadPixelsCam      badcam;
      badcam.AsciiRead("badpixels.dat");  
      
      MHPedestalCam      hpedcam;
      MCalibrationPedCam cpedcam;
      
      plist.AddToList(&geomcam);
      plist.AddToList(&pedcam);
      plist.AddToList(&hpedcam);
      plist.AddToList(&cpedcam);
      plist.AddToList(&badcam);
      
      //
      // Create and setup the eventloop
      //
      MRawEvtHeader *evtheader;
      MRawEvtData *fRawEvt;
      TH1F *hpedevt = new TH1F("hpedevt","Pedestals Drift",NEvt,0,(Float_t)NEvt/samplingRate);
      Bool_t rc;
      MEvtLoop evtloop;
      
      evtloop.SetParList(&plist);
      evtloop.SetDisplay(display);
      //tlist.SetDisplay(display);
      //
      // Execute first analysis
      //
      //if (!evtloop.Eventloop())
      //  return;
      rc = evtloop.PreProcess();
      
      //while (tlist.Process())
	//{
	//cout << "bleah" << endl;
	//}      
//tlist.SetDisplay(display);

       Int_t fNumHiGainSample=15;
	Float_t peddata[NEvt];
	Int_t evtan=0;

      if (rc){
          while (tlist.Process())
          {
            fRawEvt=(MRawEvtData *)plist.FindObject("MRawEvtData");
            MRawEvtPixelIter pixel(fRawEvt);

            Float_t ped=0;
	    evtan++;

 	    while (pixel.Next())
	     {

 	       const UInt_t idx    = pixel.GetPixelId();

               // lista di pixel da sommare
               if (idx >pixout || idx < pixin)
		  continue;

	        Byte_t *ptr = pixel.GetHiGainSamples();
                const Byte_t *end = ptr + fNumHiGainSample;
      
                UInt_t sum = 0;
	        UInt_t sqr = -1;

	        do
	          {
	            sum += *ptr;
	            //sqr += *ptr * *ptr;
	          }
	        while (++ptr != end);
      
                 ped += (Float_t)sum;

                  //cout << " Pixel " << idx << " pede: " << ped << " --- ";
                }
 
               //rc = evtloop.Process(0);
	       evtheader = (MRawEvtHeader *)plist.FindObject("MRawEvtHeader");

             ped /= pixout-pixin+1;
		Int_t evtn=evtheader->GetDAQEvtNumber();
             evtn=evtan;

//test
	//ped =  150 + sin(((Float_t)evtn/16.-TMath::Floor(evtn/16))*6.28); 
	//cout << ped << " ";
	timeped = evtn/samplingRate;
	timevt[evtn-1]=timeped;
     
        hpedevt->Fill(timeped,ped);

            //	cout <<" evt:" << evtn;
	     // Save pedestal in array for Fourier analysis
	     peddata[evtn-1]=ped;
          }
      }
      if (!evtloop.PostProcess()){
          rc = kFALSE;
          return;
      }

      TCanvas &p1 = display->AddTab("ped studies");
      p1.Divide(1,2);
      p1.cd(1);
	hpedevt->GetXaxis()->SetTitle("Time (s)");
//      hpedevt->Draw();

	TGraph *drift = new TGraph(evtan,timevt,peddata);
	drift->SetMarkerStyle(21);
	drift->SetMarkerSize(0.25);
	drift->GetXaxis()->SetTitle("time (s)");
	drift->GetYaxis()->SetTitle(" pedestal");
	drift->Draw("APL");


      // Fourier analysis
	// campionamento a 500Hz (1000 evt) => freq. max a 250Hz

	MFFT pro;
	TArrayF *pedarray = new TArrayF(NEvt,peddata);
	TArrayF *fft; // array con risultato FFT 
	fft = pro.RealFunctionFFT(pedarray);
	Int_t nsample = (Int_t)(fft->GetSize()); 
	
	const Int_t NNEvt = nsample;

	Float_t x[10000];
	Float_t prova;
	for (int i=0;i<NNEvt;i++)
	 x[i]=(Float_t)i/(NNEvt/(samplingRate/2)); 

	Float_t *ffty = fft->GetArray();
        cout << "Szie:" << fft->GetSize();
	ffty[0]=0.;
	TGraph *gr = new TGraph(nsample,x,ffty);
	p1.cd(2);
	gr->GetXaxis()->SetTitle("Frequency (Hz)");
	gr->SetMarkerStyle(20);
	gr->SetMarkerSize(0.4);
	gr->Draw("AP");		



//-------------------------------


      // Set status lines
      display->SetStatusLine1(evtloop.GetName());
      display->SetStatusLine2(rc ? "Done." : "Error!");
      // Stop automatic update
      display->StopUpdate();
      // Reallow context menus
      display->SetNoContextMenu(kFALSE);
      // Reallow user to exit window by File menu
      display->UnLock();
      // 
      // Look at one specific pixel, after all the histogram manipulations:
      //
      /*
      MHGausEvents &hpix = hpedcam.GetAverageHiGainArea(0);
      hpix.DrawClone("fourierevents");
      
      MHGausEvents &lpix = hpedcam.GetAverageHiGainArea(1);
      lpix.DrawClone("fourierevents");

      hpedcam[170].DrawClone("fourierevents");
  
      */

      MHCamera dispped0  (geomcam, "Ped;Pedestal",       "Mean per Slice");
      MHCamera dispped2  (geomcam, "Ped;PedestalRms",    "RMS per Slice");
      MHCamera dispped4  (geomcam, "Ped;Mean",           "Fitted Mean per Slice");
      MHCamera dispped6  (geomcam, "Ped;Sigma",          "Fitted Sigma per Slice");
      MHCamera dispped9  (geomcam, "Ped;DeltaPedMean",   "Rel. Diff. Mean per Slice (Fit-Calc.)");
      MHCamera dispped11 (geomcam, "Ped;DeltaRmsSigma",  "Rel. Diff. RMS per Slice (Fit-Calc.)");
  
      dispped0.SetCamContent(  pedcam, 0);
      dispped0.SetCamError(    pedcam, 1);
      dispped2.SetCamContent(  pedcam, 2);
      dispped2.SetCamError(    pedcam, 3);
      
      dispped4.SetCamContent( hpedcam, 0);
      dispped4.SetCamError(   hpedcam, 1);
      dispped6.SetCamContent( hpedcam, 2);
      dispped6.SetCamError(   hpedcam, 3);
      dispped9.SetCamContent( hpedcam, 5);
      dispped9.SetCamError(   hpedcam, 6);
      dispped11.SetCamContent(hpedcam, 8);
      dispped11.SetCamError(  hpedcam, 9);
  
      dispped0.SetYTitle("Calc. Pedestal per slice [FADC counts]");
      dispped2.SetYTitle("Calc. Pedestal RMS per slice [FADC counts]");
      dispped4.SetYTitle("Fitted Mean per slice [FADC counts]");
      dispped6.SetYTitle("Fitted Sigma per slice [FADC counts]");
      dispped9.SetYTitle("Rel. Diff. Pedestal per slice Fit-Calc [1]");
      dispped11.SetYTitle("Rel. Diff. Pedestal RMS per slice Fit-Calc [1]");


      // Histogram values
      tabname="MeanRMS ";
      tabname += samples;
      //fill.SetNameTab(tabname);
      TCanvas &b1 = display->AddTab(tabname);
      b1.Divide(4,3);

      CamDraw(b1,dispped0,1,4,*hmeaninn,*hmeanout,*hmeaninnerr,*hmeanouterr,samples,stepsize);
      CamDraw(b1,dispped2,2,4,*hrmsinn,*hrmsout,*hrmsinnerr,*hrmsouterr,samples,stepsize);  
      CamDraw(b1,dispped4,3,4,*hmuinn,*hmuout,*hmuinnerr,*hmuouterr,samples,stepsize);
      CamDraw(b1,dispped6,4,4,*hsigmainn,*hsigmaout,*hsigmainnerr,*hsigmaouterr,samples,stepsize); 
      
      tabname="MeanRrmSamples ";
      tabname += samples;
      display->SaveAsGIF(3*((samples-1)/stepsize)+2,tabname);

      // Differences
      tabname="RelDiff ";
      tabname += samples;
      TCanvas &c4 = display->AddTab(tabname);
      c4.Divide(2,3);
      
      CamDraw(c4,dispped9,1,2,*hmeandiffinn,*hmeandiffout,*hmeandiffinnerr,*hmeandiffouterr,samples,stepsize);
      CamDraw(c4,dispped11,2,2,*hrmsdiffinn,*hrmsdiffout,*hrmsdiffinnerr,*hrmsdiffouterr,samples,stepsize); 

      tabname="RelDiffSamples ";
      tabname += samples;
      display->SaveAsGIF(3*((samples-1)/stepsize)+3,tabname);

   cout << "Loop N. " <<samples << " completed..." <<endl;

    }

  /*
  TF1 *logg = new TF1("logg","[1]+TMath::Log(x-[0])",1.,30.,2);
  logg->SetParameters(1.,3.5);
  logg->SetParLimits(0,-1.,3.);
  logg->SetParLimits(1,-1.,7.);
  logg->SetLineColor(kRed);
  */

  TCanvas *canvas = new TCanvas("PedstudInner","Pedestal Studies Inner Pixels",600,900);
  canvas->Divide(2,3);
  canvas->cd(1);

  TGraphErrors *gmeaninn = new TGraphErrors(hmeaninn->GetSize(),
                                            CreateXaxis(hmeaninn->GetSize(),stepsize),hmeaninn->GetArray(),
                                            CreateXaxisErr(hmeaninnerr->GetSize(),stepsize),hmeaninnerr->GetArray());
  gmeaninn->Draw("A*");
  gmeaninn->SetTitle("Calculated Mean per Slice Inner Pixels");
  gmeaninn->GetXaxis()->SetTitle("Nr. added FADC slices");
  gmeaninn->GetYaxis()->SetTitle("Calculated Mean per slice");
  //  gmeaninn->Fit("pol0");
  //  gmeaninn->GetFunction("pol0")->SetLineColor(kGreen);
  //  //  gmeaninn->Fit(logg);

  canvas->cd(2);

  TGraphErrors *gmuinn = new TGraphErrors(hmuinn->GetSize(),
                                          CreateXaxis(hmuinn->GetSize(),stepsize),hmuinn->GetArray(),
                                          CreateXaxisErr(hmuinnerr->GetSize(),stepsize),hmuinnerr->GetArray());
  gmuinn->Draw("A*");
  gmuinn->SetTitle("Fitted Mean per Slice Inner Pixels");
  gmuinn->GetXaxis()->SetTitle("Nr. added FADC slices");
  gmuinn->GetYaxis()->SetTitle("Fitted Mean per Slice");
  //  gmuinn->Fit("pol0");
  //  gmuinn->GetFunction("pol0")->SetLineColor(kGreen);
  //gmuinn->Fit(logg);


  canvas->cd(3);

  TGraphErrors *grmsinn = new TGraphErrors(hrmsinn->GetSize(),
                                           CreateXaxis(hrmsinn->GetSize(),stepsize),hrmsinn->GetArray(),
                                           CreateXaxisErr(hrmsinnerr->GetSize(),stepsize),hrmsinnerr->GetArray());
  grmsinn->Draw("A*");
  grmsinn->SetTitle("Calculated Rms per Slice Inner Pixels");
  grmsinn->GetXaxis()->SetTitle("Nr. added FADC slices");
  grmsinn->GetYaxis()->SetTitle("Calculated Rms per Slice");
  //  //grmsinn->Fit("pol2");
  //  //grmsinn->GetFunction("pol2")->SetLineColor(kRed);
  //  grmsinn->Fit(logg);

  canvas->cd(4);

  TGraphErrors *gsigmainn = new TGraphErrors(hsigmainn->GetSize(),
                                             CreateXaxis(hsigmainn->GetSize(),stepsize),hsigmainn->GetArray(),
                                             CreateXaxisErr(hsigmainnerr->GetSize(),stepsize),hsigmainnerr->GetArray());
  gsigmainn->Draw("A*");
  gsigmainn->SetTitle("Fitted Sigma per Slice Inner Pixels");
  gsigmainn->GetXaxis()->SetTitle("Nr. added FADC slices");
  gsigmainn->GetYaxis()->SetTitle("Fitted Sigma per Slice");
  //  //  gsigmainn->Fit("pol2");
  //  //  gsigmainn->GetFunction("pol2")->SetLineColor(kRed);
  //  gsigmainn->Fit(logg);

  canvas->cd(5);

  TGraphErrors *gmeandiffinn = new TGraphErrors(hmeandiffinn->GetSize(),
                                                CreateXaxis(hmeandiffinn->GetSize(),stepsize),hmeandiffinn->GetArray(),
                                                CreateXaxisErr(hmeandiffinnerr->GetSize(),stepsize),hmeandiffinnerr->GetArray());
  gmeandiffinn->Draw("A*"); 
  gmeandiffinn->SetTitle("Rel. Difference  Mean per Slice Inner Pixels");
  gmeandiffinn->GetXaxis()->SetTitle("Nr. added FADC slices");
  gmeandiffinn->GetYaxis()->SetTitle("Rel. Difference Mean per Slice");
  //  //gmeandiffinn->Fit("pol2");
  //  //gmeandiffinn->GetFunction("pol2")->SetLineColor(kBlue);
  //  gmeandiffinn->Fit(logg);


  canvas->cd(6);

  TGraphErrors *grmsdiffinn = new TGraphErrors(hrmsdiffinn->GetSize(),
                                               CreateXaxis(hrmsdiffinn->GetSize(),stepsize),hrmsdiffinn->GetArray(),
                                               CreateXaxisErr(hrmsdiffinnerr->GetSize(),stepsize),hrmsdiffinnerr->GetArray());
  grmsdiffinn->Draw("A*");
  grmsdiffinn->SetTitle("Rel. Difference Sigma per Slice-RMS Inner Pixels");
  grmsdiffinn->GetXaxis()->SetTitle("Nr. added FADC slices");
  grmsdiffinn->GetYaxis()->SetTitle("Rel. Difference Sigma per Slice-RMS");
  //  //grmsdiffinn->Fit("pol2");
  //  //grmsdiffinn->GetFunction("pol2")->SetLineColor(kBlue);
  //  grmsdiffinn->Fit(logg);

  canvas->SaveAs("PedestalStudyInner.root");
  canvas->SaveAs("PedestalStudyInner.ps");

  TCanvas *canvas2 = new TCanvas("PedstudOut","Pedestal Studies Outer Pixels",600,900);
  canvas2->Divide(2,3);
  canvas2->cd(1);

  TGraphErrors *gmeanout = new TGraphErrors(hmeanout->GetSize(),
                                            CreateXaxis(hmeanout->GetSize(),stepsize),hmeanout->GetArray(),
                                            CreateXaxisErr(hmeanouterr->GetSize(),stepsize),hmeanouterr->GetArray());
  gmeanout->Draw("A*");
  gmeanout->SetTitle("Calculated Mean per Slice Outer Pixels");
  gmeanout->GetXaxis()->SetTitle("Nr. added FADC slices");
  gmeanout->GetYaxis()->SetTitle("Calculated Mean per Slice");
  //  gmeanout->Fit("pol0");
  //  gmeanout->GetFunction("pol0")->SetLineColor(kGreen);
  //gmeanout->Fit(logg);

  canvas2->cd(2);

  TGraphErrors *gmuout = new TGraphErrors(hmuout->GetSize(),
                                          CreateXaxis(hmuout->GetSize(),stepsize),hmuout->GetArray(),
                                          CreateXaxisErr(hmuouterr->GetSize(),stepsize),hmuouterr->GetArray());
  gmuout->Draw("A*");
  gmuout->SetTitle("Fitted Mean per Slice Outer Pixels");
  gmuout->GetXaxis()->SetTitle("Nr. added FADC slices");
  gmuout->GetYaxis()->SetTitle("Fitted Mean per Slice");
  //  gmuout->Fit("pol0");
  //  gmuout->GetFunction("pol0")->SetLineColor(kGreen);
  //gmuout->Fit(logg);

  canvas2->cd(3);

  TGraphErrors *grmsout = new TGraphErrors(hrmsout->GetSize(),
                                           CreateXaxis(hrmsout->GetSize(),stepsize),hrmsout->GetArray(),
                                           CreateXaxisErr(hrmsouterr->GetSize(),stepsize),hrmsouterr->GetArray());
  grmsout->Draw("A*");
  grmsout->SetTitle("Calculated Rms per Slice Outer Pixels");
  grmsout->GetXaxis()->SetTitle("Nr. added FADC slices");
  grmsout->GetYaxis()->SetTitle("Calculated Rms per Slice");
  //  //grmsout->Fit("pol2");
  //  //grmsout->GetFunction("pol2")->SetLineColor(kRed);
  //  grmsout->Fit(logg);

  canvas2->cd(4);

  TGraphErrors *gsigmaout = new TGraphErrors(hsigmaout->GetSize(),
                                             CreateXaxis(hsigmaout->GetSize(),stepsize),hsigmaout->GetArray(),
                                             CreateXaxisErr(hsigmaouterr->GetSize(),stepsize),hsigmaouterr->GetArray());
  gsigmaout->Draw("A*");
  gsigmaout->SetTitle("Fitted Sigma per Slice Outer Pixels");
  gsigmaout->GetXaxis()->SetTitle("Nr. added FADC slices");
  gsigmaout->GetYaxis()->SetTitle("Fitted Sigma per Slice");
  //  //gsigmaout->Fit("pol2");
  //  //gsigmaout->GetFunction("pol2")->SetLineColor(kRed);
  //  gsigmaout->Fit(logg);


  canvas2->cd(5);

  TGraphErrors *gmeandiffout = new TGraphErrors(hmeandiffout->GetSize(),
                                                CreateXaxis(hmeandiffout->GetSize(),stepsize),hmeandiffout->GetArray(),
                                                CreateXaxisErr(hmeandiffouterr->GetSize(),stepsize),hmeandiffouterr->GetArray());
  gmeandiffout->Draw("A*");
  gmeandiffout->SetTitle("Rel. Difference  Mean per Slice Outer Pixels");
  gmeandiffout->GetXaxis()->SetTitle("Nr. added FADC slices");
  gmeandiffout->GetYaxis()->SetTitle("Rel. Difference Mean per Slice");
  //  //gmeandiffout->Fit("pol2");
  //w  //gmeandiffout->GetFunction("pol2")->SetLineColor(kBlue);
  //  gmeandiffout->Fit(logg);

  canvas2->cd(6);

  TGraphErrors *grmsdiffout = new TGraphErrors(hrmsdiffout->GetSize(),
                                               CreateXaxis(hrmsdiffout->GetSize(),stepsize),hrmsdiffout->GetArray(),
                                               CreateXaxisErr(hrmsdiffouterr->GetSize(),stepsize),hrmsdiffouterr->GetArray());
  grmsdiffout->Draw("A*");
  grmsdiffout->SetTitle("Rel. Difference Sigma per Slice-RMS Outer Pixels");
  grmsdiffout->GetXaxis()->SetTitle("Nr. added FADC slices");
  grmsdiffout->GetYaxis()->SetTitle("Rel. Difference Sigma per Slice-RMS");
  //  //grmsdiffout->Fit("pol2");
  //  //grmsdiffout->GetFunction("pol2")->SetLineColor(kBlue);
  //  grmsdiffout->Fit(logg);


  canvas2->SaveAs("PedestalStudyOuter.root");
  canvas2->SaveAs("PedestalStudyOuter.ps");


}


void CamDraw(TCanvas &c, MHCamera &cam, Int_t i, Int_t j, TArrayF &a1, TArrayF &a2, 
             TArrayF &a1err, TArrayF &a2err, Int_t samp, Int_t stepsize)
{

  c.cd(i);
  MHCamera *obj1=(MHCamera*)cam.DrawCopy("hist");
  obj1->SetDirectory(NULL);
  
  c.cd(i+j);
  obj1->Draw();
  ((MHCamera*)obj1)->SetPrettyPalette();

  c.cd(i+2*j);
  TH1D *obj2 = (TH1D*)obj1->Projection();
  obj2->SetDirectory(NULL);
  
  //      obj2->Sumw2();
  obj2->Draw();
  obj2->SetBit(kCanDelete);
  
  const Double_t min   = obj2->GetBinCenter(obj2->GetXaxis()->GetFirst());
  const Double_t max   = obj2->GetBinCenter(obj2->GetXaxis()->GetLast());
  const Double_t integ = obj2->Integral("width")/2.5066283;
  const Double_t mean  = obj2->GetMean();
  const Double_t rms   = obj2->GetRMS();
  const Double_t width = max-min;
  
  if (rms == 0. || width == 0. )
    return;
  
  TArrayI s0(6);
  s0[0] = 6;
  s0[1] = 1;
  s0[2] = 2;
  s0[3] = 3;
  s0[4] = 4;
  s0[5] = 5;
  
  TArrayI inner(1);
  inner[0] = 0;
  
  TArrayI outer(1);
  outer[0] = 1;
      
  // Just to get the right (maximum) binning
  TH1D *half[2];
  half[0] = obj1->ProjectionS(s0, inner, "Inner");
  half[1] = obj1->ProjectionS(s0, outer, "Outer");

  half[0]->SetDirectory(NULL);
  half[1]->SetDirectory(NULL);
  
  for (int i=0; i<2; i++)
    {
      half[i]->SetLineColor(kRed+i);
      half[i]->SetDirectory(0);
      half[i]->SetBit(kCanDelete);
      half[i]->Draw("same");
      half[i]->Fit("gaus","Q+");

      if (i==0)
        {
          a1[(samp-1)/stepsize] = half[i]->GetFunction("gaus")->GetParameter(1);
          a1err[(samp-1)/stepsize] = half[i]->GetFunction("gaus")->GetParError(1);
          if (a1err[(samp-1)/stepsize] > 3.)
            a1err[(samp-1)/stepsize] = 1.;
        }
     if (i==1)
       {
         a2[(samp-1)/stepsize] = half[i]->GetFunction("gaus")->GetParameter(1);
         a2err[(samp-1)/stepsize] = half[i]->GetFunction("gaus")->GetParError(1);
          if (a2err[(samp-1)/stepsize] > 3.)
            a2err[(samp-1)/stepsize] = 1.;
       }
    }
  
  
}

// -----------------------------------------------------------------------------
// 
// Create the x-axis for the event graph
//
Float_t *CreateXaxis(Int_t n, Int_t step)
{

  Float_t *xaxis = new Float_t[n];

  for (Int_t i=0;i<n;i++)
    xaxis[i] = 2. + step*i;

  return xaxis;
                 
}

// -----------------------------------------------------------------------------
// 
// Create the x-axis for the event graph
//
Float_t *CreateXaxisErr(Int_t n, Int_t step)
{

  Float_t *xaxis = new Float_t[n];

  for (Int_t i=0;i<n;i++)
    xaxis[i] = step/2.;

  return xaxis;
                 
}