/*********************************/
/* Compute the hillas parameters */
/*********************************/

#include "TString.h"
#include "TArrayS.h"

#include "MParList.h"
#include "MTaskList.h"
#include "MPedestalCam.h"
#include "MBadPixelsCam.h"
#include "MReadMarsFile.h"
#include "MGeomApply.h"
#include "MPedCalcPedRun.h"
#include "MEvtLoop.h"
#include "MGeomCamMagic.h"
#include "MExtractedSignalCam.h"
#include "MCalibrationChargeCam.h"
#include "MHCalibrationChargeCam.h"
#include "MHCalibrationRelTimeCam.h"
#include "MExtractSignal.h"
#include "MCalibrationChargeCalc.h"
#include "MFCosmics.h"
#include "MContinue.h"
#include "MFillH.h"
#include "MLog.h"
#include "MCerPhotEvt.h"
#include "MPedPhotCam.h"
#include "MCalibrate.h"
#include "MPedPhotCalc.h"
#include "MHillas.h"
#include "MRawRunHeader.h"
#include "MSrcPosCam.h"
#include "MBlindPixelCalc.h"
#include "MImgCleanStd.h"
#include "MHillasSrcCalc.h"
#include "MHillasCalc.h"
#include "MWriteRootFile.h"
#include "MProgressBar.h"
#include "MArgs.h"
#include "MRunIter.h"
#include "MJPedestal.h"
#include "MJCalibration.h"

#include <iostream>
#include <fstream>
#include <stdlib.h>

using namespace std;

Bool_t readDatacards(TString& filename);
void makeHillas();

// initial and final time slices to be used in signal extraction
const Byte_t hifirst = 1;
const Byte_t hilast  = 14;
const Byte_t lofirst = 4;
const Byte_t lolast  = 13;

// declaration of variables read from datacards
TString  outname;
TString  idirname;
MRunIter caliter;
MRunIter pediter;
MRunIter datiter;
ULong_t  nmaxevents=999999999;
Short_t  calflag=1;
Float_t  lcore = 3.0;
Float_t  ltail = 1.5;
Int_t    nfiles = 0;

const TString defaultcard="input.datacard";
static void Usage()
{
  gLog <<endl;
  gLog << "Usage is:" << endl;
  gLog << "   makeHillas [-h] [-?] <datacards>" << endl << endl;
  gLog << "     <datacards>: datacards file name (dafault input.datacards)" << endl;
  gLog << "     -?/-h: This help" << endl << endl;
}

int main(int argc, char **argv)
{
  // evaluate arguments
  MArgs arg(argc, argv);
  if (arg.HasOnly("-?") || arg.HasOnly("-h") || arg.HasOnly("--help"))
    {
      Usage();
      return -1;
    }

  TString datacard  = arg.GetArgumentStr(0);
  if(!datacard.Length())
    datacard = defaultcard;

  if(!readDatacards(datacard))
    {
      cout << "Error reading datacards. Stoping" << endl;
      return -1;
    }
  makeHillas();
}

