Logo ROOT  
Reference Guide
TFoam.cxx
Go to the documentation of this file.
1// @(#)root/foam:$Id$
2// Author: S. Jadach <mailto:Stanislaw.jadach@ifj.edu.pl>, P.Sawicki <mailto:Pawel.Sawicki@ifj.edu.pl>
3
4/** \class TFoam
5
6
7TFoam is the main class of the multi-dimensional general purpose
8Monte Carlo event generator (integrator) FOAM.
9
10### FOAM Version 1.02M
11
12\authors
13 S. Jadach and P.Sawicki
14 Institute of Nuclear Physics, Cracow, Poland
15 Stanislaw. Jadach@ifj.edu.pl, Pawel.Sawicki@ifj.edu.pl
16
17### What is FOAM for?
18
19 - Suppose you want to generate randomly points (vectors) according to
20 an arbitrary probability distribution in n dimensions,
21 for which you supply your own method. FOAM can do it for you!
22 Even if your distributions has quite strong peaks and is discontinuous!
23 - FOAM generates random points with weight one or with variable weight.
24 - FOAM is capable to integrate using efficient "adaptive" MC method.
25 (The distribution does not need to be normalized to one.)
26
27### How does it work?
28
29FOAM is the simplified version of the multi-dimensional general purpose
30Monte Carlo event generator (integrator) FOAM.
31It creates hyper-rectangular "foam of cells", which is more dense around its peaks.
32See the following 2-dim. example of the map of 1000 cells for doubly peaked distribution:
33
34\image html foam_MapCamel1000.png width=400
35
36FOAM is now fully integrated with the ROOT package.
37The important bonus of the ROOT use is persistency of the FOAM objects!
38
39For more sophisticated problems full version of FOAM may be more appropriate:
40See [full version of FOAM](http://jadach.home.cern.ch/jadach/Foam/Index.html)
41
42### Simple example of the use of FOAM:
43
44Begin_Macro(source)
45../../../tutorials/foam/foam_kanwa.C
46End_Macro
47
48### Canonical nine steering parameters of FOAM
49
50
51| Name | default | Description |
52|----------|----------|--------------------------------------------------------|
53| kDim | 0 | Dimension of the integration space. Must be redefined! |
54| nCells | 1000 | No of allocated number of cells, |
55| nSampl | 200 | No. of MC events in the cell MC exploration |
56| nBin | 8 | No. of bins in edge-histogram in cell exploration |
57| OptRej | 1 | OptRej = 0, weighted; OptRej=1, wt=1 MC events |
58| OptDrive | 2 | Maximum weight reduction, =1 for variance reduction |
59| EvPerBin | 25 | Maximum number of the effective wt=1 events/bin, |
60| | | EvPerBin=0 deactivates this option |
61| Chat | 1 | =0,1,2 is the ``chat level'' in the standard output |
62| MaxWtRej | 1.1 | Maximum weight used to get w=1 MC events |
63
64The above can be redefined before calling Initialize() method,
65for instance `FoamObject->SetkDim(15)` sets dimension of the distribution to 15.
66Only `kDim` HAS TO BE redefined, the other parameters may be left at their defaults.
67`nCell` may be increased up to about million cells for wildly peaked distributions.
68Increasing `nSampl` sometimes helps, but it may cost CPU time.
69`MaxWtRej` may need to be increased for wild a distribution, while using `OptRej=0`.
70
71Past versions of FOAM: August 2003, v.1.00; September 2003 v.1.01
72Adopted starting from FOAM-2.06 by P. Sawicki
73
74Users of FOAM are kindly requested to cite the following work:
75S. Jadach, Computer Physics Communications 152 (2003) 55.
76
77*/
78
79#include "TFoam.h"
80#include "TFoamIntegrand.h"
81#include "TFoamMaxwt.h"
82#include "TFoamVect.h"
83#include "TFoamCell.h"
84#include <iostream>
85#include <iomanip>
86#include <fstream>
87#include "TH1.h"
88#include "TRefArray.h"
89#include "TObjArray.h"
90#include "TMethodCall.h"
91#include "TRandom.h"
92#include "TMath.h"
93#include "TInterpreter.h"
94
96
97//FFFFFF BoX-FORMATs for nice and flexible outputs
98#define BXOPE std::cout<<\
99"FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF"<<std::endl<<\
100"F F"<<std::endl
101#define BXTXT(text) std::cout<<\
102"F "<<std::setw(40)<< text <<" F"<<std::endl
103#define BX1I(name,numb,text) std::cout<<\
104"F "<<std::setw(10)<<name<<" = "<<std::setw(10)<<numb<<" = " <<std::setw(50)<<text<<" F"<<std::endl
105#define BX1F(name,numb,text) std::cout<<"F "<<std::setw(10)<<name<<\
106 " = "<<std::setw(15)<<std::setprecision(8)<<numb<<" = "<<std::setw(40)<<text<<" F"<<std::endl
107#define BX2F(name,numb,err,text) std::cout<<"F "<<std::setw(10)<<name<<\
108" = "<<std::setw(15)<<std::setprecision(8)<<numb<<" +- "<<std::setw(15)<<std::setprecision(8)<<err<< \
109 " = "<<std::setw(25)<<text<<" F"<<std::endl
110#define BXCLO std::cout<<\
111"F F"<<std::endl<<\
112"FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF"<<std::endl
113 //FFFFFF BoX-FORMATs ends here
114
115static const Double_t gHigh= 1.0e150;
116static const Double_t gVlow=-1.0e150;
117
118#define SW2 setprecision(7) << std::setw(12)
119
120// class to wrap a global function in a TFoamIntegrand function
122
123public:
124
126
128
130
131 // evaluate the density using the provided function pointer
132 Double_t Density (Int_t nDim, Double_t * x) override {
133 return fFunc(nDim,x);
134 }
135
136private:
137
139
140};
141
142
143////////////////////////////////////////////////////////////////////////////////
144/// Default constructor for streamer, user should not use it.
145
147 fDim(0), fNCells(0), fRNmax(0),
148 fOptDrive(0), fChat(0), fOptRej(0),
149 fNBin(0), fNSampl(0), fEvPerBin(0),
150 fMaskDiv(0), fInhiDiv(0), fOptPRD(0), fXdivPRD(0),
151 fNoAct(0), fLastCe(0), fCells(0),
152 fMCMonit(0), fMaxWtRej(0), fPrimAcu(0),
153 fHistEdg(0), fHistDbg(0), fHistWt(0),
154 fMCvect(0), fMCwt(0), fRvec(0),
155 fRho(0), fMethodCall(0), fPseRan(0),
156 fNCalls(0), fNEffev(0),
157 fSumWt(0), fSumWt2(0),
158 fSumOve(0), fNevGen(0),
159 fWtMax(0), fWtMin(0),
160 fPrime(0), fMCresult(0), fMCerror(0),
161 fAlpha(0)
162{
163}
164////////////////////////////////////////////////////////////////////////////////
165/// User constructor, to be employed by the user
166
168 fDim(0), fNCells(0), fRNmax(0),
169 fOptDrive(0), fChat(0), fOptRej(0),
170 fNBin(0), fNSampl(0), fEvPerBin(0),
171 fMaskDiv(0), fInhiDiv(0), fOptPRD(0), fXdivPRD(0),
172 fNoAct(0), fLastCe(0), fCells(0),
173 fMCMonit(0), fMaxWtRej(0), fPrimAcu(0),
174 fHistEdg(0), fHistDbg(0), fHistWt(0),
175 fMCvect(0), fMCwt(0), fRvec(0),
176 fRho(0), fMethodCall(0), fPseRan(0),
177 fNCalls(0), fNEffev(0),
178 fSumWt(0), fSumWt2(0),
179 fSumOve(0), fNevGen(0),
180 fWtMax(0), fWtMin(0),
181 fPrime(0), fMCresult(0), fMCerror(0),
182 fAlpha(0)
183{
184 if(strlen(Name) >129) {
185 Error("TFoam","Name too long %s \n",Name);
186 }
187 fName=Name; // Class name
188 fDate=" Release date: 2005.04.10"; // Release date
189 fVersion= "1.02M"; // Release version
190 fMaskDiv = 0; // Dynamic Mask for cell division, h-cubic
191 fInhiDiv = 0; // Flag allowing to inhibit cell division in certain projection/edge
192 fXdivPRD = 0; // Lists of division values encoded in one vector per direction
193 fCells = 0;
194 fAlpha = 0;
195 fPrimAcu = 0;
196 fHistEdg = 0;
197 fHistWt = 0;
198 fHistDbg = 0;
199 fDim = 0; // dimension of hyp-cubical space
200 fNCells = 1000; // Maximum number of Cells, is usually re-set
201 fNSampl = 200; // No of sampling when dividing cell
202 fOptPRD = 0; // General Option switch for PRedefined Division, for quick check
203 fOptDrive = 2; // type of Drive =1,2 for TrueVol,Sigma,WtMax
204 fChat = 1; // Chat=0,1,2 chat level in output, Chat=1 normal level
205 fOptRej = 1; // OptRej=0, wted events; OptRej=1, wt=1 events
206 /////////////////////////////////////////////////////////////////////////////
207
208 fNBin = 8; // binning of edge-histogram in cell exploration
209 fEvPerBin =25; // maximum no. of EFFECTIVE event per bin, =0 option is inactive
210 //------------------------------------------------------
211 fNCalls = 0; // No of function calls
212 fNEffev = 0; // Total no of eff. wt=1 events in build=up
213 fLastCe =-1; // Index of the last cell
214 fNoAct = 0; // No of active cells (used in MC generation)
215 fWtMin = gHigh; // Minimal weight
216 fWtMax = gVlow; // Maximal weight
217 fMaxWtRej =1.10; // Maximum weight in rejection for getting wt=1 events
218 fPseRan = 0; // Initialize private copy of random number generator
219 fMCMonit = 0; // MC efficiency monitoring
220 fRho = 0; // pointer to abstract class providing function to integrate
221 fMCvect = 0;
222 fRvec = 0;
223 fPseRan = 0; // generator of pseudorandom numbers
224 fMethodCall=0; // ROOT's pointer to global distribution function
225}
226
227////////////////////////////////////////////////////////////////////////////////
228/// Default destructor
229
231{
232 Int_t i;
233
234 if(fCells!= 0) {
235 for(i=0; i<fNCells; i++) delete fCells[i]; // TFoamCell*[]
236 delete [] fCells;
237 }
238 if (fRvec) delete [] fRvec; //double[]
239 if (fAlpha) delete [] fAlpha; //double[]
240 if (fMCvect) delete [] fMCvect; //double[]
241 if (fPrimAcu) delete [] fPrimAcu; //double[]
242 if (fMaskDiv) delete [] fMaskDiv; //int[]
243 if (fInhiDiv) delete [] fInhiDiv; //int[]
244
245 if( fXdivPRD!= 0) {
246 for(i=0; i<fDim; i++) delete fXdivPRD[i]; // TFoamVect*[]
247 delete [] fXdivPRD;
248 }
249 if (fMCMonit) delete fMCMonit;
250 if (fHistWt) delete fHistWt;
251
252 // delete histogram arrays
253 if (fHistEdg) {
254 fHistEdg->Delete();
255 delete fHistEdg;
256 }
257 if (fHistDbg) {
258 fHistDbg->Delete();
259 delete fHistDbg;
260 }
261 // delete function object if it has been created here in SetRhoInt
262 if (fRho && dynamic_cast<FoamIntegrandFunction*>(fRho) ) delete fRho;
263}
264
265
266////////////////////////////////////////////////////////////////////////////////
267/// Copy Constructor NOT IMPLEMENTED (NEVER USED)
268
269TFoam::TFoam(const TFoam &From): TObject(From)
270{
271 Error("TFoam", "COPY CONSTRUCTOR NOT IMPLEMENTED \n");
272}
273
274////////////////////////////////////////////////////////////////////////////////
275/// ### Basic initialization of FOAM invoked by the user. Mandatory!
276///
277/// This method starts the process of the cell build-up.
278/// User must invoke Initialize with two arguments or Initialize without arguments.
279/// This is done BEFORE generating first MC event and AFTER allocating FOAM object
280/// and reseting (optionally) its internal parameters/switches.
281/// The overall operational scheme of the FOAM is the following:
282///
283/// \image html foam_schema2.png width=600
284///
285/// ### This method invokes several other methods:
286///
287/// InitCells initializes memory storage for cells and begins exploration process
288/// from the root cell. The empty cells are allocated/filled using CellFill.
289/// The procedure Grow which loops over cells, picks up the cell with the biggest
290/// ``driver integral'', see Comp. Phys. Commun. 152 152 (2003) 55 for explanations,
291/// with the help of PeekMax procedure. The chosen cell is split using Divide.
292/// Subsequently, the procedure Explore called by the Divide
293/// (and by InitCells for the root cell) does the most important
294/// job in the FOAM object build-up: it performs a small MC run for each
295/// newly allocated daughter cell.
296/// Explore calculates how profitable the future split of the cell will be
297/// and defines the optimal cell division geometry with the help of Carver or Varedu
298/// procedures, for maximum weight or variance optimization respectively.
299/// All essential results of the exploration are written into
300/// the explored cell object. At the very end of the foam build-up,
301/// Finally, MakeActiveList is invoked to create a list of pointers to
302/// all active cells, for the purpose of the quick access during the MC generation.
303/// The procedure Explore employs MakeAlpha to generate random coordinates
304/// inside a given cell with the uniform distribution.
305/// The above sequence of the procedure calls is depicted in the following figure:
306///
307/// \image html foam_Initialize_schema.png width=600
308
310{
311
312
313 SetPseRan(PseRan);
314 SetRho(fun);
315 Initialize();
316}
317
318////////////////////////////////////////////////////////////////////////////////
319/// Basic initialization of FOAM invoked by the user.
320/// IMPORTANT: Random number generator and the distribution object has to be
321/// provided using SetPseRan and SetRho prior to invoking this initialiser!
322
324{
325 Bool_t addStatus = TH1::AddDirectoryStatus();
327 Int_t i;
328
329 if(fChat>0){
330 BXOPE;
331 BXTXT("****************************************");
332 BXTXT("****** TFoam::Initialize ******");
333 BXTXT("****************************************");
334 BXTXT(fName);
335 BX1F(" Version",fVersion, fDate);
336 BX1I(" kDim",fDim, " Dimension of the hyper-cubical space ");
337 BX1I(" nCells",fNCells, " Requested number of Cells (half of them active) ");
338 BX1I(" nSampl",fNSampl, " No of MC events in exploration of a cell ");
339 BX1I(" nBin",fNBin, " No of bins in histograms, MC exploration of cell ");
340 BX1I(" EvPerBin",fEvPerBin, " Maximum No effective_events/bin, MC exploration ");
341 BX1I(" OptDrive",fOptDrive, " Type of Driver =1,2 for Sigma,WtMax ");
342 BX1I(" OptRej",fOptRej, " MC rejection on/off for OptRej=0,1 ");
343 BX1F(" MaxWtRej",fMaxWtRej, " Maximum wt in rejection for wt=1 evts");
344 BXCLO;
345 }
346
347 if(fPseRan==0) Error("Initialize", "Random number generator not set \n");
348 if(fRho==0 && fMethodCall==0 ) Error("Initialize", "Distribution function not set \n");
349 if(fDim==0) Error("Initialize", "Zero dimension not allowed \n");
350
351 /////////////////////////////////////////////////////////////////////////
352 // ALLOCATE SMALL LISTS //
353 // it is done globally, not for each cell, to save on allocation time //
354 /////////////////////////////////////////////////////////////////////////
355 fRNmax= fDim+1;
356 fRvec = new Double_t[fRNmax]; // Vector of random numbers
357 if(fRvec==0) Error("Initialize", "Cannot initialize buffer fRvec \n");
358
359 if(fDim>0){
360 fAlpha = new Double_t[fDim]; // sum<1 for internal parametrization of the simplex
361 if(fAlpha==0) Error("Initialize", "Cannot initialize buffer fAlpha \n" );
362 }
363 fMCvect = new Double_t[fDim]; // vector generated in the MC run
364 if(fMCvect==0) Error("Initialize", "Cannot initialize buffer fMCvect \n" );
365
366 //====== List of directions inhibited for division
367 if(fInhiDiv == 0){
368 fInhiDiv = new Int_t[fDim];
369 for(i=0; i<fDim; i++) fInhiDiv[i]=0;
370 }
371 //====== Dynamic mask used in Explore for edge determination
372 if(fMaskDiv == 0){
373 fMaskDiv = new Int_t[fDim];
374 for(i=0; i<fDim; i++) fMaskDiv[i]=1;
375 }
376 //====== List of predefined division values in all directions (initialized as empty)
377 if(fXdivPRD == 0){
378 fXdivPRD = new TFoamVect*[fDim];
379 for(i=0; i<fDim; i++) fXdivPRD[i]=0; // Artificially extended beyond fDim
380 }
381 //====== Initialize list of histograms
382 fHistWt = new TH1D("HistWt","Histogram of MC weight",100,0.0, 1.5*fMaxWtRej); // MC weight
383 fHistEdg = new TObjArray(fDim); // Initialize list of histograms
384 TString hname;
385 TString htitle;
386 for(i=0;i<fDim;i++){
387 hname=fName+TString("_HistEdge_");
388 hname+=i;
389 htitle=TString("Edge Histogram No. ");
390 htitle+=i;
391 //std::cout<<"i= "<<i<<" hname= "<<hname<<" htitle= "<<htitle<<std::endl;
392 (*fHistEdg)[i] = new TH1D(hname.Data(),htitle.Data(),fNBin,0.0, 1.0); // Initialize histogram for each edge
393 ((TH1D*)(*fHistEdg)[i])->Sumw2();
394 }
395 //====== extra histograms for debug purposes
396 fHistDbg = new TObjArray(fDim); // Initialize list of histograms
397 for(i=0;i<fDim;i++){
398 hname=fName+TString("_HistDebug_");
399 hname+=i;
400 htitle=TString("Debug Histogram ");
401 htitle+=i;
402 (*fHistDbg)[i] = new TH1D(hname.Data(),htitle.Data(),fNBin,0.0, 1.0); // Initialize histogram for each edge
403 }
404
405 // ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| //
406 // BUILD-UP of the FOAM //
407 // ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| //
408 //
409 // Define and explore root cell(s)
410 InitCells();
411 // PrintCells(); std::cout<<" ===== after InitCells ====="<<std::endl;
412 Grow();
413 // PrintCells(); std::cout<<" ===== after Grow ====="<<std::endl;
414
415 MakeActiveList(); // Final Preparations for the M.C. generation
416
417 // Preparations for the M.C. generation
418 fSumWt = 0.0; // M.C. generation sum of Wt
419 fSumWt2 = 0.0; // M.C. generation sum of Wt**2
420 fSumOve = 0.0; // M.C. generation sum of overweighted
421 fNevGen = 0.0; // M.C. generation sum of 1d0
422 fWtMax = gVlow; // M.C. generation maximum wt
423 fWtMin = gHigh; // M.C. generation minimum wt
424 fMCresult=fCells[0]->GetIntg(); // M.C. Value of INTEGRAL,temporary assignment
425 fMCresult=fCells[0]->GetIntg(); // M.C. Value of INTEGRAL,temporary assignment
426 fMCerror =fCells[0]->GetIntg(); // M.C. Value of ERROR ,temporary assignment
427 fMCMonit = new TFoamMaxwt(5.0,1000); // monitoring M.C. efficiency
428 //
429 if(fChat>0){
430 Double_t driver = fCells[0]->GetDriv();
431 BXOPE;
432 BXTXT("*** TFoam::Initialize FINISHED!!! ***");
433 BX1I(" nCalls",fNCalls, "Total number of function calls ");
434 BX1F(" XPrime",fPrime, "Primary total integral ");
435 BX1F(" XDiver",driver, "Driver total integral ");
436 BX1F(" mcResult",fMCresult,"Estimate of the true MC Integral ");
437 BXCLO;
438 }
439 if(fChat==2) PrintCells();
440 TH1::AddDirectory(addStatus);
441} // Initialize
442
443////////////////////////////////////////////////////////////////////////////////
444/// Internal method used by Initialize.
445/// It initializes "root part" of the FOAM of the tree of cells.
446
448{
449 Int_t i;
450
451 fLastCe =-1; // Index of the last cell
452 if(fCells!= 0) {
453 for(i=0; i<fNCells; i++) delete fCells[i];
454 delete [] fCells;
455 }
456 //
457 fCells = new TFoamCell*[fNCells];
458 for(i=0;i<fNCells;i++){
459 fCells[i]= new TFoamCell(fDim); // Allocate BIG list of cells
460 fCells[i]->SetSerial(i);
461 }
462 if(fCells==0) Error("InitCells", "Cannot initialize CELLS \n" );
463
464 /////////////////////////////////////////////////////////////////////////////
465 // Single Root Hypercube //
466 /////////////////////////////////////////////////////////////////////////////
467 CellFill(1, 0); // 0-th cell ACTIVE
468
469 // Exploration of the root cell(s)
470 for(Long_t iCell=0; iCell<=fLastCe; iCell++){
471 Explore( fCells[iCell] ); // Exploration of root cell(s)
472 }
473}
474
475////////////////////////////////////////////////////////////////////////////////
476/// Internal method used by Initialize.
477/// It initializes content of the newly allocated active cell.
478
480{
481 TFoamCell *cell;
482 if (fLastCe==fNCells){
483 Error( "CellFill", "Too many cells\n");
484 }
485 fLastCe++; // 0-th cell is the first
486 if (Status==1) fNoAct++;
487
488 cell = fCells[fLastCe];
489
490 cell->Fill(Status, parent, 0, 0);
491
492 cell->SetBest( -1); // pointer for planning division of the cell
493 cell->SetXdiv(0.5); // factor for division
494 Double_t xInt2,xDri2;
495 if(parent!=0){
496 xInt2 = 0.5*parent->GetIntg();
497 xDri2 = 0.5*parent->GetDriv();
498 cell->SetIntg(xInt2);
499 cell->SetDriv(xDri2);
500 }else{
501 cell->SetIntg(0.0);
502 cell->SetDriv(0.0);
503 }
504 return fLastCe;
505}
506
507////////////////////////////////////////////////////////////////////////////////
508/// Internal method used by Initialize.
509///
510/// It explores newly defined cell with help of special short MC sampling.
511/// As a result, estimates of true and drive volume is defined/determined
512/// Average and dispersion of the weight distribution will is found along
513/// each edge and the best edge (minimum dispersion, best maximum weight)
514/// is memorized for future use.
515///
516/// The optimal division point for eventual future cell division is
517/// determined/recorded. Recorded are also minimum and maximum weight etc.
518/// The volume estimate in all (inactive) parent cells is updated.
519/// Note that links to parents and initial volume = 1/2 parent has to be
520/// already defined prior to calling this routine.
521
523{
524 Double_t wt, dx, xBest=0, yBest=0;
525 Double_t intOld, driOld;
526
527 Long_t iev;
528 Double_t nevMC;
529 Int_t i, j, k;
530 Int_t nProj, kBest;
531 Double_t ceSum[5], xproj;
532
533 TFoamVect cellSize(fDim);
534 TFoamVect cellPosi(fDim);
535
536 cell->GetHcub(cellPosi,cellSize);
537
538 TFoamCell *parent;
539
540 Double_t *xRand = new Double_t[fDim];
541
542 Double_t *volPart=0;
543
544 cell->CalcVolume();
545 dx = cell->GetVolume();
546 intOld = cell->GetIntg(); //memorize old values,
547 driOld = cell->GetDriv(); //will be needed for correcting parent cells
548
549
550 /////////////////////////////////////////////////////
551 // Special Short MC sampling to probe cell //
552 /////////////////////////////////////////////////////
553 ceSum[0]=0;
554 ceSum[1]=0;
555 ceSum[2]=0;
556 ceSum[3]=gHigh; //wtmin
557 ceSum[4]=gVlow; //wtmax
558 //
559 for(i=0;i<fDim;i++) ((TH1D *)(*fHistEdg)[i])->Reset(); // Reset histograms
560 fHistWt->Reset();
561 //
562 // ||||||||||||||||||||||||||BEGIN MC LOOP|||||||||||||||||||||||||||||
563 Double_t nevEff=0.;
564 for(iev=0;iev<fNSampl;iev++){
565 MakeAlpha(); // generate uniformly vector inside hypercube
566
567 if(fDim>0){
568 for(j=0; j<fDim; j++)
569 xRand[j]= cellPosi[j] +fAlpha[j]*(cellSize[j]);
570 }
571
572 wt=dx*Eval(xRand);
573
574 nProj = 0;
575 if(fDim>0) {
576 for(k=0; k<fDim; k++) {
577 xproj =fAlpha[k];
578 ((TH1D *)(*fHistEdg)[nProj])->Fill(xproj,wt);
579 nProj++;
580 }
581 }
582 //
583 fNCalls++;
584 ceSum[0] += wt; // sum of weights
585 ceSum[1] += wt*wt; // sum of weights squared
586 ceSum[2]++; // sum of 1
587 if (ceSum[3]>wt) ceSum[3]=wt; // minimum weight;
588 if (ceSum[4]<wt) ceSum[4]=wt; // maximum weight
589 // test MC loop exit condition
590 nevEff = ceSum[1] == 0. ? 0. : ceSum[0]*ceSum[0]/ceSum[1];
591 if( nevEff >= fNBin*fEvPerBin) break;
592 } // ||||||||||||||||||||||||||END MC LOOP|||||||||||||||||||||||||||||
593 //------------------------------------------------------------------
594 //--- predefine logics of searching for the best division edge ---
595 for(k=0; k<fDim;k++){
596 fMaskDiv[k] =1; // default is all
597 if( fInhiDiv[k]==1) fMaskDiv[k] =0; // inhibit some...
598 }
599 // Note that predefined division below overrule inhibition above
600 kBest=-1;
601 Double_t rmin,rmax,rdiv;
602 if(fOptPRD) { // quick check
603 for(k=0; k<fDim; k++) {
604 rmin= cellPosi[k];
605 rmax= cellPosi[k] +cellSize[k];
606 if( fXdivPRD[k] != 0) {
607 Int_t n= (fXdivPRD[k])->GetDim();
608 for(j=0; j<n; j++) {
609 rdiv=(*fXdivPRD[k])[j];
610 // check predefined divisions is available in this cell
611 if( (rmin +1e-99 <rdiv) && (rdiv< rmax -1e-99)) {
612 kBest=k;
613 xBest= (rdiv-cellPosi[k])/cellSize[k] ;
614 goto ee05;
615 }
616 }
617 }
618 }//k
619 }
620 ee05:
621 /////////////////////////////////////////////////////////////////////////////
622
623 fNEffev += (Long_t)nevEff;
624 nevMC = ceSum[2];
625 Double_t intTrue = ceSum[0]/(nevMC+0.000001);
626 Double_t intDriv=0.;
627 Double_t intPrim=0.;
628
629 switch(fOptDrive){
630 case 1: // VARIANCE REDUCTION
631 if(kBest == -1) Varedu(ceSum,kBest,xBest,yBest); // determine the best edge,
632 //intDriv =sqrt( ceSum[1]/nevMC -intTrue*intTrue ); // Older ansatz, numerically not bad
633 intDriv =sqrt(ceSum[1]/nevMC) -intTrue; // Foam build-up, sqrt(<w**2>) -<w>
634 intPrim =sqrt(ceSum[1]/nevMC); // MC gen. sqrt(<w**2>) =sqrt(<w>**2 +sigma**2)
635 break;
636 case 2: // WTMAX REDUCTION
637 if(kBest == -1) Carver(kBest,xBest,yBest); // determine the best edge
638 intDriv =ceSum[4] -intTrue; // Foam build-up, wtmax-<w>
639 intPrim =ceSum[4]; // MC generation, wtmax!
640 break;
641 default:
642 Error("Explore", "Wrong fOptDrive = \n" );
643 }//switch
644 //=================================================================================
645 //hist_Neff_distrib.Fill( fLastCe/2.0+0.01, nevEff+0.01); //
646 //hist_kBest_distrib.Fill( kBest+0.50, 1.0 ); // debug
647 //hist_xBest_distrib.Fill( xBest+0.01, 1.0 ); // debug
648 //=================================================================================
649 cell->SetBest(kBest);
650 cell->SetXdiv(xBest);
651 cell->SetIntg(intTrue);
652 cell->SetDriv(intDriv);
653 cell->SetPrim(intPrim);
654 // correct/update integrals in all parent cells to the top of the tree
655 Double_t parIntg, parDriv;
656 for(parent = cell->GetPare(); parent!=0; parent = parent->GetPare()){
657 parIntg = parent->GetIntg();
658 parDriv = parent->GetDriv();
659 parent->SetIntg( parIntg +intTrue -intOld );
660 parent->SetDriv( parDriv +intDriv -driOld );
661 }
662 delete [] volPart;
663 delete [] xRand;
664 //cell->Print();
665} // TFoam::Explore
666
667////////////////////////////////////////////////////////////////////////////////
668/// Internal method used by Initialize.
669///
670/// It determines the best edge candidate and the position of the cell division plane
671/// in case of the variance reduction for future cell division,
672/// using results of the MC exploration run stored in fHistEdg
673
674void TFoam::Varedu(Double_t ceSum[5], Int_t &kBest, Double_t &xBest, Double_t &yBest)
675{
676 Double_t nent = ceSum[2];
677 Double_t swAll = ceSum[0];
678 Double_t sswAll = ceSum[1];
679 Double_t ssw = sqrt(sswAll)/sqrt(nent);
680 //
681 Double_t swIn,swOut,sswIn,sswOut,xLo,xUp;
682 kBest =-1;
683 xBest =0.5;
684 yBest =1.0;
685 Double_t maxGain=0.0;
686 // Now go over all projections kProj
687 for(Int_t kProj=0; kProj<fDim; kProj++) {
688 if( fMaskDiv[kProj]) {
689 // initialize search over bins
690 Double_t sigmIn =0.0; Double_t sigmOut =0.0;
691 Double_t sswtBest = gHigh;
692 Double_t gain =0.0;
693 Double_t xMin=0.0; Double_t xMax=0.0;
694 // Double loop over all pairs jLo<jUp
695 for(Int_t jLo=1; jLo<=fNBin; jLo++) {
696 Double_t aswIn=0; Double_t asswIn=0;
697 for(Int_t jUp=jLo; jUp<=fNBin;jUp++) {
698 aswIn += ((TH1D *)(*fHistEdg)[kProj])->GetBinContent(jUp);
699 asswIn += Sqr(((TH1D *)(*fHistEdg)[kProj])->GetBinError( jUp));
700 xLo=(jLo-1.0)/fNBin;
701 xUp=(jUp*1.0)/fNBin;
702 swIn = aswIn/nent;
703 swOut = (swAll-aswIn)/nent;
704 sswIn = sqrt(asswIn) /sqrt(nent*(xUp-xLo)) *(xUp-xLo);
705 sswOut= sqrt(sswAll-asswIn)/sqrt(nent*(1.0-xUp+xLo)) *(1.0-xUp+xLo);
706 if( (sswIn+sswOut) < sswtBest) {
707 sswtBest = sswIn+sswOut;
708 gain = ssw-sswtBest;
709 sigmIn = sswIn -swIn; // Debug
710 sigmOut = sswOut-swOut; // Debug
711 xMin = xLo;
712 xMax = xUp;
713 }
714 }//jUp
715 }//jLo
716 Int_t iLo = (Int_t) (fNBin*xMin);
717 Int_t iUp = (Int_t) (fNBin*xMax);
718 //----------DEBUG printout
719 //std::cout<<"@@@@@ xMin xMax = "<<xMin <<" "<<xMax<<" iLo= "<<iLo<<" iUp= "<<iUp;
720 //std::cout<<" sswtBest/ssw= "<<sswtBest/ssw<<" Gain/ssw= "<< Gain/ssw<<std::endl;
721 //----------DEBUG auxiliary Plot
722 for(Int_t iBin=1;iBin<=fNBin;iBin++) {
723 if( ((iBin-0.5)/fNBin > xMin) && ((iBin-0.5)/fNBin < xMax) ){
724 ((TH1D *)(*fHistDbg)[kProj])->SetBinContent(iBin,sigmIn/(xMax-xMin));
725 } else {
726 ((TH1D *)(*fHistDbg)[kProj])->SetBinContent(iBin,sigmOut/(1-xMax+xMin));
727 }
728 }
729 if(gain>=maxGain) {
730 maxGain=gain;
731 kBest=kProj; // <--- !!!!! The best edge
732 xBest=xMin;
733 yBest=xMax;
734 if(iLo == 0) xBest=yBest; // The best division point
735 if(iUp == fNBin) yBest=xBest; // this is not really used
736 }
737 }
738 } //kProj
739 //----------DEBUG printout
740 //std::cout<<"@@@@@@@>>>>> kBest= "<<kBest<<" maxGain/ssw= "<< maxGain/ssw<<std::endl;
741 if( (kBest >= fDim) || (kBest<0) ) Error("Varedu", "Something wrong with kBest - kBest = %d dim = %d\n",kBest,fDim);
742} //TFoam::Varedu
743
744////////////////////////////////////////////////////////////////////////////////
745/// Internal method used by Initialize.
746///
747/// Determines the best edge-candidate and the position of the division plane
748/// for the future cell division, in the case of the optimization of the maximum weight.
749/// It exploits results of the cell MC exploration run stored in fHistEdg.
750
751void TFoam::Carver(Int_t &kBest, Double_t &xBest, Double_t &yBest)
752{
753 Int_t kProj,iBin;
754 Double_t carve,carvTot,carvMax,carvOne,binMax;
755 Int_t jLow,jUp,iLow,iUp;
756 Double_t theBin;
757 // Int_t jDivi; // TEST
758 Int_t j;
759
760 Double_t *bins = new Double_t[fNBin]; // bins of histogram for single PROJECTION
761 if(bins==0) Error("Carver", "Cannot initialize buffer Bins \n" );
762
763 kBest =-1;
764 xBest =0.5;
765 yBest =1.0;
766 carvMax = gVlow;
767 // primMax = gVlow;
768 for(kProj=0; kProj<fDim; kProj++)
769 if( fMaskDiv[kProj] ) {
770 //if( kProj==1 ){
771 //std::cout<<"==================== Carver histogram: kProj ="<<kProj<<"==================="<<std::endl;
772 //((TH1D *)(*fHistEdg)[kProj])->Print("all");
773 binMax = gVlow;
774 for(iBin=0; iBin<fNBin;iBin++){
775 bins[iBin]= ((TH1D *)(*fHistEdg)[kProj])->GetBinContent(iBin+1);
776 binMax = TMath::Max( binMax, bins[iBin]); // Maximum content/bin
777 }
778 if(binMax < 0 ) { //case of empty cell
779 delete [] bins;
780 return;
781 }
782 carvTot = 0.0;
783 for(iBin=0;iBin<fNBin;iBin++){
784 carvTot = carvTot + (binMax-bins[iBin]); // Total Carve (more stable)
785 }
786 // primTot = binMax*fNBin;
787 //std::cout <<"Carver: CarvTot "<<CarvTot<< " primTot "<<primTot<<std::endl;
788 jLow =0;
789 jUp =fNBin-1;
790 carvOne = gVlow;
791 Double_t yLevel = gVlow;
792 for(iBin=0; iBin<fNBin;iBin++) {
793 theBin = bins[iBin];
794 //----- walk to the left and find first bin > theBin
795 iLow = iBin;
796 for(j=iBin; j>-1; j-- ) {
797 if(theBin< bins[j]) break;
798 iLow = j;
799 }
800 //iLow = iBin;
801 //if(iLow>0) while( (theBin >= bins[iLow-1])&&(iLow >0) ){iLow--;} // horror!!!
802 //------ walk to the right and find first bin > theBin
803 iUp = iBin;
804 for(j=iBin; j<fNBin; j++) {
805 if(theBin< bins[j]) break;
806 iUp = j;
807 }
808 //iUp = iBin;
809 //if(iUp<fNBin-1) while( (theBin >= bins[iUp+1])&&( iUp<fNBin-1 ) ){iUp++;} // horror!!!
810 //
811 carve = (iUp-iLow+1)*(binMax-theBin);
812 if( carve > carvOne) {
813 carvOne = carve;
814 jLow = iLow;
815 jUp = iUp;
816 yLevel = theBin;
817 }
818 }//iBin
819 if( carvTot > carvMax) {
820 carvMax = carvTot;
821 //primMax = primTot;
822 //std::cout <<"Carver: primMax "<<primMax<<std::endl;
823 kBest = kProj; // Best edge
824 xBest = ((Double_t)(jLow))/fNBin;
825 yBest = ((Double_t)(jUp+1))/fNBin;
826 if(jLow == 0 ) xBest = yBest;
827 if(jUp == fNBin-1) yBest = xBest;
828 // division ratio in units of 1/fNBin, testing
829 // jDivi = jLow;
830 // if(jLow == 0 ) jDivi=jUp+1;
831 }
832 //====== extra histograms for debug purposes
833 //std::cout<<"kProj= "<<kProj<<" jLow= "<<jLow<<" jUp= "<<jUp<<std::endl;
834 for(iBin=0; iBin<fNBin; iBin++)
835 ((TH1D *)(*fHistDbg)[kProj])->SetBinContent(iBin+1,binMax);
836 for(iBin=jLow; iBin<jUp+1; iBin++)
837 ((TH1D *)(*fHistDbg)[kProj])->SetBinContent(iBin+1,yLevel);
838 }//kProj
839 if( (kBest >= fDim) || (kBest<0) ) Error("Carver", "Something wrong with kBest - kBest = %d dim = %d\n",kBest,fDim);
840 delete [] bins;
841} //TFoam::Carver
842
843////////////////////////////////////////////////////////////////////////////////
844/// Internal method used by Initialize.
845/// Provides random vector Alpha 0< Alpha(i) < 1
846
848{
849 Int_t k;
850 if(fDim<1) return;
851
852 // simply generate and load kDim uniform random numbers
853 fPseRan->RndmArray(fDim,fRvec); // kDim random numbers needed
854 for(k=0; k<fDim; k++) fAlpha[k] = fRvec[k];
855} //MakeAlpha
856
857
858////////////////////////////////////////////////////////////////////////////////
859/// Internal method used by Initialize.
860/// It grow new cells by the binary division process.
861
863{
864 Long_t iCell;
865 TFoamCell* newCell;
866
867 while ( (fLastCe+2) < fNCells ) { // this condition also checked inside Divide
868 iCell = PeekMax(); // peek up cell with maximum driver integral
869 if( (iCell<0) || (iCell>fLastCe) ) {
870 Error("Grow", "Wrong iCell \n");
871 break;
872 }
873 newCell = fCells[iCell];
874
875 if(fLastCe !=0) {
876 Int_t kEcho=10;
877 if(fLastCe>=10000) kEcho=100;
878 if( (fLastCe%kEcho)==0) {
879 if (fChat>0) {
880 if(fDim<10)
881 std::cout<<fDim<<std::flush;
882 else
883 std::cout<<"."<<std::flush;
884 if( (fLastCe%(100*kEcho))==0) std::cout<<"|"<<fLastCe<<std::endl<<std::flush;
885 }
886 }
887 }
888 if( Divide( newCell )==0) break; // and divide it into two
889 }
890 if (fChat>0) {
891 std::cout<<std::endl<<std::flush;
892 }
893 CheckAll(0); // set arg=1 for more info
894}
895
896////////////////////////////////////////////////////////////////////////////////
897/// Internal method used by Initialize.
898/// It finds cell with maximal driver integral for the purpose of the division.
899
901{
902 Long_t i;
903 Long_t iCell = -1;
904 Double_t drivMax, driv;
905
906 drivMax = gVlow;
907 for(i=0; i<=fLastCe; i++) {//without root
908 if( fCells[i]->GetStat() == 1 ) {
909 driv = TMath::Abs( fCells[i]->GetDriv());
910 //std::cout<<"PeekMax: Driv = "<<driv<<std::endl;
911 if(driv > drivMax) {
912 drivMax = driv;
913 iCell = i;
914 }
915 }
916 }
917 // std::cout << 'TFoam_PeekMax: iCell=' << iCell << std::endl;
918 if (iCell == -1)
919 std::cout << "STOP in TFoam::PeekMax: not found iCell=" << iCell << std::endl;
920 return(iCell);
921} // TFoam_PeekMax
922
923////////////////////////////////////////////////////////////////////////////////
924/// Internal method used by Initialize.
925///
926/// It divides cell iCell into two daughter cells.
927/// The iCell is retained and tagged as inactive, daughter cells are appended
928/// at the end of the buffer.
929/// New vertex is added to list of vertices.
930/// List of active cells is updated, iCell removed, two daughters added
931/// and their properties set with help of MC sampling (TFoam_Explore)
932/// Returns Code RC=-1 of buffer limit is reached, fLastCe=fnBuf.
933
935{
936 Int_t kBest;
937
938 if(fLastCe+1 >= fNCells) Error("Divide", "Buffer limit is reached, fLastCe=fnBuf \n");
939
940 cell->SetStat(0); // reset cell as inactive
941 fNoAct--;
942
943 kBest = cell->GetBest();
944 if( kBest<0 || kBest>=fDim ) Error("Divide", "Wrong kBest \n");
945
946 //////////////////////////////////////////////////////////////////
947 // define two daughter cells (active) //
948 //////////////////////////////////////////////////////////////////
949
950 Int_t d1 = CellFill(1, cell);
951 Int_t d2 = CellFill(1, cell);
952 cell->SetDau0((fCells[d1]));
953 cell->SetDau1((fCells[d2]));
954 Explore( (fCells[d1]) );
955 Explore( (fCells[d2]) );
956 return 1;
957} // TFoam_Divide
958
959
960////////////////////////////////////////////////////////////////////////////////
961/// Internal method used by Initialize.
962///
963/// It finds out number of active cells fNoAct,
964/// creates list of active cell fCellsAct and primary cumulative fPrimAcu.
965/// They are used during the MC generation to choose randomly an active cell.
966
968{
969 Long_t iCell;
971 // flush previous result
972 if(fPrimAcu != 0) delete [] fPrimAcu;
973 fCellsAct.clear();
974 fCellsAct.reserve(fNoAct);
975
976 // Count Active cells and find total Primary
977 // Fill-in tables of active cells
978
979 fPrime = 0.0;
980 for(iCell=0; iCell<=fLastCe; iCell++) {
981 if (fCells[iCell]->GetStat()==1) {
982 fPrime += fCells[iCell]->GetPrim();
983 fCellsAct.push_back(iCell);
984 }
985 }
986
987 if(fNoAct != static_cast<Int_t>(fCellsAct.size())) Error("MakeActiveList", "Wrong fNoAct \n" );
988 if(fPrime == 0.) Error("MakeActiveList", "Integrand function is zero \n" );
989
990 fPrimAcu = new Double_t[fNoAct]; // cumulative primary for MC generation
991 if( fPrimAcu==0 ) Error("MakeActiveList", "Cant allocate fCellsAct or fPrimAcu \n");
992
993 sum =0.0;
994 for(iCell=0; iCell<fNoAct; iCell++) {
995 sum = sum + ( fCells[fCellsAct[iCell]] )->GetPrim()/fPrime;
996 fPrimAcu[iCell]=sum;
997 }
998
999} //MakeActiveList
1000
1001////////////////////////////////////////////////////////////////////////////////
1002/// User may optionally reset random number generator using this method.
1003///
1004/// Usually it is done when FOAM object is restored from the disk.
1005/// IMPORTANT: this method deletes existing random number generator registered in the FOAM object.
1006/// In particular such an object is created by the streamer during the disk-read operation.
1007
1009{
1010 if(fPseRan) {
1011 Info("ResetPseRan", "Resetting random number generator \n");
1012 delete fPseRan;
1013 }
1014 SetPseRan(PseRan);
1015}
1016
1017////////////////////////////////////////////////////////////////////////////////
1018/// User may use this method to set the distribution object
1019
1021{
1022 if (fun)
1023 fRho=fun;
1024 else
1025 Error("SetRho", "Bad function \n" );
1026}
1027////////////////////////////////////////////////////////////////////////////////
1028/// User may use this method to set the distribution object as a global function pointer
1029/// (and not as an interpreted function).
1030
1032{
1033 // This is needed for both AClic and Cling
1034 if (fun) {
1035 // delete function object if it has been created here in SetRho
1036 if (fRho && dynamic_cast<FoamIntegrandFunction*>(fRho) ) delete fRho;
1037 fRho= new FoamIntegrandFunction(fun);
1038 } else
1039 Error("SetRho", "Bad function \n" );
1040}
1041
1042////////////////////////////////////////////////////////////////////////////////
1043/// User may optionally reset the distribution using this method.
1044///
1045/// Usually it is done when FOAM object is restored from the disk.
1046/// IMPORTANT: this method deletes existing distribution object registered in the FOAM object.
1047/// In particular such an object is created by the streamer diring the disk-read operation.
1048/// This method is used only in very special cases, because the distribution in most cases
1049/// should be "owned" by the FOAM object and should not be replaced by another one after initialization.
1050
1052{
1053 if(fRho) {
1054 Info("ResetRho", "!!! Resetting distribution function !!!\n");
1055 delete fRho;
1056 }
1057 SetRho(fun);
1058}
1059
1060////////////////////////////////////////////////////////////////////////////////
1061/// Internal method.
1062/// Evaluates distribution to be generated.
1063
1065{
1067
1068 if(!fRho) { //interactive mode
1069 Longptr_t paramArr[3];
1070 paramArr[0]=(Longptr_t)fDim;
1071 paramArr[1]=(Longptr_t)xRand;
1072 fMethodCall->SetParamPtrs(paramArr);
1074 } else { //compiled mode
1075 result=fRho->Density(fDim,xRand);
1076 }
1077
1078 return result;
1079}
1080
1081////////////////////////////////////////////////////////////////////////////////
1082/// Internal method.
1083/// Return randomly chosen active cell with probability equal to its
1084/// contribution into total driver integral using interpolation search.
1085
1087{
1088 Long_t lo, hi, hit;
1089 Double_t fhit, flo, fhi;
1090 Double_t random;
1091
1092 random=fPseRan->Rndm();
1093 lo = 0; hi =fNoAct-1;
1094 flo = fPrimAcu[lo]; fhi=fPrimAcu[hi];
1095 while(lo+1<hi) {
1096 hit = lo + (Int_t)( (hi-lo)*(random-flo)/(fhi-flo)+0.5);
1097 if (hit<=lo)
1098 hit = lo+1;
1099 else if(hit>=hi)
1100 hit = hi-1;
1101 fhit=fPrimAcu[hit];
1102 if (fhit>random) {
1103 hi = hit;
1104 fhi = fhit;
1105 } else {
1106 lo = hit;
1107 flo = fhit;
1108 }
1109 }
1110 if (fPrimAcu[lo]>random)
1111 pCell = fCells[fCellsAct[lo]];
1112 else
1113 pCell = fCells[fCellsAct[hi]];
1114} // TFoam::GenerCel2
1115
1116
1117////////////////////////////////////////////////////////////////////////////////
1118/// User method.
1119/// It generates randomly point/vector according to user-defined distribution.
1120/// Prior initialization with help of Initialize() is mandatory.
1121/// Generated MC point/vector is available using GetMCvect and the MC weight with GetMCwt.
1122/// MC point is generated with wt=1 or with variable weight, see OptRej switch.
1123
1125{
1126 Int_t j;
1127 Double_t wt,dx,mcwt;
1128 TFoamCell *rCell;
1129 //
1130 //********************** MC LOOP STARS HERE **********************
1131ee0:
1132 GenerCel2(rCell); // choose randomly one cell
1133
1134 MakeAlpha();
1135
1136 TFoamVect cellPosi(fDim); TFoamVect cellSize(fDim);
1137 rCell->GetHcub(cellPosi,cellSize);
1138 for(j=0; j<fDim; j++)
1139 fMCvect[j]= cellPosi[j] +fAlpha[j]*cellSize[j];
1140 dx = rCell->GetVolume(); // Cartesian volume of the Cell
1141 // weight average normalized to PRIMARY integral over the cell
1142
1143 wt=dx*Eval(fMCvect);
1144
1145 mcwt = wt / rCell->GetPrim(); // PRIMARY controls normalization
1146 fNCalls++;
1147 fMCwt = mcwt;
1148 // accumulation of statistics for the main MC weight
1149 fSumWt += mcwt; // sum of Wt
1150 fSumWt2 += mcwt*mcwt; // sum of Wt**2
1151 fNevGen++; // sum of 1d0
1152 fWtMax = TMath::Max(fWtMax, mcwt); // maximum wt
1153 fWtMin = TMath::Min(fWtMin, mcwt); // minimum wt
1154 fMCMonit->Fill(mcwt);
1155 fHistWt->Fill(mcwt,1.0); // histogram
1156 //******* Optional rejection ******
1157 if(fOptRej == 1) {
1158 Double_t random;
1159 random=fPseRan->Rndm();
1160 if( fMaxWtRej*random > fMCwt) goto ee0; // Wt=1 events, internal rejection
1161 if( fMCwt<fMaxWtRej ) {
1162 fMCwt = 1.0; // normal Wt=1 event
1163 } else {
1164 fMCwt = fMCwt/fMaxWtRej; // weight for overweighted events! kept for debug
1165 fSumOve += fMCwt-fMaxWtRej; // contribution of overweighted
1166 }
1167 }
1168 //********************** MC LOOP ENDS HERE **********************
1169} // MakeEvent
1170
1171////////////////////////////////////////////////////////////////////////////////
1172/// User may get generated MC point/vector with help of this method
1173
1175{
1176 for ( Int_t k=0 ; k<fDim ; k++) *(MCvect +k) = fMCvect[k];
1177}
1178
1179////////////////////////////////////////////////////////////////////////////////
1180/// User may get weight MC weight using this method
1181
1183{
1184 return(fMCwt);
1185}
1186////////////////////////////////////////////////////////////////////////////////
1187/// User may get weight MC weight using this method
1188
1190{
1191 mcwt=fMCwt;
1192}
1193
1194////////////////////////////////////////////////////////////////////////////////
1195/// User method which generates MC event and returns MC weight
1196
1198{
1199 MakeEvent();
1200 GetMCvect(MCvect);
1201 return(fMCwt);
1202}
1203
1204////////////////////////////////////////////////////////////////////////////////
1205/// User method.
1206/// It provides the value of the integral calculated from the averages of the MC run
1207/// May be called after (or during) the MC run.
1208
1209void TFoam::GetIntegMC(Double_t &mcResult, Double_t &mcError)
1210{
1211 Double_t mCerelat;
1212 mcResult = 0.0;
1213 mCerelat = 1.0;
1214 if (fNevGen>0) {
1215 mcResult = fPrime*fSumWt/fNevGen;
1216 mCerelat = sqrt( fSumWt2/(fSumWt*fSumWt) - 1/fNevGen);
1217 }
1218 mcError = mcResult *mCerelat;
1219}
1220
1221////////////////////////////////////////////////////////////////////////////////
1222/// User method.
1223/// It returns NORMALIZATION integral to be combined with the average weights
1224/// and content of the histograms in order to get proper absolute normalization
1225/// of the integrand and distributions.
1226/// It can be called after initialization, before or during the MC run.
1227
1228void TFoam::GetIntNorm(Double_t& IntNorm, Double_t& Errel )
1229{
1230 if(fOptRej == 1) { // Wt=1 events, internal rejection
1231 Double_t intMC,errMC;
1232 GetIntegMC(intMC,errMC);
1233 IntNorm = intMC;
1234 Errel = errMC;
1235 } else { // Wted events, NO internal rejection
1236 IntNorm = fPrime;
1237 Errel = 0;
1238 }
1239}
1240
1241////////////////////////////////////////////////////////////////////////////////
1242/// May be called optionally after the MC run.
1243/// Returns various parameters of the MC weight for efficiency evaluation
1244
1246{
1247 Double_t mCeff, wtLim;
1248 fMCMonit->GetMCeff(eps, mCeff, wtLim);
1249 wtMax = wtLim;
1250 aveWt = fSumWt/fNevGen;
1251 sigma = sqrt( fSumWt2/fNevGen -aveWt*aveWt );
1252}
1253
1254////////////////////////////////////////////////////////////////////////////////
1255/// May be called optionally by the user after the MC run.
1256/// It provides normalization and also prints some information/statistics on the MC run.
1257
1258void TFoam::Finalize(Double_t& IntNorm, Double_t& Errel)
1259{
1260 GetIntNorm(IntNorm,Errel);
1261 Double_t mcResult,mcError;
1262 GetIntegMC(mcResult,mcError);
1263 Double_t mCerelat= mcError/mcResult;
1264 //
1265 if(fChat>0) {
1266 Double_t eps = 0.0005;
1267 Double_t mCeff, mcEf2, wtMax, aveWt, sigma;
1268 GetWtParams(eps, aveWt, wtMax, sigma);
1269 mCeff=0;
1270 if(wtMax>0.0) mCeff=aveWt/wtMax;
1271 mcEf2 = sigma/aveWt;
1272 Double_t driver = fCells[0]->GetDriv();
1273 //
1274 BXOPE;
1275 BXTXT("****************************************");
1276 BXTXT("****** TFoam::Finalize ******");
1277 BXTXT("****************************************");
1278 BX1I(" NevGen",fNevGen, "Number of generated events in the MC generation ");
1279 BX1I(" nCalls",fNCalls, "Total number of function calls ");
1280 BXTXT("----------------------------------------");
1281 BX1F(" AveWt",aveWt, "Average MC weight ");
1282 BX1F(" WtMin",fWtMin, "Minimum MC weight (absolute) ");
1283 BX1F(" WtMax",fWtMax, "Maximum MC weight (absolute) ");
1284 BXTXT("----------------------------------------");
1285 BX1F(" XPrime",fPrime, "Primary total integral, R_prime ");
1286 BX1F(" XDiver",driver, "Driver total integral, R_loss ");
1287 BXTXT("----------------------------------------");
1288 BX2F(" IntMC", mcResult, mcError, "Result of the MC Integral");
1289 BX1F(" mCerelat", mCerelat, "Relative error of the MC integral ");
1290 BX1F(" <w>/WtMax",mCeff, "MC efficiency, acceptance rate");
1291 BX1F(" Sigma/<w>",mcEf2, "MC efficiency, variance/ave_wt");
1292 BX1F(" WtMax",wtMax, "WtMax(esp= 0.0005) ");
1293 BX1F(" Sigma",sigma, "variance of MC weight ");
1294 if(fOptRej==1) {
1295 Double_t avOve=fSumOve/fSumWt;
1296 BX1F("<OveW>/<W>",avOve, "Contrib. of events wt>MaxWtRej");
1297 }
1298 BXCLO;
1299 }
1300} // Finalize
1301
1302////////////////////////////////////////////////////////////////////////////////
1303/// This can be called before Initialize, after setting kDim
1304/// It defines which variables are excluded in the process of the cell division.
1305/// For example 'FoamX->SetInhiDiv(1, 1);' inhibits division of y-variable.
1306/// The resulting map of cells in 2-dim. case will look as follows:
1307///
1308/// \image html foam_Map2.png width=400
1309
1310void TFoam::SetInhiDiv(Int_t iDim, Int_t InhiDiv)
1311{
1312
1313
1314 if(fDim==0) Error("TFoam","SetInhiDiv: fDim=0 \n");
1315 if(fInhiDiv == 0) {
1316 fInhiDiv = new Int_t[ fDim ];
1317 for(Int_t i=0; i<fDim; i++) fInhiDiv[i]=0;
1318 }
1319 //
1320 if( ( 0<=iDim) && (iDim<fDim)) {
1321 fInhiDiv[iDim] = InhiDiv;
1322 } else
1323 Error("SetInhiDiv:","Wrong iDim \n");
1324}
1325
1326////////////////////////////////////////////////////////////////////////////////
1327/// This should be called before Initialize, after setting kDim
1328/// It predefines values of the cell division for certain variable iDim.
1329/// For example setting 3 predefined division lines using:
1330///
1331/// ~~~ {.cpp}
1332/// xDiv[0]=0.30; xDiv[1]=0.40; xDiv[2]=0.65;
1333/// FoamX->SetXdivPRD(0,3,xDiv);
1334/// ~~~
1335///
1336/// results in the following 2-dim. pattern of the cells:
1337///
1338/// \image html foam_Map3.png width=400
1339
1341{
1342 Int_t i;
1343
1344 if(fDim<=0) Error("SetXdivPRD", "fDim=0 \n");
1345 if( len<1 ) Error("SetXdivPRD", "len<1 \n");
1346 // allocate list of pointers, if it was not done before
1347 if(fXdivPRD == 0) {
1348 fXdivPRD = new TFoamVect*[fDim];
1349 for(i=0; i<fDim; i++) fXdivPRD[i]=0;
1350 }
1351 // set division list for direction iDim in H-cubic space!!!
1352 if( ( 0<=iDim) && (iDim<fDim)) {
1353 fOptPRD =1; // !!!!
1354 if( fXdivPRD[iDim] != 0)
1355 Error("SetXdivPRD", "Second allocation of XdivPRD not allowed \n");
1356 fXdivPRD[iDim] = new TFoamVect(len); // allocate list of division points
1357 for(i=0; i<len; i++) {
1358 (*fXdivPRD[iDim])[i]=xDiv[i]; // set list of division points
1359 }
1360 } else {
1361 Error("SetXdivPRD", "Wrong iDim \n");
1362 }
1363 // Printing predefined division points
1364 std::cout<<" SetXdivPRD, idim= "<<iDim<<" len= "<<len<<" "<<std::endl;
1365 for(i=0; i<len; i++) {
1366 if (iDim < fDim) std::cout<< (*fXdivPRD[iDim])[i] <<" ";
1367 }
1368 std::cout<<std::endl;
1369 for(i=0; i<len; i++) std::cout<< xDiv[i] <<" ";
1370 std::cout<<std::endl;
1371 //
1372}
1373
1374////////////////////////////////////////////////////////////////////////////////
1375/// User utility, miscellaneous and debug.
1376/// Checks all pointers in the tree of cells. This is useful auto diagnostic.
1377/// - level=0, no printout, failures causes STOP
1378/// - level=1, printout, failures lead to WARNINGS only
1379
1381{
1382 Int_t errors, warnings;
1383 TFoamCell *cell;
1384 Long_t iCell;
1385
1386 errors = 0; warnings = 0;
1387 if (level==1) std::cout << "///////////////////////////// FOAM_Checks /////////////////////////////////" << std::endl;
1388 for(iCell=1; iCell<=fLastCe; iCell++) {
1389 cell = fCells[iCell];
1390 // checking general rules
1391 if( ((cell->GetDau0()==0) && (cell->GetDau1()!=0) ) ||
1392 ((cell->GetDau1()==0) && (cell->GetDau0()!=0) ) ) {
1393 errors++;
1394 if (level==1) Error("CheckAll","ERROR: Cell's no %ld has only one daughter \n",iCell);
1395 }
1396 if( (cell->GetDau0()==0) && (cell->GetDau1()==0) && (cell->GetStat()==0) ) {
1397 errors++;
1398 if (level==1) Error("CheckAll","ERROR: Cell's no %ld has no daughter and is inactive \n",iCell);
1399 }
1400 if( (cell->GetDau0()!=0) && (cell->GetDau1()!=0) && (cell->GetStat()==1) ) {
1401 errors++;
1402 if (level==1) Error("CheckAll","ERROR: Cell's no %ld has two daughters and is active \n",iCell);
1403 }
1404
1405 // checking parents
1406 if( (cell->GetPare())!=fCells[0] ) { // not child of the root
1407 if ( (cell != cell->GetPare()->GetDau0()) && (cell != cell->GetPare()->GetDau1()) ) {
1408 errors++;
1409 if (level==1) Error("CheckAll","ERROR: Cell's no %ld parent not pointing to this cell\n ",iCell);
1410 }
1411 }
1412
1413 // checking daughters
1414 if(cell->GetDau0()!=0) {
1415 if(cell != (cell->GetDau0())->GetPare()) {
1416 errors++;
1417 if (level==1) Error("CheckAll","ERROR: Cell's no %ld daughter 0 not pointing to this cell \n",iCell);
1418 }
1419 }
1420 if(cell->GetDau1()!=0) {
1421 if(cell != (cell->GetDau1())->GetPare()) {
1422 errors++;
1423 if (level==1) Error("CheckAll","ERROR: Cell's no %ld daughter 1 not pointing to this cell \n",iCell);
1424 }
1425 }
1426 }// loop after cells;
1427
1428 // Check for empty cells
1429 for(iCell=0; iCell<=fLastCe; iCell++) {
1430 cell = fCells[iCell];
1431 if( (cell->GetStat()==1) && (cell->GetDriv()==0) ) {
1432 warnings++;
1433 if(level==1) Warning("CheckAll", "Warning: Cell no. %ld is active but empty \n", iCell);
1434 }
1435 }
1436 // summary
1437 if(level==1){
1438 Info("CheckAll","Check has found %d errors and %d warnings \n",errors, warnings);
1439 }
1440 if(errors>0){
1441 Info("CheckAll","Check - found total %d errors \n",errors);
1442 }
1443} // Check
1444
1445////////////////////////////////////////////////////////////////////////////////
1446/// Prints geometry of ALL cells of the FOAM
1447
1449{
1450 Long_t iCell;
1451
1452 for(iCell=0; iCell<=fLastCe; iCell++) {
1453 std::cout<<"Cell["<<iCell<<"]={ ";
1454 //std::cout<<" "<< fCells[iCell]<<" "; // extra DEBUG
1455 std::cout<<std::endl;
1456 fCells[iCell]->Print("1");
1457 std::cout<<"}"<<std::endl;
1458 }
1459}
1460
1461////////////////////////////////////////////////////////////////////////////////
1462/// Debugging tool which plots 2-dimensional cells as rectangles
1463/// in C++ format readable for root
1464
1466{
1467 std::ofstream outfile(filename, std::ios::out);
1468 Double_t x1,y1,x2,y2,x,y;
1469 Long_t iCell;
1470 Double_t offs =0.1;
1471 Double_t lpag =1-2*offs;
1472 outfile<<"{" << std::endl;
1473 outfile<<"cMap = new TCanvas(\"Map1\",\" Cell Map \",600,600);"<<std::endl;
1474 //
1475 outfile<<"TBox*a=new TBox();"<<std::endl;
1476 outfile<<"a->SetFillStyle(0);"<<std::endl; // big frame
1477 outfile<<"a->SetLineWidth(4);"<<std::endl;
1478 outfile<<"a->SetLineColor(2);"<<std::endl;
1479 outfile<<"a->DrawBox("<<offs<<","<<offs<<","<<(offs+lpag)<<","<<(offs+lpag)<<");"<<std::endl;
1480 //
1481 outfile<<"TText*t=new TText();"<<std::endl; // text for numbering
1482 outfile<<"t->SetTextColor(4);"<<std::endl;
1483 if(fLastCe<51)
1484 outfile<<"t->SetTextSize(0.025);"<<std::endl; // text for numbering
1485 else if(fLastCe<251)
1486 outfile<<"t->SetTextSize(0.015);"<<std::endl;
1487 else
1488 outfile<<"t->SetTextSize(0.008);"<<std::endl;
1489 //
1490 outfile<<"TBox*b=new TBox();"<<std::endl; // single cell
1491 outfile <<"b->SetFillStyle(0);"<<std::endl;
1492 //
1493 if(fDim==2 && fLastCe<=2000) {
1494 TFoamVect cellPosi(fDim); TFoamVect cellSize(fDim);
1495 outfile << "// =========== Rectangular cells ==========="<< std::endl;
1496 for(iCell=1; iCell<=fLastCe; iCell++) {
1497 if( fCells[iCell]->GetStat() == 1) {
1498 fCells[iCell]->GetHcub(cellPosi,cellSize);
1499 x1 = offs+lpag*( cellPosi[0]); y1 = offs+lpag*( cellPosi[1]);
1500 x2 = offs+lpag*(cellPosi[0]+cellSize[0]); y2 = offs+lpag*(cellPosi[1]+cellSize[1]);
1501 // cell rectangle
1502 if(fLastCe<=2000)
1503 outfile<<"b->DrawBox("<<x1<<","<<y1<<","<<x2<<","<<y2<<");"<<std::endl;
1504 // cell number
1505 if(fLastCe<=250) {
1506 x = offs+lpag*(cellPosi[0]+0.5*cellSize[0]); y = offs+lpag*(cellPosi[1]+0.5*cellSize[1]);
1507 outfile<<"t->DrawText("<<x<<","<<y<<","<<"\""<<iCell<<"\""<<");"<<std::endl;
1508 }
1509 }
1510 }
1511 outfile<<"// ============== End Rectangles ==========="<< std::endl;
1512 }//kDim=2
1513 //
1514 //
1515 outfile << "}" << std::endl;
1516 outfile.close();
1517}
1518
1520{
1521// Void function for backward compatibility
1522
1523 Info("LinkCells", "VOID function for backward compatibility \n");
1524 return;
1525}
1526
#define e(i)
Definition: RSha256.hxx:103
int Int_t
Definition: RtypesCore.h:45
long Longptr_t
Definition: RtypesCore.h:82
char Char_t
Definition: RtypesCore.h:37
const Bool_t kFALSE
Definition: RtypesCore.h:101
long Long_t
Definition: RtypesCore.h:54
double Double_t
Definition: RtypesCore.h:59
#define ClassImp(name)
Definition: Rtypes.h:375
#define BX1F(name, numb, text)
Definition: TFoam.cxx:105
#define BX1I(name, numb, text)
Definition: TFoam.cxx:103
#define BXCLO
Definition: TFoam.cxx:110
#define BXTXT(text)
Definition: TFoam.cxx:101
#define BX2F(name, numb, err, text)
Definition: TFoam.cxx:107
static const Double_t gVlow
Definition: TFoam.cxx:116
static const Double_t gHigh
Definition: TFoam.cxx:115
#define BXOPE
Definition: TFoam.cxx:98
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t Float_t Float_t Int_t Int_t UInt_t UInt_t Rectangle_t Int_t Int_t Window_t TString Int_t GCValues_t GetPrimarySelectionOwner GetDisplay GetScreen GetColormap GetNativeEvent const char const char dpyName wid window const char font_name cursor keysym reg const char only_if_exist regb h Point_t winding char text const char depth char const char Int_t count const char ColorStruct_t color const char filename
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t Float_t Float_t Int_t Int_t UInt_t UInt_t Rectangle_t result
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t Float_t Float_t Int_t Int_t UInt_t UInt_t Rectangle_t Int_t Int_t Window_t TString Int_t GCValues_t GetPrimarySelectionOwner GetDisplay GetScreen GetColormap GetNativeEvent const char const char dpyName wid window const char font_name cursor keysym reg const char only_if_exist regb h Point_t winding char text const char depth char const char Int_t count const char ColorStruct_t color const char Pixmap_t Pixmap_t PictureAttributes_t attr const char char ret_data h unsigned char height h Atom_t Int_t ULong_t ULong_t unsigned char prop_list Atom_t Atom_t Atom_t Time_t UChar_t len
Option_t Option_t TPoint TPoint const char x2
Option_t Option_t TPoint TPoint const char x1
Option_t Option_t TPoint TPoint const char y2
Option_t Option_t TPoint TPoint const char y1
#define hi
Definition: THbookFile.cxx:128
Double_t Density(Int_t nDim, Double_t *x) override
Definition: TFoam.cxx:132
~FoamIntegrandFunction() override
Definition: TFoam.cxx:129
FoamIntegrandFunction(FunctionPtr func)
Definition: TFoam.cxx:127
FunctionPtr fFunc
Definition: TFoam.cxx:138
Double_t(* FunctionPtr)(Int_t, Double_t *)
Definition: TFoam.cxx:125
Used by TFoam.
Definition: TFoamCell.h:12
void SetXdiv(Double_t Xdiv)
Definition: TFoamCell.h:45
TFoamCell * GetDau1() const
Definition: TFoamCell.h:62
void CalcVolume()
Calculates volume of the cell using size params which are calculated.
Definition: TFoamCell.cxx:170
Double_t GetVolume() const
Definition: TFoamCell.h:50
Double_t GetDriv() const
Definition: TFoamCell.h:52
Int_t GetStat() const
Definition: TFoamCell.h:58
TFoamCell * GetDau0() const
Definition: TFoamCell.h:61
void SetStat(Int_t Stat)
Definition: TFoamCell.h:59
Double_t GetPrim() const
Definition: TFoamCell.h:53
void Print(Option_t *option) const override
Printout of the cell geometry parameters for the debug purpose.
Definition: TFoamCell.cxx:185
Int_t GetBest() const
Definition: TFoamCell.h:43
void GetHcub(TFoamVect &, TFoamVect &) const
Provides size and position of the cell These parameter are calculated by analyzing information in all...
Definition: TFoamCell.cxx:116
TFoamCell * GetPare() const
Definition: TFoamCell.h:60
void SetSerial(Int_t Serial)
Definition: TFoamCell.h:65
Double_t GetIntg() const
Definition: TFoamCell.h:51
void Fill(Int_t, TFoamCell *, TFoamCell *, TFoamCell *)
Fills in certain data into newly allocated cell.
Definition: TFoamCell.cxx:103
void SetDau1(TFoamCell *Daug)
Definition: TFoamCell.h:64
void SetPrim(Double_t Prim)
Definition: TFoamCell.h:56
void SetDriv(Double_t Driv)
Definition: TFoamCell.h:55
void SetIntg(Double_t Intg)
Definition: TFoamCell.h:54
void SetDau0(TFoamCell *Daug)
Definition: TFoamCell.h:63
void SetBest(Int_t Best)
Definition: TFoamCell.h:44
Abstract class representing n-dimensional real positive integrand function.
Definition: TFoamIntegrand.h:9
virtual Double_t Density(Int_t ndim, Double_t *)=0
Small auxiliary class for controlling MC weight.
Definition: TFoamMaxwt.h:12
void Fill(Double_t)
Filling analyzed weight.
Definition: TFoamMaxwt.cxx:93
void GetMCeff(Double_t, Double_t &, Double_t &)
Calculates Efficiency= aveWt/wtLim for a given tolerance level epsilon<<1 using information stored in...
Definition: TFoamMaxwt.cxx:120
Auxiliary class TFoamVect of n-dimensional vector, with dynamic allocation used for the cartesian geo...
Definition: TFoamVect.h:10
TFoam is the main class of the multi-dimensional general purpose Monte Carlo event generator (integra...
Definition: TFoam.h:21
virtual void MakeActiveList()
Internal method used by Initialize.
Definition: TFoam.cxx:967
virtual void CheckAll(Int_t)
User utility, miscellaneous and debug.
Definition: TFoam.cxx:1380
TString fName
Name of a given instance of the FOAM class.
Definition: TFoam.h:24
Double_t fMCwt
MC weight.
Definition: TFoam.h:57
virtual Int_t Divide(TFoamCell *)
Internal method used by Initialize.
Definition: TFoam.cxx:934
virtual void PrintCells()
Prints geometry of ALL cells of the FOAM.
Definition: TFoam.cxx:1448
TH1D * fHistWt
Histogram of the MC wt.
Definition: TFoam.h:54
virtual void GetMCvect(Double_t *)
User may get generated MC point/vector with help of this method.
Definition: TFoam.cxx:1174
TFoamIntegrand * fRho
! Pointer to the user-defined integrand function/distribution
Definition: TFoam.h:60
TRandom * fPseRan
Pointer to user-defined generator of pseudorandom numbers.
Definition: TFoam.h:62
Double_t fMCerror
and its error
Definition: TFoam.h:72
Double_t fPrime
Primary integral R' (R=R'<wt>)
Definition: TFoam.h:70
Int_t fChat
Chat=0,1,2 chat level in output, Chat=1 normal level.
Definition: TFoam.h:32
virtual void InitCells()
Internal method used by Initialize.
Definition: TFoam.cxx:447
Double_t fNevGen
Total number of the generated MC events.
Definition: TFoam.h:68
virtual void MakeEvent()
User method.
Definition: TFoam.cxx:1124
Double_t fSumWt2
Total sum of wt and wt^2.
Definition: TFoam.h:66
Double_t fSumWt
Definition: TFoam.h:66
TFoamVect ** fXdivPRD
! Lists of division values encoded in one vector per direction
Definition: TFoam.h:42
~TFoam() override
Default destructor.
Definition: TFoam.cxx:230
virtual Long_t PeekMax()
Internal method used by Initialize.
Definition: TFoam.cxx:900
virtual void GenerCel2(TFoamCell *&)
Internal method.
Definition: TFoam.cxx:1086
virtual void Initialize()
Basic initialization of FOAM invoked by the user.
Definition: TFoam.cxx:323
virtual void GetIntNorm(Double_t &, Double_t &)
User method.
Definition: TFoam.cxx:1228
TFoamMaxwt * fMCMonit
Monitor of the MC weight for measuring MC efficiency.
Definition: TFoam.h:48
Int_t fRNmax
Maximum No. of the rand. numb. requested at once.
Definition: TFoam.h:29
Int_t * fInhiDiv
! [fDim] Flags for inhibiting cell division
Definition: TFoam.h:40
Long_t fNCalls
Total number of the function calls.
Definition: TFoam.h:64
virtual Double_t GetMCwt()
User may get weight MC weight using this method.
Definition: TFoam.cxx:1182
Double_t fWtMin
Maximum/Minimum MC weight.
Definition: TFoam.h:69
virtual void ResetRho(TFoamIntegrand *Rho)
User may optionally reset the distribution using this method.
Definition: TFoam.cxx:1051
Int_t fEvPerBin
Maximum number of effective (wt=1) events per bin.
Definition: TFoam.h:37
virtual void GetIntegMC(Double_t &, Double_t &)
User method.
Definition: TFoam.cxx:1209
Double_t fMaxWtRej
Maximum weight in rejection for getting wt=1 events.
Definition: TFoam.h:49
Int_t fNCells
Maximum number of cells.
Definition: TFoam.h:28
Double_t * fRvec
[fRNmax] random number vector from r.n. generator fDim+1 maximum elements
Definition: TFoam.h:58
virtual void RootPlot2dim(Char_t *)
Debugging tool which plots 2-dimensional cells as rectangles in C++ format readable for root.
Definition: TFoam.cxx:1465
Double_t * fAlpha
[fDim] Internal parameters of the hyper-rectangle
Definition: TFoam.h:74
Int_t fOptDrive
Optimization switch =1,2 for variance or maximum weight optimization.
Definition: TFoam.h:31
virtual void Explore(TFoamCell *Cell)
Internal method used by Initialize.
Definition: TFoam.cxx:522
virtual Double_t Eval(Double_t *)
Internal method.
Definition: TFoam.cxx:1064
virtual void MakeAlpha()
Internal method used by Initialize.
Definition: TFoam.cxx:847
virtual void ResetPseRan(TRandom *PseRan)
User may optionally reset random number generator using this method.
Definition: TFoam.cxx:1008
TMethodCall * fMethodCall
! ROOT's pointer to user-defined global distribution function
Definition: TFoam.h:61
Double_t Sqr(Double_t x) const
Definition: TFoam.h:140
Int_t fLastCe
Index of the last cell.
Definition: TFoam.h:45
Int_t fOptPRD
Option switch for predefined division, for quick check.
Definition: TFoam.h:41
TObjArray * fHistDbg
Histograms of wt, for debug.
Definition: TFoam.h:53
virtual void SetXdivPRD(Int_t, Int_t, Double_t[])
This should be called before Initialize, after setting kDim It predefines values of the cell division...
Definition: TFoam.cxx:1340
Double_t fSumOve
Total Sum of overweighted events.
Definition: TFoam.h:67
virtual void SetInhiDiv(Int_t, Int_t)
This can be called before Initialize, after setting kDim It defines which variables are excluded in t...
Definition: TFoam.cxx:1310
Double_t fWtMax
Definition: TFoam.h:69
std::vector< Long_t > fCellsAct
Index of active cells, constructed at the end of foam build-up.
Definition: TFoam.h:50
Int_t * fMaskDiv
! [fDim] Dynamic Mask for cell division
Definition: TFoam.h:39
virtual void LinkCells(void)
Definition: TFoam.cxx:1519
virtual void SetRho(TFoamIntegrand *Rho)
User may use this method to set the distribution object.
Definition: TFoam.cxx:1020
virtual void SetRhoInt(Double_t(*fun)(Int_t, Double_t *))
User may use this method to set the distribution object as a global function pointer (and not as an i...
Definition: TFoam.cxx:1031
Int_t fNoAct
Number of active cells.
Definition: TFoam.h:44
virtual void Finalize(Double_t &, Double_t &)
May be called optionally by the user after the MC run.
Definition: TFoam.cxx:1258
Long_t fNEffev
Total number of effective events (wt=1) in the foam buildup.
Definition: TFoam.h:65
Int_t fOptRej
Switch =0 for weighted events; =1 for unweighted events in MC.
Definition: TFoam.h:33
Double_t * fPrimAcu
[fNoAct] Array of cumulative probability of all active cells
Definition: TFoam.h:51
virtual Double_t MCgenerate(Double_t *MCvect)
User method which generates MC event and returns MC weight.
Definition: TFoam.cxx:1197
Double_t * fMCvect
[fDim] Generated MC vector for the outside user
Definition: TFoam.h:56
TString fVersion
Actual version of the FOAM like (1.01m)
Definition: TFoam.h:25
virtual Int_t CellFill(Int_t, TFoamCell *)
Internal method used by Initialize.
Definition: TFoam.cxx:479
Int_t fDim
Dimension of the integration/simulation space.
Definition: TFoam.h:27
Int_t fNSampl
No. of MC events, when dividing (exploring) cell.
Definition: TFoam.h:36
TString fDate
Release date of FOAM.
Definition: TFoam.h:26
TFoamCell ** fCells
[fNCells] Array of ALL cells
Definition: TFoam.h:46
virtual void SetPseRan(TRandom *PseRan)
Definition: TFoam.h:112
TObjArray * fHistEdg
Histograms of wt, one for each cell edge.
Definition: TFoam.h:52
Double_t fMCresult
True Integral R from MC series.
Definition: TFoam.h:71
virtual void GetWtParams(Double_t, Double_t &, Double_t &, Double_t &)
May be called optionally after the MC run.
Definition: TFoam.cxx:1245
Int_t fNBin
No. of bins in the edge histogram for cell MC exploration.
Definition: TFoam.h:35
virtual void Grow()
Internal method used by Initialize.
Definition: TFoam.cxx:862
TFoam()
Default constructor for streamer, user should not use it.
Definition: TFoam.cxx:146
virtual void Varedu(Double_t[], Int_t &, Double_t &, Double_t &)
Internal method used by Initialize.
Definition: TFoam.cxx:674
virtual void Carver(Int_t &, Double_t &, Double_t &)
Internal method used by Initialize.
Definition: TFoam.cxx:751
1-D histogram with a double per channel (see TH1 documentation)}
Definition: TH1.h:617
void Reset(Option_t *option="") override
Reset.
Definition: TH1.cxx:10161
static void AddDirectory(Bool_t add=kTRUE)
Sets the flag controlling the automatic add of histograms in memory.
Definition: TH1.cxx:1267
virtual Int_t Fill(Double_t x)
Increment bin with abscissa X by 1.
Definition: TH1.cxx:3338
static Bool_t AddDirectoryStatus()
Static function: cannot be inlined on Windows/NT.
Definition: TH1.cxx:735
void Execute(const char *, const char *, int *=nullptr) override
Execute method on this object with the given parameter string, e.g.
Definition: TMethodCall.h:64
void SetParamPtrs(void *paramArr, Int_t nparam=-1)
ParamArr is an array containing the function argument values.
An array of TObjects.
Definition: TObjArray.h:31
void Delete(Option_t *option="") override
Remove all objects from the array AND delete all heap based objects.
Definition: TObjArray.cxx:356
Mother of all ROOT objects.
Definition: TObject.h:41
virtual void Warning(const char *method, const char *msgfmt,...) const
Issue warning message.
Definition: TObject.cxx:949
virtual void Error(const char *method, const char *msgfmt,...) const
Issue error message.
Definition: TObject.cxx:963
virtual void Info(const char *method, const char *msgfmt,...) const
Issue info message.
Definition: TObject.cxx:937
This is the base class for the ROOT Random number generators.
Definition: TRandom.h:27
virtual void RndmArray(Int_t n, Float_t *array)
Return an array of n random numbers uniformly distributed in ]0,1].
Definition: TRandom.cxx:588
Double_t Rndm() override
Machine independent random number generator.
Definition: TRandom.cxx:552
Basic string class.
Definition: TString.h:136
const char * Data() const
Definition: TString.h:369
const Double_t sigma
Double_t y[n]
Definition: legend1.C:17
Double_t x[n]
Definition: legend1.C:17
const Int_t n
Definition: legend1.C:16
VecExpr< UnaryOp< Sqrt< T >, VecExpr< A, T, D >, T >, T, D > sqrt(const VecExpr< A, T, D > &rhs)
Short_t Max(Short_t a, Short_t b)
Returns the largest of a and b.
Definition: TMathBase.h:250
Short_t Min(Short_t a, Short_t b)
Returns the smallest of a and b.
Definition: TMathBase.h:198
Short_t Abs(Short_t d)
Returns the absolute value of parameter Short_t d.
Definition: TMathBase.h:123
const char * Name
Definition: TXMLSetup.cxx:67
static uint64_t sum(uint64_t i)
Definition: Factory.cxx:2345