Changeset 615 for BAORadio

Timestamp:

Nov 28, 2011, 3:47:28 PM (13 years ago)

Author:

torrento

Message:

Added calibcoeffNtp function

File:

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

BAORadio/AmasNancay/trunk/mergeAnaFiles.cc

@@ -63 +63 @@
 //--------------------------------------------------------------
 //Utility functions
-// template <class T>
-// void MeanSigma(TArray<T> const & a, T & mean, T & sig){
-//   if (a.NbDimensions() < 1) 
-//     throw RangeCheckError("MathArray<T>::MeanSigma(TArray<T> const & a..) Not Allocated Array a !");
-//   const T * pe;
-//   sa_size_t j,k;
-//   mean=0.;
-//   sig = 0.;
-//   T valok;
-//   if (a.AvgStep() > 0)   {  // regularly spaced elements
-//     sa_size_t step = a.AvgStep(); 
-//     sa_size_t maxx = a.Size()*step;
-//     pe = a.Data();
-//     for(k=0; k<maxx; k+=step )  { 
-//       valok = pe[k]; 
-//       mean += valok;  sig += valok*valok;
-//     }
-//   }
-//   else {    // Non regular data spacing ...
-//     int_4 ka = a.MaxSizeKA();
-//     sa_size_t step = a.Step(ka);
-//     sa_size_t gpas = a.Size(ka)*step;
-//     sa_size_t naxa = a.Size()/a.Size(ka);
-//     for(j=0; j<naxa; j++)  {
-//       pe = a.DataBlock().Begin()+a.Offset(ka,j);
-//       for(k=0; k<gpas; k+=step) { 
-//      valok = pe[k]; 
-//      mean += valok;  sig += valok*valok;
-//       }
-//     }
-//   }
-//   T dsz = (T)(a.Size());
-//   mean /= dsz;
-//   if (dsz > 1.5) {
-//     sig = sig/dsz - mean*mean;
-//     sig *= (dsz/(dsz-1));
-//     if (sig >= 0.) sig = sqrt(sig);
-//   }
-//   else sig = 0.;
-// }
-
 
 sa_size_t freqToChan(r_4 f){
@@ -120 +79 @@
   for (sa_size_t ir=0; ir<nr; ir++){
     TVector<r_4> tmp(mtx(Range(ir),Range(chLow,chHigh)),false);
-    r_8 mean, sigma;
+    double mean, sigma;
     MeanSigma(tmp,mean,sigma);
     vec(ir) = mean;
@@ -314 +273 @@
   }//eo loop on slices
 }
