/* ======================================================================== *\
!
! *
! * 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): Markus Gaug, 11/2003 <mailto:markus@ifae.es>
!
!   Copyright: MAGIC Software Development, 2000-2004
!
!
\* ======================================================================== */
#include "MAGIC.h"

void calibration()
{

  //  const TString inpath = "/mnt/Data/rootdata/CrabNebula/2004_01_27/";
  const TString inpath = "/mnt/Data/rootdata/Miscellaneous/2004_03_03/";
  //  const TString inpath = "/home/rootdata/BlindPixel/";

  MRunIter pruns;
  MRunIter cruns;

  pruns.AddRun(20132,inpath);
  cruns.AddRun(20134,inpath);
  //  cruns.AddRun(16774,inpath);

  gStyle->SetOptStat(1111);
  gStyle->SetOptFit();
  
  MStatusDisplay *display = new MStatusDisplay;
  display->SetUpdateTime(3000);
  display->Resize(850,700);
  
  MBadPixelsCam badcam;
  badcam.AsciiRead("badpixels.dat");

  /************************************/
  /* FIRST LOOP: PEDESTAL COMPUTATION */
  /************************************/
  
  MJPedestal pedloop;
  pedloop.SetInput(&pruns);
  pedloop.SetDisplay(display);
  pedloop.SetBadPixels(badcam);

  if (!pedloop.Process())
    return;

  //
  // Create a empty Parameter List and an empty Task List
  // The tasklist is identified in the eventloop by its name
  //
  MParList  plist0;
  MTaskList tlist0;
  plist0.AddToList(&tlist0);
  
  //
  // Now setup the tasks and tasklist for the pedestals:
  // ---------------------------------------------------
  //
  MReadMarsFile read("Events");
  read.DisableAutoScheme();
  static_cast<MRead&>(read).AddFiles(pruns);

  MGeomApply      geomapl;
  MExtractSignal  sigcalc;
  
  //
  // Additionally to calculating the pedestals, 
  // you can fill histograms and look at them
  //
  MFillH fill("MHPedestalCam", "MExtractedSignalCam");
  
  tlist0.AddToList(&read);
  tlist0.AddToList(&geomapl);
  tlist0.AddToList(&sigcalc);
  tlist0.AddToList(&fill);
  
  MGeomCamMagic  geomcam;
  MHPedestalCam  hpedcam;
  plist0.AddToList(&geomcam);
  plist0.AddToList(&pedloop.GetPedestalCam());
  plist0.AddToList(&badcam);
  plist0.AddToList(&hpedcam);
  
  //
  // Create and setup the eventloop
  //
  MEvtLoop evtloop0;
  
  evtloop0.SetParList(&plist0);
  evtloop0.SetDisplay(display);
  
  //
  // Execute first analysis
  //
  if (!evtloop0.Eventloop())
    return;

  tlist0.PrintStatistics();

  for (Int_t aidx=0;aidx<2;aidx++)
    {
      hpedcam.GetAverageHiGainArea(aidx).DrawClone("fourierevents");
      hpedcam.GetAverageLoGainArea(aidx).DrawClone("fourierevents");
    }
  
  for (Int_t sector=1;sector<7;sector++)
    {
      hpedcam.GetAverageHiGainSector(sector).DrawClone("fourierevents");
      hpedcam.GetAverageLoGainSector(sector).DrawClone("fourierevents");
    }
  

  // 
  // Now the short version: 
  //
  /*
  MJCalibration calloop;
  calloop.SetInput(&cruns);
  calloop.SetDisplay(display);
  calloop.SetBadPixels(pedloop.GetBadPixelsCam());
  if (!calloop.Process(pedloop.GetPedestalCam()))
    return;
    #if 0
  */
  // 
  // The longer version:
  //
  
  //
  // Create a empty Parameter List and an empty Task List 
  //
  MParList  plist;
  MTaskList tlist;
  plist.AddToList(&tlist);
  plist.AddToList(&pedloop.GetPedestalCam());
  plist.AddToList(&pedloop.GetBadPixels());    

  gLog << endl;;
  gLog << "Calculate MCalibrationCam from Runs " << cruns.GetRunsAsString() << endl;
  gLog << endl;
  
  MReadMarsFile read("Events");
  read.DisableAutoScheme();
  static_cast<MRead&>(read).AddFiles(cruns);
  
  MGeomCamMagic              geomcam;
  MExtractedSignalCam        sigcam;
  MArrivalTimeCam            timecam;
  MCalibrationRelTimeCam     relcam;
  MCalibrationChargeCam      calcam;
  MCalibrationChargePINDiode pindiode;   
  MCalibrationChargeBlindPix blindpix;   
  
  MHCalibrationRelTimeCam     histtime;
  MHCalibrationChargeCam      histcharge;
  MHCalibrationChargePINDiode histpin;
  MHCalibrationChargeBlindPix histblind;
  histcharge.SetPulserFrequency(500);
  histblind.SetSinglePheCut(600);
  //
  // As long, as we don't have digital modules,
  // we have to set the color of the pulser LED by hand
  //
  blindpix.SetColor(kCT1);
  //    pindiode.SetColor(kCT1);
  //
  // Get the previously created MPedestalCam into the new Parameter List 
  //
  plist.AddToList(&geomcam);
  plist.AddToList(&sigcam);
  plist.AddToList(&timecam);
  plist.AddToList(&relcam);
  plist.AddToList(&calcam);
  plist.AddToList(&histtime);
  plist.AddToList(&histcharge);
  //    plist.AddToList(&histpin);
  plist.AddToList(&histblind);
  
  //
  // We saw that the signal jumps between slices, 
  // thus take the sliding window
  //		
  MExtractSignal2        sigcalc2;
  MExtractPINDiode       pincalc;
  MExtractBlindPixel     blindcalc;
  sigcalc2.SetRange(2,15,6,5,14,6);
  blindcalc.SetRange(12,17);

  MArrivalTimeCalc2      timecalc;
  MCalibrationChargeCalc calcalc;
  MGeomApply             geomapl;
  
  MFillH filltime( "MHCalibrationRelTimeCam"    , "MArrivalTimeCam");
  //   MFillH fillpin  ("MHCalibrationChargePINDiode", "MExtractedSignalPINDiode");
  MFillH fillblind("MHCalibrationChargeBlindPix", "MExtractedSignalBlindPixel");
  MFillH fillcam  ("MHCalibrationChargeCam"     , "MExtractedSignalCam");
  
  //
  // Skip the HiGain vs. LoGain calibration
  // 
  calcalc.SkipHiLoGainCalibration();
  
  // 
  // Apply a filter against cosmics
  // (was directly in MCalibrationCalc in earlier versions)
  //
  MFCosmics            cosmics;
  MContinue            cont(&cosmics);
  
  tlist.AddToList(&read);
  tlist.AddToList(&geomapl);
  tlist.AddToList(&sigcalc2);
  tlist.AddToList(&blindcalc);
  //    tlist.AddToList(&pincalc);
  //
  // In case, you want to skip the cosmics rejection, 
  // uncomment the next line
  //
  tlist.AddToList(&cont);
  //
  // In case, you want to skip the somewhat lengthy calculation
  // of the arrival times using a spline, uncomment the next two lines
  //
  tlist.AddToList(&timecalc);
  tlist.AddToList(&filltime);
  //    tlist.AddToList(&fillpin);
  tlist.AddToList(&fillblind);
  tlist.AddToList(&fillcam);
  //
  tlist.AddToList(&calcalc);
  //
  // Create and setup the eventloop
  //
  MEvtLoop evtloop;
  evtloop.SetParList(&plist);
  evtloop.SetDisplay(display);
  
  //
  // Execute second analysis
  //
  if (!evtloop.Eventloop())
    return;
  
  tlist.PrintStatistics();
  
  MBadPixelsCam *badpixels = (MBadPixelsCam*)plist->FindObject("MBadPixelsCam");

  //
  // print the most important results of all pixels to a file
  //
  /*
  MLog gauglog;
  gauglog.SetOutputFile(Form("%s%s",calcam.GetName(),".txt"),1);
  calcam.SetLogStream(&gauglog);
  badpixels->Print();
  calcam.SetLogStream(&gLog);
  */
  //
  // just one example how to get the plots of individual pixels
  //
  //  histblind.DrawClone("all");
  //  histcharge[400].DrawClone("all");
  //  histcharge(5).DrawClone("all");
  //  histtime[5].DrawClone("fourierevents");
  for (Int_t aidx=0;aidx<2;aidx++)
    {
      histcharge.GetAverageHiGainArea(aidx).DrawClone("all");
      histcharge.GetAverageLoGainArea(aidx).DrawClone("all");
    }
  
  for (Int_t sector=1;sector<7;sector++)
    {
      histcharge.GetAverageHiGainSector(sector).DrawClone("all");
      histcharge.GetAverageLoGainSector(sector).DrawClone("all");
    }
  

  // 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;FFactorPhe",           "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;BlindPixConv",         "Conversion Factor to photons (Blind Pixel Method)");
  MHCamera disp10 (geomcam, "Cal;BlindPixFFactor",      "Total F-Factor (Blind Pixel Method)");
  MHCamera disp11 (geomcam, "Cal;PINDiodeConv",         "Conversion Factor tp photons (PIN Diode Method)");
  MHCamera disp12 (geomcam, "Cal;PINDiodeFFactor",      "Total F-Factor (PIN Diode Method)");
  MHCamera disp13 (geomcam, "Cal;Excluded",             "Pixels previously excluded");
  MHCamera disp14 (geomcam, "Cal;Unsuited",             "Unsuited Pixels  ");
  MHCamera disp15 (geomcam, "Cal;Unreliable",           "Unreliable Pixels");
  MHCamera disp16 (geomcam, "Cal;HiGainOscillating",    "Oscillating Pixels High Gain");
  MHCamera disp17 (geomcam, "Cal;LoGainOscillating",    "Oscillating Pixels Low Gain");
  MHCamera disp18 (geomcam, "Cal;HiGainPickup",         "Number Pickup events Hi Gain");
  MHCamera disp19 (geomcam, "Cal;LoGainPickup",         "Number Pickup events Lo Gain");
  MHCamera disp20 (geomcam, "Cal;Saturation",           "Pixels with saturated Hi Gain");
  MHCamera disp21 (geomcam, "Cal;FFactorValid",         "Pixels with valid F-Factor calibration");
  MHCamera disp22 (geomcam, "Cal;BlindPixelValid",      "Pixels with valid BlindPixel calibration");
  MHCamera disp23 (geomcam, "Cal;PINdiodeFFactorValid", "Pixels with valid PINDiode calibration");
  
  MHCamera disp24 (geomcam, "Cal;Ped",         "Pedestals");
  MHCamera disp25 (geomcam, "Cal;PedRms",      "Pedestal RMS");
  
  MHCamera disp26 (geomcam, "time;Time",        "Rel. Arrival Times");
  MHCamera disp27 (geomcam, "time;SigmaTime",   "Sigma of Rel. Arrival Times");
  MHCamera disp28 (geomcam, "time;TimeProb",    "Probability of Time Fit");
  MHCamera disp29 (geomcam, "time;NotFitValid", "Pixels with not valid fit results");
  MHCamera disp30 (geomcam, "time;Oscillating", "Oscillating Pixels");
  
  MHCamera disp31 (geomcam, "Cal;AbsTimeMean",  "Abs. Arrival Times");
  MHCamera disp32 (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, 16);
  disp9.SetCamError(  calcam, 17);
  disp10.SetCamContent(calcam,18);
  disp10.SetCamError(  calcam,19);
  
  // PIN Diode Method
  disp11.SetCamContent(calcam,21);
  disp11.SetCamError(  calcam,22);
  disp12.SetCamContent(calcam,23);
  disp12.SetCamError(  calcam,24);
  
  // Pixels with defects
  disp13.SetCamContent(calcam,26);
  disp14.SetCamContent(*badpixels,1);
  disp15.SetCamContent(*badpixels,3);
  disp16.SetCamContent(*badpixels,10);
  disp17.SetCamContent(*badpixels,11);
  disp18.SetCamContent(calcam,27);
  disp19.SetCamContent(calcam,28);
  
  // Lo Gain calibration
  disp20.SetCamContent(calcam,29);
  
  // Valid flags
  disp21.SetCamContent(calcam,15);
  disp22.SetCamContent(calcam,20);
  disp23.SetCamContent(calcam,25);
  
  // Pedestals
  disp24.SetCamContent(calcam,30);
  disp24.SetCamError(  calcam,31);
  disp25.SetCamContent(calcam,32);
  disp25.SetCamError(  calcam,33);
  
  // Relative Times
  disp26.SetCamContent(histtime,0);
  disp26.SetCamError(  histtime,1);
  disp27.SetCamContent(histtime,2);
  disp27.SetCamError(  histtime,3);
  disp28.SetCamContent(histtime,4);
  disp29.SetCamContent(histtime,5);
  disp30.SetCamContent(histtime,6);
  
  // Absolute Times
  disp31.SetCamContent(calcam,34);
  disp31.SetCamError(  calcam,35);
  disp32.SetCamContent(calcam,35);
  
  disp1.SetYTitle("Mean Charge [FADC Counts]");
  disp2.SetYTitle("\\sigma_{Charge} [FADC Counts]");
  disp3.SetYTitle("P_{Sum} [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("Conversion Factor [Phot/FADC Count]");
  disp10.SetYTitle("\\sqrt{N_{Ph}}*\\sigma_{Charge}/\\mu_{Charge} [1]");
  
  disp11.SetYTitle("Conversion Factor [Phot/FADC Count]");
  disp12.SetYTitle("\\sqrt{N_{Ph}}*\\sigma_{Charge}/\\mu_{Charge} [1]");
  
  disp13.SetYTitle("[1]");
  disp14.SetYTitle("[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("Ped [FADC Counts ]");
  disp25.SetYTitle("RMS_{Ped} [FADC Counts ]");
  
  disp26.SetYTitle("Time Offset [ns]");
  disp27.SetYTitle("Timing resolution [ns]");
  disp28.SetYTitle("P_{Time} [1]");
  
  disp29.SetYTitle("[1]");
  disp30.SetYTitle("[1]");
  
  disp31.SetYTitle("Mean Abs. Time [FADC slice]");
  disp32.SetYTitle("RMS Abs. Time [FADC slices]");
  
  gStyle->SetOptStat(1111);
  gStyle->SetOptFit();

  // Charges
  TCanvas &c1 = display->AddTab("Fit.Charge");
  c1.Divide(2, 4);
  
  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,4);
  
  CamDraw(c2, disp3,calcam,1,1,4);
  
  // Reduced Sigmas
  TCanvas &c3 = display->AddTab("Red.Sigma");
  c3.Divide(2,4);
  
  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,4);
  
  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(2, 4);
  
  CamDraw(c5, disp9 ,calcam,1,2, 2);
  CamDraw(c5, disp10,calcam,2,2, 2);
  
  // PIN Diode Method
  TCanvas &c6 = display->AddTab("PINDiode");
  c6.Divide(2,4);
  
  CamDraw(c6, disp11,calcam,1,2, 2);
  CamDraw(c6, disp12,calcam,2,2, 2);
  
  // Defects
  TCanvas &c7 = display->AddTab("Defects");
  c7.Divide(4,2);
  
  CamDraw(c7, disp13,calcam,1,4, 0);
  CamDraw(c7, disp14,calcam,2,4, 0);
  CamDraw(c7, disp18,calcam,3,4, 0);
  CamDraw(c7, disp19,calcam,4,4, 0);
  
  // BadCam
  TCanvas &c8 = display->AddTab("Defects");
  c8.Divide(3,2);
  
  CamDraw(c8, disp15,badcam,1,3, 0);
  CamDraw(c8, disp16,badcam,2,3, 0);
  CamDraw(c8, disp17,badcam,3,3, 0);
  
  // Valid flags
  TCanvas &c9 = display->AddTab("Validity");
  c9.Divide(4,2);
  
  CamDraw(c9, disp20,calcam,1,4,0);
  CamDraw(c9, disp21,calcam,2,4,0);
  CamDraw(c9, disp22,calcam,3,4,0);
  CamDraw(c9, disp23,calcam,4,4,0);
  
  // Pedestals
  TCanvas &c10 = display->AddTab("Pedestals");
  c10.Divide(2,4);
  
  CamDraw(c10,disp24,calcam,1,2,1);
  CamDraw(c10,disp25,calcam,2,2,2);
  
  // Rel. Times
  TCanvas &c11 = display->AddTab("Fitted Rel. Times");
  c11.Divide(3,4);
  
  CamDraw(c11,disp26,calcam,1,3,2);
  CamDraw(c11,disp27,calcam,2,3,2);
  CamDraw(c11,disp28,calcam,3,3,4);
  
  // Time Defects
  TCanvas &c12 = display->AddTab("Time Def.");
  c12.Divide(2,2);
  
  CamDraw(c12, disp29,calcam,1,2, 0);
  CamDraw(c12, disp30,calcam,2,2, 0);
  
  // Abs. Times
  TCanvas &c13 = display->AddTab("Abs. Times");
  c13.Divide(2,4);
  
  CamDraw(c13,disp31,calcam,1,2,2);
  CamDraw(c13,disp32,calcam,2,2,2);
#endif

}


