C++ TH1F::KolmogorovTest方法代码示例

/*************************************************************************************
 * getKS: Searches through the histograms in the plotter output, adds the MC together
 *        for each field, and compares the MC with the Data histogram using a KS test
 *  input:  the main() arguments array
 *  output: writes to stdout the (human-readable) KS statistics of pairs of histograms
 *
 *  Structure-wise: this is fine, can be implemented into class easily.
 ***********************************/
void getKS(const char* argv[]) {
  //open the input TFile
  TFile *f = new TFile(argv[2]);
  f->cd();

  //get the filesystem information from the file
  TList *alokDirs = (TList*) f->GetListOfKeys();

  //loop through the directories in the input file
  for(int idir=0; alokDirs->At(idir-1) != alokDirs->Last(); idir++) {
    TDirectory *cDir  = (TDirectory*) f->Get(alokDirs->At(idir)->GetName());
    TList *alokHistos = (TList*)      cDir->GetListOfKeys();

    // create the MC histogram and start collecting relevant MC histograms
    // from the current directory
    TList *aloh   = new TList;
    // loop through keys (histograms) in current directory
    for(int ihisto=0; alokHistos->At(ihisto) != alokHistos->Last(); ihisto++) {
      if(TString(alokHistos->At(ihisto)->GetName()).Contains("MC8TeV")) {
        TH1F *cHisto = (TH1F*) cDir->Get(alokHistos->At(ihisto)->GetName());
        aloh->Add(cHisto);
      }
    }
 
    //merge the data histograms into one histogram
    TH1F *MCHisto = (TH1F*) (aloh->Last())->Clone(TString(cDir->GetName()) + TString("MCHisto"));
    aloh->RemoveLast();
    MCHisto->Merge(aloh);

    cout<<"-------------------- "<<cDir->GetName()<<" -----------------------"<<endl;
    //now create the data histogram and run the KS test
    TH1F *DataHisto = (TH1F*) cDir->Get(alokHistos->Last()->GetName());
    cout<<"  ---> KS Test: "<<cDir->GetName()<<" has probability "<<MCHisto->KolmogorovTest(DataHisto, "D")<<"\n"<<endl;
  }
}

