/home/bes3soft/bes3soft/Boss/7.0.2/dist/7.0.2/Reconstruction/MdcPatRec/MdcxReco/MdcxReco-00-01-59/src/MdcxFindSegs.cxx

Go to the documentation of this file.

#include "MdcxReco/MdcxFindSegs.h"
#include "MdcxReco/MdcxHit.h"
#include "MdcGeom/MdcDetector.h"
#include "MdcTrkRecon/mdcWrapWire.h"
#include "MdcGeom/BesAngle.h"
#include "MdcxReco/MdcxParameters.h"

using std::cout;
using std::endl;
using std::ostream;
#include "AIDA/IHistogram1D.h"
#include "AIDA/IHistogram2D.h"
extern AIDA::IHistogram1D*  g_csmax4;
extern AIDA::IHistogram1D*  g_csmax3;

int MdcxFindSegs::wireGroups[11][288][16] = { { {0} } };

MdcxFindSegs::MdcxFindSegs(MdcDetector* g, const HepAList<MdcxHit> &hitlist, int debug){
  m_debug = debug;
  gm = g;
  dchitlist = hitlist;
  initWireGroups();
  process();
}

MdcxFindSegs::MdcxFindSegs(const HepAList<MdcxHit> &hitlist, int debug){
  m_debug = debug;
  gm = MdcDetector::instance();
  dchitlist = hitlist;
  initWireGroups();
  process();
}

MdcxFindSegs::~MdcxFindSegs() { KillList(); }

void MdcxFindSegs::initWireGroups(void) {
  //Assure initWireGroups run only once
  static bool alreadyInit = false;
  if (alreadyInit) return;

  //Make wire groups
  int lastsl= gm->nSuper();
  for (int sl = 0; sl < lastsl; sl++) { // loop over superLayers
    const MdcLayer *layArray[4];
    const MdcSuperLayer* slayer = gm->SuperLayer(sl);
    int nslay = slayer->nLayers();
    for (int i = 0; i < nslay; i++) layArray[i] = slayer->layer(i);
    int cellMax = layArray[1]->nWires();
    // loop over cells
    for (int cellIndex = 0; cellIndex < cellMax ; cellIndex++) {
      int wireIndex[4]; 
      if (slayer->slayNum() > 4) {
        wireIndex[0]= cellIndex;
        wireIndex[1]= cellIndex+1;
        wireIndex[2]= cellIndex;
        wireIndex[3]= cellIndex+1;

        wireIndex[1] = mdcWrapWire(wireIndex[1], layArray[1]->nWires());
        wireIndex[3] = mdcWrapWire(wireIndex[3], layArray[3]->nWires());
      } else {
        wireIndex[1]= cellIndex+1; 
        //take No.2 wire as referenc wire
        double phi = layArray[1]->getWire(wireIndex[1])->phiE();
        for(int i = 0; i < 4; i++ ) {
          if ( i == 1 ) continue;
          //change reference wire
          if ( i == 3 ) phi = layArray[2]->getWire(wireIndex[2])->phiE();
          BesAngle tmp(phi - layArray[i]->phiEPOffset());//yzhang temp
          if ( i==3 ) {
            wireIndex[i]=(tmp==0)?1:(int)ceil(layArray[i]->nWires()*tmp.rad()/twopi);
          }else{
            wireIndex[i]=(tmp==0)?-1:(int)floor(layArray[i]->nWires()*tmp.rad()/twopi);
          }
          wireIndex[i] = mdcWrapWire(wireIndex[i], layArray[i]->nWires());
        }
      }
      // set wireGroups
      int* w = wireGroups[sl][cellIndex];
      for (int i = 0; i < 4; i++) {
        for (int j = 0; j < 4; j++) {
          //wireGroups[sl][cellIndex][4*i+j] = wireIndex[i] + (j - 1);
          w[4*i+j] = mdcWrapWire((wireIndex[i]+j-1), layArray[i]->nWires());
        }
      }
    } //end loop cells
  } //end loop superLayers

  //Set alreadyInit flag
  alreadyInit = true;
  return;
}

