// @(#)root/hist:$Id$
// Author: Christophe.Delaere@cern.ch   21/08/2002

///////////////////////////////////////////////////////////////////////////
//
// TLimit
//
// Class to compute 95% CL limits
//
// adapted from the mclimit code from Tom Junk (CLs method)
// see http://root.cern.ch/root/doc/TomJunk.pdf
// see http://cern.ch/thomasj/searchlimits/ecl.html
// see: Tom Junk,NIM A434, p. 435-443, 1999
//
// see also the following interesting references:
//  Alex Read, "Presentation of search results: the CLs technique"
//  Journal of Physics G: Nucl. Part. Phys. 28 2693-2704 (2002).
//  http://www.iop.org/EJ/abstract/0954-3899/28/10/313/
//
// A nice article is also available in the CERN yellow report with the proceeding
// of the 2000 CERN workshop on confidence intervals.
//
//  Alex Read, "Modified Frequentist Analysis of Search Results (The CLs Method)"
//  CERN 2000-005 (30 May 2000)
//
// see note about: "Should I use TRolke, TFeldmanCousins, TLimit?"
//  in the TRolke class description.
//
///////////////////////////////////////////////////////////////////////////

#include "TLimit.h"
#include "TArrayD.h"
#include "TVectorD.h"
#include "TOrdCollection.h"
#include "TConfidenceLevel.h"
#include "TLimitDataSource.h"
#include "TRandom3.h"
#include "TH1.h"
#include "TObjArray.h"
#include "TMath.h"
#include "TIterator.h"
#include "TObjString.h"
#include "TClassTable.h"
#include "Riostream.h"

ClassImp(TLimit)

TArrayD *TLimit::fgTable = new TArrayD(0);
TOrdCollection *TLimit::fgSystNames = new TOrdCollection();

