{
//////////////////////////////////////////////////////////
//   This file has been automatically generated 
//     (Wed Dec  8 19:23:01 2004 by ROOT version4.01/02)
//   from TTree AncSvx/SVX Ancestors
//   found on file: ancsvx_softlink.root (=ancsvx_53.root)
// How to: $ root -l ancsvx_softlink.root
// root[1] AncSvx->MakeCode("charm02.C")
// Dec 2004
//////////////////////////////////////////////////////////
// The analysis code is 
// derived from Pat's Charm2.C Jun 04 HvH
//   Major changes:
// - allow for tilted planes (see 'integerization' section)
// - read in secondary vertex distribution from flat file
//   flat_xyz.out, which comes from Pythia.
// 18 jun 04 switch to all_keep = .true. in pythia. 
// paw not available, disable this part for now Dce 9 2004
//////////////////////////////////////////////////////////

//Reset ROOT and connect tree file
   gROOT->Reset();
   TFile *f = (TFile*)gROOT->GetListOfFiles()->FindObject("ancsvx_softlink.root");
   if (!f) {
      f = new TFile("ancsvx_softlink.root");
   }
   TTree *AncSvx = (TTree*)gDirectory->Get("AncSvx");

//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         X0_EVENT;
   Float_t         Y0_EVENT;
   Float_t         Z0_EVENT;
   Float_t         B_IMPACT;
   Float_t         EVENT;

   // Set branch addresses.
   AncSvx->SetBranchAddress("TRACK",&TRACK);
   AncSvx->SetBranchAddress("NFILE",&NFILE);
   AncSvx->SetBranchAddress("PTOT",&PTOT);
   AncSvx->SetBranchAddress("PTHETA",&PTHETA);
   AncSvx->SetBranchAddress("PPHI",&PPHI);
   AncSvx->SetBranchAddress("R_VERTEX",&R_VERTEX);
   AncSvx->SetBranchAddress("Z_VERTEX",&Z_VERTEX);
   AncSvx->SetBranchAddress("THET_VER",&THET_VER);
   AncSvx->SetBranchAddress("PHI_VER",&PHI_VER);
   AncSvx->SetBranchAddress("ITPARENT",&ITPARENT);
   AncSvx->SetBranchAddress("IDPARENT",&IDPARENT);
   AncSvx->SetBranchAddress("IDPART",&IDPART);
   AncSvx->SetBranchAddress("ITORIGIN",&ITORIGIN);
   AncSvx->SetBranchAddress("IDORIGIN",&IDORIGIN);
   AncSvx->SetBranchAddress("IHIT",&IHIT);
   AncSvx->SetBranchAddress("LAYER",&LAYER);
   AncSvx->SetBranchAddress("THETA",&THETA);
   AncSvx->SetBranchAddress("PHI",&PHI);
   AncSvx->SetBranchAddress("XGLOBAL",&XGLOBAL);
   AncSvx->SetBranchAddress("YGLOBAL",&YGLOBAL);
   AncSvx->SetBranchAddress("ZGLOBAL",&ZGLOBAL);
   AncSvx->SetBranchAddress("PMOMX",&PMOMX);
   AncSvx->SetBranchAddress("PMOMY",&PMOMY);
   AncSvx->SetBranchAddress("PMOMZ",&PMOMZ);
   AncSvx->SetBranchAddress("IPC123",&IPC123);
   AncSvx->SetBranchAddress("PC1THETA",&PC1THETA);
   AncSvx->SetBranchAddress("PC1PHI",&PC1PHI);
   AncSvx->SetBranchAddress("PC1RDIST",&PC1RDIST);
   AncSvx->SetBranchAddress("PC1ZDIST",&PC1ZDIST);
   AncSvx->SetBranchAddress("PC2THETA",&PC2THETA);
   AncSvx->SetBranchAddress("PC2PHI",&PC2PHI);
   AncSvx->SetBranchAddress("PC2RDIST",&PC2RDIST);
   AncSvx->SetBranchAddress("PC2ZDIST",&PC2ZDIST);
   AncSvx->SetBranchAddress("PC3THETA",&PC3THETA);
   AncSvx->SetBranchAddress("PC3PHI",&PC3PHI);
   AncSvx->SetBranchAddress("PC3RDIST",&PC3RDIST);
   AncSvx->SetBranchAddress("PC3ZDIST",&PC3ZDIST);
   AncSvx->SetBranchAddress("DELE",&DELE);
   AncSvx->SetBranchAddress("XLOCIN",&XLOCIN);
   AncSvx->SetBranchAddress("YLOCIN",&YLOCIN);
   AncSvx->SetBranchAddress("ZLOCIN",&ZLOCIN);
   AncSvx->SetBranchAddress("XLOCOUT",&XLOCOUT);
   AncSvx->SetBranchAddress("YLOCOUT",&YLOCOUT);
   AncSvx->SetBranchAddress("ZLOCOUT",&ZLOCOUT);
   AncSvx->SetBranchAddress("HITVOL0",&HITVOL0);
   AncSvx->SetBranchAddress("HITVOL1",&HITVOL1);
   AncSvx->SetBranchAddress("HITVOL2",&HITVOL2);
   AncSvx->SetBranchAddress("HITVOL3",&HITVOL3);
   AncSvx->SetBranchAddress("HITVOL4",&HITVOL4);
   AncSvx->SetBranchAddress("HITVOL5",&HITVOL5);
   AncSvx->SetBranchAddress("NHIT",&NHIT);
   AncSvx->SetBranchAddress("X0_EVENT",&X0_EVENT);
   AncSvx->SetBranchAddress("Y0_EVENT",&Y0_EVENT);
   AncSvx->SetBranchAddress("Z0_EVENT",&Z0_EVENT);
   AncSvx->SetBranchAddress("B_IMPACT",&B_IMPACT);
   AncSvx->SetBranchAddress("EVENT",&EVENT);

//     This is the loop skeleton
//       To read only selected branches, Insert statements like:
// AncSvx->SetBranchStatus("*",0);  // disable all branches
// TTreePlayer->SetBranchStatus("branchname",1);  // activate branchname
// MakeCode initially declared as  Long64_t nentries = AncSvx->GetEntries();
// this cuased a fatal error (malloc failure): 
   Int_t nentries = AncSvx->GetEntries();
   Int_t nbytes = 0;

  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 pt[nentries];
Float_t layer[nentries];
Float_t id[nentries];
Int_t ievent[10001];

  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++) {
      nbytes += AncSvx->GetEntry(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
    layer[i]=LAYER;

    ptot[i]=PTOT;
    //cout<<ptot[i]<<endl;

    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;
    pt[i]=sqrt(PMOMX**2+PMOMY**2);
    pz[i]=PMOMZ;
    id[i]=IDPART;  

    //if (layer[i]>8) {              // mirror South into North arm xxx
    //  layer[i] = layer[i] - 4;     // 9-12 -> 5-8
    //  x[i] =  -x[i];  y[i] =  -y[i];  z[i] =  -z[i]; 
    //  px[i] = -px[i]; py[i] = -py[i]; pz[i] = -pz[i];
    //  ptot[i]   = 
    //  ptheta[i] = 
    //  pphi[i]   = 
    //  theta[i]  = 
    //  phi[i]    =  
    //}
   }         // end of loop to stuff arrays
  cout <<"Total records "<<i<<", events "<<iev<<" "<<endl;
  f->Close();
//-----------------------------------------------------------------------------

const int max_tracks=200;
Float_t tracks[max_tracks][13][7];    // track#(1-2) layer(1-12) (x,y,z,px,py,pz)
Int_t idpart  [max_tracks];           // pid   (1-2)
Float_t vertex[max_tracks][7];        // track#(1-2) (xvertex,yvertex,zvertex,ptot,ptheta,pphi)

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.;

//define histograms
TH1F *g1 = new TH1F("g1","xintercept id=5",50,-0.10,0.10);
TH1F *g2 = new TH1F("g2","yintercept id=5",50,-0.10,0.10);
TH1F *g3 = new TH1F("g3","xintercept id=6",50,-0.10,0.10);
TH1F *g4 = new TH1F("g4","yintercept id=6",50,-0.10,0.10);
TH1F *g5 = new TH1F("g5","zxvertex id=5",50,-0.4,0.4);
TH1F *g6 = new TH1F("g6","zyvertex id=5",50,-0.4,0.4);
TH1F *g7 = new TH1F("g7","zxvertex id=6",50,-0.4,0.4);
TH1F *g8 = new TH1F("g8","zyvertex id=6",50,-0.4,0.4);
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",75,-.25, 0.50);
TH1F *g21 = new TH1F("g21","Si rmatch",50,-0.02,0.02);
TH1F *g22 = new TH1F("g22","zrvertex + zpythia",75, -0.25, 0.50);
TH1F *g23 = new TH1F("g23","xpythia",75, -0.25, 0.50);
TH1F *g24 = new TH1F("g24","ypythia",75, -0.25, 0.50);
TH1F *g25 = new TH1F("g25","zpythia",75, -0.25, 0.50);

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,  0.0,30.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);
TH1F *g37 = new TH1F("g37","pt",80,0.0,8.0);
TH2F *g38 = new TH2F("g38","pT vs vertex ",75, -0.25, 0.50, 80,0.0,4.0);
TH2F *g39 = new TH2F("g39","ID vs track number",5,-0.5,4.5,20,-0.5,19.5);
TH1F *g40 = new TH1F("g40","hits in this event",20,0.0,20.0);
TH1F *g41 = new TH1F("g41","particle ID",20,-0.5,19.5);
TH1F *g42 = new TH1F("g42","px",50,-10.0,10.0);
TH1F *g43 = new TH1F("g43","py",50,-5.0,5.0);
TH1F *g44 = new TH1F("g44","pz",50,-10.0,10.0);
TH1F *g45 = new TH1F("g45","ptot",50,0.0,10.0);
TH1F *g46 = new TH1F("g46","pt", 50, 0.,10.0);

