Changeset 17056 for branches/Mars_MC/fact/analysis

Timestamp:

08/29/13 18:17:12 (11 years ago)

Author:

dneise

Message:

edited both scripts, so we can use them on phido for processing of simulated data

Location:

branches/Mars_MC/fact/analysis/mc

Files:

• callisto.C (modified) (22 diffs)
• star.C (modified) (8 diffs)

branches/Mars_MC/fact/analysis/mc/callisto.C

@@ -1 @@
-int callisto(const char *seqfile="seq/2012/01/23/20120123_023.seq", const char *outpath = "output", bool use_delays=true)
+
+#include <sstream>
+#include <iostream>
+
+#include "TH1F.h"
+#include "TFile.h"
+#include "TStyle.h"
+#include "TGraph.h"
+#include "TLine.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"
+
+/* 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.
+ */
+
+using namespace std;
+
+int callisto_for_monte_carlo_simulated_data(
+            const char *drs_calib_300_path="/fhgfs/groups/app/fact/mc_test/testingdrsfile/test300samples.drs.fits",
+            const char *pedestal_file_path="/fhgfs/groups/app/fact/mc_test/ceresfitstest/mcFilesForTests/mcNsbPedestal/00000001.001_P_MonteCarlo000_Events.fits",
+            const char *data_file_path="/fhgfs/groups/app/fact/mc_test/ceresfitstest/mcFilesForTests/mcProton/00000003.387_D_MonteCarlo010_Events.fits",
+            
+            const char *root_file_output_path = "/fhgfs/groups/app/callisto_star_test/callisto_for_mc_test_output/callisto.root",
+            const char *status_display_output_path = "/fhgfs/groups/app/callisto_star_test/callisto_for_mc_test_output/callisto_status_display.root", 
+            const char *status_display_title = "callisto_status_display")
 {
+    
     // ======================================================
 
@@ -10 +87 @@
 
     // map file to use (get that from La Palma!)
-    const char *map = usemap ? "/home/fact/FACT++/FACTmap111030.txt" : NULL;
+    const char *map = usemap ? "/fhgfs/groups/app/fact/resources/monte_carlo_FACTmap.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";
+        "/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 = "/cm/shared/apps/fact/Mars_svn_LP/template-pulse.root";
 
     // ------------------------------------------------------
@@ -47 +124 @@
     // Extraction type: Extract integral and half leading edge
 
-    const int type = maximum ? (MExtralgoSpline::kIntegralRel) : (MExtralgoSpline::kIntegralFixed);
+    const MExtralgoSpline::ExtractionType_t type = maximum ? (MExtralgoSpline::kIntegralRel) : (MExtralgoSpline::kIntegralFixed);
     //const int type = MExtralgoSpline::kIntegralFixed;
 
@@ -54 +131 @@
 
     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
 
     // ======================================================
@@ -65 +138 @@
     if (map && gSystem->AccessPathName(map, kFileExists))
     {
-        gLog << "ERROR - Cannot access mapping file '" << map << "'" << endl;
+        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 << "ERROR - Sequence '" << seqfile << "' invalid!" << endl;
-        return 2;
+    TString datfile = TString(data_file_path);
+    TString drsfile = TString(drs_calib_300_path);
+    TString pedfile = TString(pedestal_file_path);
+
+    gLog.Separator("Callisto");
+    gLog << all;    
+    gLog << "Data File        : " << datfile << "\n";
+    gLog << "DRS calib     300: " << drsfile << '\n';
+
+    MDrsCalibration drscalib300;
+    if (!drscalib300.ReadFits(drsfile.Data())) {
+        gLog << err << "ERROR - Cannot access drscallib300 file '" << drsfile << "'" << endl;
+        return 5;
     }
-
-    // --------------------------------------------------------------------------------
-
-    gLog.Separator("Callisto");
-    gLog << "Calibrate data of sequence '" << seq.GetFileName() << "'" << endl;
-    gLog << endl;
+    gLog << all;
+    gLog << "Pedestal     file: " << pedfile << '\n';
+
+    gLog << "root_file_output_path: " << root_file_output_path << endl;
+    gLog << "status_display_output_path: " << status_display_output_path << endl;
+    gLog << "status_display_title: " << status_display_title << 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 << "ERROR - DRS sequence invalid!" << endl;
-        return 3;
-    }
-
-    gLog << "DRS sequence file: " << drsname.str() << '\n' << endl;
-
-    TString drsfile = seq.GetFileName(0, MSequence::kRawDrs);
-    if (drsfile.IsNull())
-    {
-        cout << "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 << "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 << "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 << "ERROR - Sequence valid but without files." << endl;
-        return 7;
-    }
-    iter.Print();
-
-    // ======================================================
-
     MStatusArray arrt, arrp;
 
@@ -143 +169 @@
     if (arrt.Read()<=0)
     {
-        cout << "ERROR - Reading LP template from " << lp_template << endl;
+        gLog << err << "ERROR - Reading LP template from " << lp_template << endl;
         return 100;
     }
@@ -150 +176 @@
     if (!lpref)
     {
-        cout << "ERROR - LP Template not found in " << lp_template << endl;
+        gLog << err << "ERROR - LP Template not found in " << lp_template << endl;
         return 101;
     }
@@ -158 +184 @@
     if (!gain)
     {
-        cout << "ERROR - Gain not found in " << lp_template << endl;
+        gLog << err << "ERROR - Gain not found in " << lp_template << endl;
         return 101;
     }
@@ -166 +192 @@
     if (arrp.Read()<=0)
     {
-        cout << "ERROR - Reading Pulse template from " << pulse_template << endl;
+        gLog << err << "ERROR - Reading Pulse template from " << pulse_template << endl;
         return 102;
     }
@@ -173 +199 @@
     if (!hpulse)
     {
-        cout << "ERROR - Pulse Template not found in " << pulse_template << endl;
+        gLog << err << "ERROR - Pulse Template not found in " << pulse_template << endl;
         return 103;
     }
     hpulse->SetDirectory(0);
-
     // ======================================================
 
@@ -184 +209 @@
     MBadPixelsCam badpixels;
     badpixels.InitSize(1440);
+    /*
     badpixels[ 424].SetUnsuitable(MBadPixelsPix::kUnsuitable);
     badpixels[ 583].SetUnsuitable(MBadPixelsPix::kUnsuitable);
@@ -190 +216 @@
     badpixels[1208].SetUnsuitable(MBadPixelsPix::kUnsuitable);
     badpixels[1399].SetUnsuitable(MBadPixelsPix::kUnsuitable);
-
+    */
     //  Twin pixel
     //     113
@@ -211 +237 @@
     hrate.DefineLabelY(0x400, "Ped");
     // hrate.DefaultLabelY("ERROR");
-
-    Bool_t isinteg =
-        (type&MExtralgoSpline::kIntegral)    ||
-        (type&MExtralgoSpline::kFixedWidth)  ||
-        (type&MExtralgoSpline::kDynWidth)
-        ? kTRUE : kFALSE;
-
     gStyle->SetOptFit(kTRUE);
 
-    // ======================================================
-
-    gLog << endl;
-    gLog.Separator("Processing DRS timing calibration run");
-
-    MTaskList tlist0;
-
-    MParList plist0;
-    plist0.AddToList(&tlist0);
-    plist0.AddToList(&drscalib1024);
-    plist0.AddToList(&badpixels);
-    plist0.AddToList(&timecam);
-
-    MEvtLoop loop0("DetermineTimeCal");
-    loop0.SetDisplay(d);
-    loop0.SetParList(&plist0);
-
-    // ------------------ Setup the tasks ---------------
-
-    MRawFitsRead read0(timfile);
-
-    MContinue cont0("MRawEvtHeader.GetTriggerID!=33792", "SelectTim");
-
-    MGeomApply apply0;
-
-    MDrsCalibApply drsapply0;
-
-    MFillH fill0("MHDrsCalibrationTime");
-    fill0.SetNameTab("DeltaT");
-
-    tlist0.AddToList(&read0);
-    tlist0.AddToList(&apply0);
-    tlist0.AddToList(&drsapply0);
-    tlist0.AddToList(&cont0);
-    tlist0.AddToList(&fill0);
-
-    if (!loop0.Eventloop(max0))
-        return 8;
-
-    if (!loop0.GetDisplay())
-        return 9;
-
-    /*
-     MHDrsCalibrationT *t = (MHDrsCalibrationT*)plist4.FindObject("MHDrsCalibrationT");
-     t->SetDisplay(d);
-     t->PlotAll();
-     */
-
-    // ======================================================
-
-    gLog << endl;
-    gLog.Separator("Processing external light pulser run");
-
-    MTaskList tlist1;
-
-    MParList plist1;
-    plist1.AddToList(&tlist1);
-    plist1.AddToList(&drscalib300);
-    plist1.AddToList(&badpixels);
-    plist1.AddToList(&timecam);
-
-    MEvtLoop loop1("DetermineCalConst");
-    loop1.SetDisplay(d);
-    loop1.SetParList(&plist1);
-
-    // ------------------ Setup the tasks ---------------
-
-    MRawFitsRead read1;
-    read1.LoadMap(map);
-    read1.AddFile(calfile);
-
-    MContinue cont1("(MRawEvtHeader.GetTriggerID&0xff00)!=0x100", "SelectCal");
-
-    MGeomApply apply1;
-
-    MDrsCalibApply drsapply1;
-
-    /*
-    MPedestalCam  fPedestalCamOut1a;
-    MPedestalCam  fPedestalCamOut1b;
-    MPedestalCam  fPedestalCamOut1c;
-
-    MPedCalcPedRun pedcalc1a;
-    MPedCalcPedRun pedcalc1b;
-    MPedCalcPedRun pedcalc1c;
-    pedcalc1a.SetPedestalsOut(&fPedestalCamOut1a);
-    pedcalc1b.SetPedestalsOut(&fPedestalCamOut1b);
-    pedcalc1c.SetPedestalsOut(&fPedestalCamOut1c);
-
-    MExtractTimeAndChargeSpline extractor1ab;
-    extractor1a.SetRange(first_slice, last_slice);
-    extractor1a.SetRiseTimeHiGain(rise_time);
-    extractor1a.SetFallTimeHiGain(fall_time);
-    extractor1a.SetChargeType(type);
-    extractor1a.SetSaturationLimit(600000);
-    extractor1a.SetNoiseCalculation(kFALSE);
-
-    pedcalc1a.SetRandomCalculation(kTRUE);
-    pedcalc1b.SetRandomCalculation(kFALSE);
-    pedcalc1a.SetExtractor(&extractor1a);
-    pedcalc1b.SetExtractor($extractor1a);
-    pedcalc1c.SetRangeFromExtractor(&extractor1a);
-    */
-
-    // ---
-
-    MExtractTimeAndChargeSpline 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);
-
-    // ---
-
-    MHCamEvent evt1f(0, "ExtCalSig", "Extracted calibration signal;;S [mV·sl]");
-    MHCamEvent evt1g(4, "ExtCalTm",  "Extracted arrival times;;T [sl]");
-    MHCamEvent evt1h(6, "CalCalTm",  "Calibrated arrival times;;T [sl]");
-
-    MHSectorVsTime hist1rmsb("ExtSigVsTm");
-    MHSectorVsTime hist1tmb("CalTmVsTm");
-    hist1rmsb.SetTitle("Extracted calibration vs event number;;S [mV·sl]");
-    hist1rmsb.SetType(0);
-    hist1tmb.SetTitle("Extracted arrival time vs event number;;T [sl]");
-    //hist1tmb.SetType(4);
-    hist1tmb.SetType(6);
-
-    MFillH fill1f(&evt1f, "MExtractedSignalCam", "FillExtSig");
-    MFillH fill1g(&evt1g, "MArrivalTimeCam",     "FillExtTm");
-    MFillH fill1h(&evt1h, "MSignalCam",          "FillCalTm");
-    MFillH fill1r(&hist1rmsb, "MExtractedSignalCam", "FillExtSigVsTm");
-    //MFillH fill1j(&hist1tmb,  "MArrivalTimeCam",     "FillExtTmVsTm");
-    MFillH fill1j(&hist1tmb,  "MSignalCam",     "FillCalTmVsTm");
-
-    fill1f.SetDrawOption("gaus");
-    fill1h.SetDrawOption("gaus");
-
-    // ---
-
-    MCalibrateDrsTimes calctm1a("CalibrateCalEvents");
-    calctm1a.SetNameUncalibrated("UncalibratedTimes");
-
-    MBadPixelsTreat treat1;
-    treat1.SetProcessPedestalRun(kFALSE);
-    treat1.SetProcessPedestalEvt(kFALSE);
-
-    // ---
-
-    MHCamEvent evt1c(6, "ExtCalTmShift", "Relative extracted arrival time of calibration pulse (w.r.t. event-median);;\\Delta T [ns]");
-    MHCamEvent evt1d(6, "CalCalTmShift", "Relative calibrated arrival time of calibration pulse (w.r.t. event-median);;\\Delta T [ns]");
-
-    evt1c.SetMedianShift();
-    evt1d.SetMedianShift();
-
-    MFillH fill1c(&evt1c, "UncalibratedTimes", "FillExtCalTm");
-    MFillH fill1d(&evt1d, "MSignalCam",        "FillCalCalTm");
-    fill1d.SetDrawOption("gaus");
-
-    // ------------------ Setup eventloop and run analysis ---------------
-
-    tlist1.AddToList(&read1);
-    tlist1.AddToList(&apply1);
-    tlist1.AddToList(&drsapply1);
-    tlist1.AddToList(&cont1);
-    tlist1.AddToList(&extractor1b);
-    if (isinteg)
-    {
-        tlist1.AddToList(&extractor1c);
-        tlist1.AddToList(&fill1a);
-    }
-    tlist1.AddToList(&calctm1a);
-    tlist1.AddToList(&treat1);
-    tlist1.AddToList(&fill1f);
-    tlist1.AddToList(&fill1g);
-    tlist1.AddToList(&fill1h);
-    tlist1.AddToList(&fill1r);
-    tlist1.AddToList(&fill1j);
-    tlist1.AddToList(&fill1c);
-    tlist1.AddToList(&fill1d);
-
-    if (!loop1.Eventloop(max1))
-        return 10;
-
-    if (!loop1.GetDisplay())
-        return 11;
-
-    if (use_delays)
-        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();
+    //~ 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());
@@ -454 +264 @@
     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;
