void SetupWriter(MWriteRootFile *write, const char *name, Bool_t ismc) const
{
    if (!write)
        return;

    write->SetName(name);
    write->AddContainer("MHillas",        "Events");
    write->AddContainer("MHillasSrc",     "Events");
    write->AddContainer("MHillasExt",     "Events");
    write->AddContainer("MPointingPos",   "Events");
    write->AddContainer("MHillasSrcAnti", "Events", kFALSE);
    write->AddContainer("MImagePar",      "Events", kFALSE);
    write->AddContainer("MNewImagePar",   "Events", kFALSE);
    write->AddContainer("MNewImagePar2",  "Events", kFALSE);
    write->AddContainer("Hadronness",     "Events", kFALSE);
    write->AddContainer("MSrcPosCam",     "Events", kFALSE);
    write->AddContainer("MSrcPosAnti",    "Events", kFALSE);
    write->AddContainer("ThetaSquared",   "Events", kFALSE);
    write->AddContainer("OpticalAxis",    "Events", kFALSE);
    write->AddContainer("Disp",           "Events", kFALSE);
    write->AddContainer("Ghostbuster",    "Events", kFALSE);
    write->AddContainer("MEnergyEst",     "Events", kFALSE);
    write->AddContainer("MTime",          "Events", kFALSE);
    write->AddContainer("MMcEvt",         "Events", kFALSE);
    write->AddContainer("MEnergyEst",     "Events", kFALSE);
    write->AddContainer("DataType",       "Events");
    write->AddContainer("Weight",         "Events");
    write->AddContainer("FileId",         "Events");
    write->AddContainer("EvtNumber",      "Events");

    write->AddContainer("MMcRunHeader",         "RunHeaders", kFALSE);
    write->AddContainer("MMcCorsikaRunHeader",  "RunHeaders", kFALSE);
    write->AddContainer("MCorsikaRunHeader",    "RunHeaders", kFALSE);

    //if (ismc)
    //    write->AddCopySource("OriginalMC");
}

