Index: trunk/src/SmearUtil.cc
===================================================================
--- trunk/src/SmearUtil.cc	(revision 69)
+++ trunk/src/SmearUtil.cc	(revision 71)
@@ -52,7 +52,7 @@
 ELG_Ncen =        0.25 ;                 // N term for central ECAL
 ELG_Ccen =        0.0055 ;                // C term for central ECAL
-ELG_Sfwd =        1.5;                    // S term for forward ECAL
-ELG_Nfwd =        0.0;                    // N term for central ECAL
-ELG_Cfwd =        0.06;                   // C term for forward ECAL
+ELG_Cfwd =        0.107 ;                 // S term for forward ECAL
+ELG_Sfwd =        2.084 ;                 // C term for forward ECAL
+ELG_Nfwd =        0.0   ;                 // N term for central ECAL
 
 HAD_Shcal =       1.5 ;                    // S term for central HCAL // hadronic calorimeter
@@ -104,4 +104,20 @@
 M_MAXITERATIONS     =   100;
 
+// Define Calorimeter Towers
+NTOWERS	= 40;
+
+const float tower_eta_edges[41] = { 
+    0., 0.087, 0.174, 0.261, 0.348, 0.435, 0.522, 0.609, 0.696, 0.783, 0.870, 0.957, 1.044, 1.131, 1.218, 1.305, 1.392, 1.479, 1.566,
+    1.653, 1.740, 1.830, 1.930, 2.043, 2.172, 2.322, 2.500, 2.650, 2.868, 2.950, 3.125, 3.300, 3.475, 3.650, 3.825, 4.000, 4.175,
+    4.350, 4.525,  4.700, 5.000}; // temporary object
+TOWER_ETA_EDGES = new float[NTOWERS+1];
+for(unsigned int i=0; i<NTOWERS+1; i++) TOWER_ETA_EDGES[i] = tower_eta_edges[i];
+
+const float tower_dphi[40] = { 
+     5, 5, 5, 5, 5,  5, 5, 5, 5, 5,  5, 5, 5, 5, 5,  5, 5, 5, 5, 10,
+    10,10,10,10,10, 10,10,10,10,10, 10,10,10,10,10, 10,10,10,20, 20 }; // temporary object
+TOWER_DPHI = new float[NTOWERS];
+for(unsigned int i=0; i<NTOWERS; i++) TOWER_DPHI[i] = tower_dphi[i];
+
 }
 
@@ -122,5 +138,5 @@
     curstring.str(temp_string);
     string varname;
-    float value;
+    float value; int ivalue;
     
     if(strstr(temp_string.c_str(),"#")) { }
@@ -129,6 +145,6 @@
     else if(strstr(temp_string.c_str(),"MAX_CALO_FWD")){curstring >> varname >> value; MAX_CALO_FWD = value;}
     else if(strstr(temp_string.c_str(),"MAX_MU")){curstring >> varname >> value; MAX_MU = value;}
-    else if(strstr(temp_string.c_str(),"TRACKING_RADIUS")){curstring >> varname >> value; TRACKING_RADIUS = (int)value;}
-    else if(strstr(temp_string.c_str(),"TRACKING_LENGTH")){curstring >> varname >> value; TRACKING_LENGTH = (int)value;}
+    else if(strstr(temp_string.c_str(),"TRACKING_RADIUS")){curstring >> varname >> ivalue; TRACKING_RADIUS = ivalue;}
+    else if(strstr(temp_string.c_str(),"TRACKING_LENGTH")){curstring >> varname >> ivalue; TRACKING_LENGTH = ivalue;}
     else if(strstr(temp_string.c_str(),"BFIELD_X")){curstring >> varname >> value; BFIELD_X = value;}
     else if(strstr(temp_string.c_str(),"BFIELD_Y")){curstring >> varname >> value; BFIELD_Y = value;}
