{
//////////////////////////////////////////////////////////
//   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][9], 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 = 3;           // 1=shades, 2=flat

//                          // for integerization:
    Float_t rseg = 0.0050;  // flat: 50 microns

//--- --- threshold cuts:
    Float_t zcut = 0.0000;

  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);
  int nbins = 200;
  TH1F *v00 = new TH1F("v00","tan(theta) at z=0.",nbins, 0.0, 1.0);
  TH1F *v01 = new TH1F("v01","tan(theta) at z=10.",nbins, 0.0, 1.0);
  TH1F *v01a= new TH1F("v01a","tan(theta) at delz=0",nbins, 0.0, 1.0);
  TH1F *v01b= new TH1F("v01a","tan(theta) at delz=1",nbins, 0.0, 1.0);
  TH1F *v02 = new TH1F("v02","original event vertex", 100, -10, 10.0);
  int nscan = 1000;
  float offset_min = -10.0;
  float offset_max =  10.0;
  TH1F *v10 = new TH1F("v10","vertex sum", nscan, offset_min, offset_max);
  TH1F *v10a= new TH1F("v10a","vertex sum for event_stop", nscan, offset_min, offset_max);
  TH1F *v10b= new TH1F("v10b","vertex sum, flattened", nscan, offset_min, offset_max);
  TH1F *v11 = new TH1F("v11","vertex delta", 100, -10.0, 10.0);
  TH1F *v11b= new TH1F("v11b","vertex delta (b)", 100, -10.0, 10.0);
  TH1F *v15 = new TH1F("v15","vertex found", 100, -10.0, 10.0);
  TH1F *v15b= new TH1F("v15b","vertex found (b)", 100, -10.0, 10.0);
  TH2F *v12 = new TH2F("v12","bin vs event", 80, 460.0, 540.0, 21, 1.0, 21.0);
  TH2F *v13 = new TH2F("v13","x vs y layer 10 event 3", 40, -20.,20., 40, -20.,20.);
  TH1F *v14 = new TH1F("v14","tracks per event", 100, 0.0, 150.0);
  TF1 *f1 = new TF1("f1","pol1",1,0);

//----------------------xxx------------------------------------------
  icount = 0; 
//  nentries=3000;     //xxx   500 records is enough for 1st event
  int event_stop=99;
  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 (EVENT != oldevent && EVENT<=event_stop) {    // new event, starting at events 2!
//    if (EVENT != oldevent ) {    // new event, starting at events 2!

      v01->Reset();              // for this event
      v10->Reset(); v10a->Reset(); v15->Reset();; v15b->Reset();
      v02->Fill(Z0_EVENT);
      //float z0_event = Z0_EVENT;

      float trial_deltaz = 99.0;

      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 );

//--- start vertexing --------------------------------------------------------------------
      int iscan      = 0;
      float offset   = 99.;
      float tantheta = 99.;
      v14->Fill(track_max);
      printf("track_max: %f \n", track_max);

      for (offset = offset_min; offset < offset_max; 
           offset += (offset_max-offset_min)/nscan) {     // loop over vertex offsets
        //printf("offset: %f \n",offset);
        v01->Reset();
        float z0_event = hit[1][9][8];
        for (itrack=1; itrack<=track_max; itrack++ ) {       // loop over tracks
          for (ilayer=9; ilayer<=12; ilayer++) {             // loop over layers
            if (hit[itrack][ilayer][3] != 0.0) {             // 3 = zglobal
	      //              tantheta = sqrt(hit[itrack][ilayer][1]**2 + 
	      //              hit[itrack][ilayer][2]**2 ) /
	      //              hit[itrack][ilayer][3];
              //v00->Fill(tantheta);
	      //offset = 1.0;
              tantheta = sqrt((hit[itrack][ilayer][1])**2 + 
                              (hit[itrack][ilayer][2])**2) /
                              (hit[itrack][ilayer][3]+offset);
	      v01->Fill(tantheta);              
              if (ilayer==10 && EVENT==3 && offset==offset_min){
                v13->Fill(hit[itrack][ilayer][1],hit[itrack][ilayer][2]);
              }
            }          // no div by 0
          }            // loop over layers 9 - 12
        
        }              // loop over tracks in this event
        
        //if (fabs(offset-z0_event) < trial_deltaz) {
        if (fabs(offset) <= 0.001 && EVENT==8) {
          printf("offset %f  z0_event %f trial_deltaz %f \n", offset, z0_event, trial_deltaz);
          trial_deltaz = fabs(offset-z0_event);
          v01a->Reset();
          for (int ib=0; ib<nbins; ib++) {
	  v01a->SetBinContent(ib,v01->GetBinContent(ib));
	  }
        }
        if (fabs(offset-1.0) <= 0.001 && EVENT==8) {
          v01b->Reset();
          for (int ib=0; ib<nbins; ib++) {
	  v01b->SetBinContent(ib,v01->GetBinContent(ib));
	  }
        }

        int nhits = 0;
        float sum = 0.;
        float avg = 0.;
        float val = 0.;
        for (int ibin=0; ibin<nbins; ibin++) {
          val = v01->GetBinContent(ibin);
          if (val != 0) {
            nhits++;
            sum = sum + val;
	    }          // nonzero bin
	  }            // loop over bins
        if (nhits != 0) { avg = sum/nhits; }
        //printf("In v01: event: %d, nhits %d, sum %f, avg %f \n", EVENT, nhits, sum, avg);
        iscan++;
        v10->SetBinContent(iscan,avg);
        if (EVENT==event_stop) {v10a->SetBinContent(iscan,avg);}
        
      }    // end vertex scan offset

      float maxbin = v10->GetMaximumBin();
      v12->Fill(maxbin, EVENT);
      float zvertex = offset_min + (maxbin-1)*(offset_max - offset_min)/nscan;
      v15->Fill(zvertex);
      printf("maxbin: %f zvertex: %f \n",maxbin, zvertex);
      float delz = z0_event - zvertex;
      v11->Fill(delz);

      v10->Fit("f1","NQ");  
      float slope = f1->GetParameter(1);
      //printf("intercept: %f, slope: %f \n", intercept, slope);   
      for (int j=0; j<nscan; j++) {
        float flatval = v10->GetBinContent(j);
        flatval = flatval - slope*( (j-1)*(offset_max - offset_min)/nscan);
        v10b->SetBinContent(j,flatval);
      }
      float maxbinb = v10b->GetMaximumBin();
      float zvertexb = -(offset_min + (maxbinb-1)*(offset_max - offset_min)/nscan);
      v15b->Fill(zvertexb);
      printf("z0_event: %f, zvertexb: %f, delzb %f.\n",z0_event,zvertexb,delzb);
      float delzb = z0_event - zvertexb;
      printf("z0_event: %f, zvertexb: %f \n", z0_event, zvertexb);
      v11b->Fill(delzb);
    }                                      // end analyze previous event
//-----------------------------------------// ------------------------------
    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;
      hit[TRACK][LAYER][8] = Z0_EVENT;
      //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