void MdcxFindSegs::process() {
  //if(3 == m_debug) cout <<" Event contains "<<dchitlist.length()<<" MdcxHit "<<endl;
  // for last super layer use 44
  int lflag[44][288] = { {0} };
  int pflag[44][288] = { {0} };
  int ix,iy,iyp,iyn,cellMax;

  //yzhang 2009-10-20 do not use poison
  int ipoison=0;
  //number of last layer
  int lastlay= gm->nLayer();

  // tag hit's lfalg
  int k=0; 
  while(dchitlist[k]){
    MdcxHit* temp = dchitlist[k]; 
    k++; 
    lflag[temp->Layer()][temp->WireNo()]=k;
  }

  // try poisoning
  if (ipoison){
    for (ix=0; ix<lastlay; ix++){
      cellMax = gm->Layer(ix)->nWires();
      // if have hits in left and right , poison this cell
      for (iy=0; iy<cellMax; iy++){
        iyp=iy+1; if (iyp == cellMax)iyp=0; 
        iyn=iy-1; if (iyn == -1)iyn=cellMax-1;
        if ((lflag[ix][iyp] != 0) && (lflag[ix][iyn] != 0)) {
          pflag[ix][iy] = 1;
          lflag[ix][iy] = 0;
          //g_poison->fill(ix,iy);
        }
      }
      //for (iy=0; iy<cellMax; iy++) {if (pflag[ix][iy] != 0) lflag[ix][iy] = 0;}
    }
  }  //FIXME (optimize ???)
  
  float csmax_4 = MdcxParameters::csmax4;
  float csmax_3 = MdcxParameters::csmax3;
  int lastsl= gm->nSuper();//lastSLayNum;

  // loop superlayer
  for (sl = 0; sl < lastsl; sl++) {
    int fsl=4*sl;
    int l0=fsl, l1=fsl+1, l2=fsl+2, l3=fsl+3; 
    int iprt = 0; 
    //initialize superlayer pointer in layArray[]
    const MdcLayer *layArray[4];
    const MdcSuperLayer* slayer = gm->SuperLayer(sl);
    int nslay = slayer->nLayers();
    for (int i = 0; i < nslay; i++) layArray[i] = slayer->layer(i);
    if(3 == m_debug){
      cout<<"slayer No.  "<<slayer->index()<<endl;//yzhang debug
    }
    //reference point (2,1)
    cellMax = layArray[1]->nWires();

    // loop over cells
    for (int cellIndex=0; cellIndex<cellMax ; cellIndex++) {
      int* w = wireGroups[sl][cellIndex];
      unsigned int sig_mark = 0;
      for (int ilayer = l0; ilayer <= l3; ilayer++) {
        for (int iwire = 3; iwire >= 0; iwire--) {
          if (lflag[ilayer][w[4*(ilayer-l0)+iwire]]) {
            sig_mark |= 0x1;
          }
          sig_mark <<= 1;
        }
      }
      sig_mark >>= 1;

      int goodSegNo = 0;
      int nPat = m_segPat.patt4_size;
      int iPat = (sig_mark & 0x0200) ? 0 : 11;
      for ( ; iPat < nPat; iPat++) {
        int pat = m_segPat.patt4[iPat];
        if ((pat & sig_mark) == pat) {
          if(3 == m_debug) {
            cout<< "pat " << std::hex << pat << std::dec << "  with wire";
            for (int tmpi = 0; tmpi < 4; tmpi++) {
              cout<<"("<<l0+tmpi<<","<<w[4*tmpi+m_segPat.wirePat4[iPat][tmpi]-1]<<")";
            }
            cout << endl;
          }
          int w0 = lflag[l0][w[0 +m_segPat.wirePat4[iPat][0]-1]] - 1;
          int w1 = lflag[l1][w[4 +m_segPat.wirePat4[iPat][1]-1]] - 1;
          int w2 = lflag[l2][w[8 +m_segPat.wirePat4[iPat][2]-1]] - 1;
          int w3 = lflag[l3][w[12+m_segPat.wirePat4[iPat][3]-1]] - 1;
          int tw[4] = {w0, w1, w2, w3};

          int namb = m_segPat.ambPat4_size[iPat];
          for (int iamb = 0; iamb < namb; iamb++) {
            int amb = m_segPat.ambigPatt4[iPat][iamb];
            fithel = trial(w0, w1, w2, w3, amb);
            if(3 == m_debug){
              cout<<"chisq "<<fithel.Chisq()
                <<" <? csmax4 "<<csmax_4<<  endl;
              if (fithel.Chisq() < csmax_4) {
                cout<<"Accept this seg"<<endl;
              }else{
                cout<<"DROP this seg"<<endl;
              }
            }
            if(g_csmax4) g_csmax4->fill(fithel.Chisq());
            if (fithel.Chisq() < csmax_4) {
              if (iprt) printseg(fithel, pat,amb);
              appendseg(fithel, pat, amb);
              goodSegNo++;
            }
          }
        }
      }
      if (goodSegNo != 0) continue;
      nPat = m_segPat.patt3_size;
      iPat = (sig_mark & 0x0200) ? 0 : 14;
      for ( ; iPat < nPat; iPat++) {
        if (iPat > nPat-3) {
          if ((iPat == nPat-2) && (sig_mark&0x2121 == 0x2121)) continue;
          if ((iPat == nPat-1) && (sig_mark&0x2122 == 0x2122)) continue;
        }
        int pat = m_segPat.patt3[iPat];
        if ((pat & sig_mark) == pat) {
          if(3 == m_debug) {
            cout<<"MdcxFindSegs: in pat "<<std::hex<<pat<<std::dec<<" with wire";
            for (int tmpi = 0; tmpi < 4; tmpi++) {
              if (m_segPat.wirePat3[iPat][tmpi] == 0) continue;
              cout<< " (" << l0+tmpi << "," << w[4*tmpi+m_segPat.wirePat3[iPat][tmpi]-1] << ")";
            }
            cout<< endl;
          }
          int wn[3];
          for (int iw = 0, iwp = 0; iwp < 4; iwp++) {
            int wireNo = m_segPat.wirePat3[iPat][iwp];
            if( wireNo == 0 ) continue;
            wn[iw++] = lflag[l0+iwp][w[4*iwp+wireNo-1]] - 1;
          }

          int namb = m_segPat.ambPat3_size[iPat];
          for (int iamb = 0; iamb < namb; iamb++) {
            int amb = m_segPat.ambigPatt3[iPat][iamb];
            fithel = trial(wn[0], wn[1], wn[2], amb);//FIXME
            if(3 == m_debug){
              cout<<"chisq "<<fithel.Chisq()
                <<" <? csmax3 "<<csmax_3<<  endl;
              if (fithel.Chisq() < csmax_3) {
                cout<<"Accept this seg"<<endl;
              }else{
                cout<<"DROP this seg"<<endl;
              }
            }
            if(g_csmax3) g_csmax3->fill(fithel.Chisq());
            if (fithel.Chisq() < csmax_3) {
              if (iprt) printseg(fithel, pat, amb);
              appendseg(fithel, pat ,amb);
            }
          }
        }
      }//end of nPat3
    }
  }

  if(3 == m_debug){
    cout << "MdcxFindSegs found " << MdcxSeglist.length() << " segs" << endl;
  }
  return;
}