-    }
+        calib[i] =1.;
 
     gROOT->SetSelectedPad(0);
@@ -472 +278 @@
     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();
+    //~ Double_t H = w*(hpulse->GetMaximumBin()+0.5)+ax->GetXmin();
 
     TLine line;
@@ -704 +510 @@
     MRawFitsRead read5;
     read5.LoadMap(map);
-    read5.AddFiles(iter);
+    read5.AddFile(datfile);
 
     MFillH fill5a(&hrate);
@@ -777 +583 @@
     //calctm4tm.SetFilter(&filtercal);
 
-    MHCamEvent evt5m(6, "ExtTm",      "Extracted arrival times of calibration pulse;;\\Delta T [ns]");
-    MHCamEvent evt5n(6, "CalTm",      "Calibrated arrival times of calibration pulse;;\\Delta T [ns]");
-    MHCamEvent evt5q(6, "ExtTmShift", "Relative extracted arrival times of calibration pulse (w.r.t. event-median);;\\Delta T [ns]");
-    MHCamEvent evt5r(6, "CalTmShift", "Relative calibrated arrival times of calibration pulse (w.r.t. event-median);;\\Delta T [ns]");
-    MHCamEvent evt5s(6, "ExtTM",      "Extracted absolute time marker position;;T [sl]");
-    MHCamEvent evt5t(6, "CalTM",      "Calibrated absolute time marker position;;T [ns]");
-    MHCamEvent evt5u(6, "ExtTMshift", "Relative extracted time marker position (w.r.t. event-median);;\\Delta T [ns]");
-    MHCamEvent evt5v(6, "CalTMshift", "Relative calibrated time marker position (w.r.t. event-median);;\\Delta T [ns]");
-    MHCamEvent evt5w(6, "ExtDiff",    "Difference between extracted arrival time of time marker and calibration pulse;;\\Delta T [ns]");
-    MHCamEvent evt5x(6, "CalDiff",    "Difference between calibrated arrival time of time marker and calibration pulse;;\\Delta T [ns]");
-
-    evt5w.SetNameSub("UncalibratedTimes");
-    evt5x.SetNameSub("MSignalCam");
-
-    evt5q.SetMedianShift();
-    evt5r.SetMedianShift();
-    evt5u.SetMedianShift();
-    evt5v.SetMedianShift();
-    //evt4w.SetMedianShift();
-    //evt4x.SetMedianShift();
-
-    MFillH fill5m(&evt5m, "UncalibratedTimes", "FillExtTm");
-    MFillH fill5n(&evt5n, "MSignalCam",        "FillCalTm");
-    MFillH fill5q(&evt5q, "UncalibratedTimes", "FillExtTmShift");
-    MFillH fill5r(&evt5r, "MSignalCam"       , "FillCalTmShift");
-    MFillH fill5s(&evt5s, "UncalTimeMarker",   "FillExtTM");
-    MFillH fill5t(&evt5t, "TimeMarker",        "FillCalTM");
-    MFillH fill5u(&evt5u, "UncalTimeMarker",   "FillExtTMshift");
-    MFillH fill5v(&evt5v, "TimeMarker",        "FillCalTMshift");
-    MFillH fill5w(&evt5w, "UncalTimeMarker",   "FillExtDiff");
-    MFillH fill5x(&evt5x, "TimeMarker",        "FillCalDiff");
-
-    fill5m.SetDrawOption("gaus");
-    fill5n.SetDrawOption("gaus");
-    fill5q.SetDrawOption("gaus");
-    fill5r.SetDrawOption("gaus");
-    //fill5s.SetDrawOption("gaus");
-    //fill5t.SetDrawOption("gaus");
-    //fill5u.SetDrawOption("gaus");
-    //fill5v.SetDrawOption("gaus");
-    //fill5w.SetDrawOption("gaus");
-    //fill5x.SetDrawOption("gaus");
-
-
     MBadPixelsTreat treat5;
     treat5.SetProcessPedestalRun(kFALSE);
     treat5.SetProcessPedestalEvt(kFALSE);
