Index: branches/Mars_McMismatchStudy/fact/analysis/callisto.C
===================================================================
--- branches/Mars_McMismatchStudy/fact/analysis/callisto.C	(revision 17960)
+++ branches/Mars_McMismatchStudy/fact/analysis/callisto.C	(revision 17962)
@@ -1,5 +1,5 @@
 #include "MLogManip.h"
 
-int callisto(const char *seqfile="seq/2012/01/23/20120123_023.seq", const char *outpath = "output")
+int callisto(const char *seqfile="seq/2012/01/23/20120123_023.seq", const char *outpath = "output", bool use_delays=false)
 {
     // ======================================================
@@ -12,182 +12,8 @@
 
     // map file to use (get that from La Palma!)
-    const char *map = usemap ? "/home/fact/FACT++/FACTmap111030.txt" : NULL;
-/*
-    Bool_t maximum = kTRUE;
-
-    const char *lp_template    = maximum ?
-        "template-lp-extractor-maximum.root" :
-        "template-lp-extractor-leading-edge.root";
-
-    const char *pulse_template = "template-pulse.root";
-*/
+//    const char *map = usemap ? "/home/fact/FACT++/FACTmap111030.txt" : NULL;
+    const char *map = usemap ? "/scratch/fact/FACTmap111030.txt" : NULL;
+
     // ------------------------------------------------------
-
-    //bool use_delays=false;
-
-    int spike_removal=3;
-
-    // The gain as extracted from our dark count spectra
-    double gain = 258;
-
-    // ------------------------------------------------------
-
-    // Extraction range in slices. It will always(!) contain the full range
-    // of integration
-    const int first_slice =  25; //  10ns
-    const int last_slice  = 225; // 125ns
-
-    // Note that rise and fall time mean different things whether you use IntegralFixed or IntegraRel:
-    //
-    //  IntegralFixed:
-    //    * fRiseTime: Number of slices left  from arrival time
-    //    * fFallTime: Number of slices right from arrival time
-    //  IntegralRel:
-    //    * fRiseTime: Number of slices left  from maximum time
-    //    * fFallTime: Number of slices right from maximum time
-    //
-/*
-    const int rise_time_cal = maximum ?  40 :  10; // was 13;   5ns
-    const int fall_time_cal = maximum ? 120 : 160; // was 23;  80ns
-
-    const int rise_time_dat = maximum ?  10 :   2; // was 13; was 10;   1ns
-    const int fall_time_dat = maximum ?  40 :  48; // was 23; was 40;  24ns
-
-    // Extraction type: Extract integral and half leading edge
-
-    const int type = maximum ? (MExtralgoSpline::kAmplitudeRel) : (MExtralgoSpline::kIntegralFixed);
-    //const int type = MExtralgoSpline::kIntegralFixed;
-
-
-    const double heighttm   = 0.5; // IntegralAbs { 1.5pe * 9.6mV/pe } / IntegralRel { 0.5 }
-*/
-    Long_t max  =    0;  // All
-    Long_t max0 =  max;  // Time marker
-    //Long_t max1 =  max;  // Light pulser
-    //Long_t max2 = 3000;  // Calibration ratio
-    //Long_t max3 =  max;  // Pedestal Rndm
-    //Long_t max4 =  max;  // Pedestal Ext
-    Long_t max5 =  max;  // Data
-
-    // ======================================================
-
-    if (map && gSystem->AccessPathName(map, kFileExists))
-    {
-        gLog << err << "ERROR - Cannot access mapping file '" << map << "'" << endl;
-        return 1;
-    }
-
-    // The sequence file which defines the files for the analysis
-    MSequence seq(seqfile);
-    if (!seq.IsValid())
-    {
-        gLog << err << "ERROR - Sequence '" << seqfile << "' invalid!" << endl;
-        return 2;
-    }
-
-    // --------------------------------------------------------------------------------
-
-    gLog.Separator("Callisto");
-    gLog << all << "Calibrate data of sequence '" << seq.GetFileName() << "'" << endl;
-    gLog << endl;
-
-    // ------------------------------------------------------
-
-    ostringstream drsname;
-    drsname << gSystem->DirName(seqfile) << "/";
-    drsname << seq.GetNight().GetNightAsInt() << "_";
-    drsname << Form("%03d", seq.GetDrsSequence()) << ".drs.seq";
-
-    MSequence drs(drsname.str().c_str());
-    if (!drs.IsValid())
-    {
-        gLog << err << "ERROR - DRS sequence invalid!" << endl;
-        return 3;
-    }
-
-    gLog << all << "DRS sequence file: " << drsname.str() << '\n' << endl;
-
-    TString drsfile = seq.GetFileName(0, MSequence::kRawDrs);
-    if (drsfile.IsNull())
-    {
-        gLog << err << "No DRS file available in sequence." << endl;
-        return 4;
-    }
-
-    TString timfile = drs.GetFileName(0, MSequence::kFitsDat);
-    TString drs1024 = drs.GetFileName(0, MSequence::kFitsDrs);
-    TString calfile = seq.GetFileName(0, MSequence::kFitsCal);
-    //TString pedfile = seq.GetFileName(0, MSequence::kFitsPed);
-
-    gLog << all;
-    gLog << "DRS calib     300: " << drsfile << '\n';
-    gLog << "DRS calib    1024: " << drs1024 << "\n\n";
-
-    MDrsCalibration drscalib300;
-    if (!drscalib300.ReadFits(drsfile.Data()))
-        return 5;
-
-    MDrsCalibration drscalib1024;
-    if (!drscalib1024.ReadFits(drs1024.Data()))
-        return 6;
-
-    gLog << all;
-    gLog << "Time calibration : " << timfile << '\n';
-    gLog << "Light Pulser file: " << calfile << '\n' << endl;
-    //gLog << "Pedestal     file: " << pedfile << '\n';
-
-    // ------------------------------------------------------
-
-    MDirIter iter;
-    if (seq.GetRuns(iter, MSequence::kFitsDat)<=0)
-    {
-        gLog << err << "ERROR - Sequence valid but without files." << endl;
-        return 7;
-    }
-    iter.Print();
-
-    // ======================================================
-
-    MStatusArray arrt, arrp;
-
-    TFile ft(lp_template);
-    if (arrt.Read()<=0)
-    {
-        gLog << err << "ERROR - Reading LP template from " << lp_template << endl;
-        return 100;
-    }
-
-    MHCamera *lpref = (MHCamera*)arrt.FindObjectInCanvas("ExtCalSig;avg", "MHCamera", "Cam");
-    if (!lpref)
-    {
-        gLog << err << "ERROR - LP Template not found in " << lp_template << endl;
-        return 101;
-    }
-    lpref->SetDirectory(0);
-
-    MHCamera *gain = (MHCamera*)arrt.FindObjectInCanvas("gain", "MHCamera", "Gain");
-    if (!gain)
-    {
-        gLog << err << "ERROR - Gain not found in " << lp_template << endl;
-        return 101;
-    }
-    gain->SetDirectory(0);
-
-    TFile fp(pulse_template);
-    if (arrp.Read()<=0)
-    {
-        gLog << err << "ERROR - Reading Pulse template from " << pulse_template << endl;
-        return 102;
-    }
-
-    TH1F *hpulse = (TH1F*)arrp.FindObjectInCanvas("hPixelEdgeMean0_0", "TH1F", "cgpPixelPulses0");
-    if (!hpulse)
-    {
-        gLog << err << "ERROR - Pulse Template not found in " << pulse_template << endl;
-        return 103;
-    }
-    hpulse->SetDirectory(0);
-
-    // ======================================================
 
     MStatusDisplay *d = new MStatusDisplay;
@@ -201,5 +27,4 @@
     badpixels[1208].SetUnsuitable(MBadPixelsPix::kUnsuitable);
     badpixels[1399].SetUnsuitable(MBadPixelsPix::kUnsuitable);
-
     //  Twin pixel
     //     113
@@ -210,5 +35,152 @@
     //    1393
 
-    MDrsCalibrationTime timecam;
+    // ------------------------------------------------------
+
+    // ------------------------------------------------------
+
+    // Calib: 51 / 90 / 197 (20% TH)
+    // Data:  52 / 64 / 104 (20% TH)
+
+    // Extraction range in slices. It will always(!) contain the full range
+    // of integration
+    const int first_slice =  25; //  10ns
+    const int last_slice  = 225; // 125ns
+
+    const double heighttm   = 0.5; // IntegralAbs { 1.5pe * 9.6mV/pe } / IntegralRel { 0.5 }
+
+    Long_t max  =    0;  // All
+    Long_t max0 =  max;  // Time marker
+    Long_t max1 =  max;  // Light pulser
+    //Long_t max2 = 3000;  // Calibration ratio
+    Long_t max3 =  max;  // Pedestal Rndm
+    Long_t max4 =  max;  // Pedestal Ext
+    Long_t max5 =  max;  // Data
+
+    // ========================= Result ==================================
+
+    //double scale = 0.1;
+    double scale = 0.1024;
+
+    // ======================================================
+
+    if (map && gSystem->AccessPathName(map, kFileExists))
+    {
+        gLog << err << "ERROR - Cannot access mapping file '" << map << "'" << endl;
+        return 1;
+    }
+
+    // The sequence file which defines the files for the analysis
+    MSequence seq(seqfile);
+    if (!seq.IsValid())
+    {
+        gLog << err << "ERROR - Sequence '" << seqfile << "' invalid!" << endl;
+        return 2;
+    }
+
+    // --------------------------------------------------------------------------------
+
+    gLog.Separator("Callisto");
+    gLog << all << "Calibrate data of sequence '" << seq.GetFileName() << "'" << endl;
+    gLog << endl;
+
+    // ------------------------------------------------------
+
+    ostringstream drsname;
+    drsname << gSystem->DirName(seqfile) << "/";
+    drsname << seq.GetNight().GetNightAsInt() << "_";
+    drsname << Form("%03d", seq.GetDrsSequence()) << ".drs.seq";
+
+    MSequence drs(drsname.str().c_str());
+    if (!drs.IsValid())
+    {
+        gLog << err << "ERROR - DRS sequence invalid!" << endl;
+        return 3;
+    }
+
+    gLog << all << "DRS sequence file: " << drsname.str() << '\n' << endl;
+
+    TString drsfile = seq.GetFileName(0, MSequence::kRawDrs);
+    if (drsfile.IsNull())
+    {
+        gLog << err << "No DRS file available in sequence." << endl;
+        return 4;
+    }
+
+    TString timfile = drs.GetFileName(0, MSequence::kFitsDat);
+    TString drs1024 = drs.GetFileName(0, MSequence::kFitsDrs);
+    TString pedfile = seq.GetFileName(0, MSequence::kFitsPed);
+    TString calfile = seq.GetFileName(0, MSequence::kFitsCal);
+
+    gLog << all;
+    gLog << "DRS calib     300: " << drsfile << '\n';
+    gLog << "DRS calib    1024: " << drs1024 << "\n\n";
+
+    MDrsCalibration drscalib300;
+    if (!drscalib300.ReadFits(drsfile.Data()))
+        return 5;
+
+    MDrsCalibration drscalib1024;
+    if (!drscalib1024.ReadFits(drs1024.Data()))
+        return 6;
+
+    gLog << all;
+    gLog << "Time calibration : " << timfile << '\n';
+    gLog << "Pedestal     file: " << pedfile << '\n';
+    gLog << "Light Pulser file: " << calfile << '\n' << endl;
+
+    // ------------------------------------------------------
+
+    MDirIter iter;
+    if (seq.GetRuns(iter, MSequence::kFitsDat)<=0)
+    {
+        gLog << err << "ERROR - Sequence valid but without files." << endl;
+        return 7;
+    }
+    iter.Print();
+
+    // ======================================================
+
+/*
+    MStatusArray arrt, arrp;
+
+    TFile ft(lp_template);
+    if (arrt.Read()<=0)
+    {
+        gLog << err << "ERROR - Reading LP template from " << lp_template << endl;
+        return 100;
+    }
+
+    MHCamera *lpref = (MHCamera*)arrt.FindObjectInCanvas("ExtCalSig;avg", "MHCamera", "Cam");
+    if (!lpref)
+    {
+        gLog << err << "ERROR - LP Template not found in " << lp_template << endl;
+        return 101;
+    }
+    lpref->SetDirectory(0);
+
+    MHCamera *gain = (MHCamera*)arrt.FindObjectInCanvas("gain", "MHCamera", "Gain");
+    if (!gain)
+    {
+        gLog << err << "ERROR - Gain not found in " << lp_template << endl;
+        return 101;
+    }
+    gain->SetDirectory(0);
+
+    TFile fp(pulse_template);
+    if (arrp.Read()<=0)
+    {
+        gLog << err << "ERROR - Reading Pulse template from " << pulse_template << endl;
+        return 102;
+    }
+
+    TH1F *hpulse = (TH1F*)arrp.FindObjectInCanvas("hPixelEdgeMean0_0", "TH1F", "cgpPixelPulses0");
+    if (!hpulse)
+    {
+        gLog << err << "ERROR - Pulse Template not found in " << pulse_template << endl;
+        return 103;
+    }
+    hpulse->SetDirectory(0);
+*/
+    // ======================================================
 
     // Plot the trigger pattern rates vs. run-number
@@ -223,9 +195,5 @@
     // hrate.DefaultLabelY("ERROR");
 
-    Bool_t isinteg =
-        (type&MExtralgoSpline::kIntegral)    ||
-        (type&MExtralgoSpline::kFixedWidth)  ||
-        (type&MExtralgoSpline::kDynWidth)
-        ? kTRUE : kFALSE;
+    MDrsCalibrationTime timecam;
 
     gStyle->SetOptFit(kTRUE);
@@ -241,5 +209,4 @@
     plist0.AddToList(&tlist0);
     plist0.AddToList(&drscalib1024);
-    plist0.AddToList(&badpixels);
     plist0.AddToList(&timecam);
 
@@ -257,5 +224,5 @@
 
     MDrsCalibApply drsapply0;
-    drsapply0.SetRemoveSpikes(spike_removal);
+    //drsapply0.SetRemoveSpikes(4);
 
     MFillH fill0("MHDrsCalibrationTime");
@@ -281,5 +248,5 @@
 
     // ======================================================
-/*
+
     gLog << endl;
     gLog.Separator("Processing external light pulser run");
@@ -308,35 +275,13 @@
 
     MDrsCalibApply drsapply1;
-    drsapply1.SetRemoveSpikes(spike_removal);
+    //drsapply1.SetRemoveSpikes(4);
+
+    MFilterData filterdata1;
 
     // ---
 
-    MExtractTimeAndChargeSpline extractor1b("ExtractPulse");
+    MExtractFACT extractor1b("ExtractPulse");
     extractor1b.SetRange(first_slice, last_slice);
-    extractor1b.SetRiseTimeHiGain(rise_time_cal);
-    extractor1b.SetFallTimeHiGain(fall_time_cal);
-    extractor1b.SetHeightTm(heighttm);
-    extractor1b.SetChargeType(type);
-    extractor1b.SetSaturationLimit(600000);
     extractor1b.SetNoiseCalculation(kFALSE);
-
-    MExtractTimeAndChargeSpline extractor1c("ExtractAmplitude");
-    extractor1c.SetRange(first_slice, last_slice);
-    extractor1c.SetChargeType(MExtralgoSpline::kAmplitude);
-    extractor1c.SetSaturationLimit(600000);
-    extractor1c.SetNoiseCalculation(kFALSE);
-    extractor1c.SetNameSignalCam("Amplitude");
-    extractor1c.SetNameTimeCam("AmplitudePos");
-
-    // ---
-
-    MHCamEvent evt1a(5, "CalRatio", "Ratio per slice between integrated signal and amplitude;; r [1/n]");
-    evt1a.SetNameSub("Amplitude", kTRUE);
-    MFillH fill1a(&evt1a, "MExtractedSignalCam", "FillRatio");
-    fill1a.SetDrawOption("gaus");
-
-    MParameterD ratio1a;
-    ratio1a.SetVal(1./(fall_time_cal+rise_time_cal));
-    fill1a.SetWeight(&ratio1a);
 
     // ---
@@ -391,10 +336,6 @@
     tlist1.AddToList(&drsapply1);
     tlist1.AddToList(&cont1);
+    tlist1.AddToList(&filterdata1);
     tlist1.AddToList(&extractor1b);
-    if (isinteg)
-    {
-        tlist1.AddToList(&extractor1c);
-        tlist1.AddToList(&fill1a);
-    }
     tlist1.AddToList(&calctm1a);
     tlist1.AddToList(&treat1);
@@ -416,123 +357,6 @@
         timecam.SetDelays(*evt1h.GetHist());
 
-    // ========================= Result ==================================
-
-    Double_t avgS = evt1f.GetHist()->GetMean();
-    Double_t medS = evt1f.GetHist()->GetMedian();
-    Double_t rmsS = evt1f.GetHist()->GetRMS();
-    Double_t maxS = evt1f.GetHist()->GetMaximum();
-
-    MArrayF der1(hpulse->GetNbinsX());
-    MArrayF der2(hpulse->GetNbinsX());
-
-    MExtralgoSpline spline(hpulse->GetArray()+1, hpulse->GetNbinsX(),
-                           der1.GetArray(), der2.GetArray());
-    spline.SetRiseFallTime(rise_time_dat, fall_time_dat);
-    spline.SetExtractionType(type);
-    spline.SetHeightTm(heighttm);
-
-    spline.Extract(hpulse->GetMaximumBin()-1);
-
-    // The pulser signal is most probably around 400mV/9.5mV
-    // IntegraFixed 2/48 corresponds to roughly 215mV*50slices
-    Double_t scale = 1./spline.GetSignal();
-
-    MArrayD calib(1440);
-    for (int i=0; i<1440; i++)
-    {
-        Double_t g = gain->GetBinContent(i+1)>0.5 ? gain->GetBinContent(i+1) : 1;
-        if (evt1f.GetHist()->GetBinContent(i+1)>0 && !badpixels[i].IsUnsuitable())
-            calib[i] = lpref->GetBinContent(i+1) / evt1f.GetHist()->GetBinContent(i+1) / g;
-    }
-
-    gROOT->SetSelectedPad(0);
-    d->AddTab("PulseTemp");
-    gPad->SetGrid();
-    hpulse->SetNameTitle("Pulse", "Single p.e. pulse template");
-    hpulse->SetDirectory(0);
-    hpulse->SetLineColor(kBlack);
-    hpulse->DrawCopy();
-
-    TAxis *ax = hpulse->GetXaxis();
-
-    Double_t w = hpulse->GetBinWidth(1);
-    Double_t T = w*(spline.GetTime()+0.5)       +ax->GetXmin();
-    Double_t H = w*(hpulse->GetMaximumBin()+0.5)+ax->GetXmin();
-
-    TLine line;
-    line.SetLineColor(kRed);
-    line.DrawLine(T-rise_time_dat*w, spline.GetHeight(),
-                  T+fall_time_dat*w, spline.GetHeight());
-    line.DrawLine(T, spline.GetHeight()/4, T, 3*spline.GetHeight()/4);
-    line.DrawLine(T-rise_time_dat*w, 0,
-                  T-rise_time_dat*w, spline.GetHeight());
-    line.DrawLine(T+fall_time_dat*w, 0,
-                  T+fall_time_dat*w, spline.GetHeight());
-
-    TGraph gg;
-    for (int ix=1; ix<=hpulse->GetNbinsX(); ix++)
-        for (int i=0; i<10; i++)
-        {
-            Double_t x = hpulse->GetBinLowEdge(ix)+i*hpulse->GetBinWidth(ix)/10.;
-            gg.SetPoint(gg.GetN(), x+w/2, spline.EvalAt(ix-1+i/10.));
-        }
-
-    gg.SetLineColor(kBlue);
-    gg.SetMarkerColor(kBlue);
-    gg.SetMarkerStyle(kFullDotMedium);
-    gg.DrawClone("L");
-
-    gROOT->SetSelectedPad(0);
-    d->AddTab("CalConst");
-    MGeomCamFACT fact;
-    MHCamera hcalco(fact);
-    hcalco.SetName("CalConst");
-    hcalco.SetTitle(Form("Relative calibration constant [%.0f/pe]", 1./scale));
-    hcalco.SetCamContent(calib);
-    hcalco.SetAllUsed();
-    //hcalco.Scale(scale);
-    hcalco.DrawCopy();
-*/
-
-    TF1 f("f", "[0]*(1-1/(1+exp(x/[1])))*exp(-x/[2])", -15, 150);
-    f.SetParameter(0, gain*0.05143);
-    f.SetParameter(1, 1.075);
-    f.SetParameter(2, 19.3);
-    f.SetNpx(2*165); //2GHz
-
-    // Convert to graph
-    TGraph g(&f);
-
-    // Convert to float
-    MArrayF x(g.GetN());
-    MArrayF y(g.GetN());
-    for (int i=0; i<g.GetN(); i++)
-    {
-        x[i] = g.GetX()[i];
-        y[i] = g.GetY()[i];
-    }
-
-    MFilterData filter;
-    filter.Filter(1, g.GetN(), y.GetArray(), y.GetArray());
-
-    // Define spline
-    MArrayF der1(g.GetN());
-    MArrayF der2(g.GetN());
-
-    MExtralgoSpline spline(y.GetArray(), y.GetSize(),
-                           der1.GetArray(), der2.GetArray());
-
-    spline.SetExtractionType(MExtralgoSpline::kAmplitudeRel);
-    spline.SetHeightTm(0.5);
-
-    // Estimate where the maximum is and extract signal
-    Long64_t maxi = TMath::LocMax(y.GetSize(), y.GetArray());
-    spline.Extract(maxi);
-
-    // Scale factor for signal extraction
-    double scale = 1./spline.GetSignal();
-
     // ======================================================
-/*
+
     gLog << endl;
     gLog.Separator("Extracting random pedestal");
@@ -563,20 +387,17 @@
 
     MDrsCalibApply drsapply3;
-    drsapply3.SetRemoveSpikes(spike_removal);
+    //drsapply3.SetRemoveSpikes(4);
+
+    MFilterData filterdata3;
 
     //---
 
-    MExtractTimeAndChargeSpline extractor3;
+    MExtractFACT extractor3;
     extractor3.SetRange(first_slice, last_slice);
-    extractor3.SetRiseTimeHiGain(rise_time_dat);
-    extractor3.SetFallTimeHiGain(fall_time_dat);
-    extractor3.SetHeightTm(heighttm);
-    extractor3.SetChargeType(type);
-    extractor3.SetSaturationLimit(600000);
     extractor3.SetNoiseCalculation(kTRUE);
 
     MCalibrateFact conv3;
     conv3.SetScale(scale);
-    conv3.SetCalibConst(calib);
+    //conv3.SetCalibConst(calib);
 
     MBadPixelsTreat treat3;
@@ -597,4 +418,5 @@
     tlist3.AddToList(&drsapply3);
     tlist3.AddToList(&cont3);
+    tlist3.AddToList(&filterdata3);
     tlist3.AddToList(&extractor3);
 //    tlist3.AddToList(&fill3a);
@@ -637,18 +459,15 @@
 
     MDrsCalibApply drsapply4;
-    drsapply4.SetRemoveSpikes(spike_removal);
-
-    MExtractTimeAndChargeSpline extractor4;
+    //drsapply4.SetRemoveSpikes(4);
+
+    MFilterData filterdata4;
+
+    MExtractFACT extractor4;
     extractor4.SetRange(first_slice, last_slice);
-    extractor4.SetRiseTimeHiGain(rise_time_dat);
-    extractor4.SetFallTimeHiGain(fall_time_dat);
-    extractor4.SetHeightTm(heighttm);
-    extractor4.SetChargeType(type);
-    extractor4.SetSaturationLimit(600000);
     extractor4.SetNoiseCalculation(kFALSE);
 
     MCalibrateFact conv4;
     conv4.SetScale(scale);
-    conv4.SetCalibConst(calib);
+    //conv4.SetCalibConst(calib);
 
     MBadPixelsTreat treat4;
@@ -668,6 +487,6 @@
     tlist4.AddToList(&drsapply4);
     tlist4.AddToList(&cont4);
+    tlist4.AddToList(&filterdata4);
     tlist4.AddToList(&extractor4);
-//    tlist4.AddToList(&fill4a);
     tlist4.AddToList(&conv4);
     tlist4.AddToList(&treat4);
@@ -679,5 +498,5 @@
     if (!loop4.GetDisplay())
         return 16;
-*/
+
     // ===================================================================
 
@@ -708,5 +527,9 @@
 
     MDrsCalibApply drsapply5;
-    drsapply5.SetRemoveSpikes(spike_removal);
+    //drsapply5.SetRemoveSpikes(4);
+
+    MTreatSaturation treatsat5;
+
+    MFilterData filterdata5;
 
     MFDataPhrase filterdat("(MRawEvtHeader.GetTriggerID&0xff00)==0",     "SelectDat");
@@ -719,32 +542,20 @@
     // ---
 
-    MTreatSaturation treatsat5;
-    MFilterData filterdata5;
-
     MExtractFACT extractor5dat;
     extractor5dat.SetRange(first_slice, last_slice);
-/*
-    MExtractTimeAndChargeSpline extractor5cal;
+    extractor5dat.SetNoiseCalculation(kFALSE);
+
+    MExtractFACT extractor5cal;
     extractor5cal.SetRange(first_slice, last_slice);
-    extractor5cal.SetRiseTimeHiGain(rise_time_cal);
-    extractor5cal.SetFallTimeHiGain(fall_time_cal);
-    extractor5cal.SetHeightTm(heighttm);
-    extractor5cal.SetChargeType(type);
-    extractor5cal.SetSaturationLimit(600000);
     extractor5cal.SetNoiseCalculation(kFALSE);
 
-    MExtractTimeAndChargeSpline extractor5tm("ExtractTM");
+    MExtractFACT extractor5tm("ExtractTM");
     extractor5tm.SetRange(last_slice, 294);
-    extractor5tm.SetRiseTimeHiGain(1);
-    extractor5tm.SetFallTimeHiGain(1);
-    extractor5tm.SetHeightTm(0.5);
-    extractor5tm.SetChargeType(MExtralgoSpline::kAmplitudeRel);
-    extractor5tm.SetSaturationLimit(600000);
     extractor5tm.SetNoiseCalculation(kFALSE);
     extractor5tm.SetNameSignalCam("TimeMarkerAmplitude");
     extractor5tm.SetNameTimeCam("TimeMarkerTime");
-*/
+
     extractor5dat.SetFilter(&filterncl);
-    //extractor5cal.SetFilter(&filtercal);
+    extractor5cal.SetFilter(&filtercal);
     //extractor4tm.SetFilter(&filtercal);
 
@@ -752,9 +563,9 @@
     MCalibrateFact conv5;
     conv5.SetScale(scale);
-    conv5.SetCalibConst(calib);
+    //conv5.SetCalibConst(calib);
 
     MCalibrateDrsTimes calctm5;
     calctm5.SetNameUncalibrated("UncalibratedTimes");
-/*
+
     MCalibrateDrsTimes calctm5tm("CalibrateTimeMarker");
     calctm5tm.SetNameArrivalTime("TimeMarkerTime");
@@ -806,10 +617,10 @@
     //fill5w.SetDrawOption("gaus");
     //fill5x.SetDrawOption("gaus");
-*/
+
 
     MBadPixelsTreat treat5;
     treat5.SetProcessPedestalRun(kFALSE);
     treat5.SetProcessPedestalEvt(kFALSE);
-/*
+
     MHSectorVsTime hist5cal("CalVsTm");
     MHSectorVsTime hist5ped("PedVsTm");
@@ -829,5 +640,5 @@
     fill5cal.SetFilter(&filtercal);
     fill5ped.SetFilter(&filterped);
-*/
+
     MHCamEvent evt5b(0, "ExtSig",   "Extracted signal;;S [mV·sl]");
     MHCamEvent evt5c(0, "CalSig",   "Calibrated and interpolated signal;;S [~phe]");
@@ -860,18 +671,10 @@
     // The second rule is for the case reading raw-files!
     MWriteRootFile write5(2, fname, "RECREATE", "Calibrated Data");
-    write5.AddContainer("MRawRunHeader",       "RunHeaders");
-    write5.AddContainer("MGeomCam",            "RunHeaders");
-    write5.AddContainer("MMcCorsikaRunHeader", "RunHeaders", kFALSE);
-    write5.AddContainer("MCorsikaRunHeader",   "RunHeaders", kFALSE);
-    write5.AddContainer("MMcRunHeader",        "RunHeaders", kFALSE);
-
-    // Common events
-    write5.AddContainer("MCorsikaEvtHeader",   "Events", kFALSE);
-    write5.AddContainer("MMcEvt",              "Events", kFALSE);
-    write5.AddContainer("IncidentAngle",       "Events", kFALSE);
-    write5.AddContainer("MPointingPos",        "Events", kFALSE);
-    write5.AddContainer("MSignalCam",          "Events");
-    write5.AddContainer("MTime",               "Events", kFALSE);
-    write5.AddContainer("MRawEvtHeader",       "Events");
+    write5.AddContainer("MRawRunHeader", "RunHeaders");
+    write5.AddContainer("MGeomCam",      "RunHeaders");
+    write5.AddContainer("MSignalCam",    "Events");
+    write5.AddContainer("MTime",         "Events");
+    write5.AddContainer("MRawEvtHeader", "Events");
+    //write.AddContainer("MTriggerPattern", "Events");
 
     // ------------------ Setup histograms and fill tasks ----------------
@@ -879,5 +682,5 @@
     MContinue test;
     test.SetFilter(&filterncl);
-/*
+
     MTaskList tlist5tm;
     tlist5tm.AddToList(&extractor5tm);
@@ -894,5 +697,5 @@
     tlist5tm.AddToList(&fill5x);
     tlist5tm.SetFilter(&filtercal);
-*/
+
     MTaskList tlist5dat;
     tlist5dat.AddToList(&fill5b);
@@ -914,15 +717,15 @@
     tlist5.AddToList(&filterdata5);
     tlist5.AddToList(&extractor5dat);
-    //tlist5.AddToList(&extractor5cal);
+    tlist5.AddToList(&extractor5cal);
     tlist5.AddToList(&calctm5);
-    //tlist5.AddToList(&tlist5tm);
+    tlist5.AddToList(&tlist5tm);
     tlist5.AddToList(&conv5);
     tlist5.AddToList(&treat5);
-    //tlist5.AddToList(&fill5ped);
-    //tlist5.AddToList(&fill5cal);
+    tlist5.AddToList(&fill5ped);
+    tlist5.AddToList(&fill5cal);
     tlist5.AddToList(&tlist5dat);
     tlist5.AddToList(&write5);
 
-    if (!loop5.Eventloop(max4))
+    if (!loop5.Eventloop(max5))
         return 18;
 
@@ -945,13 +748,2 @@
     return 0;
 }
