{
//////////////////////////////////////////////////////////
//   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.00500;  // flat: 50 microns
  Int_t irstrip = 0, jrstrip=0, nstrips, irold;
  Float_t sin_alpha = 0.447;
  Float_t multistrip, za, zb, z[10000], zcluster, zreal, delz_ncl1, delz_ncl2;
  Float_t Rseg, zrold;

TH1F *e01 = new TH1F("e01","records per event",50, 0.0, 1000.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 *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("m17","delta-r for cluster size 1",    50, -60., 140.);
TH1F *m18 = new TH1F("m18","delta-r for cluster size 2",    50, -60., 140.);
TH1F *m19 = new TH1F("m19","delta-r for cluster size 3",    50, -60., 140.);
TH1F *m20 = new TH1F("m20","delta-r for cluster size 4",    50, -60., 140.);
TH1F *m21 = new TH1F("m21","delta-r for cluster size 5",    50, -60., 140.);
TH1F *m29 = new TH1F("m29","delta-r for old integerization",50, -60., 140.); 

//----------------------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 <<"multistrip: "<< multistrip<<endl;
//----clusters
            Rseg = 10000*rseg;                           // segmentation, microns
	    za = 10000*(hit[itrack][ilayer][4] + 10.0);  //entry, microns
	    zb = 10000*(hit[itrack][ilayer][5] + 10.0);  //exit, microns
            zreal = (za+zb)/2.0;            // real coordinate
            irold = zreal/Rseg;
            zrold = (irold + 0.5) * Rseg;   // integerized
            m29->Fill(zreal - zrold);       // plot integerized error

            irstrip = za/Rseg;
            z[1] = Rseg*irstrip;
            jrstrip = zb/Rseg + 1;
	    nstrips =  jrstrip - irstrip;
            //cout <<"za, zb, irstrip, jrstrip, nstrips "<<za<<" "<<zb<<" "
	    //   <<irstrip<<" "<<jrstrip<<" "<<nstrips<<endl;

	    if (nstrips==1) {
              z[2] = z[1] + Rseg;
              //zcluster = z[1] +0.5*Rseg;
              //delz_ncl1 = zcluster - zreal;
              //m17->Fill(delz_ncl1);
	    }
	    else {
	      z[nstrips] = z[1] + (nstrips-1)*Rseg;
	      z[nstrips+1] = z[nstrips] + Rseg;
	      for (int ii=2; ii<nstrips; ii++) {
                z[ii] = z[1] + (ii-1)*Rseg;
	      }
              delz_ncl2 = (z[nstrips+1]-z[1])/2.0 - zreal;
              //m18->Fill(delz_ncl2);
	    }
            zcluster = (z[1] + z[nstrips+1])/2.0;
            if (nstrips==1) {m17->Fill(zcluster-zreal);}
            if (nstrips==2) {m18->Fill(zcluster-zreal);}
            if (nstrips==3) {m19->Fill(zcluster-zreal);}
            if (nstrips==4) {m20->Fill(zcluster-zreal);}
            if (nstrips==5) {m21->Fill(zcluster-zreal);}

	    //cout <<"nstrips,entry,exit= "<<nstrips<<" "<< za<<" "<<zb<<". z[i]: ";
	    for (int ii=1; ii<=nstrips+1; ii++) {
              //cout <<z[ii]<<"  ";
	    }
            //cout <<endl;
//---clusters

          }
        } 
        //          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);
            ra = rseg*irstrip;
            irstrip = (rb/rseg + 0.5);
            rb = rseg*irstrip;
            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;
      //cout <<"zlocout: "<<hit[TRACK][LAYER][5]<<endl;
      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