Changeset 632

Timestamp:

Jan 13, 2012, 12:19:04 PM (12 years ago)

Author:

torrento

Message:

Added: elements to ntuple; treatment of NoGain data; function meanRawDiffOnOffAllClusters

File:

• BAORadio/AmasNancay/trunk/mergeAnaFiles.cc (modified) (16 diffs)

BAORadio/AmasNancay/trunk/mergeAnaFiles.cc

@@ -59 +59 @@
   string calibBandFreq_;  //band of freq. used for calibration
   r_4 rcalibBandFreq_;    //float version
-  int maxNumberCycles_;//maximum number of cycles to be treated
+  int maxNumberCycles_;//(int)maximum number of cycles to be treated
+  string maxNumberCyclesStr_;//(string)maximum number of cycles to be treated
 } para;
 //--------------------------------------------------------------
@@ -237 +238 @@
 }
 //-----
-list<string> ListOfFileInDir(string dir, string filePettern) throw(string) {
+list<string> ListOfFileInDir(string dir, string filePattern) throw(string) {
   list<string> theList;
 
@@ -253 +254 @@
   while ( (dit = readdir(dip)) != NULL ) {
     fileName = dit->d_name;
-    found=fileName.find(filePettern);
+    found=fileName.find(filePattern);
     if (found != string::npos) {
       fullFileName = dir + "/";
@@ -1421 +1422 @@
   TMatrix<r_4> meanOFFovOFF_noCalib(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ); // OFF/Filtered_OFF
   //Tuple only for RAW things to follow 
-  static const int NINFO=11;
+  static const int NINFO=25;
   char* onffTupleName[NINFO]={"cycle"
                               ,"onoffRaw01420","onoffRaw11420"
                               ,"onoffRaw01420side","onoffRaw11420side"
                               ,"onoffRaw0f14101415","onoffRaw1f14101415"
-                              ,"offRaw0f14101415","offRaw1f14101415"
-                              ,"onRaw0f14101415","onRaw1f14101415"
+                              ,"offovoffRaw0f14101415","offovoffRaw1f14101415"
+                              ,"onovoffRaw0f14101415","onovoffRaw1f14101415"
+                              ,"onoffRaw0f14051415","onoffRaw1f14051415"
+                              ,"onoffRaw0f14051410","onoffRaw1f14051410"
+                              ,"diffonoffRaw0f14051415","diffonoffRaw1f14051415"
+                              ,"diffonoffRaw0f14051410","diffonoffRaw1f14051410"
+                              ,"diffonoffRaw0f14101415","diffonoffRaw1f14101415"
+                              ,"onRaw0f14001415","onRaw1f14001415"
+                              ,"offRaw0f14001415","offRaw1f14001415"
   };
   NTuple onoffevolution(NINFO,onffTupleName);
@@ -1444 +1452 @@
   //1400-1420Mhz
   sa_size_t ch1400 = freqToChan(1400);
-  //  sa_size_t ch1405 = freqToChan(1400);
+  sa_size_t ch1405 = freqToChan(1405);
   sa_size_t ch1410 = freqToChan(1410);
   sa_size_t ch1415 = freqToChan(1415);
-  sa_size_t ch1420 = freqToChan(1420);
+  //  sa_size_t ch1420 = freqToChan(1420);
 
   if (para.debuglev_>0){
@@ -1520 +1528 @@
     
     if (para.debuglev_>90){
-      cout << "Liste of cycles common to On & Off: <";
+      cout << dateOfRun << ": Liste of cycles common to On & Off: <";
       for (list<int>::iterator i=commonCycles.begin(); i!=commonCycles.end(); ++i){
         cout << *i << " ";
@@ -1560 +1568 @@
       TMatrix<r_4> aSpecOff(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);
       fin.GetObject(aSpecOff,ppftag);
-      
+
       //Perform the difference ON-OFF
       TMatrix<r_4> diffOnOff_noCalib = aSpecOn - aSpecOff;
@@ -1583 +1591 @@
       offOvOff.Div(aSpecOffFiltered);
       meanOFFovOFF_noCalib += offOvOff;
+
 
       totalNumberCycles++;
@@ -1595 +1604 @@
       xnt[2] = meanInRange_1420Freq_noCalib(1);
 
-
       //Follow up below the 1420.4MHz Freq.
       TVector<r_4> meanInRange_1420aFreq_noCalib(NUMBER_OF_CHANNELS);
@@ -1630 +1638 @@
       xnt[9] = meanInRange_1410a1415Freq_ON_noCalib(0);
       xnt[10] = meanInRange_1410a1415Freq_ON_noCalib(1);
+
+      //AST 09/01/12 follow 1405-1415 and 1405-1410 MHz for (ON-OFF)/OFF_filt
+      TVector<r_4> meanInRange_1405a1415Freq_noCalib(NUMBER_OF_CHANNELS);
+      meanInRange(diffOnOffOvOff_noCalib,ch1405,ch1415,meanInRange_1405a1415Freq_noCalib);
+      xnt[11] = meanInRange_1405a1415Freq_noCalib(0);
+      xnt[12] = meanInRange_1405a1415Freq_noCalib(1);
+      
+      TVector<r_4> meanInRange_1405a1410Freq_noCalib(NUMBER_OF_CHANNELS);
+      meanInRange(diffOnOffOvOff_noCalib,ch1405,ch1410,meanInRange_1405a1410Freq_noCalib);
+      xnt[13] = meanInRange_1405a1410Freq_noCalib(0);
+      xnt[14] = meanInRange_1405a1410Freq_noCalib(1);
+
+      //AST 09/01/12 follow 1405-1415, 1405-1410 and 1410-1415 MHz for (ON-OFF)
+      TVector<r_4> meanInRange_1405a1415Freq_diffOnOff_noCalib(NUMBER_OF_CHANNELS);
+      meanInRange(diffOnOff_noCalib,ch1405,ch1415,meanInRange_1405a1415Freq_diffOnOff_noCalib);
+      xnt[15] = meanInRange_1405a1415Freq_diffOnOff_noCalib(0);
+      xnt[16] = meanInRange_1405a1415Freq_diffOnOff_noCalib(1);
+      
+      TVector<r_4> meanInRange_1405a1410Freq_diffOnOff_noCalib(NUMBER_OF_CHANNELS);
+      meanInRange(diffOnOff_noCalib,ch1405,ch1410,meanInRange_1405a1410Freq_diffOnOff_noCalib);
+      xnt[17] = meanInRange_1405a1410Freq_diffOnOff_noCalib(0);
+      xnt[18] = meanInRange_1405a1410Freq_diffOnOff_noCalib(1);
+
+      TVector<r_4> meanInRange_1410a1415Freq_diffOnOff_noCalib(NUMBER_OF_CHANNELS);
+      meanInRange(diffOnOff_noCalib,ch1410,ch1415,meanInRange_1410a1415Freq_diffOnOff_noCalib);
+      xnt[19] = meanInRange_1410a1415Freq_diffOnOff_noCalib(0);
+      xnt[20] = meanInRange_1410a1415Freq_diffOnOff_noCalib(1);
+
+      //AST 10/01/12 follow 1400-1415 MHz ON, OFF (to continue plot on internal notes)
+      TVector<r_4> meanInRange_1400a1415Freq_On_noCalib(NUMBER_OF_CHANNELS);
+      meanInRange(aSpecOn,ch1400,ch1415,meanInRange_1400a1415Freq_On_noCalib);
+      xnt[21] = meanInRange_1400a1415Freq_On_noCalib(0);
+      xnt[22] = meanInRange_1400a1415Freq_On_noCalib(1);
+      TVector<r_4> meanInRange_1400a1415Freq_Off_noCalib(NUMBER_OF_CHANNELS);
+      meanInRange(aSpecOff,ch1400,ch1415,meanInRange_1400a1415Freq_Off_noCalib);
+      xnt[23] = meanInRange_1400a1415Freq_Off_noCalib(0);
+      xnt[24] = meanInRange_1400a1415Freq_Off_noCalib(1);
       
       //store infos to Ntuple
@@ -1653 +1698 @@
     stringstream tmp;
     tmp << totalNumberCycles;
-    string fileName = para.outPath_+"/rawOnOffDiff_"+StringToLower(para.sourceName_)+"-"+tmp.str()+"Cycles.ppf";
-
+    string fileName;
+    if ( para.ppfFile_ == "dataRawNoGain" ) {
+      fileName = para.outPath_+"/rawOnOffDiffNoGain_"+StringToLower(para.sourceName_)+"-"+tmp.str()+"Cycles.ppf";
+    } else {
+      fileName = para.outPath_+"/rawOnOffDiff_"+StringToLower(para.sourceName_)+"-"+tmp.str()+"Cycles.ppf";
+    }
     POutPersist fos(fileName);
     cout << "Save output in " << fileName << endl;
@@ -1803 +1852 @@
   }
 
-  
 
   {//Save the result
@@ -1828 +1876 @@
     fos << PPFNameTag(tag) << tarr;
   }
+}
+
+// AST 5/12/2011
+//-------------------------------------------------------
+// Take all statistics from cluster data to compute Raw sensitivity vs. t_int
+
+void meanRawDiffOnOffAllClusters() throw(string) {
+
+  // Sources names
+  std::vector < std::string > sources;      //we'll put all of the tokens in here 
+  std::string  temp, insources = para.sourceName_;
+  std::cout << insources << std::endl;
+  while (insources.find("-", 0) != std::string::npos) { //does the string a comma in it?
+    size_t  pos = insources.find("-", 0); //store the position of the delimiter
+    temp = insources.substr(0, pos);      //get the token
+    insources.erase(0, pos + 1);          //erase it from the source 
+    sources.push_back(temp);              //and put it into the array
+  }
+  sources.push_back(insources);           //the last token is all alone
+  std::cout << "Processing " << sources.size() << " sources" << std::endl;
+
+  // Sources max number of cycles
+  std::vector < std::string > cycles;      //we'll put all of the tokens in here
+  std::string  temp1, incycles = para.maxNumberCyclesStr_;
+  while (incycles.find("-", 0) != std::string::npos) { //does the string a hyphen in it?
+    size_t  pos = incycles.find("-", 0);               //store the position of the delimiter
+    temp1 = incycles.substr(0, pos);                   //get the token
+    incycles.erase(0, pos + 1);                        //erase it from the source 
+    cycles.push_back(temp1);               //and put it into the array
+  } 
+  cycles.push_back(incycles);              //the last token is all alone
+
+  if ( sources.size() != cycles.size() ) {
+    std::cout << "Problem: a max number of cycles for each source must be entered" << std::endl;
+    exit (1);
+  } else {
+    for (uint j=0;j<cycles.size();j++) {
+      cout<<"Source "<<sources.at(j)<<": max "<<cycles.at(j)<<" cycles"<< endl;
+    }
+  }
+
+  // Begin analysis for each source
+  list<string> listOfFiles,copylistOfFiles;
+  string directoryName;
+ 
+  //mean ON-OFF over the list of cycles
+  TMatrix<r_4> meanDiffONOFF_noCalib(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);      //(ON-OFF)/GAIN
+  TMatrix<r_4> meanDiffONOFFovOFF_noCalib(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ); //(ON-OFF)/Filtered_OFF
+  TMatrix<r_4> meanONovOFF_noCalib(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);        // ON/Filtered_OFF
+  TMatrix<r_4> meanOFFovOFF_noCalib(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);       // OFF/Filtered_OFF
+  //Tuple only for RAW things to follow
+  static const int NINFO=11;
+  char* onffTupleName[NINFO]={"srcID","date","run","cycletot","cyclesource"
+                              ,"onoffRaw0f14101415","onoffRaw1f14101415"
+                              ,"onoffRaw0f141014125","onoffRaw1f141014125"
+                              ,"onoffRaw0f14101411","onoffRaw1f14101411"
+  };
+  NTuple onoffevolution(NINFO,onffTupleName);
+  r_4 xnt[NINFO];
+
+  //Loop on Sources
+  int numberCyclesSource;  // running number of cycles for each source
+  int maxCyclesSource;     // maximum number of cycles for each source 
+  int totalNumberRuns=0;   //total number of runs (for all clusters)
+  int totalNumberCycles=0; //total number of cycles (for all clusters)
+  std::string dumm;
+
+  for (uint iSource=0;iSource<sources.size();iSource++){
+    directoryName = para.inPath_ + "/" + sources.at(iSource);
+    maxCyclesSource = atoi(cycles.at(iSource).c_str());
+    cout<<"Processing "<< maxCyclesSource<<" cycles for source "<<sources.at(iSource)<<endl;
+
+    //Make the listing of the directory
+    listOfFiles = ListOfFileInDir(directoryName,para.ppfFile_);
+    cout << "Total number of files loaded..." << listOfFiles.size() << endl;
+    //    cout << endl;
+
+    list<string>::const_iterator iFile, iFileEnd, iSpec, iSpecEnd;
+    iFileEnd = listOfFiles.end();
+
+    //1410-1415Mhz
+    sa_size_t ch1410 = freqToChan(1410);
+    sa_size_t ch1411 = freqToChan(1411);
+    sa_size_t ch14125 = freqToChan(1412.5);
+    sa_size_t ch1415 = freqToChan(1415);
+
+    if (para.debuglev_>0){
+      cout << "(1) freq. band for follow up [" <<  ch1410 << ", "<< ch1415 << "]" << endl;
+      cout << "(2) freq. band for follow up [" <<  ch1410 << ", "<< ch14125 << "]" << endl;
+      cout << "(3) freq. band for follow up [" <<  ch1410 << ", "<< ch1411 << "]" << endl;
+    }
+
+    numberCyclesSource=0; //total number of cycles in this source
+    //Loop on files
+    for (iFile = listOfFiles.begin(); iFile != iFileEnd; ++iFile) {
+      if (para.debuglev_>90){
+        cout << "load file <" << *iFile << ">" << endl;
+      }
+      
+      vector<string> tokens;
+      split(*iFile,"_",tokens);
+      string dateOfRun = tokens[1];
+      r_4 rdateOfRun = atoi(dateOfRun.c_str()) - 2011*1.e4;
+      if (para.debuglev_>90){
+        cout << "date <" << dateOfRun << ">" << endl;
+        cout << "rdate <" << rdateOfRun << ">" << endl;
+      }
+      vector<string> tokens2;
+      split(tokens[2],".",tokens2);
+      string srcLower = tokens2[0], csrcID;
+      size_t found;
+      found=srcLower.find_first_of("0123456789");
+      while (found!=string::npos){
+        csrcID.push_back(srcLower[found]);
+        found=srcLower.find_first_of("0123456789",found+1);
+      }
+      r_4 srcID = atoi(csrcID.c_str());
+      if (para.debuglev_>1000){
+        cout << "source <" << srcLower << ">" << endl;
+        cout << "srcID <" << srcID << ">" << endl;
+      }
+
+      // AST remove bad data of Abell1205
+      if ( (dateOfRun=="20110822"||dateOfRun=="20110903"||dateOfRun=="20110905"||dateOfRun=="20110912"||dateOfRun=="20110915") && srcLower == "abell1205" ) {
+        cout << ">>> Skipping bad period " << dateOfRun << endl;
+        continue;
+      }
+
+      PInPersist fin(*iFile);
+      vector<string> vec = fin.GetNameTags();
+
+      vector<string> modeList;
+      modeList.push_back("On");
+      modeList.push_back("Off");
+      vector<string>::const_iterator iMode;
+    
+      map<string, list<int> > cycleModeCollect;
+
+      for (iMode = modeList.begin(); iMode!=modeList.end(); ++iMode) {
+        list<string> listOfSpectra;
+        //Keep only required PPF objects
+        string matchstr = "specRaw"+(*iMode)+"[0-9]+"; 
+        std::remove_copy_if(
+                            vec.begin(), vec.end(), back_inserter(listOfSpectra),
+                            not1(StringMatch(matchstr))
+                            );
+      
+        listOfSpectra.sort(stringCompare);
+        iSpecEnd = listOfSpectra.end();
+        
+        matchstr = "[0-9]+";
+        //Loop of spectra matrix
+        list<int> listOfCycles;
+        for (iSpec = listOfSpectra.begin(); iSpec!=iSpecEnd;  ++iSpec){
+          int b,e;
+          regexp(iSpec->c_str(),matchstr.c_str(),&b,&e);
+          if (para.debuglev_>90){
+            cout << " spectra <" << *iSpec << ">" << endl;
+            cout << " cycle " << iSpec->substr(b) << endl;
+          }
+          listOfCycles.push_back(atoi((iSpec->substr(b)).c_str()));
+        }//end loop spectra
+        cycleModeCollect[*iMode] = listOfCycles;
+      }//end of mode
+
+      //Take the Intersection of the list Of cycles in mode Off and On    
+      list<int> commonCycles;
+      set_intersection(cycleModeCollect["On"].begin(),
+                       cycleModeCollect["On"].end(),
+                       cycleModeCollect["Off"].begin(),
+                       cycleModeCollect["Off"].end(),
+                       back_inserter(commonCycles)
+                       );
+      
+      if (para.debuglev_>90){
+        cout << "Liste of cycles common to On & Off: <";
+        for (list<int>::iterator i=commonCycles.begin(); i!=commonCycles.end(); ++i){
+          cout << *i << " ";
+        }
+        cout << ">" << endl;
+      }
+      
+      //Loop on cycles
+      for (list<int>::iterator ic=commonCycles.begin(); ic!=commonCycles.end(); ++ic){
+        
+        // AST 28.11.11 remove non-calibrated cycles for Abell1205 and Abell2440
+        if ( *ic == 1 && srcLower == "abell1205" ) {
+          if ( dateOfRun == "20110502" || dateOfRun == "20110607" || dateOfRun == "20110818" ) {
+            cout << "Skipping non-calibrated cycle " << *ic << endl;
+            continue;
+          }
+        } else if ( *ic == 1 && srcLower == "abell2440" ) {
+          if ( dateOfRun == "20110516" ) {
+            cout << "Skipping non-calibrated cycle " << *ic << endl;
+            continue;
+          }
+        } else if ( *ic == 3 && srcLower == "abell1205" ) {
+          if ( dateOfRun == "20110810" ) {
+            cout << "Skipping non-calibrated cycle " << *ic << endl;
+            continue;
+          }
+        }
+        
+        string ppftag;
+        //load ON phase
+        stringstream cycle;
+        cycle << *ic;
+        ppftag = "specRawOn"+cycle.str();
+        TMatrix<r_4> aSpecOn(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);
+        fin.GetObject(aSpecOn,ppftag);
+        
+        //Load OFF phase
+        ppftag = "specRawOff"+cycle.str();
+        TMatrix<r_4> aSpecOff(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);
+        fin.GetObject(aSpecOff,ppftag);
+        
+      //Perform the difference ON-OFF
+        TMatrix<r_4> diffOnOff_noCalib = aSpecOn - aSpecOff;
+        meanDiffONOFF_noCalib += diffOnOff_noCalib;
+        
+        TMatrix<r_4> diffOnOffOvOff_noCalib(diffOnOff_noCalib,false); //do not share data
+        TMatrix<r_4> aSpecOffFiltered(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);
+        sa_size_t halfWidth = 35; //number of freq. bin for the 1/2 width of the filtering window
+        medianFiltering(aSpecOff,halfWidth,aSpecOffFiltered);
+        diffOnOffOvOff_noCalib.Div(aSpecOffFiltered); //division in place
+        meanDiffONOFFovOFF_noCalib += diffOnOffOvOff_noCalib;
+        
+        numberCyclesSource++;
+        totalNumberCycles++;
+        
+        //Fill NTuple
+        xnt[0] = srcID;
+        xnt[1] = rdateOfRun;
+        xnt[2] = totalNumberRuns;
+        xnt[3] = totalNumberCycles;
+        xnt[4] = numberCyclesSource;
+        
+        TVector<r_4> meanInRange_1410a1415Freq_noCalib(NUMBER_OF_CHANNELS);
+        meanInRange(diffOnOffOvOff_noCalib,ch1410,ch1415,meanInRange_1410a1415Freq_noCalib);
+        TVector<r_4> meanInRange_1410a14125Freq_noCalib(NUMBER_OF_CHANNELS);
+        meanInRange(diffOnOffOvOff_noCalib,ch1410,ch14125,meanInRange_1410a14125Freq_noCalib);
+        TVector<r_4> meanInRange_1410a1411Freq_noCalib(NUMBER_OF_CHANNELS);
+        meanInRange(diffOnOffOvOff_noCalib,ch1410,ch1411,meanInRange_1410a1411Freq_noCalib);
+
+        xnt[5] = meanInRange_1410a1415Freq_noCalib(0);
+        xnt[6] = meanInRange_1410a1415Freq_noCalib(1);
+        xnt[7] = meanInRange_1410a14125Freq_noCalib(0);
+        xnt[8] = meanInRange_1410a14125Freq_noCalib(1);
+        xnt[9] = meanInRange_1410a1411Freq_noCalib(0);
+        xnt[10] = meanInRange_1410a1411Freq_noCalib(1);
+
+        onoffevolution.Fill(xnt);
+        
+        //Quit if enough
+        if ( numberCyclesSource >= maxCyclesSource ) break;
+      }//end loop on cycles
+      if (numberCyclesSource >= maxCyclesSource) break;
+      totalNumberRuns++;
+    } // end loop on files
+    cout << "End of job on source: we have treated " << numberCyclesSource << " cycles" << endl;        
+  } // end loop on sources
+  cout << "End of jobs: we have treated " << totalNumberCycles << " cycles" << endl;
+
+  //Normalisation
+  if(totalNumberCycles > 0){
+    meanDiffONOFF_noCalib       /= (r_4)totalNumberCycles;
+    meanDiffONOFFovOFF_noCalib  /= (r_4)totalNumberCycles;
+  }
+
+  {//save results
+    stringstream tmp;
+    tmp << totalNumberCycles;
+    string fileName = para.outPath_+"/rawOnOffDiff_"+StringToLower(para.sourceName_)+"-"+tmp.str()+"Cycles.ppf";
+
+    POutPersist fos(fileName);
+    cout << "Save output in " << fileName << endl;
+
+    string tag = "meanNoCalib";
+    fos << PPFNameTag(tag) << meanDiffONOFF_noCalib;
+    
+    tag = "meanOvOffNoCalib";
+    fos << PPFNameTag(tag) << meanDiffONOFFovOFF_noCalib;
+
+    tag = "onoffevol";
+    fos << PPFNameTag(tag) << onoffevolution;
+
+  }//end save
+ 
 }
 
@@ -1856 +2192 @@
     if (strcmp(arg[ka],"-h")==0) {
       cout << "Usage:  -act [meanRawDiffOnOff]|medianRawDiffOnOff|meanCalibBAODiffOnOff\n"
-           << " -mxcycles <number> (max. number of cycles to be treated)\n"
+           << " -mxcycles <number> (max. number of cycles to be treated) or <list-of-numbers> if several sources are analysed \n"
            << " -calibfreq <number> (cf. freq. used by calibration operation)\n"
            << " -calibbandfreq <number> (cf. band of freq. used by calibration operation)\n"
-           << " -src [Abell85]\n" 
+           << " -src [Abell85] or list of sources [Abell85-Abell1205-Abell2440]\n" 
            << " -inPath [.]|<top_root_dir of the ppf file>\n" 
            << " (ex. /sps/baoradio/AmasNancay/JEC/\n " 
@@ -1924 +2260 @@
   if (mxcycles != "") {
     para.maxNumberCycles_ = atoi(mxcycles.c_str());
+    para.maxNumberCyclesStr_ = mxcycles;
   } else {
     para.maxNumberCycles_ = std::numeric_limits<int>::max();
+    para.maxNumberCyclesStr_ = "500";
   }
 
@@ -1962 +2300 @@
     } else if (action == "calibCoeffNtp") {
       calibCoeffNtp();
+    } else if (action == "meanRawDiffOnOffAllClusters") {
+      meanRawDiffOnOffAllClusters();
     } else {
       msg = "Unknown action " + action;