-
-    MHSectorVsTime hist5cal("CalVsTm");
-    MHSectorVsTime hist5ped("PedVsTm");
-    hist5cal.SetTitle("Median calibrated calibration signal vs event number;;Signal [~phe]");
-    hist5ped.SetTitle("Median calibrated pedestal signal vs event number;;Signal [~phe]");
-    hist5cal.SetType(0);
-    hist5ped.SetType(0);
-    hist5cal.SetMinimum(0);
-    hist5ped.SetMinimum(0);
-    hist5cal.SetUseMedian();
-    hist5ped.SetUseMedian();
-    hist5cal.SetNameTime("MTime");
-    hist5ped.SetNameTime("MTime");
-
-    MFillH fill5cal(&hist5cal, "MSignalCam", "FillCalVsTm");
-    MFillH fill5ped(&hist5ped, "MSignalCam", "FillPedVsTm");
-    fill5cal.SetFilter(&filtercal);
-    fill5ped.SetFilter(&filterped);
 
     MHCamEvent evt5b(0, "ExtSig",   "Extracted signal;;S [mV·sl]");
@@ -868 +612 @@
     //contsw.SetInverted();
 
-    const TString fname(Form("s/([0-9]+_[0-9]+)[.]fits([.]gz)?$/%s\\/$1_C.root/",
-                             MJob::Esc(outpath).Data()));
-
     // The second rule is for the case reading raw-files!
