Changeset 8105 for trunk/MagicSoft/Mars/mtools/MRolke.cc

Timestamp:

10/17/06 17:32:52 (18 years ago)

Author:

meyer

Message:

*** empty log message ***

File:

• trunk/MagicSoft/Mars/mtools/MRolke.cc (modified) (24 diffs)

trunk/MagicSoft/Mars/mtools/MRolke.cc

@@ -1 @@
-// @(#)root/physics:$Name: not supported by cvs2svn $:$Id: MRolke.cc,v 1.1 2006-10-17 12:51:31 meyer Exp $
+// @(#)root/physics:$Name: not supported by cvs2svn $:$Id: MRolke.cc,v 1.2 2006-10-17 16:32:52 meyer Exp $
 // Author: Jan Conrad    9/2/2004
 
@@ -115 +115 @@
 MRolke::MRolke(Double_t CL, Option_t * /*option*/)
 {
-    //constructor
-    fUpperLimit  = 0.0;
-    fLowerLimit  = 0.0;
-    fCL          = CL;
-    fSwitch      = 0; // 0: unbounded likelihood
-    // 1: bounded likelihood
+   //constructor
+   fUpperLimit  = 0.0;
+   fLowerLimit  = 0.0;
+   fCL          = CL;
+   fSwitch      = 0; // 0: unbounded likelihood
+                    // 1: bounded likelihood
 }
 
@@ -132 +132 @@
 Double_t MRolke::CalculateInterval(Int_t x, Int_t y, Int_t z, Double_t bm, Double_t em,Double_t e, Int_t mid, Double_t sde, Double_t sdb, Double_t tau, Double_t b, Int_t m)
 {
-    //calculate interval
-    Int_t done = 0;
-    Double_t limit[2];
-
-    limit[1] = Interval(x,y,z,bm,em,e,mid, sde,sdb,tau,b,m);
-
-    if (limit[1] > 0) {
-        done = 1;
-    }
-
-    if (fSwitch == 0) {
-
-        Int_t trial_x = x;
-
-        while (done == 0) {
-            trial_x++;
-            limit[1] = Interval(trial_x,y,z,bm,em,e,mid, sde,sdb,tau,b,m);
-            if (limit[1] > 0) done = 1;
-        }
-    }
-
-    return limit[1];
+   //calculate interval
+   Int_t done = 0;
+   Double_t limit[2];
+
+   limit[1] = Interval(x,y,z,bm,em,e,mid, sde,sdb,tau,b,m);
+
+   if (limit[1] > 0) {
+      done = 1;
+   }
+
+   if (fSwitch == 0) {
+
+      Int_t trial_x = x;
+
+      while (done == 0) {
+         trial_x++;
+         limit[1] = Interval(trial_x,y,z,bm,em,e,mid, sde,sdb,tau,b,m);
+         if (limit[1] > 0) done = 1;
+      } 
+   }
+
+   return limit[1];
 }
 