+//AST 18/11/11 
+//-------------------------------------------------------
+//Make n-tuple with calibration coefficients 
+//
+void calibCoeffNtp() throw(string) {
+  list<string> listOfFiles;
+  string directoryName;
+  directoryName = para.inPath_ + "/" + para.sourceName_;
+
+  //Make the listing of the directory
+  listOfFiles = ListOfFileInDir(directoryName,para.ppfFile_);
+  
+  list<string>::const_iterator iFile, iFileEnd, iSpec, iSpecEnd;
+  iFileEnd = listOfFiles.end();
+
+  static const int NINFO=19;
+  char* calibTupleName[NINFO]={"date","run","cycle"
+                              ,"meancoeffoff0","meancoeffoff1"
+                              ,"meancoeffon0","meancoeffon1"
+                              ,"coeffoff0","coeffoff1"                              
+                              ,"coeffon0","coeffon1"
+                              ,"invmeancoeffoff0","invmeancoeffoff1"
+                              ,"invmeancoeffon0","invmeancoeffon1"
+                              ,"invcoeffoff0","invcoeffoff1"                              
+                              ,"invcoeffon0","invcoeffon1"
+  };
+  NTuple calibcoeffs(NINFO,calibTupleName); 
+  r_4 xnt[NINFO];
+  
+  int totalNumberRuns=0;   //total number of runs
+  int totalNumberCycles=0; //total number of cycles for normalisation
+  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];
+    
+    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;
+    }
+    
+    //
+    //Load BAO Calibration factors "per Cycle and Channels"
+    //Compute the mean per Cycle to 
+    //  fill factors "per Run and Channels" with the same cycle list
+    //
+    //
+    //TODO improve the code....
+
+    TMatrix<r_4> calibBAOfactors_Off_Cycle_Ch0;
+    TMatrix<r_4> calibBAOfactors_Off_Cycle_Ch1;
+    TMatrix<r_4> calibBAOfactors_On_Cycle_Ch0;
+    TMatrix<r_4> calibBAOfactors_On_Cycle_Ch1;
+    
+    string calibFileName;
+    ifstream ifs;
+    sa_size_t nr,nc; //values read 
+
+    //OFF Cycle per Channel
+    calibFileName = directoryName + "/" 
+      + "calib_" + dateOfRun + "_" + srcLower + "_Off_" 
+      + para.calibFreq_ +"MHz-Ch0Cycles.txt";
+    if(para.debuglev_>0) cout << "Read Calib file " << calibFileName << endl;
+    ifs.open(calibFileName.c_str());
+    if ( ! ifs.is_open() ) {
+      throw calibFileName + " cannot be opened...";
+    }   
+    calibBAOfactors_Off_Cycle_Ch0.ReadASCII(ifs,nr,nc);
+    if(para.debuglev_>9){
+      cout << "(nr,nc): "<< nr << "," << nc << endl;
+      calibBAOfactors_Off_Cycle_Ch0.Print(cout);
+      cout << endl;
+    }
+
+    TMatrix<r_4> calibBAOfactors_Off_Run_Ch0(nr,nc);
+    calibBAOfactors_Off_Run_Ch0.Column(0) = calibBAOfactors_Off_Cycle_Ch0.Column(0);
+    {//Compute the mean
+      TVector<r_4> coef(calibBAOfactors_Off_Cycle_Ch0(Range::all(),Range::last()),false);
+      double mean,sigma;
+      MeanSigma(coef,mean,sigma);
+      calibBAOfactors_Off_Run_Ch0.Column(1).SetCst(mean);
+    }
+    if(para.debuglev_>9){
+      cout << "Fill calib. with mean value " << endl; 
+      calibBAOfactors_Off_Run_Ch0.Print(cout);
+      cout << endl;
+    }
+
+    TMatrix<r_4> invcalibBAOfactors_Off_Cycle_Ch0(nr,nc);
+    invcalibBAOfactors_Off_Cycle_Ch0.Column(0) = calibBAOfactors_Off_Cycle_Ch0.Column(0);
+    {//Compute the inverse value
+      TVector<r_4> coef(calibBAOfactors_Off_Cycle_Ch0(Range::all(),Range::last()),false);
+      invcalibBAOfactors_Off_Cycle_Ch0.Column(1).SetCst(1);
+      invcalibBAOfactors_Off_Cycle_Ch0.Column(1).Div(coef);
+      if(para.debuglev_>9){
+        cout << "Fill calib. with inverse value " << endl;
+        invcalibBAOfactors_Off_Cycle_Ch0.Print(cout);
+        cout << endl;
+      }
+    }
+    TMatrix<r_4> invcalibBAOfactors_Off_Run_Ch0(nr,nc);
+    invcalibBAOfactors_Off_Run_Ch0.Column(0) = calibBAOfactors_Off_Cycle_Ch0.Column(0);
+    {//Compute the inverse mean
+      TVector<r_4> coef(invcalibBAOfactors_Off_Cycle_Ch0(Range::all(),Range::last()),false);
+      double mean,sigma;
+      MeanSigma(coef,mean,sigma);
+      invcalibBAOfactors_Off_Run_Ch0.Column(1).SetCst(mean);
+    }
+    if(para.debuglev_>9){
+      cout << "Fill calib. with inverse mean value " << endl; 
+      invcalibBAOfactors_Off_Run_Ch0.Print(cout);
+      cout << endl;
+    }
+    ifs.close();
+
+    //
+    calibFileName = directoryName + "/" 
+      + "calib_" + dateOfRun + "_" + srcLower + "_Off_" 
+      + para.calibFreq_ +"MHz-Ch1Cycles.txt";
+    if(para.debuglev_>0) cout << "Read Calib file " << calibFileName << endl;
+    ifs.open(calibFileName.c_str());
+    if ( ! ifs.is_open() ) {
+
+      throw calibFileName + " cannot be opened...";
+    }   
+    calibBAOfactors_Off_Cycle_Ch1.ReadASCII(ifs,nr,nc);
+    if(para.debuglev_>9){
+      cout << "(nr,nc): "<< nr << "," << nc << endl;
+      calibBAOfactors_Off_Cycle_Ch1.Print(cout);
+      cout << endl;
+    }
+    TMatrix<r_4> calibBAOfactors_Off_Run_Ch1(nr,nc);
+    calibBAOfactors_Off_Run_Ch1.Column(0) = calibBAOfactors_Off_Cycle_Ch1.Column(0);
+    {//Compute the mean
+      TVector<r_4> coef(calibBAOfactors_Off_Cycle_Ch1(Range::all(),Range::last()),false);
+      double mean,sigma;
+      MeanSigma(coef,mean,sigma);
+      //      cout << "Mean: " << mean << " sigma " << sigma << endl;
+      calibBAOfactors_Off_Run_Ch1.Column(1).SetCst(mean);
+    }
+    if(para.debuglev_>9){
+      cout << "Fill calib. with mean value " << endl; 
+      calibBAOfactors_Off_Run_Ch1.Print(cout);
+      cout << endl;
+    }
+    TMatrix<r_4> invcalibBAOfactors_Off_Cycle_Ch1(nr,nc);
+    invcalibBAOfactors_Off_Cycle_Ch1.Column(0) = calibBAOfactors_Off_Cycle_Ch1.Column(0);
+    {//Compute the inverse value
+      TVector<r_4> coef(calibBAOfactors_Off_Cycle_Ch1(Range::all(),Range::last()),false);
+      invcalibBAOfactors_Off_Cycle_Ch1.Column(1).SetCst(1);
+      invcalibBAOfactors_Off_Cycle_Ch1.Column(1).Div(coef);
+      if(para.debuglev_>9){
+        cout << "Fill calib. with inverse value " << endl;
+        invcalibBAOfactors_Off_Cycle_Ch1.Print(cout);
+        cout << endl;
+      }
+    }
+    TMatrix<r_4> invcalibBAOfactors_Off_Run_Ch1(nr,nc);
+    invcalibBAOfactors_Off_Run_Ch1.Column(0) = calibBAOfactors_Off_Cycle_Ch1.Column(0);
+    {//Compute the inverse mean
+      TVector<r_4> coef(invcalibBAOfactors_Off_Cycle_Ch1(Range::all(),Range::last()),false);
+      double mean,sigma;
+      MeanSigma(coef,mean,sigma);
+      invcalibBAOfactors_Off_Run_Ch1.Column(1).SetCst(mean);
+    }
+    if(para.debuglev_>9){
+      cout << "Fill calib. with inverse mean value " << endl; 
+      invcalibBAOfactors_Off_Run_Ch1.Print(cout);
+      cout << endl;
+    }
+    ifs.close();
+
+    //ON Cycle per Channel
+    calibFileName = directoryName + "/" 
+      + "calib_" + dateOfRun + "_" + srcLower + "_On_" 
+      + para.calibFreq_ +"MHz-Ch0Cycles.txt";
+    if(para.debuglev_>0) cout << "Read Calib file " << calibFileName << endl;
+    ifs.open(calibFileName.c_str());
+    if ( ! ifs.is_open() ) {
+
+      throw calibFileName + " cannot be opened...";
+    }   
+    calibBAOfactors_On_Cycle_Ch0.ReadASCII(ifs,nr,nc);
+    if(para.debuglev_>9){
+      cout << "(nr,nc): "<< nr << "," << nc << endl;
+      calibBAOfactors_On_Cycle_Ch0.Print(cout);
+      cout << endl;      
+    }
+
+    TMatrix<r_4> calibBAOfactors_On_Run_Ch0(nr,nc);
+    calibBAOfactors_On_Run_Ch0.Column(0) = calibBAOfactors_On_Cycle_Ch0.Column(0);
+    {//Compute the mean
+      TVector<r_4> coef(calibBAOfactors_On_Cycle_Ch0(Range::all(),Range::last()),false);
+      double mean,sigma;
+      MeanSigma(coef,mean,sigma);
+      //      cout << "Mean: " << mean << " sigma " << sigma << endl;
+      calibBAOfactors_On_Run_Ch0.Column(1).SetCst(mean);
+    }
+    if(para.debuglev_>9){
+      cout << "Fill calib. with mean value " << endl; 
+      calibBAOfactors_On_Run_Ch0.Print(cout);
+      cout << endl;
+    }
+
+    TMatrix<r_4> invcalibBAOfactors_On_Cycle_Ch0(nr,nc);
+    invcalibBAOfactors_On_Cycle_Ch0.Column(0) = calibBAOfactors_On_Cycle_Ch0.Column(0);
+    {//Compute the inverse value
+      TVector<r_4> coef(calibBAOfactors_On_Cycle_Ch0(Range::all(),Range::last()),false);
+      invcalibBAOfactors_On_Cycle_Ch0.Column(1).SetCst(1);
+      invcalibBAOfactors_On_Cycle_Ch0.Column(1).Div(coef);
+      if(para.debuglev_>9){
+        cout << "Fill calib. with inverse value " << endl;
+        invcalibBAOfactors_On_Cycle_Ch0.Print(cout);
+        cout << endl;
+      }
+    }
+    TMatrix<r_4> invcalibBAOfactors_On_Run_Ch0(nr,nc);
+    invcalibBAOfactors_On_Run_Ch0.Column(0) = calibBAOfactors_On_Cycle_Ch0.Column(0);
+    {//Compute the inverse mean
+      TVector<r_4> coef(invcalibBAOfactors_On_Cycle_Ch0(Range::all(),Range::last()),false);
+      double mean,sigma;
+      MeanSigma(coef,mean,sigma);
+      invcalibBAOfactors_On_Run_Ch0.Column(1).SetCst(mean);
+    }
+    if(para.debuglev_>9){
+      cout << "Fill calib. with inverse mean value " << endl; 
+      invcalibBAOfactors_On_Run_Ch0.Print(cout);
+      cout << endl;
+    }
+    ifs.close();
+    
+    calibFileName = directoryName + "/" 
+      + "calib_" + dateOfRun + "_" + srcLower + "_On_" 
+      + para.calibFreq_ +"MHz-Ch1Cycles.txt";
+    if(para.debuglev_>0) cout << "Read Calib file " << calibFileName << endl;
+    ifs.open(calibFileName.c_str());
+    if ( ! ifs.is_open() ) {
+      throw calibFileName + " cannot be opened...";
+    }   
+    calibBAOfactors_On_Cycle_Ch1.ReadASCII(ifs,nr,nc);
+    if(para.debuglev_>9){
+      cout << "(nr,nc): "<< nr << "," << nc << endl;
+      calibBAOfactors_On_Cycle_Ch1.Print(cout);
+      cout << endl;
+    }
+    TMatrix<r_4> calibBAOfactors_On_Run_Ch1(nr,nc);
+    calibBAOfactors_On_Run_Ch1.Column(0) = calibBAOfactors_On_Cycle_Ch1.Column(0);
+    {//Compute the mean
+      TVector<r_4> coef(calibBAOfactors_On_Cycle_Ch1(Range::all(),Range::last()),false);
+      double mean,sigma;
+      MeanSigma(coef,mean,sigma);
+      //      cout << "Mean: " << mean << " sigma " << sigma << endl;
+      calibBAOfactors_On_Run_Ch1.Column(1).SetCst(mean);
+    }
+    if(para.debuglev_>9){
+      cout << "Fill calib. with mean value " << endl; 
+      calibBAOfactors_On_Run_Ch1.Print(cout);
+      cout << endl;
+    }
+
+    TMatrix<r_4> invcalibBAOfactors_On_Cycle_Ch1(nr,nc);
+    invcalibBAOfactors_On_Cycle_Ch1.Column(0) = calibBAOfactors_On_Cycle_Ch1.Column(0);
+    {//Compute the inverse value
+      TVector<r_4> coef(calibBAOfactors_On_Cycle_Ch1(Range::all(),Range::last()),false);
+      invcalibBAOfactors_On_Cycle_Ch1.Column(1).SetCst(1);
+      invcalibBAOfactors_On_Cycle_Ch1.Column(1).Div(coef);
+      if(para.debuglev_>9){
+        cout << "Fill calib. with inverse value " << endl;
+        invcalibBAOfactors_On_Cycle_Ch1.Print(cout);
+        cout << endl;
+      }
+    }
+    TMatrix<r_4> invcalibBAOfactors_On_Run_Ch1(nr,nc);
+    invcalibBAOfactors_On_Run_Ch1.Column(0) = calibBAOfactors_On_Cycle_Ch1.Column(0);
+    {//Compute the inverse mean
+      TVector<r_4> coef(invcalibBAOfactors_On_Cycle_Ch1(Range::all(),Range::last()),false);
+      double mean,sigma;
+      MeanSigma(coef,mean,sigma);
+      invcalibBAOfactors_On_Run_Ch1.Column(1).SetCst(mean);
+    }
+    if(para.debuglev_>9){
+      cout << "Fill calib. with inverse mean value " << endl; 
+      invcalibBAOfactors_On_Run_Ch1.Print(cout);
+      cout << endl;
+    }
+    ifs.close();
+
+    //link <cycle> - <calibration coefficient>
+    //We cannot rely on identical cycle list of the OFF and ON calibration
+    map<int,r_4> calibBAO_Off_Run_Ch0;
+    map<int,r_4> calibBAO_Off_Run_Ch1;
+    map<int,r_4> calibBAO_On_Run_Ch0;
+    map<int,r_4> calibBAO_On_Run_Ch1;
+
+    map<int,r_4> calibBAO_Off_Cycle_Ch0;
+    map<int,r_4> calibBAO_Off_Cycle_Ch1;
+    map<int,r_4> calibBAO_On_Cycle_Ch0;
+    map<int,r_4> calibBAO_On_Cycle_Ch1;
+
+    map<int,r_4> invcalibBAO_Off_Run_Ch0;
+    map<int,r_4> invcalibBAO_Off_Run_Ch1;
+    map<int,r_4> invcalibBAO_On_Run_Ch0;
+    map<int,r_4> invcalibBAO_On_Run_Ch1;
+
+    map<int,r_4> invcalibBAO_Off_Cycle_Ch0;
+    map<int,r_4> invcalibBAO_Off_Cycle_Ch1;
+    map<int,r_4> invcalibBAO_On_Cycle_Ch0;
+    map<int,r_4> invcalibBAO_On_Cycle_Ch1;
+
+    //per Run based BAO coefficients
+    nr = calibBAOfactors_Off_Run_Ch0.NRows();
+    for (sa_size_t ir=0; ir<nr; ++ir){
+//       cout << "Calib. Off Run Ch0 cycle ["<< calibBAOfactors_Off_Run_Ch0(ir,0)<<"], val "
+//         << calibBAOfactors_Off_Run_Ch0(ir,1) << endl;
+
+      calibBAO_Off_Run_Ch0[(int)calibBAOfactors_Off_Run_Ch0(ir,0)]
+        = calibBAOfactors_Off_Run_Ch0(ir,1);
+      calibBAO_Off_Cycle_Ch0[(int)calibBAOfactors_Off_Cycle_Ch0(ir,0)]
+        = calibBAOfactors_Off_Cycle_Ch0(ir,1);
+      calibBAO_Off_Run_Ch1[(int)calibBAOfactors_Off_Run_Ch1(ir,0)]
+        = calibBAOfactors_Off_Run_Ch1(ir,1);
+      calibBAO_Off_Cycle_Ch1[(int)calibBAOfactors_Off_Cycle_Ch1(ir,0)]
+        = calibBAOfactors_Off_Cycle_Ch1(ir,1);
+
+      invcalibBAO_Off_Run_Ch0[(int)invcalibBAOfactors_Off_Run_Ch0(ir,0)]
+        = invcalibBAOfactors_Off_Run_Ch0(ir,1);
+      invcalibBAO_Off_Cycle_Ch0[(int)invcalibBAOfactors_Off_Cycle_Ch0(ir,0)]
+        = invcalibBAOfactors_Off_Cycle_Ch0(ir,1);
+      invcalibBAO_Off_Run_Ch1[(int)invcalibBAOfactors_Off_Run_Ch1(ir,0)]
+        = invcalibBAOfactors_Off_Run_Ch1(ir,1);
+      invcalibBAO_Off_Cycle_Ch1[(int)invcalibBAOfactors_Off_Cycle_Ch1(ir,0)]
+        = invcalibBAOfactors_Off_Cycle_Ch1(ir,1);
+    }//eo loop on coef Off
+    
+    nr = calibBAOfactors_On_Run_Ch0.NRows();
+    for (sa_size_t ir=0; ir<nr; ++ir){
+      calibBAO_On_Run_Ch0[(int)calibBAOfactors_On_Run_Ch0(ir,0)]
+        = calibBAOfactors_On_Run_Ch0(ir,1);
+      calibBAO_On_Cycle_Ch0[(int)calibBAOfactors_On_Cycle_Ch0(ir,0)]
+        = calibBAOfactors_On_Cycle_Ch0(ir,1);
+      calibBAO_On_Run_Ch1[(int)calibBAOfactors_On_Run_Ch1(ir,0)]
+        = calibBAOfactors_On_Run_Ch1(ir,1);
+      calibBAO_On_Cycle_Ch1[(int)calibBAOfactors_On_Cycle_Ch1(ir,0)]
+        = calibBAOfactors_On_Cycle_Ch1(ir,1);
+
+      invcalibBAO_On_Run_Ch0[(int)invcalibBAOfactors_On_Run_Ch0(ir,0)]
+        = invcalibBAOfactors_On_Run_Ch0(ir,1);
+      invcalibBAO_On_Cycle_Ch0[(int)invcalibBAOfactors_On_Cycle_Ch0(ir,0)]
+        = invcalibBAOfactors_On_Cycle_Ch0(ir,1);
+      invcalibBAO_On_Run_Ch1[(int)invcalibBAOfactors_On_Run_Ch1(ir,0)]
+        = invcalibBAOfactors_On_Run_Ch1(ir,1);
+      invcalibBAO_On_Cycle_Ch1[(int)invcalibBAOfactors_On_Cycle_Ch1(ir,0)]
+        = invcalibBAOfactors_On_Cycle_Ch1(ir,1);
+    }//eo loop on coef On
+      
+    //Loop on cyles
+    for (list<int>::iterator ic=commonCycles.begin(); ic!=commonCycles.end(); ++ic){
+
+      //Look if the cycle has been calibrated...
+      bool isCycleCalibrated = 
+        ( calibBAO_On_Run_Ch0.count(*ic)    *
+          calibBAO_On_Run_Ch1.count(*ic)    *
+          calibBAO_Off_Run_Ch0.count(*ic)   *
+          calibBAO_Off_Run_Ch1.count(*ic)   *
+          calibBAO_On_Cycle_Ch0.count(*ic)  *
+          calibBAO_On_Cycle_Ch1.count(*ic)  *
+          calibBAO_Off_Cycle_Ch0.count(*ic) *
+          calibBAO_Off_Cycle_Ch1.count(*ic) ) != 0 ? true : false;
+
+      if(para.debuglev_>9) {
+        cout << "Calibration coefficients for cycle "<<*ic << endl; 
+        if (isCycleCalibrated) {
+          cout << "Cycle calibrated " << endl;
+          cout << "Off Run Ch0 " << calibBAO_Off_Run_Ch0[*ic] << " "
+               << "Ch1 " << calibBAO_Off_Run_Ch1[*ic] << "\n"
+               << "On Run Ch0 " << calibBAO_On_Run_Ch0[*ic] << " "
+               << "Ch1 " << calibBAO_On_Run_Ch1[*ic] << "\n"
+               << "Off Cycle Ch0 " << calibBAO_Off_Cycle_Ch0[*ic] << " "
+               << "Ch1 " << calibBAO_Off_Cycle_Ch1[*ic] << "\n"
+               << "On Cycle Ch0 " << calibBAO_On_Cycle_Ch0[*ic] << " "
+               << "Ch1 " << calibBAO_On_Cycle_Ch1[*ic] << endl;
+        } else {
+          cout << "Cycle NOT calibrated " << endl;
+        }
+      }//debug
+
+      if ( ! isCycleCalibrated ) continue;
+
+      totalNumberCycles++;
+      //Fill NTuple
+      xnt[0] = rdateOfRun;
+      xnt[1] = totalNumberRuns;
+      xnt[2] = totalNumberCycles;
+      xnt[3] = calibBAO_Off_Run_Ch0[*ic];
+      xnt[4] = calibBAO_Off_Run_Ch1[*ic];
+      xnt[5] = calibBAO_On_Run_Ch0[*ic];
+      xnt[6] = calibBAO_On_Run_Ch1[*ic];
+      xnt[7] = calibBAO_Off_Cycle_Ch0[*ic];
+      xnt[8] = calibBAO_Off_Cycle_Ch1[*ic];
+      xnt[9] = calibBAO_On_Cycle_Ch0[*ic];
+      xnt[10] = calibBAO_On_Cycle_Ch1[*ic];
+      xnt[11] = invcalibBAO_Off_Run_Ch0[*ic];
+      xnt[12] = invcalibBAO_Off_Run_Ch1[*ic];
+      xnt[13] = invcalibBAO_On_Run_Ch0[*ic];
+      xnt[14] = invcalibBAO_On_Run_Ch1[*ic];
+      xnt[15] = invcalibBAO_Off_Cycle_Ch0[*ic];
+      xnt[16] = invcalibBAO_Off_Cycle_Ch1[*ic];
+      xnt[17] = invcalibBAO_On_Cycle_Ch0[*ic];
+      xnt[18] = invcalibBAO_On_Cycle_Ch1[*ic];
+
+      calibcoeffs.Fill(xnt);
+
+      //Quit if enough
+      if (totalNumberCycles >= para.maxNumberCycles_) break;    
+
+    }//eo loop on cycles
+    if (totalNumberCycles >= para.maxNumberCycles_) break;          
+    totalNumberRuns++;
+  }//eo loop on spectra in a file
+  cout << "End of jobs: we have treated " << totalNumberCycles << " cycles" << endl;
+
+  {//save the results
+    stringstream tmp;
+    tmp << totalNumberCycles;
+    string fileName = para.outPath_+"/calibcoeffTuple_"+StringToLower(para.sourceName_)+"-"+tmp.str()+"Cycles.ppf";
+
+    POutPersist fos(fileName);
+    cout << "Save output in " << fileName << endl;
+
+    string tag = "calib";
+    fos << PPFNameTag(tag) << calibcoeffs;
+
+  }//end of save
+
+}
 //-------------------------------------------------------
 //Compute the mean of Diff ON-OFF BAO-calibrated spectra and also the mean/sigma of rebinned spectra 
