Changeset 7188 for trunk/MagicSoft

trunk/MagicSoft/Mars/Changelog

-              r7183
+              r7188
 While an underscore is a placeholder for a white-space or an empty line.
                                                  -*-*- END OF LINE -*-*-
+/07/13 Abelardo Moralejo (2005/07/12)
+   * manalysis/MGeomApply.cc
+     - Now, if the MGeomApply task has a serial number >0 (indicating
+       most likely a specific telescope in a multi-telescope file), then
+       apply the geometry only to objects with the same serial number.
+       This is important for stereo setups in which the telescopes have
+       cameras with different numbers of pixels. If the MGeomApply task
+       has serial number 0 (default), it will apply the geometry to all
+       found objects as it used to do.
+   * manalysis/MMcCalibrationUpdate.cc, mpedestal/MMcPedestalCopy.cc
+     - Now the PreProcess adds the serial number also to MRawRunHeader
+       and MMcRunHeader, since in stereo MC files produced with the
+       current camera simulation at CVS there is one such container per
+       telescope. There is no difference in behaviour for single
+       telescope files. Old MC stereo files will no longer be readable.
+   * mimage/MNewImagePar.cc
+     - Made calculation of "Inner size" valid also for camera geometries
+       other than MAGIC standard one. The calculation of "inner leakage"
+       is now made only for MAGIC standard geometry (there is no easy
+       fix for that, but having fInnerSize and MHillas.fSize, the inner
+       leakage is a completely unnecessary parameter). Added a comment on
+       a possible fix for the calculation of fConc.
+   * mcalib/MMcCalibrationCalc.[h,cc]
+     - Added member fMinSize and setter for it. It is the minimum
+       (uncalibrated) size of events to be used in the calibration.
+       Previously it was fixed at 1000 ADC counts, but the value depends
+       on the extractor used, so it must be possible to change it.
+     - Add the serial number to MRawRunHeader and MMcConfigRunHeader
+       (now there is one per telescope in stereo files).
+     - Moved the creation of histograms from the constructor to the
+       PreProcess, otherwise the serial number is not added to the
+       histogram name.
+     - In PostProcess, set the average QE for area 0 (inner pixels) in
+       MCalibrationQECam object, to be used later by MCalibrateData.
+   * macros/starmcstereo.C
+     - Big update. Adapted to all changes above. Now not only relative, but
+       also absolute calibration (=SIZE in phes) is made.
+/07/13 Raquel de los Reyes (2005/07/12)
+   * merpp.cc
+     - Added the container MCameraActiveLoad to the rootified central
+       control files
+   * mreport/MReportCC.[h,cc]
+     - Changes in the InterpreteBody to write the variable UPS status and
+       the Time difference between GPS and rubidium clock.
+  * mreport/MReportCamera.cc
+     - Changes in the InterpreteCamera and InterpreteBody functions to allow
+       the data merpping from 2005_06_30.
+/07/13 Abelardo Moralejo (2005/07/12)
+   * macros/starmcstereo.C
+     - Fixed error which made some objects to be written twice in the
+       _test_ output file.
+/07/13 Patricia Liebing (2005/07/06)
+   * mhcalib/MHCalibrationChargeCam.cc:
+    - include new flag (kLoGainBlackout) into the "uncalibrated"
+      pixel classification (see docu from 28/06/2005 Markus Gaug)
+   * mbadpixels/MBadPixelsPix.[cc.h], mcalib/MCalibrationChargeCalc.cc
+     - included flag kLoGainBlackout for unsuitable pixels w/
+       corresponding docu/printout (see MHCalibrationChargeCam.cc)
+   * mhcalib/MHCalibrationPulseTimeCam.cc, MHCalibrationCam.cc
+     - add Pixel ID to Debug printout statements
+   * mjobs/MJCalibration.cc
+     - include new display for kLoGainBlackout unsuitable pixels
+/07/13 Antonio Stamerra (2005/07/04)
+   * manalysis/MMcTriggerLvl2Calc.cc
+     - the checks written on ReInit for MC have been moved to PreProcess()
+     - kFALSE substituted with kSKIP in the check for the standard MAGIC
+       geometry
+/07/13 Abelardo Moralejo (2005/07/01)
+   * macros/starmc.C
+     - Set two different cleaning tasks for the calibration loop and
+       for the second (analysis) loop. Otherwise, some pixels were
+       cleaned out in the first loop (we do not want this) if the
+       absolute method was chosen.
+/07/13 Abelardo Moralejo (2005/06/30)
+   * mmc/MTriggerDefine.h
+     - Added TRIGGER_PIXELS_4, a provisional number of trigger pixels
+       for the MGeomCamMagic1183 geometry.
+   * mgeom/MGeomCamMagic1183.[h,cc]
+     - added. Possible design for a MAGIC-II camera. Preliminary
+       version, not yet ready to use.
+/07/13 Abelardo Moralejo (2005/06/29)
+   * mgeom/MGeomCamMagic1183.[h,cc]
+     - added. Possible design for a MAGIC-II camera
+   * mgeom/Makefile, GeomLinkDef.h
+     - added new class above.
+/07/13 Markus Gaug (2005/06/28)
+   * mhcalib/MHCalibrationChargeCam.[h,cc]
+     - set the fPickupLimit and fBlackoutLimit for the low-gain arrays
+       to a 3.5 instead of the default used in MHCalibrationPix (5.0).
+     - use a new criterium to exclude bad pixels: If the high-gain was
+       saturated and the blackout-events in the low-gain exceed the
+       fNumLoGainBlackoutLimit, the pixel is declared unsuitable.
+       This excludes those pixels which have a saturating high-gain
+       channel, but the low-gain switch does not switch often enough
+       to make the distribution reliable.
+/07/13 Markus Gaug (2005/06/23)
+   * mhcalib/MHCalibrationChargeCam.cc
+     - fix a bug counting the number of saturated events. Up to now, the
+       number of saturated slices was counted (which is one for a not too
+       high number), but for some (pathological) pixels, many more slices
+       saturated and produced wrong limits.
+     - make the PickupLimit und BlackoutLimit for the low-gain explicitely
+       smaller. The too large limits caused a failed fit for some channels
+       which resulted in wrong calibration constants. (M. Gaug)
+/07/13 Markus Gaug (2005/06/22)
+   * mjobs/MJPedestal.cc
+     - change order of reading local environmnet variables in lines 1044
+       to 1068. Reason: The lines:
+            if (fExtractor==tenv.GetTask())
+               return kTRUE;
+       yielded sometimes a return before the rest of the env-variables
+       could be read in. This affected ONLY the reading of the pedestal
+       reference file and the bad pixels file name. (M. Gaug)
+     - Added "decode" task (MTriggerPatternDecode) task to the list also
+       for MC runs.
+   * mhcalib/MHCalibrationChargeCam.cc
+     - fixed the setting of unreliability for saturated events where the
+       high-gain signal was not stable over time. Now, saturation is
+       asked for before (like in the high-gain fitted-case). (M. Gaug)
+     - fixed a small bug in the precise setting of the histogram ranges.
+       May have introduced a bias of typically 1 FADC cnt in the mean
+       fitted signal (< 0.5%).
+   * mbadpixels/MBadPixelsPix.[h,cc]
+     - took out the un-used bad pixels information in function:
+       GetUnsuitableCalLevel().
+     - debugged the class-documentation
+   * mcalib/MCalibrationIntensityChargeCam.cc
+     - a very small fix in the axis title from N_phe to N_{phe}
+     - fixed one bug in the draw-functions of this class. Not (yet) used by
+       standard calibration.
+     - +some undocumented changes
+  * mjobs/MJCalibration.cc
+    - changed order in tasklist between MCalibColorSet and
+      MPedCalcFromPedRun. Caused crashes in the intensity calibrations.
+      Does not affect ordinary calibration since MPedCalcFromPedRun is
+      then not in the tasklist anyhow.
+    - changed order of MHCalibrationChargeCam and MCalibrationRelTimeCalc.
+      This wrong order caused MCalibrationRelTimeCalc to print out
+      unreliable pixels which were not set by the rel. time calibration.
+      This caused some confusion in the reading of the output because of
+      the un-logical order of flag settings and posterior output.
+    - set the error of the hi-logain intercalibration constants
+      as the real error on the mean instead of the sigma.
+    - took out the those bad pixel-informations which are not used
+      any more in the "Defect" - display. (M. Gaug)
+/07/13 Markus Gaug (2005/06/21)
+   * mhcalib/MHCalibrationChargePix.h
+     - take out one un-necessary 'virtual'-declaration for the Draw-
+       function
+   * mhcalib/MHGausEvents.cc
+     - fixed a small bug in the non-standard display with options.
+/07/13 Markus Gaug (2005/06/19)
+   * msignal/MExtractTimeAndChargeSpline.cc
+     - set fgOffsetLoGain from 1.7 to 1.39. This is the correct time
+       shift between high-gain and low-gain. Affects only very detailed
+       timing studies.
+   * msignal/MExtractTimeAndChargeDigitalFilter.cc
+     - set fOffsetLoGain from 1.7 to 1.4 which corrects now exactly for
+       the arrival time difference between low-gain and high-gain pulses.
+       This wrong number could have had some effect on image cleanings
+       using the time information. (M. Gaug)
+   * mcalib/MCalibrationBlindCam.[h,cc]
+     - remove obsolete Copy-function
+   * mcalib/MCalibConstPix.h, mcalib/MCalibConstCam.h
+     - introduce Copy-functions
+     - set class index to 1 (was 0)
+   * mcalib/MCalibrationChargePix.[h,cc]
+     - introduce Copy() function.
+/07/13 Markus Gaug (2005/06/18)
+   * mpedestal/MPedestalCam.cc
+     - modify Print() function to print out also the event-by-event
+       averaged results.
+   * mpedestal/MExtractPedestal.[h,cc]
+     - new flag fUseSpecialPixels (set to kFALSE per default)
+     - new data member: fNameRawEvtData which can be set from outside
+       (e.g. to "MRawEvtData2").
+     - with these two changes, the special pixels can be extracted.
+       Difference to previous version: if flag fUseSpecialPixels is set,
+       * the area and sector averages are not calculated.
+       * the MPedestalCam will get initialized to
+         MRawRunHeader->GetNumSpecialPixels().
+     - fix a bug in the calculation of the event-by-event averaged
+       pedestal RMS (one total number of areas or sectors was not taken
+       into account).
+     - make local variables in the calculation of rms and mean doubles
+       instead of floats.
+   * mpedestal/MPedCalcPedRun.cc
+     - do not calculate the area and sector averages in case flag
+       fUseSpecialPixels is set.
+     - check for existence of pointer fRawEvt in PostProcess.
+/07/13 Markus Gaug (2005/06/15)
+   * mcalib/MCalibrationIntensityConstCam.[h,cc]
+     - new class to store the actual applied conversion factors (now with
+       averaged photo-electrons not in the MCalibrationIntensityChargeCam
+       any more).
+   * mjobs/MJCalibration.cc
+     - set the correct (statistical) error of the inter-calibraiton factor
+       into MCalibrationChargePix::SetConversionHiLoErr() and the
+       sigma of the distribution
+       into MCalibrationChargePix::SetConversionHiLoSigma().
+     - In MCalibrationChargePix, the mean conversion error is used for the
+       statistical error of the conversion factor.
+     - In CalibrateData, the sigma is used for the statistical error of the
+       converted signal to equiv. high-gain charges.
+/07/13 Markus Gaug (2005/06/14)
+   * mcalib/MCalibrationHiLoPix.[h,cc]
+     - store also results from profile fits.
+     - increades class version
+   * mhcalib/MHCalibrationHiLoPix.[h,cc], mhcalib/MHCalibrationHiLoCam.[h,cc]
+    - added TProfiles for every pixel filled with the high-gain signal
+      vs. low-gain signal which can be fitted to extract gain and
+      offset.
+   * mhcalib/Makefile, mhcalib/HCalibLinkDef.h
+     - added MHCalibrationHiLoPix
+/07/13 Markus Gaug (2005/06/10)
+   * mcalib/MCalibrationRelTimeCalc.cc
+     - print out forgotten successfully calibrated pixels
+     - small change in spacing of one output
+   * mhcalib/MHCalibrationCam.[h,cc]
+     - introduce max. number of events, the rest gets skipped before the
+       next ReInit() is called. Default: 4096
+       Reason: The calibration causes too many un-reliable pixels if more
+               than about 5000 events are treated (@500 Hz) because of the
+               mode hopping of the VCSels. However, in the past, some
+               calibration runs have been taken (erroneously) with more
+               than 5000 events, especially the intensity scans where
+               a good precision is needed.
+    - increased version number
+    - introduce a Getter for the fAverageAreaNum array.
+    - re-normalize the area-averaged sigmas to the equivalent number per
+      pixel also for the intensity calibration.
+   * mcalib/MCalibrationChargePix.[h,cc]
+    - divide errors in stat. and syst. errors. Both can be retrieved
+      now, depending on the study, one wants to make.
+/07/13 Markus Gaug (2005/06/06)
+   * mhcalib/MHGausEvents.[h,cc]
+     - put small functions into header file to speed up calculations a bit.
+   * mcalib/MCalibrationRelTimeCalc.cc
+     - introduce flags to steer the setting of unreliability like done
+       in MCalibrationChargeCalc. Can be steered from the conf-file.
+/07/13 Thomas Bretz
+   * mbase/MTaskInteractive.[h,cc]:
+     - fixed a bug with calling ReInit... PostProcess was called instead.
+     - added the missing initialisation of fReInit in constructor
+   * sinope.cc:
+     - removed plotting the tasklist statistics twice
+     - replaced DrawClone by DrawCopy (DrawClone gave problems)
+   * mhflux/MAlphaFitter.[h,cc]:
+     - fixed the ThetaSq function. It was wrongly defined as
+       exp(-((x-c)/s)^2) instead of exp(-1/2*((x-c)/s))
+     - Set start value for corresponding fit back to same value
+       as for a standard gauss-fit.
 /07/12 Daniela Dorner

trunk/MagicSoft/Mars/NEWS

-              r7180
+              r7188
    - callisto: MCalibrationHiLoCam can now be printed from its context
      menu, eg in the TBrowser
+   - callisto: fixed logain offset (fgOffsetLoGain) from 1.7 to
+       - 1.39 (MExtractTimeAndChargeSpline)
+       - 1.40 (MExtractTimeAndChargeDigitalFilter)
+     This is important mainly for timing studies.
+   - callisto: Changed limits in MHCalibrationChargeCalc from
+       - -100.125 to -98   (fgChargeHiGainFirst)
+       - 1899.875 to 1902. (fgChargeHiGainLast)
+       - -100.25  to -99   (fgChargeLoGainFirst)
+       - 899.75   to 901.  (fgChargeLoGainLast)
+     Introduced new limits:
+       - fgNumLoGainBlackoutLimit: 0.05
+       - fgLoGainBlackoutLimit: 3.5
+       - fgLoGainPickupLimit: 3.5
+   - callisto: use a new criterium to exclude bad pixels: If the high-gain
+     was saturated and the blackout-events in the low-gain exceed the
+     fNumLoGainBlackoutLimit, the pixel is declared unsuitable.
+     This excludes those pixels which have a saturating high-gain
+     channel, but the low-gain switch does not switch often enough
+     to make the distribution reliable.
+   - callisto: fix a bug counting the number of saturated events. Up to now,
+     the number of saturated slices was counted (which is one for a not too
+     high number), but for some (pathological) pixels, many more slices
+     saturated and produced wrong limits.
+   - callisto: New options in in callisto.rc for MCalibrationRelTimeCalc:
+     + MCalibrationRelTimeCam.CheckFitResults: Yes
+     + MCalibrationRelTimeCam.CheckDeviatingBehavior: Yes
+     + MCalibrationRelTimeCam.CheckHistOverflow: Yes
+     + MCalibrationRelTimeCam.CheckOscillations: Yes
+   - callisto: introduce max. number of events for intercalibration,
+     the rest gets skipped. Default: 4096
+     The calibration causes too many un-reliable pixels if more
+     than about 5000 events are treated (@500 Hz) because of the
+     mode hopping of the VCSels. However, in the past, some
+     calibration runs have been taken (erroneously) with more
+     than 5000 events, especially the intensity scans where
+     a good precision is needed.
    - star: fixed a bug which caused MEffectiveOnTime containers not to

trunk/MagicSoft/Mars/callisto.rc

-              r7130
+              r7188
 # Use this to change the behaviour of the calibration
 # -------------------------------------------------------------------------
 # Type if you set a colour explicitely from outside (only for MC!!!)
+# Type if you set a colour explicitely from outside (TEST purposes!!)
 #MJCalibration.MCalibColorSet.ExplicitColor: green,blue,uv,ct1
 …
 #MJCalibration.MHCalibrationChargeCam.Averageing:   yes
 #MJCalibration.MHCalibrationChargeCam.HiGainNbins:  500
 #MJCalibration.MHCalibrationChargeCam.HiGainFirst:  -100.125
 #MJCalibration.MHCalibrationChargeCam.HiGainLast:   1899.875
+#MJCalibration.MHCalibrationChargeCam.HiGainFirst:  -98.
+#MJCalibration.MHCalibrationChargeCam.HiGainLast:   1902.
 #MJCalibration.MHCalibrationChargeCam.LoGainNbins:   500
 #MJCalibration.MHCalibrationChargeCam.LoGainFirst:  -100.25
 #MJCalibration.MHCalibrationChargeCam.LoGainLast:   899.75
+#MJCalibration.MHCalibrationChargeCam.LoGainFirst:  -99.
+#MJCalibration.MHCalibrationChargeCam.LoGainLast:   901.
 #MJCalibration.MHCalibrationChargeCam.TimeLowerLimit: 1.
 #MJCalibration.MHCalibrationChargeCam.TimeUpperLimit: 3.
 …
 #MJCalibration.MHCalibrationChargeCam.OverflowLimit: 0.005
 #MJCalibration.MHCalibrationChargeCam.PulserFrequency: 500
 #MJCalibration.MHCalibrationRelTimeCam.Debug:        no
 …
 # Type of used data format: raw,root,MC
 #MJCalibrateSignal.DataType: Root
 # Type if you set a colour explicitely from outside (only for MC!!!)
+# Type if you set a colour explicitely from outside (only TESTS!!!)
 #MJCalibrateSignal.MCalibColorSet.ExpicitColor: green,blue,uv,ct1
 #MJCalibrateSignal.MCalibrateData.PedestalFlag: Event
 …
 #MJCalibrateSignal.MHCalibrationChargeCam.Averageing:   yes
 #MJCalibrateSignal.MHCalibrationChargeCam.HiGainNbins:  500
+#MJCalibrateSignal.MHCalibrationChargeCam.HiGainFirst:  -100.5
+#MJCalibrateSignal.MHCalibrationChargeCam.HiGainFirst:  -98.
+#MJCalibrateSignal.MHCalibrationChargeCam.HiGainLast:   1902.
 #MJCalibrateSignal.MHCalibrationChargeCam.HiGainLast:   1899.5
 MJCalibrateSignal.MHCalibrationChargeCam.LoGainNbins:   250
 #MJCalibrateSignal.MHCalibrationChargeCam.LoGainFirst:  -100.5
 #MJCalibrateSignal.MHCalibrationChargeCam.LoGainLast:   899.5
+#MJCalibrateSignal.MHCalibrationChargeCam.LoGainFirst:  -99.
+#MJCalibrateSignal.MHCalibrationChargeCam.LoGainLast:   901.
 #MJCalibrateSignal.MHCalibrationChargeCam.TimeLowerLimit: 1.
 #MJCalibrateSignal.MHCalibrationChargeCam.TimeUpperLimit: 3.

trunk/MagicSoft/Mars/macros/starmc.C

-              r7005
+              r7188
   // differences in gain of outer pixels)
   //
   CalibrationFilename = new TString("/data1/magic/mc_data/root/Period025/gammas_nonoise/Gamma_zbin0_*.root");
+  CalibrationFilename = new TString("nonoise/Gamma_zbin0_0_7_1000to1009_w0.root");
   // File to be used in the calibration (must be a camera file without added noise)
+  Char_t* AnalysisFilename = "Gamma_zbin1_0_*.root";  // File to be analyzed
+  Char_t* AnalysisFilename = "Gamma_zbin0_0_*.root";  // File to be analyzed
 …
   // USER CHANGE: tail cuts for image analysis
+  Float_t CleanLev[2] = {7., 5.};
+  MImgCleanStd  clean(CleanLev[0], CleanLev[1]); // Applies tail cuts to image.
+  clean.SetMethod(MImgCleanStd::kAbsolute);      // In absolute units (phot or phe- as chosen below)
+  Float_t CleanLev[2] = {10., 5.}; // Tail cuts for the analysis loop
+  MImgCleanStd clean2(CleanLev[0], CleanLev[1]); // Applies tail cuts to image.
+  clean2.SetMethod(MImgCleanStd::kAbsolute);
   //  USER CHANGE: signal extraction
-  //
-  //  MExtractFixedWindowPeakSearch sigextract;
-  //  sigextract.SetWindows(6, 6, 4);
   //  MExtractTimeAndChargeDigitalFilter sigextract;
 …
   MExtractTimeAndChargeSpline sigextract;
+  sigextract.SetChargeType(MExtractTimeAndChargeSpline::kIntegral);
   sigextract.SetRiseTimeHiGain(0.5);
   sigextract.SetFallTimeHiGain(0.5);
 …
   // USER CHANGE: high to low gain ratio. DEPENDS on EXTRACTOR!!
   // One should calculate somewhere else what this factor is for each extractor!
+  Float_t hi2low = 12.; // tentative value for spline with risetime 0.5, fall time 0.5
+  //  Float_t hi2low = 10.83;
+  // value for digital filter, HG window 4, LG window 6
+  Float_t hi2low = 11.28;
+  // value for spline with risetime 0.5, fall time 0.5, low gain stretch 1.5
 …
+  MImgCleanStd  clean(0.,0.);
+  // All pixels containing any photon will be accepted. This is what we want
+  // for the calibration loop (since no noise is added)
+  clean.SetMethod(MImgCleanStd::kAbsolute);
+  // In absolute units (phot or phe- as chosen below)
   MSrcPosCam src;
   //
+  // ONLY FOR WOBBLE MODE!!
+  //   src.SetX(120.);  // units: mm
+  // ONLY FOR WOBBLE MODE!! : set the rigt source position on camera!
+  //   src.SetX(120.);   // units: mm. Value for MC w+ files
+  //   src.SetX(-120.);  // units: mm. Value for MC w- files
   src.SetReadyToSave();
 …
   //
   MProgressBar bar;
   bar.SetWindowName("Calibrating...");
 …
   tlist.RemoveFromList(&read);
+  tlist.AddToListBefore(&clean2, &clean);
+  tlist.RemoveFromList(&clean);
   //
   // Analyzed only the desired fraction of events, skip the rest:

trunk/MagicSoft/Mars/macros/starmcstereo.C

-              r3439
+              r7188
+!
 !   Author(s):  Abelardo Moralejo 1/2004 <mailto:moralejo@pd.infn.it>
-!               Thomas Bretz, 5/2002 <mailto:tbretz@astro.uni-wuerzburg.de>
+!
 !   Copyright: MAGIC Software Development, 2000-2004
 …
 /////////////////////////////////////////////////////////////////////////////
