Ignore:
Timestamp:
01/16/05 12:48:33 (20 years ago)
Author:
gaug
Message:
*** empty log message ***
Location:
trunk/MagicSoft/Mars/mcalib
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • trunk/MagicSoft/Mars/mcalib/MCalibrateData.cc

    r5853 r5854  
    280280      }
    281281
     282    fCalibConsts.Reset();
     283    fCalibFFactors.Reset();
     284
    282285    return kTRUE;
    283286}
    … …  
    378381    }
    379382   
     383    const Int_t npixels = fGeomCam->GetNumPixels();
     384
     385    if (fCalibrationMode > kNone)
     386      {
     387       
     388        if (fCalibrations->GetSize() != npixels)
     389          {
     390            *fLog << err << GetDescriptor()
     391                  << ": Size mismatch between MGeomCam and MCalibrationChargeCam ... abort!" << endl;
     392            return kFALSE;
     393          }
     394       
     395        if (fBadPixels->GetSize() != npixels)
     396          {
     397            *fLog << err << GetDescriptor()
     398                  << ": Size mismatch between MGeomCam and MBadPixelsCam ... abort!" << endl;
     399            return kFALSE;
     400          }
     401       
     402        if (fBadPixels->GetSize() != npixels)
     403          {
     404            *fLog << err << GetDescriptor()
     405                  << ": Size mismatch between MGeomCam and MBadPixelsCam ... abort!" << endl;
     406            return kFALSE;
     407          }
     408      }
     409   
     410    fCalibConsts  .Set(npixels);
     411    fCalibFFactors.Set(npixels);
     412    fHiLoConv     .Set(npixels);
     413    fHiLoConvErr  .Set(npixels);
     414   
     415    if (!UpdateConversionFactors())
     416      return kFALSE;
     417
    380418    if (TestPedestalFlag(kRun))
    381419        Calibrate(kFALSE, kTRUE);
    … …  
    386424// --------------------------------------------------------------------------
    387425//
    388 // Get conversion factor and its error from MCalibrationCam
    389 //
    390 Bool_t MCalibrateData::GetConversionFactor(UInt_t pixidx, Float_t &hiloconv, Float_t &hiloconverr,
    391                                            Float_t &calibConv, Float_t &calibConvVar,
    392                                            Float_t &calibFFactor) const
     426// Update the conversion factors and F-Factors from MCalibrationCams into
     427// the arrays. Later, the here pre-calcualted conversion factors get simply
     428// copied from memory.
     429//
     430// This function can be called from outside in case that the MCalibrationCams
     431// have been updated...
     432//
     433Bool_t MCalibrateData::UpdateConversionFactors()
    393434{
    394435    //
    … …  
    397438    const Float_t zenith = -1.;
    398439
    399     hiloconv     = 1.;
    400     hiloconverr  = 0.;
    401     calibConv    = 1.;
    402     calibConvVar = 0.;
    403     calibFFactor = 0.;
    404 
    405     Float_t calibQE       = 1.;
    406     Float_t calibQEVar    = 0.;
    407 
    408     if(fCalibrationMode!=kNone)
    409     {
    410         MCalibrationChargePix &pix = (MCalibrationChargePix&)(*fCalibrations)[pixidx];
    411 
    412         hiloconv   = pix.GetConversionHiLo   ();
    413         hiloconverr= pix.GetConversionHiLoErr();
    414 
    415         if ((*fBadPixels)[pixidx].IsUnsuitable())
    416             return kFALSE;
    417 
    418         calibConv    = pix.GetMeanConvFADC2Phe();
    419         calibConvVar = pix.GetMeanConvFADC2PheVar();
    420         calibFFactor = pix.GetMeanFFactorFADC2Phot();
    421 
    422         MCalibrationQEPix &qe = (MCalibrationQEPix&) (*fQEs)[pixidx];
    423 
    424         switch(fCalibrationMode)
     440    UInt_t skip = 0;
     441
     442    for (UInt_t pixidx=0; pixidx<fGeomCam->GetNumPixels(); pixidx++)
     443      {
     444
     445        Float_t hiloconv     = 1.;
     446        Float_t hiloconverr  = 0.;
     447        Float_t calibConv    = 1.;
     448        Float_t calibConvVar = 0.;
     449        Float_t calibFFactor = 0.;
     450
     451        Float_t calibQE       = 1.;
     452        Float_t calibQEVar    = 0.;
     453
     454        if(fCalibrationMode!=kNone)
    425455          {
    426           case kFlatCharge:
    427             {
    428               MCalibrationChargePix &avpix = (MCalibrationChargePix&)fCalibrations->GetAverageArea(0);
    429               calibConv    = avpix.GetMean() / (pix.GetMean() * fGeomCam->GetPixRatio(pixidx));
    430               calibConvVar = (avpix.GetMeanRelVar() + pix.GetMeanRelVar()) * calibConv * calibConv;
    431               if (pix.IsFFactorMethodValid())
     456            MCalibrationChargePix &pix = (MCalibrationChargePix&)(*fCalibrations)[pixidx];
     457           
     458            hiloconv   = pix.GetConversionHiLo   ();
     459            hiloconverr= pix.GetConversionHiLoErr();
     460           
     461            if ((*fBadPixels)[pixidx].IsUnsuitable())
     462              {
     463                skip++;
     464                continue;
     465              }
     466           
     467            calibConv    = pix.GetMeanConvFADC2Phe();
     468            calibConvVar = pix.GetMeanConvFADC2PheVar();
     469            calibFFactor = pix.GetMeanFFactorFADC2Phot();
     470           
     471            MCalibrationQEPix &qe = (MCalibrationQEPix&) (*fQEs)[pixidx];
     472
     473            switch(fCalibrationMode)
     474              {
     475              case kFlatCharge:
    432476                {
    433                 const Float_t convmin1 = qe.GetQECascadesFFactor(zenith)/pix.GetMeanConvFADC2Phe();
    434                 if (convmin1 > 0)
    435                   calibFFactor *= TMath::Sqrt(convmin1);
    436                 else
    437                   calibFFactor = -1.;
     477                  MCalibrationChargePix &avpix = (MCalibrationChargePix&)fCalibrations->GetAverageArea(0);
     478                  calibConv    = avpix.GetMean() / (pix.GetMean() * fGeomCam->GetPixRatio(pixidx));
     479                  calibConvVar = (avpix.GetMeanRelVar() + pix.GetMeanRelVar()) * calibConv * calibConv;
     480                  if (pix.IsFFactorMethodValid())
     481                    {
     482                      const Float_t convmin1 = qe.GetQECascadesFFactor(zenith)/pix.GetMeanConvFADC2Phe();
     483                      if (convmin1 > 0)
     484                        calibFFactor *= TMath::Sqrt(convmin1);
     485                      else
     486                        calibFFactor = -1.;
     487                    }
     488                  break;
    438489                }
    439               break;
    440             }
    441           case kBlindPixel:
    442             if (!qe.IsBlindPixelMethodValid())
    443               return kFALSE;
    444             calibQE     = qe.GetQECascadesBlindPixel   ( zenith );
    445             calibQEVar  = qe.GetQECascadesBlindPixelVar( zenith );
    446             break;
    447            
    448           case kPinDiode:
    449             if (!qe.IsPINDiodeMethodValid())
    450               return kFALSE;
    451             calibQE     = qe.GetQECascadesPINDiode   ( zenith );
    452             calibQEVar  = qe.GetQECascadesPINDiodeVar( zenith );
    453             break;
    454            
    455           case kFfactor:
    456             if (!pix.IsFFactorMethodValid())
    457               return kFALSE;
    458             calibQE     = qe.GetQECascadesFFactor   ( zenith );
    459             calibQEVar  = qe.GetQECascadesFFactorVar( zenith );
    460             break;
    461            
    462           case kCombined:
    463             if (!qe.IsCombinedMethodValid())
    464               return kFALSE;
    465             calibQE     = qe.GetQECascadesCombined   ( zenith );
    466             calibQEVar  = qe.GetQECascadesCombinedVar( zenith );
    467             break;
    468            
    469           case kDummy:
    470             hiloconv    = 1.;
    471             hiloconverr = 0.;
    472             break;
    473           } /* switch calibration mode */
    474     } /* if(fCalibrationMode!=kNone) */
    475     else
    476     {
    477       calibConv  = 1./fGeomCam->GetPixRatio(pixidx);
    478     }
    479 
    480     calibConv /= calibQE;
    481 
    482     if (calibConv != 0. && calibQE != 0.)
    483     {
    484         // Now doing:
    485         // calibConvVar  = calibConvVar/(calibConv*calibConv) + calibQEVar/(calibQE*calibQE);
    486         // calibConvVar *= (calibConv*calibConv);
    487         calibConvVar += calibQEVar*(calibConv*calibConv)/(calibQE*calibQE);
     490              case kBlindPixel:
     491                if (!qe.IsBlindPixelMethodValid())
     492                  {
     493                    skip++;
     494                    continue;
     495                  }
     496                calibQE     = qe.GetQECascadesBlindPixel   ( zenith );
     497                calibQEVar  = qe.GetQECascadesBlindPixelVar( zenith );
     498                break;
     499               
     500              case kPinDiode:
     501                if (!qe.IsPINDiodeMethodValid())
     502                  {
     503                    skip++;
     504                    continue;
     505                  }
     506                calibQE     = qe.GetQECascadesPINDiode   ( zenith );
     507                calibQEVar  = qe.GetQECascadesPINDiodeVar( zenith );
     508                break;
     509               
     510              case kFfactor:
     511                if (!pix.IsFFactorMethodValid())
     512                  {
     513                    skip++;
     514                    continue;
     515                  }
     516                calibQE     = qe.GetQECascadesFFactor   ( zenith );
     517                calibQEVar  = qe.GetQECascadesFFactorVar( zenith );
     518                break;
     519               
     520              case kCombined:
     521                if (!qe.IsCombinedMethodValid())
     522                  {
     523                    skip++;
     524                    continue;
     525                  }
     526                calibQE     = qe.GetQECascadesCombined   ( zenith );
     527                calibQEVar  = qe.GetQECascadesCombinedVar( zenith );
     528                break;
     529               
     530              case kDummy:
     531                hiloconv    = 1.;
     532                hiloconverr = 0.;
     533                break;
     534              } /* switch calibration mode */
     535          } /* if(fCalibrationMode!=kNone) */
     536        else
     537          {
     538            calibConv  = 1./fGeomCam->GetPixRatio(pixidx);
     539          }
     540       
     541        calibConv /= calibQE;
     542
     543        if (calibConv != 0. && calibQE != 0.)
     544          {
     545            // Now doing:
     546            calibConvVar  = calibConvVar/(calibConv*calibConv) + calibQEVar/(calibQE*calibQE);
     547            calibConvVar *= (calibConv*calibConv);
     548            // The above two lines had been commented by TB and replaced by the following line
     549            // (without notice to me!) but it is simply wrong!
     550            //          calibConvVar += calibQEVar*(calibConv*calibConv)/(calibQE*calibQE);
     551          }
     552
     553        calibConv *= fRenormFactor;
     554
     555        fHiLoConv     [pixidx] = hiloconv;
     556        fHiLoConvErr  [pixidx] = hiloconverr;       
     557        fCalibConsts  [pixidx] = calibConv;
     558        fCalibFFactors[pixidx] = calibFFactor;
     559
     560      } /*     for (Int_t pixidx=0; pixidx<fGeomCam->GetNumPixels(); pixidx++) */
     561   
     562    if (skip>fGeomCam->GetNumPixels()*0.9)
     563    {
     564        *fLog << warn << GetDescriptor()
     565              << ": WARNING - GetConversionFactor has skipped more than 90% of the pixels... abort." << endl;
     566        return kFALSE;
    488567    }
    489568
    … …  
    491570}
    492571
     572
     573// --------------------------------------------------------------------------
     574//
     575// Apply the conversion factors and F-Factors from the arrays to the data.
     576//
     577// The flags 'data' and 'pedestal' decide whether the signal and/or the pedetals
     578// shall be calibrated, respectively.
     579//
    493580Int_t MCalibrateData::Calibrate(Bool_t data, Bool_t pedestal) const
    494581{
    … …  
    500587    const Float_t sqrtslices = TMath::Sqrt(slices);
    501588
    502     Float_t hiloconv;
    503     Float_t hiloconverr;
    504     Float_t calibConv;
    505     Float_t calibConvErr;
    506     Float_t calibFFactor;
    507 
    508     UInt_t skip = 0;
    509589    for (UInt_t pixidx=0; pixidx<npix; pixidx++)
    510590    {
    511         if (!GetConversionFactor(pixidx, hiloconv, hiloconverr,
    512                                  calibConv, calibConvErr, calibFFactor))
    513         {
    514             skip++;
    515             continue;
    516         }
    517 
    518         calibConv *= fRenormFactor;
    519 
    520         if (data)
     591     
     592      if (data)
    521593        {
    522594            const MExtractedSignalPix &sig = (*fSignals)[pixidx];
    … …  
    527599            if (sig.IsLoGainUsed())
    528600            {
    529                 signal    = sig.GetExtractedSignalLoGain()*hiloconv;
    530                 signalErr = signal*hiloconverr;
     601                signal    = sig.GetExtractedSignalLoGain()*fHiLoConv   [pixidx];
     602                signalErr = sig.GetExtractedSignalLoGain()*fHiLoConvErr[pixidx];
    531603            }
    532604            else
    … …  
    536608            }
    537609
    538             const Float_t nphot    = signal*calibConv;
    539             const Float_t nphotErr = calibFFactor*TMath::Sqrt(TMath::Abs(nphot));
    540 
    541             //
    542             // The following part is the commented first version of the error calculation
    543             // Contact Markus Gaug for questions (or wait for the next documentation update...)
    544             //
    545             /*
    546              nphotErr = signal    > 0 ? signalErr*signalErr / (signal * signal)  : 0.
    547              + calibConv > 0 ? calibConvVar  / (calibConv * calibConv ) : 0.;
    548              nphotErr  = TMath::Sqrt(nphotErr) * nphot;
    549              */
     610            const Float_t nphot    = signal                         * fCalibConsts  [pixidx];
     611            const Float_t nphotErr = TMath::Sqrt(TMath::Abs(nphot)) * fCalibFFactors[pixidx];
    550612
    551613            MCerPhotPix *cpix = fCerPhotEvt->AddPixel(pixidx, nphot, nphotErr);
    … …  
    555617
    556618            if (sig.GetNumLoGainSaturated() > 0)
    557                 cpix->SetPixelSaturated();
    558         }
     619              cpix->SetPixelSaturated();
     620        } /* if (data) */
     621     
    559622
    560623        if (pedestal)
    561624        {
    562             /*
    563             // pedestals/(used FADC slices)   in [ADC] counts
    564             const Float_t pedes  = (*fPedestalMean)[pixidx].GetPedestal()   * slices;
    565             const Float_t pedrms = (*fPedestalRms)[pixidx].GetPedestalRms() * sqrtslices;
    566 
    567             //
    568             // pedestals/(used FADC slices)   in [number of photons]
    569             //
    570             const Float_t pedphot    = pedes  * calibConv;
    571             const Float_t pedphotrms = pedrms * calibConv;
    572 
    573             (*fPedPhot)[pixidx].Set(pedphot, pedphotrms);
    574             */
    575625            TIter NextPed(&fPedestalCams);
    576626            TIter NextPhot(&fPedPhotCams);
    … …  
    578628            MPedestalCam *pedestal = 0;
    579629            MPedPhotCam  *pedphot  = 0;
     630           
     631            const Float_t pedmeancalib = slices    *fCalibConsts[pixidx];
     632            const Float_t pedrmscalib  = sqrtslices*fCalibConsts[pixidx];
    580633
    581634            while ((pedestal=(MPedestalCam*)NextPed()) &&
    582635                   (pedphot =(MPedPhotCam*)NextPhot()))
    583636            {
    584                 // pedestals/(used FADC slices)   in [ADC] counts
    585                 const Float_t pedes  = (*pedestal)[pixidx].GetPedestal()    * slices;
    586                 const Float_t pedrms = (*pedestal)[pixidx].GetPedestalRms() * sqrtslices;
    587 
    588                 // pedestals/(used FADC slices)   in [number of photons]
    589                 const Float_t mean = pedes  * calibConv;
    590                 const Float_t rms  = pedrms * calibConv;
    591 
    592                 (*pedphot)[pixidx].Set(mean, rms);
    593                 pedphot->SetReadyToSave();
     637              // pedestals/(used FADC slices)   in [number of photons]
     638              const Float_t mean = (*pedestal)[pixidx].GetPedestal()   *pedmeancalib;
     639              const Float_t rms  = (*pedestal)[pixidx].GetPedestalRms()*pedrmscalib;
     640             
     641              (*pedphot)[pixidx].Set(mean, rms);
     642              pedphot->SetReadyToSave();
    594643            }
    595         }
    596     }
    597 
    598     if (skip>npix*0.9)
    599     {
    600         *fLog << warn << "WARNING - GetConversionFactor has skipped more than 90% of the pixels... skip." << endl;
    601         return kCONTINUE;
     644        } /* if (pedestal) */
    602645    }
    603646
    … …  
    617660Int_t MCalibrateData::Process()
    618661{
    619     /*
    620      if (fCalibrations->GetNumPixels() != (UInt_t)fSignals->GetSize())
    621      {
    622      // FIXME: MExtractedSignal must be of variable size -
    623      //        like MCerPhotEvt - because we must be able
    624      //        to reduce size by zero supression
    625      //        For the moment this check could be done in ReInit...
    626      *fLog << err << "MExtractedSignal and MCalibrationCam have different sizes... abort." << endl;
    627      return kFALSE;
    628      }
    629      */
    630 
    631662    return Calibrate(fCalibrationMode!=kSkip, TestPedestalFlag(kEvent));
    632663}

  • trunk/MagicSoft/Mars/mcalib/MCalibrateData.h

    r5853 r5854  
    2121#endif
    2222
     23#ifndef MARS_MArrayF
     24#include "MArrayF.h"
     25#endif
     26
    2327class MGeomCam;
    2428class MBadPixelsCam;
    … …  
    3438{
    3539private:
     40
    3641  MGeomCam              *fGeomCam;       //! Camera geometry container
    3742  MBadPixelsCam         *fBadPixels;     //! Bad Pixels information
    … …  
    5156  TList fPedPhotCams;                    //! List of pointers to corresponding output MPedPhotCam
    5257
    53   Int_t Calibrate(Bool_t data, Bool_t pedestal) const;
     58  MArrayF fCalibConsts;                  //! Array of calibration constants for each pixel, calculated only once!
     59  MArrayF fCalibFFactors;                //! Array of calibration F-Factors for each pixel, calculated only once! 
     60  MArrayF fHiLoConv;                     //! Array of calibration constants for each pixel, calculated only once!
     61  MArrayF fHiLoConvErr;                  //! Array of calibration F-Factors for each pixel, calculated only once! 
    5462 
    55   Bool_t GetConversionFactor(UInt_t, Float_t &, Float_t &, Float_t &, Float_t &, Float_t &) const;
     63  Int_t  Calibrate(Bool_t data, Bool_t pedestal) const;
    5664 
    5765  Int_t  PreProcess(MParList *pList);
    … …  
    93101                 const char *name=NULL, const char *title=NULL);
    94102 
     103  void   AddPedestal(const char *name="Cam");
     104  void   AddPedestal(const char *pedestal, const char *pedphot);
     105
    95106  void   EnablePedestalType(PedestalType_t i)     { fPedestalFlag |=  i;      }
    96107  void   SetPedestalFlag(PedestalType_t i=kRun)   { fPedestalFlag  =  i;      }
    … …  
    100111  void   SetSignalType      ( SignalType_t      sigtype=kPhot    ) { fSignalType     =sigtype; } 
    101112
    102   void   AddPedestal(const char *name="Cam");
    103   void   AddPedestal(const char *pedestal, const char *pedphot);
    104 
     113  Bool_t UpdateConversionFactors();
     114 
    105115  ClassDef(MCalibrateData, 1)   // Task to calibrate FADC counts into Cherenkov photons
    106116};
Note: See TracChangeset for help on using the changeset viewer.