void MdcxFindSegs::print(ostream &o, int pmax)const {
  int mcheck=pmax; if (MdcxSeglist.length() < pmax)mcheck=MdcxSeglist.length();
  o << " First " << mcheck << " Drift Chamber Segs:" << endl;
  for(int i=0; i<mcheck; i++){
    o << " Superlayer # " << MdcxSeglist[i]->SuperLayer();
    o << " # of hits " << MdcxSeglist[i]->Nhits() << endl;
  } 
} // end of print

MdcxFittedHel MdcxFindSegs::trial(int i1, int i2, int i3, int i4, int amb) {
  HepAList<MdcxHit> seg;
  MdcxHit* t1 = dchitlist[i1]; MdcxHit* t2 = dchitlist[i2];  
  MdcxHit* t3 = dchitlist[i3]; MdcxHit* t4 = dchitlist[i4];  
  seg.append(t1); seg.append(t2); seg.append(t3); seg.append(t4); 
  double dx, dy, d0; 
  if (amb==0){dx=t4->xneg()-t1->xneg(); dy=t4->yneg()-t1->yneg(); d0=-t1->d();}
  if (amb==1){dx=t4->xneg()-t1->xpos(); dy=t4->yneg()-t1->ypos(); d0= t1->d();}
  if (amb==2){dx=t4->xpos()-t1->xneg(); dy=t4->ypos()-t1->yneg(); d0=-t1->d();}
  if (amb==3){dx=t4->xpos()-t1->xpos(); dy=t4->ypos()-t1->ypos(); d0= t1->d();}
  double phi0 = atan2(dy,dx); dx=t1->x(); dy=t1->y(); 
  MdcxHel hel(d0,phi0,0.0,0.0,0.0,0.0,111,0,dx,dy);
  hel.SetTurnFlag(1);
  MdcxFittedHel temp(seg, hel, 1.0);
  if (3 == m_debug) {
    cout<<"trial4 amb "<<amb
      << ": phi0 " << phi0 << " d0 " << d0 
      << " dx " << dx << " dy " << dy << endl;
    //temp.FitPrint();
  }
  return temp;
}

