diff --git a/macros/CCEScan/pulseFit.C b/macros/CCEScan/pulseFit.C
index ee91c6e..16551bf 100755
--- a/macros/CCEScan/pulseFit.C
+++ b/macros/CCEScan/pulseFit.C
@@ -267,6 +267,7 @@
   
   double chi2ndf = fu_pulse->GetChisquare()/fu_pulse->GetNDF();
   Info("pulseFit","#chi^{2}/ndf: %6.2f/%d",fu_pulse->GetChisquare(),fu_pulse->GetNDF());
+
   
   // Make PDG-like correction of the uncertainties
   //Commented out by Michele - 8/31 2017
@@ -281,6 +282,150 @@
   ge_pulse = new TGraphErrors(vr_time_val,vr_adc_val,vr_time_err,vr_adc_err);
   ge_pulse->Fit(fu_pulse,"M0Q","0");
   Info("pulseFit","#chi^{2}/ndf: %6.2f/%d",fu_pulse->GetChisquare(),fu_pulse->GetNDF());
+
+
+
+  //MIchele - STart
+  double low_limit = -50;
+  double high_limit = 50;
+
+  //POLY3
+
+  cout << "================ poly3 ====================" << endl<<endl;
+
+  TF1 *poly3 = new TF1("poly3","[0] + [1]*x +[2]*x**2 + [3]*x**3 ",-40.,40.);
+  poly3->SetLineColor(1);
+  ge_pulse->Fit("poly3","M");
+  Info("pulseFit","#chi^{2}/ndf: %6.2f/%d",poly3->GetChisquare(),poly3->GetNDF());
+
+  double x1_poly3 = poly3->GetX(0,low_limit,0);
+  double x2_poly3 = poly3->GetX(0,0, high_limit);
+  double z1_poly3 = -99;
+  double z2_poly3 = -99;
+  bool b_poly3 = false;
+  double integral_poly3_val = -999;
+  double integral_poly3_err = 0;
+
+
+  if( TMath::Abs(poly3->Eval(x1_poly3)) > 1e-08){cout <<"No zeros on the left, the lower point is y1: "<<poly3->Eval(x1_poly3)<<endl;}
+  else{
+    z1_poly3 = x1_poly3;
+    cout <<"Zero found! z1: "<< z1_poly3<<", y1: "<<poly3->Eval(x1_poly3)<<endl;
+    if( TMath::Abs(poly3->Eval(x2_poly3)) > 1e-08){cout <<"No zero found, the lower point is y2: "<<poly3->Eval(x2_poly3)<<endl;}
+    else{
+      z2_poly3 = x2_poly3;
+      cout << "Zero found! z2: "<<z2_poly3<<", y2: "<<poly3->Eval(x2_poly3)<<endl;  
+      b_poly3 = true;
+    }
+  }  
+
+ 
+  //This has to be done before any other fit is performed
+  if(b_poly3 == true){
+
+    cout << "Poly3 behaves properly in the range ["<<low_limit<<", "<<high_limit<<"]"<<endl;
+    integral_poly3_val = poly3->Integral(z1_poly3,z2_poly3)/20.;
+    integral_poly3_err = poly3->IntegralError(z1_poly3,z2_poly3)/20.;
+    cout << "The integral between the two zeros is "<< integral_poly3_val <<"+/-"<<integral_poly3_err << endl;
+
+  }
+
+
+
+
+
+  //POLY4
+
+  cout << "============ poly4 ==============" << endl<<endl;
+
+  TF1 *poly4 = new TF1("poly4","[0] + [1]*x +[2]*x**2 + [3]*x**3 + [4]*x**4 ",-40.,40.);
+  poly4->SetLineColor(4);
+  ge_pulse->Fit("poly4","M0Q","0");
+  Info("pulseFit","#chi^{2}/ndf: %6.2f/%d",poly4->GetChisquare(),poly4->GetNDF());
+
+
+  double x1_poly4 = poly4->GetX(0,low_limit,0);
+  double x2_poly4 = poly4->GetX(0,0, high_limit);
+  double z1_poly4 = -99;
+  double z2_poly4 = -99;
+  bool b_poly4 = false;
+  double integral_poly4_val = -999;
+  double integral_poly4_err = 0;
+
+
+  if( TMath::Abs(poly4->Eval(x1_poly4)) > 1e-08){cout <<"No zeros on the left, the lower point is y1: "<<poly4->Eval(x1_poly4)<<endl;}
+  else{
+    z1_poly4 = x1_poly4;
+    cout <<"Zero found! z1: "<< z1_poly4<<", y1: "<<poly4->Eval(x1_poly4)<<endl;
+    if( TMath::Abs(poly4->Eval(x2_poly4)) > 1e-08){cout <<"No zero found, the lower point is y2: "<<poly4->Eval(x2_poly4)<<endl;}
+    else{
+      z2_poly4 = x2_poly4;
+      cout << "Zero found! z2: "<<z2_poly4<<", y2: "<<poly4->Eval(x2_poly4)<<endl;  
+      b_poly4 = true;
+    }
+  }  
+
+ 
+  //This has to be done before any other fit is performed
+  if(b_poly4 == true){
+
+    cout << "Poly4 behaves properly in the range ["<<low_limit<<", "<<high_limit<<"]"<<endl;
+    integral_poly4_val = poly4->Integral(z1_poly4,z2_poly4)/20.;
+    integral_poly4_err = poly4->IntegralError(z1_poly4,z2_poly4)/20.;
+    cout << "The integral between the two zeros is "<< integral_poly4_val <<"+/-"<<integral_poly4_err << endl;
+
+  }
+
+
+
+  //POLY5
+  cout << "================  poly5 =================" << endl<<endl;
+
+  TF1 *poly5 = new TF1("poly5","[0] + [1]*x +[2]*x**2 + [3]*x**3 + [4]*x**4 + [5]*x**5",-40.,40.);
+  poly5->SetLineColor(6);
+  ge_pulse->Fit("poly5","M");
+  Info("pulseFit","#chi^{2}/ndf: %6.2f/%d",poly5->GetChisquare(),poly5->GetNDF()); 
+
+
+
+  double x1_poly5 = poly5->GetX(0,low_limit,0);
+  double x2_poly5 = poly5->GetX(0,0, high_limit);
+  double z1_poly5 = -99;
+  double z2_poly5 = -99;
+  bool b_poly5 = false;
+  double integral_poly5_val = -999;
+  double integral_poly5_err = 0;
+
+  if( TMath::Abs(poly5->Eval(x1_poly5)) > 1e-08){cout <<"No zeros on the left, the lower point is y1: "<<poly5->Eval(x1_poly5)<<endl;}
+  else{
+    z1_poly5 = x1_poly5;
+    cout <<"Zero found! z1: "<< z1_poly5<<", y1: "<<poly5->Eval(x1_poly5)<<endl;
+    if( TMath::Abs(poly5->Eval(x2_poly5)) > 1e-08){cout <<"No zero found, the lower point is y2: "<<poly5->Eval(x2_poly5)<<endl;}
+    else{
+      z2_poly5 = x2_poly5;
+      cout << "Zero found! z2: "<<z2_poly5<<", y2: "<<poly5->Eval(x2_poly5)<<endl;  
+      b_poly5 = true;
+    }
+  }  
+
+ 
+  //This has to be done before any other fit is performed
+  if(b_poly5 == true){
+
+    cout << "Poly5 behaves properly in the range ["<<low_limit<<", "<<high_limit<<"]"<<endl;
+    integral_poly5_val = poly5->Integral(z1_poly5,z2_poly5)/20.;
+    integral_poly5_err = poly5->IntegralError(z1_poly5,z2_poly5)/20.;
+    cout << "The integral between the two zeros is "<< integral_poly5_val <<"+/-"<<integral_poly5_err << endl;
+
+  }
+  
+  //MIchele - End
+  /*
+  TF1* fu_pulse2 = STTool::halfGauss2;
+  ge_pulse->Fit(fu_pulse2,"MQ","0");
+  Info("pulseFit","#chi^{2}/ndf: %6.2f/%d",fu_pulse2->GetChisquare(),fu_pulse2->GetNDF());
+  */
+
   
   TVirtualFitter *fit_volt = TVirtualFitter::GetFitter();
   Double_t* cov_volt       = fit_volt->GetCovarianceMatrix();