@@ -151,14 +167,20 @@
     else if(strstr(temp_string.c_str(),"PT_TRACK_TAU")){curstring >> varname >> value; PT_TRACK_TAU = value;}
     else if(strstr(temp_string.c_str(),"PT_TRACKS_MIN")){curstring >> varname >> value; PT_TRACKS_MIN = value;}
-    else if(strstr(temp_string.c_str(),"TAGGING_B")){curstring >> varname >> value; TAGGING_B = (int)value;}
-    else if(strstr(temp_string.c_str(),"MISTAGGING_C")){curstring >> varname >> value; MISTAGGING_C = (int)value;}
-    else if(strstr(temp_string.c_str(),"MISTAGGING_L")){curstring >> varname >> value; MISTAGGING_L = (int)value;}
+    else if(strstr(temp_string.c_str(),"TAGGING_B")){curstring >> varname >> ivalue; TAGGING_B = ivalue;}
+    else if(strstr(temp_string.c_str(),"MISTAGGING_C")){curstring >> varname >> ivalue; MISTAGGING_C = ivalue;}
+    else if(strstr(temp_string.c_str(),"MISTAGGING_L")){curstring >> varname >> ivalue; MISTAGGING_L = ivalue;}
     else if(strstr(temp_string.c_str(),"CONERADIUS")){curstring >> varname >> value; CONERADIUS = value;}
-    else if(strstr(temp_string.c_str(),"JETALGO")){curstring >> varname >> value; JETALGO = (int)value;}
-    else if(strstr(temp_string.c_str(),"TRACKING_EFF")){curstring >> varname >> value; TRACKING_EFF = (int)value;}
+    else if(strstr(temp_string.c_str(),"JETALGO")){curstring >> varname >> ivalue; JETALGO = ivalue;}
+    else if(strstr(temp_string.c_str(),"TRACKING_EFF")){curstring >> varname >> ivalue; TRACKING_EFF = ivalue;}
     else if(strstr(temp_string.c_str(),"ELEC_pt")){curstring >> varname >> value; ELEC_pt = value;}
     else if(strstr(temp_string.c_str(),"MUON_pt")){curstring >> varname >> value; MUON_pt = value;}
     else if(strstr(temp_string.c_str(),"JET_pt")){curstring >> varname >> value; JET_pt = value;}
     else if(strstr(temp_string.c_str(),"TAUJET_pt")){curstring >> varname >> value; TAUJET_pt = value;}
+    else if(strstr(temp_string.c_str(),"NTOWERS")){curstring >> varname >> ivalue; NTOWERS = ivalue;}
+    else if(strstr(temp_string.c_str(),"TOWER_ETA_EDGES")){
+	curstring >> varname;	for(unsigned int i=0; i<NTOWERS+1; i++) {curstring >> value; TOWER_ETA_EDGES[i] = value;} }
+    else if(strstr(temp_string.c_str(),"TOWER_DPHI")){
+	curstring >> varname;	for(unsigned int i=0; i<NTOWERS; i++) {curstring >> value; TOWER_DPHI[i] = value;} }
+   
 
   }
@@ -182,9 +204,9 @@
  PROTOJET_PTMIN = 0.0; 
 