+//
+// User change.
+//
+Float_t ctx[7] = {0., 0., 0., 0., 0., 0., 0.};
+Float_t cty[7] = {-70., -40., -30., 30., 50., 60., 70.}; // in meters
+//
+// FIXME: unfortunately present version of reflector was not prepared for
+// stereo configurations and keeps no track of CT position. So the positions
+// must be set above by the user, making sure that they correspond to the
+// files one is analysing.
+//
+void starmcstereo(Int_t ct1 = 2, Int_t ct2 = 5)
+void starmcstereo(Int_t ct1 = 1, Int_t ct2 = 2)
+{
+  if (ct1 > sizeof(ctx)/sizeof(ctx[0]) ||
+      ct2 > sizeof(ctx)/sizeof(ctx[0]) )
+    {
+      cout << endl << "Wrong CT id number!" << endl;
+  // ------------- user change -----------------
+  TString* CalibrationFilename = 0;
+  // Calibration file: a file with no added noise. Comment out next line if you
+  // do not want to calibrate the data (means SIZE will be in ADC counts)
+  CalibrationFilename = new TString("nonoise/Gamma_20_0_7_200000to200009_XX_w0.root");
+  Char_t* AnalysisFilename = "Gamma_20_0_7_*_XX_w0.root";  // File to be analyzed
+  Char_t* OutFileTag      = "gammas";           // Output file tag
+  // First open input files to check that the required telescopes
+  // are in the file, and get telescope coordinates.
+  TChain *rh = new TChain("RunHeaders");
+  rh->Add(AnalysisFilename);
+  MMcCorsikaRunHeader *corsrh = new MMcCorsikaRunHeader();
+  rh->SetBranchAddress("MMcCorsikaRunHeader.", &corsrh);
+  rh->GetEvent(0);
+  // We assume that all the read files will have the same telescopes inside,
+  // so we look only into the first runheader.
+  Int_t allcts = corsrh->GetNumCT();
+  if (ct1 >  allcts || ct2 >  allcts)
+    {
+      cout << endl << "Wrong CT id number, not contained in input file!" << endl;
       return;
+    }
+  // Set telescope coordinates as read from first runheader:
+  Float_t ctx[2];
+  Float_t cty[2];
+  ctx[0] = ((*corsrh)[ct1-1])->GetCTx();
+  cty[0] = ((*corsrh)[ct1-1])->GetCTy();
+  ctx[1] = ((*corsrh)[ct2-1])->GetCTx();
+  cty[1] = ((*corsrh)[ct2-1])->GetCTy();
+  // Now find out number of pixels in each camera:
+  MMcConfigRunHeader* confrh1 = new MMcConfigRunHeader();
+  MMcConfigRunHeader* confrh2 = new MMcConfigRunHeader();
+  rh->SetBranchAddress("MMcConfigRunHeader;1.", &confrh1);
+  rh->SetBranchAddress("MMcConfigRunHeader;2.", &confrh2);
+  rh->GetEvent(0);
+  Int_t npix[2];
+  npix[0] = confrh1->GetNumPMTs();
+  npix[1] = confrh2->GetNumPMTs();
+  rh->Delete();
   Int_t CT[2] = {ct1, ct2};  // Only 2-telescope analysis for the moment
   Int_t NCTs = sizeof(CT)/sizeof(CT[0]);
+  Int_t NCTs = 2;
   // ------------- user change -----------------
+  Char_t* AnalysisFilename = "gam-yXX-00001.root";  // File to be analyzed
+  Char_t* OutFileTag      = "gammas";           // Output file tag
+  Float_t CleanLev[2] = {4., 3.}; // Tail cuts for image analysis
+  Int_t BinsHigh[2] = {5, 9}; // First and last FADC bin of the range to be integrated,
+  Int_t BinsLow[2]  = {5, 9}; // for high and low gain respectively.
+  Float_t BinsHigh[2] = {0, 79};
+  Float_t BinsLow[2]  = {0, 79};  // FADC slices (2GHz sampling)
+  Float_t CleanLev[2] = {7., 5.}; // Units: phes (absolute cleaning will be used later!)
+  // Tail cuts for the analysis loop. In the first (calibration) loop they will not
+  // be used; we run over a noiseless file and we want to accept all pixels with any
+  // number of phes.
+  MImgCleanStd**     clean = new MImgCleanStd*[NCTs];
+  MImgCleanStd* clean[0] = new MImgCleanStd(1.,1.);
+  MImgCleanStd* clean[1] = new MImgCleanStd(1.,1.);
+  // Just dummy levels. Since the calibration file will be a noiseless file, RMS of
+  // pedestal will be 0, and all events will be accepted, regardless of the cleaning level.
+  // For some reason the above lines do not work if made on a loop! (???)
+  clean[0]->SetSerialNumber(CT[0]);
+  clean[1]->SetSerialNumber(CT[1]);
+  MExtractTimeAndChargeSpline* sigextract = new MExtractTimeAndChargeSpline[NCTs];
+  MMcCalibrationUpdate*  mccalibupdate = new MMcCalibrationUpdate[NCTs];
+  MCalibrateData* calib = new MCalibrateData[NCTs];
+  MMcCalibrationCalc* mccalibcalc = new MMcCalibrationCalc[NCTs];
   // -------------------------------------------
-  //
-  // This is a demonstration program which calculates the image
-  // parameters from a Magic Monte Carlo root file. Units of Size are
-  // for the moment, FADC counts.
-  //
-  //
   // Create a empty Parameter List and an empty Task List
   // The tasklist is identified in the eventloop by its name
   //
   MParList  plist;
   MTaskList tlist;
   plist.AddToList(&tlist);
 …
   MBadPixelsCam badpix[NCTs];
+  for (Int_t ict = 0; ict < NCTs; ict++)
+  Float_t hi2lowratio = 10.0;
+  for (Int_t i = 0; i < NCTs; i++)
+    {
       TString s = "MSrcPosCam;";
       s += CT[ict];
       src[ict].SetName(s);
       src[ict].SetReadyToSave();
       plist.AddToList(&(src[ict]));
+      s += CT[i];
+      src[i].SetName(s);
+      src[i].SetReadyToSave();
+      plist.AddToList(&(src[i]));
       TString b = "MBadPixelsCam;";
+      b += CT[ict];
+      badpix[ict].SetName(b);
+      badpix[ict].SetReadyToSave();
+      plist.AddToList(&(badpix[ict]));
+    }
+      b += CT[i];
+      badpix[i].SetName(b);
+      badpix[i].InitSize(npix[i]);
+      badpix[i].SetReadyToSave();
+      plist.AddToList(&(badpix[i]));
+      sigextract[i].SetSerialNumber(CT[i]);
+      sigextract[i].SetRange(BinsHigh[0], BinsHigh[1], BinsLow[0], BinsLow[1]);
+      sigextract[i].SetRiseTimeHiGain(0.5);
+      sigextract[i].SetFallTimeHiGain(0.5);
+      sigextract[i].SetLoGainStretch(1.);
+      mccalibupdate[i].SetSerialNumber(CT[i]);
+      mccalibupdate[i].SetUserLow2HiGainFactor(hi2lowratio);
+      mccalibupdate[i].SetSignalType(MCalibrateData::kPhe);
+      calib[i].SetSerialNumber(CT[i]);
+      calib[i].SetCalibConvMinLimit(0.);
+      calib[i].SetCalibConvMaxLimit(100.); // Override limits for real data
+      calib[i].SetCalibrationMode(MCalibrateData::kFfactor);
+      // Do not change CalibrationMode (just indicates where the cal. constants will be stored)
+      calib[i].SetSignalType(mccalibupdate[i].GetSignalType());
+      mccalibcalc[i].SetSerialNumber(CT[i]);
+      mccalibcalc[i].SetMinSize(200);
+      // Minimum SIZE for an event to be used in the calculation of calibration constants.
+      // Units are ADC counts, and value depends on signal extractor!
+    }
   //
   // Now setup the tasks and tasklist:
 …
   read.DisableAutoScheme();
+  read.AddFile(AnalysisFilename);
+  if (CalibrationFilename)
+    read.AddFile(CalibrationFilename->Data());
   MGeomApply*        apply = new MGeomApply[NCTs];
   MMcPedestalCopy*   pcopy = new MMcPedestalCopy[NCTs];
-  MExtractSignal*      sigextract = new MExtractSignal[NCTs];
-  MMcCalibrationUpdate*  mccalibupdate = new MMcCalibrationUpdate[NCTs];
-  MCalibrate* calib = new MCalibrate[NCTs];
-  MImgCleanStd**     clean = new MImgCleanStd*[NCTs];
   MHillasCalc*       hcalc = new MHillasCalc[NCTs];
-  MHillasSrcCalc*    scalc = new MHillasSrcCalc[NCTs];
   TString outfile = "star_";
   outfile += CT[0];
+  if (NCTs==2)
+    {
+      outfile += "_";
+      outfile += CT[1];
+    }
+  outfile += "_";
+  outfile += CT[1];
   //
 …
   outfile = "star_";
   outfile += CT[0];
+  if (NCTs==2)
+    {
+      outfile += "_";
+      outfile += CT[1];
+    }
+  outfile += "_";
+  outfile += CT[1];
   outfile += "_";
 …
+    {
       apply[i]->SetSerialNumber(CT[i]);
       pcopy[i]->SetSerialNumber(CT[i]);
-      sigextract[i]->SetSerialNumber(CT[i]);
-      sigextract[i].SetRange(BinsHigh[0], BinsHigh[1], BinsLow[0], BinsLow[1]);
-      mccalibupdate[i]->SetSerialNumber(CT[i]);
-      calib[i]->SetSerialNumber(CT[i]);
-      clean[i] = new MImgCleanStd(CleanLev[0], CleanLev[1]);
-      clean[i]->SetSerialNumber(CT[i]);
       hcalc[i]->SetSerialNumber(CT[i]);
+      scalc[i]->SetSerialNumber(CT[i]);
+      hcalc[i].Disable(MHillasCalc::kCalcHillasSrc);
+      // Source-dependent parameters not needed in the first loop (calibration)
       write1.SetSerialNumber(CT[i]);
       write2.SetSerialNumber(CT[i]);
+      write1.AddContainer(write1.AddSerialNumber("MMcEvt"),       "Events");
+      write1.AddContainer(write1.AddSerialNumber("MHillas"),      "Events");
+      write1.AddContainer(write1.AddSerialNumber("MHillasExt"),   "Events");
+      write1.AddContainer(write1.AddSerialNumber("MHillasSrc"),   "Events");
+      write1.AddContainer(write1.AddSerialNumber("MNewImagePar"), "Events");
+      write1.AddContainer(write1.AddSerialNumber("MSrcPosCam"),   "RunHeaders");
+      write2.AddContainer(write2.AddSerialNumber("MMcEvt"),       "Events");
+      write2.AddContainer(write2.AddSerialNumber("MHillas"),      "Events");
+      write2.AddContainer(write2.AddSerialNumber("MHillasExt"),   "Events");
+      write2.AddContainer(write2.AddSerialNumber("MHillasSrc"),   "Events");
+      write2.AddContainer(write2.AddSerialNumber("MNewImagePar"), "Events");
+      write2.AddContainer(write2.AddSerialNumber("MSrcPosCam"),   "RunHeaders");
+    }
+  if (NCTs==2)
+    {
+      write1.AddContainer("MStereoPar", "Events");
+      write2.AddContainer("MStereoPar", "Events");
+    }
+      write1.AddContainer("MMcEvt",       "Events");
+      write1.AddContainer("MHillas",      "Events");
+      write1.AddContainer("MHillasExt",   "Events");
+      write1.AddContainer("MHillasSrc",   "Events");
+      write1.AddContainer("MNewImagePar", "Events");
+      write1.AddContainer("MSrcPosCam",   "Events");
+      write2.AddContainer("MMcEvt",       "Events");
+      write2.AddContainer("MHillas",      "Events");
+      write2.AddContainer("MHillasExt",   "Events");
+      write2.AddContainer("MHillasSrc",   "Events");
+      write2.AddContainer("MNewImagePar", "Events");
+      write2.AddContainer("MSrcPosCam",   "Events");
+    }
+  MStereoPar* mstereo = new MStereoPar;
+  plist.AddToList(mstereo);
+  write1.AddContainer(mstereo, "Events");
+  write2.AddContainer(mstereo, "Events");
+  // We use MWriteRootFile::AddContainer(MParContainer* ,...) instead
+  // of using the container name as above, because in the former case the
+  // serial number tag (indicating the telescope id) is added to the
+  // container name, which is fine for containers of which there is one
+  // per telescope. However, the container MStereoPar is unique, since it
+  // is filled with information coming from both telescopes.
   write1.AddContainer("MRawRunHeader", "RunHeaders");
 …
   tlist.AddToList(&read);
+  // Skip untriggered events (now camera simulation output contains by default all simulated events)
+  MContinue* trigger = new MContinue("(MMcTrig;1.fNumFirstLevel<1) && (MMcTrig;2.fNumFirstLevel<1)","Events");
+  tlist.AddToList(trigger);
   for (i = 0; i < NCTs; i++)
 …
       tlist.AddToList(clean[i]);
       tlist.AddToList(&(hcalc[i]));
+      tlist.AddToList(&(scalc[i]));
+    }
+      tlist.AddToList(&(mccalibcalc[i]));
+    }
+  MF filter1("{MMcEvt;1.fEvtNumber%2}<0.5");
+  MF filter2("{MMcEvt;1.fEvtNumber%2}>0.5");
+  //
+  // ^^^^ Filters to divide output in two: test and train samples.
+  //
+  write1.SetFilter (&filter1);
+  write2.SetFilter (&filter2);
+  //
+  // Create and set up the eventloop
+  //
+  MProgressBar bar;
+  bar.SetWindowName("Calibrating");
+  MEvtLoop evtloop;
+  evtloop.SetProgressBar(&bar);
+  evtloop.SetParList(&plist);
+  //
+  // First loop: calibration loop. Go over MC events simulated with
+  // no noise, to correlate SIZE with the number of phes and get the
+  // conversion factor (this is done by MMcCalibrationCalc).
+  //
+  if (CalibrationFilename)
+    {
+      if (!evtloop.Eventloop())
+        return;
+    }
+  tlist.PrintStatistics();
+  ///////////////////////////////////////////////////////////////////////
+  // Now prepare the second loop: go over the events you want to analyze.
+  // This time the MMcCalibrationUpdate tasks will apply the previously
+  // calculated calibration factors.
+  // First substitute the reading task:
+  MReadMarsFile read2("Events");
+  read2.AddFile(AnalysisFilename);
+  read2.DisableAutoScheme();
+  tlist.AddToListBefore(&read2, &read);
+  tlist.RemoveFromList(&read);
+  // Delete cleaning tasks and create new ones with absolute cleaning method:
+    for (Int_t i= 0; i < NCTs; i++ )
+    {
+      tlist.RemoveFromList(clean[i]);
+      delete clean[i];
+    }
+  // New cleaning tasks:
+  clean[0] = new MImgCleanStd(CleanLev[0], CleanLev[1]);
+  clean[1] = new MImgCleanStd(CleanLev[0], CleanLev[1]);
+  clean[0]->SetMethod(MImgCleanStd::kAbsolute);
+  clean[1]->SetMethod(MImgCleanStd::kAbsolute);
+  clean[0]->SetSerialNumber(CT[0]);
+  clean[1]->SetSerialNumber(CT[1]);
+  tlist.AddToListBefore(clean[0],&(hcalc[0]));
+  tlist.AddToListBefore(clean[1],&(hcalc[1]));
+  tlist.RemoveFromList(&(mccalibcalc[0]));
+  tlist.RemoveFromList(&(mccalibcalc[1])); // Remove calibration tasks from list.
+  // Now calculate also source-dependent Hillas parameters:
+  for (i = 0; i < NCTs; i++)
+    hcalc[i].Enable(MHillasCalc::kCalcHillasSrc);
+  // Add task to calculate stereo parameters:
   MStereoCalc stereocalc;
   stereocalc.SetCTids(CT[0],CT[1]);
+  //
+  // FIXME: telescope coordinates must be introduced by the user, since
+  // they are not available nor in the camera file, neither in the reflector
+  // output.
+  //
+  stereocalc.SetCT1coor(ctx[CT[0]-1],cty[CT[0]-1]);
+  stereocalc.SetCT2coor(ctx[CT[1]-1],cty[CT[1]-1]);
+  stereocalc.SetCT1coor(ctx[0],cty[0]);
+  stereocalc.SetCT2coor(ctx[1],cty[1]);
   tlist.AddToList(&stereocalc);
+  MF filter1("{MMcEvt;1.fEvtNumber%2}<0.5");
+  MF filter2("{MMcEvt;1.fEvtNumber%2}>0.5");
+  //
+  // ^^^^ Filters to divide output in two: test and train samples.
+  //
+  write1.SetFilter (&filter1);
+  write2.SetFilter (&filter2);
+  // Add writing tasks:
   tlist.AddToList(&filter1);
   tlist.AddToList(&write1);
 …
   tlist.AddToList(&write2);
+  //
+  // Create and set up the eventloop
+  //
+  MProgressBar bar;
+  MEvtLoop evtloop;
+  evtloop.SetProgressBar(&bar);
+  evtloop.SetParList(&plist);
+  //
+  // Execute your analysis
+  //
+  bar.SetWindowName("Analyzing");
   if (!evtloop.Eventloop())
     return;
+  tlist.PrintStatistics();
   for (Int_t i= 0; i < NCTs; i++ )
     delete clean[i];
+  tlist.PrintStatistics();
+  plist.FindObject("MCalibrationChargeCam;1")->Write();
+  plist.FindObject("MCalibrationChargeCam;2")->Write();
+  plist.FindObject("MCalibrationQECam;1")->Write();
+  plist.FindObject("MCalibrationQECam;2")->Write();
+  return;
+}

trunk/MagicSoft/Mars/manalysis/MGeomApply.cc

-              r5844
+              r7188
+    {
         MCamEvent *cam = dynamic_cast<MCamEvent*>(o);
+        if (cam)
+            cam->Init(geom);
+        if (!cam)
+            continue;
+        // If the MGeomApply task has a serial number >0 (indicating most likely
+        // a specific telescope in a multi-telescope file), then apply the
+        // geometry only to objects with the same serial number. This is important
+        // for stereo setups in which the telescopes have cameras with different
+        // numbers of pixels. If the MGeomApply task has serial number 0 (default),
+        // it will apply the geometry to all found objects as it used to do.
+        const TString name(o->GetName());
+        // Extract serial number from name:
+        const Int_t serial = atoi(name.Data()+name.Last(';')+1);
+        // Compare with the serial number of this task:
+        if (serial>0 && serial!=GetSerialNumber())
+            continue;
+        // Initialize object according to camera geometry:
+        cam->Init(geom);
+    }
+}
 …
     // FIXME, workaround: this call to CalcPixRatio is here just to allow
     // the use of some camera files from the 0.7 beta version in which the
+    // the use of some MC camera files from the 0.7 beta version in which the
     // array containing pixel ratios is not initialized.
     geom->CalcPixRatio();

trunk/MagicSoft/Mars/manalysis/MMcCalibrationUpdate.cc

-              r7005
+              r7188
 Bool_t MMcCalibrationUpdate::CheckRunType(MParList *pList) const
+{
     const MRawRunHeader *run = (MRawRunHeader*)pList->FindObject("MRawRunHeader");
+    const MRawRunHeader *run = (MRawRunHeader*)pList->FindObject(AddSerialNumber("MRawRunHeader"));
     if (!run)
+    {
 …
+    }
+    MMcRunHeader* mcrunh = (MMcRunHeader*) pList->FindObject("MMcRunHeader");
+    MMcRunHeader* mcrunh = (MMcRunHeader*) pList->FindObject(AddSerialNumber("MMcRunHeader"));
+    if (!mcrunh)
+    {
+        *fLog << err << AddSerialNumber("MMcRunHeader") << " not found... aborting." << endl;
+        return kFALSE;
+    }
     //
 …
     // perpendicular to the camera plane.
     //
     // FIXME! We look for AddSerialNumber("MMcConfigRunHeader") but
     // for the moment the stereo version of camera does not write one such
     // header per telescope (it should!)
+    // As it happens with most containers, we look for AddSerialNumber("MMcConfigRunHeader")
+    // because in the stereo option the camera simulation writes one such header
+    // per telescope.
     //
     MMcConfigRunHeader* mcconfig = (MMcConfigRunHeader*) pList->FindObject(AddSerialNumber("MMcConfigRunHeader"));

trunk/MagicSoft/Mars/manalysis/MMcTriggerLvl2Calc.cc

-              r3795
+              r7188
+!
 !   Author(s): Antonio Stamerra  1/2003 <mailto:antono.stamerra@pi.infn.it>
-!   Author(s): Marcos Lopez 1/2003 <mailto:marcos@gae.ucm.es>
+!
 !   Copyright: MAGIC Software Development, 2000-2003
 …
 // run type
 //
 Bool_t MMcTriggerLvl2Calc::CheckRunType(MParList *pList) const
+Bool_t MMcTriggerLvl2Calc::IsMCRun(MParList *pList) const
+{
   const MRawRunHeader *run = (MRawRunHeader*)pList->FindObject("MRawRunHeader");
 …
 Bool_t MMcTriggerLvl2Calc::ReInit(MParList *pList)
+{
-    //
-    // If it is no MC file skip this function...
-    //
-    if (!CheckRunType(pList))
+    {
-        *fLog << inf << "This is no MC file... skipping." << endl;
-        return kTRUE;
+    }
-    //
-    // Check all necessary containers
-    //
-    fMcEvt = (MMcEvt*)pList->FindObject("MMcEvt");
-     if (!fMcEvt)
+     {
-         *fLog << err << dbginf << "MMcEvt not found... exit." << endl;
-         return kFALSE;
+     }
-    fMcTrig = (MMcTrig*)pList->FindObject("MMcTrig");
-    if (!fMcTrig)
+    {
-        *fLog << err << dbginf << "MMcTrig not found... exit." << endl;
-        return kFALSE;
+    }
-    fCam = (MGeomCam*)pList->FindObject("MGeomCam");
-    if (!fCam)
+    {
-        *fLog << dbginf << "MGeomCam not found (no geometry information available)... aborting." << endl;
-        return kFALSE;
+    }
-    // Check if fCam is a Magic geometry: only MGeomCamMagic geometry and
-    // geometries with 577 pixels are now accepted by MMcTriggerLvl2
-    if (fCam->GetNumPixels()!= 577)
+      {
-        *fLog << dbginf << "MGeomCam has a wrong geometry; only MAGIC geometry (577 pixels) is accepted by now... aborting" <<endl;
-        return kFALSE;
+      }
     return kTRUE;
+}
 …
       *fLog << "Compact pixel is set with at least "<<fMMcTriggerLvl2->GetCompactNN() << " NN" <<endl;
+    //------------------------------------------------------------
+    //
+    // If it is no MC file skip this function...
+    //
+    if (!IsMCRun(pList))
+    {
+        *fLog << inf << "Reading a data file...skipping the rest of PreProcess()" << endl;
+        return kTRUE;
+    }
+    //
+    // Check all necessary containers in case of a MC run
+    //
+    fMcEvt = (MMcEvt*)pList->FindObject("MMcEvt");
+     if (!fMcEvt)
+     {
+         *fLog << err << dbginf << "MMcEvt not found... exit." << endl;
+         return kFALSE;
+     }
+    fMcTrig = (MMcTrig*)pList->FindObject("MMcTrig");
+    if (!fMcTrig)
+    {
+        *fLog << err << dbginf << "MMcTrig not found... exit." << endl;
+        return kFALSE;
+    }
+    fCam = (MGeomCam*)pList->FindObject("MGeomCam");
+    if (!fCam)
+    {
+        *fLog << err << "MGeomCam not found (no geometry information available)... aborting." << endl;
+        return kFALSE;
+    }
+    // Check if fCam is a Magic geometry: only MGeomCamMagic geometry and
+    // geometries with 577 pixels are now accepted by MMcTriggerLvl2
+    if (fCam->GetNumPixels()!= 577)
+      {
+        *fLog << warn << "MGeomCam has a wrong geometry; only MAGIC geometry (577 pixels) is accepted by now... the Task is skipped." <<endl;
+        return kSKIP;
+      }
     return kTRUE;
+}

trunk/MagicSoft/Mars/manalysis/MMcTriggerLvl2Calc.h

r3795	r7188
33	33
34	34	Bool_t ReInit(MParList *pList);
35		Bool_t ~~CheckRunType~~(MParList *pList) const;
	35	Bool_t IsMCRun(MParList *pList) const;
36	36
37	37	public:

trunk/MagicSoft/Mars/mbadpixels/MBadPixelsPix.cc

-              r5797
+              r7188
 // are coded in the following form:
 //
+//
+// * Set bits leading to an unsuitable flag:
+//
+// BIT(7 ): kLoGainSaturation     :  The Low Gain signals were saturated during calibration
+// BIT(8 ): kChargeIsPedestal     :  The calibration signal contained only pedestals - presumably dead pixel
+// BIT(12): kMeanTimeInFirstBin   :  The signal has its mean maximum in the first used FADC slice - signal extractor bad
+// BIT(13): kMeanTimeInLast2Bins  :  The signal has its mean maximum in the last two used FADC slice - signal extractor bad
+// BIT(14): kDeviatingNumPhes     :  The calculated number of phes deviates by more than +6-5.5 sigma of the phes mean of the same area idx.
+// BIT(19): kHiGainOverFlow       :  The Hi-Gain calibration histogram showed overflow in more than 0.5% of the events
+// BIT(20): kLoGainOverFlow       :  The Lo-Gain calibration histogram showed overflow in more than 0.5% of the events
+// BIT(23): kDeadPedestalRms      :  The pedestal RMS was 4.5 sigma below or 25 sigma above the average per area idx.
+// BIT(24): kFluctuatingArivalTimes: The RMS of the position of the pulse maximum is larger than 3.5 FADC slices.
+// BIT(24): kLoGainBlackout       :  A high gain saturated pixel had too many blackout events in the low gain
+//
+//
 // * Set bits leading to an unreliable flag:
 //
 …
 // BIT(5 ): kLoGainOscillating   :  The Low  Gain signals fourier transform showed abnormal behavior
 // BIT(6 ): kRelTimeOscillating  :  The High Gain arrival times fourier transform showed abnormal behavior
 // BIT(14): kDeviatingNumPhes    :  The calculated number of photo-electrons deviates too much from the mean - inconsistency
+// BIT(11): kChargeSigmaNotValid :  The sigma of the signal distribution is smaller than the pedestal RMS - presumably a pixel with a star in its FOV only during the pedestal taking
 // BIT(16): kDeviatingFFactor    :  The calculated overall F-Factor deviates too much from the mean - inconsistency
 // BIT(15): kDeviatingNumPhots   :  The calculated number of calibrated photons deviates too much from the mean - inconsistency
 //
-// * Set bits leading to an unsuitable flag:
-//
-// BIT(7 ): kLoGainSaturation    :  The Low  Gain signals were saturated during calibration
-// BIT(8 ): kChargeIsPedestal    :  The calibration signal contained only pedestals - presumably dead pixel
-// BIT(9 ): kChargeErrNotValid   :  The absolute error of the derived charge has given non-sense - presumably pedestal
-// BIT(10): kChargeRelErrNotValid:  The relative error of the derived charge was too large or too small
-// BIT(11): kChargeSigmaNotValid :  The sigma of the pedestal distribution smaller than the pedestal RMS - presumably a pixel with a star in its FOV only during the pedestal taking
-// BIT(12): kMeanTimeInFirstBin  :  The signal has its mean maximum in the first used FADC slice - signal extractor bad
-// BIT(13): kMeanTimeInLast2Bins :  The signal has its mean maximum in the last two used FADC slice - signal extractor bad
-// BIT(19): kHiGainOverFlow      :  The Hi-Gain calibration histogram showed overflow without saturating the FADC
-// BIT(20): kLoGainOverFlow      :  The Lo-Gain calibration histogram showed overflow
 //
 // * Set bits leading to not useable low-gain signal:
 //
 // BIT(17): kConversionHiLoNotValid: The calibrated Conversion between Hi-Gain and Low Gain gives absurd results
+// BIT(17): kConversionHiLoNotValid: The inter-calibration constant between Hi-Gain and Low Gain does not exist.
 //
 // These bits can be called with the enum MBadPixelsPix::UncalibratedType_t in combination

trunk/MagicSoft/Mars/mbadpixels/MBadPixelsPix.h

-              r7103
+              r7188
       kHiLoOscillating          = BIT(22),
       kDeadPedestalRms          = BIT(23),
+      kFluctuatingArrivalTimes  = BIT(24)
+      kFluctuatingArrivalTimes  = BIT(24),
+      kLoGainBlackout           = BIT(25)
     };
 …
                                      || IsUncalibrated(kLoGainSaturation   )
                                      || IsUncalibrated(kConversionHiLoNotValid)
