Logo ROOT   6.18/05
Reference Guide
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
24void 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
104TMatrixD* 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
129void 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
157void 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
262TTree* 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
320TTree* 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
382void 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
412void 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)
498void 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
601void 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
772void 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
870void 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
970void 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
1068void 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
1149void 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
1238void 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
1327void 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
1367void 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
1447void 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
1507void 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
1577void 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
1644void 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
1707void 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
1767void 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
1823void 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
1857void 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
1912void 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
2003void 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}
v
SVector< double, 2 > v
Definition: Dict.h:5
r
ROOT::R::TRInterface & r
Definition: Object.C:4
c
#define c(i)
Definition: RSha256.hxx:101
h
#define h(i)
Definition: RSha256.hxx:106
R
#define R(a, b, c, d, e, f, g, h, i)
Definition: RSha256.hxx:110
Int_t
int Int_t
Definition: RtypesCore.h:41
Bool_t
bool Bool_t
Definition: RtypesCore.h:59
Double_t
double Double_t
Definition: RtypesCore.h:55
Float_t
float Float_t
Definition: RtypesCore.h:53
N
#define N
type
int type
Definition: TGX11.cxx:120
xmin
float xmin
Definition: THbookFile.cxx:93
ymin
float ymin
Definition: THbookFile.cxx:93
xmax
float xmax
Definition: THbookFile.cxx:93
ymax
float ymax
Definition: THbookFile.cxx:93
acos
double acos(double)
cos
double cos(double)
sin
double sin(double)
TMatrixD
TMatrixT< Double_t > TMatrixD
Definition: TMatrixDfwd.h:22
gROOT
#define gROOT
Definition: TROOT.h:414
gRandom
R__EXTERN TRandom * gRandom
Definition: TRandom.h:62
Form
char * Form(const char *fmt,...)
TArrayD
Array of doubles (64 bits per element).
Definition: TArrayD.h:27
TAttFill::SetFillStyle
virtual void SetFillStyle(Style_t fstyle)
Set the fill area style.
Definition: TAttFill.h:39
TAttLine::SetLineColor
virtual void SetLineColor(Color_t lcolor)
Set the line color.
Definition: TAttLine.h:40
TAttMarker::SetMarkerColor
virtual void SetMarkerColor(Color_t mcolor=1)
Set the marker color.
Definition: TAttMarker.h:38
TAttMarker::SetMarkerStyle
virtual void SetMarkerStyle(Style_t mstyle=1)
Set the marker style.
Definition: TAttMarker.h:40
TAttMarker::SetMarkerSize
virtual void SetMarkerSize(Size_t msize=1)
Set the marker size.
Definition: TAttMarker.h:41
TCanvas
The Canvas class.
Definition: TCanvas.h:31
TDirectoryFile::Get
virtual TObject * Get(const char *namecycle)
Return pointer to object identified by namecycle.
Definition: TDirectoryFile.cxx:912
TFile
A ROOT file is a suite of consecutive data records (TKey instances) with a well defined format.
Definition: TFile.h:48
TFile::Close
virtual void Close(Option_t *option="")
Close a file.
Definition: TFile.cxx:914
TFile::Write
virtual Int_t Write(const char *name=0, Int_t opt=0, Int_t bufsiz=0)
Write memory objects to this file.
Definition: TFile.cxx:2337
TFile::Open
static TFile * Open(const char *name, Option_t *option="", const char *ftitle="", Int_t compress=ROOT::RCompressionSetting::EDefaults::kUseGeneralPurpose, Int_t netopt=0)
Create / open a file.
Definition: TFile.cxx:3980
TH1F
1-D histogram with a float per channel (see TH1 documentation)}
Definition: TH1.h:571
TH1::SetTitle
virtual void SetTitle(const char *title)
See GetStatOverflows for more information.
Definition: TH1.cxx:6309
TH1::SetTitleSize
virtual void SetTitleSize(Float_t size=0.02, Option_t *axis="X")
Set the axis' title size.
Definition: Haxis.cxx:365
TH1::GetXaxis
TAxis * GetXaxis()
Get the behaviour adopted by the object about the statoverflows. See EStatOverflows for more informat...
Definition: TH1.h:316
TH1::Fill
virtual Int_t Fill(Double_t x)
Increment bin with abscissa X by 1.
Definition: TH1.cxx:3258
TH1::GetYaxis
TAxis * GetYaxis()
Definition: TH1.h:317
TH1::Draw
virtual void Draw(Option_t *option="")
Draw this histogram with options.
Definition: TH1.cxx:2981
TH1::SetLabelSize
virtual void SetLabelSize(Float_t size=0.02, Option_t *axis="X")
Set size of axis' labels.
Definition: Haxis.cxx:285
TH2F
2-D histogram with a float per channel (see TH1 documentation)}
Definition: TH2.h:248
TLegend
This class displays a legend box (TPaveText) containing several legend entries.
Definition: TLegend.h:23
TLegend::AddEntry
TLegendEntry * AddEntry(const TObject *obj, const char *label="", Option_t *option="lpf")
Add a new entry to this legend.
Definition: TLegend.cxx:330
TLegend::Draw
virtual void Draw(Option_t *option="")
Draw this legend with its current attributes.
Definition: TLegend.cxx:423
TLegend::SetMargin
void SetMargin(Float_t margin)
Definition: TLegend.h:69
TMatrixTBase::GetNrows
Int_t GetNrows() const
Definition: TMatrixTBase.h:124
TMatrixTBase::Print
void Print(Option_t *name="") const
Print the matrix as a table of elements.
Definition: TMatrixTBase.cxx:832
TNamed::SetTitle
virtual void SetTitle(const char *title="")
Set the title of the TNamed.
Definition: TNamed.cxx:164
TNamed::GetName
virtual const char * GetName() const
Returns name of object.
Definition: TNamed.h:47
TPave::SetBorderSize
virtual void SetBorderSize(Int_t bordersize=4)
Definition: TPave.h:73
TRandom
This is the base class for the ROOT Random number generators.
Definition: TRandom.h:27
TRandom::Gaus
virtual Double_t Gaus(Double_t mean=0, Double_t sigma=1)
Samples a random number from the standard Normal (Gaussian) Distribution with the given mean and sigm...
Definition: TRandom.cxx:263
TRandom::Rndm
virtual Double_t Rndm()
Machine independent random number generator.
Definition: TRandom.cxx:541
TString
Basic string class.
Definition: TString.h:131
TString::Data
const char * Data() const
Definition: TString.h:364
TStyle
TStyle objects may be created to define special styles.
Definition: TStyle.h:27
TStyle::SetOptTitle
void SetOptTitle(Int_t tit=1)
Definition: TStyle.h:312
TStyle::SetPadTopMargin
void SetPadTopMargin(Float_t margin=0.1)
Definition: TStyle.h:336
TStyle::SetTitleX
void SetTitleX(Float_t x=0)
Definition: TStyle.h:390
TStyle::SetOptStat
void SetOptStat(Int_t stat=1)
The type of information printed in the histogram statistics box can be selected via the parameter mod...
Definition: TStyle.cxx:1444
TStyle::SetPadBottomMargin
void SetPadBottomMargin(Float_t margin=0.1)
Definition: TStyle.h:335
TStyle::SetPadRightMargin
void SetPadRightMargin(Float_t margin=0.1)
Definition: TStyle.h:338
TStyle::SetPadGridX
void SetPadGridX(Bool_t gridx)
Definition: TStyle.h:339
TStyle::SetPadLeftMargin
void SetPadLeftMargin(Float_t margin=0.1)
Definition: TStyle.h:337
TStyle::SetPadGridY
void SetPadGridY(Bool_t gridy)
Definition: TStyle.h:340
TStyle::cd
virtual void cd()
Change current style.
Definition: TStyle.cxx:438
TStyle::SetTitleW
void SetTitleW(Float_t w=0)
Definition: TStyle.h:392
TStyle::SetTitleH
void SetTitleH(Float_t h=0)
Definition: TStyle.h:393
TStyle::SetTitleY
void SetTitleY(Float_t y=0.985)
Definition: TStyle.h:391
TTree
A TTree represents a columnar dataset.
Definition: TTree.h:71
TTree::Fill
virtual Int_t Fill()
Fill all branches.
Definition: TTree.cxx:4419
TTree::Show
virtual void Show(Long64_t entry=-1, Int_t lenmax=20)
Print values of all active leaves for entry.
Definition: TTree.cxx:8981
TTree::GetMaximum
virtual Double_t GetMaximum(const char *columname)
Return maximum of column with name columname.
Definition: TTree.cxx:6013
TTree::GetEntries
virtual Long64_t GetEntries() const
Definition: TTree.h:402
TTree::GetEntry
virtual Int_t GetEntry(Long64_t entry=0, Int_t getall=0)
Read all branches of entry and return total number of bytes read.
Definition: TTree.cxx:5422
TTree::GetMinimum
virtual Double_t GetMinimum(const char *columname)
Return minimum of column with name columname.
Definition: TTree.cxx:6053
TTree::Draw
virtual void Draw(Option_t *opt)
Default Draw method for all objects.
Definition: TTree.h:371
TTree::Branch
virtual Int_t Branch(TCollection *list, Int_t bufsize=32000, Int_t splitlevel=99, const char *name="")
Create one branch for each element in the collection.
Definition: TTree.cxx:1741
TTree::Write
virtual Int_t Write(const char *name=0, Int_t option=0, Int_t bufsize=0)
Write this object to the current directory.
Definition: TTree.cxx:9370
sigma
const Double_t sigma
Definition: h1analysisProxy.h:11
y
Double_t y[n]
Definition: legend1.C:17
x
Double_t x[n]
Definition: legend1.C:17
TMath::Gaus
Double_t Gaus(Double_t x, Double_t mean=0, Double_t sigma=1, Bool_t norm=kFALSE)
Calculate a gaussian function with mean and sigma.
Definition: TMath.cxx:448
TMath::Nint
Int_t Nint(T x)
Round to nearest integer. Rounds half integers to the nearest even integer.
Definition: TMath.h:701
TMath::Max
Short_t Max(Short_t a, Short_t b)
Definition: TMathBase.h:212
TMath::Log
Double_t Log(Double_t x)
Definition: TMath.h:748
TMath::DegToRad
constexpr Double_t DegToRad()
Conversion from degree to radian:
Definition: TMath.h:82
TMath::Sqrt
Double_t Sqrt(Double_t x)
Definition: TMath.h:679
TMath::Min
Short_t Min(Short_t a, Short_t b)
Definition: TMathBase.h:180
TMath::Pi
constexpr Double_t Pi()
Definition: TMath.h:38
TMath::Sin
Double_t Sin(Double_t)
Definition: TMath.h:625
TMath::Abs
Short_t Abs(Short_t d)
Definition: TMathBase.h:120
tree
Definition: tree.py:1
m
auto * m
Definition: textangle.C:8
l
auto * l
Definition: textangle.C:4
sum
static long int sum(long int i)
Definition: Factory.cxx:2258