- DOTRIGGER=1;
 
 }
 
 void RESOLution::Logfile(string LogName) { 
+//void RESOLution::Logfile(string outputfilename) { 
 
   ofstream f_out(LogName.c_str());
@@ -266,4 +288,30 @@
         << left << setw(5) <<BFIELD_Z <<""<< right << setw(10)<<"*"<<"\n";
   f_out<<"*                                                                    *"<<"\n";
+
+
+  f_out<<"*                                                                    *"<<"\n";
+  f_out<<"#********************                                                *"<<"\n";
+  f_out<<"# Calorimetric Towers                                                *"<<"\n";
+  f_out<<"#********************                                                *"<<"\n";
+  f_out << left << setw(55) <<"* Number of calorimetric towers in eta, for eta>0: "<<""
+        << left << setw(5) << NTOWERS <<""<< right << setw(10)<<"*"<<"\n";
+  f_out << left << setw(55) <<"* Tower edges in eta, for eta>0: "<<"" << right << setw(15)<<"*"<<"\n";
+	  f_out << "*   ";
+	  for (unsigned int i=0; i<NTOWERS+1; i++) {
+		f_out << left << setw(7) << TOWER_ETA_EDGES[i];
+		if(!( (i+1) %9 )) f_out << right << setw(3) << "*" << "\n" << "*   ";
+	  }
+	  for (unsigned int i=(NTOWERS+1)%9; i<9; i++) f_out << left << setw(7) << "";
+	  f_out << right << setw(3)<<"*"<<"\n";
+  f_out << left << setw(55) <<"* Tower sizes in phi, for eta>0 [degree]:"<<"" << right << setw(15)<<"*"<<"\n";
+	  f_out << "*   ";
+	  for (unsigned int i=0; i<NTOWERS; i++) {
+		f_out << left << setw(7) << TOWER_DPHI[i];
+		if(!( (i+1) %9 )) f_out << right << setw(3) << "*" << "\n" << "*   ";
+	  }
+	  for (unsigned int i=(NTOWERS)%9; i<9; i++) f_out << left << setw(7) << "";
+	  f_out << right << setw(3)<<"*"<<"\n";
+  f_out<<"*                                                                    *"<<"\n";
+
   f_out<<"#************************************                                *"<<"\n";
   f_out<<"# Electromagnetic smearing parameters                                *"<<"\n";
@@ -446,7 +494,6 @@
   float energy = electron.E();  // before smearing
   float energyS = 0.0;		// after smearing  // \sigma/E = C + N/E + S/\sqrt{E}
-  float eta=fabs(electron.Eta()); 
- 
-  if(eta < MAX_TRACKER) { // if the electron is inside the tracker
+  
+  if(fabs(electron.Eta()) < MAX_TRACKER) { // if the electron is inside the tracker
     energyS = gRandom->Gaus(energy, sqrt( 
 					 pow(ELG_Ncen,2) + 
@@ -454,5 +501,5 @@
 					 pow(ELG_Scen*sqrt(energy),2) ));
   } 
-  if(eta > MAX_TRACKER && eta < MAX_CALO_FWD){
+  if(fabs(electron.Eta()) > MAX_TRACKER && fabs(electron.Eta()) < MAX_CALO_FWD){
     energyS = gRandom->Gaus(energy, sqrt(
 					 pow(ELG_Nfwd,2) +
@@ -609,4 +656,31 @@
 
 
+  //********** returns a segmented value for eta and phi, for calo towers *****
+void RESOLution::BinEtaPhi(const float phi, const float eta, float& iPhi, float& iEta){
+   iEta = -100; 
+   int index=-100;
+   for (unsigned int i=1; i< NTOWERS+1; i++) {
+	if(fabs(eta)>TOWER_ETA_EDGES[i-1] && fabs(eta)<TOWER_ETA_EDGES[i]) {
+		iEta = (eta>0) ? TOWER_ETA_EDGES[i-1] : -TOWER_ETA_EDGES[i];
+		index = i-1;
+		//cout << setw(15) << left << eta << "\t" << iEta << endl;
+		break;
+	}
+   }
+   if(index==-100) return; 
+   iPhi = -100;
+   float dphi = TOWER_DPHI[index]*PI/180.; 
+   for (unsigned int i=1; i < 360/TOWER_DPHI[index]; i++ ) {
+	float low = -PI+(i-1)*dphi;
+	float high= low+dphi;
+	if(phi > low && phi < high ){
+		iPhi = low; 
+		break;
+	} 
+   }
+   if (phi > PI-dphi) iPhi = PI-dphi;
+}
+
+
 //**************************** Returns the delta Phi ****************************
 float DeltaPhi(const float phi1, const float phi2) {