+                                     || IsUncalibrated(kLoGainOscillating  ) ; }
+                                     || IsUncalibrated(kLoGainOscillating  )
+                                     || IsUncalibrated(kLoGainBlackout  ); }
     Bool_t IsHiGainBad() const { return IsUnsuitable  (kUnsuitableRun      )
                                      || IsUncalibrated(kHiGainOscillating  ) ; }
     Int_t  GetUnsuitableCalLevel() const   {
+      if (!IsUnsuitable())                        return 0;
+      if (IsUncalibrated( kChargeIsPedestal    )) return 1;
+      if (IsUncalibrated( kChargeRelErrNotValid)) return 2;
+      if (IsUncalibrated( kLoGainSaturation    )) return 3;
+      if (IsUncalibrated( kMeanTimeInFirstBin  )) return 4;
+      if (IsUncalibrated( kMeanTimeInLast2Bins )) return 5;
+      if (IsUncalibrated( kDeviatingNumPhots   )) return 6;
+      if (IsUncalibrated( kHiGainOverFlow      )) return 7;
+      if (IsUncalibrated( kLoGainOverFlow      )) return 8;
+      if (IsUncalibrated( kDeadPedestalRms     )) return 9;
+      if (IsUncalibrated( kFluctuatingArrivalTimes)) return 10;
+      if (IsUncalibrated( kDeviatingNumPhes       )) return 11;
+      if (IsUncalibrated( kPreviouslyExcluded     )) return 12;
+      return 13;
+      if (!IsUnsuitable())                           return 0;
+      if (IsUncalibrated( kChargeIsPedestal    ))    return 1;
+      if (IsUncalibrated( kLoGainSaturation    ))    return 2;
+      if (IsUncalibrated( kMeanTimeInFirstBin  ))    return 3;
+      if (IsUncalibrated( kMeanTimeInLast2Bins ))    return 4;
+      if (IsUncalibrated( kHiGainOverFlow      ))    return 5;
+      if (IsUncalibrated( kLoGainOverFlow      ))    return 6;
+      if (IsUncalibrated( kDeadPedestalRms     ))    return 7;
+      if (IsUncalibrated( kFluctuatingArrivalTimes)) return 8;
+      if (IsUncalibrated( kDeviatingNumPhes    ))    return 9;
+      if (IsUncalibrated( kLoGainBlackout      ))    return 10;
+      if (IsUncalibrated( kPreviouslyExcluded  ))    return 11;
+      return 12;
+    }

trunk/MagicSoft/Mars/mbase/MParameters.h

r6915	r7188
54	54	ClassDef(MParameterI, 1) // Container to hold a generalized parameters (integer)
55	55	};
	56
56	57	/*
57	58	class MParameters : public MParContainer

trunk/MagicSoft/Mars/mcalib/CalibLinkDef.h

-              r6663
+              r7188
 #pragma link C++ class MCalibrationIntensityCam+;
 #pragma link C++ class MCalibrationIntensityChargeCam+;
+#pragma link C++ class MCalibrationIntensityConstCam+;
 #pragma link C++ class MCalibrationIntensityBlindCam+;
 #pragma link C++ class MCalibrationIntensityQECam+;
 #pragma link C++ class MCalibrationIntensityRelTimeCam+;
 #pragma link C++ class MCalibrationIntensityTestCam+;
 #pragma link C++ class MCalibrationCam+;
 #pragma link C++ class MCalibrationPix+;

trunk/MagicSoft/Mars/mcalib/MCalibConstCam.cc

-              r6457
+              r7188
 // --------------------------------------------------------------------------
 //
+// Set the size of the camera
+//
+void MCalibConstCam::InitSize(const UInt_t i)
+{
+    fArray->ExpandCreate(i);
+}
+// -------------------------------------------------------------------
+//
+// Calls TClonesArray::ExpandCreate() for:
+// - fAverageAreas
+//
+void MCalibConstCam::InitAverageAreas(const UInt_t i)
+{
+  fAverageAreas->ExpandCreate(i);
+}
+// -------------------------------------------------------------------
+//
+// Calls TClonesArray::ExpandCreate() for:
+// - fAverageSectors
+//
+void MCalibConstCam::InitAverageSectors(const UInt_t i)
+{
+  fAverageSectors->ExpandCreate(i);
+}
+// Copy 'constructor'
+//
+void MCalibConstCam::Copy(TObject &obj) const
+{
+  MParContainer::Copy(obj);
+  MCalibConstCam &cam = (MCalibConstCam&)obj;
+  Int_t n = GetSize();
+  if (n==0)
+    return;
+  cam.InitSize(n);
+  for (int i=0; i<n; i++)
+    (*this)[i].Copy(cam[i]);
+  n = GetNumAverageArea();
+  cam.InitAverageAreas(n);
+  for (int i=0; i<n; i++)
+    GetAverageArea(i).Copy(cam.GetAverageArea(i));
+  n = GetNumAverageSector();
+  cam.InitAverageSectors(n);
+  for (int i=0; i<n; i++)
+    GetAverageSector(i).Copy(cam.GetAverageSector(i));
+}
 // -------------------------------------------------------------------

trunk/MagicSoft/Mars/mcalib/MCalibConstCam.h

-              r6073
+              r7188
 #endif
+class TClonesArray;
+#ifndef ROOT_TClonesArray
+#include <TClonesArray.h>
+#endif
 class MGeomCam;
 …
   TClonesArray *fAverageSectors;  //-> Array of MCalibConstPix, one per camera sector
+  Int_t fRunNumber;               // Run number
 public:
 …
   void Clear(Option_t *o="");
+  void Copy(TObject &object) const;
   // Getters
         MCalibConstPix &GetAverageArea   ( UInt_t i );
   const MCalibConstPix &GetAverageArea   ( UInt_t i )            const;
   const Int_t           GetNumAverageArea()                      const;
         MCalibConstPix &GetAverageSector ( UInt_t i );
   const MCalibConstPix &GetAverageSector ( UInt_t i )            const;
   const Int_t           GetNumAverageSector()                    const;
   Int_t                 GetSize          ()                      const;
+        MCalibConstPix &GetAverageArea     ( UInt_t i );
+  const MCalibConstPix &GetAverageArea     ( UInt_t i )  const;
+  const Int_t           GetNumAverageArea  ()            const;
+        MCalibConstPix &GetAverageSector   ( UInt_t i );
+  const MCalibConstPix &GetAverageSector   ( UInt_t i )  const;
+  const Int_t           GetNumAverageSector()            const;
+  Int_t                 GetSize            ()            const;
         MCalibConstPix &operator[]     ( Int_t i             );
   const MCalibConstPix &operator[]     ( Int_t i             ) const;
+        MCalibConstPix &operator[]         ( Int_t i  );
+  const MCalibConstPix &operator[]         ( Int_t i  )  const;
   void  Init                           ( const MGeomCam &geom);
   void  InitSize                       ( const UInt_t i      );
   void  InitAverageAreas               ( const UInt_t i      );
   void  InitAverageSectors             ( const UInt_t i      );
+  void  InitSize                       ( const UInt_t i ) { fArray->ExpandCreate(i);          }
+  void  InitAverageAreas               ( const UInt_t i ) { fAverageAreas->ExpandCreate(i);   }
+  void  InitAverageSectors             ( const UInt_t i ) { fAverageSectors->ExpandCreate(i); }
   void Print(Option_t *o="") const;
+  // Setters
+  void SetRunNumber( const Int_t n )   {  fRunNumber = n; }
   Bool_t GetPixelContent(Double_t &val, Int_t idx, const MGeomCam &cam, Int_t type=0) const;
   void  DrawPixelContent(Int_t idx) const;
+  Bool_t GetPixelContent (Double_t &val, Int_t idx, const MGeomCam &cam, Int_t type=0) const;
+  void   DrawPixelContent(Int_t idx) const;
   ClassDef(MCalibConstCam, 0)   // Temporary Storage for calibration constants
+  ClassDef(MCalibConstCam, 1)   // Temporary Storage for calibration constants
 };

trunk/MagicSoft/Mars/mcalib/MCalibConstPix.h

-              r6073
+              r7188
+{
 private:
     Float_t fCalibConst;        // conversion factor (modified after each interlaced cal. update)
     Float_t fCalibFFactor;      // global F-Factor   (modified after each interlaced cal. update)
 public:
     MCalibConstPix();
+    // TObject
     void Clear(Option_t *o="")
+      {
+        fCalibConst   = -1.;
+        fCalibFFactor = -1.;
+      }
+    {
+        fCalibConst   = -1.;
+        fCalibFFactor = -1.;
+    }
+    void Copy(TObject &object) const
+    {
+        MCalibConstPix &pix =  (MCalibConstPix&)object;
+        pix.fCalibConst   = fCalibConst;
+        pix.fCalibFFactor = fCalibFFactor;
+    }
     // Getters
     Float_t GetCalibConst()   const { return fCalibConst;   }
 …
     void SetCalibFFactor( const Float_t f )  { fCalibFFactor = f; }
     ClassDef(MCalibConstPix, 0) // Temporay Storage Calibraion Constant of one pixel
+    ClassDef(MCalibConstPix, 1) // Temporay Storage Calibraion Constant of one pixel
 };

trunk/MagicSoft/Mars/mcalib/MCalibrationBlindCam.cc

-              r5128
+              r7188
+}
-// --------------------------------------------------------------------------
-//
-// Copy 'constructor'
-//
-void MCalibrationBlindCam::Copy(TObject& object) const
+{
-  MCalibrationBlindCam &calib = (MCalibrationBlindCam&)object;
-  MCalibrationCam::Copy(calib);
-  /*
-  const UInt_t n = GetSize();
-  if (n != 0)
+    {
-      calib.InitSize(n);
-      for (UInt_t i=0; i<n; i++)
-        (*this)[i].Copy(calib[i]);
+    }
-  */
+}
 // --------------------------------------------------------------------------

trunk/MagicSoft/Mars/mcalib/MCalibrationBlindCam.h

-              r5128
+              r7188
 public:
   MCalibrationBlindCam(Int_t nblind=1,const char *name=NULL, const char *title=NULL);
-  void Copy(TObject& object) const;
   // Inits
   void  Init ( const MGeomCam &geom ) {}
+  // Getter
   Bool_t IsFluxInsidePlexiglassAvailable () const;
 …
   Float_t GetFluxInsidePlexiglassRelVar() const;
+  // Setter
   void  SetPulserColor( const MCalibrationCam::PulserColor_t col );

trunk/MagicSoft/Mars/mcalib/MCalibrationChargePix.cc

-              r7099
+              r7188
 //  + fNumSaturated
 //
+// Class Version 4:
+//  +  Float_t fConversionHiLoSigma;             // Sigma of conversion factor betw. Hi and Lo Gain
+//  -  Float_t fMeanConvFADC2PheVar;             // Variance conversion factor (F-factor method)
+//  +  Float_t fMeanConvFADC2PheStatVar;         // Variance conversion factor, only stat. error
+//  +  Float_t fMeanConvFADC2PheSystVar;         // Variance conversion factor, only syst. error
+//  -  Float_t fMeanFFactorFADC2PhotVar;         // Variance mean F-Factor photons (F-factor method)
+//  +  Float_t fMeanFFactorFADC2PhotVar;         // Variance mean F-Factor photons, only stat. error
+//  -  Float_t fPheFFactorMethod;                // Number Phe's calculated (F-factor method)
+//  -  Float_t fPheFFactorMethodVar;             // Variance number of Phe's (F-factor method)
+//  +  Float_t fPheFFactorMethod;                // Number Phe's calculated  with F-factor method)
+//  +  Float_t fPheFFactorMethodStatVar;         // Variance number of Phe's, only stat. error
+//  +  Float_t fPheFFactorMethodSystVar;         // Variance number of Phe's, only syst. error
+//
 /////////////////////////////////////////////////////////////////////////////
 #include "MCalibrationChargePix.h"
 …
 const Float_t MCalibrationChargePix::gkFFactorErr          = 0.02;
+const Float_t MCalibrationChargePix::fgConversionHiLo           = 10.;
+const Float_t MCalibrationChargePix::fgConversionHiLoErr        = 2.5;
+const Float_t MCalibrationChargePix::fgPheFFactorMethodLimit    = 1.;
+const Float_t MCalibrationChargePix::fgConvFFactorRelErrLimit   = 0.85;
+const Float_t MCalibrationChargePix::fgConversionHiLo         = 10.;
+const Float_t MCalibrationChargePix::fgConversionHiLoErr      = 0.05;
+const Float_t MCalibrationChargePix::fgConversionHiLoSigma    = 2.5;
+const Float_t MCalibrationChargePix::fgPheFFactorMethodLimit  = 1.;
+const Float_t MCalibrationChargePix::fgConvFFactorRelErrLimit = 0.85;
 // --------------------------------------------------------------------------
 …
   fPheFFactorMethod                 =  -1.;
+  fPheFFactorMethodVar              =  -1.;
+  fPheFFactorMethodStatVar          =  -1.;
+  fPheFFactorMethodSystVar          =  -1.;
   fMeanConvFADC2Phe                 =  -1.;
+  fMeanConvFADC2PheVar              =  -1.;
+  fMeanConvFADC2PheStatVar          =  -1.;
+  fMeanConvFADC2PheSystVar          =  -1.;
   fMeanFFactorFADC2Phot             =  -1.;
   fMeanFFactorFADC2PhotVar          =  -1.;
 …
   MCalibrationPix::Clear();
+}
+// -----------------------------------------------------
+//
+void MCalibrationChargePix::Copy(TObject& object) const
+{
+  MCalibrationChargePix &pix = (MCalibrationChargePix&)object;
+  MCalibrationPix::Copy(pix);
+  //
+  // Copy the data members
+  //
+  pix.fAbsTimeMean                = fAbsTimeMean              ;
+  pix.fAbsTimeRms                 = fAbsTimeRms               ;
+  pix.fCalibFlags                 = fCalibFlags               ;
+  pix.fConversionHiLo             = fConversionHiLo           ;
+  pix.fConversionHiLoVar          = fConversionHiLoVar        ;
+  pix.fConversionHiLoSigma        = fConversionHiLoSigma      ;
+  pix.fConvFFactorRelVarLimit     = fConvFFactorRelVarLimit   ;
+  pix.fLoGainPedRmsSquare         = fLoGainPedRmsSquare       ;
+  pix.fLoGainPedRmsSquareVar      = fLoGainPedRmsSquareVar    ;
+  pix.fMeanConvFADC2Phe           = fMeanConvFADC2Phe         ;
+  pix.fMeanConvFADC2PheStatVar    = fMeanConvFADC2PheStatVar  ;
+  pix.fMeanConvFADC2PheSystVar    = fMeanConvFADC2PheSystVar  ;
+  pix.fMeanFFactorFADC2Phot       = fMeanFFactorFADC2Phot     ;
+  pix.fMeanFFactorFADC2PhotVar    = fMeanFFactorFADC2PhotVar  ;
+  pix.fPed                        = fPed                      ;
+  pix.fPedVar                     = fPedVar                   ;
+  pix.fPedRms                     = fPedRms                     ;
+  pix.fPedRmsVar                  = fPedRmsVar                  ;
+  pix.fPheFFactorMethod           = fPheFFactorMethod           ;
+  pix.fPheFFactorMethodStatVar    = fPheFFactorMethodStatVar    ;
+  pix.fPheFFactorMethodSystVar    = fPheFFactorMethodSystVar    ;
+  pix.fPheFFactorMethodLimit      = fPheFFactorMethodLimit      ;
+  pix.fRSigmaSquare               = fRSigmaSquare               ;
+  pix.fRSigmaSquareVar            = fRSigmaSquareVar            ;
+  pix.fNumSaturated               = fNumSaturated               ;
+}
 …
 void MCalibrationChargePix::SetPedestal(const Float_t ped, const Float_t pedrms, const Float_t pederr)
+{
   fPed       = ped;
   fPedRms    = pedrms;
 …
 void MCalibrationChargePix::SetPed(const Float_t ped, const Float_t pederr)
+{
   fPed       = ped;
   fPedVar    = pederr*pederr;
 …
 void MCalibrationChargePix::SetPedRMS( const Float_t pedrms, const Float_t pedrmserr)
+{
   fPedRms    = pedrms;
   fPedRmsVar = pedrmserr*pedrmserr;
+}
+// -------------------------------------------------------------------------------
+//
+// Get the conversion Error Hi-Gain to Low-Gain:
+// - If fConversionHiLoVar is smaller than 0 (i.e. has not yet been set), return -1.
+}
+// Inline Functions:
+// -----------------
+//
+// GetConversionHiLoErr:
+//   Get the conversion Error Hi-Gain to Low-Gain:
+//   If fConversionHiLoVar is smaller than 0 (i.e. has not yet
+//   been set), return -1.
 //
+Float_t MCalibrationChargePix::GetConversionHiLoErr()  const
+{
+  if (fConversionHiLoVar < 0.)
+    return -1.;
+  return TMath::Sqrt(fConversionHiLoVar);
+}
+// GetPedRms(): Get the pedestals RMS: Test bit kHiGainSaturation:
+//  If yes, return square root of fLoGainPedRmsSquare (if greater than 0,
+//  otherwise -1.), If no,  return fPedRms
+//
+// GetConvertedMean(): Get the Low Gain Mean Charge converted to High Gain
+//  amplification: Returns fLoGainMean multiplied with fConversionHiLo if
+//  IsHiGainSaturation(), else return fHiGainMean
+//
+// GetConvertedMeanErr(): Get the Error of the converted Low Gain Mean:
+//  Returns -1 if the variable fLoGainMean or fLoGainMeanVar are smaller than 0.
+//  Returns the square root of the quadratic sum of the relative variances of
+//  the fLoGainMean and fConversionHiLo, mulitplied with GetConvertedMean()
+//  in case of HiGain Saturation,
+//  else return GetMeanErr()
+//
+// GetConvertedSigma(): Get the Low Gain Sigma converted to High Gain
+//  amplification: Returns fLoGainSigma multiplied with fConversionHiLo
+//  if IsHiGainSaturation() else return fHiGainSigma
+//
+// GetConvertedSigmaErr(): Get the Error of the converted Sigma:
+//  Returns -1 if the variable fLoGainSigma or fLoGainSigmaVar are smaller than 0.
+//  if IsHiGainSaturatio()
+//  returns the square root of the quadratic sum of the relative variances of
+//  the fLoGainSigma and fConversionHiLo, mulitplied with GetConvertedSigma()
+//  else returns GetSigmaErr()
+//
+// GetConvertedRSigma(): Get the converted reduced Sigma:
+//  If fRSigmaSquare is smaller than 0 (i.e. has not yet been set), return -1.
+//  Test bit kHiGainSaturation:
+//  If yes, return square root of fRSigmaSquare, multiplied with fConversionHiLo,
+//  If no , return square root of fRSigmaSquare
+//
+// GetRSigma(): Get the reduced Sigma:
+//  If fRSigmaSquare is smaller than 0 (i.e. has not yet been set), return -1.
+//
+// GetConvertedRSigmaSquare(): Get the reduced Sigma Square:
+//  If fRSigmaSquare is smaller than 0 (i.e. has not yet been set), return -1.
+//  Test bit kHiGainSaturation:
+//  If yes, return fRSigmaSquare, multiplied with fConversionHiLo^2,
+//  If no , return fRSigmaSquare
+//
+// GetPheFFactorMethodErr(): Get the error on the number of photo-electrons
+//  (F-Factor Method): If fPheFFactorMethodVar is smaller than 0 (i.e. has
+//  not yet been set), return -1. Else returns the square root of
+//  fPheFFactorMethodVar
+//
+// GetMeanFFactorFADC2PhotErr(): Get the error on the mean total F-Factor
+//  of the signal readout (F-Factor Method): If fMeanFFactorFADC2PhotVar
+//  is smaller than 0 (i.e. has not yet been set), return -1. Else returns
+//  the square root of fMeanFFactorFADC2PhotVar
+//
+// GetPheFFactorMethodRelVar(): Get the relative variance on the number of
+//  photo-electrons (F-Factor Method): If fPheFFactorMethodVar is smaller
+//  than 0 (i.e. has not yet been set), return -1. If fPheFFactorMethod
+//  is 0, return -1. Else returns fPheFFactorMethodVar / fPheFFactorMethod^2
+//
+// GetMeanConvFADC2PheErr(): Get the error on the mean conversion factor
+//  (FFactor  Method): If fMeanConvFADC2PheVar is smaller than 0 (i.e. has
+//  not yet been set), return -1. Else returns the square root of
+//  fMeanConvFADC2PheVar
 // --------------------------------------------------------------------------
 …
-// --------------------------------------------------------------------------
-//
-// Get the pedestals RMS:
-// - Test bit kHiGainSaturation:
-//   If yes, return square root of fLoGainPedRmsSquare (if greater than 0, otherwise -1.),
-//   If no,  return fPedRms
-//
-Float_t  MCalibrationChargePix::GetPedRms()  const
+{
-  if (IsHiGainSaturation())
-    if (fLoGainPedRmsSquare < 0.)
-      return -1.;
-    else
-      return TMath::Sqrt(fLoGainPedRmsSquare);
-  return fPedRms;
+}
 // --------------------------------------------------------------------------
 …
 // --------------------------------------------------------------------------
 //
-// Get the Low Gain Mean Charge converted to High Gain amplification:
-// Returns fLoGainMean multiplied with fConversionHiLo if IsHiGainSaturation(),
-//         else return fHiGainMean
-//
-Float_t MCalibrationChargePix::GetConvertedMean()  const
+{
-  if (IsHiGainSaturation())
-    return fLoGainMean * fConversionHiLo;
-  return fHiGainMean;
+}
-// --------------------------------------------------------------------------
-//
-// Get the Error of the converted Low Gain Mean:
-//
-// Returns -1 if the variable fLoGainMean or fLoGainMeanVar are smaller than 0.
-//
-// Returns the square root of the quadratic sum of the relative variances of
-// the fLoGainMean and fConversionHiLo, mulitplied with GetConvertedMean()
-// in case of HiGain Saturation,
-// else return GetMeanErr()
-//
-Float_t MCalibrationChargePix::GetConvertedMeanErr()  const
+{
-  if (IsHiGainSaturation())
+    {
-      const Float_t logainrelvar = GetLoGainMeanRelVar();
-      if (logainrelvar < 0.)
-        return -1.;
-      return TMath::Sqrt(logainrelvar + GetConversionHiLoRelVar()) * GetConvertedMean();
+    }
-  else
-    return GetMeanErr();
+}
-// --------------------------------------------------------------------------
-//
-// Get the Low Gain Sigma converted to High Gain amplification:
-// Returns fLoGainSigma multiplied with fConversionHiLo if IsHiGainSaturation()
-// else return fHiGainSigma
-//
-Float_t MCalibrationChargePix::GetConvertedSigma()  const
+{
-  if (IsHiGainSaturation())
-    return fLoGainSigma * fConversionHiLo;
-  else
-    return fHiGainSigma;
+}
-// --------------------------------------------------------------------------
-//
-// Get the Error of the converted Sigma:
-//
-// Returns -1 if the variable fLoGainSigma or fLoGainSigmaVar are smaller than 0.
-//
-// if IsHiGainSaturatio()
-// returns the square root of the quadratic sum of the relative variances of
-// the fLoGainSigma and fConversionHiLo, mulitplied with GetConvertedSigma()
-// else returns GetSigmaErr()
-//
-Float_t MCalibrationChargePix::GetConvertedSigmaErr()  const
+{
-  if (IsHiGainSaturation())
+    {
-      if (fLoGainSigmaVar < 0.)
-        return -1.;
-      if (fLoGainSigma < 0.)
-        return -1.;
-      const Float_t sigmaRelVar =  fLoGainSigmaVar
-                                /( fLoGainSigma * fLoGainSigma );
-      return TMath::Sqrt(sigmaRelVar+GetConversionHiLoRelVar()) * GetConvertedSigma();
+    }
-  else
-    return GetSigmaErr();
+}
-// --------------------------------------------------------------------------
-//
-// Get the converted reduced Sigma:
-// - If fRSigmaSquare is smaller than 0 (i.e. has not yet been set), return -1.
-// - Test bit kHiGainSaturation:
-//   If yes, return square root of fRSigmaSquare, multiplied with fConversionHiLo,
-//   If no , return square root of fRSigmaSquare
-//
-Float_t MCalibrationChargePix::GetConvertedRSigma()  const
+{
-  if (fRSigmaSquare < 0)
-    return -1;
-  const Float_t rsigma = TMath::Sqrt(fRSigmaSquare);
-  return IsHiGainSaturation() ? rsigma*fConversionHiLo : rsigma ;
+}
-// --------------------------------------------------------------------------
-//
 // Get the error of the converted reduced Sigma:
 // - If fRSigmaSquareVar is smaller than 0 (i.e. has not yet been set), return -1.
 …
 // --------------------------------------------------------------------------
 //
-// Get the reduced Sigma:
-// - If fRSigmaSquare is smaller than 0 (i.e. has not yet been set), return -1.
-//
-Float_t MCalibrationChargePix::GetRSigma()  const
+{
-  if (fRSigmaSquare < 0)
-    return -1;
-  return TMath::Sqrt(fRSigmaSquare);
+}
-// --------------------------------------------------------------------------
-//
 // Get the error of the reduced Sigma:
 // - If fRSigmaSquareVar is smaller than 0 (i.e. has not yet been set), return -1.
 …
   return TMath::Sqrt(rsigmarelvar + meanrelvar) * GetRSigmaPerCharge();
+}
-// --------------------------------------------------------------------------
-//
-// Get the reduced Sigma Square:
-// - If fRSigmaSquare is smaller than 0 (i.e. has not yet been set), return -1.
-// - Test bit kHiGainSaturation:
-//   If yes, return fRSigmaSquare, multiplied with fConversionHiLo^2,
-//   If no , return fRSigmaSquare
-//
-Float_t MCalibrationChargePix::GetConvertedRSigmaSquare()  const
+{
-  if (fRSigmaSquare < 0)
-    return -1;
-  return IsHiGainSaturation() ? fRSigmaSquare*fConversionHiLo*fConversionHiLo : fRSigmaSquare ;
+}
 …
