NcPsiDistrib.cxx

Go to the documentation of this file.

 
#include "NcPsiDistrib.h"
#include "Riostream.h"
 
ClassImp(NcPsiDistrib); // Class implementation to enable ROOT I/O

NcPsiDistrib::NcPsiDistrib()
{
 
 fNtrials=10000;
 fNevents=100;
 fNoutcomes=100;
 fProbs=new Double_t[fNoutcomes];
 fSignal=new Int_t[fNoutcomes];
 for(Int_t i=0; i<fNoutcomes; i++)
 {
  fProbs[i]=1./fNoutcomes;
  fSignal[i]=0;
 }
 fRefPsi=0;
 fNbelow=0;
 fNabove=0;
 fPrintFreq=0;
 
 fRangeSet=kFALSE;
 fPsiHisto=new TH1D("psi","psi distribution",100,0,1);
 fPsiHisto->SetDirectory(0);
 fPsiHisto->SetXTitle("psi");
 fPsiHisto->SetFillColor(4);
 
 fSample=new NcSample();
// fSample->SetStoreMode(1);
}

NcPsiDistrib::~NcPsiDistrib()
{
 
 if(fProbs)
 {
  delete[] fProbs;
  fProbs=0;
 }
 if(fSignal)
 {
  delete[] fSignal;
  fSignal=0;
 }
 if(fPsiHisto)
 {
  delete fPsiHisto;
  fPsiHisto=0;
 }
 if(fSample)
 {
  delete fSample;
  fSample=0;
 }
}

void NcPsiDistrib::SetNtrials(Long_t n)
{
 
 fNtrials=n;
}

void NcPsiDistrib::SetNevents(Int_t n)
{
 
 fNevents=n;
}

void NcPsiDistrib::SetNoutcomes(Int_t n, Double_t* p)
{
 
 fNoutcomes=n;
 if(fProbs)
 {
  delete[] fProbs;
 }
 fProbs=new Double_t[fNoutcomes];
 
 SetProbabilities(p);
 
 if(fSignal)
 {
  delete[] fSignal;
  fSignal=new Int_t[fNoutcomes];
  for(Int_t i=0; i<fNoutcomes; i++){ fSignal[i]=0; }
 }
}

void NcPsiDistrib::SetProbabilities(Double_t* p)
{
 
 if(p)
 {
  Double_t sum=0;
  for(Int_t i=0; i<fNoutcomes; i++) { sum+=p[i]; }
  if(sum<=0)
  {
   cout << "NcPsiDistrib: Sum of probabilities is not positive" << endl;
   sum=1;
  }
  for(Int_t i=0; i<fNoutcomes; i++) { fProbs[i]=p[i]/sum; }
 }
 else
 {
  for(Int_t i=0; i<fNoutcomes; i++){ fProbs[i]=1./fNoutcomes; }
 }
}

void NcPsiDistrib::SetSignal(Int_t* s)
{
 
 if(s)
 {
  for(Int_t i=0; i<fNoutcomes; i++)
  {
   fSignal[i]=s[i];
  }
 }
 else
 {
  for(Int_t i=0; i<fNoutcomes; i++)
  {
   fSignal[i]=0;
  }
 }
}

void NcPsiDistrib::SetPsiRange(Int_t nb, Float_t low, Float_t high)
{
 
 fPsiHisto->SetBins(nb,low,high);
 fRangeSet=kTRUE;
}

void NcPsiDistrib::SetReferencePsi(Double_t ref)
{
 
 fRefPsi=ref;
}

void NcPsiDistrib::SetPrintFreq(Int_t freq)
{
 
 fPrintFreq=freq;
}

void NcPsiDistrib::Distribute(Int_t storemode)
{
 
 // Variables we need
 Int_t* data=new Int_t[fNoutcomes];
 NcMath math;
 Double_t psi=0;
 char title[100];
 
 // Reference histo with probabilities
 TH1D ref("ref","ref",fNoutcomes,0,1);
 for(Int_t i=0; i<fNoutcomes; i++)
 {
  ref.SetBinContent(i+1,fProbs[i]);
 }
 
 // Make histogram for events
 TH1D histo("histo","events",fNoutcomes,0,1);
 
 // Reset stuff
 fSample->Reset();
 fSample->SetStoreMode(storemode);
 fPsiHisto->Reset();
 sprintf(title,"psi_%d_%d",fNevents,fNoutcomes);
 fPsiHisto->SetName(title);
 fNbelow=0;
 fNabove=0;
 
 // Set range for psi histogram
 if(!fRangeSet){
  fPsiHisto->SetBins(100,0,1.1*FindMaxPsi());
 }
 
 // Loop over trials
 for(Int_t itrial=0; itrial<fNtrials; itrial++){
  if(fPrintFreq>0 && itrial%fPrintFreq==0) cout << "Trial " << itrial << endl;
 
  // Fill event histogram
  histo.Reset();
  histo.FillRandom(&ref, fNevents);
 
  // Add signal if prescribed
  for(Int_t i=0; i<fNoutcomes; i++)
  {
   histo.SetBinContent(i+1, histo.GetBinContent(i+1)+fSignal[i]);
  }
 
  // Calculate psi
  for(Int_t i=0; i<fNoutcomes; i++)
  {
   data[i]=int(histo.GetBinContent(i+1));
  }
  psi=math.PsiValue(fNoutcomes, data, fProbs);
  fSample->Enter(psi);
  fPsiHisto->Fill(psi);
 
  // Keep track of # below/above ref psi
  if(fRefPsi)
  {
   if(psi<fRefPsi) fNbelow++;
   if(psi>fRefPsi) fNabove++;
  }
 }
 
 delete[] data;
 
 // Set psi distribution histo name
 sprintf(title,"psi distribution (%d events, %d time bins): mean %.2f, sigma %.2e",fNevents,fNoutcomes,fSample->GetMean(1),fSample->GetSigma(1));
 fPsiHisto->SetTitle(title);
}

TH1D* NcPsiDistrib::GetPsiHisto()
{
 
 return fPsiHisto;
}

NcSample* NcPsiDistrib::GetPsiSample()
{
 
 return fSample;
}

Double_t NcPsiDistrib::GetFracBelow()
{
 
 return ((Double_t)fNbelow)/((Double_t)fNtrials);
}

Double_t NcPsiDistrib::GetFracAbove()
{
 
 return ((Double_t)fNabove)/((Double_t)fNtrials);
}

Float_t NcPsiDistrib::FindMaxPsi()
{
 
 // Make data array and initialise with signal
 Int_t* data=new Int_t[fNoutcomes];
 for(Int_t i=0; i<fNoutcomes; i++)
 {
  data[i]=fSignal[i];
 }
 
 // Try putting background events in each bin consecutively
 Float_t maxpsi=0, psi=0;
 NcMath m;
 for(Int_t i=0; i<fNoutcomes; i++)
 {
  data[i]+=fNevents;
  psi=m.PsiValue(fNoutcomes, data, fProbs);
  if(psi>maxpsi) maxpsi=psi;
  data[i]-=fNevents;
 }
 
 // Delete data array
 delete[] data;
 
 return maxpsi;
}