-
-int callisto(const ULong64_t seqnum, const char *outpath = "output")
-{
-    UInt_t night = seqnum/1000;
-    UInt_t num   = seqnum%1000;
-
-    TString file = Form("/scratch/fact/sequences/%04d/%02d/%02d/%06d_%03d.seq",
-                        night/10000, (night/100)%100, night%100, num);
-
-    return callisto(file.Data(), outpath);
-}
Index: branches/Mars_McMismatchStudy/fact/analysis/mc/callisto.C
===================================================================
--- branches/Mars_McMismatchStudy/fact/analysis/mc/callisto.C	(revision 17960)
+++ branches/Mars_McMismatchStudy/fact/analysis/mc/callisto.C	(revision 17962)
@@ -32,5 +32,5 @@
 #include "MHSectorVsTime.h"
 #include "MHCamEvent.h"
-#include "MExtractTimeAndChargeSpline.h"
+#include "MExtractFACT.h"
 #include "MFillH.h"
 #include "MDrsCalibApply.h"
@@ -75,8 +75,8 @@
  */
 
-int callisto(const TString drsfile="test300samples.drs.fits",
+int callisto(const TString drsfile="test300samples2.drs.fits.gz",
              const TString pedfile="00000001.001_P_MonteCarlo000_Events.fits",
              const TString datfile="00000003.387_D_MonteCarlo010_Events.fits",
-             TString outfile = "",
+             TString outpath = "",
              TString displayfile = "", TString displaytitle = "")
 {
@@ -88,9 +88,11 @@
 
     FileStat_t fstat;
-    int rc = gSystem->GetPathInfo(outfile, fstat);
+    int rc = gSystem->GetPathInfo(outpath, fstat);
     bool isdir = !rc || R_ISDIR(fstat.fMode);
 
-    const char *buf = gSystem->ConcatFileName(outfile, "callisto.root");
-    outfile = buf;
+    TString filename = datfile + "_callisto.root";
+    filename.Replace(0, filename.Last('/')+1, "");
+    const char *buf = gSystem->ConcatFileName(outpath, filename);
+    TString outfile = buf;
     delete [] buf;
 
@@ -119,116 +121,7 @@
 
     // map file to use (get that from La Palma!)
-    const char *pmap = usemap ? "TestForThomas/FACT/FACTmap111030.txt" : NULL;
-
-    Bool_t maximum = kTRUE;
-
-    //const char *lp_template    = maximum ?
-    //    "/cm/shared/apps/fact/Mars_svn_LP/template-lp-extractor-maximum.root" :
-    //    "/cm/shared/apps/fact/Mars_svn_LP/template-lp-extractor-leading-edge.root";
-
-    const char *pulse_template = "TestForThomas/FACT/template-pulse.root";
+    const char *pmap = usemap ? "/home/isdc/toscanos/FACT/Mars_svn/resources/FACTmap111030.txt" : NULL;
 
     // ------------------------------------------------------
-
-    // Calib: 51 / 90 / 197 (20% TH)
-    // Data:  52 / 64 / 104 (20% TH)
-
-    // Extraction range in slices. It will always(!) contain the full range
-    // of integration
-    const int first_slice =  20; //  10ns
-    const int last_slice  = 250; // 125ns
-
-    // Note that rise and fall time mean different things whether you use IntegralFixed or IntegraRel:
-    //
-    //  IntegralFixed:
-    //    * fRiseTime: Number of slices left  from arrival time
-    //    * fFallTime: Number of slices right from arrival time
-    //  IntegralRel:
-    //    * fRiseTime: Number of slices left  from maximum time
-    //    * fFallTime: Number of slices right from maximum time
-    //
-    const int rise_time_cal = maximum ?  40 :  10; // was 13;   5ns
-    const int fall_time_cal = maximum ? 120 : 160; // was 23;  80ns
-
-    const int rise_time_dat = maximum ?  10 :   2; // was 13; was 10;   1ns
-    const int fall_time_dat = maximum ?  40 :  48; // was 23; was 40;  24ns
-
-    // Extraction type: Extract integral and half leading edge
-
-    const MExtralgoSpline::ExtractionType_t type = maximum ? (MExtralgoSpline::kIntegralRel) : (MExtralgoSpline::kIntegralFixed);
-    //const int type = MExtralgoSpline::kIntegralFixed;
-
-
-    const double heighttm   = 0.5; // IntegralAbs { 1.5pe * 9.6mV/pe } / IntegralRel { 0.5 }
-
-    Long_t max  =    0;  // All
-    Long_t max3 =  max;  // Pedestal Rndm
-    Long_t max4 =  max;  // Pedestal Ext
-
-    // ======================================================
-
-    if (pmap && gSystem->AccessPathName(pmap, kFileExists))
-    {
-        gLog << err << "ERROR - Cannot access mapping file '" << pmap << "'" << endl;
-        return 1;
-    }
-
-    gLog.Separator("Callisto");
-    gLog << all;
-    gLog << "Data File:     " << datfile << '\n';
-    gLog << "DRS calib 300: " << drsfile << endl;;
-
-    MDrsCalibration drscalib300;
-    if (!drscalib300.ReadFits(drsfile.Data())) {
-        gLog << err << "ERROR - Cannot access drscallib300 file '" << drsfile << "'" << endl;
-        return 5;
-    }
-    gLog << all;
-    gLog << "Pedestal file: " << pedfile << '\n';
-    gLog << "Output file:   " << outfile << '\n';
-    gLog << "Display file:  " << displayfile << '\n';
-    gLog << "Display title: " << displaytitle << endl;
-
-    // ------------------------------------------------------
-    MStatusArray arrt, arrp;
-
-    // TFile ft(lp_template);
-    // if (arrt.Read()<=0)
-    // {
-    //     gLog << err << "ERROR - Reading LP template from " << lp_template << endl;
-    //     return 100;
-    // }
-
-    // MHCamera *lpref = (MHCamera*)arrt.FindObjectInCanvas("ExtCalSig;avg", "MHCamera", "Cam");
-    // if (!lpref)
-    // {
-    //     gLog << err << "ERROR - LP Template not found in " << lp_template << endl;
-    //     return 101;
-    // }
-    // lpref->SetDirectory(0);
-
-    // MHCamera *gain = (MHCamera*)arrt.FindObjectInCanvas("gain", "MHCamera", "Gain");
-    // if (!gain)
-    // {
-    //     gLog << err << "ERROR - Gain not found in " << lp_template << endl;
-    //     return 101;
-    // }
-    // gain->SetDirectory(0);
-
-    TFile fp(pulse_template);
-    if (arrp.Read()<=0)
-    {
-        gLog << err << "ERROR - Reading Pulse template from " << pulse_template << endl;
-        return 102;
-    }
-
-    TH1F *hpulse = (TH1F*)arrp.FindObjectInCanvas("hPixelEdgeMean0_0", "TH1F", "cgpPixelPulses0");
-    if (!hpulse)
-    {
-        gLog << err << "ERROR - Pulse Template not found in " << pulse_template << endl;
-        return 103;
-    }
-    hpulse->SetDirectory(0);
-    // ======================================================
 
     MStatusDisplay *d = new MStatusDisplay;
@@ -236,5 +129,4 @@
     MBadPixelsCam badpixels;
     badpixels.InitSize(1440);
-    /*
     badpixels[ 424].SetUnsuitable(MBadPixelsPix::kUnsuitable);
     badpixels[ 583].SetUnsuitable(MBadPixelsPix::kUnsuitable);
@@ -243,5 +135,4 @@
     badpixels[1208].SetUnsuitable(MBadPixelsPix::kUnsuitable);
     badpixels[1399].SetUnsuitable(MBadPixelsPix::kUnsuitable);
-    */
     //  Twin pixel
     //     113
@@ -251,4 +142,91 @@
     //    1195
     //    1393
+
+    // ------------------------------------------------------
+
+    // Calib: 51 / 90 / 197 (20% TH)
+    // Data:  52 / 64 / 104 (20% TH)
+
+    // Extraction range in slices. It will always(!) contain the full range
+    // of integration
+    const int first_slice =  25; //  10ns
+    const int last_slice  = 225; // 125ns
+
+    Long_t max  =    0;  // All
+    Long_t max3 =  max;  // Pedestal Rndm
+    Long_t max4 =  max;  // Pedestal Ext
+
+    // ======================================================
+
+    //double scale = 0.1;
+    double scale = 0.1024;
+
+    // ======================================================
+
+    if (pmap && gSystem->AccessPathName(pmap, kFileExists))
+    {
+        gLog << err << "ERROR - Cannot access mapping file '" << pmap << "'" << endl;
+        return 1;
+    }
+
+    gLog.Separator("Callisto");
+    gLog << all;
+    gLog << "Data File:     " << datfile << '\n';
+    gLog << "DRS calib 300: " << drsfile << endl;;
+
+    MDrsCalibration drscalib300;
+    if (!drscalib300.ReadFits(drsfile.Data())) {
+        gLog << err << "ERROR - Cannot access drscallib300 file '" << drsfile << "'" << endl;
+        return 5;
+    }
+    gLog << all;
+    gLog << "Pedestal file: " << pedfile << '\n';
+    gLog << "Output file:   " << outfile << '\n';
+    gLog << "Display file:  " << displayfile << '\n';
+    gLog << "Display title: " << displaytitle << endl;
+
+    // ------------------------------------------------------
+/*
+    MStatusArray arrt, arrp;
+
+    TFile ft(lp_template);
+    if (arrt.Read()<=0)
+    {
+        gLog << err << "ERROR - Reading LP template from " << lp_template << endl;
+        return 100;
+    }
+
+    MHCamera *lpref = (MHCamera*)arrt.FindObjectInCanvas("ExtCalSig;avg", "MHCamera", "Cam");
+    if (!lpref)
+    {
+        gLog << err << "ERROR - LP Template not found in " << lp_template << endl;
+        return 101;
+    }
+    lpref->SetDirectory(0);
+
+    MHCamera *gain = (MHCamera*)arrt.FindObjectInCanvas("gain", "MHCamera", "Gain");
+    if (!gain)
+    {
+        gLog << err << "ERROR - Gain not found in " << lp_template << endl;
+        return 101;
+    }
+    gain->SetDirectory(0);
+
+    TFile fp(pulse_template);
+    if (arrp.Read()<=0)
+    {
+        gLog << err << "ERROR - Reading Pulse template from " << pulse_template << endl;
+        return 102;
+    }
+
+    TH1F *hpulse = (TH1F*)arrp.FindObjectInCanvas("hPixelEdgeMean0_0", "TH1F", "cgpPixelPulses0");
+    if (!hpulse)
+    {
+        gLog << err << "ERROR - Pulse Template not found in " << pulse_template << endl;
+        return 103;
+    }
+    hpulse->SetDirectory(0);
+    */
+    // ======================================================
 
     MDrsCalibrationTime timecam;
@@ -268,5 +246,5 @@
 
     // ========================= Result ==================================
-
+/*
     //~ Double_t avgS = evt1f.GetHist()->GetMean();
     //~ Double_t medS = evt1f.GetHist()->GetMedian();
@@ -340,5 +318,5 @@
     //hcalco.Scale(scale);
     hcalco.DrawCopy();
-
+*/
     // ======================================================
 
@@ -372,34 +350,15 @@
     MDrsCalibApply drsapply3;
 
+    MFilterData filterdata3;
+
     //---
 
-    MExtractTimeAndChargeSpline extractor3;
+    MExtractFACT extractor3;
     extractor3.SetRange(first_slice, last_slice);
-    extractor3.SetRiseTimeHiGain(rise_time_dat);
-    extractor3.SetFallTimeHiGain(fall_time_dat);
-    extractor3.SetHeightTm(heighttm);
-    extractor3.SetChargeType(type);
-    extractor3.SetSaturationLimit(600000);
     extractor3.SetNoiseCalculation(kTRUE);
-
-//    MHCamEvent evt2a(0, "PedRdm", "Extracted Pedestal Signal;;S");
-
-//    MFillH fill2a(&evt2a, "MExtractedSignalCam", "FillPedRndm");
-
-    // Use this for data, but not for calibration events
-//    evt2a.SetErrorSpread(kFALSE);
-
-    /*
-     MCalibrateData conv3;
-     conv3.SetCalibrationMode(MCalibrateData::kNone);
-     conv3.SetPedestalFlag(MCalibrateData::kNo);
-     conv3.SetCalibConvMinLimit(0);
-     conv3.SetCalibConvMaxLimit(10000);
-     conv3.SetScaleFactor(scale);
-     */
 
     MCalibrateFact conv3;
     conv3.SetScale(scale);
-    conv3.SetCalibConst(calib);
+    //conv3.SetCalibConst(calib);
 
     MBadPixelsTreat treat3;
@@ -420,4 +379,5 @@
     tlist3.AddToList(&drsapply3);
     tlist3.AddToList(&cont3);
+    tlist3.AddToList(&filterdata3);
     tlist3.AddToList(&extractor3);
 //    tlist3.AddToList(&fill3a);
@@ -461,30 +421,13 @@
     MDrsCalibApply drsapply4;
 
-    MExtractTimeAndChargeSpline extractor4;
+    MFilterData filterdata4;
+
+    MExtractFACT extractor4;
     extractor4.SetRange(first_slice, last_slice);
-    extractor4.SetRiseTimeHiGain(rise_time_dat);
-    extractor4.SetFallTimeHiGain(fall_time_dat);
-    extractor4.SetHeightTm(heighttm);
-    extractor4.SetChargeType(type);
-    extractor4.SetSaturationLimit(600000);
     extractor4.SetNoiseCalculation(kFALSE);
 
-    //    MHCamEvent evt3a(0, "PedExt", "Extracted Pedestal Signal;;S");
-
-    //    MFillH fill3a(&evt3a, "MExtractedSignalCam", "FillPedExt");
-
-    // Use this for data, but not for calibration events
-//    evt3a.SetErrorSpread(kFALSE);
-/*
-    MCalibrateData conv4;
-    conv4.SetCalibrationMode(MCalibrateData::kNone);
-    conv4.SetPedestalFlag(MCalibrateData::kNo);
-    conv4.SetCalibConvMinLimit(0);
-    conv4.SetCalibConvMaxLimit(10000);
-    conv4.SetScaleFactor(scale);
-*/
     MCalibrateFact conv4;
     conv4.SetScale(scale);
-    conv4.SetCalibConst(calib);
+    //conv4.SetCalibConst(calib);
 
     MBadPixelsTreat treat4;
@@ -504,6 +447,6 @@
     tlist4.AddToList(&drsapply4);
     tlist4.AddToList(&cont4);
+    tlist4.AddToList(&filterdata4);
     tlist4.AddToList(&extractor4);
-//    tlist4.AddToList(&fill4a);
     tlist4.AddToList(&conv4);
     tlist4.AddToList(&treat4);
@@ -550,4 +493,8 @@
     MDrsCalibApply drsapply5;
 
+    MTreatSaturation treatsat5;
+
+    MFilterData filterdata5;
+
     MFDataPhrase filterdat("(MRawEvtHeader.GetTriggerID&0xff00)==0",     "SelectDat");
     MFDataPhrase filtercal("(MRawEvtHeader.GetTriggerID&0xff00)==0x100", "SelectCal");
@@ -559,29 +506,14 @@
     // ---
 
-    MExtractTimeAndChargeSpline extractor5dat;
+    MExtractFACT extractor5dat;
     extractor5dat.SetRange(first_slice, last_slice);
-    extractor5dat.SetRiseTimeHiGain(rise_time_dat);
-    extractor5dat.SetFallTimeHiGain(fall_time_dat);
-    extractor5dat.SetHeightTm(heighttm);
-    extractor5dat.SetChargeType(type);
-    extractor5dat.SetSaturationLimit(600000);
     extractor5dat.SetNoiseCalculation(kFALSE);
 
-    MExtractTimeAndChargeSpline extractor5cal;
+    MExtractFACT extractor5cal;
     extractor5cal.SetRange(first_slice, last_slice);
-    extractor5cal.SetRiseTimeHiGain(rise_time_cal);
-    extractor5cal.SetFallTimeHiGain(fall_time_cal);
-    extractor5cal.SetHeightTm(heighttm);
-    extractor5cal.SetChargeType(type);
-    extractor5cal.SetSaturationLimit(600000);
     extractor5cal.SetNoiseCalculation(kFALSE);
 
-    MExtractTimeAndChargeSpline extractor5tm("ExtractTM");
+    MExtractFACT extractor5tm("ExtractTM");
     extractor5tm.SetRange(last_slice, 294);
-    extractor5tm.SetRiseTimeHiGain(1);
-    extractor5tm.SetFallTimeHiGain(1);
-    extractor5tm.SetHeightTm(0.5);
-    extractor5tm.SetChargeType(MExtralgoSpline::kAmplitudeRel);
-    extractor5tm.SetSaturationLimit(600000);
     extractor5tm.SetNoiseCalculation(kFALSE);
     extractor5tm.SetNameSignalCam("TimeMarkerAmplitude");
@@ -593,15 +525,7 @@
 
     // ---
-/*
-    MCalibrateData conv5;
-    conv5.SetCalibrationMode(MCalibrateData::kNone);
-    conv5.SetPedestalFlag(MCalibrateData::kNo);
-    conv5.SetCalibConvMinLimit(0);
-    conv5.SetCalibConvMaxLimit(10000);
-    conv5.SetScaleFactor(scale);
-*/
     MCalibrateFact conv5;
     conv5.SetScale(scale);
-    conv5.SetCalibConst(calib);
+    //conv5.SetCalibConst(calib);
 
     MCalibrateDrsTimes calctm5;
@@ -689,4 +613,6 @@
     tlist5.AddToList(&filterped);
     tlist5.AddToList(&fill5a);
+    tlist5.AddToList(&treatsat5);
+    tlist5.AddToList(&filterdata5);
     tlist5.AddToList(&extractor5dat);
     tlist5.AddToList(&extractor5cal);
Index: branches/Mars_McMismatchStudy/fact/analysis/mc/callisto.C.GOOD
===================================================================
--- branches/Mars_McMismatchStudy/fact/analysis/mc/callisto.C.GOOD	(revision 17962)
+++ branches/Mars_McMismatchStudy/fact/analysis/mc/callisto.C.GOOD	(revision 17962)
@@ -0,0 +1,721 @@
+#include <sstream>
+#include <iostream>
+
+#include "MLog.h"
+#include "MLogManip.h"
+
+#if !defined(__CINT__) || defined(__MAKECINT__)
+
+#include "TH1F.h"
+#include "TFile.h"
+#include "TStyle.h"
+#include "TGraph.h"
+#include "TLine.h"
+
+#include "../mcore/DrsCalib.h"
+#include "MDrsCalibration.h"
+#include "MExtralgoSpline.h"
+#include "MSequence.h"
+#include "MStatusArray.h"
+#include "MHCamera.h"
+#include "MJob.h"
+#include "MWriteRootFile.h"
+#include "MHCamera.h"
+#include "MBadPixelsCam.h"
+#include "MBadPixelsPix.h"
+#include "MDirIter.h"
+#include "MTaskList.h"
+#include "MFDataPhrase.h"
+#include "MArrayF.h"
+#include "MBadPixelsTreat.h"
+#include "MCalibrateDrsTimes.h"
+#include "MHSectorVsTime.h"
+#include "MHCamEvent.h"
+#include "MExtractTimeAndChargeSpline.h"
+#include "MFillH.h"
+#include "MDrsCalibApply.h"
+#include "MGeomApply.h"
+#include "MContinue.h"
+#include "MRawFitsRead.h"
+#include "MReadMarsFile.h"
+#include "MEvtLoop.h"
+#include "MParList.h"
+#include "MStatusDisplay.h"
+#include "MDrsCalibrationTime.h"
+#include "MH3.h"
+#include "MGeomCamFACT.h"
+#include "MCalibrateFact.h"
+#include "MParameters.h"
+#include "MWriteAsciiFile.h"
+
+#endif
+
+using namespace std;
+
+/* Maybe you wanna use this macro like this:
+ * 
+ * 0.) ---- call root ----
+ *  root -b
+ * 
+ * 1.) ---- compile the stuff ----
+ *  .L fact/analysis/callisto_buildable_no_sequence_file.C++
+ *  <read a lot of warnings>
+ * 
+ * 2.) ---- you can call it then ---- 
+ *      Therefore you need to specify all the paths ... see below.
+ *   
+ * When you wanna call the stuff directly from the bash make sure to 
+ * escape the bracets and quotes correctly.
+ * 
+ * your can do:
+ *  root -b -q callisto_buildable_no_sequence_file.C++'("path1","path2",...)'
+ * or:
+ *  root -b -q callisto_buildable_no_sequence_file.C++(\"path1\",\"$HOME\",...)
+ * using bash enviroment variables like $HOME is not possible in the upper variant.
+ */
+
+int callisto(const TString drsfile="test300samples2.drs.fits.gz",
+             const TString pedfile="00000001.001_P_MonteCarlo000_Events.fits",
+             const TString datfile="00000003.387_D_MonteCarlo010_Events.fits",
+             TString outpath = "",
+             TString displayfile = "", TString displaytitle = "")
+{
+
+    // ======================================================
+
+    if (displaytitle.IsNull())
+        displaytitle = gSystem->BaseName(datfile);
+
+    FileStat_t fstat;
+    int rc = gSystem->GetPathInfo(outpath, fstat);
+    bool isdir = !rc || R_ISDIR(fstat.fMode);
+
+    TString filename = datfile + "_callisto.root";
+    filename.Replace(0, filename.Last('/')+1, "");
+    const char *buf = gSystem->ConcatFileName(outpath, filename);
+    TString outfile = buf;
+    delete [] buf;
+
+    if (displayfile.IsNull())
+    {
+        displayfile = outfile;
+        displayfile.Insert(displayfile.Last('.'), "-display");
+    }
+    else
+    {
+        if (isdir && gSystem->DirName(displayfile)==TString("."))
+        {
+            buf = gSystem->ConcatFileName(outfile, displayfile);
+            displayfile = buf;
+            delete [] buf;
+        }
+    }
+
+    // ======================================================
+
+    // true:  Display correctly mapped pixels in the camera displays
+    //        but the value-vs-index plot is in software/spiral indices
+    // false: Display pixels in hardware/linear indices,
+    //        but the order is the camera display is distorted
+    bool usemap = true;
+
+    // map file to use (get that from La Palma!)
+    const char *pmap = usemap ? "/home/isdc/toscanos/FACT/Mars_svn/resources/FACTmap111030.txt" : NULL;
+
+    Bool_t maximum = kTRUE;
+
+    //const char *lp_template    = maximum ?
+    //    "/cm/shared/apps/fact/Mars_svn_LP/template-lp-extractor-maximum.root" :
+    //    "/cm/shared/apps/fact/Mars_svn_LP/template-lp-extractor-leading-edge.root";
+
+    //    const char *lp_template    = "/home/isdc/toscanos/FACT/Mars_svn/resources/template-lp.root";
+
+    //const char *pulse_template = "/cm/shared/apps/fact/Mars_svn_LP/template-pulse.root";
+
+    const char *pulse_template = "/home/isdc/toscanos/FACT/Mars_svn/resources/template-pulse.root";
+
+    // ------------------------------------------------------
+
+    // Calib: 51 / 90 / 197 (20% TH)
+    // Data:  52 / 64 / 104 (20% TH)
+
+    // Extraction range in slices. It will always(!) contain the full range
+    // of integration
+    const int first_slice =  20; //  10ns
+    const int last_slice  = 250; // 125ns
+
+    // Note that rise and fall time mean different things whether you use IntegralFixed or IntegraRel:
+    //
+    //  IntegralFixed:
+    //    * fRiseTime: Number of slices left  from arrival time
+    //    * fFallTime: Number of slices right from arrival time
+    //  IntegralRel:
+    //    * fRiseTime: Number of slices left  from maximum time
+    //    * fFallTime: Number of slices right from maximum time
+    //
+    const int rise_time_cal = maximum ?  40 :  10; // was 13;   5ns
+    const int fall_time_cal = maximum ? 120 : 160; // was 23;  80ns
+
+    const int rise_time_dat = maximum ?  10 :   2; // was 13; was 10;   1ns
+    const int fall_time_dat = maximum ?  40 :  48; // was 23; was 40;  24ns
+
+    // Extraction type: Extract integral and half leading edge
+
+    const MExtralgoSpline::ExtractionType_t type = maximum ? (MExtralgoSpline::kIntegralRel) : (MExtralgoSpline::kIntegralFixed);
+    //const int type = MExtralgoSpline::kIntegralFixed;
+
+
+    const double heighttm   = 0.5; // IntegralAbs { 1.5pe * 9.6mV/pe } / IntegralRel { 0.5 }
+
+    Long_t max  =    0;  // All
+    Long_t max3 =  max;  // Pedestal Rndm
+    Long_t max4 =  max;  // Pedestal Ext
+
+    // ======================================================
+
+    if (pmap && gSystem->AccessPathName(pmap, kFileExists))
+    {
+        gLog << err << "ERROR - Cannot access mapping file '" << pmap << "'" << endl;
+        return 1;
+    }
+
+    gLog.Separator("Callisto");
+    gLog << all;
+    gLog << "Data File:     " << datfile << '\n';
+    gLog << "DRS calib 300: " << drsfile << endl;;
+
+    MDrsCalibration drscalib300;
+    if (!drscalib300.ReadFits(drsfile.Data())) {
+        gLog << err << "ERROR - Cannot access drscallib300 file '" << drsfile << "'" << endl;
+        return 5;
+    }
+    gLog << all;
+    gLog << "Pedestal file: " << pedfile << '\n';
+    gLog << "Output file:   " << outfile << '\n';
+    gLog << "Display file:  " << displayfile << '\n';
+    gLog << "Display title: " << displaytitle << endl;
+
+    // ------------------------------------------------------
+    MStatusArray arrt, arrp;
+
+    // TFile ft(lp_template);
+    // if (arrt.Read()<=0)
+    // {
+    //     gLog << err << "ERROR - Reading LP template from " << lp_template << endl;
+    //     return 100;
+    // }
+
+    // MHCamera *lpref = (MHCamera*)arrt.FindObjectInCanvas("ExtCalSig;avg", "MHCamera", "Cam");
+    // if (!lpref)
+    // {
+    //     gLog << err << "ERROR - LP Template not found in " << lp_template << endl;
+    //     return 101;
+    // }
+    // lpref->SetDirectory(0);
+
+    // MHCamera *gain = (MHCamera*)arrt.FindObjectInCanvas("gain", "MHCamera", "Gain");
+    // if (!gain)
+    // {
+    //     gLog << err << "ERROR - Gain not found in " << lp_template << endl;
+    //     return 101;
+    // }
+    // gain->SetDirectory(0);
+
+    TFile fp(pulse_template);
+    if (arrp.Read()<=0)
+    {
+        gLog << err << "ERROR - Reading Pulse template from " << pulse_template << endl;
+        return 102;
+    }
+
+    TH1F *hpulse = (TH1F*)arrp.FindObjectInCanvas("hPixelEdgeMean0_0", "TH1F", "cgpPixelPulses0");
+    if (!hpulse)
+    {
+        gLog << err << "ERROR - Pulse Template not found in " << pulse_template << endl;
+        return 103;
+    }
+    hpulse->SetDirectory(0);
+    // ======================================================
+
+    MStatusDisplay *d = new MStatusDisplay;
+
+    MBadPixelsCam badpixels;
+    badpixels.InitSize(1440);
+    /*
+    badpixels[ 424].SetUnsuitable(MBadPixelsPix::kUnsuitable);
+    badpixels[ 583].SetUnsuitable(MBadPixelsPix::kUnsuitable);
+    badpixels[ 830].SetUnsuitable(MBadPixelsPix::kUnsuitable);
+    badpixels[ 923].SetUnsuitable(MBadPixelsPix::kUnsuitable);
+    badpixels[1208].SetUnsuitable(MBadPixelsPix::kUnsuitable);
+    badpixels[1399].SetUnsuitable(MBadPixelsPix::kUnsuitable);
+    */
+    //  Twin pixel
+    //     113
+    //     115
+    //     354
+    //     423
+    //    1195
+    //    1393
+
+    MDrsCalibrationTime timecam;
+
+    // Plot the trigger pattern rates vs. run-number
+    MH3 hrate("MRawRunHeader.GetFileID", "MRawEvtHeader.GetTriggerID&0xff00");
+    hrate.SetWeight("1./TMath::Max(MRawRunHeader.GetRunLength,1)");
+    hrate.SetName("Rate");
+    hrate.SetTitle("Event rate [Hz];File Id;Trigger Type;");
+    hrate.InitLabels(MH3::kLabelsXY);
+    hrate.DefineLabelY(    0, "Data"); // What if TriggerID==0 already???
+    hrate.DefineLabelY(0x100, "Cal");
+    hrate.DefineLabelY(0x400, "Ped");
+    // hrate.DefaultLabelY("ERROR");
+    gStyle->SetOptFit(kTRUE);
+
+
+    // ========================= Result ==================================
+
+    //~ Double_t avgS = evt1f.GetHist()->GetMean();
+    //~ Double_t medS = evt1f.GetHist()->GetMedian();
+    //~ Double_t rmsS = evt1f.GetHist()->GetRMS();
+    //~ Double_t maxS = evt1f.GetHist()->GetMaximum();
+
+    MArrayF der1(hpulse->GetNbinsX());
+    MArrayF der2(hpulse->GetNbinsX());
+
+    MExtralgoSpline spline(hpulse->GetArray()+1, hpulse->GetNbinsX(),
+                           der1.GetArray(), der2.GetArray());
+    spline.SetRiseFallTime(rise_time_dat, fall_time_dat);
+    spline.SetExtractionType(type);
+    spline.SetHeightTm(heighttm);
+
+    spline.Extract(hpulse->GetMaximumBin()-1);
+
+    // The pulser signal is most probably around 400mV/9.5mV
+    // IntegraFixed 2/48 corresponds to roughly 215mV*50slices
+    Double_t scale = 1./spline.GetSignal();
+
+    MArrayD calib(1440);
+    for (int i=0; i<1440; i++)
+        calib[i] =1.;
+
+    gROOT->SetSelectedPad(0);
+    d->AddTab("PulseTemp");
+    gPad->SetGrid();
+    hpulse->SetNameTitle("Pulse", "Single p.e. pulse template");
+    hpulse->SetDirectory(0);
+    hpulse->SetLineColor(kBlack);
+    hpulse->DrawCopy();
+
+    TAxis *ax = hpulse->GetXaxis();
+
+    Double_t w = hpulse->GetBinWidth(1);
+    Double_t T = w*(spline.GetTime()+0.5)       +ax->GetXmin();
+    //~ Double_t H = w*(hpulse->GetMaximumBin()+0.5)+ax->GetXmin();
+
+    TLine line;
+    line.SetLineColor(kRed);
+    line.DrawLine(T-rise_time_dat*w, spline.GetHeight(),
+                  T+fall_time_dat*w, spline.GetHeight());
+    line.DrawLine(T, spline.GetHeight()/4, T, 3*spline.GetHeight()/4);
+    line.DrawLine(T-rise_time_dat*w, 0,
+                  T-rise_time_dat*w, spline.GetHeight());
+    line.DrawLine(T+fall_time_dat*w, 0,
+                  T+fall_time_dat*w, spline.GetHeight());
+
+    TGraph gg;
+    for (int ix=1; ix<=hpulse->GetNbinsX(); ix++)
+        for (int i=0; i<10; i++)
+        {
+            Double_t x = hpulse->GetBinLowEdge(ix)+i*hpulse->GetBinWidth(ix)/10.;
+            gg.SetPoint(gg.GetN(), x+w/2, spline.EvalAt(ix-1+i/10.));
+        }
+
+    gg.SetLineColor(kBlue);
+    gg.SetMarkerColor(kBlue);
+    gg.SetMarkerStyle(kFullDotMedium);
+    gg.DrawClone("L");
+
+    gROOT->SetSelectedPad(0);
+    d->AddTab("CalConst");
+    MGeomCamFACT fact;
+    MHCamera hcalco(fact);
+    hcalco.SetName("CalConst");
+    hcalco.SetTitle(Form("Relative calibration constant [%.0f/pe]", 1./scale));
+    hcalco.SetCamContent(calib);
+    hcalco.SetAllUsed();
+    //hcalco.Scale(scale);
+    hcalco.DrawCopy();
+
+    // ======================================================
+
+    gLog << endl;
+    gLog.Separator("Extracting random pedestal");
+
+    MTaskList tlist3;
+
+    MParList plist3;
+    plist3.AddToList(&tlist3);
+    plist3.AddToList(&drscalib300);
+    plist3.AddToList(&badpixels);
+    plist3.AddToList(&timecam);
+
+    MEvtLoop loop3("DetermineRndmPed");
+    loop3.SetDisplay(d);
+    loop3.SetParList(&plist3);
+
+    // ------------------ Setup the tasks ---------------
+
+    MRawFitsRead read3;
+    read3.LoadMap(pmap);
+    read3.AddFile(pedfile);
+
+    MFillH fill3a(&hrate);
+
+    MContinue cont3("(MRawEvtHeader.GetTriggerID&0xff00)!=0x400", "SelectPed");
+
+    MGeomApply apply3;
+
+    MDrsCalibApply drsapply3;
+
+    //MFilterData filterdata5;
+
+    //---
+
+    MExtractTimeAndChargeSpline extractor3;
+    extractor3.SetRange(first_slice, last_slice);
+    extractor3.SetRiseTimeHiGain(rise_time_dat);
+    extractor3.SetFallTimeHiGain(fall_time_dat);
+    extractor3.SetHeightTm(heighttm);
+    extractor3.SetChargeType(type);
+    extractor3.SetSaturationLimit(600000);
+    extractor3.SetNoiseCalculation(kTRUE);
+
+//    MHCamEvent evt2a(0, "PedRdm", "Extracted Pedestal Signal;;S");
+
+//    MFillH fill2a(&evt2a, "MExtractedSignalCam", "FillPedRndm");
+
+    // Use this for data, but not for calibration events
+//    evt2a.SetErrorSpread(kFALSE);
+
+    /*
+     MCalibrateData conv3;
+     conv3.SetCalibrationMode(MCalibrateData::kNone);
+     conv3.SetPedestalFlag(MCalibrateData::kNo);
+     conv3.SetCalibConvMinLimit(0);
+     conv3.SetCalibConvMaxLimit(10000);
+     conv3.SetScaleFactor(scale);
+     */
+
+    MCalibrateFact conv3;
+    conv3.SetScale(scale);
+    conv3.SetCalibConst(calib);
+
+    MBadPixelsTreat treat3;
+    treat3.SetProcessPedestalRun(kFALSE);
+    treat3.SetProcessPedestalEvt(kFALSE);
+    treat3.SetProcessTimes(kFALSE);
+
+    MHCamEvent evt3b(0, "PedRdm","Interpolated random pedestal;;Signal [~phe]");
+    //evt2b.SetErrorSpread(kFALSE);
+
+    MFillH fill3b(&evt3b, "MSignalCam", "FillPedRdm");
+    fill3b.SetDrawOption("gaus");
+
+    // ------------------ Setup eventloop and run analysis ---------------
+
+    tlist3.AddToList(&read3);
+    tlist3.AddToList(&apply3);
+    tlist3.AddToList(&drsapply3);
+    tlist3.AddToList(&cont3);
+    tlist3.AddToList(&extractor3);
+//    tlist3.AddToList(&fill3a);
+    tlist3.AddToList(&conv3);
+    tlist3.AddToList(&treat3);
+    tlist3.AddToList(&fill3b);
+
+    if (!loop3.Eventloop(max3))
+        return 14;
+
+    if (!loop3.GetDisplay())
+        return 15;
+
+    // ======================================================
+
+    gLog << endl;
+    gLog.Separator("Extracting pedestal");
+
+    MTaskList tlist4;
+
+    MParList plist4;
+    plist4.AddToList(&tlist4);
+    plist4.AddToList(&drscalib300);
+    plist4.AddToList(&badpixels);
+    plist4.AddToList(&timecam);
+
+    MEvtLoop loop4("DetermineExtractedPed");
+    loop4.SetDisplay(d);
+    loop4.SetParList(&plist4);
+
+    // ------------------ Setup the tasks ---------------
+
+    MRawFitsRead read4;
+    read4.LoadMap(pmap);
+    read4.AddFile(pedfile);
+
+    MContinue cont4("(MRawEvtHeader.GetTriggerID&0xff00)!=0x400", "SelectPed");
+
+    MGeomApply apply4;
+
+    MDrsCalibApply drsapply4;
+
+    MExtractTimeAndChargeSpline extractor4;
+    extractor4.SetRange(first_slice, last_slice);
+    extractor4.SetRiseTimeHiGain(rise_time_dat);
+    extractor4.SetFallTimeHiGain(fall_time_dat);
+    extractor4.SetHeightTm(heighttm);
+    extractor4.SetChargeType(type);
+    extractor4.SetSaturationLimit(600000);
+    extractor4.SetNoiseCalculation(kFALSE);
+
+    //    MHCamEvent evt3a(0, "PedExt", "Extracted Pedestal Signal;;S");
+
+    //    MFillH fill3a(&evt3a, "MExtractedSignalCam", "FillPedExt");
+
+    // Use this for data, but not for calibration events
+//    evt3a.SetErrorSpread(kFALSE);
+/*
+    MCalibrateData conv4;
+    conv4.SetCalibrationMode(MCalibrateData::kNone);
+    conv4.SetPedestalFlag(MCalibrateData::kNo);
+    conv4.SetCalibConvMinLimit(0);
+    conv4.SetCalibConvMaxLimit(10000);
+    conv4.SetScaleFactor(scale);
+*/
+    MCalibrateFact conv4;
+    conv4.SetScale(scale);
+    conv4.SetCalibConst(calib);
+
+    MBadPixelsTreat treat4;
+    treat4.SetProcessPedestalRun(kFALSE);
+    treat4.SetProcessPedestalEvt(kFALSE);
+
+    MHCamEvent evt4b(0, "PedExt","Interpolated extracted pedestal;;Signal [~phe]");
+    //evt4b.SetErrorSpread(kFALSE);
+
+    MFillH fill4b(&evt4b, "MSignalCam", "FillPedExt");
+    fill4b.SetDrawOption("gaus");
+
+    // ------------------ Setup eventloop and run analysis ---------------
+
+    tlist4.AddToList(&read4);
+    tlist4.AddToList(&apply4);
+    tlist4.AddToList(&drsapply4);
+    tlist4.AddToList(&cont4);
+    tlist4.AddToList(&extractor4);
+//    tlist4.AddToList(&fill4a);
+    tlist4.AddToList(&conv4);
+    tlist4.AddToList(&treat4);
+    tlist4.AddToList(&fill4b);
+
+    if (!loop4.Eventloop(max4))
+        return 15;
+
+    if (!loop4.GetDisplay())
+        return 16;
+
+    // ===================================================================
+
+    gLog << endl;
+    gLog.Separator("Extracting and calibration data");
+
+    MTaskList tlist5;
+
+    MParList plist5;
+    plist5.AddToList(&tlist5);
+    plist5.AddToList(&drscalib300);
+    plist5.AddToList(&badpixels);
+    plist5.AddToList(&timecam);
+
+    MEvtLoop loop5("CalibratingData");
+    loop5.SetDisplay(d);
+    loop5.SetParList(&plist5);
+
+    // ------------------ Setup the tasks ---------------
+
+    MRawFitsRead read5a;
+    MReadMarsFile read5b("Events");
+    read5a.LoadMap(pmap);
+    read5a.AddFile(datfile);
+    read5b.DisableAutoScheme();
+    read5b.AddFile(datfile);
+
+    MRead &read5 = datfile.EndsWith(".root") ? static_cast<MRead&>(read5b) : static_cast<MRead&>(read5a);
+
+    MFillH fill5a(&hrate);
+
+    MGeomApply apply5;
+
+    MDrsCalibApply drsapply5;
+
+    MTreatSaturation treatsat5;
+
+    MFDataPhrase filterdat("(MRawEvtHeader.GetTriggerID&0xff00)==0",     "SelectDat");
+    MFDataPhrase filtercal("(MRawEvtHeader.GetTriggerID&0xff00)==0x100", "SelectCal");
+    MFDataPhrase filterped("(MRawEvtHeader.GetTriggerID&0xff00)==0x400", "SelectPed");
+    MFDataPhrase filterncl("(MRawEvtHeader.GetTriggerID&0xff00)!=0x100", "SelectNonCal");
+
+    //MContinue cont4("MRawEvtHeader.GetTriggerID!=4", "SelectData");
+
+    // ---
+
+    MExtractTimeAndChargeSpline extractor5dat;
+    extractor5dat.SetRange(first_slice, last_slice);
+    extractor5dat.SetRiseTimeHiGain(rise_time_dat);
+    extractor5dat.SetFallTimeHiGain(fall_time_dat);
+    extractor5dat.SetHeightTm(heighttm);
+    extractor5dat.SetChargeType(type);
+    extractor5dat.SetSaturationLimit(600000);
+    extractor5dat.SetNoiseCalculation(kFALSE);
+
+    MExtractTimeAndChargeSpline extractor5cal;
+    extractor5cal.SetRange(first_slice, last_slice);
+    extractor5cal.SetRiseTimeHiGain(rise_time_cal);
+    extractor5cal.SetFallTimeHiGain(fall_time_cal);
+    extractor5cal.SetHeightTm(heighttm);
+    extractor5cal.SetChargeType(type);
+    extractor5cal.SetSaturationLimit(600000);
+    extractor5cal.SetNoiseCalculation(kFALSE);
+
+    MExtractTimeAndChargeSpline extractor5tm("ExtractTM");
+    extractor5tm.SetRange(last_slice, 294);
+    extractor5tm.SetRiseTimeHiGain(1);
+    extractor5tm.SetFallTimeHiGain(1);
+    extractor5tm.SetHeightTm(0.5);
+    extractor5tm.SetChargeType(MExtralgoSpline::kAmplitudeRel);
+    extractor5tm.SetSaturationLimit(600000);
+    extractor5tm.SetNoiseCalculation(kFALSE);
+    extractor5tm.SetNameSignalCam("TimeMarkerAmplitude");
+    extractor5tm.SetNameTimeCam("TimeMarkerTime");
+
+    extractor5dat.SetFilter(&filterncl);
+    extractor5cal.SetFilter(&filtercal);
+    //extractor4tm.SetFilter(&filtercal);
+
+    // ---
+/*
+    MCalibrateData conv5;
+    conv5.SetCalibrationMode(MCalibrateData::kNone);
+    conv5.SetPedestalFlag(MCalibrateData::kNo);
+    conv5.SetCalibConvMinLimit(0);
+    conv5.SetCalibConvMaxLimit(10000);
+    conv5.SetScaleFactor(scale);
+*/
+    MCalibrateFact conv5;
+    conv5.SetScale(scale);
+    conv5.SetCalibConst(calib);
+
+    MCalibrateDrsTimes calctm5;
+    calctm5.SetNameUncalibrated("UncalibratedTimes");
+
+    MCalibrateDrsTimes calctm5tm("CalibrateTimeMarker");
+    calctm5tm.SetNameArrivalTime("TimeMarkerTime");
+    calctm5tm.SetNameUncalibrated("UncalTimeMarker");
+    calctm5tm.SetNameCalibrated("TimeMarker");
+    calctm5tm.SetTimeMarker();
+    //calctm4tm.SetFilter(&filtercal);
+
+    MBadPixelsTreat treat5;
+    treat5.SetProcessPedestalRun(kFALSE);
+    treat5.SetProcessPedestalEvt(kFALSE);
+
+    MHCamEvent evt5b(0, "ExtSig",   "Extracted signal;;S [mV·sl]");
+    MHCamEvent evt5c(0, "CalSig",   "Calibrated and interpolated signal;;S [~phe]");
+    MHCamEvent evt5d(4, "ExtSigTm", "Extracted time;;T [sl]");
+    MHCamEvent evt5e(6, "CalSigTm", "Calibrated and interpolated time;;T [ns]");
+
+    MFillH fill5b(&evt5b, "MExtractedSignalCam", "FillExtSig");
+    MFillH fill5c(&evt5c, "MSignalCam",          "FillCalSig");
+    MFillH fill5d(&evt5d, "MArrivalTimeCam",     "FillExtTm");
+    MFillH fill5e(&evt5e, "MSignalCam",          "FillCalTm");
+
+    fill5c.SetDrawOption("gaus");
+    fill5d.SetDrawOption("gaus");
+    fill5e.SetDrawOption("gaus");
+
+    /*
+    fill4b.SetFilter(&filterdat);
+    fill4c.SetFilter(&filterdat);
+    fill4d.SetFilter(&filterdat);
+    fill4e.SetFilter(&filterdat);
+    */
+
+    //MFSoftwareTrigger swtrig;
+    //MContinue contsw(&swtrig, "FilterSwTrigger", "Software trigger");
+    //contsw.SetInverted();
+
+    // The second rule is for the case reading raw-files!
+    
+    MWriteRootFile write5(outfile, "RECREATE", "Calibrated Data", 2);
+    write5.AddContainer("MRawRunHeader",       "RunHeaders");
+    write5.AddContainer("MGeomCam",            "RunHeaders");
+    write5.AddContainer("MMcCorsikaRunHeader", "RunHeaders", kFALSE);
+    write5.AddContainer("MCorsikaRunHeader",   "RunHeaders", kFALSE);
+    write5.AddContainer("MMcRunHeader",        "RunHeaders", kFALSE);
+
+    // Common events
+    write5.AddContainer("MCorsikaEvtHeader",   "Events", kFALSE);
+    write5.AddContainer("MMcEvt",              "Events", kFALSE);
+    write5.AddContainer("IncidentAngle",       "Events", kFALSE);
+    write5.AddContainer("MPointingPos",        "Events", kFALSE);
+    write5.AddContainer("MSignalCam",          "Events");
+    write5.AddContainer("MTime",               "Events", kFALSE);
+    write5.AddContainer("MRawEvtHeader",       "Events");
+    //write.AddContainer("MTriggerPattern", "Events");
+    
+    // ------------------ Setup histograms and fill tasks ----------------
+
+    MContinue test;
+    test.SetFilter(&filterncl);
+
+    MTaskList tlist5tm;
+    tlist5tm.AddToList(&extractor5tm);
+    tlist5tm.AddToList(&calctm5tm);
+    tlist5tm.SetFilter(&filtercal);
+
+    MTaskList tlist5dat;
+    tlist5dat.AddToList(&fill5b);
+    tlist5dat.AddToList(&fill5c);
+    tlist5dat.AddToList(&fill5d);
+    tlist5dat.AddToList(&fill5e);
+    tlist5dat.SetFilter(&filterdat);
+
+    tlist5.AddToList(&read5);
+    tlist5.AddToList(&apply5);
+    tlist5.AddToList(&drsapply5);
+    tlist5.AddToList(&treatsat5);
+    tlist5.AddToList(&filterncl);
+    //tlist5.AddToList(&test);
+    tlist5.AddToList(&filterdat);
+    tlist5.AddToList(&filtercal);
+    tlist5.AddToList(&filterped);
+    tlist5.AddToList(&fill5a);
+    tlist5.AddToList(&extractor5dat);
+    tlist5.AddToList(&extractor5cal);
+    tlist5.AddToList(&calctm5);
+    tlist5.AddToList(&tlist5tm);
+    tlist5.AddToList(&conv5);
+    tlist5.AddToList(&treat5);
+    tlist5.AddToList(&tlist5dat);
+    tlist5.AddToList(&write5);
+
+    if (!loop5.Eventloop(max4))
+        return 18;
+
+    if (!loop5.GetDisplay())
+        return 19;
+
+    d->SetTitle(displaytitle, kFALSE);
+    d->SaveAs(displayfile);
+
+    return 0;
+}
Index: branches/Mars_McMismatchStudy/fact/analysis/mc/ganymed.C
===================================================================
--- branches/Mars_McMismatchStudy/fact/analysis/mc/ganymed.C	(revision 17960)
+++ branches/Mars_McMismatchStudy/fact/analysis/mc/ganymed.C	(revision 17962)
@@ -169,4 +169,5 @@
     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) )");
