{
//////////////////////////////////////////////////////////
//   Macro to analyze muon arm Si tracker ntuples
//   Steps : 1. Run Pisa
//           2. Run pisaRootRead
//           3. In root, run Track8.C (this macro)
//////////////////////////////////////////////////////////

//Reset ROOT and connect tree file, set up for use from WINDOWS
//   gROOT->Reset();
   #include <iostream.h>
   TFile *f = (TFile*)gROOT->GetListOfFiles()->FindObject("ancinr.root");
   if (!f) {
      f = new TFile("ancinr.root");
   }
   TTree *AncInr = (TTree*)gDirectory->Get("AncInr");

//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         NHIT;
   Float_t         Z0_EVENT;
   Float_t         B_IMPACT;
   Float_t         EVENT;

//Set branch addresses
   AncInr->SetBranchAddress("TRACK",&TRACK);
   AncInr->SetBranchAddress("NFILE",&NFILE);
   AncInr->SetBranchAddress("PTOT",&PTOT);
   AncInr->SetBranchAddress("PTHETA",&PTHETA);
   AncInr->SetBranchAddress("PPHI",&PPHI);
   AncInr->SetBranchAddress("R_VERTEX",&R_VERTEX);
   AncInr->SetBranchAddress("Z_VERTEX",&Z_VERTEX);
   AncInr->SetBranchAddress("THET_VER",&THET_VER);
   AncInr->SetBranchAddress("PHI_VER",&PHI_VER);
   AncInr->SetBranchAddress("ITPARENT",&ITPARENT);
   AncInr->SetBranchAddress("IDPARENT",&IDPARENT);
   AncInr->SetBranchAddress("IDPART",&IDPART);
   AncInr->SetBranchAddress("ITORIGIN",&ITORIGIN);
   AncInr->SetBranchAddress("IDORIGIN",&IDORIGIN);
   AncInr->SetBranchAddress("IHIT",&IHIT);
   AncInr->SetBranchAddress("LAYER",&LAYER);
   AncInr->SetBranchAddress("THETA",&THETA);
   AncInr->SetBranchAddress("PHI",&PHI);
   AncInr->SetBranchAddress("XGLOBAL",&XGLOBAL);
   AncInr->SetBranchAddress("YGLOBAL",&YGLOBAL);
   AncInr->SetBranchAddress("ZGLOBAL",&ZGLOBAL);
   AncInr->SetBranchAddress("PMOMX",&PMOMX);
   AncInr->SetBranchAddress("PMOMY",&PMOMY);
   AncInr->SetBranchAddress("PMOMZ",&PMOMZ);
   AncInr->SetBranchAddress("IPC123",&IPC123);
   AncInr->SetBranchAddress("PC1THETA",&PC1THETA);
   AncInr->SetBranchAddress("PC1PHI",&PC1PHI);
   AncInr->SetBranchAddress("PC1RDIST",&PC1RDIST);
   AncInr->SetBranchAddress("PC1ZDIST",&PC1ZDIST);
   AncInr->SetBranchAddress("PC2THETA",&PC2THETA);
   AncInr->SetBranchAddress("PC2PHI",&PC2PHI);
   AncInr->SetBranchAddress("PC2RDIST",&PC2RDIST);
   AncInr->SetBranchAddress("PC2ZDIST",&PC2ZDIST);
   AncInr->SetBranchAddress("PC3THETA",&PC3THETA);
   AncInr->SetBranchAddress("PC3PHI",&PC3PHI);
   AncInr->SetBranchAddress("PC3RDIST",&PC3RDIST);
   AncInr->SetBranchAddress("PC3ZDIST",&PC3ZDIST);
   AncInr->SetBranchAddress("NHIT",&NHIT);
   AncInr->SetBranchAddress("Z0_EVENT",&Z0_EVENT);
   AncInr->SetBranchAddress("B_IMPACT",&B_IMPACT);
   AncInr->SetBranchAddress("EVENT",&EVENT);

//     This is the loop skeleton
   Int_t nentries = AncInr->GetEntries();

// arrays of floats
Float_t event[nentries];
Float_t track[nentries];
Float_t ptot[nentries];
Float_t ptheta[nentries];
Float_t pphi[nentries];
Float_t theta[nentries];
Float_t phi[nentries];
Float_t x[nentries];
Float_t y[nentries];
Float_t z[nentries];
Float_t px[nentries];
Float_t py[nentries];
Float_t pz[nentries];
Float_t layer[nentries];
Float_t id[nentries];
Int_t ievent[10001];


   Int_t nbytes = 0;
   Float_t oldevent = 1.;
   Int_t iev = 0;
   ievent[1]=0;

// Loop over all ntuple entries
   for (Int_t i=0; i<nentries;i++) {
      nbytes += AncInr->GetEntry(i);
      //      cout<<EVENT<<" ";
// watch for new event number
      event[i]=EVENT;
      if (EVENT>oldevent) {
	oldevent=EVENT;
	iev=EVENT;
	ievent[iev]=i;
      }
// fill hit arrays
      track[i]=TRACK;
      ptot[i]=PTOT;
      ptheta[i]=PTHETA;
      pphi[i]=PPHI;
      theta[i]=THETA;
      phi[i]=PHI;
      x[i]=XGLOBAL;
      y[i]=YGLOBAL;
      z[i]=ZGLOBAL;
      px[i]=PMOMX;
      py[i]=PMOMY;
      pz[i]=PMOMZ;
      layer[i]=LAYER;
      id[i]=IDPART;
   }
//   cout<<endl;

// tracks entries below are track#(1-2) layer(1-12) (x,y,z,px,py,pz)
Float_t tracks[3][13][7];
// vertex entries below are track#(1-2) (xvertex,yvertex,zvertex,ptot,ptheta,pphi)
Float_t vertex[3][7];
Int_t idp = 0;
Int_t itrack = 0;
Int_t ilayer = 0;
Float_t pzvertex = 0.;
Float_t ptvertex = 0.;
Float_t pxvertex = 0.;
Float_t pyvertex = 0.;
Float_t degtorad = 3.14159/180.;
Int_t idpart[3];

//define histograms
TH1F *g1 = new TH1F("g1","xintercept id=5",50,-0.05,0.05);
TH1F *g2 = new TH1F("g2","yintercept id=5",50,-0.05,0.05);
TH1F *g3 = new TH1F("g3","xintercept id=6",50,-0.05,0.05);
TH1F *g4 = new TH1F("g4","yintercept id=6",50,-0.05,0.05);
TH1F *g5 = new TH1F("g5","zxvertex id=5",50,-0.2,0.2);
TH1F *g6 = new TH1F("g6","zyvertex id=5",50,-0.2,0.2);
TH1F *g7 = new TH1F("g7","zxvertex id=6",50,-0.2,0.2);
TH1F *g8 = new TH1F("g8","zyvertex id=6",50,-0.2,0.2);
TH1F *g9 = new TH1F("g9","xpair vertex",50,-0.2,0.2);
TH1F *g10 = new TH1F("g10","ypair vertex",50,-0.2,0.2);
TH1F *g11 = new TH1F("g11","pair vertex",50,-0.2,0.2);
TH1F *g12 = new TH1F("g12","pair mass",50,2.0,4.0);
TH1F *g13 = new TH1F("g13","tracker xmatch",50,-10.,10.);
TH1F *g14 = new TH1F("g14","tracker ymatch",50,-10.,10.);
TH1F *g15 = new TH1F("g15","tracker pzmatch",50,0.,2.5);
TH1F *g16 = new TH1F("g16","si dphi deg id=5",50,-1.,1.);
TH1F *g17 = new TH1F("g17","si dphi deg id=6",50,-1.,1.);
TH1F *g18 = new TH1F("g18","si dphi*pz degGeV",50,-5.,5.);
TH1F *g19 = new TH1F("g19","rintercept",50,-0.05,0.2);
TH1F *g20 = new TH1F("g20","zrvertex",50,0.2,1.);
// TH1F *g20 = new TH1F("g20","zrvertex",50,-0.1,1.);
TH1F *g21 = new TH1F("g21","Si rmatch",50,-0.02,0.02);
TH1F *g22 = new TH1F("g22","rvtx thrown",50,-0.05,0.2);
TH1F *g23 = new TH1F("g23","zvtx thrown",50,0.2,1.);
// TH1F *g23 = new TH1F("g23","zvtx thrown",50,-0.1,1.);
TH2F *g24 = new TH2F("g24","xintercept 1 vs 2",50,-0.05,0.05,50,-0.05,0.05);
TH2F *g25 = new TH2F("g25","yintercept 1 vs 2",50,-0.05,0.05,50,-0.05,0.05);
TH2F *g26 = new TH2F("g26","xint1 vs yint1",50,-0.05,0.05,50,-0.05,0.05);
TH2F *g27 = new TH2F("g27","xint2 vs yint2",50,-0.05,0.05,50,-0.05,0.05);
TH1F *g28 = new TH1F("g28","rint fit",50,-0.05,0.05);
TH1F *g29 = new TH1F("g29","rpair vertex",50,-0.1,0.1);
TH1F *g30 = new TH1F("g30","rpair intercept",50,-0.05,0.05);

//loop over events
for (Int_t ii=1; ii<iev; ii++) {
  //  cout<<ii<<" event"<<endl;

//zero out track and vertex arrays
  for (Int_t i=1; i<7; i++) {
    vertex[1][i]=0;
    vertex[2][i]=0;
    for (Int_t j=1; j<13;j++) {
      tracks[1][j][i]=0.;
      tracks[2][j][i]=0.;
    }
  }

//work on each event's hits, filling tracks and vertex arrays
  //  nbytes = AncInr->GetEntry(ievent(ii));
  //  cout <<" number of hits on tracks "<<ievent[ii+1]-ievent[ii]<<endl;

  for (Int_t i=ievent[ii]; i<ievent[ii+1];i++) {
    itrack=track[i];
    ilayer=layer[i];
    if (itrack>2) break;
    //    cout <<" track " <<itrack <<" ";
    //    cout <<" layer " <<ilayer <<endl;
    // store track info 
    idpart[itrack]=id[i];
    idp=id[i];
    tracks[itrack][ilayer][1] = x[i];
    tracks[itrack][ilayer][2] = y[i];
    tracks[itrack][ilayer][3] = z[i];
    tracks[itrack][ilayer][4] = px[i];
    tracks[itrack][ilayer][5] = py[i];
    tracks[itrack][ilayer][6] = pz[i];
    // store vertex info
    pzvertex = ptot[i]*cos(ptheta[i]*degtorad);
    ptvertex = ptot[i]*sin(ptheta[i]*degtorad);
    pxvertex = ptvertex*cos(pphi[i]*degtorad);
    pyvertex = ptvertex*sin(pphi[i]*degtorad);
    //vertex[itrack][1] = vx[ii][6];                     //vx,y,z are not defined or filled
    //vertex[itrack][2] = vy[ii][6];                     // ???
    //vertex[itrack][3] = vz[ii][6];
    vertex[itrack][4] = pxvertex;
    vertex[itrack][5] = pyvertex;
    vertex[itrack][6] = pzvertex;
  }

  Float_t xintercept[3];
  Float_t yintercept[3];
  Float_t xvtx = 0.;
  Float_t yvtx = 0.;
  Float_t xslope[3];
  Float_t yslope[3];
  Float_t xpairvtx = 0.;
  Float_t ypairvtx = 0.;
  Float_t pairvtx = 0.;
  Float_t rpairvtx = 0.;
  Float_t rpairint = 0.;
  Float_t pairmass = 0.;
  Float_t pairmom = 0.;
  Float_t pairpt = 0.;
  Float_t pairetot = 0.;
  Float_t etot[3];
// xseg also used for radius, yseg also used for phi*radius
  Float_t xseg = 0.00500;
  Float_t yseg = 0.02000;
  Float_t xmatch = 0.;
  Float_t ymatch = 0.;
  Float_t pzmatch = 0.;
  Float_t deltaphi = 0.;
  Float_t radius[5];
  Float_t phidet[5];
  Float_t rslope[3];
  Float_t rintercept[3];
  Float_t rvtx = 0.;
  Float_t rmatch = 0.;
  Float_t radvtx = 0.;
  Float_t vtxslope = 0.;
  Int_t goodtrk[3];
  Int_t xstrip = 0;
  Int_t ystrip = 0;
  Int_t rstrip = 0;
  Int_t phistrip = 0;

  Float_t fitx[4] = {13.0,21.0,29.0,37.0};
  Float_t fity[4];
  Float_t ex[4] = {.0001,.0001,.0001,.0001};
  Float_t ey[4] = {.001,.001,.001,.001};

// loop over tracks in each event
  for (Int_t i=1; i<3;i++) {
    xslope[i]=0.;
    yslope[i]=0.;
    xintercept[i]=0.;
    yintercept[i]=0.;
    rintercept[i]=0.; 
    rslope[i]=0.; 
    goodtrk[i]=0;

    // require hit in layer 3,4,5,6 (= Si),  and 11 (= MuTr entrance plane)
    if (tracks[i][3][3]!=0. && tracks[i][4][3]!=0. && tracks[i][5][3]!=0. && tracks[i][6][3]!=0. && tracks[i][11][3]!=0. && vertex[i][6]>2.5) {
      goodtrk[i]=1;
      //      cout<<"track "<<i<<" good"<<endl;
// compute x and y projected vertex
      //      cout<<tracks[i][3][1]<<" x at plane 3, z = "<<tracks[i][3][3]<<endl;
      //      cout<<tracks[i][3][2]<<" y at plane 3"<<endl;
      //      cout<<tracks[i][6][1]<<" x at plane 6, z = "<<tracks[i][6][3]<<endl;
      //      cout<<tracks[i][6][2]<<" y at plane 6"<<endl;
// integerize strip positions for planes
      //      xstrip = (tracks[i][3][1]/xseg + 0.5);
      //      ystrip = (tracks[i][3][2]/yseg + 0.5);
      //      tracks[i][3][1] = xseg*xstrip;
      //      tracks[i][3][2] = yseg*ystrip;
      //      xstrip = (tracks[i][4][1]/xseg + 0.5);
      //      ystrip = (tracks[i][4][2]/yseg + 0.5);
      //      tracks[i][4][1] = xseg*xstrip;
      //      tracks[i][4][2] = yseg*ystrip;
      //      xstrip = (tracks[i][5][1]/xseg + 0.5);
      //      ystrip = (tracks[i][5][2]/yseg + 0.5);
      //      tracks[i][5][1] = xseg*xstrip;
      //      tracks[i][5][2] = yseg*ystrip;
      //      xstrip = (tracks[i][6][1]/xseg + 0.5);
      //      ystrip = (tracks[i][6][2]/yseg + 0.5);
      //      tracks[i][6][1] = xseg*xstrip;
      //      tracks[i][6][2] = yseg*ystrip;

// calculate track slopes and intercepts at Z=0
// try using first two planes to limit multiple scattering
      //      if ((tracks[i][6][3]-tracks[i][3][3])!=0) xslope[i] = (tracks[i][6][1]-tracks[i][3][1]) / (tracks[i][6][3]-tracks[i][3][3]);
// use first and last planes t minimize resolution effects
      if ((tracks[i][4][3]-tracks[i][3][3])!=0) xslope[i] = (tracks[i][4][1]-tracks[i][3][1]) / (tracks[i][4][3]-tracks[i][3][3]);
      xintercept[i] = tracks[i][3][1] - xslope[i] * tracks[i][3][3];
      xvtx=-999.;
      if (xslope[i]!=0.) xvtx = -xintercept[i]/xslope[i];
      //      cout<<xintercept[i]<<" xintercept"<<endl;
      if (idpart[i]==5) g1->Fill(xintercept[i]);
      if (idpart[i]==6) g3->Fill(xintercept[i]);
      if (idpart[i]==5) g5->Fill(xvtx);
      if (idpart[i]==6) g7->Fill(xvtx);
      if ((tracks[i][4][3]-tracks[i][3][3])!=0) yslope[i] = (tracks[i][4][2]-tracks[i][3][2]) / (tracks[i][4][3]-tracks[i][3][3]);
      //      if ((tracks[i][6][3]-tracks[i][3][3])!=0) yslope[i] = (tracks[i][6][2]-tracks[i][3][2]) / (tracks[i][6][3]-tracks[i][3][3]);
      yintercept[i] = tracks[i][3][2] - yslope[i] * tracks[i][3][3];
      yvtx = -999.;      
      if (yslope[i]!=0.) yvtx = -yintercept[i]/yslope[i];
      //      cout<<yintercept[i]<<" yintercept"<<endl;
      if (idpart[i]==5) g2->Fill(yintercept[i]);
      if (idpart[i]==6) g4->Fill(yintercept[i]);
      if (idpart[i]==5) g6->Fill(yvtx);
      if (idpart[i]==6) g8->Fill(yvtx);
      if (goodtrk[1]==1 && goodtrk[2]==1) {
	g24->Fill(xintercept[1],xintercept[2]);
      	g25->Fill(yintercept[1],yintercept[2]);
	g26->Fill(xintercept[1],yintercept[1]);
      	g27->Fill(xintercept[2],yintercept[2]);
      }
// radial quantities
      // use xseg as radial segmentation, yseg as phi 
      if ((tracks[i][6][3]-tracks[i][3][3])!=0) {
        radius[1] = (tracks[i][3][1]**2 + tracks[i][3][2]**2)**0.5;
        radius[2] = (tracks[i][4][1]**2 + tracks[i][4][2]**2)**0.5;
        radius[3] = (tracks[i][5][1]**2 + tracks[i][5][2]**2)**0.5;
        radius[4] = (tracks[i][6][1]**2 + tracks[i][6][2]**2)**0.5;
        phidet[1] = atan(tracks[i][3][2]/tracks[i][3][1]);
        phidet[2] = atan(tracks[i][4][2]/tracks[i][4][1]);
        phidet[3] = atan(tracks[i][5][2]/tracks[i][5][1]);
        phidet[4] = atan(tracks[i][6][2]/tracks[i][6][1]);

	// integerize radius
 	rstrip = (radius[1]/xseg + 0.5);
	radius[1] = xseg*rstrip;
	rstrip = (radius[2]/xseg + 0.5);
	radius[2] = xseg*rstrip;
        rstrip = (radius[3]/xseg + 0.5);
        radius[3] = xseg*rstrip;
        rstrip = (radius[4]/xseg + 0.5);
        radius[4] = xseg*rstrip;
// integerize phi using arc length
	phistrip = (phidet[1]*radius[1]/yseg + 0.5);
	phidet[1] = yseg*phistrip/radius[1];
	phistrip = (phidet[2]*radius[2]/yseg + 0.5);
	phidet[2] = yseg*phistrip/radius[2];
        phistrip = (phidet[3]*radius[3]/yseg + 0.5);
        phidet[3] = yseg*phistrip/radius[3];
        phistrip = (phidet[4]*radius[4]/yseg + 0.5);
        phidet[4] = yseg*phistrip/radius[4];

        rslope[i] = (radius[4]-radius[1]) / (tracks[i][6][3]-tracks[i][3][3]);
        rintercept[i] = radius[1] - rslope[i] * tracks[i][3][3];
        if (rslope[i]!=0.) rvtx = -rintercept[i]/rslope[i];
        vtxslope = ((vertex[i][4]**2+vertex[i][5]**2)**0.5)/vertex[i][6];
        rmatch = radius[1] - vtxslope*tracks[i][3][3];
//        rmatch = rmatch + (radius[2] - vtxslope*tracks[i][6][3]);
//        cout<<rslope<<rvtx<<endl; 
        g19->Fill(rintercept[i]);
        g20->Fill(rvtx);
        g21->Fill(rmatch);
        radvtx = ((vertex[i][1]**2+vertex[i][2]**2)**0.5);          //xxx note vertex[][1,2,3] are not set
        g22->Fill(radvtx);                                          //xxx note radvtx is not set properly
        g23->Fill(vertex[i][3]);                                    //xxx note vertex[][3] is not set
// fit the radial intercept
//        fity[0]=radius[1];
//	fity[1]=radius[2];
//	fity[2]=radius[3];
//	fity[3]=radius[4];
//	gr = new TGraphErrors(4,fitx,fity,ex,ey);
//	gr->Fit("pol1");
//	TF1 *fit = gr->GetFunction("pol1");
//	Double_t chi2 = fit->GetChisquare();
//	Double_t p0 = fit->GetParameter(0);
//	Double_t e0 = fit->GetParError(0);
	//	cout << "intercept fit " << p0 << " chi2 "<< chi2 << endl;
//        g28->Fill(p0);
      }
// calculate track match from MuTr to last silicon plane
      xmatch = tracks[i][11][1] - (tracks[i][11][4]/tracks[i][11][6])*(tracks[i][11][3]-tracks[i][6][3]) - tracks[i][6][1];
      ymatch = tracks[i][11][2] - (tracks[i][11][5]/tracks[i][11][5])*(tracks[i][11][3]-tracks[i][6][3]) - tracks[i][6][2];
      pzmatch = tracks[i][6][6]-tracks[i][11][6];
      g13->Fill(xmatch);
      g14->Fill(ymatch);
      g15->Fill(pzmatch);
// calculate bend in phi, due to magnetic field and multiple scattering
      if (tracks[i][6][1]!=0. && tracks[i][3][1]!=0.) {
        deltaphi = phidet[4] - phidet[1];
        deltaphi = deltaphi/degtorad;
	if (idpart[i]==5) g16->Fill(deltaphi);
        if (idpart[i]==6) g17->Fill(deltaphi);
        deltaphi=deltaphi*tracks[i][6][6];
	g18->Fill(deltaphi);
      }
    }
  }

// calculate Z-vertex and kinematics for good track pairs
  if (goodtrk[1]==1&&goodtrk[2]==1) {
    cout<<"good pair \nl"
    if ((xslope[2]-xslope[1])!=0) xpairvtx = -(xintercept[2]-xintercept[1])/(xslope[2]-xslope[1]);
    if ((yslope[2]-yslope[1])!=0) ypairvtx = -(yintercept[2]-yintercept[1])/(yslope[2]-yslope[1]);
    rpairvtx = -(rintercept[1]/rslope[1] + rintercept[2]/rslope[2]) / 2. ;
    rpairint = (rintercept[1] + rintercept[2]) / 2. ;
    pairvtx = (xpairvtx+ypairvtx)/2.;
    //    pairpt = (tracks[1][3][4] + tracks[2][3][4])**2 + (tracks[1][3][5] + tracks[2][3][5])**2;
    pairpt = (xslope[1]*tracks[1][3][6] + xslope[2]*tracks[2][3][6])**2 + (yslope[1]*tracks[1][3][6] + yslope[2]*tracks[2][3][6])**2;
    pairpt = pairpt**0.5;    
    pairmom = (pairpt**2 + (tracks[1][3][6] + tracks[2][3][6])**2)**0.5;
    //    etot[1] = (0.10566**2 + tracks[1][3][4]**2 + tracks[1][3][5]**2 + tracks[1][3][6]**2)**0.5;
    //    etot[2] = (0.10566**2 + tracks[2][3][4]**2 + tracks[2][3][5]**2 + tracks[2][3][6]**2)**0.5;
    etot[1] = (0.10566**2 + (xslope[1]*tracks[1][3][6])**2 + (yslope[1]*tracks[1][3][6])**2 + tracks[1][3][6]**2)**0.5;
    etot[2] = (0.10566**2 + (xslope[2]*tracks[2][3][6])**2 + (yslope[2]*tracks[2][3][6])**2 + tracks[2][3][6]**2)**0.5;
    pairmass = ((etot[1]+etot[2])**2 - pairmom**2)**0.5;
    g9->Fill(xpairvtx);
    g10->Fill(ypairvtx);
    g11->Fill(pairvtx);
    g12->Fill(pairmass);
    g29->Fill(rpairvtx);
    g30->Fill(rpairint);
  }
}
// done
}