@@ -162 +162 @@
 Double_t MRolke::Interval(Int_t x, Int_t y, Int_t z, Double_t bm, Double_t em,Double_t e, Int_t mid, Double_t sde, Double_t sdb, Double_t tau, Double_t b, Int_t m)
 {
-    // Calculates the Confidence Interval
-
-    //Double_t dchi2 =  Chi2Percentile(1,1-fCL);
-    Double_t dchi2 = TMath::ChisquareQuantile(fCL, 1);
-
-    Double_t tempxy[2],limits[2] = {0,0};
-    Double_t slope,fmid,low,flow,high,fhigh,test,ftest,mu0,maximum,target,l,f0;
-    Double_t med = 0;
-    Double_t maxiter=1000, acc = 0.00001;
-    Int_t i;
-    Int_t bp = 0;
-
-    if ((mid != 3) && (mid != 5)) bm = (Double_t)y;
-
-    if ((mid == 3) || (mid == 5)) {
-        if (bm == 0) bm = 0.00001;
-    }
-
-    if ((mid <= 2) || (mid == 4)) bp = 1;
-
-
-    if (bp == 1 && x == 0 && bm > 0 ){
-
-        for(Int_t i = 0; i < 2; i++) {
-            x++;
-            tempxy[i] = Interval(x,y,z,bm,em,e,mid,sde,sdb,tau,b,m);
-        }
-        slope = tempxy[1] - tempxy[0];
-        limits[1] = tempxy[0] - slope;
-        limits[0] = 0.0;
-        if (limits[1] < 0) limits[1] = 0.0;
-        goto done;
-    }
-
-    if (bp != 1 && x == 0){
-
-        for(Int_t i = 0; i < 2; i++) {
-            x++;
-            tempxy[i] = Interval(x,y,z,bm,em,e,mid,sde,sdb,tau,b,m);
-      }
-        slope = tempxy[1] - tempxy[0];
-        limits[1] = tempxy[0] - slope;
-        limits[0] = 0.0;
-        if (limits[1] < 0) limits[1] = 0.0;
-        goto done;
-    }
-
-    if (bp != 1  && bm == 0){
-        for(Int_t i = 0; i < 2; i++) {
-            bm++;
-            limits[1] = Interval(x,y,z,bm,em,e,mid,sde,sdb,tau,b,m);
-            tempxy[i] = limits[1];
-        }
-        slope = tempxy[1] - tempxy[0];
-        limits[1] = tempxy[0] - slope;
-        if (limits[1] < 0) limits[1] = 0;
-        goto done;
-    }
-
-
-    if (x == 0 && bm == 0){
-        x++;
-        bm++;
-
-        limits[1] = Interval(x,y,z,bm,em,e,mid,sde,sdb,tau,b,m);
-        tempxy[0] = limits[1];
-        x  = 1;
-        bm = 2;
-        limits[1] = Interval(x,y,z,bm,em,e,mid,sde,sdb,tau,b,m);
-        tempxy[1] = limits[1];
-        x  = 2;
-        bm = 1;
-        limits[1] = Interval(x,y,z,bm,em,e,mid,sde,sdb,tau,b,m);
-        limits[1] = 3*tempxy[0] -tempxy[1] - limits[1];
-        if (limits[1] < 0) limits[1] = 0;
-        goto done;
-    }
-
-    mu0 = Likelihood(0,x,y,z,bm,em,e,mid,sde,sdb,tau,b,m,1);
-
-    maximum = Likelihood(0,x,y,z,bm,em,e,mid,sde,sdb,tau,b,m,2);
-
-    test = 0;
-
-    f0 = Likelihood(test,x,y,z,bm,em,e,mid,sde,sdb,tau,b,m,3);
-
-    if ( fSwitch == 1 ) {  // do this only for the unbounded likelihood case
-        if ( mu0 < 0 ) maximum = f0;
-    }
-
-    target = maximum - dchi2;
-
-    if (f0 > target) {
-        limits[0] = 0;
-    } else {
-        if (mu0 < 0){
-            limits[0] = 0;
-            limits[1] = 0;
-        }
-
-        low   = 0;
-        flow  = f0;
-        high  = mu0;
-        fhigh = maximum;
-
-        for(Int_t i = 0; i < maxiter; i++) {
-            l = (target-fhigh)/(flow-fhigh);
-            if (l < 0.2) l = 0.2;
-            if (l > 0.8) l = 0.8;
-
-            med = l*low + (1-l)*high;
-            if(med < 0.01){
-                limits[1]=0.0;
-                goto done;
-            }
-
-            fmid = Likelihood(med,x,y,z,bm,em,e,mid,sde,sdb,tau,b,m,3);
-
-            if (fmid > target) {
-                high  = med;
-                fhigh = fmid;
-            } else {
-                low  = med;
-                flow = fmid;
-            }
-            if ((high-low) < acc*high) break;
-      }
-        limits[0] = med;
-    }
-
-
-    if(mu0 > 0) {
-        low  = mu0;
-        flow = maximum;
-    } else {
-        low  = 0;
-        flow = f0;
-    }
-
-    test = low +1 ;
-
-    ftest = Likelihood(test,x,y,z,bm,em,e,mid,sde,sdb,tau,b,m,3);
-
-    if (ftest < target) {
-        high  = test;
-        fhigh = ftest;
-    } else {
-        slope = (ftest - flow)/(test - low);
-        high  = test + (target -ftest)/slope;
-        fhigh = Likelihood(high,x,y,z,bm,em,e,mid,sde,sdb,tau,b,m,3);
-   }
-
-    for(i = 0; i < maxiter; i++) {
-        l = (target-fhigh)/(flow-fhigh);
-        if (l < 0.2) l = 0.2;
-        if (l > 0.8) l = 0.8;
-        med  = l * low + (1.-l)*high;
-        fmid = Likelihood(med,x,y,z,bm,em,e,mid,sde,sdb,tau,b,m,3);
-
-        if (fmid < target) {
+   // Calculates the Confidence Interval
+
+   //Double_t dchi2 =  Chi2Percentile(1,1-fCL);
+   Double_t dchi2 = TMath::ChisquareQuantile(fCL, 1);
+
+   Double_t tempxy[2],limits[2] = {0,0};
+   Double_t slope,fmid,low,flow,high,fhigh,test,ftest,mu0,maximum,target,l,f0;
+   Double_t med = 0;
+   Double_t maxiter=1000, acc = 0.00001;
+   Int_t i;
+   Int_t bp = 0;
+
+   if ((mid != 3) && (mid != 5)) bm = (Double_t)y;
+
+   if ((mid == 3) || (mid == 5)) {
+      if (bm == 0) bm = 0.00001;
+   } 
+
+   if ((mid <= 2) || (mid == 4)) bp = 1;
+  
+
+   if (bp == 1 && x == 0 && bm > 0 ){
+
+      for(Int_t i = 0; i < 2; i++) {
+         x++;
+         tempxy[i] = Interval(x,y,z,bm,em,e,mid,sde,sdb,tau,b,m);
+      }
+      slope = tempxy[1] - tempxy[0];
+      limits[1] = tempxy[0] - slope;
+      limits[0] = 0.0;
+      if (limits[1] < 0) limits[1] = 0.0;
+      goto done;
+   }
+
+   if (bp != 1 && x == 0){
+
+      for(Int_t i = 0; i < 2; i++) {
+         x++;
+         tempxy[i] = Interval(x,y,z,bm,em,e,mid,sde,sdb,tau,b,m);
+      }
+      slope = tempxy[1] - tempxy[0];
+      limits[1] = tempxy[0] - slope;
+      limits[0] = 0.0;
+      if (limits[1] < 0) limits[1] = 0.0;
+      goto done;
+   }
+
+   if (bp != 1  && bm == 0){
+      for(Int_t i = 0; i < 2; i++) {
+         bm++;
+         limits[1] = Interval(x,y,z,bm,em,e,mid,sde,sdb,tau,b,m);
+         tempxy[i] = limits[1];
+      }
+      slope = tempxy[1] - tempxy[0];
+      limits[1] = tempxy[0] - slope;
+      if (limits[1] < 0) limits[1] = 0;
+      goto done;
+   }
+
+
+   if (x == 0 && bm == 0){
+      x++;
+      bm++;
+
+      limits[1] = Interval(x,y,z,bm,em,e,mid,sde,sdb,tau,b,m);
+      tempxy[0] = limits[1];
+      x  = 1;
+      bm = 2;
+      limits[1] = Interval(x,y,z,bm,em,e,mid,sde,sdb,tau,b,m);
+      tempxy[1] = limits[1];
+      x  = 2;
+      bm = 1;
+      limits[1] = Interval(x,y,z,bm,em,e,mid,sde,sdb,tau,b,m);
+      limits[1] = 3*tempxy[0] -tempxy[1] - limits[1];
+      if (limits[1] < 0) limits[1] = 0;
+      goto done;
+   }
+
+   mu0 = Likelihood(0,x,y,z,bm,em,e,mid,sde,sdb,tau,b,m,1);
+
+   maximum = Likelihood(0,x,y,z,bm,em,e,mid,sde,sdb,tau,b,m,2);
+
+   test = 0;
+
+   f0 = Likelihood(test,x,y,z,bm,em,e,mid,sde,sdb,tau,b,m,3);
+
+   if ( fSwitch == 1 ) {  // do this only for the unbounded likelihood case
+      if ( mu0 < 0 ) maximum = f0;
+   }
+
+   target = maximum - dchi2;
+
+   if (f0 > target) {
+      limits[0] = 0;
+   } else {
+      if (mu0 < 0){
+         limits[0] = 0;
+         limits[1] = 0;
+      }
+
+      low   = 0;
+      flow  = f0;
+      high  = mu0;
+      fhigh = maximum;
+
+      for(Int_t i = 0; i < maxiter; i++) {
+         l = (target-fhigh)/(flow-fhigh);
+         if (l < 0.2) l = 0.2;
+         if (l > 0.8) l = 0.8;
+
+         med = l*low + (1-l)*high;
+         if(med < 0.01){
+            limits[1]=0.0;                           
+            goto done;
+         }
+
+         fmid = Likelihood(med,x,y,z,bm,em,e,mid,sde,sdb,tau,b,m,3);
+
+         if (fmid > target) {
             high  = med;
             fhigh = fmid;
-        } else {
+         } else {
             low  = med;
             flow = fmid;
-        }
-        if (high-low < acc*high) break;
-    }
-    limits[1] = med;
+         }
+         if ((high-low) < acc*high) break;
+      }
+      limits[0] = med;
+   }
+
+
+   if(mu0 > 0) {
+      low  = mu0;
+      flow = maximum;
+   } else {
+      low  = 0;
+      flow = f0;
+   }
+
+   test = low +1 ;
+
+   ftest = Likelihood(test,x,y,z,bm,em,e,mid,sde,sdb,tau,b,m,3);
+
+   if (ftest < target) {
+      high  = test;
+      fhigh = ftest;
+   } else {
+      slope = (ftest - flow)/(test - low);
+      high  = test + (target -ftest)/slope;
+      fhigh = Likelihood(high,x,y,z,bm,em,e,mid,sde,sdb,tau,b,m,3);
+   }
+
+   for(i = 0; i < maxiter; i++) {
+      l = (target-fhigh)/(flow-fhigh);
+      if (l < 0.2) l = 0.2;
+      if (l > 0.8) l = 0.8;
+      med  = l * low + (1.-l)*high;
+      fmid = Likelihood(med,x,y,z,bm,em,e,mid,sde,sdb,tau,b,m,3);
+
+      if (fmid < target) {
+         high  = med;
+         fhigh = fmid;
+      } else {
+         low  = med;
+         flow = fmid;
+      }
+      if (high-low < acc*high) break;
+   }
+   limits[1] = med;
 
 done:
 
-    if ( (mid == 4) || (mid==5) ) {
-        limits[0] /= e;
-        limits[1] /= e;
-    }
-
-
-    fUpperLimit = limits[1];
-    fLowerLimit = TMath::Max(limits[0],0.0);
-
-
-    return limits[1];
+   if ( (mid == 4) || (mid==5) ) {
+      limits[0] /= e;
+      limits[1] /= e;
+   }
+
+
+   fUpperLimit = limits[1];
+   fLowerLimit = TMath::Max(limits[0],0.0);
+
+  
+   return limits[1];
 }
 