@@ -356 +812 @@
   sa_size_t chCalibHigh = freqToChan(para.rcalibFreq_ + (para.rcalibBandFreq_*0.5));
   //Lower and Higher freq.  just arround 1420.4MHz Freq. bin to perform mean follow up
-  sa_size_t ch1420Low  = freqToChan(1420.4-0.2);
-  sa_size_t ch1420High = freqToChan(1420.4+0.2);
-
+  sa_size_t ch1420Low;
+  sa_size_t ch1420High;
+  if (para.sourceName_ == "Abell85") {
+    ch1420Low  = freqToChan(1420.4-0.2);  // Abell85
+    ch1420High = freqToChan(1420.4+0.2);
+  } else if (para.sourceName_ == "Abell1205") {
+    ch1420Low  = freqToChan(1420.4-0.3);  // Abell1205
+    ch1420High = freqToChan(1420.4+0.2);
+  } else if (para.sourceName_ == "Abell2440") {
+    ch1420Low  = freqToChan(1420.4);      // Abell2440
+    ch1420High = freqToChan(1420.4+0.3);
+  } else {
+    ch1420Low  = freqToChan(1420.4-0.2);  // Abell85
+    ch1420High = freqToChan(1420.4+0.2);
+  }
   //Lower and Higher freq. on the sides of 1420.4Mhz Freq. bin to perform mean follow up
   sa_size_t ch1420aLow  = freqToChan(1418);