double Chi2F(const double *xx ){
 const Double_t scale = xx[0];
 fileInMC->cd();
 TTree* MyTreeMC = (TTree*) fileInMC->Get(treeNameMC.c_str());
 outFile->cd();
 TString nameDATA = Form("hDATA_%d_%d_%.5f",Data_or_MC,nIter,ScaleTrue);
 TH1F* hDATA = (TH1F*) outFile->Get(nameDATA);
 TString NameMC = Form("hMC_Chi2_%.5f",scale);
 
 //std::cerr << " NameMC = "  << NameMC.Data() << " => " << scale << std::endl;
 
 TH1F* hMC;
 if (!gROOT->FindObject(NameMC.Data())){
  hMC = new TH1F(NameMC,NameMC,numBINS,minBINS,maxBINS);
  hMC->Reset();

//  std::cerr << " AdditionalCut.Data() = " << AdditionalCut.Data() << std::endl;
  
  MyTreeMC->SetEntryList(0); 
  MyTreeMC->Draw(">> myListMCTot",(AdditionalCut + Form(" && (ET * (1+(%f)))>%f",scale,minET)).Data(),"entrylist");
  TEntryList *mylistMCTot = (TEntryList*)gDirectory->Get("myListMCTot");
  MyTreeMC->SetEntryList(mylistMCTot);
  
  TString DrawMC = Form("(%s * (1+(%f)))>>%s",variableName.c_str(),scale,NameMC.Data());
//  std::cerr << " DrawMC = " << DrawMC.Data() << std::endl;
  MyTreeMC->Draw(DrawMC);
//   MyTreeMC->Draw(DrawMC,(AdditionalCut+Form("&& (ET * (1+(%f)))>%f",scale,minET)).Data());
  hMC->Sumw2();
  hMC->Scale(hDATA->GetEffectiveEntries()/hMC->GetEffectiveEntries());
  outFile->cd();
  hMC->Write();
 }
 else {
//   std::cerr << " KM old "  << NameMC.Data() << std::endl;
  hMC = (TH1F*) outFile->Get(NameMC.Data());
 }
 
 outFile->cd();
//  hDATA.Write();
//  hMC.Write();
 
 double result = hMC->KolmogorovTest(hDATA,"M");
// double result = - hMC->KolmogorovTest(hDATA,"X");
//  double result = hMC->Chi2Test(&hDATA,"CHI2/NDF");
//=========> E' QUESTO! ==> double result = hMC->Chi2Test(hDATA,"CHI2");
 //double result = - hMC.Chi2Test(&hDATA,""); ///==== http://root.cern.ch/root/html/TH1.html#TH1:Chi2Test
 return result;
}

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

   // start with the main "subevent", as it contains most of the info
  // and always exists.
  WCSimRootTrigger* wcsimrootevent;

  TH1F *h2 = new TH1F("PMT Hits 2", "Digitized Hits", 500, 0, 3000);
  TH1F *time2 = new TH1F("Average time 2", "Average time", 600, 900, 2000);
  TH1F *pe2 = new TH1F("Q/# Digitzed PMT 2", "Q/# Digitzed PMT", 200, 0, 5);
 
  
  // Now loop over events
  for (int ev=0; ev<nevent; ev++)
  {
    // Read the event from the tree into the WCSimRootEvent instance
    wcsimT2->GetEntry(ev);      
    wcsimrootevent = wcsimrootsuperevent->GetTrigger(0);
  
    if(verbose){
      printf("********************************************************");
      printf("Evt, date %d %d\n", wcsimrootevent->GetHeader()->GetEvtNum(),
	     wcsimrootevent->GetHeader()->GetDate());
      printf("Mode %d\n", wcsimrootevent->GetMode());
      printf("Number of subevents %d\n",
	     wcsimrootsuperevent->GetNumberOfSubEvents());
      
      printf("Vtxvol %d\n", wcsimrootevent->GetVtxvol());
      printf("Vtx %f %f %f\n", wcsimrootevent->GetVtx(0),
	     wcsimrootevent->GetVtx(1),wcsimrootevent->GetVtx(2));
    }
    
    for (int index = 0 ; index < wcsimrootsuperevent->GetNumberOfEvents(); index++) 
      { 
	int ncherenkovdigihits = wcsimrootevent->GetNcherenkovdigihits();
	h2->Fill(ncherenkovdigihits);

	
	float totalq = 0.;
	float totalt = 0.;
	// Loop through elements in the TClonesArray of WCSimRootCherenkovHits
	for (int i=0; i< ncherenkovdigihits; i++)
	  {
	    TObject *Digi = (wcsimrootevent->GetCherenkovDigiHits())->At(i);
	    WCSimRootCherenkovDigiHit *wcsimrootcherenkovdigihit = 
	      dynamic_cast<WCSimRootCherenkovDigiHit*>(Digi);
	    
	    int tubeNumber     =  (wcsimrootcherenkovdigihit->GetT(),wcsimrootcherenkovdigihit->GetTubeId());
	    float q = wcsimrootcherenkovdigihit->GetQ();
	    float t = wcsimrootcherenkovdigihit->GetT();
	    totalq+=q;
	    totalt+=t;
	  }
	float av_time = totalt/ncherenkovdigihits;
	float av_q = totalq/ncherenkovdigihits;
      }
    pe2->Fill(av_q);  
    time2->Fill(av_time);
    // reinitialize super event between loops.
    wcsimrootsuperevent->ReInitialize();
  }// End of loop over events

	Double_t ks_hits = h1->KolmogorovTest(h2);
	Double_t ks_charge = pe->KolmogorovTest(pe2);
	Double_t ks_time = time->KolmogorovTest(time2);
	cout << "***********************************************************" << endl;
	cout << "ks test for # of digitized hits: " << ks_hits << endl;
	cout << "ks test for average charge: " << ks_charge << endl;
	cout << "ks test for average time: " << ks_time << endl;

  //  TCanvas c1("c1"); 
  float win_scale = 0.75;
  int n_wide(2);
  int n_high(2);
  TCanvas* c1 = new TCanvas("c1", "Test Plots", 500*n_wide*win_scale, 500*n_high*win_scale);
  c1->Draw();
  c1->Divide(2,2);
  c1->cd(1); 
  h2->SetLineColor(kRed);
  h1->Draw();
  c1->cd(1); h2->Draw("SAME");

 TLegend *leg = new TLegend(0.2,0.7,0.55,0.85, "");
 leg->SetFillColor(0);
 leg->SetBorderSize(0);
 leg->AddEntry(h1,filename, "l");
 leg->AddEntry(h2,filename2, "l");
 leg->Draw();
 
 c1->cd(2);
 pe->GetXaxis()->SetTitle("Total Charge / # digitized hits");
 pe->Draw();
 pe2->SetLineColor(kRed);
 c1->cd(2); pe2->Draw("SAME");
 c1->cd(3); 
 time->GetXaxis()->SetTitle("Total Time / # digitized hits (ns)");  
 time->Draw();
 time2->SetLineColor(kRed);
 c1->cd(3); time2->Draw("SAME");
  

}

