C++ ToyMCSampler::SetTestStatistic方法代码示例

//.........这里部分代码省略.........
   // set the test statistic 
   TestStatistic * testStat = 0;
   if (testStatType == 0) testStat = &slrts;
   if (testStatType == 1 || testStatType == 11) testStat = &ropl;
   if (testStatType == 2 || testStatType == 3 || testStatType == 4) testStat = &profll;
   if (testStatType == 5) testStat = &maxll;
   if (testStatType == 6) testStat = &nevtts;

   if (testStat == 0) { 
      Error("StandardHypoTestInvDemo","Invalid - test statistic type = %d supported values are only :\n\t\t\t 0 (SLR) , 1 (Tevatron) , 2 (PLR), 3 (PLR1), 4(MLE)",testStatType);
      return 0;
   }
  
  
   ToyMCSampler *toymcs = (ToyMCSampler*)hc->GetTestStatSampler();
   if (toymcs && (type == 0 || type == 1) ) { 
      // look if pdf is number counting or extended
      if (sbModel->GetPdf()->canBeExtended() ) { 
         if (useNumberCounting)   Warning("StandardHypoTestInvDemo","Pdf is extended: but number counting flag is set: ignore it ");
      }
      else { 
         // for not extended pdf
         if (!useNumberCounting  )  { 
            int nEvents = data->numEntries();
            Info("StandardHypoTestInvDemo","Pdf is not extended: number of events to generate taken  from observed data set is %d",nEvents);
            toymcs->SetNEventsPerToy(nEvents);
         }
         else {
            Info("StandardHypoTestInvDemo","using a number counting pdf");
            toymcs->SetNEventsPerToy(1);
         }
      }

      toymcs->SetTestStatistic(testStat);
    
      if (data->isWeighted() && !mGenerateBinned) { 
         Info("StandardHypoTestInvDemo","Data set is weighted, nentries = %d and sum of weights = %8.1f but toy generation is unbinned - it would be faster to set mGenerateBinned to true\n",data->numEntries(), data->sumEntries());
      }
      toymcs->SetGenerateBinned(mGenerateBinned);
  
      toymcs->SetUseMultiGen(mOptimize);
    
      if (mGenerateBinned &&  sbModel->GetObservables()->getSize() > 2) { 
         Warning("StandardHypoTestInvDemo","generate binned is activated but the number of ovservable is %d. Too much memory could be needed for allocating all the bins",sbModel->GetObservables()->getSize() );
      }

      // set the random seed if needed
      if (mRandomSeed >= 0) RooRandom::randomGenerator()->SetSeed(mRandomSeed); 
    
   }
  
   // specify if need to re-use same toys
   if (reuseAltToys) {
      hc->UseSameAltToys();
   }
  
   if (type == 1) { 
      HybridCalculator *hhc = dynamic_cast<HybridCalculator*> (hc);
      assert(hhc);
    
      hhc->SetToys(ntoys,ntoys/mNToysRatio); // can use less ntoys for b hypothesis 
    
      // remove global observables from ModelConfig (this is probably not needed anymore in 5.32)
      bModel->SetGlobalObservables(RooArgSet() );
      sbModel->SetGlobalObservables(RooArgSet() );
    

//void RunToyScan5(TString fileName, double startVal, double stopVal, TString outFile) {
void frequentist(TString fileName) {
  cout << "Starting frequentist " << time(NULL) << endl;
  double startVal = 0;
  double stopVal = 200;
  TString outFile = "";

  int nToys = 1 ;
  int nscanpoints = 2 ;

  /*
  gROOT->LoadMacro("RooBetaPdf.cxx+") ;
  gROOT->LoadMacro("RooRatio.cxx+") ;
  gROOT->LoadMacro("RooPosDefCorrGauss.cxx+") ;
  */

  // get relevant objects out of the "ws" file

  TFile *file = TFile::Open(fileName);
  if(!file){
    cout <<"file not found" << endl;
    return;
  } 

  RooWorkspace* w = (RooWorkspace*) file->Get("workspace");
  if(!w){
    cout <<"workspace not found" << endl;
    return;
  }

  ModelConfig* mc = (ModelConfig*) w->obj("S+B_model");
  RooAbsData* data = w->data("data");

  if( !data || !mc ){
    w->Print();
    cout << "data or ModelConfig was not found" <<endl;
    return;
  }

  RooRealVar* myPOI = (RooRealVar*) mc->GetParametersOfInterest()->first();
  myPOI->setRange(0, 1000.);

  ModelConfig* bModel = (ModelConfig*) w->obj("B_model");
  ModelConfig* sbModel = (ModelConfig*) w->obj("S+B_model");

  ProfileLikelihoodTestStat profll(*sbModel->GetPdf());
  profll.SetPrintLevel(2);
  profll.SetOneSided(1);
  TestStatistic * testStat = &profll;

  HypoTestCalculatorGeneric *  hc = 0;
  hc = new FrequentistCalculator(*data, *bModel, *sbModel);
  
  ToyMCSampler *toymcs = (ToyMCSampler*)hc->GetTestStatSampler();
  toymcs->SetMaxToys(10000);
  toymcs->SetNEventsPerToy(1);
  toymcs->SetTestStatistic(testStat);


  ((FrequentistCalculator *)hc)->SetToys(nToys,nToys);
  
  HypoTestInverter calc(*hc);
  calc.SetConfidenceLevel(0.95);
  calc.UseCLs(true);
  //calc.SetVerbose(true);
  calc.SetVerbose(2);

  cout << "About to set fixed scan " << time(NULL) << endl;
  calc.SetFixedScan(nscanpoints,startVal,stopVal);
  cout << "About to do inverter " << time(NULL) << endl;
  HypoTestInverterResult * res_toysCLs_calculator = calc.GetInterval();

  cout << "CLs = " << res_toysCLs_calculator->UpperLimit() 
	    << "   CLs_exp = " << res_toysCLs_calculator->GetExpectedUpperLimit(0) 
	    << "   CLs_exp(-1s) = " << res_toysCLs_calculator->GetExpectedUpperLimit(-1) 
	    << "   CLs_exp(+1s) = " << res_toysCLs_calculator->GetExpectedUpperLimit(1) << endl ;

  /*
  // dump results string to output file
  ofstream outStream ;
  outStream.open(outFile,ios::app) ;
  
  outStream << "CLs = " << res_toysCLs_calculator->UpperLimit() 
	    << "   CLs_exp = " << res_toysCLs_calculator->GetExpectedUpperLimit(0) 
	    << "   CLs_exp(-1s) = " << res_toysCLs_calculator->GetExpectedUpperLimit(-1) 
	    << "   CLs_exp(+1s) = " << res_toysCLs_calculator->GetExpectedUpperLimit(1) << endl ;
  
  outStream.close() ;
  */


  cout << "End of frequentist " << time(NULL) << endl;
  return ;

}

//.........这里部分代码省略.........
      
         ((HybridCalculator*)hypoCalc)->ForcePriorNuisanceAlt(*nuisPdf);
         ((HybridCalculator*)hypoCalc)->ForcePriorNuisanceNull(*nuisPdf);
   }

   // hypoCalc->ForcePriorNuisanceAlt(*sbModel->GetPriorPdf());
   // hypoCalc->ForcePriorNuisanceNull(*bModel->GetPriorPdf());

   ToyMCSampler * sampler = (ToyMCSampler *)hypoCalc->GetTestStatSampler();

   if (sampler && (calcType == 0 || calcType == 1) ) { 

      // look if pdf is number counting or extended
      if (sbModel->GetPdf()->canBeExtended() ) { 
         if (useNC)   Warning("StandardHypoTestDemo","Pdf is extended: but number counting flag is set: ignore it ");
      }
      else {
         // for not extended pdf
         if (!useNC)  { 
            int nEvents = data->numEntries();
            Info("StandardHypoTestDemo","Pdf is not extended: number of events to generate taken  from observed data set is %d",nEvents);
            sampler->SetNEventsPerToy(nEvents);
         }
         else {
            Info("StandardHypoTestDemo","using a number counting pdf");
            sampler->SetNEventsPerToy(1);
         }
      }
      
      if (data->isWeighted() && !generateBinned) { 
         Info("StandardHypoTestDemo","Data set is weighted, nentries = %d and sum of weights = %8.1f but toy generation is unbinned - it would be faster to set generateBinned to true\n",data->numEntries(), data->sumEntries());
      }
      if (generateBinned)  sampler->SetGenerateBinned(generateBinned);


      // set the test statistic
      if (testStatType == 0) sampler->SetTestStatistic(slrts); 
      if (testStatType == 1) sampler->SetTestStatistic(ropl); 
      if (testStatType == 2 || testStatType == 3) sampler->SetTestStatistic(profll); 

   }
   
   HypoTestResult *  htr = hypoCalc->GetHypoTest();
   htr->SetPValueIsRightTail(true);
   htr->SetBackgroundAsAlt(false);
   htr->Print(); // how to get meaningfull CLs at this point?

   delete sampler;
   delete slrts; 
   delete ropl; 
   delete profll;

   if (calcType != 2) {
      HypoTestPlot * plot = new HypoTestPlot(*htr,100);
      plot->SetLogYaxis(true);
      plot->Draw();
   }
   else { 
      std::cout << "Asymptotic results " << std::endl;
      
   }

   // look at expected significances 
   // found median of S+B distribution
   if (calcType != 2) { 

      SamplingDistribution * altDist = htr->GetAltDistribution();   
      HypoTestResult htExp("Expected Result");
      htExp.Append(htr);
      // find quantiles in alt (S+B) distribution 
      double p[5];
      double q[5];
      for (int i = 0; i < 5; ++i) { 
         double sig = -2  + i;
         p[i] = ROOT::Math::normal_cdf(sig,1);
      }
      std::vector<double> values = altDist->GetSamplingDistribution();
      TMath::Quantiles( values.size(), 5, &values[0], q, p, false);  

      for (int i = 0; i < 5; ++i) { 
         htExp.SetTestStatisticData( q[i] );
         double sig = -2  + i;      
         std::cout << " Expected p -value and significance at " << sig << " sigma = " 
                   << htExp.NullPValue() << " significance " << htExp.Significance() << " sigma " << std::endl; 
         
      }
   }
   else { 
      // case of asymptotic calculator 
      for (int i = 0; i < 5; ++i) { 
         double sig = -2  + i;      
         // sigma is inverted here 
         double pval = AsymptoticCalculator::GetExpectedPValues( htr->NullPValue(), htr->AlternatePValue(), -sig, false);
         std::cout << " Expected p -value and significance at " << sig << " sigma = " 
                   << pval << " significance " << ROOT::Math::normal_quantile_c(pval,1) << " sigma " << std::endl; 
         
      }
   }

}

