Index: trunk/MagicSoft/Mars/mtemp/mmpi/MSkyPlot.cc
===================================================================
--- trunk/MagicSoft/Mars/mtemp/mmpi/MSkyPlot.cc	(revision 5490)
+++ trunk/MagicSoft/Mars/mtemp/mmpi/MSkyPlot.cc	(revision 6190)
@@ -115,5 +115,4 @@
 #include "MSkyPlot.h"
 
-#include <vector>
 #include <TF1.h>
 #include <TH2.h>
@@ -135,4 +134,5 @@
 #include "MHillasExt.h"
 #include "MNewImagePar.h"
+#include "MHadronness.h"
 #include "MTime.h"
 #include "MObservatory.h"
@@ -151,4 +151,6 @@
 #include "MLogManip.h"
 
+#include <TOrdCollection.h>
+
 ClassImp(MSkyPlot);
 
@@ -160,9 +162,19 @@
 //
 MSkyPlot::MSkyPlot(const char *name, const char *title)
-    : fTime(0), fPointPos(0), fObservatory(0), fMm2Deg(-1) 
-//    fAlphaCut(12.5), BgMean(55), fMinDist(-1), fMaxDist(-1), fMinLD(-1), fMaxLD(-1)
-{
-
-*fLog << warn << "entering default constructor in MSkyPlot" << endl; 
+  : fGeomCam(NULL),   
+    fTime(NULL),      
+    fPointPos(NULL),  
+    fRepDrive(NULL),  
+    fSrcPosCam(NULL), 
+    fPntPosCam(NULL), 
+    fObservatory(NULL),
+    fHillas(NULL),    
+    fHillasExt(NULL), 
+    fHillasSrc(NULL), 
+    fNewImagePar(NULL),
+    fMm2Deg(-1)
+{
+
+  *fLog << warn << "entering default constructor in MSkyPlot" << endl; 
     //
     //   set the name and title of this object
@@ -176,4 +188,10 @@
     fHistNexcess.SetDirectory(NULL);
     fHistOn.SetDirectory(NULL);
+    fHistSignifGaus.SetDirectory(NULL);
+
+    fHistSignif.UseCurrentStyle();
+    fHistNexcess.UseCurrentStyle();
+    fHistOn.UseCurrentStyle();
+    fHistSignifGaus.UseCurrentStyle();
 
     fHistSignif.SetName("SkyPlotSignif");
@@ -210,4 +228,6 @@
     fGridFineBin = 0.01;   // degrees
 
+// elapsed time:
+    fElaTime = 0.;
 
     // some filter cuts:
@@ -217,4 +237,5 @@
     fMaxDist = 1.4;          // dist max cut (ever possible)
     fMinDist = 0.1;          // dist max cut (ever possible)
+    fHadrCut = 0.2;          // hadronness cut
 
     fNumBinsAlpha   = 36; 	 // number of bins for alpha histograms
@@ -282,11 +303,48 @@
     }
 
-    kSaveAlphaPlots=kTRUE;
-    kSaveSkyPlots=kTRUE;
-    kSaveNexPlot=kTRUE;
+    fSaveAlphaPlots=kTRUE;
+    fSaveSkyPlots=kTRUE;
+    fSaveNexPlot=kTRUE;
+    fUseRF=kFALSE;
     fAlphaHName = "alphaplot.root";
     fSkyHName   = "skyplot.root";
     fNexHName   = "Nexcess.gif";
-}
+
+    fHistAlpha = new TOrdCollection();
+    fHistAlpha->SetOwner();    
+
+}
+
+MSkyPlot::~MSkyPlot()
+{
+
+  if (fHistAlpha)
+    delete fHistAlpha;
+}
+
+// --------------------------------------------------------------------------
+//
+// Get i-th hist 
+//
+TH1D *MSkyPlot::GetAlphaPlot(Int_t i)
+{
+  if (GetSize() == 0)
+    return NULL;
+
+  return static_cast<TH1D*>(i==-1 ? fHistAlpha->At(GetSize()/2+1) : fHistAlpha->At(i));
+}
+
+// --------------------------------------------------------------------------
+//
+// Get i-th hist 
+//
+const TH1D *MSkyPlot::GetAlphaPlot(Int_t i) const 
+{
+  if (GetSize() == 0)
+    return NULL;
+
+  return static_cast<TH1D*>(i==-1 ? fHistAlpha->At(GetSize()/2+1) : fHistAlpha->At(i));
+}
+
 
 void MSkyPlot::ReadCuts(const TString parSCinit="OptSCParametersONOFFThetaRange0_1570mRad.root")
@@ -412,6 +470,4 @@
 }
 
-
-
 // --------------------------------------------------------------------------
 //
@@ -423,5 +479,7 @@
 Int_t MSkyPlot::PreProcess(MParList *plist)
 {
-*fLog << warn << "entering PreProcess in MSkyPlot::PreProcess" << endl; 
+
+  *fLog << warn << "entering PreProcess in MSkyPlot::PreProcess" << endl; 
+  
     fGeomCam = (MGeomCam*)plist->FindObject("MGeomCam");
     if (!fGeomCam)
@@ -438,5 +496,5 @@
 
 
-    fSrcPosCam = (MSrcPosCam*)plist->FindObject(AddSerialNumber("MSrcPosCam"));
+    fSrcPosCam = (MSrcPosCam*)plist->FindCreateObj(AddSerialNumber("MSrcPosCam"));
     if (!fSrcPosCam)
         *fLog << warn << "MSrcPosCam not found... no sky map." << endl;
@@ -477,5 +535,12 @@
         *fLog << err << "MNewImagePar not found... no sky map." << endl;
 
-
+    if(fUseRF)
+    {
+       fHadron = (MHadronness*)plist->FindObject(AddSerialNumber("MHadronness"));
+       if (!fHadron)
+           *fLog << err << "MHadronness not found although specified... no sky map." << endl;
+
+       *fLog << inf << "Hadronness cut set to : " << fHadrCut << endl;
+    }
 
     // FIXME: Because the pointing position could change we must check
@@ -495,7 +560,5 @@
     fHistSignif.SetBins(fNumStepsX, fMinXGrid-0.5*fBinStepGrid, fMaxXGrid+0.5*fBinStepGrid, 
 			fNumStepsY, fMinYGrid-0.5*fBinStepGrid, fMaxYGrid+0.5*fBinStepGrid);
-   // fHistSignif.SetDirectory(NULL);
     fHistSignif.SetFillStyle(4000);
-    fHistSignif.UseCurrentStyle();
 
    // fHistNexcess.SetName("SPNex2ndOrder");
@@ -503,7 +566,5 @@
     fHistNexcess.SetBins(fNumStepsX, fMinXGrid-0.5*fBinStepGrid, fMaxXGrid+0.5*fBinStepGrid, 
 			fNumStepsY, fMinYGrid-0.5*fBinStepGrid, fMaxYGrid+0.5*fBinStepGrid);
-   // fHistNexcess.SetDirectory(NULL);
     fHistNexcess.SetFillStyle(4000);
-    fHistNexcess.UseCurrentStyle();
 
    // fHistOn.SetName("SPOnCounted");
@@ -511,7 +572,5 @@
     fHistOn.SetBins(fNumStepsX, fMinXGrid-0.5*fBinStepGrid, fMaxXGrid+0.5*fBinStepGrid, 
 			fNumStepsY, fMinYGrid-0.5*fBinStepGrid, fMaxYGrid+0.5*fBinStepGrid);
-   // fHistOn.SetDirectory(NULL);
     fHistOn.SetFillStyle(4000);
-    fHistOn.UseCurrentStyle();
 
     //fHistSignifGaus.SetName("SignifDistrib");
@@ -519,17 +578,16 @@
     fHistSignifGaus.SetBins(100, -10., 10.); 
 
-    // prepare alpha plots
+    // create alpha histograms
+    for (Int_t i=0; i< fNumalphahist; i++)
+      {
     	// temp histogram for an alpha plot
-    TH1D *temp = new TH1D("alphaplot","alphaplot",fNumBinsAlpha,fAlphaLeftEdge,fAlphaRightEdge);
-    temp->SetDirectory(NULL);
-    fHistAlphaBinWidth = temp->GetBinWidth(1);
-    // create alpha histograms
-    for (Int_t i=0; i< fNumalphahist; i++) {
-    	fHistAlpha.push_back(*temp);
-	fHistAlpha[i].SetDirectory(NULL);
-    }
-//cout  << " fHistAlpha[10].GetBinContent(5) " << fHistAlpha[10].GetBinContent(5) << endl;
-
-    delete temp;
+	TH1D *temp = new TH1D("alphaplot","alphaplot",fNumBinsAlpha,fAlphaLeftEdge,fAlphaRightEdge);
+	temp->SetDirectory(NULL);
+	fHistAlpha->AddAt(temp,i);
+      }
+
+    fHistAlphaBinWidth = GetAlphaPlot()->GetBinWidth(1);
+//cout  << " (*fHistAlpha)[10].GetBinContent(5) " << (*fHistAlpha)[10].GetBinContent(5) << endl;
+
     return kTRUE;
 }