TConfidenceLevel *TLimit::ComputeLimit(TLimitDataSource * data,
                                       Int_t nmc, bool stat,
                                       TRandom * generator)
{
   // class TLimit
   // ------------
   // Algorithm to compute 95% C.L. limits using the Likelihood ratio
   // semi-bayesian method.
   // It takes signal, background and data histograms wrapped in a
   // TLimitDataSource as input and runs a set of Monte Carlo experiments in
   // order to compute the limits. If needed, inputs are fluctuated according
   // to systematics. The output is a TConfidenceLevel.
   //
   // class TLimitDataSource
   // ----------------------
   //
   // Takes the signal, background and data histograms as well as different
   // systematics sources to form the TLimit input.
   //
   //  class TConfidenceLevel
   //  ----------------------
   //
   // Final result of the TLimit algorithm. It is created just after the
   // time-consuming part and can be stored in a TFile for further processing.
   // It contains light methods to return CLs, CLb and other interesting
   // quantities.
   //
   // The actual algorithm...
   // From an input (TLimitDataSource) it produces an output TConfidenceLevel.
   // For this, nmc Monte Carlo experiments are performed.
   // As usual, the larger this number, the longer the compute time,
   // but the better the result.
   //Begin_Html
   /*
   <FONT SIZE=+0>
   <p>Supposing that there is a plotfile.root file containing 3 histograms
           (signal, background and data), you can imagine doing things like:</p>
   <p>
   <BLOCKQUOTE><PRE>
    TFile* infile=new TFile("plotfile.root","READ");
    infile->cd();
    TH1* sh=(TH1*)infile->Get("signal");
    TH1* bh=(TH1*)infile->Get("background");
    TH1* dh=(TH1*)infile->Get("data");
    TLimitDataSource* mydatasource = new TLimitDataSource(sh,bh,dh);
    TConfidenceLevel *myconfidence = TLimit::ComputeLimit(mydatasource,50000);
    std::cout &lt&lt "  CLs    : " &lt&lt myconfidence->CLs()  &lt&lt std::endl;
    std::cout &lt&lt "  CLsb   : " &lt&lt myconfidence->CLsb() &lt&lt std::endl;
    std::cout &lt&lt "  CLb    : " &lt&lt myconfidence->CLb()  &lt&lt std::endl;
    std::cout &lt&lt "&lt CLs &gt  : " &lt&lt myconfidence->GetExpectedCLs_b()  &lt&lt std::endl;
    std::cout &lt&lt "&lt CLsb &gt : " &lt&lt myconfidence->GetExpectedCLsb_b() &lt&lt std::endl;
    std::cout &lt&lt "&lt CLb &gt  : " &lt&lt myconfidence->GetExpectedCLb_b()  &lt&lt std::endl;
    delete myconfidence;
    delete mydatasource;
    infile->Close();
   </PRE></BLOCKQUOTE></p>
   <p></p>
   <p>More information can still be found on
   <a HREF="http://cern.ch/aleph-proj-alphapp/doc/tlimit.html">this</a> page.</p>
   </FONT>
   */
   //End_Html

   // The final object returned...
   TConfidenceLevel *result = new TConfidenceLevel(nmc);
   // The random generator used...
   TRandom *myrandom = generator ? generator : new TRandom3;
   // Compute some total quantities on all the channels
   Int_t maxbins = 0;
   Double_t nsig = 0;
   Double_t nbg  = 0;
   Int_t ncand   = 0;
   Int_t i;
   for (i = 0; i <= data->GetSignal()->GetLast(); i++) {
      maxbins = (((TH1 *) (data->GetSignal()->At(i)))->GetNbinsX() + 2) > maxbins ?
                 (((TH1 *) (data->GetSignal()->At(i)))->GetNbinsX() + 2) : maxbins;
      nsig   +=  ((TH1 *) (data->GetSignal()->At(i)))->Integral();
      nbg    +=  ((TH1 *) (data->GetBackground()->At(i)))->Integral();
      ncand  +=  (Int_t) ((TH1 *) (data->GetCandidates()->At(i)))->Integral();
   }
   result->SetBtot(nbg);
   result->SetStot(nsig);
   result->SetDtot(ncand);
   Double_t buffer = 0;
   fgTable->Set(maxbins * (data->GetSignal()->GetLast() + 1));
   for (Int_t channel = 0; channel <= data->GetSignal()->GetLast(); channel++)
      for (Int_t bin = 0;
           bin <= ((TH1 *) (data->GetSignal()->At(channel)))->GetNbinsX()+1;
           bin++) {
         Double_t s = (Double_t) ((TH1 *) (data->GetSignal()->At(channel)))->GetBinContent(bin);
         Double_t b = (Double_t) ((TH1 *) (data->GetBackground()->At(channel)))->GetBinContent(bin);
         Double_t d = (Double_t) ((TH1 *) (data->GetCandidates()->At(channel)))->GetBinContent(bin);
         // Compute the value of the "-2lnQ" for the actual data
         if ((b == 0) && (s > 0)) {
            std::cout << "WARNING: Ignoring bin " << bin << " of channel "
                 << channel << " which has s=" << s << " but b=" << b << std::endl;
            std::cout << "         Maybe the MC statistic has to be improved..." << std::endl;
         }
         if ((s > 0) && (b > 0))
            buffer += LogLikelihood(s, b, b, d);
         // precompute the log(1+s/b)'s in an array to speed up computation
         // background-free bins are set to have a maximum t.s. value
         // for protection (corresponding to s/b of about 5E8)
         if ((s > 0) && (b > 0))
            fgTable->AddAt(LogLikelihood(s, b, b, 1), (channel * maxbins) + bin);
         else if ((s > 0) && (b == 0))
            fgTable->AddAt(20, (channel * maxbins) + bin);
      }
   result->SetTSD(buffer);
   // accumulate MC experiments.  Hold the test statistic function fixed, but
   // fluctuate s and b within syst. errors for computing probabilities of
   // having that outcome.  (Alex Read's prescription -- errors are on the ensemble,
   // not on the observed test statistic.  This technique does not split outcomes.)
   // keep the tstats as sum log(1+s/b). convert to -2lnQ when preparing the results
   // (reason -- like to keep the < signs right)
   Double_t *tss = new Double_t[nmc];
   Double_t *tsb = new Double_t[nmc];
   Double_t *lrs = new Double_t[nmc];
   Double_t *lrb = new Double_t[nmc];
   TLimitDataSource* tmp_src = (TLimitDataSource*)(data->Clone());
   TLimitDataSource* tmp_src2 = (TLimitDataSource*)(data->Clone());
   for (i = 0; i < nmc; i++) {
      tss[i] = 0;
      tsb[i] = 0;
      lrs[i] = 0;
      lrb[i] = 0;
      // fluctuate signal and background
      TLimitDataSource*  fluctuated = Fluctuate(data, tmp_src, !i, myrandom, stat) ? tmp_src : data;
      // make an independent set of fluctuations for reweighting pseudoexperiments
      // it turns out using the same fluctuated ones for numerator and denominator
      // is biased.
      TLimitDataSource*  fluctuated2 = Fluctuate(data, tmp_src2, false, myrandom, stat) ? tmp_src2 : data;

      for (Int_t channel = 0;
           channel <= fluctuated->GetSignal()->GetLast(); channel++) {
         for (Int_t bin = 0;
              bin <=((TH1 *) (fluctuated->GetSignal()->At(channel)))->GetNbinsX()+1;
              bin++) {
            if ((Double_t) ((TH1 *) (fluctuated->GetSignal()->At(channel)))->GetBinContent(bin) != 0) {
               // s+b hypothesis
               Double_t rate = (Double_t) ((TH1 *) (fluctuated->GetSignal()->At(channel)))->GetBinContent(bin) +
                               (Double_t) ((TH1 *) (fluctuated->GetBackground()->At(channel)))->GetBinContent(bin);
               Double_t rand = myrandom->Poisson(rate);
               tss[i] += rand * fgTable->At((channel * maxbins) + bin);
               Double_t s = (Double_t) ((TH1 *) (fluctuated->GetSignal()->At(channel)))->GetBinContent(bin);
               Double_t s2= (Double_t) ((TH1 *) (fluctuated2->GetSignal()->At(channel)))->GetBinContent(bin);
               Double_t b = (Double_t) ((TH1 *) (fluctuated->GetBackground()->At(channel)))->GetBinContent(bin);
               Double_t b2= (Double_t) ((TH1 *) (fluctuated2->GetBackground()->At(channel)))->GetBinContent(bin);
               if ((s > 0) && (b2 > 0))
                  lrs[i] += LogLikelihood(s, b, b2, rand) - s - b + b2;
               else if ((s > 0) && (b2 == 0))
                  lrs[i] += 20 * rand - s;
               // b hypothesis
               rate = (Double_t) ((TH1 *) (fluctuated->GetBackground()->At(channel)))->GetBinContent(bin);
               rand = myrandom->Poisson(rate);
               tsb[i] += rand * fgTable->At((channel * maxbins) + bin);
               if ((s2 > 0) && (b > 0))
                  lrb[i] += LogLikelihood(s2, b2, b, rand) - s2 - b2 + b;
               else if ((s > 0) && (b == 0))
                  lrb[i] += 20 * rand - s;
            }
         }
      }
      lrs[i] = TMath::Exp(lrs[i] < 710 ? lrs[i] : 710);
      lrb[i] = TMath::Exp(lrb[i] < 710 ? lrb[i] : 710);
   }
   delete tmp_src;
   delete tmp_src2;
   // lrs and lrb are the LR's (no logs) = prob(s+b)/prob(b) for
   // that choice of s and b within syst. errors in the ensemble.  These are
   // the MC experiment weights for relating the s+b and b PDF's of the unsmeared
   // test statistic (in which cas one can use another test statistic if one likes).

   // Now produce the output object.
   // The final quantities are computed on-demand form the arrays tss, tsb, lrs and lrb.
   result->SetTSS(tss);
   result->SetTSB(tsb);
   result->SetLRS(lrs);
   result->SetLRB(lrb);
   if (!generator)
      delete myrandom;
   return result;
}

