Changeset 598


Ignore:
Timestamp:
Nov 10, 2011, 10:14:11 AM (13 years ago)
Author:
campagne
Message:

cleaning of the set of scripts, and improve th submit2ge options, and improve analysis (jec)

Location:
BAORadio/AmasNancay/trunk
Files:
15 deleted
9 edited

Legend:

Unmodified
Added
Removed

  • BAORadio/AmasNancay/trunk/anagainmaker.sh

    r583 r598  
    6464    Abell2440) $ECHO "INFO (${scriptName}): process ${sourceRadio}";;
    6565    Abell1205) $ECHO "INFO (${scriptName}): process ${sourceRadio}";;
     66    NGC4383) $ECHO "INFO (${scriptName}): process ${sourceRadio}";;
    6667    *) ECHO "FATAL (${scriptName}): process ${sourceRadio} not yet foreseen"
    6768    exit 1;;

  • BAORadio/AmasNancay/trunk/analyse.cc

    r595 r598  
    3737#include "fiosinit.h"
    3838#include "ppersist.h"
     39
     40
     41//~/private/work/AmasNancay/Objs/analyse -act redMeanONOFF -source Abell85 -date 20110428 -sca sca151675.sum.trans -specdir meancycle -specname mspecmtx -sigmaname sigspecmtx -npaq 25000 -numcycle 1,11 -debug 100 > & redMean.log
     42
    3943
    4044
    … …  
    256260};
    257261
     262//------------
     263//Process ON/OFF data treated with the mspec (specmfib) => no filtering at all: paquets & freq.
     264//------------
     265class ProcessONOFFReducedMeanData : public ProcessBase {
     266protected:
     267  string sigma2FileName_; //name of the file which contains the sigma2
     268  string nPaqPerWin_; // number of paquets used for mean and sigma2 computations 'cf. proc_specmfib'
     269  r_4 valNPaqPerWin_; // value
     270
     271public:
     272  ProcessONOFFReducedMeanData(){}
     273  virtual ~ProcessONOFFReducedMeanData(){}
     274 
     275  void SetSigma2FileName(const string& val) {sigma2FileName_ = val;}
     276  void SetNPaqPerWin(const string& val) {
     277    nPaqPerWin_    = val;
     278    valNPaqPerWin_ = atof(val.c_str());
     279  }
     280
     281  virtual int processCmd() throw(string);
     282};
     283
     284
     285//------------
     286//Process ON/OFF data treated with the gain (specmfib) => filtering paquets only by default
     287//------------
     288class ProcessONOFFReducedMedianData : public ProcessBase {
     289protected:
     290  string nPaqPerWin_; // number of paquets used for mean and sigma2 computations 'cf. proc_specmfib'
     291  r_4 valNPaqPerWin_; // value
     292
     293public:
     294  ProcessONOFFReducedMedianData(){}
     295  virtual ~ProcessONOFFReducedMedianData(){}
     296 
     297  void SetNPaqPerWin(const string& val) {
     298    nPaqPerWin_    = val;
     299    valNPaqPerWin_ = atof(val.c_str());
     300  }
     301
     302  virtual int processCmd() throw(string);
     303};
     304
    258305
    259306//------------
    … …  
    309356  upperFreqBin_ = c0+dc/2;
    310357}
    311 
    312 
    313 //----------------------------------------------------
    314358//----------------------------------------------------
    315359int main(int narg, char* arg[])
    … …  
    318362  //Init process types
    319363  map<string,IProcess*> process;
     364  process["redMedianONOFF"] = new ProcessONOFFReducedMedianData(); //JEC  8/11/11
     365  process["redMeanONOFF"] = new ProcessONOFFReducedMeanData(); //JEC  8/11/11
    320366  process["meanONOFF"] = new ProcessONOFFMeanData(); //JEC 27/10/11
    321367  process["rawOnOff"]  = new ProcessONOFFRawData(); 
    … …  
    352398  string numcycle;
    353399  string calibrationOpt = ""; 
     400  string nPaqPerWin = "";
    354401
    355402  string typeOfFile="medfiltmtx";
    356  
     403  string secondTypeOfFile=""; //valid for ProcessONOFFReducedMeanData JEC 8/11/11
    357404
    358405  //  bool okarg=false;
    … …  
    396443      ka+=2;
    397444    }   
     445    else if (strcmp(arg[ka],"-sigmaname")==0) {
     446      secondTypeOfFile=arg[ka+1];
     447      ka+=2;
     448    }   
     449    else if (strcmp(arg[ka],"-npaq")==0) {
     450      nPaqPerWin=arg[ka+1];
     451      ka+=2;
     452    }
    398453    else if (strcmp(arg[ka],"-freqBAOCalib")==0) {
    399454      freqBAOCalib = arg[ka+1];
    … …  
    430485       << " spectraDirectory = " << spectraDirectory << "\n"
    431486       << " spectraName = " << typeOfFile << "\n"
     487       << " sigmaName = " << secondTypeOfFile << "\n"
     488       << " nPaqPerWin = " << nPaqPerWin << "\n"
    432489       << " freqBAOCalib = " << freqBAOCalib  << "\n"
    433490       << " bandWidthBAOCalib = " << bandWidthBAOCalib << "\n"
    434        << " debuglev = "  << debuglev  << "\n"
     491       << " debug = "  << debuglev  << "\n"
    435492       << " mode = " << mode << "\n"
    436493       << " numcycle = " << numcycle << "\n"
    … …  
    517574//     }
    518575    //JEC 22/9/11 Make ON-OFF analysis WO any calibration END
     576
     577    //JEC 8/11/11 Start
     578    try {
     579      ProcessONOFFReducedMeanData* procdata=dynamic_cast<ProcessONOFFReducedMeanData*>(process[action]);
     580      if (procdata) {
     581        if (secondTypeOfFile == ""){
     582          msg = "(FATAL) missing file of the sigmas for action " + action;
     583          Usage(true);
     584          throw msg;     
     585        }
     586        if (nPaqPerWin == ""){
     587          msg = "(FATAL) missing number of paquets per window for action " + action;
     588          Usage(true);
     589          throw msg;     
     590        }
     591        procdata->SetSigma2FileName(secondTypeOfFile);
     592        procdata->SetNPaqPerWin(nPaqPerWin);
     593      }
     594    }
     595    catch(exception& e){
     596      throw e.what();
     597    }
     598    //JEC 8/11/11 End
     599
     600    //JEC 8/11/11 Start
     601    try {
     602      ProcessONOFFReducedMedianData* procdata=dynamic_cast<ProcessONOFFReducedMedianData*>(process[action]);
     603      if (procdata) {
     604        if (nPaqPerWin == ""){
     605          msg = "(FATAL) missing number of paquets per window for action " + action;
     606          Usage(true);
     607          throw msg;     
     608        }
     609        procdata->SetNPaqPerWin(nPaqPerWin);
     610      }
     611    }
     612    catch(exception& e){
     613      throw e.what();
     614    }
     615    //JEC 8/11/11 End
     616   
    519617
    520618
    … …  
    611709int Usage(bool flag) {
    612710  cout << "Analyse.cc usage...." << endl;
    613   cout << "analyse  -act <action_type>: meanONOFF,dataOnOff, rawOnOff, gain, calib\n"
     711  cout << "analyse  -act <action_type>: redMeanONOFF, meanONOFF,dataOnOff, rawOnOff, gain, calib\n"
    614712       << "         -inPath <path for input files: default='.'>\n"
    615713       << "         -outPath <path for output files: default='.'>\n"
    … …  
    619717       << "         -specdir <generic directory name of spectra fits files>\n"
    620718       << "         -specname <generic name of spectra fits files>\n"
     719       << "         -sigmaname <generic name of the sigma2 fits files>\n"
     720       << "            valid for redMeanONOFF\n"
     721       << "         -npaq <number of paquest per windows for mean&sigma2 computaion (brproc)>\n"
     722       << "            valid for redMeanONOFF\n"
    621723       << "         -freqBAOCalib <freq in MHZ> freq. of calibration BAO\n"
    622724       << "            valid for act=dataOnOff\n"
    … …  
    710812  string tag;
    711813
    712   map< pair<string, sa_size_t>, TMatrix<r_4> > spectreCollect;
     814  map< pair<string, sa_size_t>, TMatrix<r_4> > meanCollect;
     815  map< pair<string, sa_size_t>, TMatrix<r_4> > sigma2Collect;
    713816  map< pair<string, sa_size_t>, TMatrix<r_4> >::iterator iSpectre, iSpectreEnd;
    714817 
    … …  
    744847      if(debuglev_>0)cout << "compute mean for cycle " << icycle << endl;
    745848      TMatrix<r_4> spectreMean(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ); //implicit init to 0
     849      TMatrix<r_4> spectreSigma2(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ); //implicit init to 0
    746850      iFileEnd = listOfSpecFiles.end();
    747851      r_4 nSpectres  = 0;
    … …  
    750854        TMatrix<r_4> spectre(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);
    751855        aSpectrum >> spectre;
    752         spectreMean += spectre;
     856        spectreMean   += spectre;
     857        spectreSigma2 += spectre && spectre;
    753858        nSpectres++;
    754859      }// end of for files
    755       if(nSpectres>0) spectreMean /= nSpectres;
     860      if(nSpectres>0) {
     861        spectreMean /= nSpectres;
     862        spectreSigma2 /= nSpectres;
     863        spectreSigma2 -= spectreMean && spectreMean;
     864      }
    756865     
    757866      //save mean spectrum
    758       spectreCollect.insert( pair< pair<string,sa_size_t>, TMatrix<r_4> >(make_pair(mode,icycle),TMatrix<r_4>(spectreMean,false) ));
     867      meanCollect.insert( pair< pair<string,sa_size_t>, TMatrix<r_4> >(make_pair(mode,icycle),TMatrix<r_4>(spectreMean,false) ));
     868      sigma2Collect.insert( pair< pair<string,sa_size_t>, TMatrix<r_4> >(make_pair(mode,icycle),TMatrix<r_4>(spectreSigma2,false) ));
    759869    }//end of for cycles
    760870  }//end loop on mode for raw preocess
    … …  
    764874    string fileName;
    765875
    766     fileName = "./dataMean_" + dateOfRun_ + "_" + StringToLower(sourceName_) + ".ppf";
     876    fileName = "./dataMeanSigma2_" + dateOfRun_ + "_" + StringToLower(sourceName_) + ".ppf";
    767877
    768878    POutPersist fos(fileName);
    769879    id=0;
    770     iSpectreEnd = spectreCollect.end();
    771     for (iSpectre = spectreCollect.begin();
     880    iSpectreEnd = meanCollect.end();
     881    for (iSpectre = meanCollect.begin();
    772882         iSpectre != iSpectreEnd ; ++iSpectre, ++id) {
    773883      tag = "specMean";
    … …  
    780890        cout << "save tag<" << tag << ">" << endl;
    781891      }
    782 
    783892      fos << PPFNameTag(tag) << iSpectre->second;
    784893    }
     894
     895
     896    id=0;
     897    iSpectreEnd = sigma2Collect.end();
     898    for (iSpectre = sigma2Collect.begin();
     899         iSpectre != iSpectreEnd ; ++iSpectre, ++id) {
     900      tag = "specSigma2";
     901
     902      tag += (iSpectre->first).first;
     903      stringstream sid;
     904      sid << (iSpectre->first).second;
     905      tag += sid.str();
     906      if(debuglev_>9) {
     907        cout << "save tag<" << tag << ">" << endl;
     908      }
     909      fos << PPFNameTag(tag) << iSpectre->second;
     910    }
     911
    785912  }//end of save fits
    786913 
    … …  
    790917
    791918  cout << "OK ProcessONOFFMeanData finished" <<endl;
     919  return rc;
     920}
     921//----------------------------------------------
     922//----------------------------------------------
     923//JEC 8/11/11
     924//Process the files that are output of the specmfib -act = mspec (mean+sigma of signal files)
     925//Compute the reduced variables r_i = (m_i - <m_i>)/(sig_i/sqrt(Npaq.))
     926//----------------------------------------------
     927int ProcessONOFFReducedMeanData::processCmd() throw(string) {
     928  int rc = 0;
     929  try {
     930    rc = ProcessBase::processCmd();
     931  }
     932  catch (string s) {
     933    throw s;
     934  }
     935
     936  if(debuglev_>0)cout << "Process Data" << endl;
     937  vector<string> modeList;
     938  modeList.push_back("On");
     939  modeList.push_back("Off");
     940  vector<string>::const_iterator iMode;
     941 
     942  uint_4 id;
     943  string tag;
     944
     945  map<string,  TMatrix<r_4> > meanSigmaRedRatioCollect;
     946  map<string,  TMatrix<r_4> >::iterator iMeanSigmaRed, iMeanSigmaRedEnd;
     947
     948  map< pair<string, sa_size_t>, TMatrix<r_4> > reducedRatioCollect;
     949  map< pair<string, sa_size_t>, TMatrix<r_4> >::iterator iReduced, iReducedEnd;
     950
     951  map< sa_size_t, TMatrix<r_4> > meanCollect;
     952  map< sa_size_t, TMatrix<r_4> >::iterator iMean, iMeanEnd;
     953
     954  map< sa_size_t, TMatrix<r_4> > sigmaCollect;
     955  map< sa_size_t, TMatrix<r_4> >::iterator iSigma, iSigmaEnd;
     956 
     957  //
     958  //loop on modes
     959  //
     960  for (iMode = modeList.begin(); iMode != modeList.end(); ++iMode) {
     961    string mode = *iMode;
     962    if(debuglev_>0)cout << "Process Mean-mspec Mode " << mode << endl;
     963
     964
     965    string directoryName;
     966    list<string> listOfMeanFiles;
     967    list<string> listOfSigma2Files;
     968   
     969    if (listOfMeanFiles.size() != listOfSigma2Files.size()) {
     970      throw "ProcessONOFFReducedMeanData: FATAL: length (1) missmatch";
     971    }
     972   
     973    list<string>::const_iterator iFile, iFileEnd;           
     974    //Mean of all Means (one per Mode)
     975    TMatrix<r_4> spectreMean(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ); //implicit init to 0
     976    uint_4 nMeans  = 0;
     977    uint_4 nSigma2s  = 0;
     978    //
     979    //loop on cycles
     980    //
     981    for (sa_size_t icycle = ifirstCycle_; icycle <= ilastCycle_; icycle++) {
     982      directoryName = inputPath_ + "/"
     983        + sourceName_ + "/" + dateOfRun_ + StringToLower(sourceName_) + "/"
     984        + mode + "/";
     985      stringstream sicycle;
     986      sicycle << icycle;
     987      directoryName += spectraDirectory_ + sicycle.str() + "/";
     988
     989      //read directory
     990      listOfMeanFiles   = ListOfFileInDir(directoryName,typeOfFile_); //the means
     991      listOfSigma2Files = ListOfFileInDir(directoryName,sigma2FileName_); //the sigma2s
     992     
     993      //get a mean spectra per file
     994      iFileEnd = listOfMeanFiles.end();
     995      for (iFile = listOfMeanFiles.begin(); iFile != iFileEnd; ++iFile) {
     996//      if(debuglev_>10){
     997//        cout << "read file <" << *iFile << ">:";
     998//      }
     999        FitsInOutFile aSpectrum(*iFile,FitsInOutFile::Fits_RO);
     1000        TMatrix<r_4> spectre(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);
     1001        aSpectrum >> spectre;
     1002//      if(debuglev_>10){
     1003//        cout << " mean of spectre = " << Mean(spectre) << endl;
     1004//      }
     1005        //update mean
     1006        spectreMean   += spectre;
     1007        //save spectrum
     1008        meanCollect.insert( pair< sa_size_t, TMatrix<r_4> >(nMeans,TMatrix<r_4>(spectre,false) ));
     1009
     1010        nMeans++;
     1011      }// end of for files
     1012
     1013      //get a sigma2 spectra per file and store the sigma
     1014      iFileEnd = listOfSigma2Files.end();
     1015      for (iFile = listOfSigma2Files.begin(); iFile != iFileEnd; ++iFile) {
     1016//      if(debuglev_>10){
     1017//        cout << "read file <" << *iFile << ">:";
     1018//      }
     1019        FitsInOutFile aSpectrum(*iFile,FitsInOutFile::Fits_RO);
     1020        TMatrix<r_4> spectre(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);
     1021        aSpectrum >> spectre;
     1022//      if(debuglev_>10){
     1023//        cout << " mean of spectre = " << Mean(Sqrt(spectre)) << endl;
     1024//      }
     1025        //save sigma = Sqrt(sigma2)
     1026        sigmaCollect.insert( pair< sa_size_t, TMatrix<r_4> >(nSigma2s,TMatrix<r_4>(Sqrt(spectre),false) ));
     1027
     1028        nSigma2s++;
     1029      }// end of for files
     1030
     1031    }//end of for cycles
     1032
     1033   
     1034    //finalize mean of means
     1035    if(nMeans>0) {
     1036      spectreMean /= (r_4)nMeans;
     1037//       if(debuglev_>10){
     1038//      cout << "Mean mode [" << mode << "]: " << Mean(spectreMean) <<endl;
     1039//       }
     1040    } else {
     1041      throw "ProcessONOFFReducedMeanData: FATAL: nMeans=0 !!!";
     1042    }
     1043
     1044    if ( nMeans != nSigma2s ) {
     1045      cout << "ProcessONOFFReducedMeanData: nMeans, nSigma2s "
     1046           <<  nMeans << " " <<nSigma2s 
     1047           << endl;
     1048      throw "ProcessONOFFReducedMeanData: FATAL: nMeans <> nSigma2s !!!";
     1049    }
     1050   
     1051    //from here nMeans == nSigma2s
     1052    iMeanEnd   = meanCollect.end();
     1053    iSigmaEnd  = sigmaCollect.end();
     1054    TMatrix<r_4> meanReducedRatio(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);   //implicit init to 0
     1055    TMatrix<r_4> sigma2ReducedRatio(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ); //implicit init to 0
     1056    for (id=0,iMean = meanCollect.begin(), iSigma=sigmaCollect.begin();
     1057         id<nMeans;
     1058         ++id, ++iMean, ++iSigma) {
     1059      TMatrix<r_4> reducedRatio(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);
     1060//       if(debuglev_>10){
     1061//      cout << "reduced Mean [" << mode << "," << id << "]: "
     1062//           << Mean(iMean->second) <<  " "
     1063//           << Mean(spectreMean) <<  " "
     1064//           << Mean(iSigma->second) << " "
     1065//           << sqrt(valNPaqPerWin_)
     1066//           << endl;
     1067//       }
     1068     
     1069      reducedRatio = iMean->second - spectreMean;
     1070      reducedRatio.Div(iSigma->second);
     1071      reducedRatio *= sqrt(valNPaqPerWin_);
     1072
     1073      meanReducedRatio   += reducedRatio;
     1074      sigma2ReducedRatio += reducedRatio&&reducedRatio;
     1075     
     1076//       if(debuglev_>10){
     1077//      cout << "reduced Mean [" << mode << "," << id << "]: " << Mean(reducedRatio) <<endl;
     1078//       }
     1079     
     1080      reducedRatioCollect.insert( pair< pair<string,sa_size_t>, TMatrix<r_4> >(make_pair(mode,id),TMatrix<r_4>(reducedRatio,false) ));
     1081    }
     1082   
     1083    if(debuglev_>10) cout << "number of Means used: " << nMeans << "( =" << nSigma2s << ")" << endl;
     1084    meanReducedRatio   /= (r_4)nMeans;
     1085    sigma2ReducedRatio /= (r_4)nMeans;
     1086    sigma2ReducedRatio -= meanReducedRatio&&meanReducedRatio;
     1087    TMatrix<r_4> sigmaReducedRatio(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);
     1088    sigmaReducedRatio = Sqrt(sigma2ReducedRatio);
     1089   
     1090   
     1091    meanSigmaRedRatioCollect.insert( pair< string,  TMatrix<r_4> >(mode+"Mean",TMatrix<r_4>(meanReducedRatio,false)) );
     1092    meanSigmaRedRatioCollect.insert( pair< string,  TMatrix<r_4> >(mode+"Sigma",TMatrix<r_4>(sigmaReducedRatio,false)) );
     1093
     1094
     1095  }//end loop on mode for raw preocess
     1096
     1097  cout << "Save reduced variable based on Means/Sigmas" << endl;
     1098  string fileName;
     1099 
     1100  fileName = "./reducedMean_" + dateOfRun_ + "_" + StringToLower(sourceName_) + ".ppf";
     1101 
     1102  POutPersist fos(fileName);
     1103  id=0;
     1104  iReducedEnd = reducedRatioCollect.end();
     1105  for (iReduced = reducedRatioCollect.begin();
     1106       iReduced != iReducedEnd ; ++iReduced, ++id) {
     1107    tag = "redMean";
     1108   
     1109    tag += (iReduced->first).first;
     1110    stringstream sid;
     1111    sid << (iReduced->first).second;
     1112    tag += sid.str();
     1113    if(debuglev_>9) {
     1114      cout << "save tag<" << tag << ">" << endl;
     1115    }
     1116    fos << PPFNameTag(tag) << iReduced->second;
     1117  }
     1118 
     1119  iMeanSigmaRedEnd = meanSigmaRedRatioCollect.end();
     1120  for (iMeanSigmaRed = meanSigmaRedRatioCollect.begin();
     1121       iMeanSigmaRed != iMeanSigmaRedEnd;
     1122       ++iMeanSigmaRed){
     1123    tag = "redMean";
     1124    tag += iMeanSigmaRed->first;
     1125    if(debuglev_>9) {
     1126      cout << "save tag<" << tag << ">" << endl;
     1127    }
     1128    fos << PPFNameTag(tag) << iMeanSigmaRed->second;
     1129  }
     1130 
     1131  cout << "OK ProcessONOFFReducedMeanData finished" <<endl;
     1132  return rc;
     1133}
     1134//----------------------------------------------
     1135
     1136//----------------------------------------------
     1137//JEC 8/11/11
     1138//Process the files that are output of the specmfib -act = gain (median of signal files)
     1139//Compute the reduced variables r_i = (med_i - <med_i>)/(med_i/(Ln(2)*sqrt(N))
     1140//----------------------------------------------
     1141int ProcessONOFFReducedMedianData::processCmd() throw(string) {
     1142  int rc = 0;
     1143  try {
     1144    rc = ProcessBase::processCmd();
     1145  }
     1146  catch (string s) {
     1147    throw s;
     1148  }
     1149
     1150  if(debuglev_>0)cout << "Process Data" << endl;
     1151  vector<string> modeList;
     1152  modeList.push_back("On");
     1153  modeList.push_back("Off");
     1154  vector<string>::const_iterator iMode;
     1155 
     1156  uint_4 id;
     1157  string tag;
     1158
     1159  map<string,  TMatrix<r_4> > meanSigmaRedRatioCollect;
     1160  map<string,  TMatrix<r_4> >::iterator iMeanSigmaRed, iMeanSigmaRedEnd;
     1161
     1162  map< pair<string, sa_size_t>, TMatrix<r_4> > reducedRatioCollect;
     1163  map< pair<string, sa_size_t>, TMatrix<r_4> >::iterator iReduced, iReducedEnd;
     1164
     1165  map< sa_size_t, TMatrix<r_4> > medianCollect;
     1166  map< sa_size_t, TMatrix<r_4> >::iterator iMedian, iMedianEnd;
     1167
     1168  map< sa_size_t, TMatrix<r_4> > sigmaCollect;
     1169  map< sa_size_t, TMatrix<r_4> >::iterator iSigma, iSigmaEnd;
     1170 
     1171  //
     1172  //loop on modes
     1173  //
     1174  for (iMode = modeList.begin(); iMode != modeList.end(); ++iMode) {
     1175    string mode = *iMode;
     1176    if(debuglev_>0)cout << "Process Mean-mspec Mode " << mode << endl;
     1177
     1178
     1179    string directoryName;
     1180    list<string> listOfMedianFiles;
     1181   
     1182    list<string>::const_iterator iFile, iFileEnd;           
     1183    //Mean of all Medians (one per Mode)
     1184    TMatrix<r_4> spectreMean(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ); //implicit init to 0
     1185    uint_4 nMedians  = 0;
     1186
     1187    //
     1188    //loop on cycles
     1189    //
     1190    for (sa_size_t icycle = ifirstCycle_; icycle <= ilastCycle_; icycle++) {
     1191      directoryName = inputPath_ + "/"
     1192        + sourceName_ + "/" + dateOfRun_ + StringToLower(sourceName_) + "/"
     1193        + mode + "/";
     1194      stringstream sicycle;
     1195      sicycle << icycle;
     1196      directoryName += spectraDirectory_ + sicycle.str() + "/";
     1197
     1198      //read directory
     1199      listOfMedianFiles   = ListOfFileInDir(directoryName,typeOfFile_); //the means
     1200     
     1201      //get a mean spectra per file
     1202      iFileEnd = listOfMedianFiles.end();
     1203      for (iFile = listOfMedianFiles.begin(); iFile != iFileEnd; ++iFile) {
     1204//      if(debuglev_>10){
     1205//        cout << "read file <" << *iFile << ">:";
     1206//      }
     1207        FitsInOutFile aSpectrum(*iFile,FitsInOutFile::Fits_RO);
     1208        TMatrix<r_4> spectre(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);
     1209        aSpectrum >> spectre;
     1210//      if(debuglev_>10){
     1211//        cout << " mean of spectre = " << Mean(spectre) << endl;
     1212//      }
     1213        //update mean
     1214        spectreMean   += spectre;
     1215        //save spectrum
     1216        medianCollect.insert( pair< sa_size_t, TMatrix<r_4> >(nMedians,TMatrix<r_4>(spectre,false) ));
     1217        sigmaCollect.insert( pair< sa_size_t, TMatrix<r_4> >(nMedians,TMatrix<r_4>(spectre,false) ));
     1218
     1219        nMedians++;
     1220      }// end of for files
     1221
     1222    }//end of for cycles
     1223
     1224   
     1225    //finalize mean of means
     1226    if(nMedians>0) {
     1227      spectreMean /= (r_4)nMedians;
     1228//       if(debuglev_>10){
     1229//      cout << "Mean mode [" << mode << "]: " << Mean(spectreMean) <<endl;
     1230//       }
     1231    } else {
     1232      throw "ProcessONOFFReducedMeanData: FATAL: nMedians=0 !!!";
     1233    }
     1234
     1235    iMedianEnd   = medianCollect.end();
     1236    iSigmaEnd  = sigmaCollect.end();
     1237    TMatrix<r_4> meanReducedRatio(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);   //implicit init to 0
     1238    TMatrix<r_4> sigma2ReducedRatio(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ); //implicit init to 0
     1239    for (id=0,iMedian = medianCollect.begin(), iSigma=sigmaCollect.begin();
     1240         id<nMedians;
     1241         ++id, ++iMedian, ++iSigma) {
     1242      TMatrix<r_4> reducedRatio(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);
     1243//       if(debuglev_>10){
     1244//      cout << "reduced Mean [" << mode << "," << id << "]: "
     1245//           << Mean(iMedian->second) <<  " "
     1246//           << Mean(spectreMean) <<  " "
     1247//           << Mean(iSigma->second) << " "
     1248//           << sqrt(valNPaqPerWin_)
     1249//           << endl;
     1250//       }
     1251     
     1252      reducedRatio = iMedian->second - spectreMean;
     1253      reducedRatio.Div(iSigma->second);
     1254      reducedRatio *= sqrt(valNPaqPerWin_)*log(2.0);
     1255
     1256      meanReducedRatio   += reducedRatio;
     1257      sigma2ReducedRatio += reducedRatio&&reducedRatio;
     1258     
     1259//       if(debuglev_>10){
     1260//      cout << "reduced Mean [" << mode << "," << id << "]: " << Mean(reducedRatio) <<endl;
     1261//       }
     1262     
     1263      reducedRatioCollect.insert( pair< pair<string,sa_size_t>, TMatrix<r_4> >(make_pair(mode,id),TMatrix<r_4>(reducedRatio,false) ));
     1264    }
     1265   
     1266    if(debuglev_>10) cout << "number of Means used: " << nMedians << endl;
     1267    meanReducedRatio   /= (r_4)nMedians;
     1268    sigma2ReducedRatio /= (r_4)nMedians;
     1269    sigma2ReducedRatio -= meanReducedRatio&&meanReducedRatio;
     1270    TMatrix<r_4> sigmaReducedRatio(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);
     1271    sigmaReducedRatio = Sqrt(sigma2ReducedRatio);
     1272   
     1273   
     1274    meanSigmaRedRatioCollect.insert( pair< string,  TMatrix<r_4> >(mode+"Mean",TMatrix<r_4>(meanReducedRatio,false)) );
     1275    meanSigmaRedRatioCollect.insert( pair< string,  TMatrix<r_4> >(mode+"Sigma",TMatrix<r_4>(sigmaReducedRatio,false)) );
     1276
     1277
     1278  }//end loop on mode for raw preocess
     1279
     1280  cout << "Save reduced variable based on Means/Sigmas" << endl;
     1281  string fileName;
     1282 
     1283  fileName = "./reducedMedian_" + dateOfRun_ + "_" + StringToLower(sourceName_) + ".ppf";
     1284 
     1285  POutPersist fos(fileName);
     1286  id=0;
     1287  iReducedEnd = reducedRatioCollect.end();
     1288  for (iReduced = reducedRatioCollect.begin();
     1289       iReduced != iReducedEnd ; ++iReduced, ++id) {
     1290    tag = "redMedian";
     1291   
     1292    tag += (iReduced->first).first;
     1293    stringstream sid;
     1294    sid << (iReduced->first).second;
     1295    tag += sid.str();
     1296    if(debuglev_>9) {
     1297      cout << "save tag<" << tag << ">" << endl;
     1298    }
     1299    fos << PPFNameTag(tag) << iReduced->second;
     1300  }
     1301 
     1302  iMeanSigmaRedEnd = meanSigmaRedRatioCollect.end();
     1303  for (iMeanSigmaRed = meanSigmaRedRatioCollect.begin();
     1304       iMeanSigmaRed != iMeanSigmaRedEnd;
     1305       ++iMeanSigmaRed){
     1306    tag = "redMedian";
     1307    tag += iMeanSigmaRed->first;
     1308    if(debuglev_>9) {
     1309      cout << "save tag<" << tag << ">" << endl;
     1310    }
     1311    fos << PPFNameTag(tag) << iMeanSigmaRed->second;
     1312  }
     1313 
     1314  cout << "OK ProcessONOFFReducedMedianData finished" <<endl;
    7921315  return rc;
    7931316}

  • BAORadio/AmasNancay/trunk/makefile

    r595 r598  
    3838        $(CXXCOMPILE) -c  -o $(OBJ)mergeAnaFiles.o mergeAnaFiles.cc
    3939
    40 anaTools : $(EXE)anaTools
    41         echo '---anaTools made'
    42 
    43 $(EXE)anaTools : $(OBJ)anaTools.o
    44         $(CXXLINK) -o $(EXE)anaTools $(OBJ)anaTools.o  $(SOPHYAALLSLBLIST)
    45 
    46 $(OBJ)anaTools.o : anaTools.cc
    47         $(CXXCOMPILE) -c  -o $(OBJ)anaTools.o anaTools.cc
    48 

  • BAORadio/AmasNancay/trunk/submit2ge-anacalibmaker.sh

    r583 r598  
    4141jobLogName="${jobBatchName}.log.$$"
    4242
    43 qsub -P P_baoradio -l sps=1,irods=1,ct=1200,vmem=1000M,fsize=1000M -o $jobLogName -j yes -N $jobBatchName -m be -M ${LOGNAME}@lal.in2p3.fr -V <<EOF
     43qsub -P P_baoradio -l sps=1,irods=1,ct=1200,vmem=1000M,fsize=1000M -o ${PWD}/${jobLogName} -j yes -N $jobBatchName -m be -M ${LOGNAME}@lal.in2p3.fr -V <<EOF
    4444
    4545${SCRIPTPATH}/anacalibmaker.sh  -src ${sourceRadio} ${dateSelected} ${freqBAOCalib} ${bwBAOCalib} ${simulationMode}

  • BAORadio/AmasNancay/trunk/submit2ge-anagainmaker.sh

    r583 r598  
    3030jobLogName="${jobBatchName}.log.$$"
    3131
    32 qsub -P P_baoradio -l sps=1,irods=1,ct=2400,vmem=1000M,fsize=500M -o $jobLogName -j yes -N $jobBatchName -m be -M ${LOGNAME}@lal.in2p3.fr -V <<EOF
     32qsub -P P_baoradio -l sps=1,irods=1,ct=2400,vmem=1000M,fsize=500M -o ${PWD}/${jobLogName} -j yes -N $jobBatchName -m be -M ${LOGNAME}@lal.in2p3.fr -V <<EOF
    3333
    3434${SCRIPTPATH}/anagainmaker.sh -src ${sourceRadio} ${dateSelected} ${simulationMode}

  • BAORadio/AmasNancay/trunk/submit2ge-anarawonoffmaker.sh

    r595 r598  
    3838jobLogName="${jobBatchName}.log.$$"
    3939
    40 qsub -P P_baoradio -l sps=1,irods=1,ct=2400,vmem=1000M,fsize=1000M -o $jobLogName -j yes -N $jobBatchName -m be -M ${LOGNAME}@lal.in2p3.fr -V <<EOF
     40qsub -P P_baoradio -l sps=1,irods=1,ct=2400,vmem=1000M,fsize=1000M -o ${PWD}/${jobLogName} -j yes -N $jobBatchName -m be -M ${LOGNAME}@lal.in2p3.fr -V <<EOF
    4141
    4242${SCRIPTPATH}/anarawonoffmaker.sh  -src ${sourceRadio} ${dateSelected} ${simulationMode}

  • BAORadio/AmasNancay/trunk/submit2ge-hdrfitsextractor.sh

    r580 r598  
    1919jobLogName="${jobBatchName}.log.$$"
    2020
    21 qsub -P P_baoradio -l sps=1,irods=1,ct=2000,vmem=256M,fsize=700M -o $jobLogName -j yes -N $jobBatchName -m be -M ${LOGNAME}@lal.in2p3.fr -V <<EOF
     21qsub -P P_baoradio -l sps=1,irods=1,ct=2000,vmem=256M,fsize=700M -o ${PWD}/${jobLogName} -j yes -N $jobBatchName -m be -M ${LOGNAME}@lal.in2p3.fr -V <<EOF
    2222
    2323${SCRIPTPATH}/hdrfitsextractor.sh  $sourceRadio

  • BAORadio/AmasNancay/trunk/submit2ge-procspecmfib.sh

    r597 r598  
    125125cpu=`expr ${nFiles} \* 60`
    126126
    127 qsub -P P_baoradio -l sps=1,irods=1,ct=${cpu},vmem=1000M,fsize=${scratchSize}M -o $jobLogName -j yes -N $jobBatchName -m be -M ${LOGNAME}@lal.in2p3.fr -V <<EOF
     127qsub -P P_baoradio -l sps=1,irods=1,ct=${cpu},vmem=1000M,fsize=${scratchSize}M -o ${PWD}/${jobLogName} -j yes -N $jobBatchName -m be -M ${LOGNAME}@lal.in2p3.fr -V <<EOF
    128128
    129129${SCRIPTPATH}/proc_specmfib.sh ${action} -src ${sourceRadio} ${dateSelected} ${typeofproc} ${firstCycle} ${lastCycle} ${simulationMode}

  • BAORadio/AmasNancay/trunk/submit2ge-scaextractor.sh

    r580 r598  
    5959jobLogName="${jobBatchName}.log.$$"
    6060
    61 qsub -P P_baoradio -l sps=1,irods=1,ct=2000,vmem=256M,fsize=700M -o $jobLogName -j yes -N $jobBatchName -m be -M ${LOGNAME}@lal.in2p3.fr -V <<EOF
     61qsub -P P_baoradio -l sps=1,irods=1,ct=2000,vmem=256M,fsize=700M -o ${PWD}/${jobLogName} -j yes -N $jobBatchName -m be -M ${LOGNAME}@lal.in2p3.fr -V <<EOF
    6262
    6363${SCRIPTPATH}/scaextractor.sh  ${sourceRadio} ${force}
Note: See TracChangeset for help on using the changeset viewer.