void AA::selectBDsMu(AA::Selection s){
  bool ok = false;

  vector<PType> types;
  types.push_back(AA::K_PLUS);
  types.push_back(AA::K_MINUS);
  types.push_back(AA::PI_PLUS);

  //select muons
  PtlLst muons;
  PtlLst* plstMuo = 0;
  if(AA::indMuo < 0) {
    AA::selectMuo(muons);
    plstMuo = &muons;
  } else {
    plstMuo = &AA::lstCommon[AA::indMuo];
  }

  //Loop over muons
  for(PtlLstCIt p = plstMuo->begin(); p != plstMuo->end(); ++p) {
    Ptl* ptl = *p;
    const Vrt* pvp = ptl->primaryVertex();
    if(pvp == 0) return;

    //Search for Ds
    for(PtlLstCIt p3 = ptlNoV0Lst.begin(); p3 != ptlNoV0Lst.end(); ++p3){
      Ptl* ptl3 = *p3;
      if(ptl3 == ptl) continue;
      if(ptl3->pt() < 0.7) continue;
      if(ptl3->primaryVertex() != pvp) continue;
      if(ptl3->nSMT() < 2) continue;
      if(ptl3->jet() != ptl->jet()) continue;

      for(PtlLstCIt p4 = p3; p4 != ptlNoV0Lst.end(); ++p4){
	Ptl* ptl4 = *p4;
	if(ptl4 == ptl) continue;
	if(ptl4 == ptl3) continue;
	if(ptl4->pt() < 0.7) continue;
	if(ptl4->primaryVertex() != pvp) continue;
	if(ptl4->nSMT() < 2) continue;
	if(ptl4->jet() != ptl->jet()) continue;

	if(ptl3->q()*ptl4->q() > 0) continue;

	for(PtlLstCIt p5 = ptlNoV0Lst.begin(); 
	    p5 != ptlNoV0Lst.end(); ++p5){
	  Ptl* ptl5 = *p5;
	  if(ptl5 == ptl) continue;
	  if(ptl5 == ptl3) continue;
	  if(ptl5 == ptl4) continue;
	  if(ptl5->pt() < 0.7) continue;
	  if(ptl5->primaryVertex() != pvp) continue;
	  if(ptl5->nSMT() < 2) continue;
	  if(ptl5->jet() != ptl->jet()) continue;

	  //Make vertex of these 3 particles
	  TrkLst lstT;
	  lstT.push_back(ptl3);
	  lstT.push_back(ptl4);
	  lstT.push_back(ptl5);
	  Vrt vrt;
	  if(!vrt.fill(pvp->x(),&lstT)) continue;
	  if(!vrt.filter()) continue;
	  if(vrt.size() != 3) continue;

	  //Find primary vertex with minimal distance to V0 vertex
	  double lv,vlv;
	  vrt.distance(pvp,lv,vlv);
	  double lxy, vlxy;
	  vrt.distanceXY(pvp,lxy,vlxy);
	  
	  list<Chain*> chb;
	  vrt.chains(&chb);
	  list<Chain*>::iterator pc = chb.begin();
	  Chain* pc3 = *pc;
	  Chain* pc4 = *(++pc);
	  Chain* pc5 = *(++pc);
	
	  double mkk = AA::xMass(&vrt,pc3,pc4,AA::K_PLUS,AA::K_MINUS);
	  if(mkk > 1.07) continue;

	  int q = ptl3->q() + ptl4->q() + ptl5->q();
	  Ptl ds;
	  if(!ds.combine(&vrt,q)) continue;
	  if(!ds.associateToVrt(&AA::vrtPBox)) continue;
	  if(ds.primaryVertex() != pvp) continue;

	  double cxyds = ds.cosXY();
	  
	  //Mass of ds
	  double mds = ds.mass(types);
	  
	  //Combine Ds and lepton
	  lstT.clear();
	  lstT.push_back(&ds);
	  lstT.push_back(ptl);
	  Vrt vrtB;
	  if(!vrtB.fill(pvp->x(),&lstT)) continue;
	  if(!vrtB.filter()) continue;
	  if(vrtB.size() != 2) continue;
	  if(vrtB.chi2() > 16.) continue;

	  //Find primary vertex with minimal distance to V0 vertex
	  double lxyB, vlxyB;
	  vrtB.distanceXY(pvp,lxyB,vlxyB);
	  double lxyBD, vlxyBD;
	  vrt.distanceXY(&vrtB,lxyBD,vlxyBD);

	  //compute mass of l+ and Ds
	  vrtB.chains(&chb);
	  pc = chb.begin();
	  Chain* pcD = *pc;
	  Chain* pcl = *(++pc);
	  double xmb = AA::xMass(&vrtB,pcl,pcD,AA::mass[MU_PLUS],mds);

	  Ptl b;
	  q = ptl3->q() + ptl4->q() + ptl5->q() + ptl->q();
	  if(!b.combine(&vrtB,q)) continue;
	  if(!b.associateToVrt(&AA::vrtPBox)) continue;
	  if(b.primaryVertex() != pvp) continue;
	  double cxyB = b.cosXY();
	    
	  //Trigger selections
	  bool okds = true;
	  if(mds > 2.3) okds = false;

	  if(cxyds < 0.85) okds = false;
	  if(lxy*lxy < 16.*vlxy) okds = false;
	  if(lxyB*lxyB > vlxyB*16. && cxyB < 0.95) okds = false;
	  if(lxyB > lxy && lxyBD*lxyBD > vlxyBD*9.) okds = false;

	  if(xmb > 5.5) okds = false;
	  if(xmb < 1.5) okds = false;
	  //if(abs(ptl3->imp(1))>0.25 || abs(ptl3->imp(2))>0.4) okds = false;
	  //if(abs(ptl4->imp(1))>0.25 || abs(ptl4->imp(2))>0.4) okds = false;
	  //if(abs(ptl5->imp(1))>0.25 || abs(ptl5->imp(2))>0.4) okds = false;
	  //if(abs(ptl->imp(1))>0.25  || abs(ptl->imp(2))>0.4)  okds = false;

	  if(ptl->nSMT()  < 2) okds = false;
	  if(ptl->nCFT()  < 2) okds = false;
	  if(ptl3->nCFT() < 2) okds = false;
	  if(ptl4->nCFT() < 2) okds = false;
	  if(ptl5->nCFT() < 2) okds = false;

	  if(ptl->ptot() < 3.) okds = false;
	  if(ptl->muon()->chisqloc() < -0.5) okds = false;

	  if(okds) ok = true;

#ifdef AA_FILL_HIST
	  //optional filling of control histograms
	  if(okds) {
	    if(mkk < 1.006 || mkk > 1.032) continue;
	    if(ptl3->pt() < 0.7 || ptl4->pt() < 0.7) continue;
	    if(ptl5->pt() < 0.7) continue;

	    //Pt of muon with respect to D
	    double cmud = (ptl->mom(1)*ds.mom(1) + ptl->mom(2)*ds.mom(2) +
			   ptl->mom(3)*ds.mom(3))/ptl->ptot()/ds.ptot();
	    double ptmu = sqrt(abs(1.-cmud*cmud))*ptl->ptot();
	    
	    if(ptl->pt()   < 1.5) continue;
	    if(ptl->ptot() < 4.) continue;
	    if(ptl->nSMT() < 2) continue;
	    if(ptl->nCFT() < 2) continue;
	    if(ptl->muon()->chisqloc() < -0.5) continue;
	    if(ptl->muon()->chisqloc() > 25.)  continue;
	    if(abs(ptl->muon()->sctimeA()) > 7.) continue;
	    
	    if(abs(ptl3->imp(1))>0.25 || abs(ptl3->imp(2))>0.4) okds = false;
	    if(abs(ptl4->imp(1))>0.25 || abs(ptl4->imp(2))>0.4) okds = false;
	    if(abs(ptl5->imp(1))>0.25 || abs(ptl5->imp(2))>0.4) okds = false;
	    if(abs(ptl->imp(1))>0.25  || abs(ptl->imp(2))>0.4)  okds = false;
	    if(ptl3->nCFT() < 2) continue;
	    if(ptl4->nCFT() < 2) continue;
	    if(ptl5->nCFT() < 2) continue;

	    if(lxy*lxy < 25.*vlxy) continue;
	    if(cxyds < 0.95) continue;
	    if(vrtB.chi2() > 9.) continue;
	    if(ptl3->chi2Vrt() < 9. && ptl4->chi2Vrt()<9. && 
	       ptl5->chi2Vrt() < 9.) continue;
	    if(ptl5->chi2Vrt() < 2.) continue;

	    if(lxyB > lxy && lxyBD*lxyBD > vlxyBD*4.) continue;
	    if(lxyB*lxyB > vlxyB*16. && cxyB < 0.95) continue;

	    if(xmb < 2.8) continue;
	    if(xmb > 5.2) continue;
	    if(ptmu < 1.0) continue;
	    if(b.ptot() < 8.) continue;

	    if(ptl->q()*ptl5->q() < 0) {
	      HFILL(4021,mds,0.,1.);
	    } else {
	      HFILL(4022,mds,0.,1.);
	    }
	    //	    if(mds>1.6) {
	      //cout << "Selected: " << AA::runNumber << " " << AA::evtNumber 
	      //	   << " " << mds << " " << mkk << " " 
	      //<< ptl->q()*ptl5->q() << endl;
	    //}
	  }
#endif
	}
      }
    }
  }
  if(ok) AA::selectionMask |= AA::maskBDsMu;
}