bool TLimit::Fluctuate(TLimitDataSource * input, TLimitDataSource * output,
                       bool init, TRandom * generator, bool stat)
{
   // initialisation: create a sorted list of all the names of systematics
   if (init) {
      // create a "map" with the systematics names
      TIterator *errornames = input->GetErrorNames()->MakeIterator();
      TObjArray *listofnames = 0;
      delete fgSystNames;
      fgSystNames = new TOrdCollection();
      while ((listofnames = ((TObjArray *) errornames->Next()))) {
         TObjString *name = 0;
         TIterator *loniter = listofnames->MakeIterator();
         while ((name = (TObjString *) (loniter->Next())))
            if ((fgSystNames->IndexOf(name)) < 0)
               fgSystNames->AddLast(name);
      }
      fgSystNames->Sort();
   }
   // if the output is not given, create it from the input
   if (!output)
   output = (TLimitDataSource*)(input->Clone());
   // if there are no systematics, just returns the input as "fluctuated" output
   if ((fgSystNames->GetSize() <= 0)&&(!stat)) {
      return 0;
   }
   // if there are only stat, just fluctuate stats.
   if (fgSystNames->GetSize() <= 0) {
      output->SetOwner();
      for (Int_t channel = 0; channel <= input->GetSignal()->GetLast(); channel++) {
         TH1 *newsignal = (TH1*)(output->GetSignal()->At(channel));
         TH1 *oldsignal = (TH1*)(input->GetSignal()->At(channel));
         if(stat)
            for(int i=1; i<=newsignal->GetNbinsX(); i++) {
               newsignal->SetBinContent(i,oldsignal->GetBinContent(i)+generator->Gaus(0,oldsignal->GetBinError(i)));
            }
            newsignal->SetDirectory(0);
            TH1 *newbackground = (TH1*)(output->GetBackground()->At(channel));
            TH1 *oldbackground = (TH1*)(input->GetBackground()->At(channel));
            if(stat)
               for(int i=1; i<=newbackground->GetNbinsX(); i++)
                  newbackground->SetBinContent(i,oldbackground->GetBinContent(i)+generator->Gaus(0,oldbackground->GetBinError(i)));
            newbackground->SetDirectory(0);
      }
      return 1;
   }
   // Find a choice for the random variation and
   // re-toss all random numbers if any background or signal
   // goes negative.  (background = 0 is bad too, so put a little protection
   // around it -- must have at least 10% of the bg estimate).
   Bool_t retoss   = kTRUE;
   Double_t *serrf = new Double_t[(input->GetSignal()->GetLast()) + 1];
   Double_t *berrf = new Double_t[(input->GetSignal()->GetLast()) + 1];
   do {
      Double_t *toss = new Double_t[fgSystNames->GetSize()];
      for (Int_t i = 0; i < fgSystNames->GetSize(); i++)
         toss[i] = generator->Gaus(0, 1);
      retoss = kFALSE;
      for (Int_t channel = 0;
           channel <= input->GetSignal()->GetLast();
           channel++) {
         serrf[channel] = 0;
         berrf[channel] = 0;
         for (Int_t bin = 0;
              bin <((TVectorD *) (input->GetErrorOnSignal()->At(channel)))->GetNrows();
              bin++) {
            serrf[channel] += ((TVectorD *) (input->GetErrorOnSignal()->At(channel)))->operator[](bin) *
                toss[fgSystNames->BinarySearch((TObjString*) (((TObjArray *) (input->GetErrorNames()->At(channel)))->At(bin)))];
            berrf[channel] += ((TVectorD *) (input->GetErrorOnBackground()->At(channel)))->operator[](bin) *
                toss[fgSystNames->BinarySearch((TObjString*) (((TObjArray *) (input->GetErrorNames()->At(channel)))->At(bin)))];
         }
         if ((serrf[channel] < -1.0) || (berrf[channel] < -0.9)) {
            retoss = kTRUE;
            continue;
         }
      }
      delete[]toss;
   } while (retoss);
   // adjust the fluctuated signal and background counts with a legal set
   // of random fluctuations above.
   output->SetOwner();
   for (Int_t channel = 0; channel <= input->GetSignal()->GetLast();
        channel++) {
      TH1 *newsignal = (TH1*)(output->GetSignal()->At(channel));
      TH1 *oldsignal = (TH1*)(input->GetSignal()->At(channel));
      if(stat)
         for(int i=1; i<=newsignal->GetNbinsX(); i++)
            newsignal->SetBinContent(i,oldsignal->GetBinContent(i)+generator->Gaus(0,oldsignal->GetBinError(i)));
      else
         for(int i=1; i<=newsignal->GetNbinsX(); i++)
            newsignal->SetBinContent(i,oldsignal->GetBinContent(i));
      newsignal->Scale(1 + serrf[channel]);
      newsignal->SetDirectory(0);
      TH1 *newbackground = (TH1*)(output->GetBackground()->At(channel));
      TH1 *oldbackground = (TH1*)(input->GetBackground()->At(channel));
      if(stat)
         for(int i=1; i<=newbackground->GetNbinsX(); i++)
            newbackground->SetBinContent(i,oldbackground->GetBinContent(i)+generator->Gaus(0,oldbackground->GetBinError(i)));
      else
         for(int i=1; i<=newbackground->GetNbinsX(); i++)
            newbackground->SetBinContent(i,oldbackground->GetBinContent(i));
      newbackground->Scale(1 + berrf[channel]);
      newbackground->SetDirectory(0);
   }
   delete[] serrf;
   delete[] berrf;
   return 1;
}

