- Timestamp:
- Oct 5, 2011, 3:13:36 PM (13 years ago)
- Location:
- BAORadio/AmasNancay/trunk
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
BAORadio/AmasNancay/trunk/mergeRawOnOff.cc
r551 r552 55 55 //-------------------------------------------------------------- 56 56 //Utility functions 57 58 class median_of_empty_list_exception:public std::exception{ 59 virtual const char* what() const throw() { 60 return "Attempt to take the median of an empty list of numbers. " 61 "The median of an empty list is undefined."; 62 } 63 }; 64 template<class RandAccessIter> 65 double median(RandAccessIter begin, RandAccessIter end) 66 throw(median_of_empty_list_exception){ 67 if(begin == end){ throw median_of_empty_list_exception(); } 68 std::size_t size = end - begin; 69 std::size_t middleIdx = size/2; 70 RandAccessIter target = begin + middleIdx; 71 std::nth_element(begin, target, end); 72 73 if(size % 2 != 0){ //Odd number of elements 74 return *target; 75 }else{ //Even number of elements 76 double a = *target; 77 RandAccessIter targetNeighbor= target-1; 78 std::nth_element(begin, targetNeighbor, end); 79 return (a+*targetNeighbor)/2.0; 80 } 81 } 82 57 83 //-------------------------------------------------------------- 58 84 char *regexp (const char *string, const char *patrn, int *begin, int *end) { … … 291 317 292 318 if ((sa_size_t)nSpectra < TOTAL_NUM_CYCLES ) { 293 //STORE THE spectra294 // tableOfSpectra(Range::all(),Range::all(),Range(nSpectra)).Show();295 // aSpec.Show();296 319 tableOfSpectra(Range::all(),Range::all(),Range(nSpectra)) = aSpec; 297 // aSpec.Show();298 // cout << "deb 1" << endl;299 // tmp = aSpec;300 // cout << "deb 2" << endl;301 320 } else { 302 321 stringstream tmp; … … 310 329 //Resize 311 330 312 nSpectra--;//For further selection in Range313 if (para.debuglev_>90){314 cout << "Cycles < 0," << nSpectra << ">" << endl;315 }316 331 317 318 332 TMatrix<r_4> medianOfSpectra(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ); 319 333 //Compute the median for each freq. and Channel 320 334 for (sa_size_t iCh=0; iCh<NUMBER_OF_CHANNELS; iCh++){ 321 335 for (sa_size_t freq =0; freq<NUMBER_OF_FREQ; freq++){ 322 TVector<r_4> tmp0(tableOfSpectra(Range(freq),Range(iCh),Range(0,nSpectra )).CompactAllDimensions());336 TVector<r_4> tmp0(tableOfSpectra(Range(freq),Range(iCh),Range(0,nSpectra-1)).CompactAllDimensions()); 323 337 vector<r_4> tmp; 324 338 tmp0.FillTo(tmp); 325 nth_element(tmp.begin(),tmp.begin()+tmp.size()/2,tmp.end()); 326 medianOfSpectra(iCh,freq) = *(tmp.begin()+tmp.size()/2); 339 medianOfSpectra(iCh,freq) = median(tmp.begin(),tmp.end()); 327 340 } 328 341 } … … 338 351 tmp << nSpectra; 339 352 string fileName = para.outPath_+"/medianDiffOnOffRaw_"+StringToLower(para.sourceName_)+"-"+tmp.str()+"Cycles.ppf"; 340 cout << "Save median based on " << (nSpectra+1)<< " cycles " << endl;353 cout << "Save median based on " << nSpectra << " cycles " << endl; 341 354 POutPersist fos(fileName); 342 355 … … 347 360 tag = "sigmamedred"; 348 361 fos << PPFNameTag(tag) << sigmaRedMtx; 362 tag = "cycleVsfreq"; 363 364 TArray<r_4> tarr(tableOfSpectra(Range::all(),Range::all(),Range(0,nSpectra-1))); 365 fos << PPFNameTag(tag) << tarr; 349 366 } 350 367 } -
BAORadio/AmasNancay/trunk/plotDiffOnOff.pic
r550 r552 6 6 set ncycles $2 7 7 set nFreqBin $3 8 set action $4 9 10 set defatt "font=helvetica,bold,20" 11 12 13 14 if ( ${action} == "mean" ) then 15 8 16 openppf meanDiffOnOffRaw_${source}-${ncycles}Cycles.ppf 9 17 … … 13 21 objaoper mean row 0 s0 14 22 objaoper mean row 1 s1 15 graphicatt " xylimits=1250,1500,-0.01,0.01"23 graphicatt "${defatt} xylimits=1250,1500,-0.01,0.01" 16 24 newwin 1 1 17 25 plot2d s0 (n/8192)*250+1250 val n>0 "cpts blue nsta notit" 18 26 plot2d s1 (n/8192)*250+1250 val n>0 "same cpts red nsta notit" 19 settitle "<(ON-OFF)[r,c]/Gain[r]>_{r,c} $source (${ncycles} cycles) Ch 0 (blue) Ch 1 (red)" 27 settitle "$source (${ncycles} cycles) Ch 0 (blue) Ch 1 (red)" ' ' $defatt 28 setaxelabels "Freq. (MHz)" "<(ON-OFF)[r,c]/Gain[r]>_{r,c}" $defatt 29 20 30 21 31 … … 27 37 plot2d rs0 (n/${nFreqBin})*250+1250 val n>0 "marker=fcircle,5 blue nsta notit" 28 38 plot2d rs1 (n/${nFreqBin})*250+1250 val n>0 "same marker=fcircle,5 red nsta notit" 29 settitle "<(ON-OFF)[r,c]/Gain[r]>_{r,c} Rebin(${nFreqBin}) $source (${ncycles} cycles) Ch 0 (blue) Ch 1 (red)" 39 settitle "$source (${ncycles} cycles) Ch 0 (blue) Ch 1 (red)" ' ' $defatt 40 setaxelabels "Freq. (MHz)" "<(ON-OFF)[r,c]/Gain[r]>_{r,c} Rebin(${nFreqBin})" $defatt 41 30 42 31 43 ##Reduced Mean … … 36 48 plot2d rss0 (n/${nFreqBin})*250+1250 val n>0 "marker=fcircle,5 blue nsta notit" 37 49 plot2d rss1 (n/${nFreqBin})*250+1250 val n>0 "same marker=fcircle,5 red nsta notit" 38 settitle "Sigma[(ON-OFF)[r,c]/Gain[r]] Rebin(${nFreqBin}) $source (${ncycles} cycles) Ch 0 (blue) Ch 1 (red)" 50 settitle "$source (${ncycles} cycles) Ch 0 (blue) Ch 1 (red)" ' ' $defatt 51 setaxelabels "Freq. (MHz)" "Sigma[(ON-OFF)[r,c]/Gain[r]] Rebin(${nFreqBin})" $defatt 52 53 else 54 55 openppf medianDiffOnOffRaw_${source}-${ncycles}Cycles.ppf 39 56 40 57 58 ##General Median 59 del s0 s1 60 objaoper median row 0 s0 61 objaoper median row 1 s1 62 graphicatt "xylimits=1250,1500,-0.01,0.01" 63 newwin 1 1 64 plot2d s0 (n/8192)*250+1250 val n>0 "cpts blue nsta notit" 65 plot2d s1 (n/8192)*250+1250 val n>0 "same cpts red nsta notit" 66 settitle "$source (${ncycles} cycles) Ch 0 (blue) Ch 1 (red)" ' ' $defatt 67 setaxelabels "Freq. (MHz)" "Median[(ON-OFF)[r,c]/Gain[r]]" $defatt 68 69 70 ##Reduced Mean of the median 71 newwin 1 1 72 del rs0 rs1 73 objaoper meanmedred row 0 rs0 74 objaoper meanmedred row 1 rs1 75 graphicatt "xylimits=1250,1500,-0.01,0.01" 76 plot2d rs0 (n/${nFreqBin})*250+1250 val n>0 "marker=fcircle,5 blue nsta notit" 77 plot2d rs1 (n/${nFreqBin})*250+1250 val n>0 "same marker=fcircle,5 red nsta notit" 78 settitle "$source (${ncycles} cycles) Ch 0 (blue) Ch 1 (red)" ' ' $defatt 79 setaxelabels "Freq. (MHz)" "Mean Median[(ON-OFF)[r,c]/Gain[r]] Rebin(${nFreqBin})" $defatt 80 81 82 ##Reduced Sigma of the median 83 newwin 1 1 84 del rss0 rss1 85 objaoper sigmamedred row 0 rss0 86 objaoper sigmamedred row 1 rss1 87 graphicatt "xylimits=1250,1500,0,0.01" 88 plot2d rss0 (n/${nFreqBin})*250+1250 val n>0 "marker=fcircle,5 blue nsta notit" 89 plot2d rss1 (n/${nFreqBin})*250+1250 val n>0 "same marker=fcircle,5 red nsta notit" 90 settitle "$source (${ncycles} cycles) Ch 0 (blue) Ch 1 (red)" ' ' $defatt 91 setaxelabels "Freq. (MHz)" "Sigma Median[(ON-OFF)[r,c]/Gain[r]] Rebin(${nFreqBin})" $defatt 92 93 94 endif 95
Note: See TracChangeset
for help on using the changeset viewer.