@@ -541,167 +599,186 @@
 Int_t MSkyPlot::Process()
 {
-// check whether MPointingPos comtains something:
+
+//      *fLog << err << "MPointingPos ENTERING the process" << endl;
+//      *fLog << err << "MPointingPos:: fUseRF " << (int)fUseRF << endl;
+  // check whether MPointingPos comtains something:
   if (TMath::Abs(fPointPos->GetRa())<1e-3 && TMath::Abs(fPointPos->GetDec())<1e-3) 
-  {
-     *fLog << warn << "MPointingPos is not filled ... event skipped" << endl;
-     return kTRUE;
-  }
-
-    // Get RA_0, DEC_0 for the camera center (Tracking MDrive?, can be set from outside)
+    {
+      *fLog << warn << "MPointingPos is not filled ... event skipped" << endl;
+      return kTRUE;
+    }
+  
+  // Get RA_0, DEC_0 for the camera center (Tracking MDrive?, can be set from outside)
   if (fSetCenter==kTRUE)
-  {
-     if (fRepDrive)
-     {
-	 fRa0  = fRepDrive->GetRa();  // [h]
-         fDec0 = fRepDrive->GetDec(); // [deg]
-         if (fRa0 < 0. || fRa0 > 24. || fDec0 < -90. || fDec0 > 90. || (fRa0==0 && fDec0==0)) return kTRUE;  // temp!, should be changed
-     }
-     else
-     {
-	fRa0 = fPointPos->GetRa(); 
-        fDec0 = fPointPos->GetDec();
-     } 
-     *fLog << inf << "Ra (center) = " << fRa0 << ", Dec = " << fDec0 << endl;
-     fSetCenter=kFALSE;
-  }
-
-
-// some filter cuts:
+    {
+      if (fRepDrive)
+	{
+	  fRa0  = fRepDrive->GetRa();  // [h]
+	  fDec0 = fRepDrive->GetDec(); // [deg]
+	  if (fRa0 < 0. || fRa0 > 24. || fDec0 < -90. || fDec0 > 90. || (fRa0==0 && fDec0==0)) return kTRUE;  // temp!, should be changed
+	}
+      else
+	{
+	  fRa0 = fPointPos->GetRa(); 
+	  fDec0 = fPointPos->GetDec();
+	} 
+      *fLog << inf << "Ra (center) = " << fRa0 << ", Dec = " << fDec0 << endl;
+      fSetCenter=kFALSE;
+    }
+
+  // some filter cuts:
   if ( fHillas->GetSize() > fSizeMin && fHillas->GetSize() < fSizeMax && fNewImagePar->GetLeakage1() < fLeakMax)
-  {
-
-    Double_t xsource, ysource, cosgam, singam, x_0, y_0, xsourcam, ysourcam;
-    Double_t dx, dy, beta, tanbeta, alphapar, distpar;
-    Double_t logsize, lgsize, lgsize2, dist2; 
-    const Double_t log3000 = log(3000.);
-    const Float_t fDistOffset = 0.9;
-
-    // Get camera rotation angle
-    Double_t rho = 0;
-    if (fTime && fObservatory && fPointPos)
-    {
-        rho = fPointPos->RotationAngle(*fObservatory, *fTime);
-//*fLog << inf << " rho = " << rho*180./TMath::Pi() << ", Zd = " << fPointPos->GetZd() << 
-//                ", Az = " << fPointPos->GetAz() << ", Ra = " << fPointPos->GetRa() << ", Dec = " << fPointPos->GetDec() << endl;
-
-    // => coord. system: xtilde, ytilde with center in RA_0, DEC_0 
-
-// Get Rot. Angle:
-    cosgam = TMath::Cos(rho);
-    singam = TMath::Sin(rho);
-// Get x_0, y_0 for RA_0,DEC_0 of the current event
-    }
-    else 
-    {
-//	rho = fPointPos->RotationAngle(*fObservatory);
-      Double_t theta, phi, sing, cosg;
-      theta = fPointPos->GetZd()*TMath::Pi()/180.;
-      phi = fPointPos->GetAz()*TMath::Pi()/180.;
-//   printf("theta: %5.3f, phi: %5.3f\n", theta*180./4.1415, phi*180./4.1415);
-      fObservatory->RotationAngle(theta, phi, sing, cosg);        
-      cosgam = cosg;
-      singam = sing;
-    }
-   // if (fPointPos)
-   //    rho = fPointPos->RotationAngle(*fObservatory);
-
-/*
-//TEMP
-//        theta = mcevt->GetTelescopeTheta();
-    Double_t theta, phi, sing, cosg;
-    theta = fPointPos->GetZd()*TMath::Pi()/180.;
-    phi = fPointPos->GetAz()*TMath::Pi()/180.;
-//   printf("theta: %5.3f, phi: %5.3f\n", theta*180./4.1415, phi*180./4.1415);
-    fObservatory->RotationAngle(theta, phi, sing, cosg);
+    {
       
-// conclusion: diffference in rho = 7 deg
-  *fLog << "new thetaTel, phiTel, cosal, sinal, rho = " << theta << ",  "
-        << phi << ",  " << cosg << ",  " << sing << ", " << TMath::ATan2(sing,cosg)*180./TMath::Pi() << endl;
-
-    Double_t a1 =  0.876627;
-    Double_t a3 = -0.481171;
-    Double_t thetaTel=theta, phiTel=phi;
-
-    Double_t denom =  1./ sqrt( sin(thetaTel)*sin(phiTel) * sin(thetaTel)*sin(phiTel) +
-                              ( a1*cos(thetaTel)+a3*sin(thetaTel)*cos(phiTel) ) *
-                              ( a1*cos(thetaTel)+a3*sin(thetaTel)*cos(phiTel) )   );
-    Double_t cosal = - (a3 * sin(thetaTel) + a1 * cos(thetaTel) * cos(phiTel)) * denom;
-    Double_t sinal =    a1 * sin(phiTel) * denom;
-
-  *fLog << "old thetaTel, phiTel, cosal, sinal, rho = " << thetaTel << ",  "
-        << phiTel << ",  " << cosal << ",  " << sinal << ", " << TMath::ATan2(sinal,cosal)*180./TMath::Pi() << endl;
-
-
-
-// END TEMP
-*/
-
-// make the center of the plot different from the center of the camera
-/*
-    x_0 = fPntPosCam->GetX()*fMm2Deg; 
-    y_0 = fPntPosCam->GetY()*fMm2Deg; 
-*/      
-    x_0 = 0.;  
-    y_0 = 0.;
-
-    Int_t index = 0;  // index for alpha histograms
-    // loop over xtilde
-    for (Int_t gridy = 0; gridy < fNumStepsY; gridy++)   
-    {
-	ysource = fMinYGrid + gridy * fBinStepGrid;
-        // loop over ytilde
-	for (Int_t gridx = 0; gridx < fNumStepsX; gridx++)
+      Double_t xsource, ysource, cosgam, singam, x_0, y_0, xsourcam, ysourcam;
+      Double_t dx, dy, beta, tanbeta, alphapar, distpar;
+      Double_t logsize, lgsize, lgsize2, dist2, hadr; 
+      const Double_t log3000 = log(3000.);
+      const Float_t fDistOffset = 0.9;
+      
+      //Get Hadronness if available:
+      if(fUseRF) 
 	{
-
-	   xsource = fMinXGrid + gridx * fBinStepGrid;
-
- /*     derotation    : rotate  into camera coordinates */
- /*     formel: (x_cam)    (cos(gam)  -sin(gam))   (xtilde)   (x_0)
-                (     )  = (                   ) * (      ) + (   ) 
-	        (y_cam)    (sin(gam)   sin(gam))   (ytilde)   (y_0)
-*/
-           xsourcam = cosgam * xsource - singam * ysource + x_0;
-           ysourcam = singam * xsource + cosgam * ysource + y_0;
-
-
- /*    end derotatiom    */
-//           xsourcam = xsource;
-//           ysourcam = ysource;
-
-       /* calculate ALPHA und DIST according to the source position */
-           dx = fHillas->GetMeanX()*fMm2Deg - xsourcam;
-           dy = fHillas->GetMeanY()*fMm2Deg - ysourcam;
-           tanbeta = dy / dx ;
-           beta = TMath::ATan(tanbeta);
-           alphapar = (fHillas->GetDelta() - beta) * 180./ TMath::Pi();
-           distpar = sqrt( dy*dy + dx*dx );
-           if(alphapar >  90.) alphapar -= 180.;
-           if(alphapar < -90.) alphapar += 180.;
-           alphapar = TMath::Abs(alphapar);
-
-// apply cuts
-           logsize = log(fHillas->GetSize());
-           lgsize = logsize-log3000;
-           lgsize2 = lgsize*lgsize;
-           dist2 = distpar*distpar - fDistOffset*fDistOffset;
-
-           if ( (fHillas->GetLength()*fMm2Deg) < CalcLimit(fLengthUp, lgsize, lgsize2, dist2) &&
-                (fHillas->GetLength()*fMm2Deg) > CalcLimit(fLengthLo, lgsize, lgsize2, dist2))
-           if ( (fHillas->GetWidth()*fMm2Deg) < CalcLimit(fWidthUp, lgsize, lgsize2, dist2) &&
-                (fHillas->GetWidth()*fMm2Deg) > CalcLimit(fWidthLo, lgsize, lgsize2, dist2))
-           if ( distpar < CalcLimit(fDistUp, lgsize, lgsize2, dist2) &&
-                      distpar > CalcLimit(fDistLo, lgsize, lgsize2, dist2) &&
-                      distpar < fMaxDist && distpar > fMinDist)
-	   {
-// gamma candidates!
-//*fLog <<  "Length : " << fHillas->GetLength()*fMm2Deg << ", Width : " << fHillas->GetWidth()*fMm2Deg << endl;
-//*fLog <<  "distpar: " << distpar << ", Size : " << fHillas->GetSize() << endl;
-		fHistAlpha[index].Fill(alphapar);
-	   }
-	   index++;
+		hadr=fHadron->GetHadronness();
+//		cout << " will use RF " << hadr << endl;
 	}
-     }
-  }
-
+      // Get camera rotation angle
+      Double_t rho = 0;
+      if (fTime && fObservatory && fPointPos)
+	{
+	  rho = fPointPos->RotationAngle(*fObservatory, *fTime);
+	  //*fLog << inf << " rho = " << rho*180./TMath::Pi() << ", Zd = " << fPointPos->GetZd() << 
+	  //                ", Az = " << fPointPos->GetAz() << ", Ra = " << fPointPos->GetRa() << ", Dec = " << fPointPos->GetDec() << endl;
+	  
+	  // => coord. system: xtilde, ytilde with center in RA_0, DEC_0 
+	  
+	  // Get Rot. Angle:
+	  cosgam = TMath::Cos(rho);
+	  singam = TMath::Sin(rho);
+	  // Get x_0, y_0 for RA_0,DEC_0 of the current event
+	}
+      else 
+	{
+	  //	rho = fPointPos->RotationAngle(*fObservatory);
+	  Double_t theta, phi, sing, cosg;
+	  theta = fPointPos->GetZd()*TMath::Pi()/180.;
+	  phi = fPointPos->GetAz()*TMath::Pi()/180.;
+	  printf("theta: %5.3f, phi: %5.3f\n", theta*180./4.1415, phi*180./4.1415);
+	  fObservatory->RotationAngle(theta, phi, sing, cosg);        
+	  cosgam = cosg;
+	  singam = sing;
+	}
+      // if (fPointPos)
+      //    rho = fPointPos->RotationAngle(*fObservatory);
+      
+      /*
+	//TEMP
+	//        theta = mcevt->GetTelescopeTheta();
+	Double_t theta, phi, sing, cosg;
+	theta = fPointPos->GetZd()*TMath::Pi()/180.;
+	phi = fPointPos->GetAz()*TMath::Pi()/180.;
+	//   printf("theta: %5.3f, phi: %5.3f\n", theta*180./4.1415, phi*180./4.1415);
+	fObservatory->RotationAngle(theta, phi, sing, cosg);
+	
+	// conclusion: diffference in rho = 7 deg
+	//  *fLog << "new thetaTel, phiTel, cosal, sinal, rho = " << theta << ",  "
+	//        << phi << ",  " << cosg << ",  " << sing << ", " << TMath::ATan2(sing,cosg)*180./TMath::Pi() << endl;
+	
+	Double_t a1 =  0.876627;
+	Double_t a3 = -0.481171;
+	Double_t thetaTel=theta, phiTel=phi;
+	
+	Double_t denom =  1./ sqrt( sin(thetaTel)*sin(phiTel) * sin(thetaTel)*sin(phiTel) +
+	( a1*cos(thetaTel)+a3*sin(thetaTel)*cos(phiTel) ) *
+	( a1*cos(thetaTel)+a3*sin(thetaTel)*cos(phiTel) )   );
+	Double_t cosal = - (a3 * sin(thetaTel) + a1 * cos(thetaTel) * cos(phiTel)) * denom;
+	Double_t sinal =    a1 * sin(phiTel) * denom;
+	
+	//  *fLog << "old thetaTel, phiTel, cosal, sinal, rho = " << thetaTel << ",  "
+	//        << phiTel << ",  " << cosal << ",  " << sinal << ", " << TMath::ATan2(sinal,cosal)*180./TMath::Pi() << endl;
+	
+	// END TEMP
+      */
+	    
+      // make the center of the plot different from the center of the camera
+      /*
+	x_0 = fPntPosCam->GetX()*fMm2Deg; 
+	y_0 = fPntPosCam->GetY()*fMm2Deg; 
+      */      
+      x_0 = 0.;  
+      y_0 = 0.;
+      
+      Int_t index = 0;  // index for alpha histograms
+      // loop over xtilde
+      for (Int_t gridy = 0; gridy < fNumStepsY; gridy++)   
+	{
+	  ysource = fMinYGrid + gridy * fBinStepGrid;
+	  // loop over ytilde
+	  for (Int_t gridx = 0; gridx < fNumStepsX; gridx++)
+	    {
+	      
+	      xsource = fMinXGrid + gridx * fBinStepGrid;
+	      
+	      /*     derotation    : rotate  into camera coordinates */
+	      /*     formel: (x_cam)      (cos(gam)  -sin(gam))   (xtilde)   (x_0)
+		     (     )  = - (                   ) * (      ) + (   ) 
+	        (y_cam)      (sin(gam)   cos(gam))   (ytilde)   (y_0)
+	      */
+	      //xsourcam = - (cosgam * xsource - singam * ysource) + x_0;
+	      //ysourcam = - (singam * xsource + cosgam * ysource) + y_0;
+	      
+	      
+	      /*    end derotatiom    */
+	                 xsourcam = xsource;
+	                 ysourcam = ysource;
+
+	      /* calculate ALPHA und DIST according to the source position */
+	      dx = fHillas->GetMeanX()*fMm2Deg - xsourcam;
+	      dy = fHillas->GetMeanY()*fMm2Deg - ysourcam;
+	      tanbeta = dy / dx ;
+	      beta = TMath::ATan(tanbeta);
+	      alphapar = (fHillas->GetDelta() - beta) * 180./ TMath::Pi();
+	      distpar = sqrt( dy*dy + dx*dx );
+	      if(alphapar >  90.) alphapar -= 180.;
+	      if(alphapar < -90.) alphapar += 180.;
+	      alphapar = TMath::Abs(alphapar);
+	      
+	      if(fUseRF)
+	      {
+		
+//		cout << " will use RF " << hadr << endl;
+		if(hadr<fHadrCut  &&  distpar < fMaxDist && distpar > fMinDist)
+		{
+		    TH1D *hist = GetAlphaPlot(index);
+                    hist->Fill(alphapar);
+		}
+	      }
+	      else
+	      {
+	      	// apply cuts
+	      	logsize = log(fHillas->GetSize());
+	      	lgsize = logsize-log3000;
+	      	lgsize2 = lgsize*lgsize;
+	      	dist2 = distpar*distpar - fDistOffset*fDistOffset;
+	      
+	      	if ( (fHillas->GetLength()*fMm2Deg) < CalcLimit(fLengthUp, lgsize, lgsize2, dist2) &&
+		   (fHillas->GetLength()*fMm2Deg) > CalcLimit(fLengthLo, lgsize, lgsize2, dist2))
+			if ( (fHillas->GetWidth()*fMm2Deg) < CalcLimit(fWidthUp, lgsize, lgsize2, dist2) &&
+		     	(fHillas->GetWidth()*fMm2Deg) > CalcLimit(fWidthLo, lgsize, lgsize2, dist2))
+			  if ( distpar < CalcLimit(fDistUp, lgsize, lgsize2, dist2) &&
+		      	  distpar > CalcLimit(fDistLo, lgsize, lgsize2, dist2) &&
+		      	  distpar < fMaxDist && distpar > fMinDist)
+		    	{
+		      // gamma candidates!
+		      //*fLog <<  "Length : " << fHillas->GetLength()*fMm2Deg << ", Width : " << fHillas->GetWidth()*fMm2Deg << endl;
+		      //*fLog <<  "distpar: " << distpar << ", Size : " << fHillas->GetSize() << endl;
+		      	TH1D *hist = GetAlphaPlot(index);
+		      	hist->Fill(alphapar);
+		    	}
+		}
+	      	index++;
+	    }
+	}
+    }
     return kTRUE;
 }
