/* ======================================================================== *\
!
! *
! * 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-2003
!
!
\* ======================================================================== */

const TString pedfile = "../20040215_16770_P_OffCrab4_E.root";
const TString calfile = "../20040215_16771_C_OffCrab4_E.root";

//const TString pedfile = "/mnt/users/mdoro/Mars/Data/20040201_14418_P_OffMrk421-1_E.root";
//const TString calfile = "/mnt/users/mdoro/Mars/Data/20040201_1441*_C_OffMrk421-1_E.root";

//const TString pedfile = "/mnt/Data/rootdata/CrabNebula/2004_02_10/20040210_14607_P_CrabNebula_E.root";
//const TString calfile = "/mnt/Data/rootdata/CrabNebula/2004_02_10/20040210_14608_C_CrabNebula_E.root";

//const TString pedfile = "/mnt/Data/rootdata/CrabNebula/2004_01_26/20040125_10412_P_Crab-On_E.root";
//const TString calfile = "/mnt/Data/rootdata/CrabNebula/2004_01_26/20040125_1041*_C_Crab-On_E.root";

//const TString pedfile = "/mnt/Data/rootdata/Miscellaneous/2003_12_19/20031218_03522_P_Park_E.root";
//const TString calfile = "/mnt/Data/rootdata/Miscellaneous/2003_12_19/20031218_03527_C_Park_E.root";

