{//////////////////////////////////////////////////////////
//   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
//   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. 
//////////////////////////////////////////////////////////
//Reset ROOT and connect tree file
   gROOT->Reset();
   TFile *f = new TFile("ancsvx_softlink.root");  // soft link to ancsvx_nnn.root

//Declaration of leaves types
//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 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
//xxx

Float_t event[100000];
Float_t track[100000];
Float_t ptot[100000];
Float_t ptheta[100000];
Float_t pphi[100000];
Float_t theta[100000];
Float_t phi[100000];
Float_t x[100000];
Float_t y[100000];
Float_t z[100000];
Float_t px[100000];
Float_t py[100000];
Float_t pz[100000];
Float_t layer[100000];
Float_t id[100000];
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++) {
    AncSvx->GetEvent(i);
    event[i]=EVENT;           // watch for new event number:
    if (EVENT>oldevent) {
      oldevent=EVENT;
      iev=EVENT;
      ievent[iev]=i; 
      if (iev%1000==0) {cout <<"record "<<i<<", event "<<iev<<" "<<endl;}
    }
    track[i]=TRACK;           // fill hit arrays
    layer[i]=LAYER;
    ptot[i]=PTOT;
    //    cout<< "ptot: "<<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;
    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 loop to stuff arrays
  cout <<"Total records "<<i<<", events "<<iev<<" "<<endl;
  //delete f;    // close the file xxx this causes a segfault

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

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 *g16 = new TH1F("g16","si dphi deg id=5",50,-1.,1.);
TH2F *g17 = new TH2F("g17","r slope vs intercept",50,-1. ,1., 50, -1., 1.);
TH1F *g18 = new TH1F("g18","r-slope ",50,0.,0.5);
TH1F *g19 = new TH1F("g19","rintercept",50,-0.2,0.2);
TH1F *g20 = new TH1F("g20","zrvertex",75,-.25, 0.50);
TH1F *g21 = new TH1F("g21","Si rmatch",50,-0.02,0.02);
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);
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);
TH1F *g37 = new TH1F("g37","pt",80,0.0,8.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","event hit multiplicity",50,0.0,50.0);
TH1F *g41 = new TH1F("g41","particle ID",20,-0.5,19.5);

Int_t ngood = 0;
Int_t pa, pb;     // track extrapolations are done from planes a and b
pa = 10;           // 9 is the first Si plane on the North side
pb = 12;           // 11 is the 3rd (out of 4) Si plane xxx      
cout<<"Tracks made from hits in layers "<<pa<<" and "<<pb
    <<" (Si planes N-"<<pa-8<<" and N-"<<pb-8<<")."<<endl;

int nmess = 100;                          // messages every nmess events (typically 1 or 100)
//iev = 10;                              //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.;
    }
  }

  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
    if (iev<20) 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++) {
    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];

      // 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 etot[3];
  // rseg also used for radius, pseg also used for phi*radius
  //Float_t rseg   = 0.0075;     // 75 microns -> cm 
  Float_t rseg   = 0.0050;     // 50 microns for the ifvtx 
  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  

//  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)  ){
    if ((tracks[i][9][3]!=0. && tracks[i][10][3]!=0. && tracks[i][11][3]!=0. && 
         tracks[i][12][3]!=0)  ){

      goodtrk[i]=1;
      cout<<"****                              *** Track "<<i<<" good"<<endl;
      ngood++;
      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      xxx note planes 3,4,5,6 -> 5,6,7,8
                      // 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 ((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 (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);    // rstrip is Int
        rsave = radius[1];
	radius[1] = rseg*rstrip;
        //cout<<"Layer 1 old, new r: "<<rsave<<" "<<radius[1]<<endl;

	rstrip = (radius[2]/rseg + 0.5);
        rsave = radius[2];
	radius[2] = rseg*rstrip;

        rstrip = (radius[3]/rseg + 0.5);
        rsave = radius[3];
        radius[3] = rseg*rstrip;

        rstrip = (radius[4]/rseg + 0.5);
        rsave = radius[4];
        radius[4] = rseg*rstrip;

        // 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<20) cout<<"rvtx= "<<rvtx<<endl;
          g20->Fill(rvtx);              // since this is centered on 0.0, 
          float pt = sqrt(vertex[i][4]**2 + vertex[i][5]**2);
          g37->Fill(pt);
        }
        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; 
        g17->Fill(rslope[i],rintercept[i]);
        g18->Fill(rslope[i]);
        g19->Fill(rintercept[i]);
        g21->Fill(rmatch);
 

     }      // endif protect against divide by 0

      }
    }            // good track
  }              // end loop over tracks in each event
}    // end for 
cout<<"good tracks: "<<ngood<<endl;
TFile *outfile = new  TFile("beauty03.root","RECREATE");  
g7->Write(); g8->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();
g31->Write(); g32->Write(); g33->Write(); g34->Write(); g35->Write(); 
g36->Write(); g37->Write(); g39->Write(); g40->Write();
g41->Write(); 
outfile ->Close();
// done
}