Index: trunk/MagicSoft/Mars/mtemp/mifae/macros/DispSkymap.C
===================================================================
--- trunk/MagicSoft/Mars/mtemp/mifae/macros/DispSkymap.C	(revision 6043)
+++ trunk/MagicSoft/Mars/mtemp/mifae/macros/DispSkymap.C	(revision 6044)
@@ -9,5 +9,9 @@
 //************************************************************************
 
-void DispSkymap(TString filename = "hillasfile.root")
+void DispSkymap(TString onfilename  = "ONhillasfile.root",
+		TString offfilename = "OFFhillasfile.root",
+		TString outfilename = "skymaps.root",
+		Float_t HadronnessCut = 0.06, 
+		Float_t SizeCut = 1000.)
 {
   //======================================================================
@@ -20,13 +24,22 @@
   
   // Input root file/s storing the computed Hillas and Disp parameters
-  gLog << "Input file/s: " <<  filename << endl;
-  
+  gLog << "Input ON file/s: "  <<  onfilename  << endl;
+  gLog << "Input OFF file/s: " <<  offfilename << endl;
+  
+
   //----------------------------------------------------
-  MTaskList tlist;
-  MParList  plist;
-  
-  //      MReadMarsFile read("Events", filename);
-  MReadTree read("Parameters", filename);
-  read.DisableAutoScheme();
+  MTaskList tliston;
+  MParList  pliston;
+
+  MTaskList tlistoff;
+  MParList  plistoff;
+  
+  // MReadMarsFile readon("Events", onfilename);
+  MReadTree readon("Parameters", onfilename);
+  readon.DisableAutoScheme();
+
+  // MReadMarsFile readoff("Events", offfilename);
+  MReadTree readoff("Parameters", offfilename);
+  readoff.DisableAutoScheme();
   
   MGeomCamMagic geomcam;
@@ -37,5 +50,21 @@
   //      fdisp->SetCuts(0,1,7,600,0,600,0.,3000000.,0.,0.,0.,0.,0.,0.);
   //      MContinue contdisp(fdisp);
-  
+
+  // Hadroness cut
+  TString hadrcut = "MHadronness.fHadronness<";
+  hadrcut += HadronnessCut;
+  gLog << "Hadronness Cut applied = " << hadrcut << endl;
+  MF hadrfilter(hadrcut);  
+  MContinue conthadr(&hadrfilter);
+  conthadr.SetInverted(kTRUE);
+
+  // Size cut
+  TString sizecut = "MHillas.fSize>";
+  sizecut += SizeCut;
+  gLog << "Size Cut applied = " << sizecut << endl;
+  MF sizefilter(sizecut);  
+  MContinue contsize(&sizefilter);
+  contsize.SetInverted(kTRUE);
+
   // make Disp plots
   // SelectedPos = 1  means choose the right source position
@@ -43,81 +72,141 @@
   //               3               the position according to M3Long
   //               4               the position according to Asym
-  MHDisp hdisp1;
-  hdisp1.SetName("MHDispCorr");
-  hdisp1.SetSelectedPos(1);
-  MFillH filldisp1("MHDispCorr[MHDisp]", "");
-  
-  MHDisp hdisp2;
-  hdisp2.SetName("MHDispWrong");
-  hdisp2.SetSelectedPos(2);
-  MFillH filldisp2("MHDispWrong[MHDisp]", "");
-  
-  MHDisp hdisp3;
-  hdisp3.SetName("MHDispM3Long");
-  hdisp3.SetSelectedPos(3);
-  MFillH filldisp3("MHDispM3Long[MHDisp]", "");
-  
-  MHDisp hdisp4;
-  hdisp4.SetName("MHDispAsym");
-  hdisp4.SetSelectedPos(4);
-  MFillH filldisp4("MHDispAsym[MHDisp]", "");
+  MHDisp hdispon;
+  hdispon.SetName("MHDispAsymOn");
+  hdispon.SetSelectedPos(4);
+  MFillH filldispon("MHDispAsymOn[MHDisp]", "");
+
+  MHDisp hdispoff;
+  hdispoff.SetName("MHDispAsymOff");
+  hdispoff.SetSelectedPos(4);
+  MFillH filldispoff("MHDispAsymOff[MHDisp]", "");
   
   
   //*****************************
-  // entries in MParList
-  plist.AddToList(&tlist);
-  plist.AddToList(&geomcam);
-  plist.AddToList(&hdisp1);
-  plist.AddToList(&hdisp2);
-  plist.AddToList(&hdisp3);
-  plist.AddToList(&hdisp4);
+  // entries in MParList On
+  pliston.AddToList(&tliston);
+  pliston.AddToList(&geomcam);
+  pliston.AddToList(&hdispon);
+
+  // entries in MParList Off
+  plistoff.AddToList(&tlistoff);
+  plistoff.AddToList(&geomcam);
+  plistoff.AddToList(&hdispoff);
+
   
   //*****************************
-  // entries in MTaskList
-  tlist.AddToList(&read);
+  // entries in MTaskList On
+  tliston.AddToList(&readon);
+  if (fdisp != NULL)
+    tliston.AddToList(&contdisp);
+  tliston.AddToList(&conthadr);
+  tliston.AddToList(&contsize);
+  tliston.AddToList(&filldispon);
+
+  // entries in MTaskList Off
+  tlistoff.AddToList(&readoff);
   if (fdisp != NULL)
     tlist.AddToList(&contdisp);
-  tlist.AddToList(&filldisp1);
-  //  tlist.AddToList(&filldisp2);
-  tlist.AddToList(&filldisp3);
-  tlist.AddToList(&filldisp4);
+  tlistoff.AddToList(&conthadr);
+  tlistoff.AddToList(&contsize);
+  tlistoff.AddToList(&filldispoff);
   
   //*****************************
   
   //-------------------------------------------
-  // Execute event loop
-  //
-  MProgressBar bar;
-  MEvtLoop evtloop;
-  evtloop.SetParList(&plist);
-  evtloop.SetProgressBar(&bar);
-  
-  Int_t maxevents = -1;
-  if ( !evtloop.Eventloop(maxevents) )
+  // Execute event loop On
+  MProgressBar barOn;
+  MEvtLoop evtloopOn;
+  evtloopOn.SetParList(&pliston);
+  evtloopOn.SetProgressBar(&barOn);
+  
+  Int_t maxeventsOn = -1;
+  if ( !evtloopOn.Eventloop(maxeventsOn) )
     return;
   
-  tlist.PrintStatistics(0, kTRUE);
-  
-  //-------------------------------------------
-  // Display the histograms
-  //
-  //  hdisp1.DrawClone();
-  //  hdisp2.DrawClone();
-  //  hdisp3.DrawClone();
-  //  hdisp4.DrawClone();
-  
-  gLog << "Drawing DISP Skymap for the FoV (srcpos solution selected according Asym sign)..." << endl; 
-  TCanvas *c = new TCanvas("c","Disp Skymap",0,0,900,900);    
+  tliston.PrintStatistics(0, kTRUE);
+
+  //-------------------------------------------
+  // Execute event loop Off
+  MProgressBar barOff;
+  MEvtLoop evtloopOff;
+  evtloopOff.SetParList(&plistoff);
+  evtloopOff.SetProgressBar(&barOff);
+  
+  Int_t maxeventsOff = -1;
+  if ( !evtloopOff.Eventloop(maxeventsOff) )
+    return;
+  
+  tlistoff.PrintStatistics(0, kTRUE);
+  
+
+  //-------------------------------------------
+  // Get the skymaps
+  TH2F *skymapOn  = hdispon.GetSkymapXY();
+  skymapOn->SetName("fSkymapXYOn");
+  skymapOn->SetTitle("ON DATA: Disp estimated source positions Skymap");
+  skymapOn->SetTitleOffset(1.2,"Y");
+  TH2F *skymapOff = hdispoff.GetSkymapXY();
+  skymapOff->SetName("fSkymapXYOff");
+  skymapOff->SetTitle("OFF DATA: Disp estimated source positions Skymap");
+  skymapOff->SetTitleOffset(1.2,"Y");
+
+  // Normalize to the number of events and subtract Off to On skymap
+  Double_t onentries  = skymapOn->GetEntries();
+  Double_t offentries = skymapOff->GetEntries();
+  Double_t norm = onentries/offentries;
+  cout << "Number of ON events after cuts = "  << onentries  << endl;
+  cout << "Number of OFF events after cuts = " << offentries << endl;
+  cout << "Normalization factor = " << norm << endl;
+
+  TH2F *skymap = new TH2F("fSkymapXY", 
+         "ON-OFF: Disp estimated source positions Skymap", 71, -2., 2., 71, -2., 2.);
+  skymap->Add(skymapOn,skymapOff,1.,-norm);
+  skymap->SetTitleOffset(1.2,"Y");
+  
+  //-------------------------------------------
+  // Display the skymaps
+  gLog << endl;
+  gLog << "Drawing DISP Skymaps for the FoV ...... "   << endl;
+  gLog << "(srcpos solution selected according Asym sign)" << endl; 
+
+  gStyle->SetPalette(1);
+  gStyle->SetOptStat(11);
+
+  TCanvas *c = new TCanvas("c","Disp Skymaps",0,0,900,900);    
   c->SetBorderMode(0);
-  gStyle->SetPalette(1);
-  TH2F *skymap = hdisp4.GetSkymapXY();
-  skymap->SetTitleOffset(1.2,"Y");
-  //  skymap->SetStats(0);
+  c->Divide(2,2);
+
+  c->cd(1);
+  gPad->SetBorderMode(0);
+  skymapOn->DrawClone("COLZ");
+  skymapOn->SetBit(kCanDelete);
+  gPad->Modified();
+
+  c->cd(2);
+  gPad->SetBorderMode(0);
+  skymapOff->DrawClone("COLZ");
+  skymapOff->SetBit(kCanDelete);
+  gPad->Modified();
+
+  c->cd(3);
+  gPad->SetBorderMode(0);
   skymap->DrawClone("COLZ");
   skymap->SetBit(kCanDelete);
-
-
-  //-------------------------------------------
-  gLog << endl << "Disp plots were made for file '" << filename << endl; 
+  gPad->Modified();
+
+  //-------------------------------------------
+  // Save the skymaps in a .root file
+  TFile outfile(outfilename,"RECREATE");
+  skymapOn->Write();
+  skymapOff->Write();
+  skymap->Write();
+  outfile.Close();
+  gLog << endl << "Skymaps stored in file: " << outfilename <<endl;
+
+  //-------------------------------------------
+  gLog << endl << "Disp plots were made for files: " << endl;
+  gLog << "ON data:  " << onfilename  << endl; 
+  gLog << "OFF data: " << offfilename << endl; 
   gLog << "-----------------------------------------------------" << endl;