/*************************************************************/
void makeHillas()
{
#if 0
  MStatusDisplay *display = new MStatusDisplay;
  display->SetUpdateTime(3000);
  display->Resize(850,700);
  
  gStyle->SetOptStat(1111);
  gStyle->SetOptFit();
#endif
  
  /************************************/
  /* FIRST LOOP: PEDESTAL COMPUTATION */
  /************************************/
  
  MParList plist1;
  MTaskList tlist1;
  plist1.AddToList(&tlist1);
  
  // containers
  MPedestalCam   pedcam;
  MBadPixelsCam  badcam;  
  // badcam.AsciiRead("badpixels.dat");

  plist1.AddToList(&pedcam);
  plist1.AddToList(&badcam);

  //tasks
  MReadMarsFile read("Events");
  static_cast<MRead&>(read).AddFiles(pediter); 
  read.DisableAutoScheme();
  
  MGeomApply     geomapl;
  MPedCalcPedRun pedcalc;

  tlist1.AddToList(&read);
  tlist1.AddToList(&geomapl);
  tlist1.AddToList(&pedcalc);

  // Create and setup the eventloop
  MEvtLoop pedloop;
  pedloop.SetParList(&plist1);
  //pedloop.SetDisplay(display);

  cout << "*************************" << endl;
  cout << "** COMPUTING PEDESTALS **" << endl;
  cout << "*************************" << endl;

  if (!pedloop.Eventloop())
    return;
  
  tlist1.PrintStatistics();

  /*****************************/
  /* SECOND LOOP: CALIBRATION  */
  /*****************************/        

  MParList  plist2;
  MTaskList tlist2;
  plist2.AddToList(&tlist2);
  plist2.AddToList(&pedcam);
  plist2.AddToList(&badcam);

  // new containers
  MGeomCamMagic          geomcam;
  MExtractedSignalCam    sigcam;
  MCalibrationChargeCam  calcam;
  MHCalibrationChargeCam histcharge;
  MHCalibrationRelTimeCam     histtime;


  plist2.AddToList(&geomcam);
  plist2.AddToList(&sigcam);
  plist2.AddToList(&calcam);
  plist2.AddToList(&histcharge);

  //tasks
  MReadMarsFile read2("Events");
  static_cast<MRead&>(read2).AddFiles(caliter); 
  read2.DisableAutoScheme();     
  
  MExtractSignal         sigcalc;
  sigcalc.SetRange(hifirst,hilast,lofirst,lolast);
  MCalibrationChargeCalc calcalc;
  MFCosmics              cosmics;
  MContinue              cont(&cosmics);

  
  calcalc.SetPulserColor(MCalibrationCam::kCT1);
  MFillH fillcam  ("MHCalibrationChargeCam"     , "MExtractedSignalCam");

  tlist2.AddToList(&read2);
  tlist2.AddToList(&geomapl);
  tlist2.AddToList(&sigcalc);
  tlist2.AddToList(&cont);
  tlist2.AddToList(&fillcam);
  tlist2.AddToList(&calcalc);
  
  // Create and setup the eventloop
  MEvtLoop calloop;
  calloop.SetParList(&plist2);
  //calloop.SetDisplay(display);

  cout << "***************************" << endl;
  cout << "** COMPUTING CALIBRATION **" << endl;
  cout << "***************************" << endl;
  
  if (!calloop.Eventloop())
    return;
  
  tlist2.PrintStatistics();
  
  MLog gauglog;
  gauglog.SetOutputFile(Form("%s%s",calcam.GetName(),".txt"),1);
  calcam.SetLogStream(&gauglog);
  calcam.Print();
  calcam.SetLogStream(&gLog);
#if 0
  // Create histograms to display
  MHCamera disp1  (geomcam, "Cal;Charge",         "Fitted Mean Charges");
  MHCamera disp2  (geomcam, "Cal;SigmaCharge",    "Sigma of Fitted Charges");
  MHCamera disp3  (geomcam, "Cal;FitProb",        "Probability of Fit");
  MHCamera disp4  (geomcam, "Cal;RSigma",         "Reduced Sigmas");
  MHCamera disp5  (geomcam, "Cal;RSigma/Charge",  "Reduced Sigma per Charge");
  MHCamera disp6  (geomcam, "Cal;FFactorPh",      "Nr. of Photo-electrons (F-Factor Method)");
  MHCamera disp7  (geomcam, "Cal;FFactorConv",    "Conversion Factor to photons (F-Factor Method)");
  MHCamera disp8  (geomcam, "Cal;FFactorFFactor", "Total F-Factor (F-Factor Method)");
  MHCamera disp9  (geomcam, "Cal;BlindPixPh",     "Photon flux inside plexiglass (Blind Pixel Method)");
  MHCamera disp10 (geomcam, "Cal;BlindPixConv",   "Conversion Factor to photons (Blind Pixel Method)");
  MHCamera disp11 (geomcam, "Cal;BlindPixFFactor","Total F-Factor (Blind Pixel Method)");
  MHCamera disp12 (geomcam, "Cal;PINDiodePh",     "Photon flux outside plexiglass (PIN Diode Method)");
  MHCamera disp13 (geomcam, "Cal;PINDiodeConv",   "Conversion Factor tp photons (PIN Diode Method)");
  MHCamera disp14 (geomcam, "Cal;PINDiodeFFactor","Total F-Factor (PIN Diode Method)");
  MHCamera disp15 (geomcam, "Cal;Excluded",       "Pixels previously excluded");
  MHCamera disp16 (geomcam, "Cal;NotFitted",      "Pixels that could not be fitted");
  MHCamera disp17 (geomcam, "Cal;NotFitValid",    "Pixels with not valid fit results");
  MHCamera disp18 (geomcam, "Cal;HiGainOscillating",     "Oscillating Pixels HI Gain");
  MHCamera disp19 (geomcam, "Cal;LoGainOscillating",     "Oscillating Pixels LO Gain");
  MHCamera disp20 (geomcam, "Cal;HiGainPickup",          "Number Pickup events Hi Gain");
  MHCamera disp21 (geomcam, "Cal;LoGainPickup",          "Number Pickup events Lo Gain");
  MHCamera disp22 (geomcam, "Cal;Saturation",            "Pixels with saturated Hi Gain");
  MHCamera disp23 (geomcam, "Cal;FFactorValid",          "Pixels with valid F-Factor calibration");
  MHCamera disp24 (geomcam, "Cal;BlindPixelValid",       "Pixels with valid BlindPixel calibration");
  MHCamera disp25 (geomcam, "Cal;PINdiodeFFactorValid",  "Pixels with valid PINDiode calibration");
  
  MHCamera disp26 (geomcam, "Cal;Ped",         "Pedestals");
  MHCamera disp27 (geomcam, "Cal;PedRms",      "Pedestal RMS");
  
  MHCamera disp28 (geomcam, "time;Time",        "Rel. Arrival Times");
  MHCamera disp29 (geomcam, "time;SigmaTime",   "Sigma of Rel. Arrival Times");
  MHCamera disp30 (geomcam, "time;TimeProb",    "Probability of Time Fit");
  MHCamera disp31 (geomcam, "time;NotFitValid", "Pixels with not valid fit results");
  MHCamera disp32 (geomcam, "time;Oscillating", "Oscillating Pixels");
  
  MHCamera disp33 (geomcam, "Cal;AbsTimeMean", "Abs. Arrival Times");
  MHCamera disp34 (geomcam, "Cal;AbsTimeRms",  "RMS of Arrival Times");

  // Fitted charge means and sigmas
  disp1.SetCamContent(calcam,  0);
  disp1.SetCamError(  calcam,  1);
  disp2.SetCamContent(calcam,  2);
  disp2.SetCamError(  calcam,  3);

  // Fit probabilities
  disp3.SetCamContent(calcam,  4);
  
  // Reduced Sigmas and reduced sigmas per charge
  disp4.SetCamContent(calcam,  5);
  disp4.SetCamError(  calcam,  6);
  disp5.SetCamContent(calcam,  7);
  disp5.SetCamError(  calcam,  8);
  
  // F-Factor Method
  disp6.SetCamContent(calcam,  9);
  disp6.SetCamError(  calcam, 10);
  disp7.SetCamContent(calcam, 11);
  disp7.SetCamError(  calcam, 12);
  disp8.SetCamContent(calcam, 13);
  disp8.SetCamError(  calcam, 14);
  
  // Blind Pixel Method
  disp9.SetCamContent(calcam, 15);
  disp9.SetCamError(  calcam, 16);
  disp10.SetCamContent(calcam,17);
  disp10.SetCamError(  calcam,18);
  disp11.SetCamContent(calcam,19);
  disp11.SetCamError(  calcam,20);
  
  // PIN Diode Method
  disp12.SetCamContent(calcam,21);
  disp12.SetCamError(  calcam,22);
  disp13.SetCamContent(calcam,23);
  disp13.SetCamError(  calcam,24);
  disp14.SetCamContent(calcam,25);
  disp14.SetCamError(  calcam,26);
  
  // Pixels with defects
  disp15.SetCamContent(calcam,27);
  disp16.SetCamContent(calcam,28);
  disp17.SetCamContent(badcam,9);
  disp18.SetCamContent(badcam,16);
  disp19.SetCamContent(badcam,15);
  disp20.SetCamContent(calcam,29);
  disp21.SetCamContent(calcam,30);
  
  // Lo Gain calibration
  disp22.SetCamContent(calcam,31);
  
  // Valid flags
  disp23.SetCamContent(calcam,32);
  disp24.SetCamContent(calcam,33);
  disp25.SetCamContent(calcam,34);
  
  // Pedestals
  disp26.SetCamContent(calcam,35);
  disp26.SetCamError(  calcam,36);
  disp27.SetCamContent(calcam,37);
  disp27.SetCamError(  calcam,38);
  
  // Relative Times
  disp28.SetCamContent(histtime,0);
  disp28.SetCamError(  histtime,1);
  disp29.SetCamContent(histtime,2);
  disp29.SetCamError(  histtime,3);
  disp30.SetCamContent(histtime,4);
  disp31.SetCamContent(histtime,5);
  disp32.SetCamContent(histtime,6);
  
  // Absolute Times
  disp33.SetCamContent(calcam,39);
  disp33.SetCamError(  calcam,40);
  disp34.SetCamContent(calcam,41);
  
  disp1.SetYTitle("Charge [FADC units]");
  disp2.SetYTitle("\\sigma_{Charge} [FADC units]");
  disp3.SetYTitle("P_{Charge} [1]");
  
  disp4.SetYTitle("\\sqrt{\\sigma^{2}_{Charge} - RMS^{2}_{Ped}} [FADC Counts]");
  disp5.SetYTitle("Reduced Sigma / Mean Charge [1]");
  
  disp6.SetYTitle("Nr. Photo-electrons [1]");
  disp7.SetYTitle("Conversion Factor [Ph/FADC Count]");
  disp8.SetYTitle("\\sqrt{N_{Ph}}*\\sigma_{Charge}/\\mu_{Charge} [1] ");
  
  disp9.SetYTitle("Photon flux [ph/mm^2]");
  disp10.SetYTitle("Conversion Factor [Phot/FADC Count]");
  disp11.SetYTitle("\\sqrt{N_{Ph}}*\\sigma_{Charge}/\\mu_{Charge} [1]");
  
  disp12.SetYTitle("Photon flux [ph/mm^2]");
  disp13.SetYTitle("Conversion Factor [Phot/FADC Count]");
  disp14.SetYTitle("\\sqrt{N_{Ph}}*\\sigma_{Charge}/\\mu_{Charge} [1]");
  
  disp15.SetYTitle("[1]");
  disp16.SetYTitle("[1]");
  disp17.SetYTitle("[1]");
  disp18.SetYTitle("[1]");
  disp19.SetYTitle("[1]");
  disp20.SetYTitle("[1]");
  disp21.SetYTitle("[1]");
  disp22.SetYTitle("[1]");
  disp23.SetYTitle("[1]");
  disp24.SetYTitle("[1]");
  disp25.SetYTitle("[1]");

  disp26.SetYTitle("Ped [FADC Counts ]");
  disp27.SetYTitle("RMS_{Ped} [FADC Counts ]");
  
  disp28.SetYTitle("Time Offset [ns]");
  disp29.SetYTitle("Timing resolution [ns]");
  disp30.SetYTitle("P_{Time} [1]");
  
  disp31.SetYTitle("[1]");
  disp32.SetYTitle("[1]");
  
  disp33.SetYTitle("Mean Abs. Time [FADC slice]");
  disp34.SetYTitle("RMS Abs. Time [FADC slices]");
  
  gStyle->SetOptStat(1111);
  gStyle->SetOptFit();
  
  // Charges
  TCanvas &c1 = display->AddTab("Fit.Charge");
  c1.Divide(2, 3);
  
  CamDraw(c1, disp1,calcam,1, 2 , 2);
  CamDraw(c1, disp2,calcam,2, 2 , 2);

  // Fit Probability
  TCanvas &c2 = display->AddTab("Fit.Prob");
  c2.Divide(1,3);

  CamDraw(c2, disp3,calcam,1, 1 , 4);
  
  // Reduced Sigmas
  TCanvas &c3 = display->AddTab("Red.Sigma");
  c3.Divide(2,3);
  
  CamDraw(c3, disp4,calcam,1, 2 , 2);
  CamDraw(c3, disp5,calcam,2, 2 , 2);

  // F-Factor Method
  TCanvas &c4 = display->AddTab("F-Factor");
  c4.Divide(3,3);
  
  CamDraw(c4, disp6,calcam,1, 3 , 2);
  CamDraw(c4, disp7,calcam,2, 3 , 2);
  CamDraw(c4, disp8,calcam,3, 3 , 2);
  
  // Blind Pixel Method
  TCanvas &c5 = display->AddTab("BlindPix");
  c5.Divide(3, 3);
  
  CamDraw(c5, disp9,calcam,1, 3 ,  9);
  CamDraw(c5, disp10,calcam,2, 3 , 2);
  CamDraw(c5, disp11,calcam,3, 3 , 2);
  
  // PIN Diode Method
  TCanvas &c6 = display->AddTab("PINDiode");
  c6.Divide(3,3);
  
  CamDraw(c6, disp12,calcam,1, 3 , 9);
  CamDraw(c6, disp13,calcam,2, 3 , 2);
  CamDraw(c6, disp14,calcam,3, 3 , 2);
  
  // Defects
  TCanvas &c7 = display->AddTab("Defects");
  c7.Divide(4,2);
  
  CamDraw(c7, disp15,calcam,1,4, 0);
  CamDraw(c7, disp16,calcam,2,4, 0);
  CamDraw(c7, disp20,calcam,3,4, 0);
  CamDraw(c7, disp21,calcam,4,4, 0);
  
  // BadCam
  TCanvas &c8 = display->AddTab("Defects");
  c8.Divide(3,2);
  
  CamDraw(c8, disp17,badcam,1,3, 0);
  CamDraw(c8, disp18,badcam,2,3, 0);
  CamDraw(c8, disp19,badcam,3,3, 0);
  
  // Valid flags
  TCanvas &c9 = display->AddTab("Validity");
  c9.Divide(4,2);
  
  CamDraw(c9, disp22,calcam,1,4,0);
  CamDraw(c9, disp23,calcam,2,4,0);
  CamDraw(c9, disp24,calcam,3,4,0);
  CamDraw(c9, disp25,calcam,4,4,0);
  
  // Pedestals
  TCanvas &c10 = display->AddTab("Pedestals");
  c10.Divide(2,3);
  
  CamDraw(c10,disp26,calcam,1,2,1);
  CamDraw(c10,disp27,calcam,2,2,2);
  
  // Rel. Times
  TCanvas &c11 = display->AddTab("Fitted Rel. Times");
  c11.Divide(3,3);
  
  CamDraw(c11,disp28,calcam,1,3,2);
  CamDraw(c11,disp29,calcam,2,3,2);
  CamDraw(c11,disp30,calcam,3,3,4);
  
  // Time Defects
  TCanvas &c12 = display->AddTab("Time Def.");
  c12.Divide(2,2);
  
  CamDraw(c12, disp31,calcam,1,2, 0);
  CamDraw(c12, disp32,calcam,2,2, 0);
  
  // Abs. Times
  TCanvas &c13 = display->AddTab("Abs. Times");
  c13.Divide(2,3);
  
  CamDraw(c13,disp33,calcam,1,2,2);
  CamDraw(c13,disp34,calcam,2,2,2);
#endif  
  /************************************************************************/
  /*                THIRD LOOP: PEDESTAL CALIBRATION INTO PHOTONS         */
  /************************************************************************/
  MParList  plist3;
  MTaskList tlist3;
  plist3.AddToList(&tlist3);
  
  // containers
  MCerPhotEvt    nphot;
  MPedPhotCam    nphotrms;
  
  plist3.AddToList(&geomcam);
  plist3.AddToList(&pedcam);
  plist3.AddToList(&calcam);
  //  plist3.AddToList(&badcam);
  plist3.AddToList(&sigcam);
  plist3.AddToList(&nphot);
  plist3.AddToList(&nphotrms);

  
  MCalibrate::CalibrationMode_t calMode=MCalibrate::kDefault;  
  if(calflag==0)
    calMode=MCalibrate::kNone;
  

  //tasks
  MReadMarsFile read3("Events");
  static_cast<MRead&>(read3).AddFiles(pediter); 
  read3.DisableAutoScheme();
  
  MExtractSignal  extsig;
  extsig.SetRange(hifirst,hilast,lofirst,lolast);
  MCalibrate      photcalc(calMode);
  MPedPhotCalc    photrmscalc; 
  
  tlist3.AddToList(&read3);
  tlist3.AddToList(&geomapl);
  tlist3.AddToList(&extsig);
  tlist3.AddToList(&photcalc);
  tlist3.AddToList(&photrmscalc);

  // Create and setup the eventloop
  MEvtLoop evtloop3;
  evtloop3.SetParList(&plist3);
  if (!evtloop3.Eventloop())
    return;

  tlist3.PrintStatistics();
  
  /************************************************************************/
  /*                FOURTH LOOP: DATA CALIBRATION INTO PHOTONS            */
  /************************************************************************/

  MParList  plist4;
  MTaskList tlist4;
  plist4.AddToList(&tlist4);
  
  // containers 
  MHillas       hillas;
  MSrcPosCam    source;
  MRawRunHeader runhead;
  
  plist4.AddToList(&geomcam);
  plist4.AddToList(&pedcam);
  plist4.AddToList(&calcam);
  //  plist4.AddToList(&badcam);
  plist4.AddToList(&nphot);
  plist4.AddToList(&nphotrms);
  plist4.AddToList(&source);
  plist4.AddToList(&hillas);
  plist4.AddToList(&runhead);
  
  //tasks
  MReadMarsFile read4("Events");
  static_cast<MRead&>(read4).AddFiles(datiter); 
  read4.DisableAutoScheme();
  
  // set bad pixels 
  MBlindPixelCalc   blind;
  MBlindPixelCalc   blind2;
  const Short_t x[16] = {330,395,329,396,389,
                         323,388,322,384,385,
                         386,387,321,320,319,
                         394};
  const TArrayS bp(16,(Short_t*)x);
  blind.SetPixelIndices(bp);
  blind2.SetPixelIndices(bp);
  
  MImgCleanStd      clean(lcore,ltail);
  MHillasCalc       hcalc;
  MHillasSrcCalc    csrc1;
  
  MWriteRootFile write(outname,"RECREATE");

  //  write.AddContainer("MGeomCam"              , "RunHeaders");
  //  write.AddContainer("MRawRunHeader"         , "RunHeaders");
  //  write.AddContainer("MSrcPosCam"            , "RunHeaders");
  //  write.AddContainer("MCalibrationChargeCam" , "RunHeaders");
  //  write.AddContainer("MPedPhotCam"           , "RunHeaders"); // Attention, was in Events - Tree!!
  //  write.AddContainer("MPedestalCam"          , "RunHeaders");
  //  write.AddContainer("MHCalibrationRelTimeCam","RunHeaders");

  //  write.AddContainer("MCerPhotEvt"   , "Events");
  //  write.AddContainer("MRawEvtHeader" , "Events");
  //  write.AddContainer("MBadPixelsCam" , "Events");
  //  write.AddContainer("MPedPhotCam"   , "Events");

  write.AddContainer("MHillas"       , "Parameters");
  write.AddContainer("MHillasSrc"    , "Parameters");
  write.AddContainer("MHillasExt"    , "Parameters");
  write.AddContainer("MNewImagePar"  , "Parameters");
  write.AddContainer("MRawEvtHeader" , "Parameters");
  write.AddContainer("MRawRunHeader" , "Parameters");
  write.AddContainer("MConcentration" , "Parameters");
  
  tlist4.AddToList(&read4);
  tlist4.AddToList(&geomapl);
  tlist4.AddToList(&extsig);
  tlist4.AddToList(&photcalc);
  //tlist4.AddToList(&blind);
  tlist4.AddToList(&clean);
  //tlist4.AddToList(&blind2);
  tlist4.AddToList(&hcalc);
  //  tlist4.AddToList(&srcposcalc);
  tlist4.AddToList(&csrc1);
  tlist4.AddToList(&write);

  // Create and setup the eventloop
  MEvtLoop datloop;
  datloop.SetParList(&plist4);
  //  MProgressBar bar;
  //  datloop.SetProgressBar(&bar);

  cout << "*************************************************************" << endl;
  cout << "***   COMPUTING DATA USING EXTRACTED SIGNAL (IN PHOTONS)  ***" << endl;
  cout << "*************************************************************" << endl;
  
  if (!datloop.Eventloop(nmaxevents))
    return;
  
  tlist4.PrintStatistics();    
  nphotrms.Print();

}

