{
//////////////////////////////////////////////////////////
//   This file has been automatically generated 
//     (Wed Jun  2 18:07:43 2004 by ROOT version3.05/07)
//   from TTree AncSvx/SVX Ancestors
//   found on file: ancsvx_027.root
//////////////////////////////////////////////////////////
// derived from Pat's Charm2.C Jun 04 HvH
//////////////////////////////////////////////////////////
//Reset ROOT and connect tree file
   gROOT->Reset();
   TFile *f = new TFile("ancsvx_027.root");

//Declaration of leaves types
   Float_t         TRACK;
   Float_t         NFILE;
   Float_t         PTOT;
   Float_t         PTHETA;
   Float_t         PPHI;
   Float_t         R_VERTEX;
   Float_t         Z_VERTEX;
   Float_t         THET_VER;
   Float_t         PHI_VER;
   Float_t         ITPARENT;
   Float_t         IDPARENT;
   Float_t         IDPART;
   Float_t         ITORIGIN;
   Float_t         IDORIGIN;
   Float_t         IHIT;
   Float_t         LAYER;
   Float_t         THETA;
   Float_t         PHI;
   Float_t         XGLOBAL;
   Float_t         YGLOBAL;
   Float_t         ZGLOBAL;
   Float_t         PMOMX;
   Float_t         PMOMY;
   Float_t         PMOMZ;
   Float_t         IPC123;
   Float_t         PC1THETA;
   Float_t         PC1PHI;
   Float_t         PC1RDIST;
   Float_t         PC1ZDIST;
   Float_t         PC2THETA;
   Float_t         PC2PHI;
   Float_t         PC2RDIST;
   Float_t         PC2ZDIST;
   Float_t         PC3THETA;
   Float_t         PC3PHI;
   Float_t         PC3RDIST;
   Float_t         PC3ZDIST;
   Float_t         DELE;
   Float_t         XLOCIN;
   Float_t         YLOCIN;
   Float_t         ZLOCIN;
   Float_t         XLOCOUT;
   Float_t         YLOCOUT;
   Float_t         ZLOCOUT;
   Float_t         HITVOL0;
   Float_t         HITVOL1;
   Float_t         HITVOL2;
   Float_t         HITVOL3;
   Float_t         HITVOL4;
   Float_t         HITVOL5;
   Float_t         NHIT;
   Float_t         Z0_EVENT;
   Float_t         B_IMPACT;
   Float_t         EVENT;

//Get a pointer to individual leaves and set their addresses
   TObjArray *leaves = AncSvx->GetListOfLeaves();
   TLeaf *l_TRACK           = (TLeaf*)leaves->At(0);   l_TRACK->SetAddress(&TRACK);
   TLeaf *l_NFILE           = (TLeaf*)leaves->At(1);   l_NFILE->SetAddress(&NFILE);
   TLeaf *l_PTOT            = (TLeaf*)leaves->At(2);   l_PTOT->SetAddress(&PTOT);
   TLeaf *l_PTHETA          = (TLeaf*)leaves->At(3);   l_PTHETA->SetAddress(&PTHETA);
   TLeaf *l_PPHI            = (TLeaf*)leaves->At(4);   l_PPHI->SetAddress(&PPHI);
   TLeaf *l_R_VERTEX        = (TLeaf*)leaves->At(5);   l_R_VERTEX->SetAddress(&R_VERTEX);
   TLeaf *l_Z_VERTEX        = (TLeaf*)leaves->At(6);   l_Z_VERTEX->SetAddress(&Z_VERTEX);
   TLeaf *l_THET_VER        = (TLeaf*)leaves->At(7);   l_THET_VER->SetAddress(&THET_VER);
   TLeaf *l_PHI_VER         = (TLeaf*)leaves->At(8);   l_PHI_VER->SetAddress(&PHI_VER);
   TLeaf *l_ITPARENT        = (TLeaf*)leaves->At(9);   l_ITPARENT->SetAddress(&ITPARENT);
   TLeaf *l_IDPARENT        = (TLeaf*)leaves->At(10);   l_IDPARENT->SetAddress(&IDPARENT);
   TLeaf *l_IDPART          = (TLeaf*)leaves->At(11);   l_IDPART->SetAddress(&IDPART);
   TLeaf *l_ITORIGIN        = (TLeaf*)leaves->At(12);   l_ITORIGIN->SetAddress(&ITORIGIN);
   TLeaf *l_IDORIGIN        = (TLeaf*)leaves->At(13);   l_IDORIGIN->SetAddress(&IDORIGIN);
   TLeaf *l_IHIT            = (TLeaf*)leaves->At(14);   l_IHIT->SetAddress(&IHIT);
   TLeaf *l_LAYER           = (TLeaf*)leaves->At(15);   l_LAYER->SetAddress(&LAYER);
   TLeaf *l_THETA           = (TLeaf*)leaves->At(16);   l_THETA->SetAddress(&THETA);
   TLeaf *l_PHI             = (TLeaf*)leaves->At(17);   l_PHI->SetAddress(&PHI);
   TLeaf *l_XGLOBAL         = (TLeaf*)leaves->At(18);   l_XGLOBAL->SetAddress(&XGLOBAL);
   TLeaf *l_YGLOBAL         = (TLeaf*)leaves->At(19);   l_YGLOBAL->SetAddress(&YGLOBAL);
   TLeaf *l_ZGLOBAL         = (TLeaf*)leaves->At(20);   l_ZGLOBAL->SetAddress(&ZGLOBAL);
   TLeaf *l_PMOMX           = (TLeaf*)leaves->At(21);   l_PMOMX->SetAddress(&PMOMX);
   TLeaf *l_PMOMY           = (TLeaf*)leaves->At(22);   l_PMOMY->SetAddress(&PMOMY);
   TLeaf *l_PMOMZ           = (TLeaf*)leaves->At(23);   l_PMOMZ->SetAddress(&PMOMZ);
   TLeaf *l_IPC123          = (TLeaf*)leaves->At(24);   l_IPC123->SetAddress(&IPC123);
   TLeaf *l_PC1THETA        = (TLeaf*)leaves->At(25);   l_PC1THETA->SetAddress(&PC1THETA);
   TLeaf *l_PC1PHI          = (TLeaf*)leaves->At(26);   l_PC1PHI->SetAddress(&PC1PHI);
   TLeaf *l_PC1RDIST        = (TLeaf*)leaves->At(27);   l_PC1RDIST->SetAddress(&PC1RDIST);
   TLeaf *l_PC1ZDIST        = (TLeaf*)leaves->At(28);   l_PC1ZDIST->SetAddress(&PC1ZDIST);
   TLeaf *l_PC2THETA        = (TLeaf*)leaves->At(29);   l_PC2THETA->SetAddress(&PC2THETA);
   TLeaf *l_PC2PHI          = (TLeaf*)leaves->At(30);   l_PC2PHI->SetAddress(&PC2PHI);
   TLeaf *l_PC2RDIST        = (TLeaf*)leaves->At(31);   l_PC2RDIST->SetAddress(&PC2RDIST);
   TLeaf *l_PC2ZDIST        = (TLeaf*)leaves->At(32);   l_PC2ZDIST->SetAddress(&PC2ZDIST);
   TLeaf *l_PC3THETA        = (TLeaf*)leaves->At(33);   l_PC3THETA->SetAddress(&PC3THETA);
   TLeaf *l_PC3PHI          = (TLeaf*)leaves->At(34);   l_PC3PHI->SetAddress(&PC3PHI);
   TLeaf *l_PC3RDIST        = (TLeaf*)leaves->At(35);   l_PC3RDIST->SetAddress(&PC3RDIST);
   TLeaf *l_PC3ZDIST        = (TLeaf*)leaves->At(36);   l_PC3ZDIST->SetAddress(&PC3ZDIST);
   TLeaf *l_DELE            = (TLeaf*)leaves->At(37);   l_DELE->SetAddress(&DELE);
   TLeaf *l_XLOCIN          = (TLeaf*)leaves->At(38);   l_XLOCIN->SetAddress(&XLOCIN);
   TLeaf *l_YLOCIN          = (TLeaf*)leaves->At(39);   l_YLOCIN->SetAddress(&YLOCIN);
   TLeaf *l_ZLOCIN          = (TLeaf*)leaves->At(40);   l_ZLOCIN->SetAddress(&ZLOCIN);
   TLeaf *l_XLOCOUT         = (TLeaf*)leaves->At(41);   l_XLOCOUT->SetAddress(&XLOCOUT);
   TLeaf *l_YLOCOUT         = (TLeaf*)leaves->At(42);   l_YLOCOUT->SetAddress(&YLOCOUT);
   TLeaf *l_ZLOCOUT         = (TLeaf*)leaves->At(43);   l_ZLOCOUT->SetAddress(&ZLOCOUT);
   TLeaf *l_HITVOL0         = (TLeaf*)leaves->At(44);   l_HITVOL0->SetAddress(&HITVOL0);
   TLeaf *l_HITVOL1         = (TLeaf*)leaves->At(45);   l_HITVOL1->SetAddress(&HITVOL1);
   TLeaf *l_HITVOL2         = (TLeaf*)leaves->At(46);   l_HITVOL2->SetAddress(&HITVOL2);
   TLeaf *l_HITVOL3         = (TLeaf*)leaves->At(47);   l_HITVOL3->SetAddress(&HITVOL3);
   TLeaf *l_HITVOL4         = (TLeaf*)leaves->At(48);   l_HITVOL4->SetAddress(&HITVOL4);
   TLeaf *l_HITVOL5         = (TLeaf*)leaves->At(49);   l_HITVOL5->SetAddress(&HITVOL5);
   TLeaf *l_NHIT            = (TLeaf*)leaves->At(50);   l_NHIT->SetAddress(&NHIT);
   TLeaf *l_Z0_EVENT        = (TLeaf*)leaves->At(51);   l_Z0_EVENT->SetAddress(&Z0_EVENT);
   TLeaf *l_B_IMPACT        = (TLeaf*)leaves->At(52);   l_B_IMPACT->SetAddress(&B_IMPACT);
   TLeaf *l_EVENT           = (TLeaf*)leaves->At(53);   l_EVENT->SetAddress(&EVENT);

//This is the loop skeleton
//To read only selected branches, Insert a statement like:
//          leafname->GetBranch()->GetEvent(i); instead of AncSvx->GetEvent(i);
  Int_t nentries = AncSvx->GetEntries();
  cout <<" AncSvx nentries = "<<nentries<<endl;
// arrays of floats
Float_t event[nentries];
Float_t track[nentries];
Float_t ptot[nentries];
Float_t ptheta[nentries];
Float_t pphi[nentries];
Float_t theta[nentries];
Float_t phi[nentries];
Float_t x[nentries];
Float_t y[nentries];
Float_t z[nentries];
Float_t px[nentries];
Float_t py[nentries];
Float_t pz[nentries];
Float_t layer[nentries];
Float_t id[nentries];
Int_t ievent[10001];


  Int_t nbytes = 0;
  Float_t oldevent = 1.;
Int_t iev = 0;
  ievent[1]=0;
                               // loop over all records and copy info to arrays
  for (Int_t i=0; i<nentries; i++) {
    AncSvx->GetEvent(i);
    event[i]=EVENT;           // watch for new event number:
    if (EVENT>oldevent) {
      oldevent=EVENT;
      iev=EVENT;
      ievent[iev]=i; 
      if (iev%100==0) {cout <<"record "<<i<<", event "<<iev<<" "<<endl;}
    }
    track[i]=TRACK;           // fill hit arrays
    ptot[i]=PTOT;
    ptheta[i]=PTHETA;
    pphi[i]=PPHI;
    theta[i]=THETA;
    phi[i]=PHI;
    x[i]=XGLOBAL;
    y[i]=YGLOBAL;
    z[i]=ZGLOBAL;
    px[i]=PMOMX;
    py[i]=PMOMY;
    pz[i]=PMOMZ;
    layer[i]=LAYER;
    id[i]=IDPART;  
  }                           // end loop to stuff arrays
  cout <<"Total records "<<i<<", events "<<iev<<" "<<endl;

// tracks entries below are track#(1-2) layer(1-12) (x,y,z,px,py,pz)
Float_t tracks[3][13][7];
// vertex entries below are track#(1-2) (xvertex,yvertex,zvertex,ptot,ptheta,pphi)
Float_t vertex[3][7];
Int_t idp = 0;
Int_t itrack = 0;
Int_t ilayer = 0;
Float_t pzvertex = 0.;
Float_t ptvertex = 0.;
Float_t pxvertex = 0.;
Float_t pyvertex = 0.;
Float_t degtorad = 3.14159/180.;
Int_t idpart[3];

//define histograms
TH1F *g1 = new TH1F("g1","xintercept id=5",50,-0.05,0.05);
TH1F *g2 = new TH1F("g2","yintercept id=5",50,-0.05,0.05);
TH1F *g3 = new TH1F("g3","xintercept id=6",50,-0.05,0.05);
TH1F *g4 = new TH1F("g4","yintercept id=6",50,-0.05,0.05);
TH1F *g5 = new TH1F("g5","zxvertex id=5",50,-0.2,0.2);
TH1F *g6 = new TH1F("g6","zyvertex id=5",50,-0.2,0.2);
TH1F *g7 = new TH1F("g7","zxvertex id=6",50,-0.2,0.2);
TH1F *g8 = new TH1F("g8","zyvertex id=6",50,-0.2,0.2);
TH1F *g9 = new TH1F("g9","xpair vertex",50,-0.2,0.2);
TH1F *g10 = new TH1F("g10","ypair vertex",50,-0.2,0.2);
TH1F *g11 = new TH1F("g11","pair vertex",50,-0.2,0.2);
TH1F *g12 = new TH1F("g12","pair mass",50,2.0,4.0);
TH1F *g13 = new TH1F("g13","tracker xmatch",50,-10.,10.);
TH1F *g14 = new TH1F("g14","tracker ymatch",50,-10.,10.);
TH1F *g15 = new TH1F("g15","tracker pzmatch",50,0.,2.5);
TH1F *g16 = new TH1F("g16","si dphi deg id=5",50,-1.,1.);
TH1F *g17 = new TH1F("g17","si dphi deg id=6",50,-1.,1.);
TH1F *g18 = new TH1F("g18","si dphi*pz degGeV",50,-5.,5.);
TH1F *g19 = new TH1F("g19","rintercept",50,-0.05,0.2);
TH1F *g20 = new TH1F("g20","zrvertex",50,0.2,1.);
// TH1F *g20 = new TH1F("g20","zrvertex",50,-0.1,1.);
TH1F *g21 = new TH1F("g21","Si rmatch",50,-0.02,0.02);
TH1F *g22 = new TH1F("g22","rvtx thrown",50,-0.05,0.2);
TH1F *g23 = new TH1F("g23","zvtx thrown",50,0.2,1.);
// TH1F *g23 = new TH1F("g23","zvtx thrown",50,-0.1,1.);
TH2F *g24 = new TH2F("g24","xintercept 1 vs 2",50,-0.05,0.05,50,-0.05,0.05);
TH2F *g25 = new TH2F("g25","yintercept 1 vs 2",50,-0.05,0.05,50,-0.05,0.05);
TH2F *g26 = new TH2F("g26","xint1 vs yint1",50,-0.05,0.05,50,-0.05,0.05);
TH2F *g27 = new TH2F("g27","xint2 vs yint2",50,-0.05,0.05,50,-0.05,0.05);
TH1F *g28 = new TH1F("g28","rint fit",50,-0.05,0.05);
TH1F *g29 = new TH1F("g29","rpair vertex",50,-0.1,0.1);
TH1F *g30 = new TH1F("g30","rpair intercept",50,-0.05,0.05);
TH2F *g31 = new TH2F("g31","layer vs x",13,0.0,13.0, 40,-20.0,20.0);
TH2F *g32 = new TH2F("g32","layer vs y",13,0.0,13.0, 40,-20.0,20.0);
TH2F *g33 = new TH2F("g33","layer vs z",13,0.0,13.0, 40,-40.0,40.0);
TH2F *g34 = new TH2F("g34","layer vs r",13,0.0,13.0, 40,-20.0,20.0);
TH2F *g35 = new TH2F("g35","r vs z", 160,-40.0,40.0, 80,0.0,20.0);
TH2F *g36 = new TH2F("g36","r vs z for good tracks", 160,-40.0,40.0, 80,0.0,20.0);

Int_t ngood = 0;
//iev = 500;                              //xxx <10K for testing
cout <<"Now processing "<<iev<<" events:"<<endl;
for (Int_t ii=1; ii<iev; ii++) {         // loop over events
  for (Int_t i=1; i<7; i++) {            // clear track and vertex arrays
    vertex[1][i]=0;
    vertex[2][i]=0;
    for (Int_t j=1; j<13;j++) {
      tracks[1][j][i]=0.;
      tracks[2][j][i]=0.;
    }
  }

  int ntracks = ievent[ii+1]-ievent[ii];
  if (ii%100 == 0) {cout<<" event "<<ii<<", "<<" number of hits on tracks "<<ntracks;}
  if (abs(ntracks) >20) {                                                   // There are a few with huge
    cout <<" number of hits on tracks "<<ntracks <<" -> skipped "<<endl;        // neg values, followed by huge pos.
    continue;                                                        // values of ntracks.
  }                                                              // Skip this event
  else if (ii%100 == 0) {
    cout <<endl;
  }

  //work on each event's hits, filling tracks and vertex arrays
  for (Int_t i=ievent[ii]; i<ievent[ii+1]; i++) {
   itrack=track[i];
    ilayer=layer[i];
    if (itrack>2) {                      // there could be 2 muons!
      cout <<" Track " <<itrack <<" ";
      cout <<" layer " <<ilayer <<" (too many tracks) "<<endl;
	break;
    }
    // store track info 
    idpart[itrack]=id[i];
    idp=id[i];
    tracks[itrack][ilayer][1] = x[i];
    tracks[itrack][ilayer][2] = y[i];
    tracks[itrack][ilayer][3] = z[i];

    //printf("%d %d %d %d %6.1f = i,itrack,ilayer,idp,z[i]\n",i,itrack,ilayer,idp,z[i]);
    g31->Fill(layer[i],x[i]);
    g32->Fill(layer[i],y[i]);
    g33->Fill(layer[i],z[i]);
    g34->Fill(layer[i],sqrt(x[i]**2+y[i]**2));
    g35->Fill(z[i],sqrt(x[i]**2+y[i]**2));
    tracks[itrack][ilayer][4] = px[i];
    tracks[itrack][ilayer][5] = py[i];
    tracks[itrack][ilayer][6] = pz[i];
    // store vertex info
    pzvertex = ptot[i]*cos(ptheta[i]*degtorad);
    ptvertex = ptot[i]*sin(ptheta[i]*degtorad);
    pxvertex = ptvertex*cos(pphi[i]*degtorad);
    pyvertex = ptvertex*sin(pphi[i]*degtorad);
    //vertex[itrack][1] = vx[ii][6];            note vx,y,z are note defined
    //vertex[itrack][2] = vy[ii][6];            or set
    //vertex[itrack][3] = vz[ii][6];
    vertex[itrack][4] = pxvertex;
    vertex[itrack][5] = pyvertex;
    vertex[itrack][6] = pzvertex;
  }

  //-------------------------------------------------------------------------//

  Float_t xintercept[3];
  Float_t yintercept[3];
  Float_t xvtx = 0.;
  Float_t yvtx = 0.;
  Float_t xslope[3];
  Float_t yslope[3];
  Float_t xpairvtx = 0.;
  Float_t ypairvtx = 0.;
  Float_t pairvtx = 0.;
  Float_t rpairvtx = 0.;
  Float_t rpairint = 0.;
  Float_t pairmass = 0.;
  Float_t pairmom = 0.;
  Float_t pairpt = 0.;
  Float_t pairetot = 0.;
  Float_t etot[3];
  // xseg also used for radius, yseg also used for phi*radius
  Float_t xseg = 0.00500;
  Float_t yseg = 0.02000;
  Float_t xmatch = 0.;
  Float_t ymatch = 0.;
  Float_t pzmatch = 0.;
  Float_t deltaphi = 0.;
  Float_t radius[5];
  Float_t phidet[5];
  Float_t rslope[3];
  Float_t rintercept[3];
  Float_t rvtx = 0.;
  Float_t rmatch = 0.;
  Float_t radvtx = 0.;
  Float_t vtxslope = 0.;
  Int_t goodtrk[3];
  Int_t xstrip = 0;
  Int_t ystrip = 0;
  Int_t rstrip = 0;
  Int_t phistrip = 0;

  Float_t fitx[4] = {13.0,21.0,29.0,37.0};
  Float_t fity[4];
  Float_t ex[4] = {.0001,.0001,.0001,.0001};
  Float_t ey[4] = {.001,.001,.001,.001};

  // loop over tracks in each event
  for (Int_t i=1; i<3;i++) {
    xslope[i]=0.;
    yslope[i]=0.;
    xintercept[i]=0.;
    yintercept[i]=0.;
    rintercept[i]=0.; 
    rslope[i]=0.; 
    goodtrk[i]=0;

    // require hits in layers 5, 6, 7, 8 (= Si N) and  pz>2.5 GeV   or
    //                 layers 9,10,11,12 (= Si S) and  pz>2.5 GeV  
    //printf("%6.1f %6.1f %6.1f %6.1f %6.1f %6.1f %6.1f %6.1f %6.1f \n", 
    //tracks[i][5][3],tracks[i][6][3],tracks[i][7][3], 
    //tracks[i][8][3],tracks[i][9][3],tracks[i][10][3],tracks[i][11][3],tracks[i][12][3],vertex[i][6]);
    if ((tracks[i][5][3]!=0. && tracks[i][ 6][3]!=0. && tracks[i][ 7][3]!=0. && tracks[i][ 8][3]!=0. && vertex[i][6]>2.5)  ){
//      (tracks[i][9][3]!=0. && tracks[i][10][3]!=0. && tracks[i][11][3]!=0. && tracks[i][12][3]!=0. && vertex[i][6]>2.5))  {
      goodtrk[i]=1;
      //cout<<"****                              *** Track "<<i<<" good"<<endl;
      ngood++;
      g36->Fill(tracks[i][5][3],sqrt(tracks[i][5][2]**2+tracks[i][5][1]**2));
      g36->Fill(tracks[i][6][3],sqrt(tracks[i][6][2]**2+tracks[i][6][1]**2));
      g36->Fill(tracks[i][7][3],sqrt(tracks[i][7][2]**2+tracks[i][7][1]**2));
      g36->Fill(tracks[i][8][3],sqrt(tracks[i][8][2]**2+tracks[i][8][1]**2));
//xxx OK to here
      //      compute x and y projected vertex                             xxx 3,4,5,6 -> 5,6,7,8
      //      cout<<tracks[i][3][1]<<" x at plane 3, z = "<<tracks[i][3][3]<<endl;
      //      cout<<tracks[i][3][2]<<" y at plane 3"<<endl;
      //      cout<<tracks[i][6][1]<<" x at plane 6, z = "<<tracks[i][6][3]<<endl;
      //      cout<<tracks[i][6][2]<<" y at plane 6"<<endl;
      // integerize strip positions for planes
      //      xstrip = (tracks[i][3][1]/xseg + 0.5);
      //      ystrip = (tracks[i][3][2]/yseg + 0.5);
      //      tracks[i][3][1] = xseg*xstrip;
      //      tracks[i][3][2] = yseg*ystrip;
      //      xstrip = (tracks[i][4][1]/xseg + 0.5);
      //      ystrip = (tracks[i][4][2]/yseg + 0.5);
      //      tracks[i][4][1] = xseg*xstrip;
      //      tracks[i][4][2] = yseg*ystrip;
      //      xstrip = (tracks[i][5][1]/xseg + 0.5);
      //      ystrip = (tracks[i][5][2]/yseg + 0.5);
      //      tracks[i][5][1] = xseg*xstrip;
      //      tracks[i][5][2] = yseg*ystrip;
      //      xstrip = (tracks[i][6][1]/xseg + 0.5);
      //      ystrip = (tracks[i][6][2]/yseg + 0.5);
      //      tracks[i][6][1] = xseg*xstrip;
      //      tracks[i][6][2] = yseg*ystrip;

      // calculate track slopes and intercepts at Z=0
      // try using first two planes to limit multiple scattering
      //      if ((tracks[i][6][3]-tracks[i][3][3])!=0) xslope[i] = (tracks[i][6][1]-tracks[i][3][1]) / (tracks[i][6][3]-tracks[i][3][3]);

      //xxx new layer numbers 3->05 4->07  (1st,2nd -> 1st,3rd)
      // use first and last planes to minimize resolution effects
      if ((tracks[i][07][3]-tracks[i][05][3])!=0) 
         xslope[i] = (tracks[i][07][1]-tracks[i][05][1]) / (tracks[i][07][3]-tracks[i][05][3]);
      xintercept[i] = tracks[i][05][1] - xslope[i] * tracks[i][05][3];
      xvtx=-999.;
      if (xslope[i]!=0.) xvtx = -xintercept[i]/xslope[i];
      //      cout<<xintercept[i]<<" xintercept"<<endl;
      if (idpart[i]==5) g1->Fill(xintercept[i]);
      if (idpart[i]==6) g3->Fill(xintercept[i]);
      if (idpart[i]==5) g5->Fill(xvtx);
      if (idpart[i]==6) g7->Fill(xvtx);
      if ((tracks[i][07][3]-tracks[i][05][3])!=0) 
         yslope[i] = (tracks[i][07][2]-tracks[i][05][2]) / (tracks[i][07][3]-tracks[i][05][3]);
      //      if ((tracks[i][6][3]-tracks[i][3][3])!=0) yslope[i] = (tracks[i][6][2]-tracks[i][3][2]) / (tracks[i][6][3]-tracks[i][3][3]);
      yintercept[i] = tracks[i][05][2] - yslope[i] * tracks[i][05][3];
      yvtx = -999.;      
      if (yslope[i]!=0.) yvtx = -yintercept[i]/yslope[i];
      //      cout<<yintercept[i]<<" yintercept"<<endl;
      if (idpart[i]==5) g2->Fill(yintercept[i]);
      if (idpart[i]==6) g4->Fill(yintercept[i]);
      if (idpart[i]==5) g6->Fill(yvtx);
      if (idpart[i]==6) g8->Fill(yvtx);
      if (goodtrk[1]==1 && goodtrk[2]==1) {
	g24->Fill(xintercept[1],xintercept[2]);
      	g25->Fill(yintercept[1],yintercept[2]);
	g26->Fill(xintercept[1],yintercept[1]);
      	g27->Fill(xintercept[2],yintercept[2]);
      }
      // radial quantities
      // use xseg as radial segmentation, yseg as phi 
      if ((tracks[i][08][3]-tracks[i][05][3])!=0) {
        radius[1] = (tracks[i][05][1]**2 + tracks[i][05][2]**2)**0.5;
        radius[2] = (tracks[i][06][1]**2 + tracks[i][06][2]**2)**0.5;
        radius[3] = (tracks[i][07][1]**2 + tracks[i][07][2]**2)**0.5;
        radius[4] = (tracks[i][08][1]**2 + tracks[i][08][2]**2)**0.5;
        phidet[1] = atan(tracks[i][05][2]/tracks[i][05][1]);                      //xxx
        phidet[2] = atan(tracks[i][06][2]/tracks[i][06][1]);
        phidet[3] = atan(tracks[i][07][2]/tracks[i][07][1]);
        phidet[4] = atan(tracks[i][08][2]/tracks[i][08][1]);

	// integerize radius
 	rstrip = (radius[1]/xseg + 0.5);
	radius[1] = xseg*rstrip;
	rstrip = (radius[2]/xseg + 0.5);
	radius[2] = xseg*rstrip;
        rstrip = (radius[3]/xseg + 0.5);
        radius[3] = xseg*rstrip;
        rstrip = (radius[4]/xseg + 0.5);
        radius[4] = xseg*rstrip;
        // integerize phi using arc length
	phistrip = (phidet[1]*radius[1]/yseg + 0.5);
	phidet[1] = yseg*phistrip/radius[1];
	phistrip = (phidet[2]*radius[2]/yseg + 0.5);
	phidet[2] = yseg*phistrip/radius[2];
        phistrip = (phidet[3]*radius[3]/yseg + 0.5);
        phidet[3] = yseg*phistrip/radius[3];
        phistrip = (phidet[4]*radius[4]/yseg + 0.5);
        phidet[4] = yseg*phistrip/radius[4];

        rslope[i] = (radius[4]-radius[1]) / (tracks[i][08][3]-tracks[i][05][3]);
        rintercept[i] = radius[1] - rslope[i] * tracks[i][05][3];
        if (rslope[i]!=0.) rvtx = -rintercept[i]/rslope[i];
        vtxslope = ((vertex[i][4]**2+vertex[i][5]**2)**0.5)/vertex[i][6];
        rmatch = radius[1] - vtxslope*tracks[i][05][3];
        //        rmatch = rmatch + (radius[2] - vtxslope*tracks[i][6][3]);
        //        cout<<rslope<<rvtx<<endl; 
        g19->Fill(rintercept[i]);
        g20->Fill(rvtx);
        g21->Fill(rmatch);
        radvtx = ((vertex[i][1]**2+vertex[i][2]**2)**0.5);   //xxx note vertex[][1,2,3] are never set 
        g22->Fill(radvtx);                                   //xxx radvtx not properly defined
        g23->Fill(vertex[i][3]);                             //xxx vertex[][1,2,3] never set
// fit the radial intercept
//        fity[0]=radius[1];
//	fity[1]=radius[2];
//	fity[2]=radius[3];
//	fity[3]=radius[4];
//	gr = new TGraphErrors(4,fitx,fity,ex,ey);
//	gr->Fit("pol1");
//	TF1 *fit = gr->GetFunction("pol1");
//	Double_t chi2 = fit->GetChisquare();
//	Double_t p0 = fit->GetParameter(0);
//	Double_t e0 = fit->GetParError(0);
//	cout << "intercept fit " << p0 << " chi2 "<< chi2 << endl;
//        g28->Fill(p0);
      }
      // calculate track match from MuTr to last silicon plane
      //      xmatch = tracks[i][11][1] - (tracks[i][11][4]/tracks[i][11][6])*(tracks[i][11][3]-tracks[i][6][3]) - tracks[i][6][1];
      //ymatch = tracks[i][11][2] - (tracks[i][11][5]/tracks[i][11][5])*(tracks[i][11][3]-tracks[i][6][3]) - tracks[i][6][2];
      //pzmatch = tracks[i][6][6]-tracks[i][11][6];
      //g13->Fill(xmatch);
      //g14->Fill(ymatch);
      //g15->Fill(pzmatch);
      // calculate bend in phi, due to magnetic field and multiple scattering
      if (tracks[i][08][1]!=0. && tracks[i][05][1]!=0.) {
        deltaphi = phidet[4] - phidet[1];
        deltaphi = deltaphi/degtorad;
	if (idpart[i]==5) g16->Fill(deltaphi);
        if (idpart[i]==6) g17->Fill(deltaphi);
        deltaphi=deltaphi*tracks[i][6][6];
	g18->Fill(deltaphi);
      }
    }
  }

  // calculate Z-vertex and kinematics for good track pairs
  if (goodtrk[1]==1&&goodtrk[2]==1) {
    cout<<"good pair "<<endl;
    if ((xslope[2]-xslope[1])!=0) xpairvtx = -(xintercept[2]-xintercept[1])/(xslope[2]-xslope[1]);
    if ((yslope[2]-yslope[1])!=0) ypairvtx = -(yintercept[2]-yintercept[1])/(yslope[2]-yslope[1]);
    rpairvtx = -(rintercept[1]/rslope[1] + rintercept[2]/rslope[2]) / 2. ;
    rpairint = (rintercept[1] + rintercept[2]) / 2. ;
    pairvtx = (xpairvtx+ypairvtx)/2.;
    //    pairpt = (tracks[1][3][4] + tracks[2][3][4])**2 + (tracks[1][3][5] + tracks[2][3][5])**2;
    pairpt = (xslope[1]*tracks[1][05][6] + xslope[2]*tracks[2][05][6])**2 + (yslope[1]*tracks[1][05][6] + yslope[2]*tracks[2][05][6])**2;
    pairpt = pairpt**0.5;    
    pairmom = (pairpt**2 + (tracks[1][05][6] + tracks[2][05][6])**2)**0.5;
    //    etot[1] = (0.10566**2 + tracks[1][3][4]**2 + tracks[1][3][5]**2 + tracks[1][3][6]**2)**0.5;
    //    etot[2] = (0.10566**2 + tracks[2][3][4]**2 + tracks[2][3][5]**2 + tracks[2][3][6]**2)**0.5;
    etot[1] = (0.10566**2 + (xslope[1]*tracks[1][3][6])**2 + (yslope[1]*tracks[1][3][6])**2 + tracks[1][3][6]**2)**0.5;
    etot[2] = (0.10566**2 + (xslope[2]*tracks[2][3][6])**2 + (yslope[2]*tracks[2][3][6])**2 + tracks[2][3][6]**2)**0.5;
    pairmass = ((etot[1]+etot[2])**2 - pairmom**2)**0.5;
    g9->Fill(xpairvtx);
    g10->Fill(ypairvtx);
    g11->Fill(pairvtx);
    g12->Fill(pairmass);
    g29->Fill(rpairvtx);
    g30->Fill(rpairint);
  }
}
cout<<"good tracks: "<<ngood<<endl;
// done
}