+}
+// --------------------------------------------------------------------------
+//
+// Get the error on the number of photo-electrons (F-Factor Method):
+// - If fPheFFactorMethodVar is smaller than 0 (i.e. has not yet been set), return -1.
+// - Else returns the square root of fPheFFactorMethodVar
+//
+Float_t MCalibrationChargePix::GetPheFFactorMethodErr()  const
+{
+  if (fPheFFactorMethodVar < 0.)
+    return -1.;
+  return TMath::Sqrt(fPheFFactorMethodVar);
+}
+// --------------------------------------------------------------------------
+//
+// Get the error on the mean total F-Factor of the signal readout (F-Factor Method):
+// - If fMeanFFactorFADC2PhotVar is smaller than 0 (i.e. has not yet been set), return -1.
+// - Else returns the square root of fMeanFFactorFADC2PhotVar
+//
+Float_t MCalibrationChargePix::GetMeanFFactorFADC2PhotErr()  const
+{
+  if (fMeanFFactorFADC2PhotVar < 0.)
+    return -1.;
+  return TMath::Sqrt(fMeanFFactorFADC2PhotVar);
+}
+// --------------------------------------------------------------------------
+//
+// Get the relative variance on the number of photo-electrons (F-Factor Method):
+// - If fPheFFactorMethodVar is smaller than 0 (i.e. has not yet been set), return -1.
+// - If fPheFFactorMethod    is 0, return -1.
+// - Else returns fPheFFactorMethodVar / fPheFFactorMethod^2
+//
+Float_t MCalibrationChargePix::GetPheFFactorMethodRelVar()  const
+{
+  if (fPheFFactorMethodVar < 0.)
+    return -1.;
+  if (fPheFFactorMethod  == 0.)
+    return -1.;
+  return fPheFFactorMethodVar / (fPheFFactorMethod * fPheFFactorMethod);
+}
+// --------------------------------------------------------------------------
+//
+// Get the error on the mean conversion factor (FFactor  Method):
+// - If fMeanConvFADC2PheVar is smaller than 0 (i.e. has not yet been set), return -1.
+// - Else returns the square root of fMeanConvFADC2PheVar
+//
+Float_t MCalibrationChargePix::GetMeanConvFADC2PheErr()  const
+{
+  if (fMeanConvFADC2PheVar < 0.)
+    return -1.;
+  return TMath::Sqrt(fMeanConvFADC2PheVar);
+}
 // --------------------------------------------------------------------------
 …
   // Calculate the Error of Nphe
   //
+  const Float_t pheRelVar = ffactorsquareRelVar + meanSquareRelVar + rsigmaSquareRelVar;
+  fPheFFactorMethodVar    = pheRelVar * fPheFFactorMethod * fPheFFactorMethod;
+  const Float_t pheRelVar  = meanSquareRelVar + rsigmaSquareRelVar;
+  fPheFFactorMethodStatVar = pheRelVar * fPheFFactorMethod * fPheFFactorMethod;
+  fPheFFactorMethodSystVar = ffactorsquareRelVar * fPheFFactorMethod * fPheFFactorMethod;
   if (IsDebug())
 …
+    }
   if (fPheFFactorMethodVar < 0. )
+  if (fPheFFactorMethodStatVar < 0. )
     return kFALSE;
 …
   // the errors, but have to take account of this cancellation:
   //
+  Float_t convrelvar = ffactorsquareRelVar + GetMeanRelVar() + rsigmaSquareRelVar;
+  Float_t convrelvar = GetMeanRelVar() + rsigmaSquareRelVar;
+  if (IsHiGainSaturation())
+      convrelvar += GetConversionHiLoRelVar();
   const Float_t limit = IsHiGainSaturation() ? fConvFFactorRelVarLimit * 4. : fConvFFactorRelVarLimit;
 …
+    }
+  fMeanConvFADC2PheVar =  convrelvar * fMeanConvFADC2Phe * fMeanConvFADC2Phe;
+  fMeanConvFADC2PheStatVar = convrelvar * fMeanConvFADC2Phe  * fMeanConvFADC2Phe;
+  fMeanConvFADC2PheSystVar = ffactorsquareRelVar * fMeanConvFADC2Phe * fMeanConvFADC2Phe;
   SetFFactorMethodValid(kTRUE);
 …
                               + 0.25 * nphotonsrelvar;
   fMeanFFactorFADC2PhotVar    = ffactorrelvar * fMeanFFactorFADC2Phot * fMeanFFactorFADC2Phot;
+  fMeanFFactorFADC2PhotVar = ffactorrelvar * fMeanFFactorFADC2Phot * fMeanFFactorFADC2Phot;
   if (IsDebug())

trunk/MagicSoft/Mars/mcalib/MCalibrationChargePix.h

-              r7099
+              r7188
   static const Float_t fgConversionHiLo;         //! Default fConversionHiLo          (now set to: 10.)
   static const Float_t fgConversionHiLoErr;      //! Default fConversionHiLoVar       (now set to: 2.5)
+  static const Float_t fgConversionHiLoSigma;    //! Default fConversionHiLoSigma     (now set to: 2.5)
   static const Float_t fgPheFFactorMethodLimit;  //! Default fPheFFactorMethodLimit   (now set to: 5.)
   static const Float_t fgConvFFactorRelErrLimit; //! Default fConvFFactorRelErrLimit  (now set to: 0.35)
 …
   Float_t fConversionHiLo;                  // Conversion factor betw. Hi Gain and Lo Gain
   Float_t fConversionHiLoVar;               // Variance Conversion factor betw. Hi and Lo Gain
+  Float_t fConversionHiLoSigma;             // Sigma of conversion factor betw. Hi and Lo Gain
   Float_t fConvFFactorRelVarLimit;          // Limit for acceptance rel. variance Conversion FADC2Phe
   Float_t fLoGainPedRmsSquare;              // Pedestal RMS square of Low Gain
   Float_t fLoGainPedRmsSquareVar;           // Pedestal RMS square Variance of Low Gain
   Float_t fMeanConvFADC2Phe;                // Conversion factor (F-factor method)
+  Float_t fMeanConvFADC2PheVar;             // Variance conversion factor (F-factor method)
+  Float_t fMeanConvFADC2PheStatVar;         // Variance conversion factor, only stat. error
+  Float_t fMeanConvFADC2PheSystVar;         // Variance conversion factor, only syst. error
   Float_t fMeanFFactorFADC2Phot;            // Total mean F-Factor to photons (F-factor method)
   Float_t fMeanFFactorFADC2PhotVar;         // Variance mean F-Factor photons (F-factor method)
+  Float_t fMeanFFactorFADC2PhotVar;         // Variance mean F-Factor photons, only stat. error
   Float_t fPed;                             // Pedestal (from MPedestalPix) times number FADC slices
   Float_t fPedVar;                          // Variance of pedestal
   Float_t fPedRms;                          // Pedestal RMS (from MPedestalPix) times sqrt nr. FADC slices
   Float_t fPedRmsVar;                       // Pedestal RMS (from MPedestalPix) times sqrt nr. FADC slices
+  Float_t fPheFFactorMethod;                // Number Phe's calculated (F-factor method)
+  Float_t fPheFFactorMethodVar;             // Variance number of Phe's (F-factor method)
+  Float_t fPheFFactorMethod;                // Number Phe's calculated  with F-factor method)
+  Float_t fPheFFactorMethodStatVar;         // Variance number of Phe's, only stat. error
+  Float_t fPheFFactorMethodSystVar;         // Variance number of Phe's, only syst. error
   Float_t fPheFFactorMethodLimit;           // Min. number Photo-electrons for pix to be accepted.
   Float_t fRSigmaSquare;                    // Square of Reduced sigma
 …
 public:
+  MCalibrationChargePix(const char *name=NULL, const char *title=NULL);
+  MCalibrationChargePix(const char *name=NULL, const char *title=NULL);
+  // TObject
   void Clear(Option_t *o="");
+  void Copy(TObject& object) const;
   // Setter
+  void SetAbsTimeMean ( const Float_t f ) { fAbsTimeMean = f; }
+  void SetAbsTimeRms  ( const Float_t f ) { fAbsTimeRms  = f; }
+  void SetConversionHiLo    ( const Float_t c=fgConversionHiLo    )        { fConversionHiLo    = c;       }
+  void SetConversionHiLoErr ( const Float_t e=fgConversionHiLoErr )        { fConversionHiLoVar = e*e;     }
+  /*
   void SetConvFFactorRelErrLimit   ( const Float_t f=fgConvFFactorRelErrLimit) { fConvFFactorRelVarLimit = f*f;}
-  void SetFFactorMethodValid   ( const Bool_t b = kTRUE );
   void SetMeanConvFADC2Phe      ( const Float_t f)                          { fMeanConvFADC2Phe       = f; }
   void SetMeanConvFADC2PheVar   ( const Float_t f)                          { fMeanConvFADC2PheVar    = f; }
   void SetMeanFFactorFADC2Phot  ( const Float_t f)                          { fMeanFFactorFADC2Phot   = f; }
-  void SetPedestal              ( const Float_t ped, const Float_t pedrms, const Float_t pederr);
-  void SetPed                   ( const Float_t ped, const Float_t pederr);
-  void SetPedRMS              ( const Float_t pedrms, const Float_t pedrmserr);
   void SetPheFFactorMethod      ( const Float_t f)                          { fPheFFactorMethod       = f; }
   void SetPheFFactorMethodVar   ( const Float_t f)                          { fPheFFactorMethodVar    = f; }
   void SetPheFFactorMethodLimit ( const Float_t f=fgPheFFactorMethodLimit ) { fPheFFactorMethodLimit  = f; }
   void SetNumSaturated          ( const Int_t   i)                          { fNumSaturated           = i; }
+  */
+  void SetFFactorMethodValid     (const Bool_t b = kTRUE );
+  void SetPedestal               (const Float_t ped, const Float_t pedrms, const Float_t pederr);
+  void SetPed                   ( const Float_t ped, const Float_t pederr);
+  void SetPedRMS              ( const Float_t pedrms, const Float_t pedrmserr);
+  void SetAbsTimeMean            (const Float_t f)                          { fAbsTimeMean            = f; }
+  void SetAbsTimeRms             (const Float_t f)                          { fAbsTimeRms             = f; }
+  void SetConversionHiLo         (const Float_t c=fgConversionHiLo        ) { fConversionHiLo         = c; }
+  void SetConversionHiLoErr      (const Float_t e=fgConversionHiLoErr     ) { fConversionHiLoVar      = e*e;}
+  void SetConversionHiLoSigma    (const Float_t s=fgConversionHiLoSigma   ) { fConversionHiLoSigma    = s; }
+  void SetConvFFactorRelErrLimit (const Float_t f=fgConvFFactorRelErrLimit) { fConvFFactorRelVarLimit = f*f;}
+  void SetMeanConvFADC2Phe       (const Float_t f)                          { fMeanConvFADC2Phe       = f; }
+  void SetMeanConvFADC2PheVar    (const Float_t f)                          { fMeanConvFADC2PheStatVar= f; }
+  void SetMeanConvFADC2PheSystVar(const Float_t f)                          { fMeanConvFADC2PheSystVar= f; }
+  void SetMeanFFactorFADC2Phot   (const Float_t f)                          { fMeanFFactorFADC2Phot   = f; }
+  void SetNumSaturated           (const Int_t   i)                          { fNumSaturated           = i; }
+  void SetPheFFactorMethod       (const Float_t f)                          { fPheFFactorMethod       = f; }
+  void SetPheFFactorMethodVar    (const Float_t f)                          { fPheFFactorMethodStatVar= f; }
+  void SetPheFFactorMethodSystVar(const Float_t f)                          { fPheFFactorMethodSystVar= f; }
+  void SetPheFFactorMethodLimit  (const Float_t f=fgPheFFactorMethodLimit ) { fPheFFactorMethodLimit  = f; }
   // Getters
+  Float_t GetAbsTimeMean             () const { return fAbsTimeMean;             }
+  Float_t GetAbsTimeRms              () const { return fAbsTimeRms;              }
+  Float_t GetConversionHiLo          () const { return fConversionHiLo;          }
+  Float_t GetConversionHiLoErr       () const;
+  Float_t GetConvertedMean           () const;
+  Float_t GetConvertedMeanErr        () const;
+  Float_t GetConvertedSigma          () const;
+  Float_t GetConvertedSigmaErr       () const;
+  Float_t GetConvertedRSigma         () const;
+  Float_t GetConvertedRSigmaErr      () const;
+  Float_t GetConvertedRSigmaSquare   () const;
+  Float_t GetMeanConvFADC2Phe        () const { return fMeanConvFADC2Phe;        }
+  Float_t GetMeanConvFADC2PheErr     () const;
+  Float_t GetMeanConvFADC2PheVar     () const { return fMeanConvFADC2PheVar;     }
+  Float_t GetMeanFFactorFADC2Phot    () const { return fMeanFFactorFADC2Phot;    }
+  Float_t GetMeanFFactorFADC2PhotErr () const;
+  Float_t GetMeanFFactorFADC2PhotVar () const { return fMeanFFactorFADC2PhotVar; }
+  Int_t   GetNumSaturated            () const { return fNumSaturated;            }
+  Float_t GetPed                     () const { return fPed;                     }
+  Float_t GetPedErr                  () const { return TMath::Sqrt(fPedVar);     }
+  Float_t GetPedRms                  () const;
+  Float_t GetPedRmsErr               () const;
+  Float_t GetPheFFactorMethod        () const { return fPheFFactorMethod;        }
+  Float_t GetPheFFactorMethodErr     () const;
+  Float_t GetPheFFactorMethodVar     () const { return fPheFFactorMethodVar;     }
+  Float_t GetPheFFactorMethodRelVar  () const;
+  Float_t GetRSigma                  () const;
+  Float_t GetRSigmaErr               () const;
+  Float_t GetRSigmaRelVar            () const;
+  Float_t GetRSigmaPerCharge         () const;
+  Float_t GetRSigmaPerChargeErr      () const;
+  Float_t GetAbsTimeMean        () const { return fAbsTimeMean;    }
+  Float_t GetAbsTimeRms         () const { return fAbsTimeRms;     }
+  Float_t GetConversionHiLo     () const { return fConversionHiLo; }
+  Float_t GetConversionHiLoErr  () const { return fConversionHiLoVar<0 ? -1 : TMath::Sqrt(fConversionHiLoVar);  }
+  Float_t GetConversionHiLoSigma() const { return fConversionHiLoSigma; }
+  Float_t GetConvertedMean      () const { return IsHiGainSaturation() ? fLoGainMean * fConversionHiLo : fHiGainMean; }
+  Float_t GetConvertedMeanErr   () const
+  {
+      if (!IsHiGainSaturation())
+          return GetMeanErr();
+      const Float_t logainrelvar = GetLoGainMeanRelVar();
+      return logainrelvar<0 ? -1 : TMath::Sqrt(logainrelvar + GetConversionHiLoRelVar()) * GetConvertedMean();
+  }
+  Float_t GetConvertedSigma() const { return IsHiGainSaturation() ? fLoGainSigma * fConversionHiLo : fHiGainSigma; }
+  Float_t GetConvertedSigmaErr() const
+  {
+    if (!IsHiGainSaturation())
+        return GetSigmaErr();
+  Bool_t IsFFactorMethodValid        () const;
+    if (fLoGainSigmaVar<0 || fLoGainSigma<0)
+        return -1.;
+    const Float_t sigmaRelVar = fLoGainSigmaVar/(fLoGainSigma*fLoGainSigma);
+    return TMath::Sqrt(sigmaRelVar+GetConversionHiLoRelVar()) * GetConvertedSigma();
+  }
+  Float_t GetConvertedRSigma() const
+  {
+      if (fRSigmaSquare < 0)
+          return -1;
+      const Float_t rsigma = TMath::Sqrt(fRSigmaSquare);
+      return IsHiGainSaturation() ? rsigma*fConversionHiLo : rsigma ;
+  }
+  Float_t GetConvertedRSigmaErr() const;
+  Float_t GetConvertedRSigmaSquare() const
+  {
+      if (fRSigmaSquare < 0)
+          return -1;
+      return IsHiGainSaturation() ? fRSigmaSquare*fConversionHiLo*fConversionHiLo : fRSigmaSquare ;
+  }
+  Float_t GetMeanConvFADC2Phe() const { return fMeanConvFADC2Phe; }
+  Float_t GetMeanConvFADC2PheErr  () const { return fMeanConvFADC2PheStatVar<0 ? -1 : TMath::Sqrt(fMeanConvFADC2PheStatVar); }
+  Float_t GetMeanConvFADC2PheSystErr() const { return fMeanConvFADC2PheSystVar<0 ? -1 : TMath::Sqrt(fMeanConvFADC2PheSystVar); }
+  Float_t GetMeanConvFADC2PheTotErr() const
+  {
+      if (fMeanConvFADC2PheSystVar<0 || fMeanConvFADC2PheStatVar<0)
+          return -1.;
+      return TMath::Sqrt(fMeanConvFADC2PheSystVar+fMeanConvFADC2PheStatVar);
+  }
+  Float_t GetFFactorFADC2Phe        () const { return gkFFactor;   }
+  Float_t GetMeanConvFADC2PheVar    () const { return fMeanConvFADC2PheStatVar; }
+  Float_t GetMeanConvFADC2PheSystVar() const { return fMeanConvFADC2PheSystVar; }
+  Float_t GetMeanFFactorFADC2Phot   () const { return fMeanFFactorFADC2Phot;    }
+  Float_t GetMeanFFactorFADC2PhotErr() const { return fMeanFFactorFADC2PhotVar<0 ? -1. : TMath::Sqrt(fMeanFFactorFADC2PhotVar); }
+  Float_t GetMeanFFactorFADC2PhotVar() const { return fMeanFFactorFADC2PhotVar; }
+  Int_t   GetNumSaturated           () const { return fNumSaturated;            }
+  Float_t GetPed                    () const { return fPed;                     }
+  Float_t GetPedErr                 () const { return TMath::Sqrt(fPedVar);     }
+  Float_t GetPedRms                 () const
+  {
+    if (!IsHiGainSaturation())
+        return fPedRms;
+    return fLoGainPedRmsSquare<0 ? -1 : TMath::Sqrt(fLoGainPedRmsSquare);
+  }
+  Float_t GetPedRmsErr              () const;
+  Float_t GetPheFFactorMethod       () const { return fPheFFactorMethod; }
+  Float_t GetPheFFactorMethodErr    () const { return fPheFFactorMethodStatVar<0 ? -1 : TMath::Sqrt(fPheFFactorMethodStatVar); }
+  Float_t GetPheFFactorMethodSystErr() const { return fPheFFactorMethodSystVar<0 ? -1 : TMath::Sqrt(fPheFFactorMethodSystVar); }
+  Float_t GetPheFFactorMethodTotErr () const
+  {
+      if (fPheFFactorMethodStatVar<0 || fPheFFactorMethodSystVar<0)
+          return -1.;
+      return TMath::Sqrt(fPheFFactorMethodStatVar+fPheFFactorMethodSystVar);
+  }
+  Float_t GetPheFFactorMethodVar    () const { return fPheFFactorMethodStatVar; }
+  Float_t GetPheFFactorMethodSystVar() const { return fPheFFactorMethodSystVar; }
+  Float_t GetPheFFactorMethodRelVar () const { return fPheFFactorMethodStatVar<=0 ? -1 : fPheFFactorMethodStatVar / (fPheFFactorMethod * fPheFFactorMethod); }
+  Float_t GetPheFFactorMethodRelSystVar() const { return fPheFFactorMethodSystVar<=0 ? -1. : fPheFFactorMethodSystVar / (fPheFFactorMethod * fPheFFactorMethod); }
+  Float_t GetRSigma                 () const { return fRSigmaSquare<0 ? -1 : TMath::Sqrt(fRSigmaSquare); }
+  Float_t GetRSigmaErr         () const;
+  Float_t GetRSigmaRelVar      () const;
+  Float_t GetRSigmaPerCharge   () const;
+  Float_t GetRSigmaPerChargeErr() const;
+  Bool_t IsFFactorMethodValid() const;
   // Calculations

trunk/MagicSoft/Mars/mcalib/MCalibrationHiLoPix.cc

-              r5749
+              r7188
 !   Author(s): Markus Gaug   02/2004 <mailto:markus@ifae.es>
+!
 !   Copyright: MAGIC Software Development, 2000-2004
+!   Copyright: MAGIC Software Development, 2000-2005
+!
 \* ======================================================================== */
 /////////////////////////////////////////////////////////////////////////////
 //                                                                         //
 // MCalibrationHiLoPix                                                     //
 //                                                                         //
 // Storage container for high-gain vs. low-gain charge calibration results //
 // of one Pixel (PMT).                                                     //
 // The following "calibration" constants can be retrieved:                 //
+//
+// MCalibrationHiLoPix
+//
+// Storage container for high-gain vs. low-gain charge calibration results
+// of one Pixel (PMT).
+// The following "calibration" constants can be retrieved:
 // - GetHiLoRatio(): The mean conversion High-gain vs. Low-gain
 //   with which the low-gain result has to be multiplied
+// - GetHiLoSigma(): The Gauss sigma of histogrammed High-gain vs. Low-gain ratios
+// - GetHiLoSigma(): The Gauss sigma of histogrammed High-gain vs.
+//   Low-gain ratios
 //
+// See also: MHCalibrationHiLoPix, MHCalibrationHiLoCam              //
+//                                                                         //
+// See also: MHCalibrationHiLoPix, MHCalibrationHiLoCam
+//
+// Class Version 2:
+//  + Float_t fOffsetPerSlice;  // Offset from fit (per FADC slice)
+//  + Float_t fGainRatio;       // Ratio of gains from fit
+//
 /////////////////////////////////////////////////////////////////////////////
 #include "MCalibrationHiLoPix.h"
 …
 //
 MCalibrationHiLoPix::MCalibrationHiLoPix(const char *name, const char *title)
+  : fOffsetPerSlice(-9999.), fGainRatio(-1.)
+{
 …
+}

trunk/MagicSoft/Mars/mcalib/MCalibrationHiLoPix.h

-              r5946
+              r7188
+{
 private:
+    Float_t fOffsetPerSlice;            // Offset from fit (per FADC slice)
+    Float_t fGainRatio;                 // Ratio of gains from fit
 public:
+  MCalibrationHiLoPix(const char *name=NULL, const char *title=NULL);
+  MCalibrationHiLoPix(const char *name=NULL, const char *title=NULL);
+  ~MCalibrationHiLoPix() {}
+  // Setter
+  void SetGainRatio     (const Float_t f) { fGainRatio      = f; }
+  void SetOffsetPerSlice(const Float_t f) { fOffsetPerSlice = f; }
+  // Getter
   Float_t GetHiLoChargeRatio()         const { return GetHiGainMean();     }
   Float_t GetHiLoChargeRatioErr()      const { return GetHiGainMeanErr();  }
 …
   Float_t GetHiLoTimeDiffSigmaErr()    const { return GetLoGainSigmaErr(); }
   Float_t GetHiLoTimeDiffProb()        const { return GetLoGainProb();     }
+  Float_t GetGainRatio     ()          const { return fGainRatio;          }
+  Float_t GetOffsetPerSlice()          const { return fOffsetPerSlice;     }
   ClassDef(MCalibrationHiLoPix, 1)      // Container HiLo conversion Calibration Results Pixel
+  ClassDef(MCalibrationHiLoPix, 2)      // Container HiLo conversion Calibration Results Pixel
 };

trunk/MagicSoft/Mars/mcalib/MCalibrationIntensityChargeCam.cc

-              r6412
+              r7188
           varerr[i] = pix.GetRSigmaErr();
+        }
       if (option.Contains("abstime"))
+      if (option.Contains("abstimemean"))
+        {
           var   [i] = pix.GetAbsTimeMean();
           varerr[i] = pix.GetAbsTimeRms();
+        }
+      if (option.Contains("abstimerms"))
+        {
+          var   [i] = pix.GetAbsTimeRms();
+          varerr[i] = pix.GetAbsTimeRms()/2.;
+        }
       if (option.Contains("blackout"))
+        {
 …
           var   [i] = pix.GetRSigmaPerCharge();
           varerr[i] = pix.GetRSigmaPerChargeErr();
+        }
+      if (option.Contains("conversionfactor"))
+        {
+          const MCalibrationChargePix &apix = (MCalibrationChargePix&)cam->GetAverageArea(0);
+          const Float_t mean = pix.GetConvertedMean();
+          const Float_t phe  = apix.GetPheFFactorMethod();
+          var[i]    = phe/mean;
+          varerr[i] = TMath::Sqrt(apix.GetPheFFactorMethodErr()*apix.GetPheFFactorMethodErr()/mean/mean
+                                  + phe*phe/mean/mean/mean/mean*pix.GetConvertedMeanErr()*pix.GetConvertedMeanErr());
+        }
+  }
 …
           if (option.Contains("rsigma"))