+    //    energyest.SetRule("1737 * (MHillas.fSize / 300.)^0.810");
 
     // Cuts before writing
Index: branches/Mars_McMismatchStudy/fact/analysis/mc/ganymed.C.ORIGINAL
===================================================================
--- branches/Mars_McMismatchStudy/fact/analysis/mc/ganymed.C.ORIGINAL	(revision 17962)
+++ branches/Mars_McMismatchStudy/fact/analysis/mc/ganymed.C.ORIGINAL	(revision 17962)
@@ -0,0 +1,525 @@
+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;
+}
Index: branches/Mars_McMismatchStudy/fact/analysis/mc/star.C.NEWCLEANING
===================================================================
--- branches/Mars_McMismatchStudy/fact/analysis/mc/star.C.NEWCLEANING	(revision 17962)
+++ branches/Mars_McMismatchStudy/fact/analysis/mc/star.C.NEWCLEANING	(revision 17962)
@@ -0,0 +1,420 @@
+#include <sstream>
+#include <iostream>
+
+#include "MLog.h"
+#include "MLogManip.h"
+
+#if !defined(__CINT__) || defined(__MAKECINT__)
+
+#include "TH1F.h"
+#include "TFile.h"
+#include "TStyle.h"
+#include "TGraph.h"
+#include "TLine.h"
+
+#include "../mcore/DrsCalib.h"
+#include "MDrsCalibration.h"
+//#include "MLogManip.h"
+#include "MExtralgoSpline.h"
+#include "MSequence.h"
+#include "MStatusArray.h"
+#include "MHCamera.h"
+#include "MJob.h"
+#include "MWriteRootFile.h"
+#include "MHCamera.h"
+#include "MBadPixelsCam.h"
+#include "MBadPixelsPix.h"
+#include "MDirIter.h"
+#include "MTaskList.h"
+#include "MFDataPhrase.h"
+#include "MArrayF.h"
+#include "MBadPixelsTreat.h"
+#include "MCalibrateDrsTimes.h"
+#include "MHSectorVsTime.h"
+#include "MHCamEvent.h"
+#include "MExtractTimeAndChargeSpline.h"
+#include "MFillH.h"
+#include "MDrsCalibApply.h"
+#include "MGeomApply.h"
+#include "MContinue.h"
+#include "MRawFitsRead.h"
+#include "MEvtLoop.h"
+#include "MParList.h"
+#include "MStatusDisplay.h"
+#include "MDrsCalibrationTime.h"
+#include "MH3.h"
+#include "MGeomCamFACT.h"
+#include "MCalibrateFact.h"
+#include "MParameters.h"
+#include "MWriteAsciiFile.h"
+#include "MFMagicCuts.h"
+
+#include "MMuonSetup.h"
+#include "MReadMarsFile.h"
+#include "MHillasCalc.h"
+#include "MHn.h"
+#include "MMuonSearchParCalc.h"
+#include "MMuonCalibParCalc.h"
+#include "MBinning.h"
+#include "MImgCleanTime.h"
+
+
+#endif
+
+using namespace std;
+
+int star(const TString callistofile, TString outfile,
+         TString displayfile, TString displaytitle,
+         const TString txtfile,
+         Double_t lvl1=4.0, Double_t lvl2=2.5, Double_t deltat = 17.5);
+
+
+int star(const TString infile = "callisto.root", TString outfile = "",
+         Double_t lvl1=4.0,  Double_t lvl2=2.5, Double_t deltat = 17.5)
+{
+    return star(infile, outfile, "", "", "", lvl1, lvl2, deltat);
+}
+
+int star(const TString callistofile, TString outfile,
+         TString displayfile,  TString displaytitle,
+         const TString txtfile,
+         Double_t lvl1, Double_t lvl2, Double_t deltat)
+{
+
+    //Double_t lvl1=7.8, 
+    //Double_t lvl2=3.9,
+    //    lvl1 = 2.5;
+    //    lvl2 = 0.5;
+    // ------------------------------------------------------
+
+    if (displaytitle.IsNull())
+        displaytitle = gSystem->BaseName(callistofile);
+
+    FileStat_t fstat;
+    int rc = gSystem->GetPathInfo(outfile, fstat);
+    bool isdir = !rc || R_ISDIR(fstat.fMode);
+
+    char *buf = gSystem->ConcatFileName(outfile, "star.root");
+    outfile = buf;
+    delete [] buf;
+
+    if (displayfile.IsNull())
+    {
+        displayfile = outfile;
+        displayfile.Insert(displayfile.Last('.'), "-display");
+    }
+    else
+    {
+        if (isdir && gSystem->DirName(displayfile)==TString("."))
+        {
+            buf = gSystem->ConcatFileName(outfile, displayfile);
+            displayfile = buf;
+            delete [] buf;
+        }
+    }
+
+    // ------------------------------------------------------
+
+    gLog.Separator("Star");
+    gLog << all << "Calculate image parameters of sequence ";
+    gLog << endl;
+    gLog << "Callisto file: " << callistofile << endl;
+    gLog << "Output file:   " << outfile << endl;
+    gLog << "Display file:  " << displayfile << endl;
+    gLog << "Display title: " << displaytitle << endl;
+    gLog << endl;
+    gLog << "lvl1: "   << lvl1 <<endl;
+    gLog << "lvl2: "   << lvl2 <<endl;
+    gLog << "deltat: " << deltat <<endl;
+    
+    
+    gLog << endl;
+
+    // ------------------------------------------------------
+
+    gLog.Separator();
+
+    // --------------------------------------------------------
+
+    MBinning bins1(  36,    -90,    90, "BinningAlpha");
+    MBinning bins3(  67, -0.005, 0.665, "BinningTheta", "asin");
+    MBinning bins4(  25,      0,   2.5, "BinningDist");
+    MBinning binsM1(100,      0,     5, "BinningMuonRadius");
+    MBinning binsM2( 60,      0,   0.3, "BinningMuonDeviation");
+    MBinning binsM3(150,      0,    60, "BinningMuonTime");
+    MBinning binsM4(300,      0,    30, "BinningMuonTimeRms");
+    MBinning binsM5( 20,      0,   360, "BinningMuonArcPhi");
+    MBinning binsM6( 50,      0,   0.5, "BinningMuonArcWidth");
+    MBinning binsM7( 30,      5,  5000, "BinningMuonSize", "log");
+    MBinning binsM8(100,      0,     5, "BinningMuonRelTimeSigma");
+    MBinning binsM9(100,      0,     5, "BinningRadius");
+
+    // ---------------------------------------------------------
+
+    // Filter to start muon analysis
+    //MFDataPhrase fmuon1("MHillas.fSize>150 && MNewImagePar.fConcCOG<0.1", "MuonPreCut");
+    // Filter to calculate further muon parameters
+//    MFDataPhrase fmuon2("(MMuonSearchPar.fRadius*MGeomCam.fConvMm2Deg>0.6) && (MMuonSearchPar.fRadius*MGeomCam.fConvMm2Deg<1.35) &&"
+//                        "(MMuonSearchPar.fDeviation*MGeomCam.fConvMm2Deg<0.152)", "MuonSearchCut");
+    MFDataPhrase fmuon2("MMuonSearchPar.fDeviation*MGeomCam.fConvMm2Deg>0.015 &&"
+                        "MMuonSearchPar.fDeviation*MGeomCam.fConvMm2Deg>(MMuonSearchPar.fRadius*MGeomCam.fConvMm2Deg-0.1)*0.15/0.4",
+                        "MuonSearchCut");
+
+    // Filter to fill the MHMuonPar
+    MFDataPhrase fmuon3("MMuonCalibPar.fArcPhi>190 && MMuonCalibPar.fRelTimeSigma<3.2",
+                        "MuonFinalCut");
+    // Filter to write Muons to Muon tree
+    //MFDataMember fmuon4("MMuonCalibPar.fArcPhi", '>', -0.5, "MuonWriteCut");
+
+    // ---------------------------------------------------------
+
+    MStatusDisplay *d = new MStatusDisplay;
+
+    MTaskList tlist;
+
+    MParList plist2;
+    plist2.AddToList(&tlist);
+    plist2.AddToList(&bins1);
+    plist2.AddToList(&bins3);
+    plist2.AddToList(&bins4);
+    plist2.AddToList(&binsM1);
+    plist2.AddToList(&binsM2);
+    plist2.AddToList(&binsM3);
+    plist2.AddToList(&binsM4);
+    plist2.AddToList(&binsM5);
+    plist2.AddToList(&binsM6);
+    plist2.AddToList(&binsM7);
+    plist2.AddToList(&binsM8);
+    plist2.AddToList(&binsM9);
+
+    MMuonSetup muonsetup;
+    plist2.AddToList(&muonsetup);
+
+    MEvtLoop loop;
+    loop.SetDisplay(d);
+    loop.SetParList(&plist2);
+
+    MReadMarsFile read("Events");
+    read.DisableAutoScheme();
+    read.AddFile(callistofile);
+
+    MContinue cont("MRawEvtHeader.GetTriggerID!=4", "SelectData");
+
+    MGeomApply apply;
+
+    MImgCleanTime clean;
+    clean.SetMinCount(0);
+    clean.SetMinSize(25);
+    clean.SetDeltaT(1.25*17.5*0.1111);
+
+    MHillasCalc hcalc;
+    hcalc.Disable(MHillasCalc::kCalcConc);
+
+    MH3 hrate("MRawRunHeader.GetFileID"/*, "MRawEvtHeader.GetTriggerID"*/);
+    hrate.SetWeight("1./TMath::Max(MRawRunHeader.GetRunLength,1)");
+    hrate.SetName("Rate");
+    hrate.SetTitle("Event rate after cleaning;File Id;Event Rate [Hz];");
+    hrate.InitLabels(MH3::kLabelsX);
+    hrate.GetHist().SetMinimum(0);
+
+    MFillH fillR(&hrate, "", "FillRate");
+
+    // ------------------ Setup histograms and fill tasks ----------------
+    MHCamEvent evtC1(0, "Cleaned",  "Average signal after Cleaning;;S [phe]");
+    MHCamEvent evtC2(0, "UsedPix",  "Fraction of Events in which Pixels are used;;Fraction");
+    MHCamEvent evtC3(8, "PulsePos", "Pulse Position of signal after cleaning;;T [ns]");
+    evtC1.SetErrorSpread(kFALSE);
+    evtC2.SetErrorSpread(kFALSE);
+    evtC2.SetThreshold(0);
+
+    MFillH fillC1(&evtC1, "MSignalCam", "FillSignalCam");
+    MFillH fillC2(&evtC2, "MSignalCam", "FillCntUsedPixels");
+    MFillH fillC3(&evtC3, "MSignalCam", "FillPulsePos");
+
+    MFillH fillD1("MHHillas",      "MHillas",      "FillHillas");
+    MFillH fillD2("MHHillasExt",   "",             "FillHillasExt");
+    MFillH fillD3("MHHillasSrc",   "MHillasSrc",   "FillHillasSrc");
+    MFillH fillD4("MHImagePar",    "MImagePar",    "FillImagePar");
+    MFillH fillD5("MHNewImagePar", "MNewImagePar", "FillNewImagePar");
+    MFillH fillD6("MHVsSize",      "MHillasSrc",   "FillVsSize");
+
+    MH3 halpha("MHillasSrc.fDist*MGeomCam.fConvMm2Deg", "fabs(MHillasSrc.fAlpha)");
+    halpha.SetName("AvsD;Dist;Alpha");
+
+    MFillH filla(&halpha);
+    filla.SetNameTab("AvsD");
+    filla.SetDrawOption("colz");
+
+    // ----------------------- Muon Analysis ----------------------
+    //writem.SetFilter(&fmuon4);
+
+    MHn mhn1;
+    mhn1.AddHist("MMuonSearchPar.fRadius*MGeomCam.fConvMm2Deg");
+    mhn1.InitName("MuonRadius");
+    mhn1.InitTitle("MuonRadius");
+
+    mhn1.AddHist("MMuonSearchPar.fDeviation*MGeomCam.fConvMm2Deg");
+    mhn1.InitName("MuonDeviation");
+    mhn1.InitTitle("MuonDeviation");
+
+    mhn1.AddHist("MMuonSearchPar.fTime");
+    mhn1.InitName("MuonTime");
+    mhn1.InitTitle("MuonTime");
+
+    mhn1.AddHist("MMuonSearchPar.fTimeRms");
+    mhn1.InitName("MuonTimeRms");
+    mhn1.InitTitle("MuonTimeRms");
+
+    MHn mhn3;
+    mhn3.AddHist("MMuonSearchPar.fRadius*MGeomCam.fConvMm2Deg","MMuonSearchPar.fDeviation*MGeomCam.fConvMm2Deg");
+    mhn3.InitName("MuonRadius;MuonRadius;MuonDeviation");
+    mhn3.InitTitle("MuonRadius");
+    mhn3.SetDrawOption("colz");
+
+    MHn mhn2;
+    mhn2.AddHist("MMuonCalibPar.fArcPhi");
+    mhn2.InitName("MuonArcPhi");
+    mhn2.InitTitle("MuonArcPhi");
+
+    mhn2.AddHist("MMuonCalibPar.fArcWidth");
+    mhn2.InitName("MuonArcWidth");
+    mhn2.InitTitle("MuonArcWidth");
+
+    mhn2.AddHist("MMuonCalibPar.fMuonSize");
+    mhn2.InitName("MuonSize");
+    mhn2.InitTitle("MuonSize");
+
+    mhn2.AddHist("MMuonCalibPar.fRelTimeSigma");
+    mhn2.InitName("MuonRelTimeSigma");
+    mhn2.InitTitle("MuonRelTimeSigma");
+
+    MHn mhn4;
+    mhn4.AddHist("MMuonCalibPar.fArcPhi","MMuonCalibPar.fArcWidth");
+    mhn4.InitName("MuonArcPhi;MuonArcPhi;MuonArcWidth");
+    mhn4.InitTitle("MuonArcPhi");
+    mhn4.SetDrawOption("colz");
+
+    MFillH fillmhn1(&mhn1, "", "FillMuonSearchPar");
+    MFillH fillmhn2(&mhn2, "", "FillMuonCalibPar");
+    MFillH fillmhn3(&mhn3, "", "FillMuonDevVsRadius");
+    MFillH fillmhn4(&mhn4, "", "FillMuonWidthVsPhi");
+    fillmhn1.SetNameTab("MuonS");
+    fillmhn3.SetNameTab("DevVsRad");
+    fillmhn2.SetNameTab("MuonC");
+    fillmhn4.SetNameTab("WidthVsPhi");
+    fillmhn2.SetFilter(&fmuon2);
+    fillmhn4.SetFilter(&fmuon2);
+
+    MMuonSearchParCalc muscalc;
+    //muscalc.SetFilter(&fmuon1);
+
+    MMuonCalibParCalc mcalc;
+    mcalc.SetFilter(&fmuon2);
+
+    MFillH fillmuon("MHSingleMuon", "", "FillMuon");
+    MFillH fillmpar("MHMuonPar",    "", "FillMuonPar");
+    fillmuon.SetFilter(&fmuon2);
+    fillmpar.SetFilter(&fmuon3);
+    fillmuon.SetBit(MFillH::kDoNotDisplay);
+
+    // ---------------------------------------------------------------
+
+    MWriteAsciiFile writeascii(txtfile, "MHillas");
+    writeascii.AddColumns("MHillasExt");
+    writeascii.AddColumns("MHillasSrc");
+
+
+    MWriteRootFile write(outfile, "RECREATE", "Image parameters", 2);
+    write.AddContainer("MTime",               "Events", kFALSE);
+    write.AddContainer("MHillas",             "Events");
+    write.AddContainer("MHillasExt",          "Events");
+    write.AddContainer("MHillasSrc",          "Events");
+    write.AddContainer("MImagePar",           "Events");
+    write.AddContainer("MNewImagePar",        "Events");
+    write.AddContainer("MRawEvtHeader",       "Events");
+    write.AddContainer("ThetaSquared",        "Events");
+    write.AddContainer("Disp",                "Events");
+    write.AddContainer("MCorsikaEvtHeader",   "Events", kFALSE);
+    write.AddContainer("MMcEvt",              "Events", kFALSE);
+    write.AddContainer("IncidentAngle",       "Events", kFALSE);
+    write.AddContainer("MPointingPos",        "Events", kFALSE);
+
+    write.AddContainer("MRawRunHeader",       "RunHeaders");
+    write.AddContainer("MGeomCam",            "RunHeaders");
+    write.AddContainer("MMcCorsikaRunHeader", "RunHeaders", kFALSE);
+    write.AddContainer("MCorsikaRunHeader",   "RunHeaders", kFALSE);
+    write.AddContainer("MMcRunHeader",        "RunHeaders", kFALSE);
+
+    MFDataPhrase fmuonhn("MMuonCalibPar.fRelTimeSigma>=0",
+                         "MuonHistCut");
+    fillmhn2.SetFilter(&fmuonhn);
+
+    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;
+    MFMagicCuts cuts;
+    cuts.SetHadronnessCut(MFMagicCuts::kNoCut);
+    cuts.SetVariables(param);
+
+
+    tlist.AddToList(&read);
+    tlist.AddToList(&cont);
+    tlist.AddToList(&apply);
+    tlist.AddToList(&clean);
+    tlist.AddToList(&hcalc);
+    tlist.AddToList(&cuts);
+    tlist.AddToList(&fillC1);
+    tlist.AddToList(&fillC2);
+    tlist.AddToList(&fillC3);
+    tlist.AddToList(&fillR);
+    tlist.AddToList(&fillD1);
+    tlist.AddToList(&fillD2);
+    tlist.AddToList(&fillD3);
+    tlist.AddToList(&fillD4);
+    tlist.AddToList(&fillD5);
+    tlist.AddToList(&fillD6);
+    tlist.AddToList(&filla);
+    //tlist.AddToList(&fmuon1);
+    tlist.AddToList(&muscalc);
+    tlist.AddToList(&fillmhn1);
+    tlist.AddToList(&fillmhn3);
+    tlist.AddToList(&fmuon2);
+    tlist.AddToList(&fillmuon);
+    tlist.AddToList(&mcalc);
+    tlist.AddToList(&fmuonhn);
+    tlist.AddToList(&fillmhn2);
+    tlist.AddToList(&fillmhn4);
+    tlist.AddToList(&fmuon3);
+    tlist.AddToList(&fillmpar);
+    //tlist.AddToList(&fmuon4);
+    //tlist.AddToList(&writem);
+    tlist.AddToList(&write);
+    if (!txtfile.IsNull())
+        tlist.AddToList(&writeascii);
+
+    if (!loop.Eventloop())
+        return 3;
+
+    if (!loop.GetDisplay())
+        return 4;
+
+    d->SetTitle(displaytitle, kFALSE);
+    d->SaveAs(displayfile);
+
+    return 0;
+}
+
Index: branches/Mars_McMismatchStudy/fact/analysis/star.C
===================================================================
--- branches/Mars_McMismatchStudy/fact/analysis/star.C	(revision 17960)
+++ branches/Mars_McMismatchStudy/fact/analysis/star.C	(revision 17962)
@@ -1,2 +1,3 @@
+//Star version from revision 17885 (used in the processing with the old image cleaning)
 #include "MLogManip.h"
 