@@ -397 +865 @@
     string srcLower = tokens2[0];
 
+
+    
     PInPersist fin(*iFile);
     vector<string> vec = fin.GetNameTags();
@@ -612 +1082 @@
     nr = calibBAOfactors_Off_Run_Ch0.NRows();
     for (sa_size_t ir=0; ir<nr; ++ir){
-      cout << "Calib. Off Run Ch0 cycle ["<< calibBAOfactors_Off_Run_Ch0(ir,0)<<"], val "
-           << calibBAOfactors_Off_Run_Ch0(ir,1) << endl;
+//       cout << "Calib. Off Run Ch0 cycle ["<< calibBAOfactors_Off_Run_Ch0(ir,0)<<"], val "
+//         << calibBAOfactors_Off_Run_Ch0(ir,1) << endl;
 
       calibBAO_Off_Run_Ch0[(int)calibBAOfactors_Off_Run_Ch0(ir,0)]
@@ -664 +1134 @@
                << "Ch1 " << calibBAO_On_Cycle_Ch1[*ic] << endl;
         } else {
-          cout << "Cycle << " << *ic <<" NOT calibrated for file " << *iFile << endl;
+          cout << "Cycle NOT calibrated " << endl;
         }
       }//debug
@@ -708 +1178 @@
       //JEC 29/10/11 add ON-OFF/OFF
       TMatrix<r_4> diffOnOffOvOff_noCalib(diffOnOff_noCalib,false); //do not share data
