void UASM(){ //AnalyticSphereMethod *ASM){
  const Int_t MAXSOLN=2048*8;
  AnalyticSphereMethod *ASM=new AnalyticSphereMethod();
  Double_t aveX = 0.0;                                                                
  Double_t aveY = 0.0;                                                                
  Double_t aveZ = 0.0;   
  Float_t phisolns[MAXSOLN]={0};
  Float_t thesolns[MAXSOLN]={0};
  Float_t rMag;
  Int_t solcnt = 0;                                                                   
  Int_t nentry=0;
  TH1F *h1phi=new TH1F("h1phi","h1phi",360,-180,180);
  TH1F *h1the=new TH1F("h1the","h1the",360,-180,180);
  ASMvtx.r=-NNN; ASMvtx.phi=-NNN; ASMvtx.the=-NNN; ASMvtx.phiTH1=-NNN; ASMvtx.theTH1=-NNN;
   for(int irx1=0; irx1<MCH16; irx1++) {         
    if(!Rx[irx1].LVtx) continue;
    ASM->SetAnt(0,Double_t(Rx[irx1].x),Double_t(Rx[irx1].y),Double_t(Rx[irx1].z),Double_t(Rx[irx1].HitTime));        
    for(int irx2=irx1+1; irx2<MCH16; irx2++) {                                        
      if(!Rx[irx2].LVtx) continue;                                                    
      ASM->SetAnt(1,Double_t(Rx[irx2].x),Double_t(Rx[irx2].y),Double_t(Rx[irx2].z),Double_t(Rx[irx2].HitTime));  
      for(int irx3=irx2+1; irx3<MCH16; irx3++) { 
	if(!Rx[irx3].LVtx) continue; 
	ASM->SetAnt(2,Double_t(Rx[irx3].x),Double_t(Rx[irx3].y),Double_t(Rx[irx3].z),Double_t(Rx[irx3].HitTime));
	for(int irx4=irx3+1; irx4<MCH16; irx4++) {                                    
	  if(!Rx[irx4].LVtx) continue;       
	  if(solcnt-1>=MAXSOLN){
	    cout<<"ASM MAXSOLN!!!!"<<endl;
	    continue;
	  }
	  nentry++;
	  ASM->SetAnt(3,Double_t(Rx[irx4].x),Double_t(Rx[irx4].y),Double_t(Rx[irx4].z),Double_t(Rx[irx4].HitTime));        
  	  if(ASM->Solve()) {                                                          
	    aveX += ASM->GetX();                                                     
	    aveY += ASM->GetY();
	    aveZ += ASM->GetZ();        
	    rMag=sqrt(ASM->GetX()*ASM->GetX() + ASM->GetY()*ASM->GetY() + ASM->GetZ()*ASM->GetZ());
	    phisolns[solcnt]=atan2(ASM->GetY(),ASM->GetX())*TMath::RadToDeg();
	    thesolns[solcnt]=acos(ASM->GetZ()/rMag)*TMath::RadToDeg();
	    h1phi->Fill(phisolns[solcnt]);
	    h1the->Fill(thesolns[solcnt]);
	    solcnt++;                                                                 
	    ASMsolnCh[irx1]++; ASMsolnCh[irx2]++; ASMsolnCh[irx3]++; ASMsolnCh[irx4]++;
	  }    
	}
      }
    }
   }
   if(solcnt>0) {                                                              
     aveX /= solcnt;                                                          
     aveY /= solcnt;                                                         
     aveZ /= solcnt;
     ASMvtx.r=sqrt(aveX*aveX + aveY*aveY + aveZ*aveZ);
     ASMvtx.phi=TMath::RadToDeg()*atan2(aveY,aveX);
     ASMvtx.the=TMath::RadToDeg()*acos(aveZ/ASMvtx.r);
     ASMvtx.dphi=TMath::RMS(solcnt,phisolns);
     ASMvtx.dthe=TMath::RMS(solcnt,thesolns);
     ASMvtx.ncomb=solcnt;
     //     h1phi->Fit("gaus");     h1the->Fit("gaus");     ASMvtx.phiTH1=h1phi->GetFunction("gaus")->GetParameter(1);     ASMvtx.theTH1=h1the->GetFunction("gaus")->GetParameter(1);
   }
   //   printf("ASM.h nfitvtx=%d nentry=%d nsoln=%d xyz <ASMphi>=%g/<ASMthe>=%g <ASMphisolns>=%g / <ASMthesolns>=%g ASM.phiTH1=%g ASM.theTH1=%g\n",nFitVtx,nentry,solcnt,ASMvtx.phi,ASMvtx.the,TMath::Mean(solcnt,phisolns),TMath::RMS(solcnt,thesolns),ASMvtx.phiTH1,ASMvtx.theTH1);
   printf("ASM.h nfitvtx=%d nentry=%d nsoln=%d xyz <ASMphi>=%g/<ASMthe>=%g <ASMphisolns>=%g / <ASMthesolns>=%g\n",nFitVtx,nentry,solcnt,ASMvtx.phi,ASMvtx.the,TMath::Mean(solcnt,phisolns),TMath::RMS(solcnt,thesolns));
   delete ASM;
   delete h1phi;
   delete h1the;
}