@@ -724,8 +801,8 @@
     numdegfreed = (fAlphaBgUp - fAlphaBgLow) / fHistAlphaBinWidth - 2.; 
 
-    vector <TH1D>::iterator alpha_iterator;
     int index2 = 0;  // index of the TH2F histograms
 
-    alpha_iterator = fHistAlpha.begin();
+    TOrdCollectionIter Next(fHistAlpha);
+    TH1D *alpha_iterator = NULL;
 
     fHistNexcess.Reset();
@@ -734,8 +811,8 @@
     fHistSignifGaus.Reset();
     
-    while( alpha_iterator != fHistAlpha.end() ) {
+    while ( (alpha_iterator = (TH1D*)Next())) {
 	 // find the bin in the 2dim histogram
-         nrow    = index2/fHistOn.GetNbinsX() + 1;                    // row of the histogram (y)
-         ncolumn = index2%fHistOn.GetNbinsX()+1;   // column of the histogram (x)
+         nrow    = index2/fHistOn.GetNbinsX() + 1;              // row of the histogram (y)
+         ncolumn = index2%fHistOn.GetNbinsX()+1;   		// column of the histogram (x)
          //ncolumn = TMath::Nint(fmod(index2,fHistOn.GetNbinsX()))+1;   // column of the histogram (x)
 
@@ -747,58 +824,60 @@
 
          // fit parabola to a background region
-         (*alpha_iterator).Fit(fitbgpar,"EQRN");  // NWR OK?????????????????????????
+         alpha_iterator->Fit(fitbgpar,"EQRN");  // NWR OK?????????????????????????
          // get Chi2
          chisquarefit = fitbgpar->GetChisquare();
          if (chisquarefit/numdegfreed<2. && chisquarefit/numdegfreed>0.01);
          else  *fLog << warn << "Fit is bad, chi2/ndf = " << chisquarefit/numdegfreed << endl;
-
+	 
          // estimate Noff from the fit:
        	 Noff_fit = (1./3. * (fitbgpar->GetParameter(0)) * TMath::Power(fAlphaONMax,3.) +
-       	            (fitbgpar->GetParameter(1)) * fAlphaONMax) /  fHistAlphaBinWidth;
-
+		     (fitbgpar->GetParameter(1)) * fAlphaONMax) /  fHistAlphaBinWidth;
+	 
 	 Nex = Non - Noff_fit;
 
          fHistNexcess.SetBinContent(ncolumn, nrow, Nex);	// fill in the fHistNexcess
-
+	 
          if (Noff_fit<0.) Sign = 0.;
-//         else Sign = LiMa17(Non,Noff_fit,1.);
+	 //         else Sign = LiMa17(Non,Noff_fit,1.);
          else Sign = MMath::SignificanceLiMaSigned(Non, Noff_fit, 1.);
-
+	 
          fHistSignif.SetBinContent(ncolumn, nrow, Sign);	// fill in the fHistSignif
          fHistSignifGaus.Fill(Sign);
-
-         alpha_iterator++;
+	 
          index2++;
     }