-      TMatrix<r_4> aSpecOffFiltered(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);
+      TMatrix<r_4> aSpecOffFitltered(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);
+      medianFiltering(aSpecOff,halfWidth,aSpecOffFitltered);
       
-      diffOnOffOvOff_noCalib.Div(aSpecOffFiltered); //division in place
+      diffOnOffOvOff_noCalib.Div(aSpecOffFitltered); //division in place
       meanDiffONOFFovOFF_noCalib += diffOnOffOvOff_noCalib;
 
@@ -721 +1191 @@
       meanDiffONOFF_perCycleCalib += diffOnOff_perCycleCalib;
 
+      //
       totalNumberCycles++;
-
       //Fill NTuple
       xnt[0] = totalNumberCycles;
@@ -796 +1266 @@
       //JEC 18/11/11 follow up the 1400-1420MHz OFF only
       TMatrix<r_4> aSpecOffovOff(aSpecOff,false);
-      aSpecOffovOff.Div(aSpecOffFiltered);
+      aSpecOffovOff.Div(aSpecOffFitltered);
 
       TVector<r_4> meanInRange_1410a1415Freq_OFF_noCalib(NUMBER_OF_CHANNELS);
-      //      meanInRange(aSpecOff,ch1410,ch1415,meanInRange_1410a1415Freq_OFF_noCalib);
       meanInRange(aSpecOffovOff,ch1410,ch1415,meanInRange_1410a1415Freq_OFF_noCalib);
 