//.........这里部分代码省略.........
         }         
      }
   }


   SimpleLikelihoodRatioTestStat slrts(*sbModel->GetPdf(),*bModel->GetPdf());
   if (sbModel->GetSnapshot()) slrts.SetNullParameters(*sbModel->GetSnapshot());
   if (bModel->GetSnapshot()) slrts.SetAltParameters(*bModel->GetSnapshot());

   // ratio of profile likelihood - need to pass snapshot for the alt
   RatioOfProfiledLikelihoodsTestStat 
      ropl(*sbModel->GetPdf(), *bModel->GetPdf(), bModel->GetSnapshot());
   ropl.SetSubtractMLE(false);
   
   //MyProfileLikelihoodTestStat profll(*sbModel->GetPdf());
   ProfileLikelihoodTestStat profll(*sbModel->GetPdf());
   if (testStatType == 3) profll.SetOneSided(1);
   if (optimize) profll.SetReuseNLL(true);

   TestStatistic * testStat = &slrts;
   if (testStatType == 1) testStat = &ropl;
   if (testStatType == 2 || testStatType == 3) testStat = &profll;
  
   
   HypoTestCalculatorGeneric *  hc = 0;
   if (type == 0) hc = new FrequentistCalculator(*data, *bModel, *sbModel);
   else hc = new HybridCalculator(*data, *bModel, *sbModel);

   ToyMCSampler *toymcs = (ToyMCSampler*)hc->GetTestStatSampler();
   //=== DEBUG
   ///// toymcs->SetWS( w ) ;
   //=== DEBUG
   toymcs->SetNEventsPerToy(1);
   toymcs->SetTestStatistic(testStat);
   if (optimize) toymcs->SetUseMultiGen(true);


   if (type == 1) { 
      HybridCalculator *hhc = (HybridCalculator*) hc;
      hhc->SetToys(ntoys,ntoys); 

      // check for nuisance prior pdf 
      if (bModel->GetPriorPdf() && sbModel->GetPriorPdf() ) {
         hhc->ForcePriorNuisanceAlt(*bModel->GetPriorPdf());
         hhc->ForcePriorNuisanceNull(*sbModel->GetPriorPdf());
      }
      else {
         if (bModel->GetNuisanceParameters() || sbModel->GetNuisanceParameters() ) {
            Error("RA2bHypoTestInvDemo","Cannnot run Hybrid calculator because no prior on the nuisance parameter is specified");
            return 0;
         }
      }
   } 
   else 
      ((FrequentistCalculator*) hc)->SetToys(ntoys,ntoys); 

   // Get the result
   RooMsgService::instance().getStream(1).removeTopic(RooFit::NumIntegration);


   TStopwatch tw; tw.Start(); 
   const RooArgSet * poiSet = sbModel->GetParametersOfInterest();
   RooRealVar *poi = (RooRealVar*)poiSet->first();

   // fit the data first
   sbModel->GetPdf()->fitTo(*data);