@@ -351 +351 @@
 Double_t MRolke::Likelihood(Double_t mu, Int_t x, Int_t y, Int_t z, Double_t bm,Double_t em, Double_t e, Int_t mid, Double_t sde, Double_t sdb, Double_t tau, Double_t b, Int_t m, Int_t what)
 {
-    // Chooses between the different profile likelihood functions to use for the
-    // different models.
-    // Returns evaluation of the profile likelihood functions.
-
-    switch (mid) {
-    case 1: return EvalLikeMod1(mu,x,y,z,e,tau,b,m,what);
-    case 2: return EvalLikeMod2(mu,x,y,em,e,sde,tau,b,what);
-    case 3: return EvalLikeMod3(mu,x,bm,em,e,sde,sdb,b,what);
-    case 4: return EvalLikeMod4(mu,x,y,tau,b,what);
-    case 5: return EvalLikeMod5(mu,x,bm,sdb,b,what);
-    case 6: return EvalLikeMod6(mu,x,z,e,b,m,what);
-    case 7: return EvalLikeMod7(mu,x,em,e,sde,b,what);
-    }
-
-    return 0;
+   // Chooses between the different profile likelihood functions to use for the
+   // different models.
+   // Returns evaluation of the profile likelihood functions.
+
+   switch (mid) {
+      case 1: return EvalLikeMod1(mu,x,y,z,e,tau,b,m,what);
+      case 2: return EvalLikeMod2(mu,x,y,em,e,sde,tau,b,what);
+      case 3: return EvalLikeMod3(mu,x,bm,em,e,sde,sdb,b,what);
+      case 4: return EvalLikeMod4(mu,x,y,tau,b,what);
+      case 5: return EvalLikeMod5(mu,x,bm,sdb,b,what);
+      case 6: return EvalLikeMod6(mu,x,z,e,b,m,what);
+      case 7: return EvalLikeMod7(mu,x,em,e,sde,b,what);
+   }
+
+   return 0;
 }
 