@@ -806 +1275 @@
 
       TMatrix<r_4> aSpecOnovOff(aSpecOn,false);
-      aSpecOnovOff.Div(aSpecOffFiltered);
+      aSpecOnovOff.Div(aSpecOffFitltered);
 
       TVector<r_4> meanInRange_1410a1415Freq_ON_noCalib(NUMBER_OF_CHANNELS);
-      //meanInRange(aSpecOn,ch1410,ch1415,meanInRange_1410a1415Freq_ON_noCalib);
       meanInRange(aSpecOnovOff,ch1410,ch1415,meanInRange_1410a1415Freq_ON_noCalib);
 
@@ -884 +1352 @@
     
     tag = "onoffevol";
-    fos << PPFNameTag(tag) << onoffevolution;
- 
+    fos << PPFNameTag(tag) << onoffevolution;    
   }//end of save
-
-
 }
 //JEC 14/11/11 New meanRawDiffOnOffCycles START
@@ -911 +1376 @@
   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;
@@ -951 +1417 @@
   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];
-    if (para.debuglev_>90){
-      cout << "date <" << dateOfRun << ">" << endl;
-    }
-    vector<string> tokens2;
-    split(tokens[2],".",tokens2);
-    string srcLower = tokens2[0];
-
+    }
     PInPersist fin(*iFile);
     vector<string> vec = fin.GetNameTags();
