Changeset 5162
- Timestamp:
- 10/01/04 21:24:31 (20 years ago)
- Location:
- trunk/MagicSoft/Mars/msignal
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/MagicSoft/Mars/msignal/MExtractTimeAndChargeDigitalFilter.cc
r5155 r5162 56 56 #include "MPedestalPix.h" 57 57 58 #include "TFile.h" // remove!! 59 #include "TH1F.h" // remove!! 58 #include "TFile.h" 59 #include "TH1F.h" 60 #include "TH2F.h" 60 61 #include "TString.h" 62 #include "TMatrix.h" 61 63 62 64 #include <fstream> … … 74 76 const Int_t MExtractTimeAndChargeDigitalFilter::fgBinningResolutionHiGain = 10; 75 77 const Int_t MExtractTimeAndChargeDigitalFilter::fgBinningResolutionLoGain = 10; 78 const Int_t MExtractTimeAndChargeDigitalFilter::fgSignalStartBinHiGain = 4; 79 const Int_t MExtractTimeAndChargeDigitalFilter::fgSignalStartBinLoGain = 4; 76 80 // -------------------------------------------------------------------------- 77 81 // … … 94 98 SetWindowSize(); 95 99 SetBinningResolution(); 100 SetSignalStartBin(); 96 101 97 102 ReadWeightsFile(""); … … 165 170 } 166 171 167 //----------------------------------------------------------------------------168 //169 // Read a pre-defined weights file into the class.170 // This is mandatory for the extraction171 //172 // If filenname is empty, then all weights will be set to 1.173 //174 Bool_t MExtractTimeAndChargeDigitalFilter::ReadWeightsFile(TString filename)175 {176 177 fAmpWeightsHiGain .Set(fBinningResolutionHiGain*fWindowSizeHiGain);178 fAmpWeightsLoGain .Set(fBinningResolutionLoGain*fWindowSizeLoGain);179 fTimeWeightsHiGain.Set(fBinningResolutionHiGain*fWindowSizeHiGain);180 fTimeWeightsLoGain.Set(fBinningResolutionLoGain*fWindowSizeLoGain);181 182 if (filename.IsNull())183 {184 for (UInt_t i=0; i<fAmpWeightsHiGain.GetSize(); i++)185 {186 fAmpWeightsHiGain [i] = 1.;187 fTimeWeightsHiGain[i] = 1.;188 }189 for (UInt_t i=0; i<fAmpWeightsLoGain.GetSize(); i++)190 {191 fAmpWeightsLoGain [i] = 1.;192 fTimeWeightsLoGain[i] = 1.;193 }194 return kTRUE;195 }196 197 ifstream fin(filename.Data());198 199 if (!fin)200 {201 *fLog << err << GetDescriptor()202 << ": No weights file found: " << filename << endl;203 return kFALSE;204 }205 206 Int_t len = 0;207 Int_t cnt = 0;208 Bool_t hi = kFALSE;209 Bool_t lo = kFALSE;210 211 TString str;212 213 while (1)214 {215 216 str.ReadLine(fin);217 if (!fin)218 break;219 220 221 if (str.Contains("# High Gain Weights:"))222 {223 str.ReplaceAll("# High Gain Weights:","");224 sscanf(str.Data(),"%2i%2i",&fWindowSizeHiGain,&fBinningResolutionHiGain);225 *fLog << inf << "Found number of High Gain slices: " << fWindowSizeHiGain226 << " and High Gain resolution: " << fBinningResolutionHiGain << endl;227 len = fBinningResolutionHiGain*fWindowSizeHiGain;228 fAmpWeightsHiGain .Set(len);229 fTimeWeightsHiGain.Set(len);230 hi = kTRUE;231 continue;232 }233 234 if (str.Contains("# Low Gain Weights:"))235 {236 str.ReplaceAll("# Low Gain Weights:","");237 sscanf(str.Data(),"%2i%2i",&fWindowSizeLoGain,&fBinningResolutionLoGain);238 *fLog << inf << "Found number of Low Gain slices: " << fWindowSizeLoGain239 << " and Low Gain resolution: " << fBinningResolutionLoGain << endl;240 len = fBinningResolutionLoGain*fWindowSizeHiGain;241 fAmpWeightsLoGain .Set(len);242 fTimeWeightsLoGain.Set(len);243 lo = kTRUE;244 continue;245 }246 247 if (str.Contains("#"))248 continue;249 250 if (len == 0)251 continue;252 253 sscanf(str.Data(),"\t%f\t%f",lo ? &fAmpWeightsLoGain [cnt] : &fAmpWeightsHiGain [cnt],254 lo ? &fTimeWeightsLoGain[cnt] : &fTimeWeightsHiGain[cnt]);255 256 if (++cnt == len)257 {258 len = 0;259 cnt = 0;260 }261 }262 263 if (cnt != len)264 {265 *fLog << err << GetDescriptor()266 << ": Size mismatch in weights file " << filename << endl;267 return kFALSE;268 }269 270 if (!hi)271 {272 *fLog << err << GetDescriptor()273 << ": No correct header found in weights file " << filename << endl;274 return kFALSE;275 }276 277 return kTRUE;278 }279 280 //----------------------------------------------------------------------------281 //282 // Read a pre-defined weights file into the class.283 // This is mandatory for the extraction284 //285 Bool_t MExtractTimeAndChargeDigitalFilter::WriteWeightsFile(TString filename)286 {287 288 289 Float_t amp[60];290 Float_t tim[60];291 292 ofstream fn(filename.Data());293 294 fn << "# High Gain Weights: 6 10" << endl;295 fn << "# (Amplitude) (Time) " << endl;296 297 for (UInt_t i=0; i<60; i++)298 {299 fn << "\t" << amp[i] << "\t" << tim[i] << endl;300 }301 302 fn << "# Low Gain Weights: 6 10" << endl;303 fn << "# (Amplitude) (Time) " << endl;304 305 for (UInt_t i=0; i<60; i++)306 {307 fn << "\t" << amp[i] << "\t" << tim[i] << endl;308 }309 return kTRUE;310 }311 312 172 // -------------------------------------------------------------------------- 313 173 // … … 337 197 return kTRUE; 338 198 } 199 339 200 340 201 void MExtractTimeAndChargeDigitalFilter::FindTimeAndChargeHiGain(Byte_t *ptr, Byte_t *logain, Float_t &sum, Float_t &dsum, … … 653 514 return kTRUE; 654 515 } 516 517 //---------------------------------------------------------------------------- 518 // 519 // Read a pre-defined weights file into the class. 520 // This is mandatory for the extraction 521 // 522 // If filenname is empty, then all weights will be set to 1. 523 // 524 Bool_t MExtractTimeAndChargeDigitalFilter::ReadWeightsFile(TString filename) 525 { 526 527 fAmpWeightsHiGain .Set(fBinningResolutionHiGain*fWindowSizeHiGain); 528 fAmpWeightsLoGain .Set(fBinningResolutionLoGain*fWindowSizeLoGain); 529 fTimeWeightsHiGain.Set(fBinningResolutionHiGain*fWindowSizeHiGain); 530 fTimeWeightsLoGain.Set(fBinningResolutionLoGain*fWindowSizeLoGain); 531 532 if (filename.IsNull()) 533 { 534 for (UInt_t i=0; i<fAmpWeightsHiGain.GetSize(); i++) 535 { 536 fAmpWeightsHiGain [i] = 1.; 537 fTimeWeightsHiGain[i] = 1.; 538 } 539 for (UInt_t i=0; i<fAmpWeightsLoGain.GetSize(); i++) 540 { 541 fAmpWeightsLoGain [i] = 1.; 542 fTimeWeightsLoGain[i] = 1.; 543 } 544 return kTRUE; 545 } 546 547 ifstream fin(filename.Data()); 548 549 if (!fin) 550 { 551 *fLog << err << GetDescriptor() 552 << ": No weights file found: " << filename << endl; 553 return kFALSE; 554 } 555 556 Int_t len = 0; 557 Int_t cnt = 0; 558 Bool_t hi = kFALSE; 559 Bool_t lo = kFALSE; 560 561 TString str; 562 563 while (1) 564 { 565 566 str.ReadLine(fin); 567 if (!fin) 568 break; 569 570 571 if (str.Contains("# High Gain Weights:")) 572 { 573 str.ReplaceAll("# High Gain Weights:",""); 574 sscanf(str.Data(),"%2i%2i",&fWindowSizeHiGain,&fBinningResolutionHiGain); 575 *fLog << inf << "Found number of High Gain slices: " << fWindowSizeHiGain 576 << " and High Gain resolution: " << fBinningResolutionHiGain << endl; 577 len = fBinningResolutionHiGain*fWindowSizeHiGain; 578 fAmpWeightsHiGain .Set(len); 579 fTimeWeightsHiGain.Set(len); 580 hi = kTRUE; 581 continue; 582 } 583 584 if (str.Contains("# Low Gain Weights:")) 585 { 586 str.ReplaceAll("# Low Gain Weights:",""); 587 sscanf(str.Data(),"%2i%2i",&fWindowSizeLoGain,&fBinningResolutionLoGain); 588 *fLog << inf << "Found number of Low Gain slices: " << fWindowSizeLoGain 589 << " and Low Gain resolution: " << fBinningResolutionLoGain << endl; 590 len = fBinningResolutionLoGain*fWindowSizeHiGain; 591 fAmpWeightsLoGain .Set(len); 592 fTimeWeightsLoGain.Set(len); 593 lo = kTRUE; 594 continue; 595 } 596 597 if (str.Contains("#")) 598 continue; 599 600 if (len == 0) 601 continue; 602 603 sscanf(str.Data(),"\t%f\t%f",lo ? &fAmpWeightsLoGain [cnt] : &fAmpWeightsHiGain [cnt], 604 lo ? &fTimeWeightsLoGain[cnt] : &fTimeWeightsHiGain[cnt]); 605 606 if (++cnt == len) 607 { 608 len = 0; 609 cnt = 0; 610 } 611 } 612 613 if (cnt != len) 614 { 615 *fLog << err << GetDescriptor() 616 << ": Size mismatch in weights file " << filename << endl; 617 return kFALSE; 618 } 619 620 if (!hi) 621 { 622 *fLog << err << GetDescriptor() 623 << ": No correct header found in weights file " << filename << endl; 624 return kFALSE; 625 } 626 627 return kTRUE; 628 } 629 630 //---------------------------------------------------------------------------- 631 // 632 // Create the weights file 633 // Beware that the shape-histogram has to contain the pulse starting at bin 1 634 // 635 Bool_t MExtractTimeAndChargeDigitalFilter::WriteWeightsFile(TString filename, TH1F *shapehi, TH2F *autocorrhi, 636 TH1F *shapelo, TH2F *autocorrlo ) 637 { 638 639 const Int_t nbinshi = shapehi->GetNbinsX(); 640 Float_t binwidth = shapehi->GetBinWidth(1); 641 642 TH1F *derivativehi = new TH1F(Form("%s%s",shapehi->GetName(),"_der"), 643 Form("%s%s",shapehi->GetTitle()," derivative"), 644 nbinshi, 645 shapehi->GetBinLowEdge(1), 646 shapehi->GetBinLowEdge(nbinshi)+binwidth); 647 648 // 649 // Calculate the derivative of shapehi 650 // 651 for (Int_t i = 1; i<nbinshi+1;i++) 652 { 653 derivativehi->SetBinContent(i, 654 ((shapehi->GetBinContent(i+1)-shapehi->GetBinContent(i-1))/2./binwidth)); 655 derivativehi->SetBinError(i, 656 (sqrt(shapehi->GetBinError(i+1)*shapehi->GetBinError(i+1) 657 +shapehi->GetBinError(i-1)*shapehi->GetBinError(i-1))/2./binwidth)); 658 } 659 660 // 661 // normalize the shapehi, such that the integral for fWindowSize slices is one! 662 // 663 Float_t sum = 0; 664 Int_t lasttemp = fBinningResolutionHiGain * (fSignalStartBinHiGain + fWindowSizeHiGain); 665 lasttemp = lasttemp > nbinshi ? nbinshi : lasttemp; 666 667 for (Int_t i=fBinningResolutionHiGain*fSignalStartBinHiGain; i<lasttemp; i++) { 668 sum += shapehi->GetBinContent(i); 669 } 670 sum /= fBinningResolutionHiGain; 671 672 shapehi->Scale(1./sum); 673 derivativehi->Scale(1./sum); 674 675 // 676 // read in the noise auto-correlation function: 677 // 678 TMatrix Bhi(fWindowSizeHiGain,fWindowSizeHiGain); 679 680 for (Int_t i=0; i<fWindowSizeHiGain; i++){ 681 for (Int_t j=0; j<fWindowSizeHiGain; j++){ 682 Bhi[i][j]=autocorrhi->GetBinContent(i+1+fSignalStartBinHiGain,j+1+fSignalStartBinHiGain); 683 } 684 } 685 Bhi.Invert(); 686 687 const Int_t nsizehi = fWindowSizeHiGain*fBinningResolutionHiGain; 688 fAmpWeightsHiGain.Set(nsizehi); 689 fTimeWeightsHiGain.Set(nsizehi); 690 691 // 692 // Loop over relative time in one BinningResolution interval 693 // 694 Int_t start = fBinningResolutionHiGain*fSignalStartBinHiGain + fBinningResolutionHalfHiGain; 695 696 for (Int_t i = -fBinningResolutionHalfHiGain+1; i<=fBinningResolutionHalfHiGain; i++) 697 { 698 699 TMatrix g(fWindowSizeHiGain,1); 700 TMatrix gT(1,fWindowSizeHiGain); 701 TMatrix d(fWindowSizeHiGain,1); 702 TMatrix dT(1,fWindowSizeHiGain); 703 704 for (Int_t count=0; count < fWindowSizeHiGain; count++){ 705 706 g[count][0]=shapehi->GetBinContent(start 707 +fBinningResolutionHiGain*count+i); 708 gT[0][count]=shapehi->GetBinContent(start 709 +fBinningResolutionHiGain*count+i); 710 d[count][0]=derivativehi->GetBinContent(start 711 +fBinningResolutionHiGain*count+i); 712 dT[0][count]=derivativehi->GetBinContent(start 713 +fBinningResolutionHiGain*count+i); 714 } 715 716 TMatrix m_denom = (gT*(Bhi*g))*(dT*(Bhi*d)) - (dT*(Bhi*g))*(dT*(Bhi*g)); 717 Float_t denom = m_denom[0][0]; // ROOT thinks, m_denom is still a matrix 718 719 TMatrix m_first = dT*(Bhi*d); // ROOT thinks, m_first is still a matrix 720 Float_t first = m_first[0][0]/denom; 721 722 TMatrix m_last = gT*(Bhi*d); // ROOT thinks, m_last is still a matrix 723 Float_t last = m_last[0][0]/denom; 724 725 TMatrix m1 = gT*Bhi; 726 m1 *= first; 727 728 TMatrix m2 = dT*Bhi; 729 m2 *=last; 730 731 TMatrix w_amp = m1 - m2; 732 733 TMatrix m_first1 = gT*(Bhi*g); 734 Float_t first1 = m_first1[0][0]/denom; 735 736 TMatrix m_last1 = gT*(Bhi*d); 737 Float_t last1 = m_last1 [0][0]/denom; 738 739 TMatrix m11 = dT*Bhi; 740 m11 *=first1; 741 742 TMatrix m21 = gT*Bhi; 743 m21 *=last1; 744 745 TMatrix w_time= m11 - m21; 746 747 for (Int_t count=0; count < fWindowSizeHiGain; count++) 748 { 749 const Int_t idx = i+fBinningResolutionHalfHiGain+fBinningResolutionHiGain*count-1; 750 fAmpWeightsHiGain [idx] = w_amp [0][count]; 751 fTimeWeightsHiGain[idx] = w_time[0][count]; 752 } 753 754 } // end loop over i 755 756 // 757 // Low Gain histograms 758 // 759 if (shapelo) 760 { 761 const Int_t nbinslo = shapelo->GetNbinsX(); 762 binwidth = shapelo->GetBinWidth(1); 763 764 TH1F *derivativelo = new TH1F(Form("%s%s",shapelo->GetName(),"_der"), 765 Form("%s%s",shapelo->GetTitle()," derivative"), 766 nbinslo, 767 shapelo->GetBinLowEdge(1), 768 shapelo->GetBinLowEdge(nbinslo)+binwidth); 769 770 // 771 // Calculate the derivative of shapelo 772 // 773 for (Int_t i = 1; i<nbinslo+1;i++) 774 { 775 derivativelo->SetBinContent(i, 776 ((shapelo->GetBinContent(i+1)-shapelo->GetBinContent(i-1))/2./binwidth)); 777 derivativelo->SetBinError(i, 778 (sqrt(shapelo->GetBinError(i+1)*shapelo->GetBinError(i+1) 779 +shapelo->GetBinError(i-1)*shapelo->GetBinError(i-1))/2./binwidth)); 780 } 781 782 // 783 // normalize the shapelo, such that the integral for fWindowSize slices is one! 784 // 785 sum = 0; 786 lasttemp = fBinningResolutionLoGain * (fSignalStartBinLoGain + fWindowSizeLoGain); 787 lasttemp = lasttemp > nbinslo ? nbinslo : lasttemp; 788 789 for (Int_t i=fBinningResolutionLoGain*fSignalStartBinLoGain; i<lasttemp; i++) 790 sum += shapelo->GetBinContent(i); 791 792 sum /= fBinningResolutionLoGain; 793 794 shapelo->Scale(1./sum); 795 derivativelo->Scale(1./sum); 796 797 // 798 // read in the noise auto-correlation function: 799 // 800 TMatrix Blo(fWindowSizeLoGain,fWindowSizeLoGain); 801 802 for (Int_t i=0; i<fWindowSizeLoGain; i++){ 803 for (Int_t j=0; j<fWindowSizeLoGain; j++){ 804 Blo[i][j]=autocorrlo->GetBinContent(i+1+fSignalStartBinLoGain,j+1+fSignalStartBinLoGain); 805 } 806 } 807 Blo.Invert(); 808 809 const Int_t nsizelo = fWindowSizeLoGain*fBinningResolutionLoGain; 810 fAmpWeightsLoGain.Set(nsizelo); 811 fTimeWeightsLoGain.Set(nsizelo); 812 813 // 814 // Loop over relative time in one BinningResolution interval 815 // 816 Int_t start = fBinningResolutionLoGain*fSignalStartBinLoGain + fBinningResolutionHalfLoGain; 817 818 for (Int_t i = -fBinningResolutionHalfLoGain+1; i<=fBinningResolutionHalfLoGain; i++) 819 { 820 821 TMatrix g(fWindowSizeLoGain,1); 822 TMatrix gT(1,fWindowSizeLoGain); 823 TMatrix d(fWindowSizeLoGain,1); 824 TMatrix dT(1,fWindowSizeLoGain); 825 826 for (Int_t count=0; count < fWindowSizeLoGain; count++){ 827 828 g[count][0] = shapelo->GetBinContent(start 829 +fBinningResolutionLoGain*count+i); 830 gT[0][count]= shapelo->GetBinContent(start 831 +fBinningResolutionLoGain*count+i); 832 d[count][0] = derivativelo->GetBinContent(start 833 +fBinningResolutionLoGain*count+i); 834 dT[0][count]= derivativelo->GetBinContent(start 835 +fBinningResolutionLoGain*count+i); 836 } 837 838 TMatrix m_denom = (gT*(Blo*g))*(dT*(Blo*d)) - (dT*(Blo*g))*(dT*(Blo*g)); 839 Float_t denom = m_denom[0][0]; // ROOT thinks, m_denom is still a matrix 840 841 TMatrix m_first = dT*(Blo*d); // ROOT thinks, m_first is still a matrix 842 Float_t first = m_first[0][0]/denom; 843 844 TMatrix m_last = gT*(Blo*d); // ROOT thinks, m_last is still a matrix 845 Float_t last = m_last[0][0]/denom; 846 847 TMatrix m1 = gT*Blo; 848 m1 *= first; 849 850 TMatrix m2 = dT*Blo; 851 m2 *=last; 852 853 TMatrix w_amp = m1 - m2; 854 855 TMatrix m_first1 = gT*(Blo*g); 856 Float_t first1 = m_first1[0][0]/denom; 857 858 TMatrix m_last1 = gT*(Blo*d); 859 Float_t last1 = m_last1 [0][0]/denom; 860 861 TMatrix m11 = dT*Blo; 862 m11 *=first1; 863 864 TMatrix m21 = gT*Blo; 865 m21 *=last1; 866 867 TMatrix w_time= m11 - m21; 868 869 for (Int_t count=0; count < fWindowSizeLoGain; count++) 870 { 871 const Int_t idx = i+fBinningResolutionHalfLoGain+fBinningResolutionLoGain*count-1; 872 fAmpWeightsLoGain [idx] = w_amp [0][count]; 873 fTimeWeightsLoGain[idx] = w_time[0][count]; 874 } 875 876 } // end loop over i 877 } 878 879 880 ofstream fn(filename.Data()); 881 882 fn << "# High Gain Weights: " << fWindowSizeHiGain << " " << fBinningResolutionHiGain << endl; 883 fn << "# (Amplitude) (Time) " << endl; 884 885 for (Int_t i=0; i<nsizehi; i++) 886 fn << "\t" << fAmpWeightsHiGain[i] << "\t" << fTimeWeightsHiGain[i] << endl; 887 888 fn << "# Low Gain Weights: " << fWindowSizeLoGain << " " << fBinningResolutionLoGain << endl; 889 fn << "# (Amplitude) (Time) " << endl; 890 891 for (Int_t i=0; i<nsizehi; i++) 892 fn << "\t" << fAmpWeightsLoGain[i] << "\t" << fTimeWeightsLoGain[i] << endl; 893 894 return kTRUE; 895 } -
trunk/MagicSoft/Mars/msignal/MExtractTimeAndChargeDigitalFilter.h
r5155 r5162 10 10 #endif 11 11 12 class TH1F; 13 class TH2F; 12 14 class MPedestalPix; 13 15 class MExtractTimeAndChargeDigitalFilter : public MExtractTimeAndCharge … … 23 25 static const Int_t fgBinningResolutionHiGain; 24 26 static const Int_t fgBinningResolutionLoGain; 27 static const Int_t fgSignalStartBinHiGain; 28 static const Int_t fgSignalStartBinLoGain; 25 29 26 MArrayF fHiGainSignal; 27 MArrayF fLoGainSignal; 30 MArrayF fHiGainSignal; //! Need fast access to the signals in a float way 31 MArrayF fLoGainSignal; //! Store them in separate arrays 28 32 29 33 Float_t fTimeShiftHiGain; 30 34 Float_t fTimeShiftLoGain; 31 35 36 Int_t fSignalStartBinHiGain; 37 Int_t fSignalStartBinLoGain; 38 32 39 Int_t fWindowSizeHiGain; 33 40 Int_t fWindowSizeLoGain; … … 61 68 MExtractTimeAndChargeDigitalFilter(const char *name=NULL, const char *title=NULL); 62 69 63 Bool_t WriteWeightsFile(TString filename="cosmic_weights.dat"); 70 Bool_t WriteWeightsFile(TString filename, 71 TH1F *shapehi, TH2F *autocorrhi, 72 TH1F *shapelo=NULL, TH2F *autocorrlo=NULL ); 64 73 65 74 Bool_t ReadWeightsFile(TString filename="cosmic_weights.dat"); … … 74 83 } 75 84 85 86 void SetSignalStartBin( const Int_t sh=fgSignalStartBinHiGain, const Int_t sl=fgSignalStartBinLoGain) { 87 fSignalStartBinHiGain = sh; 88 fSignalStartBinLoGain = sl; 89 } 76 90 77 91 ClassDef(MExtractTimeAndChargeDigitalFilter, 0) // Hendrik's digital filter
Note:
See TracChangeset
for help on using the changeset viewer.