Index: trunk/MagicSoft/Mars/mtemp/mifae/programs/pedvsslices.cc
===================================================================
--- trunk/MagicSoft/Mars/mtemp/mifae/programs/pedvsslices.cc	(revision 4016)
+++ trunk/MagicSoft/Mars/mtemp/mifae/programs/pedvsslices.cc	(revision 4016)
@@ -0,0 +1,253 @@
+#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");
+
+}
