{
//////////////////////////////////////////////////////////
//   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[8000][13][8], 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;
  Int_t itrack, lnhit, ilayer, nlhit, ixyz, icount;
  int sili_endcap_type = 2;           // 1=shades, 2=flat

//                             // 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, nstrips, irold;
  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;

  TH1F *e00 = new TH1F("e00","tracks per event",50, 0.0, 1000.0);
  TH1F *e01 = new TH1F("e01","records per event",50, 0.0, 1000.0);
  TH1F *e02 = new TH1F("e02","phi hit a",50, -1.0, 7.0);
  TH1F *e03 = new TH1F("e03","phi hit b",50, -1.0, 7.0);
  TH1F *e04 = new TH1F("e04","phifirst-philast",  70, -7.0, 7.0);
  TH1F *e04b = new TH1F("e04b","phifirst-philast",  70, -0.01, 0.01);
  TH2F *e05 = new TH2F("e05","phifirst vs. philast",  50, -1.0, 7.0, 50, -1.0, 7.0);

  TH1F *e101 = new TH1F("e101","re: records per event",50, 0.0, 1000.0);
  TH1F *e102 = new TH1F("e102","re: phi hit a",50, -1.0, 7.0);
  TH1F *e103 = new TH1F("e103","re: phi hit b",50, -1.0, 7.0);
  TH1F *e104 = new TH1F("e104","re: phifirst-philast",  70, -7.0, 7.0);
  TH1F *e104b = new TH1F("e104b","re: phifirst-philast",  70, -0.01, 0.01);
  TH2F *e105 = new TH2F("e105","re: phifirst vs. philast",  50, -1.0, 7.0, 50, -1.0, 7.0);

  TH1F *m01 = new TH1F("m01","all tracks, # of 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",50, -50.0, 50.0);
  TH1F *m11 = new TH1F("m11","delz vertex, perfect tracking",50, -0.1, 0.3);
  TH1F *m12 = new TH1F("m12" ,"delz vertex, integerized",50, -50.0, 50.0);
  TH1F *m12a= new TH1F("m12a","delz vertex, integerized",50, -50.0, 50.0);
  TH1F *m13 = new TH1F("m13" ,"delz vertex, integerized",50, -0.1, 0.3);
  TH1F *m13a= new TH1F("m13a","delz vertex, clustered",50, -0.1, 0.3);
  TH1F *m14 = new TH1F("m14","event z vertex",50, -20.0, 20.0);
  TH1F *m22 = new TH1F("m22","delta-r 1st layer local, due to integerization",50, -100, 100.);
  TH1F *m23 = new TH1F("m23","delta-r last layer local, due to integerization",50, -100, 100);

  TH2F *m102 = new TH2F("m102","re: r vs z first hit", 40, 10.0, 30.0, 20, 0.0, 20.0);
  TH2F *m103 = new TH2F("m103","re: y vs x first hit", 40, -20.0, 20.0, 40, -20.0, 20.0);
  TH2F *m104 = new TH2F("m104","re: r vs z last hit", 40, 20.0, 40.0, 20, 0.0, 20.0);
  TH2F *m105 = new TH2F("m105","re: y vs x last hit", 40, -20.0, 20.0, 40, -20.0, 20.0);
  TH2F *m106 = new TH2F("m106","re: y vs x at z=0", 40, -50.0, 50.0, 40, -50.0, 50.0);
  TH1F *m107 = new TH1F("m107","re: r in plane of station 1", 40, 0., 200.0);
  TH2F *m108 = new TH2F("m108","re: r vs z first and last hits", 40, 0.0, 40.0, 20, 0.0, 20.0);
  TH1F *m109 = new TH1F("m109","re: x last hit for y<1cm",40, -20.0, 20.0);

  TH1F *m110 = new TH1F("m110","re: delz vertex, perfect tracking",50, -50.0, 50.0);
  TH1F *m111 = new TH1F("m111","re: delz vertex, perfect tracking",50, -0.1, 0.3);
  TH1F *m112 = new TH1F("m112","re: delz vertex, integerized",50, -50.0, 50.0);
  TH1F *m113 = new TH1F("m113","re: delz vertex, integerized",50, -0.1, 0.3);
  TH1F *m114 = new TH1F("m114","re: event z vertex",50, -20.0, 20.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 *m17a = new TH1F("m17a","delta-r for cluster size 1",    50, -60., 140.);
  TH1F *m18a = new TH1F("m18a","delta-r for cluster size 2",    50, -60., 140.);
  TH1F *m19a = new TH1F("m19a","delta-r for cluster size 3",    50, -60., 140.);
  TH1F *m20a = new TH1F("m20a","delta-r for cluster size 4",    50, -60., 140.);
  TH1F *m21a = new TH1F("m21a","delta-r for cluster size 5",    50, -60., 140.);

  TH1F *m29 = new TH1F("m29","delta-r for old integerization",50, -60., 140.); 

  TH1F *m30 = new TH1F("m30","raw number of strips traversed",20, 0.25, 10.25); 
  TH1F *m31 = new TH1F("m31","path length (um) for n=1",50, 295., 305.); 
  TH1F *m32 = new TH1F("m32","first layer path length all strips",50, 0., 350.); 
  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); 

  TH1F *m30a = new TH1F("m30a","number of strips traversed after threshold cut",20, 0.25, 10.25); 
  TH1F *m31a = new TH1F("m31a","path length (um) for n=1",50, 295., 305.); 
  TH1F *m32a = new TH1F("m32a","last layer path length all strips",50, 0., 350.); 
  TH2F *m33a = new TH2F("m33a","middle path lengths (um) for n>=3",50, 0., 350., 10, 0, 10); 
  TH2F *m34a = new TH2F("m34a","first path length (um) for n>=2",50, 0., 350., 10, 0, 10); 
  TH2F *m35a = new TH2F("m35a","last path length (um) for n>=2",50, 0., 350., 10, 0, 10); 
  TH1F *m36  = new TH1F("m36","total path, original ",50, 0., 350.); 
  TH2F *m37  = new TH2F("m37","nstrips vs. total path, summed parts ",50, 0., 350.,5, 0.5, 5.5); 

  TH1F *m40 = new TH1F("m40","raw number of strips traversed after cut",20, 0.25, 10.25); 
  TH1F *m41 = new TH1F("m41","path length (um) for n=1 after cut",50, 295., 305.); 
  TH1F *m42 = new TH1F("m42","path length (um) for n>1 after cut",50, 0., 350.); 
  TH2F *m43 = new TH2F("m43","middle path lengths (um) for n>=3 after cut",50, 0., 350., 10, 0, 10); 
  TH2F *m44 = new TH2F("m44","first path length (um) for n>=2 after cut",50, 0., 350., 10, 0, 10); 
  TH2F *m45 = new TH2F("m45","last path length (um) for n>=2 after cut",50, 0., 350., 10, 0, 10); 

  TH1F *m47 = new TH1F("m47","delta-r for cluster size 1 (after thresh)",    50, -60., 140.);
  TH1F *m48 = new TH1F("m48","delta-r for cluster size 2 (after thresh)",    50, -60., 140.);
  TH1F *m49 = new TH1F("m49","delta-r for cluster size 3 (after thresh)",    50, -60., 140.);
  TH1F *m50 = new TH1F("m50","delta-r for cluster size 4 (after thresh)",    50, -60., 140.);
  TH1F *m51 = new TH1F("m51","delta-r for cluster size 5 (after thresh)",    50, -60., 140.);

  TH1F *m99 = new TH1F("m99","xlocout",    50, -10., 10.);
  TH1F *m98 = new TH1F("m98","z12-z9",    50, 12., 22.);
  TH1F *mxxx = new TH1F("mxxx","out-in",    60, -20.0, 20.0);
  TH1F *myyy = new TH1F("myyy","out-in",    50, -0.0100, 0.0100);

//----------------------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
    rglobal = sqrt(XGLOBAL**2 + YGLOBAL**2);
    //if (LAYER==9){
    //  cout<<"1) EVENT, oldevent: "<<EVENT<<" "<<oldevent<<", track: "<<TRACK<<", layer: "<<LAYER
    //      <<", z: "<<ZGLOBAL<<", Rglobal: "<<rglobal<<endl;
    // }
    //cout<<"1a) EVENT, oldevent= "<<EVENT<<" "<<oldevent<<endl;
    if (EVENT != oldevent) {                    // new event, starting at events 2!
      cout<<"2) EVENT, oldevent= "<<EVENT<<" "<<oldevent<<endl;
      oldevent = EVENT;
      cout<<"New event seen: "<<EVENT<<endl;
      cout<<"#of records for event "<<EVENT-1<<" was "<<mult0<<endl;
      e01->Fill( float(mult0) );
      e00->Fill( track_max );

      //track_max=6;    //xxx for test
      
      for (itrack=1; itrack<=track_max; itrack++ ) {
        //cout<<"event"<<EVENT<<" itrack "<<itrack<<" pphi:"<<PPHI<<endl;
        m98->Fill(hit[itrack][12][3]-hit[itrack][9][3]);
        nlhit = 0;
        for (ilayer=9; ilayer<=12; ilayer++) {     // count layers hit
          delr[ilayer] = -999.;
          if (hit[itrack][ilayer][1]!=0.) { 
            xglobal = hit[itrack][ilayer][1];
            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;
          }          // layer is hit
        }            // loop over layers 9 - 12
 
        //          cout <<" itrack, nlhit= "<<itrack<<" "<<nlhit<<endl;
        m01->Fill(float(nlhit));

//--- start  tracking---------------------------------------------------------------------
        if (nlhit>=3) {                            // 3 or more layers hit
          xfirst[0] = 9.0;
          xlast [0] = 12.0;
                           
          if (hit[itrack][9][1]==0){               // 10, 11, 12                               ||
            xfirst[0] = 10.0;
            for (ixyz=1; ixyz<=5; ixyz++) {
              xfirst[ixyz] = hit[itrack][10][ixyz];
              xlast [ixyz] = hit[itrack][12][ixyz];
            }
          }
          else if (hit[itrack][12][1]==0){            // 9, 10, 11                             ||
            xlast[0] = 11.0;
            for (ixyz=1; ixyz<=5; ixyz++) {
              xfirst[ixyz] = hit[itrack][ 9][ixyz];
              xlast [ixyz] = hit[itrack][11][ixyz];
            }
          }
          else {                                       // 9, 10, 11, 12                        ||
            for (ixyz=1; ixyz<=5; ixyz++) {
              xfirst[ixyz] = hit[itrack][ 9][ixyz];
              xlast [ixyz] = hit[itrack][12][ixyz];
            }
          }             // xfirst, xlast is first, last hit coordinates -----------------------^^

          rfirst = sqrt(xfirst[1]**2+xfirst[2]**2);        // plot hits in first and last layers
          rfreal = rfirst;
          phifirst = atan2(xfirst[1],xfirst[2]);
          m02->Fill(xfirst[3],rfirst);
          m08->Fill(xfirst[3],rfirst);
          m03->Fill(xfirst[1],xfirst[2]);
          e02->Fill(phifirst);

          rlast = sqrt(xlast[1]**2+xlast[2]**2);
          rlreal = rlast;
          philast = atan2(xlast[1],xlast[2]);
          m04->Fill(xlast[3],rlast);
          m08->Fill(xlast[3],rlast);
          m05->Fill(xlast[1],xlast[2]);
          e03->Fill(philast);
 
          e04->Fill(phifirst-philast);
          e04b->Fill(phifirst-philast);
          e05->Fill(phifirst,philast);
          if (abs(xlast[2])<1.0) {m09->Fill(xlast[1]);}    // slice in y

          zs1 = 0.0;                 // plane at z=0
          m14->Fill(zvertex);
          if ( (xlast[3]-xfirst[3]) != 0.0) {
            xs1 = (zs1-xfirst[3])*(xlast[1]-xfirst[1])/(xlast[3]-xfirst[3]) +xfirst[1];
            ys1 = (zs1-xfirst[3])*(xlast[2]-xfirst[2])/(xlast[3]-xfirst[3]) +xfirst[2];
            rs1 = sqrt(xs1**2+ys1**2);
            m06->Fill(xs1,ys1);            //x-y at z=0
            m07->Fill(rs1);
	  }
          zv = -999.;
	    if ( (rfirst-rlast) != 0.0) {
            zv = (rfirst*xlast[3] - rlast*xfirst[3]) / (rfirst-rlast) - zvertex;
          }   // z vertex
          m10->Fill(zv);
          if (fabs(zv)<0.1) {m11->Fill(zv);}

//------- integerize radii, and adjust z accordingly-------------------------------

          irstrip = (rfirst/rseg + 0.5);
          rfirst = rseg*irstrip;
          m22->Fill(1e4*(rfirst-rfreal));
          irstrip = (rlast/rseg + 0.5);
          rlast = rseg*irstrip;
          m23->Fill(1e4*(rlast-rlreal));
          zv = -999.;
          if ( (rfirst-rlast) != 0.0) {
            zv = (rfirst*xlast[3] - rlast*xfirst[3]) / (rfirst-rlast) - zvertex;
          }   // z vertex
          m12->Fill(zv);
          if (fabs(zv)<0.1) {
            m13->Fill(zv);
          }

//------- clustering first layer hit: --------------

          rfirst = sqrt(xfirst[1]**2+xfirst[2]**2);        // restore rfirst, last

          for (int ii=0; ii<=99; ii++) { z[ii]=0.; zpath[ii]=0.; }  //clear

          rfin  = rfirst - 0.5*(xfirst[5] - xfirst[4]);    // there is some ambiguity
          rfout = rfirst + 0.5*(xfirst[5] - xfirst[4]);    // about the order of in/out
          if (rfin>rfout) {
            rtmp  = rfin;
            rfin  = rfout;
            rfout = rtmp;
	  }
          //rfin  = 1e4*r in;                                //
          //rfout = 1e4*rfout;                               // cm to microns
          //cout<<"xfirst 0,4,5 = "<<xfirst[0]<<" "<<xfirst[4]<<" "<<xfirst[5]<<" [5]-[4]: "
          //    <<xfirst[5]-xfirst[4]<<endl;

          //mxxx->Fill(rfin);  //xxx
          path_org = sqrt((rfout-rfin)**2 + 0.03**2);
          m36->Fill(1e4*path_org);

          irstrip = rfin/rseg;             // integerize the entry point
          z[1] = rseg*irstrip;             // cell boundary below the entry point
          jrstrip = rfout/rseg + 1;         // integerize the exit point
	  nstrips =  jrstrip - irstrip;    //
          if (nstrips<1 || nstrips>99) {
            cout<<"rfin, rfout, n ,i , j strips= "<<rfin<<" "<<rfout<<" "<<nstrips<<" "<<irstrip<<" "
            <<jrstrip<<endl;    
          }
          m30->Fill(float(nstrips));    // this one is also filled for last layer
          //cout <<"rfin, rfout, irstrip, jrstrip, nstrips "<<rfin<<" "<<rfout<<" "
	  //   <<irstrip<<" "<<jrstrip<<" "<<nstrips<<endl;

          path_total = 0.;
	  if (nstrips==1) {
            z[2] = z[1] + rseg;            // cell boundary above the exit point
            zpath[1] = sqrt((rfout-rfin)**2+(0.0300)**2);   // single-cell path length
            m31->Fill(1e4*zpath[1]);
            path_total = zpath[1];
	  }
	  else {                                  // 2 or more strips
	    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;
	    }
            sinalpha = (rfout-rfin)/sqrt((rfout-rfin)**2+ 0.0300**2);
            zpath[1]       = (z[2]-rfin      )/sinalpha;
            zpath[nstrips] = (rfout-z[nstrips])/sinalpha;
            path_total = path_total + zpath[1] + zpath[nstrips];
            for (int ii=2; ii<nstrips; ii++) {
              zpath[ii] = rseg/sinalpha;
              path_total = path_total + zpath[ii];
            }
	  }                                      // 2 or more strips hit

	  if (nstrips>=2) {
            for (int ii=2; ii<=nstrips-1; ii++) {
              m33->Fill(1e4*zpath[ii], float(nstrips));
            }
            m34->Fill(1e4*zpath[1]      , float(nstrips));
            m35->Fill(1e4*zpath[nstrips], float(nstrips));
          }

          m37->Fill(1e4*path_total,nstrips);
          if (nstrips<1) {
            cout<<"ERROR 1: nstrips<1 "<<nstrips<<endl;
	  }
//--- first layer threshold cut: ----------------------

          path_total = 0.;
          if (zpath[1]<zcut) {        // remove 1st strip if below threshold
            //cout<<"a-nstrips, zpath[1,2] = "<<nstrips<<" "<<zpath[1]<<" "<<zpath[2]<<endl;
            nstrips = nstrips -1;
            for (int ii=1; ii<=nstrips+1; ii++) {
              z[ii] = z[ii+1];
              zpath[ii] = zpath[ii+1];
            }                         // end remove 1st strip
            //cout<<"b-nstrips, zpath[1,2] = "<<nstrips<<" "<<zpath[1]<<" "<<zpath[2]<<endl;
          }
          if (zpath[nstrips]<zcut){   // remove last strip if below threshold
            nstrips = nstrips -1;
          }
          m30a->Fill(float(nstrips));    // after threshold cuts, also filled for last layer

          for (int ii=1; ii<=nstrips; ii++) {
            m32->Fill(1e4*zpath[ii]);
          } 
          if (nstrips<1) {
            cout<<"ERROR 2: nstrips<1 "<<nstrips<<endl;
	  }

          rfcluster = (z[1] + z[nstrips+1])/2.0;
          if (nstrips==1) {m17->Fill(1e4*(rfcluster-rfreal));}
          if (nstrips==2) {m18->Fill(1e4*(rfcluster-rfreal));}
          if (nstrips==3) {m19->Fill(1e4*(rfcluster-rfreal));}
          if (nstrips==4) {m20->Fill(1e4*(rfcluster-rfreal));}
          if (nstrips==5) {m21->Fill(1e4*(rfcluster-rfreal));}

//------- clustering last layer hit: --------------

          rlast  = sqrt(xlast [1]**2+xlast [2]**2);      // restore last

          for (int ii=0; ii<=99; ii++) { z[ii]=0.; zpath[ii]=0.; }  //clear

          rlin  = rlast - 0.5*(xlast[5] - xlast[4]);    // there is some ambiguity
          rlout = rlast + 0.5*(xlast[5] - xlast[4]);    // about the order of in/out
          if (rlin>rlout) {
            rtmp  = rlin;
            rlin  = rlout;
            rlout = rtmp;
	  }
          mxxx->Fill(rlin);  //xxx

          irstrip = rlin/rseg;             // integerize the entry point
          z[1] = rseg*irstrip;             // cell boundary below the entry point
          jrstrip = rlout/rseg + 1;         // integerize the exit point
	  nstrips =  jrstrip - irstrip;    //
          if (nstrips<1 || nstrips>99) {
            cout<<"rfin, rfout, n ,i , j strips= "<<rfin<<" "<<rfout<<" "<<nstrips<<" "<<irstrip<<" "
            <<jrstrip<<endl;    
          }
          m30->Fill(float(nstrips));    // this one is also filled for first layer
          //cout <<"rfin, rfout, irstrip, jrstrip, nstrips "<<rfin<<" "<<rfout<<" "
	  //   <<irstrip<<" "<<jrstrip<<" "<<nstrips<<endl;

	  if (nstrips==1) {
            z[2] = z[1] + rseg;            // cell boundary above the exit point
            zpath[1] = sqrt((rlout-rlin)**2+(0.0300)**2);   // single-cell path length
            m31a->Fill(1e4*zpath[1]);
	  }
	  else {                                  // 2 or more strips
	    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;
	    }
            sinalpha = (rlout-rlin)/sqrt((rlout-rlin)**2+ 0.0300**2);
            zpath[1]       = (z[2]-rlin      )/sinalpha;
            zpath[nstrips] = (rlout-z[nstrips])/sinalpha;
            for (int ii=2; ii<nstrips; ii++) {
              zpath[ii] = rseg/sinalpha;
            }
	  }                                      // 2 or more strips hit

	  if (nstrips>=2) {
            for (int ii=2; ii<=nstrips-1; ii++) {
              m33a->Fill(1e4*zpath[ii], float(nstrips));
            }
            m34a->Fill(1e4*zpath[1]      , float(nstrips));
            m35a->Fill(1e4*zpath[nstrips], float(nstrips));
          }

          if (nstrips<1) {
            //cout<<"ERROR 3: nstrips<1 "<<nstrips<<endl;
	  }
//--- threshold cut last plane: ----------------------

          if (zpath[1]<zcut) {        // remove 1st strip if below threshold
            nstrips = nstrips -1;
            for (int ii=1; ii<=nstrips+1; ii++) {
              z[ii] = z[ii+1];
              zpath[ii] = zpath[ii+1];
            }                         // end remove 1st strip
          }
          if (zpath[nstrips]<zcut){   // remove last strip if below threshold
            nstrips = nstrips -1;
          }
          m30a->Fill(float(nstrips));    // after threshold cuts, also filled for first layer

          for (int ii=1; ii<=nstrips; ii++) {
            m32a->Fill(1e4*zpath[ii]);
          } 
          if (nstrips<1) {
            //cout<<"ERROR 4: nstrips<1 "<<nstrips<<endl;
	  }

          rlcluster = (z[1] + z[nstrips+1])/2.0;
          if (nstrips==1) {m17a->Fill(1e4*(rlcluster-rlreal));}
          if (nstrips==2) {m18a->Fill(1e4*(rlcluster-rlreal));}
          if (nstrips==3) {m19a->Fill(1e4*(rlcluster-rlreal));}
          if (nstrips==4) {m20a->Fill(1e4*(rlcluster-rlreal));}
          if (nstrips==5) {m21a->Fill(1e4*(rlcluster-rlreal));}


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


//--- track with first, last clusters -----------------------------------
          // note this works with flat planes only for now: 
          if ( (rfcluster-rlcluster) != 0.0) {
            zv = (rfcluster*xlast[3] - rlcluster*xfirst[3]) / (rfcluster-rlcluster) - zvertex;
          }   // z vertex
          m12a->Fill(zv);
          if (fabs(zv)<0.1) {
            m13a->Fill(zv);
          }

//xxx


        }               // >= 3 planes hit
      }                 // loop over all tracks
      mult0 = 1;                   // reset hit multiplicity counter
      //      hit[][][]
    }
    else {                                 // we're in the same event
      //cout<<"3) EVENT, oldevent= "<<EVENT<<" "<<oldevent<<endl;
      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;
      hit[TRACK][LAYER][6] = XLOCIN;
      hit[TRACK][LAYER][7] = XLOCOUT;
      m99->Fill(XLOCOUT);
      //cout <<"zlocout: "<<hit[TRACK][LAYER][5]<<endl;
      zvertex   = Z_VERTEX;
      track_max = TRACK;
      //      rglobal = sqrt(XGLOBAL**2 + YGLOBAL**2);
      //if (LAYER==9) {
      //cout<<"b) EVENT: "<<EVENT<<", track: "<<TRACK<<", layer: "<<LAYER
      //   <<".  zloc in, out= "<<ZLOCIN<<" "<<ZLOCOUT <<"  ZLOCOUT-ZLOCIN"<<ZLOCOUT-ZLOCIN <<endl;
      //}
      if (LAYER==12 && ZGLOBAL>36.5 && sili_endcap_type==2){
	cout<<"***********************************************************"<<endl
            <<"***   error: cone data file, but flat switch set        ***"<<endl
            <<"***          set sili_endcap_type to 1                  ***"<<endl
            <<"***********************************************************"<<endl<<endl;
      }
      float dist = sqrt(ZGLOBAL**2 +YGLOBAL**2 + XGLOBAL**2);

      if (LAYER==9 && dist<19.5 && sili_endcap_type==1){
	cout<<"***********************************************************"<<endl
            <<"***   error: flat data file, but cone switch set        ***"<<endl
            <<"***          set sili_endcap_type to 2                  ***"<<endl
            <<"***********************************************************"<<endl<<endl;
      } 
    }
  }               // loop over records

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