Ignore:
Timestamp:
02/16/04 11:18:25 (21 years ago)
Author:
tbretz
Message:
*** empty log message ***
Location:
trunk/MagicSoft/Mars/mcalib
Files:
3 edited

Legend:

Unmodified
Added
Removed

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

    r3116 r3183  
    1717!
    1818!   Author(s): Javier Lopez 12/2003 <mailto:jlopez@ifae.es>
    19 !   Modified by: Javier Rico  01/2004 <mailto:jrico@ifae.es>
     19!   Author(s): Javier Rico  01/2004 <mailto:jrico@ifae.es>
    2020!
    2121!   Copyright: MAGIC Software Development, 2000-2004

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

    r3140 r3183  
    1717!
    1818!   Author(s): Javier Lopez    12/2003 <mailto:jlopez@ifae.es>
    19 !   Modified by: Javier Rico   01/2004 <mailto:jrico@ifae.es>
    20 !             Wolfgang Wittek 02/2004 <mailto:wittek@mppmu.mpg.de>
     19!   Author(s): Javier Rico     01/2004 <mailto:jrico@ifae.es>
     20!   Author(s): Wolfgang Wittek 02/2004 <mailto:wittek@mppmu.mpg.de>
    2121!
    2222!   Copyright: MAGIC Software Development, 2000-2004

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

    r3004 r3183  
    1818!   Author(s): Abelardo Moralejo, 12/2003 <mailto:moralejo@pd.infn.it>
    1919!
    20 !   Copyright: MAGIC Software Development, 2000-2003
     20!   Copyright: MAGIC Software Development, 2000-2004
    2121!
    2222!
    … …  
    6565MMcCalibrationCalc::MMcCalibrationCalc(const char *name, const char *title)
    6666{
    67   fName  = name  ? name  : "MMcCalibrationCalc";
    68   fTitle = title ? title : "Calculate and write conversion factors into MCalibrationCam Container";
    69 
    70   fADC2Phot = 0.;
    71   fEvents = 0;
    72 
    73   fHistRatio = new TH1F(AddSerialNumber("HistRatio"), "log10(fPassPhotCone/fSize)", 1500, -3., 3.);
    74   fHistRatio->GetXaxis()->SetTitle("log_{10}(fPassPhotCone / fSize) (in photons/ADC count)");
     67    fName  = name  ? name  : "MMcCalibrationCalc";
     68    fTitle = title ? title : "Calculate and write conversion factors into MCalibrationCam Container";
     69
     70    fHistRatio = new TH1F(AddSerialNumber("HistRatio"), "log10(fPassPhotCone/fSize)", 1500, -3., 3.);
     71    fHistRatio->SetXTitle("log_{10}(fPassPhotCone / fSize) [phot/ADC count]");
    7572}
    7673
    … …  
    8380Bool_t MMcCalibrationCalc::CheckRunType(MParList *pList) const
    8481{
    85   const MRawRunHeader *run = (MRawRunHeader*)pList->FindObject("MRawRunHeader");
    86   if (!run)
    87     {
    88       *fLog << warn << dbginf << "Warning - cannot check file type, MRawRunHeader not found." << endl;
    89       return kTRUE;
    90     }
    91 
    92   return  run->GetRunType() == kRTMonteCarlo;
     82    const MRawRunHeader *run = (MRawRunHeader*)pList->FindObject("MRawRunHeader");
     83    if (!run)
     84    {
     85        *fLog << warn << "Warning - cannot check file type, MRawRunHeader not found." << endl;
     86        return kTRUE;
     87    }
     88
     89    return  run->GetRunType() == kRTMonteCarlo;
    9390}
    9491
    … …  
    9996Int_t MMcCalibrationCalc::PreProcess(MParList *pList)
    10097{
    101 
    102   fCalCam = (MCalibrationCam*) pList->FindObject(AddSerialNumber("MCalibrationCam"));
    103 
    104   if ( !fCalCam )
    105     {
    106       *fLog << err << dbginf << "Cannot find " << AddSerialNumber("MCalibrationCam") << "... aborting." << endl;
    107       return kFALSE;
    108     }
    109 
    110   fHillas = (MHillas*) pList->FindCreateObj(AddSerialNumber("MHillas"));
    111   if ( !fHillas)
    112     {
    113       *fLog << err << dbginf << "Cannot find " << AddSerialNumber("MHillas") << "... aborting." << endl;
    114       return kFALSE;
    115     }
    116 
    117   fNew = (MNewImagePar*) pList->FindCreateObj(AddSerialNumber("MNewImagePar"));
    118   if ( !fNew)
    119     {
    120       *fLog << err << dbginf << "Cannot find " << AddSerialNumber("MNewImagePar") << "... aborting." << endl;
    121       return kFALSE;
    122     }
    123 
    124   fMcEvt = (MMcEvt*) pList->FindCreateObj(AddSerialNumber("MMcEvt"));
    125   if ( !fMcEvt)
    126     {
    127       *fLog << err << dbginf << "Cannot find " << AddSerialNumber("MMcEvt") << "... aborting." << endl;
    128       return kFALSE;
    129     }
    130 
    131   return kTRUE;
    132 
     98    fHistRatio->Reset();
     99    fADC2Phot = 0;
     100
     101    fCalCam = (MCalibrationCam*) pList->FindObject(AddSerialNumber("MCalibrationCam"));
     102    if (!fCalCam)
     103    {
     104        *fLog << err << AddSerialNumber("MCalibrationCam") << "not found... aborting." << endl;
     105        return kFALSE;
     106    }
     107
     108    fHillas = (MHillas*) pList->FindObject(AddSerialNumber("MHillas"));
     109    if ( !fHillas)
     110    {
     111        *fLog << err << AddSerialNumber("MHillas") << "not found... aborting." << endl;
     112        return kFALSE;
     113    }
     114
     115    fNew = (MNewImagePar*)pList->FindObject(AddSerialNumber("MNewImagePar"));
     116    if (!fNew)
     117    {
     118        *fLog << err << AddSerialNumber("MNewImagePar") << "not found... aborting." << endl;
     119        return kFALSE;
     120    }
     121
     122    fMcEvt = (MMcEvt*) pList->FindObject(AddSerialNumber("MMcEvt"));
     123    if (!fMcEvt)
     124    {
     125        *fLog << err << AddSerialNumber("MMcEvt") << "not found... aborting." << endl;
     126        return kFALSE;
     127    }
     128
     129    return kTRUE;
    133130}
    134131
    … …  
    144141  //
    145142  if (!CheckRunType(pList))
    146     {
    147       *fLog << err << dbginf << "This is no MC file... aborting." << endl;
     143  {
     144      *fLog << err << "MMcCalibrationCalc can only used with MC files... aborting." << endl;
    148145      return kFALSE;
    149     }
     146  }
    150147
    151148  //
    152149  // Now check the existence of all necessary containers.
    153150  //
    154 
    155151  fGeom = (MGeomCam*) pList->FindObject(AddSerialNumber("MGeomCam"));
    156   if ( ! fGeom )
    157     {
    158       *fLog << err << dbginf << "Cannot find " << AddSerialNumber("MGeomCam") << "... aborting." << endl;
     152  if (!fGeom)
     153  {
     154      *fLog << err << AddSerialNumber("MGeomCam") << " mot found... aborting." << endl;
    159155      return kFALSE;
    160     }
     156  }
    161157
    162158  fHeaderFadc = (MMcFadcHeader*)pList->FindObject(AddSerialNumber("MMcFadcHeader"));
    163159  if (!fHeaderFadc)
    164     {
    165       *fLog << err << dbginf << AddSerialNumber("MMcFadcHeader") << " not found... aborting." << endl;
     160  {
     161      *fLog << err << AddSerialNumber("MMcFadcHeader") << " not found... aborting." << endl;
    166162      return kFALSE;
    167     }
     163  }
    168164
    169165  for (UInt_t ipix = 0; ipix < fGeom->GetNumPixels(); ipix++)
    170     {
     166  {
    171167      if (fHeaderFadc->GetPedestalRmsHigh(ipix) > 0 ||
    172           fHeaderFadc->GetPedestalRmsLow(ipix) > 0 )
    173         {
    174           *fLog << err << endl << endl << dbginf << "You are trying to calibrate the data using a Camera file produced with added noise. Please use a noiseless file for calibration. Aborting..." << endl << endl;
     168          fHeaderFadc->GetPedestalRmsLow(ipix)  > 0 )
     169      {
     170          *fLog << err << "Trying to calibrate the data using a Camera file produced with added noise." << endl;
     171          *fLog << "Please use a noiseless file for calibration... aborting." << endl << endl;
    175172          return kFALSE;
    176         }
    177     }
     173      }
     174  }
    178175
    179176  return kTRUE;
    … …  
    187184Int_t MMcCalibrationCalc::Process()
    188185{
    189 
    190   //
    191   // Exclude events with some saturated pixel
    192   //
    193   if ( fNew->GetNumSaturatedPixels() > 0 )
     186    //
     187    // Exclude events with some saturated pixel
     188    //
     189    if (fNew->GetNumSaturatedPixels()>0)
     190        return kTRUE;
     191
     192    //
     193    // Exclude events with low Size (larger fluctuations)
     194    // FIXME? The present cut (1000 "inner-pixel-counts") is somehow
     195    // arbitrary. Might it be optimized?
     196    //
     197    if (fHillas->GetSize()<1000)
     198        return kTRUE;
     199
     200    fADC2Phot += fMcEvt->GetPassPhotCone()/fHillas->GetSize();
     201
     202    fHistRatio->Fill(TMath::Log10(fMcEvt->GetPassPhotCone()/fHillas->GetSize()));
     203
    194204    return kTRUE;
    195 
    196   //
    197   // Exclude events with low Size (larger fluctuations)
    198   // FIXME? The present cut (1000 "inner-pixel-counts") is somehow arbitrary.
    199   // Might it be optimized?
    200   //
    201   if ( fHillas->GetSize() < 1000 )
     205}
     206
     207// --------------------------------------------------------------------------
     208//
     209// Fill the MCalibrationCam object
     210//
     211Int_t MMcCalibrationCalc::PostProcess()
     212{
     213    const Stat_t n = fHistRatio->GetEntries();
     214    if (n<1)
     215    {
     216        *fLog << err << "No events read... aborting." << endl;
     217        return kFALSE;
     218    }
     219
     220    fADC2Phot /= n;
     221
     222    //
     223    // For the calibration we no longer use the mean,
     224    // but the peak of the distribution:
     225    //
     226    const Int_t reach = 2;
     227
     228    Stat_t summax = 0;
     229    Int_t  mode   = 0;
     230
     231    // FIXME: Is this necessary? We could use GetMaximumBin instead..
     232    for (Int_t ibin = 1+reach; ibin <= fHistRatio->GetNbinsX()-reach; ibin++)
     233    {
     234        const Stat_t sum = fHistRatio->Integral(ibin-reach, ibin+reach);
     235
     236        if (sum <= summax)
     237            continue;
     238
     239        summax = sum;
     240        mode = ibin;
     241    }
     242
     243    fADC2Phot = TMath::Power(10, fHistRatio->GetBinCenter(mode));
     244
     245    const Int_t num = fCalCam->GetSize();
     246    for (int i=0; i<num; i++)
     247    {
     248        MCalibrationPix &calpix = (*fCalCam)[i];
     249
     250        const Float_t factor = fADC2Phot*calpix.GetMeanConversionBlindPixelMethod();
     251
     252        calpix.SetConversionBlindPixelMethod(factor, 0., 0.);
     253    }
     254
    202255    return kTRUE;
    203 
    204   fADC2Phot += fMcEvt->GetPassPhotCone() / fHillas->GetSize();
    205   fEvents ++;
    206 
    207   fHistRatio->Fill(log10(fMcEvt->GetPassPhotCone()/fHillas->GetSize()));
    208 
    209   return kTRUE;
    210 }
    211 
    212 // --------------------------------------------------------------------------
    213 //
    214 // Fill the MCalibrationCam object
    215 //
    216 Int_t MMcCalibrationCalc::PostProcess()
    217 {
    218 
    219   if (fEvents > 0)
    220     fADC2Phot /= fEvents;
    221   else
    222     {
    223       *fLog << err << dbginf << "No events were read! Aborting." << endl;
    224       return kFALSE;
    225     }
    226 
    227   //
    228   // For the calibration we no longer use the mean, but thepeak of the distribution:
    229   //
    230 
    231   Float_t summax = 0.;
    232   Int_t mode = 0;
    233   Int_t reach = 2;
    234   for (Int_t ibin = 1+reach; ibin <= fHistRatio->GetNbinsX()-reach; ibin++)
    235     {
    236       Float_t sum = 0;
    237       for(Int_t k = ibin-reach; k <= ibin+reach; k++)
    238         sum += fHistRatio->GetBinContent(k);
    239       if (sum > summax)
    240         {
    241           summax = sum;
    242           mode = ibin;
    243         }
    244     }
    245 
    246   fADC2Phot = pow(10., fHistRatio->GetXaxis()->GetBinCenter(mode));
    247 
    248   const int num = fCalCam->GetSize();
    249  
    250   for (int i=0; i<num; i++)
    251     {
    252       MCalibrationPix &calpix = (*fCalCam)[i];
    253 
    254       Float_t factor =  fADC2Phot*calpix.GetMeanConversionBlindPixelMethod();
    255 
    256       calpix.SetConversionBlindPixelMethod(factor, 0., 0.);
    257      
    258     }
    259 
    260   return kTRUE;
    261 
    262 }
     256}
Note: See TracChangeset for help on using the changeset viewer.