Changeset 5854 for trunk/MagicSoft
- Timestamp:
- 01/16/05 12:48:33 (20 years ago)
- Location:
- trunk/MagicSoft/Mars
- Files:
-
- 3 edited
-
Changelog (modified) (1 diff)
-
mcalib/MCalibrateData.cc (modified) (11 diffs)
-
mcalib/MCalibrateData.h (modified) (5 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/MagicSoft/Mars/Changelog
r5853 r5854 31 31 ( SetSignalType(MCalibrateData::kPhe) ). 32 32 Default is old version in photons 33 - speed up Process by storing pre-calculated calibration constants 34 in arrays (needed 40% of CPU time of the eventloop before, now: ) 33 35 34 36 -
trunk/MagicSoft/Mars/mcalib/MCalibrateData.cc
r5853 r5854 280 280 } 281 281 282 fCalibConsts.Reset(); 283 fCalibFFactors.Reset(); 284 282 285 return kTRUE; 283 286 } … … 378 381 } 379 382 383 const Int_t npixels = fGeomCam->GetNumPixels(); 384 385 if (fCalibrationMode > kNone) 386 { 387 388 if (fCalibrations->GetSize() != npixels) 389 { 390 *fLog << err << GetDescriptor() 391 << ": Size mismatch between MGeomCam and MCalibrationChargeCam ... abort!" << endl; 392 return kFALSE; 393 } 394 395 if (fBadPixels->GetSize() != npixels) 396 { 397 *fLog << err << GetDescriptor() 398 << ": Size mismatch between MGeomCam and MBadPixelsCam ... abort!" << endl; 399 return kFALSE; 400 } 401 402 if (fBadPixels->GetSize() != npixels) 403 { 404 *fLog << err << GetDescriptor() 405 << ": Size mismatch between MGeomCam and MBadPixelsCam ... abort!" << endl; 406 return kFALSE; 407 } 408 } 409 410 fCalibConsts .Set(npixels); 411 fCalibFFactors.Set(npixels); 412 fHiLoConv .Set(npixels); 413 fHiLoConvErr .Set(npixels); 414 415 if (!UpdateConversionFactors()) 416 return kFALSE; 417 380 418 if (TestPedestalFlag(kRun)) 381 419 Calibrate(kFALSE, kTRUE); … … 386 424 // -------------------------------------------------------------------------- 387 425 // 388 // Get conversion factor and its error from MCalibrationCam 389 // 390 Bool_t MCalibrateData::GetConversionFactor(UInt_t pixidx, Float_t &hiloconv, Float_t &hiloconverr, 391 Float_t &calibConv, Float_t &calibConvVar, 392 Float_t &calibFFactor) const 426 // Update the conversion factors and F-Factors from MCalibrationCams into 427 // the arrays. Later, the here pre-calcualted conversion factors get simply 428 // copied from memory. 429 // 430 // This function can be called from outside in case that the MCalibrationCams 431 // have been updated... 432 // 433 Bool_t MCalibrateData::UpdateConversionFactors() 393 434 { 394 435 // … … 397 438 const Float_t zenith = -1.; 398 439 399 hiloconv = 1.; 400 hiloconverr = 0.; 401 calibConv = 1.; 402 calibConvVar = 0.; 403 calibFFactor = 0.; 404 405 Float_t calibQE = 1.; 406 Float_t calibQEVar = 0.; 407 408 if(fCalibrationMode!=kNone) 409 { 410 MCalibrationChargePix &pix = (MCalibrationChargePix&)(*fCalibrations)[pixidx]; 411 412 hiloconv = pix.GetConversionHiLo (); 413 hiloconverr= pix.GetConversionHiLoErr(); 414 415 if ((*fBadPixels)[pixidx].IsUnsuitable()) 416 return kFALSE; 417 418 calibConv = pix.GetMeanConvFADC2Phe(); 419 calibConvVar = pix.GetMeanConvFADC2PheVar(); 420 calibFFactor = pix.GetMeanFFactorFADC2Phot(); 421 422 MCalibrationQEPix &qe = (MCalibrationQEPix&) (*fQEs)[pixidx]; 423 424 switch(fCalibrationMode) 440 UInt_t skip = 0; 441 442 for (UInt_t pixidx=0; pixidx<fGeomCam->GetNumPixels(); pixidx++) 443 { 444 445 Float_t hiloconv = 1.; 446 Float_t hiloconverr = 0.; 447 Float_t calibConv = 1.; 448 Float_t calibConvVar = 0.; 449 Float_t calibFFactor = 0.; 450 451 Float_t calibQE = 1.; 452 Float_t calibQEVar = 0.; 453 454 if(fCalibrationMode!=kNone) 425 455 { 426 case kFlatCharge: 427 { 428 MCalibrationChargePix &avpix = (MCalibrationChargePix&)fCalibrations->GetAverageArea(0); 429 calibConv = avpix.GetMean() / (pix.GetMean() * fGeomCam->GetPixRatio(pixidx)); 430 calibConvVar = (avpix.GetMeanRelVar() + pix.GetMeanRelVar()) * calibConv * calibConv; 431 if (pix.IsFFactorMethodValid()) 456 MCalibrationChargePix &pix = (MCalibrationChargePix&)(*fCalibrations)[pixidx]; 457 458 hiloconv = pix.GetConversionHiLo (); 459 hiloconverr= pix.GetConversionHiLoErr(); 460 461 if ((*fBadPixels)[pixidx].IsUnsuitable()) 462 { 463 skip++; 464 continue; 465 } 466 467 calibConv = pix.GetMeanConvFADC2Phe(); 468 calibConvVar = pix.GetMeanConvFADC2PheVar(); 469 calibFFactor = pix.GetMeanFFactorFADC2Phot(); 470 471 MCalibrationQEPix &qe = (MCalibrationQEPix&) (*fQEs)[pixidx]; 472 473 switch(fCalibrationMode) 474 { 475 case kFlatCharge: 432 476 { 433 const Float_t convmin1 = qe.GetQECascadesFFactor(zenith)/pix.GetMeanConvFADC2Phe(); 434 if (convmin1 > 0) 435 calibFFactor *= TMath::Sqrt(convmin1); 436 else 437 calibFFactor = -1.; 477 MCalibrationChargePix &avpix = (MCalibrationChargePix&)fCalibrations->GetAverageArea(0); 478 calibConv = avpix.GetMean() / (pix.GetMean() * fGeomCam->GetPixRatio(pixidx)); 479 calibConvVar = (avpix.GetMeanRelVar() + pix.GetMeanRelVar()) * calibConv * calibConv; 480 if (pix.IsFFactorMethodValid()) 481 { 482 const Float_t convmin1 = qe.GetQECascadesFFactor(zenith)/pix.GetMeanConvFADC2Phe(); 483 if (convmin1 > 0) 484 calibFFactor *= TMath::Sqrt(convmin1); 485 else 486 calibFFactor = -1.; 487 } 488 break; 438 489 } 439 break; 440 } 441 case kBlindPixel: 442 if (!qe.IsBlindPixelMethodValid()) 443 return kFALSE; 444 calibQE = qe.GetQECascadesBlindPixel ( zenith ); 445 calibQEVar = qe.GetQECascadesBlindPixelVar( zenith ); 446 break; 447 448 case kPinDiode: 449 if (!qe.IsPINDiodeMethodValid()) 450 return kFALSE; 451 calibQE = qe.GetQECascadesPINDiode ( zenith ); 452 calibQEVar = qe.GetQECascadesPINDiodeVar( zenith ); 453 break; 454 455 case kFfactor: 456 if (!pix.IsFFactorMethodValid()) 457 return kFALSE; 458 calibQE = qe.GetQECascadesFFactor ( zenith ); 459 calibQEVar = qe.GetQECascadesFFactorVar( zenith ); 460 break; 461 462 case kCombined: 463 if (!qe.IsCombinedMethodValid()) 464 return kFALSE; 465 calibQE = qe.GetQECascadesCombined ( zenith ); 466 calibQEVar = qe.GetQECascadesCombinedVar( zenith ); 467 break; 468 469 case kDummy: 470 hiloconv = 1.; 471 hiloconverr = 0.; 472 break; 473 } /* switch calibration mode */ 474 } /* if(fCalibrationMode!=kNone) */ 475 else 476 { 477 calibConv = 1./fGeomCam->GetPixRatio(pixidx); 478 } 479 480 calibConv /= calibQE; 481 482 if (calibConv != 0. && calibQE != 0.) 483 { 484 // Now doing: 485 // calibConvVar = calibConvVar/(calibConv*calibConv) + calibQEVar/(calibQE*calibQE); 486 // calibConvVar *= (calibConv*calibConv); 487 calibConvVar += calibQEVar*(calibConv*calibConv)/(calibQE*calibQE); 490 case kBlindPixel: 491 if (!qe.IsBlindPixelMethodValid()) 492 { 493 skip++; 494 continue; 495 } 496 calibQE = qe.GetQECascadesBlindPixel ( zenith ); 497 calibQEVar = qe.GetQECascadesBlindPixelVar( zenith ); 498 break; 499 500 case kPinDiode: 501 if (!qe.IsPINDiodeMethodValid()) 502 { 503 skip++; 504 continue; 505 } 506 calibQE = qe.GetQECascadesPINDiode ( zenith ); 507 calibQEVar = qe.GetQECascadesPINDiodeVar( zenith ); 508 break; 509 510 case kFfactor: 511 if (!pix.IsFFactorMethodValid()) 512 { 513 skip++; 514 continue; 515 } 516 calibQE = qe.GetQECascadesFFactor ( zenith ); 517 calibQEVar = qe.GetQECascadesFFactorVar( zenith ); 518 break; 519 520 case kCombined: 521 if (!qe.IsCombinedMethodValid()) 522 { 523 skip++; 524 continue; 525 } 526 calibQE = qe.GetQECascadesCombined ( zenith ); 527 calibQEVar = qe.GetQECascadesCombinedVar( zenith ); 528 break; 529 530 case kDummy: 531 hiloconv = 1.; 532 hiloconverr = 0.; 533 break; 534 } /* switch calibration mode */ 535 } /* if(fCalibrationMode!=kNone) */ 536 else 537 { 538 calibConv = 1./fGeomCam->GetPixRatio(pixidx); 539 } 540 541 calibConv /= calibQE; 542 543 if (calibConv != 0. && calibQE != 0.) 544 { 545 // Now doing: 546 calibConvVar = calibConvVar/(calibConv*calibConv) + calibQEVar/(calibQE*calibQE); 547 calibConvVar *= (calibConv*calibConv); 548 // The above two lines had been commented by TB and replaced by the following line 549 // (without notice to me!) but it is simply wrong! 550 // calibConvVar += calibQEVar*(calibConv*calibConv)/(calibQE*calibQE); 551 } 552 553 calibConv *= fRenormFactor; 554 555 fHiLoConv [pixidx] = hiloconv; 556 fHiLoConvErr [pixidx] = hiloconverr; 557 fCalibConsts [pixidx] = calibConv; 558 fCalibFFactors[pixidx] = calibFFactor; 559 560 } /* for (Int_t pixidx=0; pixidx<fGeomCam->GetNumPixels(); pixidx++) */ 561 562 if (skip>fGeomCam->GetNumPixels()*0.9) 563 { 564 *fLog << warn << GetDescriptor() 565 << ": WARNING - GetConversionFactor has skipped more than 90% of the pixels... abort." << endl; 566 return kFALSE; 488 567 } 489 568 … … 491 570 } 492 571 572 573 // -------------------------------------------------------------------------- 574 // 575 // Apply the conversion factors and F-Factors from the arrays to the data. 576 // 577 // The flags 'data' and 'pedestal' decide whether the signal and/or the pedetals 578 // shall be calibrated, respectively. 579 // 493 580 Int_t MCalibrateData::Calibrate(Bool_t data, Bool_t pedestal) const 494 581 { … … 500 587 const Float_t sqrtslices = TMath::Sqrt(slices); 501 588 502 Float_t hiloconv;503 Float_t hiloconverr;504 Float_t calibConv;505 Float_t calibConvErr;506 Float_t calibFFactor;507 508 UInt_t skip = 0;509 589 for (UInt_t pixidx=0; pixidx<npix; pixidx++) 510 590 { 511 if (!GetConversionFactor(pixidx, hiloconv, hiloconverr, 512 calibConv, calibConvErr, calibFFactor)) 513 { 514 skip++; 515 continue; 516 } 517 518 calibConv *= fRenormFactor; 519 520 if (data) 591 592 if (data) 521 593 { 522 594 const MExtractedSignalPix &sig = (*fSignals)[pixidx]; … … 527 599 if (sig.IsLoGainUsed()) 528 600 { 529 signal = sig.GetExtractedSignalLoGain()* hiloconv;530 signalErr = sig nal*hiloconverr;601 signal = sig.GetExtractedSignalLoGain()*fHiLoConv [pixidx]; 602 signalErr = sig.GetExtractedSignalLoGain()*fHiLoConvErr[pixidx]; 531 603 } 532 604 else … … 536 608 } 537 609 538 const Float_t nphot = signal*calibConv; 539 const Float_t nphotErr = calibFFactor*TMath::Sqrt(TMath::Abs(nphot)); 540 541 // 542 // The following part is the commented first version of the error calculation 543 // Contact Markus Gaug for questions (or wait for the next documentation update...) 544 // 545 /* 546 nphotErr = signal > 0 ? signalErr*signalErr / (signal * signal) : 0. 547 + calibConv > 0 ? calibConvVar / (calibConv * calibConv ) : 0.; 548 nphotErr = TMath::Sqrt(nphotErr) * nphot; 549 */ 610 const Float_t nphot = signal * fCalibConsts [pixidx]; 611 const Float_t nphotErr = TMath::Sqrt(TMath::Abs(nphot)) * fCalibFFactors[pixidx]; 550 612 551 613 MCerPhotPix *cpix = fCerPhotEvt->AddPixel(pixidx, nphot, nphotErr); … … 555 617 556 618 if (sig.GetNumLoGainSaturated() > 0) 557 cpix->SetPixelSaturated(); 558 } 619 cpix->SetPixelSaturated(); 620 } /* if (data) */ 621 559 622 560 623 if (pedestal) 561 624 { 562 /*563 // pedestals/(used FADC slices) in [ADC] counts564 const Float_t pedes = (*fPedestalMean)[pixidx].GetPedestal() * slices;565 const Float_t pedrms = (*fPedestalRms)[pixidx].GetPedestalRms() * sqrtslices;566 567 //568 // pedestals/(used FADC slices) in [number of photons]569 //570 const Float_t pedphot = pedes * calibConv;571 const Float_t pedphotrms = pedrms * calibConv;572 573 (*fPedPhot)[pixidx].Set(pedphot, pedphotrms);574 */575 625 TIter NextPed(&fPedestalCams); 576 626 TIter NextPhot(&fPedPhotCams); … … 578 628 MPedestalCam *pedestal = 0; 579 629 MPedPhotCam *pedphot = 0; 630 631 const Float_t pedmeancalib = slices *fCalibConsts[pixidx]; 632 const Float_t pedrmscalib = sqrtslices*fCalibConsts[pixidx]; 580 633 581 634 while ((pedestal=(MPedestalCam*)NextPed()) && 582 635 (pedphot =(MPedPhotCam*)NextPhot())) 583 636 { 584 // pedestals/(used FADC slices) in [ADC] counts 585 const Float_t pedes = (*pedestal)[pixidx].GetPedestal() * slices; 586 const Float_t pedrms = (*pedestal)[pixidx].GetPedestalRms() * sqrtslices; 587 588 // pedestals/(used FADC slices) in [number of photons] 589 const Float_t mean = pedes * calibConv; 590 const Float_t rms = pedrms * calibConv; 591 592 (*pedphot)[pixidx].Set(mean, rms); 593 pedphot->SetReadyToSave(); 637 // pedestals/(used FADC slices) in [number of photons] 638 const Float_t mean = (*pedestal)[pixidx].GetPedestal() *pedmeancalib; 639 const Float_t rms = (*pedestal)[pixidx].GetPedestalRms()*pedrmscalib; 640 641 (*pedphot)[pixidx].Set(mean, rms); 642 pedphot->SetReadyToSave(); 594 643 } 595 } 596 } 597 598 if (skip>npix*0.9) 599 { 600 *fLog << warn << "WARNING - GetConversionFactor has skipped more than 90% of the pixels... skip." << endl; 601 return kCONTINUE; 644 } /* if (pedestal) */ 602 645 } 603 646 … … 617 660 Int_t MCalibrateData::Process() 618 661 { 619 /*620 if (fCalibrations->GetNumPixels() != (UInt_t)fSignals->GetSize())621 {622 // FIXME: MExtractedSignal must be of variable size -623 // like MCerPhotEvt - because we must be able624 // to reduce size by zero supression625 // For the moment this check could be done in ReInit...626 *fLog << err << "MExtractedSignal and MCalibrationCam have different sizes... abort." << endl;627 return kFALSE;628 }629 */630 631 662 return Calibrate(fCalibrationMode!=kSkip, TestPedestalFlag(kEvent)); 632 663 } -
trunk/MagicSoft/Mars/mcalib/MCalibrateData.h
r5853 r5854 21 21 #endif 22 22 23 #ifndef MARS_MArrayF 24 #include "MArrayF.h" 25 #endif 26 23 27 class MGeomCam; 24 28 class MBadPixelsCam; … … 34 38 { 35 39 private: 40 36 41 MGeomCam *fGeomCam; //! Camera geometry container 37 42 MBadPixelsCam *fBadPixels; //! Bad Pixels information … … 51 56 TList fPedPhotCams; //! List of pointers to corresponding output MPedPhotCam 52 57 53 Int_t Calibrate(Bool_t data, Bool_t pedestal) const; 58 MArrayF fCalibConsts; //! Array of calibration constants for each pixel, calculated only once! 59 MArrayF fCalibFFactors; //! Array of calibration F-Factors for each pixel, calculated only once! 60 MArrayF fHiLoConv; //! Array of calibration constants for each pixel, calculated only once! 61 MArrayF fHiLoConvErr; //! Array of calibration F-Factors for each pixel, calculated only once! 54 62 55 Bool_t GetConversionFactor(UInt_t, Float_t &, Float_t &, Float_t &, Float_t &, Float_t &) const;63 Int_t Calibrate(Bool_t data, Bool_t pedestal) const; 56 64 57 65 Int_t PreProcess(MParList *pList); … … 93 101 const char *name=NULL, const char *title=NULL); 94 102 103 void AddPedestal(const char *name="Cam"); 104 void AddPedestal(const char *pedestal, const char *pedphot); 105 95 106 void EnablePedestalType(PedestalType_t i) { fPedestalFlag |= i; } 96 107 void SetPedestalFlag(PedestalType_t i=kRun) { fPedestalFlag = i; } … … 100 111 void SetSignalType ( SignalType_t sigtype=kPhot ) { fSignalType =sigtype; } 101 112 102 void AddPedestal(const char *name="Cam"); 103 void AddPedestal(const char *pedestal, const char *pedphot); 104 113 Bool_t UpdateConversionFactors(); 114 105 115 ClassDef(MCalibrateData, 1) // Task to calibrate FADC counts into Cherenkov photons 106 116 };
Note:
See TracChangeset
for help on using the changeset viewer.