@@ -371 +371 @@
 Double_t MRolke::EvalLikeMod1(Double_t mu, Int_t x, Int_t y, Int_t z, Double_t e, Double_t tau, Double_t b, Int_t m, Int_t what)
 {
-    // Calculates the Profile Likelihood for MODEL 1:
-    //  Poisson background/ Binomial Efficiency
-    // what = 1: Maximum likelihood estimate is returned
-    // what = 2: Profile Likelihood of Maxmimum Likelihood estimate is returned.
-    // what = 3: Profile Likelihood of Test hypothesis is returned
-    // otherwise parameters as described in the beginning of the class)
-
-    Double_t f  = 0;
-    Double_t zm = Double_t(z)/m;
-
-    if (what == 1) {
-        f = (x-y/tau)/zm;
-    }
-
-    if (what == 2) {
-        mu = (x-y/tau)/zm;
-        b  = y/tau;
-        Double_t e = zm;
-        f = LikeMod1(mu,b,e,x,y,z,tau,m);
-    }
-
-    if (what == 3) {
-        if (mu == 0){
-            b = (x+y)/(1.0+tau);
-            e = zm;
-            f = LikeMod1(mu,b,e,x,y,z,tau,m);
-        } else {
-            MRolke g;
-            g.ProfLikeMod1(mu,b,e,x,y,z,tau,m);
-            f = LikeMod1(mu,b,e,x,y,z,tau,m);
-        }
-    }
-
-    return f;
+   // Calculates the Profile Likelihood for MODEL 1:
+   //  Poisson background/ Binomial Efficiency
+   // what = 1: Maximum likelihood estimate is returned
+   // what = 2: Profile Likelihood of Maxmimum Likelihood estimate is returned.
+   // what = 3: Profile Likelihood of Test hypothesis is returned
+   // otherwise parameters as described in the beginning of the class)
+
+   Double_t f  = 0;
+   Double_t zm = Double_t(z)/m;
+
+   if (what == 1) {
+      f = (x-y/tau)/zm;
+   }
+
+   if (what == 2) {
+      mu = (x-y/tau)/zm;
+      b  = y/tau;
+      Double_t e = zm;
+      f = LikeMod1(mu,b,e,x,y,z,tau,m);
+   }
+
+   if (what == 3) {
+      if (mu == 0){
+         b = (x+y)/(1.0+tau);
+         e = zm;
+         f = LikeMod1(mu,b,e,x,y,z,tau,m);
+      } else {
+         MRolke g;
+         g.ProfLikeMod1(mu,b,e,x,y,z,tau,m);
+         f = LikeMod1(mu,b,e,x,y,z,tau,m);
+      }
+   }
+
+   return f;
 }
 
@@ -410 +410 @@
 Double_t MRolke::LikeMod1(Double_t mu,Double_t b, Double_t e, Int_t x, Int_t y, Int_t z, Double_t tau, Int_t m)
 {
-    // Profile Likelihood function for MODEL 1:
-    // Poisson background/ Binomial Efficiency
-
-    return 2*(x*TMath::Log(e*mu+b)-(e*mu +b)-LogFactorial(x)+y*TMath::Log(tau*b)-tau*b-LogFactorial(y) + TMath::Log(TMath::Factorial(m)) - TMath::Log(TMath::Factorial(m-z)) - TMath::Log(TMath::Factorial(z))+ z * TMath::Log(e) + (m-z)*TMath::Log(1-e));
+   // Profile Likelihood function for MODEL 1:
+   // Poisson background/ Binomial Efficiency
+
+   return 2*(x*TMath::Log(e*mu+b)-(e*mu +b)-LogFactorial(x)+y*TMath::Log(tau*b)-tau*b-LogFactorial(y) + TMath::Log(TMath::Factorial(m)) - TMath::Log(TMath::Factorial(m-z)) - TMath::Log(TMath::Factorial(z))+ z * TMath::Log(e) + (m-z)*TMath::Log(1-e));
 }
 
@@ -419 +419 @@
 void MRolke::ProfLikeMod1(Double_t mu,Double_t &b,Double_t &e,Int_t x,Int_t y, Int_t z,Double_t tau,Int_t m)
 {
-    // Void needed to calculate estimates of efficiency and background for model 1
-
-    Double_t med = 0.0,fmid;
-    Int_t maxiter =1000;
-    Double_t acc = 0.00001;
-    Double_t emin = ((m+mu*tau)-TMath::Sqrt((m+mu*tau)*(m+mu*tau)-4 * mu* tau * z))/2/mu/tau;
-
-    Double_t low  = TMath::Max(1e-10,emin+1e-10);
-    Double_t high = 1 - 1e-10;
-
-    for(Int_t i = 0; i < maxiter; i++) {
-        med = (low+high)/2.;
-
-        fmid = LikeGradMod1(med,mu,x,y,z,tau,m);
-
-        if(high < 0.5) acc = 0.00001*high;
-        else           acc = 0.00001*(1-high);
-
-        if ((high - low) < acc*high) break;
-
-        if(fmid > 0) low  = med;
-        else         high = med;
-    }
-
-    e = med;
-    Double_t eta = Double_t(z)/e -Double_t(m-z)/(1-e);
-
-    b = Double_t(y)/(tau -eta/mu);
+   // Void needed to calculate estimates of efficiency and background for model 1
+
+   Double_t med = 0.0,fmid;
+   Int_t maxiter =1000;
+   Double_t acc = 0.00001;
+   Double_t emin = ((m+mu*tau)-TMath::Sqrt((m+mu*tau)*(m+mu*tau)-4 * mu* tau * z))/2/mu/tau;
+
+   Double_t low  = TMath::Max(1e-10,emin+1e-10);
+   Double_t high = 1 - 1e-10;
+
+   for(Int_t i = 0; i < maxiter; i++) {
+      med = (low+high)/2.;
+
+      fmid = LikeGradMod1(med,mu,x,y,z,tau,m);
+
+      if(high < 0.5) acc = 0.00001*high;
+      else           acc = 0.00001*(1-high);
+
+      if ((high - low) < acc*high) break;
+
+      if(fmid > 0) low  = med;
+      else         high = med;
+   }
+
+   e = med;
+   Double_t eta = Double_t(z)/e -Double_t(m-z)/(1-e);
+
+   b = Double_t(y)/(tau -eta/mu);
 }
 
@@ -453 +453 @@
 {
    //gradient model
-    Double_t eta, etaprime, bprime,f;
-    eta = static_cast<double>(z)/e - static_cast<double>(m-z)/(1.0 - e);
-    etaprime = (-1) * (static_cast<double>(m-z)/((1.0 - e)*(1.0 - e)) + static_cast<double>(z)/(e*e));
-    Double_t b = y/(tau - eta/mu);
-    bprime = (b*b * etaprime)/mu/y;
-    f =  (mu + bprime) * (x/(e * mu + b) - 1)+(y/b - tau) * bprime + eta;
-    return f;
+   Double_t eta, etaprime, bprime,f;
+   eta = static_cast<double>(z)/e - static_cast<double>(m-z)/(1.0 - e);
+   etaprime = (-1) * (static_cast<double>(m-z)/((1.0 - e)*(1.0 - e)) + static_cast<double>(z)/(e*e));
+   Double_t b = y/(tau - eta/mu);
+   bprime = (b*b * etaprime)/mu/y;
+   f =  (mu + bprime) * (x/(e * mu + b) - 1)+(y/b - tau) * bprime + eta;
+   return f;
 }
 