TConfidenceLevel *TLimit::ComputeLimit(TH1* s, TH1* b, TH1* d,
                                       Int_t nmc, bool stat,
                                       TRandom * generator)
{
   // Compute limit.

   TLimitDataSource* lds = new TLimitDataSource(s,b,d);
   TConfidenceLevel* out = ComputeLimit(lds,nmc,stat,generator);
   delete lds;
   return out;
}

TConfidenceLevel *TLimit::ComputeLimit(TH1* s, TH1* b, TH1* d,
                                       TVectorD* se, TVectorD* be, TObjArray* l,
                                       Int_t nmc, bool stat,
                                       TRandom * generator)
{
   // Compute limit.

   TLimitDataSource* lds = new TLimitDataSource(s,b,d,se,be,l);
   TConfidenceLevel* out = ComputeLimit(lds,nmc,stat,generator);
   delete lds;
   return out;
}

TConfidenceLevel *TLimit::ComputeLimit(Double_t s, Double_t b, Int_t d,
                                       Int_t nmc,
                                       bool stat,
                                       TRandom * generator)
{
   // Compute limit.

   TH1D* sh = new TH1D("__sh","__sh",1,0,2);
   sh->Fill(1,s);
   TH1D* bh = new TH1D("__bh","__bh",1,0,2);
   bh->Fill(1,b);
   TH1D* dh = new TH1D("__dh","__dh",1,0,2);
   dh->Fill(1,d);
   TLimitDataSource* lds = new TLimitDataSource(sh,bh,dh);
   TConfidenceLevel* out = ComputeLimit(lds,nmc,stat,generator);
   delete lds;
   delete sh;
   delete bh;
   delete dh;
   return out;
}

