Changeset 2014 in Sophya for trunk/ArchTOIPipe/ProcWSophya


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/ProcWSophya
Files:
2 added
2 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
Note: See TracChangeset for help on using the changeset viewer.