{
//////////////////////////////////////////////////////////
//   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 =0;
  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=  TMath::Pi();
  Float_t degtorad = pi/180.;
Float_t r_global, s_global, the_global, phi_global, del_the, del_phi, dr_the, dr_phi, resid;

TH1F *e01 = new TH1F("e01","event multplicity",50, 0.0, 30000.0);
TH1F *e02 = new TH1F("e02","PTOT",50, 0.0, 20.0);
TH1F *e03 = new TH1F("e03","PTHETA",90, 0.0, 180.0);
TH1F *e04 = new TH1F("e04","PPHI",90, -180.0, 180.0);
TH1F *e05 = new TH1F("e05","eta",80, -4.0, 4.0);
TH2F *e06 = new TH2F("e06","eta vs theta",90, 0.0, 180.0,  90, -4.0, 4.0);
TH2F *e07 = new TH2F("e07","eta vs phi",100, -200.0, 200.0,  90, -4.0, 4.0);
TH2F *e08 = new TH2F("e08","eta vs phi",100, -200.0, 200.0,  90, -4.0, 4.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 *r01 = new TH1F("r01","theta from x,y,zglobal", 50, .0, 1.5);
TH1F *r02 = new TH1F("r02","phi from x,y,zglobal", 50, -1.0, 1.0);
TH1F *r03 = new TH1F("r03","del_theta", 50, -0.002, 0.002);
TH1F *r04 = new TH1F("r04","del_phi  ", 50, -0.005, 0.015);
TH1F *r05 = new TH1F("r05","del due to del_theta  ", 50, -0.05, 0.05);
TH1F *r06 = new TH1F("r06","del due to del_phi  ", 50, -0.05, 0.05);
TH1F *r07 = new TH1F("r07","residual  ", 50, -0.05, 0.05);

TH1F *r091 = new TH1F("r091","layer 9 del due to del_theta  ", 50, -0.05, 0.05);
TH1F *r092 = new TH1F("r092","layer 9 del due to del_phi  ", 50, -0.05, 0.05);
TH1F *r101 = new TH1F("r101","layer 10 del due to del_theta  ", 50, -0.05, 0.05);
TH1F *r102 = new TH1F("r102","layer 10 del due to del_phi  ", 50, -0.05, 0.05);
TH1F *r111 = new TH1F("r111","layer 11 del due to del_theta  ", 50, -0.05, 0.05);
TH1F *r112 = new TH1F("r112","layer 11 del due to del_phi  ", 50, -0.05, 0.05);
TH1F *r121 = new TH1F("r121","layer 12 del due to del_theta  ", 50, -0.05, 0.05);
TH1F *r122 = new TH1F("r122","layer 12 del due to del_phi  ", 50, -0.05, 0.05);


  for (Long64_t i=0; i<nentries;i++) {            // master loop over records
    oldz = ZGLOBAL;                               // save previous z
    nbytes += AncSvx->GetEntry(i);                // pull in the record
//    cout<<" event, track, layer, z = "<<EVENT<<" "<<TRACK<<" "<<LAYER<<"       "<<ZGLOBAL<<endl;
    if (EVENT != oldevent) {                      // new event
      oldevent = EVENT;
      cout<<"New event seen: "<<EVENT<<endl;
      if (oldevent != 1) {        // do end-of-event stats
        cout<<"#of records for event "<<EVENT-1<<" was "<<mult0<<endl;
        e01->Fill( float(mult0) );
      }
      mult0 = 1;                   // reset multiplicity counter

    }
    else {                                        // we're in the same event
      if (TRACK != oldtrack) {                    // here, don't double-count tracks
        oldtrack = TRACK;
        mult0++;
        //cout <<"mult0 = "<<mult0<<endl;

        z_global = ZGLOBAL - Z_VERTEX;
        r_global = sqrt(XGLOBAL*XGLOBAL+YGLOBAL*YGLOBAL);
        s_global = sqrt(XGLOBAL*XGLOBAL+YGLOBAL*YGLOBAL+z_global*z_global);
        the_global = TMath::ATan(r_global/z_global);   // in radians
        phi_global = TMath::ATan(YGLOBAL/XGLOBAL);    // in radians
       
        del_the = PTHETA*degtorad - the_global;
        del_phi =   PPHI*degtorad - phi_global;
       
        dr_the = del_the * s_global;
        dr_phi = del_phi * s_global * sin(PTHETA*degtorad);
         
        resid = sqrt(dr_the*dr_the+dr_phi*dr_phi);
        
        switch (LAYER) {
	  case 9:
            r091->Fill(dr_the);
            r092->Fill(dr_phi);
            break;
	  case 10:
            r101->Fill(dr_the);
            r102->Fill(dr_phi);
            break;
	  case 11:
            r111->Fill(dr_the);
            r112->Fill(dr_phi);
            break;
	  case 12:
            r121->Fill(dr_the);
            r122->Fill(dr_phi);
            break;
	}


        if (LAYER==8) {
//        cout<<"track, layer: "<<TRACK<<", "<<LAYER<<" x,y,z local= "
//                              <<XLOCIN<<" "<<YLOCIN<<" "<<ZLOCIN<<endl;
          g01->Fill(XLOCIN);
          g02->Fill(YLOCIN);
          g03->Fill(ZLOCIN);
          g04->Fill(ZLOCIN,XLOCIN);
          g05->Fill(XGLOBAL,YGLOBAL);
          g05a->Fill(XGLOBAL,YGLOBAL);
          g06->Fill(ZGLOBAL,sqrt(XGLOBAL**2+YGLOBAL**2));
        }
        else if (LAYER==12) {
          g11->Fill(XLOCIN);
          g12->Fill(YLOCIN);
          g13->Fill(ZLOCIN);
          g14->Fill(ZLOCIN,XLOCIN);
          g15->Fill(XGLOBAL,YGLOBAL);
          g15a->Fill(XGLOBAL,YGLOBAL);
          g16->Fill(ZGLOBAL,sqrt(XGLOBAL**2+YGLOBAL**2));
        }

        if (R_VERTEX == 0) {       // primary particle
	  e02->Fill(PTOT);
	  e03->Fill(PTHETA);
	  e04->Fill(PPHI);
          sign = PTHETA/fabs(PTHETA);
          eta = sign*log( fabs(tan(pi*PTHETA/360.0)) );
          e05->Fill(eta);
          e06->Fill(PTHETA,eta);
          e07->Fill(PPHI,eta);
          if (fabs(PPHI)<45.0) {e08->Fill(PPHI,eta);}
          //cout <<ITPARENT<<" "<<IDPARENT<<" "<<IDPART<<" "<<ITORIGIN<<" "<<IDORIGIN<<endl;
        }

      }                        // unique track
      else {                   // same track - successive 'hits' in same layer





        if (LAYER==8) {
        }
        else if (LAYER==12) {
          delz = ZGLOBAL-oldz;
          g27->Fill(delz);
          //cout<<"2) "<<EVENT<<" "<<TRACK<<" "<<LAYER<<" "<<ZGLOBAL<<" "<<oldz<<" "<<delz<<endl;
          //cout<<"3) x,y,z "<<XGLOBAL<<" "<<YGLOBAL<<" "<<ZGLOBAL<<" the, phi "<<PTHETA<<" "<<PPHI<<endl;
          r_global = sqrt(XGLOBAL*XGLOBAL+YGLOBAL*YGLOBAL);
          s_global = sqrt(XGLOBAL*XGLOBAL+YGLOBAL*YGLOBAL+ZGLOBAL*ZGLOBAL);
          the_global = TMath::ATan(r_global/ZGLOBAL);   // in radians
          phi_global = TMath::ATan(YGLOBAL/XGLOBAL);    // in radians
         
          del_the = PTHETA*degtorad - the_global;
          del_phi =   PPHI*degtorad - phi_global;
         
          dr_the = del_the * s_global;
          dr_phi = del_phi * s_global * sin(PTHETA*degtorad);
           
          resid = sqrt(dr_the*dr_the+dr_phi*dr_phi);
           
	  cout <<" PTHETA="<<PTHETA<<" PPHI="<<PPHI;
          cout <<" del_the, del_phi="<<del_the<<" "<<del_phi;
          cout <<" dr_the, dr_phi, resid= "<<dr_the<<" "<<dr_phi<<" "<<resid<<endl;
          
          r01->Fill(the_global);
          r02->Fill(phi_global);
          r03->Fill(del_the);
          r04->Fill(del_phi);
          r05->Fill(dr_the);
          r06->Fill(dr_phi);
          r07->Fill(resid);
        
          g21->Fill(XLOCIN);
          g22->Fill(YLOCIN);
          g23->Fill(ZLOCIN);
          g24->Fill(ZLOCIN,XLOCIN);
          g25 ->Fill(XGLOBAL,YGLOBAL);
          g25a->Fill(XGLOBAL,YGLOBAL);
          phi = TMath::ATan(XGLOBAL/YGLOBAL);              // range +- pi/2
          phi2 = float ( int(phi*10000) % 2618  )/10000.0;   // 0.2618 is 2pi / 24
          r   = sqrt(XGLOBAL**2+YGLOBAL**2);
          x2 = -r*sin(phi2);
          y2 = -r*cos(phi2);
          g25d->Fill(x2,y2);
          if (fabs(delz)>0.) {             // or use fabs()
            g25b->Fill(phi ,r);
            g25c->Fill(phi2,r);
            g25e->Fill(XGLOBAL,YGLOBAL);
          }
          if (abs(XGLOBAL)<0.3 && YGLOBAL< -10) {
            g28->Fill(delz);
	  }
          g26->Fill(ZGLOBAL,sqrt(r));
        }    // layer 12
      }      // same track
    }        // same event
  }          // loop over records
      cout<<"End of file seen: "<<EVENT<<endl;
      cout<<"#of records for event "<<EVENT<<" was "<<mult0<<endl;
      e01->Fill( float(mult0) );
}            // end of macro