-
-// fit with gaus 
+    
+    // fit with gaus 
     fHistSignifGaus.Fit("gaus","N");
-
-
-//temp
-/*
-   Double_t decl, hang;
-   MStarCamTrans mstarc(*fGeomCam, *fObservatory);
-   mstarc.LocToCel(fRepDrive->GetNominalZd(),fRepDrive->GetNominalAz(),decl, hang);
-
-*fLog << warn << "MDrive, RA, DEC = " << fRepDrive->GetRa() << ", " << fRepDrive->GetDec() << endl;
-*fLog << warn << "MStarCamTrans, H angle , DEC = " << hang << ", " << decl << endl;
-*/
-//endtemp
-
-
-// save alpha plots:
-//  TString stri1 = "alphaplots.root";
-  if(kSaveAlphaPlots==kTRUE) SaveAlphaPlots(fAlphaHName);
-
-// save sky plots:
-//  TString stri2 = "skyplots.root";
-  if(kSaveSkyPlots==kTRUE) SaveSkyPlots(fSkyHName);
-
-// save nex plot:
-  if(kSaveNexPlot==kTRUE) SaveNexPlot(fNexHName);
-
-  fHistAlpha.clear();
-  return kTRUE;
+    
+    
+    //temp
+    /*
+      Double_t decl, hang;
+      MStarCamTrans mstarc(*fGeomCam, *fObservatory);
+      mstarc.LocToCel(fRepDrive->GetNominalZd(),fRepDrive->GetNominalAz(),decl, hang);
+      
+      *fLog << warn << "MDrive, RA, DEC = " << fRepDrive->GetRa() << ", " << fRepDrive->GetDec() << endl;
+      *fLog << warn << "MStarCamTrans, H angle , DEC = " << hang << ", " << decl << endl;
+      */
+    //endtemp
+    
+    
+    // save alpha plots:
+    //  TString stri1 = "alphaplots.root";
+    if(fSaveAlphaPlots==kTRUE) SaveAlphaPlots(fAlphaHName);
+    
+    // save sky plots:
+    //  TString stri2 = "skyplots.root";
+    if(fSaveSkyPlots==kTRUE) SaveSkyPlots(fSkyHName);
+    
+    // save nex plot:
+    if(fSaveNexPlot==kTRUE) SaveNexPlot(fNexHName);
+    
+    fHistAlpha->Clear();
+
+    delete fitbgpar;
+    
+    return kTRUE;
 }
 