//.........这里部分代码省略.........
    double NewMin = min;
    double NewMax = max;

    int rebinning = Nbins; 

    if (RMS>0 && hd->GetType() ) 
      {
	dNdX = 100. / ( 10 * RMS);
	NewMin = Mean - 10 * RMS;
	NewMax = Mean + 10 * RMS;
      }
    
    if ((dX * dNdX)>0  && hd->GetType() ) 
      rebinning = (int)(double(Nbins) / (dX * dNdX));
    
    if ( rebinning > 1 && hd->GetType() ) 
      {
	href->Rebin(rebinning);
	hnew->Rebin(rebinning);
      }

    if ( hd->GetType() == 1 )
      { 
	href->GetXaxis()->SetRangeUser(0.0, NewMax);
	hnew->GetXaxis()->SetRangeUser(0.0, NewMax);
      }
    else if ( hd->GetType() == 2 )
      {
	href->GetXaxis()->SetRangeUser(NewMin, NewMax);
	hnew->GetXaxis()->SetRangeUser(NewMin, NewMax);
      }

    // perform statistical tests
    double ks_score = hnew->KolmogorovTest(href,"D");
    double chi2_score = hnew->Chi2Test(href, "p");
    //double result = KS_TEST ? ks_score : chi2_score;
    double result = (ks_score>chi2_score) ? ks_score : chi2_score;
    
      href->SetNormFactor(Nevents_new);     	
      hnew->SetNormFactor(Nevents_new);
    //hnew->SetNormFactor(1);

    // ensure that the peaks of both histograms will be shown by making a dummy histogram
    float Nentries_ref = href->GetEntries();
    float Nentries_new = hnew->GetEntries();
    float XaxisMin_ref = 0, XaxisMax_ref = 0, YaxisMin_ref = 0, YaxisMax_ref = 0;
    float XaxisMin_new = 0, XaxisMax_new = 0, YaxisMin_new = 0, YaxisMax_new = 0;
    if (Nentries_ref>0) YaxisMax_ref = (href->GetMaximum()+TMath::Sqrt(href->GetMaximum()))*(Nentries_new/Nentries_ref);
    if (Nentries_new>0) YaxisMax_new = (hnew->GetMaximum()+TMath::Sqrt(hnew->GetMaximum()));

    XaxisMin_ref = href->GetXaxis()->GetXmin()>NewMin  ? href->GetXaxis()->GetXmin() : NewMin;
    XaxisMax_ref = href->GetXaxis()->GetXmax()<=NewMax ? href->GetXaxis()->GetXmax() : NewMax;
    YaxisMax_ref = (YaxisMax_ref>=YaxisMax_new) ? YaxisMax_ref : YaxisMax_new;

    if (TMath::Abs(XaxisMin_ref - XaxisMax_ref)<1E-6)
      {
	XaxisMin_ref = 0;
	XaxisMax_ref = 1;
      }
    
    TH1F *hdumb = new TH1F("hdumb","", rebinning, XaxisMin_ref, XaxisMax_ref);
    hdumb->SetMinimum(1E-1); //--For Rick
    hdumb->SetMaximum(1.05*YaxisMax_ref);
    //    if (href->GetMaximum() < hnew->GetMaximum())
    //  href->SetAxisRange(0, 1.1 * hnew->GetMaximum(), "Y");
        