void HypoTestInvDemo(const char * fileName ="GausModel_b.root",
                     const char * wsName = "w",
                     const char * modelSBName = "model_sb",
                     const char * modelBName = "model_b",
                     const char * dataName = "data_obs",                  
                     int type = 0,  // calculator type 
                     int testStatType = 0, // test stat type
                     int npoints = 10,   
                     int ntoys=1000,
                     bool useCls = true )
{ 
   /*
    type = 0 Freq calculator 
    type = 1 Hybrid 

    testStatType = 0 LEP
                 = 1 Tevatron 
                 = 2 PL


   */

   if (fileName==0) { 
      std::cout << "give input filename " << std::endl;
      return; 
   }
   TFile * file = new TFile(fileName); 

   RooWorkspace * w = dynamic_cast<RooWorkspace*>( file->Get(wsName) );
   if (!w) {      
      return; 
   }
   w->Print();


   RooAbsData * data = w->data(dataName); 
   if (!data) { 
      Error("HypoTestDemo","Not existing data %s",dataName);
   }

   
   // get models from WS
  // get the modelConfig out of the file
  ModelConfig* bModel = (ModelConfig*) w->obj(modelBName);
  ModelConfig* sbModel = (ModelConfig*) w->obj(modelSBName);


   SimpleLikelihoodRatioTestStat slrts(*bModel->GetPdf(),*sbModel->GetPdf());
   slrts.SetNullParameters(*bModel->GetSnapshot());
   slrts.SetAltParameters(*sbModel->GetSnapshot());

   RatioOfProfiledLikelihoodsTestStat 
   ropl(*bModel->GetPdf(), *sbModel->GetPdf(), sbModel->GetSnapshot());
   ropl.SetSubtractMLE(false);
   
   ProfileLikelihoodTestStat profll(*sbModel->GetPdf());
   profll.SetOneSided(0);

   TestStatistic * testStat = &slrts;
   if (testStatType == 1) testStat = &ropl;
   if (testStatType == 2) testStat = &profll;
  
   
   HypoTestCalculatorGeneric *  hc = 0;
   if (type == 0) hc = new FrequentistCalculator(*data, *sbModel, *bModel);
   else new HybridCalculator(*data, *sbModel, *bModel);

   ToyMCSampler *toymcs = (ToyMCSampler*)hc->GetTestStatSampler();
   //toymcs->SetNEventsPerToy(1);
   toymcs->SetTestStatistic(testStat);


    if (type == 1) { 
      HybridCalculator *hhc = (HybridCalculator*) hc;
      hhc->SetToys(ntoys,ntoys); 
      // hhc->ForcePriorNuisanceAlt(*pdfNuis);
      // hhc->ForcePriorNuisanceNull(*pdfNuis);
   } 
   else 
      ((FrequentistCalculator*) hc)->SetToys(ntoys,ntoys); 

  // Get the result
   RooMsgService::instance().getStream(1).removeTopic(RooFit::NumIntegration);


   TStopwatch tw; tw.Start(); 
   const RooArgSet * poiSet = sbModel->GetParametersOfInterest();
   RooRealVar *poi = (RooRealVar*)poiSet->first();

   // fit the data first
   sbModel->GetPdf()->fitTo(*data);
   double poihat  = poi->getVal();
   //poi->setVal(30);
   //poi->setError(10);


   HypoTestInverter calc(*hc);
   // GENA: for two-sided interval
   //calc.SetConfidenceLevel(0.95);
   // GENA: for 95% upper limit
//.........这里部分代码省略.........