@@ -843,4 +922,6 @@
     gStyle->SetPadLeftMargin(0.13);
 
+    Char_t timet[100];
+
     TH2D tmp = fHistNexcess;
     tmp.SetMaximum(470); 
@@ -861,4 +942,10 @@
     gPad->SetTheta(20);
     fHistNexcess.Draw("lego2");
+    TLatex tu(0.5,0.8,"elapsed time:");
+    tu.Draw();
+    sprintf(timet,"%.1f min",fElaTime);
+    TLatex tut(0.5,0.7,timet);
+    tut.Draw();
+
     can.Print(nexplotname);  
 //    can.Print("test.root");
@@ -867,137 +954,137 @@
 void MSkyPlot::SaveSkyPlots(TString skyplotfilename)
 {
-   TFile rootfile(skyplotfilename, "RECREATE",
-                   "sky plots in some variations");
-     fHistSignif.Write();
-     fHistNexcess.Write();
-     fHistOn.Write();
-     fHistSignif.Write();
-
-// from Wolfgang: 
-   //--------------------------------------------------------------------
-    // Dec-Hour lines
-    Double_t rho, sinrho, cosrho;
-    Double_t theta, phi, sing, cosg;
-// do I need it?
-    if (fTime && fObservatory && fPointPos)
-    {
-        rho = fPointPos->RotationAngle(*fObservatory, *fTime);
-        sinrho=TMath::Sin(rho);
-        cosrho=TMath::Cos(rho);
-    }
-
-       theta = fPointPos->GetZd()*TMath::Pi()/180.;
-       phi = fPointPos->GetAz()*TMath::Pi()/180.;
-//   printf("theta: %5.3f, phi: %5.3f\n", theta*180./4.1415, phi*180./4.1415);
-       fObservatory->RotationAngle(theta, phi, sing, cosg);
-
-*fLog << "1: sinrho, cosrho = " << sinrho << ", " << cosrho << endl;
-*fLog << "2: sinrho, cosrho = " << sing << ", " << cosg << endl;
-       sinrho=sing;
-       cosrho=cosg;
-
-    Double_t fDistCam = fGeomCam->GetCameraDist() * 1000.0;     //  [mm]
-    Double_t gridbinning = fGridBinning;
-    Double_t gridfinebin = fGridFineBin;
-    Int_t    numgridlines = (Int_t)(4.0/gridbinning);
-    Double_t aberr = 1.07;
-    Double_t mmtodeg = 180.0 / TMath::Pi() / fDistCam ;
-
-    Double_t declin, hangle;  // [deg], [h]
-    MStarCamTrans mstarc(*fGeomCam, *fObservatory);
-    if (fRepDrive) mstarc.LocToCel(fRepDrive->GetNominalZd(),fRepDrive->GetNominalAz(),declin, hangle);
-    else mstarc.LocToCel(theta*180./TMath::Pi(),phi*180./TMath::Pi(),declin, hangle); // NOT GOOD!
-
-//    std::vector <TGraph> graphvec;
-    TLatex *pix;
-
-    Char_t tit[100];
-    Double_t xtxt;
-    Double_t ytxt;
-
-    Double_t xlo = -534.0 * mmtodeg;
-    Double_t xup =  534.0 * mmtodeg;
-
-    Double_t ylo = -534.0 * mmtodeg;
-    Double_t yup =  534.0 * mmtodeg;
-
-    Double_t xx, yy;
-
-
-    // direction for the center of the camera
-    Double_t dec0 = declin;
-    Double_t h0   = hangle*360./24.; //deg
-    Double_t RaHOffset = fRepDrive->GetRa() - h0;
+
+  TFile rootfile(skyplotfilename, "RECREATE",
+		 "sky plots in some variations");
+  fHistSignif.Write();
+  fHistNexcess.Write();
+  fHistOn.Write();
+  fHistSignif.Write();
+     
+  // from Wolfgang: 
+  //--------------------------------------------------------------------
+  // Dec-Hour lines
+  Double_t rho, sinrho, cosrho;
+  Double_t theta, phi, sing, cosg;
+  // do I need it?
+  if (fTime && fObservatory && fPointPos)
+    {
+      rho = fPointPos->RotationAngle(*fObservatory, *fTime);
+      sinrho=TMath::Sin(rho);
+      cosrho=TMath::Cos(rho);
+    }
+  
+  theta = fPointPos->GetZd()*TMath::Pi()/180.;
+  phi = fPointPos->GetAz()*TMath::Pi()/180.;
+  //   printf("theta: %5.3f, phi: %5.3f\n", theta*180./4.1415, phi*180./4.1415);
+  fObservatory->RotationAngle(theta, phi, sing, cosg);
+  
+  *fLog << "1: sinrho, cosrho = " << sinrho << ", " << cosrho << endl;
+  *fLog << "2: sinrho, cosrho = " << sing << ", " << cosg << endl;
+  sinrho=sing;
+  cosrho=cosg;
+
+  Double_t fDistCam = fGeomCam->GetCameraDist() * 1000.0;     //  [mm]
+  Double_t gridbinning = fGridBinning;
+  Double_t gridfinebin = fGridFineBin;
+  Int_t    numgridlines = (Int_t)(4.0/gridbinning);
+  Double_t aberr = 1.07;
+  Double_t mmtodeg = 180.0 / TMath::Pi() / fDistCam ;
+  
+  Double_t declin, hangle;  // [deg], [h]
+  MStarCamTrans mstarc(*fGeomCam, *fObservatory);
+  if (fRepDrive) mstarc.LocToCel(fRepDrive->GetNominalZd(),fRepDrive->GetNominalAz(),declin, hangle);
+  else mstarc.LocToCel(theta*180./TMath::Pi(),phi*180./TMath::Pi(),declin, hangle); // NOT GOOD!
+  
+  TLatex *pix;
+  
+  Char_t tit[100];
+  Double_t xtxt;
+  Double_t ytxt;
+  
+  Double_t xlo = -534.0 * mmtodeg;
+  Double_t xup =  534.0 * mmtodeg;
+  
+  Double_t ylo = -534.0 * mmtodeg;
+  Double_t yup =  534.0 * mmtodeg;
+  
+  Double_t xx, yy;
+  
+  
+  // direction for the center of the camera
+  Double_t dec0 = declin;
+  Double_t h0   = hangle*360./24.; //deg
+  //    Double_t RaHOffset = fRepDrive->GetRa() - h0;
     //trans.LocToCel(theta0, phi0, dec0, h0);
-
-    gStyle->SetOptFit(0);
-    gStyle->SetOptStat(0);
-    gStyle->SetPalette(1);
-    TCanvas *c1 = new TCanvas("SPlotsRaDecLines","SPlotsRaDecLines", 400, 0, 700, 600);
-    c1->Divide(2,2);
-    c1->cd(1);
-    fHistSignif.Draw("colz");
-    c1->cd(2);
-    fHistNexcess.Draw("colz");
-    c1->cd(3);
-    fHistOn.Draw("colz");
-    c1->cd(4);
-    gPad->SetLogy();
-    fHistSignifGaus.Draw(); 
-
-    //-----   lines for fixed dec   ------------------------------------
-
-    const Int_t Ndec = numgridlines;
-    Double_t dec[Ndec];
-    Double_t ddec = gridbinning;
-
-    Int_t nh = (Int_t)(4.0/gridfinebin);
-    const Int_t Nh   = nh+1;
-    Double_t h[Nh];
-    Double_t dh = gridfinebin/cos(dec0/180.0*3.1415926);
-    if ( dh > 360.0/((Double_t)(Nh-1)) )
-      dh = 360.0/((Double_t)(Nh-1));
-
-// start to copy
-    for (Int_t j=0; j<Ndec; j++)
+  
+  gStyle->SetOptFit(0);
+  gStyle->SetOptStat(0);
+  gStyle->SetPalette(1);
+  TCanvas *c1 = new TCanvas("SPlotsRaDecLines","SPlotsRaDecLines", 400, 0, 700, 600);
+  c1->Divide(2,2);
+  c1->cd(1);
+  fHistSignif.Draw("colz");
+  c1->cd(2);
+  fHistNexcess.Draw("colz");
+  c1->cd(3);
+  fHistOn.Draw("colz");
+  c1->cd(4);
+  gPad->SetLogy();
+  fHistSignifGaus.Draw(); 
+  
+  //-----   lines for fixed dec   ------------------------------------
+  
+  const Int_t Ndec = numgridlines;
+  Double_t dec[Ndec];
+  Double_t ddec = gridbinning;
+  
+  Int_t nh = (Int_t)(4.0/gridfinebin);
+  const Int_t Nh   = nh+1;
+  Double_t h[Nh];
+  Double_t dh = gridfinebin/cos(dec0/180.0*3.1415926);
+  if ( dh > 360.0/((Double_t)(Nh-1)) )
+    dh = 360.0/((Double_t)(Nh-1));
+  
+  // start to copy
+  for (Int_t j=0; j<Ndec; j++)
     {
       dec[j] = dec0 + ((Double_t)j)*ddec
-                        - ((Double_t)(Ndec/2)-1.0)*ddec;
-    }
-
-    for (Int_t k=0; k<Nh; k++)
+	- ((Double_t)(Ndec/2)-1.0)*ddec;
+    }
+  
+  for (Int_t k=0; k<Nh; k++)
     {
       h[k] = h0 + ((Double_t)k)*dh - ((Double_t)(Nh/2)-1.0)*dh;
     }
