Logo ROOT   6.30.04
Reference Guide
 All Namespaces Files Pages
createData.C
Go to the documentation of this file.
1 /// \file
2 /// \ingroup tutorial_tmva
3 /// Plot the variables.
4 ///
5 /// \macro_code
6 /// \author Andreas Hoecker
7 
8 
9 #include "TROOT.h"
10 #include "TMath.h"
11 #include "TTree.h"
12 #include "TArrayD.h"
13 #include "TStyle.h"
14 #include "TFile.h"
15 #include "TRandom.h"
16 #include "Riostream.h"
17 #include "TCanvas.h"
18 #include "TMatrixD.h"
19 #include "TH2F.h"
20 #include "TLegend.h"
21 #include "TBranch.h"
22 #include <vector>
23 
24 void plot( TString fname = "data.root", TString var0="var0", TString var1="var1" )
25 {
26  TFile* dataFile = TFile::Open( fname );
27 
28  if (!dataFile) {
29  cout << "ERROR: cannot open file: " << fname << endl;
30  return;
31  }
32 
33  TTree *treeS = (TTree*)dataFile->Get("TreeS");
34  TTree *treeB = (TTree*)dataFile->Get("TreeB");
35 
36  TCanvas* c = new TCanvas( "c", "", 0, 0, 550, 550 );
37 
38  TStyle *TMVAStyle = gROOT->GetStyle("Plain"); // our style is based on Plain
39  TMVAStyle->SetOptStat(0);
40  TMVAStyle->SetPadTopMargin(0.02);
41  TMVAStyle->SetPadBottomMargin(0.16);
42  TMVAStyle->SetPadRightMargin(0.03);
43  TMVAStyle->SetPadLeftMargin(0.15);
44  TMVAStyle->SetPadGridX(0);
45  TMVAStyle->SetPadGridY(0);
46 
47  TMVAStyle->SetOptTitle(0);
48  TMVAStyle->SetTitleW(.4);
49  TMVAStyle->SetTitleH(.10);
50  TMVAStyle->SetTitleX(.5);
51  TMVAStyle->SetTitleY(.9);
52  TMVAStyle->SetMarkerStyle(20);
53  TMVAStyle->SetMarkerSize(1.6);
54  TMVAStyle->cd();
55 
56 
57  Float_t xmin = TMath::Min( treeS->GetMinimum( var0 ), treeB->GetMinimum( var0 ) );
58  Float_t xmax = TMath::Max( treeS->GetMaximum( var0 ), treeB->GetMaximum( var0 ) );
59  Float_t ymin = TMath::Min( treeS->GetMinimum( var1 ), treeB->GetMinimum( var1 ) );
60  Float_t ymax = TMath::Max( treeS->GetMaximum( var1 ), treeB->GetMaximum( var1 ) );
61 
62  Int_t nbin = 500;
63  TH2F* frameS = new TH2F( "DataS", "DataS", nbin, xmin, xmax, nbin, ymin, ymax );
64  TH2F* frameB = new TH2F( "DataB", "DataB", nbin, xmin, xmax, nbin, ymin, ymax );
65 
66  // project trees
67  treeS->Draw( Form("%s:%s>>DataS",var1.Data(),var0.Data()), "", "0" );
68  treeB->Draw( Form("%s:%s>>DataB",var1.Data(),var0.Data()
69 ), "", "0" );
70 
71  // set style
72  frameS->SetMarkerSize( 0.1 );
73  frameS->SetMarkerColor( 4 );
74 
75  frameB->SetMarkerSize( 0.1 );
76  frameB->SetMarkerColor( 2 );
77 
78  // legend
79  frameS->SetTitle( var1+" versus "+var0+" for signal and background" );
80  frameS->GetXaxis()->SetTitle( var0 );
81  frameS->GetYaxis()->SetTitle( var1 );
82 
83  frameS->SetLabelSize( 0.04, "X" );
84  frameS->SetLabelSize( 0.04, "Y" );
85  frameS->SetTitleSize( 0.05, "X" );
86  frameS->SetTitleSize( 0.05, "Y" );
87 
88  // and plot
89  frameS->Draw();
90  frameB->Draw( "same" );
91 
92  // Draw legend
93  TLegend *legend = new TLegend( 1 - c->GetRightMargin() - 0.32, 1 - c->GetTopMargin() - 0.12,
94  1 - c->GetRightMargin(), 1 - c->GetTopMargin() );
95  legend->SetFillStyle( 1 );
96  legend->AddEntry(frameS,"Signal","p");
97  legend->AddEntry(frameB,"Background","p");
98  legend->Draw("same");
99  legend->SetBorderSize(1);
100  legend->SetMargin( 0.3 );
101 
102 }
103 
104 TMatrixD* produceSqrtMat( const TMatrixD& covMat )
105 {
106  Double_t sum = 0;
107  Int_t size = covMat.GetNrows();;
108  TMatrixD* sqrtMat = new TMatrixD( size, size );
109 
110  for (Int_t i=0; i< size; i++) {
111 
112  sum = 0;
113  for (Int_t j=0;j< i; j++) sum += (*sqrtMat)(i,j) * (*sqrtMat)(i,j);
114 
115  (*sqrtMat)(i,i) = TMath::Sqrt(TMath::Abs(covMat(i,i) - sum));
116 
117  for (Int_t k=i+1 ;k<size; k++) {
118 
119  sum = 0;
120  for (Int_t l=0; l<i; l++) sum += (*sqrtMat)(k,l) * (*sqrtMat)(i,l);
121 
122  (*sqrtMat)(k,i) = (covMat(k,i) - sum) / (*sqrtMat)(i,i);
123 
124  }
125  }
126  return sqrtMat;
127 }
128 
129 void getGaussRnd( TArrayD& v, const TMatrixD& sqrtMat, TRandom& R )
130 {
131  // generate "size" correlated Gaussian random numbers
132 
133  // sanity check
134  const Int_t size = sqrtMat.GetNrows();
135  if (size != v.GetSize())
136  cout << "<getGaussRnd> too short input vector: " << size << " " << v.GetSize() << endl;
137 
138  Double_t* tmpVec = new Double_t[size];
139 
140  for (Int_t i=0; i<size; i++) {
141  Double_t x, y, z;
142  y = R.Rndm();
143  z = R.Rndm();
144  x = 2*TMath::Pi()*z;
145  tmpVec[i] = TMath::Sin(x) * TMath::Sqrt(-2.0*TMath::Log(y));
146  }
147 
148  for (Int_t i=0; i<size; i++) {
149  v[i] = 0;
150  for (Int_t j=0; j<=i; j++) v[i] += sqrtMat(i,j) * tmpVec[j];
151  }
152 
153  delete[] tmpVec;
154 }
155 
156 // create the data
157 void create_lin_Nvar_withFriend(Int_t N = 2000)
158 {
159  const Int_t nvar = 4;
160  const Int_t nvar2 = 1;
161  Float_t xvar[nvar];
162 
163  // output flie
164  TFile* dataFile = TFile::Open( "data.root", "RECREATE" );
165 
166  // create signal and background trees
167  TTree* treeS = new TTree( "TreeS", "TreeS", 1 );
168  TTree* treeB = new TTree( "TreeB", "TreeB", 1 );
169  for (Int_t ivar=0; ivar<nvar-nvar2; ivar++) {
170  cout << "Creating branch var" << ivar+1 << " in signal tree" << endl;
171  treeS->Branch( TString(Form( "var%i", ivar+1 )).Data(), &xvar[ivar], TString(Form( "var%i/F", ivar+1 )).Data() );
172  treeB->Branch( TString(Form( "var%i", ivar+1 )).Data(), &xvar[ivar], TString(Form( "var%i/F", ivar+1 )).Data() );
173  }
174  TTree* treeSF = new TTree( "TreeSF", "TreeS", 1 );
175  TTree* treeBF = new TTree( "TreeBF", "TreeB", 1 );
176  for (Int_t ivar=nvar-nvar2; ivar<nvar; ivar++) {
177  treeSF->Branch( TString(Form( "var%i", ivar+1 )).Data(), &xvar[ivar], TString(Form( "var%i/F", ivar+1 )).Data() );
178  treeBF->Branch( TString(Form( "var%i", ivar+1 )).Data(), &xvar[ivar], TString(Form( "var%i/F", ivar+1 )).Data() );
179  }
180 
181 
182  TRandom R( 100 );
183  Float_t xS[nvar] = { 0.2, 0.3, 0.5, 0.9 };
184  Float_t xB[nvar] = { -0.2, -0.3, -0.5, -0.6 };
185  Float_t dx[nvar] = { 1.0, 1.0, 1.0, 1.0 };
186  TArrayD* v = new TArrayD( nvar );
187  Float_t rho[20];
188  rho[1*2] = 0.4;
189  rho[1*3] = 0.6;
190  rho[1*4] = 0.9;
191  rho[2*3] = 0.7;
192  rho[2*4] = 0.8;
193  rho[3*4] = 0.93;
194 
195  // create covariance matrix
196  TMatrixD* covMatS = new TMatrixD( nvar, nvar );
197  TMatrixD* covMatB = new TMatrixD( nvar, nvar );
198  for (Int_t ivar=0; ivar<nvar; ivar++) {
199  (*covMatS)(ivar,ivar) = dx[ivar]*dx[ivar];
200  (*covMatB)(ivar,ivar) = dx[ivar]*dx[ivar];
201  for (Int_t jvar=ivar+1; jvar<nvar; jvar++) {
202  (*covMatS)(ivar,jvar) = rho[(ivar+1)*(jvar+1)]*dx[ivar]*dx[jvar];
203  (*covMatS)(jvar,ivar) = (*covMatS)(ivar,jvar);
204 
205  (*covMatB)(ivar,jvar) = rho[(ivar+1)*(jvar+1)]*dx[ivar]*dx[jvar];
206  (*covMatB)(jvar,ivar) = (*covMatB)(ivar,jvar);
207  }
208  }
209 
210  cout << "signal covariance matrix: " << endl;
211  covMatS->Print();
212  cout << "background covariance matrix: " << endl;
213  covMatB->Print();
214 
215  // produce the square-root matrix
216  TMatrixD* sqrtMatS = produceSqrtMat( *covMatS );
217  TMatrixD* sqrtMatB = produceSqrtMat( *covMatB );
218 
219  // loop over species
220  for (Int_t itype=0; itype<2; itype++) {
221 
222  Float_t* x;
223  TMatrixD* m;
224  if (itype == 0) { x = xS; m = sqrtMatS; cout << "- produce signal" << endl; }
225  else { x = xB; m = sqrtMatB; cout << "- produce background" << endl; }
226 
227  // event loop
228  TTree* tree = (itype==0) ? treeS : treeB;
229  TTree* treeF = (itype==0) ? treeSF : treeBF;
230  for (Int_t i=0; i<N; i++) {
231 
232  if (i%1000 == 0) cout << "... event: " << i << " (" << N << ")" << endl;
233  getGaussRnd( *v, *m, R );
234 
235  for (Int_t ivar=0; ivar<nvar; ivar++) xvar[ivar] = (*v)[ivar] + x[ivar];
236 
237  tree->Fill();
238  treeF->Fill();
239  }
240  }
241 
242 // treeS->AddFriend(treeSF);
243 // treeB->AddFriend(treeBF);
244 
245  // write trees
246  treeS->Write();
247  treeB->Write();
248  treeSF->Write();
249  treeBF->Write();
250 
251  treeS->Show(0);
252  treeB->Show(1);
253 
254  dataFile->Close();
255  cout << "created data file: " << dataFile->GetName() << endl;
256 
257 
258 }
259 
260 
261 // create the tree
262 TTree* makeTree_lin_Nvar( TString treeName, TString treeTitle, Float_t* x, Float_t* dx, const Int_t nvar, Int_t N )
263 {
264  Float_t xvar[nvar];
265 
266  // create tree
267  TTree* tree = new TTree(treeName, treeTitle, 1);
268 
269  for (Int_t ivar=0; ivar<nvar; ivar++) {
270  tree->Branch( TString(Form( "var%i", ivar+1 )).Data(), &xvar[ivar], TString(Form( "var%i/F", ivar+1 )).Data() );
271  }
272 
273  TRandom R( 100 );
274  TArrayD* v = new TArrayD( nvar );
275  Float_t rho[20];
276  rho[1*2] = 0.4;
277  rho[1*3] = 0.6;
278  rho[1*4] = 0.9;
279  rho[2*3] = 0.7;
280  rho[2*4] = 0.8;
281  rho[3*4] = 0.93;
282 
283  // create covariance matrix
284  TMatrixD* covMat = new TMatrixD( nvar, nvar );
285  for (Int_t ivar=0; ivar<nvar; ivar++) {
286  (*covMat)(ivar,ivar) = dx[ivar]*dx[ivar];
287  for (Int_t jvar=ivar+1; jvar<nvar; jvar++) {
288  (*covMat)(ivar,jvar) = rho[(ivar+1)*(jvar+1)]*dx[ivar]*dx[jvar];
289  (*covMat)(jvar,ivar) = (*covMat)(ivar,jvar);
290  }
291  }
292  //cout << "covariance matrix: " << endl;
293  //covMat->Print();
294 
295  // produce the square-root matrix
296  TMatrixD* sqrtMat = produceSqrtMat( *covMat );
297 
298  // event loop
299  for (Int_t i=0; i<N; i++) {
300 
301  if (i%1000 == 0) cout << "... event: " << i << " (" << N << ")" << endl;
302  getGaussRnd( *v, *sqrtMat, R );
303 
304  for (Int_t ivar=0; ivar<nvar; ivar++) xvar[ivar] = (*v)[ivar] + x[ivar];
305 
306  tree->Fill();
307  }
308 
309  // write trees
310 // tree->Write();
311 
312  tree->Show(0);
313 
314  cout << "created tree: " << tree->GetName() << endl;
315  return tree;
316 }
317 
318 
319 // create the data
320 TTree* makeTree_circ(TString treeName, TString treeTitle, Int_t nvar = 2, Int_t N = 6000, Float_t radius = 1.0, Bool_t distort = false)
321 {
322  Int_t Nn = 0;
323  Float_t xvar[nvar]; //variable array size does not work in interactive mode
324 
325  // create signal and background trees
326  TTree* tree = new TTree( treeName, treeTitle, 1 );
327  for (Int_t ivar=0; ivar<nvar; ivar++) {
328  tree->Branch( TString(Form( "var%i", ivar+1 )).Data(), &xvar[ivar], TString(Form( "var%i/F", ivar+1 )).Data() );
329  }
330 
331  TRandom R( 100 );
332  //Float_t phimin = -30, phimax = 130;
333  Float_t phimin = -70, phimax = 130;
334  Float_t phisig = 5;
335  Float_t rsig = 0.1;
336  Float_t fnmin = -(radius+4.0*rsig);
337  Float_t fnmax = +(radius+4.0*rsig);
338  Float_t dfn = fnmax-fnmin;
339 
340  // event loop
341  for (Int_t i=0; i<N; i++) {
342  Double_t r1=R.Rndm(),r2=R.Rndm(), r3;
343  r3= r1>r2? r1 :r2;
344  Float_t phi;
345  if (distort) phi = r3*(phimax - phimin) + phimin;
346  else phi = R.Rndm()*(phimax - phimin) + phimin;
347  phi += R.Gaus()*phisig;
348 
349  Float_t r = radius;
350  r += R.Gaus()*rsig;
351 
352  xvar[0] = r*cos(TMath::DegToRad()*phi);
353  xvar[1] = r*sin(TMath::DegToRad()*phi);
354 
355  for( Int_t j = 2; j<nvar; ++j )
356  xvar[j] = dfn*R.Rndm()+fnmin;
357 
358  tree->Fill();
359  }
360 
361  for (Int_t i=0; i<Nn; i++) {
362 
363  xvar[0] = dfn*R.Rndm()+fnmin;
364  xvar[1] = dfn*R.Rndm()+fnmin;
365 
366  for( Int_t j = 2; j<nvar; ++j )
367  xvar[j] = dfn*R.Rndm()+fnmin;
368 
369 
370  tree->Fill();
371  }
372 
373  tree->Show(0);
374  // write trees
375  cout << "created tree: " << tree->GetName() << endl;
376  return tree;
377 }
378 
379 
380 
381 // create the data
382 void create_lin_Nvar_2(Int_t N = 50000)
383 {
384  const int nvar = 4;
385 
386  // output flie
387  TFile* dataFile = TFile::Open( "data.root", "RECREATE" );
388 
389 
390  Float_t xS[nvar] = { 0.2, 0.3, 0.5, 0.9 };
391  Float_t xB[nvar] = { -0.2, -0.3, -0.5, -0.6 };
392  Float_t dx[nvar] = { 1.0, 1.0, 1.0, 1.0 };
393 
394  // create signal and background trees
395  TTree* treeS = makeTree_lin_Nvar( "TreeS", "Signal tree", xS, dx, nvar, N );
396  TTree* treeB = makeTree_lin_Nvar( "TreeB", "Background tree", xB, dx, nvar, N );
397 
398  treeS->Write();
399  treeB->Write();
400 
401  treeS->Show(0);
402  treeB->Show(0);
403 
404  dataFile->Close();
405  cout << "created data file: " << dataFile->GetName() << endl;
406 }
407 
408 
409 
410 
411 // create the data
412 void create_lin_Nvar(Int_t N = 50000)
413 {
414  const Int_t nvar = 4;
415  Float_t xvar[nvar];
416 
417  // output flie
418  TFile* dataFile = TFile::Open( "data.root", "RECREATE" );
419 
420  // create signal and background trees
421  TTree* treeS = new TTree( "TreeS", "TreeS", 1 );
422  TTree* treeB = new TTree( "TreeB", "TreeB", 1 );
423  for (Int_t ivar=0; ivar<nvar; ivar++) {
424  treeS->Branch( TString(Form( "var%i", ivar+1 )).Data(), &xvar[ivar], TString(Form( "var%i/F", ivar+1 )).Data() );
425  treeB->Branch( TString(Form( "var%i", ivar+1 )).Data(), &xvar[ivar], TString(Form( "var%i/F", ivar+1 )).Data() );
426  }
427 
428  TRandom R( 100 );
429  Float_t xS[nvar] = { 0.2, 0.3, 0.5, 0.9 };
430  Float_t xB[nvar] = { -0.2, -0.3, -0.5, -0.6 };
431  Float_t dx[nvar] = { 1.0, 1.0, 1.0, 1.0 };
432  TArrayD* v = new TArrayD( nvar );
433  Float_t rho[20];
434  rho[1*2] = 0.4;
435  rho[1*3] = 0.6;
436  rho[1*4] = 0.9;
437  rho[2*3] = 0.7;
438  rho[2*4] = 0.8;
439  rho[3*4] = 0.93;
440 
441  // create covariance matrix
442  TMatrixD* covMatS = new TMatrixD( nvar, nvar );
443  TMatrixD* covMatB = new TMatrixD( nvar, nvar );
444  for (Int_t ivar=0; ivar<nvar; ivar++) {
445  (*covMatS)(ivar,ivar) = dx[ivar]*dx[ivar];
446  (*covMatB)(ivar,ivar) = dx[ivar]*dx[ivar];
447  for (Int_t jvar=ivar+1; jvar<nvar; jvar++) {
448  (*covMatS)(ivar,jvar) = rho[(ivar+1)*(jvar+1)]*dx[ivar]*dx[jvar];
449  (*covMatS)(jvar,ivar) = (*covMatS)(ivar,jvar);
450 
451  (*covMatB)(ivar,jvar) = rho[(ivar+1)*(jvar+1)]*dx[ivar]*dx[jvar];
452  (*covMatB)(jvar,ivar) = (*covMatB)(ivar,jvar);
453  }
454  }
455  cout << "signal covariance matrix: " << endl;
456  covMatS->Print();
457  cout << "background covariance matrix: " << endl;
458  covMatB->Print();
459 
460  // produce the square-root matrix
461  TMatrixD* sqrtMatS = produceSqrtMat( *covMatS );
462  TMatrixD* sqrtMatB = produceSqrtMat( *covMatB );
463 
464  // loop over species
465  for (Int_t itype=0; itype<2; itype++) {
466 
467  Float_t* x;
468  TMatrixD* m;
469  if (itype == 0) { x = xS; m = sqrtMatS; cout << "- produce signal" << endl; }
470  else { x = xB; m = sqrtMatB; cout << "- produce background" << endl; }
471 
472  // event loop
473  TTree* tree = (itype==0) ? treeS : treeB;
474  for (Int_t i=0; i<N; i++) {
475 
476  if (i%1000 == 0) cout << "... event: " << i << " (" << N << ")" << endl;
477  getGaussRnd( *v, *m, R );
478 
479  for (Int_t ivar=0; ivar<nvar; ivar++) xvar[ivar] = (*v)[ivar] + x[ivar];
480 
481  tree->Fill();
482  }
483  }
484 
485  // write trees
486  treeS->Write();
487  treeB->Write();
488 
489  treeS->Show(0);
490  treeB->Show(1);
491 
492  dataFile->Close();
493  cout << "created data file: " << dataFile->GetName() << endl;
494 }
495 
496 // create the category data
497 // type = 1 (offset) or 2 (variable = -99)
498 void create_lin_Nvar_categories(Int_t N = 10000, Int_t type = 2)
499 {
500  const Int_t nvar = 4;
501  Float_t xvar[nvar];
502  Float_t eta;
503 
504  // output flie
505  TFile* dataFile = TFile::Open( "data.root", "RECREATE" );
506 
507  // create signal and background trees
508  TTree* treeS = new TTree( "TreeS", "TreeS", 1 );
509  TTree* treeB = new TTree( "TreeB", "TreeB", 1 );
510  for (Int_t ivar=0; ivar<nvar; ivar++) {
511  treeS->Branch( TString(Form( "var%i", ivar+1 )).Data(), &xvar[ivar], TString(Form( "var%i/F", ivar+1 )).Data() );
512  treeB->Branch( TString(Form( "var%i", ivar+1 )).Data(), &xvar[ivar], TString(Form( "var%i/F", ivar+1 )).Data() );
513  }
514 
515  // add category variable
516  treeS->Branch( "eta", &eta, "eta/F" );
517  treeB->Branch( "eta", &eta, "eta/F" );
518 
519  TRandom R( 100 );
520  Float_t xS[nvar] = { 0.2, 0.3, 0.5, 0.9 };
521  Float_t xB[nvar] = { -0.2, -0.3, -0.5, -0.6 };
522  Float_t dx[nvar] = { 1.0, 1.0, 1.0, 1.0 };
523  TArrayD* v = new TArrayD( nvar );
524  Float_t rho[20];
525  rho[1*2] = 0.0;
526  rho[1*3] = 0.0;
527  rho[1*4] = 0.0;
528  rho[2*3] = 0.0;
529  rho[2*4] = 0.0;
530  rho[3*4] = 0.0;
531  if (type != 1) {
532  rho[1*2] = 0.6;
533  rho[1*3] = 0.7;
534  rho[1*4] = 0.9;
535  rho[2*3] = 0.8;
536  rho[2*4] = 0.9;
537  rho[3*4] = 0.93;
538  }
539 
540  // create covariance matrix
541  TMatrixD* covMatS = new TMatrixD( nvar, nvar );
542  TMatrixD* covMatB = new TMatrixD( nvar, nvar );
543  for (Int_t ivar=0; ivar<nvar; ivar++) {
544  (*covMatS)(ivar,ivar) = dx[ivar]*dx[ivar];
545  (*covMatB)(ivar,ivar) = dx[ivar]*dx[ivar];
546  for (Int_t jvar=ivar+1; jvar<nvar; jvar++) {
547  (*covMatS)(ivar,jvar) = rho[(ivar+1)*(jvar+1)]*dx[ivar]*dx[jvar];
548  (*covMatS)(jvar,ivar) = (*covMatS)(ivar,jvar);
549 
550  (*covMatB)(ivar,jvar) = rho[(ivar+1)*(jvar+1)]*dx[ivar]*dx[jvar];
551  (*covMatB)(jvar,ivar) = (*covMatB)(ivar,jvar);
552  }
553  }
554  cout << "signal covariance matrix: " << endl;
555  covMatS->Print();
556  cout << "background covariance matrix: " << endl;
557  covMatB->Print();
558 
559  // produce the square-root matrix
560  TMatrixD* sqrtMatS = produceSqrtMat( *covMatS );
561  TMatrixD* sqrtMatB = produceSqrtMat( *covMatB );
562 
563  // loop over species
564  for (Int_t itype=0; itype<2; itype++) {
565 
566  Float_t* x;
567  TMatrixD* m;
568  if (itype == 0) { x = xS; m = sqrtMatS; cout << "- produce signal" << endl; }
569  else { x = xB; m = sqrtMatB; cout << "- produce background" << endl; }
570 
571  // event loop
572  TTree* tree = (itype==0) ? treeS : treeB;
573  for (Int_t i=0; i<N; i++) {
574 
575  if (i%1000 == 0) cout << "... event: " << i << " (" << N << ")" << endl;
576  getGaussRnd( *v, *m, R );
577 
578  eta = 2.5*2*(R.Rndm() - 0.5);
579  Float_t offset = 0;
580  if (type == 1) offset = TMath::Abs(eta) > 1.3 ? 0.8 : -0.8;
581  for (Int_t ivar=0; ivar<nvar; ivar++) xvar[ivar] = (*v)[ivar] + x[ivar] + offset;
582  if (type != 1 && TMath::Abs(eta) > 1.3) xvar[nvar-1] = -5;
583 
584  tree->Fill();
585  }
586  }
587 
588  // write trees
589  treeS->Write();
590  treeB->Write();
591 
592  treeS->Show(0);
593  treeB->Show(1);
594 
595  dataFile->Close();
596  cout << "created data file: " << dataFile->GetName() << endl;
597 }
598 
599 
600 // create the data
601 void create_lin_Nvar_weighted(Int_t N = 10000, int WeightedSignal=0, int WeightedBkg=1, Float_t BackgroundContamination=0, Int_t seed=100)
602 {
603  const Int_t nvar = 4;
604  Float_t xvar[nvar];
605  Float_t weight;
606 
607 
608  cout << endl << endl << endl;
609  cout << "please use .L createData.C++ if you want to run this MC geneation" <<endl;
610  cout << "otherwise you will wait for ages!!! " << endl;
611  cout << endl << endl << endl;
612 
613 
614  // output flie
615  TString fileName;
616  if (BackgroundContamination) fileName = Form("linCorGauss%d_weighted+background.root",seed);
617  else fileName = Form("linCorGauss%d_weighted.root",seed);
618 
619  TFile* dataFile = TFile::Open( fileName.Data(), "RECREATE" );
620 
621  // create signal and background trees
622  TTree* treeS = new TTree( "TreeS", "TreeS", 1 );
623  TTree* treeB = new TTree( "TreeB", "TreeB", 1 );
624  for (Int_t ivar=0; ivar<nvar; ivar++) {
625  treeS->Branch( TString(Form( "var%i", ivar+1 )).Data(), &xvar[ivar], TString(Form( "var%i/F", ivar+1 )).Data() );
626  treeB->Branch( TString(Form( "var%i", ivar+1 )).Data(), &xvar[ivar], TString(Form( "var%i/F", ivar+1 )).Data() );
627  }
628  if (WeightedSignal||BackgroundContamination>0||1) treeS->Branch( "weight", &weight,"weight/F" );
629  if (WeightedBkg) treeB->Branch( "weight", &weight,"weight/F" );
630 
631  TRandom R( seed );
632  Float_t xS[nvar] = { 0.2, 0.3, 0.4, 0.8 };
633  Float_t xB[nvar] = { -0.2, -0.3, -0.4, -0.5 };
634  Float_t dx[nvar] = { 1.0, 1.0, 1.0, 1.0 };
635  TArrayD* v = new TArrayD( nvar );
636  Float_t rho[20];
637  rho[1*2] = 0.4;
638  rho[1*3] = 0.6;
639  rho[1*4] = 0.9;
640  rho[2*3] = 0.7;
641  rho[2*4] = 0.8;
642  rho[3*4] = 0.93;
643 
644  // create covariance matrix
645  TMatrixD* covMatS = new TMatrixD( nvar, nvar );
646  TMatrixD* covMatB = new TMatrixD( nvar, nvar );
647  for (Int_t ivar=0; ivar<nvar; ivar++) {
648  (*covMatS)(ivar,ivar) = dx[ivar]*dx[ivar];
649  (*covMatB)(ivar,ivar) = dx[ivar]*dx[ivar];
650  for (Int_t jvar=ivar+1; jvar<nvar; jvar++) {
651  (*covMatS)(ivar,jvar) = rho[(ivar+1)*(jvar+1)]*dx[ivar]*dx[jvar];
652  (*covMatS)(jvar,ivar) = (*covMatS)(ivar,jvar);
653 
654  (*covMatB)(ivar,jvar) = rho[(ivar+1)*(jvar+1)]*dx[ivar]*dx[jvar];
655  (*covMatB)(jvar,ivar) = (*covMatB)(ivar,jvar);
656  }
657  }
658  cout << "signal covariance matrix: " << endl;
659  covMatS->Print();
660  cout << "background covariance matrix: " << endl;
661  covMatB->Print();
662 
663  // produce the square-root matrix
664  TMatrixD* sqrtMatS = produceSqrtMat( *covMatS );
665  TMatrixD* sqrtMatB = produceSqrtMat( *covMatB );
666 
667  // loop over species
668  for (Int_t itype=0; itype<2; itype++) {
669 
670  Float_t* x;
671  TMatrixD* m;
672  if (itype == 0) { x = xS; m = sqrtMatS; cout << "- produce signal" << endl; }
673  else { x = xB; m = sqrtMatB; cout << "- produce background" << endl; }
674 
675  // event loop
676  TTree* tree = (itype==0) ? treeS : treeB;
677  Int_t i=0;
678  do {
679  getGaussRnd( *v, *m, R );
680 
681  for (Int_t ivar=0; ivar<nvar; ivar++) xvar[ivar] = (*v)[ivar] + x[ivar];
682  // for (Int_t ivar=0; ivar<nvar; ivar++) xvar[ivar] = R.Uniform()*10.-5.;
683 
684  // weight = 0.5 / (TMath::Gaus( (xvar[nvar-1]-x[nvar-1]), 0, 1.1) );
685  // weight = TMath::Gaus(0.675,0,1) / (TMath::Gaus( (xvar[nvar-1]-x[nvar-1]), 0, 1.) );
686  weight = 0.8 / (TMath::Gaus( ((*v)[nvar-1]), 0, 1.09) );
687  Double_t tmp=R.Uniform()/0.00034;
688  if (itype==0 && !WeightedSignal) {
689  weight = 1;
690  tree->Fill();
691  i++;
692  } else if (itype==1 && !WeightedBkg) {
693  weight = 1;
694  tree->Fill();
695  i++;
696  }
697  else {
698  if (tmp < weight){
699  weight = 1./weight;
700  tree->Fill();
701  if (i%10 == 0) cout << "... event: " << i << " (" << N << ")" << endl;
702  i++;
703  }
704  }
705  } while (i<N);
706  }
707 
708 
709  if (BackgroundContamination > 0){ // add "background contamination" in the Signal (which later is again "subtracted" with
710  // using (statistically indepentent) background events with negative weight)
711  Float_t* x=xB;
712  TMatrixD* m = sqrtMatB;
713  TTree* tree = treeS;
714  for (Int_t i=0; i<N*BackgroundContamination*2; i++) {
715  if (i%1000 == 0) cout << "... event: " << i << " (" << N << ")" << endl;
716  getGaussRnd( *v, *m, R );
717  for (Int_t ivar=0; ivar<nvar; ivar++) xvar[ivar] = (*v)[ivar] + x[ivar];
718 
719  // add weights
720  if (i%2) weight = 1;
721  else weight = -1;
722 
723  tree->Fill();
724  }
725  }
726 
727 
728 
729  // write trees
730  treeS->Write();
731  treeB->Write();
732 
733  treeS->Show(0);
734  treeB->Show(1);
735 
736  TH1F *h[4];
737  TH1F *hw[4];
738  for (Int_t i=0;i<4;i++){
739  char buffer[5];
740  sprintf(buffer,"h%d",i);
741  h[i]= new TH1F(buffer,"",100,-5,5);
742  sprintf(buffer,"hw%d",i);
743  hw[i] = new TH1F(buffer,"",100,-5,5);
744  hw[i]->SetLineColor(3);
745  }
746 
747  for (int ie=0;ie<treeS->GetEntries();ie++){
748  treeS->GetEntry(ie);
749  for (Int_t i=0;i<4;i++){
750  h[i]->Fill(xvar[i]);
751  hw[i]->Fill(xvar[i],weight);
752  }
753  }
754 
755  TCanvas *c = new TCanvas("c","",800,800);
756  c->Divide(2,2);
757 
758  for (Int_t i=0;i<4;i++){
759  c->cd(i+1);
760  h[i]->Draw();
761  hw[i]->Draw("same");
762  }
763 
764 
765  // dataFile->Close();
766  cout << "created data file: " << dataFile->GetName() << endl;
767 }
768 
769 
770 
771 // create the data
772 void create_lin_Nvar_Arr(Int_t N = 1000)
773 {
774  const Int_t nvar = 4;
775  std::vector<float>* xvar[nvar];
776 
777  // output flie
778  TFile* dataFile = TFile::Open( "data.root", "RECREATE" );
779 
780  // create signal and background trees
781  TTree* treeS = new TTree( "TreeS", "TreeS", 1 );
782  TTree* treeB = new TTree( "TreeB", "TreeB", 1 );
783  for (Int_t ivar=0; ivar<nvar; ivar++) {
784  xvar[ivar] = new std::vector<float>();
785  treeS->Branch( TString(Form( "var%i", ivar+1 )).Data(), "vector<float>", &xvar[ivar], 64000, 1 );
786  treeB->Branch( TString(Form( "var%i", ivar+1 )).Data(), "vector<float>", &xvar[ivar], 64000, 1 );
787  }
788 
789  TRandom R( 100 );
790  Float_t xS[nvar] = { 0.2, 0.3, 0.5, 0.9 };
791  Float_t xB[nvar] = { -0.2, -0.3, -0.5, -0.6 };
792  Float_t dx[nvar] = { 1.0, 1.0, 1.0, 1.0 };
793  TArrayD* v = new TArrayD( nvar );
794  Float_t rho[20];
795  rho[1*2] = 0.4;
796  rho[1*3] = 0.6;
797  rho[1*4] = 0.9;
798  rho[2*3] = 0.7;
799  rho[2*4] = 0.8;
800  rho[3*4] = 0.93;
801 
802  // create covariance matrix
803  TMatrixD* covMatS = new TMatrixD( nvar, nvar );
804  TMatrixD* covMatB = new TMatrixD( nvar, nvar );
805  for (Int_t ivar=0; ivar<nvar; ivar++) {
806  (*covMatS)(ivar,ivar) = dx[ivar]*dx[ivar];
807  (*covMatB)(ivar,ivar) = dx[ivar]*dx[ivar];
808  for (Int_t jvar=ivar+1; jvar<nvar; jvar++) {
809  (*covMatS)(ivar,jvar) = rho[(ivar+1)*(jvar+1)]*dx[ivar]*dx[jvar];
810  (*covMatS)(jvar,ivar) = (*covMatS)(ivar,jvar);
811 
812  (*covMatB)(ivar,jvar) = rho[(ivar+1)*(jvar+1)]*dx[ivar]*dx[jvar];
813  (*covMatB)(jvar,ivar) = (*covMatB)(ivar,jvar);
814  }
815  }
816  cout << "signal covariance matrix: " << endl;
817  covMatS->Print();
818  cout << "background covariance matrix: " << endl;
819  covMatB->Print();
820 
821  // produce the square-root matrix
822  TMatrixD* sqrtMatS = produceSqrtMat( *covMatS );
823  TMatrixD* sqrtMatB = produceSqrtMat( *covMatB );
824 
825  // loop over species
826  for (Int_t itype=0; itype<2; itype++) {
827 
828  Float_t* x;
829  TMatrixD* m;
830  if (itype == 0) { x = xS; m = sqrtMatS; cout << "- produce signal" << endl; }
831  else { x = xB; m = sqrtMatB; cout << "- produce background" << endl; }
832 
833  // event loop
834  TTree* tree = (itype==0) ? treeS : treeB;
835  for (Int_t i=0; i<N; i++) {
836 
837  if (i%100 == 0) cout << "... event: " << i << " (" << N << ")" << endl;
838 
839  Int_t aSize = (Int_t)(gRandom->Rndm()*10); // size of array varies between events
840  for (Int_t ivar=0; ivar<nvar; ivar++) {
841  xvar[ivar]->clear();
842  xvar[ivar]->reserve(aSize);
843  }
844  for(Int_t iA = 0; iA<aSize; iA++) {
845  //for (Int_t ivar=0; ivar<nvar; ivar++) {
846  getGaussRnd( *v, *m, R );
847  for (Int_t ivar=0; ivar<nvar; ivar++) xvar[ivar]->push_back((*v)[ivar] + x[ivar]);
848  //}
849  }
850  tree->Fill();
851  }
852  }
853 
854  // write trees
855  treeS->Write();
856  treeB->Write();
857 
858  treeS->Show(0);
859  treeB->Show(1);
860 
861  dataFile->Close();
862  cout << "created data file: " << dataFile->GetName() << endl;
863 
864  //plot();
865 }
866 
867 
868 
869 // create the data
870 void create_lin_Nvar_double()
871 {
872  Int_t N = 10000;
873  const Int_t nvar = 4;
874  Double_t xvar[nvar];
875  Double_t xvarD[nvar];
876  Float_t xvarF[nvar];
877 
878  // output flie
879  TFile* dataFile = TFile::Open( "data.root", "RECREATE" );
880 
881  // create signal and background trees
882  TTree* treeS = new TTree( "TreeS", "TreeS", 1 );
883  TTree* treeB = new TTree( "TreeB", "TreeB", 1 );
884  for (Int_t ivar=0; ivar<nvar; ivar++) {
885  if (ivar<2) {
886  treeS->Branch( TString(Form( "var%i", ivar+1 )).Data(), &xvarD[ivar], TString(Form( "var%i/D", ivar+1 )).Data() );
887  treeB->Branch( TString(Form( "var%i", ivar+1 )).Data(), &xvarD[ivar], TString(Form( "var%i/D", ivar+1 )).Data() );
888  }
889  else {
890  treeS->Branch( TString(Form( "var%i", ivar+1 )).Data(), &xvarF[ivar], TString(Form( "var%i/F", ivar+1 )).Data() );
891  treeB->Branch( TString(Form( "var%i", ivar+1 )).Data(), &xvarF[ivar], TString(Form( "var%i/F", ivar+1 )).Data() );
892  }
893  }
894 
895  TRandom R( 100 );
896  Double_t xS[nvar] = { 0.2, 0.3, 0.5, 0.6 };
897  Double_t xB[nvar] = { -0.2, -0.3, -0.5, -0.6 };
898  Double_t dx[nvar] = { 1.0, 1.0, 1.0, 1.0 };
899  TArrayD* v = new TArrayD( nvar );
900  Double_t rho[20];
901  rho[1*2] = 0.4;
902  rho[1*3] = 0.6;
903  rho[1*4] = 0.9;
904  rho[2*3] = 0.7;
905  rho[2*4] = 0.8;
906  rho[3*4] = 0.93;
907 
908  // create covariance matrix
909  TMatrixD* covMatS = new TMatrixD( nvar, nvar );
910  TMatrixD* covMatB = new TMatrixD( nvar, nvar );
911  for (Int_t ivar=0; ivar<nvar; ivar++) {
912  (*covMatS)(ivar,ivar) = dx[ivar]*dx[ivar];
913  (*covMatB)(ivar,ivar) = dx[ivar]*dx[ivar];
914  for (Int_t jvar=ivar+1; jvar<nvar; jvar++) {
915  (*covMatS)(ivar,jvar) = rho[(ivar+1)*(jvar+1)]*dx[ivar]*dx[jvar];
916  (*covMatS)(jvar,ivar) = (*covMatS)(ivar,jvar);
917 
918  (*covMatB)(ivar,jvar) = rho[(ivar+1)*(jvar+1)]*dx[ivar]*dx[jvar];
919  (*covMatB)(jvar,ivar) = (*covMatB)(ivar,jvar);
920  }
921  }
922  cout << "signal covariance matrix: " << endl;
923  covMatS->Print();
924  cout << "background covariance matrix: " << endl;
925  covMatB->Print();
926 
927  // produce the square-root matrix
928  TMatrixD* sqrtMatS = produceSqrtMat( *covMatS );
929  TMatrixD* sqrtMatB = produceSqrtMat( *covMatB );
930 
931  // loop over species
932  for (Int_t itype=0; itype<2; itype++) {
933 
934  Double_t* x;
935  TMatrixD* m;
936  if (itype == 0) { x = xS; m = sqrtMatS; cout << "- produce signal" << endl; }
937  else { x = xB; m = sqrtMatB; cout << "- produce background" << endl; }
938 
939  // event loop
940  TTree* tree = (itype==0) ? treeS : treeB;
941  for (Int_t i=0; i<N; i++) {
942 
943  if (i%1000 == 0) cout << "... event: " << i << " (" << N << ")" << endl;
944  getGaussRnd( *v, *m, R );
945 
946  for (Int_t ivar=0; ivar<nvar; ivar++) xvar[ivar] = (*v)[ivar] + x[ivar];
947  for (Int_t ivar=0; ivar<nvar; ivar++) {
948  if (ivar<2) xvarD[ivar] = xvar[ivar];
949  else xvarF[ivar] = xvar[ivar];
950  }
951 
952  tree->Fill();
953  }
954  }
955 
956  // write trees
957  treeS->Write();
958  treeB->Write();
959 
960  treeS->Show(0);
961  treeB->Show(1);
962 
963  dataFile->Close();
964  cout << "created data file: " << dataFile->GetName() << endl;
965 
966  plot();
967 }
968 
969 // create the data
970 void create_lin_Nvar_discrete()
971 {
972  Int_t N = 10000;
973  const Int_t nvar = 4;
974  Float_t xvar[nvar];
975  Int_t xvarI[2];
976 
977  // output flie
978  TFile* dataFile = TFile::Open( "data.root", "RECREATE" );
979 
980  // create signal and background trees
981  TTree* treeS = new TTree( "TreeS", "TreeS", 1 );
982  TTree* treeB = new TTree( "TreeB", "TreeB", 1 );
983  for (Int_t ivar=0; ivar<nvar-2; ivar++) {
984  treeS->Branch( TString(Form( "var%i", ivar+1 )).Data(), &xvar[ivar], TString(Form( "var%i/F", ivar+1 )).Data() );
985  treeB->Branch( TString(Form( "var%i", ivar+1 )).Data(), &xvar[ivar], TString(Form( "var%i/F", ivar+1 )).Data() );
986  }
987  for (Int_t ivar=0; ivar<2; ivar++) {
988  treeS->Branch( TString(Form( "var%i", ivar+nvar-2+1 )).Data(), &xvarI[ivar], TString(Form( "var%i/I", ivar+nvar-2+1 )).Data() );
989  treeB->Branch( TString(Form( "var%i", ivar+nvar-2+1 )).Data(), &xvarI[ivar], TString(Form( "var%i/I", ivar+nvar-2+1 )).Data() );
990  }
991 
992  TRandom R( 100 );
993  Float_t xS[nvar] = { 0.2, 0.3, 1, 2 };
994  Float_t xB[nvar] = { -0.2, -0.3, 0, 0 };
995  Float_t dx[nvar] = { 1.0, 1.0, 1, 2 };
996  TArrayD* v = new TArrayD( nvar );
997  Float_t rho[20];
998  rho[1*2] = 0.4;
999  rho[1*3] = 0.6;
1000  rho[1*4] = 0.9;
1001  rho[2*3] = 0.7;
1002  rho[2*4] = 0.8;
1003  rho[3*4] = 0.93;
1004  // no correlations
1005  for (int i=0; i<20; i++) rho[i] = 0;
1006 
1007  // create covariance matrix
1008  TMatrixD* covMatS = new TMatrixD( nvar, nvar );
1009  TMatrixD* covMatB = new TMatrixD( nvar, nvar );
1010  for (Int_t ivar=0; ivar<nvar; ivar++) {
1011  (*covMatS)(ivar,ivar) = dx[ivar]*dx[ivar];
1012  (*covMatB)(ivar,ivar) = dx[ivar]*dx[ivar];
1013  for (Int_t jvar=ivar+1; jvar<nvar; jvar++) {
1014  (*covMatS)(ivar,jvar) = rho[(ivar+1)*(jvar+1)]*dx[ivar]*dx[jvar];
1015  (*covMatS)(jvar,ivar) = (*covMatS)(ivar,jvar);
1016 
1017  (*covMatB)(ivar,jvar) = rho[(ivar+1)*(jvar+1)]*dx[ivar]*dx[jvar];
1018  (*covMatB)(jvar,ivar) = (*covMatB)(ivar,jvar);
1019  }
1020  }
1021  cout << "signal covariance matrix: " << endl;
1022  covMatS->Print();
1023  cout << "background covariance matrix: " << endl;
1024  covMatB->Print();
1025 
1026  // produce the square-root matrix
1027  TMatrixD* sqrtMatS = produceSqrtMat( *covMatS );
1028  TMatrixD* sqrtMatB = produceSqrtMat( *covMatB );
1029 
1030  // loop over species
1031  for (Int_t itype=0; itype<2; itype++) {
1032 
1033  Float_t* x;
1034  TMatrixD* m;
1035  if (itype == 0) { x = xS; m = sqrtMatS; cout << "- produce signal" << endl; }
1036  else { x = xB; m = sqrtMatB; cout << "- produce background" << endl; }
1037 
1038  // event loop
1039  TTree* tree = (itype==0) ? treeS : treeB;
1040  for (Int_t i=0; i<N; i++) {
1041 
1042  if (i%1000 == 0) cout << "... event: " << i << " (" << N << ")" << endl;
1043  getGaussRnd( *v, *m, R );
1044 
1045  for (Int_t ivar=0; ivar<nvar; ivar++) xvar[ivar] = (*v)[ivar] + x[ivar];
1046 
1047  xvarI[0] = TMath::Nint(xvar[nvar-2]);
1048  xvarI[1] = TMath::Nint(xvar[nvar-1]);
1049 
1050  tree->Fill();
1051  }
1052  }
1053 
1054  // write trees
1055  treeS->Write();
1056  treeB->Write();
1057 
1058  treeS->Show(0);
1059  treeB->Show(1);
1060 
1061  dataFile->Close();
1062  cout << "created data file: " << dataFile->GetName() << endl;
1063 
1064  plot();
1065 }
1066 
1067 // create the data
1068 void create_ManyVars()
1069 {
1070  Int_t N = 20000;
1071  const Int_t nvar = 20;
1072  Float_t xvar[nvar];
1073 
1074  // output flie
1075  TFile* dataFile = TFile::Open( "data.root", "RECREATE" );
1076 
1077  // create signal and background trees
1078  TTree* treeS = new TTree( "TreeS", "TreeS", 1 );
1079  TTree* treeB = new TTree( "TreeB", "TreeB", 1 );
1080  for (Int_t ivar=0; ivar<nvar; ivar++) {
1081  treeS->Branch( TString(Form( "var%i", ivar )).Data(), &xvar[ivar], TString(Form( "var%i/F", ivar )).Data() );
1082  treeB->Branch( TString(Form( "var%i", ivar )).Data(), &xvar[ivar], TString(Form( "var%i/F", ivar )).Data() );
1083  }
1084 
1085  Float_t xS[nvar];
1086  Float_t xB[nvar];
1087  Float_t dx[nvar];
1088  for (Int_t ivar=0; ivar<nvar; ivar++) {
1089  xS[ivar] = 0 + ivar*0.05;
1090  xB[ivar] = 0 - ivar*0.05;
1091  dx[ivar] = 1;
1092  }
1093 
1094  xS[0] = 0.2;
1095  xB[0] = -0.2;
1096  dx[0] = 1.0;
1097  xS[1] = 0.3;
1098  xB[1] = -0.3;
1099  dx[1] = 1.0;
1100  xS[2] = 0.4;
1101  xB[2] = -0.4;
1102  dx[2] = 1.0 ;
1103  xS[3] = 0.8 ;
1104  xB[3] = -0.5 ;
1105  dx[3] = 1.0 ;
1106 // TArrayD* v = new TArrayD( nvar );
1107  Float_t rho[20];
1108  rho[1*2] = 0.4;
1109  rho[1*3] = 0.6;
1110  rho[1*4] = 0.9;
1111  rho[2*3] = 0.7;
1112  rho[2*4] = 0.8;
1113  rho[3*4] = 0.93;
1114 
1115  TRandom R( 100 );
1116 
1117  // loop over species
1118  for (Int_t itype=0; itype<2; itype++) {
1119 
1120  Float_t* x = (itype == 0) ? xS : xB;
1121 
1122  // event loop
1123  TTree* tree = (itype == 0) ? treeS : treeB;
1124  for (Int_t i=0; i<N; i++) {
1125 
1126  if (i%1000 == 0) cout << "... event: " << i << " (" << N << ")" << endl;
1127  for (Int_t ivar=0; ivar<nvar; ivar++) {
1128  if (ivar == 1500 && itype!=10) xvar[ivar] = 1;
1129  else xvar[ivar] = x[ivar] + R.Gaus()*dx[ivar];
1130  }
1131 
1132  tree->Fill();
1133  }
1134  }
1135 
1136  // write trees
1137  treeS->Write();
1138  treeB->Write();
1139 
1140  treeS->Show(0);
1141  treeB->Show(1);
1142 
1143  dataFile->Close();
1144  plot();
1145  cout << "created data file: " << dataFile->GetName() << endl;
1146 }
1147 
1148 // create the data
1149 void create_lin_NvarObsolete()
1150 {
1151  Int_t N = 20000;
1152  const Int_t nvar = 20;
1153  Float_t xvar[nvar];
1154 
1155  // output flie
1156  TFile* dataFile = TFile::Open( "data.root", "RECREATE" );
1157 
1158  // create signal and background trees
1159  TTree* treeS = new TTree( "TreeS", "TreeS", 1 );
1160  TTree* treeB = new TTree( "TreeB", "TreeB", 1 );
1161  for (Int_t ivar=0; ivar<nvar; ivar++) {
1162  treeS->Branch( TString(Form( "var%i", ivar )).Data(), &xvar[ivar], TString(Form( "var%i/F", ivar )).Data() );
1163  treeB->Branch( TString(Form( "var%i", ivar )).Data(), &xvar[ivar], TString(Form( "var%i/F", ivar )).Data() );
1164  }
1165 
1166  TRandom R( 100 );
1167  Float_t xS[nvar] = { 0.5, 0.5, 0.0, 0.0, 0.0, 0.0 };
1168  Float_t xB[nvar] = { -0.5, -0.5, -0.0, -0.0, -0.0, -0.0 };
1169  Float_t dx[nvar] = { 1.0, 1.0, 1.0, 1.0, 1.0, 1.0 };
1170  TArrayD* v = new TArrayD( nvar );
1171  Float_t rho[50];
1172  for (Int_t i=0; i<50; i++) rho[i] = 0;
1173  rho[1*2] = 0.3;
1174  rho[1*3] = 0.0;
1175  rho[1*4] = 0.0;
1176  rho[2*3] = 0.0;
1177  rho[2*4] = 0.0;
1178  rho[3*4] = 0.0;
1179 
1180  // create covariance matrix
1181  TMatrixD* covMatS = new TMatrixD( nvar, nvar );
1182  TMatrixD* covMatB = new TMatrixD( nvar, nvar );
1183  for (Int_t ivar=0; ivar<nvar; ivar++) {
1184  (*covMatS)(ivar,ivar) = dx[ivar]*dx[ivar];
1185  (*covMatB)(ivar,ivar) = dx[ivar]*dx[ivar];
1186  for (Int_t jvar=ivar+1; jvar<nvar; jvar++) {
1187  (*covMatS)(ivar,jvar) = rho[(ivar+1)*(jvar+1)]*dx[ivar]*dx[jvar];
1188  (*covMatS)(jvar,ivar) = (*covMatS)(ivar,jvar);
1189 
1190  (*covMatB)(ivar,jvar) = rho[(ivar+1)*(jvar+1)]*dx[ivar]*dx[jvar];
1191  (*covMatB)(jvar,ivar) = (*covMatB)(ivar,jvar);
1192  }
1193  }
1194  cout << "signal covariance matrix: " << endl;
1195  covMatS->Print();
1196  cout << "background covariance matrix: " << endl;
1197  covMatB->Print();
1198 
1199  // produce the square-root matrix
1200  TMatrixD* sqrtMatS = produceSqrtMat( *covMatS );
1201  TMatrixD* sqrtMatB = produceSqrtMat( *covMatB );
1202 
1203  // loop over species
1204  for (Int_t itype=0; itype<2; itype++) {
1205 
1206  Float_t* x;
1207  TMatrixD* m;
1208  if (itype == 0) { x = xS; m = sqrtMatS; cout << "- produce signal" << endl; }
1209  else { x = xB; m = sqrtMatB; cout << "- produce background" << endl; }
1210 
1211  // event loop
1212  TTree* tree = (itype==0) ? treeS : treeB;
1213  for (Int_t i=0; i<N; i++) {
1214 
1215  if (i%1000 == 0) cout << "... event: " << i << " (" << N << ")" << endl;
1216  getGaussRnd( *v, *m, R );
1217 
1218  for (Int_t ivar=0; ivar<nvar; ivar++) xvar[ivar] = (*v)[ivar] + x[ivar];
1219 
1220  tree->Fill();
1221  }
1222  }
1223 
1224  // write trees
1225  treeS->Write();
1226  treeB->Write();
1227 
1228  treeS->Show(0);
1229  treeB->Show(1);
1230 
1231  dataFile->Close();
1232  cout << "created data file: " << dataFile->GetName() << endl;
1233 
1234  plot();
1235 }
1236 
1237 // create the data
1238 void create_lin(Int_t N = 2000)
1239 {
1240  const Int_t nvar = 2;
1241  Double_t xvar[nvar];
1242  Float_t weight;
1243 
1244  // output flie
1245  TFile* dataFile = TFile::Open( "data.root", "RECREATE" );
1246 
1247  // create signal and background trees
1248  TTree* treeS = new TTree( "TreeS", "TreeS", 1 );
1249  TTree* treeB = new TTree( "TreeB", "TreeB", 1 );
1250  for (Int_t ivar=0; ivar<nvar; ivar++) {
1251  treeS->Branch( TString(Form( "var%i", ivar )).Data(), &xvar[ivar], TString(Form( "var%i/D", ivar )).Data() );
1252  treeB->Branch( TString(Form( "var%i", ivar )).Data(), &xvar[ivar], TString(Form( "var%i/D", ivar )).Data() );
1253  }
1254  treeS->Branch( "weight", &weight, "weight/F" );
1255  treeB->Branch( "weight", &weight, "weight/F" );
1256 
1257  TRandom R( 100 );
1258  Float_t xS[nvar] = { 0.0, 0.0 };
1259  Float_t xB[nvar] = { -0.0, -0.0 };
1260  Float_t dx[nvar] = { 1.0, 1.0 };
1261  TArrayD* v = new TArrayD( 2 );
1262  Float_t rhoS = 0.21;
1263  Float_t rhoB = 0.0;
1264 
1265  // create covariance matrix
1266  TMatrixD* covMatS = new TMatrixD( nvar, nvar );
1267  TMatrixD* covMatB = new TMatrixD( nvar, nvar );
1268  for (Int_t ivar=0; ivar<nvar; ivar++) {
1269  (*covMatS)(ivar,ivar) = dx[ivar]*dx[ivar];
1270  (*covMatB)(ivar,ivar) = dx[ivar]*dx[ivar];
1271  for (Int_t jvar=ivar+1; jvar<nvar; jvar++) {
1272  (*covMatS)(ivar,jvar) = rhoS*dx[ivar]*dx[jvar];
1273  (*covMatS)(jvar,ivar) = (*covMatS)(ivar,jvar);
1274 
1275  (*covMatB)(ivar,jvar) = rhoB*dx[ivar]*dx[jvar];
1276  (*covMatB)(jvar,ivar) = (*covMatB)(ivar,jvar);
1277  }
1278  }
1279  cout << "signal covariance matrix: " << endl;
1280  covMatS->Print();
1281  cout << "background covariance matrix: " << endl;
1282  covMatB->Print();
1283 
1284  // produce the square-root matrix
1285  TMatrixD* sqrtMatS = produceSqrtMat( *covMatS );
1286  TMatrixD* sqrtMatB = produceSqrtMat( *covMatB );
1287 
1288  // loop over species
1289  for (Int_t itype=0; itype<2; itype++) {
1290 
1291  Float_t* x;
1292  TMatrixD* m;
1293  if (itype == 0) { x = xS; m = sqrtMatS; cout << "- produce signal" << endl; }
1294  else { x = xB; m = sqrtMatB; cout << "- produce background" << endl; }
1295 
1296  // event loop
1297  TTree* tree = (itype==0) ? treeS : treeB;
1298  for (Int_t i=0; i<N; i++) {
1299 
1300  if (i%1000 == 0) cout << "... event: " << i << " (" << N << ")" << endl;
1301  getGaussRnd( *v, *m, R );
1302  for (Int_t ivar=0; ivar<nvar; ivar++) xvar[ivar] = (*v)[ivar] + x[ivar];
1303 
1304  // add weights
1305  if (itype == 0) weight = 1.0; // this is the signal weight
1306  else weight = 2.0; // this is the background weight
1307 
1308  tree->Fill();
1309  }
1310  }
1311 
1312  // write trees
1313  treeS->Write();
1314  treeB->Write();
1315 
1316  treeS->Show(0);
1317  treeB->Show(1);
1318 
1319  dataFile->Close();
1320  cout << "created data file: " << dataFile->GetName() << endl;
1321 
1322  plot();
1323 }
1324 
1325 // create the data
1326 
1327 void create_fullcirc(Int_t nmax = 20000, Bool_t distort=false)
1328 {
1329  TFile* dataFile = TFile::Open( "circledata.root", "RECREATE" );
1330  int nvar = 2;
1331  int nsig = 0, nbgd=0;
1332  Float_t weight=1;
1333  Float_t xvar[100];
1334  // create signal and background trees
1335  TTree* treeS = new TTree( "TreeS", "TreeS", 1 );
1336  TTree* treeB = new TTree( "TreeB", "TreeB", 1 );
1337  for (Int_t ivar=0; ivar<nvar; ivar++) {
1338  treeS->Branch( TString(Form( "var%i", ivar )).Data(), &xvar[ivar], TString(Form( "var%i/F", ivar)).Data() );
1339  treeB->Branch( TString(Form( "var%i", ivar )).Data(), &xvar[ivar], TString(Form( "var%i/F", ivar)).Data() );
1340  }
1341  treeS->Branch("weight", &weight, "weight/F");
1342  treeB->Branch("weight", &weight, "weight/F");
1343 
1344  TRandom R( 100 );
1345  do {
1346  for (Int_t ivar=0; ivar<nvar; ivar++) { xvar[ivar]=2.*R.Rndm()-1.;}
1347  Float_t xout = xvar[0]*xvar[0]+xvar[1]*xvar[1];
1348  if (nsig<10) cout << "xout = " << xout<<endl;
1349  if (xout < 0.3 || (xout >0.3 && xout<0.5 && R.Rndm() > xout)) {
1350  if (distort && xvar[0] < 0 && R.Rndm()>0.1) continue;
1351  treeS->Fill();
1352  nsig++;
1353  }
1354  else {
1355  if (distort && xvar[0] > 0 && R.Rndm()>0.1) continue;
1356  treeB->Fill();
1357  nbgd++;
1358  }
1359  } while ( nsig < nmax || nbgd < nmax);
1360 
1361  dataFile->Write();
1362  dataFile->Close();
1363 
1364 }
1365 
1366 // create the data
1367 void create_circ(Int_t N = 6000, Bool_t distort = false)
1368 {
1369  Int_t Nn = 0;
1370  const Int_t nvar = 2;
1371  Float_t xvar[nvar];
1372 
1373  // output flie
1374  TFile* dataFile = TFile::Open( "data.root", "RECREATE" );
1375 
1376  // create signal and background trees
1377  TTree* treeS = new TTree( "TreeS", "TreeS", 1 );
1378  TTree* treeB = new TTree( "TreeB", "TreeB", 1 );
1379  for (Int_t ivar=0; ivar<nvar; ivar++) {
1380  treeS->Branch( TString(Form( "var%i", ivar )).Data(), &xvar[ivar], TString(Form( "var%i/F", ivar )).Data() );
1381  treeB->Branch( TString(Form( "var%i", ivar )).Data(), &xvar[ivar], TString(Form( "var%i/F", ivar )).Data() );
1382  }
1383 // TTree *treeB = treeS->CloneTree();
1384 // for (Int_t ivar=0; ivar<nvar; ivar++) {
1385 // treeS->SetBranchAddress( Form( "var%i", ivar ), &xvar[ivar] );
1386 // treeB->SetBranchAddress( Form( "var%i", ivar ), &xvar[ivar] );
1387 // }
1388 // treeB->SetName ( "TreeB" );
1389 // treeB->SetTitle( "TreeB" );
1390 
1391  TRandom R( 100 );
1392  //Float_t phimin = -30, phimax = 130;
1393  Float_t phimin = -70, phimax = 130;
1394  Float_t phisig = 5;
1395  Float_t rS = 1.1;
1396  Float_t rB = 0.75;
1397  Float_t rsig = 0.1;
1398  Float_t fnmin = -(rS+4.0*rsig);
1399  Float_t fnmax = +(rS+4.0*rsig);
1400  Float_t dfn = fnmax-fnmin;
1401  // loop over species
1402  for (Int_t itype=0; itype<2; itype++) {
1403 
1404  // event loop
1405  TTree* tree = (itype==0) ? treeS : treeB;
1406  for (Int_t i=0; i<N; i++) {
1407  Double_t r1=R.Rndm(),r2=R.Rndm(), r3;
1408  if (itype==0) r3= r1>r2? r1 :r2;
1409  else r3= r2;
1410  Float_t phi;
1411  if (distort) phi = r3*(phimax - phimin) + phimin;
1412  else phi = R.Rndm()*(phimax - phimin) + phimin;
1413  phi += R.Gaus()*phisig;
1414 
1415  Float_t r = (itype==0) ? rS : rB;
1416  r += R.Gaus()*rsig;
1417 
1418  xvar[0] = r*cos(TMath::DegToRad()*phi);
1419  xvar[1] = r*sin(TMath::DegToRad()*phi);
1420 
1421  tree->Fill();
1422  }
1423 
1424  for (Int_t i=0; i<Nn; i++) {
1425 
1426  xvar[0] = dfn*R.Rndm()+fnmin;
1427  xvar[1] = dfn*R.Rndm()+fnmin;
1428 
1429  tree->Fill();
1430  }
1431  }
1432 
1433  // write trees
1434  treeS->Write();
1435  treeB->Write();
1436 
1437  treeS->Show(0);
1438  treeB->Show(1);
1439 
1440  dataFile->Close();
1441  cout << "created data file: " << dataFile->GetName() << endl;
1442 
1443  plot();
1444 }
1445 
1446 
1447 void create_schachbrett(Int_t nEvents = 20000) {
1448 
1449  const Int_t nvar = 2;
1450  Float_t xvar[nvar];
1451 
1452  // output flie
1453  TFile* dataFile = TFile::Open( "data.root", "RECREATE" );
1454 
1455  // create signal and background trees
1456  TTree* treeS = new TTree( "TreeS", "TreeS", 1 );
1457  TTree* treeB = new TTree( "TreeB", "TreeB", 1 );
1458  for (Int_t ivar=0; ivar<nvar; ivar++) {
1459  treeS->Branch( TString(Form( "var%i", ivar )).Data(), &xvar[ivar], TString(Form( "var%i/F", ivar )).Data() );
1460  treeB->Branch( TString(Form( "var%i", ivar )).Data(), &xvar[ivar], TString(Form( "var%i/F", ivar )).Data() );
1461  }
1462 
1463  Int_t nSeed = 12345;
1464  TRandom *m_rand = new TRandom(nSeed);
1465  Double_t sigma=0.3;
1466  Double_t meanX;
1467  Double_t meanY;
1468  Int_t xtype=1, ytype=1;
1469  Int_t iev=0;
1470  Int_t m_nDim = 2; // actually the boundary, there is a "bump" for every interger value
1471  // between in the Inteval [-m_nDim,m_nDim]
1472  while (iev < nEvents){
1473  xtype=1;
1474  for (Int_t i=-m_nDim; i <= m_nDim; i++){
1475  ytype = 1;
1476  for (Int_t j=-m_nDim; j <= m_nDim; j++){
1477  meanX=Double_t(i);
1478  meanY=Double_t(j);
1479  xvar[0]=m_rand->Gaus(meanY,sigma);
1480  xvar[1]=m_rand->Gaus(meanX,sigma);
1481  Int_t type = xtype*ytype;
1482  TTree* tree = (type==1) ? treeS : treeB;
1483  tree->Fill();
1484  iev++;
1485  ytype *= -1;
1486  }
1487  xtype *= -1;
1488  }
1489  }
1490 
1491 
1492  // write trees
1493  treeS->Write();
1494  treeB->Write();
1495 
1496  treeS->Show(0);
1497  treeB->Show(1);
1498 
1499  dataFile->Close();
1500  cout << "created data file: " << dataFile->GetName() << endl;
1501 
1502  plot();
1503 
1504 }
1505 
1506 
1507 void create_schachbrett_5D(Int_t nEvents = 200000) {
1508  const Int_t nvar = 5;
1509  Float_t xvar[nvar];
1510 
1511  // output flie
1512  TFile* dataFile = TFile::Open( "data.root", "RECREATE" );
1513 
1514  // create signal and background trees
1515  TTree* treeS = new TTree( "TreeS", "TreeS", 1 );
1516  TTree* treeB = new TTree( "TreeB", "TreeB", 1 );
1517  for (Int_t ivar=0; ivar<nvar; ivar++) {
1518  treeS->Branch( TString(Form( "var%i", ivar )).Data(), &xvar[ivar], TString(Form( "var%i/F", ivar )).Data() );
1519  treeB->Branch( TString(Form( "var%i", ivar )).Data(), &xvar[ivar], TString(Form( "var%i/F", ivar )).Data() );
1520  }
1521 
1522  Int_t nSeed = 12345;
1523  TRandom *m_rand = new TRandom(nSeed);
1524  Double_t sigma=0.3;
1525  Int_t itype[nvar];
1526  Int_t iev=0;
1527  Int_t m_nDim = 2; // actually the boundary, there is a "bump" for every interger value
1528  // between in the Inteval [-m_nDim,m_nDim]
1529 
1530  int idx[nvar];
1531  while (iev < nEvents){
1532  itype[0]=1;
1533  for (idx[0]=-m_nDim; idx[0] <= m_nDim; idx[0]++){
1534  itype[1]=1;
1535  for (idx[1]=-m_nDim; idx[1] <= m_nDim; idx[1]++){
1536  itype[2]=1;
1537  for (idx[2]=-m_nDim; idx[2] <= m_nDim; idx[2]++){
1538  itype[3]=1;
1539  for (idx[3]=-m_nDim; idx[3] <= m_nDim; idx[3]++){
1540  itype[4]=1;
1541  for (idx[4]=-m_nDim; idx[4] <= m_nDim; idx[4]++){
1542  Int_t type = itype[0];
1543  for (Int_t i=0;i<nvar;i++){
1544  xvar[i]=m_rand->Gaus(Double_t(idx[i]),sigma);
1545  if (i>0) type *= itype[i];
1546  }
1547  TTree* tree = (type==1) ? treeS : treeB;
1548  tree->Fill();
1549  iev++;
1550  itype[4] *= -1;
1551  }
1552  itype[3] *= -1;
1553  }
1554  itype[2] *= -1;
1555  }
1556  itype[1] *= -1;
1557  }
1558  itype[0] *= -1;
1559  }
1560  }
1561 
1562  // write trees
1563  treeS->Write();
1564  treeB->Write();
1565 
1566  treeS->Show(0);
1567  treeB->Show(1);
1568 
1569  dataFile->Close();
1570  cout << "created data file: " << dataFile->GetName() << endl;
1571 
1572  plot();
1573 
1574 }
1575 
1576 
1577 void create_schachbrett_4D(Int_t nEvents = 200000) {
1578 
1579  const Int_t nvar = 4;
1580  Float_t xvar[nvar];
1581 
1582  // output flie
1583  TFile* dataFile = TFile::Open( "data.root", "RECREATE" );
1584 
1585  // create signal and background trees
1586  TTree* treeS = new TTree( "TreeS", "TreeS", 1 );
1587  TTree* treeB = new TTree( "TreeB", "TreeB", 1 );
1588  for (Int_t ivar=0; ivar<nvar; ivar++) {
1589  treeS->Branch( TString(Form( "var%i", ivar )).Data(), &xvar[ivar], TString(Form( "var%i/F", ivar )).Data() );
1590  treeB->Branch( TString(Form( "var%i", ivar )).Data(), &xvar[ivar], TString(Form( "var%i/F", ivar )).Data() );
1591  }
1592 
1593  Int_t nSeed = 12345;
1594  TRandom *m_rand = new TRandom(nSeed);
1595  Double_t sigma=0.3;
1596  Int_t itype[nvar];
1597  Int_t iev=0;
1598  Int_t m_nDim = 2; // actually the boundary, there is a "bump" for every interger value
1599  // between in the Inteval [-m_nDim,m_nDim]
1600 
1601  int idx[nvar];
1602  while (iev < nEvents){
1603  itype[0]=1;
1604  for (idx[0]=-m_nDim; idx[0] <= m_nDim; idx[0]++){
1605  itype[1]=1;
1606  for (idx[1]=-m_nDim; idx[1] <= m_nDim; idx[1]++){
1607  itype[2]=1;
1608  for (idx[2]=-m_nDim; idx[2] <= m_nDim; idx[2]++){
1609  itype[3]=1;
1610  for (idx[3]=-m_nDim; idx[3] <= m_nDim; idx[3]++){
1611  Int_t type = itype[0];
1612  for (Int_t i=0;i<nvar;i++){
1613  xvar[i]=m_rand->Gaus(Double_t(idx[i]),sigma);
1614  if (i>0) type *= itype[i];
1615  }
1616  TTree* tree = (type==1) ? treeS : treeB;
1617  tree->Fill();
1618  iev++;
1619  itype[3] *= -1;
1620  }
1621  itype[2] *= -1;
1622  }
1623  itype[1] *= -1;
1624  }
1625  itype[0] *= -1;
1626  }
1627  }
1628 
1629  // write trees
1630  treeS->Write();
1631  treeB->Write();
1632 
1633  treeS->Show(0);
1634  treeB->Show(1);
1635 
1636  dataFile->Close();
1637  cout << "created data file: " << dataFile->GetName() << endl;
1638 
1639  plot();
1640 
1641 }
1642 
1643 
1644 void create_schachbrett_3D(Int_t nEvents = 20000) {
1645 
1646  const Int_t nvar = 3;
1647  Float_t xvar[nvar];
1648 
1649  // output flie
1650  TFile* dataFile = TFile::Open( "data.root", "RECREATE" );
1651 
1652  // create signal and background trees
1653  TTree* treeS = new TTree( "TreeS", "TreeS", 1 );
1654  TTree* treeB = new TTree( "TreeB", "TreeB", 1 );
1655  for (Int_t ivar=0; ivar<nvar; ivar++) {
1656  treeS->Branch( TString(Form( "var%i", ivar )).Data(), &xvar[ivar], TString(Form( "var%i/F", ivar )).Data() );
1657  treeB->Branch( TString(Form( "var%i", ivar )).Data(), &xvar[ivar], TString(Form( "var%i/F", ivar )).Data() );
1658  }
1659 
1660  Int_t nSeed = 12345;
1661  TRandom *m_rand = new TRandom(nSeed);
1662  Double_t sigma=0.3;
1663  Int_t itype[nvar];
1664  Int_t iev=0;
1665  Int_t m_nDim = 2; // actually the boundary, there is a "bump" for every interger value
1666  // between in the Inteval [-m_nDim,m_nDim]
1667 
1668  int idx[nvar];
1669  while (iev < nEvents){
1670  itype[0]=1;
1671  for (idx[0]=-m_nDim; idx[0] <= m_nDim; idx[0]++){
1672  itype[1]=1;
1673  for (idx[1]=-m_nDim; idx[1] <= m_nDim; idx[1]++){
1674  itype[2]=1;
1675  for (idx[2]=-m_nDim; idx[2] <= m_nDim; idx[2]++){
1676  Int_t type = itype[0];
1677  for (Int_t i=0;i<nvar;i++){
1678  xvar[i]=m_rand->Gaus(Double_t(idx[i]),sigma);
1679  if (i>0) type *= itype[i];
1680  }
1681  TTree* tree = (type==1) ? treeS : treeB;
1682  tree->Fill();
1683  iev++;
1684  itype[2] *= -1;
1685  }
1686  itype[1] *= -1;
1687  }
1688  itype[0] *= -1;
1689  }
1690  }
1691 
1692  // write trees
1693  treeS->Write();
1694  treeB->Write();
1695 
1696  treeS->Show(0);
1697  treeB->Show(1);
1698 
1699  dataFile->Close();
1700  cout << "created data file: " << dataFile->GetName() << endl;
1701 
1702  plot();
1703 
1704 }
1705 
1706 
1707 void create_schachbrett_2D(Int_t nEvents = 100000, Int_t nbumps=2) {
1708 
1709  const Int_t nvar = 2;
1710  Float_t xvar[nvar];
1711 
1712  // output flie
1713  TFile* dataFile = TFile::Open( "data.root", "RECREATE" );
1714 
1715  // create signal and background trees
1716  TTree* treeS = new TTree( "TreeS", "TreeS", 1 );
1717  TTree* treeB = new TTree( "TreeB", "TreeB", 1 );
1718  for (Int_t ivar=0; ivar<nvar; ivar++) {
1719  treeS->Branch( TString(Form( "var%i", ivar )).Data(), &xvar[ivar], TString(Form( "var%i/F", ivar )).Data() );
1720  treeB->Branch( TString(Form( "var%i", ivar )).Data(), &xvar[ivar], TString(Form( "var%i/F", ivar )).Data() );
1721  }
1722 
1723  Int_t nSeed = 345;
1724  TRandom *m_rand = new TRandom(nSeed);
1725  Double_t sigma=0.35;
1726  Int_t itype[nvar];
1727  Int_t iev=0;
1728  Int_t m_nDim = nbumps; // actually the boundary, there is a "bump" for every interger value
1729  // between in the Inteval [-m_nDim,m_nDim]
1730 
1731  int idx[nvar];
1732  while (iev < nEvents){
1733  itype[0]=1;
1734  for (idx[0]=-m_nDim; idx[0] <= m_nDim; idx[0]++){
1735  itype[1]=1;
1736  for (idx[1]=-m_nDim; idx[1] <= m_nDim; idx[1]++){
1737  Int_t type = itype[0];
1738  for (Int_t i=0;i<nvar;i++){
1739  xvar[i]=m_rand->Gaus(Double_t(idx[i]),sigma);
1740  if (i>0) type *= itype[i];
1741  }
1742  TTree* tree = (type==1) ? treeS : treeB;
1743  tree->Fill();
1744  iev++;
1745  itype[1] *= -1;
1746  }
1747  itype[0] *= -1;
1748  }
1749  }
1750 
1751  // write trees
1752  treeS->Write();
1753  treeB->Write();
1754 
1755  treeS->Show(0);
1756  treeB->Show(1);
1757 
1758  dataFile->Close();
1759  cout << "created data file: " << dataFile->GetName() << endl;
1760 
1761  plot();
1762 
1763 }
1764 
1765 
1766 
1767 void create_3Bumps(Int_t nEvents = 5000) {
1768  // signal is clustered around (1,0) and (-1,0) where one is two times(1,0)
1769  // bkg (0,0)
1770 
1771 
1772 
1773  const Int_t nvar = 2;
1774  Float_t xvar[nvar];
1775 
1776  // output flie
1777  TString filename = "data_3Bumps.root";
1778  TFile* dataFile = TFile::Open( filename, "RECREATE" );
1779 
1780  // create signal and background trees
1781  TTree* treeS = new TTree( "TreeS", "TreeS", 1 );
1782  TTree* treeB = new TTree( "TreeB", "TreeB", 1 );
1783  for (Int_t ivar=0; ivar<nvar; ivar++) {
1784  treeS->Branch( TString(Form( "var%i", ivar )).Data(), &xvar[ivar], TString(Form( "var%i/F", ivar )).Data() );
1785  treeB->Branch( TString(Form( "var%i", ivar )).Data(), &xvar[ivar], TString(Form( "var%i/F", ivar )).Data() );
1786  }
1787 
1788  Int_t nSeed = 12345;
1789  TRandom *m_rand = new TRandom(nSeed);
1790  Double_t sigma=0.2;
1791  Int_t type;
1792  Int_t iev=0;
1793  Double_t Centers[nvar][6] = {{-1,0,0,0,1,1},{0,0,0,0,0,0}}; //
1794 
1795 
1796  while (iev < nEvents){
1797  for (int idx=0; idx<6; idx++){
1798  if (idx==1 || idx==2 || idx==3) type = 0;
1799  else type=1;
1800  for (Int_t ivar=0;ivar<nvar;ivar++){
1801  xvar[ivar]=m_rand->Gaus(Centers[ivar][idx],sigma);
1802  }
1803  TTree* tree = (type==1) ? treeS : treeB;
1804  tree->Fill();
1805  iev++;
1806  }
1807  }
1808 
1809  // write trees
1810  treeS->Write();
1811  treeB->Write();
1812 
1813  treeS->Show(0);
1814  treeB->Show(1);
1815 
1816  dataFile->Close();
1817  cout << "created data file: " << dataFile->GetName() << endl;
1818 
1819  plot(filename);
1820 
1821 }
1822 
1823 void createOnionData(Int_t nmax = 50000){
1824  // output file
1825  TFile* dataFile = TFile::Open( "oniondata.root", "RECREATE" );
1826  int nvar = 4;
1827  int nsig = 0, nbgd=0;
1828  Float_t xvar[100];
1829  // create signal and background trees
1830  TTree* treeS = new TTree( "TreeS", "TreeS", 1 );
1831  TTree* treeB = new TTree( "TreeB", "TreeB", 1 );
1832  for (Int_t ivar=0; ivar<nvar; ivar++) {
1833  treeS->Branch( TString(Form( "var%i", ivar+1 )).Data(), &xvar[ivar], TString(Form( "var%i/F", ivar+1 )).Data() );
1834  treeB->Branch( TString(Form( "var%i", ivar+1 )).Data(), &xvar[ivar], TString(Form( "var%i/F", ivar+1 )).Data() );
1835  }
1836 
1837  TRandom R( 100 );
1838  do {
1839  for (Int_t ivar=0; ivar<nvar; ivar++) { xvar[ivar]=R.Rndm();}
1840  Float_t xout = sin(2.*acos(-1.)*(xvar[0]*xvar[1]*xvar[2]*xvar[3]+xvar[0]*xvar[1]));
1841  if (nsig<100) cout << "xout = " << xout<<endl;
1842  Int_t i = (Int_t) ((1.+xout)*4.99);
1843  if (i%2 == 0 && nsig < nmax) {
1844  treeS->Fill();
1845  nsig++;
1846  }
1847  if (i%2 != 0 && nbgd < nmax){
1848  treeB->Fill();
1849  nbgd++;
1850  }
1851  } while ( nsig < nmax || nbgd < nmax);
1852 
1853  dataFile->Write();
1854  dataFile->Close();
1855 }
1856 
1857 void create_multiclassdata(Int_t nmax = 20000)
1858 {
1859  TFile* dataFile = TFile::Open( "data.root", "RECREATE" );
1860  int ncls = 3;
1861  int nvar = 4;
1862  int ndat = 0;
1863  Int_t cls;
1864  Float_t thecls;
1865  Float_t weight=1;
1866  Float_t xcls[100];
1867  Float_t xmean[3][4] = {
1868  { 0. , 0.3, 0.5, 0.9 },
1869  { -0.2 , -0.3, 0.5, 0.4 },
1870  { 0.2 , 0.1, -0.1, 0.7 }} ;
1871 
1872  Float_t xvar[100];
1873  // create tree using class flag stored in int variable cls
1874  TTree* treeR = new TTree( "TreeR", "TreeR", 1 );
1875  for (Int_t ivar=0; ivar<nvar; ivar++) {
1876  treeR->Branch( TString(Form( "var%i", ivar )).Data(), &xvar[ivar], TString(Form( "var%i/F", ivar)).Data() );
1877  }
1878  for (Int_t icls=0; icls<ncls; icls++) {
1879  treeR->Branch(TString(Form( "cls%i", icls )).Data(), &xcls[icls], TString(Form( "cls%i/F", icls)).Data() );
1880  }
1881 
1882  treeR->Branch("cls", &thecls, "cls/F");
1883  treeR->Branch("weight", &weight, "weight/F");
1884 
1885  TRandom R( 100 );
1886  do {
1887  for (Int_t icls=0; icls<ncls; icls++) xcls[icls]=0.;
1888  cls = R.Integer(ncls);
1889  thecls = cls;
1890  xcls[cls]=1.;
1891  for (Int_t ivar=0; ivar<nvar; ivar++) {
1892  xvar[ivar]=R.Gaus(xmean[cls][ivar],1.);
1893  }
1894 
1895  if (ndat<30) cout << "cls=" << cls <<" xvar = " << xvar[0]<<" " <<xvar[1]<<" " << xvar[2]<<" " <<xvar[3]<<endl;
1896 
1897  treeR->Fill();
1898  ndat++;
1899  } while ( ndat < nmax );
1900 
1901  dataFile->Write();
1902  dataFile->Close();
1903 
1904 }
1905 
1906 
1907 
1908 
1909 
1910 
1911 // create the data
1912 void create_array_with_different_lengths(Int_t N = 100)
1913 {
1914  const Int_t nvar = 4;
1915  Int_t nvarCurrent = 4;
1916  Float_t xvar[nvar];
1917 
1918  // output flie
1919  TFile* dataFile = TFile::Open( "data.root", "RECREATE" );
1920 
1921  // create signal and background trees
1922  TTree* treeS = new TTree( "TreeS", "TreeS", 1 );
1923  TTree* treeB = new TTree( "TreeB", "TreeB", 1 );
1924  treeS->Branch( "arrSize", &nvarCurrent, "arrSize/I" );
1925  treeS->Branch( "arr", xvar, "arr[arrSize]/F" );
1926  treeB->Branch( "arrSize", &nvarCurrent, "arrSize/I" );
1927  treeB->Branch( "arr", xvar, "arr[arrSize]/F" );
1928 
1929  TRandom R( 100 );
1930  Float_t xS[nvar] = { 0.2, 0.3, 0.5, 0.9 };
1931  Float_t xB[nvar] = { -0.2, -0.3, -0.5, -0.6 };
1932  Float_t dx[nvar] = { 1.0, 1.0, 1.0, 1.0 };
1933  TArrayD* v = new TArrayD( nvar );
1934  Float_t rho[20];
1935  rho[1*2] = 0.4;
1936  rho[1*3] = 0.6;
1937  rho[1*4] = 0.9;
1938  rho[2*3] = 0.7;
1939  rho[2*4] = 0.8;
1940  rho[3*4] = 0.93;
1941 
1942  // create covariance matrix
1943  TMatrixD* covMatS = new TMatrixD( nvar, nvar );
1944  TMatrixD* covMatB = new TMatrixD( nvar, nvar );
1945  for (Int_t ivar=0; ivar<nvar; ivar++) {
1946  (*covMatS)(ivar,ivar) = dx[ivar]*dx[ivar];
1947  (*covMatB)(ivar,ivar) = dx[ivar]*dx[ivar];
1948  for (Int_t jvar=ivar+1; jvar<nvar; jvar++) {
1949  (*covMatS)(ivar,jvar) = rho[(ivar+1)*(jvar+1)]*dx[ivar]*dx[jvar];
1950  (*covMatS)(jvar,ivar) = (*covMatS)(ivar,jvar);
1951 
1952  (*covMatB)(ivar,jvar) = rho[(ivar+1)*(jvar+1)]*dx[ivar]*dx[jvar];
1953  (*covMatB)(jvar,ivar) = (*covMatB)(ivar,jvar);
1954  }
1955  }
1956  cout << "signal covariance matrix: " << endl;
1957  covMatS->Print();
1958  cout << "background covariance matrix: " << endl;
1959  covMatB->Print();
1960 
1961  // produce the square-root matrix
1962  TMatrixD* sqrtMatS = produceSqrtMat( *covMatS );
1963  TMatrixD* sqrtMatB = produceSqrtMat( *covMatB );
1964 
1965  // loop over species
1966  for (Int_t itype=0; itype<2; itype++) {
1967 
1968  Float_t* x;
1969  TMatrixD* m;
1970  if (itype == 0) { x = xS; m = sqrtMatS; cout << "- produce signal" << endl; }
1971  else { x = xB; m = sqrtMatB; cout << "- produce background" << endl; }
1972 
1973  // event loop
1974  TTree* tree = (itype==0) ? treeS : treeB;
1975  for (Int_t i=0; i<N; i++) {
1976 
1977  if (i%1000 == 0) cout << "... event: " << i << " (" << N << ")" << endl;
1978  getGaussRnd( *v, *m, R );
1979 
1980  for (Int_t ivar=0; ivar<nvar; ivar++) xvar[ivar] = (*v)[ivar] + x[ivar];
1981 
1982 
1983  nvarCurrent = (i%4)+1;
1984 
1985  tree->Fill();
1986  }
1987  }
1988 
1989  // write trees
1990  treeS->Write();
1991  treeB->Write();
1992 
1993  treeS->Show(0);
1994  treeB->Show(1);
1995 
1996  dataFile->Close();
1997  cout << "created data file: " << dataFile->GetName() << endl;
1998 }
1999 
2000 
2001 
2002 // create the data
2003 void create_MultipleBackground(Int_t N = 50000)
2004 {
2005  const int nvar = 4;
2006 
2007  // output flie
2008  TFile* dataFile = TFile::Open( "tmva_example_multiple_background.root", "RECREATE" );
2009 
2010 
2011  Float_t xS[nvar] = { 0.2, 0.3, 0.5, 0.9 };
2012  Float_t xB0[nvar] = { -0.2, -0.3, -0.5, -0.6 };
2013  Float_t xB1[nvar] = { -0.2, 0.3, 0.5, -0.6 };
2014  Float_t dx0[nvar] = { 1.0, 1.0, 1.0, 1.0 };
2015  Float_t dx1[nvar] = { -1.0, -1.0, -1.0, -1.0 };
2016 
2017  // create signal and background trees
2018  TTree* treeS = makeTree_lin_Nvar( "TreeS", "Signal tree", xS, dx0, nvar, N );
2019  TTree* treeB0 = makeTree_lin_Nvar( "TreeB0", "Background 0", xB0, dx0, nvar, N );
2020  TTree* treeB1 = makeTree_lin_Nvar( "TreeB1", "Background 1", xB1, dx1, nvar, N );
2021  TTree* treeB2 = makeTree_circ( "TreeB2", "Background 2", nvar, N, 1.5, true);
2022 
2023  treeS->Write();
2024  treeB0->Write();
2025  treeB1->Write();
2026  treeB2->Write();
2027 
2028  //treeS->Show(0);
2029  //treeB0->Show(0);
2030  //treeB1->Show(0);
2031  //treeB2->Show(0);
2032 
2033  dataFile->Close();
2034  cout << "created data file: " << dataFile->GetName() << endl;
2035 }