Changeset 3636
- Timestamp:
- 04/03/04 17:27:50 (21 years ago)
- Location:
- trunk/MagicSoft/Mars
- Files:
-
- 17 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/MagicSoft/Mars/Changelog
r3635 r3636 18 18 19 19 -*-*- END OF LINE -*-*- 20 21 2004/04/02: Markus Gaug 22 * mcalib/MHGausEvents.[h,cc] 23 * mcalib/MHCalibrationChargePix.[h,cc] 24 * mcalib/MHCalibrationChargeHiGainPix.[h,cc] 25 * mcalib/MHCalibrationChargeLoGainPix.[h,cc] 26 * mcalib/MHCalibrationRelTimePix.[h,cc] 27 - put fPixId, fPickup, fPickupLimt, CountPickup(), RepeatFit() and 28 ChangeHistId() into MHGausEvents (before in the derived classes) 29 - put fChargeNbins, fChargeFirst, fChargeLast, 30 - put fRelTimeNbins, fRelTimeFirst, fRelTimeLast together 31 into MHGausEvents as fNbins, fFirst and fLast 32 33 * mcalib/MHCalibrationRelTimePix.[h,cc] 34 - remove Renormalization to time slices. Need to think about 35 more direct way to implement 36 37 * mcalib/MHCalibrationCam.[h,cc] 38 * mcalib/MHCalibrationChargeCam.[h,cc] 39 * mcalib/MHCalibrationRelTimeCam.[h,cc] 40 - put most of the functionality into the base class MHCalibrationCam 41 - derived classes overload the functions SetupHists, ReInitHists, 42 FillHists, FinalizeHists and FinalizeBadPixels. 43 - functions FitHiGainArrays, FitLoGainArrays, FitHiGainHists, 44 FitLoGainHists and InitHists can be used from base class. 20 45 21 46 -
trunk/MagicSoft/Mars/mcalib/MHCalibrationCam.cc
r3631 r3636 53 53 #include "MLogManip.h" 54 54 55 #include "MCalibrationPix.h" 56 #include "MCalibrationCam.h" 57 55 58 #include "MHGausEvents.h" 56 59 60 #include "MBadPixelsPix.h" 61 #include "MBadPixelsCam.h" 62 57 63 #include "MGeomCam.h" 64 #include "MGeomPix.h" 65 66 #include "MParList.h" 58 67 59 68 ClassImp(MHCalibrationCam); … … 61 70 using namespace std; 62 71 72 const Int_t MHCalibrationCam::fgAverageNbins = 2000; 63 73 const Int_t MHCalibrationCam::fgPulserFrequency = 500; 64 74 // -------------------------------------------------------------------------- 65 75 // 66 76 // Default Constructor. 77 // 78 // Sets: 79 // - all pointers to NULL 67 80 // 68 81 // Initializes and sets owner of: … … 75 88 // 76 89 MHCalibrationCam::MHCalibrationCam(const char *name, const char *title) 90 : fGeom(NULL), fBadPixels(NULL), fCam(NULL) 77 91 { 78 92 fName = name ? name : "MHCalibrationCam"; … … 97 111 fAverageLoGainSectors->SetOwner(); 98 112 113 SetAverageNbins(); 99 114 SetPulserFrequency(); 100 115 } … … 306 321 // -------------------------------------------------------------------------- 307 322 // 323 // Gets the pointers to: 324 // - MGeomCam 325 // 326 // Calls SetupHists(const MParList *pList) 327 // 328 // Calls Delete-Function of: 329 // - MHCalibrationCam::fHiGainArray, MHCalibrationCam::fLoGainArray 330 // - MHCalibrationCam::fAverageHiGainAreas, MHCalibrationCam::fAverageLoGainAreas 331 // - MHCalibrationCam::fAverageHiGainSectors, MHCalibrationCam::fAverageLoGainSectors 332 // 333 Bool_t MHCalibrationCam::SetupFill(const MParList *pList) 334 { 335 336 fGeom = (MGeomCam*)pList->FindObject("MGeomCam"); 337 if (!fGeom) 338 { 339 *fLog << err << "MGeomCam not found... aborting." << endl; 340 return kFALSE; 341 } 342 343 fHiGainArray->Delete(); 344 fLoGainArray->Delete(); 345 346 fAverageHiGainAreas->Delete(); 347 fAverageLoGainAreas->Delete(); 348 349 fAverageHiGainSectors->Delete(); 350 fAverageLoGainSectors->Delete(); 351 352 return SetupHists(pList); 353 } 354 355 356 Bool_t MHCalibrationCam::SetupHists(const MParList *pList) 357 { 358 return kTRUE; 359 } 360 361 // -------------------------------------------------------------------------- 362 // 363 // Gets or creates the pointers to: 364 // - MBadPixelsCam 365 // 366 // Searches pointer to: 367 // - MArrivalTimeCam 368 // 369 // Initializes, if empty to MArrivalTimeCam::GetSize() for: 370 // - MHCalibrationCam::fHiGainArray, MHCalibrationCam::fLoGainArray 371 // 372 // Initializes, if empty to MGeomCam::GetNumAreas() for: 373 // - MHCalibrationCam::fAverageHiGainAreas, MHCalibrationCam::fAverageLoGainAreas 374 // 375 // Initializes, if empty to MGeomCam::GetNumSectors() for: 376 // - MHCalibrationCam::fAverageHiGainSectors, MHCalibrationCam::fAverageLoGainSectors 377 // 378 // Initializes TArray's to MGeomCam::GetNumAreas and MGeomCam::GetNumSectors, respectively 379 // Fills with number of valid pixels (if !MBadPixelsPix::IsBad()): 380 // - MHCalibrationCam::fAverageAreaNum[area index] 381 // - MHCalibrationCam::fAverageSectorNum[area index] 382 // 383 // Calls InitializeHists() for every entry in: 384 // - MHCalibrationCam::fHiGainArray 385 // - MHCalibrationCam::fAverageHiGainAreas 386 // - MHCalibrationCam::fAverageHiGainSectors 387 // 388 // Sets Titles and Names for the Histograms 389 // - MHCalibrationCam::fAverageHiGainAreas 390 // - MHCalibrationCam::fAverageHiGainSectors 391 // 392 Bool_t MHCalibrationCam::ReInit(MParList *pList) 393 { 394 395 fBadPixels = (MBadPixelsCam*)pList->FindCreateObj("MBadPixelsCam"); 396 if (!fBadPixels) 397 return kFALSE; 398 399 400 const Int_t npixels = fGeom->GetNumPixels(); 401 const Int_t nsectors = fGeom->GetNumSectors(); 402 const Int_t nareas = fGeom->GetNumAreas(); 403 404 // 405 // The function TArrayF::Set() already sets all entries to 0. 406 // 407 fAverageAreaNum. Set(nareas); 408 fAverageAreaSat. Set(nareas); 409 fAverageAreaSigma. Set(nareas); 410 fAverageAreaSigmaErr. Set(nareas); 411 fAverageAreaRelSigma. Set(nareas); 412 fAverageAreaRelSigmaErr.Set(nareas); 413 fAverageSectorNum. Set(nsectors); 414 415 for (Int_t i=0; i<npixels; i++) 416 { 417 418 if ((*fBadPixels)[i].IsBad()) 419 continue; 420 421 fAverageAreaNum [(*fGeom)[i].GetAidx() ]++; 422 fAverageSectorNum[(*fGeom)[i].GetSector()]++; 423 } 424 425 return ReInitHists(pList); 426 } 427 428 429 Bool_t MHCalibrationCam::ReInitHists(MParList *pList) 430 { 431 return kTRUE; 432 } 433 434 435 436 //-------------------------------------------------------------------------------- 437 // 438 // Retrieves from MGeomCam: 439 // - number of pixels 440 // - number of pixel areas 441 // - number of sectors 442 // 443 // For all TObjArray's (including the averaged ones), the following steps are performed: 444 // 445 // 1) Test size and return kFALSE if not matching 446 // 2) 447 // 448 Bool_t MHCalibrationCam::Fill(const MParContainer *par, const Stat_t w) 449 { 450 451 const Int_t npixels = fGeom->GetNumPixels(); 452 const Int_t nareas = fGeom->GetNumAreas(); 453 const Int_t nsectors = fGeom->GetNumSectors(); 454 455 if (fHiGainArray->GetEntries() != npixels) 456 { 457 gLog << err << "ERROR - Size mismatch... abort." << endl; 458 return kFALSE; 459 } 460 461 if (fLoGainArray->GetEntries() != npixels) 462 { 463 gLog << err << "ERROR - Size mismatch... abort." << endl; 464 return kFALSE; 465 } 466 467 if (fAverageHiGainAreas->GetEntries() != nareas) 468 { 469 *fLog << err << "ERROR - Size mismatch in number of areas ... abort." << endl; 470 return kFALSE; 471 } 472 473 if (fAverageLoGainAreas->GetEntries() != nareas) 474 { 475 *fLog << err << "ERROR - Size mismatch in number of areas ... abort." << endl; 476 return kFALSE; 477 } 478 479 if (fAverageHiGainSectors->GetEntries() != nsectors) 480 { 481 *fLog << err << "ERROR - Size mismatch in number of sectors ... abort." << endl; 482 return kFALSE; 483 } 484 485 if (fAverageLoGainSectors->GetEntries() != nsectors) 486 { 487 *fLog << err << "ERROR - Size mismatch in number of sectors ... abort." << endl; 488 return kFALSE; 489 } 490 491 return FillHists(par, w); 492 } 493 494 Bool_t MHCalibrationCam::FillHists(const MParContainer *par, const Stat_t w) 495 { 496 return kTRUE; 497 } 498 499 // -------------------------------------------------------------------------- 500 // 501 // 1) FinalizeHists() 502 // 2) FinalizeBadPixels() 503 // 3) CalcAverageSigma() 504 // 505 Bool_t MHCalibrationCam::Finalize() 506 { 507 if (!FinalizeHists()) 508 return kFALSE; 509 510 FinalizeBadPixels(); 511 CalcAverageSigma(); 512 513 return kTRUE; 514 } 515 516 Bool_t MHCalibrationCam::FinalizeHists() 517 { 518 return kTRUE; 519 } 520 521 void MHCalibrationCam::FinalizeBadPixels() 522 { 523 } 524 525 526 // ------------------------------------------------------------- 527 // 528 // If MBadPixelsPix::IsBad(): 529 // - calls MHGausEvents::SetExcluded() 530 // 531 // Calls: 532 // - MHGausEvents::InitBins() 533 // - MHGausEvents::ChangeHistId(i) 534 // - MHGausEvents::SetEventFrequency(fPulserFrequency) 535 // 536 void MHCalibrationCam::InitHists(MHGausEvents &hist, MBadPixelsPix &bad, const Int_t i) 537 { 538 539 if (bad.IsBad()) 540 hist.SetExcluded(); 541 542 hist.InitBins(); 543 hist.ChangeHistId(i); 544 hist.SetEventFrequency(fPulserFrequency); 545 546 } 547 548 void MHCalibrationCam::FitHiGainArrays(MCalibrationCam &calcam, MBadPixelsCam &badcam, 549 MBadPixelsPix::UncalibratedType_t fittyp, 550 MBadPixelsPix::UncalibratedType_t osctyp) 551 { 552 553 for (Int_t i=0; i<fHiGainArray->GetSize(); i++) 554 { 555 556 MHGausEvents &hist = (*this)[i]; 557 558 if (hist.IsExcluded()) 559 continue; 560 561 MCalibrationPix &pix = calcam[i]; 562 MBadPixelsPix &bad = badcam[i]; 563 564 FitHiGainHists(hist,pix,bad,fittyp,osctyp); 565 566 } 567 568 for (Int_t j=0; j<fAverageHiGainAreas->GetSize(); j++) 569 { 570 571 MHGausEvents &hist = GetAverageHiGainArea(j); 572 MCalibrationPix &pix = calcam.GetAverageArea(j); 573 MBadPixelsPix &bad = calcam.GetAverageBadArea(j); 574 575 FitHiGainHists(hist,pix,bad,fittyp,osctyp); 576 } 577 578 579 for (Int_t j=0; j<fAverageHiGainSectors->GetSize(); j++) 580 { 581 582 MHGausEvents &hist = GetAverageHiGainSector(j); 583 MCalibrationPix &pix = calcam.GetAverageSector(j); 584 MBadPixelsPix &bad = calcam.GetAverageBadSector(j); 585 586 FitHiGainHists(hist,pix,bad,fittyp,osctyp); 587 } 588 589 } 590 591 void MHCalibrationCam::FitLoGainArrays(MCalibrationCam &calcam, MBadPixelsCam &badcam, 592 MBadPixelsPix::UncalibratedType_t fittyp, 593 MBadPixelsPix::UncalibratedType_t osctyp) 594 { 595 596 for (Int_t i=0; i<fLoGainArray->GetSize(); i++) 597 { 598 599 MHGausEvents &hist = (*this)(i); 600 601 if (hist.IsExcluded()) 602 continue; 603 604 MCalibrationPix &pix = calcam[i]; 605 MBadPixelsPix &bad = badcam[i]; 606 607 FitLoGainHists(hist,pix,bad,fittyp,osctyp); 608 609 } 610 611 for (Int_t j=0; j<fAverageLoGainAreas->GetSize(); j++) 612 { 613 614 MHGausEvents &hist = GetAverageLoGainArea(j); 615 MCalibrationPix &pix = calcam.GetAverageArea(j); 616 MBadPixelsPix &bad = calcam.GetAverageBadArea(j); 617 618 FitLoGainHists(hist,pix,bad,fittyp,osctyp); 619 } 620 621 622 for (Int_t j=0; j<fAverageLoGainSectors->GetSize(); j++) 623 { 624 625 MHGausEvents &hist = GetAverageLoGainSector(j); 626 MCalibrationPix &pix = calcam.GetAverageSector(j); 627 MBadPixelsPix &bad = calcam.GetAverageBadSector(j); 628 629 FitLoGainHists(hist,pix,bad,fittyp,osctyp); 630 } 631 } 632 633 //------------------------------------------------------------ 634 // 635 // For all averaged areas, the fitted sigma is multiplied with the square root of 636 // the number involved pixels 637 // 638 void MHCalibrationCam::CalcAverageSigma() 639 { 640 641 for (Int_t j=0; j<fGeom->GetNumAreas(); j++) 642 { 643 644 MCalibrationPix &pix = (*fCam).GetAverageArea(j); 645 646 fAverageAreaSigma[j] = pix.GetSigma () * TMath::Sqrt((Float_t)fAverageAreaNum[j]); 647 fAverageAreaSigmaErr[j] = pix.GetSigmaErr () * TMath::Sqrt((Float_t)fAverageAreaNum[j]); 648 649 pix.SetSigma (fAverageAreaSigma[j]); 650 pix.SetSigmaErr(fAverageAreaSigmaErr[j]); 651 652 fAverageAreaRelSigma[j] = fAverageAreaSigma[j] / pix.GetMean(); 653 654 Float_t relsigmaerr = fAverageAreaSigmaErr[j]*fAverageAreaSigmaErr[j] 655 / (fAverageAreaSigma[j] *fAverageAreaSigma[j] ); 656 relsigmaerr += pix.GetMeanErr()*pix.GetMeanErr() 657 / (pix.GetMean() *pix.GetMean() ); 658 relsigmaerr *= fAverageAreaRelSigma[j]; 659 fAverageAreaRelSigmaErr[j] = TMath::Sqrt(relsigmaerr); 660 } 661 } 662 663 // --------------------------------------------------------------------------- 664 // 665 // Returns if the histogram is empty and sets the following flag: 666 // - MBadPixelsPix::SetUnsuitable(MBadPixelsPix::kUnsuitableRun) 667 // 668 // Fits the histograms with a Gaussian, in case of failure 669 // calls MHGausEvents::RepeatFit(), in case of repeated failure 670 // calls MHGausEvents::BypassFit() and sets the following flags: 671 // - MBadPixelsPix::SetUncalibrated( MBadPixelsPix::UncalibratedType_t fittyp ) 672 // - MBadPixelsPix::SetUnsuitable( MBadPixelsPix::kUnreliableRun ) 673 // 674 // Creates the fourier spectrum and tests MHGausEvents::IsFourierSpectrumOK(). 675 // In case no, sets the following flags: 676 // - MBadPixelsPix::SetUncalibrated( MBadPixelsPix::UncalibratedType_t osctyp ) 677 // - MBadPixelsPix::SetUnsuitable( MBadPixelsPix::kUnreliableRun ) 678 // 679 // Retrieves the results and store them in MCalibrationPix 680 // 681 void MHCalibrationCam::FitHiGainHists(MHGausEvents &hist, 682 MCalibrationPix &pix, 683 MBadPixelsPix &bad, 684 MBadPixelsPix::UncalibratedType_t fittyp, 685 MBadPixelsPix::UncalibratedType_t osctyp) 686 { 687 688 if (hist.IsEmpty()) 689 return; 690 691 // 692 // 2) Fit the Hi Gain histograms with a Gaussian 693 // 694 if (!hist.FitGaus()) 695 // 696 // 3) In case of failure set the bit Fitted to false and take histogram means and RMS 697 // 698 if (!hist.RepeatFit()) 699 { 700 hist.BypassFit(); 701 bad.SetUncalibrated( fittyp ); 702 } 703 704 // 705 // 4) Check for oscillations 706 // 707 hist.CreateFourierSpectrum(); 708 709 if (!hist.IsFourierSpectrumOK()) 710 bad.SetUncalibrated( osctyp ); 711 712 // 713 // 5) Retrieve the results and store them in this class 714 // 715 pix.SetHiGainMean ( hist.GetMean() ); 716 pix.SetHiGainMeanErr ( hist.GetMeanErr() ); 717 pix.SetHiGainSigma ( hist.GetSigma() ); 718 pix.SetHiGainSigmaErr ( hist.GetSigmaErr() ); 719 pix.SetHiGainProb ( hist.GetProb() ); 720 pix.SetHiGainNumPickup ( hist.GetPickup() ); 721 722 } 723 724 725 // --------------------------------------------------------------------------- 726 // 727 // Returns if the histogram is empty and sets the following flag: 728 // - MBadPixelsPix::SetUnsuitable(MBadPixelsPix::kUnsuitableRun) 729 // 730 // Fits the histograms with a Gaussian, in case of failure 731 // calls MHGausEvents::RepeatFit(), in case of repeated failure 732 // calls MHGausEvents::BypassFit() and sets the following flags: 733 // - MBadPixelsPix::SetUncalibrated( MBadPixelsPix::UncalibratedType_t fittyp ) 734 // - MBadPixelsPix::SetUnsuitable( MBadPixelsPix::kUnreliableRun ) 735 // 736 // Creates the fourier spectrum and tests MHGausEvents::IsFourierSpectrumOK(). 737 // In case no, sets the following flags: 738 // - MBadPixelsPix::SetUncalibrated( MBadPixelsPix::UncalibratedType_t osctyp ) 739 // - MBadPixelsPix::SetUnsuitable( MBadPixelsPix::kUnreliableRun ) 740 // 741 // Retrieves the results and store them in MCalibrationPix 742 // 743 void MHCalibrationCam::FitLoGainHists(MHGausEvents &hist, 744 MCalibrationPix &pix, 745 MBadPixelsPix &bad, 746 MBadPixelsPix::UncalibratedType_t fittyp, 747 MBadPixelsPix::UncalibratedType_t osctyp) 748 { 749 750 if (hist.IsEmpty()) 751 return; 752 753 754 // 755 // 2) Fit the Hi Gain histograms with a Gaussian 756 // 757 if (!hist.FitGaus()) 758 // 759 // 3) In case of failure set the bit Fitted to false and take histogram means and RMS 760 // 761 if (!hist.RepeatFit()) 762 { 763 hist.BypassFit(); 764 bad.SetUncalibrated( fittyp ); 765 } 766 767 // 768 // 4) Check for oscillations 769 // 770 hist.CreateFourierSpectrum(); 771 772 if (!hist.IsFourierSpectrumOK()) 773 bad.SetUncalibrated( osctyp ); 774 775 // 776 // 5) Retrieve the results and store them in this class 777 // 778 pix.SetLoGainMean ( hist.GetMean() ); 779 pix.SetLoGainMeanErr ( hist.GetMeanErr() ); 780 pix.SetLoGainSigma ( hist.GetSigma() ); 781 pix.SetLoGainSigmaErr ( hist.GetSigmaErr() ); 782 pix.SetLoGainProb ( hist.GetProb() ); 783 pix.SetLoGainNumPickup ( hist.GetPickup() ); 784 785 } 786 787 788 789 // -------------------------------------------------------------------------- 790 // 308 791 // Dummy, needed by MCamEvent 309 792 // … … 343 826 for (Int_t i=0; i<nareas;i++) 344 827 { 828 345 829 pad->cd(2*(i+1)-1); 346 830 GetAverageHiGainArea(i).Draw(opt); -
trunk/MagicSoft/Mars/mcalib/MHCalibrationCam.h
r3631 r3636 21 21 #endif 22 22 23 #ifndef MARS_MBadPixelsPix 24 #include "MBadPixelsPix.h" 25 #endif 26 23 27 class TText; 24 28 class TArrayI; 25 29 class TArrayF; 26 30 class MHGausEvents; 31 class MGeomCam; 32 class MCalibrationCam; 33 class MCalibrationPix; 34 class MBadPixelsCam; 35 class MBadPixelsPix; 27 36 class MHCalibrationCam : public MH, public MCamEvent 28 37 { 29 38 39 private: 40 41 static const Int_t fgAverageNbins; //! The default for fAverageNbins (now set to: 2000) 42 static const Int_t fgPulserFrequency; //! The default for fPulserFrequency (now set to: 500) 43 30 44 protected: 31 45 32 static const Int_t fgPulserFrequency; // The default for fPulserFrequency (now set to: 500)46 Int_t fAverageNbins; // Number of bins for the average histograms 33 47 Int_t fPulserFrequency; // Light pulser frequency 34 48 … … 40 54 TObjArray *fAverageLoGainSectors; //-> Array of calibration pixels, one per camera sector 41 55 56 MGeomCam *fGeom; //! Camera geometry 57 MBadPixelsCam *fBadPixels; //! Bad Pixels storage container 58 MCalibrationCam *fCam; //! Calibration Cam with the results 59 42 60 TArrayI fAverageAreaNum; // Number of pixels in average pixels per area 43 61 TArrayI fAverageAreaSat; // Number of saturated slices in average pixels per area … … 48 66 TArrayI fAverageSectorNum; // Number of pixels in average pixels per sector 49 67 68 virtual Bool_t SetupHists(const MParList *pList); 69 virtual Bool_t ReInitHists(MParList *pList); 70 virtual Bool_t FillHists(const MParContainer *par, const Stat_t w=1); 71 virtual Bool_t FinalizeHists(); 72 virtual void FinalizeBadPixels(); 73 74 void InitHists(MHGausEvents &hist, MBadPixelsPix &bad, const Int_t i); 75 76 void FitHiGainArrays(MCalibrationCam &calcam, MBadPixelsCam &badcam, 77 MBadPixelsPix::UncalibratedType_t fittyp, 78 MBadPixelsPix::UncalibratedType_t osctyp); 79 80 void FitLoGainArrays(MCalibrationCam &calcam, MBadPixelsCam &badcam, 81 MBadPixelsPix::UncalibratedType_t fittyp, 82 MBadPixelsPix::UncalibratedType_t osctyp); 83 84 void FitHiGainHists(MHGausEvents &hist, 85 MCalibrationPix &pix, 86 MBadPixelsPix &bad, 87 MBadPixelsPix::UncalibratedType_t fittyp, 88 MBadPixelsPix::UncalibratedType_t osctyp); 89 90 void FitLoGainHists(MHGausEvents &hist, 91 MCalibrationPix &pix, 92 MBadPixelsPix &bad, 93 MBadPixelsPix::UncalibratedType_t fittyp, 94 MBadPixelsPix::UncalibratedType_t osctyp); 95 96 void CalcAverageSigma(); 97 50 98 void DrawAverageSigma(Bool_t sat, Bool_t inner, 51 99 Float_t sigma, Float_t sigmaerr, … … 57 105 ~MHCalibrationCam(); 58 106 107 void SetAverageNbins( const Int_t bins=fgAverageNbins ) { fAverageNbins = bins; } 59 108 void SetPulserFrequency(const Int_t f=fgPulserFrequency) { fPulserFrequency = f; } 60 109 … … 76 125 MHGausEvents &GetAverageLoGainSector(UInt_t i); 77 126 const MHGausEvents &GetAverageLoGainSector(UInt_t i) const; 127 128 virtual Bool_t SetupFill(const MParList *pList); 129 virtual Bool_t ReInit ( MParList *pList); 130 virtual Bool_t Fill (const MParContainer *par, const Stat_t w=1); 131 virtual Bool_t Finalize ( ); 78 132 79 133 // Clone -
trunk/MagicSoft/Mars/mcalib/MHCalibrationChargeCam.cc
r3631 r3636 127 127 #include "MGeomPix.h" 128 128 129 #include "MHGausEvents.h" 130 129 131 #include "MBadPixelsCam.h" 130 132 #include "MBadPixelsPix.h" … … 140 142 using namespace std; 141 143 142 const Int_t MHCalibrationChargeCam::fgAverageNbins = 4000;143 144 const Float_t MHCalibrationChargeCam::fgNumHiGainSaturationLimit = 0.01; 144 145 const Float_t MHCalibrationChargeCam::fgNumLoGainSaturationLimit = 0.005; … … 150 151 // - all pointers to NULL 151 152 // 152 // Initializes and sets owner of (done in MHCalibrationCam):153 // - fHiGainArray, fLoGainArray154 // - fAverageHiGainAreas, fAverageLoGainAreas155 // - fAverageHiGainSectors, fAverageLoGainSectors156 //157 153 // Initializes: 158 154 // - fNumHiGainSaturationLimit to fgNumHiGainSaturationLimit 159 155 // - fNumLoGainSaturationLimit to fgNumLoGainSaturationLimit 160 // - fPulserFrequency to fgPulserFrequency (in MHCalibrationCam)161 156 // 162 157 MHCalibrationChargeCam::MHCalibrationChargeCam(const char *name, const char *title) 163 : f Cam(NULL), fRawEvt(NULL), fGeom(NULL), fBadPixels(NULL)158 : fRawEvt(NULL) 164 159 { 165 160 fName = name ? name : "MHCalibrationChargeCam"; … … 174 169 // Gets the pointers to: 175 170 // - MRawEvtData 176 // - MGeomCam 177 // 178 // Calls Delete-Function of: 179 // - MHCalibrationCam::fHiGainArray, MHCalibrationCam::fLoGainArray 180 // - MHCalibrationCam::fAverageHiGainAreas, MHCalibrationCam::fAverageLoGainAreas 181 // - MHCalibrationCam::fAverageHiGainSectors, MHCalibrationCam::fAverageLoGainSectors 182 // 183 Bool_t MHCalibrationChargeCam::SetupFill(const MParList *pList) 171 // 172 Bool_t MHCalibrationChargeCam::SetupHists(const MParList *pList) 184 173 { 185 174 … … 191 180 } 192 181 193 fGeom = (MGeomCam*)pList->FindObject("MGeomCam");194 if (!fGeom)195 {196 *fLog << err << "MGeomCam not found... aborting." << endl;197 return kFALSE;198 }199 200 fHiGainArray->Delete();201 fLoGainArray->Delete();202 203 fAverageHiGainAreas->Delete();204 fAverageLoGainAreas->Delete();205 206 fAverageHiGainSectors->Delete();207 fAverageLoGainSectors->Delete();208 209 182 return kTRUE; 210 183 } … … 213 186 // 214 187 // Gets or creates the pointers to: 215 // - MBadPixelsCam216 188 // - MCalibrationChargeCam 217 189 // … … 219 191 // - MExtractedSignalCam 220 192 // 221 // Initializes, if empty to MExtractedSignalCam::GetSize() for: 222 // - MHCalibrationCam::fHiGainArray, MHCalibrationCam::fLoGainArray 223 // 224 // Initializes, if empty to MGeomCam::GetNumAreas() for: 225 // - MHCalibrationCam::fAverageHiGainAreas, MHCalibrationCam::fAverageLoGainAreas 226 // 227 // Initializes, if empty to MGeomCam::GetNumSectors() for: 228 // - MHCalibrationCam::fAverageHiGainSectors, MHCalibrationCam::fAverageLoGainSectors 229 // 230 // Initializes TArrays to MGeomCam::GetNumAreas and MGeomCam::GetNumSectors, respectively 231 // Fills with number of valid pixels (if !MBadPixelsPix::IsBad()): 232 // - MHCalibrationCam::fAverageAreaNum[area index] 233 // - MHCalibrationCam::fAverageSectorNum[area index] 234 // 235 // Calls InitializeHiGainHists() and InitializeLoGainHists() for every entry in: 193 // Calls InitializeHists() for every entry in: 236 194 // - MHCalibrationCam::fHiGainArray, MHCalibrationCam::fLoGainArray 237 195 // - MHCalibrationCam::fAverageHiGainAreas, MHCalibrationCam::fAverageLoGainAreas … … 239 197 // 240 198 // Sets Titles and Names for the Charge Histograms and 241 // Sets number of bins to fAverageNbins for:199 // Sets number of bins to MHCalibrationCam::fAverageNbins for: 242 200 // - MHCalibrationCam::fAverageHiGainAreas, MHCalibrationCam::fAverageLoGainAreas 243 201 // - MHCalibrationCam::fAverageHiGainSectors, MHCalibrationCam::fAverageLoGainSectors 244 202 // 245 Bool_t MHCalibrationChargeCam::ReInit (MParList *pList)203 Bool_t MHCalibrationChargeCam::ReInitHists(MParList *pList) 246 204 { 247 205 248 fBadPixels = (MBadPixelsCam*)pList->FindCreateObj("MBadPixelsCam"); 249 if (!fBadPixels) 250 return kFALSE; 251 252 253 fCam = (MCalibrationChargeCam*)pList->FindCreateObj("MCalibrationChargeCam"); 206 fCam = (MCalibrationCam*)pList->FindCreateObj("MCalibrationChargeCam"); 254 207 if (!fCam) 255 208 return kFALSE; … … 262 215 } 263 216 264 const Int_t n = signal->GetSize();217 const Int_t npixels = fGeom->GetNumPixels(); 265 218 const Int_t nsectors = fGeom->GetNumSectors(); 266 219 const Int_t nareas = fGeom->GetNumAreas(); 267 220 268 //269 // The function TArrayF::Set() already sets all entries to 0.270 //271 fAverageAreaNum. Set(nareas);272 fAverageAreaSat. Set(nareas);273 fAverageAreaSigma. Set(nareas);274 fAverageAreaSigmaErr. Set(nareas);275 fAverageAreaRelSigma. Set(nareas);276 fAverageAreaRelSigmaErr.Set(nareas);277 fAverageSectorNum. Set(nsectors);278 279 for (Int_t i=0; i<n; i++)280 {281 if ((*fBadPixels)[i].IsBad())282 continue;283 fAverageAreaNum [(*fGeom)[i].GetAidx() ]++;284 fAverageSectorNum[(*fGeom)[i].GetSector()]++;285 }286 287 221 if (fHiGainArray->GetEntries()==0) 288 222 { 289 fHiGainArray->Expand(n );290 for (Int_t i=0; i<n ; i++)223 fHiGainArray->Expand(npixels); 224 for (Int_t i=0; i<npixels; i++) 291 225 { 292 226 (*fHiGainArray)[i] = new MHCalibrationChargeHiGainPix; 293 Init ializeHiGainHists((MHCalibrationChargePix&)(*this)[i],(*fBadPixels)[i],i);227 InitHists((*this)[i],(*fBadPixels)[i],i); 294 228 } 295 229 } 296 230 297 298 231 if (fLoGainArray->GetEntries()==0) 299 232 { 300 fLoGainArray->Expand(n );301 302 for (Int_t i=0; i<n ; i++)233 fLoGainArray->Expand(npixels); 234 235 for (Int_t i=0; i<npixels; i++) 303 236 { 304 237 (*fLoGainArray)[i] = new MHCalibrationChargeLoGainPix; 305 Init ializeLoGainHists((MHCalibrationChargePix&)(*this)(i),(*fBadPixels)[i],i);238 InitHists((*this)(i),(*fBadPixels)[i],i); 306 239 } 307 240 … … 318 251 "Average HiGain FADC sums area idx "); 319 252 253 InitHists(GetAverageHiGainArea(j),fCam->GetAverageBadArea(j),j); 254 320 255 MHCalibrationChargePix &hist = (MHCalibrationChargePix&)GetAverageHiGainArea(j); 321 256 322 InitializeHiGainHists(hist,fCam->GetAverageBadArea(j),j);323 324 257 hist.GetHGausHist()->SetTitle("Summed FADC slices average HiGain Area Idx "); 325 hist.Set ChargeNbins(fAverageNbins);258 hist.SetNbins(fAverageNbins); 326 259 hist.GetHAbsTime()->SetTitle("Absolute Arrival Time average HiGain Area Idx "); 327 260 } … … 340 273 MHCalibrationChargePix &hist = (MHCalibrationChargePix&)GetAverageLoGainArea(j); 341 274 342 Init ializeLoGainHists(hist,fCam->GetAverageBadArea(j),j);275 InitHists(hist,fCam->GetAverageBadArea(j),j); 343 276 344 277 hist.GetHGausHist()->SetTitle("Summed FADC slices average LoGain Area Idx "); 345 hist.Set ChargeNbins(fAverageNbins);278 hist.SetNbins(fAverageNbins); 346 279 hist.GetHAbsTime()->SetTitle("Absolute Arrival Time average LoGain Area Idx "); 347 280 … … 361 294 MHCalibrationChargePix &hist = (MHCalibrationChargePix&)GetAverageHiGainSector(j); 362 295 363 Init ializeHiGainHists(hist,fCam->GetAverageBadSector(j),j);296 InitHists(hist,fCam->GetAverageBadSector(j),j); 364 297 365 298 hist.GetHGausHist()->SetTitle("Summed FADC slices average HiGain Sector "); 366 hist.Set ChargeNbins(fAverageNbins);299 hist.SetNbins(fAverageNbins); 367 300 hist.GetHAbsTime()->SetTitle("Absolute Arrival Time average HiGain Sector "); 368 301 … … 382 315 MHCalibrationChargePix &hist = (MHCalibrationChargePix&)GetAverageLoGainSector(j); 383 316 384 Init ializeLoGainHists(hist,fCam->GetAverageBadSector(j),j);317 InitHists(hist,fCam->GetAverageBadSector(j),j); 385 318 386 319 hist.GetHGausHist()->SetTitle("Summed FADC slices average LoGain Sector "); 387 hist.Set ChargeNbins(fAverageNbins);320 hist.SetNbins(fAverageNbins); 388 321 hist.GetHAbsTime()->SetTitle("Absolute Arrival Time average LoGain Sector "); 389 322 … … 394 327 } 395 328 396 // -------------------------------------------------------------397 //398 // If MBadPixelsPix::IsBad():399 // - calls MHCalibrationChargePix::SetExcluded()400 //401 // Calls:402 // - MHCalibrationChargePix::Init()403 // - MHCalibrationChargePix::ChangeHistId(i)404 // - MHGausEvents::SetEventFrequency(fPulserFrequency)405 //406 void MHCalibrationChargeCam::InitializeHiGainHists(MHCalibrationChargePix &hist, MBadPixelsPix &bad, const Int_t i)407 {408 409 if (bad.IsBad())410 {411 *fLog << warn << "Excluded pixel: " << i << " from calibration " << endl;412 hist.SetExcluded();413 }414 415 hist.Init();416 hist.ChangeHistId(i);417 hist.SetEventFrequency(fPulserFrequency);418 419 }420 421 // -------------------------------------------------------------422 //423 // If MBadPixelsPix::IsBad():424 // - calls MHCalibrationChargePix::SetExcluded()425 //426 // Calls:427 // - MHCalibrationChargePix::Init()428 // - MHCalibrationChargePix::ChangeHistId(i)429 // - MHGausEvents::SetEventFrequency(fPulserFrequency)430 //431 void MHCalibrationChargeCam::InitializeLoGainHists(MHCalibrationChargePix &hist,MBadPixelsPix &bad, const Int_t i)432 {433 434 if (bad.IsBad())435 hist.SetExcluded();436 437 hist.Init();438 hist.ChangeHistId(i);439 hist.SetEventFrequency(fPulserFrequency);440 }441 442 329 443 330 // -------------------------------------------------------------------------- 444 331 // 445 332 // Retrieves from MExtractedSignalCam: 333 // - first used LoGain FADC slice 334 // 335 // Retrieves from MGeomCam: 446 336 // - number of pixels 447 // - first used LoGain FADC slice448 //449 // Retrieves from MGeomCam:450 337 // - number of pixel areas 451 338 // - number of sectors … … 453 340 // For all TObjArray's (including the averaged ones), the following steps are performed: 454 341 // 455 // 1) Test size and return kFALSE if not matching 456 // 457 // 2) Fill Charges histograms (MHGausEvents::FillHistAndArray()) with: 342 // 1) Fill Charges histograms (MHGausEvents::FillHistAndArray()) with: 458 343 // - MExtractedSignalPix::GetExtractedSignalHiGain(); 459 344 // - MExtractedSignalPix::GetExtractedSignalLoGain(); 460 345 // 461 // 3) Set number of saturated slices (MHCalibrationChargePix::SetSaturated()) with:346 // 2) Set number of saturated slices (MHCalibrationChargePix::SetSaturated()) with: 462 347 // - MExtractedSignalPix::GetNumHiGainSaturated(); 463 348 // - MExtractedSignalPix::GetNumLoGainSaturated(); 464 349 // 465 // 4) Fill AbsTime histograms (MHCalibrationChargePix::FillAbsTime()) with:350 // 3) Fill AbsTime histograms (MHCalibrationChargePix::FillAbsTime()) with: 466 351 // - MRawEvtPixelIter::GetIdxMaxHiGainSample(); 467 352 // - MRawEvtPixelIter::GetIdxMaxLoGainSample(first slice); 468 353 // 469 Bool_t MHCalibrationChargeCam::Fill (const MParContainer *par, const Stat_t w)354 Bool_t MHCalibrationChargeCam::FillHists(const MParContainer *par, const Stat_t w) 470 355 { 471 356 … … 477 362 } 478 363 479 const Int_t npixels = signal->GetSize(); 480 const Int_t lofirst = signal->GetFirstUsedSliceLoGain(); 364 const Int_t npixels = fGeom->GetNumPixels(); 481 365 const Int_t nareas = fGeom->GetNumAreas(); 482 366 const Int_t nsectors = fGeom->GetNumSectors(); 483 484 if (fHiGainArray->GetEntries() != npixels) 485 { 486 *fLog << err << "ERROR - Size mismatch in number of pixels ... abort." << endl; 487 return kFALSE; 488 } 489 490 if (fLoGainArray->GetEntries() != npixels) 491 { 492 *fLog << err << "ERROR - Size mismatch in number of pixels ... abort." << endl; 493 return kFALSE; 494 } 495 496 if (fAverageHiGainAreas->GetEntries() != nareas) 497 { 498 *fLog << err << "ERROR - Size mismatch in number of areas ... abort." << endl; 499 return kFALSE; 500 } 501 502 if (fAverageLoGainAreas->GetEntries() != nareas) 503 { 504 *fLog << err << "ERROR - Size mismatch in number of areas ... abort." << endl; 505 return kFALSE; 506 } 507 508 if (fAverageHiGainSectors->GetEntries() != nsectors) 509 { 510 *fLog << err << "ERROR - Size mismatch in number of sectors ... abort." << endl; 511 return kFALSE; 512 } 513 514 if (fAverageLoGainSectors->GetEntries() != nsectors) 515 { 516 *fLog << err << "ERROR - Size mismatch in number of sectors ... abort." << endl; 517 return kFALSE; 518 } 367 const Int_t lofirst = signal->GetFirstUsedSliceLoGain(); 519 368 520 369 Float_t sumhiarea [nareas], sumloarea [nareas], timehiarea [nareas], timeloarea [nareas]; … … 526 375 { 527 376 sumhiarea [j] = sumloarea [j] = timehiarea [j] = timeloarea [j] = 0.; 377 sathiarea [j] = satloarea [j] = 0; 378 } 379 380 for (UInt_t j=0; j<nsectors; j++) 381 { 528 382 sumhisector[j] = sumlosector[j] = timehisector[j] = timelosector[j] = 0.; 529 sathiarea [j] = satloarea [j] = 0;530 383 sathisector[j] = satlosector[j] = 0; 531 384 } 532 385 386 533 387 for (Int_t i=0; i<npixels; i++) 534 388 { … … 596 450 } 597 451 598 599 452 for (UInt_t j=0; j<nareas; j++) 600 453 { … … 642 495 // For all TObjArray's (including the averaged ones), the following steps are performed: 643 496 // 644 // 1) Return if the pixel is excluded 645 // 646 // 2) Call to FinalizeHiGainHists() and FinalizeLoGainHists() 647 // 648 // For all averaged areas, the fitted sigma is multiplied with the square root of 649 // the number involved pixels 650 // 651 Bool_t MHCalibrationChargeCam::Finalize() 497 // 1) Returns if the pixel is excluded. 498 // 2) Tests saturation. In case yes, set the flag: MCalibrationPix::SetHiGainSaturation() 499 // or the flag: MBadPixelsPix::SetUncalibrated( MBadPixelsPix::kLoGainSaturated ) 500 // 3) Store the absolute arrival times in the MCalibrationChargePix's. If flag 501 // MCalibrationPix::IsHiGainSaturation() is set, the Low-Gain arrival times are stored, 502 // otherwise the Hi-Gain ones. 503 // 4) Calls to MHCalibrationCam::FitHiGainArrays() and MCalibrationCam::FitLoGainArrays() 504 // and sets the flags (if occurring): 505 // - MBadPixelsPix::SetUncalibrated( MBadPixelsPix::kHiGainNotFitted ) 506 // - MBadPixelsPix::SetUncalibrated( MBadPixelsPix::kLoGainNotFitted ) 507 // - MBadPixelsPix::SetUncalibrated( MBadPixelsPix::kHiGainOscillating ) 508 // - MBadPixelsPix::SetUncalibrated( MBadPixelsPix::kLoGainOscillating ) 509 // 5) Retrieves the results and store them in MCalibrationChargePix 510 // 511 Bool_t MHCalibrationChargeCam::FinalizeHists() 652 512 { 653 513 … … 656 516 657 517 MHCalibrationChargePix &histhi = (MHCalibrationChargePix&)(*this)[i]; 658 518 MCalibrationChargePix &pix = (MCalibrationChargePix&)(*fCam)[i]; 519 659 520 if (histhi.IsExcluded()) 660 521 continue; 661 522 662 MCalibrationChargePix &pix = (*fCam)[i]; 523 if (histhi.GetSaturated() > fNumHiGainSaturationLimit*histhi.GetHGausHist()->GetEntries()) 524 { 525 pix.SetHiGainSaturation(); 526 histhi.CreateFourierSpectrum(); 527 continue; 528 } 529 530 pix.SetAbsTimeMean ( histhi.GetAbsTimeMean()); 531 pix.SetAbsTimeRms ( histhi.GetAbsTimeRms() ); 532 } 533 534 for (Int_t i=0; i<fLoGainArray->GetSize(); i++) 535 { 536 537 MHCalibrationChargePix &histlo = (MHCalibrationChargePix&)(*this)(i); 663 538 MBadPixelsPix &bad = (*fBadPixels)[i]; 664 665 FinalizeHiGainHists(histhi,pix,bad);666 667 }668 669 for (Int_t i=0; i<fLoGainArray->GetSize(); i++)670 {671 672 MHCalibrationChargePix &histlo = (MHCalibrationChargePix&)(*this)(i);673 539 674 540 if (histlo.IsExcluded()) 675 541 continue; 676 542 677 MCalibrationChargePix &pix = (*fCam)[i]; 678 MBadPixelsPix &bad = (*fBadPixels)[i]; 679 680 FinalizeLoGainHists(histlo,pix,bad); 681 682 } 683 543 if (histlo.GetSaturated() > fNumLoGainSaturationLimit*histlo.GetHGausHist()->GetEntries()) 544 { 545 *fLog << warn << "Saturated Lo Gain histogram in pixel: " << i << endl; 546 bad.SetUncalibrated( MBadPixelsPix::kLoGainSaturation ); 547 histlo.CreateFourierSpectrum(); 548 continue; 549 } 550 551 MCalibrationChargePix &pix = (MCalibrationChargePix&)(*fCam)[i]; 552 553 if (pix.IsHiGainSaturation()) 554 { 555 pix.SetAbsTimeMean ( histlo.GetAbsTimeMean()); 556 pix.SetAbsTimeRms ( histlo.GetAbsTimeRms() ); 557 } 558 } 559 684 560 for (Int_t j=0; j<fAverageHiGainAreas->GetSize(); j++) 685 561 { 686 562 687 563 MHCalibrationChargePix &histhi = (MHCalibrationChargePix&)GetAverageHiGainArea(j); 688 MCalibrationChargePix &pix = fCam->GetAverageArea(j); 689 MBadPixelsPix &bad = fCam->GetAverageBadArea(j); 690 691 FinalizeHiGainHists(histhi,pix,bad); 564 MCalibrationChargePix &pix = (MCalibrationChargePix&)fCam->GetAverageArea(j); 565 566 if (histhi.GetSaturated() > fNumHiGainSaturationLimit*histhi.GetHGausHist()->GetEntries()) 567 { 568 pix.SetHiGainSaturation(); 569 histhi.CreateFourierSpectrum(); 570 continue; 571 } 572 573 pix.SetAbsTimeMean ( histhi.GetAbsTimeMean()); 574 pix.SetAbsTimeRms ( histhi.GetAbsTimeRms() ); 692 575 } 693 576 … … 696 579 697 580 MHCalibrationChargePix &histlo = (MHCalibrationChargePix&)GetAverageLoGainArea(j); 698 MCalibrationChargePix &pix = fCam->GetAverageArea(j); 699 MBadPixelsPix &bad = fCam->GetAverageBadArea(j); 700 701 FinalizeLoGainHists(histlo,pix,bad); 581 MCalibrationChargePix &pix = (MCalibrationChargePix&)fCam->GetAverageArea(j); 582 583 if (histlo.GetSaturated() > fNumLoGainSaturationLimit*histlo.GetHGausHist()->GetEntries()) 584 { 585 *fLog << warn << "Saturated Lo Gain histogram in area idx: " << j << endl; 586 histlo.CreateFourierSpectrum(); 587 continue; 588 } 589 590 if (pix.IsHiGainSaturation()) 591 { 592 pix.SetAbsTimeMean ( histlo.GetAbsTimeMean()); 593 pix.SetAbsTimeRms ( histlo.GetAbsTimeRms() ); 594 } 702 595 } 703 596 … … 706 599 707 600 MHCalibrationChargePix &histhi = (MHCalibrationChargePix&)GetAverageHiGainSector(j); 708 MCalibrationChargePix &pix = fCam->GetAverageSector(j); 601 MCalibrationChargePix &pix = (MCalibrationChargePix&)fCam->GetAverageSector(j); 602 603 if (histhi.GetSaturated() > fNumHiGainSaturationLimit*histhi.GetHGausHist()->GetEntries()) 604 { 605 pix.SetHiGainSaturation(); 606 histhi.CreateFourierSpectrum(); 607 continue; 608 } 609 610 pix.SetAbsTimeMean ( histhi.GetAbsTimeMean()); 611 pix.SetAbsTimeRms ( histhi.GetAbsTimeRms() ); 612 } 613 614 for (Int_t j=0; j<fAverageLoGainSectors->GetSize(); j++) 615 { 616 617 MHCalibrationChargePix &histlo = (MHCalibrationChargePix&)GetAverageLoGainSector(j); 618 MCalibrationChargePix &pix = (MCalibrationChargePix&)fCam->GetAverageSector(j); 709 619 MBadPixelsPix &bad = fCam->GetAverageBadSector(j); 710 620 711 FinalizeHiGainHists(histhi,pix,bad); 712 } 713 714 for (Int_t j=0; j<fAverageLoGainSectors->GetSize(); j++) 715 { 716 717 MHCalibrationChargePix &histlo = (MHCalibrationChargePix&)GetAverageLoGainSector(j); 718 MCalibrationChargePix &pix = fCam->GetAverageSector(j); 719 MBadPixelsPix &bad = fCam->GetAverageBadSector(j); 720 721 FinalizeLoGainHists(histlo,pix,bad); 722 } 723 724 for (Int_t j=0; j<fGeom->GetNumAreas();j++) 725 { 726 727 MCalibrationChargePix &pix = fCam->GetAverageArea(j); 621 if (histlo.GetSaturated() > fNumLoGainSaturationLimit*histlo.GetHGausHist()->GetEntries()) 622 { 623 *fLog << warn << "Saturated Lo Gain histogram in sector: " << j << endl; 624 bad.SetUncalibrated( MBadPixelsPix::kLoGainSaturation ); 625 histlo.CreateFourierSpectrum(); 626 continue; 627 } 728 628 729 629 if (pix.IsHiGainSaturation()) 730 fAverageAreaSat[j]++;731 732 fAverageAreaSigma[j] = pix.GetSigmaCharge () * TMath::Sqrt((Float_t)fAverageAreaNum[j]);733 fAverageAreaSigmaErr[j] = pix.GetSigmaChargeErr () * TMath::Sqrt((Float_t)fAverageAreaNum[j]);734 735 pix.SetSigmaCharge (fAverageAreaSigma[j]);736 pix.SetSigmaChargeErr(fAverageAreaSigmaErr[j]);737 738 fAverageAreaRelSigma[j] = fAverageAreaSigma[j] / pix.GetMeanCharge();739 740 Float_t relsigmaerr = fAverageAreaSigmaErr[j]*fAverageAreaSigmaErr[j]741 / (fAverageAreaSigma[j] *fAverageAreaSigma[j]);742 relsigmaerr += pix.GetMeanChargeErr()*pix.GetMeanChargeErr()743 / (pix.GetMeanCharge() *pix.GetMeanCharge() );744 relsigmaerr *= fAverageAreaRelSigma[j];745 fAverageAreaRelSigmaErr[j] = TMath::Sqrt(relsigmaerr);746 747 }748 630 { 631 pix.SetAbsTimeMean ( histlo.GetAbsTimeMean()); 632 pix.SetAbsTimeRms ( histlo.GetAbsTimeRms() ); 633 } 634 } 635 636 // 637 // Perform the fitting for the High Gain (done in MHCalibrationCam) 638 // 639 FitHiGainArrays((MCalibrationCam&)(*fCam),(*fBadPixels), 640 MBadPixelsPix::kHiGainNotFitted, 641 MBadPixelsPix::kHiGainOscillating); 642 // 643 // Perform the fitting for the Low Gain (done in MHCalibrationCam) 644 // 645 FitLoGainArrays((MCalibrationCam&)(*fCam),(*fBadPixels), 646 MBadPixelsPix::kLoGainNotFitted, 647 MBadPixelsPix::kLoGainOscillating); 648 749 649 return kTRUE; 750 650 } 751 651 752 // --------------------------------------------------------------------------- 753 // 754 // Returns if the histogram is empty and sets the following flag: 755 // - MBadPixelsPix::SetUnsuitable(MBadPixelsPix::kUnsuitableRun) 756 // 757 // Returns if the number of saturated slices (MHCalibrationChargePix::GetSaturated()) is 758 // higher than the allowed limit (fNumHiGainSaturationLimit*histogram-entries), creates 759 // the fourier spectrum and sets the following flags: 760 // - MCalibrationChargePix::SetHiGainSaturation() 761 // 762 // Fits the histograms with a Gaussian, in case of failure 763 // calls MHCalibrationChargePix::RepeatFit(), in case of repeated failure 764 // calls MHGausEvents::BypassFit() and sets the following flags: 765 // - MBadPixelsPix::SetUncalibrated( MBadPixelsPix::kHiGainNotFitted ) 766 // - MBadPixelsPix::SetUnsuitable( MBadPixelsPix::kUnreliableRun ) 767 // 768 // Counts the number of pickup events 769 // 770 // Creates the fourier spectrum and tests MHGausEvents::IsFourierSpectrumOK(). 771 // In case no, sets the following flags: 772 // - MBadPixelsPix::SetUncalibrated( MBadPixelsPix::kHiGainOscillating ) 773 // - MBadPixelsPix::SetUnsuitable( MBadPixelsPix::kUnreliableRun ) 774 // 775 // Retrieves the results and store them in MCalibrationChargePix 776 // 777 void MHCalibrationChargeCam::FinalizeHiGainHists(MHCalibrationChargePix &hist, 778 MCalibrationChargePix &pix, 779 MBadPixelsPix &bad) 652 // -------------------------------------------------------------------------- 653 // 654 // Takes the decisions under which conditions a pixel is declared: 655 // MBadPixelsPix::kUnsuitableRun or MBadPixelsPix::kUnreliableRun, namely: 656 // * if MBadPixelsPix::kHiGainNotFitted and !MCalibrationPix::IsHiGainSaturation() 657 // sets MBadPixelsPix::kUnreliableRun 658 // * if MBadPixelsPix::kHiGainOscillating and !MCalibrationPix::IsHiGainSaturation() 659 // sets MBadPixelsPix::kUnreliableRun 660 // * if MBadPixelsPix::kLoGainNotFitted and MCalibrationPix::IsLoGainSaturation() 661 // sets MBadPixelsPix::kUnreliableRun 662 // * if MBadPixelsPix::kLoGainOscillating and MCalibrationPix::IsLoGainSaturation() 663 // sets MBadPixelsPix::kUnreliableRun 664 // * if MBadPixelsPix::kLoGainSaturation 665 // sets MBadPixelsPix::kUnsuitableRun 666 // 667 void MHCalibrationChargeCam::FinalizeBadPixels() 780 668 { 781 782 if (hist.IsEmpty()) 783 { 784 *fLog << warn << "Empty Hi Gain histogram in pixel: " << pix.GetPixId() << endl; 785 bad.SetUnsuitable(MBadPixelsPix::kUnsuitableRun); 786 return; 787 } 788 789 if (hist.GetSaturated() > fNumHiGainSaturationLimit*hist.GetHGausHist()->GetEntries()) 790 { 791 *fLog << warn << "Saturated Hi Gain histogram in pixel: " << pix.GetPixId() << endl; 792 pix.SetHiGainSaturation(); 793 hist.CreateFourierSpectrum(); 794 return; 795 } 796 797 // 798 // 2) Fit the Hi Gain histograms with a Gaussian 799 // 800 if (!hist.FitGaus()) 801 // 802 // 3) In case of failure set the bit Fitted to false and take histogram means and RMS 803 // 804 if (!hist.RepeatFit()) 805 { 806 hist.BypassFit(); 807 *fLog << warn << "Hi Gain could not be fitted in pixel: " << pix.GetPixId() << endl; 808 bad.SetUncalibrated( MBadPixelsPix::kHiGainNotFitted ); 809 bad.SetUnsuitable( MBadPixelsPix::kUnreliableRun ); 810 } 811 812 // 813 // 4) Check for pickup 814 // 815 hist.CountPickup(); 816 817 // 818 // 5) Check for oscillations 819 // 820 hist.CreateFourierSpectrum(); 821 822 if (!hist.IsFourierSpectrumOK()) 823 { 824 *fLog << warn << "Oscillating Hi Gain in pixel: " << pix.GetPixId() << endl; 825 bad.SetUncalibrated( MBadPixelsPix::kHiGainOscillating ); 826 bad.SetUnsuitable( MBadPixelsPix::kUnreliableRun ); 827 } 828 829 // 830 // 6) Retrieve the results and store them in this class 831 // 832 pix.SetHiGainMeanCharge( hist.GetMean() ); 833 pix.SetHiGainMeanChargeErr( hist.GetMeanErr() ); 834 pix.SetHiGainSigmaCharge( hist.GetSigma() ); 835 pix.SetHiGainSigmaChargeErr( hist.GetSigmaErr() ); 836 pix.SetHiGainChargeProb ( hist.GetProb() ); 837 838 pix.SetAbsTimeMean ( hist.GetAbsTimeMean()); 839 pix.SetAbsTimeRms ( hist.GetAbsTimeRms() ); 840 841 pix.SetHiGainNumPickup ( hist.GetPickup() ); 842 669 670 for (Int_t i=0; i<fBadPixels->GetSize(); i++) 671 { 672 673 MBadPixelsPix &bad = (*fBadPixels)[i]; 674 MCalibrationPix &pix = (*fCam)[i]; 675 676 if (bad.IsUncalibrated( MBadPixelsPix::kHiGainNotFitted )) 677 if (!pix.IsHiGainSaturation()) 678 bad.SetUnsuitable( MBadPixelsPix::kUnreliableRun ); 679 680 if (bad.IsUncalibrated( MBadPixelsPix::kHiGainOscillating )) 681 bad.SetUnsuitable( MBadPixelsPix::kUnreliableRun ); 682 683 if (bad.IsUncalibrated( MBadPixelsPix::kLoGainNotFitted )) 684 if (pix.IsHiGainSaturation()) 685 bad.SetUnsuitable( MBadPixelsPix::kUnreliableRun ); 686 687 if (bad.IsUncalibrated( MBadPixelsPix::kLoGainOscillating )) 688 if (pix.IsHiGainSaturation()) 689 bad.SetUnsuitable( MBadPixelsPix::kUnreliableRun ); 690 691 if (bad.IsUncalibrated( MBadPixelsPix::kLoGainSaturation )) 692 bad.SetUnsuitable( MBadPixelsPix::kUnsuitableRun ); 693 } 843 694 } 844 845 // ---------------------------------------------------------------------------846 //847 // Returns if the histogram is empty. If yes and the flag MCalibrationChargePix::IsHiGainSaturation()848 // is set, sets the following flag:849 // - MBadPixelsPix::SetUnsuitable(MBadPixelsPix::kUnsuitableRun)850 //851 // Returns if the number of saturated slices (MHCalibrationChargePix::GetSaturated()) is852 // higher than the allowed limit (fNumLoGainSaturationLimit*histogram-entries), creates853 // the fourier spectrum and sets the following flags:854 // - MBadPixelsPix::SetUncalibrated( MBadPixelsPix::kLoGainSaturation );855 // - MBadPixelsPix::SetUnsuitable( MBadPixelsPix::kUnsuitableRun );856 //857 // Fits the histograms with a Gaussian, in case of failure calls MHCalibrationChargePix::RepeatFit(),858 // in case of repeated failure calls MHGausEvents::BypassFit() and sets the following flags:859 // - MBadPixelsPix::SetUncalibrated( MBadPixelsPix::kLoGainNotFitted )860 // If the flag MCalibrationChargePix::IsHiGainSaturation() is set, sets additionally:861 // - MBadPixelsPix::SetUnsuitable( MBadPixelsPix::kUnreliableRun )862 //863 // Counts the number of pickup events864 //865 // Creates the fourier spectrum and tests MHGausEvents::IsFourierSpectrumOK().866 // In case no, sets the following flags:867 // - MBadPixelsPix::SetUncalibrated( MBadPixelsPix::kLoGainOscillating )868 // If the flag MCalibrationChargePix::IsHiGainSaturation() is set, sets additionally:869 // - MBadPixelsPix::SetUnsuitable( MBadPixelsPix::kUnreliableRun )870 //871 // Retrieves the results and store them in MCalibrationChargePix872 //873 void MHCalibrationChargeCam::FinalizeLoGainHists(MHCalibrationChargePix &hist,874 MCalibrationChargePix &pix,875 MBadPixelsPix &bad)876 {877 878 if (hist.IsEmpty())879 {880 if (pix.IsHiGainSaturation())881 bad.SetUnsuitable(MBadPixelsPix::kUnsuitableRun);882 return;883 }884 885 if (hist.GetSaturated() > fNumLoGainSaturationLimit*hist.GetHGausHist()->GetEntries())886 {887 *fLog << warn << "Saturated Lo Gain histogram in pixel: " << pix.GetPixId() << endl;888 bad.SetUncalibrated( MBadPixelsPix::kLoGainSaturation );889 bad.SetUnsuitable( MBadPixelsPix::kUnsuitableRun );890 hist.CreateFourierSpectrum();891 return;892 }893 894 //895 // 2) Fit the Lo Gain histograms with a Gaussian896 //897 if (!hist.FitGaus())898 //899 // 3) In case of failure set the bit kFitted to false and take histogram means and RMS900 //901 if (!hist.RepeatFit())902 {903 hist.BypassFit();904 bad.SetUncalibrated( MBadPixelsPix::kLoGainNotFitted );905 if (pix.IsHiGainSaturation())906 bad.SetUnsuitable(MBadPixelsPix::kUnreliableRun);907 }908 909 //910 // 4) Check for pickup911 //912 hist.CountPickup();913 914 //915 // 5) Check for oscillations916 //917 hist.CreateFourierSpectrum();918 919 if (!hist.IsFourierSpectrumOK())920 {921 bad.SetUncalibrated( MBadPixelsPix::kLoGainOscillating );922 if (pix.IsHiGainSaturation())923 bad.SetUnsuitable(MBadPixelsPix::kUnreliableRun);924 }925 //926 // 6) Retrieve the results and store them in this class927 //928 pix.SetLoGainMeanCharge( hist.GetMean() );929 pix.SetLoGainMeanChargeErr( hist.GetMeanErr() );930 pix.SetLoGainSigmaCharge( hist.GetSigma() );931 pix.SetLoGainSigmaChargeErr( hist.GetSigmaErr() );932 pix.SetLoGainChargeProb ( hist.GetProb() );933 934 if (pix.IsHiGainSaturation())935 {936 pix.SetAbsTimeMean ( hist.GetAbsTimeMean());937 pix.SetAbsTimeRms ( hist.GetAbsTimeRms() );938 }939 940 pix.SetLoGainNumPickup ( hist.GetPickup() );941 942 }943 695 944 696 // -------------------------------------------------------------------------- -
trunk/MagicSoft/Mars/mcalib/MHCalibrationChargeCam.h
r3631 r3636 7 7 8 8 class MRawEvtData; 9 class MGeomCam;10 class MBadPixelsCam;11 9 class MBadPixelsPix; 12 10 class MCalibrationChargeCam; … … 17 15 private: 18 16 19 static const Int_t fgAverageNbins; // The default for fAverageNbins20 17 static const Float_t fgNumHiGainSaturationLimit; // The default for fNumHiGainSaturationLimit 21 18 static const Float_t fgNumLoGainSaturationLimit; // The default for fNumLoGainSaturationLimit 22 19 23 Int_t fAverageNbins; // Number of bins for the average histograms24 20 Float_t fNumHiGainSaturationLimit; // Rel. amount sat. higain FADC slices until pixel is called saturated 25 21 Float_t fNumLoGainSaturationLimit; // Rel. amount sat. logain FADC slices until pixel is called saturated 26 22 27 MCalibrationChargeCam *fCam; //! Calibration Cam with the results 28 MRawEvtData *fRawEvt; //! Raw event data 29 MGeomCam *fGeom; //! Camera geometry 30 MBadPixelsCam *fBadPixels; //! Bad Pixels storage container 23 MRawEvtData *fRawEvt; //! Raw event data 31 24 32 void InitializeHiGainHists(MHCalibrationChargePix &hist, MBadPixelsPix &bad, const Int_t i); 33 void InitializeLoGainHists(MHCalibrationChargePix &hist, MBadPixelsPix &bad, const Int_t i); 34 25 Bool_t SetupHists(const MParList *pList); 26 Bool_t ReInitHists(MParList *pList); 27 Bool_t FillHists(const MParContainer *par, const Stat_t w=1); 28 Bool_t FinalizeHists(); 29 void FinalizeBadPixels(); 30 35 31 void FinalizeHiGainHists(MHCalibrationChargePix &hist, MCalibrationChargePix &pix, MBadPixelsPix &bad); 36 32 void FinalizeLoGainHists(MHCalibrationChargePix &hist, MCalibrationChargePix &pix, MBadPixelsPix &bad); … … 41 37 ~MHCalibrationChargeCam() {} 42 38 43 44 void SetAverageNbins( const Int_t bins=fgAverageNbins ) { fAverageNbins = bins; }45 39 void SetNumLoGainSaturationLimit( const Float_t lim=fgNumLoGainSaturationLimit) { fNumLoGainSaturationLimit = lim; } 46 40 void SetNumHiGainSaturationLimit( const Float_t lim=fgNumHiGainSaturationLimit) { fNumHiGainSaturationLimit = lim; } … … 49 43 Float_t GetNumLoGainSaturationLimit() const { return fNumLoGainSaturationLimit; } 50 44 51 Bool_t SetupFill(const MParList *pList);52 Bool_t ReInit ( MParList *pList);53 Bool_t Fill (const MParContainer *par, const Stat_t w=1);54 Bool_t Finalize ( );55 56 45 Bool_t GetPixelContent ( Double_t &val, Int_t idx, const MGeomCam &cam, Int_t type=0) const; 57 46 void DrawPixelContent( Int_t num ) const; 58 47 59 ClassDef(MHCalibrationChargeCam, 1) // Histogram class for camera c alibration48 ClassDef(MHCalibrationChargeCam, 1) // Histogram class for camera charge calibration 60 49 }; 61 50 -
trunk/MagicSoft/Mars/mcalib/MHCalibrationChargeHiGainPix.cc
r3613 r3636 60 60 // 61 61 // Sets: 62 // - the default number for f ChargeNbins(fgChargeNbins)63 // - the default number for f ChargeFirst(fgChargeFirst)64 // - the default number for f ChargeLast(fgChargeLast)62 // - the default number for fNbins (fgChargeNbins) 63 // - the default number for fFirst (fgChargeFirst) 64 // - the default number for fLast (fgChargeLast) 65 65 // - the default number for fAbsTimeNbins (fgAbstTimeNbins) 66 66 // - the default number for fAbsTimeFirst (fgAbsTimeFirst) … … 79 79 fTitle = title ? title : "Fill the FADC sums of the HiGainPix events and perform the fits Pixel "; 80 80 81 Set ChargeNbins();82 Set ChargeFirst();83 Set ChargeLast();81 SetNbins ( fgChargeNbins ); 82 SetFirst ( fgChargeFirst ); 83 SetLast ( fgChargeLast ); 84 84 85 85 SetAbsTimeNbins(); … … 87 87 SetAbsTimeLast(); 88 88 89 fHGausHist.SetName ("HCalibrationChargeHiGainPix");89 fHGausHist.SetName ("HCalibrationChargeHiGainPix"); 90 90 fHGausHist.SetTitle("Distribution of Summed Hi Gain FADC slices Pixel "); 91 91 92 fHAbsTime.SetName ("HAbsTimeHiGainPix");92 fHAbsTime.SetName ("HAbsTimeHiGainPix"); 93 93 fHAbsTime.SetTitle("Distribution of Absolute Arrival Times Hi Gain Pixel "); 94 94 } -
trunk/MagicSoft/Mars/mcalib/MHCalibrationChargeHiGainPix.h
r3613 r3636 12 12 private: 13 13 14 static const Int_t fgChargeNbins; // Default for fChargeNBins(now set to: 2000 )15 static const Axis_t fgChargeFirst; // Default for fChargeFirst(now set to: -0.5 )16 static const Axis_t fgChargeLast; // Default for fChargeLast(now set to: 1999.5)17 static const Int_t fgAbsTimeNbins; // Default for fAbsTimeNbins (now set to: 20 )18 static const Axis_t fgAbsTimeFirst; // Default for fAbsTimeFirst (now set to: -0.5 )19 static const Axis_t fgAbsTimeLast; // Default for fAbsTimeLast (now set to: 19.5 )14 static const Int_t fgChargeNbins; // Default for MHGausEvents::fNbins (now set to: 2000 ) 15 static const Axis_t fgChargeFirst; // Default for MHGausEvents::fFirst (now set to: -0.5 ) 16 static const Axis_t fgChargeLast; // Default for MHGausEvents::fLast (now set to: 1999.5) 17 static const Int_t fgAbsTimeNbins; // Default for fAbsTimeNbins (now set to: 20 ) 18 static const Axis_t fgAbsTimeFirst; // Default for fAbsTimeFirst (now set to: -0.5 ) 19 static const Axis_t fgAbsTimeLast; // Default for fAbsTimeLast (now set to: 19.5 ) 20 20 21 21 public: … … 25 25 26 26 // Setters 27 void SetChargeNbins(const Int_t bins =fgChargeNbins) { fChargeNbins = bins; }28 void SetChargeFirst(const Axis_t first=fgChargeFirst) { fChargeFirst = first; }29 void SetChargeLast (const Axis_t last =fgChargeLast) { fChargeLast = last; }30 31 27 void SetAbsTimeNbins(const Int_t bins =fgAbsTimeNbins) { fAbsTimeNbins = bins; } 32 28 void SetAbsTimeFirst(const Axis_t first=fgAbsTimeFirst) { fAbsTimeFirst = first; } -
trunk/MagicSoft/Mars/mcalib/MHCalibrationChargeLoGainPix.cc
r3614 r3636 60 60 // 61 61 // Sets: 62 // - the default number for f ChargeNbins(fgChargeNbins)63 // - the default number for f ChargeFirst(fgChargeFirst)64 // - the default number for f ChargeLast(fgChargeLast)62 // - the default number for fNbins (fgChargeNbins) 63 // - the default number for fFirst (fgChargeFirst) 64 // - the default number for fLast (fgChargeLast) 65 65 // - the default number for fAbsTimeNbins (fgAbstTimeNbins) 66 66 // - the default number for fAbsTimeFirst (fgAbsTimeFirst) … … 79 79 fTitle = title ? title : "Fill the FADC sums of the Low Gain events and perform the fits Pixel "; 80 80 81 Set ChargeNbins();82 Set ChargeFirst();83 Set ChargeLast();81 SetNbins ( fgChargeNbins ); 82 SetFirst ( fgChargeFirst ); 83 SetLast ( fgChargeLast ); 84 84 85 85 SetAbsTimeNbins(); … … 87 87 SetAbsTimeLast(); 88 88 89 fHGausHist.SetName ("HCalibrationChargeLoGainPix");89 fHGausHist.SetName ("HCalibrationChargeLoGainPix"); 90 90 fHGausHist.SetTitle("Distribution of Summed Lo Gain FADC slices Pixel "); 91 91 92 fHAbsTime.SetName ("HAbsTimeLoGainPix");92 fHAbsTime.SetName ("HAbsTimeLoGainPix"); 93 93 fHAbsTime.SetTitle("Distribution of Absolute Arrival Times Lo Gain Pixel "); 94 94 } -
trunk/MagicSoft/Mars/mcalib/MHCalibrationChargeLoGainPix.h
r3614 r3636 11 11 private: 12 12 13 static const Int_t fgChargeNbins; // Default for fChargeNBins(now set to: 200 )14 static const Axis_t fgChargeFirst; // Default for fChargeFirst(now set to: -0.5 )15 static const Axis_t fgChargeLast; // Default for fChargeLast(now set to: 199.5 )16 static const Int_t fgAbsTimeNbins; // Default for fAbsTimeNbins (now set to: 15 )17 static const Axis_t fgAbsTimeFirst; // Default for fAbsTimeFirst (now set to: -0.5 )18 static const Axis_t fgAbsTimeLast; // Default for fAbsTimeLast (now set to: 14.5 )13 static const Int_t fgChargeNbins; // Default for MHGausEvents::fNbins (now set to: 200 ) 14 static const Axis_t fgChargeFirst; // Default for MHGausEvents::fFirst (now set to: -0.5 ) 15 static const Axis_t fgChargeLast; // Default for MHGausEvents::fLast (now set to: 199.5 ) 16 static const Int_t fgAbsTimeNbins; // Default for fAbsTimeNbins (now set to: 15 ) 17 static const Axis_t fgAbsTimeFirst; // Default for fAbsTimeFirst (now set to: -0.5 ) 18 static const Axis_t fgAbsTimeLast; // Default for fAbsTimeLast (now set to: 14.5 ) 19 19 20 20 public: … … 24 24 25 25 // Setters 26 void SetChargeNbins(const Int_t bins =fgChargeNbins) { fChargeNbins = bins; }27 void SetChargeFirst(const Axis_t first=fgChargeFirst) { fChargeFirst = first; }28 void SetChargeLast (const Axis_t last =fgChargeLast) { fChargeLast = last; }29 30 26 void SetAbsTimeNbins(const Int_t bins =fgAbsTimeNbins) { fAbsTimeNbins = bins; } 31 27 void SetAbsTimeFirst(const Axis_t first=fgAbsTimeFirst) { fAbsTimeFirst = first; } -
trunk/MagicSoft/Mars/mcalib/MHCalibrationChargePix.cc
r3625 r3636 56 56 const Axis_t MHCalibrationChargePix::fgAbsTimeFirst = -0.5; 57 57 const Axis_t MHCalibrationChargePix::fgAbsTimeLast = 14.5; 58 const Float_t MHCalibrationChargePix::fgPickupLimit = 5.;59 58 // -------------------------------------------------------------------------- 60 59 // … … 62 61 // 63 62 // Sets: 64 // - the default number for f ChargeNbins(fgChargeNbins)65 // - the default number for f ChargeFirst(fgChargeFirst)66 // - the default number for f ChargeLast(fgChargeLast)63 // - the default number for fNbins (fgChargeNbins) 64 // - the default number for fFirst (fgChargeFirst) 65 // - the default number for fLast (fgChargeLast) 67 66 // - the default number for fAbsTimeNbins (fgAbstTimeNbins) 68 67 // - the default number for fAbsTimeFirst (fgAbsTimeFirst) 69 68 // - the default number for fAbsTimeLast (fgAbsTimeLast) 70 // - the default number for fPickupLimit (fgPickupLimit)71 69 // 72 70 // - the default name of the fHGausHist ("HCalibrationCharge") … … 88 86 // 89 87 MHCalibrationChargePix::MHCalibrationChargePix(const char *name, const char *title) 90 : fHAbsTime() , fFlags(0)88 : fHAbsTime() 91 89 { 92 90 … … 94 92 fTitle = title ? title : "Statistics of the FADC sums of calibration events"; 95 93 96 Set ChargeNbins();97 Set ChargeFirst();98 Set ChargeLast();94 SetNbins ( fgChargeNbins ); 95 SetFirst ( fgChargeFirst ); 96 SetLast ( fgChargeLast ); 99 97 100 98 SetAbsTimeNbins(); 101 99 SetAbsTimeFirst(); 102 100 SetAbsTimeLast(); 103 104 SetPickupLimit();105 101 106 102 fHGausHist.SetName("HCalibrationCharge"); … … 124 120 // 125 121 // Sets: 126 // - fHGausHist.SetBins(f ChargeNbins,fChargeFirst,fChargeLast);122 // - fHGausHist.SetBins(fNbins,fFirst,fLast); 127 123 // - fHAbsTime.SetBins(fAbsTimeNbins,fAbsTimeFirst,fAbsTimeLast); 128 124 // 129 void MHCalibrationChargePix::Init() 130 { 131 132 fHGausHist.SetBins(fChargeNbins,fChargeFirst,fChargeLast); 125 void MHCalibrationChargePix::InitBins() 126 { 127 MHGausEvents::InitBins(); 133 128 fHAbsTime.SetBins(fAbsTimeNbins,fAbsTimeFirst,fAbsTimeLast); 134 129 } … … 137 132 // 138 133 // Sets: 139 // - all variables to 0., except fPixId to -1 140 // - all flags to kFALSE 134 // - all variables to 0. 141 135 // 142 136 // Executes: … … 146 140 { 147 141 148 fPixId = -1;149 142 fSaturated = 0.; 150 fPickup = 0.; 151 152 SetExcluded ( kFALSE ); 153 143 154 144 MHGausEvents::Clear(); 155 145 return; … … 166 156 // -------------------------------------------------------------------------- 167 157 // 168 // Set fPixId to id158 // Calls MHGausEvents::ChangeHistId() 169 159 // 170 160 // Add id to names and titles of: 171 // - this172 // - fHGausHist173 161 // - fHAbsTime 174 162 // … … 176 164 { 177 165 178 fPixId = id; 179 180 fHGausHist.SetName(Form("%s%d", fHGausHist.GetName(), id)); 181 fHGausHist.SetTitle(Form("%s%d", fHGausHist.GetTitle(), id)); 182 183 fHAbsTime.SetName(Form("%s%d", fHAbsTime.GetName(), id)); 166 MHGausEvents::ChangeHistId(id); 167 168 fHAbsTime.SetName (Form("%s%d", fHAbsTime.GetName(), id)); 184 169 fHAbsTime.SetTitle(Form("%s%d", fHAbsTime.GetTitle(), id)); 185 170 186 fName = Form("%s%d", fName.Data(), id);187 fTitle = Form("%s%d", fTitle.Data(), id);188 }189 190 // --------------------------------------------------------------------------191 //192 // Set Excluded bit from outside193 //194 void MHCalibrationChargePix::SetExcluded(const Bool_t b)195 {196 b ? SETBIT(fFlags,kExcluded) : CLRBIT(fFlags,kExcluded);197 }198 199 const Bool_t MHCalibrationChargePix::IsExcluded() const200 {201 return TESTBIT(fFlags,kExcluded);202 171 } 203 172 … … 317 286 318 287 } 319 320 // -----------------------------------------------------------------------------321 //322 // Repeats the Gauss fit in a smaller range, defined by:323 //324 // min = GetMean() - fPickupLimit * GetSigma();325 // max = GetMean() + fPickupLimit * GetSigma();326 //327 Bool_t MHCalibrationChargePix::RepeatFit(const Option_t *option)328 {329 330 //331 // Get new fitting ranges332 //333 Axis_t rmin = GetMean() - fPickupLimit * GetSigma();334 Axis_t rmax = GetMean() + fPickupLimit * GetSigma();335 336 GetFGausFit()->SetRange(rmin,rmax);337 338 GetHGausHist()->Fit(GetFGausFit(),option);339 340 SetMean ( GetFGausFit()->GetParameter(1) );341 SetMeanErr ( GetFGausFit()->GetParameter(2) );342 SetSigma ( GetFGausFit()->GetParError(1) );343 SetSigmaErr ( GetFGausFit()->GetParError(2) );344 SetProb ( GetFGausFit()->GetProb() );345 346 //347 // The fit result is accepted under condition:348 // 1) The results are not nan's349 // 2) The NDF is not smaller than fNDFLimit (5)350 // 3) The Probability is greater than fProbLimit (default 0.001 == 99.9%)351 //352 if ( TMath::IsNaN ( GetMean() )353 || TMath::IsNaN ( GetMeanErr() )354 || TMath::IsNaN ( GetProb() )355 || TMath::IsNaN ( GetSigma() )356 || TMath::IsNaN ( GetSigmaErr() )357 || GetFGausFit()->GetNDF() < fNDFLimit358 || GetProb() < fProbLimit )359 return kFALSE;360 361 SetGausFitOK(kTRUE);362 return kTRUE;363 364 }365 366 367 368 void MHCalibrationChargePix::CountPickup()369 {370 fPickup = (Int_t)GetHGausHist()->Integral(GetHGausHist()->GetXaxis()->FindBin(GetMean()+fPickupLimit*GetSigma()),371 GetHGausHist()->GetXaxis()->GetLast(),372 "width");373 }374 375 376 377 378 379 -
trunk/MagicSoft/Mars/mcalib/MHCalibrationChargePix.h
r3631 r3636 12 12 private: 13 13 14 static const Int_t fgChargeNbins; // Default for fChargeNBins (now set to: 2000 ) 15 static const Axis_t fgChargeFirst; // Default for fChargeFirst (now set to: -0.5 ) 16 static const Axis_t fgChargeLast; // Default for fChargeLast (now set to: 1999.5) 17 static const Int_t fgAbsTimeNbins; // Default for fAbsTimeNbins (now set to: 15 ) 18 static const Axis_t fgAbsTimeFirst; // Default for fAbsTimeFirst (now set to: -0.5 ) 19 static const Axis_t fgAbsTimeLast; // Default for fAbsTimeLast (now set to: 14.5 ) 20 static const Float_t fgPickupLimit; // Default for fPickupLimit (now set to: 5. ) 14 static const Int_t fgChargeNbins; // Default for MHGausEvents::fNBins (now set to: 2000 ) 15 static const Axis_t fgChargeFirst; // Default for MHGausEvents::fFirst (now set to: -0.5 ) 16 static const Axis_t fgChargeLast; // Default for MHGausEvents::fLast (now set to: 1999.5) 17 static const Int_t fgAbsTimeNbins; // Default for fAbsTimeNbins (now set to: 15 ) 18 static const Axis_t fgAbsTimeFirst; // Default for fAbsTimeFirst (now set to: -0.5 ) 19 static const Axis_t fgAbsTimeLast; // Default for fAbsTimeLast (now set to: 14.5 ) 21 20 22 21 protected: 23 22 24 Int_t fPixId; // The pixel ID25 26 23 TH1F fHAbsTime; // Histogram containing the absolute arrival times 27 24 28 Int_t fChargeNbins; // Number of bins used for the fHGausHist29 Axis_t fChargeFirst; // Lower bound bin used for the fHGausHist30 Axis_t fChargeLast; // Upper bound bin used for the fHGausHist31 32 25 Int_t fAbsTimeNbins; // Number of bins used for the fHAbsTime 33 26 Axis_t fAbsTimeFirst; // Lower bound bin used for the fHAbsTime 34 27 Axis_t fAbsTimeLast; // Upper bound bin used for the fHAbsTime 35 28 36 Float_t fPickupLimit; // Upper number of sigmas from the fitted mean above which events are considered as pickup37 38 29 Float_t fSaturated; // Number of events classified as saturated 39 Float_t fPickup; // Number of events classified as pick-up40 41 Byte_t fFlags; // Bit-field for the flags42 enum { kExcluded }; // Possible bits to be set43 30 44 31 public: … … 49 36 virtual void Clear(Option_t *o=""); 50 37 virtual void Reset(); 51 virtual void Init(); 52 virtual void ChangeHistId(Int_t i); 38 virtual void InitBins(); 53 39 54 40 // Setters 55 virtual void SetChargeNbins(const Int_t bins =fgChargeNbins) { fChargeNbins = bins; }56 virtual void SetChargeFirst(const Axis_t first=fgChargeFirst) { fChargeFirst = first; }57 virtual void SetChargeLast( const Axis_t last =fgChargeLast) { fChargeLast = last; }58 59 41 virtual void SetAbsTimeNbins(const Int_t bins =fgAbsTimeNbins) { fAbsTimeNbins = bins; } 60 42 virtual void SetAbsTimeFirst(const Axis_t first=fgAbsTimeFirst) { fAbsTimeFirst = first; } 61 43 virtual void SetAbsTimeLast( const Axis_t last =fgAbsTimeLast) { fAbsTimeLast = last; } 62 44 63 virtual void SetPickupLimit( const Float_t lim =fgPickupLimit) { fPickupLimit = lim; }64 65 45 void SetSaturated ( const Float_t f ) { fSaturated += f; } 66 void SetExcluded ( const Bool_t b=kTRUE );67 46 68 47 // Getters … … 70 49 const TH1F *GetHAbsTime() const { return &fHAbsTime; } 71 50 72 const Float_t GetIntegral() const;73 74 51 const Float_t GetAbsTimeMean( ) const; 75 52 const Float_t GetAbsTimeRms() const; 76 53 const Float_t GetIntegral() const; 77 54 const Float_t GetSaturated() const { return fSaturated; } 78 const Float_t GetPickup() const { return fPickup; }79 80 const Bool_t IsExcluded() const;81 55 82 56 // Fill histos 83 57 Bool_t FillAbsTime(const Float_t t); 84 58 85 // Fits86 Bool_t RepeatFit(const Option_t *option="RQ0");87 88 59 // Draws 89 60 virtual void Draw(Option_t *opt=""); 90 61 91 62 // Miscelleaneous 92 void C ountPickup();93 63 void ChangeHistId(Int_t id); 64 94 65 ClassDef(MHCalibrationChargePix, 1) // Base Histogram class for a Charge Calibration Pixel 95 66 }; -
trunk/MagicSoft/Mars/mcalib/MHCalibrationRelTimeCam.cc
r3625 r3636 22 22 ! 23 23 \* ======================================================================== */ 24 25 24 ///////////////////////////////////////////////////////////////////////////// 26 25 // // 27 // MHCalibrationRelTimeCam 26 // MHCalibrationRelTimeCam // 28 27 // // 29 // Hold the CalibrationRelTime information for all pixels in the camera // 30 // // 28 // Fills the extracted relative arrival times of MArrivalTimeCam into 29 // the MHGausEvents-classes MHCalibrationRelTimePix for every: 30 // 31 // - pixel, stored in the TObjArray's MHCalibrationCam::fHiGainArray 32 // 33 // - average pixel per area index (e.g. inner and outer for the MAGIC camera), 34 // stored in the TObjArray's MHCalibrationCam::fAverageHiGainAreas 35 // 36 // - average pixel per camera sector (e.g. sectors 1-6 for the MAGIC camera), 37 // stored in the TObjArray's MHCalibrationCam::fAverageHiGainSectors 38 // 39 // Every time is filled into a histogram and an array, in order to perform 40 // a Fourier analysis (see MHGausEvents). The signals are moreover averaged on an 41 // event-by-event basis and written into the corresponding average pixels. 42 // 43 // The histograms are fitted to a Gaussian, mean and sigma with its errors 44 // and the fit probability are extracted. If none of these values are NaN's and 45 // if the probability is bigger than fProbLimit (default: 0.5%), the fit is valid. 46 // Otherwise, the fit is repeated within ranges of the previous mean +- 5 sigma. 47 // In case this does not make the fit valid, the histogram means and RMS's are 48 // taken directly and the following flags are set: 49 // - MBadPixelsPix::SetUncalibrated( MBadPixelsPix::kRelTimeNotFitted ) and 50 // - MBadPixelsPix::SetUnsuitable( MBadPixelsPix::kUnreliableRun ) 51 // 52 // Outliers of more than MHCalibrationChargePix::fPickupLimit (default: 5) sigmas 53 // from the mean are counted as PickUp events (stored in MHCalibrationRelTimePix::fPickup) 54 // 55 // The class also fills arrays with the signal vs. event number, creates a fourier 56 // spectrum and investigates if the projected fourier components follow an exponential 57 // distribution. In case that the probability of the exponential fit is less than 58 // fProbLimit (default: 0.5%), the following flags are set: 59 // - MBadPixelsPix::SetUncalibrated( MBadPixelsPix::kRelTimeOscillating ) and 60 // - MBadPixelsPix::SetUnsuitable( MBadPixelsPix::kUnreliableRun ) 61 // 62 // This same procedure is performed for the average pixels. 63 // 64 // The following results are written into MCalibrationRelTimeCam: 65 // 66 // - MCalibrationPix::SetMean() 67 // - MCalibrationPix::SetMeanErr() 68 // - MCalibrationPix::SetSigma() 69 // - MCalibrationPix::SetSigmaErr() 70 // - MCalibrationPix::SetProb() 71 // - MCalibrationPix::SetNumPickup() 72 // 73 // For all averaged areas, the fitted sigma is multiplied with the square root of 74 // the number involved pixels in order to be able to compare it to the average of 75 // sigmas in the camera. 76 // 31 77 ///////////////////////////////////////////////////////////////////////////// 32 78 #include "MHCalibrationRelTimeCam.h" 79 #include "MHCalibrationRelTimePix.h" 33 80 34 81 #include "MLog.h" … … 37 84 #include "MParList.h" 38 85 39 #include "MHCalibrationRelTimePix.h" 86 #include "MCalibrationRelTimeCam.h" 87 #include "MCalibrationPix.h" 40 88 41 89 #include "MArrivalTimeCam.h" 42 90 #include "MArrivalTimePix.h" 43 91 92 #include "MGeomCam.h" 93 #include "MGeomPix.h" 94 95 #include "MBadPixelsCam.h" 96 #include "MBadPixelsPix.h" 97 44 98 ClassImp(MHCalibrationRelTimeCam); 45 99 46 100 using namespace std; 47 const Float_t MHCalibrationRelTimeCam::fgTimeSliceWidth = 3.3;48 const Int_t MHCalibrationRelTimeCam::fgPulserFrequency = 500;49 101 // -------------------------------------------------------------------------- 50 102 // 51 // Default constructor. Creates a MHCalibrationRelTimePix object for each pixel103 // Default Constructor. 52 104 // 53 105 MHCalibrationRelTimeCam::MHCalibrationRelTimeCam(const char *name, const char *title) … … 55 107 56 108 fName = name ? name : "MHCalibrationRelTimeCam"; 57 fTitle = title ? title : "Histogram container for the relative time calibration of the camera"; 58 59 // 60 // loop over all Pixels and create two histograms 61 // one for the Low and one for the High gain 62 // connect all the histogram with the container fHist 63 // 64 fArray = new TObjArray; 65 fArray->SetOwner(); 66 67 68 SetTimeSliceWidth(); 69 SetPulserFrequency(); 109 fTitle = title ? title : "Histogram class for the relative time calibration of the camera"; 110 70 111 } 71 112 72 113 // -------------------------------------------------------------------------- 73 114 // 74 // Delete the array conatining the pixel pedest information 75 // 76 MHCalibrationRelTimeCam::~MHCalibrationRelTimeCam() 77 { 78 delete fArray; 79 } 80 81 // -------------------------------------- 82 // 83 void MHCalibrationRelTimeCam::Clear(Option_t *o) 84 { 85 86 fArray->ForEach(TObject, Clear)(); 87 return; 88 } 89 90 // -------------------------------------------------------------------------- 91 // 92 // Get i-th pixel (pixel number) 93 // 94 MHCalibrationRelTimePix &MHCalibrationRelTimeCam::operator[](UInt_t i) 95 { 96 return *(MHCalibrationRelTimePix*)(fArray->At(i)); 97 } 98 99 // -------------------------------------------------------------------------- 100 // 101 // Get i-th pixel (pixel number) 102 // 103 const MHCalibrationRelTimePix &MHCalibrationRelTimeCam::operator[](UInt_t i) const 104 { 105 return *(MHCalibrationRelTimePix*)(fArray->At(i)); 106 } 107 108 109 // -------------------------------------------------------------------------- 110 // 111 // Our own clone function is necessary since root 3.01/06 or Mars 0.4 112 // I don't know the reason 113 // 114 TObject *MHCalibrationRelTimeCam::Clone(const char *) const 115 { 116 117 const Int_t n = fArray->GetSize(); 118 119 // 120 // FIXME, this might be done faster and more elegant, by direct copy. 121 // 122 MHCalibrationRelTimeCam *cam = new MHCalibrationRelTimeCam; 123 124 cam->fArray->Expand(n); 125 126 for (int i=0; i<n; i++) 127 { 128 delete (*cam->fArray)[i]; 129 (*cam->fArray)[i] = (*fArray)[i]->Clone(); 130 } 131 132 return cam; 133 } 134 135 136 137 // -------------------------------------------------------------------------- 138 // 139 // 140 Bool_t MHCalibrationRelTimeCam::SetupFill(const MParList *pList) 141 { 115 // Gets or creates the pointers to: 116 // - MCalibrationRelTimeCam 117 // 118 // Searches pointer to: 119 // - MArrivalTimeCam 120 // 121 // Initializes, if empty to MArrivalTimeCam::GetSize() for: 122 // - MHCalibrationCam::fHiGainArray, MHCalibrationCam::fLoGainArray 123 // 124 // Initializes, if empty to MGeomCam::GetNumAreas() for: 125 // - MHCalibrationCam::fAverageHiGainAreas, MHCalibrationCam::fAverageLoGainAreas 126 // 127 // Initializes, if empty to MGeomCam::GetNumSectors() for: 128 // - MHCalibrationCam::fAverageHiGainSectors, MHCalibrationCam::fAverageLoGainSectors 129 // 130 // Calls InitializeHists() for every entry in: 131 // - MHCalibrationCam::fHiGainArray 132 // - MHCalibrationCam::fAverageHiGainAreas 133 // - MHCalibrationCam::fAverageHiGainSectors 134 // 135 // Sets Titles and Names for the Histograms 136 // - MHCalibrationCam::fAverageHiGainAreas 137 // - MHCalibrationCam::fAverageHiGainSectors 138 // 139 // Sets number of bins to MHCalibrationCam::fAverageNbins for: 140 // - MHCalibrationCam::fAverageHiGainAreas, MHCalibrationCam::fAverageLoGainAreas 141 // - MHCalibrationCam::fAverageHiGainSectors, MHCalibrationCam::fAverageLoGainSectors 142 // 143 Bool_t MHCalibrationRelTimeCam::ReInitHists(MParList *pList) 144 { 145 146 fCam = (MCalibrationCam*)pList->FindCreateObj("MCalibrationRelTimeCam"); 147 if (!fCam) 148 return kFALSE; 149 150 MArrivalTimeCam *signal = (MArrivalTimeCam*)pList->FindObject("MArrivalTimeCam"); 151 if (!signal) 152 { 153 *fLog << err << "MArrivalTimeCam not found... abort." << endl; 154 return kFALSE; 155 } 156 157 const Int_t npixels = fGeom->GetNumPixels(); 158 const Int_t nsectors = fGeom->GetNumSectors(); 159 const Int_t nareas = fGeom->GetNumAreas(); 160 161 if (fHiGainArray->GetEntries()==0) 162 { 163 fHiGainArray->Expand(npixels); 164 for (Int_t i=0; i<npixels; i++) 165 { 166 (*fHiGainArray)[i] = new MHCalibrationRelTimePix; 167 InitHists((*this)[i],(*fBadPixels)[i],i); 168 } 169 } 170 171 if (fLoGainArray->GetEntries()==0) 172 { 173 fLoGainArray->Expand(npixels); 174 for (Int_t i=0; i<npixels; i++) 175 { 176 (*fLoGainArray)[i] = new MHCalibrationRelTimePix; 177 InitHists((*this)(i),(*fBadPixels)[i],i); 178 } 179 } 180 181 if (fAverageHiGainAreas->GetEntries()==0) 182 { 183 fAverageHiGainAreas->Expand(nareas); 184 185 for (Int_t j=0; j<nareas; j++) 186 { 187 (*fAverageHiGainAreas)[j] = 188 new MHCalibrationRelTimePix("AverageHiGainArea", 189 "Average Relative Arrival Times area idx "); 190 191 InitHists(GetAverageHiGainArea(j),fCam->GetAverageBadArea(j),j); 192 193 GetAverageHiGainArea(j).GetHGausHist()->SetTitle("Relative Arrival Times average Area Idx "); 194 GetAverageHiGainArea(j).SetNbins(fAverageNbins); 195 } 196 } 197 198 if (fAverageLoGainAreas->GetEntries()==0) 199 { 200 fAverageLoGainAreas->Expand(nareas); 201 202 for (Int_t j=0; j<nareas; j++) 203 { 204 (*fAverageLoGainAreas)[j] = 205 new MHCalibrationRelTimePix("AverageLoGainArea", 206 "Relative Arrival Times average Area idx "); 207 208 InitHists(GetAverageLoGainArea(j),fCam->GetAverageBadArea(j),j); 209 210 GetAverageLoGainArea(j).GetHGausHist()->SetTitle("Relative Arrival Times average Area Idx "); 211 GetAverageLoGainArea(j).SetNbins(fAverageNbins); 212 } 213 } 214 215 if (fAverageHiGainSectors->GetEntries()==0) 216 { 217 fAverageHiGainSectors->Expand(nsectors); 218 219 for (Int_t j=0; j<nsectors; j++) 220 { 221 (*fAverageHiGainSectors)[j] = 222 new MHCalibrationRelTimePix("AverageHiGainSector", 223 "Relative Arrival Times average sector "); 224 225 GetAverageHiGainSector(j).GetHGausHist()->SetTitle("Relative Arrival Times average Sector "); 226 GetAverageHiGainSector(j).SetNbins(fAverageNbins); 227 } 228 } 229 230 if (fAverageLoGainSectors->GetEntries()==0) 231 { 232 fAverageLoGainSectors->Expand(nsectors); 233 234 for (Int_t j=0; j<nsectors; j++) 235 { 236 (*fAverageLoGainSectors)[j] = 237 new MHCalibrationRelTimePix("AverageLoGainSector", 238 "Relative Arrival Times average sector "); 239 240 InitHists(GetAverageLoGainSector(j),fCam->GetAverageBadSector(j),j); 241 242 GetAverageLoGainSector(j).GetHGausHist()->SetTitle("Relative Arrival Times average Sector "); 243 GetAverageLoGainSector(j).SetNbins(fAverageNbins); 244 } 245 } 142 246 143 247 return kTRUE; 144 145 } 146 147 // -------------------------------------------------------------------------- 148 // 149 Bool_t MHCalibrationRelTimeCam::Fill(const MParContainer *par, const Stat_t w) 248 } 249 250 251 //-------------------------------------------------------------------------------- 252 // 253 // Retrieves pointer to MArrivalTimeCam: 254 // 255 // Retrieves from MGeomCam: 256 // - number of pixels 257 // - number of pixel areas 258 // - number of sectors 259 // 260 // Fill RelTimes histograms (MHGausEvents::FillHistAndArray()) with: 261 // - MArrivalTimePix::GetArrivalTime(pixid) - MArrivalTimePix::GetArrivalTime(1); 262 // (i.e. the time difference between pixel i and pixel 1 (hardware number: 2) 263 // 264 Bool_t MHCalibrationRelTimeCam::FillHists(const MParContainer *par, const Stat_t w) 150 265 { 151 266 … … 157 272 } 158 273 159 const Int_t n = arrtime->GetSize(); 160 161 if (fArray->GetEntries()==0) 162 { 163 fArray->Expand(n); 274 const Int_t npixels = fGeom->GetNumPixels(); 275 const Int_t nareas = fGeom->GetNumAreas(); 276 const Int_t nsectors = fGeom->GetNumSectors(); 277 278 Float_t sumareahi [nareas], sumarealo [nareas]; 279 Float_t sumsectorhi[nsectors], sumsectorlo[nsectors]; 280 Int_t numareahi [nareas], numarealo [nareas]; 281 Int_t numsectorhi[nsectors], numsectorlo[nsectors]; 282 283 for (UInt_t j=0; j<nareas; j++) 284 { 285 sumareahi[j] = sumarealo[j] = 0.; 286 numareahi[j] = numarealo[j] = 0; 287 } 288 289 for (UInt_t j=0; j<nsectors; j++) 290 { 291 sumsectorhi[j] = sumsectorlo[j] = 0.; 292 numsectorhi[j] = numsectorlo[j] = 0; 293 } 294 295 const MArrivalTimePix &refpix = (*arrtime)[1]; 296 const Float_t reftime = refpix.IsLoGainUsed() 297 ? refpix.GetArrivalTimeLoGain() : refpix.GetArrivalTimeHiGain(); 298 299 for (Int_t i=0; i<npixels; i++) 300 { 301 302 MHGausEvents &histhi = (*this)[i]; 303 MHGausEvents &histlo = (*this)(i); 304 305 if (histhi.IsExcluded()) 306 continue; 307 308 const MArrivalTimePix &pix = (*arrtime)[i]; 309 const Float_t time = pix.IsLoGainUsed() 310 ? pix.GetArrivalTimeLoGain() : pix.GetArrivalTimeHiGain(); 311 const Float_t reltime = time - reftime; 164 312 165 for (Int_t i=0; i<n; i++) 313 const Int_t aidx = (*fGeom)[i].GetAidx(); 314 const Int_t sector = (*fGeom)[i].GetSector(); 315 316 if (pix.IsLoGainUsed()) 317 { 318 histlo.FillHistAndArray(reltime); 319 sumarealo [aidx] += reltime; 320 numarealo [aidx] ++; 321 sumsectorlo[sector] += reltime; 322 numsectorlo[sector] ++; 323 } 324 else 166 325 { 167 (*fArray)[i] = new MHCalibrationRelTimePix;168 MHCalibrationRelTimePix &hist = (*this)[i];169 hist.ChangeHistId(i);170 hist.Init();171 hist.SetEventFrequency(fPulserFrequency);326 histhi.FillHistAndArray(reltime) ; 327 sumareahi [aidx] += reltime; 328 numareahi [aidx] ++; 329 sumsectorhi[sector] += reltime; 330 numsectorhi[sector] ++; 172 331 } 173 332 } 174 333 175 if (fArray->GetEntries() != n)176 { 177 gLog << err << "ERROR - Size mismatch... abort." << endl;178 return kFALSE;179 } 180 181 const MArrivalTimePix &refpix = (*arrtime)[1];182 const Float_t reftime = refpix.IsLoGainUsed() ? refpix.GetArrivalTimeLoGain() : refpix.GetArrivalTimeHiGain();183 184 for ( Int_t i=0; i<n; i++)185 { 186 const MArrivalTimePix &pix = (*arrtime)[i];187 const Float_t time = pix.IsLoGainUsed() ? pix.GetArrivalTimeLoGain() : pix.GetArrivalTimeHiGain(); 188 189 MHCalibrationRelTimePix &hist = (*this)[i];190 hist.FillHistAndArray(time - reftime);191 } 192 334 for (UInt_t j=0; j<nareas; j++) 335 { 336 MHGausEvents &histhi = GetAverageHiGainArea(j); 337 histhi.FillHistAndArray(numareahi[j] == 0 ? 0. : sumareahi[j]/numareahi[j]); 338 339 MHGausEvents &histlo = GetAverageLoGainArea(j); 340 histlo.FillHistAndArray(numarealo[j] == 0 ? 0. : sumarealo[j]/numarealo[j]); 341 } 342 343 for (UInt_t j=0; j<nsectors; j++) 344 { 345 MHGausEvents &histhi = GetAverageHiGainSector(j); 346 histhi.FillHistAndArray(numsectorhi[j] == 0 ? 0. : sumsectorhi[j]/numsectorhi[j]); 347 348 MHGausEvents &histlo = GetAverageLoGainSector(j); 349 histlo.FillHistAndArray(numsectorlo[j] == 0 ? 0. : sumsectorlo[j]/numsectorlo[j]); 350 } 351 193 352 return kTRUE; 194 353 } 195 354 196 Bool_t MHCalibrationRelTimeCam::Finalize() 197 { 198 199 for (Int_t i=0; i<fArray->GetSize(); i++) 355 // -------------------------------------------------------------------------- 356 // 357 // Calls: 358 // - MHCalibrationCam::FitHiGainArrays() with flags: 359 // MBadPixelsPix::kRelTimeNotFitted and MBadPixelsPix::kRelTimeOscillating 360 // - MHCalibrationCam::FitLoGainArrays() with flags: 361 // MBadPixelsPix::kRelTimeNotFitted and MBadPixelsPix::kRelTimeOscillating 362 // 363 Bool_t MHCalibrationRelTimeCam::FinalizeHists() 364 { 365 366 FitHiGainArrays((MCalibrationCam&)(*fCam),*fBadPixels, 367 MBadPixelsPix::kRelTimeNotFitted, 368 MBadPixelsPix::kRelTimeOscillating); 369 370 FitLoGainArrays((MCalibrationCam&)(*fCam),*fBadPixels, 371 MBadPixelsPix::kRelTimeNotFitted, 372 MBadPixelsPix::kRelTimeOscillating); 373 374 return kTRUE; 375 } 376 377 // -------------------------------------------------------------------------- 378 // 379 // Sets all pixels to MBadPixelsPix::kUnreliableRun, if following flags are set: 380 // - MBadPixelsPix::kRelTimeNotFitted 381 // - MBadPixelsPix::kRelTimeOscillating 382 // 383 void MHCalibrationRelTimeCam::FinalizeBadPixels() 384 { 385 386 for (Int_t i=0; i<fBadPixels->GetSize(); i++) 200 387 { 201 388 202 MHCalibrationRelTimePix &hist = (*this)[i]; 389 MBadPixelsPix &bad = (*fBadPixels)[i]; 390 391 if (bad.IsUncalibrated( MBadPixelsPix::kRelTimeNotFitted )) 392 bad.SetUnsuitable( MBadPixelsPix::kUnreliableRun ); 393 394 if (bad.IsUncalibrated( MBadPixelsPix::kRelTimeOscillating)) 395 bad.SetUnsuitable( MBadPixelsPix::kUnreliableRun ); 203 396 204 // 205 // 1) Return if the charge distribution is already succesfully fitted 206 // or if the histogram is empty 207 // 208 if (hist.IsGausFitOK() || hist.IsEmpty()) 209 continue; 210 211 // 212 // 2) Fit the Hi Gain histograms with a Gaussian 213 // 214 hist.FitGaus(); 215 216 // 217 // 3) If fit does not succeed , bypass the fit and take the histogram means and sigmas 218 // 219 if (!hist.IsGausFitOK()) 220 hist.BypassFit(); 221 222 // 223 // 4) Create the fourier transform of the arrays 224 // 225 hist.CreateFourierSpectrum(); 226 227 // 228 // 5) Renormalize to the real time in ns. 229 // 230 hist.Renorm(fTimeSliceWidth); 231 232 } 233 return kTRUE; 397 } 234 398 } 235 399 … … 250 414 // ============================================= 251 415 // 252 // 4: Returned probability of Gauss fit to Charge distribution416 // 4: Returned probability of Gauss fit to RelTime distribution 253 417 // 254 418 // Localized defects: … … 261 425 { 262 426 263 if (f Array->GetSize() <= idx)427 if (fHiGainArray->GetSize() <= idx) 264 428 return kFALSE; 265 429 430 MHCalibrationRelTimePix &hist = (MHCalibrationRelTimePix&)(*this)[idx]; 431 266 432 switch (type) 267 433 { 268 434 case 0: 269 val = (*this)[idx].GetMean();435 val = hist.GetMean(); 270 436 break; 271 437 case 1: 272 val = (*this)[idx].GetMeanErr();438 val = hist.GetMeanErr(); 273 439 break; 274 440 case 2: 275 val = (*this)[idx].GetSigma();441 val = hist.GetSigma(); 276 442 break; 277 443 case 3: 278 val = (*this)[idx].GetSigmaErr();444 val = hist.GetSigmaErr(); 279 445 break; 280 446 case 4: 281 val = (*this)[idx].GetProb();447 val = hist.GetProb(); 282 448 break; 283 449 case 5: 284 if (! (*this)[idx].IsGausFitOK())450 if (!hist.IsGausFitOK()) 285 451 val = 1.; 286 452 break; 287 453 case 6: 288 if (! (*this)[idx].IsFourierSpectrumOK())454 if (!hist.IsFourierSpectrumOK()) 289 455 val = 1.; 290 456 break; … … 297 463 void MHCalibrationRelTimeCam::DrawPixelContent(Int_t idx) const 298 464 { 299 (*this)[idx].DrawClone(); 300 } 465 const MHCalibrationRelTimePix &hist = (MHCalibrationRelTimePix&)(*this)[idx]; 466 hist.DrawClone(); 467 } -
trunk/MagicSoft/Mars/mcalib/MHCalibrationRelTimeCam.h
r3625 r3636 2 2 #define MARS_MHCalibrationRelTimeCam 3 3 4 #ifndef ROOT_TObjArray5 #include <TObjArray.h>4 #ifndef MARS_MHCalibrationCam 5 #include "MHCalibrationCam.h" 6 6 #endif 7 7 8 #ifndef MARS_MH 9 #include "MH.h" 10 #endif 11 #ifndef MARS_MCamEvent 12 #include "MCamEvent.h" 13 #endif 14 15 class MHCalibrationRelTimePix; 16 class MHCalibrationRelTimeCam : public MH, public MCamEvent 8 class MGeomCam; 9 class MHCalibrationRelTimeCam : public MHCalibrationCam 17 10 { 18 11 19 12 private: 20 13 21 static const Float_t fgTimeSliceWidth; // Default for fTimeSliceWidth 22 static const Int_t fgPulserFrequency; // Default for fPulserFrequency 23 24 Float_t fTimeSliceWidth; // FADC slice time width 25 Int_t fPulserFrequency; // The pulser frequency 14 Bool_t ReInitHists(MParList *pList); 15 Bool_t FillHists(const MParContainer *par, const Stat_t w=1); 16 Bool_t FinalizeHists(); 17 void FinalizeBadPixels(); 26 18 27 TObjArray *fArray; //-> List of MHCalibrationRelTimePix's28 29 19 public: 30 20 31 21 MHCalibrationRelTimeCam(const char *name=NULL, const char *title=NULL); 32 ~MHCalibrationRelTimeCam() ;22 ~MHCalibrationRelTimeCam() {} 33 23 34 void Clear(Option_t *o="");35 36 MHCalibrationRelTimePix &operator[](UInt_t i);37 const MHCalibrationRelTimePix &operator[](UInt_t i) const;38 39 Bool_t SetupFill(const MParList *pList);40 Bool_t Fill(const MParContainer *par, const Stat_t w=1);41 Bool_t Finalize();42 43 // Setters44 void SetTimeSliceWidth( const Float_t width=fgTimeSliceWidth) { fTimeSliceWidth = width; }45 void SetPulserFrequency( const Int_t f=fgPulserFrequency) { fPulserFrequency = f; }46 47 TObject *Clone(const char *) const;48 49 24 Bool_t GetPixelContent(Double_t &val, Int_t idx, const MGeomCam &cam, Int_t type=0) const; 50 25 void DrawPixelContent(Int_t idx) const; -
trunk/MagicSoft/Mars/mcalib/MHCalibrationRelTimePix.cc
r3625 r3636 40 40 using namespace std; 41 41 // 42 const Int_t MHCalibrationRelTimePix::fg ChargeNbins = 900;43 const Axis_t MHCalibrationRelTimePix::fg ChargeFirst = -13.;44 const Axis_t MHCalibrationRelTimePix::fg ChargeLast = 13.;42 const Int_t MHCalibrationRelTimePix::fgRelTimeNbins = 900; 43 const Axis_t MHCalibrationRelTimePix::fgRelTimeFirst = -13.; 44 const Axis_t MHCalibrationRelTimePix::fgRelTimeLast = 13.; 45 45 // -------------------------------------------------------------------------- 46 46 // … … 48 48 // 49 49 // Sets: 50 // - the default number for f ChargeNbins (fgChargeNbins)51 // - the default number for f ChargeFirst (fgChargeFirst)52 // - the default number for f ChargeLast (fgChargeLast)50 // - the default number for fNbins (fgRelTimeNbins) 51 // - the default number for fFirst (fgRelTimeFirst) 52 // - the default number for fLast (fgRelTimeLast) 53 53 // 54 54 // - the default name of the fHGausHist ("HCalibrationRelTime") … … 64 64 // 65 65 MHCalibrationRelTimePix::MHCalibrationRelTimePix(const char *name, const char *title) 66 : fPixId(-1)67 66 { 68 67 … … 70 69 fTitle = title ? title : "Histogrammed Calibration Relative Arrival Time events"; 71 70 72 Set ChargeNbins();73 Set ChargeFirst();74 Set ChargeLast();71 SetNbins ( fgRelTimeNbins ); 72 SetFirst ( fgRelTimeFirst ); 73 SetLast ( fgRelTimeLast ); 75 74 76 75 // Create a large number of bins, later we will rebin … … 96 95 97 96 fPixId = -1; 97 98 98 MHGausEvents::Clear(); 99 99 return; … … 109 109 } 110 110 111 // --------------------------------------------------------------------------112 //113 // Sets:114 // - fHGausHist.SetBins(fChargeNbins,fChargeFirst,fChargeLast);115 //116 void MHCalibrationRelTimePix::Init()117 {118 111 119 fHGausHist.SetBins(fChargeNbins,fChargeFirst,fChargeLast);120 121 }122 123 // --------------------------------------------------------------------------124 //125 // - Set fPixId to id126 //127 // Add id to names and titles of:128 // - fHGausHist129 //130 void MHCalibrationRelTimePix::ChangeHistId(Int_t id)131 {132 133 fPixId = id;134 135 fHGausHist.SetName( Form("%s%d", fHGausHist.GetName(), id));136 fHGausHist.SetTitle( Form("%s%d", fHGausHist.GetTitle(), id));137 138 }139 140 141 // ----------------------------------------------------------------------142 //143 // Renorm the results from FADC slices to times.144 // The parameters slicewidth is the inverse of the FADC frequency and has the unit ns.145 //146 void MHCalibrationRelTimePix::Renorm(const Float_t slicewidth)147 {148 149 SetMean( GetMean() * slicewidth );150 SetMeanErr( GetMeanErr() * slicewidth );151 SetSigma( GetSigma() * slicewidth );152 SetSigmaErr( GetSigmaErr()* slicewidth );153 154 }155 -
trunk/MagicSoft/Mars/mcalib/MHCalibrationRelTimePix.h
r3625 r3636 11 11 private: 12 12 13 static const Int_t fgChargeNbins; // Default for fChargeNBins (now set to: 900 ) 14 static const Axis_t fgChargeFirst; // Default for fChargeFirst (now set to: -13.5 ) 15 static const Axis_t fgChargeLast; // Default for fChargeLast (now set to: 13.5 ) 16 17 Int_t fChargeNbins; // Number of bins used for the fHGausHist 18 Axis_t fChargeFirst; // Lower bound bin used for the fHGausHist 19 Axis_t fChargeLast; // Upper bound bin used for the fHGausHist 20 21 Int_t fPixId; // The pixel ID 13 static const Int_t fgRelTimeNbins; //! Default for MHGausEvents::fNBins (now set to: 900 ) 14 static const Axis_t fgRelTimeFirst; //! Default for MHGausEvents::fFirst (now set to: -13.5 ) 15 static const Axis_t fgRelTimeLast; //! Default for MHGausEvents::fLast (now set to: 13.5 ) 22 16 23 17 public: … … 26 20 ~MHCalibrationRelTimePix() {} 27 21 28 29 22 void Clear(Option_t *o=""); 30 23 void Reset(); 31 void Init();32 33 // Setters34 void SetChargeNbins(const Int_t bins =fgChargeNbins) { fChargeNbins = bins; }35 void SetChargeFirst(const Axis_t first=fgChargeFirst) { fChargeFirst = first; }36 void SetChargeLast( const Axis_t last =fgChargeLast) { fChargeLast = last; }37 38 // Others39 void ChangeHistId(Int_t i);40 void Renorm(const Float_t slicewidth);41 24 42 25 ClassDef(MHCalibrationRelTimePix, 1) // Histogram class for Relative Time Calibration pixel -
trunk/MagicSoft/Mars/mcalib/MHGausEvents.cc
r3625 r3636 40 40 // It is created with a default name and title and resides in directory NULL. 41 41 // - Any class deriving from MHGausEvents needs to apply a binning to fHGausHist 42 // (e.g. with the function TH1F::SetBins(..) ) 42 // (e.g. by setting the variables fNbins, fFirst, fLast and calling the function 43 // InitBins() or by directly calling fHGausHist.SetBins(..) ) 43 44 // - The histogram is filled with the functions FillHist() or FillHistAndArray(). 44 45 // - The histogram can be fitted with the function FitGaus(). This involves stripping … … 48 49 // The NDF is compared to fNDFLimit and a check is made whether results are NaNs. 49 50 // Accordingly, a flag IsGausFitOK() is set. 51 // - One can repeat the fit within a given amount of sigmas from the previous mean 52 // (specified by the variable fPickupLimit) with the function RepeatFit() 53 // - One can completely skip the fitting to set mean, sigma and its errors directly 54 // from the histograms with the function BypassFit() 50 55 // 51 56 // 2) The TArrayF fEvents: … … 96 101 const Float_t MHGausEvents::fgProbLimit = 0.001; 97 102 const Int_t MHGausEvents::fgNDFLimit = 2; 103 const Float_t MHGausEvents::fgPickupLimit = 5.; 98 104 const Int_t MHGausEvents::fgPowerProbabilityBins = 20; 99 105 const Int_t MHGausEvents::fgBinsAfterStripping = 40; … … 105 111 // - the default Probability Limit for fProbLimit (fgProbLimit) 106 112 // - the default NDF Limit for fNDFLimit (fgNDFLimit) 113 // - the default number for fPickupLimit (fgPickupLimit) 107 114 // - the default number of bins after stripping for fBinsAfterStipping (fgBinsAfterStipping) 108 115 // - the default name of the fHGausHist ("HGausHist") 109 116 // - the default title of the fHGausHist ("Histogram of Events with Gaussian Distribution") 110 117 // - the default directory of the fHGausHist (NULL) 118 // - the default for fNbins (100) 119 // - the default for fFirst (0.) 120 // - the default for fLast (100.) 111 121 // 112 122 // Initializes: … … 121 131 fPowerSpectrum(NULL), 122 132 fGraphEvents(NULL), fGraphPowerSpectrum(NULL), 133 fNbins(100), fFirst(0.), fLast(100.), fPixId(-1), 123 134 fHGausHist(), fEvents(0), 124 135 fFGausFit(NULL), fFExpFit(NULL) … … 133 144 SetProbLimit(); 134 145 SetNDFLimit(); 146 SetPickupLimit(); 135 147 SetBinsAfterStripping(); 136 148 … … 183 195 // -------------------------------------------------------------------------- 184 196 // 197 // Default Clear(), can be overloaded. 198 // 185 199 // Sets: 186 200 // - all other pointers to NULL 187 // - all variables to 0., except f EventFrequency201 // - all variables to 0., except fPixId to -1 and keep fEventFrequency 188 202 // - all flags to kFALSE 189 203 // … … 197 211 SetExpFitOK ( kFALSE ); 198 212 SetFourierSpectrumOK( kFALSE ); 213 SetExcluded ( kFALSE ); 199 214 200 215 fMean = 0.; … … 205 220 206 221 fCurrentSize = 0; 222 fPixId = -1; 207 223 208 224 if (fHPowerProbability) … … 247 263 248 264 // -------------------------------------------------------------------------- 265 // 266 // Default Reset(), can be overloaded. 249 267 // 250 268 // Executes: … … 264 282 // -------------------------------------------------------------------------- 265 283 // 284 // Default InitBins, can be overloaded. 285 // 286 // Executes: 287 // - fHGausHist.SetBins(fNbins,fFirst,fLast) 288 // 289 void MHGausEvents::InitBins() 290 { 291 fHGausHist.SetBins(fNbins,fFirst,fLast); 292 } 293 294 // -------------------------------------------------------------------------- 295 // 266 296 // Executes: 267 297 // - FillArray() … … 303 333 } 304 334 335 // -------------------------------------------------------------------------- 336 // 337 // Set Excluded bit from outside 338 // 339 void MHGausEvents::SetExcluded(const Bool_t b) 340 { 341 b ? SETBIT(fFlags,kExcluded) : CLRBIT(fFlags,kExcluded); 342 } 343 305 344 306 345 const Double_t MHGausEvents::GetChiSquare() const … … 366 405 { 367 406 return TESTBIT(fFlags,kExpFitOK); 407 } 408 409 const Bool_t MHGausEvents::IsExcluded() const 410 { 411 return TESTBIT(fFlags,kExcluded); 368 412 } 369 413 … … 409 453 if (fFExpFit) 410 454 return; 455 456 if (fEvents.GetSize() < 8) 457 { 458 *fLog << warn << "Cannot create Fourier spectrum in pixel: " << fPixId 459 << ". Number of events smaller than 8 " << endl; 460 return; 461 } 462 411 463 412 464 // … … 424 476 // This cuts only the non-used zero's, but MFFT will later cut the rest 425 477 MArray::StripZeros(fEvents); 478 479 if (fEvents.GetSize() < 8) 480 { 481 *fLog << warn << "Cannot create Fourier spectrum. " << endl; 482 *fLog << warn << "Number of events (after stripping of zeros) is smaller than 8 " 483 << "in pixel: " << fPixId << endl; 484 return; 485 } 426 486 427 487 MFFT fourier; … … 506 566 if (!fFGausFit) 507 567 { 508 *fLog << warn << dbginf << "WARNING: Could not create fit function for Gauss fit" << endl; 568 *fLog << warn << dbginf << "WARNING: Could not create fit function for Gauss fit " 569 << "in pixel: " << fPixId << endl; 509 570 return kFALSE; 510 571 } … … 528 589 // The fit result is accepted under condition: 529 590 // 1) The results are not nan's 530 // 2) The NDF is not smaller than fNDFLimit ( 5)531 // 3) The Probability is greater than fProbLimit (default 0.001 == 99.9%)591 // 2) The NDF is not smaller than fNDFLimit (default: fgNDFLimit) 592 // 3) The Probability is greater than fProbLimit (default: fgProbLimit) 532 593 // 533 594 if ( TMath::IsNaN(fMean) … … 545 606 546 607 // ----------------------------------------------------------------------------- 608 // 609 // If flag IsGausFitOK() is set (histogram already successfully fitted), 610 // returns kTRUE 611 // 612 // If both fMean and fSigma are still zero, call FitGaus() 613 // 614 // Repeats the Gauss fit in a smaller range, defined by: 615 // 616 // min = GetMean() - fPickupLimit * GetSigma(); 617 // max = GetMean() + fPickupLimit * GetSigma(); 618 // 619 // The fit results are retrieved and stored in class-own variables. 620 // 621 // A flag IsGausFitOK() is set according to whether the fit probability 622 // is smaller or bigger than fProbLimit, whether the NDF is bigger than 623 // fNDFLimit and whether results are NaNs. 624 // 625 Bool_t MHGausEvents::RepeatFit(const Option_t *option) 626 { 627 628 if (IsGausFitOK()) 629 return kTRUE; 630 631 if ((fMean == 0.) && (fSigma == 0.)) 632 return FitGaus(); 633 634 // 635 // Get new fitting ranges 636 // 637 Axis_t rmin = fMean - fPickupLimit * fSigma; 638 Axis_t rmax = fMean + fPickupLimit * fSigma; 639 640 Axis_t hmin = fHGausHist.GetBinCenter(fHGausHist.GetXaxis()->GetFirst()); 641 Axis_t hmax = fHGausHist.GetBinCenter(fHGausHist.GetXaxis()->GetLast()) ; 642 643 fFGausFit->SetRange(hmin < rmin ? rmin : hmin , hmax > rmax ? rmax : hmax); 644 645 fHGausHist.Fit(fFGausFit,option); 646 647 fMean = fFGausFit->GetParameter(1); 648 fMeanErr = fFGausFit->GetParameter(2); 649 fSigma = fFGausFit->GetParError(1) ; 650 fSigmaErr = fFGausFit->GetParError(2) ; 651 fProb = fFGausFit->GetProb() ; 652 653 // 654 // The fit result is accepted under condition: 655 // 1) The results are not nan's 656 // 2) The NDF is not smaller than fNDFLimit (default: fgNDFLimit) 657 // 3) The Probability is greater than fProbLimit (default: fgProbLimit) 658 // 659 if ( TMath::IsNaN ( fMean ) 660 || TMath::IsNaN ( fMeanErr ) 661 || TMath::IsNaN ( fProb ) 662 || TMath::IsNaN ( fSigma ) 663 || TMath::IsNaN ( fSigmaErr ) 664 || fFGausFit->GetNDF() < fNDFLimit 665 || fProb < fProbLimit ) 666 return kFALSE; 667 668 SetGausFitOK(kTRUE); 669 return kTRUE; 670 671 } 672 673 // ----------------------------------------------------------------------------- 547 674 // 548 675 // Bypasses the Gauss fit by taking mean and RMS from the histogram … … 559 686 if (entries <= 0.) 560 687 { 561 *fLog << warn << GetDescriptor() << ": Cannot bypass fit. Number of entries smaller or equal 0" << endl; 688 *fLog << warn << GetDescriptor() 689 << ": Cannot bypass fit. Number of entries smaller or equal 0 in pixel: " << fPixId << endl; 562 690 return; 563 691 } 564 692 565 SetMean ( fHGausHist.GetMean());566 SetMeanErr ( fHGausHist.GetRMS() / TMath::Sqrt(entries));567 SetSigma ( fHGausHist.GetRMS() );568 SetSigmaErr ( fHGausHist.GetRMS() / TMath::Sqrt(entries) / 2. );693 fMean = fHGausHist.GetMean(); 694 fMeanErr = fHGausHist.GetRMS() / TMath::Sqrt(entries); 695 fSigma = fHGausHist.GetRMS() ; 696 fSigmaErr = fHGausHist.GetRMS() / TMath::Sqrt(entries) / 2.; 569 697 } 570 698 … … 577 705 578 706 *fLog << all << endl; 579 *fLog << all << "Results of the Gauss Fit : "<< endl;707 *fLog << all << "Results of the Gauss Fit in pixel: " << fPixId << endl; 580 708 *fLog << all << "Mean: " << GetMean() << endl; 581 709 *fLog << all << "Sigma: " << GetSigma() << endl; … … 585 713 *fLog << all << endl; 586 714 715 } 716 717 // -------------------------------------------------------------------------- 718 // 719 // - Set fPixId to id 720 // 721 // Add id to names and titles of: 722 // - fHGausHist 723 // 724 void MHGausEvents::ChangeHistId(Int_t id) 725 { 726 727 fPixId = id; 728 729 fHGausHist.SetName( Form("%s%d", fHGausHist.GetName(), id)); 730 fHGausHist.SetTitle( Form("%s%d", fHGausHist.GetTitle(), id)); 731 732 fName = Form("%s%d", fName.Data(), id); 733 fTitle = Form("%s%d", fTitle.Data(), id); 734 587 735 } 588 736 … … 766 914 } 767 915 768 916 const Double_t MHGausEvents::GetPickup() const 917 { 918 919 if ((fMean == 0.) && (fSigma == 0.)) 920 return -1.; 921 922 const Int_t first = fHGausHist.GetXaxis()->FindBin(fMean+fPickupLimit*fSigma); 923 const Int_t last = fHGausHist.GetXaxis()->GetLast(); 924 925 if (first >= last) 926 return 0.; 927 928 return fHGausHist.Integral(first, last, "width"); 929 } 930 931 -
trunk/MagicSoft/Mars/mcalib/MHGausEvents.h
r3625 r3636 20 20 private: 21 21 22 const static Float_t fgProbLimit; // Default for fProbLimit (now 0.001) 23 const static Int_t fgNDFLimit; // Default for fNDFLimit (now 2) 24 const static Int_t fgPowerProbabilityBins; // Default for fPowerProbabilityBins (now 20) 25 const static Int_t fgBinsAfterStripping; // Default for fBinsAfterStripping (now 40) 26 22 const static Float_t fgProbLimit; //! Default for fProbLimit (now set to: 0.001) 23 const static Int_t fgNDFLimit; //! Default for fNDFLimit (now set to: 2) 24 const static Float_t fgPickupLimit; //! Default for fPickupLimit (now set to: 5. ) 25 const static Int_t fgPowerProbabilityBins; //! Default for fPowerProbabilityBins (now set to: 20) 26 const static Int_t fgBinsAfterStripping; //! Default for fBinsAfterStripping (now set to: 40) 27 27 28 Int_t fPowerProbabilityBins; // Bins for the projected power spectrum 28 29 Int_t fBinsAfterStripping; // Bins for the Gauss Histogram after stripping off the zeros at both ends … … 41 42 Double_t fProb; // Probability of the Gauss fit (derived from Chi-Square and NDF 42 43 43 enum { kGausFitOK, kExpFitOK, kFourierSpectrumOK }; // Bits to hold information about fit results 44 enum { kGausFitOK, kExpFitOK, kFourierSpectrumOK, 45 kExcluded }; // Bits to hold information about fit results 44 46 45 Byte_t fFlags; // 47 Byte_t fFlags; // Byte to hold the bits fit result bits 46 48 47 Int_t fCurrentSize; // 48 49 Int_t fCurrentSize; // Current size of the array fEvents 50 49 51 protected: 50 52 53 Int_t fNbins; // Number histogram bins for fHGausHist (used by Init()) 54 Axis_t fFirst; // Lower histogram edge for fHGausHist (used by Init()) 55 Axis_t fLast; // Upper histogram edge for fHGausHist (used by Init()) 56 Int_t fPixId; // Pixel ID 57 51 58 TH1F fHGausHist; // Histogram which should hold the Gaussian distribution 52 59 TArrayF fEvents; // Array which holds the entries of GausHist … … 57 64 Float_t fProbLimit; // Probability limit for judgement if fit is OK 58 65 Int_t fNDFLimit; // NDF limit for judgement if fit is OK 66 Float_t fPickupLimit; // Upper number of sigmas from the mean until events are considered as pickup 59 67 60 68 Float_t *CreateEventXaxis(Int_t n); // Create an x-axis for the Event TGraphs … … 75 83 virtual void Clear(Option_t *o=""); 76 84 virtual void Reset(); 77 85 virtual void InitBins(); 86 78 87 // Setters 79 88 void SetEventFrequency(const Float_t f) { fEventFrequency = f; } … … 84 93 void SetSigmaErr( const Double_t d ) { fSigmaErr = d; } 85 94 void SetProb ( const Double_t d ) { fProb = d; } 95 void SetPixId ( const Int_t i ) { fPixId = i; } 96 97 void SetNbins ( const Int_t i ) { fNbins = i; } 98 void SetFirst ( const Double_t d ) { fFirst = d; } 99 void SetLast ( const Double_t d ) { fLast = d; } 86 100 101 void SetNDFLimit( const Int_t lim=fgNDFLimit ) { fNDFLimit = lim; } 102 void SetPickupLimit( const Float_t lim =fgPickupLimit) { fPickupLimit = lim; } 87 103 void SetProbLimit( const Float_t lim=fgProbLimit ) { fProbLimit = lim; } 88 void SetNDFLimit( const Int_t lim=fgNDFLimit ) { fNDFLimit = lim; }89 104 90 105 // Setters ONLY for MC: 91 void SetGausFitOK( const Bool_t b ); 92 void SetExpFitOK( const Bool_t b ); 93 void SetFourierSpectrumOK( const Bool_t b ); 106 void SetGausFitOK( const Bool_t b=kTRUE ); 107 void SetExpFitOK( const Bool_t b=kTRUE ); 108 void SetFourierSpectrumOK( const Bool_t b=kTRUE ); 109 void SetExcluded ( const Bool_t b=kTRUE ); 94 110 95 111 // Getters … … 101 117 const Double_t GetProb() const { return fProb; } 102 118 const Int_t GetNdf() const; 119 const Double_t GetPickup() const; 120 const Int_t GetPixId() const { return fPixId; } 103 121 104 122 const Double_t GetSlope() const; … … 136 154 const Bool_t IsEmpty() const; 137 155 const Bool_t IsFourierSpectrumOK() const; 156 const Bool_t IsExcluded () const; 138 157 139 158 // Fill … … 145 164 Bool_t FitGaus(Option_t *option="RQ0", 146 165 const Double_t xmin=0., const Double_t xmax=0.); // Fit the histogram HGausHist with a Gaussian 147 void BypassFit(); // Take mean and RMS from the histogram148 166 Bool_t RepeatFit(const Option_t *option="RQ0"); // Repeat fit within limits defined by fPickupLimit 167 void BypassFit(); // Take mean and RMS from the histogram 149 168 150 169 // Draws … … 155 174 156 175 // Miscelleaneous 176 virtual void ChangeHistId(Int_t id); // Changes names and titles of the histogram 157 177 void CreateFourierSpectrum(); // Create the fourier spectrum out of fEvents 158 178 void CreateGraphEvents(); // Create the TGraph fGraphEvents of fEvents
Note:
See TracChangeset
for help on using the changeset viewer.