@@ -465 +465 @@
 Double_t MRolke::EvalLikeMod2(Double_t mu, Int_t x, Int_t y, Double_t em, Double_t e,Double_t sde, Double_t tau, Double_t b, Int_t what)
 {
-    // Calculates the Profile Likelihood for MODEL 2:
-    //  Poisson background/ Gauss Efficiency
-    // what = 1: Maximum likelihood estimate is returned
-    // what = 2: Profile Likelihood of Maxmimum Likelihood estimate is returned.
-    // what = 3: Profile Likelihood of Test hypothesis is returned
-    // otherwise parameters as described in the beginning of the class)
-
-    Double_t v =  sde*sde;
-    Double_t coef[4],roots[3];
-    Double_t f = 0;
-
-    if (what == 1) {
-        f = (x-y/tau)/em;
-    }
-
-    if (what == 2) {
-        mu = (x-y/tau)/em;
-        b = y/tau;
-        e = em;
-        f = LikeMod2(mu,b,e,x,y,em,tau,v);
-    }
-
-    if (what == 3) {
-        if (mu == 0 ) {
-            b = (x+y)/(1+tau);
-            f = LikeMod2(mu,b,e,x,y,em,tau,v);
+   // Calculates the Profile Likelihood for MODEL 2:
+   //  Poisson background/ Gauss Efficiency
+   // what = 1: Maximum likelihood estimate is returned
+   // what = 2: Profile Likelihood of Maxmimum Likelihood estimate is returned.
+   // what = 3: Profile Likelihood of Test hypothesis is returned
+   // otherwise parameters as described in the beginning of the class)
+
+   Double_t v =  sde*sde;
+   Double_t coef[4],roots[3];
+   Double_t f = 0;
+
+   if (what == 1) {
+      f = (x-y/tau)/em;
+   }
+
+   if (what == 2) {
+      mu = (x-y/tau)/em;
+      b = y/tau;
+      e = em;
+      f = LikeMod2(mu,b,e,x,y,em,tau,v);
+   }
+
+   if (what == 3) {
+      if (mu == 0 ) {
+         b = (x+y)/(1+tau);
+         f = LikeMod2(mu,b,e,x,y,em,tau,v);
       } else {
-            coef[3] = mu;
-            coef[2] = mu*mu*v-2*em*mu-mu*mu*v*tau;
-            coef[1] = ( - x)*mu*v - mu*mu*mu*v*v*tau - mu*mu*v*em + em*mu*mu*v*tau + em*em*mu - y*mu*v;
-            coef[0] = x*mu*mu*v*v*tau + x*em*mu*v - y*mu*mu*v*v + y*em*mu*v;
-
-            TMath::RootsCubic(coef,roots[0],roots[1],roots[2]);
-
-            e = roots[1];
-            b = y/(tau + (em - e)/mu/v);
-            f = LikeMod2(mu,b,e,x,y,em,tau,v);
-      }
-    }
-
-    return f;
+         coef[3] = mu;
+         coef[2] = mu*mu*v-2*em*mu-mu*mu*v*tau;
+         coef[1] = ( - x)*mu*v - mu*mu*mu*v*v*tau - mu*mu*v*em + em*mu*mu*v*tau + em*em*mu - y*mu*v;
+         coef[0] = x*mu*mu*v*v*tau + x*em*mu*v - y*mu*mu*v*v + y*em*mu*v;
+
+         TMath::RootsCubic(coef,roots[0],roots[1],roots[2]);
+
+         e = roots[1];
+         b = y/(tau + (em - e)/mu/v);
+         f = LikeMod2(mu,b,e,x,y,em,tau,v);
+      }
+   }
+
+   return f;
 }
 
@@ -512 +512 @@
 {
    // Profile Likelihood function for MODEL 2:
-    // Poisson background/Gauss Efficiency
-
-    return 2*(x*TMath::Log(e*mu+b)-(e*mu+b)-LogFactorial(x)+y*TMath::Log(tau*b)-tau*b-LogFactorial(y)-0.9189385-TMath::Log(v)/2-(em-e)*(em-e)/v/2);
+   // Poisson background/Gauss Efficiency
+
+   return 2*(x*TMath::Log(e*mu+b)-(e*mu+b)-LogFactorial(x)+y*TMath::Log(tau*b)-tau*b-LogFactorial(y)-0.9189385-TMath::Log(v)/2-(em-e)*(em-e)/v/2);
 }
 
