#define qcd_udsg_macro_cxx
#include "qcd_udsg_macro.h"
#include "TH2.h"
#include "TStyle.h"
#include "TCanvas.h"

void qcd_udsg_macro::Loop()
{
//   In a Root session, you can do:
//      Root > .L qcd_udsg_macro.C
//      Root > qcd_udsg_macro t
//      Root > t.GetEntry(12); // Fill t data members with entry number 12
//      Root > t.Show();       // Show values of entry 12
//      Root > t.Show(16);     // Read and show values of entry 16
//      Root > t.Loop();       // Loop on all entries
//

//     This is the loop skeleton
//       To read only selected branches, Insert statements like:
// METHOD1:
//    fChain->SetBranchStatus("*",0);  // disable all branches
//    fChain->SetBranchStatus("branchname",1);  // activate branchname
// METHOD2: replace line
//    fChain->GetEntry(i);  // read all branches
//by  b_branchname->GetEntry(i); //read only this branch
   if (fChain == 0) return;

//
// This macro loops over all values of eta for the calorimeter 
//  (-4.6 to 4.6, in 74 steps, mostly of 0.1 eta), and 
//  recalculates the missing energy for an 
//  event, with respect to a user-supplied vertex.
//
   TH1F *Hist1 = new TH1F("Hist1","MEx",100,-75,75);
   TH1F *Hist2 = new TH1F("Hist2","MEy",100,-75,75);
   TH1F *Hist3 = new TH1F("Hist3","RETx",100,-75,75);
   TH1F *Hist4 = new TH1F("Hist4","RETy",100,-75,75);
   TH1F *Hist5 = new TH1F("Hist5","MEx - RETx",100,-75,75);
   TH1F *Hist6 = new TH1F("Hist6","MEy - RETy",100,-75,75);
   TH1F *Hist7 = new TH1F("Hist7","MET",50,0,75);
   TH1F *Hist8 = new TH1F("Hist8","RET",50,0,75);
   TH1F *Hist9 = new TH1F("Hist9","(MET - RET)/MET",200,-2,2);
   TH1F *Hist10 = new TH1F("Hist10","MET - RET",100,-75,75);
   TH1F *Hist11 = new TH1F("Hist11","NEW_MET",50,0,75);
   TH1F *Hist12 = new TH1F("Hist12","(MET - NEW_MET)/MET",200,-2,2);
   TH1F *Hist13 = new TH1F("Hist13","MET - NEW_MET",100,-75,75);
   TH1F *Hist14 = new TH1F("Hist14","MET_Zvtx",100,-75,75);
   TH1F *Hist15 = new TH1F("Hist15","VTX_Z",100,-75,75);

   TH1F *Hist21 = new TH1F("Hist21","MET_Zvtx",100,-750,750);
   TH1F *Hist22 = new TH1F("Hist22","MET - RET",100,-20,20);    
   TH1F *Hist23 = new TH1F("Hist23","MET - RET",100,-20,20);
   TH1F *Hist24 = new TH1F("Hist24","(MET - RET)/MET",200,-2,2);
   TH1F *Hist25 = new TH1F("Hist25","(MET - RET)/MET",200,-2,2);
   TH1F *Hist26 = new TH1F("Hist26","MET - RET",100,-20,20);
   TH1F *Hist27 = new TH1F("Hist27","MET - RET",100,-20,20);
   TH1F *Hist28 = new TH1F("Hist28","MET - RET",100,-20,20);
   TH1F *Hist29 = new TH1F("Hist29","MET - RET",100,-20,20);
   TH1F *Hist30 = new TH1F("Hist30","(MET - RET)/MET",200,-2,2);
   TH1F *Hist31 = new TH1F("Hist31","(MET - RET)/MET",200,-2,2);
   TH1F *Hist32 = new TH1F("Hist32","(MET - RET)/MET",200,-2,2);
   TH1F *Hist33 = new TH1F("Hist33","(MET - RET)/MET",200,-2,2);
   TH1F *Hist34 = new TH1F("Hist34","NEW_MET - RET",100,-20,20);
   TH1F *Hist35 = new TH1F("Hist35","COMBTET",50,0,75);
   TH1F *Hist36 = new TH1F("Hist36","MET - COMBTET",100,-20,20);
   TH1F *Hist37 = new TH1F("Hist37","NEW_MET - COMBTET",100,-20,20);
   TH1F *Hist38 = new TH1F("Hist38","(NEW_MET - RET)/NEW_MET",200,-2,2);
   TH1F *Hist39 = new TH1F("Hist39","(MET - COMB_MET)/MET",200,-2,2);
   TH1F *Hist40 = new TH1F("Hist40","(NEW_MET - COMB_MET)/NEW_MET",200,-2,2);
   TH2F *Hist41 = new 
         TH2F("Hist41","MET - RET vs Z vertex",80,-20,20,100,-75,75);
   TH2F *Hist42 = new 
         TH2F("Hist42","(MET - RET)/MET vs Z vertex",80,-2,2,100,-75,75);
   TH1F *Hist43 = new TH1F("Hist43","MET - NEW_MET",100,-20,20);
   TH1F *Hist44 = new TH1F("Hist44","NRET_EMX",100,-20,20);
   TH1F *Hist45 = new TH1F("Hist45","NRET_HDX",100,-20,20);
   TH1F *Hist46 = new TH1F("Hist46","NRET_EMY",100,-20,20);
   TH1F *Hist47 = new TH1F("Hist47","NRET_HDY",100,-20,20);
   TH1F *Hist48 = new TH1F("Hist48","COMBVERTx",100,-75,75);
   TH1F *Hist49 = new TH1F("Hist49","COMBVERTy",100,-75,75);
   TH1F *Hist50 = new TH1F("Hist50","THETPRIME",100,-6,6);
   TH1F *Hist51 = new TH1F("Hist51","THETNOT",100,-6,6);
   TH1F *Hist52 = new TH1F("Hist52","THETPRIME-THENOT",100,-4,4);
   TH1F *Hist53 = new TH1F("Hist53","HADTHETPRIME",100,-6,6);
   TH1F *Hist54 = new TH1F("Hist54","HADTHETNOT",100,-6,6);
   TH1F *Hist55 = new TH1F("Hist55","HADTHETPRIME-HADTHENOT",100,-4,4);
   TH1F *Hist56 = new TH1F("Hist56","COMBVERTEMx",100,-75,75);
   TH1F *Hist57 = new TH1F("Hist57","COMBVERTEMy",100,-75,75);
   TH1F *Hist58 = new TH1F("Hist58","COMBVERTHDx",100,-75,75);
   TH1F *Hist59 = new TH1F("Hist59","COMBVERTHDy",100,-75,75);
   TH1F *Hist60 = new TH1F("Hist60","NRETEMx",100,-20,20);
   TH1F *Hist61 = new TH1F("Hist61","NRETEMy",100,-20,20);
   TH1F *Hist62 = new TH1F("Hist62","NRETHDx",100,-20,20);
   TH1F *Hist63 = new TH1F("Hist63","NRETHDy",100,-20,20);
   TH1F *Hist64 = new TH1F("Hist64","TRICOMBTET",50,0,75);
   TH1F *Hist65 = new TH1F("Hist65","MET - TRICOMBTET",100,-20,20);
   TH1F *Hist66 = new TH1F("Hist66","(MET - TRICOMB_MET)/MET",200,-2,2);
   TH1F *Hist67 = new TH1F("Hist67","NEW_MET - TRICOMBTET",100,-20,20);
   TH1F *Hist68 = new 
         TH1F("Hist68","(NEW_MET - TRICOMB_MET)/NEW_MET",200,-2,2);
   TH1F *Hist69 = new TH1F("Hist69","COMBTET - TRICOMBTET",100,-20,20);
   TH1F *Hist70 = new 
         TH1F("Hist70","(COMBTET - TRICOMB_MET)/COMBTET",200,-2,2);


   Int_t nentries = Int_t(fChain->GetEntries());

   Int_t nbytes = 0, nb = 0;

   Double_t TOTRETx = 0;			
   Double_t TOTRETy = 0;			
   Double_t TOTEMRETx = 0;			
   Double_t TOTEMRETy = 0;			
   Double_t TOTHDRETx = 0;			
   Double_t TOTHDRETy = 0;			
   Double_t TOTRET = 0;
   Double_t NEWVERTx = 0;
   Double_t NEWVERTy = 0; 
   Double_t NEWVERTEMx = 0;
   Double_t NEWVERTEMy = 0; 
   Double_t NEWVERTHDx = 0;
   Double_t NEWVERTHDy = 0; 
   Double_t NEWTET = 0;
   Double_t THETPRIME = 0;
   Double_t THETNOT = 0;
   Double_t HADTHETPRIME = 0;
   Double_t HADTHETNOT = 0;
   Double_t VTX_Z = 0;
   Double_t eta = 0;
   Double_t COMBVERTx =0;
   Double_t COMBVERTy = 0;
   Double_t COMBVERTEMx =0;
   Double_t COMBVERTEMy = 0;
   Double_t COMBVERTHDx =0;
   Double_t COMBVERTHDy = 0;
   Double_t COMBTET = 0;
   Double_t VERTx =0;
   Double_t VERTy = 0;
   Double_t VERTEMx =0;
   Double_t VERTEMy = 0;
   Double_t VERTHDx =0;
   Double_t VERTHDy = 0;
   Double_t VERTET = 0;
   Double_t AVE_EMR = 0;
   Double_t AVE_EMX = 0;
   Double_t AVE_HADR = 0;
   Double_t AVE_HADX = 0;
   Double_t NRET_EMX = 0;
   Double_t NRET_EMY = 0;
   Double_t NRET_HDX = 0;
   Double_t NRET_HDY = 0;

   Double_t X_EM[74];
   Double_t X_HAD[74];
   Double_t EM_R[74];   
   Double_t HAD_R[74]; 
// 
   EM_R[0] = 7.2;   //guesstimate (These are the average heights for the
   EM_R[1] = 7.2;   //guesstimate    electromagnetic part of the calorimeter.)
   EM_R[2] = 7.6;
   EM_R[3] = 10.4;
   EM_R[4] = 13.2;
   EM_R[5] = 15.5;
   EM_R[6] = 17.2;
   EM_R[7] = 19.0;
   EM_R[8] = 21.0;
   EM_R[9] = 23.2;
   EM_R[10] = 25.7;
   EM_R[11] = 28.4;
   EM_R[12] = 31.4;
   EM_R[13] = 34.8;
   EM_R[14] = 38.6;
   EM_R[15] = 42.7;
   EM_R[16] = 47.4;
   EM_R[17] = 52.5;
   EM_R[18] = 58.3;
   EM_R[19] = 64.8;
   EM_R[20] = 72.1;
   EM_R[21] = 80.4;
   EM_R[22] = 88.0;
   EM_R[23] = 92.6;  //guesstimate
   EM_R[24] = 90.2;  //guesstimate
   EM_R[25] = 90.2;  //guesstimate
   EM_R[26] = 92.6;
   EM_R[27] = 92.6;
   EM_R[28] = 92.6;
   EM_R[29] = 92.6;
   EM_R[30] = 92.6;
   EM_R[31] = 92.6;
   EM_R[32] = 92.6;
   EM_R[33] = 92.6;
   EM_R[34] = 92.6;
   EM_R[35] = 92.6;
   EM_R[36] = 92.6;
   EM_R[37] = 92.6;
   EM_R[38] = 92.6;
   EM_R[39] = 92.6;
   EM_R[40] = 92.6;
   EM_R[41] = 92.6;
   EM_R[42] = 92.6;
   EM_R[43] = 92.6;
   EM_R[44] = 92.6;
   EM_R[45] = 92.6;
   EM_R[46] = 92.6;
   EM_R[47] = 92.6;
   EM_R[48] = 90.6;  //guesstimate
   EM_R[49] = 90.2;  //guesstimate
   EM_R[50] = 92.6;  //guesstimate
   EM_R[51] = 88.0;
   EM_R[52] = 80.4;
   EM_R[53] = 72.1;
   EM_R[54] = 64.8;
   EM_R[55] = 58.3;
   EM_R[56] = 52.5;
   EM_R[57] = 47.4;
   EM_R[58] = 42.7;
   EM_R[59] = 38.6;
   EM_R[60] = 34.8;
   EM_R[61] = 31.4;
   EM_R[62] = 28.4;
   EM_R[63] = 25.7;
   EM_R[64] = 23.2;
   EM_R[65] = 21.0;
   EM_R[66] = 19.0;
   EM_R[67] = 17.2;
   EM_R[68] = 15.5;
   EM_R[69] = 13.2;
   EM_R[70] = 10.4;
   EM_R[71] = 7.6;
   EM_R[72] = 7.2;  //guesstimate
   EM_R[73] = 7.2;  //guesstimate
//
   HAD_R[0] = 5.8;   //guesstimate (These are the average heights for the
   HAD_R[1] = 7.8;   	     //  hadronic part of the calorimeter.)
   HAD_R[2] = 11.3;
   HAD_R[3] = 15.7;
   HAD_R[4] = 20.1;
   HAD_R[5] = 23.5;
   HAD_R[6] = 26.0;
   HAD_R[7] = 28.7;
   HAD_R[8] = 31.8;
   HAD_R[9] = 35.1;
   HAD_R[10] = 38.8;
   HAD_R[11] = 43.0;
   HAD_R[12] = 47.6;
   HAD_R[13] = 52.6;
   HAD_R[14] = 58.3;
   HAD_R[15] = 61.7;
   HAD_R[16] = 67.5;
   HAD_R[17] = 78.4;
   HAD_R[18] = 85.3;
   HAD_R[19] = 93.6;
   HAD_R[20] = 101.8;
   HAD_R[21] = 115.9;
   HAD_R[22] = 123.7;
   HAD_R[23] = 153.4;
   HAD_R[24] = 144.9;  //guesstimate
   HAD_R[25] = 156.1;
   HAD_R[26] = 163.8;
   HAD_R[27] = 155.1;
   HAD_R[28] = 163.6;
   HAD_R[29] = 167.2;
   HAD_R[30] = 158.6; //guesstimate
   HAD_R[31] = 158.6;
   HAD_R[32] = 158.6;
   HAD_R[33] = 158.6;
   HAD_R[34] = 158.6;
   HAD_R[35] = 158.6;
   HAD_R[36] = 158.6;
   HAD_R[37] = 158.6;
   HAD_R[38] = 158.6;
   HAD_R[39] = 158.6;
   HAD_R[40] = 158.6;
   HAD_R[41] = 158.6;
   HAD_R[42] = 158.6;
   HAD_R[43] = 158.6; //guesstimate
   HAD_R[44] = 167.2;
   HAD_R[45] = 163.6;
   HAD_R[46] = 155.1;
   HAD_R[47] = 163.8;
   HAD_R[48] = 156.1;
   HAD_R[49] = 144.9;  //guesstimate
   HAD_R[50] = 153.4;
   HAD_R[51] = 123.7;
   HAD_R[52] = 115.9;
   HAD_R[53] = 101.8;
   HAD_R[54] = 93.6;
   HAD_R[55] = 85.3;
   HAD_R[56] = 78.4;
   HAD_R[57] = 67.5;
   HAD_R[58] = 61.7;
   HAD_R[59] = 58.3;
   HAD_R[60] = 52.6;
   HAD_R[61] = 47.6;
   HAD_R[62] = 43.0;
   HAD_R[63] = 38.8;
   HAD_R[64] = 35.1;
   HAD_R[65] = 31.8;
   HAD_R[66] = 28.7;
   HAD_R[67] = 26.0;
   HAD_R[68] = 23.5;
   HAD_R[69] = 20.1;
   HAD_R[70] = 15.7;
   HAD_R[71] = 11.3;
   HAD_R[72] = 7.8;
   HAD_R[73] = 5.8;  //guesstimate
//

   for (Int_t jentry=0; jentry<nentries;jentry++) {

      Int_t ientry = LoadTree(jentry); //in case of a TChain, ientry is the 							//entry number in the current file
      nb = fChain->GetEntry(jentry);   nbytes += nb;
      TOTRETx = 0;				
      TOTRETy = 0;				
      TOTEMRETx = 0;				
      TOTEMRETy = 0;				
      TOTHDRETx = 0;				
      TOTHDRETy = 0;				
      TOTRET = 0;				
      NEWVERTx = 0;
      NEWVERTy = 0;
      NEWVERTEMx = 0;
      NEWVERTEMy = 0;
      NEWVERTHDx = 0;
      NEWVERTHDy = 0;
      NEWTET = 0;
      THETPRIME = 0;
      THETNOT = 0;
      HADTHETPRIME = 0;
      HADTHETNOT = 0;
      AVE_EMR = 0;
      AVE_EMX = 0;
      AVE_HADR = 0;
      AVE_HADX = 0;
      NRET_EMX = 0;
      NRET_EMY = 0;
      NRET_HDX = 0;
      NRET_HDY = 0;
      COMBVERTx =0;
      COMBVERTy = 0;
      COMBVERTEMx =0;
      COMBVERTEMy = 0;
      COMBVERTHDx =0;
      COMBVERTHDy = 0;
      COMBTET = 0;
      VERTx =0;
      VERTy = 0;
      VERTEMx =0;
      VERTEMy = 0;
      VERTHDx =0;
      VERTHDy = 0;
      VERTET = 0;

// 
      VTX_Z = MET_Zvtx;  //change this line to include the new vertex
//
//      cout << " Please input the z position of the vertex you wish
//                to use, in cm.: " << endl;
//      cin >> VTX_Z;
//
     
      for (Int_t i = 0; i < 74; i++) {
        eta = 0;
        X_EM[i] = 0;
        X_HAD[i] = 0;
        THETPRIME = 0;
        THETNOT = 0;
        HADTHETPRIME = 0;
        HADTHETNOT = 0;
        TOTEMRETx += N_NRETEMx[i];		        
        TOTEMRETy += N_NRETEMy[i];		        
        TOTHDRETx += N_NRETHDx[i];
        TOTHDRETy += N_NRETHDy[i];
//
	if (i==0) eta = -4.5;
	if (i==1) eta = -4.25;
	if (i==2) eta = -3.9;
	if (i==3) eta = -3.55;
	if (i==4) eta = -3.3;
	if (i==73) eta = -4.5;
	if (i==72) eta = -4.25;
	if (i==71) eta = -3.9;
	if (i==70) eta = -3.55;
	if (i==69) eta = -3.3;

	if (i>4 && i<37) eta = i/10. - 3.65;
	if (i<69 && i>36) eta = i/10. - 3.55;
//
	X_EM[i] = EM_R[i]/(tan(2. * atan(exp((-1.)*eta)) ));
//
	THETPRIME = atan2( EM_R[i] , (X_EM[i] - VTX_Z) );
	 	// VTX_Z is the event vertex supplied by the user
	if (X_EM[i] != 0) THETNOT = atan2( EM_R[i] , X_EM[i] );

//
	X_HAD[i] = HAD_R[i]/(tan(2. * atan(exp((-1.)*eta)) ));
//
	HADTHETPRIME = atan2( HAD_R[i] , (X_HAD[i] - VTX_Z) );
	 	// VTX_Z is the event vertex supplied by the user
	if (X_HAD[i] != 0) HADTHETNOT = atan2( HAD_R[i] , X_HAD[i] );

	NEWVERTEMx += ( N_NRETEMx[i] * (sin(THETPRIME)/sin(THETNOT)) );
	NEWVERTEMy += ( N_NRETEMy[i] * (sin(THETPRIME)/sin(THETNOT)) );
	NEWVERTHDx += ( N_NRETHDx[i] * (sin(HADTHETPRIME)/sin(HADTHETNOT)) );
	NEWVERTHDy += ( N_NRETHDy[i] * (sin(HADTHETPRIME)/sin(HADTHETNOT)) );

      }  // if (i<74) 				

      TOTRETx = TOTEMRETx + TOTHDRETx;
      TOTRETy = TOTEMRETy + TOTHDRETy;

      TOTRET = sqrt(TOTRETx*TOTRETx + TOTRETy*TOTRETy);	

      NEWVERTx = NEWVERTEMx + NEWVERTHDx;
      NEWVERTy = NEWVERTEMy + NEWVERTHDy;

      NEWTET = sqrt(NEWVERTx*NEWVERTx + NEWVERTy*NEWVERTy);

      for (Int_t j = 0; j < 74; j++) {
        eta = 0;
        X_EM[j] = 0;
        X_HAD[j] = 0;
        X_EM[j+1] = 0;
        X_HAD[j+1] = 0;
	AVE_EMR = 0;
        AVE_EMX = 0;
	AVE_HADR = 0;
        AVE_HADX = 0;
        NRET_EMX = 0;
        NRET_EMY = 0;
        NRET_HDX = 0;
        NRET_HDY = 0;
        THETPRIME = 0;
        THETNOT = 0;
        HADTHETPRIME = 0;
        HADTHETNOT = 0;
//
	if (j==0) eta = -4.5;
	if (j==1) eta = -4.25;
	if (j==2) eta = -3.9;
	if (j==3) eta = -3.55;
	if (j==4) eta = -3.3;
	if (j==73) eta = -4.5;
	if (j==72) eta = -4.25;
	if (j==71) eta = -3.9;
	if (j==70) eta = -3.55;
	if (j==69) eta = -3.3;
	if (j>4 && j<37) eta = j/10. - 3.65;
	if (j<69 && j>36) eta = j/10. - 3.55;

//
	X_EM[j] = EM_R[j]/(tan(2. * atan(exp((-1.)*eta)) ));
	X_HAD[j] = HAD_R[j]/(tan(2. * atan(exp((-1.)*eta)) ));
        if (j==36) {
	  X_EM[j+1] = X_EM[j];
	  X_HAD[j+1] = X_HAD[j];
	}
	else {
	  X_EM[j+1] = EM_R[j+1]/(tan(2. * atan(exp( (-1.)*(eta+.1) )) ));
	  X_HAD[j+1] = HAD_R[j+1]/(tan(2. * atan(exp( (-1.)*(eta+.1) )) ));
	}
//
        NRET_EMX = N_NRETEMx[j] + N_NRETEMx[j+1];
        NRET_EMY = N_NRETEMy[j] + N_NRETEMy[j+1];
        NRET_HDX = N_NRETHDx[j] + N_NRETHDx[j+1];
        NRET_HDY = N_NRETHDy[j] + N_NRETHDy[j+1];
//
        AVE_EMR = ( EM_R[j] + EM_R[j+1] )/2;
        AVE_EMX = ( X_EM[j] + X_EM[j+1] )/2;
        AVE_HADR = ( HAD_R[j] + HAD_R[j+1] )/2;
        AVE_HADX = ( X_HAD[j] + X_HAD[j+1] )/2;
//
	THETPRIME = atan2( AVE_EMR , (AVE_EMX - VTX_Z) );
	 	// VTX_Z is the event vertex supplied by the user
	if (AVE_EMX != 0) THETNOT = atan2( AVE_EMR , AVE_EMX );
//
	HADTHETPRIME = atan2( AVE_HADR , (AVE_HADX - VTX_Z) );
	 	// VTX_Z is the event vertex supplied by the user
	if (AVE_HADX != 0) HADTHETNOT = atan2( AVE_HADR , AVE_HADX );
//
	if (sin(THETNOT) != 0)
	  COMBVERTEMx += ( NRET_EMX * (sin(THETPRIME)/sin(THETNOT)) );
	if (sin(THETNOT) != 0)
	  COMBVERTEMy += ( NRET_EMY * (sin(THETPRIME)/sin(THETNOT)) );
	if (sin(HADTHETNOT) != 0)
	  COMBVERTHDx += ( NRET_HDX * (sin(HADTHETPRIME)/sin(HADTHETNOT)) );
	if (sin(HADTHETNOT) != 0)
	  COMBVERTHDy += ( NRET_HDY * (sin(HADTHETPRIME)/sin(HADTHETNOT)) );
//
        Hist44->Fill( NRET_EMX );
        Hist45->Fill( NRET_HDX );
        Hist46->Fill( NRET_EMY );
        Hist47->Fill( NRET_HDY );
        Hist50->Fill( THETPRIME );
        Hist51->Fill( THETNOT );
        Hist52->Fill( THETPRIME - THETNOT );
        Hist53->Fill( HADTHETPRIME );
        Hist54->Fill( HADTHETNOT );
        Hist55->Fill( HADTHETPRIME - HADTHETNOT );
        Hist60->Fill( N_NRETEMx[j] );
        Hist61->Fill( N_NRETEMy[j] );
        Hist62->Fill( N_NRETHDx[j] );
        Hist63->Fill( N_NRETHDy[j] );
                                                       
        j++;  //go through j's 2 at a time
//
      }  // if (j<74) 

     for (Int_t k = 0; k < 74; k++) {
        eta = 0;
        X_EM[k] = 0;
        X_HAD[k] = 0;
        X_EM[k+1] = 0;
        X_HAD[k+1] = 0;
//        X_EM[k+2] = 0;
//        X_HAD[k+2] = 0;
	AVE_EMR = 0;
        AVE_EMX = 0;
	AVE_HADR = 0;
        AVE_HADX = 0;
        NRET_EMX = 0;
        NRET_EMY = 0;
        NRET_HDX = 0;
        NRET_HDY = 0;
        THETPRIME = 0;
        THETNOT = 0;
        HADTHETPRIME = 0;
        HADTHETNOT = 0;
//
	if (k==0) eta = -4.5;
	if (k==1) eta = -4.25;
	if (k==2) eta = -3.9;
	if (k==3) eta = -3.55;
	if (k==4) eta = -3.3;
	if (k==73) eta = -4.5;
	if (k==72) eta = -4.25;
	if (k==71) eta = -3.9;
	if (k==70) eta = -3.55;
	if (k==69) eta = -3.3;
	if (k>4 && k<37) eta = k/10. - 3.65;
	if (k<69 && k>36) eta = k/10. - 3.55;

//
	X_EM[k] = EM_R[k]/(tan(2. * atan(exp((-1.)*eta)) ));
	X_HAD[k] = HAD_R[k]/(tan(2. * atan(exp((-1.)*eta)) ));
        if (k==36) {
	  X_EM[k+1] = EM_R[k+1]/(tan(2. * atan(exp( (-1.)*(eta+.2) )) ));
	  X_HAD[k+1] = HAD_R[k+1]/(tan(2. * atan(exp( (-1.)*(eta+.2) )) ));
	}
	else {
	  X_EM[k+1] = EM_R[k+1]/(tan(2. * atan(exp( (-1.)*(eta+.1) )) ));
	  X_HAD[k+1] = HAD_R[k+1]/(tan(2. * atan(exp( (-1.)*(eta+.1) )) ));
	}
//
	if (k<72) {
        NRET_EMX = N_NRETEMx[k] + N_NRETEMx[k+1] + N_NRETEMx[k+2];
        NRET_EMY = N_NRETEMy[k] + N_NRETEMy[k+1] + N_NRETEMy[k+2];
        NRET_HDX = N_NRETHDx[k] + N_NRETHDx[k+1] + N_NRETHDx[k+2];
        NRET_HDY = N_NRETHDy[k] + N_NRETHDy[k+1] + N_NRETHDy[k+2];
	}
        else {
          NRET_EMX = N_NRETEMx[k] + N_NRETEMx[k+1];
          NRET_EMY = N_NRETEMy[k] + N_NRETEMy[k+1];
          NRET_HDX = N_NRETHDx[k] + N_NRETHDx[k+1];
          NRET_HDY = N_NRETHDy[k] + N_NRETHDy[k+1];
	}

//
//        AVE_EMR = ( EM_R[k] + EM_R[k+1] + EM_R[k+2] )/3;
//        AVE_EMX = ( X_EM[k] + X_EM[k+1] + X_EM[k+2] )/3;
//        AVE_HADR = ( HAD_R[k] + HAD_R[k+1] + HAD_R[k+2] )/3;
//        AVE_HADX = ( X_HAD[k] + X_HAD[k+1] + X_HAD[k+2] )/3;
        AVE_EMR = EM_R[k+1] ;
        AVE_EMX = X_EM[k+1] ;
        AVE_HADR = HAD_R[k+1] ;
        AVE_HADX = X_HAD[k+1] ;
//

	THETPRIME = atan2( AVE_EMR , (AVE_EMX - VTX_Z) );
	 	// VTX_Z is the event vertex supplied by the user
	if (AVE_EMX != 0) THETNOT = atan2( AVE_EMR , AVE_EMX );
//
	HADTHETPRIME = atan2( AVE_HADR , (AVE_HADX - VTX_Z) );
	 	// VTX_Z is the event vertex supplied by the user
	if (AVE_HADX != 0) HADTHETNOT = atan2( AVE_HADR , AVE_HADX );
//
	if (sin(THETNOT) != 0)
	  VERTEMx += ( NRET_EMX * (sin(THETPRIME)/sin(THETNOT)) );
	if (sin(THETNOT) != 0)
	  VERTEMy += ( NRET_EMY * (sin(THETPRIME)/sin(THETNOT)) );
	if (sin(HADTHETNOT) != 0)
	  VERTHDx += ( NRET_HDX * (sin(HADTHETPRIME)/sin(HADTHETNOT)) );
	if (sin(HADTHETNOT) != 0)
	  VERTHDy += ( NRET_HDY * (sin(HADTHETPRIME)/sin(HADTHETNOT)) );
//
        k++;  //go through k's 2 at a time
        k++;  //go through k's 3 at a time
//
      }  // if (k<74) 


      COMBVERTx = COMBVERTEMx + COMBVERTHDx;
      COMBVERTy = COMBVERTEMy + COMBVERTHDy;

      COMBTET = sqrt(COMBVERTx*COMBVERTx + COMBVERTy*COMBVERTy);

      VERTx = VERTEMx + VERTHDx;
      VERTy = VERTEMy + VERTHDy;

      VERTET = sqrt(VERTx*VERTx + VERTy*VERTy);

      Hist1->Fill(MET_MEx);
      Hist2->Fill(MET_MEy);
      Hist3->Fill(TOTRETx);
      Hist4->Fill(TOTRETy);
      Hist5->Fill(MET_MEx - TOTRETx);
      Hist6->Fill(MET_MEy - TOTRETy);
      Hist7->Fill(MET_MET);
      Hist8->Fill(TOTRET);
      if (MET_MET != 0)
         Hist9->Fill( 1 - TOTRET/MET_MET );
      Hist10->Fill( MET_MET - TOTRET );
      Hist11->Fill(NEWTET);
      if (MET_MET != 0)
         Hist12->Fill( 1 - NEWTET/MET_MET );
//      Hist12->SetLineColor(51);
      Hist12->SetLineColor(77);
      Hist13->Fill( MET_MET - NEWTET );
      Hist14->Fill( MET_Zvtx );
      Hist21->Fill( MET_Zvtx );
      Hist15->Fill( VTX_Z );
      Hist34->Fill( NEWTET - TOTRET );
      if (NEWTET != 0)
         Hist38->Fill( 1 - TOTRET/NEWTET );
      Hist35->Fill( COMBTET );
      Hist36->Fill( MET_MET - COMBTET );
      if (MET_MET != 0)
         Hist39->Fill( 1 - COMBTET/MET_MET );
      Hist37->Fill( NEWTET - COMBTET );
      if (NEWTET != 0)
         Hist40->Fill( 1 - COMBTET/NEWTET );
      Hist43->Fill( NEWTET - MET_MET );
      Hist41->Fill(MET_MET - TOTRET,MET_Zvtx);
      if (MET_MET != 0)
         Hist42->Fill(1 - TOTRET/MET_MET,MET_Zvtx);
      Hist48->Fill( COMBVERTx );
      Hist49->Fill( COMBVERTy );
      Hist56->Fill( COMBVERTEMx );
      Hist57->Fill( COMBVERTEMy );
      Hist58->Fill( COMBVERTHDx );
      Hist59->Fill( COMBVERTHDy );
      Hist64->Fill( VERTET );
      Hist65->Fill( MET_MET - VERTET );
      if (MET_MET != 0)
         Hist66->Fill( 1 - VERTET/MET_MET );
      Hist67->Fill( NEWTET - VERTET );
      if (NEWTET != 0)
         Hist68->Fill( 1 - VERTET/NEWTET );
      Hist69->Fill( COMBTET - VERTET );
      if (COMBTET != 0)
         Hist70->Fill( 1 - VERTET/COMBTET );

      // if (Cut(ientry) < 0) continue;
    } // if (jentry<nentries)

   Hist12->Draw();
   Hist9->Draw("same");

}  // Loop (Main program)