@@ -303,5 +304,5 @@
 
     TString file = Form("/scratch/fact/sequences/%04d/%02d/%02d/%06d_%03d.seq",
-                        night/10000, (night/100)%100, night%100, num);
+                        night/10000, (night/100)%100, night%100, night, num);
 
     return star(file.Data(), lvl1, lvl2, inpath, outpath);
Index: branches/Mars_McMismatchStudy/fact/analysis/star.C.ORIGINAL
===================================================================
--- branches/Mars_McMismatchStudy/fact/analysis/star.C.ORIGINAL	(revision 17962)
+++ branches/Mars_McMismatchStudy/fact/analysis/star.C.ORIGINAL	(revision 17962)
@@ -0,0 +1,297 @@
+#include "MLogManip.h"
+
+int star(const char *seqfile="sequences/20111205_013.seq",
+         const char *inpath = "callisto_new", const char *outpath = "callisto_new")
+{
+    double deltat = 17.5;
+
+    // The sequence file which defines the files for the analysis
+    MSequence seq(seqfile);
+    if (!seq.IsValid())
+    {
+        gLog << err << "ERROR - Sequence invalid!" << endl;
+        return 1;
+    }
+
+    // ------------------------------------------------------
+
+    gLog.Separator("Star");
+    gLog << all << "Calculate image parameters of sequence ";
+    gLog << seq.GetFileName() << endl;
+    gLog << endl;
+
+    // ------------------------------------------------------
+
+    gLog << "Inpath:   " << inpath << endl;
+    gLog << "Outpath:  " << outpath << endl;
+
+    const TString rule(Form("s/([0-9]+_[0-9]+)_C[.]root?$/%s\\/$1_I.root/",
+                            MJob::Esc(outpath).Data()));
+    gLog << "Rule:     " << rule << endl;
+
+    MDirIter iter;
+    if (seq.GetRuns(iter, MSequence::kFactCal, inpath)<=0)
+    {
+        gLog << err << "ERROR - Sequence valid but without files." << endl;
+        return 2;
+    }
+
+    iter.Print();
+
+    gLog.Separator();
+
+    // --------------------------------------------------------
+
+    MBinning bins1(  36,    -90,    90, "BinningAlpha");
+    MBinning bins3(  67, -0.005, 0.665, "BinningTheta", "asin");
+    MBinning bins4(  25,      0,   2.5, "BinningDist");
+
+    MBinning binsM1(100,      0,     5, "BinningMuonRadius");
+    MBinning binsM2( 60,      0,   0.3, "BinningMuonDeviation");
+    MBinning binsM3(150,      0,    60, "BinningMuonTime");
+    MBinning binsM4(300,      0,    30, "BinningMuonTimeRms");
+    MBinning binsM5( 20,      0,   360, "BinningMuonArcPhi");
+    MBinning binsM6( 50,      0,   0.5, "BinningMuonArcWidth");
+    MBinning binsM7( 30,      5,  5000, "BinningMuonSize", "log");
+    MBinning binsM8(100,      0,     5, "BinningMuonRelTimeSigma");
+
+    MBinning binsM9(100,      0,     5, "BinningRadius");
+
+    // ---------------------------------------------------------
+
+    // Filter to start muon analysis
+    //MFDataPhrase fmuon1("MHillas.fSize>150 && MNewImagePar.fConcCOG<0.1", "MuonPreCut");
+    // Filter to calculate further muon parameters
+//    MFDataPhrase fmuon2("(MMuonSearchPar.fRadius*MGeomCam.fConvMm2Deg>0.6) && (MMuonSearchPar.fRadius*MGeomCam.fConvMm2Deg<1.35) &&"
+//                        "(MMuonSearchPar.fDeviation*MGeomCam.fConvMm2Deg<0.152)", "MuonSearchCut");
+    MFDataPhrase fmuon2("MMuonSearchPar.fDeviation*MGeomCam.fConvMm2Deg>0.015 &&"
+                        "MMuonSearchPar.fDeviation*MGeomCam.fConvMm2Deg>(MMuonSearchPar.fRadius*MGeomCam.fConvMm2Deg-0.1)*0.15/0.4",
+                        "MuonSearchCut");
+
+    // Filter to fill the MHMuonPar
+    MFDataPhrase fmuon3("MMuonCalibPar.fArcPhi>190 && MMuonCalibPar.fRelTimeSigma<3.2",
+                        "MuonFinalCut");
+    // Filter to write Muons to Muon tree
+    //MFDataMember fmuon4("MMuonCalibPar.fArcPhi", '>', -0.5, "MuonWriteCut");
+
+    // ---------------------------------------------------------
+
+    MStatusDisplay *d = new MStatusDisplay;
+
+    MTaskList tlist;
+
+    MParList plist2;
+    plist2.AddToList(&tlist);
+    plist2.AddToList(&bins1);
+    plist2.AddToList(&bins3);
+    plist2.AddToList(&bins4);
+    plist2.AddToList(&binsM1);
+    plist2.AddToList(&binsM2);
+    plist2.AddToList(&binsM3);
+    plist2.AddToList(&binsM4);
+    plist2.AddToList(&binsM5);
+    plist2.AddToList(&binsM6);
+    plist2.AddToList(&binsM7);
+    plist2.AddToList(&binsM8);
+    plist2.AddToList(&binsM9);
+
+    MMuonSetup muonsetup;
+    plist2.AddToList(&muonsetup);
+
+    MEvtLoop loop;
+    loop.SetDisplay(d);
+    loop.SetParList(&plist2);
+
+    MReadMarsFile read("Events");
+    read.DisableAutoScheme();
+    read.AddFiles(iter);
+
+    MContinue cont("MRawEvtHeader.GetTriggerID!=4", "SelectData");
+    //MContinue cont("(MRawEvtHeader.GetTriggerID&0xff00)!=0",     "SelectDat");
+    //MContinue cont("(MRawEvtHeader.GetTriggerID&0xff00)!=0x400", "SelectPed");
+
+    MGeomApply apply;
+
+    MImgCleanTime clean;
+    clean.SetMinCount(0);
+    clean.SetMinSize(25);
+    clean.SetDeltaT(1.25*17.5*0.1111);
+
+    MHillasCalc hcalc;
+    hcalc.Disable(MHillasCalc::kCalcConc);
+
+    MH3 hrate("MRawRunHeader.GetFileID"/*, "MRawEvtHeader.GetTriggerID"*/);
+    hrate.SetWeight("1./TMath::Max(MRawRunHeader.GetRunLength,1)");
+    hrate.SetName("Rate");
+    hrate.SetTitle("Event rate after cleaning;File Id;Event Rate [Hz];");
+    hrate.InitLabels(MH3::kLabelsX);
+    hrate.GetHist().SetMinimum(0);
+
+    MFillH fillR(&hrate, "", "FillRate");
+
+    // ------------------ Setup histograms and fill tasks ----------------
+    MHCamEvent evtC1(0, "Cleaned",  "Average signal after Cleaning;;S [phe]");
+    MHCamEvent evtC2(0, "UsedPix",  "Fraction of Events in which Pixels are used;;Fraction");
+    MHCamEvent evtC3(8, "PulsePos", "Pulse Position of signal after cleaning;;T [ns]");
+    evtC1.SetErrorSpread(kFALSE);
+    evtC2.SetErrorSpread(kFALSE);
+    evtC2.SetThreshold(0);
+
+    MFillH fillC1(&evtC1, "MSignalCam", "FillSignalCam");
+    MFillH fillC2(&evtC2, "MSignalCam", "FillCntUsedPixels");
+    MFillH fillC3(&evtC3, "MSignalCam", "FillPulsePos");
+
+    MFillH fillD1("MHHillas",      "MHillas",      "FillHillas");
+    MFillH fillD2("MHHillasExt",   "",             "FillHillasExt");
+    MFillH fillD3("MHHillasSrc",   "MHillasSrc",   "FillHillasSrc");
+    MFillH fillD4("MHImagePar",    "MImagePar",    "FillImagePar");
+    MFillH fillD5("MHNewImagePar", "MNewImagePar", "FillNewImagePar");
+    MFillH fillD6("MHVsSize",      "MHillasSrc",   "FillVsSize");
+
+    MH3 halpha("MHillasSrc.fDist*MGeomCam.fConvMm2Deg", "fabs(MHillasSrc.fAlpha)");
+    halpha.SetName("AvsD;Dist;Alpha");
+
+    MFillH filla(&halpha);
+    filla.SetNameTab("AvsD");
+    filla.SetDrawOption("colz");
+
+    // ----------------------- Muon Analysis ----------------------
+    //writem.SetFilter(&fmuon4);
+
+    MHn mhn1;
+    mhn1.AddHist("MMuonSearchPar.fRadius*MGeomCam.fConvMm2Deg");
+    mhn1.InitName("MuonRadius");
+    mhn1.InitTitle("MuonRadius");
+
+    mhn1.AddHist("MMuonSearchPar.fDeviation*MGeomCam.fConvMm2Deg");
+    mhn1.InitName("MuonDeviation");
+    mhn1.InitTitle("MuonDeviation");
+
+    mhn1.AddHist("MMuonSearchPar.fTime");
+    mhn1.InitName("MuonTime");
+    mhn1.InitTitle("MuonTime");
+
+    mhn1.AddHist("MMuonSearchPar.fTimeRms");
+    mhn1.InitName("MuonTimeRms");
+    mhn1.InitTitle("MuonTimeRms");
+
+    MHn mhn3;
+    mhn3.AddHist("MMuonSearchPar.fRadius*MGeomCam.fConvMm2Deg","MMuonSearchPar.fDeviation*MGeomCam.fConvMm2Deg");
+    mhn3.InitName("MuonRadius;MuonRadius;MuonDeviation");
+    mhn3.InitTitle("MuonRadius");
+    mhn3.SetDrawOption("colz");
+
+    MHn mhn2;
+    mhn2.AddHist("MMuonCalibPar.fArcPhi");
+    mhn2.InitName("MuonArcPhi");
+    mhn2.InitTitle("MuonArcPhi");
+
+    mhn2.AddHist("MMuonCalibPar.fArcWidth");
+    mhn2.InitName("MuonArcWidth");
+    mhn2.InitTitle("MuonArcWidth");
+
+    mhn2.AddHist("MMuonCalibPar.fMuonSize");
+    mhn2.InitName("MuonSize");
+    mhn2.InitTitle("MuonSize");
+
+    mhn2.AddHist("MMuonCalibPar.fRelTimeSigma");
+    mhn2.InitName("MuonRelTimeSigma");
+    mhn2.InitTitle("MuonRelTimeSigma");
+
+    MHn mhn4;
+    mhn4.AddHist("MMuonCalibPar.fArcPhi","MMuonCalibPar.fArcWidth");
+    mhn4.InitName("MuonArcPhi;MuonArcPhi;MuonArcWidth");
+    mhn4.InitTitle("MuonArcPhi");
+    mhn4.SetDrawOption("colz");
+
+    MFillH fillmhn1(&mhn1, "", "FillMuonSearchPar");
+    MFillH fillmhn2(&mhn2, "", "FillMuonCalibPar");
+    MFillH fillmhn3(&mhn3, "", "FillMuonDevVsRadius");
+    MFillH fillmhn4(&mhn4, "", "FillMuonWidthVsPhi");
+    fillmhn1.SetNameTab("MuonS");
+    fillmhn3.SetNameTab("DevVsRad");
+    fillmhn2.SetNameTab("MuonC");
+    fillmhn4.SetNameTab("WidthVsPhi");
+    fillmhn2.SetFilter(&fmuon2);
+    fillmhn4.SetFilter(&fmuon2);
+
+    MMuonSearchParCalc muscalc;
+    //muscalc.SetFilter(&fmuon1);
+
+    MMuonCalibParCalc mcalc;
+    mcalc.SetFilter(&fmuon2);
+
+    MFillH fillmuon("MHSingleMuon", "", "FillMuon");
+    MFillH fillmpar("MHMuonPar",    "", "FillMuonPar");
+    fillmuon.SetFilter(&fmuon2);
+    fillmpar.SetFilter(&fmuon3);
+    fillmuon.SetBit(MFillH::kDoNotDisplay);
+
+    // ---------------------------------------------------------------
+
+    MWriteRootFile write(2, rule, "RECREATE", "Image parameters"); // EffectiveOnTime
+    write.AddContainer("MTime",           "Events");
+    write.AddContainer("MHillas",         "Events");
+    write.AddContainer("MHillasExt",      "Events");
+    write.AddContainer("MHillasSrc",      "Events");
+    write.AddContainer("MImagePar",       "Events");
+    write.AddContainer("MNewImagePar",    "Events");
+    write.AddContainer("MRawEvtHeader",   "Events");
+    write.AddContainer("MRawRunHeader",   "RunHeaders");
+    write.AddContainer("MGeomCam",        "RunHeaders");
+
+    MFDataPhrase fmuonhn("MMuonCalibPar.fRelTimeSigma>=0",
+                         "MuonHistCut");
+    fillmhn2.SetFilter(&fmuonhn);
+
+    tlist.AddToList(&read);
+    tlist.AddToList(&cont);
+    tlist.AddToList(&apply);
+    tlist.AddToList(&clean);
+    tlist.AddToList(&hcalc);
+    tlist.AddToList(&fillC1);
+    tlist.AddToList(&fillC2);
+    tlist.AddToList(&fillC3);
+    tlist.AddToList(&fillR);
+    tlist.AddToList(&fillD1);
+    tlist.AddToList(&fillD2);
+    tlist.AddToList(&fillD3);
+    tlist.AddToList(&fillD4);
+    tlist.AddToList(&fillD5);
+    tlist.AddToList(&fillD6);
+    tlist.AddToList(&filla);
+    //tlist.AddToList(&fmuon1);
+    tlist.AddToList(&muscalc);
+    tlist.AddToList(&fillmhn1);
+    tlist.AddToList(&fillmhn3);
+    tlist.AddToList(&fmuon2);
+    tlist.AddToList(&fillmuon);
+    tlist.AddToList(&mcalc);
+    tlist.AddToList(&fmuonhn);
+    tlist.AddToList(&fillmhn2);
+    tlist.AddToList(&fillmhn4);
+    tlist.AddToList(&fmuon3);
+    tlist.AddToList(&fillmpar);
+    //tlist.AddToList(&fmuon4);
+    //tlist.AddToList(&writem);
+    tlist.AddToList(&write);
+
+    if (!loop.Eventloop())
+        return 3;
+
+    if (!loop.GetDisplay())
+        return 4;
+
+    TString title = "--  Image parameters #";
+    title += seq.GetSequence();
+    title += "  --";
+    d->SetTitle(title, kFALSE);
+
+    TString path;
+    path += Form("%s/20%06d_%03d-images.root", outpath,
+                 seq.GetSequence()/1000, seq.GetSequence()%1000);
+
+    d->SaveAs(path);
+
+    return 0;
+}