#if 0
void CamDraw(TCanvas &c, MHCamera &cam, MCamEvent &evt, Int_t i, Int_t j, Int_t fit)
{

  c.cd(i);
  gPad->SetBorderMode(0);
  MHCamera *obj1=(MHCamera*)cam.DrawCopy("hist");
  //  obj1->AddNotify(evt);
  
  c.cd(i+j);
  gPad->SetBorderMode(0);
  obj1->Draw();
  ((MHCamera*)obj1)->SetPrettyPalette();

  if (fit != 0)
    {
      c.cd(i+2*j);
      gPad->SetBorderMode(0);
      TH1D *obj2 = (TH1D*)obj1->Projection(obj1.GetName());
      
//      obj2->Sumw2();
      obj2->Draw();
      obj2->SetBit(kCanDelete);

      const Double_t min   = obj2->GetBinCenter(obj2->GetXaxis()->GetFirst());
      const Double_t max   = obj2->GetBinCenter(obj2->GetXaxis()->GetLast());
      const Double_t integ = obj2->Integral("width")/2.5066283;
      const Double_t mean  = obj2->GetMean();
      const Double_t rms   = obj2->GetRMS();
      const Double_t width = max-min;

      if (rms == 0. || width == 0. )
        return;
      
      switch (fit)
        {
        case 1:
          TF1 *sgaus = new TF1("sgaus","gaus(0)",min,max);
          sgaus->SetBit(kCanDelete);
          sgaus->SetParNames("Area","#mu","#sigma");
          sgaus->SetParameters(integ/rms,mean,rms);
          sgaus->SetParLimits(0,0.,integ);
          sgaus->SetParLimits(1,min,max);
          sgaus->SetParLimits(2,0,width/1.5);
          obj2->Fit("sgaus","QLR");
          obj2->GetFunction("sgaus")->SetLineColor(kYellow);
          break;

        case 2:
          TString dgausform = "([0]-[3])/[2]*exp(-0.5*(x-[1])*(x-[1])/[2]/[2])";
          dgausform += "+[3]/[5]*exp(-0.5*(x-[4])*(x-[4])/[5]/[5])";
          TF1 *dgaus = new TF1("dgaus",dgausform.Data(),min,max);
          dgaus->SetBit(kCanDelete);
          dgaus->SetParNames("A_{tot}","#mu_{1}","#sigma_{1}","A_{2}","#mu_{2}","#sigma_{2}");
          dgaus->SetParameters(integ,(min+mean)/2.,width/4.,
                               integ/width/2.,(max+mean)/2.,width/4.);
          // The left-sided Gauss 
          dgaus->SetParLimits(0,integ-1.5,integ+1.5);
          dgaus->SetParLimits(1,min+(width/10.),mean);
          dgaus->SetParLimits(2,0,width/2.);
          // The right-sided Gauss 
          dgaus->SetParLimits(3,0,integ);
          dgaus->SetParLimits(4,mean,max-(width/10.));
          dgaus->SetParLimits(5,0,width/2.);
          obj2->Fit("dgaus","QLRM");
          obj2->GetFunction("dgaus")->SetLineColor(kYellow);
          break;
          
        case 3:
          TString tgausform = "([0]-[3]-[6])/[2]*exp(-0.5*(x-[1])*(x-[1])/[2]/[2])";
          tgausform += "+[3]/[5]*exp(-0.5*(x-[4])*(x-[4])/[5]/[5])";
          tgausform += "+[6]/[8]*exp(-0.5*(x-[7])*(x-[7])/[8]/[8])";
          TF1 *tgaus = new TF1("tgaus",tgausform.Data(),min,max);
          tgaus->SetBit(kCanDelete);
          tgaus->SetParNames("A_{tot}","#mu_{1}","#sigma_{1}",
                             "A_{2}","#mu_{2}","#sigma_{2}",
                             "A_{3}","#mu_{3}","#sigma_{3}");
          tgaus->SetParameters(integ,(min+mean)/2,width/4.,
                               integ/width/3.,(max+mean)/2.,width/4.,
                               integ/width/3.,mean,width/2.);
          // The left-sided Gauss 
          tgaus->SetParLimits(0,integ-1.5,integ+1.5);
          tgaus->SetParLimits(1,min+(width/10.),mean);
          tgaus->SetParLimits(2,width/15.,width/2.);
          // The right-sided Gauss 
          tgaus->SetParLimits(3,0.,integ);
          tgaus->SetParLimits(4,mean,max-(width/10.));
          tgaus->SetParLimits(5,width/15.,width/2.);
          // The Gauss describing the outliers
          tgaus->SetParLimits(6,0.,integ);
          tgaus->SetParLimits(7,min,max);
          tgaus->SetParLimits(8,width/4.,width/1.5);
          obj2->Fit("tgaus","QLRM");
          obj2->GetFunction("tgaus")->SetLineColor(kYellow);
          break;
        case 4:
          obj2->Fit("pol0","Q");
          obj2->GetFunction("pol0")->SetLineColor(kYellow);
          break;
        case 9:
          break;
        default:
          obj2->Fit("gaus","Q");
          obj2->GetFunction("gaus")->SetLineColor(kYellow);
          break;
        }
      
        TArrayI s0(3);
        s0[0] = 6;
        s0[1] = 1;
        s0[2] = 2;

        TArrayI s1(3);
        s1[0] = 3;
        s1[1] = 4;
        s1[2] = 5;

        TArrayI inner(1);
        inner[0] = 0;

        TArrayI outer(1);
        outer[0] = 1;

        // Just to get the right (maximum) binning
        TH1D *half[4];
        half[0] = obj1->ProjectionS(s0, inner, "Sector 6-1-2 Inner");
        half[1] = obj1->ProjectionS(s1, inner, "Sector 3-4-5 Inner");
        half[2] = obj1->ProjectionS(s0, outer, "Sector 6-1-2 Outer");
        half[3] = obj1->ProjectionS(s1, outer, "Sector 3-4-5 Outer");

        for (int i=0; i<4; i++)      
        {
            half[i]->SetLineColor(kRed+i);
            half[i]->SetDirectory(0);
            half[i]->SetBit(kCanDelete);
            half[i]->Draw("same");
        }

      gPad->Modified();
      gPad->Update();
      
    }
}
#endif