@@ -521 +521 @@
 Double_t MRolke::EvalLikeMod3(Double_t mu, Int_t x, Double_t bm, Double_t em, Double_t e, Double_t sde, Double_t sdb, Double_t b, Int_t what)
 {
-    // Calculates the Profile Likelihood for MODEL 3:
-    // Gauss  background/ Gauss Efficiency
-    // what = 1: Maximum likelihood estimate is returned
-    // what = 2: Profile Likelihood of Maxmimum Likelihood estimate is returned.
-    // what = 3: Profile Likelihood of Test hypothesis is returned
-    // otherwise parameters as described in the beginning of the class)
-
-    Double_t f = 0.;
-    Double_t  v = sde*sde;
-    Double_t  u = sdb*sdb;
-
-    if (what == 1) {
-        f = (x-bm)/em;
-    }
-
-
-    if (what == 2) {
-        mu = (x-bm)/em;
-        b  = bm;
-        e  = em;
-        f  = LikeMod3(mu,b,e,x,bm,em,u,v);
-    }
-
-
-    if(what == 3) {
-        if(mu == 0.0){
-            b = ((bm-u)+TMath::Sqrt((bm-u)*(bm-u)+4*x*u))/2.;
-            e = em;
+   // Calculates the Profile Likelihood for MODEL 3:
+   // Gauss  background/ Gauss Efficiency
+   // what = 1: Maximum likelihood estimate is returned
+   // what = 2: Profile Likelihood of Maxmimum Likelihood estimate is returned.
+   // what = 3: Profile Likelihood of Test hypothesis is returned
+   // otherwise parameters as described in the beginning of the class)
+
+   Double_t f = 0.;
+   Double_t  v = sde*sde;
+   Double_t  u = sdb*sdb;
+
+   if (what == 1) {
+      f = (x-bm)/em;
+   }
+
+
+   if (what == 2) {
+      mu = (x-bm)/em;
+      b  = bm;
+      e  = em;
+      f  = LikeMod3(mu,b,e,x,bm,em,u,v);
+   }
+
+
+   if(what == 3) {
+      if(mu == 0.0){
+         b = ((bm-u)+TMath::Sqrt((bm-u)*(bm-u)+4*x*u))/2.;
+         e = em;
          f = LikeMod3(mu,b,e,x,bm,em,u,v);
-        } else {
-            Double_t temp[3];
-            temp[0] = mu*mu*v+u;
-            temp[1] = mu*mu*mu*v*v+mu*v*u-mu*mu*v*em+mu*v*bm-2*u*em;
-            temp[2] = mu*mu*v*v*bm-mu*v*u*em-mu*v*bm*em+u*em*em-mu*mu*v*v*x;
-            e = (-temp[1]+TMath::Sqrt(temp[1]*temp[1]-4*temp[0]*temp[2]))/2/temp[0];
+      } else {
+         Double_t temp[3];
+         temp[0] = mu*mu*v+u;
+         temp[1] = mu*mu*mu*v*v+mu*v*u-mu*mu*v*em+mu*v*bm-2*u*em;
+         temp[2] = mu*mu*v*v*bm-mu*v*u*em-mu*v*bm*em+u*em*em-mu*mu*v*v*x;
+         e = (-temp[1]+TMath::Sqrt(temp[1]*temp[1]-4*temp[0]*temp[2]))/2/temp[0];
          b = bm-(u*(em-e))/v/mu;
          f = LikeMod3(mu,b,e,x,bm,em,u,v);
-        }
-    }
-
-    return f;
+      }
+   }
+
+   return f;
 }
 
@@ -567 +567 @@
 Double_t MRolke::LikeMod3(Double_t mu,Double_t b,Double_t e,Int_t x,Double_t bm,Double_t em,Double_t u,Double_t v)
 {
-    // Profile Likelihood function for MODEL 3:
-    // Gauss background/Gauss Efficiency
-
-    return 2*(x * TMath::Log(e*mu+b)-(e*mu+b)-LogFactorial(x)-1.837877-TMath::Log(u)/2-(bm-b)*(bm-b)/u/2-TMath::Log(v)/2-(em-e)*(em-e)/v/2);
+   // Profile Likelihood function for MODEL 3:
+   // Gauss background/Gauss Efficiency
+
+   return 2*(x * TMath::Log(e*mu+b)-(e*mu+b)-LogFactorial(x)-1.837877-TMath::Log(u)/2-(bm-b)*(bm-b)/u/2-TMath::Log(v)/2-(em-e)*(em-e)/v/2);
 }
 
@@ -576 +576 @@
 Double_t MRolke::EvalLikeMod4(Double_t mu, Int_t x, Int_t y, Double_t tau, Double_t b, Int_t what)
 {
-    // Calculates the Profile Likelihood for MODEL 4:
-    // Poiss  background/Efficiency known
-    // what = 1: Maximum likelihood estimate is returned
-    // what = 2: Profile Likelihood of Maxmimum Likelihood estimate is returned.
-    // what = 3: Profile Likelihood of Test hypothesis is returned
-    // otherwise parameters as described in the beginning of the class)
-
-    Double_t f = 0.0;
-
-    if (what == 1) f = x-y/tau;
-    if (what == 2) {
-        mu = x-y/tau;
-        b  = Double_t(y)/tau;
-        f  = LikeMod4(mu,b,x,y,tau);
-    }
-    if (what == 3) {
-        if (mu == 0.0) {
-            b = Double_t(x+y)/(1+tau);
-            f = LikeMod4(mu,b,x,y,tau);
-        } else {
-            b = (x+y-(1+tau)*mu+sqrt((x+y-(1+tau)*mu)*(x+y-(1+tau)*mu)+4*(1+tau)*y*mu))/2/(1+tau);
-            f = LikeMod4(mu,b,x,y,tau);
-        }
-    }
-    return f;
+   // Calculates the Profile Likelihood for MODEL 4:
+   // Poiss  background/Efficiency known
+   // what = 1: Maximum likelihood estimate is returned
+   // what = 2: Profile Likelihood of Maxmimum Likelihood estimate is returned.
+   // what = 3: Profile Likelihood of Test hypothesis is returned
+   // otherwise parameters as described in the beginning of the class)
+
+   Double_t f = 0.0;
+
+   if (what == 1) f = x-y/tau;
+   if (what == 2) {
+      mu = x-y/tau;
+      b  = Double_t(y)/tau;
+      f  = LikeMod4(mu,b,x,y,tau);
+   }
+   if (what == 3) {
+      if (mu == 0.0) {
+         b = Double_t(x+y)/(1+tau);
+         f = LikeMod4(mu,b,x,y,tau);
+      } else {
+         b = (x+y-(1+tau)*mu+sqrt((x+y-(1+tau)*mu)*(x+y-(1+tau)*mu)+4*(1+tau)*y*mu))/2/(1+tau);
+         f = LikeMod4(mu,b,x,y,tau);
+      }
+   }
+   return f;
 }
 
@@ -606 +606 @@
 Double_t MRolke::LikeMod4(Double_t mu,Double_t b,Int_t x,Int_t y,Double_t tau)
 {
-    // Profile Likelihood function for MODEL 4:
-    // Poiss background/Efficiency known
-
-    return 2*(x*TMath::Log(mu+b)-(mu+b)-LogFactorial(x)+y*TMath::Log(tau*b)-tau*b-LogFactorial(y) );
+   // Profile Likelihood function for MODEL 4:
+   // Poiss background/Efficiency known
+
+   return 2*(x*TMath::Log(mu+b)-(mu+b)-LogFactorial(x)+y*TMath::Log(tau*b)-tau*b-LogFactorial(y) );
 }
 