//.........这里部分代码省略.........
  hsig->SetYTitle("Iso (GeV)");

  hsig->SetLineColor(2);
  TH1F* hspike = new TH1F("hspike","",120,-1,11);
  hspike->SetLineColor(4);

  TCanvas* c1 = new TCanvas("c1","",500,500);

  for(int ieta=0; ieta< 1; ieta++){
    for(int ipt=0; ipt < nPtBin; ipt++){


      hsig->Reset();
      hspike->Reset();

      // first obtain fit for the random cone corrected values
      hfit_sig  = (TH1F*)hdata_S[ieta][ipt]->Clone();
      hfit_sig  -> SetName("hfit_sig");
      hfit_sig  -> Rebin(4);

      double sigPar[20] = {hfit_sig->GetMaximum(), 1., 0.5, 0.3,
			   1.,-.3,-1., 0.01, 0.5, 0.01, 1., 1.};

      TF1 *fsig = new TF1("fsig", exp_conv, -1., 11., 12);
      fsig->SetParameters(sigPar);
      fsig->SetNpx(2500);
      fsig->SetLineColor(2);
      hfit_sig->Fit(fsig,"","",-1,5.0);

      c1->Print(Form("hfit_sig_pt%d.gif",(int)fBinsPt[ipt]));
      
      fsig->SetParameter(0,2);
      fsig->SetParameter(1,fsig->GetParameter(1)*8.39614e-01/6.83080e-01);//correction from RC
      fsig->SetParameter(2,4.83182e-01);
      fsig->SetParameter(3,fsig->GetParameter(3)*2.33769e-01/2.26323e-01);
    
      cout << "fsig Integral = " << fsig->Integral(-1,11) << endl;
      for(int i=0;i<4;i++)
	cout << "fsig par " << i << "= " << fsig->GetParameter(i) << endl;

      // second obtain fit for the spikes
      hfit_spike  = (TH1F*)hTemplate_S[ieta][ipt]->Clone();
      hfit_spike  -> SetName("hfit_spike");
      hfit_spike  -> Rebin(4);

      TF1 *fspike = new TF1("fspike", exp_conv, -1., 11., 12);
      fspike->SetParameters(sigPar);
      fspike->SetParameter(0,hfit_spike->GetMaximum());
      fspike->SetLineColor(4);
      fspike->SetNpx(2500);
      hfit_spike->Fit(fspike,"","",-1,5.0);

      c1->Print(Form("hfit_spike_pt%d.gif",(int)fBinsPt[ipt]));
      
      fspike->SetParameter(0,2.0);
      cout << "fspike Integral = " << fspike->Integral(-1,11) << endl;
      float scale = (float)fsig->Integral(-1,11)/(float)fspike->Integral(-1,11);
      fspike->SetParameter(0,2.0*scale);
      cout << "now fspike Integral = " << fspike->Integral(-1,11) << endl;

      for(int i=0;i<4;i++)
	cout << "fspike par " << i << "= " << fspike->GetParameter(i) << endl;

    
      hsig->FillRandom("fsig",10000);
      hspike->FillRandom("fspike",10000);

      cout << "KS test for sideband data and sideband MC = " << 
	hsig->KolmogorovTest(hspike,"X") << endl;
      cout << "KS test 2 for sideband data and sideband MC = " << 
	hsig->KolmogorovTest(hspike) << endl;
     
      hsig->Reset();
      hsig->SetMaximum(1.2);
      hsig->Draw();
//       hspike->Draw("same");
      fsig->Draw("same");
      fspike->Draw("same");

      TLegend* leg = new TLegend(0.35,0.70,0.55,0.85);
      leg->SetHeader(Form("p_{T}[%d,%d]",(int)fBinsPt[ipt],(int)fBinsPt[ipt+1]));
      leg->SetFillColor(0);
      leg->SetFillStyle(0);
      leg->SetTextSize(0.04);
      leg->SetBorderSize(0);
      leg->AddEntry(fsig,"Random-cone-corrected MC","f");
      leg->AddEntry(fspike,"Spike","f");   
      leg->Draw("same");
      
      c1->Print(Form("CorrRCSpike_pt%d.eps",(int)fBinsPt[ipt]));
      c1->Print(Form("CorrRCSpike_pt%d.gif",(int)fBinsPt[ipt]));
      c1->Print(Form("CorrRCSpike_pt%d.C",(int)fBinsPt[ipt]));
      
    } 
  }

    
   

}