@@ -305,20 +450,64 @@
                                        2.0*cov_volt[2]/(tau_val*fu_pulse->GetParameter(0))
                                        );
                                        
+
+  //Michele 3/9/2017
+  //Estimate the  systematics
+  double int_sys = -999.0;
+  double b_int_sys = true; 
+
+  if(b_poly3 == true){
+
+    int_sys = TMath::Abs(integral_poly3_val - int_val)/2.;
+
+  }else{
+    
+    if(b_poly4 == true){
+
+      int_sys = TMath::Abs(integral_poly4_val - int_val)/2.;
+      
+    }else{
+
+      if(b_poly5 == true){
+
+	    int_sys = TMath::Abs(integral_poly5_val - int_val)/2.;
+      }else{
+	b_int_sys = false;
+      }
+
+    }
+    
+  }
+
+  //Michele - end
+
+
+
   
-  mg_pulse->Add(ge_pulse,"pe");
+  mg_pulse->Add(ge_pulse,"PE");
   
   
   
   Info("pulseFit","pulse Height: %8.4f+/-%8.4f",max_val,max_err);
   Info("pulseFit","         tau: %8.4f+/-%8.4f ns",tau_val,tau_err);
-  Info("pulseFit","   charge eq: %8.4f+/-%8.4f",int_val,int_err);
-  
+
+  //Michele
+  if(b_int_sys == true){
+    Info("pulseFit","   charge eq: %8.4f +/- %8.4f(stat) +/- %8.4f (sys)",int_val,int_err, int_sys);
+  }else{
+
+    Info("pulseFit","   charge eq: %8.4f +/- %8.4f(stat) +/- No sys eval.",int_val,int_err);
+  }
+  //
   // Plot pulse
   TCanvas* c_pul = new TCanvas("c_pul","Pulse Canvas",800,600);
   
   mg_pulse->Draw("A");
   fu_pulse->Draw("Same");
+  //fu_pulse2->Draw("Same");
+  poly3->Draw("Same");
+  poly4->Draw("Same");
+  poly5->Draw("Same");
   mg_pulse->Draw("");
   mg_pulse->GetXaxis()->SetLimits(-30.0,30.0);
   mg_pulse->GetXaxis()->SetTitle("#delta#it{t} [ns]");
@@ -332,6 +521,8 @@
   // 0,1: Maximal signal height
   v_res(0) = max_val;
   v_res(1) = max_err;
+  //v_res(1) = max_sys;
+
   // 2,3: tau (width of pulse)
   v_res(2) = tau_val;
   v_res(3) = tau_err;