Changeset 4601 for trunk/MagicSoft/Mars/mpedestal
- Timestamp:
- 08/12/04 16:41:37 (21 years ago)
- Location:
- trunk/MagicSoft/Mars/mpedestal
- Files:
-
- 4 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/MagicSoft/Mars/mpedestal/MPedCalcFromLoGain.cc
r4565 r4601 114 114 // fLoGainWindowSize = fgLoGainWindowSize = 0 115 115 // 116 // Variables: 117 // fgHiGainFirst; First FADC slice Hi-Gain (currently set to: 3) 118 // fgHiGainLast: Last FADC slice Hi-Gain (currently set to: 14) 119 // fgLoGainFirst: First FADC slice Lo-Gain (currently set to: 3) 120 // fgLoGainLast: Last FADC slice Lo-Gain (currently set to: 14) 121 // fgHiGainWindowSize: The extraction window Hi-Gain 122 // fgLoGainWindowSize: The extraction window Lo-Gain 123 // fgMaxHiGainVar: The maximum difference between the highest and lowest slice 124 // in the high gain window allowed in order to use low gain 125 // 116 126 // Input Containers: 117 127 // MRawEvtData … … 201 211 fRunHeader = NULL; 202 212 fPedestals = NULL; 213 214 // If the size is yet set, set the size 215 if (fSumx.GetSize()>0) 216 { 217 // Reset contents of arrays. 218 fSumx.Reset(); 219 fSumx2.Reset(); 220 fSumAB0.Reset(); 221 fSumAB1.Reset(); 222 fNumEventsUsed.Reset(); 223 fTotalCounter.Reset(); 224 } 203 225 } 204 226 … … 297 319 // - MPedestalCam 298 320 // 299 Int_t MPedCalcFromLoGain::PreProcess( MParList *pList ) 300 { 301 Clear(); 302 303 fRawEvt = (MRawEvtData*)pList->FindObject("MRawEvtData"); 304 if (!fRawEvt) 305 { 306 *fLog << err << "MRawEvtData not found... aborting." << endl; 307 return kFALSE; 308 } 309 310 fRunHeader = (MRawRunHeader*)pList->FindObject(AddSerialNumber("MRawRunHeader")); 311 if (!fRunHeader) 312 { 313 *fLog << err << AddSerialNumber("MRawRunHeader") << " not found... aborting." << endl; 314 return kFALSE; 315 } 316 317 fGeom = (MGeomCam*)pList->FindObject("MGeomCam"); 318 if (!fGeom) 319 { 320 *fLog << err << "MGeomCam not found... aborting." << endl; 321 return kFALSE; 322 } 323 324 *fLog << " searching for the container " << fPedContainerName << endl; 325 326 fPedestals = (MPedestalCam*)pList->FindCreateObj( AddSerialNumber("MPedestalCam"),fPedContainerName); 327 if (!fPedestals) 328 { 329 *fLog << err << fPedContainerName << " can not be created" << endl; 330 return kFALSE; 331 } 332 333 334 return kTRUE; 321 Int_t MPedCalcFromLoGain::PreProcess(MParList *pList) 322 { 323 Clear(); 324 325 fRawEvt = (MRawEvtData*)pList->FindObject("MRawEvtData"); 326 if (!fRawEvt) 327 { 328 *fLog << err << "MRawEvtData not found... aborting." << endl; 329 return kFALSE; 330 } 331 332 fRunHeader = (MRawRunHeader*)pList->FindObject(AddSerialNumber("MRawRunHeader")); 333 if (!fRunHeader) 334 { 335 *fLog << err << AddSerialNumber("MRawRunHeader") << " not found... aborting." << endl; 336 return kFALSE; 337 } 338 339 fGeom = (MGeomCam*)pList->FindObject("MGeomCam"); 340 if (!fGeom) 341 { 342 *fLog << err << "MGeomCam not found... aborting." << endl; 343 return kFALSE; 344 } 345 346 fPedestals = (MPedestalCam*)pList->FindCreateObj("MPedestalCam", AddSerialNumber(fPedContainerName)); 347 if (!fPedestals) 348 return kFALSE; 349 350 return kTRUE; 335 351 } 336 352 … … 353 369 Bool_t MPedCalcFromLoGain::ReInit(MParList *pList) 354 370 { 355 356 371 Int_t lastdesired = (Int_t)fLoGainLast; 357 372 Int_t lastavailable = (Int_t)fRunHeader->GetNumSamplesLoGain()-1; … … 413 428 414 429 415 Int_t npixels = fPedestals->GetSize(); 416 430 // If the size is not yet set, set the size 417 431 if (fSumx.GetSize()==0) 418 { 432 { 433 const Int_t npixels = fPedestals->GetSize(); 434 419 435 fSumx. Set(npixels); 420 436 fSumx2.Set(npixels); … … 423 439 fNumEventsUsed.Set(npixels); 424 440 fTotalCounter.Set(npixels); 425 441 442 // Reset contents of arrays. 426 443 fSumx.Reset(); 427 444 fSumx2.Reset(); … … 430 447 fNumEventsUsed.Reset(); 431 448 fTotalCounter.Reset(); 432 433 434 if (fWindowSizeHiGain == 0 && fWindowSizeLoGain ==0)435 436 *fLog << err << GetDescriptor() 449 } 450 451 if (fWindowSizeHiGain==0 && fWindowSizeLoGain==0) 452 { 453 *fLog << err << GetDescriptor() 437 454 << ": Number of extracted Slices is 0, abort ... " << endl; 438 455 return kFALSE; 439 } 440 441 456 } 457 442 458 *fLog << endl; 443 459 *fLog << inf << GetDescriptor() << ": Taking " << (int)fWindowSizeHiGain … … 449 465 } 450 466 467 void MPedCalcFromLoGain::Calc(ULong_t n, UInt_t idx) 468 { 469 const ULong_t nsamplestot = n*fWindowSizeLoGain; 470 471 const Float_t sum = fSumx.At(idx); 472 const Float_t sum2 = fSumx2.At(idx); 473 const Float_t ped = sum/nsamplestot; 474 475 // 1. Calculate the Variance of the sums: 476 Float_t var = (sum2-sum*sum/n)/(n-1.); 477 478 // 2. Scale the variance to the number of slices: 479 var /= (Float_t)(fWindowSizeLoGain); 480 481 // 3. Calculate the RMS from the Variance: 482 const Float_t rms = var<0 ? 0 : TMath::Sqrt(var); 483 484 // 4. Calculate the amplitude of the 150MHz "AB" noise 485 const Float_t abOffs = (fSumAB0[idx] - fSumAB1[idx]) / nsamplestot; 486 487 (*fPedestals)[idx].Set(ped, rms, abOffs, n); 488 489 fTotalCounter[idx]++; 490 } 491 451 492 // -------------------------------------------------------------------------- 452 493 // … … 457 498 Int_t MPedCalcFromLoGain::Process() 458 499 { 459 MRawEvtPixelIter pixel(fRawEvt); 460 461 while (pixel.Next()) { 500 MRawEvtPixelIter pixel(fRawEvt); 501 502 while (pixel.Next()) 503 { 504 const UInt_t idx = pixel.GetPixelId(); 505 506 Byte_t *ptr = pixel.GetHiGainSamples() + fHiGainFirst; 507 Byte_t *end = ptr + fWindowSizeHiGain; 508 509 UInt_t sum = 0; 510 UInt_t sqr = 0; 511 512 UInt_t max = 0; 513 UInt_t min = 255; 462 514 463 const UInt_t idx = pixel.GetPixelId(); 464 465 Byte_t *ptr = pixel.GetHiGainSamples() + fHiGainFirst; 466 Byte_t *end = ptr + fWindowSizeHiGain; 467 468 UInt_t sum = 0; 469 UInt_t sqr = 0; 470 471 UInt_t max = 0; 472 UInt_t min = 255; 473 474 // Find the maximum and minimum signal per slice in the high gain window 475 do { 476 if (*ptr > max) { 477 max = *ptr; 478 } 479 if (*ptr < min) { 480 min = *ptr; 481 } 482 } while (++ptr != end); 483 484 // If the maximum in the high gain window is smaller than 485 if ((max-min < fMaxHiGainVar) && (max < 255)) { 486 487 Byte_t *firstSlice = pixel.GetLoGainSamples() + fLoGainFirst; 488 Byte_t *lastSlice = firstSlice + fWindowSizeLoGain; 489 490 Byte_t *slice = firstSlice; 491 do { 492 sum += *slice; 493 sqr += *slice * *slice; 494 } while (++slice != lastSlice); 495 496 const Float_t msum = (Float_t)sum; 497 const Float_t sqrsum = msum*msum; 498 499 fSumx[idx] += msum; 500 fSumx2[idx] += sqrsum; 501 fNumEventsUsed[idx]++; 502 503 // Calculate the amplitude of the 150MHz "AB" noise 504 505 Int_t abFlag = (fRunHeader->GetNumSamplesHiGain() 506 + fLoGainFirst 507 + pixel.HasABFlag()) & 0x1; 508 for (Int_t islice=0; islice<fWindowSizeLoGain; islice+=2) { 509 Int_t sliceAB0 = islice + abFlag; 510 Int_t sliceAB1 = islice + 1 - abFlag; 511 fSumAB0[idx] += firstSlice[sliceAB0]; 512 fSumAB1[idx] += firstSlice[sliceAB1]; 513 } 514 515 if (fPedestalUpdate && (fNumEventsUsed[idx] >= fNumEventsDump)) { 516 517 const ULong_t n = fNumEventsDump; 518 519 const ULong_t nsamplestot = n*fWindowSizeLoGain; 520 521 const Float_t sum = fSumx.At(idx); 522 const Float_t sum2 = fSumx2.At(idx); 523 const Float_t ped = sum/(nsamplestot); 524 525 // 1. Calculate the Variance of the sums: 526 Float_t var = (sum2-sum*sum/n)/(n-1.); 527 528 // 2. Scale the variance to the number of slices: 529 var /= (Float_t)(fWindowSizeLoGain); 530 531 // 3. Calculate the RMS from the Variance: 532 Float_t rms = var < 0 ? 0. : TMath::Sqrt(var); 533 534 // 4. Calculate the amplitude of the 150MHz "AB" noise 535 Float_t abOffs = (fSumAB0[idx] - fSumAB1[idx]) / nsamplestot; 536 537 (*fPedestals)[idx].Set(ped, rms, abOffs, n); 538 539 fTotalCounter[idx]++; 540 fNumEventsUsed[idx]=0; 541 fSumx[idx]=0; 542 fSumx2[idx]=0; 543 fSumAB0[idx]=0; 544 fSumAB1[idx]=0; 545 } 546 } 547 } 548 549 fPedestals->SetReadyToSave(); 550 551 return kTRUE; 515 // Find the maximum and minimum signal per slice in the high gain window 516 do { 517 if (*ptr > max) { 518 max = *ptr; 519 } 520 if (*ptr < min) { 521 min = *ptr; 522 } 523 } while (++ptr != end); 524 525 // If the maximum in the high gain window is smaller than 526 if (max-min>=fMaxHiGainVar || max>=255) 527 continue; 528 529 Byte_t *firstSlice = pixel.GetLoGainSamples() + fLoGainFirst; 530 Byte_t *lastSlice = firstSlice + fWindowSizeLoGain; 531 532 Byte_t *slice = firstSlice; 533 do { 534 sum += *slice; 535 sqr += *slice * *slice; 536 } while (++slice != lastSlice); 537 538 const Float_t msum = (Float_t)sum; 539 const Float_t sqrsum = msum*msum; 540 541 fSumx[idx] += msum; 542 fSumx2[idx] += sqrsum; 543 fNumEventsUsed[idx]++; 544 545 // Calculate the amplitude of the 150MHz "AB" noise 546 547 Int_t abFlag = (fRunHeader->GetNumSamplesHiGain() 548 + fLoGainFirst 549 + pixel.HasABFlag()) & 0x1; 550 for (Int_t islice=0; islice<fWindowSizeLoGain; islice+=2) 551 { 552 Int_t sliceAB0 = islice + abFlag; 553 Int_t sliceAB1 = islice - abFlag + 1; 554 fSumAB0[idx] += firstSlice[sliceAB0]; 555 fSumAB1[idx] += firstSlice[sliceAB1]; 556 } 557 558 if (!fPedestalUpdate || fNumEventsUsed[idx]<fNumEventsDump) 559 continue; 560 561 Calc(fNumEventsDump, idx); 562 563 fNumEventsUsed[idx]=0; 564 fSumx[idx]=0; 565 fSumx2[idx]=0; 566 fSumAB0[idx]=0; 567 fSumAB1[idx]=0; 568 } 569 570 if (fPedestalUpdate) 571 fPedestals->SetReadyToSave(); 572 573 return kTRUE; 552 574 } 553 575 … … 558 580 Int_t MPedCalcFromLoGain::PostProcess() 559 581 { 560 561 // Compute pedestals and rms from the whole run 562 563 if (!fPedestalUpdate) { 564 565 MRawEvtPixelIter pixel(fRawEvt); 566 567 while (pixel.Next()) { 568 569 const Int_t idx = pixel.GetPixelId(); 570 571 const ULong_t n = fNumEventsUsed[idx]; 572 573 if (n < 2) 574 continue; 575 576 const ULong_t nsamplestot = n*fWindowSizeLoGain; 577 578 const Float_t sum = fSumx.At(idx); 579 const Float_t sum2 = fSumx2.At(idx); 580 const Float_t ped = sum/(nsamplestot); 581 582 // 1. Calculate the Variance of the sums: 583 Float_t var = (sum2-sum*sum/n)/(n-1.); 584 585 // 2. Scale the variance to the number of slices: 586 var /= (Float_t)(fWindowSizeLoGain); 587 588 // 3. Calculate the RMS from the Variance: 589 Float_t rms = var < 0 ? 0. : TMath::Sqrt(var); 590 591 // 4. Calculate the amplitude of the 150MHz "AB" noise 592 Float_t abOffs = (fSumAB0[idx] - fSumAB1[idx]) / nsamplestot; 593 594 (*fPedestals)[idx].Set(ped, rms, abOffs, n); 595 596 fTotalCounter[idx]++; 582 // Compute pedestals and rms from the whole run 583 if (fPedestalUpdate) 584 return kTRUE; 585 586 *fLog << flush << inf << "Calculating Pedestals..." << flush; 587 588 const Int_t npix = fGeom->GetNumPixels(); 589 for (Int_t idx=0; idx<npix; idx++) 590 { 591 const ULong_t n = fNumEventsUsed[idx]; 592 if (n>1) 593 Calc(n, idx); 597 594 } 598 595 599 596 fPedestals->SetReadyToSave(); 600 } 601 602 fSumx.Reset(); 603 fSumx2.Reset(); 604 fSumAB0.Reset(); 605 fSumAB1.Reset(); 606 fNumEventsUsed.Reset(); 607 fTotalCounter.Reset(); 608 609 return kTRUE; 610 } 597 return kTRUE; 598 } 599 600 Int_t MPedCalcFromLoGain::ReadEnv(const TEnv &env, TString prefix, Bool_t print) 601 { 602 if (MExtractor::ReadEnv(env, prefix, print)==kERROR) 603 return kERROR; 604 605 Bool_t rc=kFALSE; 606 607 Byte_t hw = fWindowSizeHiGain; 608 Byte_t lw = fWindowSizeLoGain; 609 if (IsEnvDefined(env, prefix, "WindowSizeHiGain", print)) 610 { 611 hw = GetEnvValue(env, prefix, "WindowSizeHiGain", hw); 612 rc = kTRUE; 613 } 614 if (IsEnvDefined(env, prefix, "WindowSizeLoGain", print)) 615 { 616 lw = GetEnvValue(env, prefix, "WindowSizeLoGain", lw); 617 rc = kTRUE; 618 } 619 620 if (rc) 621 SetWindowSize(hw, lw); 622 623 Int_t num = fNumEventsDump; 624 if (IsEnvDefined(env, prefix, "NumEventsDump", print)) 625 { 626 num = GetEnvValue(env, prefix, "NumEventsDump", num); 627 rc = kTRUE; 628 } 629 SetDumpEvents(num); 630 631 Byte_t max = fMaxHiGainVar; 632 if (IsEnvDefined(env, prefix, "MaxHiGainVar", print)) 633 { 634 max = GetEnvValue(env, prefix, "MaxHiGainVar", max); 635 rc = kTRUE; 636 } 637 SetMaxHiGainVar(max); 638 639 Bool_t upd = fPedestalUpdate; 640 if (IsEnvDefined(env, prefix, "PedestalUpdate", print)) 641 { 642 upd = GetEnvValue(env, prefix, "PedestalUpdate", upd); 643 rc = kTRUE; 644 } 645 SetPedestalUpdate(upd); 646 647 return rc; 648 } -
trunk/MagicSoft/Mars/mpedestal/MPedCalcFromLoGain.h
r4565 r4601 19 19 class MPedCalcFromLoGain : public MExtractor 20 20 { 21 static const Byte_t fgHiGainFirst; // First FADC slice Hi-Gain (currently set to: 3) 22 static const Byte_t fgHiGainLast; // Last FADC slice Hi-Gain (currently set to: 14) 23 static const Byte_t fgLoGainFirst; // First FADC slice Lo-Gain (currently set to: 3) 24 static const Byte_t fgLoGainLast; // Last FADC slice Lo-Gain (currently set to: 14) 25 static const Byte_t fgHiGainWindowSize; // The extraction window Hi-Gain 26 static const Byte_t fgLoGainWindowSize; // The extraction window Lo-Gain 27 static const Byte_t fgMaxHiGainVar; // The maximum difference between the highest and lowest slice 28 // in the high gain window allowed in order to use low gain 21 static const Byte_t fgHiGainFirst; // First FADC slice Hi-Gain (currently set to: 3) 22 static const Byte_t fgHiGainLast; // Last FADC slice Hi-Gain (currently set to: 14) 23 static const Byte_t fgLoGainFirst; // First FADC slice Lo-Gain (currently set to: 3) 24 static const Byte_t fgLoGainLast; // Last FADC slice Lo-Gain (currently set to: 14) 25 static const Byte_t fgHiGainWindowSize; // The extraction window Hi-Gain 26 static const Byte_t fgLoGainWindowSize; // The extraction window Lo-Gain 27 static const Byte_t fgMaxHiGainVar; // The maximum difference between the highest and lowest slice 29 28 30 Int_t fNumEventsDump;// Number of event after which MPedestalCam gets updated29 Int_t fNumEventsDump; // Number of event after which MPedestalCam gets updated 31 30 32 Byte_t fMaxHiGainVar;33 Byte_t fWindowSizeHiGain;// Number of Hi Gain slices in window34 Byte_t fWindowSizeLoGain;// Number of Lo Gain slices in window31 Byte_t fMaxHiGainVar; 32 Byte_t fWindowSizeHiGain; // Number of Hi Gain slices in window 33 Byte_t fWindowSizeLoGain; // Number of Lo Gain slices in window 35 34 36 Bool_t fPedestalUpdate;35 Bool_t fPedestalUpdate; 37 36 38 MGeomCam *fGeom; // Camera geometry 39 TString fPedContainerName; // name of the 'MPedestalCam' container 40 41 TArrayI fNumEventsUsed; // Number of events used for pedestal calc for each pixel 42 TArrayI fTotalCounter; // Counter for dumping values to Pedestal Container 43 TArrayD fSumx; // sum of values 44 TArrayD fSumx2; // sum of squared values 45 TArrayD fSumAB0; // sum of ABFlag=0 slices 46 TArrayD fSumAB1; // sum of ABFlag=1 slices 37 MGeomCam *fGeom; // Camera geometry 38 TString fPedContainerName; // name of the 'MPedestalCam' container 47 39 48 49 Int_t PreProcess (MParList *pList); 50 Bool_t ReInit (MParList *pList); 51 Int_t Process (); 52 Int_t PostProcess(); 53 40 TArrayI fNumEventsUsed; // Number of events used for pedestal calc for each pixel 41 TArrayI fTotalCounter; // Counter for dumping values to Pedestal Container 42 TArrayD fSumx; // sum of values 43 TArrayD fSumx2; // sum of squared values 44 TArrayD fSumAB0; // sum of ABFlag=0 slices 45 TArrayD fSumAB1; // sum of ABFlag=1 slices 46 47 // MParContainer 48 Int_t PreProcess (MParList *pList); 49 Bool_t ReInit (MParList *pList); 50 Int_t Process (); 51 Int_t PostProcess(); 52 Int_t ReadEnv(const TEnv &env, TString prefix, Bool_t print); 53 54 // Calculation 55 void Calc(ULong_t n, UInt_t idx); 56 54 57 public: 55 MPedCalcFromLoGain(const char *name=NULL, const char *title=NULL); 56 57 void Clear(const Option_t *o=""); 58 void SetRange(Byte_t hifirst=0, Byte_t hilast=0, Byte_t lofirst=0, Byte_t lolast=0); 59 void SetWindowSize(Byte_t windowh=0, Byte_t windowl=0); 60 void SetMaxHiGainVar(Byte_t maxvar=0); 61 void SetDumpEvents(UInt_t dumpevents = 0) {fNumEventsDump = dumpevents;} 62 void SetPedestalUpdate(Bool_t pedupdate) {fPedestalUpdate = pedupdate;} 63 64 void SetPedContainerName(const char *name) { fPedContainerName = name; } 58 MPedCalcFromLoGain(const char *name=NULL, const char *title=NULL); 65 59 66 TArrayI *GetNumEventsUsed() {return &fNumEventsUsed;}; 67 68 ClassDef(MPedCalcFromLoGain, 0) // Task to calculate pedestals from pedestal runs raw data 60 // TObject 61 void Clear(const Option_t *o=""); 62 63 // Setter 64 void SetRange(Byte_t hifirst=0, Byte_t hilast=0, Byte_t lofirst=0, Byte_t lolast=0); 65 void SetWindowSize(Byte_t windowh=0, Byte_t windowl=0); 66 void SetMaxHiGainVar(Byte_t maxvar=0); 67 void SetDumpEvents(UInt_t dumpevents = 0) {fNumEventsDump = dumpevents;} 68 void SetPedestalUpdate(Bool_t pedupdate) {fPedestalUpdate = pedupdate;} 69 70 void SetPedContainerName(const char *name) { fPedContainerName = name; } 71 72 // Getter 73 TArrayI *GetNumEventsUsed() { return &fNumEventsUsed; } 74 75 ClassDef(MPedCalcFromLoGain, 0) // Task to calculate pedestals from pedestal runs raw data 69 76 }; 70 77 -
trunk/MagicSoft/Mars/mpedestal/MPedCalcPedRun.cc
r4556 r4601 443 443 while (pixel.Next()) 444 444 { 445 446 445 const UInt_t idx = pixel.GetPixelId(); 447 446 const UInt_t aidx = (*fGeom)[idx].GetAidx(); … … 620 619 return kTRUE; 621 620 } 621 622 Int_t MPedCalcPedRun::ReadEnv(const TEnv &env, TString prefix, Bool_t print) 623 { 624 if (MExtractor::ReadEnv(env, prefix, print)==kERROR) 625 return kERROR; 626 627 Byte_t hw = fWindowSizeHiGain; 628 Byte_t lw = fWindowSizeLoGain; 629 630 Bool_t rc = kFALSE; 631 632 if (IsEnvDefined(env, prefix, "WindowSizeHiGain", print)) 633 { 634 hw = GetEnvValue(env, prefix, "WindowSizeHiGain", hw); 635 rc=kTRUE; 636 } 637 638 if (IsEnvDefined(env, prefix, "WindowSizeLoGain", print)) 639 { 640 lw = GetEnvValue(env, prefix, "WindowSizeLoGain", lw); 641 rc=kTRUE; 642 } 643 644 if (rc) 645 SetWindowSize(hw, lw); 646 647 return rc; 648 } -
trunk/MagicSoft/Mars/mpedestal/MPedCalcPedRun.h
r4556 r4601 44 44 Int_t Process (); 45 45 Int_t PostProcess(); 46 Int_t ReadEnv(const TEnv &env, TString prefix, Bool_t print); 46 47 47 48 public:
Note:
See TracChangeset
for help on using the changeset viewer.