//.........这里部分代码省略.........
	PtGood_RANK_matchedClosest.Fill( sum2pho.pt(),ww );
	EtaGood_RANK_matchedClosest.Fill( etaMaxSC ,ww);
	NvtGood_RANK_matchedClosest.Fill( nvtx_,ww );
	if (!isData) NpuGood_RANK_matchedClosest.Fill(npu,ww);
      }
      


      //*** BDT
      int TMVAind = -1;
      float TMVAmax = -1000.;

      for ( int uu = 0; uu < nvtx_; uu++){
			
	diphopt   = vAna.diphopt(uu);
	logsumpt2 = vAna.logsumpt2(uu);
	nch       = vAna.nch(uu);
	nchthr    = vAna.nchthr(uu);
	ptbal     = vAna.ptbal(uu);
	ptasym    = vAna.ptasym(uu);
	ptmax     = vAna.ptmax(uu);
	ptmax3    = vAna.ptmax3(uu);
	sumtwd    = vAna.sumtwd(uu);

	htemp->Reset();
	for (unsigned int kk = 0; kk < PVtracks->size(); ++kk){
	  if (PVtracks_PVindex->at(kk) == uu){
	    float tkpt = sqrt(PVtracks->at(kk).perp2()) ;
	    htemp ->Fill( (tkpt < 10. ? tkpt : 10.) );
	  }
	}

	if ( htemp->GetEntries()!=0 ) 
	  dk = htemp->KolmogorovTest(hTrackPt);

	//--- Evaluate TMVA
	Double_t mva = tmvaReader_->EvaluateMVA( "BDTG" ); 
	
	BDToutput.Fill( mva );
	if ( fabs( TrueVertex_Z - PV_z->at(uu) ) < 1.)  
	  BDToutput_sig.Fill( mva );
	else 
	  BDToutput_bkg.Fill( mva );

	// take the highest score
	if ( mva > TMVAmax) {
	  TMVAmax = mva;
	  TMVAind = uu;
	}
	
      } // end loop over vertices


      //-- matching 1cm
      if ( fabs( TrueVertex_Z - PV_z->at(TMVAind) ) < 1.) {
	PtGood_BDT.Fill( sum2pho.pt(),ww );
	EtaGood_BDT.Fill( etaMaxSC ,ww);
	NvtGood_BDT.Fill( nvtx_ ,ww);
	if (!isData) NpuGood_BDT.Fill(npu,ww);
      }
      //-- matching closest vtx
      if ( iClosest == TMVAind){
	PtGood_BDT_matchedClosest.Fill( sum2pho.pt(),ww );
	EtaGood_BDT_matchedClosest.Fill( etaMaxSC ,ww);
	NvtGood_BDT_matchedClosest.Fill( nvtx_,ww );
	if (!isData) NpuGood_BDT_matchedClosest.Fill(npu,ww);

//.........这里部分代码省略.........
	  fitFunc = new TF1 ("fitFunc","gaus(0)",histoNoise->GetXaxis()->GetXmin(),histoNoise->GetXaxis()->GetXmax());
	}
	fitFunc->SetRange(histoNoise->GetXaxis()->GetXmin(),histoNoise->GetXaxis()->GetXmax());
	fitFunc->SetParameters(histoNoiseStrip->Integral(),noiseMean_,noiseRMS_);
	result = histoNoiseStrip->Fit(fitFunc,"QSR");

	if(result.Get()){
	    fitStatus_     = result->Status();
	    fitGausNormalization_  = fitFunc->GetParameter(0);
	    fitGausMean_   = fitFunc->GetParameter(1);
	    fitGausSigma_  = fitFunc->GetParameter(2);
	    fitGausNormalizationError_  = fitFunc->GetParError(0);
	    fitGausMeanError_  = fitFunc->GetParError(1);
	    fitGausSigmaError_ = fitFunc->GetParError(2);
	    fitChi2_           = result->Chi2();
	    fitChi2Probab_     = result->Prob();

	    noiseIntegral3SigmaFromFit_ = (histoNoiseStrip->Integral(histoNoiseStrip->FindBin(noiseMean_+fitGausSigma_*3),histoNoiseStrip->GetNbinsX()+1) + histoNoiseStrip->Integral(0,histoNoiseStrip->FindBin(noiseMean_-fitGausSigma_*3)))/histoNoiseStrip->Integral();
	    noiseIntegral4SigmaFromFit_ = (histoNoiseStrip->Integral(histoNoiseStrip->FindBin(noiseMean_+fitGausSigma_*4),histoNoiseStrip->GetNbinsX()+1) + histoNoiseStrip->Integral(0,histoNoiseStrip->FindBin(noiseMean_-fitGausSigma_*4)))/histoNoiseStrip->Integral();
	    noiseIntegral5SigmaFromFit_ = (histoNoiseStrip->Integral(histoNoiseStrip->FindBin(noiseMean_+fitGausSigma_*5),histoNoiseStrip->GetNbinsX()+1) + histoNoiseStrip->Integral(0,histoNoiseStrip->FindBin(noiseMean_-fitGausSigma_*5)))/histoNoiseStrip->Integral();
	    
	    jBValue_   = (histoNoiseStrip->Integral()/6)*(noiseSkewness_*noiseSkewness_+(noiseKurtosis_*noiseKurtosis_)/4);	  
	    jBProbab_  = ROOT::Math::chisquared_cdf_c(jBValue_,2);

	    if(randomHisto == 0 or randomHisto == NULL)
	      randomHisto = (TH1F*) histoNoiseStrip->Clone("randomHisto");	    	    
	    randomHisto->Reset();
	    randomHisto->SetDirectory(0);     
	
      
	    if(integral != 0){	      
	      if(generateRandomDistribution){
		randomHisto->FillRandom("fitFunc",histoNoiseStrip->Integral());	    
		kSValue_  = histoNoiseStrip->KolmogorovTest(randomHisto,"MN");
		kSProbab_ = histoNoiseStrip->KolmogorovTest(randomHisto,"N");	    
		aDValue_  = histoNoiseStrip->AndersonDarlingTest(randomHisto,"T");
		aDProbab_ = histoNoiseStrip->AndersonDarlingTest(randomHisto);
	      }
	      else{
		
		randomHisto->Add(fitFunc);		
		kSValue_  = histoNoiseStrip->KolmogorovTest(randomHisto,"MN"); 
		kSProbab_ = histoNoiseStrip->KolmogorovTest(randomHisto,"N");
		// AD test
		ROOT::Fit::BinData data1;
		ROOT::Fit::BinData data2;
		ROOT::Fit::FillData(data1,histoNoiseStrip,0);
		data2.Initialize(randomHisto->GetNbinsX()+1,1);
		for(int ibin = 0; ibin < randomHisto->GetNbinsX(); ibin++){ 
		  if(histoNoiseStrip->GetBinContent(ibin+1) != 0 or randomHisto->GetBinContent(ibin+1) >= 1)
		    data2.Add(randomHisto->GetBinCenter(ibin+1),randomHisto->GetBinContent(ibin+1),randomHisto->GetBinError(ibin+1));
		}
	  
		double probab;
		double value;
		ROOT::Math::GoFTest::AndersonDarling2SamplesTest(data1,data2,probab,value);
		aDValue_ = value;
		aDProbab_ = probab;
	      }
	    }
	}
	else
	  noFitResult++;
	
	if(not histoBranches){
	  noiseDistribution_.clear();