ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/benhoob/HWW/createData.C
Revision: 1.1
Committed: Mon Feb 14 12:39:14 2011 UTC (14 years, 3 months ago) by benhoob
Content type: text/plain
Branch: MAIN
CVS Tags: HEAD
Log Message: 
Initial commit

File Contents

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