Ignore:
Timestamp:
10/25/06 19:36:26 (18 years ago)
Author:
tbretz
Message:
*** empty log message ***
Location:
trunk/MagicSoft/Mars/msignal
Files:
6 edited

Legend:

Unmodified
Added
Removed

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

    r8154 r8165  
    11/* ======================================================================== *\
    2 ! $Name: not supported by cvs2svn $:$Id: MExtractTimeAndCharge.cc,v 1.55 2006-10-24 08:24:52 tbretz Exp $
     2! $Name: not supported by cvs2svn $:$Id: MExtractTimeAndCharge.cc,v 1.56 2006-10-25 18:36:26 tbretz Exp $
    33! --------------------------------------------------------------------------
    44!
    … …  
    6969
    7070#include <TRandom.h>
     71#include <TVector3.h>
     72
     73#include "MMath.h"
    7174
    7275#include "MLog.h"
    … …  
    7578#include "MParList.h"
    7679
    77 //#include "MArrayB.h"
    78 //#include "MArrayF.h"
    79 
    80 //#include "MRawEvtData.h"
    8180#include "MRawEvtPixelIter.h"
    82 //#include "MRawRunHeader.h"
    83 
    84 //#include "MPedestalCam.h"
    85 //#include "MPedestalPix.h"
    86 #include "MPedestalSubtractedEvt.h"
     81#include "MRawRunHeader.h"
    8782
    8883#include "MArrivalTimeCam.h"
    … …  
    9287#include "MExtractedSignalPix.h"
    9388
     89#include "MPedestalSubtractedEvt.h"
     90
    9491ClassImp(MExtractTimeAndCharge);
    9592
    9693using namespace std;
    9794
    98 const Float_t MExtractTimeAndCharge::fgLoGainStartShift = -6.0;
     95const Float_t MExtractTimeAndCharge::fgLoGainStartShift = -2.5;
    9996const Byte_t  MExtractTimeAndCharge::fgLoGainSwitch     =  120;
    10097
    … …  
    105102// Sets:
    106103// - fWindowSizeHiGain and fWindowSizeLoGain to 0
    107 // - fLoGainStartShift to fgLoGainStartShift+fgOffsetLoGain
     104// - fLoGainStartShift to fgLoGainStartShift
    108105// - fLoGainSwitch     to fgLoGainSwitch
    109106//
    110107MExtractTimeAndCharge::MExtractTimeAndCharge(const char *name, const char *title)
    111     : /*fLoGainFirstSave(0),*/ fWindowSizeHiGain(0), fWindowSizeLoGain(0)
     108    : fWindowSizeHiGain(0), fWindowSizeLoGain(0)
    112109{
    113110    fName  = name  ? name  : "MExtractTimeAndCharge";
    … …  
    133130Int_t MExtractTimeAndCharge::PreProcess(MParList *pList)
    134131{
    135 
    136   if (!MExtractTime::PreProcess(pList))
    137     return kFALSE;
    138  
    139   fSignals = (MExtractedSignalCam*)pList->FindCreateObj("MExtractedSignalCam",AddSerialNumber(fNameSignalCam));
    140   if (!fSignals)
    141       return kFALSE;
    142 
    143   *fLog << flush << inf;
    144   return kTRUE;
     132    if (!MExtractTime::PreProcess(pList))
     133        return kFALSE;
     134
     135    fSignals = (MExtractedSignalCam*)pList->FindCreateObj("MExtractedSignalCam",AddSerialNumber(fNameSignalCam));
     136    if (!fSignals)
     137        return kFALSE;
     138
     139    *fLog << flush << inf;
     140    return kTRUE;
    145141}
    146142
    … …  
    157153Bool_t MExtractTimeAndCharge::ReInit(MParList *pList)
    158154{
    159 
    160   if (!MExtractTime::ReInit(pList))
    161     return kFALSE;
    162 
    163   if (!InitArrays())
    164     return kFALSE;
    165  
    166   if (fSignals)
    167     fSignals->SetUsedFADCSlices(fHiGainFirst, fHiGainLast/*+fHiLoLast*/, fNumHiGainSamples,
    168                                 fLoGainFirst, fLoGainLast, fNumLoGainSamples);
    169 
    170   return kTRUE;
     155    if (!MExtractTime::ReInit(pList))
     156        return kFALSE;
     157
     158    if (!InitArrays(fRunHeader->GetNumSamplesHiGain()+fRunHeader->GetNumSamplesLoGain()))
     159        return kFALSE;
     160
     161    if (fSignals)
     162        fSignals->SetUsedFADCSlices(fHiGainFirst, fHiGainLast, fNumHiGainSamples,
     163                                    fLoGainFirst, fLoGainLast, fNumLoGainSamples);
     164    return kTRUE;
    171165}
    172166
    … …  
    188182        const Int_t pixidx = pixel.GetPixelId();
    189183
    190         //
    191         // Preparations for the pedestal subtraction (with AB-noise correction)
    192         //
    193         Float_t *sig = fSignal->GetSamples(pixidx);
    194 
    195         // ====================================================================
    196         //    MOVE THIS TO MExtractTimeAndCharge
    197         // ====================================================================
    198         // number of hi- and lo-gains
     184        const Float_t *sig = fSignal->GetSamples(pixidx);
     185
    199186        const Int_t numh = pixel.GetNumHiGainSamples();
    200         const Int_t numl = pixel.GetNumLoGainSamples();
    201 /*
    202         // COPY HERE PRODUCING ARRAY WITH SAMPLES
    203         static MArrayB sample;
    204         sample.Set(numh+numl);
    205 
    206         if (numl>0)
    207         {
    208             memcpy(sample.GetArray(),      pixel.GetHiGainSamples(), numh);
    209             memcpy(sample.GetArray()+numh, pixel.GetLoGainSamples(), numl);
    210         }
    211 
    212         // get pedestal information
    213         const MPedestalPix &pedpix = (*fPedestals)[pixidx];
    214 
    215         // pedestal information
    216         const Int_t   ab  = pixel.HasABFlag() ? 1 : 0;
    217         const Float_t ped = pedpix.GetPedestal();
    218 
    219         // destination array
    220         static MArrayF sig;
    221         sig.Set(numh+numl);
    222 
    223         // start of destination array, end of hi-gain destination array
    224         // and start of hi-gain samples
    225         Float_t *beg = sig.GetArray();
    226         Float_t *end = beg + sig.GetSize();
    227 
    228         // determine with which pedestal (+/- AB offset) to start
    229         const Bool_t  noswap  = (ab&1)==((beg-beg)&1);
    230         const Float_t offh    = noswap ? pedpix.GetPedestalABoffset() : -pedpix.GetPedestalABoffset();
    231         const Float_t mean[2] = { ped + offh, ped - offh };
    232 
    233         //const Int_t rangeh = fHiGainLast - fHiGainFirst + 1  + fHiLoLast;
    234         //const Int_t rangel = fLoGainLast - fLoGainFirst + 1;
    235 
    236         const Byte_t *src = numl>0 ? sample.GetArray() : pixel.GetHiGainSamples();
    237         //const Byte_t *src = sample.GetArray();
    238 
    239         // Copy hi-gains into array and substract pedestal
    240         // FIXME: Shell we really subtract the pedestal from saturating slices???
    241         for (Float_t *ptr=beg; ptr<end; ptr++)
    242             *ptr = (Float_t)*src++ - mean[(ptr-beg)&1];
    243 
    244         // Determin saturation of hi-gains
    245         Byte_t *p0 = fSignal->GetSamplesRaw(pixidx); //numl>0 ? sample.GetArray() : pixel.GetHiGainSamples();
    246 //        Byte_t *p0 = sample.GetArray();
    247 
    248         Byte_t maxcont  =  0;
    249         Int_t  maxposhi = -1;
    250 
    251         Byte_t *sathi0 = 0; // first saturating hi-gain slice
    252         Byte_t *sathi1 = 0; // last saturating hi-gain slice
    253         for (Byte_t *ptr=p0+fHiGainFirst; ptr<p0+fHiGainLast+fHiLoLast+1; ptr++)
    254         {
    255             // Get hi-gain maxbincontent
    256             // Put this into its own function and loop?
    257             if (*ptr>maxcont)
    258             {
    259                 // ORGINAL:
    260                 //Int_t range = (fHiGainLast - fHiGainFirst + 1 + fHiLoLast) - fWindowSizeHiGain + 1;
    261                 // maxpos>1 && maxpos<=range
    262 
    263                 // This is a strange workaround to fit the previous
    264                 // Spline setup: TO BE CHECKED!
    265                 const Int_t pos = ptr-p0;
    266                 if (pos<fHiGainLast+1)
    267                 {
    268                     maxposhi = pos-fHiGainFirst;
    269 
    270                     if (maxposhi>1 && maxposhi<fHiGainLast-fHiGainFirst*//*+1 -fWindowSizeHiGain+1*//*+fHiLoLast*//*)
    271                         maxcont = *ptr;
    272                 }
    273             }
    274 
    275             // FIXME: Do we also have to really COUNT the number
    276             // of saturating slices, eg. for crosschecks?
    277             if (*ptr>=fSaturationLimit)
    278             {
    279                 sathi1 = ptr;
    280                 if (!sathi0)
    281                     sathi0 = ptr;
    282             }
    283         }
    284 */
    285         Int_t sathi0 = fHiGainFirst;          // First slice to extract and first saturating slice
    286         Int_t sathi1 = fHiGainLast/*+fHiLoLast*/; // Last  slice to extract and last saturating slice
     187
     188        Int_t sathi0 = fHiGainFirst;  // First slice to extract and first saturating slice
     189        Int_t sathi1 = fHiGainLast;   // Last  slice to extract and last saturating slice
    287190
    288191        Int_t maxcont;
    … …  
    290193        // Would it be better to take lastsat-firstsat?
    291194        Int_t numsathi = fSignal->GetSaturation(pixidx, fSaturationLimit, sathi0, sathi1);
    292 /*
    293         // sathi2 is the number (not the index) of first saturating slice
    294 //        const Byte_t sathi2   = sathi0<0 ? 0 : sathi0+1;
    295 
    296         // Number (not index) of first saturating slice
    297  //       const Byte_t sathi2 = sathi0==0 ? 0 : sathi0-p0+1;
    298 //        const Byte_t numsathi = sathi0==0 ? 0 : sathi1-sathi0+1;
    299 
    300 //        Int_t maxb = maxcont;
    301 
    302         // ====================================================================
    303         // FIXME: Move range out of the extractors...
    304 
    305         // Initialize maxcont for the case, it gets not set by the derived class
    306         Float_t sumhi2 =0., deltasumhi2 =-1; // Set hi-gain of MExtractedSignalPix valid
    307         Float_t timehi2=0., deltatimehi2=-1; // Set hi-gain of MArrivalTimePix valid
    308 */
     195
    309196        Float_t sumhi =0., deltasumhi =-1; // Set hi-gain of MExtractedSignalPix valid
    310197        Float_t timehi=0., deltatimehi=-1; // Set hi-gain of MArrivalTimePix valid
    311         if (numsathi<1)
     198
     199        // Do not even try to extract the hi-gain if we have
     200        // more than one saturating slice
     201        if (numsathi<2)
    312202        {
    313 
    314         const Int_t rangehi = fHiGainLast - fHiGainFirst + 1/* + fHiLoLast*/;
    315         FindTimeAndChargeHiGain2(sig/*.GetArray()*/+fHiGainFirst, rangehi,
    316                                  sumhi, deltasumhi, timehi, deltatimehi,
    317                                  numsathi, maxposhi);
    318 
    319         timehi += fHiGainFirst;
     203            const Int_t rangehi = fHiGainLast - fHiGainFirst + 1;
     204            FindTimeAndChargeHiGain2(sig+fHiGainFirst, rangehi,
     205                                     sumhi, deltasumhi, timehi, deltatimehi,
     206                                     numsathi, maxposhi);
     207
     208
     209            // Make sure that in cases the time couldn't be correctly determined
     210            // more meaningfull default values are assigned
     211            //if (timehi<fHiGainFirst || timehi>=fHiGainLast-1)
     212            if (deltatimehi>-0.5 && (timehi<0 || timehi>=fHiGainLast-fHiGainFirst+1))
     213                timehi = gRandom->Uniform(fHiGainLast-fHiGainFirst+1);
     214
     215            timehi += fHiGainFirst;
    320216        }
    321 /*
    322         Float_t sumhi = sumhi2;
    323         Float_t deltasumhi = deltasumhi2;
    324         Float_t timehi = timehi2;
    325         Float_t deltatimehi=deltatimehi2;
    326 //        Byte_t sathi = sathi2;
    327 
    328 
    329         //
    330         // Find signal in hi- and lo-gain
    331         //
    332         Float_t sumhi =0., deltasumhi =-1; // Set hi-gain of MExtractedSignalPix valid
    333         Float_t timehi=0., deltatimehi=-1; // Set hi-gain of MArrivalTimePix valid
    334         Byte_t sathi=0;
    335 
    336         // Initialize maxcont for the case, it gets not set by the derived class
    337         maxcont = fLoGainSwitch + 1;
    338 
    339         //const Int_t pixidx = pixel.GetPixelId();
    340         //const MPedestalPix  &ped = (*fPedestals)[pixidx];
    341         const Bool_t higainabflag = pixel.HasABFlag();
    342         FindTimeAndChargeHiGain(pixel.GetHiGainSamples()+fHiGainFirst, pixel.GetLoGainSamples(),
    343                                 sumhi, deltasumhi, timehi, deltatimehi,
    344                                 sathi, pedpix, higainabflag);
    345         timehi += fHiGainFirst;
    346 
    347   */
    348         // Make sure that in cases the time couldn't be correctly determined
    349         // more meaningfull default values are assigned
    350         if (timehi<=fHiGainFirst || timehi>=fHiGainLast)
    351             timehi = gRandom->Uniform(fHiGainLast-fHiGainFirst)+fHiGainFirst;
    352 
    353         Float_t sumlo =0., deltasumlo =-1.; // invalidate logain of MExtractedSignalPix
    354         Float_t timelo=0., deltatimelo=-1;  // invalidate logain of MArrivalTimePix
    355         Byte_t satlo=0;
     217
     218        // If we have saturating slices try to get a better estimate
     219        // of the arrival time than timehi or sathi0. This is
     220        // usefull to know where to start lo-gain extraction.
     221        if (numsathi>1)
     222        {
     223            const Int_t p = sathi0>1 ? sathi0-2 : sathi0-1;
     224            if (sathi0>0)
     225            {
     226                // Find the place at which the signal is maxcont/2
     227                const TVector3 vx(sig[p], sig[p+1], sig[p+2]);
     228                const TVector3 vy(p, p+1, p+2);
     229                timehi=MMath::InterpolParabLin(vx, vy, maxcont/2);
     230            }
     231            else
     232                timehi=0;
     233        }
     234
     235        Float_t sumlo =0, deltasumlo =-1;  // invalidate logain of MExtractedSignalPix
     236        Float_t timelo=0, deltatimelo=-1;  // invalidate logain of MArrivalTimePix
    356237        Int_t numsatlo=0;
    357238
    … …  
    363244        // MEANS: Hi gain has saturated, but the signal is to dim
    364245        // to extract the lo-gain properly
     246        // This could happen because the maxcont was not extracted from
     247        // all slices!
    365248        // THIS produces pulse positions ~= -1
    366249        // The signal might be handled in MCalibrateData, but hwat's about
    … …  
    376259
    377260        // FIXME: What to do with the pixel if it saturates too early???
    378         if (pixel.HasLoGain() && (maxcont > fLoGainSwitch /*|| sathi>0*/) )
     261        if (pixel.HasLoGain() && maxcont>fLoGainSwitch)
    379262        {
     263            Int_t first = numh+fLoGainFirst;
     264            Int_t last  = numh+fLoGainLast;
     265
    380266            // To determin the window in which the lo-gain is extracted
    381267            // clearly more information about the relation between the
    382268            // extraction window and the reslting time is necessary.
    383             //fLoGainStartShift = -6.0;
    384 
    385             const Byte_t fLoGainFirstSave = fLoGainFirst;
    386 
    387             // sathi is the number (not index!) of the first saturating slice
    388             // 0 indicates that no saturation was found
    389             // FIXME: Is 0.5 should be expressed by the rise time of
    390             //        the hi-gain signal!
    391             const Float_t pos = numsathi==0 ? timehi : sathi0-0.5;
    392 
    393             const Int_t lostart = TMath::FloorNint(pos+6);
    394 
    395             const Int_t newfirst = TMath::FloorNint(pos+fLoGainStartShift);
    396             fLoGainFirst = newfirst>fLoGainFirstSave ? newfirst : fLoGainFirstSave;
    397 
    398 //            if (fLoGainLast-fLoGainFirst >= fWindowSizeLoGain)
    399             {
    400 /*
    401                 Float_t sumlo2 =0., deltasumlo2 =-1.; // invalidate logain of MExtractedSignalPix
    402                 Float_t timelo2=0., deltatimelo2=-1.;  // invalidate logain of MArrivalTimePix
    403 
    404                 // Determin saturation in lo-gains
    405                 Byte_t *satlo0 = 0; // first saturating hi-gain slice
    406                 Byte_t *satlo1 = 0; // last saturating hi-gain slice
    407                 Int_t maxposlo = -1;
    408                 Byte_t maxlo = 0;
    409 
    410                 Byte_t *p0 = fSignal->GetSamplesRaw(pixidx);
    411                 for (Byte_t *ptr=p0+numh+fLoGainFirst; ptr<p0+numh+fLoGainLast+1; ptr++)
    412                 {
    413                     if (*ptr>maxlo)
    414                     {
    415                         // This is a strange workaround to fit the previous
    416                         // Spline setup: TO BE CHECKED!
    417                         const Int_t ipos = ptr-p0-numh;
    418                         if (ipos>=fLoGainFirst && ipos<=fLoGainLast)
    419                         {
    420                             maxposlo = ipos-fLoGainFirst;
    421                             maxlo = *ptr;
    422                         }
    423                     }
    424                     if (*ptr>=fSaturationLimit)
    425                     {
    426                         satlo1 = ptr;
    427                         if (!satlo0)
    428                             satlo0 = ptr;
    429                     }
    430                     // FIXME: Do we also have to really COUNT the number
    431                     // of saturating slices, eg. for crosschecks?
    432                 }
    433 
    434                 // Number of saturating lo-gain slices (this is not necessary
    435                 // identical to a counter counting the number of saturating
    436                 // slices!)
    437                 const Byte_t numsatlo = satlo0==0 ? 0 : satlo1-satlo0+1;
    438 */
    439 
    440                 Int_t satlo0 = numh+fLoGainFirst;   // First slice to extract and first saturating slice
    441                 Int_t satlo1 = numh+fLoGainLast;    // Last  slice to extract and last saturating slice
    442 
    443                 // Would it be better to take lastsat-firstsat?
    444                 Int_t maxlo;
    445                 Int_t maxposlo = fSignal->GetMax(pixidx, satlo0, satlo1, maxlo);
    446                 numsatlo = fSignal->GetSaturation(pixidx, fSaturationLimit, satlo0, satlo1);
    447 
    448                 //                const Byte_t satlo2 = numsatlo;
    449 
    450                 const Int_t rangelo = fLoGainLast - fLoGainFirst + 1;
    451                 FindTimeAndChargeLoGain2(sig/*.GetArray()*/+numh+fLoGainFirst, rangelo,
    452                                          sumlo, deltasumlo, timelo, deltatimelo,
    453                                          numsatlo, maxposlo);
    454                 timelo += fLoGainFirst;
    455 
    456 //                sumlo =sumlo2, deltasumlo=deltasumlo2; // invalidate logain of MExtractedSignalPix
    457 //                timelo=timelo2, deltatimelo=deltatimelo2;  // invalidate logain of MArrivalTimePix
    458 //                satlo = satlo2;
    459 
    460                 /*
    461                 // THIS IS NOW THE JOB OF THE EXTRACTOR!
    462                 //deltasumlo  = 0; // make logain of MExtractedSignalPix valid
    463                 //deltatimelo = 0; // make logain of MArrivalTimePix valid
    464 
    465                 const Bool_t logainabflag = (higainabflag + pixel.GetNumHiGainSamples()) & 0x1;
    466                 FindTimeAndChargeLoGain(pixel.GetLoGainSamples()+fLoGainFirst,
    467                                         sumlo, deltasumlo, timelo, deltatimelo,
    468                                         satlo, pedpix, logainabflag);
    469                 timelo += fLoGainFirst;
    470                 */
    471 
    472 
    473                 // If more than 6 lo-gain slices have saturated this is
    474                 // no physical event. We just assume that something with
    475                 // the extraction is wrong
    476                 if (numsatlo>6)
    477                     deltasumlo=deltatimelo=-1;
    478             }
    479             fLoGainFirst = fLoGainFirstSave;
     269            //
     270            // numh + fLoGainStartShift == 14 / fLoGainStartShift=-1
     271            //
     272            // The lo-gain is expected to have its raising edge
     273            // at timehi+numh+fOffsetLoGain. We start to search for the
     274            // lo-gain fLoGainStartShift slices earlier.
     275            //
     276            // Instead of fLoGainStartShift the extractor should now how many
     277            // slices before the expected raising edge the start of the
     278            // search must be placed and from there we can step 1.5 slices
     279            // back to be on the safe side.
     280            first = TMath::FloorNint(timehi+numh+fOffsetLoGain+fLoGainStartShift);
     281            if (first<0)
     282                first = 0;
     283            if (first>last)
     284                first=last;
     285
     286            /*
     287            // Currently we have to stick to this check because at least
     288            // the spline has arrays of this size...
     289            if (first>last)
     290                first = last;
     291            if (first<numh+fLoGainFirst)
     292                first = numh+fLoGainFirst;
     293             */
     294            Int_t satlo0 = first;   // First slice to extract and first saturating slice
     295            Int_t satlo1 = last;    // Last  slice to extract and last saturating slice
     296
     297            // Would it be better to take lastsat-firstsat?
     298            Int_t maxlo;
     299            Int_t maxposlo = fSignal->GetMax(pixidx, satlo0, satlo1, maxlo);
     300            numsatlo = fSignal->GetSaturation(pixidx, fSaturationLimit, satlo0, satlo1);
     301
     302            FindTimeAndChargeLoGain2(sig+first, last-first+1,
     303                                     sumlo, deltasumlo, timelo, deltatimelo,
     304                                     numsatlo, maxposlo);
    480305
    481306            // Make sure that in cases the time couldn't be correctly determined
    482307            // more meaningfull default values are assigned
    483             if (timelo<=fLoGainFirst || timelo>=fLoGainLast)
    484                 timelo = gRandom->Uniform(fLoGainLast-fLoGainFirst)+fLoGainFirst;
    485        }
    486 
     308            //if (timehi<fHiGainFirst || timehi>=fHiGainLast-1)
     309            if (deltatimelo>-0.5 && (timelo<0 || timelo>=last-first+1))
     310                timelo = gRandom->Uniform(last-first+1);
     311
     312            timelo += first-numh;
     313
     314            // If more than 6 lo-gain slices have saturated this is
     315            // no physical event. We just assume that something with
     316            // the extraction is wrong
     317            if (numsatlo>6)
     318                deltasumlo=deltatimelo=-1;
     319
     320            //if (TMath::Abs(timelo-fOffsetLoGain - timehi)>1.0)
     321            //    deltatimelo = -1;
     322        }
     323
     324        // Now store the result in the corresponding containers
    487325        MExtractedSignalPix &pix = (*fSignals)[pixidx];
    488326        MArrivalTimePix     &tix = (*fArrTime)[pixidx];
    489327        pix.SetExtractedSignal(sumhi, deltasumhi, sumlo, deltasumlo);
    490328        pix.SetGainSaturation(numsathi, numsatlo);
    491 //        pix.SetGainSaturation(sathi, satlo);
    492329
    493330        tix.SetArrivalTime(timehi, deltatimehi, timelo-fOffsetLoGain, deltatimelo);
    494331        tix.SetGainSaturation(numsathi, numsatlo);
    495 //        tix.SetGainSaturation(sathi, satlo);
    496332    } /* while (pixel.Next()) */
    497333
    … …  
    509345Int_t MExtractTimeAndCharge::ReadEnv(const TEnv &env, TString prefix, Bool_t print)
    510346{
    511 
    512347    Bool_t rc = MExtractTime::ReadEnv(env, prefix, print);
    513 
    514     if (rc)
    515       SetLoGainStartShift();
    516348
    517349    if (IsEnvDefined(env, prefix, "LoGainStartShift", print))
    … …  
    533365void MExtractTimeAndCharge::Print(Option_t *o) const
    534366{
    535   MExtractTime::Print(o);
    536 //  if (IsA()==Class())
    537 //    *fLog << GetDescriptor() << ":" << endl;
    538 
    539   *fLog << dec;
    540   //*fLog << " Taking " << fWindowSizeHiGain
    541   //      << " HiGain samples from slice " << (Int_t)fHiGainFirst
    542   //      << " to " << (Int_t)(fHiGainLast/*+fHiLoLast*/) << " incl" << endl;
    543   //*fLog << " Taking " << fWindowSizeLoGain
    544   //      << " LoGain samples from slice " << (Int_t)fLoGainFirst
    545   //      << " to " << (Int_t)fLoGainLast << " incl" << endl;
    546   //
    547   *fLog << " LoGainStartShift:   " << fLoGainStartShift << endl;
    548   *fLog << " LoGainSwitch:       " << (Int_t)fLoGainSwitch << endl;
    549 }
     367    MExtractTime::Print(o);
     368
     369    *fLog << dec;
     370    *fLog << " LoGainStartShift:   " << fLoGainStartShift << endl;
     371    *fLog << " LoGainSwitch:       " << (Int_t)fLoGainSwitch << endl;
     372}

  • trunk/MagicSoft/Mars/msignal/MExtractTimeAndCharge.h

    r8154 r8165  
    77
    88class MPedestalPix;
     9
    910class MExtractTimeAndCharge : public MExtractTime
    1011{
    1112private:
    12 
    1313  static const Float_t fgLoGainStartShift; //! Default for fLoGainStartShift (now set to: -2.8)
    1414  static const Byte_t  fgLoGainSwitch;     //! Default for fLoGainSwitch     (now set to: 100)
    1515 
    16 private:
    17 
    18 //  Byte_t  fLoGainFirstSave;       //! Temporary variable to store the original position of low-gain start slice
    1916  Float_t fLoGainStartShift;      // Shift to start searching the low-gain signal obtained from the high-gain times.
    20 
    2117  Byte_t  fLoGainSwitch;          // Limit for max. bin content before the low-gain gets extracted
    2218
    2319protected:
    24 
    2520  Int_t  fWindowSizeHiGain;       //  Window Size High-Gain
    2621  Int_t  fWindowSizeLoGain;       //  Window Size Low-Gain
    … …  
    3126
    3227public:
    33 
    3428  MExtractTimeAndCharge(const char *name=NULL, const char *title=NULL);
    3529 
    … …  
    3832  Byte_t GetLoGainSwitch      () const { return fLoGainSwitch;     }
    3933
    40   void Print(Option_t *o="") const; //*MENU*
    41  
    4234  void SetLoGainStartShift( const Float_t f=fgLoGainStartShift ) { fLoGainStartShift = f + fOffsetLoGain;  }
    4335  void SetLoGainSwitch    ( const Byte_t  i=fgLoGainSwitch     ) { fLoGainSwitch     = i; }
    4436
    45   virtual void SetWindowSize(Int_t windowh, Int_t windowl) { fWindowSizeHiGain = windowh;
    46                                                            fWindowSizeLoGain = windowl;  }
     37  virtual void SetWindowSize(Int_t windowh, Int_t windowl) { fWindowSizeHiGain = windowh; fWindowSizeLoGain = windowl;  }
    4738
    48   Bool_t ReInit(MParList *pList);
    49  
    50   virtual Bool_t InitArrays() { return kTRUE; }
    51 /*
    52   virtual void FindTimeAndChargeHiGain(Byte_t *firstused, Byte_t *logain, Float_t &sum, Float_t &dsum,
    53                                        Float_t &time, Float_t &dtime,
    54                                        Byte_t &sat, const MPedestalPix &ped, const Bool_t abflag) { }
    55  
    56   virtual void FindTimeAndChargeLoGain(Byte_t *firstused, Float_t &sum,  Float_t &dsum,
    57                                        Float_t &time, Float_t &dtime,
    58                                        Byte_t &sat, const MPedestalPix &ped, const Bool_t abflag) { }
    59   */
     39  virtual Bool_t InitArrays(Int_t n) { return kTRUE; }
     40
    6041  virtual void FindTimeAndChargeHiGain2(const Float_t *firstused, Int_t num, Float_t &sum, Float_t &dsum,
    6142                                       Float_t &time, Float_t &dtime,
    62                                        Byte_t sat, Int_t maxpos) { }
     43                                       Byte_t sat, Int_t maxpos) const { }
    6344 
    6445  virtual void FindTimeAndChargeLoGain2(const Float_t *firstused, Int_t num, Float_t &sum,  Float_t &dsum,
    6546                                       Float_t &time, Float_t &dtime,
    66                                        Byte_t sat, Int_t maxpos) { }
     47                                       Byte_t sat, Int_t maxpos) const { }
    6748
    6849  // For MExtractPedestal
     50  Bool_t ReInit(MParList *pList);
     51
     52  // TObject
     53  void Print(Option_t *o="") const; //*MENU*
    6954
    7055  ClassDef(MExtractTimeAndCharge, 2)   // Time And Charge Extractor Base Class

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

    r8154 r8165  
    11/* ======================================================================== *\
    2 ! $Name: not supported by cvs2svn $:$Id: MExtractTimeAndChargeDigitalFilter.cc,v 1.73 2006-10-24 08:24:52 tbretz Exp $
     2! $Name: not supported by cvs2svn $:$Id: MExtractTimeAndChargeDigitalFilter.cc,v 1.74 2006-10-25 18:36:26 tbretz Exp $
    33! --------------------------------------------------------------------------
    44!
    … …  
    6262#include <fstream>
    6363
    64 #include <TH1.h>
    65 #include <TH2.h>
    66 #include <TMatrix.h>
     64#include <TRandom.h>
    6765
    6866#include "MLog.h"
    … …  
    7674#include "MExtralgoDigitalFilter.h"
    7775
    78 #include "MPedestalPix.h"
    79 
    8076ClassImp(MExtractTimeAndChargeDigitalFilter);
    8177
    8278using namespace std;
    8379
    84 const Byte_t MExtractTimeAndChargeDigitalFilter::fgHiGainFirst             =  0;
    85 const Byte_t MExtractTimeAndChargeDigitalFilter::fgHiGainLast              = 15;
    86 const Byte_t MExtractTimeAndChargeDigitalFilter::fgLoGainFirst             =  1;
    87 const Byte_t MExtractTimeAndChargeDigitalFilter::fgLoGainLast              = 14;
    88 //const Int_t  MExtractTimeAndChargeDigitalFilter::fgWindowSizeHiGain        =  4;
    89 //const Int_t  MExtractTimeAndChargeDigitalFilter::fgWindowSizeLoGain        =  6;
    90 const Int_t  MExtractTimeAndChargeDigitalFilter::fgBinningResolutionHiGain = 10;
    91 const Int_t  MExtractTimeAndChargeDigitalFilter::fgBinningResolutionLoGain = 10;
    92 const Int_t  MExtractTimeAndChargeDigitalFilter::fgSignalStartBinHiGain    =  4;
    93 const Int_t  MExtractTimeAndChargeDigitalFilter::fgSignalStartBinLoGain    =  4;
    94 const Float_t MExtractTimeAndChargeDigitalFilter::fgOffsetLoGain           =  0.95;
    95 //const Float_t MExtractTimeAndChargeDigitalFilter::fgLoGainStartShift       = -1.8;
     80const Byte_t  MExtractTimeAndChargeDigitalFilter::fgHiGainFirst             =  0;
     81const Byte_t  MExtractTimeAndChargeDigitalFilter::fgHiGainLast              = 15;
     82const Byte_t  MExtractTimeAndChargeDigitalFilter::fgLoGainFirst             =  1;
     83const Byte_t  MExtractTimeAndChargeDigitalFilter::fgLoGainLast              = 14;
     84const Int_t   MExtractTimeAndChargeDigitalFilter::fgBinningResolutionHiGain = 10;
     85const Int_t   MExtractTimeAndChargeDigitalFilter::fgBinningResolutionLoGain = 10;
     86const Int_t   MExtractTimeAndChargeDigitalFilter::fgSignalStartBinHiGain    =  4;
     87const Int_t   MExtractTimeAndChargeDigitalFilter::fgSignalStartBinLoGain    =  4;
     88const Float_t MExtractTimeAndChargeDigitalFilter::fgOffsetLoGain            =  0.95;
    9689
    9790// --------------------------------------------------------------------------
    … …  
    109102    : fBinningResolutionHiGain(fgBinningResolutionHiGain),
    110103    fBinningResolutionLoGain(fgBinningResolutionLoGain),
    111     fAutomaticWeights(kTRUE), fRandomIter(0)
     104    fAutomaticWeights(kTRUE)
    112105{
    113106    fName  = name  ? name  : "MExtractTimeAndChargeDigitalFilter";
    … …  
    116109    SetRange(fgHiGainFirst, fgHiGainLast, fgLoGainFirst, fgLoGainLast);
    117110    SetWindowSize(3, 5);
    118 //    SetBinningResolution();
    119 //    SetSignalStartBin();
    120 
    121     SetOffsetLoGain(fgOffsetLoGain);          // Order between both
    122 //    SetLoGainStartShift(fgLoGainStartShift);  // is important
     111    SetOffsetLoGain(fgOffsetLoGain);
    123112}
    124113
    … …  
    189178  fSqrtHiGainSamples = TMath::Sqrt(fNumHiGainSamples);
    190179  fSqrtLoGainSamples = TMath::Sqrt(fNumLoGainSamples);
    191 
    192180}
    193181
    … …  
    256244// Gets called in the ReInit() and initialized the arrays
    257245//
    258 Bool_t MExtractTimeAndChargeDigitalFilter::InitArrays()
     246Bool_t MExtractTimeAndChargeDigitalFilter::InitArrays(Int_t n)
    259247{
    260248    if (!fRunHeader)
    261249        return kFALSE;
    262 
    263     //const Int_t rangehi = (Int_t)(fHiGainLast - fHiGainFirst + 1 + fHiLoLast);
    264     //const Int_t rangelo = (Int_t)(fLoGainLast - fLoGainFirst + 1);
    265 
    266     //cout << "ARRAYS INITIALIZED TO : " << rangehi << " " << rangelo << endl;
    267 
    268     //fHiGainSignal.Set(rangehi);
    269     //fLoGainSignal.Set(rangelo);
    270250
    271251    return GetWeights();
    … …  
    296276                                                                  Float_t &sum, Float_t &dsum,
    297277                                                                  Float_t &time, Float_t &dtime,
    298                                                                   Byte_t sat, Int_t maxpos)
     278                                                                  Byte_t sat, Int_t maxpos) const
    299279{
    300280    // Do some handling if maxpos is last slice!
    … …  
    308288    if (IsNoiseCalculation())
    309289    {
    310         if (fRandomIter == fBinningResolutionHiGain)
    311             fRandomIter = 0;
    312 
    313         sum = df.ExtractNoise(fRandomIter);
    314         fRandomIter++;
     290        sum = df.ExtractNoise(gRandom->Integer(fBinningResolutionHiGain));
    315291        return;
    316292    }
    … …  
    324300                                                                  Float_t &sum, Float_t &dsum,
    325301                                                                  Float_t &time, Float_t &dtime,
    326                                                                   Byte_t sat, Int_t maxpos)
     302                                                                  Byte_t sat, Int_t maxpos) const
    327303{
    328304    MExtralgoDigitalFilter df(fBinningResolutionLoGain, fWindowSizeLoGain,
    … …  
    335311    if (IsNoiseCalculation())
    336312    {
    337         if (fRandomIter == fBinningResolutionLoGain)
    338             fRandomIter = 0;
    339 
    340         sum = df.ExtractNoise(fRandomIter);
     313        sum = df.ExtractNoise(gRandom->Integer(fBinningResolutionHiGain));
    341314        return;
    342315    }
    … …  
    347320}
    348321
    349 /*
     322
    350323// --------------------------------------------------------------------------
    351324//
    352 // Apply the digital filter algorithm to the high-gain slices.
    353 //
    354 void MExtractTimeAndChargeDigitalFilter::FindTimeAndChargeHiGain(Byte_t *first, Byte_t *logain, Float_t &sum, Float_t &dsum,
    355                                                                  Float_t &time, Float_t &dtime,
    356                                                                  Byte_t &sat, const MPedestalPix &ped, const Bool_t abflag)
    357 {
    358     Int_t range = fHiGainLast - fHiGainFirst + 1;
    359     const Byte_t *end = first + range;
    360     Byte_t       *p  = first;
    361 
    362     const Float_t pedes  = ped.GetPedestal();
    363     const Float_t ABoffs = ped.GetPedestalABoffset();
    364 
    365     const Float_t pedmean[2] = { pedes + ABoffs, pedes - ABoffs };
    366 
    367     Byte_t maxcont = 0;
    368     Int_t  maxpos  = 0;  // obsolete for Digital Filter (used in spline!)
    369     sat = 0;
    370 
    371     //
    372     // Check for saturation in all other slices
    373     //
    374     Int_t ids = fHiGainFirst;
    375     Float_t *sample = fHiGainSignal.GetArray();
    376 
    377     while (p<end)
    378     {
    379         *sample++ = (Float_t)*p - pedmean[(ids++ + abflag) & 0x1];
    380 
    381         if (*p > maxcont)
    382         {
    383             maxpos = ids-fHiGainFirst-1;
    384             // range-fWindowSizeHiGain+1 == fHiLoLast isn't it?
    385             //if (maxpos > 1 && maxpos < (range - fWindowSizeHiGain + 1))
    386             if (maxpos > 1 && maxpos < (range-1*//* - fWindowSizeHiGain + 1*//*))
    387                 maxcont = *p;
    388         }
    389 
    390         if (*p++ >= fSaturationLimit)
    391             if (!sat)
    392                 sat = ids-fHiGainFirst;
    393     }
    394 
    395     if (fHiLoLast != 0)
    396     {
    397         end = logain + fHiLoLast;
    398 
    399         while (logain<end)
    400         {
    401             *sample++ = (Float_t)*logain - pedmean[(ids++ + abflag) & 0x1];
    402 
    403             if (*logain > maxcont)
    404             {
    405                 maxpos = ids-fHiGainFirst-1;
    406 
    407                 // FIXME!!! MUST BE THERE!
    408                 // range-fWindowSizeHiGain+1 == fHiLoLast isn't it?
    409                 //if (maxpos > 1 && maxpos < (range - fWindowSizeHiGain + 1))
    410                 //    maxcont = *logain;
    411             }
    412 
    413             if (*logain++ >= fSaturationLimit)
    414                 if (!sat)
    415                     sat = ids-fHiGainFirst;
    416 
    417             range++;
    418         }
    419     }
    420 
    421     FindTimeAndChargeHiGain2(fHiGainSignal.GetArray(), range,
    422                              sum, dsum, time, dtime,
    423                              sat, 0);
    424 }
    425 
    426 // --------------------------------------------------------------------------
    427 //
    428 // Apply the digital filter algorithm to the low-gain slices.
    429 //
    430 void MExtractTimeAndChargeDigitalFilter::FindTimeAndChargeLoGain(Byte_t *ptr, Float_t &sum, Float_t &dsum,
    431                                                                  Float_t &time, Float_t &dtime,
    432                                                                  Byte_t &sat, const MPedestalPix &ped, const Bool_t abflag)
    433 {
    434     const Int_t range = fLoGainLast - fLoGainFirst + 1;
    435 
    436     const Byte_t *end = ptr + range;
    437     Byte_t *p     = ptr;
    438 
    439     //
    440     // Prepare the low-gain pedestal
    441     //
    442     const Float_t pedes  = ped.GetPedestal();
    443     const Float_t ABoffs = ped.GetPedestalABoffset();
    444 
    445     const Float_t pedmean[2] = { pedes + ABoffs, pedes - ABoffs };
    446 
    447     //
    448     // Check for saturation in all other slices
    449     //
    450     Float_t *sample = fLoGainSignal.GetArray();
    451     Int_t    ids    = fLoGainFirst;
    452 
    453     while (p<end)
    454     {
    455         *sample++ = (Float_t)*p - pedmean[(ids++ + abflag) & 0x1];
    456 
    457         if (*p++ >= fSaturationLimit)
    458             sat++;
    459     }
    460 
    461     FindTimeAndChargeLoGain2(fLoGainSignal.GetArray(), range,
    462                              sum, dsum, time, dtime,
    463                              sat, 0);
    464 
    465 }
    466 */
    467 // --------------------------------------------------------------------------
    468 //
    469325// Read the setup from a TEnv, eg:
    470 //   MJPedestal.MExtractor.WindowSizeHiGain: 6
    471 //   MJPedestal.MExtractor.WindowSizeLoGain: 6
    472 //   MJPedestal.MExtractor.BinningResolutionHiGain: 10
    473 //   MJPedestal.MExtractor.BinningResolutionLoGain: 10
    474326//   MJPedestal.MExtractor.WeightsFile: filename
    475327//   MJPedestal.MExtractor.AutomaticWeights: off
    … …  
    485337      rc = kTRUE;
    486338  }
    487    /*
    488   Byte_t hw = fWindowSizeHiGain;
    489   Byte_t lw = fWindowSizeLoGain;
    490 
    491   if (IsEnvDefined(env, prefix, "WindowSizeHiGain", print))
    492     {
    493       hw = GetEnvValue(env, prefix, "WindowSizeHiGain", hw);
    494       rc = kTRUE;
    495     }
    496   if (IsEnvDefined(env, prefix, "WindowSizeLoGain", print))
    497     {
    498       lw = GetEnvValue(env, prefix, "WindowSizeLoGain", lw);
    499       rc = kTRUE;
    500     }
    501  
    502   if (rc)
    503     SetWindowSize(hw, lw);
    504 
    505   Bool_t rc2 = kFALSE;
    506   Int_t brh = fBinningResolutionHiGain;
    507   Int_t brl = fBinningResolutionLoGain;
    508  
    509   if (IsEnvDefined(env, prefix, "BinningResolutionHiGain", print))
    510     {
    511       brh = GetEnvValue(env, prefix, brh);
    512       rc2 = kTRUE;
    513     }
    514   if (IsEnvDefined(env, prefix, "BinningResolutionLoGain", print))
    515     {
    516       brl = GetEnvValue(env, prefix, brl);
    517       rc2 = kTRUE;
    518     }
    519  
    520   if (rc2)
    521     {
    522       SetBinningResolution(brh, brl);
    523       rc = kTRUE;
    524     }
    525    */
     339
    526340  if (IsEnvDefined(env, prefix, "WeightsFile", print))
    527341    {
    … …  
    859673    return ReadWeightsFile(name, path);
    860674}
    861 /*
     675
    862676//----------------------------------------------------------------------------
    863677//
    864 // Create the weights file
    865 // Beware that the shape-histogram has to contain the pulse starting at bin 1
    866 //
    867 Bool_t MExtractTimeAndChargeDigitalFilter::WriteWeightsFile(TString filename, TH1F *shapehi, TH2F *autocorrhi,
    868                                                             TH1F *shapelo, TH2F *autocorrlo )
    869 {
    870 
    871   const Int_t nbinshi = shapehi->GetNbinsX();
    872   Float_t binwidth    = shapehi->GetBinWidth(1);
    873 
    874   TH1F *derivativehi = new TH1F(Form("%s%s",shapehi->GetName(),"_der"),
    875                                 Form("%s%s",shapehi->GetTitle()," derivative"),
    876                                 nbinshi,
    877                                 shapehi->GetBinLowEdge(1),
    878                                 shapehi->GetBinLowEdge(nbinshi)+binwidth);
    879 
    880   //
    881   // Calculate the derivative of shapehi
    882   //
    883   for (Int_t i = 1; i<nbinshi+1;i++)
    884     {
    885       derivativehi->SetBinContent(i,
    886                                 ((shapehi->GetBinContent(i+1)-shapehi->GetBinContent(i-1))/2./binwidth));
    887       derivativehi->SetBinError(i,
    888                               (sqrt(shapehi->GetBinError(i+1)*shapehi->GetBinError(i+1)
    889                                     +shapehi->GetBinError(i-1)*shapehi->GetBinError(i-1))/2./binwidth));
    890     }
    891  
    892   //
    893   // normalize the shapehi, such that the integral for fWindowSize slices is one!
    894   //
    895   Float_t sum    = 0;
    896   Int_t lasttemp = fBinningResolutionHiGain * (fSignalStartBinHiGain + fWindowSizeHiGain);
    897   lasttemp       = lasttemp > nbinshi ? nbinshi : lasttemp;
    898  
    899   for (Int_t i=fBinningResolutionHiGain*fSignalStartBinHiGain; i<lasttemp; i++) {
    900     sum += shapehi->GetBinContent(i);
    901   }
    902   sum /= fBinningResolutionHiGain;
    903 
    904   shapehi->Scale(1./sum);
    905   derivativehi->Scale(1./sum);
    906 
    907   //
    908   // read in the noise auto-correlation function:
    909   //
    910   TMatrix Bhi(fWindowSizeHiGain,fWindowSizeHiGain);
    911 
    912   for (Int_t i=0; i<fWindowSizeHiGain; i++){
    913     for (Int_t j=0; j<fWindowSizeHiGain; j++){
    914       Bhi[i][j]=autocorrhi->GetBinContent(i+1,j+1); //+fSignalStartBinHiGain +fSignalStartBinHiGain
    915     }
    916   } 
    917   Bhi.Invert();
    918 
    919   const Int_t nsizehi = fWindowSizeHiGain*fBinningResolutionHiGain;
    920   fAmpWeightsHiGain.Set(nsizehi);
    921   fTimeWeightsHiGain.Set(nsizehi); 
    922  
    923   //
    924   // Loop over relative time in one BinningResolution interval
    925   //
    926   Int_t start = fBinningResolutionHiGain*(fSignalStartBinHiGain + 1);
    927 
    928   for (Int_t i = -fBinningResolutionHiGain/2+1; i<=fBinningResolutionHiGain/2; i++)
    929     {
    930  
    931       TMatrix g(fWindowSizeHiGain,1);
    932       TMatrix gT(1,fWindowSizeHiGain);
    933       TMatrix d(fWindowSizeHiGain,1);
    934       TMatrix dT(1,fWindowSizeHiGain);
    935    
    936       for (Int_t count=0; count < fWindowSizeHiGain; count++){
    937      
    938         g[count][0]=shapehi->GetBinContent(start
    939                                          +fBinningResolutionHiGain*count+i);
    940         gT[0][count]=shapehi->GetBinContent(start
    941                                           +fBinningResolutionHiGain*count+i);
    942         d[count][0]=derivativehi->GetBinContent(start
    943                                               +fBinningResolutionHiGain*count+i);
    944         dT[0][count]=derivativehi->GetBinContent(start
    945                                                +fBinningResolutionHiGain*count+i);
    946       }
    947    
    948       TMatrix m_denom = (gT*(Bhi*g))*(dT*(Bhi*d)) - (dT*(Bhi*g))*(dT*(Bhi*g));
    949       Float_t   denom = m_denom[0][0];  // ROOT thinks, m_denom is still a matrix
    950      
    951       TMatrix m_first = dT*(Bhi*d);       // ROOT thinks, m_first is still a matrix
    952       Float_t   first = m_first[0][0]/denom;
    953      
    954       TMatrix m_last  = gT*(Bhi*d);       // ROOT thinks, m_last  is still a matrix
    955       Float_t   last  = m_last[0][0]/denom;
    956      
    957       TMatrix m1 = gT*Bhi;
    958       m1 *= first;
    959      
    960       TMatrix m2 = dT*Bhi;
    961       m2 *=last;
    962      
    963       TMatrix w_amp = m1 - m2;
    964      
    965       TMatrix m_first1 = gT*(Bhi*g);
    966       Float_t   first1 = m_first1[0][0]/denom;
    967      
    968       TMatrix m_last1  = gT*(Bhi*d);
    969       Float_t   last1  = m_last1 [0][0]/denom;
    970      
    971       TMatrix m11 = dT*Bhi;
    972       m11 *=first1;
    973      
    974       TMatrix m21 = gT*Bhi;
    975       m21 *=last1;
    976      
    977       TMatrix w_time= m11 - m21;
    978      
    979       for (Int_t count=0; count < fWindowSizeHiGain; count++)
    980         {
    981           const Int_t idx = i+fBinningResolutionHiGain/2+fBinningResolutionHiGain*count-1;
    982           fAmpWeightsHiGain [idx] = w_amp [0][count];
    983           fTimeWeightsHiGain[idx] = w_time[0][count];
    984         }
    985      
    986     } // end loop over i
    987 
    988   //
    989   // Low Gain histograms
    990   //
    991   TH1F *derivativelo = NULL;
    992   if (shapelo)
    993     {
    994       const Int_t nbinslo  = shapelo->GetNbinsX();
    995       binwidth = shapelo->GetBinWidth(1);
    996      
    997       derivativelo = new TH1F(Form("%s%s",shapelo->GetName(),"_der"),
    998                               Form("%s%s",shapelo->GetTitle()," derivative"),
    999                               nbinslo,
    1000                               shapelo->GetBinLowEdge(1),
    1001                               shapelo->GetBinLowEdge(nbinslo)+binwidth);
    1002      
    1003       //
    1004       // Calculate the derivative of shapelo
    1005       //
    1006       for (Int_t i = 1; i<nbinslo+1;i++)
    1007         {
    1008           derivativelo->SetBinContent(i,
    1009                                       ((shapelo->GetBinContent(i+1)-shapelo->GetBinContent(i-1))/2./binwidth));
    1010           derivativelo->SetBinError(i,
    1011                                     (sqrt(shapelo->GetBinError(i+1)*shapelo->GetBinError(i+1)
    1012                                           +shapelo->GetBinError(i-1)*shapelo->GetBinError(i-1))/2./binwidth));
    1013         }
    1014      
    1015       //
    1016       // normalize the shapelo, such that the integral for fWindowSize slices is one!
    1017       //
    1018       sum      = 0;
    1019       lasttemp = fBinningResolutionLoGain * (fSignalStartBinLoGain + fWindowSizeLoGain);
    1020       lasttemp = lasttemp > nbinslo ? nbinslo : lasttemp;
    1021      
    1022       for (Int_t i=fBinningResolutionLoGain*fSignalStartBinLoGain; i<lasttemp; i++)
    1023         sum += shapelo->GetBinContent(i);
    1024 
    1025       sum /= fBinningResolutionLoGain;
    1026      
    1027       shapelo->Scale(1./sum);
    1028       derivativelo->Scale(1./sum);
    1029      
    1030       //
    1031       // read in the noise auto-correlation function:
    1032       //
    1033       TMatrix Blo(fWindowSizeLoGain,fWindowSizeLoGain);
    1034      
    1035       for (Int_t i=0; i<fWindowSizeLoGain; i++){
    1036         for (Int_t j=0; j<fWindowSizeLoGain; j++){
    1037           Blo[i][j]=autocorrlo->GetBinContent(i+1+fSignalStartBinLoGain,j+1+fSignalStartBinLoGain);
    1038         }
    1039       } 
    1040       Blo.Invert();
    1041 
    1042       const Int_t nsizelo = fWindowSizeLoGain*fBinningResolutionLoGain;
    1043       fAmpWeightsLoGain.Set(nsizelo);
    1044       fTimeWeightsLoGain.Set(nsizelo); 
    1045  
    1046       //
    1047       // Loop over relative time in one BinningResolution interval
    1048       //
    1049       Int_t start = fBinningResolutionLoGain*fSignalStartBinLoGain + fBinningResolutionLoGain/2;
    1050      
    1051       for (Int_t i = -fBinningResolutionLoGain/2+1; i<=fBinningResolutionLoGain/2; i++)
    1052         {
    1053          
    1054           TMatrix g(fWindowSizeLoGain,1);
    1055           TMatrix gT(1,fWindowSizeLoGain);
    1056           TMatrix d(fWindowSizeLoGain,1);
    1057           TMatrix dT(1,fWindowSizeLoGain);
    1058          
    1059           for (Int_t count=0; count < fWindowSizeLoGain; count++){
    1060            
    1061             g[count][0] = shapelo->GetBinContent(start
    1062                                              +fBinningResolutionLoGain*count+i);
    1063             gT[0][count]= shapelo->GetBinContent(start
    1064                                               +fBinningResolutionLoGain*count+i);
    1065             d[count][0] = derivativelo->GetBinContent(start
    1066                                                   +fBinningResolutionLoGain*count+i);
    1067             dT[0][count]= derivativelo->GetBinContent(start
    1068                                                    +fBinningResolutionLoGain*count+i);
    1069           }
    1070          
    1071           TMatrix m_denom = (gT*(Blo*g))*(dT*(Blo*d)) - (dT*(Blo*g))*(dT*(Blo*g));
    1072           Float_t   denom = m_denom[0][0];  // ROOT thinks, m_denom is still a matrix
    1073          
    1074           TMatrix m_first = dT*(Blo*d);       // ROOT thinks, m_first is still a matrix
    1075           Float_t   first = m_first[0][0]/denom;
    1076          
    1077           TMatrix m_last  = gT*(Blo*d);       // ROOT thinks, m_last  is still a matrix
    1078           Float_t   last  = m_last[0][0]/denom;
    1079          
    1080           TMatrix m1 = gT*Blo;
    1081           m1 *= first;
    1082          
    1083           TMatrix m2 = dT*Blo;
    1084           m2 *=last;
    1085          
    1086           TMatrix w_amp = m1 - m2;
    1087          
    1088           TMatrix m_first1 = gT*(Blo*g);
    1089           Float_t   first1 = m_first1[0][0]/denom;
    1090          
    1091           TMatrix m_last1  = gT*(Blo*d);
    1092           Float_t   last1  = m_last1 [0][0]/denom;
    1093          
    1094           TMatrix m11 = dT*Blo;
    1095           m11 *=first1;
    1096          
    1097           TMatrix m21 = gT*Blo;
    1098           m21 *=last1;
    1099          
    1100           TMatrix w_time= m11 - m21;
    1101          
    1102           for (Int_t count=0; count < fWindowSizeLoGain; count++)
    1103             {
    1104               const Int_t idx = i+fBinningResolutionLoGain/2+fBinningResolutionLoGain*count-1;
    1105               fAmpWeightsLoGain [idx] = w_amp [0][count];
    1106               fTimeWeightsLoGain[idx] = w_time[0][count];
    1107             }
    1108      
    1109         } // end loop over i
    1110     }
    1111  
    1112   ofstream fn(filename.Data());
    1113 
    1114   fn << "# High Gain Weights: " << fWindowSizeHiGain << " " << fBinningResolutionHiGain << endl;
    1115   fn << "# (Amplitude)  (Time) " << endl;
    1116 
    1117   for (Int_t i=0; i<nsizehi; i++)
    1118     fn << "\t" << fAmpWeightsHiGain[i] << "\t" << fTimeWeightsHiGain[i] << endl;
    1119 
    1120   fn << "# Low Gain Weights: " << fWindowSizeLoGain << " " << fBinningResolutionLoGain << endl;
    1121   fn << "# (Amplitude)  (Time) " << endl;
    1122 
    1123   for (Int_t i=0; i<nsizehi; i++)
    1124     fn << "\t" << fAmpWeightsLoGain[i] << "\t" << fTimeWeightsLoGain[i] << endl;
    1125 
    1126   delete derivativehi;
    1127   if (derivativelo)
    1128     delete derivativelo;
    1129 
    1130   return kTRUE;
    1131 }
    1132 */
     678// Print the setup of the digital filter extraction used. Use
     679//  the option "weights" if you want to print also all weights.
     680//
    1133681void MExtractTimeAndChargeDigitalFilter::Print(Option_t *o) const
    1134682{
    … …  
    1137685
    1138686    MExtractTimeAndCharge::Print(o);
    1139     *fLog << " Window Size HiGain: " << fWindowSizeHiGain        << "  LoGain: " << fWindowSizeLoGain << endl;
    1140     *fLog << " Binning Res HiGain: " << fBinningResolutionHiGain << "  LoGain: " << fBinningResolutionHiGain << endl;
     687    *fLog << " Window Size HiGain: " << setw(2) << fWindowSizeHiGain        << "  LoGain: " << setw(2) << fWindowSizeLoGain << endl;
     688    *fLog << " Binning Res HiGain: " << setw(2) << fBinningResolutionHiGain << "  LoGain: " << setw(2) << fBinningResolutionHiGain << endl;
    1141689    *fLog << " Weights File desired: " << (fNameWeightsFile.IsNull()?"-":fNameWeightsFile.Data()) << endl;
    1142690    if (!fNameWeightsFileSet.IsNull())

  • trunk/MagicSoft/Mars/msignal/MExtractTimeAndChargeDigitalFilter.h

    r8154 r8165  
    1010#endif
    1111
    12 #ifndef MARS_MArrayI
    13 #include "MArrayI.h"
    14 #endif
    15 
    16 class TH1F;
    17 class TH2F;
    18 class MPedestalPix;
    1912class MCalibrationPattern;
    2013
    … …  
    2215{
    2316private:
    24     static const Byte_t fgHiGainFirst;             //! Default for fHiGainFirst       (now set to: 0)
    25     static const Byte_t fgHiGainLast;              //! Default for fHiGainLast        (now set to:14)
    26     static const Byte_t fgLoGainFirst;             //! Default for fLoGainFirst       (now set to: 3)
    27     static const Byte_t fgLoGainLast;              //! Default for fLoGainLast        (now set to:14)
    28 //    static const Int_t  fgWindowSizeHiGain;        //! Default for fWindowSizeHiGain  (now set to: 6)
    29 //    static const Int_t  fgWindowSizeLoGain;        //! Default for fWindowSizeLoGain  (now set to: 6)
    30     static const Int_t  fgBinningResolutionHiGain; //! Default for fBinningResolutionHiGain (now set to: 10)
    31     static const Int_t  fgBinningResolutionLoGain; //! Default for fBinningResolutionLoGain (now set to: 10)
    32     static const Int_t  fgSignalStartBinHiGain;    //! Default for fSignalStartBinHiGain (now set to: 4)
    33     static const Int_t  fgSignalStartBinLoGain;    //! Default for fSignalStartBinLoGain (now set to: 4)
    34     static const TString fgNameWeightsFile;        //! "cosmics_weights.dat"
    35     static const Float_t fgOffsetLoGain;           //! Default for fOffsetLoGain (now set to 1.7)
    36 //    static const Float_t fgLoGainStartShift;       //! Default for fLoGainStartShift (now set to -1.8)
     17    static const Byte_t  fgHiGainFirst;             //! Default for fHiGainFirst       (now set to: 0)
     18    static const Byte_t  fgHiGainLast;              //! Default for fHiGainLast        (now set to:14)
     19    static const Byte_t  fgLoGainFirst;             //! Default for fLoGainFirst       (now set to: 3)
     20    static const Byte_t  fgLoGainLast;              //! Default for fLoGainLast        (now set to:14)
     21    static const Int_t   fgBinningResolutionHiGain; //! Default for fBinningResolutionHiGain (now set to: 10)
     22    static const Int_t   fgBinningResolutionLoGain; //! Default for fBinningResolutionLoGain (now set to: 10)
     23    static const Int_t   fgSignalStartBinHiGain;    //! Default for fSignalStartBinHiGain (now set to: 4)
     24    static const Int_t   fgSignalStartBinLoGain;    //! Default for fSignalStartBinLoGain (now set to: 4)
     25    static const TString fgNameWeightsFile;         //! "cosmics_weights.dat"
     26    static const Float_t fgOffsetLoGain;            //! Default for fOffsetLoGain (now set to 1.7)
    3727
    38     MCalibrationPattern  *fCalibPattern;          //! Calibration DM pattern
     28    MCalibrationPattern  *fCalibPattern;            //! Calibration DM pattern
    3929
    40 //    MArrayF fHiGainSignal;                        //! Need fast access to the signals in a float way
    41 //    MArrayF fLoGainSignal;                        //! Store them in separate arrays
     30    Int_t   fBinningResolutionHiGain;               //  Number of weights per bin High-Gain
     31    Int_t   fBinningResolutionLoGain;               //  Number of weights per bin Low-Gain
    4232
    43 //    Int_t   fSignalStartBinHiGain;                //! Start bin from when on to apply weights
    44 //    Int_t   fSignalStartBinLoGain;                //! Start bin from when on to apply weights
     33    MArrayF fAmpWeightsHiGain;                      //! Amplitude weights High-Gain (from weights file)
     34    MArrayF fTimeWeightsHiGain;                     //! Time weights High-Gain (from weights file)
     35    MArrayF fAmpWeightsLoGain;                      //! Amplitude weights Low-Gain (from weights file)
     36    MArrayF fTimeWeightsLoGain;                     //! Time weights Low-Gain (from weights file)
    4537
    46     Int_t   fBinningResolutionHiGain;             //  Number of weights per bin High-Gain
    47     Int_t   fBinningResolutionLoGain;             //  Number of weights per bin Low-Gain
     38    MArrayF fPulseHiGain;                           //! Pulse Shape Hi-Gain (for chisq)
     39    MArrayF fPulseLoGain;                           //! Pulse Shape Lo-Gain (for chisq)
    4840
    49     MArrayF fAmpWeightsHiGain;                    //! Amplitude weights High-Gain (from weights file)
    50     MArrayF fTimeWeightsHiGain;                   //! Time weights High-Gain (from weights file)
    51     MArrayF fAmpWeightsLoGain;                    //! Amplitude weights Low-Gain (from weights file)
    52     MArrayF fTimeWeightsLoGain;                   //! Time weights Low-Gain (from weights file)
    53 
    54     MArrayF fPulseHiGain;                    //!
    55     MArrayF fPulseLoGain;                    //!
    56 
    57     TString fNameWeightsFile;                     // Name of the weights file
    58     Bool_t  fAutomaticWeights;                    // Flag whether weight should be determined automatically
    59     TString fNameWeightsFileSet;                  //! Flag if weights have alreayd been set
    60     Int_t   fRandomIter;                          //! Counter used to randomize weights for noise calculation
     41    TString fNameWeightsFile;                       // Name of the weights file
     42    Bool_t  fAutomaticWeights;                      // Flag whether weight should be determined automatically
     43    TString fNameWeightsFileSet;                    //! Flag if weights have alreayd been set
    6144
    6245    // MExtractTimeAndChargeDigitalFilter
    … …  
    6952
    7053    // MExtractTimeAndCharge
    71     Bool_t InitArrays();
     54    Bool_t InitArrays(Int_t n);
    7255
    7356    // MTask
    … …  
    8265    MExtractTimeAndChargeDigitalFilter(const char *name=NULL, const char *title=NULL);
    8366    ~MExtractTimeAndChargeDigitalFilter() { }
    84 /*
    85     Bool_t WriteWeightsFile(TString filename,
    86                             TH1F *shapehi, TH2F *autocorrhi,
    87                             TH1F *shapelo=NULL, TH2F *autocorrlo=NULL );
    8867
    89   */
    9068    void SetNameWeightsFile(TString s="")
    9169    {
    … …  
    9977    void SetBinningResolution(const Int_t rh=fgBinningResolutionHiGain, const Int_t rl=fgBinningResolutionLoGain)
    10078    {
    101         fBinningResolutionHiGain = rh & ~1;
    102         fBinningResolutionLoGain = rl & ~1;
    103     }
    104 /*
    105     void SetSignalStartBin( const Int_t sh=fgSignalStartBinHiGain, const Int_t sl=fgSignalStartBinLoGain)
    106     {
    107         fSignalStartBinHiGain = sh;
    108         fSignalStartBinLoGain = sl;
     79        fBinningResolutionHiGain = rh;
     80        fBinningResolutionLoGain = rl;
    10981    }
    11082
    111  */
    11283    void SetWindowSize( Int_t windowh, Int_t windowl);
    11384    const char* GetNameWeightsFile() const  { return fNameWeightsFile.Data(); }
    11485
    11586    void Print(Option_t *o="") const; //*MENU*
    116 /*
    117     void FindTimeAndChargeHiGain(Byte_t *firstused, Byte_t *logain, Float_t &sum, Float_t &dsum,
    118                                  Float_t &time, Float_t &dtime,
    119                                  Byte_t &sat, const MPedestalPix &ped, const Bool_t abflag);
    120     void FindTimeAndChargeLoGain(Byte_t *firstused, Float_t &sum,  Float_t &dsum,
    121                                  Float_t &time, Float_t &dtime,
    122                                  Byte_t &sat, const MPedestalPix &ped, const Bool_t abflag);
    123  */
     87
    12488    void FindTimeAndChargeHiGain2(const Float_t *firstused, Int_t num, Float_t &sum, Float_t &dsum,
    12589                                  Float_t &time, Float_t &dtime,
    126                                   Byte_t sat, Int_t maxpos);
     90                                  Byte_t sat, Int_t maxpos) const;
    12791
    12892    void FindTimeAndChargeLoGain2(const Float_t *firstused, Int_t num, Float_t &sum,  Float_t &dsum,
    12993                                  Float_t &time, Float_t &dtime,
    130                                   Byte_t sat, Int_t maxpos);
     94                                  Byte_t sat, Int_t maxpos) const;
    13195
    13296    ClassDef(MExtractTimeAndChargeDigitalFilter, 2)   // Hendrik's digital filter

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

    r8158 r8165  
    11/* ======================================================================== *\
    2 ! $Name: not supported by cvs2svn $:$Id: MExtractTimeAndChargeSpline.cc,v 1.62 2006-10-24 12:39:00 tbretz Exp $
     2! $Name: not supported by cvs2svn $:$Id: MExtractTimeAndChargeSpline.cc,v 1.63 2006-10-25 18:36:26 tbretz Exp $
    33! --------------------------------------------------------------------------
    44!
    … …  
    165165const Float_t MExtractTimeAndChargeSpline::fgLoGainStretch    = 1.5;
    166166const Float_t MExtractTimeAndChargeSpline::fgOffsetLoGain     = 1.3;
    167 //const Float_t MExtractTimeAndChargeSpline::fgLoGainStartShift = -2.4;
    168167
    169168// --------------------------------------------------------------------------
    … …  
    191190  SetLoGainStretch();
    192191  SetOffsetLoGain(fgOffsetLoGain);
    193 //  SetLoGainStartShift(fgLoGainStartShift);
    194192
    195193  SetRiseTimeHiGain();
    … …  
    225223void MExtractTimeAndChargeSpline::SetChargeType( ExtractionType_t typ )
    226224{
    227 
    228225  fExtractionType = typ;
    229226
    230   if (fExtractionType == kAmplitude)
    231     {
    232       fNumHiGainSamples = 1.;
    233       fNumLoGainSamples = fLoGainLast ? 1. : 0.;
    234       fSqrtHiGainSamples = 1.;
    235       fSqrtLoGainSamples = 1.;
    236       fWindowSizeHiGain  = 1;
    237       fWindowSizeLoGain  = 1;
    238       fRiseTimeHiGain    = 0.5;
    239 
     227  InitArrays(fHiGainFirstDeriv.GetSize());
     228
     229  switch (fExtractionType)
     230  {
    240231      SetResolutionPerPheHiGain(0.053);
    241232      SetResolutionPerPheLoGain(0.016);
    242      
    243233      return;
    244     }
    245 
    246   if (fExtractionType == kIntegral)
    247     {
    248 
    249       fNumHiGainSamples  = fRiseTimeHiGain + fFallTimeHiGain;
    250       fNumLoGainSamples  = fLoGainLast ? fRiseTimeLoGain + fFallTimeLoGain : 0.;
    251 
    252       fSqrtHiGainSamples = TMath::Sqrt(fNumHiGainSamples);
    253       fSqrtLoGainSamples = TMath::Sqrt(fNumLoGainSamples);
    254       fWindowSizeHiGain  = TMath::Nint(fRiseTimeHiGain + fFallTimeHiGain);
    255       // to ensure that for the case: 1.5, the window size becomes: 2 (at any compiler)
    256       fWindowSizeLoGain  = TMath::Nint(TMath::Ceil((fRiseTimeLoGain + fFallTimeLoGain)*fLoGainStretch));
    257     }
    258 
    259     switch (fWindowSizeHiGain)
    260     {
    261     case 1:
    262       SetResolutionPerPheHiGain(0.041);
    263       break;
    264     case 2:
    265       SetResolutionPerPheHiGain(0.064);
    266       break;
    267     case 3:
    268     case 4:
    269       SetResolutionPerPheHiGain(0.050);
    270       break;
    271     case 5:
    272     case 6:
    273       SetResolutionPerPheHiGain(0.030);
    274       break;
    275     default:
    276         *fLog << warn << GetDescriptor() << ": Could not set the high-gain extractor resolution per phe for window size "
    277               << fWindowSizeHiGain << endl;
    278         break;
    279     }
    280 
    281     switch (fWindowSizeLoGain)
    282     {
    283     case 1:
    284     case 2:
    285       SetResolutionPerPheLoGain(0.005);
    286       break;
    287     case 3:
    288     case 4:
    289       SetResolutionPerPheLoGain(0.017);
    290       break;
    291     case 5:
    292     case 6:
    293     case 7:
    294       SetResolutionPerPheLoGain(0.005);
    295       break;
    296     case 8:
    297     case 9:
    298       SetResolutionPerPheLoGain(0.005);
    299       break;
    300     default:
    301         *fLog << warn << "Could not set the low-gain extractor resolution per phe for window size "
    302               << fWindowSizeLoGain << endl;
    303         break;
    304      }
     234
     235  case kIntegral:
     236      switch (fWindowSizeHiGain)
     237      {
     238      case 1:
     239          SetResolutionPerPheHiGain(0.041);
     240          break;
     241      case 2:
     242          SetResolutionPerPheHiGain(0.064);
     243          break;
     244      case 3:
     245      case 4:
     246          SetResolutionPerPheHiGain(0.050);
     247          break;
     248      case 5:
     249      case 6:
     250          SetResolutionPerPheHiGain(0.030);
     251          break;
     252      default:
     253          *fLog << warn << GetDescriptor() << ": Could not set the high-gain extractor resolution per phe for window size "
     254              << fWindowSizeHiGain << endl;
     255          break;
     256      }
     257
     258      switch (fWindowSizeLoGain)
     259      {
     260      case 1:
     261      case 2:
     262          SetResolutionPerPheLoGain(0.005);
     263          break;
     264      case 3:
     265      case 4:
     266          SetResolutionPerPheLoGain(0.017);
     267          break;
     268      case 5:
     269      case 6:
     270      case 7:
     271          SetResolutionPerPheLoGain(0.005);
     272          break;
     273      case 8:
     274      case 9:
     275          SetResolutionPerPheLoGain(0.005);
     276          break;
     277      default:
     278          *fLog << warn << "Could not set the low-gain extractor resolution per phe for window size "
     279              << fWindowSizeLoGain << endl;
     280          break;
     281      }
     282  }
    305283}
    306284
    … …  
    311289// Gets called in the ReInit() and initialized the arrays
    312290//
    313 Bool_t MExtractTimeAndChargeSpline::InitArrays()
    314 {
    315 
    316   Int_t range = fHiGainLast - fHiGainFirst + 1;// + fHiLoLast;
    317 
    318 //  fHiGainSignal     .Set(range);
    319   fHiGainFirstDeriv .Set(range);
    320   fHiGainSecondDeriv.Set(range);
    321 
    322   range = fLoGainLast - fLoGainFirst + 1;
    323 
    324 //  fLoGainSignal     .Set(range);
    325   fLoGainFirstDeriv .Set(range);
    326   fLoGainSecondDeriv.Set(range);
    327 
    328 //  fHiGainSignal     .Reset();
    329   fHiGainFirstDeriv .Reset();
    330   fHiGainSecondDeriv.Reset();
    331 
    332 //  fLoGainSignal     .Reset();
    333   fLoGainFirstDeriv .Reset();
    334   fLoGainSecondDeriv.Reset();
    335  
    336   if (fExtractionType == kAmplitude)
    337     {
    338       fNumHiGainSamples = 1.;
    339       fNumLoGainSamples = fLoGainLast ? 1. : 0.;
    340       fSqrtHiGainSamples = 1.;
    341       fSqrtLoGainSamples = 1.;
    342       fWindowSizeHiGain  = 1;
    343       fWindowSizeLoGain  = 1;
    344       fRiseTimeHiGain    = 0.5;
    345     }
    346 
    347   fRiseTimeLoGain    = fRiseTimeHiGain * fLoGainStretch;
    348   fFallTimeLoGain    = fFallTimeHiGain * fLoGainStretch;     
    349 
    350   if (fExtractionType == kIntegral)
    351     {
    352 
    353       fNumHiGainSamples  = fRiseTimeHiGain + fFallTimeHiGain;
    354       fNumLoGainSamples  = fLoGainLast ? fRiseTimeLoGain + fFallTimeLoGain : 0.;
    355       //      fNumLoGainSamples  *= 0.75;     
    356 
    357       fSqrtHiGainSamples = TMath::Sqrt(fNumHiGainSamples);
    358       fSqrtLoGainSamples = TMath::Sqrt(fNumLoGainSamples);
    359       fWindowSizeHiGain  = (Int_t)(fRiseTimeHiGain + fFallTimeHiGain);
    360       fWindowSizeLoGain  = (Int_t)(fRiseTimeLoGain + fFallTimeLoGain);
    361     }
    362 
    363   return kTRUE;
    364 
     291Bool_t MExtractTimeAndChargeSpline::InitArrays(Int_t n)
     292{
     293    // Initialize arrays to the maximum number of entries necessary
     294    fHiGainFirstDeriv .Set(n);
     295    fHiGainSecondDeriv.Set(n);
     296
     297    fLoGainFirstDeriv .Set(n);
     298    fLoGainSecondDeriv.Set(n);
     299
     300    fRiseTimeLoGain = fRiseTimeHiGain * fLoGainStretch;
     301    fFallTimeLoGain = fFallTimeHiGain * fLoGainStretch;
     302
     303    switch (fExtractionType)
     304    {
     305    case kAmplitude:
     306        fNumHiGainSamples  = 1.;
     307        fNumLoGainSamples  = fLoGainLast ? 1. : 0.;
     308        fSqrtHiGainSamples = 1.;
     309        fSqrtLoGainSamples = 1.;
     310        fWindowSizeHiGain  = 1;
     311        fWindowSizeLoGain  = 1;
     312        fRiseTimeHiGain    = 0.5;
     313        break;
     314
     315    case kIntegral:
     316        fNumHiGainSamples  = fRiseTimeHiGain + fFallTimeHiGain;
     317        fNumLoGainSamples  = fLoGainLast ? fRiseTimeLoGain + fFallTimeLoGain : 0.;
     318        fSqrtHiGainSamples = TMath::Sqrt(fNumHiGainSamples);
     319        fSqrtLoGainSamples = TMath::Sqrt(fNumLoGainSamples);
     320        fWindowSizeHiGain  = TMath::CeilNint(fRiseTimeHiGain + fFallTimeHiGain);
     321        fWindowSizeLoGain  = TMath::CeilNint(fRiseTimeLoGain + fFallTimeLoGain);
     322        break;
     323    }
     324
     325    return kTRUE;
    365326}
    366327
    … …  
    368329                                                           Float_t &sum, Float_t &dsum,
    369330                                                           Float_t &time, Float_t &dtime,
    370                                                            Byte_t sat, Int_t maxpos)
    371 {
    372     //    const Int_t range = fHiGainLast - fHiGainFirst + 1 + fHiLoLast;
    373 
     331                                                           Byte_t sat, Int_t maxpos) const
     332{
    374333    // Do some handling if maxpos is last slice!
    375 
    376334    MExtralgoSpline s(ptr, num, fHiGainFirstDeriv.GetArray(), fHiGainSecondDeriv.GetArray());
    377335
    378336    s.SetRiseFallTime(fRiseTimeHiGain, fFallTimeHiGain);
    379 //    s.SetResolution(fResolution);
    380337
    381338    if (IsNoiseCalculation())
    382339    {
    383 //        if (fRandomIter == int(1./fResolution))
    384 //            fRandomIter = 0;
    385 
    386         sum = s.ExtractNoise(/*fRandomIter*/);
    387 //        fRandomIter++;
     340        sum = s.ExtractNoise();
    388341        return;
    389342    }
    … …  
    397350                                                           Float_t &sum,  Float_t &dsum,
    398351                                                           Float_t &time, Float_t &dtime,
    399                                                            Byte_t sat, Int_t maxpos)
     352                                                           Byte_t sat, Int_t maxpos) const
    400353{
    401354    MExtralgoSpline s(ptr, num, fLoGainFirstDeriv.GetArray(), fLoGainSecondDeriv.GetArray());
    402355
    403356    s.SetRiseFallTime(fRiseTimeLoGain, fFallTimeLoGain);
    404 //    s.SetResolution(fResolution);
    405357
    406358    if (IsNoiseCalculation())
    407359    {
    408 //        if (fRandomIter == int(1./fResolution))
    409 //            fRandomIter = 0;
    410 
    411         sum = s.ExtractNoise(/*fRandomIter*/);
     360        sum = s.ExtractNoise();
    412361        return;
    413362    }
    … …  
    417366    s.GetSignal(sum, dsum);
    418367}
    419 
    420 // --------------------------------------------------------------------------
    421 //
    422 // Calculates the arrival time and charge for each pixel
    423 //
    424 /*
    425 void MExtractTimeAndChargeSpline::FindTimeAndChargeHiGain(Byte_t *first, Byte_t *logain, Float_t &sum, Float_t &dsum,
    426                                                           Float_t &time, Float_t &dtime,
    427                                                           Byte_t &sat, const MPedestalPix &ped, const Bool_t abflag)
    428 {
    429     Int_t range = fHiGainLast - fHiGainFirst + 1;
    430     const Byte_t *end = first + range;
    431     Byte_t       *p  = first;
    432 
    433     const Float_t pedes  = ped.GetPedestal();
    434     const Float_t ABoffs = ped.GetPedestalABoffset();
    435 
    436     const Float_t pedmean[2] = { pedes + ABoffs, pedes - ABoffs };
    437 
    438     Byte_t maxcont = 0;
    439     Int_t  maxpos  = 0;
    440     sat = 0;
    441 
    442     //
    443     // Check for saturation in all other slices
    444     //
    445     Int_t ids = fHiGainFirst;
    446     Float_t *sample = fHiGainSignal.GetArray();
    447     while (p<end)
    448     {
    449         *sample++ = (Float_t)*p - pedmean[(ids++ + abflag) & 0x1];
    450 
    451         if (*p > maxcont)
    452         {
    453             maxpos = ids-fHiGainFirst-1;
    454             // range-fWindowSizeHiGain+1 == fHiLoLast isn't it?
    455             //if (maxpos > 1 && maxpos < (range - fWindowSizeHiGain + 1))
    456             if (maxpos > 1 && maxpos < (range-1*//* - fWindowSizeHiGain + 1*//*))
    457                 maxcont = *p;
    458         }
    459 
    460         if (*p++ >= fSaturationLimit)
    461             if (!sat)
    462                 sat = ids-fHiGainFirst;
    463     }
    464 
    465     if (fHiLoLast != 0)
    466     {
    467         end = logain + fHiLoLast;
    468 
    469         while (logain<end)
    470         {
    471             *sample++ = (Float_t)*logain - pedmean[(ids++ + abflag) & 0x1];
    472 
    473             if (*logain > maxcont)
    474             {
    475                 maxpos = ids-fHiGainFirst-1;
    476 
    477                 // FIXME!!! MUST BE THERE!
    478                 // range-fWindowSizeHiGain+1 == fHiLoLast isn't it?
    479                 //if (maxpos > 1 && maxpos < (range - fWindowSizeHiGain + 1))
    480                 //    maxcont = *logain;
    481             }
    482 
    483             if (*logain++ >= fSaturationLimit)
    484                 if (!sat)
    485                     sat = ids-fHiGainFirst;
    486 
    487             range++;
    488         }
    489     }
    490 
    491     FindTimeAndChargeHiGain2(fHiGainSignal.GetArray(), range,
    492                              sum, dsum, time, dtime,
    493                              sat, maxpos);
    494 }
    495 
    496 // --------------------------------------------------------------------------
    497 //
    498 // Calculates the arrival time and charge for each pixel
    499 //
    500 void MExtractTimeAndChargeSpline::FindTimeAndChargeLoGain(Byte_t *first, Float_t &sum, Float_t &dsum,
    501                                                           Float_t &time, Float_t &dtime,
    502                                                           Byte_t &sat, const MPedestalPix &ped, const Bool_t abflag)
    503 {
    504     Int_t range = fLoGainLast - fLoGainFirst + 1;
    505     const Byte_t *end = first + range;
    506     Byte_t       *p  = first;
    507 
    508     const Float_t pedes  = ped.GetPedestal();
    509     const Float_t ABoffs = ped.GetPedestalABoffset();
    510 
    511     const Float_t pedmean[2] = { pedes + ABoffs, pedes - ABoffs };
    512 
    513     Int_t  maxpos = 0;
    514     Int_t  max    = -9999;
    515 
    516     //
    517     // Check for saturation in all other slices
    518     //
    519     Int_t    ids    = fLoGainFirst;
    520     Float_t *sample = fLoGainSignal.GetArray();
    521     while (p<end)
    522     {
    523         *sample++ = (Float_t)*p - pedmean[(ids++ + abflag) & 0x1];
    524 
    525         if (*p > max)
    526         {
    527             maxpos = ids-fLoGainFirst-1;
    528             max    = *p;
    529         }
    530 
    531         if (*p++ >= fSaturationLimit)
    532             sat++;
    533     }
    534 
    535     FindTimeAndChargeLoGain2(fLoGainSignal.GetArray(), range,
    536                              sum, dsum, time, dtime,
    537                              sat, maxpos);
    538 }
    539 */
    540368
    541369// --------------------------------------------------------------------------

  • trunk/MagicSoft/Mars/msignal/MExtractTimeAndChargeSpline.h

    r8154 r8165  
    4242    Float_t fLoGainStretch;                  // The stretch of the low-gain w.r.t. the high-gain pulse
    4343
    44     Int_t   fRandomIter;                     //! Counter used to randomize weights for noise calculation
     44//    Int_t   fRandomIter;                     //! Counter used to randomize weights for noise calculation
    4545
    4646    Int_t   ReadEnv(const TEnv &env, TString prefix, Bool_t print);
    47     Bool_t  InitArrays();
     47    Bool_t  InitArrays(Int_t n);
    4848
    4949public:
    … …  
    9494    void FindTimeAndChargeHiGain2(const Float_t *firstused, Int_t num, Float_t &sum, Float_t &dsum,
    9595                                  Float_t &time, Float_t &dtime,
    96                                   Byte_t sat, Int_t maxpos);
     96                                  Byte_t sat, Int_t maxpos) const;
    9797
    9898    void FindTimeAndChargeLoGain2(const Float_t *firstused, Int_t num, Float_t &sum,  Float_t &dsum,
    9999                                  Float_t &time, Float_t &dtime,
    100                                   Byte_t sat, Int_t maxpos);
     100                                  Byte_t sat, Int_t maxpos) const;
    101101
    102102    ClassDef(MExtractTimeAndChargeSpline, 4)   // Task to Extract Arrival Times and Charges using a Cubic Spline
Note: See TracChangeset for help on using the changeset viewer.