Logo ROOT  

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