+              pvar = pix.GetRSigma();
+          if (option.Contains("abstimemean"))
+              pvar = pix.GetAbsTimeMean();
+          if (option.Contains("abstimerms"))
+              pvar = pix.GetAbsTimeRms();
+          if (option.Contains("conversionhilo"))
+              pvar = pix.GetConversionHiLo();
+          if (option.Contains("convertedmean"))
+              pvar = pix.GetConvertedMean();
+          if (option.Contains("convertedsigma"))
+              pvar = pix.GetConvertedSigma();
+          if (option.Contains("convertedrsigma"))
+              pvar = pix.GetConvertedRSigma();
+          if (option.Contains("meanconvfadc2phe"))
+              pvar = pix.GetMeanConvFADC2Phe();
+          if (option.Contains("meanffactorfadc2phot"))
+              pvar = pix.GetMeanFFactorFADC2Phot();
+          if (option.Contains("ped"))
+              pvar = pix.GetPed();
+          if (option.Contains("pedrms"))
+              pvar = pix.GetPedRms();
+          if (option.Contains("pheffactormethod"))
+              pvar = pix.GetPheFFactorMethod();
+          if (option.Contains("rsigmapercharge"))
+              pvar = pix.GetRSigmaPerCharge();
+          if (option.Contains("conversionfactor"))
+          {
+              const MCalibrationChargePix &apix = (MCalibrationChargePix&)cam->GetAverageArea(aidx);
+              pvar = apix.GetPheFFactorMethod()/pix.GetConvertedMean();
+          }
+          variab  += pvar;
+          variab2 += pvar*pvar;
+          num++;
+          camcharge.Fill(j,pvar);
+          camcharge.SetUsed(j);
+        }
+      if (num > 1)
+        {
+          variab  /= num;
+          variance = (variab2 - variab*variab*num) / (num-1);
+          vararea[i]    = variab;
+          varareaerr[i] = variance>0 ? TMath::Sqrt(variance/num) : 999999999.;
+          //
+          // Make also a Gauss-fit to the distributions. The RMS can be determined by
+          // outlier, thus we look at the sigma and the RMS and take the smaller one, afterwards.
+          //
+          h = camcharge.ProjectionS(TArrayI(),TArrayI(1,&aidx),"_py",750);
+          h->SetDirectory(NULL);
+          h->Fit("gaus","QL");
+          TF1 *fit = h->GetFunction("gaus");
+          Float_t ci2   = fit->GetChisquare();
+          Float_t sigma = fit->GetParameter(2);
+          if (ci2 > 500. || sigma > varareaerr[i])
+            {
+              h->Fit("gaus","QLM");
+              fit = h->GetFunction("gaus");
+              ci2   = fit->GetChisquare();
+              sigma = fit->GetParameter(2);
+            }
+          const Float_t mean  = fit->GetParameter(1);
+          const Float_t ndf   = fit->GetNDF();
+          *fLog << inf << "Camera Nr: " << i << endl;
+          *fLog << inf << option.Data() << " area idx: " << aidx << " Results: " << endl;
+          *fLog << inf << "Mean: " << Form("%4.3f",mean)
+                << "+-" << Form("%4.3f",fit->GetParError(1))
+                << "  Sigma: " << Form("%4.3f",sigma) << "+-" << Form("%4.3f",fit->GetParError(2))
+                << "  Chisquare: " << Form("%4.3f",ci2) << "  NDF  : " << ndf << endl;
+          delete h;
+          gROOT->GetListOfFunctions()->Remove(fit);
+          if (sigma<varareaerr[i] && ndf>2 && ci2<500.)
+            {
+              vararea   [i] = mean;
+              varareaerr[i] = sigma/TMath::Sqrt((Float_t)num);
+            }
+        }
+      else
+        {
+          vararea[i]    = -1.;
+          varareaerr[i] = 0.;
+        }
+      nr[i] = i;
+      nrerr[i] = 0.;
+    }
+  TGraphErrors *gr = new TGraphErrors(size,
+                                     nr.GetArray(),vararea.GetArray(),
+                                     nrerr.GetArray(),varareaerr.GetArray());
+  gr->SetTitle(Form("%s Area %3i Average",option.Data(),aidx));
+  gr->GetXaxis()->SetTitle("Camera Nr.");
+  //  gr->GetYaxis()->SetTitle("<Q> [1]");
+  return gr;
+}
+// -------------------------------------------------------------------
+//
+// Returns a TGraphErrors with the mean effective number of photon
+// vs. the calibration camera number. With the string 'method', different
+// calibration methods can be called.
+//
+TGraphErrors *MCalibrationIntensityChargeCam::GetPhotVsTime( const Option_t *method )
+{
+  const Int_t size = GetSize();
+  if (size == 0)
+    return NULL;
+  TString option(method);
+  TArrayF photarr(size);
+  TArrayF photarrerr(size);
+  TArrayF nr(size);
+  TArrayF nrerr(size);
+  for (Int_t i=0;i<GetSize();i++)
+    {
+      //
+      // Get the calibration cam from the intensity cam
+      //
+      MCalibrationChargeCam *cam = (MCalibrationChargeCam*)GetCam(i);
+      //
+      // Get the calibration pix from the calibration cam
+      //
+      Float_t phot    = 0.;
+      Float_t photerr = 0.;
+      if (option.Contains("BlindPixel"))
+        {
+          phot    = cam->GetNumPhotonsBlindPixelMethod();
+          photerr = cam->GetNumPhotonsBlindPixelMethodErr();
+        }
+      if (option.Contains("FFactor"))
+        {
+          phot    = cam->GetNumPhotonsFFactorMethod();
+          photerr = cam->GetNumPhotonsFFactorMethodErr();
+        }
+      if (option.Contains("PINDiode"))
+        {
+          phot    = cam->GetNumPhotonsPINDiodeMethod();
+          photerr = cam->GetNumPhotonsPINDiodeMethodErr();
+        }
+      photarr[i]       = phot;
+      photarrerr[i]    = photerr;
+      nr[i] = i;
+      nrerr[i] = 0.;
+    }
+  TGraphErrors *gr = new TGraphErrors(size,
+                                     nr.GetArray(),photarr.GetArray(),
+                                     nrerr.GetArray(),photarrerr.GetArray());
+  gr->SetTitle("Photons Average");
+  gr->GetXaxis()->SetTitle("Camera Nr.");
+  gr->GetYaxis()->SetTitle("<N_{phot}> [1]");
+  return gr;
+}
+// -------------------------------------------------------------------
+//
+// Returns a TGraphErrors with the mean effective number of photo-electrons per
+// area index 'aidx' vs. the calibration camera number
+//
+TGraphErrors *MCalibrationIntensityChargeCam::GetPhePerAreaVsTime( const Int_t aidx, const MGeomCam &geom)
+{
+  const Int_t size = GetSize();
+  if (size == 0)
+    return NULL;
+  TArrayF phearea(size);
+  TArrayF pheareaerr(size);
+  TArrayF time(size);
+  TArrayF timeerr(size);
+  for (Int_t i=0;i<GetSize();i++)
+    {
+      //
+      // Get the calibration cam from the intensity cam
+      //
+      MCalibrationChargeCam *cam = (MCalibrationChargeCam*)GetCam(i);
+      //
+      // Get the calibration pix from the calibration cam
+      //
+      const MCalibrationChargePix &apix = (MCalibrationChargePix&)cam->GetAverageArea(aidx);
+      const Float_t phe          = apix.GetPheFFactorMethod();
+      const Float_t pheerr       = apix.GetPheFFactorMethodErr();
+      phearea[i]       = phe;
+      pheareaerr[i]    = pheerr;
+      time[i] = i;
+      timeerr[i] = 0.;
+    }
+  TGraphErrors *gr = new TGraphErrors(size,
+                                     time.GetArray(),phearea.GetArray(),
+                                     timeerr.GetArray(),pheareaerr.GetArray());
+  gr->SetTitle(Form("Phes Area %d Average",aidx));
+  gr->GetXaxis()->SetTitle("Camera Nr.");
+  gr->GetYaxis()->SetTitle("<N_{phe}> [1]");
+  return gr;
+}
+// -------------------------------------------------------------------
+//
+// Returns a TGraphErrors with the event-by-event averaged charge per
+// area index 'aidx' vs. the calibration camera number
+//
+TGraphErrors *MCalibrationIntensityChargeCam::GetChargePerAreaVsTime( const Int_t aidx, const MGeomCam &geom)
+{
+  const Int_t size = GetSize();
+  if (size == 0)
+    return NULL;
+  TArrayF chargearea(size);
+  TArrayF chargeareaerr(size);
+  TArrayF nr(size);
+  TArrayF nrerr(size);
+  for (Int_t i=0;i<GetSize();i++)
+    {
+      //
+      // Get the calibration cam from the intensity cam
+      //
+      MCalibrationChargeCam *cam = (MCalibrationChargeCam*)GetCam(i);
+      //
+      // Get the calibration pix from the calibration cam
+      //
+      const MCalibrationChargePix &apix = (MCalibrationChargePix&)cam->GetAverageArea(aidx);
+      const Float_t charge          = apix.GetConvertedMean();
+      const Float_t chargeerr       = apix.GetConvertedSigma();
+      chargearea[i]       = charge;
+      chargeareaerr[i]    = chargeerr;
+      nr[i]    = i;
+      nrerr[i] = 0.;
+    }
+  TGraphErrors *gr = new TGraphErrors(size,
+                                     nr.GetArray(),chargearea.GetArray(),
+                                     nrerr.GetArray(),chargeareaerr.GetArray());
+  gr->SetTitle(Form("Averaged Charges Area Idx %d",aidx));
+  gr->GetXaxis()->SetTitle("Camera Nr.");
+  gr->GetYaxis()->SetTitle("<Q> [FADC cnts]");
+  return gr;
+}
+TH1F *MCalibrationIntensityChargeCam::GetVarFluctuations( const Int_t aidx, const MGeomCam &geom, const Option_t *varname )
+{
+  const Int_t size = GetSize();
+  if (size == 0)
+    return NULL;
+  TString option(varname);
+  option.ToLower();
+  TH1F *hist = new TH1F("hist",Form("%s - Rel. Fluctuations %s Pixel",option.Data(),aidx ? "Outer" : "Inner"),
+,0.,100.);
+  hist->SetXTitle("Relative Fluctuation [%]");
+  hist->SetYTitle("Nr. channels [1]");
+  hist->SetFillColor(kRed+aidx);
+  MCalibrationChargeCam *cam = (MCalibrationChargeCam*)GetCam();
+  //
+  // Loop over pixels
+  //
+  for (Int_t npix=0;npix<cam->GetSize();npix++)
+    {
+      if (geom[npix].GetAidx() != aidx)
+        continue;
+      Double_t variab   = 0.;
+      Double_t variab2  = 0.;
+      Double_t variance = 0.;
+      Int_t    num      = 0;
+      Float_t  pvar     = 0.;
+      Float_t  relrms   = 99.9;
+      //
+      // Loop over the Cams for each pixel
+      //
+      for (Int_t i=0; i<GetSize(); i++)
+        {
+          MCalibrationChargeCam *cam = (MCalibrationChargeCam*)GetCam(i);
+          //
+          // Get the calibration pix from the calibration cam
+          //
+          MCalibrationChargePix &pix = (MCalibrationChargePix&)(*cam)[npix];
+          //
+          // Don't use bad pixels
+          //
+          if (!pix.IsFFactorMethodValid())
+            continue;
+          if (option.Contains("rsigma"))
             pvar = pix.GetRSigma();
           if (option.Contains("abstime"))
+          if (option.Contains("abstimemean"))
             pvar = pix.GetAbsTimeMean();
+          if (option.Contains("abstimerms"))
+            pvar = pix.GetAbsTimeRms();
           if (option.Contains("conversionhilo"))
             pvar = pix.GetConversionHiLo();
 …
           if (option.Contains("rsigmapercharge"))
             pvar = pix.GetRSigmaPerCharge();
+          if (option.Contains("conversionfactor"))
+          {
+              const MCalibrationChargePix &apix = (MCalibrationChargePix&)cam->GetAverageArea(0);
+              pvar = apix.GetPheFFactorMethod()/pix.GetConvertedMean();
+          }
           variab  += pvar;
           variab2 += pvar*pvar;
           num++;
-          camcharge.Fill(j,pvar);
-          camcharge.SetUsed(j);
+        }
-      if (num > 1)
+        {
-          variab  /= num;
-          variance = (variab2 - variab*variab*num) / (num-1);
-          vararea[i] = variab;
-          if (variance > 0.)
-            varareaerr[i] = TMath::Sqrt(variance);
-          else
-            varareaerr[i] = 999999999.;
-          //
-          // Make also a Gauss-fit to the distributions. The RMS can be determined by
-          // outlier, thus we look at the sigma and the RMS and take the smaller one, afterwards.
-          //
-          h = camcharge.ProjectionS(TArrayI(),TArrayI(1,&aidx),"_py",750);
-          h->SetDirectory(NULL);
-          h->Fit("gaus","QL");
-          TF1 *fit = h->GetFunction("gaus");
-          Float_t ci2   = fit->GetChisquare();
-          Float_t sigma = fit->GetParameter(2);
-          if (ci2 > 500. || sigma > varareaerr[i])
+            {
-              h->Fit("gaus","QLM");
-              fit = h->GetFunction("gaus");
-              ci2   = fit->GetChisquare();
-              sigma = fit->GetParameter(2);
+            }
-          const Float_t mean  = fit->GetParameter(1);
-          const Float_t ndf   = fit->GetNDF();
-          *fLog << inf << "Camera Nr: " << i << endl;
-          *fLog << inf << option.Data() << " area idx: " << aidx << " Results: " << endl;
-          *fLog << inf << "Mean: " << Form("%4.3f",mean)
-                << "+-" << Form("%4.3f",fit->GetParError(1))
-                << "  Sigma: " << Form("%4.3f",sigma) << "+-" << Form("%4.3f",fit->GetParError(2))
-                << "  Chisquare: " << Form("%4.3f",fit->GetChisquare()) << "  NDF  : " << ndf << endl;
-          delete h;
-          gROOT->GetListOfFunctions()->Remove(fit);
-          if (sigma < varareaerr[i] && ndf > 2)
+            {
-              vararea   [i] = mean;
-              varareaerr[i] = sigma;
+            }
+        }
-      else
+        {
-          vararea[i]    = -1.;
-          varareaerr[i] = 0.;
+        }
-      nr[i] = i;
-      nrerr[i] = 0.;
+    }
-  TGraphErrors *gr = new TGraphErrors(size,
-                                     nr.GetArray(),vararea.GetArray(),
-                                     nrerr.GetArray(),varareaerr.GetArray());
-  gr->SetTitle(Form("%s Area %3i Average",option.Data(),aidx));
-  gr->GetXaxis()->SetTitle("Camera Nr.");
-  //  gr->GetYaxis()->SetTitle("<Q> [1]");
-  return gr;
+}
-// -------------------------------------------------------------------
-//
-// Returns a TGraphErrors with the mean effective number of photon
-// vs. the calibration camera number. With the string 'method', different
-// calibration methods can be called.
-//
-TGraphErrors *MCalibrationIntensityChargeCam::GetPhotVsTime( const Option_t *method )
+{
-  const Int_t size = GetSize();
-  if (size == 0)
-    return NULL;
-  TString option(method);
-  TArrayF photarr(size);
-  TArrayF photarrerr(size);
-  TArrayF nr(size);
-  TArrayF nrerr(size);
-  for (Int_t i=0;i<GetSize();i++)
+    {
-      //
-      // Get the calibration cam from the intensity cam
-      //
-      MCalibrationChargeCam *cam = (MCalibrationChargeCam*)GetCam(i);
-      //
-      // Get the calibration pix from the calibration cam
-      //
-      Float_t phot    = 0.;
-      Float_t photerr = 0.;
-      if (option.Contains("BlindPixel"))
+        {
-          phot    = cam->GetNumPhotonsBlindPixelMethod();
-          photerr = cam->GetNumPhotonsBlindPixelMethodErr();
+        }
-      if (option.Contains("FFactor"))
+        {
-          phot    = cam->GetNumPhotonsFFactorMethod();
-          photerr = cam->GetNumPhotonsFFactorMethodErr();
+        }
-      if (option.Contains("PINDiode"))
+        {
-          phot    = cam->GetNumPhotonsPINDiodeMethod();
-          photerr = cam->GetNumPhotonsPINDiodeMethodErr();
+        }
-      photarr[i]       = phot;
-      photarrerr[i]    = photerr;
-      nr[i] = i;
-      nrerr[i] = 0.;
+    }
-  TGraphErrors *gr = new TGraphErrors(size,
-                                     nr.GetArray(),photarr.GetArray(),
-                                     nrerr.GetArray(),photarrerr.GetArray());
-  gr->SetTitle("Photons Average");
-  gr->GetXaxis()->SetTitle("Camera Nr.");
-  gr->GetYaxis()->SetTitle("<N_phot> [1]");
-  return gr;
+}
-// -------------------------------------------------------------------
-//
-// Returns a TGraphErrors with the mean effective number of photo-electrons per
-// area index 'aidx' vs. the calibration camera number
-//
-TGraphErrors *MCalibrationIntensityChargeCam::GetPhePerAreaVsTime( const Int_t aidx, const MGeomCam &geom)
+{
-  const Int_t size = GetSize();
-  if (size == 0)
-    return NULL;
-  TArrayF phearea(size);
-  TArrayF pheareaerr(size);
-  TArrayF time(size);
-  TArrayF timeerr(size);
-  for (Int_t i=0;i<GetSize();i++)
+    {
-      //
-      // Get the calibration cam from the intensity cam
-      //
-      MCalibrationChargeCam *cam = (MCalibrationChargeCam*)GetCam(i);
-      //
-      // Get the calibration pix from the calibration cam
-      //
-      const MCalibrationChargePix &apix = (MCalibrationChargePix&)cam->GetAverageArea(aidx);
-      const Float_t phe          = apix.GetPheFFactorMethod();
-      const Float_t pheerr       = apix.GetPheFFactorMethodErr();
-      phearea[i]       = phe;
-      pheareaerr[i]    = pheerr;
-      time[i] = i;
-      timeerr[i] = 0.;
+    }
-  TGraphErrors *gr = new TGraphErrors(size,
-                                     time.GetArray(),phearea.GetArray(),
-                                     timeerr.GetArray(),pheareaerr.GetArray());
-  gr->SetTitle(Form("Phes Area %d Average",aidx));
-  gr->GetXaxis()->SetTitle("Camera Nr.");
-  gr->GetYaxis()->SetTitle("<N_phes> [1]");
-  return gr;
+}
-// -------------------------------------------------------------------
-//
-// Returns a TGraphErrors with the event-by-event averaged charge per
-// area index 'aidx' vs. the calibration camera number
-//
-TGraphErrors *MCalibrationIntensityChargeCam::GetChargePerAreaVsTime( const Int_t aidx, const MGeomCam &geom)
+{
-  const Int_t size = GetSize();
-  if (size == 0)
-    return NULL;
-  TArrayF chargearea(size);
-  TArrayF chargeareaerr(size);
-  TArrayF nr(size);
-  TArrayF nrerr(size);
-  for (Int_t i=0;i<GetSize();i++)
+    {
-      //
-      // Get the calibration cam from the intensity cam
-      //
-      MCalibrationChargeCam *cam = (MCalibrationChargeCam*)GetCam(i);
-      //
-      // Get the calibration pix from the calibration cam
-      //
-      const MCalibrationChargePix &apix = (MCalibrationChargePix&)cam->GetAverageArea(aidx);
-      const Float_t charge          = apix.GetConvertedMean();
-      const Float_t chargeerr       = apix.GetConvertedSigma();
-      chargearea[i]       = charge;
-      chargeareaerr[i]    = chargeerr;
-      nr[i]    = i;
-      nrerr[i] = 0.;
+    }
-  TGraphErrors *gr = new TGraphErrors(size,
-                                     nr.GetArray(),chargearea.GetArray(),
-                                     nrerr.GetArray(),chargeareaerr.GetArray());
-  gr->SetTitle(Form("Averaged Charges Area Idx %d",aidx));
-  gr->GetXaxis()->SetTitle("Camera Nr.");
-  gr->GetYaxis()->SetTitle("<Q> [FADC cnts]");
-  return gr;
+}
-TH1F *MCalibrationIntensityChargeCam::GetVarFluctuations( const Int_t aidx, const MGeomCam &geom, const Option_t *varname )
+{
-  const Int_t size = GetSize();
-  if (size == 0)
-    return NULL;
-  TString option(varname);
-  TH1F *hist = new TH1F("hist",Form("%s - Rel. Fluctuations %s Pixel",option.Data(),aidx ? "Outer" : "Inner"),
-,0.,100.);
-  hist->SetXTitle("Relative Fluctuation [%]");
-  hist->SetYTitle("Nr. channels [1]");
-  hist->SetFillColor(kRed+aidx);
-  MCalibrationChargeCam *cam = (MCalibrationChargeCam*)GetCam();
-  //
-  // Loop over pixels
-  //
-  for (Int_t npix=0;npix<cam->GetSize();npix++)
+    {
-      if (geom[npix].GetAidx() != aidx)
-        continue;
-      Double_t variab   = 0.;
-      Double_t variab2  = 0.;
-      Double_t variance = 0.;
-      Int_t    num      = 0;
-      Float_t  pvar     = 0.;
-      Float_t  relrms   = 99.9;
-      //
-      // Loop over the Cams for each pixel
-      //
-      for (Int_t i=0; i<GetSize(); i++)
+        {
-          MCalibrationChargeCam *cam = (MCalibrationChargeCam*)GetCam(i);
-          //
-          // Get the calibration pix from the calibration cam
-          //
-          MCalibrationChargePix &pix = (MCalibrationChargePix&)(*cam)[npix];
-          //
-          // Don't use bad pixels
-          //
-          if (!pix.IsFFactorMethodValid())
-            continue;
-          if (option.Contains("RSigma"))
-            pvar = pix.GetRSigma();
-          if (option.Contains("AbsTime"))
-            pvar = pix.GetAbsTimeMean();
-          if (option.Contains("ConversionHiLo"))
-            pvar = pix.GetConversionHiLo();
-          if (option.Contains("ConvertedMean"))
-            pvar = pix.GetConvertedMean();
-          if (option.Contains("ConvertedSigma"))
-            pvar = pix.GetConvertedSigma();
-          if (option.Contains("ConvertedRSigma"))
-            pvar = pix.GetConvertedRSigma();
-          if (option.Contains("MeanConvFADC2Phe"))
-            pvar = pix.GetMeanConvFADC2Phe();
-          if (option.Contains("MeanFFactorFADC2Phot"))
-            pvar = pix.GetMeanFFactorFADC2Phot();
-          if (option.Contains("Ped"))
-            pvar = pix.GetPed();
-          if (option.Contains("PedRms"))
-            pvar = pix.GetPedRms();
-          if (option.Contains("PheFFactorMethod"))
-            pvar = pix.GetPheFFactorMethod();
-          if (option.Contains("RSigmaPerCharge"))
-            pvar = pix.GetRSigmaPerCharge();
-          variab  += pvar;
-          variab2 += pvar*pvar;
-          num++;
+        }

trunk/MagicSoft/Mars/mcalib/MCalibrationRelTimeCalc.cc

-              r6963
+              r7188
 //                  - FinalizeBadPixels()
 //
+// Class Version 2:
+//  + Byte_t fCheckFlags; // Bit-field to hold the possible check flags
+//
+//
 //  Input Containers:
 //   MCalibrationRelTimeCam
 …
 #include "MBadPixelsPix.h"
 ClassImp(MCalibrationRelTimeCalc);
 …
 const Float_t MCalibrationRelTimeCalc::fgRelTimeResolutionLimit = 1.0;
 // --------------------------------------------------------------------------
 //
 …
   fName  = name  ? name  : "MCalibrationRelTimeCalc";
   fTitle = title ? title : "Task to finalize the relative time calibration";
+  SetCheckFitResults       ( kFALSE );
+  SetCheckDeviatingBehavior( kFALSE );
+  SetCheckHistOverflow     ( kFALSE );
+  SetCheckOscillations     ( kFALSE );
   SetRelTimeResolutionLimit();
   SetOutputPath();
   SetOutputFile("");
   Clear();
+}
 …
   PrintUncalibrated(MBadPixelsPix::kDeviatingTimeResolution,
                     Form("%s%2.1f%s","Time resolution less than ",fRelTimeResolutionLimit," FADC slices from Mean:   "));
+                    Form("%s%2.1f%s","Time resol. less than ",fRelTimeResolutionLimit," FADC sl. from Mean:     "));
   PrintUncalibrated(MBadPixelsPix::kRelTimeOscillating,
                     "Pixels with changing Rel. Times over time:             ");
 …
       MCalibrationRelTimePix &pix = (MCalibrationRelTimePix&)(*relcam)[i];
+      if (bad.IsUncalibrated( MBadPixelsPix::kDeviatingTimeResolution))
+        bad.SetUnsuitable(   MBadPixelsPix::kUnreliableRun   );
+      if (bad.IsUncalibrated( MBadPixelsPix::kRelTimeNotFitted))
+        bad.SetUnsuitable(   MBadPixelsPix::kUnreliableRun   );
+      if (bad.IsUncalibrated( MBadPixelsPix::kRelTimeOscillating))
+        bad.SetUnsuitable(   MBadPixelsPix::kUnreliableRun   );
+      if (bad.IsUnsuitable(   MBadPixelsPix::kUnsuitableRun    ))
+        pix.SetExcluded();
+    }
+      if (IsCheckDeviatingBehavior())
+          if (bad.IsUncalibrated(MBadPixelsPix::kDeviatingTimeResolution))
+              bad.SetUnsuitable(MBadPixelsPix::kUnreliableRun);
+      if (IsCheckFitResults())
+          if (bad.IsUncalibrated(MBadPixelsPix::kRelTimeNotFitted))
+              bad.SetUnsuitable(MBadPixelsPix::kUnreliableRun);
+      if (IsCheckOscillations())
+          if (bad.IsUncalibrated(MBadPixelsPix::kRelTimeOscillating))
+              bad.SetUnsuitable(MBadPixelsPix::kUnreliableRun);
+      if (bad.IsUnsuitable(MBadPixelsPix::kUnsuitableRun))
+          pix.SetExcluded();
+    }
+}
 …
   if (fGeom->InheritsFrom("MGeomCamMagic"))
     *fLog << " " << setw(7) << "Uncalibrated Pixels:            "
           << Form("%s%3i%s%3i","Inner: ",counts[0]," Outer: ",counts[1]) << endl;
+    *fLog << " " << setw(7) << "Uncalibrated Pixels:            Inner: "
+          << Form("%3i",counts[0]) << " Outer: " << Form("%3i",counts[1]) << endl;
   counts.Reset();
 …
     relcam->SetNumUnreliable(counts[aidx], aidx);
   *fLog << " " << setw(7) << "Unreliable Pixels:              "
         << Form("%s%3i%s%3i","Inner: ",counts[0]," Outer: ",counts[1]) << endl;
+  *fLog << " " << setw(7) << "Unreliable Pixels:              Inner: "
+        << Form("%3i",counts[0]) << " Outer: " << Form("%3i",counts[1]) << endl;
+}
 …
 void MCalibrationRelTimeCalc::SetOutputPath(TString path)
+{
   fOutputPath = path;
   if (fOutputPath.EndsWith("/"))
     fOutputPath = fOutputPath(0, fOutputPath.Length()-1);
+    fOutputPath = path;
+    if (fOutputPath.EndsWith("/"))
+        fOutputPath = fOutputPath(0, fOutputPath.Length()-1);
+}
 …
 const char* MCalibrationRelTimeCalc::GetOutputFile()
+{
+  return Form("%s/%s", (const char*)fOutputPath, (const char*)fOutputFile);
+}
+    return Form("%s/%s", (const char*)fOutputPath, (const char*)fOutputFile);
+}
+// --------------------------------------------------------------------------
+//
+// MCalibrationRelTimeCam.CheckFitResults: Yes
+// MCalibrationRelTimeCam.CheckDeviatingBehavior: Yes
+// MCalibrationRelTimeCam.CheckHistOverflow: Yes
+// MCalibrationRelTimeCam.CheckOscillations: Yes
+//
+Int_t MCalibrationRelTimeCalc::ReadEnv(const TEnv &env, TString prefix, Bool_t print)
+{
+    Bool_t rc = kFALSE;
+    if (IsEnvDefined(env, prefix, "CheckFitResults", print))
+    {
+        SetCheckFitResults(GetEnvValue(env, prefix, "CheckFitResults", IsCheckFitResults()));
+        rc = kTRUE;
+    }
+    if (IsEnvDefined(env, prefix, "CheckDeviatingBehavior", print))
+    {
+        SetCheckDeviatingBehavior(GetEnvValue(env, prefix, "CheckDeviatingBehavior", IsCheckDeviatingBehavior()));
+        rc = kTRUE;
+    }
+    if (IsEnvDefined(env, prefix, "CheckHistOverflow", print))
+    {
+        SetCheckHistOverflow(GetEnvValue(env, prefix, "CheckHistOverflow", IsCheckHistOverflow()));
+        rc = kTRUE;
+    }
+    if (IsEnvDefined(env, prefix, "CheckOscillations", print))
+    {
+        SetCheckOscillations(GetEnvValue(env, prefix, "CheckOscillations", IsCheckOscillations()));
+        rc = kTRUE;
+    }
+    return rc;
+}