@@ -993 +1447 @@
         regexp(iSpec->c_str(),matchstr.c_str(),&b,&e);
         if (para.debuglev_>90){
-          cout << " spectra <" << *iSpec << ">" << endl;
+          cout << " spactra <" << *iSpec << ">" << endl;
           cout << " cycle " << iSpec->substr(b) << endl;
         }
@@ -1021 +1475 @@
     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
@@ -1054 +1490 @@
       //Perform the difference ON-OFF
       TMatrix<r_4> diffOnOff_noCalib = aSpecOn - aSpecOff;
-
       meanDiffONOFF_noCalib += diffOnOff_noCalib;
       
       //JEC 29/10/11 add ON-OFF/OFF
       TMatrix<r_4> diffOnOffOvOff_noCalib(diffOnOff_noCalib,false); //do not share data
-      TMatrix<r_4> aSpecOffFiltered(NUMBER_OF_CHANNELS,NUMBER_OF_FREQ);
+      TMatrix<r_4> aSpecOffFitltered(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);
+      medianFiltering(aSpecOff,halfWidth,aSpecOffFitltered);
       
-      diffOnOffOvOff_noCalib.Div(aSpecOffFiltered); //division in place
+      diffOnOffOvOff_noCalib.Div(aSpecOffFitltered); //division in place
       meanDiffONOFFovOFF_noCalib += diffOnOffOvOff_noCalib;