@@ -615 +615 @@
 Double_t MRolke::EvalLikeMod5(Double_t mu, Int_t x, Double_t bm, Double_t sdb, Double_t b, Int_t what)
 {
-    // Calculates the Profile Likelihood for MODEL 5:
-    // Gauss  background/Efficiency known
-    // what = 1: Maximum likelihood estimate is returned
-    // what = 2: Profile Likelihood of Maxmimum Likelihood estimate is returned.
-    // what = 3: Profile Likelihood of Test hypothesis is returned
-    // otherwise parameters as described in the beginning of the class)
-
-    Double_t u=sdb*sdb;
-    Double_t f = 0;
-
-    if(what == 1) {
+   // Calculates the Profile Likelihood for MODEL 5:
+   // Gauss  background/Efficiency known
+   // what = 1: Maximum likelihood estimate is returned
+   // what = 2: Profile Likelihood of Maxmimum Likelihood estimate is returned.
+   // what = 3: Profile Likelihood of Test hypothesis is returned
+   // otherwise parameters as described in the beginning of the class)
+
+   Double_t u=sdb*sdb;
+   Double_t f = 0;
+
+   if(what == 1) {
       f = x - bm;
-    }
-    if(what == 2) {
-        mu = x-bm;
-        b  = bm;
-        f  = LikeMod5(mu,b,x,bm,u);
-    }
-
-    if (what == 3) {
-        b = ((bm-u-mu)+TMath::Sqrt((bm-u-mu)*(bm-u-mu)-4*(mu*u-mu*bm-u*x)))/2;
-        f = LikeMod5(mu,b,x,bm,u);
-    }
-    return f;
+   }
+   if(what == 2) {
+      mu = x-bm;
+      b  = bm;
+      f  = LikeMod5(mu,b,x,bm,u);
+   }
+
+   if (what == 3) {
+      b = ((bm-u-mu)+TMath::Sqrt((bm-u-mu)*(bm-u-mu)-4*(mu*u-mu*bm-u*x)))/2;
+      f = LikeMod5(mu,b,x,bm,u);
+   }
+   return f;
 }
 
@@ -644 +644 @@
 Double_t MRolke::LikeMod5(Double_t mu,Double_t b,Int_t x,Double_t bm,Double_t u)
 {
-    // Profile Likelihood function for MODEL 5:
-    // Gauss background/Efficiency known
-
-    return 2*(x*TMath::Log(mu+b)-(mu + b)-LogFactorial(x)-0.9189385-TMath::Log(u)/2-((bm-b)*(bm-b))/u/2);
+   // Profile Likelihood function for MODEL 5:
+   // Gauss background/Efficiency known
+
+   return 2*(x*TMath::Log(mu+b)-(mu + b)-LogFactorial(x)-0.9189385-TMath::Log(u)/2-((bm-b)*(bm-b))/u/2);
 }
 
@@ -653 +653 @@
 Double_t MRolke::EvalLikeMod6(Double_t mu, Int_t x, Int_t z, Double_t e, Double_t b, Int_t m, Int_t what)
 {
-    // Calculates the Profile Likelihood for MODEL 6:
-    // Gauss  known/Efficiency binomial
-    // what = 1: Maximum likelihood estimate is returned
-    // what = 2: Profile Likelihood of Maxmimum Likelihood estimate is returned.
-    // what = 3: Profile Likelihood of Test hypothesis is returned
-    // otherwise parameters as described in the beginning of the class)
-
-    Double_t coef[4],roots[3];
-    Double_t f = 0.;
-    Double_t zm = Double_t(z)/m;
-
-    if(what==1){
-        f = (x-b)/zm;
-    }
-
-    if(what==2){
-        mu = (x-b)/zm;
-        e  = zm;
-        f  = LikeMod6(mu,b,e,x,z,m);
-    }
-    if(what == 3){
-        if(mu==0){
-            e = zm;
-        } else {
-            coef[3] = mu*mu;
-            coef[2] = mu * b - mu * x - mu*mu - mu * m;
-            coef[1] = mu * x - mu * b + mu * z - m * b;
-            coef[0] = b * z;
-            TMath::RootsCubic(coef,roots[0],roots[1],roots[2]);
-            e = roots[1];
-        }
-        f =LikeMod6(mu,b,e,x,z,m);
-    }
-    return f;
+   // Calculates the Profile Likelihood for MODEL 6:
+   // Gauss  known/Efficiency binomial
+   // what = 1: Maximum likelihood estimate is returned
+   // what = 2: Profile Likelihood of Maxmimum Likelihood estimate is returned.
+   // what = 3: Profile Likelihood of Test hypothesis is returned
+   // otherwise parameters as described in the beginning of the class)
+
+   Double_t coef[4],roots[3];
+   Double_t f = 0.;
+   Double_t zm = Double_t(z)/m;
+
+   if(what==1){
+      f = (x-b)/zm;
+   }
+
+   if(what==2){
+      mu = (x-b)/zm;
+      e  = zm;
+      f  = LikeMod6(mu,b,e,x,z,m);
+   }
+   if(what == 3){
+      if(mu==0){
+         e = zm;
+      } else {
+         coef[3] = mu*mu;
+         coef[2] = mu * b - mu * x - mu*mu - mu * m;
+         coef[1] = mu * x - mu * b + mu * z - m * b;
+         coef[0] = b * z;
+         TMath::RootsCubic(coef,roots[0],roots[1],roots[2]);
+         e = roots[1];
+      }
+      f =LikeMod6(mu,b,e,x,z,m);
+   }
+   return f;
 }
 
