doc/v626/createData_8C.html

#include "TROOT.h"

#include "TMath.h"

#include "TTree.h"

#include "TArrayD.h"

#include "TStyle.h"

#include "TFile.h"

#include "TRandom.h"

#include "Riostream.h"

#include "TCanvas.h"

#include "TMatrixD.h"

#include "TH2F.h"

#include "TLegend.h"

#include "TBranch.h"

#include <vector>


void plot( TString fname = "data.root", TString var0="var0", TString var1="var1" )

{

   TFile* dataFile = TFile::Open( fname );


   if (!dataFile) {

      cout << "ERROR: cannot open file: " << fname << endl;

      return;

   }


   TTree *treeS = (TTree*)dataFile->Get("TreeS");

   TTree *treeB = (TTree*)dataFile->Get("TreeB");


   TCanvas* c = new TCanvas( "c", "", 0, 0, 550, 550 );


   TStyle *TMVAStyle = gROOT->GetStyle("Plain"); // our style is based on Plain

   TMVAStyle->SetOptStat(0);

   TMVAStyle->SetPadTopMargin(0.02);

   TMVAStyle->SetPadBottomMargin(0.16);

   TMVAStyle->SetPadRightMargin(0.03);

   TMVAStyle->SetPadLeftMargin(0.15);

   TMVAStyle->SetPadGridX(0);

   TMVAStyle->SetPadGridY(0);


   TMVAStyle->SetOptTitle(0);

   TMVAStyle->SetTitleW(.4);

   TMVAStyle->SetTitleH(.10);

   TMVAStyle->SetTitleX(.5);

   TMVAStyle->SetTitleY(.9);

   TMVAStyle->SetMarkerStyle(20);

   TMVAStyle->SetMarkerSize(1.6);

   TMVAStyle->cd();


   Float_t xmin = TMath::Min( treeS->GetMinimum( var0 ), treeB->GetMinimum( var0 ) );

   Float_t xmax = TMath::Max( treeS->GetMaximum( var0 ), treeB->GetMaximum( var0 ) );

   Float_t ymin = TMath::Min( treeS->GetMinimum( var1 ), treeB->GetMinimum( var1 ) );

   Float_t ymax = TMath::Max( treeS->GetMaximum( var1 ), treeB->GetMaximum( var1 ) );


   Int_t nbin = 500;

   TH2F* frameS = new TH2F( "DataS", "DataS", nbin, xmin, xmax, nbin, ymin, ymax );

   TH2F* frameB = new TH2F( "DataB", "DataB", nbin, xmin, xmax, nbin, ymin, ymax );


   // project trees

   treeS->Draw( Form("%s:%s>>DataS",var1.Data(),var0.Data()), "", "0" );

   treeB->Draw( Form("%s:%s>>DataB",var1.Data(),var0.Data()

), "", "0" );


   // set style

   frameS->SetMarkerSize( 0.1 );

   frameS->SetMarkerColor( 4 );


   frameB->SetMarkerSize( 0.1 );

   frameB->SetMarkerColor( 2 );


   // legend

   frameS->SetTitle( var1+" versus "+var0+" for signal and background" );

   frameS->GetXaxis()->SetTitle( var0 );

   frameS->GetYaxis()->SetTitle( var1 );


   frameS->SetLabelSize( 0.04, "X" );

   frameS->SetLabelSize( 0.04, "Y" );

   frameS->SetTitleSize( 0.05, "X" );

   frameS->SetTitleSize( 0.05, "Y" );


   // and plot

   frameS->Draw();

   frameB->Draw( "same" );


   // Draw legend

   TLegend *legend = new TLegend( 1 - c->GetRightMargin() - 0.32, 1 - c->GetTopMargin() - 0.12,

                                  1 - c->GetRightMargin(), 1 - c->GetTopMargin() );

   legend->SetFillStyle( 1 );

   legend->AddEntry(frameS,"Signal","p");

   legend->AddEntry(frameB,"Background","p");

   legend->Draw("same");

   legend->SetBorderSize(1);

   legend->SetMargin( 0.3 );


}


TMatrixD* produceSqrtMat( const TMatrixD& covMat )

{

   Double_t sum = 0;

   Int_t size = covMat.GetNrows();;

   TMatrixD* sqrtMat = new TMatrixD( size, size );


   for (Int_t i=0; i< size; i++) {


      sum = 0;

      for (Int_t j=0;j< i; j++) sum += (*sqrtMat)(i,j) * (*sqrtMat)(i,j);


      (*sqrtMat)(i,i) = TMath::Sqrt(TMath::Abs(covMat(i,i) - sum));


      for (Int_t k=i+1 ;k<size; k++) {


         sum = 0;

         for (Int_t l=0; l<i; l++) sum += (*sqrtMat)(k,l) * (*sqrtMat)(i,l);


         (*sqrtMat)(k,i) = (covMat(k,i) - sum) / (*sqrtMat)(i,i);


      }

   }

   return sqrtMat;

}


void getGaussRnd( TArrayD& v, const TMatrixD& sqrtMat, TRandom& R )

{

   // generate "size" correlated Gaussian random numbers


   // sanity check

   const Int_t size = sqrtMat.GetNrows();

   if (size != v.GetSize())

      cout << "<getGaussRnd> too short input vector: " << size << " " << v.GetSize() << endl;


   Double_t* tmpVec = new Double_t[size];


   for (Int_t i=0; i<size; i++) {

      Double_t x, y, z;

      y = R.Rndm();

      z = R.Rndm();

      x = 2*TMath::Pi()*z;

      tmpVec[i] = TMath::Sin(x) * TMath::Sqrt(-2.0*TMath::Log(y));

   }


   for (Int_t i=0; i<size; i++) {

      v[i] = 0;

      for (Int_t j=0; j<=i; j++) v[i] += sqrtMat(i,j) * tmpVec[j];

   }


   delete[] tmpVec;

}


// create the data

void create_lin_Nvar_withFriend(Int_t N = 2000)

{

   const Int_t nvar  = 4;

   const Int_t nvar2 = 1;

   Float_t xvar[nvar];


   // output flie

   TFile* dataFile = TFile::Open( "data.root", "RECREATE" );


   // create signal and background trees

   TTree* treeS = new TTree( "TreeS", "TreeS", 1 );

   TTree* treeB = new TTree( "TreeB", "TreeB", 1 );

   for (Int_t ivar=0; ivar<nvar-nvar2; ivar++) {

     cout << "Creating branch var" << ivar+1 << " in signal tree" << endl;

      treeS->Branch( TString(Form( "var%i", ivar+1 )).Data(), &xvar[ivar], TString(Form( "var%i/F", ivar+1 )).Data() );

      treeB->Branch( TString(Form( "var%i", ivar+1 )).Data(), &xvar[ivar], TString(Form( "var%i/F", ivar+1 )).Data() );

   }

   TTree* treeSF = new TTree( "TreeSF", "TreeS", 1 );

   TTree* treeBF = new TTree( "TreeBF", "TreeB", 1 );

   for (Int_t ivar=nvar-nvar2; ivar<nvar; ivar++) {

      treeSF->Branch( TString(Form( "var%i", ivar+1 )).Data(), &xvar[ivar], TString(Form( "var%i/F", ivar+1 )).Data() );

      treeBF->Branch( TString(Form( "var%i", ivar+1 )).Data(), &xvar[ivar], TString(Form( "var%i/F", ivar+1 )).Data() );

   }


   TRandom R( 100 );

   Float_t xS[nvar] = {  0.2,  0.3,  0.5,  0.9 };

   Float_t xB[nvar] = { -0.2, -0.3, -0.5, -0.6 };

   Float_t dx[nvar] = {  1.0,  1.0, 1.0, 1.0 };

   TArrayD* v = new TArrayD( nvar );

   Float_t rho[20];

   rho[1*2] = 0.4;

   rho[1*3] = 0.6;

   rho[1*4] = 0.9;

   rho[2*3] = 0.7;

   rho[2*4] = 0.8;

   rho[3*4] = 0.93;


   // create covariance matrix

   TMatrixD* covMatS = new TMatrixD( nvar, nvar );

   TMatrixD* covMatB = new TMatrixD( nvar, nvar );

   for (Int_t ivar=0; ivar<nvar; ivar++) {

      (*covMatS)(ivar,ivar) = dx[ivar]*dx[ivar];

      (*covMatB)(ivar,ivar) = dx[ivar]*dx[ivar];

      for (Int_t jvar=ivar+1; jvar<nvar; jvar++) {

         (*covMatS)(ivar,jvar) = rho[(ivar+1)*(jvar+1)]*dx[ivar]*dx[jvar];

         (*covMatS)(jvar,ivar) = (*covMatS)(ivar,jvar);


         (*covMatB)(ivar,jvar) = rho[(ivar+1)*(jvar+1)]*dx[ivar]*dx[jvar];

         (*covMatB)(jvar,ivar) = (*covMatB)(ivar,jvar);

      }

   }


   cout << "signal covariance matrix: " << endl;

   covMatS->Print();

   cout << "background covariance matrix: " << endl;

   covMatB->Print();


   // produce the square-root matrix

   TMatrixD* sqrtMatS = produceSqrtMat( *covMatS );

   TMatrixD* sqrtMatB = produceSqrtMat( *covMatB );


   // loop over species

   for (Int_t itype=0; itype<2; itype++) {


      Float_t*  x;

      TMatrixD* m;

      if (itype == 0) { x = xS; m = sqrtMatS; cout << "- produce signal" << endl; }

      else            { x = xB; m = sqrtMatB; cout << "- produce background" << endl; }


      // event loop

      TTree* tree  = (itype==0) ? treeS : treeB;

      TTree* treeF = (itype==0) ? treeSF : treeBF;

      for (Int_t i=0; i<N; i++) {


         if (i%1000 == 0) cout << "... event: " << i << " (" << N << ")" << endl;

         getGaussRnd( *v, *m, R );


         for (Int_t ivar=0; ivar<nvar; ivar++) xvar[ivar] = (*v)[ivar] + x[ivar];


         tree->Fill();

         treeF->Fill();

      }

   }


//    treeS->AddFriend(treeSF);

//    treeB->AddFriend(treeBF);


   // write trees

   treeS->Write();

   treeB->Write();

   treeSF->Write();

   treeBF->Write();


   treeS->Show(0);

   treeB->Show(1);


   dataFile->Close();

   cout << "created data file: " << dataFile->GetName() << endl;


}


// create the tree

TTree* makeTree_lin_Nvar( TString treeName, TString treeTitle, Float_t* x, Float_t* dx, const Int_t nvar, Int_t N )

{

   Float_t xvar[nvar];


   // create tree

   TTree* tree = new TTree(treeName, treeTitle, 1);


   for (Int_t ivar=0; ivar<nvar; ivar++) {

      tree->Branch( TString(Form( "var%i", ivar+1 )).Data(), &xvar[ivar], TString(Form( "var%i/F", ivar+1 )).Data() );

   }


   TRandom R( 100 );

   TArrayD* v = new TArrayD( nvar );

   Float_t rho[20];

   rho[1*2] = 0.4;

   rho[1*3] = 0.6;

   rho[1*4] = 0.9;

   rho[2*3] = 0.7;

   rho[2*4] = 0.8;

   rho[3*4] = 0.93;


   // create covariance matrix

   TMatrixD* covMat = new TMatrixD( nvar, nvar );

   for (Int_t ivar=0; ivar<nvar; ivar++) {

      (*covMat)(ivar,ivar) = dx[ivar]*dx[ivar];

      for (Int_t jvar=ivar+1; jvar<nvar; jvar++) {

         (*covMat)(ivar,jvar) = rho[(ivar+1)*(jvar+1)]*dx[ivar]*dx[jvar];

         (*covMat)(jvar,ivar) = (*covMat)(ivar,jvar);

      }

   }

   //cout << "covariance matrix: " << endl;

   //covMat->Print();


   // produce the square-root matrix

   TMatrixD* sqrtMat = produceSqrtMat( *covMat );


   // event loop

   for (Int_t i=0; i<N; i++) {


      if (i%1000 == 0) cout << "... event: " << i << " (" << N << ")" << endl;

      getGaussRnd( *v, *sqrtMat, R );


      for (Int_t ivar=0; ivar<nvar; ivar++) xvar[ivar] = (*v)[ivar] + x[ivar];


      tree->Fill();

   }


   // write trees

//   tree->Write();


   tree->Show(0);


   cout << "created tree: " << tree->GetName() << endl;

   return tree;

}


// create the data

TTree* makeTree_circ(TString treeName, TString treeTitle, Int_t nvar = 2, Int_t N  = 6000, Float_t radius = 1.0, Bool_t distort = false)

{

   Int_t Nn = 0;

   Float_t xvar[nvar]; //variable array size does not work in interactive mode


   // create signal and background trees

   TTree* tree = new TTree( treeName, treeTitle, 1 );

   for (Int_t ivar=0; ivar<nvar; ivar++) {

      tree->Branch( TString(Form( "var%i", ivar+1 )).Data(), &xvar[ivar], TString(Form( "var%i/F", ivar+1 )).Data() );

   }


   TRandom R( 100 );

   //Float_t phimin = -30, phimax = 130;

   Float_t phimin = -70, phimax = 130;

   Float_t phisig = 5;

   Float_t rsig = 0.1;

   Float_t fnmin = -(radius+4.0*rsig);

   Float_t fnmax = +(radius+4.0*rsig);

   Float_t dfn = fnmax-fnmin;


   // event loop

   for (Int_t i=0; i<N; i++) {

      Double_t r1=R.Rndm(),r2=R.Rndm(), r3;

      r3= r1>r2? r1 :r2;

      Float_t phi;

      if (distort) phi = r3*(phimax - phimin) + phimin;

      else  phi = R.Rndm()*(phimax - phimin) + phimin;

      phi += R.Gaus()*phisig;


      Float_t r = radius;

      r += R.Gaus()*rsig;


      xvar[0] = r*cos(TMath::DegToRad()*phi);

      xvar[1] = r*sin(TMath::DegToRad()*phi);


      for( Int_t j = 2; j<nvar; ++j )

    xvar[j] = dfn*R.Rndm()+fnmin;


      tree->Fill();

   }


   for (Int_t i=0; i<Nn; i++) {


      xvar[0] = dfn*R.Rndm()+fnmin;

      xvar[1] = dfn*R.Rndm()+fnmin;


      for( Int_t j = 2; j<nvar; ++j )

    xvar[j] = dfn*R.Rndm()+fnmin;


      tree->Fill();

   }


   tree->Show(0);

   // write trees

   cout << "created tree: " << tree->GetName() << endl;

   return tree;

}


// create the data

void create_lin_Nvar_2(Int_t N = 50000)

{

   const int nvar = 4;


   // output flie

   TFile* dataFile = TFile::Open( "data.root", "RECREATE" );


   Float_t xS[nvar] = {  0.2,  0.3,  0.5,  0.9 };

   Float_t xB[nvar] = { -0.2, -0.3, -0.5, -0.6 };

   Float_t dx[nvar] = {  1.0,  1.0, 1.0, 1.0 };


   // create signal and background trees

   TTree* treeS = makeTree_lin_Nvar( "TreeS", "Signal tree", xS, dx, nvar, N );

   TTree* treeB = makeTree_lin_Nvar( "TreeB", "Background tree", xB, dx, nvar, N );


   treeS->Write();

   treeB->Write();


   treeS->Show(0);

   treeB->Show(0);


   dataFile->Close();

   cout << "created data file: " << dataFile->GetName() << endl;

}


// create the data

void create_lin_Nvar(Int_t N = 50000)

{

   const Int_t nvar = 4;

   Float_t xvar[nvar];


   // output flie

   TFile* dataFile = TFile::Open( "data.root", "RECREATE" );


   // create signal and background trees

   TTree* treeS = new TTree( "TreeS", "TreeS", 1 );

   TTree* treeB = new TTree( "TreeB", "TreeB", 1 );

   for (Int_t ivar=0; ivar<nvar; ivar++) {

      treeS->Branch( TString(Form( "var%i", ivar+1 )).Data(), &xvar[ivar], TString(Form( "var%i/F", ivar+1 )).Data() );

      treeB->Branch( TString(Form( "var%i", ivar+1 )).Data(), &xvar[ivar], TString(Form( "var%i/F", ivar+1 )).Data() );

   }


   TRandom R( 100 );

   Float_t xS[nvar] = {  0.2,  0.3,  0.5,  0.9 };

   Float_t xB[nvar] = { -0.2, -0.3, -0.5, -0.6 };

   Float_t dx[nvar] = {  1.0,  1.0, 1.0, 1.0 };

   TArrayD* v = new TArrayD( nvar );

   Float_t rho[20];

   rho[1*2] = 0.4;

   rho[1*3] = 0.6;

   rho[1*4] = 0.9;

   rho[2*3] = 0.7;

   rho[2*4] = 0.8;

   rho[3*4] = 0.93;


   // create covariance matrix

   TMatrixD* covMatS = new TMatrixD( nvar, nvar );

   TMatrixD* covMatB = new TMatrixD( nvar, nvar );

   for (Int_t ivar=0; ivar<nvar; ivar++) {

      (*covMatS)(ivar,ivar) = dx[ivar]*dx[ivar];

      (*covMatB)(ivar,ivar) = dx[ivar]*dx[ivar];

      for (Int_t jvar=ivar+1; jvar<nvar; jvar++) {

         (*covMatS)(ivar,jvar) = rho[(ivar+1)*(jvar+1)]*dx[ivar]*dx[jvar];

         (*covMatS)(jvar,ivar) = (*covMatS)(ivar,jvar);


         (*covMatB)(ivar,jvar) = rho[(ivar+1)*(jvar+1)]*dx[ivar]*dx[jvar];

         (*covMatB)(jvar,ivar) = (*covMatB)(ivar,jvar);

      }

   }

   cout << "signal covariance matrix: " << endl;

   covMatS->Print();

   cout << "background covariance matrix: " << endl;

   covMatB->Print();


   // produce the square-root matrix

   TMatrixD* sqrtMatS = produceSqrtMat( *covMatS );

   TMatrixD* sqrtMatB = produceSqrtMat( *covMatB );


   // loop over species

   for (Int_t itype=0; itype<2; itype++) {


      Float_t*  x;

      TMatrixD* m;

      if (itype == 0) { x = xS; m = sqrtMatS; cout << "- produce signal" << endl; }

      else            { x = xB; m = sqrtMatB; cout << "- produce background" << endl; }


      // event loop

      TTree* tree = (itype==0) ? treeS : treeB;

      for (Int_t i=0; i<N; i++) {


         if (i%1000 == 0) cout << "... event: " << i << " (" << N << ")" << endl;

         getGaussRnd( *v, *m, R );


         for (Int_t ivar=0; ivar<nvar; ivar++) xvar[ivar] = (*v)[ivar] + x[ivar];


         tree->Fill();

      }

   }


   // write trees

   treeS->Write();

   treeB->Write();


   treeS->Show(0);

   treeB->Show(1);


   dataFile->Close();

   cout << "created data file: " << dataFile->GetName() << endl;

}


// create the category data

// type = 1 (offset) or 2 (variable = -99)

void create_lin_Nvar_categories(Int_t N = 10000, Int_t type = 2)

{

   const Int_t nvar = 4;

   Float_t xvar[nvar];

   Float_t eta;


   // output flie

   TFile* dataFile = TFile::Open( "data.root", "RECREATE" );


   // create signal and background trees

   TTree* treeS = new TTree( "TreeS", "TreeS", 1 );

   TTree* treeB = new TTree( "TreeB", "TreeB", 1 );

   for (Int_t ivar=0; ivar<nvar; ivar++) {

      treeS->Branch( TString(Form( "var%i", ivar+1 )).Data(), &xvar[ivar], TString(Form( "var%i/F", ivar+1 )).Data() );

      treeB->Branch( TString(Form( "var%i", ivar+1 )).Data(), &xvar[ivar], TString(Form( "var%i/F", ivar+1 )).Data() );

   }


   // add category variable

   treeS->Branch( "eta", &eta, "eta/F" );

   treeB->Branch( "eta", &eta, "eta/F" );


   TRandom R( 100 );

   Float_t xS[nvar] = {  0.2,  0.3,  0.5,  0.9 };

   Float_t xB[nvar] = { -0.2, -0.3, -0.5, -0.6 };

   Float_t dx[nvar] = {  1.0,  1.0, 1.0, 1.0 };

   TArrayD* v = new TArrayD( nvar );

   Float_t rho[20];

   rho[1*2] = 0.0;

   rho[1*3] = 0.0;

   rho[1*4] = 0.0;

   rho[2*3] = 0.0;

   rho[2*4] = 0.0;

   rho[3*4] = 0.0;

   if (type != 1) {

      rho[1*2] = 0.6;

      rho[1*3] = 0.7;

      rho[1*4] = 0.9;

      rho[2*3] = 0.8;

      rho[2*4] = 0.9;

      rho[3*4] = 0.93;

   }


   // create covariance matrix

   TMatrixD* covMatS = new TMatrixD( nvar, nvar );

   TMatrixD* covMatB = new TMatrixD( nvar, nvar );

   for (Int_t ivar=0; ivar<nvar; ivar++) {

      (*covMatS)(ivar,ivar) = dx[ivar]*dx[ivar];

      (*covMatB)(ivar,ivar) = dx[ivar]*dx[ivar];

      for (Int_t jvar=ivar+1; jvar<nvar; jvar++) {

         (*covMatS)(ivar,jvar) = rho[(ivar+1)*(jvar+1)]*dx[ivar]*dx[jvar];

         (*covMatS)(jvar,ivar) = (*covMatS)(ivar,jvar);


         (*covMatB)(ivar,jvar) = rho[(ivar+1)*(jvar+1)]*dx[ivar]*dx[jvar];

         (*covMatB)(jvar,ivar) = (*covMatB)(ivar,jvar);

      }

   }

   cout << "signal covariance matrix: " << endl;

   covMatS->Print();

   cout << "background covariance matrix: " << endl;

   covMatB->Print();


   // produce the square-root matrix

   TMatrixD* sqrtMatS = produceSqrtMat( *covMatS );

   TMatrixD* sqrtMatB = produceSqrtMat( *covMatB );


   // loop over species

   for (Int_t itype=0; itype<2; itype++) {


      Float_t*  x;

      TMatrixD* m;

      if (itype == 0) { x = xS; m = sqrtMatS; cout << "- produce signal" << endl; }

      else            { x = xB; m = sqrtMatB; cout << "- produce background" << endl; }


      // event loop

      TTree* tree = (itype==0) ? treeS : treeB;

      for (Int_t i=0; i<N; i++) {


         if (i%1000 == 0) cout << "... event: " << i << " (" << N << ")" << endl;

         getGaussRnd( *v, *m, R );


         eta = 2.5*2*(R.Rndm() - 0.5);

         Float_t offset = 0;

         if (type == 1) offset = TMath::Abs(eta) > 1.3 ? 0.8 : -0.8;

         for (Int_t ivar=0; ivar<nvar; ivar++) xvar[ivar] = (*v)[ivar] + x[ivar] + offset;

         if (type != 1 && TMath::Abs(eta) > 1.3) xvar[nvar-1] = -5;


         tree->Fill();

      }

   }


   // write trees

   treeS->Write();

   treeB->Write();


   treeS->Show(0);

   treeB->Show(1);


   dataFile->Close();

   cout << "created data file: " << dataFile->GetName() << endl;

}


// create the data

void create_lin_Nvar_weighted(Int_t N = 10000, int WeightedSignal=0, int WeightedBkg=1, Float_t BackgroundContamination=0, Int_t seed=100)

{

   const Int_t nvar = 4;

   Float_t xvar[nvar];

   Float_t weight;


   cout << endl << endl << endl;

   cout << "please use .L createData.C++ if you want to run this MC geneation" <<endl;

   cout << "otherwise you will wait for ages!!! " << endl;

   cout << endl << endl << endl;


   // output flie

   TString fileName;

   if (BackgroundContamination) fileName = Form("linCorGauss%d_weighted+background.root",seed);

   else                         fileName = Form("linCorGauss%d_weighted.root",seed);


   TFile* dataFile = TFile::Open( fileName.Data(), "RECREATE" );


   // create signal and background trees

   TTree* treeS = new TTree( "TreeS", "TreeS", 1 );

   TTree* treeB = new TTree( "TreeB", "TreeB", 1 );

   for (Int_t ivar=0; ivar<nvar; ivar++) {

      treeS->Branch( TString(Form( "var%i", ivar+1 )).Data(), &xvar[ivar], TString(Form( "var%i/F", ivar+1 )).Data() );

      treeB->Branch( TString(Form( "var%i", ivar+1 )).Data(), &xvar[ivar], TString(Form( "var%i/F", ivar+1 )).Data() );

   }

   if (WeightedSignal||BackgroundContamination>0||1) treeS->Branch( "weight", &weight,"weight/F" );

   if (WeightedBkg)    treeB->Branch( "weight", &weight,"weight/F" );


   TRandom R( seed );

   Float_t xS[nvar] = {  0.2,  0.3,  0.4,  0.8 };

   Float_t xB[nvar] = { -0.2, -0.3, -0.4, -0.5 };

   Float_t dx[nvar] = {  1.0,  1.0, 1.0, 1.0 };

   TArrayD* v = new TArrayD( nvar );

   Float_t rho[20];

   rho[1*2] = 0.4;

   rho[1*3] = 0.6;

   rho[1*4] = 0.9;

   rho[2*3] = 0.7;

   rho[2*4] = 0.8;

   rho[3*4] = 0.93;


   // create covariance matrix

   TMatrixD* covMatS = new TMatrixD( nvar, nvar );

   TMatrixD* covMatB = new TMatrixD( nvar, nvar );

   for (Int_t ivar=0; ivar<nvar; ivar++) {

      (*covMatS)(ivar,ivar) = dx[ivar]*dx[ivar];

      (*covMatB)(ivar,ivar) = dx[ivar]*dx[ivar];

      for (Int_t jvar=ivar+1; jvar<nvar; jvar++) {

         (*covMatS)(ivar,jvar) = rho[(ivar+1)*(jvar+1)]*dx[ivar]*dx[jvar];

         (*covMatS)(jvar,ivar) = (*covMatS)(ivar,jvar);


         (*covMatB)(ivar,jvar) = rho[(ivar+1)*(jvar+1)]*dx[ivar]*dx[jvar];

         (*covMatB)(jvar,ivar) = (*covMatB)(ivar,jvar);

      }

   }

   cout << "signal covariance matrix: " << endl;

   covMatS->Print();

   cout << "background covariance matrix: " << endl;

   covMatB->Print();


   // produce the square-root matrix

   TMatrixD* sqrtMatS = produceSqrtMat( *covMatS );

   TMatrixD* sqrtMatB = produceSqrtMat( *covMatB );


   // loop over species

   for (Int_t itype=0; itype<2; itype++) {


      Float_t*  x;

      TMatrixD* m;

      if (itype == 0) { x = xS; m = sqrtMatS; cout << "- produce signal" << endl; }

      else            { x = xB; m = sqrtMatB; cout << "- produce background" << endl; }


      // event loop

      TTree* tree = (itype==0) ? treeS : treeB;

      Int_t i=0;

      do {

         getGaussRnd( *v, *m, R );


         for (Int_t ivar=0; ivar<nvar; ivar++) xvar[ivar] = (*v)[ivar] + x[ivar];

         //         for (Int_t ivar=0; ivar<nvar; ivar++) xvar[ivar] = R.Uniform()*10.-5.;


         //         weight = 0.5 / (TMath::Gaus( (xvar[nvar-1]-x[nvar-1]), 0, 1.1) );

         // weight = TMath::Gaus(0.675,0,1) / (TMath::Gaus( (xvar[nvar-1]-x[nvar-1]), 0, 1.) );

         weight = 0.8 / (TMath::Gaus( ((*v)[nvar-1]), 0, 1.09) );

         Double_t tmp=R.Uniform()/0.00034;

         if (itype==0 && !WeightedSignal) {

            weight = 1;

            tree->Fill();

            i++;

         } else if (itype==1 && !WeightedBkg) {

            weight = 1;

            tree->Fill();

            i++;

         }

         else {

            if (tmp < weight){

               weight = 1./weight;

               tree->Fill();

               if (i%10 == 0) cout << "... event: " << i << " (" << N << ")" << endl;

               i++;

            }

         }

      } while (i<N);

   }


   if (BackgroundContamination > 0){  // add "background contamination" in the Signal (which later is again "subtracted" with

            // using (statistically indepentent) background events with negative weight)

      Float_t*  x=xB;

      TMatrixD* m = sqrtMatB;

      TTree* tree = treeS;

      for (Int_t i=0; i<N*BackgroundContamination*2; i++) {

         if (i%1000 == 0) cout << "... event: " << i << " (" << N << ")" << endl;

         getGaussRnd( *v, *m, R );

         for (Int_t ivar=0; ivar<nvar; ivar++) xvar[ivar] = (*v)[ivar] + x[ivar];


         // add weights

         if (i%2) weight = 1;

         else weight = -1;


         tree->Fill();

      }

   }


   // write trees

   treeS->Write();

   treeB->Write();


   treeS->Show(0);

   treeB->Show(1);


   TH1F *h[4];

   TH1F *hw[4];

   for (Int_t  i=0;i<4;i++){

      char buffer[5];

      sprintf(buffer,"h%d",i);

      h[i]= new TH1F(buffer,"",100,-5,5);

      sprintf(buffer,"hw%d",i);

      hw[i] = new TH1F(buffer,"",100,-5,5);

      hw[i]->SetLineColor(3);

   }


   for (int ie=0;ie<treeS->GetEntries();ie++){

      treeS->GetEntry(ie);

      for (Int_t  i=0;i<4;i++){

         h[i]->Fill(xvar[i]);

         hw[i]->Fill(xvar[i],weight);

      }

   }


   TCanvas *c = new TCanvas("c","",800,800);

   c->Divide(2,2);


   for (Int_t  i=0;i<4;i++){

      c->cd(i+1);

      h[i]->Draw();

      hw[i]->Draw("same");

   }


   //   dataFile->Close();

   cout << "created data file: " << dataFile->GetName() << endl;

}


// create the data

void create_lin_Nvar_Arr(Int_t N = 1000)

{

   const Int_t nvar = 4;

   std::vector<float>* xvar[nvar];


   // output flie

   TFile* dataFile = TFile::Open( "data.root", "RECREATE" );


   // create signal and background trees

   TTree* treeS = new TTree( "TreeS", "TreeS", 1 );

   TTree* treeB = new TTree( "TreeB", "TreeB", 1 );

   for (Int_t ivar=0; ivar<nvar; ivar++) {

      xvar[ivar] = new std::vector<float>();

      treeS->Branch( TString(Form( "var%i", ivar+1 )).Data(), "vector<float>", &xvar[ivar], 64000, 1 );

      treeB->Branch( TString(Form( "var%i", ivar+1 )).Data(), "vector<float>", &xvar[ivar], 64000, 1 );

   }


   TRandom R( 100 );

   Float_t xS[nvar] = {  0.2,  0.3,  0.5,  0.9 };

   Float_t xB[nvar] = { -0.2, -0.3, -0.5, -0.6 };

   Float_t dx[nvar] = {  1.0,  1.0, 1.0, 1.0 };

   TArrayD* v = new TArrayD( nvar );

   Float_t rho[20];

   rho[1*2] = 0.4;

   rho[1*3] = 0.6;

   rho[1*4] = 0.9;

   rho[2*3] = 0.7;

   rho[2*4] = 0.8;

   rho[3*4] = 0.93;


   // create covariance matrix

   TMatrixD* covMatS = new TMatrixD( nvar, nvar );

   TMatrixD* covMatB = new TMatrixD( nvar, nvar );

   for (Int_t ivar=0; ivar<nvar; ivar++) {

      (*covMatS)(ivar,ivar) = dx[ivar]*dx[ivar];

      (*covMatB)(ivar,ivar) = dx[ivar]*dx[ivar];

      for (Int_t jvar=ivar+1; jvar<nvar; jvar++) {

         (*covMatS)(ivar,jvar) = rho[(ivar+1)*(jvar+1)]*dx[ivar]*dx[jvar];

         (*covMatS)(jvar,ivar) = (*covMatS)(ivar,jvar);


         (*covMatB)(ivar,jvar) = rho[(ivar+1)*(jvar+1)]*dx[ivar]*dx[jvar];

         (*covMatB)(jvar,ivar) = (*covMatB)(ivar,jvar);

      }

   }

   cout << "signal covariance matrix: " << endl;

   covMatS->Print();

   cout << "background covariance matrix: " << endl;

   covMatB->Print();


   // produce the square-root matrix

   TMatrixD* sqrtMatS = produceSqrtMat( *covMatS );

   TMatrixD* sqrtMatB = produceSqrtMat( *covMatB );


   // loop over species

   for (Int_t itype=0; itype<2; itype++) {


      Float_t*  x;

      TMatrixD* m;

      if (itype == 0) { x = xS; m = sqrtMatS; cout << "- produce signal" << endl; }

      else            { x = xB; m = sqrtMatB; cout << "- produce background" << endl; }


      // event loop

      TTree* tree = (itype==0) ? treeS : treeB;

      for (Int_t i=0; i<N; i++) {


         if (i%100 == 0) cout << "... event: " << i << " (" << N << ")" << endl;


         Int_t aSize = (Int_t)(gRandom->Rndm()*10); // size of array varies between events

         for (Int_t ivar=0; ivar<nvar; ivar++) {

            xvar[ivar]->clear();

            xvar[ivar]->reserve(aSize);

         }

         for(Int_t iA = 0; iA<aSize; iA++) {

            //for (Int_t ivar=0; ivar<nvar; ivar++) {

               getGaussRnd( *v, *m, R );

               for (Int_t ivar=0; ivar<nvar; ivar++) xvar[ivar]->push_back((*v)[ivar] + x[ivar]);

               //}

         }

         tree->Fill();

      }

   }


   // write trees

   treeS->Write();

   treeB->Write();


   treeS->Show(0);

   treeB->Show(1);


   dataFile->Close();

   cout << "created data file: " << dataFile->GetName() << endl;


   //plot();

}


// create the data

void create_lin_Nvar_double()

{

   Int_t N = 10000;

   const Int_t nvar = 4;

   Double_t xvar[nvar];

   Double_t xvarD[nvar];

   Float_t  xvarF[nvar];


   // output flie

   TFile* dataFile = TFile::Open( "data.root", "RECREATE" );


   // create signal and background trees

   TTree* treeS = new TTree( "TreeS", "TreeS", 1 );

   TTree* treeB = new TTree( "TreeB", "TreeB", 1 );

   for (Int_t ivar=0; ivar<nvar; ivar++) {

      if (ivar<2) {

         treeS->Branch( TString(Form( "var%i", ivar+1 )).Data(), &xvarD[ivar], TString(Form( "var%i/D", ivar+1 )).Data() );

         treeB->Branch( TString(Form( "var%i", ivar+1 )).Data(), &xvarD[ivar], TString(Form( "var%i/D", ivar+1 )).Data() );

      }

      else {

         treeS->Branch( TString(Form( "var%i", ivar+1 )).Data(), &xvarF[ivar], TString(Form( "var%i/F", ivar+1 )).Data() );

         treeB->Branch( TString(Form( "var%i", ivar+1 )).Data(), &xvarF[ivar], TString(Form( "var%i/F", ivar+1 )).Data() );

      }

   }


   TRandom R( 100 );

   Double_t xS[nvar] = {  0.2,  0.3,  0.5,  0.6 };

   Double_t xB[nvar] = { -0.2, -0.3, -0.5, -0.6 };

   Double_t dx[nvar] = {  1.0,  1.0, 1.0, 1.0 };

   TArrayD* v = new TArrayD( nvar );

   Double_t rho[20];

   rho[1*2] = 0.4;

   rho[1*3] = 0.6;

   rho[1*4] = 0.9;

   rho[2*3] = 0.7;

   rho[2*4] = 0.8;

   rho[3*4] = 0.93;


   // create covariance matrix

   TMatrixD* covMatS = new TMatrixD( nvar, nvar );

   TMatrixD* covMatB = new TMatrixD( nvar, nvar );

   for (Int_t ivar=0; ivar<nvar; ivar++) {

      (*covMatS)(ivar,ivar) = dx[ivar]*dx[ivar];

      (*covMatB)(ivar,ivar) = dx[ivar]*dx[ivar];

      for (Int_t jvar=ivar+1; jvar<nvar; jvar++) {

         (*covMatS)(ivar,jvar) = rho[(ivar+1)*(jvar+1)]*dx[ivar]*dx[jvar];

         (*covMatS)(jvar,ivar) = (*covMatS)(ivar,jvar);


         (*covMatB)(ivar,jvar) = rho[(ivar+1)*(jvar+1)]*dx[ivar]*dx[jvar];

         (*covMatB)(jvar,ivar) = (*covMatB)(ivar,jvar);

      }

   }

   cout << "signal covariance matrix: " << endl;

   covMatS->Print();

   cout << "background covariance matrix: " << endl;

   covMatB->Print();


   // produce the square-root matrix

   TMatrixD* sqrtMatS = produceSqrtMat( *covMatS );

   TMatrixD* sqrtMatB = produceSqrtMat( *covMatB );


   // loop over species

   for (Int_t itype=0; itype<2; itype++) {


      Double_t*  x;

      TMatrixD* m;

      if (itype == 0) { x = xS; m = sqrtMatS; cout << "- produce signal" << endl; }

      else            { x = xB; m = sqrtMatB; cout << "- produce background" << endl; }


      // event loop

      TTree* tree = (itype==0) ? treeS : treeB;

      for (Int_t i=0; i<N; i++) {


         if (i%1000 == 0) cout << "... event: " << i << " (" << N << ")" << endl;

         getGaussRnd( *v, *m, R );


         for (Int_t ivar=0; ivar<nvar; ivar++) xvar[ivar] = (*v)[ivar] + x[ivar];

         for (Int_t ivar=0; ivar<nvar; ivar++) {

            if (ivar<2) xvarD[ivar] = xvar[ivar];

            else        xvarF[ivar] = xvar[ivar];

         }


         tree->Fill();

      }

   }


   // write trees

   treeS->Write();

   treeB->Write();


   treeS->Show(0);

   treeB->Show(1);


   dataFile->Close();

   cout << "created data file: " << dataFile->GetName() << endl;


   plot();

}


// create the data

void create_lin_Nvar_discrete()

{

   Int_t N = 10000;

   const Int_t nvar = 4;

   Float_t xvar[nvar];

   Int_t   xvarI[2];


   // output flie

   TFile* dataFile = TFile::Open( "data.root", "RECREATE" );


   // create signal and background trees

   TTree* treeS = new TTree( "TreeS", "TreeS", 1 );

   TTree* treeB = new TTree( "TreeB", "TreeB", 1 );

   for (Int_t ivar=0; ivar<nvar-2; ivar++) {

      treeS->Branch( TString(Form( "var%i", ivar+1 )).Data(), &xvar[ivar], TString(Form( "var%i/F", ivar+1 )).Data() );

      treeB->Branch( TString(Form( "var%i", ivar+1 )).Data(), &xvar[ivar], TString(Form( "var%i/F", ivar+1 )).Data() );

   }

   for (Int_t ivar=0; ivar<2; ivar++) {

      treeS->Branch( TString(Form( "var%i", ivar+nvar-2+1 )).Data(), &xvarI[ivar], TString(Form( "var%i/I", ivar+nvar-2+1 )).Data() );

      treeB->Branch( TString(Form( "var%i", ivar+nvar-2+1 )).Data(), &xvarI[ivar], TString(Form( "var%i/I", ivar+nvar-2+1 )).Data() );

   }


   TRandom R( 100 );

   Float_t xS[nvar] = {  0.2,  0.3,  1,  2 };

   Float_t xB[nvar] = { -0.2, -0.3,  0,  0 };

   Float_t dx[nvar] = {  1.0,  1.0, 1, 2 };

   TArrayD* v = new TArrayD( nvar );

   Float_t rho[20];

   rho[1*2] = 0.4;

   rho[1*3] = 0.6;

   rho[1*4] = 0.9;

   rho[2*3] = 0.7;

   rho[2*4] = 0.8;

   rho[3*4] = 0.93;

   // no correlations

   for (int i=0; i<20; i++) rho[i] = 0;


   // create covariance matrix

   TMatrixD* covMatS = new TMatrixD( nvar, nvar );

   TMatrixD* covMatB = new TMatrixD( nvar, nvar );

   for (Int_t ivar=0; ivar<nvar; ivar++) {

      (*covMatS)(ivar,ivar) = dx[ivar]*dx[ivar];

      (*covMatB)(ivar,ivar) = dx[ivar]*dx[ivar];

      for (Int_t jvar=ivar+1; jvar<nvar; jvar++) {

         (*covMatS)(ivar,jvar) = rho[(ivar+1)*(jvar+1)]*dx[ivar]*dx[jvar];

         (*covMatS)(jvar,ivar) = (*covMatS)(ivar,jvar);


         (*covMatB)(ivar,jvar) = rho[(ivar+1)*(jvar+1)]*dx[ivar]*dx[jvar];

         (*covMatB)(jvar,ivar) = (*covMatB)(ivar,jvar);

      }

   }

   cout << "signal covariance matrix: " << endl;

   covMatS->Print();

   cout << "background covariance matrix: " << endl;

   covMatB->Print();


   // produce the square-root matrix

   TMatrixD* sqrtMatS = produceSqrtMat( *covMatS );

   TMatrixD* sqrtMatB = produceSqrtMat( *covMatB );


   // loop over species

   for (Int_t itype=0; itype<2; itype++) {


      Float_t*  x;

      TMatrixD* m;

      if (itype == 0) { x = xS; m = sqrtMatS; cout << "- produce signal" << endl; }

      else            { x = xB; m = sqrtMatB; cout << "- produce background" << endl; }


      // event loop

      TTree* tree = (itype==0) ? treeS : treeB;

      for (Int_t i=0; i<N; i++) {


         if (i%1000 == 0) cout << "... event: " << i << " (" << N << ")" << endl;

         getGaussRnd( *v, *m, R );


         for (Int_t ivar=0; ivar<nvar; ivar++) xvar[ivar] = (*v)[ivar] + x[ivar];


         xvarI[0] =  TMath::Nint(xvar[nvar-2]);

         xvarI[1] =  TMath::Nint(xvar[nvar-1]);


         tree->Fill();

      }

   }


   // write trees

   treeS->Write();

   treeB->Write();


   treeS->Show(0);

   treeB->Show(1);


   dataFile->Close();

   cout << "created data file: " << dataFile->GetName() << endl;


   plot();

}


// create the data

void create_ManyVars()

{

   Int_t N = 20000;

   const Int_t nvar = 20;

   Float_t xvar[nvar];


   // output flie

   TFile* dataFile = TFile::Open( "data.root", "RECREATE" );


   // create signal and background trees

   TTree* treeS = new TTree( "TreeS", "TreeS", 1 );

   TTree* treeB = new TTree( "TreeB", "TreeB", 1 );

   for (Int_t ivar=0; ivar<nvar; ivar++) {

      treeS->Branch( TString(Form( "var%i", ivar )).Data(), &xvar[ivar], TString(Form( "var%i/F", ivar )).Data() );

      treeB->Branch( TString(Form( "var%i", ivar )).Data(), &xvar[ivar], TString(Form( "var%i/F", ivar )).Data() );

   }


   Float_t xS[nvar];

   Float_t xB[nvar];

   Float_t dx[nvar];

   for (Int_t ivar=0; ivar<nvar; ivar++) {

      xS[ivar] = 0 + ivar*0.05;

      xB[ivar] = 0 - ivar*0.05;

      dx[ivar] = 1;

   }


   xS[0] =   0.2;

   xB[0] =  -0.2;

   dx[0] =   1.0;

   xS[1] =   0.3;

   xB[1] =  -0.3;

   dx[1] =   1.0;

   xS[2] =   0.4;

   xB[2] =  -0.4;

   dx[2] =  1.0 ;

   xS[3] =   0.8 ;

   xB[3] =  -0.5 ;

   dx[3] =   1.0 ;

//   TArrayD* v = new TArrayD( nvar );

   Float_t rho[20];

   rho[1*2] = 0.4;

   rho[1*3] = 0.6;

   rho[1*4] = 0.9;

   rho[2*3] = 0.7;

   rho[2*4] = 0.8;

   rho[3*4] = 0.93;


   TRandom R( 100 );


   // loop over species

   for (Int_t itype=0; itype<2; itype++) {


      Float_t* x = (itype == 0) ? xS : xB;


      // event loop

      TTree* tree = (itype == 0) ? treeS : treeB;

      for (Int_t i=0; i<N; i++) {


         if (i%1000 == 0) cout << "... event: " << i << " (" << N << ")" << endl;

         for (Int_t ivar=0; ivar<nvar; ivar++) {

            if (ivar == 1500 && itype!=10) xvar[ivar] = 1;

            else                           xvar[ivar] = x[ivar] + R.Gaus()*dx[ivar];

         }


         tree->Fill();

      }

   }


   // write trees

   treeS->Write();

   treeB->Write();


   treeS->Show(0);

   treeB->Show(1);


   dataFile->Close();

   plot();

   cout << "created data file: " << dataFile->GetName() << endl;

}


// create the data

void create_lin_NvarObsolete()

{

   Int_t N = 20000;

   const Int_t nvar = 20;

   Float_t xvar[nvar];


   // output flie

   TFile* dataFile = TFile::Open( "data.root", "RECREATE" );


   // create signal and background trees

   TTree* treeS = new TTree( "TreeS", "TreeS", 1 );

   TTree* treeB = new TTree( "TreeB", "TreeB", 1 );

   for (Int_t ivar=0; ivar<nvar; ivar++) {

      treeS->Branch( TString(Form( "var%i", ivar )).Data(), &xvar[ivar], TString(Form( "var%i/F", ivar )).Data() );

      treeB->Branch( TString(Form( "var%i", ivar )).Data(), &xvar[ivar], TString(Form( "var%i/F", ivar )).Data() );

   }


   TRandom R( 100 );

   Float_t xS[nvar] = {  0.5,  0.5,  0.0,  0.0,  0.0,  0.0 };

   Float_t xB[nvar] = { -0.5, -0.5, -0.0, -0.0, -0.0, -0.0 };

   Float_t dx[nvar] = {  1.0,  1.0, 1.0, 1.0, 1.0, 1.0 };

   TArrayD* v = new TArrayD( nvar );

   Float_t rho[50];

   for (Int_t i=0; i<50; i++) rho[i] = 0;

   rho[1*2] = 0.3;

   rho[1*3] = 0.0;

   rho[1*4] = 0.0;

   rho[2*3] = 0.0;

   rho[2*4] = 0.0;

   rho[3*4] = 0.0;


   // create covariance matrix

   TMatrixD* covMatS = new TMatrixD( nvar, nvar );

   TMatrixD* covMatB = new TMatrixD( nvar, nvar );

   for (Int_t ivar=0; ivar<nvar; ivar++) {

      (*covMatS)(ivar,ivar) = dx[ivar]*dx[ivar];

      (*covMatB)(ivar,ivar) = dx[ivar]*dx[ivar];

      for (Int_t jvar=ivar+1; jvar<nvar; jvar++) {

         (*covMatS)(ivar,jvar) = rho[(ivar+1)*(jvar+1)]*dx[ivar]*dx[jvar];

         (*covMatS)(jvar,ivar) = (*covMatS)(ivar,jvar);


         (*covMatB)(ivar,jvar) = rho[(ivar+1)*(jvar+1)]*dx[ivar]*dx[jvar];

         (*covMatB)(jvar,ivar) = (*covMatB)(ivar,jvar);

      }

   }

   cout << "signal covariance matrix: " << endl;

   covMatS->Print();

   cout << "background covariance matrix: " << endl;

   covMatB->Print();


   // produce the square-root matrix

   TMatrixD* sqrtMatS = produceSqrtMat( *covMatS );

   TMatrixD* sqrtMatB = produceSqrtMat( *covMatB );


   // loop over species

   for (Int_t itype=0; itype<2; itype++) {


      Float_t*  x;

      TMatrixD* m;

      if (itype == 0) { x = xS; m = sqrtMatS; cout << "- produce signal" << endl; }

      else            { x = xB; m = sqrtMatB; cout << "- produce background" << endl; }


      // event loop

      TTree* tree = (itype==0) ? treeS : treeB;

      for (Int_t i=0; i<N; i++) {


         if (i%1000 == 0) cout << "... event: " << i << " (" << N << ")" << endl;

         getGaussRnd( *v, *m, R );


         for (Int_t ivar=0; ivar<nvar; ivar++) xvar[ivar] = (*v)[ivar] + x[ivar];


         tree->Fill();

      }

   }


   // write trees

   treeS->Write();

   treeB->Write();


   treeS->Show(0);

   treeB->Show(1);


   dataFile->Close();

   cout << "created data file: " << dataFile->GetName() << endl;


   plot();

}


// create the data

void create_lin(Int_t N = 2000)

{

   const Int_t nvar = 2;

   Double_t xvar[nvar];

   Float_t weight;


   // output flie

   TFile* dataFile = TFile::Open( "data.root", "RECREATE" );


   // create signal and background trees

   TTree* treeS = new TTree( "TreeS", "TreeS", 1 );

   TTree* treeB = new TTree( "TreeB", "TreeB", 1 );

   for (Int_t ivar=0; ivar<nvar; ivar++) {

      treeS->Branch( TString(Form( "var%i", ivar )).Data(), &xvar[ivar], TString(Form( "var%i/D", ivar )).Data() );

      treeB->Branch( TString(Form( "var%i", ivar )).Data(), &xvar[ivar], TString(Form( "var%i/D", ivar )).Data() );

   }

   treeS->Branch( "weight", &weight, "weight/F" );

   treeB->Branch( "weight", &weight, "weight/F" );


   TRandom R( 100 );

   Float_t xS[nvar] = {  0.0,  0.0 };

   Float_t xB[nvar] = { -0.0, -0.0 };

   Float_t dx[nvar] = {  1.0,  1.0 };

   TArrayD* v = new TArrayD( 2 );

   Float_t rhoS =  0.21;

   Float_t rhoB =  0.0;


   // create covariance matrix

   TMatrixD* covMatS = new TMatrixD( nvar, nvar );

   TMatrixD* covMatB = new TMatrixD( nvar, nvar );

   for (Int_t ivar=0; ivar<nvar; ivar++) {

      (*covMatS)(ivar,ivar) = dx[ivar]*dx[ivar];

      (*covMatB)(ivar,ivar) = dx[ivar]*dx[ivar];

      for (Int_t jvar=ivar+1; jvar<nvar; jvar++) {

         (*covMatS)(ivar,jvar) = rhoS*dx[ivar]*dx[jvar];

         (*covMatS)(jvar,ivar) = (*covMatS)(ivar,jvar);


         (*covMatB)(ivar,jvar) = rhoB*dx[ivar]*dx[jvar];

         (*covMatB)(jvar,ivar) = (*covMatB)(ivar,jvar);

      }

   }

   cout << "signal covariance matrix: " << endl;

   covMatS->Print();

   cout << "background covariance matrix: " << endl;

   covMatB->Print();


   // produce the square-root matrix

   TMatrixD* sqrtMatS = produceSqrtMat( *covMatS );

   TMatrixD* sqrtMatB = produceSqrtMat( *covMatB );


   // loop over species

   for (Int_t itype=0; itype<2; itype++) {


      Float_t*  x;

      TMatrixD* m;

      if (itype == 0) { x = xS; m = sqrtMatS; cout << "- produce signal" << endl; }

      else            { x = xB; m = sqrtMatB; cout << "- produce background" << endl; }


      // event loop

      TTree* tree = (itype==0) ? treeS : treeB;

      for (Int_t i=0; i<N; i++) {


         if (i%1000 == 0) cout << "... event: " << i << " (" << N << ")" << endl;

         getGaussRnd( *v, *m, R );

         for (Int_t ivar=0; ivar<nvar; ivar++) xvar[ivar] = (*v)[ivar] + x[ivar];


         // add weights

         if (itype == 0) weight = 1.0; // this is the signal weight

         else            weight = 2.0; // this is the background weight


         tree->Fill();

      }

   }


   // write trees

   treeS->Write();

   treeB->Write();


   treeS->Show(0);

   treeB->Show(1);


   dataFile->Close();

   cout << "created data file: " << dataFile->GetName() << endl;


   plot();

}


// create the data


void create_fullcirc(Int_t nmax  = 20000,  Bool_t distort=false)

{

  TFile* dataFile = TFile::Open( "circledata.root", "RECREATE" );

   int nvar = 2;

   int nsig = 0, nbgd=0;

   Float_t weight=1;

   Float_t xvar[100];

   // create signal and background trees

   TTree* treeS = new TTree( "TreeS", "TreeS", 1 );

   TTree* treeB = new TTree( "TreeB", "TreeB", 1 );

   for (Int_t ivar=0; ivar<nvar; ivar++) {

      treeS->Branch( TString(Form( "var%i", ivar )).Data(), &xvar[ivar], TString(Form( "var%i/F", ivar)).Data() );

      treeB->Branch( TString(Form( "var%i", ivar )).Data(), &xvar[ivar], TString(Form( "var%i/F", ivar)).Data() );

   }

   treeS->Branch("weight", &weight, "weight/F");

   treeB->Branch("weight", &weight, "weight/F");


   TRandom R( 100 );

   do {

      for (Int_t ivar=0; ivar<nvar; ivar++) { xvar[ivar]=2.*R.Rndm()-1.;}

      Float_t xout = xvar[0]*xvar[0]+xvar[1]*xvar[1];

      if (nsig<10) cout << "xout = " << xout<<endl;

      if (xout < 0.3  || (xout >0.3 && xout<0.5 && R.Rndm() > xout)) {

         if (distort && xvar[0] < 0 && R.Rndm()>0.1) continue;

         treeS->Fill();

         nsig++;

      }

      else {

         if (distort && xvar[0] > 0 && R.Rndm()>0.1) continue;

         treeB->Fill();

         nbgd++;

      }

   } while ( nsig < nmax || nbgd < nmax);


   dataFile->Write();

   dataFile->Close();


}


// create the data

void create_circ(Int_t N  = 6000, Bool_t distort = false)

{

   Int_t Nn = 0;

   const Int_t nvar = 2;

   Float_t xvar[nvar];


   // output flie

   TFile* dataFile = TFile::Open( "data.root", "RECREATE" );


   // create signal and background trees

   TTree* treeS = new TTree( "TreeS", "TreeS", 1 );

   TTree* treeB = new TTree( "TreeB", "TreeB", 1 );

   for (Int_t ivar=0; ivar<nvar; ivar++) {

      treeS->Branch( TString(Form( "var%i", ivar )).Data(), &xvar[ivar], TString(Form( "var%i/F", ivar )).Data() );

      treeB->Branch( TString(Form( "var%i", ivar )).Data(), &xvar[ivar], TString(Form( "var%i/F", ivar )).Data() );

   }

//    TTree *treeB  = treeS->CloneTree();

//    for (Int_t ivar=0; ivar<nvar; ivar++) {

//       treeS->SetBranchAddress( Form( "var%i", ivar ), &xvar[ivar] );

//       treeB->SetBranchAddress( Form( "var%i", ivar ), &xvar[ivar] );

//    }

//    treeB->SetName ( "TreeB" );

//    treeB->SetTitle( "TreeB" );


   TRandom R( 100 );

   //Float_t phimin = -30, phimax = 130;

   Float_t phimin = -70, phimax = 130;

   Float_t phisig = 5;

   Float_t rS = 1.1;

   Float_t rB = 0.75;

   Float_t rsig = 0.1;

   Float_t fnmin = -(rS+4.0*rsig);

   Float_t fnmax = +(rS+4.0*rsig);

   Float_t dfn = fnmax-fnmin;

   // loop over species

   for (Int_t itype=0; itype<2; itype++) {


      // event loop

      TTree* tree = (itype==0) ? treeS : treeB;

      for (Int_t i=0; i<N; i++) {

    Double_t r1=R.Rndm(),r2=R.Rndm(), r3;

    if (itype==0) r3= r1>r2? r1 :r2;

    else r3= r2;

    Float_t phi;

    if (distort) phi = r3*(phimax - phimin) + phimin;

    else  phi = R.Rndm()*(phimax - phimin) + phimin;

         phi += R.Gaus()*phisig;


         Float_t r = (itype==0) ? rS : rB;

         r += R.Gaus()*rsig;


         xvar[0] = r*cos(TMath::DegToRad()*phi);

         xvar[1] = r*sin(TMath::DegToRad()*phi);


         tree->Fill();

      }


      for (Int_t i=0; i<Nn; i++) {


         xvar[0] = dfn*R.Rndm()+fnmin;

         xvar[1] = dfn*R.Rndm()+fnmin;


         tree->Fill();

      }

   }


   // write trees

   treeS->Write();

   treeB->Write();


   treeS->Show(0);

   treeB->Show(1);


   dataFile->Close();

   cout << "created data file: " << dataFile->GetName() << endl;


   plot();

}


void create_schachbrett(Int_t nEvents = 20000) {


   const Int_t nvar = 2;

   Float_t xvar[nvar];


   // output flie

   TFile* dataFile = TFile::Open( "data.root", "RECREATE" );


   // create signal and background trees

   TTree* treeS = new TTree( "TreeS", "TreeS", 1 );

   TTree* treeB = new TTree( "TreeB", "TreeB", 1 );

   for (Int_t ivar=0; ivar<nvar; ivar++) {

      treeS->Branch( TString(Form( "var%i", ivar )).Data(), &xvar[ivar], TString(Form( "var%i/F", ivar )).Data() );

      treeB->Branch( TString(Form( "var%i", ivar )).Data(), &xvar[ivar], TString(Form( "var%i/F", ivar )).Data() );

   }


   Int_t   nSeed   = 12345;

   TRandom *m_rand = new TRandom(nSeed);

   Double_t sigma=0.3;

   Double_t meanX;

   Double_t meanY;

   Int_t xtype=1, ytype=1;

   Int_t iev=0;

   Int_t m_nDim = 2; // actually the boundary, there is a "bump" for every interger value

                     // between in the Inteval [-m_nDim,m_nDim]

   while (iev < nEvents){

      xtype=1;

      for (Int_t i=-m_nDim; i <=  m_nDim; i++){

         ytype  =  1;

         for (Int_t j=-m_nDim; j <=  m_nDim; j++){

            meanX=Double_t(i);

            meanY=Double_t(j);

            xvar[0]=m_rand->Gaus(meanY,sigma);

            xvar[1]=m_rand->Gaus(meanX,sigma);

            Int_t type   = xtype*ytype;

            TTree* tree = (type==1) ? treeS : treeB;

            tree->Fill();

            iev++;

            ytype *= -1;

         }

         xtype *= -1;

      }

   }


   // write trees

   treeS->Write();

   treeB->Write();


   treeS->Show(0);

   treeB->Show(1);


   dataFile->Close();

   cout << "created data file: " << dataFile->GetName() << endl;


   plot();


}


void create_schachbrett_5D(Int_t nEvents = 200000) {

   const Int_t nvar = 5;

   Float_t xvar[nvar];


   // output flie

   TFile* dataFile = TFile::Open( "data.root", "RECREATE" );


   // create signal and background trees

   TTree* treeS = new TTree( "TreeS", "TreeS", 1 );

   TTree* treeB = new TTree( "TreeB", "TreeB", 1 );

   for (Int_t ivar=0; ivar<nvar; ivar++) {

      treeS->Branch( TString(Form( "var%i", ivar )).Data(), &xvar[ivar], TString(Form( "var%i/F", ivar )).Data() );

      treeB->Branch( TString(Form( "var%i", ivar )).Data(), &xvar[ivar], TString(Form( "var%i/F", ivar )).Data() );

   }


   Int_t   nSeed   = 12345;

   TRandom *m_rand = new TRandom(nSeed);

   Double_t sigma=0.3;

   Int_t itype[nvar];

   Int_t iev=0;

   Int_t m_nDim = 2; // actually the boundary, there is a "bump" for every interger value

                     // between in the Inteval [-m_nDim,m_nDim]


   int idx[nvar];

   while (iev < nEvents){

      itype[0]=1;

      for (idx[0]=-m_nDim; idx[0] <=  m_nDim; idx[0]++){

         itype[1]=1;

         for (idx[1]=-m_nDim; idx[1] <=  m_nDim; idx[1]++){

            itype[2]=1;

            for (idx[2]=-m_nDim; idx[2] <=  m_nDim; idx[2]++){

               itype[3]=1;

               for (idx[3]=-m_nDim; idx[3] <=  m_nDim; idx[3]++){

                  itype[4]=1;

                  for (idx[4]=-m_nDim; idx[4] <=  m_nDim; idx[4]++){

                     Int_t type   = itype[0];

                     for (Int_t i=0;i<nvar;i++){

                        xvar[i]=m_rand->Gaus(Double_t(idx[i]),sigma);

                        if (i>0) type *= itype[i];

                     }

                     TTree* tree = (type==1) ? treeS : treeB;

                     tree->Fill();

                     iev++;

                     itype[4] *= -1;

                  }

                  itype[3] *= -1;

               }

               itype[2] *= -1;

            }

            itype[1] *= -1;

         }

         itype[0] *= -1;

      }

   }


   // write trees

   treeS->Write();

   treeB->Write();


   treeS->Show(0);

   treeB->Show(1);


   dataFile->Close();

   cout << "created data file: " << dataFile->GetName() << endl;


   plot();


}


void create_schachbrett_4D(Int_t nEvents = 200000) {


   const Int_t nvar = 4;

   Float_t xvar[nvar];


   // output flie

   TFile* dataFile = TFile::Open( "data.root", "RECREATE" );


   // create signal and background trees

   TTree* treeS = new TTree( "TreeS", "TreeS", 1 );

   TTree* treeB = new TTree( "TreeB", "TreeB", 1 );

   for (Int_t ivar=0; ivar<nvar; ivar++) {

      treeS->Branch( TString(Form( "var%i", ivar )).Data(), &xvar[ivar], TString(Form( "var%i/F", ivar )).Data() );

      treeB->Branch( TString(Form( "var%i", ivar )).Data(), &xvar[ivar], TString(Form( "var%i/F", ivar )).Data() );

   }


   Int_t   nSeed   = 12345;

   TRandom *m_rand = new TRandom(nSeed);

   Double_t sigma=0.3;

   Int_t itype[nvar];

   Int_t iev=0;

   Int_t m_nDim = 2; // actually the boundary, there is a "bump" for every interger value

                     // between in the Inteval [-m_nDim,m_nDim]


   int idx[nvar];

   while (iev < nEvents){

      itype[0]=1;

      for (idx[0]=-m_nDim; idx[0] <=  m_nDim; idx[0]++){

         itype[1]=1;

         for (idx[1]=-m_nDim; idx[1] <=  m_nDim; idx[1]++){

            itype[2]=1;

            for (idx[2]=-m_nDim; idx[2] <=  m_nDim; idx[2]++){

               itype[3]=1;

               for (idx[3]=-m_nDim; idx[3] <=  m_nDim; idx[3]++){

                  Int_t type   = itype[0];

                  for (Int_t i=0;i<nvar;i++){

                     xvar[i]=m_rand->Gaus(Double_t(idx[i]),sigma);

                     if (i>0) type *= itype[i];

                  }

                  TTree* tree = (type==1) ? treeS : treeB;

                  tree->Fill();

                  iev++;

                  itype[3] *= -1;

               }

               itype[2] *= -1;

            }

            itype[1] *= -1;

         }

         itype[0] *= -1;

      }

   }


   // write trees

   treeS->Write();

   treeB->Write();


   treeS->Show(0);

   treeB->Show(1);


   dataFile->Close();

   cout << "created data file: " << dataFile->GetName() << endl;


   plot();


}


void create_schachbrett_3D(Int_t nEvents = 20000) {


   const Int_t nvar = 3;

   Float_t xvar[nvar];


   // output flie

   TFile* dataFile = TFile::Open( "data.root", "RECREATE" );


   // create signal and background trees

   TTree* treeS = new TTree( "TreeS", "TreeS", 1 );

   TTree* treeB = new TTree( "TreeB", "TreeB", 1 );

   for (Int_t ivar=0; ivar<nvar; ivar++) {

      treeS->Branch( TString(Form( "var%i", ivar )).Data(), &xvar[ivar], TString(Form( "var%i/F", ivar )).Data() );

      treeB->Branch( TString(Form( "var%i", ivar )).Data(), &xvar[ivar], TString(Form( "var%i/F", ivar )).Data() );

   }


   Int_t   nSeed   = 12345;

   TRandom *m_rand = new TRandom(nSeed);

   Double_t sigma=0.3;

   Int_t itype[nvar];

   Int_t iev=0;

   Int_t m_nDim = 2; // actually the boundary, there is a "bump" for every interger value

                     // between in the Inteval [-m_nDim,m_nDim]


   int idx[nvar];

   while (iev < nEvents){

      itype[0]=1;

      for (idx[0]=-m_nDim; idx[0] <=  m_nDim; idx[0]++){

         itype[1]=1;

         for (idx[1]=-m_nDim; idx[1] <=  m_nDim; idx[1]++){

            itype[2]=1;

            for (idx[2]=-m_nDim; idx[2] <=  m_nDim; idx[2]++){

               Int_t type   = itype[0];

               for (Int_t i=0;i<nvar;i++){

                  xvar[i]=m_rand->Gaus(Double_t(idx[i]),sigma);

                  if (i>0) type *= itype[i];

               }

               TTree* tree = (type==1) ? treeS : treeB;

               tree->Fill();

               iev++;

               itype[2] *= -1;

            }

            itype[1] *= -1;

         }

         itype[0] *= -1;

      }

   }


   // write trees

   treeS->Write();

   treeB->Write();


   treeS->Show(0);

   treeB->Show(1);


   dataFile->Close();

   cout << "created data file: " << dataFile->GetName() << endl;


   plot();


}


void create_schachbrett_2D(Int_t nEvents = 100000, Int_t nbumps=2) {


   const Int_t nvar = 2;

   Float_t xvar[nvar];


   // output flie

   TFile* dataFile = TFile::Open( "data.root", "RECREATE" );


   // create signal and background trees

   TTree* treeS = new TTree( "TreeS", "TreeS", 1 );

   TTree* treeB = new TTree( "TreeB", "TreeB", 1 );

   for (Int_t ivar=0; ivar<nvar; ivar++) {

      treeS->Branch( TString(Form( "var%i", ivar )).Data(), &xvar[ivar], TString(Form( "var%i/F", ivar )).Data() );

      treeB->Branch( TString(Form( "var%i", ivar )).Data(), &xvar[ivar], TString(Form( "var%i/F", ivar )).Data() );

   }


   Int_t   nSeed   = 345;

   TRandom *m_rand = new TRandom(nSeed);

   Double_t sigma=0.35;

   Int_t itype[nvar];

   Int_t iev=0;

   Int_t m_nDim = nbumps; // actually the boundary, there is a "bump" for every interger value

                     // between in the Inteval [-m_nDim,m_nDim]


   int idx[nvar];

   while (iev < nEvents){

      itype[0]=1;

      for (idx[0]=-m_nDim; idx[0] <=  m_nDim; idx[0]++){

         itype[1]=1;

         for (idx[1]=-m_nDim; idx[1] <=  m_nDim; idx[1]++){

            Int_t type   = itype[0];

            for (Int_t i=0;i<nvar;i++){

               xvar[i]=m_rand->Gaus(Double_t(idx[i]),sigma);

               if (i>0) type *= itype[i];

            }

            TTree* tree = (type==1) ? treeS : treeB;

            tree->Fill();

            iev++;

            itype[1] *= -1;

         }

         itype[0] *= -1;

      }

   }


   // write trees

   treeS->Write();

   treeB->Write();


   treeS->Show(0);

   treeB->Show(1);


   dataFile->Close();

   cout << "created data file: " << dataFile->GetName() << endl;


   plot();


}


void create_3Bumps(Int_t nEvents = 5000) {

   // signal is clustered around (1,0) and (-1,0) where one is two times(1,0)

   // bkg                        (0,0)


   const Int_t nvar = 2;

   Float_t xvar[nvar];


   // output flie

   TString filename = "data_3Bumps.root";

   TFile* dataFile = TFile::Open( filename, "RECREATE" );


   // create signal and background trees

   TTree* treeS = new TTree( "TreeS", "TreeS", 1 );

   TTree* treeB = new TTree( "TreeB", "TreeB", 1 );

   for (Int_t ivar=0; ivar<nvar; ivar++) {

      treeS->Branch( TString(Form( "var%i", ivar )).Data(), &xvar[ivar], TString(Form( "var%i/F", ivar )).Data() );

      treeB->Branch( TString(Form( "var%i", ivar )).Data(), &xvar[ivar], TString(Form( "var%i/F", ivar )).Data() );

   }


   Int_t   nSeed   = 12345;

   TRandom *m_rand = new TRandom(nSeed);

   Double_t sigma=0.2;

   Int_t type;

   Int_t iev=0;

   Double_t Centers[nvar][6] = {{-1,0,0,0,1,1},{0,0,0,0,0,0}}; //


   while (iev < nEvents){

      for (int idx=0; idx<6; idx++){

         if (idx==1 || idx==2 || idx==3) type = 0;

         else type=1;

         for (Int_t ivar=0;ivar<nvar;ivar++){

            xvar[ivar]=m_rand->Gaus(Centers[ivar][idx],sigma);

         }

         TTree* tree = (type==1) ? treeS : treeB;

         tree->Fill();

         iev++;

      }

   }


   // write trees

   treeS->Write();

   treeB->Write();


   treeS->Show(0);

   treeB->Show(1);


   dataFile->Close();

   cout << "created data file: " << dataFile->GetName() << endl;


   plot(filename);


}


void createOnionData(Int_t nmax = 50000){

   // output file

   TFile* dataFile = TFile::Open( "oniondata.root", "RECREATE" );

   int nvar = 4;

   int nsig = 0, nbgd=0;

   Float_t xvar[100];

   // create signal and background trees

   TTree* treeS = new TTree( "TreeS", "TreeS", 1 );

   TTree* treeB = new TTree( "TreeB", "TreeB", 1 );

   for (Int_t ivar=0; ivar<nvar; ivar++) {

      treeS->Branch( TString(Form( "var%i", ivar+1 )).Data(), &xvar[ivar], TString(Form( "var%i/F", ivar+1 )).Data() );

      treeB->Branch( TString(Form( "var%i", ivar+1 )).Data(), &xvar[ivar], TString(Form( "var%i/F", ivar+1 )).Data() );

   }


   TRandom R( 100 );

   do {

      for (Int_t ivar=0; ivar<nvar; ivar++) { xvar[ivar]=R.Rndm();}

      Float_t xout = sin(2.*acos(-1.)*(xvar[0]*xvar[1]*xvar[2]*xvar[3]+xvar[0]*xvar[1]));

      if (nsig<100) cout << "xout = " << xout<<endl;

      Int_t i = (Int_t) ((1.+xout)*4.99);

      if (i%2 == 0 && nsig < nmax) {

    treeS->Fill();

    nsig++;

      }

      if (i%2 != 0 && nbgd < nmax){

    treeB->Fill();

    nbgd++;

      }

   } while ( nsig < nmax || nbgd < nmax);


   dataFile->Write();

   dataFile->Close();

}


void create_multiclassdata(Int_t nmax  = 20000)

{

   TFile* dataFile = TFile::Open( "data.root", "RECREATE" );

   int ncls = 3;

   int nvar = 4;

   int ndat = 0;

   Int_t cls;

   Float_t thecls;

   Float_t weight=1;

   Float_t xcls[100];

   Float_t xmean[3][4] = {

      { 0.   ,  0.3,  0.5, 0.9 },

      { -0.2 , -0.3,  0.5, 0.4 },

      { 0.2  ,  0.1, -0.1, 0.7 }} ;


   Float_t xvar[100];

   // create tree using class flag stored in int variable cls

   TTree* treeR = new TTree( "TreeR", "TreeR", 1 );

   for (Int_t ivar=0; ivar<nvar; ivar++) {

      treeR->Branch( TString(Form( "var%i", ivar )).Data(), &xvar[ivar], TString(Form( "var%i/F", ivar)).Data() );

   }

   for (Int_t icls=0; icls<ncls; icls++) {

      treeR->Branch(TString(Form( "cls%i", icls )).Data(), &xcls[icls], TString(Form( "cls%i/F", icls)).Data() );

   }


   treeR->Branch("cls", &thecls, "cls/F");

   treeR->Branch("weight", &weight, "weight/F");


   TRandom R( 100 );

   do {

      for (Int_t icls=0; icls<ncls; icls++) xcls[icls]=0.;

      cls = R.Integer(ncls);

      thecls = cls;

      xcls[cls]=1.;

      for (Int_t ivar=0; ivar<nvar; ivar++) {

         xvar[ivar]=R.Gaus(xmean[cls][ivar],1.);

      }


      if (ndat<30) cout << "cls=" << cls <<" xvar = " << xvar[0]<<" " <<xvar[1]<<" " << xvar[2]<<" " <<xvar[3]<<endl;


      treeR->Fill();

      ndat++;

   } while ( ndat < nmax );


   dataFile->Write();

   dataFile->Close();


}


// create the data

void create_array_with_different_lengths(Int_t N = 100)

{

   const Int_t nvar = 4;

   Int_t nvarCurrent = 4;

   Float_t xvar[nvar];


   // output flie

   TFile* dataFile = TFile::Open( "data.root", "RECREATE" );


   // create signal and background trees

   TTree* treeS = new TTree( "TreeS", "TreeS", 1 );

   TTree* treeB = new TTree( "TreeB", "TreeB", 1 );

   treeS->Branch( "arrSize", &nvarCurrent, "arrSize/I" );

   treeS->Branch( "arr", xvar, "arr[arrSize]/F" );

   treeB->Branch( "arrSize", &nvarCurrent, "arrSize/I" );

   treeB->Branch( "arr", xvar, "arr[arrSize]/F" );


   TRandom R( 100 );

   Float_t xS[nvar] = {  0.2,  0.3,  0.5,  0.9 };

   Float_t xB[nvar] = { -0.2, -0.3, -0.5, -0.6 };

   Float_t dx[nvar] = {  1.0,  1.0, 1.0, 1.0 };

   TArrayD* v = new TArrayD( nvar );

   Float_t rho[20];

   rho[1*2] = 0.4;

   rho[1*3] = 0.6;

   rho[1*4] = 0.9;

   rho[2*3] = 0.7;

   rho[2*4] = 0.8;

   rho[3*4] = 0.93;


   // create covariance matrix

   TMatrixD* covMatS = new TMatrixD( nvar, nvar );

   TMatrixD* covMatB = new TMatrixD( nvar, nvar );

   for (Int_t ivar=0; ivar<nvar; ivar++) {

      (*covMatS)(ivar,ivar) = dx[ivar]*dx[ivar];

      (*covMatB)(ivar,ivar) = dx[ivar]*dx[ivar];

      for (Int_t jvar=ivar+1; jvar<nvar; jvar++) {

         (*covMatS)(ivar,jvar) = rho[(ivar+1)*(jvar+1)]*dx[ivar]*dx[jvar];

         (*covMatS)(jvar,ivar) = (*covMatS)(ivar,jvar);


         (*covMatB)(ivar,jvar) = rho[(ivar+1)*(jvar+1)]*dx[ivar]*dx[jvar];

         (*covMatB)(jvar,ivar) = (*covMatB)(ivar,jvar);

      }

   }

   cout << "signal covariance matrix: " << endl;

   covMatS->Print();

   cout << "background covariance matrix: " << endl;

   covMatB->Print();


   // produce the square-root matrix

   TMatrixD* sqrtMatS = produceSqrtMat( *covMatS );

   TMatrixD* sqrtMatB = produceSqrtMat( *covMatB );


   // loop over species

   for (Int_t itype=0; itype<2; itype++) {


      Float_t*  x;

      TMatrixD* m;

      if (itype == 0) { x = xS; m = sqrtMatS; cout << "- produce signal" << endl; }

      else            { x = xB; m = sqrtMatB; cout << "- produce background" << endl; }


      // event loop

      TTree* tree = (itype==0) ? treeS : treeB;

      for (Int_t i=0; i<N; i++) {


         if (i%1000 == 0) cout << "... event: " << i << " (" << N << ")" << endl;

         getGaussRnd( *v, *m, R );


         for (Int_t ivar=0; ivar<nvar; ivar++) xvar[ivar] = (*v)[ivar] + x[ivar];


    nvarCurrent = (i%4)+1;


         tree->Fill();

      }

   }


   // write trees

   treeS->Write();

   treeB->Write();


   treeS->Show(0);

   treeB->Show(1);


   dataFile->Close();

   cout << "created data file: " << dataFile->GetName() << endl;

}


// create the data

void create_MultipleBackground(Int_t N = 50000)

{

   const int nvar = 4;


   // output flie

   TFile* dataFile = TFile::Open( "tmva_example_multiple_background.root", "RECREATE" );


   Float_t xS[nvar] = {  0.2,  0.3,  0.5,  0.9 };

   Float_t xB0[nvar] = { -0.2, -0.3, -0.5, -0.6 };

   Float_t xB1[nvar] = { -0.2, 0.3, 0.5, -0.6 };

   Float_t dx0[nvar] = {  1.0,  1.0, 1.0, 1.0 };

   Float_t dx1[nvar] = {  -1.0,  -1.0, -1.0, -1.0 };


   // create signal and background trees

   TTree* treeS = makeTree_lin_Nvar( "TreeS", "Signal tree", xS, dx0, nvar, N );

   TTree* treeB0 = makeTree_lin_Nvar( "TreeB0", "Background 0", xB0, dx0, nvar, N );

   TTree* treeB1 = makeTree_lin_Nvar( "TreeB1", "Background 1", xB1, dx1, nvar, N );

   TTree* treeB2 = makeTree_circ( "TreeB2", "Background 2", nvar, N, 1.5, true);


   treeS->Write();

   treeB0->Write();

   treeB1->Write();

   treeB2->Write();


   //treeS->Show(0);

   //treeB0->Show(0);

   //treeB1->Show(0);

   //treeB2->Show(0);


   dataFile->Close();

   cout << "created data file: " << dataFile->GetName() << endl;

}

r
ROOT::R::TRInterface & r
Definition Object.C:4

c
#define c(i)
Definition RSha256.hxx:101

h
#define h(i)
Definition RSha256.hxx:106

R
#define R(a, b, c, d, e, f, g, h, i)
Definition RSha256.hxx:110

Riostream.h

size
size_t size(const MatrixT &matrix)
retrieve the size of a square matrix

Int_t
int Int_t
Definition RtypesCore.h:45

Bool_t
bool Bool_t
Definition RtypesCore.h:63

Double_t
double Double_t
Definition RtypesCore.h:59

Float_t
float Float_t
Definition RtypesCore.h:57

TArrayD.h

TBranch.h

TCanvas.h

TFile.h

N
#define N

type
int type
Definition TGX11.cxx:121

TH2F.h

xmin
float xmin
Definition THbookFile.cxx:95

ymin
float ymin
Definition THbookFile.cxx:95

xmax
float xmax
Definition THbookFile.cxx:95

ymax
float ymax
Definition THbookFile.cxx:95

TLegend.h

TMath.h

TMatrixD.h

TMatrixD
TMatrixT< Double_t > TMatrixD
Definition TMatrixDfwd.h:23

TROOT.h

gROOT
#define gROOT
Definition TROOT.h:404

TRandom.h

gRandom
R__EXTERN TRandom * gRandom
Definition TRandom.h:62

Form
char * Form(const char *fmt,...)

TStyle.h

TTree.h

R

TArrayD
Array of doubles (64 bits per element).
Definition TArrayD.h:27

TAttFill::SetFillStyle
virtual void SetFillStyle(Style_t fstyle)
Set the fill area style.
Definition TAttFill.h:39

TAttLine::SetLineColor
virtual void SetLineColor(Color_t lcolor)
Set the line color.
Definition TAttLine.h:40

TAttMarker::SetMarkerColor
virtual void SetMarkerColor(Color_t mcolor=1)
Set the marker color.
Definition TAttMarker.h:38

TAttMarker::SetMarkerStyle
virtual void SetMarkerStyle(Style_t mstyle=1)
Set the marker style.
Definition TAttMarker.h:40

TAttMarker::SetMarkerSize
virtual void SetMarkerSize(Size_t msize=1)
Set the marker size.
Definition TAttMarker.h:41

TCanvas
The Canvas class.
Definition TCanvas.h:23

TDirectoryFile::Get
TObject * Get(const char *namecycle) override
Return pointer to object identified by namecycle.
Definition TDirectoryFile.cxx:937

TFile
A ROOT file is a suite of consecutive data records (TKey instances) with a well defined format.
Definition TFile.h:54

TFile::Open
static TFile * Open(const char *name, Option_t *option="", const char *ftitle="", Int_t compress=ROOT::RCompressionSetting::EDefaults::kUseCompiledDefault, Int_t netopt=0)
Create / open a file.
Definition TFile.cxx:4025

TFile::Write
Int_t Write(const char *name=nullptr, Int_t opt=0, Int_t bufsiz=0) override
Write memory objects to this file.
Definition TFile.cxx:2374

TFile::Close
void Close(Option_t *option="") override
Close a file.
Definition TFile.cxx:899

TH1F
1-D histogram with a float per channel (see TH1 documentation)}
Definition TH1.h:575

TH1::SetTitle
virtual void SetTitle(const char *title)
See GetStatOverflows for more information.
Definition TH1.cxx:6667

TH1::SetTitleSize
virtual void SetTitleSize(Float_t size=0.02, Option_t *axis="X")
Set the axis' title size.
Definition Haxis.cxx:365

TH1::GetXaxis
TAxis * GetXaxis()
Get the behaviour adopted by the object about the statoverflows. See EStatOverflows for more informat...
Definition TH1.h:320

TH1::Fill
virtual Int_t Fill(Double_t x)
Increment bin with abscissa X by 1.
Definition TH1.cxx:3351

TH1::GetYaxis
TAxis * GetYaxis()
Definition TH1.h:321

TH1::Draw
virtual void Draw(Option_t *option="")
Draw this histogram with options.
Definition TH1.cxx:3074

TH1::SetLabelSize
virtual void SetLabelSize(Float_t size=0.02, Option_t *axis="X")
Set size of axis' labels.
Definition Haxis.cxx:285

TH2F
2-D histogram with a float per channel (see TH1 documentation)}
Definition TH2.h:251

TLegend
This class displays a legend box (TPaveText) containing several legend entries.
Definition TLegend.h:23

TLegend::AddEntry
TLegendEntry * AddEntry(const TObject *obj, const char *label="", Option_t *option="lpf")
Add a new entry to this legend.
Definition TLegend.cxx:330

TLegend::Draw
virtual void Draw(Option_t *option="")
Draw this legend with its current attributes.
Definition TLegend.cxx:423

TLegend::SetMargin
void SetMargin(Float_t margin)
Definition TLegend.h:69

TMatrixTBase::GetNrows
Int_t GetNrows() const
Definition TMatrixTBase.h:123

TMatrixTBase::Print
void Print(Option_t *name="") const
Print the matrix as a table of elements.
Definition TMatrixTBase.cxx:636

TMatrixT< Double_t >

TNamed::SetTitle
virtual void SetTitle(const char *title="")
Set the title of the TNamed.
Definition TNamed.cxx:164

TNamed::GetName
virtual const char * GetName() const
Returns name of object.
Definition TNamed.h:47

TPave::SetBorderSize
virtual void SetBorderSize(Int_t bordersize=4)
Definition TPave.h:73

TRandom
This is the base class for the ROOT Random number generators.
Definition TRandom.h:27

TRandom::Gaus
virtual Double_t Gaus(Double_t mean=0, Double_t sigma=1)
Samples a random number from the standard Normal (Gaussian) Distribution with the given mean and sigm...
Definition TRandom.cxx:274

TRandom::Rndm
virtual Double_t Rndm()
Machine independent random number generator.
Definition TRandom.cxx:552

TString
Basic string class.
Definition TString.h:136

TString::Data
const char * Data() const
Definition TString.h:369

TStyle
TStyle objects may be created to define special styles.
Definition TStyle.h:29

TStyle::SetOptTitle
void SetOptTitle(Int_t tit=1)
Definition TStyle.h:318

TStyle::SetPadTopMargin
void SetPadTopMargin(Float_t margin=0.1)
Definition TStyle.h:342

TStyle::SetTitleX
void SetTitleX(Float_t x=0)
Definition TStyle.h:396

TStyle::SetOptStat
void SetOptStat(Int_t stat=1)
The type of information printed in the histogram statistics box can be selected via the parameter mod...
Definition TStyle.cxx:1589

TStyle::SetPadBottomMargin
void SetPadBottomMargin(Float_t margin=0.1)
Definition TStyle.h:341

TStyle::SetPadRightMargin
void SetPadRightMargin(Float_t margin=0.1)
Definition TStyle.h:344

TStyle::SetPadGridX
void SetPadGridX(Bool_t gridx)
Definition TStyle.h:345

TStyle::SetPadLeftMargin
void SetPadLeftMargin(Float_t margin=0.1)
Definition TStyle.h:343

TStyle::SetPadGridY
void SetPadGridY(Bool_t gridy)
Definition TStyle.h:346

TStyle::cd
virtual void cd()
Change current style.
Definition TStyle.cxx:527

TStyle::SetTitleW
void SetTitleW(Float_t w=0)
Definition TStyle.h:398

TStyle::SetTitleH
void SetTitleH(Float_t h=0)
Definition TStyle.h:399

TStyle::SetTitleY
void SetTitleY(Float_t y=0.985)
Definition TStyle.h:397

TTree
A TTree represents a columnar dataset.
Definition TTree.h:79

TTree::Fill
virtual Int_t Fill()
Fill all branches.
Definition TTree.cxx:4594

TTree::Show
virtual void Show(Long64_t entry=-1, Int_t lenmax=20)
Print values of all active leaves for entry.
Definition TTree.cxx:9325

TTree::GetEntry
virtual Int_t GetEntry(Long64_t entry, Int_t getall=0)
Read all branches of entry and return total number of bytes read.
Definition TTree.cxx:5622

TTree::GetMaximum
virtual Double_t GetMaximum(const char *columname)
Return maximum of column with name columname.
Definition TTree.cxx:6219

TTree::GetEntries
virtual Long64_t GetEntries() const
Definition TTree.h:460

TTree::Branch
TBranch * Branch(const char *name, T *obj, Int_t bufsize=32000, Int_t splitlevel=99)
Add a new branch, and infer the data type from the type of obj being passed.
Definition TTree.h:350

TTree::GetMinimum
virtual Double_t GetMinimum(const char *columname)
Return minimum of column with name columname.
Definition TTree.cxx:6259

TTree::Draw
virtual void Draw(Option_t *opt)
Default Draw method for all objects.
Definition TTree.h:428

TTree::Write
virtual Int_t Write(const char *name=0, Int_t option=0, Int_t bufsize=0)
Write this object to the current directory.
Definition TTree.cxx:9693

sigma
const Double_t sigma
Definition h1analysisProxy.h:11

y
Double_t y[n]
Definition legend1.C:17

x
Double_t x[n]
Definition legend1.C:17

TMath::Gaus
Double_t Gaus(Double_t x, Double_t mean=0, Double_t sigma=1, Bool_t norm=kFALSE)
Calculate a gaussian function with mean and sigma.
Definition TMath.cxx:448

TMath::Nint
Int_t Nint(T x)
Round to nearest integer. Rounds half integers to the nearest even integer.
Definition TMath.h:663

TMath::Max
Short_t Max(Short_t a, Short_t b)
Definition TMathBase.h:208

TMath::Log
Double_t Log(Double_t x)
Definition TMath.h:710

TMath::DegToRad
constexpr Double_t DegToRad()
Conversion from degree to radian:
Definition TMath.h:81

TMath::Sqrt
Double_t Sqrt(Double_t x)
Definition TMath.h:641

TMath::Min
Short_t Min(Short_t a, Short_t b)
Definition TMathBase.h:176

TMath::Pi
constexpr Double_t Pi()
Definition TMath.h:37

TMath::Sin
Double_t Sin(Double_t)
Definition TMath.h:589

TMath::Abs
Short_t Abs(Short_t d)
Definition TMathBase.h:120

tree
Definition tree.py:1

v
@ v
Definition rootcling_impl.cxx:3670

m
auto * m
Definition textangle.C:8

l
auto * l
Definition textangle.C:4

sum
static uint64_t sum(uint64_t i)
Definition Factory.cxx:2345