trunk/MagicSoft/Mars/mcalib/MCalibrationRelTimeCalc.h

-              r5679
+              r7188
   MGeomCam                   *fGeom;             //! Camera geometry
+  // enums
+  enum  Check_t
+  {
+      kCheckHistOverflow,
+      kCheckDeviatingBehavior,
+      kCheckOscillations,
+      kCheckFitResults
+  };                                             // Possible Checks
+  Byte_t fCheckFlags;                            // Bit-field to hold the possible check flags
   enum  { kDebug };                              //  Possible flags
 …
   void   PrintUncalibrated( MBadPixelsPix::UncalibratedType_t typ, const char *text) const;
+  // Query checks
+  Bool_t IsCheckDeviatingBehavior() const { return TESTBIT(fCheckFlags,kCheckDeviatingBehavior); }
+  Bool_t IsCheckHistOverflow     () const { return TESTBIT(fCheckFlags,kCheckHistOverflow);      }
+  Bool_t IsCheckOscillations     () const { return TESTBIT(fCheckFlags,kCheckOscillations);      }
+  Bool_t IsCheckFitResults       () const { return TESTBIT(fCheckFlags,kCheckFitResults);        }
+  // MTask
   Bool_t ReInit     (MParList *pList);
   Int_t  Process    () { return kTRUE; }
   Int_t  PostProcess();
+  // MParContainer
+  Int_t  ReadEnv(const TEnv &env, TString prefix, Bool_t print);
 public:
   MCalibrationRelTimeCalc(const char *name=NULL, const char *title=NULL);
+  // TObject
   void Clear(const Option_t *o="");
+  // MCalibrationRelTimeCalc
   Int_t Finalize();
+  // Getter
   Bool_t IsDebug() const  {  return TESTBIT(fFlags,kDebug); }
   void SetDebug                 ( const Bool_t  b=kTRUE ) { b ? SETBIT(fFlags,kDebug) : CLRBIT(fFlags,kDebug); }
+  // Setter
   void SetOutputPath            ( TString path="."                         );
   void SetOutputFile            ( TString file="TimeCalibStat.txt"         ) { fOutputFile          = file; }
   void SetRelTimeResolutionLimit( const Float_t f=fgRelTimeResolutionLimit ) { fRelTimeResolutionLimit = f; }
+  ClassDef(MCalibrationRelTimeCalc, 1)   // Task finalizing the relative time Calibration
+  // Checks
+  void SetCheckFitResults(const Bool_t b=kTRUE)        { b ? SETBIT(fCheckFlags,kCheckFitResults) : CLRBIT(fCheckFlags,kCheckFitResults); }
+  void SetCheckDeviatingBehavior(const Bool_t b=kTRUE) { b ? SETBIT(fCheckFlags,kCheckDeviatingBehavior) : CLRBIT(fCheckFlags,kCheckDeviatingBehavior); }
+  void SetCheckHistOverflow(const Bool_t b=kTRUE)      { b ? SETBIT(fCheckFlags,kCheckHistOverflow) : CLRBIT(fCheckFlags,kCheckHistOverflow); }
+  void SetCheckOscillations(const Bool_t b=kTRUE)      { b ? SETBIT(fCheckFlags,kCheckOscillations) : CLRBIT(fCheckFlags,kCheckOscillations); }
+  void SetDebug(const Bool_t b=kTRUE)                  { b ? SETBIT(fFlags, kDebug) : CLRBIT(fFlags, kDebug); }
+  ClassDef(MCalibrationRelTimeCalc, 2)   // Task finalizing the relative time Calibration
 };

trunk/MagicSoft/Mars/mcalib/MMcCalibrationCalc.cc

-              r7005
+              r7188
 using namespace std;
+MMcCalibrationCalc::MMcCalibrationCalc(const char *name, const char *title)
+// --------------------------------------------------------------------------
+//
+// Constructor. Default value for fMinSize is 1000 ADC counts. This must be
+// set in general by the user (SetMinSize), since it depends among other things
+// on the signal extractor used.
+//
+MMcCalibrationCalc::MMcCalibrationCalc(const char *name, const char *title): fMinSize(1000)
+{
     fName  = name  ? name  : "MMcCalibrationCalc";
     fTitle = title ? title : "Calculate and write conversion factors into MCalibrationChargeCam and MCalibrationQECam containers";
-    fHistADC2PhotEl = new TH1F(AddSerialNumber("ADC2PhotEl"), "log10(fPhotElfromShower/fSize)", 1500, -3., 3.);
-    fHistADC2PhotEl->SetXTitle("log_{10}(fPhotElfromShower / fSize) [photel/ADC count]");
-    fHistPhot2PhotEl = new TH1F(AddSerialNumber("Phot2PhotEl"), "Photon conversion efficiency", 1000, 0., 1.);
-    fHistPhot2PhotEl->SetXTitle("Overall photon conversion efficiency [photoelectron/photon]");
+}
 …
 Bool_t MMcCalibrationCalc::CheckRunType(MParList *pList) const
+{
     const MRawRunHeader *run = (MRawRunHeader*)pList->FindObject("MRawRunHeader");
+    const MRawRunHeader *run = (MRawRunHeader*)pList->FindObject(AddSerialNumber("MRawRunHeader"));
     if (!run)
+    {
+        *fLog << warn << "Warning - cannot check file type, MRawRunHeader not found." << endl;
+        return kTRUE;
+      *fLog << warn << "Warning - cannot check file type, " << AddSerialNumber("MRawRunHeader")
+            << " not found." << endl;
+      return kTRUE;
+    }
 …
 // --------------------------------------------------------------------------
 //
 // Make sure, that there is an MCalibrationCam Object in the Parameter List.
+// Look for all necessary containers and create histograms
 //
 Int_t MMcCalibrationCalc::PreProcess(MParList *pList)
+{
-    fHistADC2PhotEl->Reset();
-    fHistPhot2PhotEl->Reset();
     fADC2PhotEl = 0;
     fPhot2PhotEl = 0;
 …
         return kFALSE;
+    }
+    //
+    // Create histograms:
+    //
+    fHistADC2PhotEl = new TH1F(AddSerialNumber("ADC2PhotEl"), "log10(fPhotElfromShower/fSize)", 1500, -3., 3.);
+    fHistADC2PhotEl->SetXTitle("log_{10}(fPhotElfromShower / fSize) [photel/ADC count]");
+    fHistPhot2PhotEl = new TH1F(AddSerialNumber("Phot2PhotEl"), "Photon conversion efficiency", 1000, 0., 1.);
+    fHistPhot2PhotEl->SetXTitle("Overall photon conversion efficiency [photoelectron/photon]");
     return kTRUE;
 …
   // perpendicular to the camera plane.
   //
+  // FIXME! We should look for AddSerialNumber("MMcConfigRunHeader") but
+  // for the moment the stereo version of camera does not write one such
+  // header per telescope (it should!)
+  //
+  MMcConfigRunHeader* mcconfig = (MMcConfigRunHeader*) pList->FindObject("MMcConfigRunHeader");
+  MMcConfigRunHeader* mcconfig = (MMcConfigRunHeader*) pList->FindObject(AddSerialNumber("MMcConfigRunHeader"));
   if (!mcconfig)
+    {
       *fLog << err << "MMcConfigRunHeader" << " not found... aborting." << endl;
+      *fLog << err << AddSerialNumber("MMcConfigRunHeader") << " not found... aborting." << endl;
       return kFALSE;
+    }
 …
     //
     // Exclude events with low Size (larger fluctuations)
-    // FIXME? The present cut (1000 "inner-pixel-counts") is somehow
-    // arbitrary. Might it be optimized?
     //
     if (size < 1000)
+    if (size < fMinSize)
         return kTRUE;
 …
+    }
+    fPhot2PhotEl = fHistPhot2PhotEl->GetMean();   // Average quantum efficiency
+    fPhot2PhotEl = fHistPhot2PhotEl->GetMean();
+    // Average quantum efficiency. For now we will set this value for all pixels, although
+    // even in MC there may be differences from pixel to pixel, if the .dat file containing
+    // QE vs lambda, input of the camera simulation, has different QE curves for different
+    // pixels. FIXME?
+    MCalibrationQEPix &avqepix = (MCalibrationQEPix&)(fQECam->GetAverageArea(0));
+    avqepix.SetAverageQE(fPhot2PhotEl); // Needed by MCalibrateData!
     //
 …
         // FIXME: we are now assuming that all inner pixels have the same gain, and all
         // outer pixels have the same gain (different from inner ones though). This can
         // only be like this in camera 0.7, but may change in future versions of camera.
+        // only be like this in camera 0.7, but might change in future versions of camera.
         //
 …
+    }
     return kTRUE;
+}

trunk/MagicSoft/Mars/mcalib/MMcCalibrationCalc.h

-              r4710
+              r7188
                                          // outer pixels w.r.t inner ones
+    Float_t fMinSize;
+    // Minimum SIZE (before calibration, ADC counts) an event must have to be considered in the
+    // calculation of the calibration constants.
     TH1F*   fHistADC2PhotEl;
     TH1F*   fHistPhot2PhotEl; // Histograms for monitoring the calibration.
 …
     TH1F*   GetHistPhot2PhotEl() { return fHistPhot2PhotEl; }
+    void SetMinSize(Float_t x) { fMinSize = x; }
     ClassDef(MMcCalibrationCalc, 0)   // Task which obtains, for MC files, the calibration factor from ADC counts to photons.
 };

trunk/MagicSoft/Mars/mcalib/Makefile

r6662	r7188
40	40	MCalibrateRelTimes.cc \
41	41	MCalibrationIntensityCam.cc \
	42	MCalibrationIntensityConstCam.cc \
42	43	MCalibrationIntensityChargeCam.cc \
43	44	MCalibrationIntensityBlindCam.cc \

trunk/MagicSoft/Mars/merpp.cc

r7116	r7188
332	332	w->AddContainer("MCameraCalibration", "Camera");
333	333	w->AddContainer("MCameraCooling", "Camera");
	334	w->AddContainer("MCameraActiveLoad", "Camera");
334	335	w->AddContainer("MCameraHV", "Camera");
335	336	w->AddContainer("MCameraLV", "Camera");

trunk/MagicSoft/Mars/mhcalib/HCalibLinkDef.h

r6598	r7188
15	15	#pragma link C++ class MHCalibrationTestCam+;
16	16	#pragma link C++ class MHCalibrationHiLoCam+;
	17	#pragma link C++ class MHCalibrationHiLoPix+;
17	18	#pragma link C++ class MHCalibrationPulseTimeCam+;
18	19	#pragma link C++ class MHPedestalCam+;

trunk/MagicSoft/Mars/mhcalib/MHCalibrationCam.cc

-              r7095
+              r7188
 //  + Bool_t   fIsLoGainFitRanges;            // Are low-gain fit ranges defined?
 //
+// Class Version 5:
+//  + Int_t   fMaxNumEvts;                    // Max Number of events
 //
 /////////////////////////////////////////////////////////////////////////////
 …
 const Float_t MHCalibrationCam::fgProbLimit        = 0.0001;
 const Float_t MHCalibrationCam::fgOverflowLimit    = 0.005;
+const Int_t   MHCalibrationCam::fgMaxNumEvts       = 4096;
 const TString MHCalibrationCam::gsHistName   = "Hist";
 …
     fHistName(gsHistName),fHistTitle(gsHistTitle),
     fHistXTitle(gsHistXTitle),fHistYTitle(gsHistYTitle),
+    fCurrentNumEvts(0),
     fColor(MCalibrationCam::kNONE), fIntensBad(NULL),
     fBadPixels(NULL), fIntensCam(NULL), fCam(NULL), fGeom(NULL),
 …
     SetProbLimit();
     SetOverflowLimit();
+    SetMaxNumEvts();
     SetAverageing  (kTRUE);
 …
+{
   SetIsReset();
+  fCurrentNumEvts = 0;
   ResetHistTitles();
 …
+  }
+  fCurrentNumEvts = 0;
   return SetupHists(pList);
+}
 …
           MBadPixelsPix &bad = fIntensBad ? (*fIntensBad)[i] : (*fBadPixels)[i];
           if (bad.IsBad())
+          if (bad.IsUnsuitable(MBadPixelsPix::kUnsuitableRun))
             continue;
 …
         SetLoGain(fRunHeader->GetNumSamplesLoGain());
+    }
+  fCurrentNumEvts = 0;
   if (!ReInitHists(pList))
 …
 Bool_t MHCalibrationCam::Fill(const MParContainer *par, const Stat_t w)
+{
+  if (fCurrentNumEvts >= fMaxNumEvts)
+    return kTRUE;
+  fCurrentNumEvts++;
   SetIsReset(kFALSE);
 …
   if (GetNumExecutions() < 2)
     return kTRUE;
+  *fLog << inf << GetDescriptor() << ": Number of event used to fill histograms == " << fCurrentNumEvts << endl;
   if (fHiGainArray->GetSize() == 0 && fLoGainArray->GetSize() == 0)
 …
 void MHCalibrationCam::CalcAverageSigma()
+{
-  if (!fCam)
-    return;
   if (!IsAverageing())
     return;
+  MCalibrationCam *cam = fIntensCam ? fIntensCam->GetCam() : fCam;
+  if (!cam)
+    return;
   for (UInt_t j=0; j<fGeom->GetNumAreas(); j++)
+    {
       MCalibrationPix &pix    = fCam->GetAverageArea(j);
+      MCalibrationPix &pix    = cam->GetAverageArea(j);
       const Float_t numsqr    = TMath::Sqrt((Float_t)fAverageAreaNum[j]);
 …
       fAverageAreaRelSigmaVar[j] += pix.GetMeanRelVar();
       fAverageAreaRelSigmaVar[j] *= fAverageAreaRelSigma[j];
+    }
+  for (UInt_t j=0; j<fGeom->GetNumSectors(); j++)
+    {
+      MCalibrationPix &pix    = cam->GetAverageSector(j);
+      const Float_t numsqr    = TMath::Sqrt((Float_t)fAverageSectorNum[j]);
+      pix.SetSigma   (pix.GetSigma() * numsqr);
+      pix.SetSigmaVar(pix.GetSigmaErr() * pix.GetSigmaErr() * numsqr);
+    }
+}
 …
+    {
       *fLog << dbginf << GetDescriptor() << ": ID " << GetName()
+            << " "<<pix.GetPixId()
             << " HiGainSaturation: "   << pix.IsHiGainSaturation()
             << " HiGainMean: "         << hist.GetMean    ()
 …
+    {
       *fLog << dbginf << GetDescriptor() << "ID: " << hist.GetName()
+            << " HiGainSaturation: "   << pix.IsHiGainSaturation()
+            << " "<<pix.GetPixId()
+            << " HiGainSaturation: "   << pix.IsHiGainSaturation()
             << " LoGainMean: "         << hist.GetMean    ()
             << " LoGainMeanErr: "      << hist.GetMeanErr ()
 …
       rc = kTRUE;
+    }
+  if (IsEnvDefined(env, prefix, "MaxNumEvts", print))
+    {
+      SetMaxNumEvts(GetEnvValue(env, prefix, "MaxNumEvts", fMaxNumEvts));
+      rc = kTRUE;
+    }
   if (IsEnvDefined(env, prefix, "OverflowLimit", print))

trunk/MagicSoft/Mars/mhcalib/MHCalibrationCam.h

-              r7095
+              r7188
 private:
+  static const Double_t fgLowerFitLimitHiGain; //! The default for fLowerFitLimitHiGain (now at: 0)
+  static const Double_t fgUpperFitLimitHiGain; //! The default for fUpperFitLimitHiGain (now at: 0)
+  static const Double_t fgLowerFitLimitLoGain; //! The default for fLowerFitLimitLoGain (now at: 0)
+  static const Double_t fgUpperFitLimitLoGain; //! The default for fUpperFitLimitLoGain (now at: 0)
+  static const Int_t   fgPulserFrequency;  //! The default for fPulserFrequency (now set to: 500)
+  static const Float_t fgProbLimit;        //! The default for fProbLimit (now set to: 0.0001)
+  static const Float_t fgOverflowLimit;    //! The default for fOverflowLimit (now at: 0.005)
+  static const Double_t fgLowerFitLimitHiGain; //! Default for fLowerFitLimitHiGain
+  static const Double_t fgUpperFitLimitHiGain; //! Default for fUpperFitLimitHiGain
+  static const Double_t fgLowerFitLimitLoGain; //! Default for fLowerFitLimitLoGain
+  static const Double_t fgUpperFitLimitLoGain; //! Default for fUpperFitLimitLoGain
+  static const Int_t   fgPulserFrequency;  //! Default for fPulserFrequency
+  static const Float_t fgProbLimit;        //! Default for fProbLimit
+  static const Float_t fgOverflowLimit;    //! Default for fOverflowLimit
+  static const Int_t   fgMaxNumEvts;       //! Default for fMaxNumEvts
   static const TString gsHistName;         //! Default Histogram names
 …
   Float_t fNumHiGainSaturationLimit;      // Rel. amount sat. higain FADC slices until pixel is called saturated
   Float_t fNumLoGainSaturationLimit;      // Rel. amount sat. logain FADC slices until pixel is called saturated
+  Int_t   fMaxNumEvts;                    // Max Number of events
+  Int_t   fCurrentNumEvts;                //! Current number of events
   MArrayI fRunNumbers;                    // Numbers of runs used
 …
   void   DrawPixelContent( Int_t num )  const  {}
+  const MArrayI           &GetAverageAreaNum     ()          const { return fAverageAreaNum; }
   const Int_t              GetAverageAreas       ()          const;
         MHCalibrationPix  &GetAverageHiGainArea  (UInt_t i);
 …
         MHCalibrationPix  &GetAverageLoGainSector(UInt_t i);
   const MHCalibrationPix  &GetAverageLoGainSector(UInt_t i)  const;
+  const Int_t              GetAverageSectors     ()          const;
+  const MArrayI           &GetAverageSectorNum    ()         const { return fAverageSectorNum; }
+  const Int_t              GetAverageSectors      ()         const;
   const MCalibrationCam::PulserColor_t GetColor   ()     const  { return fColor;                    }
   const Float_t        GetNumHiGainSaturationLimit()     const  { return fNumHiGainSaturationLimit; }
 …
   void SetFirst                   ( const Axis_t f )       { fFirst   = f; }
   void SetLast                    ( const Axis_t f )       { fLast    = f; }
+  void SetMaxNumEvts              ( const Int_t  i=fgMaxNumEvts ) { fMaxNumEvts  = i; }
   void SetProbLimit               ( const Float_t f=fgProbLimit) { fProbLimit = f; }
 …
   void SetPulserFrequency      ( const Int_t   i=fgPulserFrequency )   { fPulserFrequency  = i; }
   ClassDef(MHCalibrationCam, 5) // Base Histogram class for Calibration Camera
+  ClassDef(MHCalibrationCam, 6) // Base Histogram class for Calibration Camera
 };

trunk/MagicSoft/Mars/mhcalib/MHCalibrationChargeCam.cc

-              r7126
+              r7188
 // from the mean are counted as Pickup events (stored in MHCalibrationPix::fPickup)
 //
 // Unless more than fNumHiGainSaturationLimit (default: 1%) of the overall FADC
+// If more than fNumHiGainSaturationLimit (default: 15%) of the overall FADC
 // slices show saturation, the following flag is set:
 // - MCalibrationChargePix::SetHiGainSaturation();
 …
 // sigmas in the camera.
 //
+// Class Version 2:
+//  + Float_t fNumLoGainBlackoutLimit; // Rel. amount blackout logain events until pixel is declared unsuitable
+//
 /////////////////////////////////////////////////////////////////////////////
 #include "MHCalibrationChargeCam.h"
 …
 const Int_t   MHCalibrationChargeCam::fgChargeHiGainNbins =  500;
 const Axis_t  MHCalibrationChargeCam::fgChargeHiGainFirst = -100.125;
 const Axis_t  MHCalibrationChargeCam::fgChargeHiGainLast  =  1899.875;
+const Axis_t  MHCalibrationChargeCam::fgChargeHiGainFirst = -98.;
+const Axis_t  MHCalibrationChargeCam::fgChargeHiGainLast  =  1902.;
 const Int_t   MHCalibrationChargeCam::fgChargeLoGainNbins =  500;
 const Axis_t  MHCalibrationChargeCam::fgChargeLoGainFirst = -100.25;
 const Axis_t  MHCalibrationChargeCam::fgChargeLoGainLast  =  899.75;
+const Axis_t  MHCalibrationChargeCam::fgChargeLoGainFirst = -99.;
+const Axis_t  MHCalibrationChargeCam::fgChargeLoGainLast  =  901.;
 const Float_t MHCalibrationChargeCam::fgProbLimit         = 0.00000001;
 const TString MHCalibrationChargeCam::gsHistName          = "Charge";
 …
 const Float_t MHCalibrationChargeCam::fgNumHiGainSaturationLimit = 0.15;
 const Float_t MHCalibrationChargeCam::fgNumLoGainSaturationLimit = 0.005;
+const Float_t MHCalibrationChargeCam::fgNumLoGainBlackoutLimit   = 0.05;
+const Float_t MHCalibrationChargeCam::fgLoGainBlackoutLimit      = 3.5;
+const Float_t MHCalibrationChargeCam::fgLoGainPickupLimit        = 3.5;
 const Float_t MHCalibrationChargeCam::fgTimeLowerLimit           = 1.;
 const Float_t MHCalibrationChargeCam::fgTimeUpperLimit           = 3.;
 …
 // - fTimeLowerLimit           to fgTimeLowerLimit
 // - fTimeUpperLimit           to fgTimeUpperLimit
+// - fNumLoGainBlackoutLimit   to fgNumLoGainBlackoutLimit
 //
 // - fNbins to fgChargeHiGainNbins
 …
   SetNumHiGainSaturationLimit(fgNumHiGainSaturationLimit);
   SetNumLoGainSaturationLimit(fgNumLoGainSaturationLimit);
+  SetNumLoGainBlackoutLimit  (fgNumLoGainBlackoutLimit);
   SetTimeLowerLimit();
 …
           pix.SetAbsTimeFirst(-0.5);
           pix.SetAbsTimeLast(logainsamples-0.5);
+          pix.SetPickupLimit(fgLoGainPickupLimit);
+          pix.SetBlackoutLimit(fgLoGainBlackoutLimit);
           InitHists(pix,(*badcam)[i],i);
 …
       const Float_t sumhi = pix.GetExtractedSignalHiGain();
       const Int_t   sathi = (Int_t)pix.GetNumHiGainSaturated();
+      const Bool_t  sathi = pix.GetNumHiGainSaturated()>0;
       if (IsOscillations())
 …
         histhi.FillHist(sumhi);
       histhi.AddSaturated(sathi);
+      histhi.AddSaturated(sathi);
       const Int_t aidx   = (*fGeom)[i].GetAidx();
 …
       fSumhiarea[aidx]  += sumhi;
-      fSathiarea[aidx]  += sathi;
       fSumhisector[sector]  += sumhi;
+      fSathisector[sector]  += sathi;
+      if (sathi)
+      {
+          fSathiarea[aidx]++;
+          fSathisector[sector]++;
+      }
       if (IsLoGain())
 …
             continue;
+          }
+        MCalibrationChargePix &pix = (MCalibrationChargePix&)(*chargecam)[i] ;
+        if (!pix.IsHiGainSaturation())
+            continue;
         h = histlo.GetHGausHist();
 …
+          }
+        MCalibrationChargePix  &pix    = (MCalibrationChargePix&)(*chargecam)[i] ;
+        if (pix.IsHiGainSaturation())
+          FinalizeAbsTimes(histlo, pix, bad, fFirstLoGain, fLastLoGain);
+        FinalizeAbsTimes(histlo, pix, bad, fFirstLoGain, fLastLoGain);
+      }
 …
                     MBadPixelsPix::kLoGainOscillating);
+  //
+  // Check for pixels which have the high-gain saturated, but the low-gain
+  // switch not applied in sufficient cases. Have to exclude these pixels,
+  // although they not abnormal. They simply cannot be calibrated...
+  //
+  for (Int_t i=0; i<fLoGainArray->GetSize(); i++)
+  {
+      MHCalibrationChargePix &histlo = (MHCalibrationChargePix&)(*this)(i);
+      if (histlo.IsExcluded())
+          continue;
+      MCalibrationChargePix  &pix = (MCalibrationChargePix&)(*chargecam)[i] ;
+      if (!pix.IsHiGainSaturation())
+          continue;
+      //
+      // Now,treat only low-gain calibrated pixels:
+      //
+      const Double_t lim = fNumLoGainBlackoutLimit*histlo.GetHGausHist()->GetEntries();
+      if (histlo.GetBlackout() <= lim)
+          continue;
+      MBadPixelsPix &bad = (*badcam)[i];
+      bad.SetUnsuitable(MBadPixelsPix::kUnsuitableRun);
+      bad.SetUncalibrated(MBadPixelsPix::kLoGainBlackout);
+  }
   return kTRUE;
 …
       MCalibrationPix  &pix    = (*chargecam)[i];
       if (bad.IsUncalibrated( MBadPixelsPix::kHiGainNotFitted ))
         if (!pix.IsHiGainSaturation())
           bad.SetUnsuitable(   MBadPixelsPix::kUnreliableRun    );
+      if (bad.IsUncalibrated(MBadPixelsPix::kHiGainNotFitted))
+          if (!pix.IsHiGainSaturation())
+              bad.SetUnsuitable(MBadPixelsPix::kUnreliableRun);
       if (bad.IsUncalibrated( MBadPixelsPix::kLoGainNotFitted ))
         if (pix.IsHiGainSaturation())
           bad.SetUnsuitable(   MBadPixelsPix::kUnreliableRun    );
       if (bad.IsUncalibrated( MBadPixelsPix::kLoGainSaturation ))
           bad.SetUnsuitable(   MBadPixelsPix::kUnsuitableRun    );
+      if (bad.IsUncalibrated(MBadPixelsPix::kLoGainNotFitted))
+          if (pix.IsHiGainSaturation())
+              bad.SetUnsuitable(MBadPixelsPix::kUnreliableRun);
+      if (bad.IsUncalibrated(MBadPixelsPix::kLoGainSaturation))
+          bad.SetUnsuitable(MBadPixelsPix::kUnsuitableRun);
       if (IsOscillations())