+
 
       //JEC 15/11/11 add ON/OFF and OFF/OFF
       TMatrix<r_4> onOvOff(aSpecOn,false);
-      onOvOff.Div(aSpecOffFiltered);
+      onOvOff.Div(aSpecOffFitltered);
       meanONovOFF_noCalib += onOvOff;
       
       TMatrix<r_4> offOvOff(aSpecOff,false);
-      offOvOff.Div(aSpecOffFiltered);
+      offOvOff.Div(aSpecOffFitltered);
       meanOFFovOFF_noCalib += offOvOff;
 
       totalNumberCycles++;
+
 
       //Fill NTuple
       xnt[0] = totalNumberCycles;
  
+
       //Follow up arround the 1420.4MHz Freq.
       TVector<r_4> meanInRange_1420Freq_noCalib(NUMBER_OF_CHANNELS);
@@ -1097 +1535 @@
 
 
-      //JEC 25/10/11 follow 1410-1415 MHz bande protege et n'inclue pas le 1420.4Mhz de la Galaxie
-      TVector<r_4> meanInRange_1410a1415Freq_noCalib(NUMBER_OF_CHANNELS);
-      meanInRange(diffOnOffOvOff_noCalib,ch1410,ch1415,meanInRange_1410a1415Freq_noCalib);
-      xnt[5] = meanInRange_1410a1415Freq_noCalib(0);
-      xnt[6] = meanInRange_1410a1415Freq_noCalib(1);
-
-      //JEC 18/11/11 follow up the 1410-1415MHz OFF only
+      //JEC 25/10/11 follow 1400-1420 MHz bande protege et n'inclue pas le 1420.4Mhz de la Galaxie
+      TVector<r_4> meanInRange_1400a1420Freq_noCalib(NUMBER_OF_CHANNELS);
+      meanInRange(diffOnOffOvOff_noCalib,ch1400,ch1420,meanInRange_1400a1420Freq_noCalib);
+      xnt[5] = meanInRange_1400a1420Freq_noCalib(0);
+      xnt[6] = meanInRange_1400a1420Freq_noCalib(1);
+
+      //JEC 18/11/11 follow up the 1400-1420MHz OFF only
       TMatrix<r_4> aSpecOffovOff(aSpecOff,false);
-      aSpecOffovOff.Div(aSpecOffFiltered);
+      aSpecOffovOff.Div(aSpecOffFitltered);
 
       TVector<r_4> meanInRange_1410a1415Freq_OFF_noCalib(NUMBER_OF_CHANNELS);
@@ -1114 +1552 @@
 
       TMatrix<r_4> aSpecOnovOff(aSpecOn,false);
-      aSpecOnovOff.Div(aSpecOffFiltered);
+      aSpecOnovOff.Div(aSpecOffFitltered);
 
       TVector<r_4> meanInRange_1410a1415Freq_ON_noCalib(NUMBER_OF_CHANNELS);
@@ -1162 +1600 @@
 
     tag = "onoffevol";
-    fos << PPFNameTag(tag) << onoffevolution;
+    fos << PPFNameTag(tag) << onoffevolution;    
 
   }//end save
-
-
 }
 //JEC 14/11/11 New meanRawDiffOnOffCycles END
@@ -1526 +1962 @@
     } else if (action == "meanCalibBAODiffOnOff") {
       meanCalibBAODiffOnOffCycles();
+    } else if (action == "calibCoeffNtp") {
+      calibCoeffNtp();
     } else {
       msg = "Unknown action " + action;