Changeset 591


Ignore:
Timestamp:
Oct 26, 2011, 11:36:56 AM (13 years ago)
Author:
campagne
Message:

continue improve analysis (jec)

Location:
BAORadio/AmasNancay/trunk
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • BAORadio/AmasNancay/trunk/etude_mergeAna.pic

    r588 r591  
    1818
    1919echo $nbins1420
    20 
    2120
    2221
    … …  
    2625clearscript IntAt1420
    2726clearscript sigmaRaw1420Side
     27clearscript sigmaRaw1400a1420
    2828############################################################################
    2929defscript evolAtCalibFreq
    … …  
    150150
    151151
    152 
    153 
    154 #col 13 = onoffRaw01420side
    155 #col 14 = onoffRaw11420side
     152#col 13 = onoffRaw01420side   Ch 0
     153#col 14 = onoffRaw11420side   Ch 1
    156154
    157155ntcol2var onoffevol 13 linCh0All
    … …  
    161159line2vec vecCh1All $linCh1All
    162160
    163 
     161#packing values only valid for 500 cycles...
    164162set packVal "1 10 25 50 100"
    165163set sigmaCh0 ""
    … …  
    188186  mean1 = ${outvec.sum}/${outvec.size}
    189187  sigma1 = sqrt(${outvec.sumsq}/${outvec.size}-${mean1}*${mean1})
    190   errsig1 = ${sigma0}/sqrt(2*${outvec.size})   
     188  errsig1 = ${sigma1}/sqrt(2*${outvec.size})   
    191189  set sigmaCh1 "${sigmaCh1} ${sigma1}" 
    192190  set errsigCh1 "${errsigCh1} ${errsig1}"
    … …  
    222220
    223221endscript
     222##################################################
     223defscript sigmaRaw1400a1420
     224
     225c++compile rebining
     226c++link rebining.so dorebin
     227
     228
     229#col 19 = onoffRaw0f14001420   Ch 0
     230#col 20 = onoffRaw1f14001420  Ch 1
     231
     232ntcol2var onoffevol 19 linCh0All
     233line2vec vecCh0All $linCh0All
     234
     235ntcol2var onoffevol 20 linCh1All
     236line2vec vecCh1All $linCh1All
     237
     238#packing values only valid for 500 cycles...
     239set packVal "1 10 25 50 100"
     240set sigmaCh0 ""
     241set errsigCh0 ""
     242set sigmaCh1 ""
     243set errsigCh1 ""
     244
     245
     246foreach ipack ( $packVal )
     247
     248#Packing per ipack
     249  del invec
     250  del outvec
     251  cp vecCh0All invec
     252  call dorebin $ipack
     253  mean0 = ${outvec.sum}/${outvec.size}
     254  sigma0 = sqrt(${outvec.sumsq}/${outvec.size}-${mean0}*${mean0})
     255  errsig0 = ${sigma0}/sqrt(2*${outvec.size})   
     256  set sigmaCh0 "${sigmaCh0} ${sigma0}" 
     257  set errsigCh0 "${errsigCh0} ${errsig0}"       
     258#
     259  del invec
     260  del outvec
     261  cp vecCh1All invec
     262  call dorebin $ipack
     263  mean1 = ${outvec.sum}/${outvec.size}
     264  sigma1 = sqrt(${outvec.sumsq}/${outvec.size}-${mean1}*${mean1})
     265  errsig1 = ${sigma1}/sqrt(2*${outvec.size})   
     266  set sigmaCh1 "${sigmaCh1} ${sigma1}" 
     267  set errsigCh1 "${errsigCh1} ${errsig1}"
     268end
     269
     270
     271newnt nt0 x y ey
     272newnt nt1 x y ey
     273
     274set x0 ( $packVal )
     275set y0 ( $sigmaCh0 )
     276set ey0 ( $errsigCh0 )
     277
     278set x1 ( $packVal )
     279set y1 ( $sigmaCh1 )
     280set ey1 ( $errsigCh1 )
     281       
     282for i 0:$#x0
     283 line2nt nt0 $x0[i] $y0[i] $ey0[i]
     284 line2nt nt1 $x1[i] $y1[i] $ey1[i]
     285end
     286
     287
     288newwin 1 1
     289plot2de nt0 x y 0 ey 1 "blue marker=fcircle,9 notit nsta"     
     290plot2de nt1 x y 0 ey 1 "same red  marker=fcircle,9 notit nsta"
     291n = ${#x0}-1
     292func $y0[0]/sqrt(x) $x0[0] $x0[n]  100 "same"
     293settitle "Sigma [1400,1420]MHz ${source} Ch 0 (blue) Ch 1 (red)" ' ' $defatt
     294setaxelabels "num of cycles" "Sigma (a.u)" $axedefatt
     295
     296
     297
     298endscript
    224299###########################################
     300
    225301defscript evolAt1420SideFreq
    226302
    … …  
    361437# Calibration coeff. per cycles
    362438#
     439
     440graphicatt "xylimits=1250,1500,-0.01,0.01"
     441del sonoffovoff0 sonoffovoff1
     442objaoper meanOvOffNoCalib row 0 sonoffovoff0
     443objaoper meanOvOffNoCalib row 1 sonoffovoff1
     444newwin 1 1
     445plot2d sonoffovoff0 (n/8192)*250+1250 val n>0 "blue cpts notit nsta"
     446plot2d sonoffovoff1 (n/8192)*250+1250 val n>0 "red same cpts notit nsta"
     447settitle "Raw (ON-OFF)/OFF  ${source} Ch 0 (blue) Ch 1 (red) ${ncycles}cycles" ' ' $defatt
     448setaxelabels "Freq. (MHz)" "I (a.u)" $axedefatt
    363449
    364450

  • BAORadio/AmasNancay/trunk/mergeAnaFiles.cc

    r587 r591  
    111111
    112112//-------------
     113//JEC 25/10/11 Perform a median filtering with a sliding window of "halfwidth" half width
     114//             It takes care of the edges and is based on the median function (above)
     115void medianFiltering(const TMatrix<r_4> mtx,
     116                     sa_size_t halfwidth,
     117                     TMatrix<r_4>& vec) {
     118 
     119  sa_size_t nr = mtx.NRows();
     120  sa_size_t nc = mtx.NCols();
     121  sa_size_t chMin = 0;
     122  sa_size_t chMax = nc-1;
     123 
     124  for (sa_size_t ir=0; ir<nr; ir++){
     125    for (sa_size_t ic=0; ic<nc; ic++) {
     126      sa_size_t chLow = ic-halfwidth;
     127      chLow = (chLow >= chMin) ? chLow : chMin;
     128      sa_size_t chHigh = ic+halfwidth;
     129      chHigh = ( chHigh <= chMax ) ? chHigh : chMax;
     130      TVector<r_4> tmp(mtx(Range(ir),Range(chLow,chHigh)),false);
     131      vector<r_4> val;
     132      tmp.FillTo(val);
     133      vec(ir,ic) = median(val.begin(),val.end());
     134    }
     135  }
     136}
     137//-------------
    113138void split(const string& str, const string& delimiters , vector<string>& tokens) {
    114139    // Skip delimiters at beginning.
    … …  
    264289
    265290  //mean ON-OFF over the list of cycles
     291  TMatrix<r_4> meanDiffONOFFovOFF_noCalib(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);          //set to 0
    266292  TMatrix<r_4> meanDiffONOFF_noCalib(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);          //set to 0
    267293  TMatrix<r_4> meanDiffONOFF_perRunCalib(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);      //set to 0
    268294  TMatrix<r_4> meanDiffONOFF_perCycleCalib(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);    //set to 0
    269   static const int NINFO=19;
     295  static const int NINFO=21;
    270296  char* onffTupleName[NINFO]={"cycle"
    271                           ,"onoffRaw0","onoffRaw1"
    272                           ,"onoffRun0","onoffRun1"
    273                           ,"onoffCycle0","onoffCycle1"
    274                           ,"onoffRaw01420","onoffRaw11420"
    275                           ,"onoffRun01420","onoffRun11420"
    276                           ,"onoffCycle01420","onoffCycle11420"
    277                           ,"onoffRaw01420side","onoffRaw11420side"
    278                           ,"onoffRun01420side","onoffRun11420side"
    279                           ,"onoffCycle01420side","onoffCycle11420side"
     297                              ,"onoffRaw0","onoffRaw1"
     298                              ,"onoffRun0","onoffRun1"
     299                              ,"onoffCycle0","onoffCycle1"
     300                              ,"onoffRaw01420","onoffRaw11420"
     301                              ,"onoffRun01420","onoffRun11420"
     302                              ,"onoffCycle01420","onoffCycle11420"
     303                              ,"onoffRaw01420side","onoffRaw11420side"
     304                              ,"onoffRun01420side","onoffRun11420side"
     305                              ,"onoffCycle01420side","onoffCycle11420side"
     306                              ,"onoffRaw0f14001420","onoffRaw1f14001420"
    280307  };
    281308  NTuple onoffevolution(NINFO,onffTupleName);
    … …  
    294321  sa_size_t ch1420bLow  = freqToChan(1422);
    295322  sa_size_t ch1420bHigh = freqToChan(1423);
    296 
     323 
     324  //25/10/11 follow 1400-1420Mhz
     325  sa_size_t ch1420 = freqToChan(1420);
     326  sa_size_t ch1400 = freqToChan(1400);
    297327
    298328  if (para.debuglev_>0){
    … …  
    631661      meanDiffONOFF_noCalib += diffOnOff_noCalib;
    632662     
     663      //JEC 29/10/11 add ON-OFF/OFF
     664      TMatrix<r_4> diffOnOffOvOff_noCalib(diffOnOff_noCalib,false); //do not share data
     665      TMatrix<r_4> aSpecOffFitltered(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);
     666      sa_size_t halfWidth = 35; //number of freq. bin for the 1/2 width of the filtering window
     667      medianFiltering(aSpecOff,halfWidth,aSpecOffFitltered);
     668     
     669      diffOnOffOvOff_noCalib.Div(aSpecOffFitltered); //division in place
     670      meanDiffONOFFovOFF_noCalib += diffOnOffOvOff_noCalib;
     671
    633672      TMatrix<r_4> diffOnOff_perRunCalib = aSpecOn_BAOCalibRun - aSpecOff_BAOCalibRun;
    634673      meanDiffONOFF_perRunCalib += diffOnOff_perRunCalib;
    … …  
    703742      xnt[18] = (meanInRange_1420aFreq_perCycle(1) + meanInRange_1420bFreq_perCycle(1))/2.;
    704743
     744
     745      //JEC 25/10/11 follow 1400-1420 MHz bande protege et n'inclue pas le 1420.4Mhz de la Galaxie
     746      TVector<r_4> meanInRange_1400a1420Freq_noCalib(NUMBER_OF_CHANNELS);
     747      meanInRange(diffOnOff_noCalib,ch1400,ch1420,meanInRange_1400a1420Freq_noCalib);
     748      xnt[19] = meanInRange_1400a1420Freq_noCalib(0);
     749      xnt[20] = meanInRange_1400a1420Freq_noCalib(1);
     750     
     751
    705752     
    706753      //store infos to Ntuple
    … …  
    717764  //Normalisation
    718765  if(totalNumberCycles > 0){
     766    //JEC 29/10 add ON-OFF/OFF
     767    meanDiffONOFFovOFF_noCalib  /= (r_4)totalNumberCycles;
    719768    meanDiffONOFF_noCalib       /= (r_4)totalNumberCycles;
    720769    meanDiffONOFF_perRunCalib   /= (r_4)totalNumberCycles;
    … …  
    745794    string tag = "meanNoCalib";
    746795    fos << PPFNameTag(tag) << meanDiffONOFF_noCalib;
     796   
     797    //JEC 29/10/11
     798    tag = "meanOvOffNoCalib";
     799    fos << PPFNameTag(tag) << meanDiffONOFFovOFF_noCalib;
     800
    747801    tag = "meanPerRunCalib";
    748802    fos << PPFNameTag(tag) << meanDiffONOFF_perRunCalib;
Note: See TracChangeset for help on using the changeset viewer.