int ganymed(const char *dataset, const char *outfile, Float_t ra, Float_t dec, Bool_t ismc=false)
{
    MDirIter files;
    if (files.ReadList(dataset)<=0)
        return 2;
    files.Sort();

    // =============================================================

    // Crab (05h34'32" 22d00'52")
    MPointingPos source("MSourcePos");
    //source.SetSkyPosition(5.5755, 22.0144); // ra[h], dec[deg]
    //source.SetSkyPosition(12.356089, 30.176897); // ra[h], dec[deg]//1218
    //source.SetSkyPosition(16.897867,39.760201);// mrk501
    //source.SetSkyPosition(11.0742, 38.2089); // ra[h], dec[deg]//mrk421

    if (ismc)
        source.SetLocalPosition(ra, dec); // RA is ZD, DEC, is AZ!!!!
    else
        source.SetSkyPosition(ra, dec); // ra[h], dec[deg]//mrk421

    UInt_t fNumOffSourcePos = 5;

    // =============================================================

    // FIXME: If fPathIn read cuts and energy estimator from file!
    MContinue contq("MImagePar.fNumIslands>3.5"
                    "|| MNewImagePar.fNumUsedPixels<5.5"
                    "|| MNewImagePar.fLeakage1>0.1",
                    "CutQ");

    MContinue cont0("MHillas.fSize<0"
                    "|| MHillasExt.fSlopeSpreadWeighted>0.68"
                    "|| MHillasExt.fSlopeSpreadWeighted<0.18"
                    "|| log10(MHillas.GetArea*MGeomCam.fConvMm2Deg^2)>(log10(MHillas.fSize)-2)*1.1-1.55"
                    "|| MNewImagePar.fConcCore < 0.13"
                    "|| MNewImagePar.fConcCOG  < 0.15"
                    , "Cut0");

    // Cuts before plots
    MContinue cont2("0"
                    "|| MHillasExt.fM3Long*sign(MHillasSrc.fCosDeltaAlpha)*MGeomCam.fConvMm2Deg<-0.07"
                    "|| log10(MNewImagePar.fConc1)<-log10(MHillas.fSize)*0.57+0.13"
                    //"|| log10(MNewImagePar.fConc1)>-log10(MHillas.fSize)*0.57+0.5"
                    //"|| MHillasSrc.fDist*MGeomCam.fConvMm2Deg>log10(MHillas.fSize)*0.42+0.29"
                    "|| MHillas.GetArea*MGeomCam.fConvMm2Deg^2<0.01"
                    "|| MHillasSrc.fDist*MGeomCam.fConvMm2Deg>log10(MHillas.fSize)*0.5+0.095"
                    "|| MHillasSrc.fDist*MGeomCam.fConvMm2Deg<0.4"
                    "|| MHillasSrc.fDist*MGeomCam.fConvMm2Deg>1.2",
                    "Cut2");

    TArrayD param(12);
    // Parametrization of Disp
    param[0]  =  1.15136;     // constant
    param[8]  =  0.0681437;   // slope
    param[9]  =  2.62932;     // leak
    param[10] =  1.51279;     // size-offset
    param[11] =  0.0507821;   // size-slope
    // Parametrization for sign of disp (m3long, slope)
    param[5]  = -0.07;
    param[6]  =  7.2;
    param[7]  =  0.5;
    // ThetaSq-Cut
    param[1]  =  0.11;  // 0.215
    // Area-Cut
    param[2]  =  0.215468;
    param[3]  =  5.63973;
    param[4]  =  0.0836169;

    // --------------------------------------------------------------------------------

    // Setup fitter and histograms
    MAlphaFitter fit;
    fit.SetBackgroundFitMin(0.09);
    fit.SetBackgroundFitMax(0.25);
    fit.SetPolynomOrder(1);
    fit.SetSignalFunction(MAlphaFitter::kThetaSq);
    fit.SetScaleUser(1./fNumOffSourcePos); // includes fit.SetScaleMode(MAlphaFitter::kUserScale);

    // --------------------------------------------------------------------------------

    MBinning bins3( 67, -0.005, 0.665, "BinningTheta",  "asin");
    MBinning binsD( 25,   0,      2.5, "BinningDist");
    MBinning binsC( 25,   0,      2.5, "BinningDistC");
    MBinning binsR( 60,   0,      300, "BinningRate");

    // --------------------------------------------------------------------------------

    MParList plist;
    plist.AddToList(&source);
    plist.AddToList(&fit);

    // Initialize default binnings
    plist.AddToList(&bins3);
    plist.AddToList(&binsD);
    plist.AddToList(&binsC);
    plist.AddToList(&binsR);

    MParameterI par("DataType");
    MParameterD weight("Weight");
    plist.AddToList(&par);
    plist.AddToList(&weight);

    // Setup Tasklist
    MTaskList tlist;
    plist.AddToList(&tlist);

    // La Palma Magic1
    MObservatory obs;
    plist.AddToList(&obs);

    // --------------------------------------------------------------------------------

    MParameterD scale;
    scale.SetVal(1./fNumOffSourcePos);

    MHThetaSq *halphaoff = plist.FindCreateObj("MHThetaSq", "HistOff");
    halphaoff->ForceUsingSize();

    MFillH falpha(halphaoff, "", "FillHist");
    MFillH fenergy1("HistEOff [MHThetaSq]", "", "FillHistE");

    MFMagicCuts cuts;
    cuts.SetHadronnessCut(MFMagicCuts::kNoCut);
    cuts.SetVariables(param);

    // Cuts before thetasq
    MContinue cont1(&cuts, "Cut1");

    // Energy estimation
    MEnergyEstimate energyest;
    energyest.SetRule("MHillas.fSize^0.815*pow(10, 0.0199*(MHillasSrc.fDist*1.17e-2 - (-6.95))^2*(MHillasSrc.fDist*1.17e-2>(-6.95)) + 0.5997 *(MHillasSrc.fDist*1.17e-2 - 0.7568)^2 * (MHillasSrc.fDist*1.17e-2>0.7568) )");

    // Cuts before writing
    MContinue cont3("",    "Cut3");
    contq.SetAllowEmpty();
    cont0.SetAllowEmpty();
    cont1.SetAllowEmpty();
    cont2.SetAllowEmpty();
    cont3.SetAllowEmpty();

    cont1.SetInverted();

    // Filter for VsSize
    MFDataPhrase ftheta("ThetaSquared.fVal < ThetaSquaredCut.fVal", "CutT");

    // ------------- Loop Off Data --------------------
    MRead *readoffdata = NULL;
    if (ismc)
    {
        MReadMarsFile *read = new MReadMarsFile("Events");
        read->DisableAutoScheme();
        readoffdata = read;
    }
    else
    {
        MReadReports *read = new MReadReports;
        read->AddTree("Events", "MTime.", MReadReports::kMaster);
        read->AddTree("Drive",MReadReports::kRequired);
        read->AddTree("Rates",MReadReports::kRequired);
        readoffdata = read;
    }

    readoffdata->AddFiles(files);

    TString fname0 = Form("%s-summary.root",  outfile);
    TString fname1 = Form("%s-analysis.root", outfile);

    MWriteRootFile dummy0(fname0.Data(), "RECREATE");
    MWriteRootFile dummy1(fname1.Data(), "RECREATE");

    MWriteRootFile *write0 = &dummy0;
    MWriteRootFile *write1 = &dummy1;
    SetupWriter(write0, "WriteAfterCut0", ismc);
    SetupWriter(write1, "WriteAfterCut3", ismc);

    MWriteRootFile write2(fname1.Data(), "RECRATE");
    write2.AddContainer("MReportRates", "Rates");
    write2.AddContainer("MTimeRates",   "Rates");

    MParameterCalc setevtnum("MRawEvtHeader.fDAQEvtNumber", "SetEvtNumber");
    setevtnum.SetNameParameter("EvtNumber");

    MParameterCalc setrunnum("MRawRunHeader.GetFileID", "SetFileId");
    setrunnum.SetNameParameter("FileId");

    MFillH fill1a("MHHillasOffPre  [MHHillas]",      "MHillas",      "FillHillasPre");
    MFillH fill2a("MHHillasOffPost [MHHillas]",      "MHillas",      "FillHillasPost");
    MFillH fill3a("MHVsSizeOffPost [MHVsSize]",      "MHillasSrc",   "FillVsSizePost");
    MFillH fill3c("MHVsSizeOffTheta [MHVsSize]",     "MHillasSrc",   "FillVsSizeTheta");
    MFillH fill4a("MHHilExtOffPost [MHHillasExt]",   "MHillasSrc",   "FillHilExtPost");
    MFillH fill5a("MHHilSrcOffPost [MHHillasSrc]",   "MHillasSrc",   "FillHilSrcPost");
    MFillH fill6a("MHImgParOffPost [MHImagePar]",    "MImagePar",    "FillImgParPost");
    MFillH fill7a("MHNewParOffPost [MHNewImagePar]", "MNewImagePar", "FillNewParPost");
    //MFillH fill9a("MHEffOffTime    [MHEffectiveOnTime]", "MTime",        "FillEffOnTime");
    fill1a.SetNameTab("PreCut");
    fill2a.SetNameTab("PostCut");
    fill3a.SetNameTab("VsSize");
    fill3c.SetNameTab("CutT");
    fill4a.SetNameTab("HilExt");
    fill5a.SetNameTab("HilSrc");
    fill6a.SetNameTab("ImgPar");
    fill7a.SetNameTab("NewPar");
    //fill9a.SetNameTab("EffOffT");

    fill3c.SetFilter(&ftheta);

    //MFDataMember fbin("Bin.fVal", '>', 0);
    //fill9a.SetFilter(&fbin);

    MTaskList tlist2;
    tlist2.SetNumPasses(fNumOffSourcePos);
    fill2a.SetWeight(&scale);
    fill3a.SetWeight(&scale);
    fill3c.SetWeight(&scale);
    fill4a.SetWeight(&scale);
    fill5a.SetWeight(&scale);
    fill6a.SetWeight(&scale);
    fill7a.SetWeight(&scale);

    // How to get source position from off- and on-data?
    MSrcPosCalc scalc;
    scalc.SetMode(MSrcPosCalc::kWobble);
    scalc.SetCallback(&tlist2);

    MSrcPosCorrect scor;

    MHillasCalc hcalc;
    MHillasCalc hcalc2("MHillasCalcAnti");
    hcalc.SetFlags(MHillasCalc::kCalcHillasSrc);
    hcalc2.SetFlags(MHillasCalc::kCalcHillasSrc);
    hcalc2.SetNameHillasSrc("MHillasSrcAnti");
    hcalc2.SetNameSrcPosCam("MSrcPosAnti");

    MSrcPosRndm srcrndm;

    tlist2.AddToList(&scalc);
    //tlist2.AddToList(&scor);
    tlist2.AddToList(&srcrndm);
    tlist2.AddToList(&hcalc);
    tlist2.AddToList(&hcalc2);
    tlist2.AddToList(&setrunnum);
    tlist2.AddToList(&setevtnum);
    tlist2.AddToList(&cuts);
    tlist2.AddToList(&energyest);
    tlist2.AddToList(write0);
    tlist2.AddToList(&cont0);  // Post/Pre shown cuts
    tlist2.AddToList(&cont1);  // Calculate ThetaSq
    tlist2.AddToList(&fill2a);
    tlist2.AddToList(&ftheta); // Calculate theta Cut
    if (!ismc)
    {
        tlist2.AddToList(&falpha);
        tlist2.AddToList(&fenergy1);
    }
    tlist2.AddToList(&fill3c); // With theta cut
    tlist2.AddToList(&cont2);  
    tlist2.AddToList(&fill3a);
    tlist2.AddToList(&fill4a);
    tlist2.AddToList(&fill5a);
    tlist2.AddToList(&fill6a);
    tlist2.AddToList(&fill7a);
    tlist2.AddToList(&cont3);
    tlist2.AddToList(write1);

    MPointingPosCalc pcalc;
    MPointingDevCalc devcalc;

    tlist.AddToList(readoffdata);
    if (!ismc)
        tlist.AddToList(&pcalc, "Drive");
    tlist.AddToList(&contq,  "Events");
    tlist.AddToList(&fill1a, "Events");
    tlist.AddToList(&tlist2, "Events");

    // by setting it here it is distributed to all consecutive tasks
    tlist.SetAccelerator(MTask::kAccDontReset|MTask::kAccDontTime);

    par.SetVal(0);
    weight.SetVal(-1./fNumOffSourcePos);

    MStatusDisplay *d = new MStatusDisplay;

    // Create and setup the eventloop
    MEvtLoop evtloop("Off");
    evtloop.SetParList(&plist);
    evtloop.SetDisplay(d);

    // Execute first analysis
    if (!evtloop.Eventloop())
        return 2;

    if (!evtloop.GetDisplay())
        return 2;

    // ------------- Loop On Data --------------------
    MRead *readondata = NULL;
    if (ismc)
    {
        MReadMarsFile *read = new MReadMarsFile("Events");
        read->DisableAutoScheme();
        readondata = read;
    }
    else
    {
        MReadReports *read = new MReadReports;
        read->AddTree("Events", "MTime.", MReadReports::kMaster);
        read->AddTree("Drive",        MReadReports::kRequired);
        read->AddTree("Rates",        MReadReports::kRequired);
        read->AddTree("Humidity",     MReadReports::kRequired);
        read->AddTree("Temperatures", MReadReports::kRequired);
        readondata = read;
    }
    readondata->AddFiles(files);

    scalc.SetMode(MSrcPosCalc::kDefault);
    scalc.SetNumRandomOffPositions(fNumOffSourcePos);

    MFillH fill1b("MHHillasOnPre  [MHHillas]",      "MHillas",      "FillHillasPre");
    MFillH fill2b("MHHillasOnPost [MHHillas]",      "MHillas",      "FillHillasPost");
    MFillH fill3b("MHVsSizeOnPost [MHVsSize]",      "MHillasSrc",   "FillVsSizePost");
    MFillH fill3d("MHVsSizeOnTheta [MHVsSize]",     "MHillasSrc",   "FillVsSizeTheta");
    MFillH fill4b("MHHilExtOnPost [MHHillasExt]",   "MHillasSrc",   "FillHilExtPost");
    MFillH fill5b("MHHilSrcOnPost [MHHillasSrc]",   "MHillasSrc",   "FillHilSrcPost");
    MFillH fill6b("MHImgParOnPost [MHImagePar]",    "MImagePar",    "FillImgParPost");
    MFillH fill7b("MHNewParOnPost [MHNewImagePar]", "MNewImagePar", "FillNewParPost");
    //MFillH fill9b("MHEffOnTime    [MHEffectiveOnTime]", "MTime",        "FillEffOnTime");
    fill1b.SetNameTab("PreCut");
    fill2b.SetNameTab("PostCut");
    fill3b.SetNameTab("VsSize");
    fill3d.SetNameTab("CutT");
    fill4b.SetNameTab("HilExt");
    fill5b.SetNameTab("HilSrc");
    fill6b.SetNameTab("ImgPar");
    fill7b.SetNameTab("NewPar");
    //fill9b.SetNameTab("EffOnT");
    fill1b.SetDrawOption("same");
    fill2b.SetDrawOption("same");
    fill3b.SetDrawOption("same");
    fill3d.SetDrawOption("same");
    fill4b.SetDrawOption("same");
    fill5b.SetDrawOption("same");
    fill6b.SetDrawOption("same");
    fill7b.SetDrawOption("same");
    //fill9b.SetFilter(&fbin);

    fill3d.SetFilter(&ftheta);

    MHThetaSq *halphaon = plist.FindCreateObj("MHThetaSq", "Hist");
    halphaon->ForceUsingSize();

    MFillH falpha2(halphaon, "", "FillHist");
    MFillH fenergy2("HistE [MHThetaSq]", "", "FillHistE");

    tlist2.SetNumPasses();

    tlist.Replace(readondata);

    MHSrcPosCam hsrcpos(kTRUE);

    MFillH fillsrc(&hsrcpos, "MSrcPosCam", "FillSrcPosCam");
    fillsrc.SetNameTab("SrcPos");

    tlist2.AddToListBefore(&fillsrc, &hcalc);

    MH3 hvs1("MPointingPos.fZd");
    MH3 hvs2("MPointingPos.fZd");
    MHVsTime hvstM("MTimeDrive.GetMjd");
    MHVsTime hvst0("MPointingPos.fZd");
    MHVsTime hvst1("MReportRates.fTriggerRate");
    MH3      hvst2("MPointingPos.fZd", "MReportRates.fTriggerRate");
    MHVsTime hvst("MReportTemperatures.GetTempMean");
    MHVsTime hvsh("MReportHumidity.GetMean");
    hvs1.SetName("ObsTime;Theta");
    hvs2.SetName("OnTime;Theta");
    hvstM.SetName("MJD");
    hvst0.SetName("ZdTime");
    hvst1.SetName("Rate");
    hvst2.SetName("RateZd;Theta;Rate");
    hvst.SetName("CamTemp");
    hvsh.SetName("CamHum");
    hvs1.SetTitle("Observation time vs. Zenith Angle;\\Theta [\\circ];T_{on} [s]");
    hvs2.SetTitle("On-Time vs. Zenith Angle;\\Theta [\\circ];T_{on} [s]");
    hvstM.SetTitle("MJD vs time;Time;MJD");
    hvst0.SetTitle("Zenith distance vs time;Time;Zd [\\circ]");
    hvst1.SetTitle("Camera trigger rate vs time;Time;Rate [Hz]");
    hvst2.SetTitle("Camera trigger rate vs Zd; Zd [\\circ];Rate [Hz]");
    hvst.SetTitle("Average sensor compartment temperature;Time;T [\\circ C]");
    hvsh.SetTitle("Average camera humidity;Time;Hum [%]");

    MFillH fillvs1(&hvs1, "", "FillObsTime");
    MFillH fillvs2(&hvs2, "", "FillOnTime");
    MFillH fillvstM(&hvstM, "MTime", "FillMjdVsTime");
    MFillH fillvst0(&hvst0, "MTimeDrive", "FillZdVsTime");
    MFillH fillvst1(&hvst1, "MTimeRates", "FillRate");
    MFillH fillvst2(&hvst2, "", "FillRateVsZd");
    MFillH fillvst(&hvst, "MTimeTemperatures", "FillCamTemp");
    MFillH fillvsh(&hvsh, "MTimeHumidity",     "FillCamHum");
    fillvs1.SetNameTab("ObsTime");
    fillvs2.SetNameTab("OnTime");
    fillvst1.SetNameTab("Rate");
    fillvst2.SetNameTab("RateZd");
    fillvst.SetNameTab("CamTemp");
    fillvsh.SetNameTab("CamHum");

    hvs1.SetWeight("MReportRates.fElapsedTime");
    hvs2.SetWeight("MReportRates.fElapsedOnTime");

    hvst1.SetMinimum(0);
    hvst2.SetDrawOption("colz");

    MH3 hrate("MRawRunHeader.GetFileID"/*, "MRawEvtHeader.GetTriggerID"*/);
    hrate.SetWeight("1./TMath::Max(MRawRunHeader.GetRunLength,1)");
    hrate.SetName("RateVsRun");
    hrate.SetTitle("Event rate after quality cuts;File Id;Event Rate [Hz];");
    hrate.InitLabels(MH3::kLabelsX);
    hrate.GetHist().SetMinimum(0);

    MFillH fillr(&hrate, "", "FillRateVsRun");
    tlist.AddToListAfter(&fillr, &contq, "Events");

    tlist.Replace(&fill1b);

    tlist2.Replace(&fill2b);
    tlist2.Replace(&fill3b);
    tlist2.Replace(&fill3d);
    tlist2.Replace(&fill4b);
    tlist2.Replace(&fill5b);
    tlist2.Replace(&fill6b);
    tlist2.Replace(&fill7b);
    tlist2.Replace(&falpha2);
    tlist2.Replace(&fenergy2);
    if (!ismc)
    {
        tlist.AddToList(&fillvs1,  "Rates");
        tlist.AddToList(&fillvs2,  "Rates");
        tlist.AddToList(&fillvst1, "Rates");
        tlist.AddToList(&fillvst2, "Rates");
        tlist.AddToList(&write2,   "Rates");
        tlist.AddToList(&fillvst,  "Temperatures");
        tlist.AddToList(&fillvsh,  "Humidity");
        tlist.AddToList(&fillvst0, "Drive");
        tlist.AddToList(&fillvstM, "Drive");
    }

    // by setting it here it is distributed to all consecutive tasks
    tlist.SetAccelerator(MTask::kAccDontReset|MTask::kAccDontTime);

    par.SetVal(1);
    weight.SetVal(1);

    // Execute first analysis
    if (!evtloop.Eventloop())
        return 2;

    if (!evtloop.GetDisplay())
        return 2;

    delete readondata;
    delete readoffdata;

    TObjArray cont;
    //cont.Add(const_cast<TEnv*>(GetEnv()));
    //cont.Add(const_cast<MSequence*>(&fSequence));

    //TNamed cmdline("CommandLine", fCommandLine.Data());
    //cont.Add(&cmdline);

    if (d)
    {
        TString title = "--  Analysis #";
        title += dataset;
        title += "  --";
        d->SetTitle(title, kFALSE);

        cont.Add(d);
        TString path;
        path += fname1;

        d->Write();
    }

//    if (!WriteContainer(cont, fname1, "RECREATE"))
//        return;

    return 0;
}