-    MWriteRootFile write5(2, fname, "RECREATE", "Calibrated Data");
+    
+    MWriteRootFile write5(root_file_output_path, "RECREATE", "Calibrated Data", 2);
     write5.AddContainer("MRawRunHeader", "RunHeaders");
     write5.AddContainer("MGeomCam",      "RunHeaders");
@@ -879 +621 @@
     write5.AddContainer("MRawEvtHeader", "Events");
     //write.AddContainer("MTriggerPattern", "Events");
-
+    
+    
     // ------------------ Setup histograms and fill tasks ----------------
 
@@ -888 +631 @@
     tlist5tm.AddToList(&extractor5tm);
     tlist5tm.AddToList(&calctm5tm);
-    tlist5tm.AddToList(&fill5m);
-    tlist5tm.AddToList(&fill5n);
-    tlist5tm.AddToList(&fill5q);
-    tlist5tm.AddToList(&fill5r);
-    //tlist5tm.AddToList(&fill5s);
-    //tlist5tm.AddToList(&fill5t);
-    tlist5tm.AddToList(&fill5u);
-    tlist5tm.AddToList(&fill5v);
-    tlist5tm.AddToList(&fill5w);
-    tlist5tm.AddToList(&fill5x);
     tlist5tm.SetFilter(&filtercal);
 