+        {
+          if (bad.IsUncalibrated( MBadPixelsPix::kHiGainOscillating ))
+            bad.SetUnsuitable(   MBadPixelsPix::kUnreliableRun    );
+          if (bad.IsUncalibrated( MBadPixelsPix::kLoGainOscillating ))
+            if (pix.IsHiGainSaturation())
+              bad.SetUnsuitable(   MBadPixelsPix::kUnreliableRun    );
+            if (bad.IsUncalibrated(MBadPixelsPix::kHiGainOscillating))
+                if (!pix.IsHiGainSaturation())
+                    bad.SetUnsuitable(MBadPixelsPix::kUnreliableRun);
+            if (bad.IsUncalibrated(MBadPixelsPix::kLoGainOscillating))
+                if (pix.IsHiGainSaturation())
+                    bad.SetUnsuitable(MBadPixelsPix::kUnreliableRun);
+        }
+    }
+}
 …
+  }
+  if (IsEnvDefined(env, prefix, "NumLoGainBlackoutLimit", print))
+  {
+      SetNumLoGainBlackoutLimit(GetEnvValue(env, prefix, "NumLoGainBlackoutLimit", fNumLoGainBlackoutLimit));
+      rc = kTRUE;
+  }
   TEnv refenv(fReferenceFile);

trunk/MagicSoft/Mars/mhcalib/MHCalibrationChargeCam.h

-              r7043
+              r7188
   static const TString gsAbsHistYTitle;              //! Default Histogram y-axis titles abs.times
+  static const Float_t fgNumHiGainSaturationLimit;   //! The default for fNumHiGainSaturationLimit
+  static const Float_t fgNumLoGainSaturationLimit;   //! The default for fNumLoGainSaturationLimit
+  static const Float_t fgNumHiGainSaturationLimit;   //! Default for fNumHiGainSaturationLimit
+  static const Float_t fgNumLoGainSaturationLimit;   //! Default for fNumLoGainSaturationLimit
+  static const Float_t fgNumLoGainBlackoutLimit;     //! Default for fNumLoGainBlackoutLimit (now at: 0.05)
+  static const Float_t fgLoGainBlackoutLimit;        //! Default for low-gain blackout limit (now at: 3.5)
+  static const Float_t fgLoGainPickupLimit;          //! Default for low-gian pickup limit   (now at: 3.5)
   static const Float_t fgTimeLowerLimit;             //! Default for fTimeLowerLimit
 …
   Axis_t  fLoGainFirst;                              // Lower histogram limit low gain
   Axis_t  fLoGainLast;                               // Upper histogram limit low gain
+  Float_t fNumLoGainBlackoutLimit;                   // Rel. amount blackout logain events until pixel is declared unsuitable
   TString fAbsHistName;                              // Histogram names abs.times
 …
   void SetLoGainLast    ( const Axis_t f )    { fLoGainLast    = f; }
+  void SetNumLoGainBlackoutLimit( const Float_t f=fgNumLoGainBlackoutLimit ) { fNumLoGainBlackoutLimit = f; }
   void SetReferenceFile ( const TString ref=fgReferenceFile ) { fReferenceFile  = ref; }
   void SetTimeLowerLimit( const Float_t f=fgTimeLowerLimit )  { fTimeLowerLimit = f; }
   void SetTimeUpperLimit( const Float_t f=fgTimeUpperLimit )  { fTimeUpperLimit = f; }
+  // MHCamEvent
   Bool_t GetPixelContent ( Double_t &val, Int_t idx, const MGeomCam &cam, Int_t type=0) const { return kTRUE; }
   void   DrawPixelContent( Int_t num )  const;
   ClassDef(MHCalibrationChargeCam, 1)   // Histogram class for Charge Camera Calibration
+  ClassDef(MHCalibrationChargeCam, 2)   // Histogram class for Charge Camera Calibration
 };

trunk/MagicSoft/Mars/mhcalib/MHCalibrationChargePix.h

r5137	r7188
42	42
43	43	// Draws
44		v~~irtual v~~oid Draw(Option_t *opt="");
	44	void Draw(Option_t *opt="");
45	45
46	46	ClassDef(MHCalibrationChargePix, 1) // Base Histogram class for Charge Pixel Calibration

trunk/MagicSoft/Mars/mhcalib/MHGausEvents.cc

-              r6800
+              r7188
   option.ToLower();
   Int_t win = 1;
+  Int_t win   = 1;
   Int_t nofit = 0;
   if (option.Contains("events"))
+    {
+      option.ReplaceAll("events","");
+      win += 1;
+    }
+    win += 1;
   if (option.Contains("fourier"))
+    {
+      option.ReplaceAll("fourier","");
+      win += 2;
+    }
+    win += 2;
   if (IsEmpty())
 …
+}
-const Double_t MHGausEvents::GetChiSquare()  const
+{
-  return ( fFGausFit ? fFGausFit->GetChisquare() : 0.);
+}
-const Double_t MHGausEvents::GetExpChiSquare()  const
+{
-  return ( fFExpFit ? fFExpFit->GetChisquare() : 0.);
+}
-const Int_t MHGausEvents::GetExpNdf()  const
+{
-  return ( fFExpFit ? fFExpFit->GetNDF() : 0);
+}
-const Double_t MHGausEvents::GetExpProb()  const
+{
-  return ( fFExpFit ? fFExpFit->GetProb() : 0.);
+}
-const Int_t MHGausEvents::GetNdf() const
+{
-  return ( fFGausFit ? fFGausFit->GetNDF() : 0);
+}
-const Double_t MHGausEvents::GetOffset()  const
+{
-  return ( fFExpFit ? fFExpFit->GetParameter(0) : 0.);
+}
-// --------------------------------------------------------------------------
-//
-// If fFExpFit exists, returns fit parameter 1 (Slope of Exponential fit),
-// otherwise 0.
-//
-const Double_t MHGausEvents::GetSlope()  const
+{
-  return ( fFExpFit ? fFExpFit->GetParameter(1) : 0.);
+}
 // --------------------------------------------------------------------------
 //
 …
   //  att.Copy(fHGausHist.GetXaxis());
+}
-// --------------------------------------------------------------------------
-//
-// Return kFALSE if number of entries is 0
-//
-const Bool_t MHGausEvents::IsEmpty() const
+{
-  return !(fHGausHist.GetEntries());
+}
-// --------------------------------------------------------------------------
-//
-// Return kTRUE  if number of entries is bin content of fNbins+1
-//
-const Bool_t MHGausEvents::IsOnlyOverflow() const
+{
-  return fHGausHist.GetEntries() == fHGausHist.GetBinContent(fNbins+1);
+}
-// --------------------------------------------------------------------------
-//
-// Return kTRUE  if number of entries is bin content of 0
-//
-const Bool_t MHGausEvents::IsOnlyUnderflow() const
+{
-  return fHGausHist.GetEntries() == fHGausHist.GetBinContent(0);
+}
-const Bool_t MHGausEvents::IsExcluded() const
+{
-  return TESTBIT(fFlags,kExcluded);
+}
-const Bool_t MHGausEvents::IsExpFitOK() const
+{
-  return TESTBIT(fFlags,kExpFitOK);
+}
-const Bool_t MHGausEvents::IsFourierSpectrumOK() const
+{
-  return TESTBIT(fFlags,kFourierSpectrumOK);
+}
-const Bool_t MHGausEvents::IsGausFitOK() const
+{
-  return TESTBIT(fFlags,kGausFitOK);
+}
 // -----------------------------------------------------------------------------------
 …
+}
-// --------------------------------------------------------------------------
-//
-// Set Excluded bit from outside
-//
-void MHGausEvents::SetExcluded(const Bool_t b)
+{
-    b ? SETBIT(fFlags,kExcluded) : CLRBIT(fFlags,kExcluded);
+}
-// -------------------------------------------------------------------
-//
-// The flag setters are to be used ONLY for Monte-Carlo!!
-//
-void  MHGausEvents::SetExpFitOK(const Bool_t b)
+{
-  b ? SETBIT(fFlags,kExpFitOK) : CLRBIT(fFlags,kExpFitOK);
+}
-// -------------------------------------------------------------------
-//
-// The flag setters are to be used ONLY for Monte-Carlo!!
-//
-void  MHGausEvents::SetFourierSpectrumOK(const Bool_t b)
+{
-  b ? SETBIT(fFlags,kFourierSpectrumOK) : CLRBIT(fFlags,kFourierSpectrumOK);
+}
-// -------------------------------------------------------------------
-//
-// The flag setters are to be used ONLY for Monte-Carlo!!
-//
-void  MHGausEvents::SetGausFitOK(const Bool_t b)
+{
-  b ? SETBIT(fFlags,kGausFitOK) : CLRBIT(fFlags,kGausFitOK);
+}
 // ----------------------------------------------------------------------------
 //

trunk/MagicSoft/Mars/mhcalib/MHGausEvents.h

-              r6800
+              r7188
 #ifndef ROOT_TH1
 #include <TH1.h>
+#endif
+#ifndef ROOF_TF1
+#include <TF1.h>
 #endif
 …
 class TVirtualPad;
 class TGraph;
-class MArrayF;
 class TH1F;
 class TH1I;
 class TF1;
 class MHGausEvents : public MH
+{
 …
   const static Int_t    fgPowerProbabilityBins; //! Default for fPowerProbabilityBins (now set to: 20)
   Float_t *CreateEventXaxis(Int_t n);  // Create an x-axis for the Event TGraphs
   Float_t *CreatePSDXaxis(Int_t n);    // Create an x-axis for the PSD TGraphs
+  Float_t *CreateEventXaxis(Int_t n);           // Create an x-axis for the Event TGraphs
+  Float_t *CreatePSDXaxis  (Int_t n);           // Create an x-axis for the PSD TGraphs
 protected:
-  Float_t  fEventFrequency;            // Event frequency in Hertz (to be set)
   Int_t    fBinsAfterStripping;        // Bins for the Gauss Histogram after stripping off the zeros at both ends
   UInt_t   fCurrentSize;               // Current size of the array fEvents
+  Float_t  fEventFrequency;            // Event frequency in Hertz (to be set)
   Byte_t   fFlags;                     // Bit field for the fit result bits
   Int_t    fPowerProbabilityBins;      // Bins for the projected power spectrum
 …
   // Draws
   void Draw(Option_t *option="");                  // *MENU*
   void DrawEvents(Option_t *option="");              // *MENU*
   void DrawPowerSpectrum(Option_t *option="");         // *MENU*
   void DrawPowerProjection(Option_t *option="");       // *MENU*
+  void Draw               ( Option_t *option="" ); // *MENU*
+  void DrawEvents         ( Option_t *option="" ); // *MENU*
+  void DrawPowerSpectrum  ( Option_t *option="" ); // *MENU*
+  void DrawPowerProjection( Option_t *option="" ); // *MENU*
   // Fill
   void   FillArray       ( const Float_t f );
   Bool_t FillHist        ( const Float_t f );
   Bool_t FillHistAndArray( const Float_t f );
+  void   FillArray        ( const Float_t f );
+  Bool_t FillHist         ( const Float_t f );
+  Bool_t FillHistAndArray ( const Float_t f );
   // Fits
   Bool_t FitGaus(  Option_t *option="RQ0",
                    const Double_t xmin=0.,
                    const Double_t xmax=0.);           // *MENU*
+  Bool_t FitGaus          ( Option_t *option="RQ0",
+                            const Double_t xmin=0.,
+                            const Double_t xmax=0.); // *MENU*
   // Inits
 …
   // Getters
   const Double_t GetChiSquare()          const;
   const Double_t GetExpChiSquare()       const;
   const Int_t    GetExpNdf()             const;
   const Double_t GetExpProb()            const;
+  const Double_t GetChiSquare()          const { return ( fFGausFit ? fFGausFit->GetChisquare() : 0.); }
+  const Double_t GetExpChiSquare()       const { return ( fFExpFit  ? fFExpFit->GetChisquare()  : 0.); }
+  const Int_t    GetExpNdf()             const { return ( fFExpFit  ? fFExpFit->GetNDF()        : 0 ); }
+  const Double_t GetExpProb()            const { return ( fFExpFit  ? fFExpFit->GetProb()       : 0.); }
         MArrayF *GetEvents()                   { return &fEvents;            }
   const MArrayF *GetEvents()             const { return &fEvents;            }
   const Float_t  GetEventFrequency () const { return fEventFrequency; }
   TF1     *GetFExpFit()                  { return fFExpFit;            }
+  const Float_t  GetEventFrequency ()    const { return fEventFrequency;     }
+        TF1     *GetFExpFit()                  { return fFExpFit;            }
   const TF1     *GetFExpFit()            const { return fFExpFit;            }
         TF1     *GetFGausFit()                 { return fFGausFit;           }
 …
   const Double_t GetMean()               const { return fMean;               }
   const Double_t GetMeanErr()            const { return fMeanErr;            }
   const Int_t    GetNdf()                const;
   const Int_t    GetNbins()               const { return fNbins;            }
   const Double_t GetOffset()             const;
+  const Int_t    GetNdf()                const { return ( fFGausFit ? fFGausFit->GetNDF()       : 0);  }
+  const Int_t    GetNbins()              const { return fNbins;              }
+  const Double_t GetOffset()             const { return ( fFExpFit  ? fFExpFit->GetParameter(0) : 0.); }
         MArrayF *GetPowerSpectrum()            { return fPowerSpectrum;      }
   const MArrayF *GetPowerSpectrum()      const { return fPowerSpectrum;      }
 …
   const Double_t GetSigma()              const { return fSigma;              }
   const Double_t GetSigmaErr()           const { return fSigmaErr;           }
   const Double_t GetSlope()              const;
+  const Double_t GetSlope()              const { return ( fFExpFit  ? fFExpFit->GetParameter(1) : 0.); }
   const Bool_t IsExcluded()              const;
   const Bool_t IsExpFitOK()              const;
   const Bool_t IsEmpty()                 const;
   const Bool_t IsFourierSpectrumOK()     const;
   const Bool_t IsGausFitOK()             const;
   const Bool_t IsOnlyOverflow()          const;
   const Bool_t IsOnlyUnderflow()         const;
+  const Bool_t   IsExcluded()            const { return TESTBIT(fFlags,kExcluded);          }
+  const Bool_t   IsExpFitOK()            const { return TESTBIT(fFlags,kExpFitOK);          }
+  const Bool_t   IsEmpty()               const { return !(fHGausHist.GetEntries());         }
+  const Bool_t   IsFourierSpectrumOK()   const { return TESTBIT(fFlags,kFourierSpectrumOK); }
+  const Bool_t   IsGausFitOK()           const { return TESTBIT(fFlags,kGausFitOK);         }
+  const Bool_t   IsOnlyOverflow()        const { return fHGausHist.GetEntries() == fHGausHist.GetBinContent(fNbins+1); }
+  const Bool_t   IsOnlyUnderflow()       const { return fHGausHist.GetEntries() == fHGausHist.GetBinContent(0);        }
   // Prints
 …
   // Setters
   void  SetEventFrequency   ( const Float_t  f                   ) { fEventFrequency = f;   }
   void  SetExcluded         ( const Bool_t   b=kTRUE             );
   void  SetExpFitOK         ( const Bool_t   b=kTRUE             );
   void  SetFourierSpectrumOK( const Bool_t   b=kTRUE             );
   void  SetGausFitOK        ( const Bool_t   b=kTRUE             );
+  void  SetExcluded         ( const Bool_t   b=kTRUE             ) { b ? SETBIT(fFlags,kExcluded)  : CLRBIT(fFlags,kExcluded); }
+  void  SetExpFitOK         ( const Bool_t   b=kTRUE             ) { b ? SETBIT(fFlags,kExpFitOK)  : CLRBIT(fFlags,kExpFitOK); }
+  void  SetFourierSpectrumOK( const Bool_t   b=kTRUE             ) { b ? SETBIT(fFlags,kFourierSpectrumOK) : CLRBIT(fFlags,kFourierSpectrumOK); }
+  void  SetGausFitOK        ( const Bool_t   b=kTRUE             ) { b ? SETBIT(fFlags,kGausFitOK) : CLRBIT(fFlags,kGausFitOK);}
   void  SetLast             ( const Double_t d                   ) { fLast           = d;   }
   void  SetFirst            ( const Double_t d                   ) { fFirst          = d;   }

trunk/MagicSoft/Mars/mhcalib/Makefile

r6598	r7188
41	41	MHCalibrationRelTimeCam.cc \
42	42	MHCalibrationHiLoCam.cc \
	43	MHCalibrationHiLoPix.cc \
43	44	MHCalibrationTestCam.cc \
44	45	MHCalibrationPulseTimeCam.cc \

trunk/MagicSoft/Mars/mhflux/MHFalseSource.h

-              r7178
+              r7188
     TObject *GetCatalog();
+    void MakeSymmetric(TH1 *h);
 private:
 …
     void ProjectOn(const TH3D &src, TH2D *h, TH2D *all);
     void ProjectOnOff(TH2D *h, TH2D *all);
-    void MakeSymmetric(TH1 *h);
 public:

trunk/MagicSoft/Mars/mimage/MNewImagePar.cc

-              r7176
+              r7188
     Float_t maxpix2 = 0;                                 // [#phot]
     /*
     MSignalPix *pix = 0;
     TIter Next(evt);
     while ((pix=(MSignalPix*)Next()))
+    */
+    const Bool_t ismagiclike =
+        geom.GetNumPixels() == 577 &&
+        geom.GetNumAreas()  == 2   &&
+        geom.GetPixRatio(396) > geom.GetPixRatio(397);
     UInt_t npix = evt.GetNumPixels();
     for (UInt_t i=0; i<npix; i++)
 …
            edgepix2 += nphot;
+        //
+        // Now we are working on absolute values of nphot, which
+        // must take pixel size into account
+        //
+        nphot *= geom.GetPixRatio(i/*pixid*/);
+        // count inner pixels: To dependent on MAGIC Camera --> FIXME
+        if (i<397){
+        const Double_t ratio = geom.GetPixRatio(i);
+        if (TMath::Nint(ratio)==1) // Means this is a small (= inner pixel)
+        {
             fInnerSize += nphot;
+            if(i>270){
+                edgepixin2 += nphot;
+                if(i>330)
+                    edgepixin1 += nphot;
+            // Do calculation of "inner leakage" only for MAGIC-like geometry,
+            // i.e., 577 pixels, pixels of 2 different areas, inner part
+            // from pixel 0 to 396:
+            if (ismagiclike)
+            {
+                if(i > 270) // last two "rings" of inner pixels
+                {
+                    edgepixin2 += nphot;
+                    if(i > 330)  // last "ring" of inner pixels
+                        edgepixin1 += nphot;
+                }
+            }
+        }
+        // Compute Concentration 1-2
+        //
+        // Now convert nphot from absolute number of photons or phe to signal
+        // density (divide by pixel area), to find the pixel with highest signal
+        // density:
+        //
+        nphot *= ratio;
+        // Look for signal density in two highest pixels:
         if (nphot>maxpix1)
+        {
 …
     fInnerLeakage1 = edgepixin1 / fInnerSize;
     fInnerLeakage2 = edgepixin2 / fInnerSize;
     fLeakage1 = edgepix1 / hillas.GetSize();
     fLeakage2 = edgepix2 / hillas.GetSize();
+    // FIXME?: in case the pixel with highest signal density is an outer pixel,
+    // the value of fConc (ratio of signal in two highest pixels to SIZE) should
+    // rather be 2*fConc1, under the simplest assumption that the light density
+    // inside the outer (large) pixel is uniform.
     fConc  = (maxpix1+maxpix2)/hillas.GetSize();         // [ratio]
     fConc1 = maxpix1/hillas.GetSize();                   // [ratio]

trunk/MagicSoft/Mars/mjobs/MJCalibration.cc

-              r7126
+              r7188
       gPad->SetBorderMode(0);
       if (geomcam.InheritsFrom("MGeomCamMagic"))
         DisplayDoubleProject(&disp35, "noisy", "dead");
+        DisplayDoubleProject(&disp37, "noisy", "dead");
       //
 …
       // for the datacheck, fix the ranges!!
       //
       const Double_t max = 12.;
+      const Double_t max = 10.;
       obj8->SetMinimum(0.);
       obj8->SetMaximum(max);
 …
       t1->SetTextColor(gStyle->GetColorPalette(Int_t(1./max*numcol)));
       t1->SetTextAlign(12);
-      TText *t2 = pave->AddText(Form("%s%3i%s","Signal Rel. error too large:                          ",
-                                       CountBadPixels(&disp24,2)," pixels"));
-      t2->SetTextColor(gStyle->GetColorPalette(Int_t(2./max*numcol)));
-      t2->SetTextAlign(12);
       TText *t4 = pave->AddText(Form("Low Gain Saturation:                                   %3i pixels",
                                        CountBadPixels(&disp24,3)));
       t4->SetTextColor(gStyle->GetColorPalette(Int_t(3./max*numcol)));
+                                       CountBadPixels(&disp24,2)));
+      t4->SetTextColor(gStyle->GetColorPalette(Int_t(2./max*numcol)));
       t4->SetTextAlign(12);
       TText *t5 = pave->AddText(Form("Mean Arr. Time In First Extraction Bin:      %3i pixels",
                                        CountBadPixels(&disp24,4)));
       t5->SetTextColor(gStyle->GetColorPalette(Int_t(4./max*numcol)));
+                                       CountBadPixels(&disp24,3)));
+      t5->SetTextColor(gStyle->GetColorPalette(Int_t(3./max*numcol)));
       t5->SetTextAlign(12);
       TText *t6 = pave->AddText(Form("Mean Arr. Time In Last 2 Extraction Bins:  %3i pixels",
                                        CountBadPixels(&disp24,5)));
       t6->SetTextColor(gStyle->GetColorPalette(Int_t(5./max*numcol)));
+                                       CountBadPixels(&disp24,4)));
+      t6->SetTextColor(gStyle->GetColorPalette(Int_t(4./max*numcol)));
       t6->SetTextAlign(12);
-      TText *t9 = pave->AddText(Form("Deviating Number of Photons:                    %3i pixels",
-                                       CountBadPixels(&disp24,6)));
-      t9->SetTextColor(gStyle->GetColorPalette(Int_t(6./max*numcol)));
-      t9->SetTextAlign(12);
       TText *t10= pave->AddText(Form("High-Gain Histogram Overflow:                  %3i pixels",
                                        CountBadPixels(&disp24,7 )));
       t10->SetTextColor(gStyle->GetColorPalette(Int_t(7./max*numcol)));
+                                       CountBadPixels(&disp24,5 )));
+      t10->SetTextColor(gStyle->GetColorPalette(Int_t(5./max*numcol)));
       t10->SetTextAlign(12);
       TText *t11= pave->AddText(Form("Low-Gain Histogram Overflow:                   %3i pixels",
                                        CountBadPixels(&disp24,8 )));
       t11->SetTextColor(gStyle->GetColorPalette(Int_t(8./max*numcol)));
+                                       CountBadPixels(&disp24,6 )));
+      t11->SetTextColor(gStyle->GetColorPalette(Int_t(6./max*numcol)));
       t11->SetTextAlign(12);
       TText *t12= pave->AddText(Form("Presumably dead from Ped. Rms:              %3i pixels",
                                        CountBadPixels(&disp24,9 )));
       t12->SetTextColor(gStyle->GetColorPalette(Int_t(9./max*numcol)));
+                                       CountBadPixels(&disp24,7 )));
+      t12->SetTextColor(gStyle->GetColorPalette(Int_t(7./max*numcol)));
       t12->SetTextAlign(12);
       TText *t13= pave->AddText(Form("Fluctuating Pulse Arrival Times:                 %3i pixels",
                                        CountBadPixels(&disp24,10)));
       t13->SetTextColor(gStyle->GetColorPalette(Int_t(10./max*numcol)));
+                                       CountBadPixels(&disp24,8 )));
+      t13->SetTextColor(gStyle->GetColorPalette(Int_t(8./max*numcol)));
       t13->SetTextAlign(12);
       TText *t17 = pave->AddText(Form("Deviating Number of Photo-electrons:       %3i pixels",
                                        CountBadPixels(&disp24,11)));
       t17->SetTextColor(gStyle->GetColorPalette(Int_t(11./max*numcol)));
+                                       CountBadPixels(&disp24,9 )));
+      t17->SetTextColor(gStyle->GetColorPalette(Int_t(9./max*numcol)));
       t17->SetTextAlign(12);
       TText *t14= pave->AddText(Form("Previously Excluded:                                   %3i pixels",
                                        CountBadPixels(&disp24,12)));
       t14->SetTextColor(gStyle->GetColorPalette(Int_t(12./max*numcol)));
+                                       CountBadPixels(&disp24,10)));
+      t14->SetTextColor(gStyle->GetColorPalette(Int_t(10./max*numcol)));
       t14->SetTextAlign(12);
+      TText *t15= pave->AddText(Form("Too many Low-Gain Blackout Events              %3i pixels",
+                                       CountBadPixels(&disp24,25 )));
+      t15->SetTextColor(gStyle->GetColorPalette(Int_t(11./max*numcol)));
+      t15->SetTextAlign(12);
       pave->Draw();
 …
       const MCalibrationHiLoPix &pix = (MCalibrationHiLoPix&)hilocam[i];
+      const Float_t ratio = pix.GetHiLoChargeRatio();
+      const Float_t sigma = pix.GetHiLoChargeRatioSigma();
+      const Float_t ratio  = pix.GetHiLoChargeRatio();
+      const Float_t raterr = pix.GetHiLoChargeRatioErr();
+      const Float_t sigma  = pix.GetHiLoChargeRatioSigma();
       if (ratio < 0.)
 …
       cpix.SetConversionHiLo(ratio);
+      cpix.SetConversionHiLoErr(sigma);
+      cpix.SetConversionHiLoErr(raterr);
+      cpix.SetConversionHiLoSigma(sigma);
+    }
 …
       tlist.AddToList(&steer);
-    tlist.AddToList(&fillcam);
     if (IsRelTimes())
+    {
 …
         tlist.AddToList(&timecalc);
+    }
+    tlist.AddToList(&fillcam);
     if (IsUseBlindPixel())
 …
     if (fIsPixelCheck)
