{
//////////////////////////////////////////////////////////
//   This file has been automatically generated 
//     (Fri Apr  8 17:55:21 2005 by ROOT version4.01/02)
//   from TTree AncSvx/SVX Ancestors
//   found on file: ancsvx_057.root
//   do:
//  root -l ancsvx_057.root
//  .ls shows an Ntuple: AncSvx;1
//  AncSvx->MakeCode("macroname.C");
//
//////////////////////////////////////////////////////////


//Reset ROOT and connect tree file
   gROOT->Reset();
   TFile *f = (TFile*)gROOT->GetListOfFiles()->FindObject("ancsvx_softlink.root");
   if (!f) {
     cout<<" using ancsvx_softlink.root "<<endl;
      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

   Long64_t nentries = AncSvx->GetEntries();

  Int_t oldevent = 1;
  Int_t mult0    = 0;
  Int_t mult1    = 0;
  Int_t oldtrack = 0;
  Float_t oldz   = 0;
  Float_t delz   = 0;
  Int_t nbytes   = 0;
  Float_t phi    = 0.0;
  Float_t phi2   = 0.0;
  Float_t r      = 0.0;
  Float_t eta, x2, y2, sign;
  Float_t pi=  3.14159;
  Float_t hit[8000][13][7], track_max, xa[4], xb[4], ra, rb, xs1, ys1, zs1, rs1;
  Float_t scale, zv, rglobal, zvertex, phia, phib;
  Int_t itrack, lnhit, ilayer, nlhit, ixyz, icount;
                             // for integerization:
  Float_t rseg   = 0.00447;  // =50*cos(26.6 degrees)
  Int_t irstrip = 0;
  Float_t sin_alpha = 0.447;
  Float_t multistrip;

TH1F *e01 = new TH1F("e01","records per event",50, 0.0, 50.0);
TH1F *e02 = new TH1F("e02","phi track 1",50, -1.0, 7.0);
TH1F *e03 = new TH1F("e03","phi track 2",50, -1.0, 7.0);
TH1F *e04 = new TH1F("e04","phia-phib",  70, -7.0, 7.0);
TH1F *e04b = new TH1F("e04b","phia-phib",  70, -0.01, 0.01);
TH2F *e05 = new TH2F("e05","phia vs. phib",  50, -1.0, 7.0, 50, -1.0, 7.0);

TH1F *g01 = new TH1F("g01","xlocal layer 8",50,-1.0, 1.0);
TH1F *g02 = new TH1F("g02","ylocal layer 8",50,-0.025, 0.025);
TH1F *g03 = new TH1F("g03","zlocal layer 8",50,-5.0, 5.0);
TH2F *g04 = new TH2F("g04","z-x local layer 8",50,-5.0, 5.0, 50, -5.0, 5.0);
TH2F *g05 = new TH2F("g05","x-y global layer 8", 40, -20.0, 20.0, 40, -20.0, 20.0);
TH2F *g05a = new TH2F("g05a","x-y global layer 8", 24, -12.0, 12.0, 24, -12.0, 12.0);
TH2F *g06 = new TH2F("g06","r-z global layer 8", 15, -25.0, -10.0, 15, 0.0, 15.0);

TH1F *g11 = new TH1F("g11","xlocal layer 12",50,-1.0, 1.0);
TH1F *g12 = new TH1F("g12","ylocal layer 12",50,-0.025, 0.025);
TH1F *g13 = new TH1F("g13","zlocal layer 12",50,-5.0, 5.0);
TH2F *g14 = new TH2F("g14","z-x local layer 12",50,-5.0, 5.0, 50, -2.0, 2.0);
TH2F *g15 = new TH2F("g15","x-y global layer 12 unique", 100, -5.0, 5.0, 100, -18.0, 0.0);
TH2F *g15a = new TH2F("g15a","x-y global layer 12 unique", 40, -20.0, 20.0, 40, -20.0, 20.0);
TH2F *g16 = new TH2F("g16","r-z global layer 12", 30,  25.0, 40.0, 20, 0.0, 20.0);

TH1F *g21 = new TH1F("g21","xlocal layer 12",50,-1.0, 1.0);
TH1F *g22 = new TH1F("g22","ylocal layer 12",50,-0.025, 0.025);
TH1F *g23 = new TH1F("g23","zlocal layer 12",50,-5.0, 5.0);
TH2F *g24 = new TH2F("g24","z-x local layer 12",50,-5.0, 5.0, 50, -2.0, 2.0);
TH2F *g25  = new TH2F("g25" ,"x-y global layer 12", 40,-20.0, 20.0, 40, -20.0, 20.0);
TH2F *g25d  = new TH2F("g25d" ,"x-y global layer 12 folded", 40,-20.0, 20.0, 40, -20.0, 20.0);
TH2F *g25a = new TH2F("g25a","x-y global layer 12", 60, -2.5, 2.5, 48, -18.0, -2.0);
TH2F *g25b = new TH2F("g25b","phi-r global layer 12", 100, -0.04, 0.16, 48, 2.0, 18.0);
TH2F *g25c = new TH2F("g25c","phi-r global layer 12 for overlaps", 100, -0.40, 0.40, 48, 2.0, 18.0);
TH2F *g25e = new TH2F("g25e","x-y global layer 12 for overlaps", 100, -5.0, 5.0, 100, -18.0, 0.0);
TH2F *g26 = new TH2F("g26","r-z global layer 12", 30,  25.0, 40.0, 20, -20.0, 0.0);
TH1F *g27 = new TH1F("g27","all delz layer 12", 50,-1.0, 1.0);
TH1F *g28 = new TH1F("g28","delz layer 12 |X|<.3", 50,-1.0, 1.0);

TH1F *m01 = new TH1F("m01","all tracks, north layers hit", 10, -0.5, 9.5);
TH2F *m02 = new TH2F("m02","r vs z first hit", 40, 10.0, 30.0, 20, 0.0, 20.0);
TH2F *m03 = new TH2F("m03","y vs x first hit", 40, -20.0, 20.0, 40, -20.0, 20.0);
TH2F *m04 = new TH2F("m04","r vs z last hit", 40, 20.0, 40.0, 20, 0.0, 20.0);
TH2F *m05 = new TH2F("m05","y vs x last hit", 40, -20.0, 20.0, 40, -20.0, 20.0);
TH2F *m06 = new TH2F("m06","y vs x at z=0", 40, -50.0, 50.0, 40, -50.0, 50.0);
TH1F *m07 = new TH1F("m07","r in plane of station 1", 40, 0., 200.0);
TH2F *m08 = new TH2F("m08","r vs z first and last hits", 40, 0.0, 40.0, 20, 0.0, 20.0);
TH1F *m09 = new TH1F("m09","x last hit for y<1cm",40, -20.0, 20.0);
TH1F *m10 = new TH1F("m10","delz vertex, perfect tracking",40, -40.0, 40.0);
TH1F *m11 = new TH1F("m11","delz vertex, perfect tracking",40, -0.1, 0.1);
TH1F *m12 = new TH1F("m12","delz vertex, integerized",40, -40.0, 40.0);
TH1F *m13 = new TH1F("m13","delz vertex, integerized",40, -0.1, 0.1);
TH1F *m14 = new TH1F("m14","event z vertex",40, -40.0, 40.0);

TH1F *m15 = new TH1F("m15","zlocin",40, -5.0, 5.0);
TH1F *m16 = new TH1F("m16","zlocout - zlocin",40, 0.0, 0.04);
TH1F *m17 = new TH1F("m16","zlocout - zlocin",40, 0.0, 0.04);



//----------------------xxx------------------------------------------
icount = 0;
//nentries=1000; 
  for (Long64_t i=0; i<nentries;i++) {            // master loop over records
    nbytes += AncSvx->GetEntry(i);                // pull in the record
    //    cout<<" event, track, layer= "<<EVENT<<" "<<TRACK<<" "<<LAYER<<endl;
            rglobal = sqrt(XGLOBAL**2 + YGLOBAL**2);
      cout<<"EVENT: "<<EVENT<<", track: "<<TRACK<<", layer: "<<LAYER
      	  <<", z: "<<ZGLOBAL<<", Rglobal: "<<rglobal<<endl;
    if (EVENT != oldevent) {                    // new event, starting at events 2!
      oldevent = EVENT;
      cout<<"New event seen: "<<EVENT<<endl;
      cout<<"#of records for event "<<EVENT-1<<" was "<<mult0<<endl;
      e01->Fill( float(mult0) );
      cout<<"track seen in prev event: "<<track_max<<endl;
                                                // loop over tracks in the event
      for (itrack=1; itrack<=track_max; itrack++ ) {
        //cout<<"event"<<EVENT<<" itrack "<<itrack<<" pphi:"<<PPHI<<endl;
        nlhit = 0;
        for (ilayer=9; ilayer<=12; ilayer++) {     // count layers hit
          if (hit[itrack][ilayer][1]!=0.) { 
          nlhit++;
            multistrip = hit[itrack][ilayer][5] - hit[itrack][ilayer][4];
            m15->Fill(hit[itrack][ilayer][5]);        // zlocin
            m16->Fill(fabs(multistrip));                    // zlocout - zlocin
	  }
        } 
        //          cout <<" itrack, nlhit= "<<itrack<<" "<<nlhit<<endl;
        m01->Fill(float(nlhit));
        if (nlhit>=3) {                            // 3 or more layers hit
          if (hit[itrack][9][1]==0){               // 10, 11, 12
            for (ixyz=1; ixyz<=3; ixyz++) {
              xa[ixyz] = hit[itrack][10][ixyz];
              xb[ixyz] = hit[itrack][12][ixyz];
            }
          }
          else if (hit[itrack][12][1]==0){            // 9, 10, 11
            for (ixyz=1; ixyz<=3; ixyz++) {
              xa[ixyz] = hit[itrack][9][ixyz];
              xb[ixyz] = hit[itrack][11][ixyz];
            }
          }
          else {                                       // 9, 10, 11, 12
            for (ixyz=1; ixyz<=3; ixyz++) {
              xa[ixyz] = hit[itrack][9][ixyz];
              xb[ixyz] = hit[itrack][12][ixyz];}
          }             // xa, xb is first, last hit
          ra = sqrt(xa[1]**2+xa[2]**2);        /// plot hits in first and last layers
          phia = atan(xa[2]/xa[1]);
          if (xa[1]<0) {phia = phia+3.14159;}
          if (xa[1]>0 && xa[2]<0) {phia = phia+2.0*3.14159;}
          m02->Fill(xa[3],ra);
          m08->Fill(xa[3],ra);
          m03->Fill(xa[1],xa[2]);
          e02->Fill(phia);
          rb = sqrt(xb[1]**2+xb[2]**2);
          phib = atan(xb[2]/xb[1]);
          if (xb[1]<0) {phib = phib+3.14159;}
          if (xb[1]>0 && xb[2]<0) {phib = phib + 2.0*3.14159;}
          m04->Fill(xb[3],rb);
          m08->Fill(xb[3],rb);
          m05->Fill(xb[1],xb[2]);
          e03->Fill(phib);
          e04->Fill(phia-phib);
          e04b->Fill(phia-phib);
          if ( phia == phib ){cout<<"phia,b identical"<<endl;}
          e05->Fill(phia,phib);
          if (abs(xb[2])<1.0) {m09->Fill(xb[1]);}    // slice in y


          zs1 = 0.0;                 // plane at z=0
          m14->Fill(zvertex);
          if ( (xb[3]-xa[3]) != 0.0) {
            xs1 = (zs1-xa[3])*(xb[1]-xa[1])/(xb[3]-xa[3]) +xa[1];
            ys1 = (zs1-xa[3])*(xb[2]-xa[2])/(xb[3]-xa[3]) +xa[2];
            rs1 = sqrt(xs1**2+ys1**2);
            m06->Fill(xs1,ys1);            //x-y at station 1
            m07->Fill(rs1);
            zv = (ra*xb[3] - rb*xa[3]) / (ra-rb) - zvertex;   // z vertex
            m10->Fill(zv);
            if (fabs(zv)<0.1) {m11->Fill(zv);}

            // integerize radii, and adjust z accordingly

            irstrip = (ra/rseg + 0.5);
            rsave = ra;
            ra = rseg*irstrip;
            xa[3] = xa[3] + (ra-rsave)*sin_alpha;
            irstrip = (rb/rseg + 0.5);
            rsave = rb;
            rb = rseg*irstrip;
            xb[3] = xb[3] + (rb-rsave)*sin_alpha;
            zv = (ra*xb[3] - rb*xa[3]) / (ra-rb) - zvertex;   // z vertex
            m12->Fill(zv);
            if (fabs(zv)<0.1) {m13->Fill(zv);}


            if (icount<100 && rs1<10.0) {
              icount++;
              cout<<icount<<": track  "<<itrack<<endl;
              cout<<"xa[1,2,3]= "<<xa[1]<<" "<<xa[2]<<" "<<xa[3]<<endl;
              cout<<"xb[1,2,3]= "<<xb[1]<<" "<<xb[2]<<" "<<xb[3]<<endl;
              cout<<"xs1,2,3= "<<xs1<<" "<<ys1<<" "<<zs1<<endl;
            }             // show a few events
          }
        }               // >= 3 planes hit
      }                 // loop over all tracks
                
      mult0 = 1;                   // reset hit multiplicity counter
      //      hit[][][]
    }
    else {                                 // we're in the same event
      mult0++;                             // hit count
      hit[TRACK][LAYER][1] = XGLOBAL;
      hit[TRACK][LAYER][2] = YGLOBAL;
      hit[TRACK][LAYER][3] = ZGLOBAL;
      hit[TRACK][LAYER][4] = ZLOCIN;
      hit[TRACK][LAYER][5] = ZLOCOUT;
      zvertex   = Z_VERTEX;
      track_max = TRACK;
      //      rglobal = sqrt(XGLOBAL**2 + YGLOBAL**2);
      //cout<<"EVENT: "<<EVENT<<", track: "<<TRACK<<", layer: "<<LAYER
      //	  <<", z: "<<ZGLOBAL<<", rglobal: "<<rglobal<<endl;
    }       
  }               // loop over records

//------------------------------------------------------------------------------------
}            // end of macro