@@ -922 +655 @@
     tlist5.AddToList(&conv5);
     tlist5.AddToList(&treat5);
-    tlist5.AddToList(&fill5ped);
-    tlist5.AddToList(&fill5cal);
     tlist5.AddToList(&tlist5dat);
     tlist5.AddToList(&write5);
@@ -933 +664 @@
         return 19;
 
-    TString title = "--  Calibrated signal #";
-    title += seq.GetSequence();
-    title += " (";
-    title += drsfile;
-    title += ")  --";
-    d->SetTitle(title, kFALSE);
-
-    TString path;
-    path += Form("%s/20%6d_%03d-calibration.root", outpath,
-                 seq.GetSequence()/1000, seq.GetSequence()%1000);
-
-    d->SaveAs(path);
+    d->SetTitle(status_display_title, kFALSE);
+    d->SaveAs(status_display_output_path);
 
     return 0;

branches/Mars_MC/fact/analysis/mc/star.C

@@ -1 @@
-int star(const char *seqfile="sequences/20111205_013.seq", Double_t lvl1=7.8, Double_t lvl2=3.9, const char *inpath = "output", const char *outpath = "output")
+#include <sstream>
+#include <iostream>
+
+#include "TH1F.h"
+#include "TFile.h"
+#include "TStyle.h"
+#include "TGraph.h"
+#include "TLine.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 "MMuonSetup.h"
+#include "MReadMarsFile.h"
+#include "MHillasCalc.h"
+#include "MHn.h"
+#include "MMuonSearchParCalc.h"
+#include "MMuonCalibParCalc.h"
+#include "MBinning.h"
+#include "MImgCleanStd.h"
+
+
+using namespace std;
+
+int star_for_monte_carlo_simulated_data(
+    const char *mars_data_file_path = "/fhgfs/groups/app/callisto_star_test/callisto_for_mc_test_output/callisto.root",
+    const char *root_output_file_path = "/fhgfs/groups/app/callisto_star_test/callisto_for_mc_test_output/star.root",
+    const char *status_display_output_path = "/fhgfs/groups/app/callisto_star_test/callisto_for_mc_test_output/star_status_display.root", 
+    const char *status_display_title = "star_status_display_title",
+    const char *hillas_txt_file_path = "/fhgfs/groups/app/callisto_star_test/star_hillas.txt",
+    Double_t lvl1=4.0, 
+    Double_t lvl2=2.5,
+    Double_t deltat = 17.5)
 {
-    double deltat = 17.5;
+    //Double_t lvl1=7.8, 
+    //Double_t lvl2=3.9,
     //    lvl1 = 2.5;
     //    lvl2 = 0.5;
-
-    // The sequence file which defines the files for the analysis
-    MSequence seq(seqfile);
-    if (!seq.IsValid())
-    {
-        gLog << "ERROR - Sequence invalid!" << endl;
-        return 1;
-    }
-
     // ------------------------------------------------------
 
     gLog.Separator("Star");
-    gLog << "Calculate image parameters of sequence ";
-    gLog << seq.GetFileName() << endl;
+    gLog << all << "Calculate image parameters of sequence ";
     gLog << endl;
+    gLog << "mars_data_file_path: " << mars_data_file_path << endl;
+    gLog << "root_output_file_path: " << root_output_file_path << endl;
+    gLog << "status_display_output_path: " << status_display_output_path << endl;
+    gLog << "status_display_title: " << status_display_title << endl;
+    gLog << endl;
+    gLog << "lvl1: " << lvl1 <<endl;
+    gLog << "lvl2: " << lvl2 <<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 << "ERROR - Sequence valid but without files." << endl;
-        return 2;
-    }
-
-    iter.Print();
 
     gLog.Separator();
