#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 in 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 *Hist16 = new TH1F("Hist16","THETPRIME",100,-6,6);
   TH1F *Hist17 = new TH1F("Hist17","THETNOT",100,-6,6);
   TH1F *Hist18 = new TH1F("Hist18","THETPRIME-THENOT",100,-6,6);
   TH1F *Hist19 = new TH1F("Hist19","X[i]",100,-200,200);
   TH1F *Hist20 = new TH1F("Hist20","R[i]",100,0,100);
   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);


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

   Int_t nbytes = 0, nb = 0;

   Double_t TOTRETx = 0;			
   Double_t TOTRETy = 0;			
   Double_t TOTRET = 0;
   Double_t NEWVERTx = 0;
   Double_t NEWVERTy = 0; 
   Double_t NEWTET = 0;
   Double_t THETPRIME = 0;
   Double_t THETNOT = 0;
   Double_t VTX_Z = 0;
   Double_t eta = 0;
   Double_t COMBVERTx =0;
   Double_t COMBVERTy = 0;
   Double_t COMBTET = 0;
   Double_t NRET_X = 0;
   Double_t NRET_Y = 0;

   Double_t X[74];
   Double_t R[74];   
// 
   R[0] = 6.8;   //guesstimate (These are the average heights for the 
   R[1] = 8.2;   //guesstimate    calorimeter, taken as a whole.)
   R[2] = 10.6;
   R[3] = 14.6;
   R[4] = 18.6;
   R[5] = 21.8;
   R[6] = 24.0;
   R[7] = 26.6;
   R[8] = 29.4;
   R[9] = 32.6;
   R[10] = 36.0;
   R[11] = 39.8;
   R[12] = 44.0;
   R[13] = 48.8;
   R[14] = 54.0;
   R[15] = 57.0;
   R[16] = 62.2;
   R[17] = 72.6;
   R[18] = 78.8;
   R[19] = 86.5;
   R[20] = 96.2;
   R[21] = 107.0;
   R[22] = 113.0;
   R[23] = 117.7;  //guesstimate
   R[24] = 122.0;  //guesstimate
   R[25] = 126.2;
   R[26] = 129.8;
   R[27] = 137.8;
   R[28] = 146.3;
   R[29] = 149.9;
   R[30] = 141.1; //guesstimate
   R[31] = 141.1;
   R[32] = 141.1;
   R[33] = 141.1;
   R[34] = 141.1;
   R[35] = 141.1;
   R[36] = 141.1;
   R[37] = 141.1;
   R[38] = 141.1;
   R[39] = 141.1;
   R[40] = 141.1;
   R[41] = 141.1;
   R[42] = 141.1;
   R[43] = 141.1; //guesstimate
   R[44] = 149.9;
   R[45] = 146.3;
   R[46] = 137.8;
   R[47] = 129.8;
   R[48] = 126.2;
   R[49] = 122.0;  //guesstimate
   R[50] = 117.7;  //guesstimate
   R[51] = 113.0;
   R[52] = 107.0;
   R[53] = 96.2;
   R[54] = 86.5;
   R[55] = 78.8;
   R[56] = 72.6;
   R[57] = 62.2;
   R[58] = 57.0;
   R[59] = 54.0;
   R[60] = 48.8;
   R[61] = 44.0;
   R[62] = 39.8;
   R[63] = 36.0;
   R[64] = 32.6;
   R[65] = 29.4;
   R[66] = 26.6;
   R[67] = 24.0;
   R[68] = 21.8;
   R[69] = 18.6;
   R[70] = 14.6;
   R[71] = 10.6;
   R[72] = 8.2;  //guesstimate
   R[73] = 6.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;				
      TOTRET = 0;				
      NEWVERTx = 0;
      NEWVERTy = 0;
      NEWTET = 0;
      THETPRIME = 0;
      THETNOT = 0;
      COMBVERTx =0;
      COMBVERTy = 0;
      COMBTET = 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[i] = 0;
        THETPRIME = 0;
        THETNOT = 0;
        TOTRETx += N_NRETEMx[i] + N_NRETHDx[i];		        
        TOTRETy += N_NRETEMy[i] + 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[i] = R[i]/(tan(2. * atan(exp((-1.)*eta)) ));
//
	THETPRIME = atan2( R[i] , (X[i] - VTX_Z) );
	 	// VTX_Z is the event vertex supplied by the user
	if (X[i] != 0) THETNOT = atan2( R[i] , X[i] );


	NEWVERTx += N_NRETEMx[i] * (sin(THETPRIME)/sin(THETNOT));
	NEWVERTx += N_NRETHDx[i] * (sin(THETPRIME)/sin(THETNOT));
	NEWVERTy += N_NRETEMy[i] * (sin(THETPRIME)/sin(THETNOT));
	NEWVERTy += N_NRETHDy[i] * (sin(THETPRIME)/sin(THETNOT));

        Hist16->Fill( THETPRIME );
        Hist17->Fill( THETNOT );
        Hist18->Fill( THETPRIME - THETNOT );
        Hist19->Fill( X[i] );
        Hist20->Fill( R[i] );

      }  // if (i<74) 				

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

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

      for (Int_t j = 0; j < 74; j++) {
        eta = 0;
        X[j] = 0;
	Double_t AVE_R = 0;
        Double_t AVE_X = 0;
        NRET_X = 0;
        NRET_Y = 0;
        THETPRIME = 0;
        THETNOT = 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[j] = R[j]/(tan(2. * atan(exp((-1.)*eta)) ));
//
        NRET_X = N_NRETEMx[j] + N_NRETEMx[j+1] + N_NRETHDx[j] + N_NRETHDx[j+1];
        NRET_Y = N_NRETEMy[j] + N_NRETEMy[j+1] + N_NRETHDy[j] + N_NRETHDy[j+1];
//
        AVE_R = (R[j] + R[j+1])/2;
        AVE_X = (X[j] + X[j+1])/2;
//
	THETPRIME = atan2( AVE_R , (AVE_X - VTX_Z) );
	 	// VTX_Z is the event vertex supplied by the user
	if (AVE_X != 0) THETNOT = atan2( AVE_R , AVE_X );
//
	COMBVERTx += NRET_X * (sin(THETPRIME)/sin(THETNOT));
	COMBVERTy += NRET_Y * (sin(THETPRIME)/sin(THETNOT));

        j++;  //go through j's 2 at a time
//
      }


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

      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 );
      Hist41->Fill(MET_MET - TOTRET,MET_Zvtx);
      if (MET_MET != 0)
         Hist42->Fill(1 - TOTRET/MET_MET,MET_Zvtx);

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

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

}  // Loop (Main program)