Changeset 616 for BAORadio


Ignore:
Timestamp:
Nov 28, 2011, 3:59:31 PM (13 years ago)
Author:
torrento
Message:

Added calibcoeffNtp function (dismiss previous one)

File:
1 edited

Legend:

Unmodified
Added
Removed

  • BAORadio/AmasNancay/trunk/mergeAnaFiles.cc

    r615 r616  
    6363//--------------------------------------------------------------
    6464//Utility functions
     65// template <class T>
     66// void MeanSigma(TArray<T> const & a, T & mean, T & sig){
     67//   if (a.NbDimensions() < 1)
     68//     throw RangeCheckError("MathArray<T>::MeanSigma(TArray<T> const & a..) Not Allocated Array a !");
     69//   const T * pe;
     70//   sa_size_t j,k;
     71//   mean=0.;
     72//   sig = 0.;
     73//   T valok;
     74//   if (a.AvgStep() > 0)   {  // regularly spaced elements
     75//     sa_size_t step = a.AvgStep();
     76//     sa_size_t maxx = a.Size()*step;
     77//     pe = a.Data();
     78//     for(k=0; k<maxx; k+=step )  {
     79//       valok = pe[k];
     80//       mean += valok;  sig += valok*valok;
     81//     }
     82//   }
     83//   else {    // Non regular data spacing ...
     84//     int_4 ka = a.MaxSizeKA();
     85//     sa_size_t step = a.Step(ka);
     86//     sa_size_t gpas = a.Size(ka)*step;
     87//     sa_size_t naxa = a.Size()/a.Size(ka);
     88//     for(j=0; j<naxa; j++)  {
     89//       pe = a.DataBlock().Begin()+a.Offset(ka,j);
     90//       for(k=0; k<gpas; k+=step) {
     91//      valok = pe[k];
     92//      mean += valok;  sig += valok*valok;
     93//       }
     94//     }
     95//   }
     96//   T dsz = (T)(a.Size());
     97//   mean /= dsz;
     98//   if (dsz > 1.5) {
     99//     sig = sig/dsz - mean*mean;
     100//     sig *= (dsz/(dsz-1));
     101//     if (sig >= 0.) sig = sqrt(sig);
     102//   }
     103//   else sig = 0.;
     104// }
     105
    65106
    66107sa_size_t freqToChan(r_4 f){
    … …  
    79120  for (sa_size_t ir=0; ir<nr; ir++){
    80121    TVector<r_4> tmp(mtx(Range(ir),Range(chLow,chHigh)),false);
    81     double mean, sigma;
     122    r_8 mean, sigma;
    82123    MeanSigma(tmp,mean,sigma);
    83124    vec(ir) = mean;
    … …  
    273314  }//eo loop on slices
    274315}
     316
    275317//AST 18/11/11
    276318//-------------------------------------------------------
    … …  
    449491    }
    450492    ifs.close();
    451 
    452493    //
    453494    calibFileName = directoryName + "/"
    … …  
    621662    }
    622663    ifs.close();
    623 
    624664    //link <cycle> - <calibration coefficient>
    625665    //We cannot rely on identical cycle list of the OFF and ON calibration
    … …  
    770810
    771811}
     812
    772813//-------------------------------------------------------
    773814//Compute the mean of Diff ON-OFF BAO-calibrated spectra and also the mean/sigma of rebinned spectra
    … …  
    812853  sa_size_t chCalibHigh = freqToChan(para.rcalibFreq_ + (para.rcalibBandFreq_*0.5));
    813854  //Lower and Higher freq.  just arround 1420.4MHz Freq. bin to perform mean follow up
    814   sa_size_t ch1420Low;
    815   sa_size_t ch1420High;
    816   if (para.sourceName_ == "Abell85") {
    817     ch1420Low  = freqToChan(1420.4-0.2);  // Abell85
    818     ch1420High = freqToChan(1420.4+0.2);
    819   } else if (para.sourceName_ == "Abell1205") {
    820     ch1420Low  = freqToChan(1420.4-0.3);  // Abell1205
    821     ch1420High = freqToChan(1420.4+0.2);
    822   } else if (para.sourceName_ == "Abell2440") {
    823     ch1420Low  = freqToChan(1420.4);      // Abell2440
    824     ch1420High = freqToChan(1420.4+0.3);
    825   } else {
    826     ch1420Low  = freqToChan(1420.4-0.2);  // Abell85
    827     ch1420High = freqToChan(1420.4+0.2);
    828   }
     855  sa_size_t ch1420Low  = freqToChan(1420.4-0.2);
     856  sa_size_t ch1420High = freqToChan(1420.4+0.2);
     857
    829858  //Lower and Higher freq. on the sides of 1420.4Mhz Freq. bin to perform mean follow up
    830859  sa_size_t ch1420aLow  = freqToChan(1418);
    … …  
    865894    string srcLower = tokens2[0];
    866895
    867 
    868    
    869896    PInPersist fin(*iFile);
    870897    vector<string> vec = fin.GetNameTags();
    … …  
    10821109    nr = calibBAOfactors_Off_Run_Ch0.NRows();
    10831110    for (sa_size_t ir=0; ir<nr; ++ir){
    1084 //       cout << "Calib. Off Run Ch0 cycle ["<< calibBAOfactors_Off_Run_Ch0(ir,0)<<"], val "
    1085 //         << calibBAOfactors_Off_Run_Ch0(ir,1) << endl;
     1111      cout << "Calib. Off Run Ch0 cycle ["<< calibBAOfactors_Off_Run_Ch0(ir,0)<<"], val "
     1112           << calibBAOfactors_Off_Run_Ch0(ir,1) << endl;
    10861113
    10871114      calibBAO_Off_Run_Ch0[(int)calibBAOfactors_Off_Run_Ch0(ir,0)]
    … …  
    11341161               << "Ch1 " << calibBAO_On_Cycle_Ch1[*ic] << endl;
    11351162        } else {
    1136           cout << "Cycle NOT calibrated " << endl;
     1163          cout << "Cycle << " << *ic <<" NOT calibrated for file " << *iFile << endl;
    11371164        }
    11381165      }//debug
    … …  
    11781205      //JEC 29/10/11 add ON-OFF/OFF
    11791206      TMatrix<r_4> diffOnOffOvOff_noCalib(diffOnOff_noCalib,false); //do not share data
    1180       TMatrix<r_4> aSpecOffFitltered(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);
     1207      TMatrix<r_4> aSpecOffFiltered(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);
    11811208      sa_size_t halfWidth = 35; //number of freq. bin for the 1/2 width of the filtering window
    1182       medianFiltering(aSpecOff,halfWidth,aSpecOffFitltered);
     1209      medianFiltering(aSpecOff,halfWidth,aSpecOffFiltered);
    11831210     
    1184       diffOnOffOvOff_noCalib.Div(aSpecOffFitltered); //division in place
     1211      diffOnOffOvOff_noCalib.Div(aSpecOffFiltered); //division in place
    11851212      meanDiffONOFFovOFF_noCalib += diffOnOffOvOff_noCalib;
    11861213
    … …  
    11911218      meanDiffONOFF_perCycleCalib += diffOnOff_perCycleCalib;
    11921219
    1193       //
    11941220      totalNumberCycles++;
     1221
    11951222      //Fill NTuple
    11961223      xnt[0] = totalNumberCycles;
    … …  
    12661293      //JEC 18/11/11 follow up the 1400-1420MHz OFF only
    12671294      TMatrix<r_4> aSpecOffovOff(aSpecOff,false);
    1268       aSpecOffovOff.Div(aSpecOffFitltered);
     1295      aSpecOffovOff.Div(aSpecOffFiltered);
    12691296
    12701297      TVector<r_4> meanInRange_1410a1415Freq_OFF_noCalib(NUMBER_OF_CHANNELS);
     1298      //      meanInRange(aSpecOff,ch1410,ch1415,meanInRange_1410a1415Freq_OFF_noCalib);
    12711299      meanInRange(aSpecOffovOff,ch1410,ch1415,meanInRange_1410a1415Freq_OFF_noCalib);
    12721300
    … …  
    12751303
    12761304      TMatrix<r_4> aSpecOnovOff(aSpecOn,false);
    1277       aSpecOnovOff.Div(aSpecOffFitltered);
     1305      aSpecOnovOff.Div(aSpecOffFiltered);
    12781306
    12791307      TVector<r_4> meanInRange_1410a1415Freq_ON_noCalib(NUMBER_OF_CHANNELS);
     1308      //meanInRange(aSpecOn,ch1410,ch1415,meanInRange_1410a1415Freq_ON_noCalib);
    12801309      meanInRange(aSpecOnovOff,ch1410,ch1415,meanInRange_1410a1415Freq_ON_noCalib);
    12811310
    … …  
    13521381   
    13531382    tag = "onoffevol";
    1354     fos << PPFNameTag(tag) << onoffevolution;   
     1383    fos << PPFNameTag(tag) << onoffevolution;
     1384 
    13551385  }//end of save
     1386
     1387
    13561388}
    13571389//JEC 14/11/11 New meanRawDiffOnOffCycles START
    … …  
    13761408  TMatrix<r_4> meanONovOFF_noCalib(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ); //  ON/Filtered_OFF
    13771409  TMatrix<r_4> meanOFFovOFF_noCalib(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ); // OFF/Filtered_OFF
    1378 
    13791410  //Tuple only for RAW things to follow
    13801411  static const int NINFO=11;
    … …  
    14171448  for (iFile = listOfFiles.begin(); iFile != iFileEnd; ++iFile) {
    14181449    if (para.debuglev_>90){
    1419     }
     1450      cout << "load file <" << *iFile << ">" << endl;
     1451    }
     1452
     1453    vector<string> tokens;
     1454    split(*iFile,"_",tokens);
     1455    string dateOfRun = tokens[1];
     1456    if (para.debuglev_>90){
     1457      cout << "date <" << dateOfRun << ">" << endl;
     1458    }
     1459    vector<string> tokens2;
     1460    split(tokens[2],".",tokens2);
     1461    string srcLower = tokens2[0];
     1462
    14201463    PInPersist fin(*iFile);
    14211464    vector<string> vec = fin.GetNameTags();
    … …  
    14471490        regexp(iSpec->c_str(),matchstr.c_str(),&b,&e);
    14481491        if (para.debuglev_>90){
    1449           cout << " spactra <" << *iSpec << ">" << endl;
     1492          cout << " spectra <" << *iSpec << ">" << endl;
    14501493          cout << " cycle " << iSpec->substr(b) << endl;
    14511494        }
    … …  
    14751518    for (list<int>::iterator ic=commonCycles.begin(); ic!=commonCycles.end(); ++ic){
    14761519     
     1520      // AST 28.11.11 remove non-calibrated cycles for Abell1205 and Abell2440
     1521      if ( *ic == 1 && srcLower == "abell1205" ) {
     1522          if ( dateOfRun == "20110502" || dateOfRun == "20110607" || dateOfRun == "20110818" ) {
     1523            cout << "Skipping non-calibrated cycle " << *ic << endl;
     1524            continue;
     1525          }
     1526      } else if ( *ic == 1 && srcLower == "abell2440" ) {
     1527          if ( dateOfRun == "20110516" ) {
     1528            cout << "Skipping non-calibrated cycle " << *ic << endl;
     1529            continue;
     1530          }
     1531      } else if ( *ic == 3 && srcLower == "abell1205" ) {
     1532          if ( dateOfRun == "20110810" ) {
     1533            cout << "Skipping non-calibrated cycle " << *ic << endl;
     1534            continue;
     1535          }
     1536      }
     1537
    14771538      string ppftag;
    14781539      //load ON phase
    … …  
    14901551      //Perform the difference ON-OFF
    14911552      TMatrix<r_4> diffOnOff_noCalib = aSpecOn - aSpecOff;
     1553
    14921554      meanDiffONOFF_noCalib += diffOnOff_noCalib;
    14931555     
    14941556      //JEC 29/10/11 add ON-OFF/OFF
    14951557      TMatrix<r_4> diffOnOffOvOff_noCalib(diffOnOff_noCalib,false); //do not share data
    1496       TMatrix<r_4> aSpecOffFitltered(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);
     1558      TMatrix<r_4> aSpecOffFiltered(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);
    14971559      sa_size_t halfWidth = 35; //number of freq. bin for the 1/2 width of the filtering window
    1498       medianFiltering(aSpecOff,halfWidth,aSpecOffFitltered);
     1560      medianFiltering(aSpecOff,halfWidth,aSpecOffFiltered);
    14991561     
    1500       diffOnOffOvOff_noCalib.Div(aSpecOffFitltered); //division in place
     1562      diffOnOffOvOff_noCalib.Div(aSpecOffFiltered); //division in place
    15011563      meanDiffONOFFovOFF_noCalib += diffOnOffOvOff_noCalib;
    1502 
    15031564
    15041565      //JEC 15/11/11 add ON/OFF and OFF/OFF
    15051566      TMatrix<r_4> onOvOff(aSpecOn,false);
    1506       onOvOff.Div(aSpecOffFitltered);
     1567      onOvOff.Div(aSpecOffFiltered);
    15071568      meanONovOFF_noCalib += onOvOff;
    15081569     
    15091570      TMatrix<r_4> offOvOff(aSpecOff,false);
    1510       offOvOff.Div(aSpecOffFitltered);
     1571      offOvOff.Div(aSpecOffFiltered);
    15111572      meanOFFovOFF_noCalib += offOvOff;
    15121573
    15131574      totalNumberCycles++;
    1514 
    15151575
    15161576      //Fill NTuple
    15171577      xnt[0] = totalNumberCycles;
    15181578 
    1519 
    15201579      //Follow up arround the 1420.4MHz Freq.
    15211580      TVector<r_4> meanInRange_1420Freq_noCalib(NUMBER_OF_CHANNELS);
    … …  
    15351594
    15361595
    1537       //JEC 25/10/11 follow 1400-1420 MHz bande protege et n'inclue pas le 1420.4Mhz de la Galaxie
    1538       TVector<r_4> meanInRange_1400a1420Freq_noCalib(NUMBER_OF_CHANNELS);
    1539       meanInRange(diffOnOffOvOff_noCalib,ch1400,ch1420,meanInRange_1400a1420Freq_noCalib);
    1540       xnt[5] = meanInRange_1400a1420Freq_noCalib(0);
    1541       xnt[6] = meanInRange_1400a1420Freq_noCalib(1);
    1542 
    1543       //JEC 18/11/11 follow up the 1400-1420MHz OFF only
     1596      //JEC 25/10/11 follow 1410-1415 MHz bande protege et n'inclue pas le 1420.4Mhz de la Galaxie
     1597      TVector<r_4> meanInRange_1410a1415Freq_noCalib(NUMBER_OF_CHANNELS);
     1598      meanInRange(diffOnOffOvOff_noCalib,ch1410,ch1415,meanInRange_1410a1415Freq_noCalib);
     1599      xnt[5] = meanInRange_1410a1415Freq_noCalib(0);
     1600      xnt[6] = meanInRange_1410a1415Freq_noCalib(1);
     1601
     1602      //JEC 18/11/11 follow up the 1410-1415MHz OFF only
    15441603      TMatrix<r_4> aSpecOffovOff(aSpecOff,false);
    1545       aSpecOffovOff.Div(aSpecOffFitltered);
     1604      aSpecOffovOff.Div(aSpecOffFiltered);
    15461605
    15471606      TVector<r_4> meanInRange_1410a1415Freq_OFF_noCalib(NUMBER_OF_CHANNELS);
    … …  
    15521611
    15531612      TMatrix<r_4> aSpecOnovOff(aSpecOn,false);
    1554       aSpecOnovOff.Div(aSpecOffFitltered);
     1613      aSpecOnovOff.Div(aSpecOffFiltered);
    15551614
    15561615      TVector<r_4> meanInRange_1410a1415Freq_ON_noCalib(NUMBER_OF_CHANNELS);
    … …  
    16001659
    16011660    tag = "onoffevol";
    1602     fos << PPFNameTag(tag) << onoffevolution;   
     1661    fos << PPFNameTag(tag) << onoffevolution;
    16031662
    16041663  }//end save
     1664
     1665
    16051666}
    16061667//JEC 14/11/11 New meanRawDiffOnOffCycles END
    … …  
    19622023    } else if (action == "meanCalibBAODiffOnOff") {
    19632024      meanCalibBAODiffOnOffCycles();
    1964     } else if (action == "calibCoeffNtp") {
    1965       calibCoeffNtp();
    19662025    } else {
    19672026      msg = "Unknown action " + action;
Note: See TracChangeset for help on using the changeset viewer.