Index: trunk/MagicSoft/Mars/macros/pointspreadfunction.C
===================================================================
--- trunk/MagicSoft/Mars/macros/pointspreadfunction.C	(revision 3110)
+++ trunk/MagicSoft/Mars/macros/pointspreadfunction.C	(revision 3110)
@@ -0,0 +1,824 @@
+/* ======================================================================== *\
+!
+! *
+! * 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): Javier Lopez, 12/2003 <mailto:jlopez@ifae.es>
+!   Author(s): Alex Armada,   1/2004 <mailto:armada@ifae.es>
+!
+!   Copyright: MAGIC Software Development, 2000-2003
+!
+!
+\* ======================================================================== */
+
+
+//------------------------------------------------------------------------- //
+//                                                                          //
+//  This macro fits the dc signal of a star using a two dimension gaussian  //
+//  for each dc measurement. Then the values of parameters of the fit are   //
+//  stored in histograms and shown at the end of the macro.                 //
+//                                                                          //
+//  USAGE:                                                                  //
+//  It has two arguments,                                                   //
+//  1- The first one is the dc file with the tracked star                   //
+//  2- The second one is a continuos light file used to intercalibrate      //
+//     the gain of the photomultipliers.                                    //
+//     (It's possible not to use the calibration file and then the gain     //
+//      of the pmts are supouse to be the same for all of them.             //
+//                                                                          //
+//--------------------------------------------------------------------------// 
+
+const Int_t numPixels = 577;
+Int_t nPixelsFitted; // 396[inners] 36[3rings] 60[4rings] 90[5rings]
+Bool_t isPixelsFitted[numPixels];
+Float_t z[numPixels],x[numPixels],y[numPixels],errorz[numPixels];
+Float_t chisquare;
+
+//______________________________________________________________________________
+void fcn(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)
+{
+    Int_t i;
+
+//calculate chisquare
+    Double_t chisq = 0;
+    Double_t delta;
+    nPixelsFitted=0;
+    for (i=1;i<numPixels; i++) {
+	if (isPixelsFitted[i])
+	{
+	    if (errorz[i] != 0.0)
+	    {
+		delta  = (z[i]-func(x[i],y[i],par))/errorz[i];
+		chisq += delta*delta;
+		nPixelsFitted++;
+	    }
+	    else
+		cout << "This should never happen errorz[" << i << "] " << errorz[i] << endl;
+	}
+    }
+    f = chisq;
+    chisquare = chisq;
+}
+
+//______________________________________________________________________________
+Double_t func(float x,float y,Double_t *par)
+{
+    Double_t value=par[0]*TMath::exp(-(x-par[1])*(x-par[1])/(2*par[2]*par[2]))*TMath::exp(-(y-par[3])*(y-par[3])/(2*par[4]*par[4]));
+    return value;
+}
+
+Bool_t HandleInput()
+{
+    TTimer timer("gSystem->ProcessEvents();", 50, kFALSE);
+    
+    while (1)
+    {
+        //
+        // While reading the input process gui events asynchronously
+        //
+        timer.TurnOn();
+        TString input = Getline("Type 'q' to exit, <return> to go on: ");
+        timer.TurnOff();
+	
+        if (input=="q\n")
+            return kFALSE;
+	
+        if (input=="\n")
+            return kTRUE;
+    };
+    
+    return kFALSE;
+}
+
+void pointspreadfunction(TString fname, TString clname = "NULL", Int_t userNumLines = 1000)
+{
+
+// First run over continuos light files to have a DC calibration
+
+
+    Float_t currmean[numPixels];
+    Float_t currsquaremean[numPixels];
+    Float_t currrms[numPixels];
+    Float_t meanofcurrmean = 0;
+
+    for (Int_t swpix=0; swpix<numPixels; swpix++)
+    {
+ 	currmean[swpix] = 0.;
+ 	currsquaremean[swpix] = 0.;
+  	currrms[swpix] = 0.;
+    }
+    
+    Int_t numLines=0;
+
+    //containers
+    MGeomCamMagic geomcam;
+    MCameraDC     curr;
+    MCameraDC     currbis;
+    
+    const Float_t conv4mm2deg = geomcam.GetConvMm2Deg();
+	 
+    if (clname != "NULL")
+    {
+	MParList plist0;
+	
+	MTaskList     tlist0;
+	plist0.AddToList(&tlist0);
+	
+	plist0.AddToList(&geomcam);
+	plist0.AddToList(&curr);
+	
+	MReportFileRead read0(clname);
+	read0.SetHasNoHeader();
+	read0.AddToList("MReportCurrents");
+	
+	tlist0.AddToList(&read0);
+	
+	MEvtLoop evtloop0;
+	evtloop0.SetParList(&plist0);
+	
+	
+	if (!evtloop0.PreProcess())
+	    return;
+	
+	while (tlist0.Process())
+	{
+	    for (Int_t swpix=0; swpix<numPixels; swpix++)
+	    {
+		meanofcurrmean += curr[swpix];
+		currmean[swpix] += curr[swpix];
+		currsquaremean[swpix] += curr[swpix]*curr[swpix];
+	    }
+	    numLines++;
+	}
+	
+	evtloop0.PostProcess();
+	
+	meanofcurrmean /= (numLines*numPixels);
+	for (Int_t swpix=0; swpix<numPixels; swpix++)
+	{
+	    
+	    currmean[swpix] /= numLines;
+	    currsquaremean[swpix] /= numLines;
+	    currrms[swpix] = sqrt(fabs(currsquaremean[swpix] - currmean[swpix]*currmean[swpix]));
+	    
+	    curr[swpix] = currmean[swpix];
+	    currbis[swpix] = currrms[swpix];
+	    
+	    currmean[swpix] /= meanofcurrmean;	
+	    currrms[swpix] /= meanofcurrmean;	
+	    
+// 	cout << "swpixel " << swpix << " dc fCalib " << currmean[swpix] << " +- " << currrms[swpix] << endl;
+	}
+	
+	
+	//    curr.Print();
+/*     MHCamera display0(geomcam);
+       display0.SetPrettyPalette();
+       display0.Draw();
+       
+       display0.SetCamContent(currbis);
+       cout << "PSF>> DC mean values drawn" << endl;
+       // Remove the comments if you want to go through the file
+       // event-by-event:
+       if (!HandleInput())
+       break;
+*/
+    }
+    else
+    {
+	for (Int_t swpix=0; swpix<numPixels; swpix++)
+	{
+	    currmean[swpix] = 1.;
+	    currrms[swpix] = 0.2;
+	}
+	
+    }
+
+// Now we can run over the dc data to extract the psf rms. 
+	
+    const Int_t maxNumLines = 10000;
+	
+    Double_t ux[maxNumLines];
+    Double_t uy[maxNumLines];
+    Double_t sx[maxNumLines];
+    Double_t sy[maxNumLines];
+    Double_t chisqu[maxNumLines];
+    Double_t time[maxNumLines];
+
+    MParList plist;
+
+    MGeomCamMagic geomcam;
+    MCameraDC     curr;
+    MTaskList     tlist;
+
+    plist.AddToList(&geomcam);
+    plist.AddToList(&curr);
+    plist.AddToList(&tlist);
+
+    MReportFileRead read(fname);
+    read.SetHasNoHeader();
+    read.AddToList("MReportCurrents");
+
+    tlist.AddToList(&read);
+    
+    MEvtLoop evtloop;
+    evtloop.SetParList(&plist);
+    
+    if (!evtloop.PreProcess())
+        return;
+    
+    MHCamera display(geomcam);
+    display.SetPrettyPalette();
+//    display.SetInvDeepBlueSeaPalette()
+    display.Draw();
+    gPad->SetLogy();
+    gPad->cd(1);
+
+    Double_t amin,edm,errdef;
+    Int_t nvpar,nparx,icstat;
+
+    Double_t max,maxerror;
+    Double_t xmean,xsigma,ymean,ysigma;
+    Double_t xmeanerror,xsigmaerror,ymeanerror,ysigmaerror;
+
+//    TVector xellipsecenter;
+//    TVector yellipsecenter;
+
+    TEllipse ellipse;
+    ellipse.SetFillStyle(4000);
+    ellipse.SetLineWidth(2);
+    ellipse.SetLineColor(2);
+
+    ellipse.Draw();
+
+    Int_t nbinsxy = 80;
+    Float_t minxy = -600*conv4mm2deg;
+    Float_t maxxy = 600*conv4mm2deg;
+    Float_t fromdegtobin = (maxxy-minxy)/nbinsxy;
+
+    TH2D psfhist("psfhist","",nbinsxy,minxy,maxxy,nbinsxy,minxy,maxxy);
+    psfhist->GetXaxis()->SetTitle("[deg]");
+    psfhist->GetYaxis()->SetTitle("[deg]");
+    psfhist->GetZaxis()->SetTitle("DC [uA]");
+
+    TCanvas *psfcanvas = new TCanvas("psfcanvas","Point Spread Funtion 2D",200,20,900,700);
+
+    tlist.Process();
+
+    Int_t numLines=0;
+    Float_t minDCStar = 6.0;
+
+    Int_t numPixelsInStar = 0;
+    Float_t maxDC = 0;
+    Int_t swpixelmaxDC;
+    
+    Bool_t isPixelsFittedTmp[numPixels];
+    
+    for(Int_t swpixel=1; swpixel<numPixels; swpixel++)
+	isPixelsFittedTmp[swpixel] = kFALSE;
+    
+    for(Int_t swpixel=1; swpixel<numPixels; swpixel++)
+    {
+	if(curr[swpixel] > maxDC)
+	{
+	    swpixelmaxDC = swpixel;
+	    maxDC = curr[swpixelmaxDC];
+	}
+	
+	if(curr[swpixel]>minDCStar)
+	{
+	    numPixelsInStar++;
+	    isPixelsFitted[swpixel] = kTRUE;
+	}
+	else
+	    isPixelsFitted[swpixel] = kFALSE;
+    }
+    
+    if (numPixelsInStar == 0)
+    {
+	cout << "PSF>> Warning: none pixel over minDCStar(" << minDCStar  << ')' << endl; 
+	return;
+    }
+    
+//1st neighboor ring
+    for(Int_t swpixel=1; swpixel<numPixels; swpixel++)
+	if (isPixelsFitted[swpixel])
+	    for(Int_t next=0; next<geomcam[swpixel].GetNumNeighbors(); next++)
+		isPixelsFittedTmp[geomcam[swpixel].GetNeighbor(next)] = kTRUE;
+    
+    for(Int_t swpixel=1; swpixel<numPixels; swpixel++)
+	if (isPixelsFittedTmp[swpixel])
+	    isPixelsFitted[swpixel] = kTRUE;
+    
+//2on neighboor ring
+    for(Int_t swpixel=1; swpixel<numPixels; swpixel++)
+	if (isPixelsFitted[swpixel])
+	    for(Int_t next=0; next<geomcam[swpixel].GetNumNeighbors(); next++)
+		isPixelsFittedTmp[geomcam[swpixel].GetNeighbor(next)] = kTRUE;
+    
+    for(Int_t swpixel=1; swpixel<numPixels; swpixel++)
+	if (isPixelsFittedTmp[swpixel])
+	    isPixelsFitted[swpixel] = kTRUE;
+    
+
+//3rt neighboor ring
+    for(Int_t swpixel=1; swpixel<numPixels; swpixel++)
+	if (isPixelsFitted[swpixel])
+	    for(Int_t next=0; next<geomcam[swpixel].GetNumNeighbors(); next++)
+		isPixelsFittedTmp[geomcam[swpixel].GetNeighbor(next)] = kTRUE;
+    
+    for(Int_t swpixel=1; swpixel<numPixels; swpixel++)
+	if (isPixelsFittedTmp[swpixel])
+	    isPixelsFitted[swpixel] = kTRUE;
+	      
+   
+    for(Int_t swpixel=1; swpixel<numPixels; swpixel++)
+    {
+//		cout << "isPixelsFitted[" << swpixel << "] ";
+//		cout << (Int_t)isPixelsFitted[swpixel] << endl;
+	curr[swpixel] = (Float_t)isPixelsFitted[swpixel]; 
+    }
+    
+/*    MHCamera display0(geomcam);
+    display0.SetPrettyPalette();
+    display0.Draw();
+    
+    display0.SetCamContent(curr);
+    cout << "PSF>> Fitted pixels drawn" << endl;
+    // Remove the comments if you want to go through the file
+    // event-by-event:
+    if (!HandleInput())
+	break;
+*/
+    TMinuit *gMinuit = new TMinuit(7);  //initialize TMinuit with a maximum of 5 params
+    gMinuit->SetFCN(fcn);
+
+    Double_t arglist[10];
+    Int_t ierflg = 0;
+
+    arglist[0] = 1;
+    gMinuit->mnexcm("SET ERR", arglist ,1,ierflg);
+//    arglist[0] = -1;
+    arglist[0] = 0;
+    gMinuit->mnexcm("SET PRI", arglist ,1,ierflg);
+
+// Set starting values and step sizes for parameters
+    Double_t vstart[5];
+    Double_t step[5];
+    Double_t lowlimit[5] = {minDCStar, -2., 0.05, -2, 0.05};
+    Double_t uplimit[5] = {30., 2., 1., 2., 1.};
+
+    vstart[0] = maxDC;
+    vstart[1] = geomcam[swpixelmaxDC].GetX()*conv4mm2deg;
+    vstart[2] = 30*TMath::Sqrt(numPixelsInStar/2)*conv4mm2deg;
+    vstart[3] = geomcam[swpixelmaxDC].GetY()*conv4mm2deg;
+    vstart[4] = 30*TMath::Sqrt(numPixelsInStar/2)*conv4mm2deg;
+    
+    for(Int_t i=0; i<5; i++)
+    {
+	if (vstart[i] != 0)
+	    step[i] = TMath::Abs(vstart[i]/sqrt(2));
+	else
+	    step[i] = 1;
+
+	cout << "DGB vstat[" << i << "] " << vstart[i] << " step[" << i << "] " << step[i] <<endl;
+    }
+//    cout << endl;
+  
+    gMinuit->mnparm(0, "max", vstart[0], step[0], lowlimit[0], uplimit[0],ierflg);
+    gMinuit->mnparm(1, "xmean", vstart[1], step[1], lowlimit[1], uplimit[1],ierflg);
+    gMinuit->mnparm(2, "xsigma", vstart[2], step[2], lowlimit[2], uplimit[2],ierflg);
+    gMinuit->mnparm(3, "ymean", vstart[3], step[3], lowlimit[3], uplimit[3],ierflg);
+    gMinuit->mnparm(4, "ysigma", vstart[4], step[4], lowlimit[4], uplimit[4],ierflg);
+
+    numLines++;
+
+    while (tlist.Process() && numLines < maxNumLines)
+	{
+//        curr.Print();
+
+	for (Int_t swpixel=1; swpixel<577; swpixel++)
+	{
+//	    cout << swpixel << '\t' << curr[swpixel] << '\t' << geomcam[swpixel].GetX() << '\t' << geomcam[swpixel].GetY() << endl;
+	    
+	    x[swpixel] = geomcam[swpixel].GetX()*conv4mm2deg;
+	    y[swpixel] = geomcam[swpixel].GetY()*conv4mm2deg;
+	    z[swpixel] = curr[swpixel]/currmean[swpixel];
+	    errorz[swpixel] = TMath::Sqrt((curr[swpixel]*currrms[swpixel]/(currmean[swpixel]*currmean[swpixel]))*(curr[swpixel]*currrms[swpixel]/(currmean[swpixel]*currmean[swpixel]))+(0.1)/currmean[swpixel]*(0.1)/currmean[swpixel]);
+
+
+	    psfhist->SetBinContent((Int_t)((x[swpixel]+600*conv4mm2deg)/fromdegtobin),(Int_t)((y[swpixel]+600*conv4mm2deg)/fromdegtobin),z[swpixel]);
+	}
+	
+	psfcanvas->cd(1);
+	psfhist->Draw("lego2");
+	
+// Now ready for minimization step
+	arglist[0] = 500;
+	arglist[1] = 1.;
+	gMinuit->mnexcm("MIGRAD", arglist ,2,ierflg);
+
+// Print results
+/*	gMinuit->mnstat(amin,edm,errdef,nvpar,nparx,icstat);
+	gMinuit->mnprin(3,amin);
+*/
+	gMinuit->GetParameter(0,max,maxerror);
+	gMinuit->GetParameter(1,xmean,xmeanerror);
+	gMinuit->GetParameter(2,xsigma,xsigmaerror);
+	gMinuit->GetParameter(3,ymean,ymeanerror);
+	gMinuit->GetParameter(4,ysigma,ysigmaerror);
+
+/*	cout << endl;
+	cout << "numLine " << numLines << endl; 
+	cout << "max \t" << max <<  " +- " << maxerror << endl;
+	cout << "xmean \t" << xmean << " +- " << xmeanerror << endl;
+	cout << "xsigma \t" << TMath::Abs(xsigma) << " +- " << xsigmaerror << endl;
+	cout << "ymean \t" << ymean << " +- " << ymeanerror << endl;
+	cout << "ysigma \t" << TMath::Abs(ysigma) << " +- " << ysigmaerror << endl;
+	cout << "chisquare/ndof \t" << chisquare/(nPixelsFitted-5) << endl;
+*/
+	
+	chisqu[numLines]=chisquare/(nPixelsFitted-5);
+
+        if(chisqu[numLines]<100.)
+	{
+	    ux[numLines]=xmean;
+	    uy[numLines]=ymean;
+	    sx[numLines]=TMath::Abs(xsigma);
+	    sy[numLines]=TMath::Abs(ysigma);
+	    time[numLines]=numLines;
+
+	    display.SetCamContent(curr);
+	    gPad->cd(1);
+//	ellipse.DrawEllipse(xmean,ymean,TMath::Abs(xsigma),TMath::Abs(ysigma),0,360,0);
+	    ellipse.SetX1(xmean/conv4mm2deg);
+	    ellipse.SetY1(ymean/conv4mm2deg);
+	    ellipse.SetR1(TMath::Abs(xsigma)/conv4mm2deg);
+	    ellipse.SetR2(TMath::Abs(ysigma)/conv4mm2deg);
+
+	    gPad->Modified();
+	    gPad->Update();
+
+/*	if (numLines%10==0)
+	{
+	    char imagepath[60];
+	    sprintf(&imagepath,"/home/Javi/mnt_users/Data/CaCo/img/ZetaTau_%04d.gif",numLines/10);
+	    cout << imagepath << endl;
+	    gPad->SaveAs(imagepath);
+	}
+*/
+	    numLines++;
+	    // Remove the comments if you want to go through the file
+	    //event-by-event:
+	    if (userNumLines>0)
+	    {
+		if (numLines>userNumLines)
+		    break;
+	    }
+	    else
+	    {
+		if (!HandleInput())
+		    break;
+	    }
+	}
+    }
+
+    
+    evtloop.PostProcess();
+
+//Draw the ditributions of the sigmas the point spread function
+
+    cout<<"Number of lines "<<numLines<<endl;
+
+    gROOT->Reset();
+    gStyle->SetCanvasColor(0);
+    gStyle->SetCanvasBorderMode(0);
+    gStyle->SetPadBorderMode(0);
+    gStyle->SetFrameBorderMode(0);
+    gStyle->SetOptStat(00000000);
+
+    c1 = new TCanvas("c1","Time evolution & distributions",0,0,1200,850);
+    c1->Divide(3,2);
+
+    c1->cd(1);
+
+    TMath math;
+   
+    Double_t minmeanx, maxmeanx;
+    minmeanx = ux[math.LocMin(numLines,ux)];
+    maxmeanx = ux[math.LocMax(numLines,ux)];
+
+    Double_t minmeany, maxmeany;
+    minmeany = uy[math.LocMin(numLines,uy)];
+    maxmeany = uy[math.LocMax(numLines,uy)];
+
+    Double_t minmean, maxmean;
+    minmean = math.Min(minmeanx,minmeany);
+    maxmean = math.Max(maxmeanx,maxmeany);
+
+    Double_t diff;
+    diff = maxmean - minmean;
+    diff = 0.1*diff;
+    minmean = minmean - diff;
+    maxmean = maxmean + diff;
+
+    Double_t mintime, maxtime;
+    mintime = time[math.LocMin(numLines,time)];
+    maxtime = time[math.LocMax(numLines,time)];
+
+    TH2D *h1 = new TH2D("h1",fname,1,mintime-1,maxtime+1,1,minmean,maxmean);
+    h1->GetXaxis()->SetTitle("Event number");
+    h1->GetYaxis()->SetTitle("mean position (deg)");
+    h1->Draw();
+
+    TGraph *grtimeevolmeanx = new TGraph(numLines,time,ux);
+    grtimeevolmeanx->SetMarkerColor(3);
+    grtimeevolmeanx->SetMarkerStyle(20);
+    grtimeevolmeanx->SetMarkerSize (0.4);
+    grtimeevolmeanx->Draw("P");
+
+    TGraph *grtimeevolmeany = new TGraph(numLines,time,uy);
+    grtimeevolmeany->SetMarkerColor(6);
+    grtimeevolmeany->SetMarkerStyle(24);
+    grtimeevolmeany->SetMarkerSize (0.4);
+    grtimeevolmeany->Draw("P");
+
+    legmeanxy = new TLegend(0.8,0.85,0.95,0.95);
+    legmeanxy.SetTextSize(0.03);
+    legmeanxy.AddEntry(grtimeevolmeanx,"mean x","P");
+    legmeanxy.AddEntry(grtimeevolmeany,"mean y","P");
+    legmeanxy.Draw();
+
+    c1->cd(2);
+
+    TMath math;
+   
+    Double_t minsigmax, maxsigmax;
+    minsigmax = sx[math.LocMin(numLines,sx)];
+    maxsigmax = sx[math.LocMax(numLines,sx)];
+
+    Double_t minsigmay, maxsigmay;
+    minsigmay = sy[math.LocMin(numLines,sy)];
+    maxsigmay = sy[math.LocMax(numLines,sy)];
+
+    Double_t minsigma, maxsigma;
+    minsigma = math.Min(minsigmax,minsigmay);
+    maxsigma = math.Max(maxsigmax,maxsigmay);
+
+    diff = maxsigma - minsigma;
+    diff = 0.1*diff;
+    minsigma = minsigma - diff;
+    maxsigma = maxsigma + diff;
+
+    TH2D *h2 = new TH2D("h2","",1,mintime-1,maxtime+1,1,minsigma,maxsigma);
+    h2->GetXaxis()->SetTitle("Event number");
+    h2->GetYaxis()->SetTitle("PSF Rms (deg)");
+    h2->Draw();
+
+    TGraph* grtimeevolsigmax= new TGraph(numLines,time,sx);
+    grtimeevolsigmax->SetMarkerColor(3);
+    grtimeevolsigmax->SetMarkerStyle(20);
+    grtimeevolsigmax->SetMarkerSize (0.4);
+    grtimeevolsigmax->Draw("P");
+
+    TGraph* grtimeevolsigmay= new TGraph(numLines,time,sy);
+    grtimeevolsigmay->SetMarkerColor(6);
+    grtimeevolsigmay->SetMarkerStyle(24);
+    grtimeevolsigmay->SetMarkerSize (0.4);
+    grtimeevolsigmay->Draw("P");
+
+    legsigmaxy = new TLegend(0.8,0.85,0.95,0.95);
+    legsigmaxy.SetTextSize(0.03);
+    legsigmaxy.AddEntry(grtimeevolsigmax,"sigma x","P");
+    legsigmaxy.AddEntry(grtimeevolsigmay,"sigma y","P");
+    legsigmaxy.Draw();
+   
+    c1->cd(3);
+
+    Double_t minchisqu, maxchisqu;
+
+    minchisqu = chisqu[math.LocMin(numLines,chisqu)];
+    maxchisqu = chisqu[math.LocMax(numLines,chisqu)];
+
+    diff = maxchisqu - minchisqu;
+    diff = 0.1*diff;
+    minchisqu = minchisqu - diff;
+    maxchisqu = maxchisqu + diff;
+
+    TH2D *h3 = new TH2D("h3","",1,mintime-1,maxtime+1,1,minchisqu,maxchisqu);
+    h3->GetXaxis()->SetTitle("Event number");
+    h3->Draw();
+
+    TGraph * grtimeevolchisqu = new TGraph(numLines,time,chisqu);
+    grtimeevolchisqu->SetMarkerColor(215);
+    grtimeevolchisqu->SetMarkerStyle(20);
+    grtimeevolchisqu->SetMarkerSize(0.4);
+    grtimeevolchisqu->Draw("P");
+
+    legchisqu = new TLegend(0.55,0.90,0.95,0.95);
+    legchisqu.SetTextSize(0.03);
+    legchisqu.AddEntry(grtimeevolchisqu,"chi square / ndof","P");
+    legchisqu.Draw();
+
+//***************************************
+
+    const Int_t nbins = 100; 
+  
+    TH1D *xsigmahist = new TH1D("xsigmahist","",nbins,minsigma,maxsigma);
+    TH1D *ysigmahist = new TH1D("ysigmahist","",nbins,minsigma,maxsigma);
+    TH1D *xmeanhist = new TH1D("xmeanhist","",nbins,minmean,maxmean);
+    TH1D *ymeanhist = new TH1D("ymeanhist","",nbins,minmean,maxmean);
+    TH1D *chisquhist = new TH1D("chisquhist","",nbins,minchisqu,maxchisqu);
+
+    for (Int_t i=0; i<numLines; i++)
+    {
+	xsigmahist->Fill(TMath::Abs(sx[i]));
+	ysigmahist->Fill(TMath::Abs(sy[i]));
+	xmeanhist->Fill(ux[i]);
+	ymeanhist->Fill(uy[i]);
+	chisquhist->Fill(chisqu[i]);
+	
+    }
+
+    c1->cd(5);
+
+//    TCanvas *sigmacanvas = new TCanvas("sigmacanvas","Point Spread Funtion RMS",100,10,800,600);
+
+    TMath math;
+    Double_t maxsigma;
+    Int_t binmaxx, binmaxy;
+    xsigmahist->SetLineColor(3);
+    xsigmahist->SetLineWidth(2);
+    xsigmahist->SetXTitle("RMS [deg]");
+    binmaxx = xsigmahist->GetMaximumBin();
+    binmaxx = xsigmahist->GetBinContent(binmaxx);
+       
+    ysigmahist->SetLineColor(6);
+    ysigmahist->SetLineWidth(2);
+    binmaxy = ysigmahist->GetMaximumBin();
+    binmaxy = ysigmahist->GetBinContent(binmaxy);
+   
+    maxsigma = math.Max(binmaxx,binmaxy);
+    maxsigma += 1;
+
+    xsigmahist->SetMaximum(maxsigma);
+    ysigmahist->SetMaximum(maxsigma);
+    xsigmahist->DrawCopy();
+    ysigmahist->DrawCopy("Same");
+
+    TLegend *legendsigma = new TLegend(.3,.8,.95,.95);
+    legendsigma->SetTextSize(0.03);
+    char xsigmatitle[100];
+    char ysigmatitle[100];
+    sprintf(xsigmatitle,"PSF Rms on X axis -- %5.2f +/- %5.2f deg",xsigmahist->GetMean(),xsigmahist->GetRMS());
+    sprintf(ysigmatitle,"PSF Rms on Y axis -- %5.2f +/- %5.2f deg",ysigmahist->GetMean(),ysigmahist->GetRMS());
+    legendsigma->AddEntry(xsigmahist,xsigmatitle,"F");
+    legendsigma->AddEntry(ysigmahist,ysigmatitle,"F");
+    legendsigma->Draw();
+
+    c1->cd(4);
+   
+    Double_t maxmean;
+
+    xmeanhist->SetLineColor(3);
+    xmeanhist->SetLineWidth(2);
+    xmeanhist->SetXTitle("mean [deg]");
+    binmaxx = xmeanhist->GetMaximumBin();
+    binmaxx = xmeanhist->GetBinContent(binmaxx);
+
+    ymeanhist->SetLineColor(6);
+    ymeanhist->SetLineWidth(2);
+    binmaxy = ymeanhist->GetMaximumBin();
+    binmaxy = ymeanhist->GetBinContent(binmaxy);
+
+    maxmean = math.Max(binmaxx,binmaxy);
+    maxmean += 1;
+
+    xmeanhist->SetMaximum(maxmean);
+    ymeanhist->SetMaximum(maxmean);
+    xmeanhist->DrawCopy();
+    ymeanhist->DrawCopy("Same");
+
+    TLegend *legendmean = new TLegend(.35,.8,.95,.95);
+    legendmean->SetTextSize(0.03);
+    char xmeantitle[100];
+    char ymeantitle[100];
+    sprintf(xmeantitle,"mean on X axis -- %5.2f +/- %5.2f deg",xmeanhist->GetMean(),xmeanhist->GetRMS());
+    sprintf(ymeantitle,"mean on Y axis -- %5.2f +/- %5.2f deg",ymeanhist->GetMean(),ymeanhist->GetRMS());
+    legendmean->AddEntry(xmeanhist,xmeantitle,"F");
+    legendmean->AddEntry(ymeanhist,ymeantitle,"F");
+    legendmean->Draw();
+
+    //meancanvas->Modified();
+    //meancanvas->Update();
+
+    c1->cd(6);
+
+    chisquhist->SetLineColor(215);
+    chisquhist->SetLineWidth(2);
+    chisquhist->SetXTitle("chi square / ndof");
+    TAxis * axis = chisquhist->GetXaxis();
+    axis->SetLabelSize(0.025);
+    chisquhist->DrawCopy();
+    
+    TLegend *legendchisqu = new TLegend(.4,.85,.95,.95);
+    legendchisqu->SetTextSize(0.03);
+    char chisqutitle[100];
+    sprintf(chisqutitle,"chi square / ndof -- %5.2f +/- %5.2f ",chisquhist->GetMean(),chisquhist->GetRMS());
+    legendchisqu->AddEntry(chisquhist,chisqutitle,"F");
+    legendchisqu->Draw();
+    
+    
+    return;
+      
+}
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+