Int_t pa, pb;     // track extrapolations are done from planes a and b
pa = 9;           // 9 is innermost on the N
pb = 12;           // there are 2x4 layers, south to North 5->12      
cout<<"Tracks will be made from hits in layers "<<pa<<" and "<<pb
    <<" (Si planes "<<pa-4<<" and "<<pb-4<<")."<<endl;

float xpythia, ypythia, zpythia;
//FILE* vertex_file = fopen("flat_xyz.out","r");
//cout <<"vertex_file="<<vertex_file<<endl;
//if (vertex_file == NULL) {  cout<<" error opening flat_xyz.out"<<endl; }
//else { 
//  cout<<" success opening flat_xyz.out"<<endl;
//  //fscanf(vertex_file,"%f%f%f",&xpythia,&ypythia,&zpythia);
//  //cout <<xpythia<<" "<<ypythia<<" "<<zpythia<<" "<<endl;
//}

int nmess = 100;                          // messages every nmess events (typically 1 or 100)
iev = 1000;                              //xxx <10K uncomment 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<=6; i++) {            // clear track and vertex arrays
    vertex[1][i]=0;
    vertex[2][i]=0;
    for (Int_t j=1; j<=12;j++) {
      tracks[1][j][i]=0.;
      tracks[2][j][i]=0.;
    }
  }
  xpythia = 0.0; ypythia = 0.0; zpythia = 0.0; 
  //  fscanf(vertex_file,"%f%f%f",&xpythia,&ypythia,&zpythia);        // read a vertex out of the other file
  //cout <<xpythia<<" "<<ypythia<<" "<<zpythia<<" "<<endl;
  g23->Fill(xpythia);
  g24->Fill(ypythia);
  g25->Fill(zpythia);
  int nhits = ievent[ii+1]-ievent[ii];
  g40->Fill(float(nhits));

  if (ii%nmess == 0) {cout<<" event "<<ii<<"/"<<iev<<", "<<" number of hits in this event "<<nhits;}
  if (abs(nhits) >200) {                                             // There are a few with huge
    cout <<" number of hits on tracks "<<nhits <<" -> skipped "<<endl; // neg values, followed by huge pos.
    continue;                                                          // values of nhits.
  }                                                                    // Skip this event
  else if (ii%nmess == 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++) {
    if (i == ievent[ii]) {            // once per event
      g42->Fill(px[i]);
      g43->Fill(py[i]);
      g44->Fill(pz[i]);
      g45->Fill(ptot[i]);
      g46->Fill(pz[i]);
    }
    itrack=track[i];
    ilayer=layer[i];
    //cout<<" itrack, ilayer= "<<itrack<<"  "<<ilayer<<endl;
    if (itrack>=200) {
      cout <<" Track " <<itrack <<" ";
      cout <<" layer " <<ilayer <<" (too many tracks) "<<endl;
      //break;   should be 'continue' ?
    }
    else {                                      // store track info 
      idpart[itrack]=id[i];
      g41->Fill(float(id[i]));

      tracks[itrack][ilayer][1] = x[i];
      tracks[itrack][ilayer][2] = y[i];
      tracks[itrack][ilayer][3] = z[i];
      tracks[itrack][ilayer][4] = px[i];
      tracks[itrack][ilayer][5] = py[i];
      tracks[itrack][ilayer][6] = pz[i];  
      //cout <<">>>itrack "<<itrack<<", layer "<<ilayer<<", tracks[itrack][ilayer][1,2,3]= "<<
      //tracks[itrack][ilayer][1]<<" "<<tracks[itrack][ilayer][2]<<" "<<tracks[itrack][ilayer][3]<<"."<<endl;

    // 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;

      //printf("%d %d %d %d %6.1f = i,itrack,ilayer,idpart[i],z[i]\n",i,itrack,ilayer,idpart[i],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));
    }      // <200 track in the event
  }

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

  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];
  // rseg also used for radius, pseg also used for phi*radius
  //Float_t rseg = 0.00500;
  Float_t rseg   = 0.00447;  // =50*cos(26.6 degrees)
  Float_t sin_alpha = 0.447;
  Float_t rsave;
  Float_t pseg = 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 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};
  Int_t id_track = 0;

  // loop over tracks in each event
  for (Int_t i=1; i<=1;i++) {               // one only!
    id_track = idpart[i];    
    g39->Fill(float(i),float(id_track));            // plot particle ID
    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  


    //cout <<"tracks[i][9,10,11,12]= "<<tracks[i][9][3]<<" "<<tracks[i][10][3]<<" "<<tracks[i][11][3]<<" "<<tracks[i][12][3]
    //<<" vertex[i][6]= "<<vertex[i][6]<<endl;
    if ((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;
      g36->Fill(tracks[i][ 9][3],sqrt(tracks[i][ 9][2]**2+tracks[i][ 9][1]**2));
      g36->Fill(tracks[i][10][3],sqrt(tracks[i][10][2]**2+tracks[i][10][1]**2));
      g36->Fill(tracks[i][11][3],sqrt(tracks[i][11][2]**2+tracks[i][11][1]**2));
      g36->Fill(tracks[i][12][3],sqrt(tracks[i][12][2]**2+tracks[i][12][1]**2));
      //      compute x and y projected vertex. note planes S->N are 5,6,7,8 and  9,10,11,12 
                      // calculate track slopes and intercepts at Z=0
                      // pick two planes a and b to limit the effects of multiple scattering
                      // Variables pa and pb set above ouside the event loop
      if ((tracks[i][pb][3]-tracks[i][pa][3])!=0) 
         xslope[i] = (tracks[i][pb][1]-tracks[i][pa][1]) / (tracks[i][pb][3]-tracks[i][pa][3]);
      xintercept[i] = tracks[i][pa][1] - xslope[i] * tracks[i][pa][3];
      xvtx=-999.;
      if (xslope[i]!=0.) xvtx = -xintercept[i]/xslope[i];
      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][pb][3]-tracks[i][pa][3])!=0) 
         yslope[i] = (tracks[i][pb][2]-tracks[i][pa][2]) / (tracks[i][pb][3]-tracks[i][pa][3]);
      yintercept[i] = tracks[i][pa][2] - yslope[i] * tracks[i][pa][3];
      yvtx = -999.;      
      if (yslope[i]!=0.) yvtx = -yintercept[i]/yslope[i];
      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) {
	g26->Fill(xintercept[1],yintercept[1]);
      	g27->Fill(xintercept[2],yintercept[2]);
      }
      // radial quantities
      // use rseg as radial segmentation, pseg as phi 
      if ((tracks[i][pb][3]-tracks[i][pa][3])!=0) {
        radius[1] = sqrt(tracks[i][09][1]**2 + tracks[i][09][2]**2);
        radius[2] = sqrt(tracks[i][10][1]**2 + tracks[i][10][2]**2);
        radius[3] = sqrt(tracks[i][11][1]**2 + tracks[i][11][2]**2);
        radius[4] = sqrt(tracks[i][12][1]**2 + tracks[i][12][2]**2);
        phidet[1] = atan(tracks[i][09][2]/tracks[i][09][1]);
        phidet[2] = atan(tracks[i][10][2]/tracks[i][10][1]);
        phidet[3] = atan(tracks[i][11][2]/tracks[i][11][1]);
        phidet[4] = atan(tracks[i][12][2]/tracks[i][12][1]);

	// integerize radius, and adjust z accordingly
 	rstrip = (radius[1]/rseg + 0.5);
        rsave = radius[1];
	radius[1] = rseg*rstrip;
        //cout<<"old, new r: "<<rsave<<" "<<radius[1]<<endl;
        //cout<<"old, new z: "<<tracks[i][9][3]<<" ";
        tracks[i][9][3] = tracks[i][9][3]+(radius[1]-rsave)*sin_alpha;
	//cout <<tracks[i][9][3]<<endl;

	rstrip = (radius[2]/rseg + 0.5);
        rsave = radius[2];
	radius[2] = rseg*rstrip;
        tracks[i][10][3] = tracks[i][10][3]+(radius[2]-rsave)*sin_alpha;

        rstrip = (radius[3]/rseg + 0.5);
        rsave = radius[3];
        radius[3] = rseg*rstrip;
        tracks[i][11][3] = tracks[i][11][3]+(radius[3]-rsave)*sin_alpha;

        rstrip = (radius[4]/rseg + 0.5);
        rsave = radius[4];
        radius[4] = rseg*rstrip;
        tracks[i][12][3] = tracks[i][12][3]+(radius[4]-rsave)*sin_alpha;

        // integerize phi using arc length
	phistrip = (phidet[1]*radius[1]/pseg + 0.5);
	phidet[1] = pseg*phistrip/radius[1];
	phistrip = (phidet[2]*radius[2]/pseg + 0.5);
	phidet[2] = pseg*phistrip/radius[2];
        phistrip = (phidet[3]*radius[3]/pseg + 0.5);
        phidet[3] = pseg*phistrip/radius[3];
        phistrip = (phidet[4]*radius[4]/pseg + 0.5);
        phidet[4] = pseg*phistrip/radius[4];
                                                            // rz:    
        rslope[i] = (radius[pb-8]-radius[pa-8]) / (tracks[i][pb][3]-tracks[i][pa][3]);
        rintercept[i] = radius[pa-8] - rslope[i] * tracks[i][pa][3];
        if (rslope[i]!=0.) {
          rvtx = -rintercept[i]/rslope[i];
          if (ii<50) cout<<"rvtx= "<<rvtx<<endl;
          g20->Fill(rvtx);              // since this is centered on 0.0, 
          g22->Fill(rvtx+zpythia);      // place fills inside the 'if'!
          g23->Fill(xpythia);
          g24->Fill(ypythia);
          g25->Fill(zpythia);
          float p_t = sqrt(vertex[i][4]**2 + vertex[i][5]**2);
          g37->Fill(p_t);
          g38->Fill(rvtx+zpythia,p_t);
        }
        vtxslope = ((vertex[i][4]**2+vertex[i][5]**2)**0.5)/vertex[i][6];
        rmatch = radius[1] - vtxslope*tracks[i][pa][3];
        //        rmatch = rmatch + (radius[2] - vtxslope*tracks[i][6][3]);
        //        cout<<rslope<<rvtx<<endl; 
        g19->Fill(rintercept[i]);
        g21->Fill(rmatch);

      }      // endif protect against divide by 0

      // 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);
      }
    }            // good track
  }              // end loop over tracks in each event

  // 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][pa][6] + xslope[2]*tracks[2][pa][6])**2 + (yslope[1]*tracks[1][pa][6] + yslope[2]*tracks[2][pa][6])**2;
    pairpt = pairpt**0.5;    
    pairmom = (pairpt**2 + (tracks[1][pa][6] + tracks[2][pa][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);
  }
}
TFile *outfile = new  TFile("beauty03.root","RECREATE");  
g1->Write(); g2->Write(); g3->Write(); g4->Write(); g5->Write();
g6->Write(); g7->Write(); g8->Write(); g9->Write(); g10->Write();
g11->Write(); g12->Write(); g13->Write(); g14->Write(); g15->Write();
g16->Write(); g17->Write(); g18->Write(); g19->Write(); g20->Write();
g21->Write(); g22->Write(); g23->Write(); g24->Write(); g25->Write();
g26->Write(); g27->Write(); g28->Write(); g29->Write(); g30->Write(); 
g31->Write(); g32->Write(); g33->Write(); g34->Write(); g35->Write(); 
g36->Write(); g37->Write(); g38->Write(); g39->Write(); g40->Write();
g41->Write(); g42->Write(); g43->Write(); g44->Write(); g45->Write(); 
g46->Write(); //g37->Write(); g38->Write(); g39->Write(); g40->Write();
outfile ->Close();
// done

}           // charm02.C script