-
-    Double_t xh[Nh];
-    Double_t yh[Nh];
-
-    for (Int_t j=0; j<Ndec; j++)
+  
+  Double_t xh[Nh];
+  Double_t yh[Nh];
+  
+  for (Int_t j=0; j<Ndec; j++)
     {
       if (fabs(dec[j]) > 90.0) continue;
-
+      
       for (Int_t k=0; k<Nh; k++)
-      {
-        Double_t hh0 = h0   *24.0/360.0;
-        Double_t hx = h[k]*24.0/360.0;
-        mstarc.Cel0CelToCam(dec0, hh0, dec[j], hx,
-                           xx, yy);
-        xx = xx * mmtodeg * aberr;
-        yy = yy * mmtodeg * aberr;
-//        xh[k] = xx * mmtodeg * aberr;
-//        yh[k] = yy * mmtodeg * aberr;
-        xh[k] = cosg * xx + sing * yy;
-        yh[k] =-sing * xx + cosg * yy;
-//        xh[k] = cosrho * xx + sinrho * yy;
-//        yh[k] =-sinrho * xx + cosrho * yy;
-
-
-        //gLog << "dec0, h0 = " << dec0 << ",  " << h0
-        //     << " : dec, h, xh, yh = " << dec[j] << ",  "
-        //     << h[k] << ";   " << xh[k] << ",  " << yh[k] << endl;
-      }
-
-//      c1->cd(2);
+	{
+	  Double_t hh0 = h0   *24.0/360.0;
+	  Double_t hx = h[k]*24.0/360.0;
+	  mstarc.Cel0CelToCam(dec0, hh0, dec[j], hx,
+			      xx, yy);
+	  xx = xx * mmtodeg * aberr;
+	  yy = yy * mmtodeg * aberr;
+	  //        xh[k] = xx * mmtodeg * aberr;
+	  //        yh[k] = yy * mmtodeg * aberr;
+	  xh[k] = cosg * xx + sing * yy;
+	  yh[k] =-sing * xx + cosg * yy;
+	  //        xh[k] = cosrho * xx + sinrho * yy;
+	  //        yh[k] =-sinrho * xx + cosrho * yy;
+	  
+	  
+	  //gLog << "dec0, h0 = " << dec0 << ",  " << h0
+	  //     << " : dec, h, xh, yh = " << dec[j] << ",  "
+	  //     << h[k] << ";   " << xh[k] << ",  " << yh[k] << endl;
+	}
+      
+      //      c1->cd(2);
       TGraph * graph1 = new TGraph(Nh,xh,yh);
       //graph1->SetLineColor(j+1);
@@ -1017,7 +1104,7 @@
       graph1->Draw("L");
       //delete graph1;
-//      graphvec.push_back(*graph1);
-//      graphvec[j].Draw("L");
-
+      //      graphvec.push_back(*graph1);
+      //      graphvec[j].Draw("L");
+      
       sprintf(tit,"Dec = %6.2f", dec[j]);
       xtxt = xlo + (xup-xlo)*0.1;
@@ -1028,63 +1115,62 @@
       //pix->Draw("");
       //delete pix;
-
-    }
-//stop copy
-    //-----   lines for fixed hour angle   ------------------------------------
-
-    Int_t ndec1 = (Int_t)(4.0/gridfinebin);
-    const Int_t Ndec1 = ndec1;
-    Double_t dec1[Ndec1];
-    Double_t ddec1 = gridfinebin;
-
-    const Int_t Nh1   = numgridlines;
-    Double_t h1[Nh1];
-    Double_t dh1 = gridbinning/cos(dec0/180.0*3.1415926);
-    if ( dh1 > 360.0/((Double_t)Nh1) )
-      dh1 = 360.0/((Double_t)Nh1);
-
-// start copy
-    for (Int_t j=0; j<Ndec1; j++)
+      
+    }
+  //stop copy
+  //-----   lines for fixed hour angle   ------------------------------------
+  
+  Int_t ndec1 = (Int_t)(4.0/gridfinebin);
+  const Int_t Ndec1 = ndec1;
+  Double_t dec1[Ndec1];
+  Double_t ddec1 = gridfinebin;
+  
+  const Int_t Nh1   = numgridlines;
+  Double_t h1[Nh1];
+  Double_t dh1 = gridbinning/cos(dec0/180.0*3.1415926);
+  if ( dh1 > 360.0/((Double_t)Nh1) )
+    dh1 = 360.0/((Double_t)Nh1);
+  
+  // start copy
+  for (Int_t j=0; j<Ndec1; j++)
     {
       dec1[j] = dec0 + ((Double_t)j)*ddec1
-                        - ((Double_t)(Ndec1/2)-1.0)*ddec1;
-    }
-
-    for (Int_t k=0; k<Nh1; k++)
+	- ((Double_t)(Ndec1/2)-1.0)*ddec1;
+    }
+  
+  for (Int_t k=0; k<Nh1; k++)
     {
       h1[k] = h0 + ((Double_t)k)*dh1 - ((Double_t)(Nh1/2)-1.0)*dh1;
     }