+    {
         chargecam[fCheckedPixId].DrawClone("");
         chargecam(fCheckedPixId).DrawClone("");
+        chargecam[fCheckedPixId].DrawClone("datacheck");
+        chargecam(fCheckedPixId).DrawClone("datacheck");
         if (IsRelTimes())

trunk/MagicSoft/Mars/mjobs/MJPedestal.cc

-              r7130
+              r7188
     fDeadPixelCheck = GetEnv("DeadPixelsCheck", fDeadPixelCheck);
+    fBadPixelsFile = GetEnv("BadPixelsFile",fBadPixelsFile.Data());
+    fReferenceFile = GetEnv("ReferenceFile",fReferenceFile.Data());
+    ReadReferenceFile();
+    // ------------- Do not put simple resource below --------------
     // Setup an environment task
     MTaskEnv tenv("ExtractSignal");
 …
     // ..and store it
     SetExtractor((MExtractor*)tenv.GetTask());
-    fBadPixelsFile = GetEnv("BadPixelsFile",fBadPixelsFile.Data());
-    fReferenceFile = GetEnv("ReferenceFile",fReferenceFile.Data());
-    ReadReferenceFile();
     return kTRUE;
 …
     fcalib.SetInverted();
+    tlist.AddToList(&decode);
     if (fIsPulsePosCheck)
+    {
         fillpul.SetFilter(&fcalib);
-        tlist.AddToList(&decode);
         tlist.AddToList(&fcalib);
         tlist.AddToList(&fillpul);

trunk/MagicSoft/Mars/mmain/MEventDisplay.cc

-              r7176
+              r7188
     // Now read event...
     //
+    // FIXME: Can we safly replace ReadinEvent() by RedinEvent(1)?
     ReadinEvent();
     MReadTree *reader = (MReadTree*)fEvtLoop->FindTask("MRead");

trunk/MagicSoft/Mars/mpedestal/MExtractPedestal.cc

-              r7070
+              r7188
 const TString MExtractPedestal::fgNamePedestalCam = "MPedestalCam";
+const TString MExtractPedestal::fgNameRawEvtData  = "MRawEvtData";
 const UInt_t  MExtractPedestal::fgNumDump = 500;
 …
 MExtractPedestal::MExtractPedestal(const char *name, const char *title)
     : fGeom(NULL), fPedestalsIn(NULL), fPedestalsInter(NULL), fPedestalsOut(NULL),
       fExtractor(NULL), fExtractWinFirst(0), fExtractWinSize(0)
+      fExtractor(NULL), fExtractWinFirst(0), fExtractWinSize(0), fUseSpecialPixels(kFALSE)
+{
     fName  = name  ? name  : "MExtractPedestal";
 …
     SetNamePedestalCamIn();
     SetNamePedestalCamOut();
+    SetNamePedestalCamInter();
+    SetNameRawEvtData();
     SetNumDump();
 …
   Clear();
+  fRawEvt = (MRawEvtData*)pList->FindObject(AddSerialNumber(fNameRawEvtData));
+  if (!fRawEvt)
+      if (!fUseSpecialPixels)
+      {
+          *fLog << err << AddSerialNumber(fNameRawEvtData) << " not found... aborting." << endl;
+          return kFALSE;
+      }
   fRawEvt = (MRawEvtData*)pList->FindObject(AddSerialNumber("MRawEvtData"));
   if (!fRawEvt)
 …
 Int_t MExtractPedestal::Process()
+{
+    //
+    // Necessary check for extraction of special pixels
+    // together with data which does not yet have them
+    //
+    if (fSumx.GetSize()==0)
+        return kTRUE;
     if (fExtractor)
         fExtractor->SetNoiseCalculation(fRandomCalculation);
 …
 Bool_t MExtractPedestal::ReInit(MParList *pList)
+{
+  // If the size is not yet set, set the size
+  if (fSumx.GetSize()==0)
+  {
+      const Int_t npixels  = fPedestalsOut->GetSize();
+      const Int_t areas    = fPedestalsOut->GetNumAverageArea();
+      const Int_t sectors  = fPedestalsOut->GetNumAverageSector();
+      fSumx.  Set(npixels);
+      fSumx2. Set(npixels);
+      fSumAB0.Set(npixels);
+      fSumAB1.Set(npixels);
+      fAreaSumx.  Set(areas);
+      fAreaSumx2. Set(areas);
+      fAreaSumAB0.Set(areas);
+      fAreaSumAB1.Set(areas);
+      fAreaFilled.Set(areas);
+      fAreaValid .Set(areas);
+      fSectorSumx.  Set(sectors);
+      fSectorSumx2. Set(sectors);
+      fSectorSumAB0.Set(sectors);
+      fSectorSumAB1.Set(sectors);
+      fSectorFilled.Set(sectors);
+      fSectorValid .Set(sectors);
+      for (Int_t i=0; i<npixels; i++)
+      {
+          const UInt_t aidx   = (*fGeom)[i].GetAidx();
+          const UInt_t sector = (*fGeom)[i].GetSector();
+          fAreaValid  [aidx]  ++;
+          fSectorValid[sector]++;
+      }
+  }
+  if (fExtractor)
+  {
+      if (!((MTask*)fExtractor)->ReInit(pList))
+          return kFALSE;
+      SetExtractWindow(fExtractor->GetHiGainFirst(), (Int_t)TMath::Nint(fExtractor->GetNumHiGainSamples()));
+  }
+  return kTRUE;
+    // Necessary check for special pixels which might not yet have existed
+    if (!fRawEvt)
+    {
+        if (fRunHeader->GetFormatVersion() > 3)
+            return kTRUE;
+        *fLog << err << "ERROR - " << fNameRawEvtData << " [MRawEvtData] has not ";
+        *fLog << "been found and format version > 3... abort." << endl;
+        return kFALSE;
+    }
+    // If the size is not yet set, set the size
+    if (fSumx.GetSize()==0)
+    {
+        // Initialize the normal pixels (size of MPedestalCam already set by MGeomApply)
+        const Int_t npixels  = fPedestalsOut->GetSize();
+        fSumx.  Set(npixels);
+        fSumx2. Set(npixels);
+        fSumAB0.Set(npixels);
+        fSumAB1.Set(npixels);
+        if (fUseSpecialPixels)
+        {
+            // Initialize size of MPedestalCam in case of special pixels (not done by MGeomApply)
+            const UShort_t nspecial = fRunHeader->GetNumSpecialPixels();
+            if (nspecial == 0)
+            {
+                *fLog << warn << "WARNING - Number of special pixels is 0." << endl;
+                return kTRUE;
+            }
+            fPedestalsOut->InitSize((UInt_t)nspecial);
+        }
+        else
+        {
+            // Initialize the averaged areas and sectors (do not exist for special pixels)
+            const Int_t areas    = fPedestalsOut->GetNumAverageArea();
+            const Int_t sectors  = fPedestalsOut->GetNumAverageSector();
+            fAreaSumx.  Set(areas);
+            fAreaSumx2. Set(areas);
+            fAreaSumAB0.Set(areas);
+            fAreaSumAB1.Set(areas);
+            fAreaFilled.Set(areas);
+            fAreaValid .Set(areas);
+            fSectorSumx.  Set(sectors);
+            fSectorSumx2. Set(sectors);
+            fSectorSumAB0.Set(sectors);
+            fSectorSumAB1.Set(sectors);
+            fSectorFilled.Set(sectors);
+            fSectorValid .Set(sectors);
+            for (Int_t i=0; i<npixels; i++)
+            {
+                const UInt_t aidx   = (*fGeom)[i].GetAidx();
+                const UInt_t sector = (*fGeom)[i].GetSector();
+                fAreaValid  [aidx]  ++;
+                fSectorValid[sector]++;
+            }
+        }
+    }
+    if (fExtractor)
+    {
+        if (!((MTask*)fExtractor)->ReInit(pList))
+            return kFALSE;
+        SetExtractWindow(fExtractor->GetHiGainFirst(), (Int_t)TMath::Nint(fExtractor->GetNumHiGainSamples()));
+    }
+    return kTRUE;
+}
 Int_t MExtractPedestal::PostProcess()
+{
   fPedestalsIn = NULL;
   return fExtractor ? fExtractor->CallPostProcess() : kTRUE;
+    fPedestalsIn = NULL;
+    return fExtractor ? fExtractor->CallPostProcess() : kTRUE;
+}
 …
+    }
+    // find resource for fUseSpecialPixels
+    if (IsEnvDefined(env, prefix, "UseSpecialPixels", print))
+    {
+        SetUseSpecialPixels(GetEnvValue(env, prefix, "UseSpecialPixels", fUseSpecialPixels));
+        rc = kTRUE;
+    }
     // find resource for numeventsdump
     if (IsEnvDefined(env, prefix, "NumEventsDump", print))
 …
 void MExtractPedestal::CalcPixResults(const UInt_t nevts, const UInt_t pixid)
+{
     const Float_t sum  = fSumx[pixid];
     const Float_t sum2 = fSumx2[pixid];
+    const Double_t sum  = fSumx[pixid];
+    const Double_t sum2 = fSumx2[pixid];
     // 1. Calculate the mean of the sums:
     Float_t ped        = sum/nevts;
+    Double_t ped = sum/nevts;
     // 2. Calculate the Variance of the sums:
     Float_t var = (sum2-sum*sum/nevts)/(nevts-1.);
+    Double_t var = (sum2-sum*sum/nevts)/(nevts-1.);
     // 3. Calculate the amplitude of the 150MHz "AB" noise
     Float_t abOffs = (fSumAB0[pixid] - fSumAB1[pixid]) / nevts;
+    Double_t abOffs = (fSumAB0[pixid] - fSumAB1[pixid]) / nevts;
     // 4. Scale the mean, variance and AB-noise to the number of slices:
 …
     // 5. Calculate the RMS from the Variance:
     const Float_t rms = var<0 ? 0 : TMath::Sqrt(var);
+    const Double_t rms = var<0 ? 0 : TMath::Sqrt(var);
     (*fPedestalsOut)[pixid].Set(ped, rms, abOffs, nevts);
 …
 void MExtractPedestal::CalcAreaResults(const UInt_t nevts, const UInt_t napix, const UInt_t aidx)
+{
     const Float_t sum  = fAreaSumx[aidx];
     const Float_t sum2 = fAreaSumx2[aidx];
+    const Double_t sum  = fAreaSumx[aidx];
+    const Double_t sum2 = fAreaSumx2[aidx];
     // 1. Calculate the mean of the sums:
     Float_t ped        = sum/nevts;
+    Double_t ped = sum/nevts;
     // 2. Calculate the Variance of the sums:
     Float_t var = (sum2-sum*sum/nevts)/(nevts-1.);
+    Double_t var = (sum2/napix-sum*sum/nevts)/(nevts-1.);
     // 3. Calculate the amplitude of the 150MHz "AB" noise
     Float_t abOffs = (fAreaSumAB0[aidx] - fAreaSumAB1[aidx]) / nevts;
+    Double_t abOffs = (fAreaSumAB0[aidx] - fAreaSumAB1[aidx]) / nevts;
     // 4. Scale the mean, variance and AB-noise to the number of slices:
 …
     // 6. Calculate the RMS from the Variance:
     const Float_t rms = var<0 ? 0 : TMath::Sqrt(var);
+    const Double_t rms = var<0 ? 0 : TMath::Sqrt(var);
     fPedestalsOut->GetAverageArea(aidx).Set(ped, rms, abOffs, nevts);
 …
 void MExtractPedestal::CalcSectorResults(const UInt_t nevts, const UInt_t nspix, const UInt_t sector)
+{
     const Float_t sum  = fSectorSumx[sector];
     const Float_t sum2 = fSectorSumx2[sector];
+    const Double_t sum  = fSectorSumx[sector];
+    const Double_t sum2 = fSectorSumx2[sector];
     // 1. Calculate the mean of the sums:
     Float_t ped        = sum/nevts;
+    Double_t ped        = sum/nevts;
     // 2. Calculate the Variance of the sums:
     Float_t var = (sum2-sum*sum/nevts)/(nevts-1.);
+    Double_t var = (sum2/nspix-sum*sum/nevts)/(nevts-1.);
     // 3. Calculate the amplitude of the 150MHz "AB" noise
     Float_t abOffs = (fSectorSumAB0[sector] - fSectorSumAB1[sector]) / nevts;
+    Double_t abOffs = (fSectorSumAB0[sector] - fSectorSumAB1[sector]) / nevts;
     // 4. Scale the mean, variance and AB-noise to the number of slices:
 …
     // 6. Calculate the RMS from the Variance:
     const Float_t rms = var<0 ? 0 : TMath::Sqrt(var);
+    const Double_t rms = var<0 ? 0 : TMath::Sqrt(var);
     fPedestalsOut->GetAverageSector(sector).Set(ped, rms, abOffs, nevts);
 …
     *fLog << "Intermediate Storage is       " << (fIntermediateStorage?"on":"off") << endl;
     *fLog << "Pedestal Update is            " << (fPedestalUpdate?"on":"off") << endl;
+    *fLog << "Special pixel mode            " << (fUseSpecialPixels?"on":"off") << endl;
     if (fPedestalUpdate)
+    {

trunk/MagicSoft/Mars/mpedestal/MExtractPedestal.h

-              r6913
+              r7188
 private:
   static const TString  fgNamePedestalCam;  //! "MPedestalCam"
+  static const TString  fgNameRawEvtData;   //! "MRawEvtData"
   static const UInt_t   fgNumDump;          //!
   TString fNamePedestalCamIn;        // Name of the incoming 'MPedestalCam' container
   TString fNamePedestalCamOut;       // Name of the outgoing 'MPedestalCam' container
+  TString fNamePedestalCamInter;     // Name of the intermediate 'MPedestalCam' container
+  TString fNamePedestalCamInter;     // Name of the intermediate 'MPedestalCam' container
+  TString fNameRawEvtData;           // Name of MRawEvtData
   Bool_t  fRandomCalculation;        // Is pedestalextraction by extractor random?
 …
   Bool_t  fPedestalUpdate;           // Flag if the pedestal shall be updated after every fNumEventsDump
+  Bool_t  fUseSpecialPixels;         // Flag if the special pixels shall be treated
   MArrayD fSumx;                     // sum of values
 …
   // names
+  void SetNamePedestalCamIn   (const char *name=fgNamePedestalCam.Data()) { fNamePedestalCamIn    = name; }
+  void SetNamePedestalCamInter(const char *name=fgNamePedestalCam.Data()) { fNamePedestalCamInter = name; }
+  void SetNamePedestalCamOut  (const char *name=fgNamePedestalCam.Data()) { fNamePedestalCamOut   = name; }
+  void SetNamePedestalCamIn   (const char *name=fgNamePedestalCam) { fNamePedestalCamIn    = name; }
+  void SetNamePedestalCamInter(const char *name=fgNamePedestalCam) { fNamePedestalCamInter = name; }
+  void SetNamePedestalCamOut  (const char *name=fgNamePedestalCam) { fNamePedestalCamOut   = name; }
+  void SetNameRawEvtData      (const char *name=fgNameRawEvtData)  { fNameRawEvtData       = name; }
   // pointers
 …
   // flags
   void SetIntermediateStorage (Bool_t b=kTRUE) { fIntermediateStorage = b; }
+  void SetUseSpecialPixels    (Bool_t b=kTRUE) { fUseSpecialPixels    = b; }
   void SetPedestalUpdate      (Bool_t b=kTRUE) { fPedestalUpdate      = b; }
   void SetRandomCalculation   (Bool_t b=kTRUE) { fRandomCalculation   = b; }

trunk/MagicSoft/Mars/mpedestal/MMcPedestalCopy.cc

-              r3803
+              r7188
 Bool_t MMcPedestalCopy::CheckRunType(MParList *pList) const
+{
     const MRawRunHeader *run = (MRawRunHeader*)pList->FindObject("MRawRunHeader");
+    const MRawRunHeader *run = (MRawRunHeader*)pList->FindObject(AddSerialNumber("MRawRunHeader"));
     if (!run)
+    {
+        *fLog << warn << dbginf << "Warning - cannot check file type, MRawRunHeader not found." << endl;
+      *fLog << warn << dbginf << "Warning - cannot check file type, "
+            << AddSerialNumber("MRawRunHeader") << " not found." << endl;
         return kTRUE;
+    }
 …
+      }
     MMcRunHeader *mcrun = (MMcRunHeader*)pList->FindObject("MMcRunHeader");
+    MMcRunHeader *mcrun = (MMcRunHeader*)pList->FindObject(AddSerialNumber("MMcRunHeader"));
     if (!mcrun)
+        *fLog << warn << dbginf << "MMcRunHeader not found... assuming camera<0.7" << endl;
+      *fLog << warn << dbginf << AddSerialNumber("MMcRunHeader")
+            << " not found... assuming camera<0.7" << endl;
     const int num = pedcam->GetSize();

trunk/MagicSoft/Mars/mpedestal/MPedCalcPedRun.cc

-              r7130
+              r7188
   //
+  // In the other case, produce the MPedestalCam to be use the subsequent (non-pedestal) events
+  // In the other case, produce the MPedestalCam to be use the subsequent
+  // (non-pedestal) events
   //
   *fLog << inf << "Finalizing pedestal calculations..." << flush;
 …
 Int_t MPedCalcPedRun::Calc()
+{
   if (fIsNotPedRun && !IsPedBitSet())
     return kTRUE;
 …
   while (pixel.Next())
+  {
+      const UInt_t idx    = pixel.GetPixelId();
+      const UInt_t aidx   = (*fGeom)[idx].GetAidx();
+      const UInt_t sector = (*fGeom)[idx].GetSector();
+      const UInt_t idx = pixel.GetPixelId();
       Float_t sum = 0.;
 …
         CalcSums(pixel, sum, ab0, ab1);
+      fSumx[idx]           += sum;
+      fSumx[idx] += sum;
+      if (fIntermediateStorage)
+        (*fPedestalsInter)[idx].Set(sum,0.,0.,fUsedEvents);
+      const Float_t sqrsum = sum*sum;
+      fSumx2[idx]  += sqrsum;
+      fSumAB0[idx] += ab0;
+      fSumAB1[idx] += ab1;
+      if (fUseSpecialPixels)
+        continue;
+      const UInt_t aidx   = (*fGeom)[idx].GetAidx();
+      const UInt_t sector = (*fGeom)[idx].GetSector();
       fAreaSumx[aidx]      += sum;
       fSectorSumx[sector]  += sum;
+      if (fIntermediateStorage)
+        (*fPedestalsInter)[idx].Set(sum,0.,0.,fUsedEvents);
+      const Float_t sqrsum   = sum*sum;
+      fSumx2[idx]           += sqrsum;
+      fAreaSumx2[aidx]      += sqrsum;
+      fSectorSumx2[sector]  += sqrsum;
+      fSumAB0[idx]         += ab0;
+      fSumAB1[idx]         += ab1;
+      fAreaSumx2[aidx]     += sqrsum;
+      fSectorSumx2[sector] += sqrsum;
       fAreaSumAB0[aidx]    += ab0;
 …
+}
 // --------------------------------------------------------------------------
 //
 …
 Int_t MPedCalcPedRun::Finalize()
+{
+  if (fUsedEvents == 0)
+    if (fUsedEvents == 0)
+        return kTRUE;
+    //
+    // Necessary check for extraction of special pixels
+    // together with data which does not yet have them
+    //
+    if (fSumx.GetSize()==0)
+        return kTRUE;
+    MRawEvtPixelIter pixel(fRawEvt);
+    while (pixel.Next())
+        CalcPixResults(fUsedEvents, pixel.GetPixelId());
+    if (!fUseSpecialPixels)
+    {
+        //
+        // Loop over the (two) area indices to get the averaged pedestal per aidx
+        //
+        for (UInt_t aidx=0; aidx<fAreaValid.GetSize(); aidx++)
+            if (fAreaValid[aidx]>0)
+                CalcAreaResults(fUsedEvents, fAreaValid[aidx], aidx);
+        //
+        // Loop over the (six) sector indices to get the averaged pedestal per sector
+        //
+        for (UInt_t sector=0; sector<fSectorValid.GetSize(); sector++)
+            if (fSectorValid[sector]>0)
+                CalcSectorResults(fUsedEvents, fSectorValid[sector], sector);
+    }
+    fPedestalsOut->SetTotalEntries(fUsedEvents*fExtractWinSize);
+    fPedestalsOut->SetReadyToSave();
     return kTRUE;
-  MRawEvtPixelIter pixel(fRawEvt);
-  while (pixel.Next())
-    CalcPixResults(fUsedEvents, pixel.GetPixelId());
-  //
-  // Loop over the (two) area indices to get the averaged pedestal per aidx
-  //
-  for (UInt_t aidx=0; aidx<fAreaValid.GetSize(); aidx++)
-    if (fAreaValid[aidx]>0)
-      CalcAreaResults(fUsedEvents, fAreaValid[aidx], aidx);
-  //
-  // Loop over the (six) sector indices to get the averaged pedestal per sector
-  //
-  for (UInt_t sector=0; sector<fSectorValid.GetSize(); sector++)
-    if (fSectorValid[sector]>0)
-      CalcSectorResults(fUsedEvents, fSectorValid[sector], sector);
-  fPedestalsOut->SetTotalEntries(fUsedEvents*fExtractWinSize);
-  fPedestalsOut->SetReadyToSave();
-  return kTRUE;
+}
 Int_t MPedCalcPedRun::PostProcess()
+{
+  if (!Finalize())
+    return kFALSE;
+  return MExtractPedestal::PostProcess();
+    if (!fRawEvt)
+        return kTRUE;
+    if (!Finalize())
+        return kFALSE;
+    return MExtractPedestal::PostProcess();
+}
 …
 void MPedCalcPedRun::Print(Option_t *o) const
+{
     MExtractPedestal::Print(o);

trunk/MagicSoft/Mars/mpedestal/MPedestalCam.cc

-              r5558
+              r7188
+}
+void MPedestalCam::PrintArr(const TCollection &list) const
+{
+    Int_t id = 0;
+    TIter Next(&list);
+    MPedestalPix *pix = 0;
+    while ((pix=(MPedestalPix*)Next()))
+        *fLog << Form("Nr.: %4i  Pedestal: %5.2f    RMS: %5.2f    ABOffset: %5.2f ",
+                      id++, pix->GetPedestal(), pix->GetPedestalRms(), pix->GetPedestalABoffset()) << endl;
+}
 void MPedestalCam::Print(Option_t *o) const
+{
     *fLog << all << GetDescriptor() << ":" << endl;
+    int id = 0;
+    TIter Next(fArray);
+    MPedestalPix *pix;
+    while ((pix=(MPedestalPix*)Next()))
+    {
+        id++;
+        if (!pix->IsValid())
+            continue;
+        *fLog << id-1 << ": ";
+        *fLog << pix->GetPedestal() << " " << pix->GetPedestalRms() << endl;
+    }
+    *fLog << "Pixels:" << endl;
+    *fLog << endl;
+    PrintArr(*fArray);
+    *fLog << endl;
+    *fLog << "Event-by-event averaged areas:" << endl;
+    *fLog << endl;
+    PrintArr(*fAverageAreas);
+    *fLog << endl;
+    *fLog << "Event-by-event averaged sectors:" << endl;
+    *fLog << endl;
+    PrintArr(*fAverageSectors);
+}

trunk/MagicSoft/Mars/mpedestal/MPedestalCam.h

r5558	r7188
24	24
25	25	UInt_t fTotalEntries; // Total number of times, the Process was executed (to estimate the error of pedestal)
	26
	27	void PrintArr(const TCollection &list) const;
26	28
27	29	public:

trunk/MagicSoft/Mars/mreport/MReportCC.cc

-              r7026
+              r7188
 // From here maily weather station data is decoded such as
 // temperature, humidity, wind-speed and solar radiation
+//
+// Class Version 2:
+// ----------------
+//  +  Float_t fUPSStatus; // arbitrary units (still not properly defined)
+//  +  Float_t fDifRubGPS; // [us] Difference between the Rubidium clock time and the time provided by the GPS receiver
+//
 //
 //////////////////////////////////////////////////////////////////////////////
 …
     Int_t len;
     const Int_t n=sscanf(str.Data(),
                          "%*f %*f %*f %*f %f %f %f %f %*f %*f %n",
+                         "%*f %*f %*f %*f %f %f %f %f %f %f %n",
                          &fTemperature, &fSolarRadiation, &fWindSpeed,
                          &fHumidity, &len);
     if (n!=4)
+                         &fHumidity, &fUPSStatus, &fDifRubGPS, &len);
+    if (n!=6)
+    {
         *fLog << warn << "WARNING - Wrong number of arguments." << endl;

trunk/MagicSoft/Mars/mreport/MReportCC.h

-              r7027
+              r7188
     Float_t fSolarRadiation; // [W/m^2] IR-Radiation
+    Float_t fUPSStatus;      // arbitrary units (still not properly defined)
+    Float_t fDifRubGPS;      // [us] Difference between the Rubidium clock time and the time provided by the GPS receiver
     Int_t InterpreteBody(TString &str, Int_t ver);
 …
     void Print(Option_t *opt) const;
     ClassDef(MReportCC, 1) // Class for CC-REPORT information
+    ClassDef(MReportCC, 2) // Class for CC-REPORT information
 };

trunk/MagicSoft/Mars/mreport/MReportCamera.cc

r6962	r7188
568	568	return kCONTINUE;
569	569
570		if (ver > gkActiveLoadControlVersNum)
	570	if (ver >= gkActiveLoadControlVersNum)
571	571	{
572	572	if (!InterpreteActiveLoad(str))

trunk/MagicSoft/Mars/msignal/MExtractTimeAndChargeSpline.cc

r7099	r7188
159	159	const Float_t MExtractTimeAndChargeSpline::fgFallTimeHiGain = 0.5;
160	160	const Float_t MExtractTimeAndChargeSpline::fgLoGainStretch = 1.5;
161		const Float_t MExtractTimeAndChargeSpline::fgOffsetLoGain = 1.7; // 5 ns
	161	const Float_t MExtractTimeAndChargeSpline::fgOffsetLoGain = 1.39; // 5 ns
162	162	const Float_t MExtractTimeAndChargeSpline::fgLoGainStartShift = -1.8;
163	163

trunk/MagicSoft/include-Classes/MMcFormat/MTriggerDefine.h

-              r5253
+              r7188
 //      In this file are the fundamental definitions for the class MCTrigger
 //
+// Number of pixels in the trigger region (used by camera simulation,
+// see camera.cxx.
 //
 #define TRIGGER_PIXELS_1      397
 #define TRIGGER_PIXELS_2      397
 #define TRIGGER_PIXELS_3      1657
+#define TRIGGER_PIXELS_4      547  // For MGeomCamMagic1183, PRELIMINARY!
 #define TRIGGER_PIXELS_5      397
 #define TRIGGER_PIXELS_6      1657

Context Navigation

Legend:

Download in other formats: