Main Page | Modules | Namespace List | Class Hierarchy | Alphabetical List | Class List | Directories | File List | Namespace Members | Class Members | File Members | Related Pages

NueExpGenerator Class Reference

#include <NueExpGenerator.h>

List of all members.

Public Member Functions

 NueExpGenerator ()
void Run (string errfile, string data, string output)
void SetOutputFile (string name)
Position * GenerateExperimentSet (double *num, double *err)
double CalculateDeltaLog (double nexp, double nobs, double b=-1)
void WriteToFile (Position *pos)
void SetSeed (double in)
void SetNumberOfExp (int num)
void SetScale (double in)
void ReScale (double *num)
void SetMinTotal (double num)
void CleanPos (Position *pos)
void SetOffset (int in)
void SetReadNumber (int in)
bool DetermineCuts (Position *pos)
void SetMethod (int met)
void SetFitMethod (int met)
void SetObservation (int obs)
double CalculateChi2withBestFit (double b, double s, double k, double errBg, double errK, int n)

Private Attributes

string outFile
vector< string > files
double fSeed
TRandom3 fRand
int fNumberOfExp
double fErrors [6]
double fNumbers [6]
vector< Position * > posList
double fMinTotal
double fOscPar [6]
int fOffSet
int fReadUntil
double fScale
int fMethod
int fFitMethod
double fObservation

Constructor & Destructor Documentation

NueExpGenerator::NueExpGenerator (   ) 
 

Definition at line 7 of file NueExpGenerator.cxx.

References fFitMethod, fMethod, fMinTotal, fNumberOfExp, fObservation, fOffSet, fReadUntil, fScale, fSeed, and outFile.

00007                                 { 
00008   outFile = "default.root";
00009   fSeed = -10;
00010   fNumberOfExp = 10000;
00011   fScale = 1.0;
00012   fMinTotal = 0.0;
00013 
00014   fOffSet = -1;
00015   fReadUntil = -1;
00016   fMethod = 0;
00017   fFitMethod = 0;
00018   fObservation = 35;
00019 }

Member Function Documentation

double NueExpGenerator::CalculateChi2withBestFit (  double  b,
double  s,
double  k,
double  errBg,
double  errK,
int  n
) 
 

Definition at line 247 of file NueExpGenerator.cxx.

References s().

Referenced by GenerateExperimentSet().

00248 {
00249     //First calculate the best fit to the particular point mu given these values of b,s,k
00250     // These are the same formulae used in the analytic feldman cousins
00251 
00252     double sol_A = 1 + errK/errBg*s*s;
00253     double sol_B = errBg + s*k - 2*errK*s*s*b/errBg - b + s*s*errK; 
00254     double sol_C = -n*errBg + s*k*errBg - s*s*errK*b - s*k*b + s*s*b*b*errK/errBg;
00255 
00256     double rad = sol_B*sol_B-4*sol_A*sol_C;
00257     if(rad < 0) std::cout<<"negative radical - not sure about this...."<<std::endl;
00258 
00259    double betaHat = (-sol_B + TMath::Sqrt(rad))/(2*sol_A);
00260    double kHat = k - s*errK/errBg*(b-betaHat);
00261 
00262    double nexp = betaHat + s*kHat;
00263    double chisq = 2*(nexp - n + n*TMath::Log(n/nexp)) 
00264                              + (betaHat-b)*(betaHat-b)/errBg + (kHat - k)*(kHat-k)/errK;
00265 
00266    double bBest = 0.0;
00267      // Working out  chi^2 BF 
00268      double chiBF = 0.0;
00269     
00270      if(n > b){
00271           chiBF = 0.0;
00272      }else{
00273        double  bBest = 0.5*(b - errBg + TMath::Sqrt( (b-errBg)*(b - errBg) + 4*n*errBg));
00274         // then we have no signal at best fit point
00275 
00276         chiBF = 2*(bBest - n + n*TMath::Log(n/bBest)) + (b-bBest)*(b-bBest)/errBg;
00277     }
00278 
00279    if(chiBF > chisq-1e-10 && chiBF > 0) 
00280        cout<<n<<"  "<<chisq<<"  "<<chiBF<<"  "<<bBest<<"  "<<betaHat<<"  "<<nexp<<"  "<<kHat<<"  "<<k<<"  "<<errBg<<"  "<<errK<<"  "<<s<<endl;
00281            
00282     return (chisq - chiBF);
00283 }

double NueExpGenerator::CalculateDeltaLog (  double  nexp,
double  nobs,
double  b = -1
) 
 

Definition at line 216 of file NueExpGenerator.cxx.

Referenced by GenerateExperimentSet().

00217 {   
00218   //I assume that if the observation has a value greater than the minimum number
00219   // of events, than it can be predicted perfectly
00220   double bound = fMinTotal;
00221   if(b > -0.5) bound = b; 
00222 
00223   double chi2 = 0;
00224   
00225   if(nobs < 1e-2){
00226      chi2 = 2*nexp;  //protection against nobs = 0
00227      if(nobs < bound)
00228      chi2 -= 2*(bound);
00229   }
00230   else{
00231      chi2 = 2*(nexp - nobs + nobs*TMath::Log(nobs/nexp));
00232   
00233      if(nobs < bound)
00234        chi2 -= 2*(bound - nobs + nobs*TMath::Log(nobs/bound));
00235    }
00236 
00237 //   if(chi2 < -1e-10) cout<<"really: "<< 2*(nexp - nobs + nobs*TMath::Log(nobs/nexp))<<"  "<<nexp<<" "<<fMinTotal<<"  "<<chi2<<endl;
00238 
00239   if(nexp+1e-4 < fMinTotal){
00240 //    cout<<"The expected number is less then the minimum number of events - "
00241 //        <<" something is messed up, please check: min: "<<fMinTotal<<" nexp: "<<nexp<<endl;
00242   }
00243   
00244   return chi2;
00245 }

void NueExpGenerator::CleanPos (  Position *  pos  ) 
 

Definition at line 59 of file NueExpGenerator.cxx.

References Position::deltaLog, Position::fDirectory, and Position::nevent.

Referenced by Run().

00060 {
00061    delete pos->nevent;
00062    delete pos->deltaLog;
00063    delete pos->fDirectory;
00064 //   delete pos;
00065 }

bool NueExpGenerator::DetermineCuts (  Position *  pos  ) 
 

Definition at line 67 of file NueExpGenerator.cxx.

References count, Position::deltaLog, Position::ninety, Position::sixtyeight, Position::threesigma, and total().

Referenced by Run().

00068 {
00069   TH1D* hist = pos->deltaLog;
00070 
00071   double total = hist->GetSum();
00072   double count = 0.0;
00073   bool found68 = false;
00074   bool found90 = false;
00075   bool found3s = false;
00076 
00077   for(int i = 0; i < hist->GetNbinsX(); i++){
00078     count += hist->GetBinContent(i);
00079     if(count/total > 0.68 && !found68){ found68 = true; pos->sixtyeight = hist->GetBinCenter(i);}
00080     if(count/total > 0.90 && !found90){ found90 = true; pos->ninety = hist->GetBinCenter(i);}
00081     if(count/total > 0.9973 && !found3s){ found3s = true; pos->threesigma = hist->GetBinCenter(i);}
00082 
00083     if(count/total > 0.9973) break;
00084   }
00085 
00086   return found90;
00087 }

Position * NueExpGenerator::GenerateExperimentSet (  double *  num,
double *  err
) 
 

Definition at line 96 of file NueExpGenerator.cxx.

References CalculateChi2withBestFit(), CalculateDeltaLog(), Position::deltaLog, Position::fDirectory, Position::fEventNo, fFitMethod, fMethod, fNumbers, fOscPar, fRand, Position::id, ID, max, and Position::nevent.

Referenced by Run().

00097 {
00098    float nominal = 0.0;
00099 //   int nevent = 0;
00100    double totBG = 0.0;
00101 
00102    Position* pos = new Position;
00103       
00104    for(int i = 0; i < 5; i++){
00105       pos->fEventNo[i] = num[i];      
00106       nominal += num[i];
00107       if(i != 2) totBG += num[i];
00108    }
00109    fNumbers[5] = nominal;
00110 
00111    char ID[100];
00112    int one  = (int) rint(fOscPar[0]*1e3);
00113    int two  =  (int) rint(1e2*fOscPar[1]);
00114    int three = (int) rint(1e4*fOscPar[2]);
00115    if(fOscPar[3] > 0) 
00116     sprintf(ID, "Pos_%03d_%03d_%02d_N", one, two, three);
00117    else
00118     sprintf(ID, "Pos_%03d_%03d_%02d_I", one, two, three);
00119 
00120 //   cout<<int(fOscPar[0]*1e3)<<"  "<<rint(1e2*fOscPar[1])<<"  "<<rint(1e4*fOscPar[2])<<"  "
00121 //       <<fOscPar[0]*1e3<<"  "<<1e2*fOscPar[1]<<"  "<<1e4*fOscPar[2]<<endl;
00122 //   cout<<"Name: "<<ID<<endl;
00123    pos->fDirectory = new TDirectory(ID, ID);
00124    pos->fDirectory->cd();
00125 
00126    TString n1 = "nevent";
00127    TString n2 = "DeltaLogLikely";
00128 
00129    pos->id = string(ID);
00130    int max = 199;
00131    if(nominal > 60) max = int(3*nominal);
00132    pos->nevent = new TH1D(n1, n1, max+1, -1, max);
00133    pos->deltaLog = new TH1D(n2, n2, 50000, 0, 1000);
00134 
00135    TString n3[6];
00136    n3[0] = n1+"nc";  n3[1] = n1+"cc";  n3[2] = n1+"nue";
00137    n3[3] = n1+"tau"; n3[4] = n1+"bnue"; n3[5] = n1 + "tot";
00138 
00139 //   for(int i = 0; i < 6; i++){
00140 //     pos->sysevents[i] = new TH1D(n3[i], n3[i], max+1, -1, max);
00141 //     statevents[i] = new TH1D(n1, n1, max+1, -1, max);
00142 //   }
00143 
00144 
00145    gDirectory->cd("/");
00146 
00147 
00148    //fMethod 0, 1, 2
00149    // 0  shift the observation but not expectation
00150    // 1  shift the observation and expectation together
00151    // 2  Poissson(obs) and Gauss(exp) 
00152 
00153 //   cout<<err[5]<<"  "<<totBG<<"  "<<err[2]<<"   "<<num[2]<<endl; 
00154    for(int k = 0; k < fNumberOfExp; k++)
00155    {
00156      double shift = fRand.Gaus(0, err[5]);
00157      if(shift < -1){       cout<<"Wow "<<shift<<" "<<err[5]<<endl; shift = -1; }
00158 
00159      double shiftedBg = totBG*(1+shift); 
00160      
00161      double bgSeed = shiftedBg;
00162      if(fMethod >= 2) bgSeed = totBG;
00163 
00164      int bgObs =  fRand.Poisson(bgSeed);
00165 
00166      shift = fRand.Gaus(0, err[ClassType::nue]);
00167      if(shift < -1) {     cout<<"Wow "<<shift<<" "<<err[ClassType::nue]<<endl; shift = -1; }
00168  
00169      double shiftedSig = num[ClassType::nue]*(1+shift);
00170      
00171      double sigSeed = shiftedSig;
00172      if(fMethod >= 2) sigSeed = num[ClassType::nue];
00173   
00174      int sigObs = fRand.Poisson(sigSeed);
00175 
00176      int nObs = sigObs + bgObs;
00177      
00178      double nexp = nominal;
00179      double bexp, sexp;
00180      bexp = sexp = 0.0;
00181      if(fMethod == 1 || fMethod >= 2){
00182           nexp = shiftedBg + shiftedSig;
00183           bexp = shiftedBg;
00184           sexp = shiftedSig;
00185      }
00186 
00187 //     cout<<fMethod<<"  "<<nexp<<"  "<<nominal<<endl;
00188 
00189      pos->nevent->Fill(nObs);
00190      if(k%5000000 == 0 && k > 0) cout<<k<<"/"<<fNumberOfExp<<endl;
00191 
00192      double dll = 0;
00193  
00194      // Next we implement the special condition for n == nObs
00195      if(nObs == fObservation){
00196          bexp = totBG;
00197          sexp = num[2];
00198          nexp = bexp + sexp;
00199      }   
00200 
00201      if(fFitMethod == 0) dll = CalculateDeltaLog(nexp, nObs, bexp);
00202      else{     
00203         double errBg = err[5]*err[5]*bexp*bexp;;
00204         double errK = err[2]*err[2]*sexp*sexp/(num[2]*num[2]);
00205        dll = CalculateChi2withBestFit(bexp, num[2], sexp/num[2], errBg,  errK, nObs);
00206      }
00207 //     CalculateDeltaLog(nexp, nObs, shiftedBg);
00208 
00209      pos->deltaLog->Fill(dll);
00210    }
00211 
00212 //   cout<<"I'm going home"<<endl;
00213    return pos;
00214 }

void NueExpGenerator::ReScale (  double *  num  ) 
 

Definition at line 89 of file NueExpGenerator.cxx.

References fScale.

Referenced by Run().

00090 {
00091   for(int i = 0; i < 5; i++){
00092      num[i] = fScale*num[i];
00093   }
00094 }

void NueExpGenerator::Run (  string  errfile,
string  data,
string  output
) 
 

Definition at line 21 of file NueExpGenerator.cxx.

References NueGenConfig::AddInputFile(), NueGenConfig::CheckConfig(), CleanPos(), count, DetermineCuts(), fErrors, fMinTotal, fNumberOfExp, fNumbers, fOffSet, fOscPar, fRand, fReadUntil, fScale, fSeed, GenerateExperimentSet(), NueGenConfig::GetErrors(), NueGenConfig::GetOscPar(), NueGenConfig::LoadNextNumberSet(), ReScale(), NueGenConfig::SetOffset(), SetOutputFile(), NueGenConfig::SetReadNumber(), and WriteToFile().

00022 {
00023   SetOutputFile(output);
00024   NueGenConfig* gc = new NueGenConfig(errinput);
00025   gc->AddInputFile(data);
00026   if(!gc->CheckConfig()) return;
00027 
00028   if(fOffSet >= 0) gc->SetOffset(fOffSet);
00029   if(fReadUntil > 0) gc->SetReadNumber(fReadUntil);
00030 
00031   double* temp = gc->GetErrors();
00032   for(int i = 0; i < 6; i++) fErrors[i] = temp[i];
00033 
00034   if(fSeed > -1) fRand.SetSeed((int) fSeed);
00035 
00036   int count = 0;
00037 
00038   while(gc->LoadNextNumberSet(fNumbers)){
00039      if(fScale != 1){
00040          ReScale(fNumbers);
00041          if(count == 0)   fMinTotal *= fScale;
00042      }
00043      double *temp = gc->GetOscPar();
00044      for(int i = 0; i < 4; i++) fOscPar[i] = temp[i];
00045      Position* pos = GenerateExperimentSet(fNumbers, fErrors);
00046      DetermineCuts(pos);
00047      WriteToFile(pos); 
00048      CleanPos(pos);
00049      delete pos;
00050      int temp2 = 1;
00051      if(fNumberOfExp < 1000*10000) temp2 = int(1000*10000.0/fNumberOfExp);
00052      if(count%(temp2) == 0) cout<<count<<" positions generated."<<endl;
00053      count++;
00054   }
00055 
00056   WriteToFile(0);
00057 }

void NueExpGenerator::SetFitMethod (  int  met  )  [inline]
 

Definition at line 64 of file NueExpGenerator.h.

00064 {fFitMethod = met;};

void NueExpGenerator::SetMethod (  int  met  )  [inline]
 

Definition at line 63 of file NueExpGenerator.h.

00063 {fMethod = met;};

void NueExpGenerator::SetMinTotal (  double  num  )  [inline]
 

Definition at line 55 of file NueExpGenerator.h.

00055 {fMinTotal = num;};

void NueExpGenerator::SetNumberOfExp (  int  num  )  [inline]
 

Definition at line 52 of file NueExpGenerator.h.

00052 {fNumberOfExp = num;};

void NueExpGenerator::SetObservation (  int  obs  )  [inline]
 

Definition at line 65 of file NueExpGenerator.h.

00065 { fObservation = obs; };

void NueExpGenerator::SetOffset (  int  in  )  [inline]
 

Definition at line 59 of file NueExpGenerator.h.

00059 { fOffSet = in; };

void NueExpGenerator::SetOutputFile (  string  name  )  [inline]
 

Definition at line 46 of file NueExpGenerator.h.

Referenced by Run().

00046 { outFile = name;};

void NueExpGenerator::SetReadNumber (  int  in  )  [inline]
 

Definition at line 60 of file NueExpGenerator.h.

00060 {fReadUntil = in; };

void NueExpGenerator::SetScale (  double  in  )  [inline]
 

Definition at line 53 of file NueExpGenerator.h.

00053 {fScale = in;};

void NueExpGenerator::SetSeed (  double  in  )  [inline]
 

Definition at line 51 of file NueExpGenerator.h.

00051 {fSeed = in;};

void NueExpGenerator::WriteToFile (  Position *  pos  ) 
 

Definition at line 286 of file NueExpGenerator.cxx.

References count, Position::deltaLog, fErrors, fMinTotal, fNumberOfExp, fNumbers, fOscPar, Position::id, Position::nevent, Position::ninety, outFile, Position::sixtyeight, and Position::threesigma.

Referenced by Run().

00287 {
00288   static TFile* file = 0;
00289   static TTree* tree = 0;
00290   static char name[256];
00291   static int count = 0;
00292   static double sixty;
00293   static double ninety;
00294   static double threesig;
00295                                                                                 
00296   if(file == 0){
00297      file = new TFile(outFile.c_str(),"RECREATE");
00298      tree = new TTree("PositionTree","PositionTree");
00299      tree->Branch("NC",&fNumbers[0],"NC/D");
00300      tree->Branch("NuMuCC",&fNumbers[1],"NuMuCC/D");
00301      tree->Branch("NueCC",&fNumbers[2], "NueCC/D");
00302      tree->Branch("NuTauCC",&fNumbers[3],"NuTauCC/D");
00303      tree->Branch("BNueCC",&fNumbers[4],"BNueCC/D");
00304      tree->Branch("Total", &fNumbers[5], "Total/D");
00305 
00306 //     tree->Branch("NCErr",&fErrors[0],"NCErr/D");
00307 //     tree->Branch("NuMuCCErr",&fErrors[1],"NuMuCCErr/D");
00308 //     tree->Branch("NueCCErr",&fErrors[2], "NueCCErr/D");
00309 //     tree->Branch("NuTauCCErr",&fErrors[3],"NuTauCCErr/D");
00310 //     tree->Branch("BNueCCErr",&fErrors[4],"BNueCCErr/D");
00311 
00312      tree->Branch("Th13",&fOscPar[0],"Th13/D");
00313      tree->Branch("DeltaCP",&fOscPar[1],"DeltaCP/D");
00314      tree->Branch("DeltaM32",&fOscPar[2], "DeltaM32/D");
00315      tree->Branch("Hierarchy",&fOscPar[3],"Hierarchy/D");
00316      tree->Branch("ParamID", name, "ParamID/C");
00317 
00318      tree->Branch("SixtyEightCut", &sixty, "SixtyEightCut/D");
00319      tree->Branch("NinetyCut", &ninety, "NinetyCut/D");
00320      tree->Branch("ThreeSigma",&threesig, "ThreeSigma/D");
00321 //     tree->SetAutoDelete(kFALSE);
00322 
00323   }
00324   file->cd();
00325  
00326   if(pos != 0){   
00327      string fnh_name = pos->id;
00328      sprintf(name, "%s", fnh_name.c_str());
00329      TDirectory *filedir = file->mkdir(fnh_name.c_str());
00330      filedir->cd();
00331      pos->nevent->Write();
00332      pos->deltaLog->Write();     
00333 //     for(int i = 0; i < 6; i++) pos->sysevents[i]->Write();
00334      file->cd();
00335      sixty = pos->sixtyeight;
00336      ninety = pos->ninety;
00337      threesig = pos->threesigma;
00338 //     cout<<sixty<<"  "<<ninety<<endl;
00339      tree->Fill();
00340      count++;
00341   }
00342   else{ 
00343     cout<<"Writing the trees"<<endl;
00344     tree->Write();
00345 
00346     TTree* tree2 = 0;
00347     tree2 = new TTree("GenerationTree","GenerationTree");
00348     tree2->Branch("MinTotal",&fMinTotal,"MinTotal/D");
00349     tree2->Branch("NumExperiment",&fNumberOfExp,"NumExperiment/I");
00350     tree2->Branch("NCErr",&fErrors[0],"NCErr/D");
00351     tree2->Branch("NuMuCCErr",&fErrors[1],"NuMuCCErr/D");
00352     tree2->Branch("NueCCErr",&fErrors[2], "NueCCErr/D");
00353     tree2->Branch("NuTauCCErr",&fErrors[3],"NuTauCCErr/D");
00354     tree2->Branch("BNueCCErr",&fErrors[4],"BNueCCErr/D");
00355     tree2->Branch("TotalBGErr",&fErrors[5],"TotalBGErr/D");
00356     tree2->Fill();
00357 
00358     tree2->Write(); 
00359     file->Close();
00360   }
00361 }

Member Data Documentation

double NueExpGenerator::fErrors[6] [private]
 

Definition at line 76 of file NueExpGenerator.h.

Referenced by Run(), and WriteToFile().

int NueExpGenerator::fFitMethod [private]
 

Definition at line 88 of file NueExpGenerator.h.

Referenced by GenerateExperimentSet(), and NueExpGenerator().

vector<string> NueExpGenerator::files [private]
 

Definition at line 71 of file NueExpGenerator.h.

int NueExpGenerator::fMethod [private]
 

Definition at line 87 of file NueExpGenerator.h.

Referenced by GenerateExperimentSet(), and NueExpGenerator().

double NueExpGenerator::fMinTotal [private]
 

Definition at line 79 of file NueExpGenerator.h.

Referenced by NueExpGenerator(), Run(), and WriteToFile().

int NueExpGenerator::fNumberOfExp [private]
 

Definition at line 75 of file NueExpGenerator.h.

Referenced by NueExpGenerator(), Run(), and WriteToFile().

double NueExpGenerator::fNumbers[6] [private]
 

Definition at line 77 of file NueExpGenerator.h.

Referenced by GenerateExperimentSet(), Run(), and WriteToFile().

double NueExpGenerator::fObservation [private]
 

Definition at line 89 of file NueExpGenerator.h.

Referenced by NueExpGenerator().

int NueExpGenerator::fOffSet [private]
 

Definition at line 83 of file NueExpGenerator.h.

Referenced by NueExpGenerator(), and Run().

double NueExpGenerator::fOscPar[6] [private]
 

Definition at line 81 of file NueExpGenerator.h.

Referenced by GenerateExperimentSet(), Run(), and WriteToFile().

TRandom3 NueExpGenerator::fRand [private]
 

Definition at line 74 of file NueExpGenerator.h.

Referenced by GenerateExperimentSet(), and Run().

int NueExpGenerator::fReadUntil [private]
 

Definition at line 84 of file NueExpGenerator.h.

Referenced by NueExpGenerator(), and Run().

double NueExpGenerator::fScale [private]
 

Definition at line 86 of file NueExpGenerator.h.

Referenced by NueExpGenerator(), ReScale(), and Run().

double NueExpGenerator::fSeed [private]
 

Definition at line 72 of file NueExpGenerator.h.

Referenced by NueExpGenerator(), and Run().

string NueExpGenerator::outFile [private]
 

Definition at line 70 of file NueExpGenerator.h.

Referenced by NueExpGenerator(), and WriteToFile().

vector<Position*> NueExpGenerator::posList [private]
 

Definition at line 78 of file NueExpGenerator.h.

The documentation for this class was generated from the following files:
Generated on Mon Feb 15 11:09:47 2010 for loon by  doxygen 1.3.9.1