-
-
-    Double_t xd[Ndec1];
-    Double_t yd[Ndec1];
-
-    for (Int_t k=0; k<Nh1; k++)
+  
+  Double_t xd[Ndec1];
+  Double_t yd[Ndec1];
+  
+  for (Int_t k=0; k<Nh1; k++)
     {
       Int_t count = 0;
       for (Int_t j=0; j<Ndec1; j++)
-      {
-        if (fabs(dec1[j]) > 90.0) continue;
-
-        Double_t hh0 = h0   *24.0/360.0;
-        Double_t hhx = h1[k]*24.0/360.0;
-        mstarc.Cel0CelToCam(dec0, hh0, dec1[j], hhx,
+	{
+	  if (fabs(dec1[j]) > 90.0) continue;
+	  
+	  Double_t hh0 = h0   *24.0/360.0;
+	  Double_t hhx = h1[k]*24.0/360.0;
+	  mstarc.Cel0CelToCam(dec0, hh0, dec1[j], hhx,
                            xx, yy);
-//        xd[count] = xx * mmtodeg * aberr;
-//        yd[count] = yy * mmtodeg * aberr;
-
-        xx = xx * mmtodeg * aberr;
-        yy = yy * mmtodeg * aberr;
-        xd[count] = cosg * xx + sing * yy;
-        yd[count] =-sing * xx + cosg * yy;
-
-        //gLog << "dec0, h0 = " << dec0 << ",  " << h0
-        //     << " : dec1, h1, xd, yd = " << dec1[j] << ",  "
-        //     << h1[k] << ";   " << xd[count] << ",  " << yd[count] << endl;
-
-        count++;
-      }
-
-//      c1->cd(2);
+	  //        xd[count] = xx * mmtodeg * aberr;
+	  //        yd[count] = yy * mmtodeg * aberr;
+	  
+	  xx = xx * mmtodeg * aberr;
+	  yy = yy * mmtodeg * aberr;
+	  xd[count] = cosg * xx + sing * yy;
+	  yd[count] =-sing * xx + cosg * yy;
+
+	  //gLog << "dec0, h0 = " << dec0 << ",  " << h0
+	  //     << " : dec1, h1, xd, yd = " << dec1[j] << ",  "
+	  //     << h1[k] << ";   " << xd[count] << ",  " << yd[count] << endl;
+	  
+	  count++;
+	}
+      
+      //      c1->cd(2);
       TGraph * graph1 = new TGraph(count,xd,yd);
       //graph1->SetLineColor(k+1);
@@ -1111,42 +1197,41 @@
       //pix->Draw("");
       //delete pix;
-
-    }
-
-//    c1->cd(2);
-    sprintf(tit,"Dec0 = %6.2f [deg]   Ra0 = %6.2f [h]", dec0, fRa0);
-    xtxt = xlo + (xup-xlo)*0.05 + 0.80;
-    ytxt = ylo + (yup-ylo)*0.75;
-    pix = new TLatex(xtxt, ytxt, tit);
-    pix->SetTextColor(1);
-    pix->SetTextSize(.05);
-    c1->cd(1);
-    pix->Draw("");
-    c1->cd(2);
-    pix->Draw("");
-    c1->cd(3);
-    pix->Draw("");
-    //delete pix;
-
-    c1->Write();
-    // we suppose that the {skyplotfilename} ends with .root
-    Int_t sizeofout = skyplotfilename.Sizeof();
-    TString outps = skyplotfilename.Remove(sizeofout-5,5) + "ps";
-    c1->Print(outps);  // temporary!!!!!
-
-    TCanvas *c2 = new TCanvas("SkyPlotsWithRaDecLines","SkyPlotsWithRaDecLines", 0, 0, 300, 600);
-    c2->Divide(1,2);
-    c2->SetBorderMode(0);
-    c2->cd(1);
-    fHistSignif.Draw("colz");
-    c2->cd(2);
-    fHistNexcess.Draw("colz");
-    c2->Write();
-
-   rootfile.Close();
-   delete c1;
-   delete c2;
-   delete pix;
-
+      
+    }
+  
+  //    c1->cd(2);
+  sprintf(tit,"Dec0 = %6.2f [deg]   Ra0 = %6.2f [h]", dec0, fRa0);
+  xtxt = xlo + (xup-xlo)*0.05 + 0.80;
+  ytxt = ylo + (yup-ylo)*0.75;
+  pix = new TLatex(xtxt, ytxt, tit);
+  pix->SetTextColor(1);
+  pix->SetTextSize(.05);
+  c1->cd(1);
+  pix->Draw("");
+  c1->cd(2);
+  pix->Draw("");
+  c1->cd(3);
+  pix->Draw("");
+  //delete pix;
+  
+  c1->Write();
+  // we suppose that the {skyplotfilename} ends with .root
+  Int_t sizeofout = skyplotfilename.Sizeof();
+  TString outps = skyplotfilename.Remove(sizeofout-5,5) + "ps";
+  c1->Print(outps);  // temporary!!!!!
+
+  TCanvas *c2 = new TCanvas("SkyPlotsWithRaDecLines","SkyPlotsWithRaDecLines", 0, 0, 300, 600);
+  c2->Divide(1,2);
+  c2->SetBorderMode(0);
+  c2->cd(1);
+  fHistSignif.Draw("colz");
+  c2->cd(2);
+  fHistNexcess.Draw("colz");
+  c2->Write();
+  
+  rootfile.Close();
+  delete c1;
+  delete c2;
+  delete pix;
 
 }
