Changeset 2014 in Sophya for trunk


Ignore:
Timestamp:
May 28, 2002, 7:30:45 PM (23 years ago)
Author:
ansari
Message:

Amelioration de calcul d'offset (possibilite de coupure sur coordonnees
galactique)
Nouveau processeur nettoyeur (SimpleCleaner)
Ajout programme de traitement aksj02.cc , donnees Archeops KS1/KS3

Reza 28/5/2002

Location:
trunk/ArchTOIPipe
Files:
3 added
3 edited

Legend:

Unmodified
Added
Removed

  • trunk/ArchTOIPipe/ProcWSophya/simoffset.cc

    r2008 r2014  
    1818  : poly((degpol > 1)?degpol:1)
    1919{
     20  SetParams(mwsz, nptfit, degpol);
     21  totnscount = 0;
     22  totnbblock = 0;
     23  bmincut = bmaxcut = -9999.;
     24  bgalcut = false;
     25  ns_flgcut = ns_bgalcut = 0;
     26  SavePolyNTuple();
     27}
     28
     29SimpleOffsetEstimator::~SimpleOffsetEstimator()
     30{
     31}
     32
     33void SimpleOffsetEstimator::SetParams(int mwsz, int nptfit, int degpol)
     34{
    2035  mWSz = (mwsz > 8) ? mwsz : 8;
    2136  nPtFit = (nptfit > degpol+2) ? nptfit : degpol+2;
    2237  if (nPtFit%2 == 0) nPtFit++;
    2338  degPol = (degpol > 1)?degpol:1;
    24   totnscount = 0;
    25   totnbblock = 0;
    26   SavePolyNTuple();
    27 }
    28 
    29 SimpleOffsetEstimator::~SimpleOffsetEstimator()
    30 {
     39 
     40}
     41
     42void SimpleOffsetEstimator::SetBGalCut(double bmin, double bmax)
     43{
     44  bgalcut = true;
     45  bmincut = bmin;
     46  bmaxcut = bmax;
    3147}
    3248
    … …  
    3652     << " SimpleOffsetEstimator::PrintStatus() - MeanWSize= " << mWSz << " NPtFit="
    3753     << nPtFit << " DegPoly=" << degPol << " poly.Degre()=" << poly.Degre() << endl;
     54  if (bgalcut)
     55    os << " bGalCut = Yes , bGalMin = " << bmincut << " bGalMax= " << bmaxcut << endl;
     56  else os << " bGalCut = No " << endl;
    3857  TOIProcessor::PrintStatus(os);
    39   os << " ProcessedSampleCount=" << ProcessedSampleCount() << endl;
     58  os << " ProcessedSampleCount=" << ProcessedSampleCount()
     59     << " NS_FlgCut= " << ns_flgcut << " NS_BGalCut= " << ns_bgalcut << endl;
    4060  os << " ------------------------------------------------------ " << endl; 
    4161}
    … …  
    4565  cout << "SimpleOffsetEstimator::init" << endl;
    4666  declareInput("in");
     67  declareInput("bgal");
    4768  declareOutput("offset");
    4869  declareOutput("out");
    4970  declareOutput("incopie");
    50   declareOutput("poly_a0");
    51   declareOutput("poly_a1");
    52   declareOutput("poly_a2");
    53   declareOutput("poly_sn0");
    54   declareOutput("mean_y");
    55   declareOutput("sig_y");
    56   declareOutput("mean_x");
     71  declareOutput("bgalcopie");
     72  //  declareOutput("mean_y");
     73  //  declareOutput("sig_y");
     74  //  declareOutput("mean_x");
    5775  name = "SimpleOffsetEstimator";
    5876}
    … …  
    6684  bool fgout = checkOutputTOIIndex(1);
    6785  bool fgincopie = checkOutputTOIIndex(2);
    68   bool fga0 = checkOutputTOIIndex(3);
    69   bool fga1 = checkOutputTOIIndex(4);
    70   bool fga2 = checkOutputTOIIndex(5);
    71   bool fgsn0 = checkOutputTOIIndex(6);
    72   bool fgmeany = checkOutputTOIIndex(7);
    73   bool fgsigy = checkOutputTOIIndex(8);
    74   bool fgmeanx = checkOutputTOIIndex(9);
    75  
     86  bool fgbgalcopie = checkOutputTOIIndex(3);
     87  //  bool fgmeany = checkOutputTOIIndex(4);
     88  //  bool fgsigy = checkOutputTOIIndex(5);
     89  //  bool fgmeanx = checkOutputTOIIndex(6);
     90
    7691  if (!checkInputTOIIndex(0)) {
    7792    cerr << " SimpleOffsetEstimator::run() - Input TOI (in) not connected! "
    7893         << endl;
    7994    throw ParmError("SimpleOffsetEstimator::run() Input TOI (in) not connected!");
     95  }
     96  if (bgalcut && !checkInputTOIIndex(1)) {
     97    cerr << " SimpleOffsetEstimator::run() - Input TOI bgal not connected! "
     98         << endl;
     99    throw ParmError("SimpleOffsetEstimator::run() Input TOI bgal not connected!");
    80100  }
    81101  if (!fgoffset && !fgout) {
    … …  
    88108
    89109  // NTuple pour sauvegarde des coeff de poly 
    90   char * nomsnt[] = {"sncur", "sn0", "meanx", "meany", "sigy", "a0", "a1", "a2", "ycur"};
    91   XNTuple xntp(0, 9, 0, 0, nomsnt);
     110  char * nomsnt[] = {"sncur", "sn0", "meanx", "meany", "sigy",
     111                       "a0", "a1", "a2", "ycur", "nok", "nbblkok"};
     112  XNTuple xntp(0, 11, 0, 0, nomsnt);
    92113
    93114  try {
    … …  
    103124    Vector vinc(wsize);
    104125    TVector<uint_8> vfgc(wsize);
     126    Vector bgal(wsize);
    105127
    106128   
    … …  
    131153    double y_last = 0.;
    132154
     155    double last_meanok = 0.;
     156    double last_sigok = 1.;
     157    double meanx_forlast = 0.;
     158    bool fg_last_meansig = false;
     159 
    133160    // Boucle sur les sampleNum
    134161    // 1ere partie, on traite par paquets de wsize
    … …  
    142169
    143170      getData(0, k, wsize, vin.Data(), vfg.Data());
     171      if (bgalcut) {
     172        getData(1, k, wsize, bgal.Data());
     173        if (fgbgalcopie) putData(3, k, wsize, bgal.Data());
     174      }
     175
     176      fg_last_meansig = false;
     177
     178      if (kb == 0) {
     179        last_meanok = vin.Sum()/wsize;
     180        last_sigok = vin.SumX2()/wsize - last_meanok*last_meanok;
     181      }
    144182
    145183      // Calcul moyenne et sigma du bloc
    146184      nok = 0.;  meanx = 0.;
    147185      mean = 0.;  sig = 0.;
    148      
     186      meanx_forlast = 0.;
    149187      for(j=0; j<wsize; j++) {
    150         if ( vfg(j) ) continue;
     188        meanx_forlast += k+j;
     189        if ( vfg(j) ) { ns_flgcut++;  continue; }
     190        if (bgalcut && (bgal(j) > bmincut) && (bgal(j) < bmaxcut) )
     191          { ns_bgalcut++;  continue; }
    151192        mean += vin(j);
    152193        sig += vin(j)*vin(j);
    … …  
    172213        meanx /= nok;
    173214        sig = sig/nok-mean*mean;
    174         if (sig < 1.e-6*mean) sig = 1.e-6*mean;
     215        if (sig < 1.e-10*mean) sig = 1.e-10*mean;
     216        if (sig < 1.e-39) sig = 1.e-39;
    175217        int kk = nbblkok%nPtFit;
    176218        nbblkok++;
    … …  
    178220        YErr(kk) = sig;
    179221        X(kk) = meanx;
    180       }
    181      
    182       if ((nok>3.5) && (nbblkok >= nPtFit) ) {
     222        last_meanok = mean;
     223        last_sigok = sig;
     224      }
     225      else {
     226        int kk = nbblkok%nPtFit;
     227        Y(kk) = last_meanok;
     228        YErr(kk) = last_sigok*10.;
     229        X(kk) = meanx_forlast/wsize;   
     230        fg_last_meansig = true;
     231      }
     232     
     233      if (((nok>3.5) || fg_last_meansig) && (nbblkok >= nPtFit) ) {
    183234        // On force le fit a garder une certaine continuite
    184235        Y(nPtFit) = y_last;
    185236        double smin, smax;
    186237        YErr(Range(0,0,nPtFit)).MinMax(smin, smax);
     238        if (smax < 1.e-39)   smax = 1.e-39;
    187239        YErr(nPtFit) = smax/10.;
    188240        X(nPtFit) = sn_last;
    … …  
    190242        X0 = X;
    191243        X0 -= sn0;
    192         poly.Fit(X0,Y,YErr,degPol,errCoef);
     244        Vector polsave(degPol);
     245        for(int jj=0; jj<=poly.Degre(); jj++) polsave(jj) = poly[jj];
     246        try {
     247          poly.Fit(X0,Y,YErr,degPol,errCoef);
     248        }
     249        catch (CaughtSignalExc& excsig) {
     250          cout << " -- simoffset.cc/ catched CaughtSignalExc - Msg= "
     251               << excsig.Msg() << " kb=" << kb << " k=" << k << endl;
     252          cout << " X0= " << X0 << " Y=" << Y << " YErr=" << YErr << endl;
     253          for(int jj=0; jj<=poly.Degre(); jj++) poly[jj] = polsave(jj);
     254        }
    193255      }
    194256      else {
    … …  
    208270      }
    209271      if (ntpoly) {    // Remplissage du XNTuple de controle
    210         float xnt[10];
     272        float xnt[12];
    211273        xnt[0] = k;
    212274        xnt[1] = sn0;
    … …  
    218280        xnt[7] = poly[2];
    219281        xnt[8] = poly(k-sn0);
     282        xnt[9] = nok;
     283        xnt[10] = nbblkok;
    220284        xntp.Fill(NULL, xnt, NULL, NULL);
    221285      }
    … …  
    244308        vinc.ReSize(wszt);
    245309        vfgc.ReSize(wszt);
     310        bgal.ReSize(wszt);
    246311        int jb = 0;
    247312        for(int kbsc=kbs; kbsc<nblk; kbsc++) {
    … …  
    270335        putData(1, k, wsize, vinc.Data(), vfgc.Data());
    271336      }
    272      
    273       if (fga0) {
    274         vout = poly[0];
    275         putData(3, k, wsize, vout.Data());
    276       }
    277       if (fga1) {
    278         vout = poly[1];
    279         putData(4, k, wsize, vout.Data());
    280       }
    281       if (fga2) {
    282         vout = poly[2];
    283         putData(5, k, wsize, vout.Data());
    284       }
    285       if (fgsn0) {
    286         vout = sn0;
    287         putData(6, k, wsize, vout.Data());
    288       }
     337      /*
    289338      if (fgmeany) {
    290339        vout = mean;
    291         putData(7, k, wsize, vout.Data());
     340        putData(4, k, wsize, vout.Data());
    292341      }
    293342      if (fgsigy) {
    294343        vout = sig;
    295         putData(8, k, wsize, vout.Data());
     344        putData(5, k, wsize, vout.Data());
    296345      }
    297346
    298347      if (fgmeanx) {
    299348        vout = meanx;
    300         putData(9, k, wsize, vout.Data());
    301       }
     349        putData(6, k, wsize, vout.Data());
     350      }
     351      */
    302352
    303353      klast+=wsize;
    … …  
    316366      k = klast+1;
    317367      getData(0, k, wsize, vin.Data(), vfg.Data());
     368
     369      if (bgalcut) {
     370        getData(1, k, wsize, bgal.Data());
     371        if (fgbgalcopie) putData(3, k, wsize, bgal.Data());
     372      }
     373
    318374      for(j=0; j<wsize; j++)
    319375        voff(j) = poly(k+j-sn0);
    … …  
    324380        putData(1, k, wsize, vin.Data(), vfg.Data());
    325381      }
    326      
    327       if (fga0) {
    328         vout = poly[0];
    329         putData(3, k, wsize, vout.Data());
    330       }
    331       if (fga1) {
    332         vout = poly[1];
    333         putData(4, k, wsize, vout.Data());
    334       }
    335       if (fga2) {
    336         vout = poly[2];
    337         putData(5, k, wsize, vout.Data());
    338       }
    339       if (fgsn0) {
    340         vout = sn0;
    341         putData(6, k, wsize, vout.Data());
    342       }
     382
     383      /*     
    343384      if (fgmeany) {
    344385        vout = mean;
    345         putData(7, k, wsize, vout.Data());
     386        putData(4, k, wsize, vout.Data());
    346387      }
    347388      if (fgsigy) {
    348389        vout = sig;
    349         putData(8, k, wsize, vout.Data());
     390        putData(5, k, wsize, vout.Data());
    350391      }
    351392     
    352393      if (fgmeanx) {
    353394        vout = meanx;
    354         putData(9, k, wsize, vout.Data());
    355       }
     395        putData(6, k, wsize, vout.Data());
     396      }
     397      */
    356398     
    357399      klast+=wsize;

  • trunk/ArchTOIPipe/ProcWSophya/simoffset.h

    r2004 r2014  
    1818// Structure generale :
    1919//
    20 //             -------------------------
    21 //  toi in --> |                       | ---> out  (toi = in_offset)
    22 //             | SimpleOffsetEstimator | ---> offset (toi)
    23 //            |                       | ---> Autres toi optionnel
    24 //             |                       |
    25 //             -------------------------
     20//               -------------------------
     21//    toi in --> |                       | ---> out  (toi = in_offset)
     22//               | SimpleOffsetEstimator | ---> offset (toi)
     23//  bgal(opt)--> |                       | ---> Autres toi optionnel
     24//               |                       |
     25//               -------------------------
    2626
    2727class SimpleOffsetEstimator : public TOIProcessor {
    … …  
    3030  virtual  ~SimpleOffsetEstimator();
    3131
     32  void     SetParams(int mwsz=256, int nptfit=5, int degpol=2);
    3233  virtual void  init(); 
    3334  virtual void  run();
     35
     36  virtual void  SetBGalCut(double bmin, double bmax);
    3437
    3538  inline int_8  ProcessedSampleCount() const { return totnscount; }
    … …  
    4750  bool ntpoly;
    4851  string ntpolyname;
     52  double bmincut, bmaxcut;
     53  bool bgalcut;
     54  int_4 ns_flgcut;
     55  int_4 ns_bgalcut;
    4956};
    5057

  • trunk/ArchTOIPipe/TestPipes/simofftst.cc

    r2008 r2014  
    33//                               Christophe Magneville
    44//                               Reza Ansari
    5 // $Id: simofftst.cc,v 1.2 2002-05-16 13:13:00 ansari Exp $
     5// $Id: simofftst.cc,v 1.3 2002-05-28 17:30:45 ansari Exp $
    66
    77/*   Test de processeurs ds simtoipr.cc   - Reza Avril 2001
    … …  
    1111            -intoi boloMuV_26 -wtoi 8192 -wdegli 1024 \
    1212            inputbolo.fits filt.fits xx.ppf
    13 # Avec Filtre de Fourier
     13# Avec Filtre
    1414csh> simofftst -start 104385384 -end 104399964 -range -500,500 \
    1515            -intoi boloMuV_26 -wtoi 8192 -wdegli 1024 \
    16             -wfft 4096 -keepfft 5 -killfreq 15,4,2.5 \
    1716            inputbolo.fits filtfft.fits spectre.ppf
    18 # Autre exemple b545k02
    19 cool> ./simofftst -start 104389122 -end 104649122 -range -1000.,1000. -intoi boloMu
    20 V_23 -wtoi 8192 -wdegli 1024 -wfft 4096 -keepfft 5 -killfreq 15,4,2.5 Cmv/b545k02
    21 2_16.00_16.49.fits b545.fits b545.ppf
     17csh> $exedir/simofftst -start $sn0 -end $sn1  -intoi $cleanname \
     18            -wtoi 8192 -wdegli 1024 -degli 3.,2.,3,2,0 \
     19            -soff 128,9,3 -gfilt 16,2.5 -prstat \
     20            -soffnt $outppf -bgalcut -20.,20. -bgal $glat \
     21            -bgalfile $point $cleanfile $outfile
    2222
    2323*/
    … …  
    3838#include "timing.h"
    3939#include "histinit.h"
     40#include "psighand.h"
    4041#include <stdexcept>
    4142
    … …  
    5253         << "         [-degli ns,ns2,mxpt,mnpt,wrec] [-soff wsz,nptfit,degpol] \n"
    5354         << "         [-soffnt PPFName] [-gfilt wsz,sigma] \n"
     55         << "         [-bgalcut bmin,bmax] [-bgal toiname] [-bgalfile pointFitsName] \n"
    5456         << "         [-nooutflg] [-nowrtms] [-nowrticd] [-prstat] [-useseqbuff] \n"
    5557         << "         inFitsName outFitsName \n"
    … …  
    7779         << endl;
    7880  }
     81  if (fgerr) exit(1);
    7982}
    8083
    … …  
    124127  double gfilt_sigma = 2.;
    125128
     129  // Coordinate file for galactic cut
     130  string pointfile = "pointgal.fits";
     131  double bmin = -20;
     132  double bmax = 20.;
     133  string bgaltoi = "bgal";
     134  bool bgalcut = false;
     135
    126136  // File names
    127137  string infile;
    … …  
    183193      ia++;
    184194    }   
     195    else if (strcmp(arg[ia],"-bgalcut") == 0) {
     196      if (ia == narg-1) Usage(true);
     197      sscanf(arg[ia+1],"%lf,%lf", &bmin, &bmax);
     198      bgalcut = true;
     199      ia++;
     200    }
     201    else if (strcmp(arg[ia],"-bgalfile") == 0) {
     202      if (ia == narg-1) Usage(true);
     203      pointfile = arg[ia+1];
     204      ia++;
     205    }
     206    else if (strcmp(arg[ia],"-bgal") == 0) {
     207      if (ia == narg-1) Usage(true);
     208      bgaltoi = arg[ia+1];
     209      ia++;
     210    }
    185211    else if (strcmp(arg[ia],"-wfft") == 0) {
    186212      if (ia == narg-1) Usage(true); 
    … …  
    226252  cout << " Initializing SOPHYA ... " << endl;
    227253  SophyaInit();
     254  SophyaConfigureSignalhandling(true);
     255
    228256  InitTim();
    229257
    … …  
    263291    double G_a = 1./(G_sigma*sqrt(M_PI*2.));
    264292    SimpleFilter filt(gfilt_wsz, SimpleFilter::GaussFilter, G_a, G_sigma);
    265    
     293
     294         
     295    FITSTOIReader* rbgal = NULL;
     296    if (bgalcut) {  // if Galactic cut
     297      cout << "> Creating bgal FITSTOIReader object - InFile=" << pointfile
     298           << " bgaltoiname= " << bgaltoi << endl;
     299      cout << "  offes.SetBGalCut( " << bmin << "," << bmax << ")" << endl;
     300      rbgal = new FITSTOIReader(pointfile);
     301      offes.SetBGalCut(bmin, bmax);     
     302    }
     303
    266304    cout << "> Creating FITSTOIWriter OutFitsName= " << outfile << endl;
    267305    FITSTOIWriter w(outfile);
    … …  
    281319      plombier.Connect(offes, "out", w, "degoff", "", 0, fg_wrtflag);
    282320     
     321    if (bgalcut) {
     322      inname = "bgal";
     323      plombier.Connect(*rbgal, bgaltoi, offes, inname);     
     324    }
    283325
    284326    if (fg_wrticd) {
    … …  
    315357    mgr->joinAll();
    316358    PrtTim("End threads");
     359    if (bgalcut)  delete rbgal;
    317360
    318361  }
Note: See TracChangeset for help on using the changeset viewer.