MdcxFittedHel MdcxFindSegs::trial(int i1, int i2, int i3, int amb) {
  HepAList<MdcxHit> seg;
  MdcxHit* t1 = dchitlist[i1];
  MdcxHit* t2 = dchitlist[i2];
  MdcxHit* t3 = dchitlist[i3];
  seg.append(t1); seg.append(t2); seg.append(t3);
  double dx,dy,d0; 
  if (amb==0){dx=t3->xneg()-t1->xneg(); dy=t3->yneg()-t1->yneg(); d0=-t1->d();}
  if (amb==1){dx=t3->xneg()-t1->xpos(); dy=t3->yneg()-t1->ypos(); d0= t1->d();}
  if (amb==2){dx=t3->xpos()-t1->xneg(); dy=t3->ypos()-t1->yneg(); d0=-t1->d();}
  if (amb==3){dx=t3->xpos()-t1->xpos(); dy=t3->ypos()-t1->ypos(); d0= t1->d();}

  double phi0 = atan2(dy,dx); dx=t1->x(); dy=t1->y(); 
  MdcxHel hel(d0,phi0,0.0,0.0,0.0,0.0,111,0,dx,dy);
  hel.SetTurnFlag(1);
  MdcxFittedHel temp(seg,hel,1.0);

  if (3 == m_debug) {
    cout<<"trial3 amb "<<amb
      << ": phi0 " << phi0 << " d0 " << d0 
      << " dx " << dx << " dy " << dy << endl;
  }
  return temp;
}

void MdcxFindSegs::appendseg(MdcxFittedHel& fithel,int pat,int amb){
  MdcxSeg *tempseg = new MdcxSeg(fithel, pat, amb);
  MdcxSeglist.append(tempseg);
}

void MdcxFindSegs::printseg(MdcxFittedHel& fithel, int pat, int amb, int subtry){
  cout << "Seg: pat " << pat << " amb " << amb << " subtry " << subtry 
    << " sl " << sl  << " fail "
    << fithel.Fail() << " chi2 " << fithel.Chisq() << endl;
}

---