@@ -1157,5 +1242,4 @@
                    "all the alpha plots");
 
-    vector <TH1D>::iterator alpha_iterator;
     int index = 0;  // index of the TH2F histograms
     Char_t strtitle[100];
@@ -1164,8 +1248,8 @@
     Int_t nrow, ncolumn, non;
 
-
-    alpha_iterator = fHistAlpha.begin();
-
-    while( alpha_iterator != fHistAlpha.end() ) {
+    TH1D *alpha_iterator = NULL;
+    TOrdCollectionIter Next(fHistAlpha);
+
+    while( (alpha_iterator = (TH1D*)Next())) {
 
           nrow    = index/fHistOn.GetNbinsX() + 1;                    // row of the histogram (y)
@@ -1182,8 +1266,7 @@
 			xpos, ypos, signif, nex, non);
            
-	  (*alpha_iterator).SetName(strname);
-	  (*alpha_iterator).SetTitle(strtitle);
-          (*alpha_iterator).Write();
-	  alpha_iterator++;
+	  alpha_iterator->SetName(strname);
+	  alpha_iterator->SetTitle(strtitle);
+          alpha_iterator->Write();
           index++;
     }
@@ -1192,140 +1275,5 @@
 }
 
-// --------------------------------------------------------------------------
-//
-// Draw the histogram
-//
-void MSkyPlot::Draw(Option_t *opt)
-{
-/*
-    TVirtualPad *pad = gPad ? gPad : MakeDefCanvas(this);
-    pad->SetBorderMode(0);
-
-    AppendPad("");
-
-    pad->Divide(1, 2, 0, 0.03);
-
-    TObject *catalog = GetCatalog();
-
-    // Initialize upper part
-    pad->cd(1);
-    gPad->SetBorderMode(0);
-    gPad->Divide(3, 1);
-*/
-/*
-    // PAD #1
-    pad->GetPad(1)->cd(1);
-    gPad->SetBorderMode(0);
-    fHist.GetZaxis()->SetRangeUser(0,fAlphaONMax);
-    TH1 *h3 = fHist.Project3D("yx_on");
-    fHist.GetZaxis()->SetRange(0,9999);
-    h3->SetDirectory(NULL);
-    h3->SetXTitle(fHist.GetXaxis()->GetTitle());
-    h3->SetYTitle(fHist.GetYaxis()->GetTitle());
-    h3->Draw("colz");
-    h3->SetBit(kCanDelete);
-    catalog->Draw("mirror same");
-
-    // PAD #2
-    pad->GetPad(1)->cd(2);
-    gPad->SetBorderMode(0);
-    fHist.GetZaxis()->SetRange(0,0);
-    TH1 *h4 = fHist.Project3D("yx_sig"); // Do this to get the correct binning....
-    fHist.GetZaxis()->SetRange(0,9999);
-    h4->SetTitle("Significance");
-    h4->SetDirectory(NULL);
-    h4->SetXTitle(fHist.GetXaxis()->GetTitle());
-    h4->SetYTitle(fHist.GetYaxis()->GetTitle());
-    h4->Draw("colz");
-    h4->SetBit(kCanDelete);
-    catalog->Draw("mirror same");
-
-    // PAD #3
-    pad->GetPad(1)->cd(3);
-    gPad->SetBorderMode(0);
-    fHist.GetZaxis()->SetRangeUser(fBgMean-fAlphaCut/2, fBgMean+fAlphaCut/2);
-    TH1 *h2 = fHist.Project3D("yx_off");
-    h2->SetDirectory(NULL);
-    h2->SetXTitle(fHist.GetXaxis()->GetTitle());
-    h2->SetYTitle(fHist.GetYaxis()->GetTitle());
-    h2->Draw("colz");
-    h2->SetBit(kCanDelete);
-    catalog->Draw("mirror same");
-
-
-    // Initialize lower part
-    pad->cd(2);
-    gPad->SetBorderMode(0);
-    gPad->Divide(3, 1);
-
-    // PAD #4
-    pad->GetPad(2)->cd(1);
-    gPad->SetBorderMode(0);
-    TH1 *h1 = fHist.ProjectionZ("Alpha");
-    h1->SetDirectory(NULL);
-    h1->SetTitle("Distribution of \\alpha");
-    h1->SetXTitle(fHist.GetZaxis()->GetTitle());
-    h1->SetYTitle("Counts");
-    h1->Draw(opt);
-    h1->SetBit(kCanDelete);
-
-    // PAD #5
-    pad->GetPad(2)->cd(2);
-    gPad->SetBorderMode(0);
-    TH1 *h5 = (TH1*)h3->Clone("Alpha_yx_diff");
-    h5->Add(h2, -1);
-    h5->SetTitle("Difference of on- and off-distribution");
-    h5->SetDirectory(NULL);
-    h5->Draw("colz");
-    h5->SetBit(kCanDelete);
-    catalog->Draw("mirror same");
-
-    // PAD #6
-    pad->GetPad(2)->cd(3);
-    gPad->SetBorderMode(0);
-    TH1 *h0 = fHist.Project3D("yx_all");
-    h0->SetDirectory(NULL);
-    h0->SetXTitle(fHist.GetXaxis()->GetTitle());
-    h0->SetYTitle(fHist.GetYaxis()->GetTitle());
-    h0->Draw("colz");
-    h0->SetBit(kCanDelete);
-    catalog->Draw("mirror same");
-*/
-}
-
-// --------------------------------------------------------------------------
-//
-// Everything which is in the main pad belongs to this class!
-//
-Int_t MSkyPlot::DistancetoPrimitive(Int_t px, Int_t py)
-{
-    return 0;
-}
-
-// --------------------------------------------------------------------------
-//
-// Set all sub-pads to Modified()
-//
-/*
-void MSkyPlot::Modified()
-{
-    if (!gPad)
-        return;
-
-    TVirtualPad *padsave = gPad;
-    padsave->Modified();
-    padsave->GetPad(1)->cd(1);
-    gPad->Modified();
-    padsave->GetPad(1)->cd(3);
-    gPad->Modified();
-    padsave->GetPad(2)->cd(1);
-    gPad->Modified();
-    padsave->GetPad(2)->cd(2);
-    gPad->Modified();
-    padsave->GetPad(2)->cd(3);
-    gPad->Modified();
-    gPad->cd();
-}
-*/
+
 // --------------------------------------------------------------------------
 //
