Changeset 4752


Ignore:
Timestamp:
08/27/04 11:36:24 (20 years ago)
Author:
tbretz
Message:
*** empty log message ***
Location:
trunk/MagicSoft/Mars
Files:
2 added
12 edited

Legend:

Unmodified
Added
Removed

  • trunk/MagicSoft/Mars/Changelog

    r4748 r4752  
    1919
    2020                                                 -*-*- END OF LINE -*-*-
     21 2004/08/27: Thomas Bretz
     22
     23   * Makefile:
     24     - added comments how to link statically
     25
     26   * callisto.cc:
     27     - fixed some output
     28
     29   * mbadpixels/Makefile:
     30     - added a comment
     31
     32   * mbase/BaseLinkDef.h, mbase/Makefile:
     33     - added MArrayI
     34
     35   * mbase/MArrayI.[h,cc]:
     36     - added
     37
     38   * mbase/MArrayD.cc:
     39     - fixed some comments
     40
     41   * mcalib/MCalibrateData.[h,cc]:
     42     - unified CalibratePedestal and CalibrateData. Calling GetConvFactor twice
     43       took a lot of time.
     44
     45   * mjobs/MJCalibrateSignal.cc, mjobs/MJPedestal.cc:
     46     - added two empty lines to output if finished
     47
     48   * mpedestal/MPedPhotCam.cc:
     49     - use faster MArrays in ReCalc
     50     - accelerated GetPixelContent
     51
     52   * msignal/MExtractTimeFastSpline.cc:
     53     - accelerated a bit by defining
     54          Float_t higainklo = fHiGainSecondDeriv[klo];
     55          Float_t higainkhi = fHiGainSecondDeriv[khi];
     56       instead of accesing the arrays many times inside the loops.
     57       Somebody should do the same for logain.
     58
     59
     60
    2161 2004/08/26: Wolfgang Wittek
    2262

  • trunk/MagicSoft/Mars/Makefile

    r4722 r4752  
    106106        @echo " Linking $@ ..."
    107107        $(CXX) $(CXXFLAGS) $(ROOTGLIBS) $(SOLIB) $@.o $(MARS_LIB) -o $@
     108
     109# Use this to link the programs statically - for gprof
     110#$(PROGRAMS): $(OBJS) $(HEADERS) $(PROGRAMS:=.o)
     111#       @echo " Linking $@ ..."
     112#       $(CXX) $(CXXFLAGS) $(ROOTGLIBS) $(OBJS) $(SUBDIRS:=/*.o) $@.o $(MARS_LIB) -o $@
    108113else
    109114$(DYLIB): $(LIBRARIES) $(OBJS) $(HEADERS)

  • trunk/MagicSoft/Mars/callisto.cc

    r4748 r4752  
    325325            return 1;
    326326        }
    327 
    328         gLog << flush << all << endl << endl;
    329327    }
    330328
    … …  
    380378            return 1;
    381379        }
    382 
    383         gLog << flush << all << endl << endl;
    384380    }
    385381

  • trunk/MagicSoft/Mars/mbase/BaseLinkDef.h

    r4747 r4752  
    6060#pragma link C++ class MArrayS;
    6161#pragma link C++ class MArrayD;
     62#pragma link C++ class MArrayI;
    6263
    6364#pragma link C++ class MTime+;

  • trunk/MagicSoft/Mars/mbase/MArrayD.cc

    r4747 r4752  
    2626//////////////////////////////////////////////////////////////////////////////
    2727//
    28 // MArrayS
     28// MArrayD
    2929//
    30 // Array of UShort_t. It is almost the same than TArrayD but derives from
     30// Array of Double_t. It is almost the same than TArrayD but derives from
    3131// TObject
    3232//

  • trunk/MagicSoft/Mars/mbase/Makefile

    r4747 r4752  
    5252           MArrayS.cc \
    5353           MArrayD.cc \
     54           MArrayI.cc \
    5455           MTime.cc \
    5556           MClone.cc \

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

    r4639 r4752  
    265265    }
    266266
     267    // Sizes might have changed
     268    if (fPedestalFlag && (Int_t)fPedestal->GetSize() != fSignals->GetSize())
     269    {
     270        *fLog << err << "Size mismatch of MPedestalCam and MCalibrationCam... abort." << endl;
     271        return kFALSE;
     272    }
     273
    267274    if(fCalibrationMode == kBlindPixel && !fQEs->IsBlindPixelMethodValid())
    268275    {
    … …  
    315322
    316323    if (TestPedestalFlag(kRun))
    317         if (!CalibratePedestal())
    318           return kFALSE;
     324        Calibrate(kFALSE, kTRUE);
    319325
    320326    return kTRUE;
    … …  
    323329// --------------------------------------------------------------------------
    324330//
    325 // Calibrate the Pedestal values
    326 //
    327 Bool_t MCalibrateData::CalibratePedestal()
    328 {
    329     //
    330     // Fill MPedPhot container using the informations from
    331     // MPedestalCam, MExtractedSignalCam and MCalibrationCam
    332     //
    333     const Float_t slices     = fSignals->GetNumUsedHiGainFADCSlices();
    334     const Float_t sqrtslices = TMath::Sqrt(slices);
    335 
    336     // is pixid equal to pixidx ?
    337     if ((Int_t)fPedestal->GetSize() != fSignals->GetSize())
    338     {
    339         *fLog << err << "Sizes of MPedestalCam and MCalibrationCam are different" << endl;
    340         return kFALSE;
    341     }
    342 
    343     const Int_t n = fPedestal->GetSize();
    344     for (Int_t pixid=0; pixid<n; pixid++)
    345     {
    346         const MPedestalPix &ped = (*fPedestal)[pixid];
    347 
    348         // pedestals/(used FADC slices)   in [ADC] counts
    349         const Float_t pedes  = ped.GetPedestal()    * slices;
    350         const Float_t pedrms = ped.GetPedestalRms() * sqrtslices;
    351 
    352         //
    353         // get phe/ADC conversion factor
    354         //
    355         Float_t hiloconv;
    356         Float_t hiloconverr;
    357         Float_t calibConv;
    358         Float_t calibConvVar;
    359         Float_t calibFFactor;
    360         if (!GetConversionFactor(pixid, hiloconv, hiloconverr,
    361                                  calibConv, calibConvVar, calibFFactor))
    362             continue;
    363 
    364         //
    365         // pedestals/(used FADC slices)   in [number of photons]
    366         //
    367         const Float_t pedphot    = pedes  * calibConv;
    368         const Float_t pedphotrms = pedrms * calibConv;
    369 
    370         (*fPedPhot)[pixid].Set(pedphot, pedphotrms);
    371     }
    372 
    373     fPedPhot->SetReadyToSave();
    374 
    375     return kTRUE;
    376 }
    377 
    378 // --------------------------------------------------------------------------
    379 //
    380331// Get conversion factor and its error from MCalibrationCam
    381332//
    382333Bool_t MCalibrateData::GetConversionFactor(UInt_t pixidx, Float_t &hiloconv, Float_t &hiloconverr,
    383                                            Float_t &calibConv, Float_t &calibConvVar, Float_t &calibFFactor)
     334                                           Float_t &calibConv, Float_t &calibConvVar,
     335                                           Float_t &calibFFactor) const
    384336{
    385337    //
    … …  
    486438}
    487439
    488 void MCalibrateData::CalibrateData()
    489 {
    490     const UInt_t npix = fSignals->GetSize();
     440void MCalibrateData::Calibrate(Bool_t data, Bool_t pedestal) const
     441{
     442    if (!data && !pedestal)
     443        return;
     444
     445    const UInt_t  npix       = fSignals->GetSize();
     446    const Float_t slices     = fSignals->GetNumUsedHiGainFADCSlices();
     447    const Float_t sqrtslices = TMath::Sqrt(slices);
    491448
    492449    Float_t hiloconv;
    493450    Float_t hiloconverr;
    494     Float_t calibrationConversionFactor;
    495     Float_t calibrationConversionFactorErr;
     451    Float_t calibConv;
     452    Float_t calibConvErr;
    496453    Float_t calibFFactor;
    497454
    … …  
    499456    {
    500457        if (!GetConversionFactor(pixidx, hiloconv, hiloconverr,
    501                                  calibrationConversionFactor, calibrationConversionFactorErr, calibFFactor))
     458                                 calibConv, calibConvErr, calibFFactor))
    502459            continue;
    503460
    504         const MExtractedSignalPix &sig = (*fSignals)[pixidx];
    505 
    506         Float_t signal = 0;
    507         Float_t signalErr = 0.;
    508 
    509         if (sig.IsLoGainUsed())
    510         {
    511             signal    = sig.GetExtractedSignalLoGain()*hiloconv;
    512             signalErr = signal*hiloconverr;
    513         }
    514         else
    515         {
    516             if (sig.GetExtractedSignalHiGain() <= 9999.)
    517                 signal = sig.GetExtractedSignalHiGain();
    518         }
    519 
    520         const Float_t nphot    = signal*calibrationConversionFactor;
    521         const Float_t nphotErr = calibFFactor*TMath::Sqrt(TMath::Abs(nphot));
    522 
    523         //
    524         // The following part is the commented first version of the error calculation
    525         // Contact Markus Gaug for questions (or wait for the next documentation update...)
    526         //
    527         /*
    528          nphotErr = signal    > 0 ? signalErr*signalErr / (signal * signal)  : 0.
    529          + calibConv > 0 ? calibConvVar  / (calibConv * calibConv ) : 0.;
    530          nphotErr  = TMath::Sqrt(nphotErr) * nphot;
    531          */
    532 
    533         MCerPhotPix *cpix = fCerPhotEvt->AddPixel(pixidx, nphot, nphotErr);
    534 
    535         if (sig.GetNumHiGainSaturated() > 0)
    536             cpix->SetPixelHGSaturated();
    537 
    538         if (sig.GetNumLoGainSaturated() > 0)
    539             cpix->SetPixelSaturated();
    540     }
    541 
    542     fCerPhotEvt->FixSize();
    543     fCerPhotEvt->SetReadyToSave();
     461        if (data)
     462        {
     463            const MExtractedSignalPix &sig = (*fSignals)[pixidx];
     464
     465            Float_t signal = 0;
     466            Float_t signalErr = 0.;
     467
     468            if (sig.IsLoGainUsed())
     469            {
     470                signal    = sig.GetExtractedSignalLoGain()*hiloconv;
     471                signalErr = signal*hiloconverr;
     472            }
     473            else
     474            {
     475                if (sig.GetExtractedSignalHiGain() <= 9999.)
     476                    signal = sig.GetExtractedSignalHiGain();
     477            }
     478
     479            const Float_t nphot    = signal*calibConv;
     480            const Float_t nphotErr = calibFFactor*TMath::Sqrt(TMath::Abs(nphot));
     481
     482            //
     483            // The following part is the commented first version of the error calculation
     484            // Contact Markus Gaug for questions (or wait for the next documentation update...)
     485            //
     486            /*
     487             nphotErr = signal    > 0 ? signalErr*signalErr / (signal * signal)  : 0.
     488             + calibConv > 0 ? calibConvVar  / (calibConv * calibConv ) : 0.;
     489             nphotErr  = TMath::Sqrt(nphotErr) * nphot;
     490             */
     491
     492            MCerPhotPix *cpix = fCerPhotEvt->AddPixel(pixidx, nphot, nphotErr);
     493
     494            if (sig.GetNumHiGainSaturated() > 0)
     495                cpix->SetPixelHGSaturated();
     496
     497            if (sig.GetNumLoGainSaturated() > 0)
     498                cpix->SetPixelSaturated();
     499        }
     500        if (pedestal)
     501        {
     502            const MPedestalPix &ped = (*fPedestal)[pixidx];
     503
     504            // pedestals/(used FADC slices)   in [ADC] counts
     505            const Float_t pedes  = ped.GetPedestal()    * slices;
     506            const Float_t pedrms = ped.GetPedestalRms() * sqrtslices;
     507
     508            //
     509            // pedestals/(used FADC slices)   in [number of photons]
     510            //
     511            const Float_t pedphot    = pedes  * calibConv;
     512            const Float_t pedphotrms = pedrms * calibConv;
     513
     514            (*fPedPhot)[pixidx].Set(pedphot, pedphotrms);
     515        }
     516    }
     517
     518    if (pedestal)
     519        fPedPhot->SetReadyToSave();
     520
     521    if (data)
     522    {
     523        fCerPhotEvt->FixSize();
     524        fCerPhotEvt->SetReadyToSave();
     525    }
    544526}
    545527
    … …  
    563545     */
    564546
    565     if (fCalibrationMode!=kSkip)
    566         CalibrateData();
    567 
    568     if (TestPedestalFlag(kEvent))
    569         if (!CalibratePedestal())
    570             return kFALSE;
    571 
    572   return kTRUE;
     547    Calibrate(fCalibrationMode!=kSkip, TestPedestalFlag(kEvent));
     548    return kTRUE;
    573549}
    574550

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

    r4632 r4752  
    5353  TString  fNamePedPhotContainer;           // name of fPedPhot
    5454
    55   Bool_t CalibratePedestal();
    56   void   CalibrateData();
     55  void Calibrate(Bool_t data, Bool_t pedestal) const;
    5756 
    58   Bool_t GetConversionFactor(UInt_t, Float_t &, Float_t &, Float_t &, Float_t &, Float_t &);
     57  Bool_t GetConversionFactor(UInt_t, Float_t &, Float_t &, Float_t &, Float_t &, Float_t &) const;
    5958 
    6059  Int_t  PreProcess(MParList *pList);

  • trunk/MagicSoft/Mars/mjobs/MJCalibrateSignal.cc

    r4733 r4752  
    355355
    356356    *fLog << all << GetDescriptor() << ": Done." << endl;
     357    *fLog << endl << endl;
    357358
    358359    return kTRUE;

  • trunk/MagicSoft/Mars/mjobs/MJPedestal.cc

    r4733 r4752  
    650650
    651651    *fLog << all << GetDescriptor() << ": Done." << endl;
     652    *fLog << endl << endl;
    652653
    653654    return kTRUE;

  • trunk/MagicSoft/Mars/mpedestal/MPedPhotCam.cc

    r4609 r4752  
    3939#include "MPedPhotCam.h"
    4040
    41 #include <TArrayI.h>
    42 #include <TArrayD.h>
    4341#include <TClonesArray.h>
     42
     43#include "MArrayI.h"
     44#include "MArrayD.h"
    4445
    4546#include "MLog.h"
    … …  
    233234    const Int_t na = GetNumAreas();
    234235
    235     TArrayI acnt(na);
    236     TArrayI scnt(ns);
    237     TArrayD asumx(na);
    238     TArrayD ssumx(ns);
    239     TArrayD asumr(na);
    240     TArrayD ssumr(ns);
     236    // Using MArray instead of TArray because they don't do range checks
     237    MArrayI acnt(na);
     238    MArrayI scnt(ns);
     239    MArrayD asumx(na);
     240    MArrayD ssumx(ns);
     241    MArrayD asumr(na);
     242    MArrayD ssumr(ns);
    241243
    242244    for (int i=0; i<np; i++)
    … …  
    246248
    247249        // Create sums for areas and sectors
     250        const MPedPhotPix &pix = (*this)[i];
     251
     252        const UInt_t  ne   = pix.GetNumEvents();
     253        const Float_t mean = ne*pix.GetMean();
     254        const Float_t rms  = ne*pix.GetRms();
     255
    248256        const UInt_t aidx = geom[i].GetAidx();
     257        asumx[aidx] += mean;
     258        asumr[aidx] += rms;
     259        acnt[aidx]  += ne;
     260
    249261        const UInt_t sect = geom[i].GetSector();
    250 
    251         const MPedPhotPix &pix = (*this)[i];
    252 
    253         const UInt_t  ne   = pix.GetNumEvents();
    254         const Float_t mean = pix.GetMean();
    255         const Float_t rms  = pix.GetRms();
    256 
    257         asumx[aidx] += ne*mean;
    258         asumr[aidx] += ne*rms;
    259         acnt[aidx]  += ne;
    260 
    261         ssumx[sect] += ne*mean;
    262         ssumr[sect] += ne*rms;
     262        ssumx[sect] += mean;
     263        ssumr[sect] += rms;
    263264        scnt[sect]  += ne;
    264265    }
    … …  
    306307Bool_t MPedPhotCam::GetPixelContent(Double_t &val, Int_t idx, const MGeomCam &cam, Int_t type) const
    307308{
    308 
    309     const Float_t ped      = (*this)[idx].GetMean();
    310     const Float_t rms      = (*this)[idx].GetRms();
    311     const UInt_t  nevt     = (*this)[idx].GetNumEvents();
    312 
    313 
    314     Float_t pederr;
    315     Float_t rmserr;
    316 
    317     if (nevt>0)
    318     {
    319         pederr = rms/TMath::Sqrt((Float_t)nevt);
    320         rmserr = rms/TMath::Sqrt((Float_t)nevt)/2.;
    321     }
    322     else
    323     {
    324         pederr = -1;
    325         rmserr = -1;
    326     }
    327 
    328309    switch (type)
    329310    {
    330311    case 0:
    331         val = ped;// (*this)[idx].GetMean();
     312        val = (*this)[idx].GetMean();
    332313        break;
    333314    case 1:
    334         val = rms; // (*this)[idx].GetRms();
     315        val = (*this)[idx].GetRms();
    335316        break;
    336317    case 2:
    337         val = pederr; // new
     318        val = (*this)[idx].GetNumEvents()>0 ? (*this)[idx].GetRms()/TMath::Sqrt((Float_t)(*this)[idx].GetNumEvents()) : -1;
    338319        break;
    339320    case 3:
    340       val = rmserr; // new
     321        val = (*this)[idx].GetNumEvents()>0 ? (*this)[idx].GetRms()/TMath::Sqrt((Float_t)(*this)[idx].GetNumEvents())/2. : -1;
    341322        break;
    342323    default:

  • trunk/MagicSoft/Mars/msignal/MExtractTimeFastSpline.cc

    r4732 r4752  
    200200      pp = fHiGainSecondDeriv[i-1] + 4.;
    201201      fHiGainSecondDeriv[i] = -1.0/pp;
    202       fHiGainFirstDeriv [i] = *(p+1) - 2.* *(p) + *(p-1);
    203       fHiGainFirstDeriv [i] = (6.0*fHiGainFirstDeriv[i]-fHiGainFirstDeriv[i-1])/pp;
     202      const Double_t deriv = *(p+1) - 2.* *(p) + *(p-1);
     203      fHiGainFirstDeriv [i] = (6.0*deriv-fHiGainFirstDeriv[i-1])/pp;
    204204    }
    205205
    … …  
    232232  // interval maxpos+1.
    233233  //
     234
     235  Float_t higainklo = fHiGainSecondDeriv[klo];
     236  Float_t higainkhi = fHiGainSecondDeriv[khi];
    234237  while (x<upper-0.3)
    235238    {
    … …  
    241244      y = a*klocont
    242245        + b*khicont
    243         + (a*a*a-a)*fHiGainSecondDeriv[klo]
    244         + (b*b*b-b)*fHiGainSecondDeriv[khi];
     246        + (a*a*a-a)*higainklo
     247        + (b*b*b-b)*higainkhi;
    245248     
    246249      if (y > abmax)
    … …  
    265268      khicont = (Float_t)*(first+khi);
    266269
     270      higainklo = fHiGainSecondDeriv[klo];
     271      higainkhi = fHiGainSecondDeriv[khi];
    267272      while (x<upper-0.3)
    268273        {
    … …  
    274279          y = a* klocont
    275280            + b* khicont
    276             + (a*a*a-a)*fHiGainSecondDeriv[klo]
    277             + (b*b*b-b)*fHiGainSecondDeriv[khi];
     281            + (a*a*a-a)*higainklo
     282            + (b*b*b-b)*higainkhi;
    278283         
    279284          if (y > abmax)
    … …  
    295300  step  = 0.04; // step size of 83 ps
    296301
     302  higainklo = fHiGainSecondDeriv[klo];
     303  higainkhi = fHiGainSecondDeriv[khi];
    297304  while (x<up)
    298305    {
    … …  
    304311      y = a* klocont
    305312        + b* khicont
    306         + (a*a*a-a)*fHiGainSecondDeriv[klo]
    307         + (b*b*b-b)*fHiGainSecondDeriv[khi];
     313        + (a*a*a-a)*higainklo
     314        + (b*b*b-b)*higainkhi;
    308315     
    309316      if (y > abmax)
    … …  
    329336  b     = x - lower;
    330337
     338  higainklo = fHiGainSecondDeriv[klo];
     339  higainkhi = fHiGainSecondDeriv[khi];
    331340  while (x>lo)
    332341    {
    … …  
    338347      y = a* klocont
    339348        + b* khicont
    340         + (a*a*a-a)*fHiGainSecondDeriv[klo]
    341         + (b*b*b-b)*fHiGainSecondDeriv[khi];
     349        + (a*a*a-a)*higainklo
     350        + (b*b*b-b)*higainkhi;
    342351     
    343352      if (y > abmax)
Note: See TracChangeset for help on using the changeset viewer.