Bool_t readDatacards(TString& filename)
{
  ifstream ifun(filename.Data());
  if(!ifun)
    {
      cout << "File " << filename << " not found" << endl;
      return kFALSE;
    }

  TString word;
  
  while(ifun >> word)
    {
      // skip comments
      if(word[0]=='/' && word[1]=='/')
	{
	  while(ifun.get()!='\n'); // skip line
	  continue;
	}

      // number of events
      if(strcmp(word.Data(),"NEVENTS")==0)
	ifun >> nmaxevents;


      // input file directory
      if(strcmp(word.Data(),"IDIR")==0)
	{
	  if(idirname.Length())
	    cout << "readDataCards Warning: overriding input directory file name" << endl;
	  ifun >> idirname;
	}

      // pedestal runs
      if(strcmp(word.Data(),"PRUNS")==0)
	{
	  if(pediter.GetNumRuns())
	    cout << "readDataCards Warning: adding pedestal runs to the existing list" << endl;
	  ifun >> word;
	  pediter.AddRuns(word.Data(),idirname.Data());
	}

      // calibration runs
      if(strcmp(word.Data(),"CRUNS")==0)
	{
	  if(caliter.GetNumRuns())
	    cout << "readDataCards Warning: adding calibration runs to the existing list" << endl;
	  ifun >> word;
	  caliter.AddRuns(word.Data(),idirname.Data());
	}

      // data runs
      if(strcmp(word.Data(),"DRUNS")==0)
	{
	  if(datiter.GetNumRuns())
	    cout << "readDataCards Warning: adding data runs to the existing list" << endl;
	  ifun >> word;
	  datiter.AddRuns(word.Data(),idirname.Data());
	}
      
      // output file name
      if(strcmp(word.Data(),"OUTFILE")==0)
	{
	  if(outname.Length())
	    cout << "readDataCards Warning: overriding output file name" << endl;
	  ifun >> outname;
	}

      // calibration flag
      if(strcmp(word.Data(),"CALFLAG")==0)
	ifun >> calflag;

      // cleaning level
      if(strcmp(word.Data(),"CLEANLEVEL")==0)
	{
	  ifun >> lcore;
	  ifun >> ltail;
	}
    }

  pediter.Reset();
  caliter.Reset();
  datiter.Reset();
  TString pfile;

  // Dump read values
  cout << "************************************************" << endl;
  cout << "* Datacards read from file " << filename << endl;
  cout << "************************************************" << endl;
  cout << "Pedestal file (s): "  << endl;
  while(!(pfile=pediter.Next()).IsNull())
    cout << pfile << endl;
  cout << "Calibration file (s): "  << endl;
  while(!(pfile=caliter.Next()).IsNull())
    cout << pfile << endl;
  cout << "Data file (s): "  << endl;
  while(!(pfile=datiter.Next()).IsNull())
    cout << pfile << endl;
  cout << "Maximum number of events: " << nmaxevents << endl;
  cout << "Output file name: " << outname << endl;
  cout << "Calibration flag: " << calflag << endl;
  cout << "Cleaning level: ("<<lcore<<","<<ltail<<")" << endl;
  cout << "***********" << endl << endl;

  if(!pediter.GetNumEntries())
    {
      cout << "No pedestal file name specified" << endl;
      return kFALSE;
    }
  if(!caliter.GetNumEntries())
    {
      cout << "No calibration file name specified" << endl;
      return kFALSE;
    }
  if(!datiter.GetNumEntries())
    {
      cout << "No data file name specified" << endl;
      return kFALSE;
    }
  if(!outname.Length())
    {
      cout << "No output file name specified" << endl;
      return kFALSE;
    }


  return kTRUE;
}