void calibration(TString pedname=pedfile,
             TString calname=calfile)
{

    gStyle->SetOptStat(1111);
    gStyle->SetOptFit();

    MStatusDisplay *display = new MStatusDisplay;
    display->SetUpdateTime(500);
    display->Resize(850,700);

    //
    // Create a empty Parameter List and an empty Task List
    // The tasklist is identified in the eventloop by its name
    //
    MParList  plist;

    MTaskList tlist;
    plist.AddToList(&tlist);

    //
    // Now setup the tasks and tasklist for the pedestals:
    // ---------------------------------------------------
    //

    MReadMarsFile read("Events", pedname);
    read.DisableAutoScheme();

    MGeomApply      geomapl;
    MExtractSignal  sigcalc;

    //
    // Set the extraction range higher:
    //		
    //sigcalc.SetRange(1,14,1,14);

    MPedCalcPedRun pedcalc;

    //
    // Additionally to calculating the pedestals, 
    // you can fill histograms and look at them
    //
    MFillH fill("MHPedestalCam", "MExtractedSignalCam");

    tlist.AddToList(&read);
    tlist.AddToList(&geomapl);
    tlist.AddToList(&sigcalc);
    tlist.AddToList(&pedcalc);
    tlist.AddToList(&fill);

    MGeomCamMagic  geomcam;
    MPedestalCam   pedcam;
    MHPedestalCam  hpedcam;
    plist.AddToList(&geomcam);
    plist.AddToList(&pedcam);
    plist.AddToList(&hpedcam);

    //
    // Create and setup the eventloop
    //
    MEvtLoop evtloop;
    evtloop.SetParList(&plist);
    evtloop.SetDisplay(display);
 
    //
    // Execute first analysis
    //
    if (!evtloop.Eventloop())
        return;

    tlist.PrintStatistics();

    // 
    // Look at one specific pixel, after all the histogram manipulations:
    //
    hpedcam[9].DrawClone("fourierevents");


    MHCamera dispped0  (geomcam, "Ped;Pedestal",               "Mean per Slice");
    MHCamera dispped1  (geomcam, "Ped;PedestalErr",            "Mean Error per Slice");
    MHCamera dispped2  (geomcam, "Ped;PedestalRms",            "RMS per Slice");
    MHCamera dispped3  (geomcam, "Ped;PedestalRmsErr",         "RMS Error per Slice");

    MHCamera dispped4  (geomcam, "Ped;Mean",                   "Fitted Mean per Slice");
    MHCamera dispped5  (geomcam, "Ped;MeanErr",                "Fitted Error of Mean per Slice");
    MHCamera dispped6  (geomcam, "Ped;Sigma",                  "Fitted Sigma per Slice");
    MHCamera dispped7  (geomcam, "Ped;SigmaErr",               "Fitted Error of Sigma per Slice");
    MHCamera dispped8  (geomcam, "Ped;Prob",                   "Probability of Fit");
    MHCamera dispped9  (geomcam, "Ped;DeltaPedestalMean",      "Rel. Diff. Mean per Slice (Calc.-Fitte)");
    MHCamera dispped10 (geomcam, "Ped;DeltaPedestalMeanError", "Rel. Diff. Mean Error per Slice (Calc.-Fitted)");
    MHCamera dispped11  (geomcam, "Ped;DeltaRmsSigma",         "Rel. Diff. RMS per Slice (Calc.-Fitted)");
    MHCamera dispped12  (geomcam, "Ped;DeltaRmsSigmaError",    "Rel. Diff. RMS Error per Slice (Calc.-Fitted)");
    MHCamera dispped13  (geomcam, "Ped;FitOK",                 "Gaus Fit not OK");
    MHCamera dispped14  (geomcam, "Ped;FourierOK",             "Fourier Analysis not OK");

    dispped0.SetCamContent(  pedcam, 0);
    dispped0.SetCamError(    pedcam, 1);
    dispped1.SetCamContent(  pedcam, 1);
    dispped2.SetCamContent(  pedcam, 2);
    dispped2.SetCamError(    pedcam, 3);
    dispped3.SetCamContent(  pedcam, 3);

    dispped4.SetCamContent( hpedcam, 0);
    dispped4.SetCamError(   hpedcam, 1);
    dispped5.SetCamContent( hpedcam, 1);
    dispped6.SetCamContent( hpedcam, 2);
    dispped6.SetCamError(   hpedcam, 3);
    dispped7.SetCamContent( hpedcam, 3);
    dispped8.SetCamContent( hpedcam, 4);
    dispped9.SetCamContent( hpedcam, 5);
    dispped9.SetCamError(   hpedcam, 6);
    dispped10.SetCamContent(hpedcam, 7);
    dispped11.SetCamContent(hpedcam, 8);
    dispped11.SetCamError(  hpedcam, 9);
    dispped12.SetCamContent(hpedcam, 10);
    dispped13.SetCamContent(hpedcam, 11);
    dispped14.SetCamContent(hpedcam, 12);

    dispped0.SetYTitle("Calc. Pedestal per slice [FADC counts]");
    dispped1.SetYTitle("Calc. Pedestal Error per slice [FADC counts]");
    dispped2.SetYTitle("Calc. Pedestal RMS per slice [FADC counts]");
    dispped3.SetYTitle("Calc. Pedestal RMS Error per slice [FADC counts]");
    dispped4.SetYTitle("Fitted Mean per slice [FADC counts]");
    dispped5.SetYTitle("Error of Fitted Mean per slice [FADC counts]");
    dispped6.SetYTitle("Fitted Sigma per slice [FADC counts]");
    dispped7.SetYTitle("Error of Fitted Sigma per slice [FADC counts]");
    dispped8.SetYTitle("Fit Probability [1]");
    dispped9.SetYTitle("Rel. Diff. Pedestal Calc.-Fitted per slice [1]");
    dispped10.SetYTitle("Rel. Diff. Pedestal Error Calc.-Fitted per slice [1]");
    dispped11.SetYTitle("Rel. Diff. Pedestal RMS Calc.-Fitted per slice [1]");
    dispped12.SetYTitle("Rel. Diff. Pedestal RMS Error Calc.-Fitted per slice [1]");
    dispped13.SetYTitle("[1]");
    dispped14.SetYTitle("[1]");
    
    // Histogram values
    TCanvas &b1 = display->AddTab("Ped.Calc."); 
    b1.Divide(4,3);

    CamDraw(b1,dispped0,pedcam,1,4,1);
    CamDraw(b1,dispped1,pedcam,2,4,2);
    CamDraw(b1,dispped2,pedcam,3,4,2);
    CamDraw(b1,dispped3,pedcam,4,4,2);

    // Fitted values 
    TCanvas &b2 = display->AddTab("Ped.Fit"); 
    b2.Divide(4,3);

    CamDraw(b2,dispped4,hpedcam,1,4,1);
    CamDraw(b2,dispped5,hpedcam,2,4,2);
    CamDraw(b2,dispped6,hpedcam,3,4,2);
    CamDraw(b2,dispped7,hpedcam,4,4,2);


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

    CamDraw(b3,dispped8,hpedcam,1,1,3);

    // Differences
    TCanvas &c4 = display->AddTab("Rel.Diff.Calc.-Fit");
    c4.Divide(4,3);

    CamDraw(c4,dispped9,hpedcam,1,4,1);
    CamDraw(c4,dispped10,hpedcam,2,4,1);
    CamDraw(c4,dispped11,hpedcam,3,4,1);
    CamDraw(c4,dispped12,hpedcam,4,4,1);

    // Defects
    TCanvas &c5 = display->AddTab("Defects");
    c5.Divide(2,2);

    CamDraw(c5,dispped13,hpedcam,1,2,0);
    CamDraw(c5,dispped14,hpedcam,2,2,0);

    //
    // Create a empty Parameter List and an empty Task List
    //
    MParList  plist2;
    MTaskList tlist2;
    plist2.AddToList(&tlist2);

    MExtractedSignalCam     sigcam;
    MCalibrationChargeCam   calcam;
    MHCalibrationRelTimeCam timecam;
    //
    // Get the previously created MPedestalCam into the new Parameter List 
    //
    plist2.AddToList(&geomcam);
    plist2.AddToList(&pedcam);
    plist2.AddToList(&sigcam);
    plist2.AddToList(&calcam);
    plist2.AddToList(&timecam);
    //
    // Get the MAGIC geometry
    //
    tlist2.AddToList(&geomapl);
    //
    // Now setup the new tasks and tasklist for the calibration
    // ---------------------------------------------------
    //

    MReadMarsFile read2("Events", calname);
    read2.DisableAutoScheme();

    //
    // We saw that the signal jumps between slices, 
    // thus take the sliding window
    //		
    MExtractPINDiode       pincalc;
    MExtractBlindPixel     blindcalc;
    MExtractSignal2        sigcalc2;
    MArrivalTimeCalc2      timecalc;
    MCalibrationChargeCalc calcalc;

    MFillH filltime("MHCalibrationRelTimeCam", "MArrivalTime");
    MFillH fillpin("MHCalibrationChargePINDiode", "MExtractedSignalPINDiode");

    //
    // Set the range (other than default) 
    // of FADC slices for the blind pixel
    //
    // calcalc.SetBlindPixelRange(10,25);

    //
    // Set the cut upon which a superposition of the blind pixel 
    // FADC slices will be filled into the SinglePHE histogram
    //
    //    calcalc.SetBlindPixelSinglePheCut(500);

    //
    // Skip the HiGain vs. LoGain calibration
    // 
    calcalc.SkipHiLoGainCalibration();

    //
    // As long, as we don't have digital modules, 
    // we have to set the color of the pulser LED by hand
    //
//    calcalc.SetPulserColor(MCalibrationCalc::kECT1);

    //
    // In case, we want to exclude a pre-defined list of bad pixels:
    // (This is a preliminary feature)
    //
    // calcalc.ExcludePixelsFromAsciiFile("badpixels.dat");
    
    //
    // In case, you want to skip the quality checks
    // (NOT RECOMMENDED!!!)
    //
    // calcalc.SkipQualityChecks();

    //
    // In case, we want to apply another fit function to the 
    // blind pixel 
    //
    //    MCalibrationBlindPix *bp = calcam.GetBlindPixel();
    //    bp->ChangeFitFunc(MHCalibrationBlindPixel::kEPoisson5);

    // 
    // Apply a filter against cosmics
    // (was directly in MCalibrationCalc in earlier versions)
    //
    MFCosmics            cosmics;
    MContinue            cont(&cosmics);

    tlist2.AddToList(&read2);
    tlist2.AddToList(&blindcalc);
    tlist2.AddToList(&pincalc);
    tlist2.AddToList(&sigcalc2);
    //
    // In case, you want to skip the cosmics rejection, 
    // uncomment the next line
    //
    tlist2.AddToList(&cont);
    //
    // In case, you want to skip the somewhat lengthy calculation
    // of the arrival times using a spline, uncomment the next two lines
    //
    tlist2.AddToList(&timecalc);
    tlist2.AddToList(&filltime);
    tlist2.AddToList(&fillpin);
    //
    tlist2.AddToList(&calcalc);
    //
    // Create and setup the eventloop
    //
    MEvtLoop evtloop2;
    evtloop2.SetParList(&plist2);
    evtloop2.SetDisplay(display);
    
    //
    // Execute second analysis
    //
    if (!evtloop2.Eventloop())
        return;

    tlist2.PrintStatistics();

    //
    // print the most important results of all pixels
    //
    calcam.Print();

    //
    // just one example how to get the plots of individual pixels
    //
    calcam[564].DrawClone();

    // 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 Phe's (F-Factor Method)");
    MHCamera disp7  (geomcam, "Cal;FFactorConv",    "Conversion Factor (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 (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 (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;Oscillating",    "Oscillating Pixels");
    MHCamera disp19 (geomcam, "Cal;Saturation",     "Pixels with saturated Hi Gain");

    MHCamera disp20 (geomcam, "Cal;Ped",         "Pedestals");
    MHCamera disp21 (geomcam, "Cal;PedRms",      "Pedestal RMS");

    MHCamera disp22 (geomcam, "time;Time",        "Rel. Arrival Times");
    MHCamera disp23 (geomcam, "time;SigmaTime",   "Sigma of Rel. Arrival Times");
    MHCamera disp24 (geomcam, "time;TimeProb",    "Probability of Time Fit");
    MHCamera disp25 (geomcam, "time;NotFitValid", "Pixels with not valid fit results");
    MHCamera disp26 (geomcam, "time;Oscillating", "Oscillating Pixels");

    MHCamera disp27 (geomcam, "Cal;AbsTimeMean", "Abs. Arrival Times");
    MHCamera disp28 (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(calcam,29);
    disp18.SetCamContent(calcam,30);

    // Lo Gain calibration
    disp19.SetCamContent(calcam,31);


    // Pedestals
    disp20.SetCamContent(calcam,35);
    disp20.SetCamError(  calcam,36);
    disp21.SetCamContent(calcam,37);
    disp21.SetCamError(  calcam,38);


    // Relative Times
    disp22.SetCamContent(timecam,0);
    disp22.SetCamError(  timecam,1);
    disp23.SetCamContent(timecam,2);
    disp23.SetCamError(  timecam,3);
    disp24.SetCamContent(timecam,4);
    disp25.SetCamContent(timecam,5);
    disp26.SetCamContent(timecam,6);

    // Absolute Times
    disp27.SetCamContent(calcam,39);
    disp27.SetCamError(  calcam,40);
    disp28.SetCamContent(calcam,41);
    disp28.SetCamError(  calcam,42);


    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 [PhE/FADC Count]");
    disp8.SetYTitle("\\sqrt{N_{PhE}}*\\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("Ped [FADC Counts ]");
    disp21.SetYTitle("RMS_{Ped} [FADC Counts ]");

    disp22.SetYTitle("Time Offset [ns]");
    disp23.SetYTitle("Timing resolution [ns]");
    disp24.SetYTitle("P_{Time} [1]");

    disp25.SetYTitle("[1]");
    disp26.SetYTitle("[1]");

    disp27.SetYTitle("Mean Abs. Time [FADC slice]");
    disp28.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, disp17,calcam,3,4, 0);
    CamDraw(c7, disp18,calcam,4,4,0);

    // Lo Gain Calibration
    TCanvas &c8 = display->AddTab("LowGain");
    c8.Divide(1,3);

    CamDraw(c8, disp19,calcam,1,4,0);


    // Pedestals
    TCanvas &c9 = display->AddTab("Pedestals");
    c9.Divide(2,3);

    CamDraw(c9,disp20,calcam,1,2,1);
    CamDraw(c9,disp21,calcam,2,2,2);

    // Rel. Times
    TCanvas &c10 = display->AddTab("Fitted Rel. Times");
    c10.Divide(3,3);

    CamDraw(c10,disp22,calcam,1,3,2);
    CamDraw(c10,disp23,calcam,2,3,2);
    CamDraw(c10,disp24,calcam,3,3,4);

    // Time Defects
    TCanvas &c11 = display->AddTab("Time Def.");
    c11.Divide(2,2);

    CamDraw(c11, disp25,calcam,1,2, 0);
    CamDraw(c11, disp26,calcam,2,2, 0);

    // Abs. Times
    TCanvas &c12 = display->AddTab("Abs. Times");
    c12.Divide(2,3);

    CamDraw(c12,disp27,calcam,1,2,2);
    CamDraw(c12,disp28,calcam,2,2,2);
}

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();
      
      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;
        }
      
      gPad->Modified();
      gPad->Update();
      
    }
}
