source: trunk/MagicSoft/Mars/mtemp/mifae/macros/DispSkymap.C@ 6723

//************************************************************************
//
// Authors : Eva Domingo,     1/2005 <mailto:domingo@ifae.es>
//
//
// Macro for generating DISP Skymap of the FoV
// -------------------------------------------
//
//************************************************************************

void DispSkymap(TString onfilename  = "ONhillasfile.root",
                TString offfilename = "OFFhillasfile.root",
                TString outfilename = "skymaps.root",
                Float_t HadronnessCut = 0.06, 
                Float_t SizeCut = 1000.)
{
  //======================================================================
  // Make Disp plots
  //======================================================================
  
  gLog << "-----------------------------------------------------" << endl; 
  gLog << "Make Disp related plots" << endl;
  gLog << "-----------------------------------------------------" << endl; 
  
  // Input root file/s storing the computed Hillas and Disp parameters
  gLog << "Input ON file/s: "  <<  onfilename  << endl;
  gLog << "Input OFF file/s: " <<  offfilename << endl;
  

  //----------------------------------------------------
  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;

  // set cuts to select an event sample to apply Disp
  MFDisp *fdisp = NULL;
  //      fdisp = new MFDisp;
  //      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
  //               2                   wrong
  //               3               the position according to M3Long
  //               4               the position according to Asym
  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 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 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);
  tlistoff.AddToList(&conthadr);
  tlistoff.AddToList(&contsize);
  tlistoff.AddToList(&filldispoff);
  
  //*****************************
  
  //-------------------------------------------
  // Execute event loop On
  MProgressBar barOn;
  MEvtLoop evtloopOn;
  evtloopOn.SetParList(&pliston);
  evtloopOn.SetProgressBar(&barOn);
  
  Int_t maxeventsOn = -1;
  if ( !evtloopOn.Eventloop(maxeventsOn) )
    return;
  
  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);
  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);
  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; 

}