@@ -692 +692 @@
 Double_t MRolke::LikeMod6(Double_t mu,Double_t b,Double_t e,Int_t x,Int_t z,Int_t m)
 {
-    // Profile Likelihood function for MODEL 6:
-    // background known/ Efficiency binomial
-
-    Double_t f = 0.0;
-
-    if (z > 100 || m > 100) {
-        f = 2*(x*TMath::Log(e*mu+b)-(e*mu+b)-LogFactorial(x)+(m*TMath::Log(m)  - m)-(z*TMath::Log(z) - z)  - ((m-z)*TMath::Log(m-z) - m + z)+z*TMath::Log(e)+(m-z)*TMath::Log(1-e));
-    } else {
-        f = 2*(x*TMath::Log(e*mu+b)-(e*mu+b)-LogFactorial(x)+TMath::Log(TMath::Factorial(m))-TMath::Log(TMath::Factorial(z))-TMath::Log(TMath::Factorial(m-z))+z*TMath::Log(e)+(m-z)*TMath::Log(1-e));
-    }
-    return f;
+   // Profile Likelihood function for MODEL 6:
+   // background known/ Efficiency binomial
+
+   Double_t f = 0.0;
+
+   if (z > 100 || m > 100) {
+      f = 2*(x*TMath::Log(e*mu+b)-(e*mu+b)-LogFactorial(x)+(m*TMath::Log(m)  - m)-(z*TMath::Log(z) - z)  - ((m-z)*TMath::Log(m-z) - m + z)+z*TMath::Log(e)+(m-z)*TMath::Log(1-e));
+   } else {
+      f = 2*(x*TMath::Log(e*mu+b)-(e*mu+b)-LogFactorial(x)+TMath::Log(TMath::Factorial(m))-TMath::Log(TMath::Factorial(z))-TMath::Log(TMath::Factorial(m-z))+z*TMath::Log(e)+(m-z)*TMath::Log(1-e));
+   }
+   return f;
 }
 
@@ -709 +709 @@
 Double_t MRolke::EvalLikeMod7(Double_t mu, Int_t x, Double_t em, Double_t e, Double_t sde, Double_t b, Int_t what)
 {
-    // Calculates the Profile Likelihood for MODEL 7:
-    // background known/Efficiency Gauss
-    // what = 1: Maximum likelihood estimate is returned
-    // what = 2: Profile Likelihood of Maxmimum Likelihood estimate is returned.
-    // what = 3: Profile Likelihood of Test hypothesis is returned
-    // otherwise parameters as described in the beginning of the class)
-
-    Double_t v=sde*sde;
-    Double_t f = 0.;
-
-    if(what ==  1) {
-        f = (x-b)/em;
-    }
-
-    if(what == 2) {
-        mu = (x-b)/em;
-        e  = em;
-        f  = LikeMod7(mu, b, e, x, em, v);
-    }
-
-    if(what == 3) {
-        if(mu==0) {
-            e = em;
-        } else {
-            e = ( -(mu*em-b-mu*mu*v)-TMath::Sqrt((mu*em-b-mu*mu*v)*(mu*em-b-mu*mu*v)+4*mu*(x*mu*v-mu*b*v + b * em)))/( - mu)/2;
-        }
-        f = LikeMod7(mu, b, e, x, em, v);
-    }
-
-    return f;
+   // Calculates the Profile Likelihood for MODEL 7:
+   // background known/Efficiency Gauss
+   // what = 1: Maximum likelihood estimate is returned
+   // what = 2: Profile Likelihood of Maxmimum Likelihood estimate is returned.
+   // what = 3: Profile Likelihood of Test hypothesis is returned
+   // otherwise parameters as described in the beginning of the class)
+
+   Double_t v=sde*sde;
+   Double_t f = 0.;
+
+   if(what ==  1) {
+      f = (x-b)/em;
+   }
+
+   if(what == 2) {
+      mu = (x-b)/em;
+      e  = em;
+      f  = LikeMod7(mu, b, e, x, em, v);
+   }
+
+   if(what == 3) {
+      if(mu==0) {
+         e = em;
+      } else {
+         e = ( -(mu*em-b-mu*mu*v)-TMath::Sqrt((mu*em-b-mu*mu*v)*(mu*em-b-mu*mu*v)+4*mu*(x*mu*v-mu*b*v + b * em)))/( - mu)/2;
+      }
+      f = LikeMod7(mu, b, e, x, em, v);
+   }
+
+   return f;
 }
 
@@ -744 +744 @@
 Double_t MRolke::LikeMod7(Double_t mu,Double_t b,Double_t e,Int_t x,Double_t em,Double_t v)
 {
-    // Profile Likelihood function for MODEL 6:
-    // background known/ Efficiency binomial
-
-    return 2*(x*TMath::Log(e*mu+b)-(e*mu + b)-LogFactorial(x)-0.9189385-TMath::Log(v)/2-(em-e)*(em-e)/v/2);
+   // Profile Likelihood function for MODEL 6:
+   // background known/ Efficiency binomial
+
+   return 2*(x*TMath::Log(e*mu+b)-(e*mu + b)-LogFactorial(x)-0.9189385-TMath::Log(v)/2-(em-e)*(em-e)/v/2);
 }
 
@@ -753 +753 @@
 Double_t MRolke::EvalPolynomial(Double_t x, const Int_t  coef[], Int_t N)
 {
-    // evaluate polynomial
-
-    const Int_t   *p;
-    p = coef;
-    Double_t ans = *p++;
-    Int_t i = N;
-
-    do
-        ans = ans * x  +  *p++;
-    while( --i );
-
-    return ans;
+  // evaluate polynomial
+
+   const Int_t   *p;
+   p = coef;
+   Double_t ans = *p++;
+   Int_t i = N;
+
+   do
+      ans = ans * x  +  *p++;
+   while( --i );
+
+   return ans;
 }
 
@@ -770 +770 @@
 Double_t MRolke::EvalMonomial(Double_t x, const Int_t coef[], Int_t N)
 {
-    // evaluate mononomial
-
-    Double_t ans;
-    const Int_t   *p;
-
-    p   = coef;
-    ans = x + *p++;
-    Int_t i = N-1;
-
-    do
-        ans = ans * x  + *p++;
-    while( --i );
-
-    return ans;
+   // evaluate mononomial
+
+   Double_t ans;
+   const Int_t   *p;
+
+   p   = coef;
+   ans = x + *p++;
+   Int_t i = N-1;
+
+   do
+      ans = ans * x  + *p++;
+   while( --i );
+
+   return ans;
 }
 //--------------------------------------------------------------------------
@@ -793 +793 @@
 Double_t MRolke::LogFactorial(Int_t n)
 {
-    if (n>69)
-        return n*TMath::Log(n)-n + 0.5*TMath::Log(TMath::TwoPi()*n);
-    else
-        return TMath::Log(TMath::Factorial(n));
-}
+   if (n>69)
+       return n*TMath::Log(n)-n + 0.5*TMath::Log(TMath::TwoPi()*n);
+   else
+       return TMath::Log(TMath::Factorial(n));
+}