TConfidenceLevel *TLimit::ComputeLimit(Double_t s, Double_t b, Int_t d,
                                       TVectorD* se, TVectorD* be, TObjArray* l,
                                       Int_t nmc,
                                       bool stat,
                                       TRandom * generator)
{
   // Compute limit.

   TH1D* sh = new TH1D("__sh","__sh",1,0,2);
   sh->Fill(1,s);
   TH1D* bh = new TH1D("__bh","__bh",1,0,2);
   bh->Fill(1,b);
   TH1D* dh = new TH1D("__dh","__dh",1,0,2);
   dh->Fill(1,d);
   TLimitDataSource* lds = new TLimitDataSource(sh,bh,dh,se,be,l);
   TConfidenceLevel* out = ComputeLimit(lds,nmc,stat,generator);
   delete lds;
   delete sh;
   delete bh;
   delete dh;
   return out;
}

Double_t TLimit::LogLikelihood(Double_t s, Double_t b, Double_t b2, Double_t d)
{
   // Compute LogLikelihood (static function)

   return d*TMath::Log((s+b)/b2);
}
 TLimit.cxx:1
 TLimit.cxx:2
 TLimit.cxx:3
 TLimit.cxx:4
 TLimit.cxx:5
 TLimit.cxx:6
 TLimit.cxx:7
 TLimit.cxx:8
 TLimit.cxx:9
 TLimit.cxx:10
 TLimit.cxx:11
 TLimit.cxx:12
 TLimit.cxx:13
 TLimit.cxx:14
 TLimit.cxx:15
 TLimit.cxx:16
 TLimit.cxx:17
 TLimit.cxx:18
 TLimit.cxx:19
 TLimit.cxx:20
 TLimit.cxx:21
 TLimit.cxx:22
 TLimit.cxx:23
 TLimit.cxx:24
 TLimit.cxx:25
 TLimit.cxx:26
 TLimit.cxx:27
 TLimit.cxx:28
 TLimit.cxx:29
 TLimit.cxx:30
 TLimit.cxx:31
 TLimit.cxx:32
 TLimit.cxx:33
 TLimit.cxx:34
 TLimit.cxx:35
 TLimit.cxx:36
 TLimit.cxx:37
 TLimit.cxx:38
 TLimit.cxx:39
 TLimit.cxx:40
 TLimit.cxx:41
 TLimit.cxx:42
 TLimit.cxx:43
 TLimit.cxx:44
 TLimit.cxx:45
 TLimit.cxx:46
 TLimit.cxx:47
 TLimit.cxx:48
 TLimit.cxx:49
 TLimit.cxx:50
 TLimit.cxx:51
 TLimit.cxx:52
 TLimit.cxx:53
 TLimit.cxx:54
 TLimit.cxx:55
 TLimit.cxx:56
 TLimit.cxx:57
 TLimit.cxx:58
 TLimit.cxx:59
 TLimit.cxx:60
 TLimit.cxx:61
 TLimit.cxx:62
 TLimit.cxx:63
 TLimit.cxx:64
 TLimit.cxx:65
 TLimit.cxx:66
 TLimit.cxx:67
 TLimit.cxx:68
 TLimit.cxx:69
 TLimit.cxx:70
 TLimit.cxx:71
 TLimit.cxx:72
 TLimit.cxx:73
 TLimit.cxx:74
 TLimit.cxx:75
 TLimit.cxx:76
 TLimit.cxx:77
 TLimit.cxx:78
 TLimit.cxx:79
 TLimit.cxx:80
 TLimit.cxx:81
 TLimit.cxx:82
 TLimit.cxx:83
 TLimit.cxx:84
 TLimit.cxx:85
 TLimit.cxx:86
 TLimit.cxx:87
 TLimit.cxx:88
 TLimit.cxx:89
 TLimit.cxx:90
 TLimit.cxx:91
 TLimit.cxx:92
 TLimit.cxx:93
 TLimit.cxx:94
 TLimit.cxx:95
 TLimit.cxx:96
 TLimit.cxx:97
 TLimit.cxx:98
 TLimit.cxx:99
 TLimit.cxx:100
 TLimit.cxx:101
 TLimit.cxx:102
 TLimit.cxx:103
 TLimit.cxx:104
 TLimit.cxx:105
 TLimit.cxx:106
 TLimit.cxx:107
 TLimit.cxx:108
 TLimit.cxx:109
 TLimit.cxx:110
 TLimit.cxx:111
 TLimit.cxx:112
 TLimit.cxx:113
 TLimit.cxx:114
 TLimit.cxx:115
 TLimit.cxx:116
 TLimit.cxx:117
 TLimit.cxx:118
 TLimit.cxx:119
 TLimit.cxx:120
 TLimit.cxx:121
 TLimit.cxx:122
 TLimit.cxx:123
 TLimit.cxx:124
 TLimit.cxx:125
 TLimit.cxx:126
 TLimit.cxx:127
 TLimit.cxx:128
 TLimit.cxx:129
 TLimit.cxx:130
 TLimit.cxx:131
 TLimit.cxx:132
 TLimit.cxx:133
 TLimit.cxx:134
 TLimit.cxx:135
 TLimit.cxx:136
 TLimit.cxx:137
 TLimit.cxx:138
 TLimit.cxx:139
 TLimit.cxx:140
 TLimit.cxx:141
 TLimit.cxx:142
 TLimit.cxx:143
 TLimit.cxx:144
 TLimit.cxx:145
 TLimit.cxx:146
 TLimit.cxx:147
 TLimit.cxx:148
 TLimit.cxx:149
 TLimit.cxx:150
 TLimit.cxx:151
 TLimit.cxx:152
 TLimit.cxx:153
 TLimit.cxx:154
 TLimit.cxx:155
 TLimit.cxx:156
 TLimit.cxx:157
 TLimit.cxx:158
 TLimit.cxx:159
 TLimit.cxx:160
 TLimit.cxx:161
 TLimit.cxx:162
 TLimit.cxx:163
 TLimit.cxx:164
 TLimit.cxx:165
 TLimit.cxx:166
 TLimit.cxx:167
 TLimit.cxx:168
 TLimit.cxx:169
 TLimit.cxx:170
 TLimit.cxx:171
 TLimit.cxx:172
 TLimit.cxx:173
 TLimit.cxx:174
 TLimit.cxx:175
 TLimit.cxx:176
 TLimit.cxx:177
 TLimit.cxx:178
 TLimit.cxx:179
 TLimit.cxx:180
 TLimit.cxx:181
 TLimit.cxx:182
 TLimit.cxx:183
 TLimit.cxx:184
 TLimit.cxx:185
 TLimit.cxx:186
 TLimit.cxx:187
 TLimit.cxx:188
 TLimit.cxx:189
 TLimit.cxx:190
 TLimit.cxx:191
 TLimit.cxx:192
 TLimit.cxx:193
 TLimit.cxx:194
 TLimit.cxx:195
 TLimit.cxx:196
 TLimit.cxx:197
 TLimit.cxx:198
 TLimit.cxx:199
 TLimit.cxx:200
 TLimit.cxx:201
 TLimit.cxx:202
 TLimit.cxx:203
 TLimit.cxx:204
 TLimit.cxx:205
 TLimit.cxx:206
 TLimit.cxx:207
 TLimit.cxx:208
 TLimit.cxx:209
 TLimit.cxx:210
 TLimit.cxx:211
 TLimit.cxx:212
 TLimit.cxx:213
 TLimit.cxx:214
 TLimit.cxx:215
 TLimit.cxx:216
 TLimit.cxx:217
 TLimit.cxx:218
 TLimit.cxx:219
 TLimit.cxx:220
 TLimit.cxx:221
 TLimit.cxx:222
 TLimit.cxx:223
 TLimit.cxx:224
 TLimit.cxx:225
 TLimit.cxx:226
 TLimit.cxx:227
 TLimit.cxx:228
 TLimit.cxx:229
 TLimit.cxx:230
 TLimit.cxx:231
 TLimit.cxx:232
 TLimit.cxx:233
 TLimit.cxx:234
 TLimit.cxx:235
 TLimit.cxx:236
 TLimit.cxx:237
 TLimit.cxx:238
 TLimit.cxx:239
 TLimit.cxx:240
 TLimit.cxx:241
 TLimit.cxx:242
 TLimit.cxx:243
 TLimit.cxx:244
 TLimit.cxx:245
 TLimit.cxx:246
 TLimit.cxx:247
 TLimit.cxx:248
 TLimit.cxx:249
 TLimit.cxx:250
 TLimit.cxx:251
 TLimit.cxx:252
 TLimit.cxx:253
 TLimit.cxx:254
 TLimit.cxx:255
 TLimit.cxx:256
 TLimit.cxx:257
 TLimit.cxx:258
 TLimit.cxx:259
 TLimit.cxx:260
 TLimit.cxx:261
 TLimit.cxx:262
 TLimit.cxx:263
 TLimit.cxx:264
 TLimit.cxx:265
 TLimit.cxx:266
 TLimit.cxx:267
 TLimit.cxx:268
 TLimit.cxx:269
 TLimit.cxx:270
 TLimit.cxx:271
 TLimit.cxx:272
 TLimit.cxx:273
 TLimit.cxx:274
 TLimit.cxx:275
 TLimit.cxx:276
 TLimit.cxx:277
 TLimit.cxx:278
 TLimit.cxx:279
 TLimit.cxx:280
 TLimit.cxx:281
 TLimit.cxx:282
 TLimit.cxx:283
 TLimit.cxx:284
 TLimit.cxx:285
 TLimit.cxx:286
 TLimit.cxx:287
 TLimit.cxx:288
 TLimit.cxx:289
 TLimit.cxx:290
 TLimit.cxx:291
 TLimit.cxx:292
 TLimit.cxx:293
 TLimit.cxx:294
 TLimit.cxx:295
 TLimit.cxx:296
 TLimit.cxx:297
 TLimit.cxx:298
 TLimit.cxx:299
 TLimit.cxx:300
 TLimit.cxx:301
 TLimit.cxx:302
 TLimit.cxx:303
 TLimit.cxx:304
 TLimit.cxx:305
 TLimit.cxx:306
 TLimit.cxx:307
 TLimit.cxx:308
 TLimit.cxx:309
 TLimit.cxx:310
 TLimit.cxx:311
 TLimit.cxx:312
 TLimit.cxx:313
 TLimit.cxx:314
 TLimit.cxx:315
 TLimit.cxx:316
 TLimit.cxx:317
 TLimit.cxx:318
 TLimit.cxx:319
 TLimit.cxx:320
 TLimit.cxx:321
 TLimit.cxx:322
 TLimit.cxx:323
 TLimit.cxx:324
 TLimit.cxx:325
 TLimit.cxx:326
 TLimit.cxx:327
 TLimit.cxx:328
 TLimit.cxx:329
 TLimit.cxx:330
 TLimit.cxx:331
 TLimit.cxx:332
 TLimit.cxx:333
 TLimit.cxx:334
 TLimit.cxx:335
 TLimit.cxx:336
 TLimit.cxx:337
 TLimit.cxx:338
 TLimit.cxx:339
 TLimit.cxx:340
 TLimit.cxx:341
 TLimit.cxx:342
 TLimit.cxx:343
 TLimit.cxx:344
 TLimit.cxx:345
 TLimit.cxx:346
 TLimit.cxx:347
 TLimit.cxx:348
 TLimit.cxx:349
 TLimit.cxx:350
 TLimit.cxx:351
 TLimit.cxx:352
 TLimit.cxx:353
 TLimit.cxx:354
 TLimit.cxx:355
 TLimit.cxx:356
 TLimit.cxx:357
 TLimit.cxx:358
 TLimit.cxx:359
 TLimit.cxx:360
 TLimit.cxx:361
 TLimit.cxx:362
 TLimit.cxx:363
 TLimit.cxx:364
 TLimit.cxx:365
 TLimit.cxx:366
 TLimit.cxx:367
 TLimit.cxx:368
 TLimit.cxx:369
 TLimit.cxx:370
 TLimit.cxx:371
 TLimit.cxx:372
 TLimit.cxx:373
 TLimit.cxx:374
 TLimit.cxx:375
 TLimit.cxx:376
 TLimit.cxx:377
 TLimit.cxx:378
 TLimit.cxx:379
 TLimit.cxx:380
 TLimit.cxx:381
 TLimit.cxx:382
 TLimit.cxx:383
 TLimit.cxx:384
 TLimit.cxx:385
 TLimit.cxx:386
 TLimit.cxx:387
 TLimit.cxx:388
 TLimit.cxx:389
 TLimit.cxx:390
 TLimit.cxx:391
 TLimit.cxx:392
 TLimit.cxx:393
 TLimit.cxx:394
 TLimit.cxx:395
 TLimit.cxx:396
 TLimit.cxx:397
 TLimit.cxx:398
 TLimit.cxx:399
 TLimit.cxx:400
 TLimit.cxx:401
 TLimit.cxx:402
 TLimit.cxx:403
 TLimit.cxx:404
 TLimit.cxx:405
 TLimit.cxx:406
 TLimit.cxx:407
 TLimit.cxx:408
 TLimit.cxx:409
 TLimit.cxx:410
 TLimit.cxx:411
 TLimit.cxx:412
 TLimit.cxx:413
 TLimit.cxx:414
 TLimit.cxx:415
 TLimit.cxx:416
 TLimit.cxx:417
 TLimit.cxx:418
 TLimit.cxx:419