{
//////////////////////////////////////////////////////////
//   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");
//
//////////////////////////////////////////////////////////
//  #include <iostream>;
//  using namespace std;

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

   Int_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[40], rfirst, rlast, rtmp, xs1, ys1, zs1, rs1, track_max;
  Float_t xfirst[6], xlast[6];      // 0=layer#, 1,2,3 = x,y,zglobal, 5,6=zlocal
  Float_t scale, zv, rglobal, zvertex, phifirst, philast, phi, xglobal, path_org;
  Float_t delr_shade, tanshade, cosshade, path_total, sstot, mytot, alf, bet;
  Int_t itrack, lnhit, ilayer, nlhit, ixyz, icount;
  Int_t sili_endcap_type = 0;           // type to be determined by data
  Float_t sumx[13]  = {0.0};
  Float_t sumxx[13] = {0.0};
  Float_t nx[13]    = {0.0};
  Float_t ss[13]    = {0.0};
  Float_t mm[13]    = {0.0};
  Float_t my[13]    = {0.0};
  Float_t sumy[13]  = {0.0};
//                             // for integerization:
    Float_t rseg = 0.0050;  // flat: 50 microns
//  Float_t rseg = 0.0075;  // flat: 75 microns
//  Float_t rseg = 0.0080;  // flat: 80 microns
//  Float_t rseg = 0.0120;  // flat: 120 microns
//  Float_t rseg = 0.0200;  // flat: 200 microns

//--- --- threshold cuts:
    Float_t zcut = 0.0000;
//  Float_t zcut = 0.0025;
//  Float_t zcut = 0.0034;
//  Float_t zcut = 0.0050;
//  Float_t zcut = 0.0075;
//  Float_t zcut = 0.0125; 
//
  Int_t irstrip = 0, jrstrip=0, irold, int_ang, int_r, int_y,
        iystrip, jystrip, nystrips, ixstrip, jxstrip, nxstrips, nele, ntmp,
        int_lr;
  Int_t  nstrip[11] = {0};
  Float_t sin_alpha = 0.447;
  Float_t multistrip, z[100], rfcluster, rlcluster, zreal, delz_ncl1, zpath[100];
  Float_t zrold, sinalpha, ztmp, delr[13], delrfirst, delrlast, xin, xout, xreal;
  Float_t vpanel[13][4], realr1[13], realr2[13];
  Float_t vlength, vpfirst[3], vplast[3], rfreal, rlreal, rfin, rfout, rlin, rlout,
          ang, r2, ytmp, zthick, zlow, y[100], yin, yout, xfin, xfout, x[100],
          ryseg, strips[11];
  Float_t  estrip[11] = {0.0};

  TH1F *e00 = new TH1F("e00","tracks per event",50, 0.0, 1000.0);
  TH1F *d01 = new TH1F("d01","z layer 12",50, 30.0, 40.0);
  TH2F *d02 = new TH2F("d02","midpoint layer 9 ",120, -15., 15., 120, -15., 15.);
  TH2F *d03 = new TH2F("d03","midpoint layer 10",120, -15., 15., 120, -15., 15.);
  TH2F *d04 = new TH2F("d04","midpoint layer 11",120, -15., 15., 120, -15., 15.);
  TH2F *d05 = new TH2F("d05","midpoint layer 12",120, -15., 15., 120, -15., 15.);

  TH1F *d06 = new TH1F("d06","inner midpoint angle",240, -3.927, 3.927);
  TH1F *d07 = new TH1F("d07","outer midpoint angle",240, -3.927, 3.927);
  TH1F *d08 = new TH1F("d08","inner midpoint int",240, -5, 55);
  TH1F *d09 = new TH1F("d09","outer midpoint int",240, -5, 55);
  
  TH1F *d10 = new TH1F("d10","midpoint layer 9 ", 250, 3.1, 17.1);
  TH1F *d11 = new TH1F("d11","midpoint layer 10", 50, 3.1, 17.1);
  TH1F *d12 = new TH1F("d12","midpoint layer 11", 50, 3.1, 17.1);
  TH1F *d13 = new TH1F("d13","midpoint layer 12", 50, 3.1, 17.1);
  TH1F *d14 = new TH1F("d14","int layer xx", 100, 0., 25.);
  
  TH2F *d15 = new TH2F("d15"," z local - r=1", 100, -5., 5., 100, -5., 5.);
  TH2F *d16 = new TH2F("d16"," z local - r=2", 100, -5., 5., 100, -5., 5.);
  TH2F *d17 = new TH2F("d17"," z local - r=2", 100, -5., 5., 100, -5., 5.);
  TH2F *d18 = new TH2F("d18"," radius index, left-right index",5,0.25,2.75,5,0.25,2.75);

  TH1F *m23 = new TH1F("m23","electrons for ~ normal tracks",100, 0., 66000); 
  TH1F *m24 = new TH1F("m24","DEL-E for ~ normal tracks",100, 0., 0.3); 
  TH1F *m25 = new TH1F("m25","iystrip zlow= 3.3",50, 0., 2000.); 
  TH1F *m26 = new TH1F("m26","iystrip zlow= 2.0",50, 0., 2000.); 
  TH1F *m27 = new TH1F("m27","iystrip zlow =3.95",50, 0., 2000.); 
  TH1F *m36 = new TH1F("m36","total path, original ",50, 0., 350.); 
  TH2F *m28 = new TH2F("m28","x,y pixel numbers first plane",50, -100., 100., 50, 650, 1650);
  TH1F *m29 = new TH1F("m29","nystrips ",40, -1., 19.); 
  TH1F *m31 = new TH1F("m31","path length (um) for n=1",50, 299., 309.); 
  TH2F *m33 = new TH2F("m33","middle path lengths (um) for n>=3",50, 0., 350., 10, 0, 10); 
  TH2F *m34 = new TH2F("m34","first path length (um) for n>=2",50, 0., 350., 10, 0, 10); 
  TH2F *m35 = new TH2F("m35","last path length (um) for n>=2",50, 0., 350., 10, 0, 10); 
  TH2F *m37  = new TH2F("m37","nystrips vs. total path, summed parts",50,0,350.,5,0.5,5.5);

  TFile *ntuple_file = new TFile("from_anc.root","RECREATE");
  TNtuple *digi = new TNtuple("digi","digitization, from ancsvx",
  "dummy:event:track:layer:xglobal:yglobal:zglobal:zlocin:zlocout:xlocin:xlocout:dele:iang:ir:ilr:nstrips:istrip1:istrip2:istrip3:istrip4:istrip5:istrip6:istrip7:istrip8:istrip9:istrip10:mstrips:estrip1:estrip2:estrip3:estrip4:estrip5:estrip6:estrip7:estrip8:estrip9:estrip10:a37:a38:a39");            // must be one string

//----------------------xxx------------------------------------------
  icount = 0; 
//  nentries=3000;     //xxx   500 records is enough for 1st event

  cout <<" ***** nentries **** "<<nentries<<endl;
  for (int i=0; i<nentries;i++) {            // master loop over records
    nbytes += AncSvx->GetEntry(i);                // pull in the record
    if (EVENT != oldevent) {                    // new event, starting at events 2!
      cout<<"2) EVENT, oldeenvt= "<<EVENT<<" "<<oldevent<<endl;
      oldevent = EVENT;
      cout<<"New event seen: "<<EVENT<<endl;
      cout<<"#of records for event "<<EVENT-1<<" was "<<mult0<<endl;
      //-------------------------------- determine detector type ----------------------------
      if (sili_endcap_type==0) {
        cout <<"   event "<<EVENT<<": determining detector type"<<endl;
        for (ilayer=9; ilayer<=12; ilayer++) {
          ss[ilayer] = sqrt(nx[ilayer]*sumxx[ilayer] - sumx[ilayer]*sumx[ilayer])/nx[ilayer];
          mm[ilayer] = sumx[ilayer]/nx[ilayer];
          my[ilayer] = sumy[ilayer]/nx[ilayer];
        }
        sstot = (ss[9] + ss[10] + ss[11] + ss[12])/4.0;
        mytot = (my[9] + my[10] + my[11] + my[12])/4.0;
        cout<<"   sigma ZGLOBAL for layers 9-12: "<<sstot<<endl;
        cout<<"   mean  thickness   layers 9-12: "<<mytot<<endl;
        if (sstot>0.3 && mytot>0.0275) sili_endcap_type = 1;
        if (sstot<0.3 && mytot>0.0275) sili_endcap_type = 2;
        if (sstot>0.3 && mytot<0.0275) sili_endcap_type = 3;
        if (sstot<0.3 && mytot<0.0275) sili_endcap_type = 4;
        cout<<"   ==> analyzing for type "<<sili_endcap_type<<endl;
                                   zthick = 0.0300;
	if (sili_endcap_type>=3) { zthick = 0.0250; }
        ryseg = 0.0050;        // 50 micron strip width default
        //ryseg = 0.0075;
	//        e01->Fill( float(mult0) );
        //e00->Fill( track_max );
      }
      if (sili_endcap_type>=1 && sili_endcap_type<=4) {
      }   // detector type recognized
    }
//-------------------------------------------------------------------------------------------
//                           PREPARE THE NEXT EVENT
//-------------------------------------------------------------------------------------------
    else if ( LAYER>=5 && LAYER<=12 ) {                      // we're in the same event
      mult0++;                             // hit count
      //------------------------determine detector type (use 1st event): --------------------
      if (sili_endcap_type==0){
        int il = 0 + LAYER ;        // convert to integer
        sumx [il] = sumx [il] + ZGLOBAL;
        sumxx[il] = sumxx[il] + ZGLOBAL*ZGLOBAL;
        nx   [il] = nx   [il] + 1.0;
        sumy [il] = sumy [il] + YLOCOUT-YLOCIN;
           
        if (il==12) {d01->Fill(ZGLOBAL);}
      }
      else {             //------------- event 2 and up --------------------------------

	//---- digitize: calculate center of local silicon wafer ----------------------------
        if ( sili_endcap_type ==1 || sili_endcap_type ==2 ) {      
          alf = atan2(XGLOBAL,YGLOBAL);
          //cout<<"x y global, alpha: "<<XGLOBAL<<" "<< YGLOBAL<<" "<< alf<<endl;
          bet = asin(XLOCIN/sqrt(XGLOBAL*XGLOBAL+YGLOBAL*YGLOBAL));
          //cout<<"xloc, rglob, beta: "<<XLOCIN<<" "
          //          <<sqrt(XGLOBAL*XGLOBAL+YGLOBAL*YGLOBAL)<<" "<<bet<<endl;
          r2 = sqrt(XGLOBAL*XGLOBAL+YGLOBAL*YGLOBAL+XLOCIN*XLOCIN);
          //  cout<<"r2: "<<r2<<endl;
          x2 = (r2-ZLOCIN) * sin(alf+bet);
          y2 = (r2-ZLOCIN) * cos(alf+bet);
          if (LAYER== 9) {d02->Fill(x2,y2);}
          if (LAYER==10) {d03->Fill(x2,y2);}
          if (LAYER==11) {d04->Fill(x2,y2);}
          if (LAYER==12) {d05->Fill(x2,y2);}
          //---integerize these positions: ---------------------------------------------
          if ((r2-ZLOCIN)<10.0) {d06->Fill(alf+bet);}
          if ((r2-ZLOCIN)>10.0) {d07->Fill(alf+bet);}
          ang = (alf+bet)/(7.5*3.14159/180.0) +0.5 + 24.0;
          int_ang = TMath::Nint(ang);
          if (int_ang>48) {int_ang=int_ang-48;}
          if (int_ang< 1) {int_ang=int_ang+48;}          //   <============= angle index
          //cout <<"ang, int_ang: "<<ang<<" "<<int_ang<<endl;
          if ((r2-ZLOCIN)<10.0) {
            d08->Fill(float(int_ang));
            int_r = 1;
          }
          if ((r2-ZLOCIN)>10.0) {
            d09->Fill(float(int_ang));
            int_r = 2;                                   //   <============ radius index
          }
          int_lr = 1;                                    //   <======== left-right index
          if (XLOCIN<0.) {int_lr = 2;}
          d18->Fill(float(int_lr), float(int_r));
        }                                                // DOEendcaps, type 1 and 2   

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

        else if ( sili_endcap_type ==3 || sili_endcap_type ==4 ) {
                                    delz = 0.60;         // straight
          if (sili_endcap_type==3) {delz = 0.56;}        // tilted
	  x2 = XGLOBAL - XLOCIN;
          y2 = YGLOBAL + ZLOCIN;
          //cout<<"y2: "<<y2<<endl;
          if (LAYER== 9) {d10->Fill(y2);}
          if (LAYER==10) {d11->Fill(y2);}
          if (LAYER==11) {d12->Fill(y2);}
          if (LAYER==12) {d13->Fill(y2);}
          //---integerize these positions: ---------------------------------------------
          ytmp =  (y2-3.82)/delz + 1.0;      // 3.5+0.32 see svx.f line 1630
	  cout<<"ytmp: "<<ytmp<<endl;
          int_y = TMath::Nint(ytmp);       // int_y is the 8-chip module #
          d14->Fill(float(int_y));
        }             // LDRD endcaps, type 3 and 4

        //--------- now find the strip number, using ZLOCIN ----------------------------
        if ( sili_endcap_type ==1 || sili_endcap_type ==2 ) {
          if (int_r==1) { 
            zlow = 3.3;              // = silen /2 in svx.f 
            d15->Fill(XLOCIN,ZLOCIN);
	  }
          else if (int_r==2 && LAYER==10) { 
            if ( sili_endcap_type==1 ) {
              zlow = 2.0;            // = silen7 /2 in svx.f
            }
            else if ( sili_endcap_type==2 ) {
              zlow = 3.95;           // same as layers 11, 12
	    }
            d16->Fill(XLOCIN,ZLOCIN);
          }
          else if (int_r==2 && LAYER>=11) { 
            zlow = 3.95;             // = silen4 /2 
            d17->Fill(XLOCIN,ZLOCIN);
          }

          for (int ii=0; ii<=99; ii++) { y[ii]=0.; zpath[ii]=0.; }  //clear
          //         [4]=YLOCIN  [5]=YLOCOUT  (Z is locally 'radial')

          yin  = ZLOCIN  + zlow + ryseg;    // strip #1 closest to beam
          yout = ZLOCOUT + zlow + ryseg;    // start counting at 1
          ytmp = -100.;                     //
          if (yin>yout) {                   // order yin, yout:
            ytmp = yin; yin = yout; yout = ytmp;
	  }

          path_org = sqrt((yout-yin)**2 + zthick**2);
          m36->Fill(1e4*path_org);          // should be >~=300

          iystrip = yin/ryseg;              // integerize the entry point
          if(iystrip<1) {cout<<"iystrip= "<<iystrip<<endl;}
          if (zlow>3.0&&zlow<3.5 ) { m25->Fill(float(iystrip)); }     //=3.3  type 1,2
          if (zlow<2.1 )           { m26->Fill(float(iystrip)); }     //=2.0  type 1 
          if (zlow>3.94)           { m27->Fill(float(iystrip)); }     //=3.95 type 1,2
          y[1] = ryseg*iystrip;             // cell boundary below the entry point
          jystrip = yout/ryseg + 1;         // integerize the exit point
	  nystrips =  jystrip - iystrip;    //
          m29->Fill(float(nystrips));

          path_total = 0.;
	  if (nystrips==1) {
            y[2] = y[1] + ryseg;            // cell boundary above the exit point
            zpath[1] = sqrt((yout-yin)**2+zthick**2);   // single-cell path length
            m31->Fill(1e4*zpath[1]);
            path_total = zpath[1];
            m24->Fill(DELE);
            nele = path_total * DELE * 733e4;
            m23->Fill(float(nele));
	  }
	  else {                                  // 2 or more strips
	    y[nystrips] = y[1] + (nystrips-1)*ryseg;
	    y[nystrips+1] = y[nystrips] + ryseg;
	    for (int ii=2; ii<nystrips; ii++) {
             y[ii] = y[1] + (ii-1)*ryseg;
	    }
            sinalpha = (yout-yin)/sqrt((yout-yin)**2+ zthick**2);
            zpath[1]       = (y[2]-yin      )/sinalpha;
            zpath[nystrips] = (yout-y[nystrips])/sinalpha;
            path_total = path_total + zpath[1] + zpath[nystrips];
            for (int ii=2; ii<nystrips; ii++) {
              zpath[ii] = ryseg/sinalpha;
              path_total = path_total + zpath[ii];
            }
	  }                                      // 2 or more strips hit
	  if (nystrips>=2) {
            for (int ii=2; ii<=nystrips-1; ii++) {
              m33->Fill(1e4*zpath[ii], float(nystrips));
            }
            m34->Fill(1e4*zpath[1]      , float(nystrips));
            m35->Fill(1e4*zpath[nystrips], float(nystrips));
          }

          m37->Fill(1e4*path_total,nystrips);
          if (nystrips<1) {
            cout<<"ERROR 1: nystrips<1 "<<nystrips<<endl;
	  }

          ntmp = TMath::Min(nystrips,10);         // store strips hit 
          estrip[0] = float(ntmp);                // store 10 max
          nstrip[0] = ntmp;                       // store 10 max
          for (int ii=1; ii<=ntmp; ii++) {
            nstrip[ii] = ii + iystrip-1;         // iystrip is 1st strip # <===========
            estrip[ii] = zpath[ii] * DELE * 733e4;  //           <=====================
          }

//--- done! now store all: ----------------------

          hit[1] = EVENT;
          hit[ 2] = TRACK; 
          hit[ 3] = LAYER;        // layer
 
          hit[ 4] = XGLOBAL;
          hit[ 5] = YGLOBAL;
          hit[ 6] = ZGLOBAL;
          hit[ 7] = ZLOCIN;
          hit[ 8] = ZLOCOUT;
          hit[ 9] = XLOCIN;
          hit[10] = XLOCOUT;
	  hit[11] = DELE;               // total energy loss
          hit[12] = int_ang;            // angle index of si panel
	  hit[13] = int_r;              // radius index of si panel
	  hit[14] = int_lr;             // radius index of si panel
	  for (ii=0; ii<=10; ii++) {    // loop over hit cluster
            hit[15+ii] = nstrip[ii];    // silicon strip number
            hit[26+ii] = estrip[ii];    // # electrons in the signal
	  }

        }                               // strip # type 1, 2
        if (LAYER>=9) { digi->Fill(hit); }               // fill the ntuple

      }           // event 1 / 2 and up
    }             // same event / new event
  }               // loop over records

ntuple_file->Write();                  // save the ntuple to disk
//------------------------------------------------------------------------------------
}            // end of macro