void CamDraw(TCanvas &c, MHCamera &cam, TObject &evt, Int_t i, Int_t j, Int_t fit)
{
  
  TArrayI s0(6);
  s0[0] = 1;
  s0[1] = 2;
  s0[2] = 3;
  s0[3] = 4;
  s0[4] = 5;
  s0[5] = 6;

  TArrayI s1(3);
  s1[0] = 6;
  s1[1] = 1;
  s1[2] = 2;
  
  TArrayI s2(3);
  s2[0] = 3;
  s2[1] = 4;
  s2[2] = 5;
  
  TArrayI inner(1);
  inner[0] = 0;
  
  TArrayI outer(1);
  outer[0] = 1;

  c.cd(i);
  gPad->SetBorderMode(0);
  gPad->SetTicks();
  cam.GetXaxis()->SetLabelOffset(0.005);
  cam.GetXaxis()->SetLabelSize(0.06);  
  cam.GetYaxis()->SetLabelOffset(0.005);
  cam.GetYaxis()->SetLabelSize(0.06);  
  cam.GetXaxis()->SetTitleOffset(0.85);
  cam.GetXaxis()->SetTitleSize(0.06);  
  cam.GetYaxis()->SetTitleOffset(0.7);
  cam.GetYaxis()->SetTitleSize(0.06);  
  MHCamera *obj1 = (MHCamera*)cam.DrawCopy("hist");
  obj1->SetDirectory(NULL);


  c.cd(i+j);
  //  obj1->AddNotify(&evt);
  obj1->SetPrettyPalette();
  obj1->Draw();

  if (fit != 0)
    {
      c.cd(i+2*j);
      gPad->SetBorderMode(0);
      gPad->SetTicks();
      TProfile *obj2 = obj1->RadialProfile(Form("%s%s",obj1->GetName(),"_rad"));
      obj2->SetDirectory(NULL);
      obj2->GetXaxis()->SetLabelOffset(0.005);
      obj2->GetXaxis()->SetLabelSize(0.06);  
      obj2->GetYaxis()->SetLabelOffset(0.005);
      obj2->GetYaxis()->SetLabelSize(0.06);  
      obj2->GetXaxis()->SetTitleOffset(0.85);
      obj2->GetXaxis()->SetTitleSize(0.06);  
      obj2->GetYaxis()->SetTitleOffset(0.7);
      obj2->GetYaxis()->SetTitleSize(0.06);  
      obj2->Draw();
      obj2->SetBit(kCanDelete);
      
      TProfile *hprof[2];
      hprof[0] = obj1->RadialProfileS(s0, inner,Form("%s%s",obj1->GetName(), "Inner"));
      hprof[1] = obj1->RadialProfileS(s0, outer,Form("%s%s",obj1->GetName(), "Outer"));

      
      for (Int_t k=0; k<2; k++)      
        {
          Double_t min = cam.GetGeomCam().GetMinRadius(k);
          Double_t max = cam.GetGeomCam().GetMaxRadius(k);

          hprof[k]->SetLineColor(kRed+k);
          hprof[k]->SetDirectory(0);
          hprof[k]->SetBit(kCanDelete);
          hprof[k]->Draw("same");
          hprof[k]->Fit("pol1","Q","",min,max);
          hprof[k]->GetFunction("pol1")->SetLineColor(kRed+k);
          hprof[k]->GetFunction("pol1")->SetLineWidth(1);
        }
      
      gPad->Modified();
      gPad->Update();

      c.cd(i+3*j);
      gPad->SetBorderMode(0);
      gPad->SetTicks();
      TH1D *obj3 = (TH1D*)obj1->Projection(Form("%s%s",obj1->GetName(),"_py"));
      obj3->SetDirectory(NULL);
//      obj3->Sumw2();
      obj3->GetXaxis()->SetLabelOffset(0.005);
      obj3->GetXaxis()->SetLabelSize(0.06);  
      obj3->GetYaxis()->SetLabelOffset(0.005);
      obj3->GetYaxis()->SetLabelSize(0.06);  
      obj3->GetXaxis()->SetTitleOffset(0.85);
      obj3->GetXaxis()->SetTitleSize(0.06);  
      obj3->GetYaxis()->SetTitleOffset(0.7);
      obj3->GetYaxis()->SetTitleSize(0.06);  
      obj3->Draw();
      obj3->SetBit(kCanDelete);

      gPad->Modified();
      gPad->Update();

      const Double_t min   = obj3->GetBinCenter(obj3->GetXaxis()->GetFirst());
      const Double_t max   = obj3->GetBinCenter(obj3->GetXaxis()->GetLast());
      const Double_t integ = obj3->Integral("width")/2.5066283;
      const Double_t mean  = obj3->GetMean();
      const Double_t rms   = obj3->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);
          obj3->Fit("sgaus","QLR");
          obj3->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.);
          obj3->Fit("dgaus","QLRM");
          obj3->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);
          obj3->Fit("tgaus","QLRM");
          obj3->GetFunction("tgaus")->SetLineColor(kYellow);
          break;
        case 4:
          obj3->Fit("pol0","Q");
          obj3->GetFunction("pol0")->SetLineColor(kYellow);
          break;
        case 9:
          break;
        default:
          obj3->Fit("gaus","Q");
          obj3->GetFunction("gaus")->SetLineColor(kYellow);
          break;
        }
      


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

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

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