@@ -45 +94 @@
     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");
@@ -54 +102 @@
     MBinning binsM7( 30,      5,  5000, "BinningMuonSize", "log");
     MBinning binsM8(100,      0,     5, "BinningMuonRelTimeSigma");
-
     MBinning binsM9(100,      0,     5, "BinningRadius");
 
@@ -104 +151 @@
     MReadMarsFile read("Events");
     read.DisableAutoScheme();
-    read.AddFiles(iter);
+    read.AddFile(mars_data_file_path);
 
     MContinue cont("MRawEvtHeader.GetTriggerID!=4", "SelectData");
@@ -228 +275 @@
     // ---------------------------------------------------------------
 
-    MWriteRootFile write(2, rule, "RECREATE", "Image parameters"); // EffectiveOnTime
+    MWriteAsciiFile write_ascii_hillas(hillas_txt_file_path, "MHillas");
+    write_ascii_hillas.AddColumns("MHillasExt");
+    write_ascii_hillas.AddColumns("MHillasSrc");
+
+
+    MWriteRootFile write(root_output_file_path, "RECREATE", "Image parameters", 2);
     write.AddContainer("MTime",           "Events");
     write.AddContainer("MHillas",         "Events");
@@ -238 +290 @@
     write.AddContainer("MRawRunHeader",   "RunHeaders");
     write.AddContainer("MGeomCam",        "RunHeaders");
-
     MFDataPhrase fmuonhn("MMuonCalibPar.fRelTimeSigma>=0",
                          "MuonHistCut");
@@ -274 +325 @@
     //tlist.AddToList(&writem);
     tlist.AddToList(&write);
+    tlist.AddToList(&write_ascii_hillas);
 
     if (!loop.Eventloop())
@@ -281 +333 @@
         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);
+    d->SetTitle(status_display_title, kFALSE);
+    d->SaveAs(status_display_output_path);
 
     return 0;