ROOT   Reference Guide
Searching...
No Matches
TFoam.cxx
Go to the documentation of this file.
1// @(#)root/foam:$Id$
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
14 Institute of Nuclear Physics, Cracow, Poland
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:
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.
67nCell may be increased up to about million cells for wildly peaked distributions.
68Increasing nSampl sometimes helps, but it may cost CPU time.
69MaxWtRej 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 "TObjArray.h"
89#include "TMethodCall.h"
90#include "TRandom.h"
91#include "TMath.h"
92#include "TInterpreter.h"
93
95
96//FFFFFF BoX-FORMATs for nice and flexible outputs
97#define BXOPE std::cout<<\
98"FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF"<<std::endl<<\
99"F F"<<std::endl
100#define BXTXT(text) std::cout<<\
101"F "<<std::setw(40)<< text <<" F"<<std::endl
102#define BX1I(name,numb,text) std::cout<<\
103"F "<<std::setw(10)<<name<<" = "<<std::setw(10)<<numb<<" = " <<std::setw(50)<<text<<" F"<<std::endl
104#define BX1F(name,numb,text) std::cout<<"F "<<std::setw(10)<<name<<\
105 " = "<<std::setw(15)<<std::setprecision(8)<<numb<<" = "<<std::setw(40)<<text<<" F"<<std::endl
106#define BX2F(name,numb,err,text) std::cout<<"F "<<std::setw(10)<<name<<\
107" = "<<std::setw(15)<<std::setprecision(8)<<numb<<" +- "<<std::setw(15)<<std::setprecision(8)<<err<< \
108 " = "<<std::setw(25)<<text<<" F"<<std::endl
109#define BXCLO std::cout<<\
110"F F"<<std::endl<<\
111"FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF"<<std::endl
112 //FFFFFF BoX-FORMATs ends here
113
114static const Double_t gHigh= 1.0e150;
115static const Double_t gVlow=-1.0e150;
116
117#define SW2 setprecision(7) << std::setw(12)
118
119// class to wrap a global function in a TFoamIntegrand function
121
122public:
123
125
127
129
130 // evaluate the density using the provided function pointer
131 Double_t Density (Int_t nDim, Double_t * x) override {
132 return fFunc(nDim,x);
133 }
134
135private:
136
138
139};
140
141
142////////////////////////////////////////////////////////////////////////////////
143/// Default constructor for streamer, user should not use it.
144
146 fDim(0), fNCells(0), fRNmax(0),
147 fOptDrive(0), fChat(0), fOptRej(0),
148 fNBin(0), fNSampl(0), fEvPerBin(0),
150 fNoAct(0), fLastCe(0),
151 fMCMonit(nullptr), fMaxWtRej(0), fPrimAcu(nullptr),
152 fHistEdg(nullptr), fHistDbg(nullptr), fHistWt(nullptr),
153 fMCvect(nullptr), fMCwt(0), fRvec(nullptr),
154 fRho(nullptr), fMethodCall(nullptr), fPseRan(nullptr),
155 fNCalls(0), fNEffev(0),
156 fSumWt(0), fSumWt2(0),
157 fSumOve(0), fNevGen(0),
158 fWtMax(0), fWtMin(0),
159 fPrime(0), fMCresult(0), fMCerror(0),
160 fAlpha(nullptr)
161{
162}
163////////////////////////////////////////////////////////////////////////////////
164/// User constructor, to be employed by the user
165
166TFoam::TFoam(const Char_t* Name) :
167 fDim(0), fNCells(0), fRNmax(0),
168 fOptDrive(0), fChat(0), fOptRej(0),
169 fNBin(0), fNSampl(0), fEvPerBin(0),
171 fNoAct(0), fLastCe(0),
172 fMCMonit(nullptr), fMaxWtRej(0), fPrimAcu(nullptr),
173 fHistEdg(nullptr), fHistDbg(nullptr), fHistWt(nullptr),
174 fMCvect(nullptr), fMCwt(0), fRvec(nullptr),
175 fRho(nullptr), fMethodCall(nullptr), fPseRan(nullptr),
176 fNCalls(0), fNEffev(0),
177 fSumWt(0), fSumWt2(0),
178 fSumOve(0), fNevGen(0),
179 fWtMax(0), fWtMin(0),
180 fPrime(0), fMCresult(0), fMCerror(0),
181 fAlpha(nullptr)
182{
183 if(strlen(Name) >129) {
184 Error("TFoam","Name too long %s \n",Name);
185 }
186 fName=Name; // Class name
187 fDate=" Release date: 2005.04.10"; // Release date
188 fVersion= "1.02M"; // Release version
190 fInhiDiv = nullptr; // Flag allowing to inhibit cell division in certain projection/edge
191 fXdivPRD = nullptr; // Lists of division values encoded in one vector per direction
192 fAlpha = nullptr;
193 fPrimAcu = nullptr;
194 fHistEdg = nullptr;
195 fHistWt = nullptr;
196 fHistDbg = nullptr;
197 fDim = 0; // dimension of hyp-cubical space
198 fNCells = 1000; // Maximum number of Cells, is usually re-set
199 fNSampl = 200; // No of sampling when dividing cell
200 fOptPRD = 0; // General Option switch for PRedefined Division, for quick check
201 fOptDrive = 2; // type of Drive =1,2 for TrueVol,Sigma,WtMax
202 fChat = 1; // Chat=0,1,2 chat level in output, Chat=1 normal level
203 fOptRej = 1; // OptRej=0, wted events; OptRej=1, wt=1 events
204 /////////////////////////////////////////////////////////////////////////////
205
206 fNBin = 8; // binning of edge-histogram in cell exploration
207 fEvPerBin =25; // maximum no. of EFFECTIVE event per bin, =0 option is inactive
208 //------------------------------------------------------
209 fNCalls = 0; // No of function calls
210 fNEffev = 0; // Total no of eff. wt=1 events in build=up
211 fLastCe =-1; // Index of the last cell
212 fNoAct = 0; // No of active cells (used in MC generation)
213 fWtMin = gHigh; // Minimal weight
214 fWtMax = gVlow; // Maximal weight
215 fMaxWtRej =1.10; // Maximum weight in rejection for getting wt=1 events
216 fPseRan = nullptr; // Initialize private copy of random number generator
217 fMCMonit = nullptr; // MC efficiency monitoring
218 fRho = nullptr; // pointer to abstract class providing function to integrate
219 fMCvect = nullptr;
220 fRvec = nullptr;
221 fPseRan = nullptr; // generator of pseudorandom numbers
222 fMethodCall=nullptr; // ROOT's pointer to global distribution function
223}
224
225////////////////////////////////////////////////////////////////////////////////
226/// Default destructor
227
229{
230 Int_t i;
231
232 if(fCells!= nullptr) {
233 for(i=0; i<fNCells; i++) delete fCells[i]; // TFoamCell*[]
234 delete [] fCells;
235 }
236 if (fRvec) delete [] fRvec; //double[]
237 if (fAlpha) delete [] fAlpha; //double[]
238 if (fMCvect) delete [] fMCvect; //double[]
239 if (fPrimAcu) delete [] fPrimAcu; //double[]
241 if (fInhiDiv) delete [] fInhiDiv; //int[]
242
243 if( fXdivPRD!= nullptr) {
244 for(i=0; i<fDim; i++) delete fXdivPRD[i]; // TFoamVect*[]
245 delete [] fXdivPRD;
246 }
247 if (fMCMonit) delete fMCMonit;
248 if (fHistWt) delete fHistWt;
249
250 // delete histogram arrays
251 if (fHistEdg) {
252 fHistEdg->Delete();
253 delete fHistEdg;
254 }
255 if (fHistDbg) {
256 fHistDbg->Delete();
257 delete fHistDbg;
258 }
259 // delete function object if it has been created here in SetRhoInt
260 if (fRho && dynamic_cast<FoamIntegrandFunction*>(fRho) ) delete fRho;
261}
262
263
264////////////////////////////////////////////////////////////////////////////////
265/// Copy Constructor NOT IMPLEMENTED (NEVER USED)
266
267TFoam::TFoam(const TFoam &From): TObject(From)
268{
269 Error("TFoam", "COPY CONSTRUCTOR NOT IMPLEMENTED \n");
270}
271
272////////////////////////////////////////////////////////////////////////////////
273/// ### Basic initialization of FOAM invoked by the user. Mandatory!
274///
275/// This method starts the process of the cell build-up.
276/// User must invoke Initialize with two arguments or Initialize without arguments.
277/// This is done BEFORE generating first MC event and AFTER allocating FOAM object
278/// and resetting (optionally) its internal parameters/switches.
279/// The overall operational scheme of the FOAM is the following:
280///
281/// \image html foam_schema2.png width=600
282///
283/// ### This method invokes several other methods:
284///
285/// InitCells initializes memory storage for cells and begins exploration process
286/// from the root cell. The empty cells are allocated/filled using CellFill.
287/// The procedure Grow which loops over cells, picks up the cell with the biggest
288/// driver integral'', see Computer Physics Communications 152 152 (2003) 55 for explanations,
289/// with the help of PeekMax procedure. The chosen cell is split using Divide.
290/// Subsequently, the procedure Explore called by the Divide
291/// (and by InitCells for the root cell) does the most important
292/// job in the FOAM object build-up: it performs a small MC run for each
293/// newly allocated daughter cell.
294/// Explore calculates how profitable the future split of the cell will be
295/// and defines the optimal cell division geometry with the help of Carver or Varedu
296/// procedures, for maximum weight or variance optimization respectively.
297/// All essential results of the exploration are written into
298/// the explored cell object. At the very end of the foam build-up,
299/// Finally, MakeActiveList is invoked to create a list of pointers to
300/// all active cells, for the purpose of the quick access during the MC generation.
301/// The procedure Explore employs MakeAlpha to generate random coordinates
302/// inside a given cell with the uniform distribution.
303/// The above sequence of the procedure calls is depicted in the following figure:
304///
305/// \image html foam_Initialize_schema.png width=600
306
308{
309
310
311 SetPseRan(PseRan);
312 SetRho(fun);
313 Initialize();
314}
315
316////////////////////////////////////////////////////////////////////////////////
317/// Basic initialization of FOAM invoked by the user.
318/// IMPORTANT: Random number generator and the distribution object has to be
319/// provided using SetPseRan and SetRho prior to invoking this initialiser!
320
322{
325 Int_t i;
326
327 if(fChat>0){
328 BXOPE;
329 BXTXT("****************************************");
330 BXTXT("****** TFoam::Initialize ******");
331 BXTXT("****************************************");
332 BXTXT(fName);
333 BX1F(" Version",fVersion, fDate);
334 BX1I(" kDim",fDim, " Dimension of the hyper-cubical space ");
335 BX1I(" nCells",fNCells, " Requested number of Cells (half of them active) ");
336 BX1I(" nSampl",fNSampl, " No of MC events in exploration of a cell ");
337 BX1I(" nBin",fNBin, " No of bins in histograms, MC exploration of cell ");
338 BX1I(" EvPerBin",fEvPerBin, " Maximum No effective_events/bin, MC exploration ");
339 BX1I(" OptDrive",fOptDrive, " Type of Driver =1,2 for Sigma,WtMax ");
340 BX1I(" OptRej",fOptRej, " MC rejection on/off for OptRej=0,1 ");
341 BX1F(" MaxWtRej",fMaxWtRej, " Maximum wt in rejection for wt=1 evts");
342 BXCLO;
343 }
344
345 if(fPseRan==nullptr) Error("Initialize", "Random number generator not set \n");
346 if(fRho==nullptr && fMethodCall==nullptr ) Error("Initialize", "Distribution function not set \n");
347 if(fDim==0) Error("Initialize", "Zero dimension not allowed \n");
348
349 /////////////////////////////////////////////////////////////////////////
350 // ALLOCATE SMALL LISTS //
351 // it is done globally, not for each cell, to save on allocation time //
352 /////////////////////////////////////////////////////////////////////////
353 fRNmax= fDim+1;
354 fRvec = new Double_t[fRNmax]; // Vector of random numbers
355 if(fRvec==nullptr) Error("Initialize", "Cannot initialize buffer fRvec \n");
356
357 if(fDim>0){
358 fAlpha = new Double_t[fDim]; // sum<1 for internal parametrization of the simplex
359 if(fAlpha==nullptr) Error("Initialize", "Cannot initialize buffer fAlpha \n" );
360 }
361 fMCvect = new Double_t[fDim]; // vector generated in the MC run
362 if(fMCvect==nullptr) Error("Initialize", "Cannot initialize buffer fMCvect \n" );
363
364 //====== List of directions inhibited for division
365 if(fInhiDiv == nullptr){
366 fInhiDiv = new Int_t[fDim];
367 for(i=0; i<fDim; i++) fInhiDiv[i]=0;
368 }
369 //====== Dynamic mask used in Explore for edge determination
373 }
374 //====== List of predefined division values in all directions (initialized as empty)
375 if(fXdivPRD == nullptr){
376 fXdivPRD = new TFoamVect*[fDim];
377 for(i=0; i<fDim; i++) fXdivPRD[i]=nullptr; // Artificially extended beyond fDim
378 }
379 //====== Initialize list of histograms
380 fHistWt = new TH1D("HistWt","Histogram of MC weight",100,0.0, 1.5*fMaxWtRej); // MC weight
381 fHistEdg = new TObjArray(fDim); // Initialize list of histograms
382 TString hname;
383 TString htitle;
384 for(i=0;i<fDim;i++){
385 hname=fName+TString("_HistEdge_");
386 hname+=i;
387 htitle=TString("Edge Histogram No. ");
388 htitle+=i;
389 //std::cout<<"i= "<<i<<" hname= "<<hname<<" htitle= "<<htitle<<std::endl;
390 (*fHistEdg)[i] = new TH1D(hname.Data(),htitle.Data(),fNBin,0.0, 1.0); // Initialize histogram for each edge
391 ((TH1D*)(*fHistEdg)[i])->Sumw2();
392 }
393 //====== extra histograms for debug purposes
394 fHistDbg = new TObjArray(fDim); // Initialize list of histograms
395 for(i=0;i<fDim;i++){
396 hname=fName+TString("_HistDebug_");
397 hname+=i;
398 htitle=TString("Debug Histogram ");
399 htitle+=i;
400 (*fHistDbg)[i] = new TH1D(hname.Data(),htitle.Data(),fNBin,0.0, 1.0); // Initialize histogram for each edge
401 }
402
403 // ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| //
404 // BUILD-UP of the FOAM //
405 // ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| //
406 //
407 // Define and explore root cell(s)
408 InitCells();
409 // PrintCells(); std::cout<<" ===== after InitCells ====="<<std::endl;
410 Grow();
411 // PrintCells(); std::cout<<" ===== after Grow ====="<<std::endl;
412
413 MakeActiveList(); // Final Preparations for the M.C. generation
414
415 // Preparations for the M.C. generation
416 fSumWt = 0.0; // M.C. generation sum of Wt
417 fSumWt2 = 0.0; // M.C. generation sum of Wt**2
418 fSumOve = 0.0; // M.C. generation sum of overweighted
419 fNevGen = 0.0; // M.C. generation sum of 1d0
420 fWtMax = gVlow; // M.C. generation maximum wt
421 fWtMin = gHigh; // M.C. generation minimum wt
422 fMCresult= getCell(0)->GetIntg(); // M.C. Value of INTEGRAL,temporary assignment
423 fMCresult= getCell(0)->GetIntg(); // M.C. Value of INTEGRAL,temporary assignment
424 fMCerror = getCell(0)->GetIntg(); // M.C. Value of ERROR ,temporary assignment
425 fMCMonit = new TFoamMaxwt(5.0,1000); // monitoring M.C. efficiency
426 //
427 if(fChat>0){
428 Double_t driver = getCell(0)->GetDriv();
429 BXOPE;
430 BXTXT("*** TFoam::Initialize FINISHED!!! ***");
431 BX1I(" nCalls",fNCalls, "Total number of function calls ");
432 BX1F(" XPrime",fPrime, "Primary total integral ");
433 BX1F(" XDiver",driver, "Driver total integral ");
434 BX1F(" mcResult",fMCresult,"Estimate of the true MC Integral ");
435 BXCLO;
436 }
437 if(fChat==2) PrintCells();
439} // Initialize
440
441////////////////////////////////////////////////////////////////////////////////
442/// Internal method used by Initialize.
443/// It initializes "root part" of the FOAM of the tree of cells.
444
446{
447 Int_t i;
448
449 fLastCe =-1; // Index of the last cell
450 if(fCells!= nullptr) {
451 for(i=0; i<fNCells; i++) delete fCells[i];
452 delete [] fCells;
453 }
454 //
455 fCells = new TFoamCell*[fNCells];
456 for(i=0;i<fNCells;i++){
457 fCells[i]= new TFoamCell(fDim); // Allocate BIG list of cells
458 fCells[i]->SetSerial(i);
459 }
460 if(fCells==nullptr) Error("InitCells", "Cannot initialize CELLS \n" );
461
462 /////////////////////////////////////////////////////////////////////////////
463 // Single Root Hypercube //
464 /////////////////////////////////////////////////////////////////////////////
465 CellFill(1, nullptr); // 0-th cell ACTIVE
466
467 // Exploration of the root cell(s)
468 for(Long_t iCell=0; iCell<=fLastCe; iCell++){
469 Explore( getCell(iCell) ); // Exploration of root cell(s)
470 }
471}
472
473////////////////////////////////////////////////////////////////////////////////
474/// Internal method used by Initialize.
475/// It initializes content of the newly allocated active cell.
476
478{
479 TFoamCell *cell;
480 if (fLastCe==fNCells){
481 Error( "CellFill", "Too many cells\n");
482 }
483 fLastCe++; // 0-th cell is the first
484 if (Status==1) fNoAct++;
485
486 cell = getCell(fLastCe);
487
488 cell->Fill(Status, parent, nullptr, nullptr);
489
490 cell->SetBest( -1); // pointer for planning division of the cell
491 cell->SetXdiv(0.5); // factor for division
492 Double_t xInt2,xDri2;
493 if(parent!=nullptr){
494 xInt2 = 0.5*parent->GetIntg();
495 xDri2 = 0.5*parent->GetDriv();
496 cell->SetIntg(xInt2);
497 cell->SetDriv(xDri2);
498 }else{
499 cell->SetIntg(0.0);
500 cell->SetDriv(0.0);
501 }
502 return fLastCe;
503}
504
505////////////////////////////////////////////////////////////////////////////////
506/// Internal method used by Initialize.
507///
508/// It explores newly defined cell with help of special short MC sampling.
509/// As a result, estimates of true and drive volume is defined/determined
510/// Average and dispersion of the weight distribution will is found along
511/// each edge and the best edge (minimum dispersion, best maximum weight)
512/// is memorized for future use.
513///
514/// The optimal division point for eventual future cell division is
515/// determined/recorded. Recorded are also minimum and maximum weight etc.
516/// The volume estimate in all (inactive) parent cells is updated.
517/// Note that links to parents and initial volume = 1/2 parent has to be
518/// already defined prior to calling this routine.
519
521{
522 Double_t wt, dx, xBest=0, yBest=0;
523 Double_t intOld, driOld;
524
525 Long_t iev;
526 Double_t nevMC;
527 Int_t i, j, k;
528 Int_t nProj, kBest;
529 Double_t ceSum[5], xproj;
530
531 TFoamVect cellSize(fDim);
532 TFoamVect cellPosi(fDim);
533
534 cell->GetHcub(cellPosi,cellSize);
535
536 TFoamCell *parent;
537
538 Double_t *xRand = new Double_t[fDim];
539
540 Double_t *volPart=nullptr;
541
542 cell->CalcVolume();
543 dx = cell->GetVolume();
544 intOld = cell->GetIntg(); //memorize old values,
545 driOld = cell->GetDriv(); //will be needed for correcting parent cells
546
547
548 /////////////////////////////////////////////////////
549 // Special Short MC sampling to probe cell //
550 /////////////////////////////////////////////////////
551 ceSum[0]=0;
552 ceSum[1]=0;
553 ceSum[2]=0;
554 ceSum[3]=gHigh; //wtmin
555 ceSum[4]=gVlow; //wtmax
556 //
557 for(i=0;i<fDim;i++) ((TH1D *)(*fHistEdg)[i])->Reset(); // Reset histograms
558 fHistWt->Reset();
559 //
560 // ||||||||||||||||||||||||||BEGIN MC LOOP|||||||||||||||||||||||||||||
561 Double_t nevEff=0.;
562 for(iev=0;iev<fNSampl;iev++){
563 MakeAlpha(); // generate uniformly vector inside hypercube
564
565 if(fDim>0){
566 for(j=0; j<fDim; j++)
567 xRand[j]= cellPosi[j] +fAlpha[j]*(cellSize[j]);
568 }
569
570 wt=dx*Eval(xRand);
571
572 nProj = 0;
573 if(fDim>0) {
574 for(k=0; k<fDim; k++) {
575 xproj =fAlpha[k];
576 ((TH1D *)(*fHistEdg)[nProj])->Fill(xproj,wt);
577 nProj++;
578 }
579 }
580 //
581 fNCalls++;
582 ceSum[0] += wt; // sum of weights
583 ceSum[1] += wt*wt; // sum of weights squared
584 ceSum[2]++; // sum of 1
585 if (ceSum[3]>wt) ceSum[3]=wt; // minimum weight;
586 if (ceSum[4]<wt) ceSum[4]=wt; // maximum weight
587 // test MC loop exit condition
588 nevEff = ceSum[1] == 0. ? 0. : ceSum[0]*ceSum[0]/ceSum[1];
589 if( nevEff >= fNBin*fEvPerBin) break;
590 } // ||||||||||||||||||||||||||END MC LOOP|||||||||||||||||||||||||||||
591 //------------------------------------------------------------------
592 //--- predefine logics of searching for the best division edge ---
593 for(k=0; k<fDim;k++){
594 fMaskDiv[k] =1; // default is all
595 if( fInhiDiv[k]==1) fMaskDiv[k] =0; // inhibit some...
596 }
597 // Note that predefined division below overrule inhibition above
598 kBest=-1;
599 Double_t rmin,rmax,rdiv;
600 if(fOptPRD) { // quick check
601 for(k=0; k<fDim; k++) {
602 rmin= cellPosi[k];
603 rmax= cellPosi[k] +cellSize[k];
604 if( fXdivPRD[k] != nullptr) {
605 Int_t n= (fXdivPRD[k])->GetDim();
606 for(j=0; j<n; j++) {
607 rdiv=(*fXdivPRD[k])[j];
608 // check predefined divisions is available in this cell
609 if( (rmin +1e-99 <rdiv) && (rdiv< rmax -1e-99)) {
610 kBest=k;
611 xBest= (rdiv-cellPosi[k])/cellSize[k] ;
612 goto ee05;
613 }
614 }
615 }
616 }//k
617 }
618 ee05:
619 /////////////////////////////////////////////////////////////////////////////
620
621 fNEffev += (Long_t)nevEff;
622 nevMC = ceSum[2];
623 Double_t intTrue = ceSum[0]/(nevMC+0.000001);
624 Double_t intDriv=0.;
625 Double_t intPrim=0.;
626
627 switch(fOptDrive){
628 case 1: // VARIANCE REDUCTION
629 if(kBest == -1) Varedu(ceSum,kBest,xBest,yBest); // determine the best edge,
630 //intDriv =sqrt( ceSum[1]/nevMC -intTrue*intTrue ); // Older ansatz, numerically not bad
631 intDriv =sqrt(ceSum[1]/nevMC) -intTrue; // Foam build-up, sqrt(<w**2>) -<w>
632 intPrim =sqrt(ceSum[1]/nevMC); // MC gen. sqrt(<w**2>) =sqrt(<w>**2 +sigma**2)
633 break;
634 case 2: // WTMAX REDUCTION
635 if(kBest == -1) Carver(kBest,xBest,yBest); // determine the best edge
636 intDriv =ceSum[4] -intTrue; // Foam build-up, wtmax-<w>
637 intPrim =ceSum[4]; // MC generation, wtmax!
638 break;
639 default:
640 Error("Explore", "Wrong fOptDrive = \n" );
641 }//switch
642 //=================================================================================
643 //hist_Neff_distrib.Fill( fLastCe/2.0+0.01, nevEff+0.01); //
644 //hist_kBest_distrib.Fill( kBest+0.50, 1.0 ); // debug
645 //hist_xBest_distrib.Fill( xBest+0.01, 1.0 ); // debug
646 //=================================================================================
647 cell->SetBest(kBest);
648 cell->SetXdiv(xBest);
649 cell->SetIntg(intTrue);
650 cell->SetDriv(intDriv);
651 cell->SetPrim(intPrim);
652 // correct/update integrals in all parent cells to the top of the tree
653 Double_t parIntg, parDriv;
654 for(parent = cell->GetPare(); parent!=nullptr; parent = parent->GetPare()){
655 parIntg = parent->GetIntg();
656 parDriv = parent->GetDriv();
657 parent->SetIntg( parIntg +intTrue -intOld );
658 parent->SetDriv( parDriv +intDriv -driOld );
659 }
660 delete [] volPart;
661 delete [] xRand;
662 //cell->Print();
663} // TFoam::Explore
664
665////////////////////////////////////////////////////////////////////////////////
666/// Internal method used by Initialize.
667///
668/// It determines the best edge candidate and the position of the cell division plane
669/// in case of the variance reduction for future cell division,
670/// using results of the MC exploration run stored in fHistEdg
671
672void TFoam::Varedu(Double_t ceSum[5], Int_t &kBest, Double_t &xBest, Double_t &yBest)
673{
674 Double_t nent = ceSum[2];
675 Double_t swAll = ceSum[0];
676 Double_t sswAll = ceSum[1];
677 Double_t ssw = sqrt(sswAll)/sqrt(nent);
678 //
679 Double_t swIn,swOut,sswIn,sswOut,xLo,xUp;
680 kBest =-1;
681 xBest =0.5;
682 yBest =1.0;
683 Double_t maxGain=0.0;
684 // Now go over all projections kProj
685 for(Int_t kProj=0; kProj<fDim; kProj++) {
687 // initialize search over bins
688 Double_t sigmIn =0.0; Double_t sigmOut =0.0;
689 Double_t sswtBest = gHigh;
690 Double_t gain =0.0;
691 Double_t xMin=0.0; Double_t xMax=0.0;
692 // Double loop over all pairs jLo<jUp
693 for(Int_t jLo=1; jLo<=fNBin; jLo++) {
694 Double_t aswIn=0; Double_t asswIn=0;
695 for(Int_t jUp=jLo; jUp<=fNBin;jUp++) {
696 aswIn += ((TH1D *)(*fHistEdg)[kProj])->GetBinContent(jUp);
697 asswIn += Sqr(((TH1D *)(*fHistEdg)[kProj])->GetBinError( jUp));
698 xLo=(jLo-1.0)/fNBin;
699 xUp=(jUp*1.0)/fNBin;
700 swIn = aswIn/nent;
701 swOut = (swAll-aswIn)/nent;
702 sswIn = sqrt(asswIn) /sqrt(nent*(xUp-xLo)) *(xUp-xLo);
703 sswOut= sqrt(sswAll-asswIn)/sqrt(nent*(1.0-xUp+xLo)) *(1.0-xUp+xLo);
704 if( (sswIn+sswOut) < sswtBest) {
705 sswtBest = sswIn+sswOut;
706 gain = ssw-sswtBest;
707 sigmIn = sswIn -swIn; // Debug
708 sigmOut = sswOut-swOut; // Debug
709 xMin = xLo;
710 xMax = xUp;
711 }
712 }//jUp
713 }//jLo
714 Int_t iLo = (Int_t) (fNBin*xMin);
715 Int_t iUp = (Int_t) (fNBin*xMax);
716 //----------DEBUG printout
717 //std::cout<<"@@@@@ xMin xMax = "<<xMin <<" "<<xMax<<" iLo= "<<iLo<<" iUp= "<<iUp;
718 //std::cout<<" sswtBest/ssw= "<<sswtBest/ssw<<" Gain/ssw= "<< Gain/ssw<<std::endl;
719 //----------DEBUG auxiliary Plot
720 for(Int_t iBin=1;iBin<=fNBin;iBin++) {
721 if( ((iBin-0.5)/fNBin > xMin) && ((iBin-0.5)/fNBin < xMax) ){
722 ((TH1D *)(*fHistDbg)[kProj])->SetBinContent(iBin,sigmIn/(xMax-xMin));
723 } else {
724 ((TH1D *)(*fHistDbg)[kProj])->SetBinContent(iBin,sigmOut/(1-xMax+xMin));
725 }
726 }
727 if(gain>=maxGain) {
728 maxGain=gain;
729 kBest=kProj; // <--- !!!!! The best edge
730 xBest=xMin;
731 yBest=xMax;
732 if(iLo == 0) xBest=yBest; // The best division point
733 if(iUp == fNBin) yBest=xBest; // this is not really used
734 }
735 }
736 } //kProj
737 //----------DEBUG printout
738 //std::cout<<"@@@@@@@>>>>> kBest= "<<kBest<<" maxGain/ssw= "<< maxGain/ssw<<std::endl;
739 if( (kBest >= fDim) || (kBest<0) ) Error("Varedu", "Something wrong with kBest - kBest = %d dim = %d\n",kBest,fDim);
740} //TFoam::Varedu
741
742////////////////////////////////////////////////////////////////////////////////
743/// Internal method used by Initialize.
744///
745/// Determines the best edge-candidate and the position of the division plane
746/// for the future cell division, in the case of the optimization of the maximum weight.
747/// It exploits results of the cell MC exploration run stored in fHistEdg.
748
749void TFoam::Carver(Int_t &kBest, Double_t &xBest, Double_t &yBest)
750{
751 Int_t kProj,iBin;
752 Double_t carve,carvTot,carvMax,carvOne,binMax;
753 Int_t jLow,jUp,iLow,iUp;
754 Double_t theBin;
755 // Int_t jDivi; // TEST
756 Int_t j;
757
758 Double_t *bins = new Double_t[fNBin]; // bins of histogram for single PROJECTION
759 if(bins==nullptr) Error("Carver", "Cannot initialize buffer Bins \n" );
760
761 kBest =-1;
762 xBest =0.5;
763 yBest =1.0;
764 carvMax = gVlow;
765 // primMax = gVlow;
766 for(kProj=0; kProj<fDim; kProj++)
768 //if( kProj==1 ){
769 //std::cout<<"==================== Carver histogram: kProj ="<<kProj<<"==================="<<std::endl;
770 //((TH1D *)(*fHistEdg)[kProj])->Print("all");
771 binMax = gVlow;
772 for(iBin=0; iBin<fNBin;iBin++){
773 bins[iBin]= ((TH1D *)(*fHistEdg)[kProj])->GetBinContent(iBin+1);
774 binMax = TMath::Max( binMax, bins[iBin]); // Maximum content/bin
775 }
776 if(binMax < 0 ) { //case of empty cell
777 delete [] bins;
778 return;
779 }
780 carvTot = 0.0;
781 for(iBin=0;iBin<fNBin;iBin++){
782 carvTot = carvTot + (binMax-bins[iBin]); // Total Carve (more stable)
783 }
784 // primTot = binMax*fNBin;
785 //std::cout <<"Carver: CarvTot "<<CarvTot<< " primTot "<<primTot<<std::endl;
786 jLow =0;
787 jUp =fNBin-1;
788 carvOne = gVlow;
789 Double_t yLevel = gVlow;
790 for(iBin=0; iBin<fNBin;iBin++) {
791 theBin = bins[iBin];
792 //----- walk to the left and find first bin > theBin
793 iLow = iBin;
794 for(j=iBin; j>-1; j-- ) {
795 if(theBin< bins[j]) break;
796 iLow = j;
797 }
798 //iLow = iBin;
799 //if(iLow>0) while( (theBin >= bins[iLow-1])&&(iLow >0) ){iLow--;} // horror!!!
800 //------ walk to the right and find first bin > theBin
801 iUp = iBin;
802 for(j=iBin; j<fNBin; j++) {
803 if(theBin< bins[j]) break;
804 iUp = j;
805 }
806 //iUp = iBin;
807 //if(iUp<fNBin-1) while( (theBin >= bins[iUp+1])&&( iUp<fNBin-1 ) ){iUp++;} // horror!!!
808 //
809 carve = (iUp-iLow+1)*(binMax-theBin);
810 if( carve > carvOne) {
811 carvOne = carve;
812 jLow = iLow;
813 jUp = iUp;
814 yLevel = theBin;
815 }
816 }//iBin
817 if( carvTot > carvMax) {
818 carvMax = carvTot;
819 //primMax = primTot;
820 //std::cout <<"Carver: primMax "<<primMax<<std::endl;
821 kBest = kProj; // Best edge
822 xBest = ((Double_t)(jLow))/fNBin;
823 yBest = ((Double_t)(jUp+1))/fNBin;
824 if(jLow == 0 ) xBest = yBest;
825 if(jUp == fNBin-1) yBest = xBest;
826 // division ratio in units of 1/fNBin, testing
827 // jDivi = jLow;
828 // if(jLow == 0 ) jDivi=jUp+1;
829 }
830 //====== extra histograms for debug purposes
831 //std::cout<<"kProj= "<<kProj<<" jLow= "<<jLow<<" jUp= "<<jUp<<std::endl;
832 for(iBin=0; iBin<fNBin; iBin++)
833 ((TH1D *)(*fHistDbg)[kProj])->SetBinContent(iBin+1,binMax);
834 for(iBin=jLow; iBin<jUp+1; iBin++)
835 ((TH1D *)(*fHistDbg)[kProj])->SetBinContent(iBin+1,yLevel);
836 }//kProj
837 if( (kBest >= fDim) || (kBest<0) ) Error("Carver", "Something wrong with kBest - kBest = %d dim = %d\n",kBest,fDim);
838 delete [] bins;
839} //TFoam::Carver
840
841////////////////////////////////////////////////////////////////////////////////
842/// Internal method used by Initialize.
843/// Provides random vector Alpha 0< Alpha(i) < 1
844
846{
847 Int_t k;
848 if(fDim<1) return;
849
850 // simply generate and load kDim uniform random numbers
851 fPseRan->RndmArray(fDim,fRvec); // kDim random numbers needed
852 for(k=0; k<fDim; k++) fAlpha[k] = fRvec[k];
853} //MakeAlpha
854
855
856////////////////////////////////////////////////////////////////////////////////
857/// Internal method used by Initialize.
858/// It grow new cells by the binary division process.
859
861{
862 Long_t iCell;
863 TFoamCell* newCell;
864
865 while ( (fLastCe+2) < fNCells ) { // this condition also checked inside Divide
866 iCell = PeekMax(); // peek up cell with maximum driver integral
867 if( (iCell<0) || (iCell>fLastCe) ) {
868 Error("Grow", "Wrong iCell \n");
869 break;
870 }
871 newCell = getCell(iCell);
872
873 if(fLastCe !=0) {
874 Int_t kEcho=10;
875 if(fLastCe>=10000) kEcho=100;
876 if( (fLastCe%kEcho)==0) {
877 if (fChat>0) {
878 if(fDim<10)
879 std::cout<<fDim<<std::flush;
880 else
881 std::cout<<"."<<std::flush;
882 if( (fLastCe%(100*kEcho))==0) std::cout<<"|"<<fLastCe<<std::endl<<std::flush;
883 }
884 }
885 }
886 if( Divide( newCell )==0) break; // and divide it into two
887 }
888 if (fChat>0) {
889 std::cout<<std::endl<<std::flush;
890 }
892}
893
894////////////////////////////////////////////////////////////////////////////////
895/// Internal method used by Initialize.
896/// It finds cell with maximal driver integral for the purpose of the division.
897
899{
900 Long_t i;
901 Long_t iCell = -1;
902 Double_t drivMax, driv;
903
904 drivMax = gVlow;
905 for(i=0; i<=fLastCe; i++) {//without root
906 if( getCell(i)->GetStat() == 1 ) {
907 driv = TMath::Abs( getCell(i)->GetDriv());
908 //std::cout<<"PeekMax: Driv = "<<driv<<std::endl;
909 if(driv > drivMax) {
910 drivMax = driv;
911 iCell = i;
912 }
913 }
914 }
915 // std::cout << 'TFoam_PeekMax: iCell=' << iCell << std::endl;
916 if (iCell == -1)
917 std::cout << "STOP in TFoam::PeekMax: not found iCell=" << iCell << std::endl;
918 return(iCell);
919} // TFoam_PeekMax
920
921////////////////////////////////////////////////////////////////////////////////
922/// Internal method used by Initialize.
923///
924/// It divides cell iCell into two daughter cells.
925/// The iCell is retained and tagged as inactive, daughter cells are appended
926/// at the end of the buffer.
927/// New vertex is added to list of vertices.
928/// List of active cells is updated, iCell removed, two daughters added
929/// and their properties set with help of MC sampling (TFoam_Explore)
930/// Returns Code RC=-1 of buffer limit is reached, fLastCe=fnBuf.
931
933{
934 Int_t kBest;
935
936 if(fLastCe+1 >= fNCells) Error("Divide", "Buffer limit is reached, fLastCe=fnBuf \n");
937
938 cell->SetStat(0); // reset cell as inactive
939 fNoAct--;
940
941 kBest = cell->GetBest();
942 if( kBest<0 || kBest>=fDim ) Error("Divide", "Wrong kBest \n");
943
944 //////////////////////////////////////////////////////////////////
945 // define two daughter cells (active) //
946 //////////////////////////////////////////////////////////////////
947
948 Int_t d1 = CellFill(1, cell);
949 Int_t d2 = CellFill(1, cell);
950 cell->SetDau0((getCell(d1)));
951 cell->SetDau1((getCell(d2)));
952 Explore( (getCell(d1)) );
953 Explore( (getCell(d2)) );
954 return 1;
955} // TFoam_Divide
956
957
958////////////////////////////////////////////////////////////////////////////////
959/// Internal method used by Initialize.
960///
961/// It finds out number of active cells fNoAct,
962/// creates list of active cell fCellsAct and primary cumulative fPrimAcu.
963/// They are used during the MC generation to choose randomly an active cell.
964
966{
967 Long_t iCell;
969 // flush previous result
970 if(fPrimAcu != nullptr) delete [] fPrimAcu;
971 fCellsAct.clear();
972 fCellsAct.reserve(fNoAct);
973
974 // Count Active cells and find total Primary
975 // Fill-in tables of active cells
976
977 fPrime = 0.0;
978 for(iCell=0; iCell<=fLastCe; iCell++) {
979 if (getCell(iCell)->GetStat()==1) {
980 fPrime += getCell(iCell)->GetPrim();
981 fCellsAct.push_back(iCell);
982 }
983 }
984
985 if(fNoAct != static_cast<Int_t>(fCellsAct.size())) Error("MakeActiveList", "Wrong fNoAct \n" );
986 if(fPrime == 0.) Error("MakeActiveList", "Integrand function is zero \n" );
987
988 fPrimAcu = new Double_t[fNoAct]; // cumulative primary for MC generation
989 if( fPrimAcu==nullptr ) Error("MakeActiveList", "Cant allocate fCellsAct or fPrimAcu \n");
990
991 sum =0.0;
992 for(iCell=0; iCell<fNoAct; iCell++) {
993 sum = sum + ( getCell(fCellsAct[iCell]) )->GetPrim()/fPrime;
994 fPrimAcu[iCell]=sum;
995 }
996
997} //MakeActiveList
998
999////////////////////////////////////////////////////////////////////////////////
1000/// User may optionally reset random number generator using this method.
1001///
1002/// Usually it is done when FOAM object is restored from the disk.
1003/// IMPORTANT: this method deletes existing random number generator registered in the FOAM object.
1004/// In particular such an object is created by the streamer during the disk-read operation.
1005
1007{
1008 if(fPseRan) {
1009 Info("ResetPseRan", "Resetting random number generator \n");
1010 delete fPseRan;
1011 }
1012 SetPseRan(PseRan);
1013}
1014
1015////////////////////////////////////////////////////////////////////////////////
1016/// User may use this method to set the distribution object
1017
1019{
1020 if (fun)
1021 fRho=fun;
1022 else
1023 Error("SetRho", "Bad function \n" );
1024}
1025////////////////////////////////////////////////////////////////////////////////
1026/// User may use this method to set the distribution object as a global function pointer
1027/// (and not as an interpreted function).
1028
1030{
1031 // This is needed for both AClic and Cling
1032 if (fun) {
1033 // delete function object if it has been created here in SetRho
1034 if (fRho && dynamic_cast<FoamIntegrandFunction*>(fRho) ) delete fRho;
1035 fRho= new FoamIntegrandFunction(fun);
1036 } else
1037 Error("SetRho", "Bad function \n" );
1038}
1039
1040////////////////////////////////////////////////////////////////////////////////
1041/// User may optionally reset the distribution using this method.
1042///
1043/// Usually it is done when FOAM object is restored from the disk.
1044/// IMPORTANT: this method deletes existing distribution object registered in the FOAM object.
1045/// In particular such an object is created by the streamer diring the disk-read operation.
1046/// This method is used only in very special cases, because the distribution in most cases
1047/// should be "owned" by the FOAM object and should not be replaced by another one after initialization.
1048
1050{
1051 if(fRho) {
1052 Info("ResetRho", "!!! Resetting distribution function !!!\n");
1053 delete fRho;
1054 }
1055 SetRho(fun);
1056}
1057
1058////////////////////////////////////////////////////////////////////////////////
1059/// Internal method.
1060/// Evaluates distribution to be generated.
1061
1063{
1065
1066 if(!fRho) { //interactive mode
1067 Longptr_t paramArr[3];
1068 paramArr[0]=(Longptr_t)fDim;
1069 paramArr[1]=(Longptr_t)xRand;
1070 fMethodCall->SetParamPtrs(paramArr);
1072 } else { //compiled mode
1073 result=fRho->Density(fDim,xRand);
1074 }
1075
1076 return result;
1077}
1078
1079////////////////////////////////////////////////////////////////////////////////
1080/// Internal method.
1081/// Return randomly chosen active cell with probability equal to its
1082/// contribution into total driver integral using interpolation search.
1083
1085{
1086 Long_t lo, hi, hit;
1087 Double_t fhit, flo, fhi;
1088 Double_t random;
1089
1090 random=fPseRan->Rndm();
1091 lo = 0; hi =fNoAct-1;
1092 flo = fPrimAcu[lo]; fhi=fPrimAcu[hi];
1093 while(lo+1<hi) {
1094 hit = lo + (Int_t)( (hi-lo)*(random-flo)/(fhi-flo)+0.5);
1095 if (hit<=lo)
1096 hit = lo+1;
1097 else if(hit>=hi)
1098 hit = hi-1;
1099 fhit=fPrimAcu[hit];
1100 if (fhit>random) {
1101 hi = hit;
1102 fhi = fhit;
1103 } else {
1104 lo = hit;
1105 flo = fhit;
1106 }
1107 }
1108 if (fPrimAcu[lo]>random)
1109 pCell = getCell(fCellsAct[lo]);
1110 else
1111 pCell = getCell(fCellsAct[hi]);
1112} // TFoam::GenerCel2
1113
1114
1115////////////////////////////////////////////////////////////////////////////////
1116/// User method.
1117/// It generates randomly point/vector according to user-defined distribution.
1118/// Prior initialization with help of Initialize() is mandatory.
1119/// Generated MC point/vector is available using GetMCvect and the MC weight with GetMCwt.
1120/// MC point is generated with wt=1 or with variable weight, see OptRej switch.
1121
1123{
1124 Int_t j;
1125 Double_t wt,dx,mcwt;
1126 TFoamCell *rCell;
1127 //
1128 //********************** MC LOOP STARS HERE **********************
1129ee0:
1130 GenerCel2(rCell); // choose randomly one cell
1131
1132 MakeAlpha();
1133
1134 TFoamVect cellPosi(fDim); TFoamVect cellSize(fDim);
1135 rCell->GetHcub(cellPosi,cellSize);
1136 for(j=0; j<fDim; j++)
1137 fMCvect[j]= cellPosi[j] +fAlpha[j]*cellSize[j];
1138 dx = rCell->GetVolume(); // Cartesian volume of the Cell
1139 // weight average normalized to PRIMARY integral over the cell
1140
1141 wt=dx*Eval(fMCvect);
1142
1143 mcwt = wt / rCell->GetPrim(); // PRIMARY controls normalization
1144 fNCalls++;
1145 fMCwt = mcwt;
1146 // accumulation of statistics for the main MC weight
1147 fSumWt += mcwt; // sum of Wt
1148 fSumWt2 += mcwt*mcwt; // sum of Wt**2
1149 fNevGen++; // sum of 1d0
1150 fWtMax = TMath::Max(fWtMax, mcwt); // maximum wt
1151 fWtMin = TMath::Min(fWtMin, mcwt); // minimum wt
1152 fMCMonit->Fill(mcwt);
1153 fHistWt->Fill(mcwt,1.0); // histogram
1154 //******* Optional rejection ******
1155 if(fOptRej == 1) {
1156 Double_t random;
1157 random=fPseRan->Rndm();
1158 if( fMaxWtRej*random > fMCwt) goto ee0; // Wt=1 events, internal rejection
1159 if( fMCwt<fMaxWtRej ) {
1160 fMCwt = 1.0; // normal Wt=1 event
1161 } else {
1162 fMCwt = fMCwt/fMaxWtRej; // weight for overweighted events! kept for debug
1163 fSumOve += fMCwt-fMaxWtRej; // contribution of overweighted
1164 }
1165 }
1166 //********************** MC LOOP ENDS HERE **********************
1167} // MakeEvent
1168
1169////////////////////////////////////////////////////////////////////////////////
1170/// User may get generated MC point/vector with help of this method
1171
1173{
1174 for ( Int_t k=0 ; k<fDim ; k++) *(MCvect +k) = fMCvect[k];
1175}
1176
1177////////////////////////////////////////////////////////////////////////////////
1178/// User may get weight MC weight using this method
1179
1181{
1182 return(fMCwt);
1183}
1184////////////////////////////////////////////////////////////////////////////////
1185/// User may get weight MC weight using this method
1186
1188{
1189 mcwt=fMCwt;
1190}
1191
1192////////////////////////////////////////////////////////////////////////////////
1193/// User method which generates MC event and returns MC weight
1194
1196{
1197 MakeEvent();
1198 GetMCvect(MCvect);
1199 return(fMCwt);
1200}
1201
1202////////////////////////////////////////////////////////////////////////////////
1203/// User method.
1204/// It provides the value of the integral calculated from the averages of the MC run
1205/// May be called after (or during) the MC run.
1206
1207void TFoam::GetIntegMC(Double_t &mcResult, Double_t &mcError)
1208{
1209 Double_t mCerelat;
1210 mcResult = 0.0;
1211 mCerelat = 1.0;
1212 if (fNevGen>0) {
1213 mcResult = fPrime*fSumWt/fNevGen;
1214 mCerelat = sqrt( fSumWt2/(fSumWt*fSumWt) - 1/fNevGen);
1215 }
1216 mcError = mcResult *mCerelat;
1217}
1218
1219////////////////////////////////////////////////////////////////////////////////
1220/// User method.
1221/// It returns NORMALIZATION integral to be combined with the average weights
1222/// and content of the histograms in order to get proper absolute normalization
1223/// of the integrand and distributions.
1224/// It can be called after initialization, before or during the MC run.
1225/// \param IntNorm variable where the integral will be stored
1226/// \param ErrInt variable where the absolute error of the integral will be stored if internal rejection is active, otherwise set to 0.
1227
1228void TFoam::GetIntNorm(Double_t& IntNorm, Double_t& ErrInt )
1229{
1230 if(fOptRej == 1) { // Wt=1 events, internal rejection
1231 Double_t intMC,errMC;
1232 GetIntegMC(intMC,errMC);
1233 IntNorm = intMC;
1234 ErrInt = errMC;
1235 } else { // Wted events, NO internal rejection
1236 IntNorm = fPrime;
1237 ErrInt = 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 = getCell(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 == nullptr) {
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 == nullptr) {
1348 fXdivPRD = new TFoamVect*[fDim];
1349 for(i=0; i<fDim; i++) fXdivPRD[i]=nullptr;
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] != nullptr)
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 = getCell(iCell);
1390 // checking general rules
1391 if( ((cell->GetDau0()==nullptr) && (cell->GetDau1()!=nullptr) ) ||
1392 ((cell->GetDau1()==nullptr) && (cell->GetDau0()!=nullptr) ) ) {
1393 errors++;
1394 if (level==1) Error("CheckAll","ERROR: Cell's no %ld has only one daughter \n",iCell);
1395 }
1396 if( (cell->GetDau0()==nullptr) && (cell->GetDau1()==nullptr) && (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()!=nullptr) && (cell->GetDau1()!=nullptr) && (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())!=getCell(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()!=nullptr) {
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()!=nullptr) {
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 = getCell(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<<" "<< getCell(iCell)<<" "; // extra DEBUG
1455 std::cout<<std::endl;
1456 getCell(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( getCell(iCell)->GetStat() == 1) {
1498 getCell(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
1527TFoamCell* TFoam::getCell(std::size_t i) const
1528{
1529 // Ensure that the fCells member of the cells is filled
1530 if(fCells[i]->GetCells() == nullptr) {
1531 for (Int_t j=0; j < fNCells; ++j) {
1532 fCells[j]->SetCells(fCells);
1533 }
1534 }
1535
1536 return fCells[i];
1537}
#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
long Long_t
Definition RtypesCore.h:54
constexpr Bool_t kFALSE
Definition RtypesCore.h:101
double Double_t
Definition RtypesCore.h:59
#define ClassImp(name)
Definition Rtypes.h:377
#define BX1F(name, numb, text)
Definition TFoam.cxx:104
#define BX1I(name, numb, text)
Definition TFoam.cxx:102
#define BXCLO
Definition TFoam.cxx:109
#define BXTXT(text)
Definition TFoam.cxx:100
#define BX2F(name, numb, err, text)
Definition TFoam.cxx:106
static const Double_t gVlow
Definition TFoam.cxx:115
static const Double_t gHigh
Definition TFoam.cxx:114
#define BXOPE
Definition TFoam.cxx:97
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
Double_t Density(Int_t nDim, Double_t *x) override
Definition TFoam.cxx:131
~FoamIntegrandFunction() override
Definition TFoam.cxx:128
FoamIntegrandFunction(FunctionPtr func)
Definition TFoam.cxx:126
FunctionPtr fFunc
Definition TFoam.cxx:137
Double_t(* FunctionPtr)(Int_t, Double_t *)
Definition TFoam.cxx:124
Used by TFoam.
Definition TFoamCell.h:12
void SetXdiv(Double_t Xdiv)
Definition TFoamCell.h:59
TFoamCell * GetDau1() const
Definition TFoamCell.h:78
void CalcVolume()
Calculates volume of the cell using size params which are calculated.
Double_t GetVolume() const
Definition TFoamCell.h:64
Double_t GetDriv() const
Definition TFoamCell.h:66
Int_t GetStat() const
Definition TFoamCell.h:72
TFoamCell * GetDau0() const
Definition TFoamCell.h:77
void SetStat(Int_t Stat)
Definition TFoamCell.h:73
Double_t GetPrim() const
Definition TFoamCell.h:67
void Print(Option_t *option) const override
Printout of the cell geometry parameters for the debug purpose.
void SetCells(TFoamCell **cells)
Definition TFoamCell.h:84
Int_t GetBest() const
Definition TFoamCell.h:57
void GetHcub(TFoamVect &, TFoamVect &) const
Provides size and position of the cell These parameter are calculated by analyzing information in all...
Definition TFoamCell.cxx:71
TFoamCell * GetPare() const
Definition TFoamCell.h:76
void SetSerial(Int_t Serial)
Definition TFoamCell.h:82
Double_t GetIntg() const
Definition TFoamCell.h:65
void Fill(Int_t, TFoamCell *, TFoamCell *, TFoamCell *)
Fills in certain data into newly allocated cell.
Definition TFoamCell.cxx:58
void SetDau1(TFoamCell *Daug)
Definition TFoamCell.h:81
void SetPrim(Double_t Prim)
Definition TFoamCell.h:70
void SetDriv(Double_t Driv)
Definition TFoamCell.h:69
void SetIntg(Double_t Intg)
Definition TFoamCell.h:68
void SetDau0(TFoamCell *Daug)
Definition TFoamCell.h:80
void SetBest(Int_t Best)
Definition TFoamCell.h:58
Abstract class representing n-dimensional real positive integrand function.
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.
void GetMCeff(Double_t, Double_t &, Double_t &)
Calculates Efficiency= aveWt/wtLim for a given tolerance level epsilon<<1 using information stored in...
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
TFoamCell * getCell(std::size_t i) const
Definition TFoam.cxx:1527
virtual void MakeActiveList()
Internal method used by Initialize.
Definition TFoam.cxx:965
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:932
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:1172
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:445
Double_t fNevGen
Total number of the generated MC events.
Definition TFoam.h:68
virtual void MakeEvent()
User method.
Definition TFoam.cxx:1122
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:228
virtual Long_t PeekMax()
Internal method used by Initialize.
Definition TFoam.cxx:898
virtual void GenerCel2(TFoamCell *&)
Internal method.
Definition TFoam.cxx:1084
virtual void Initialize()
Basic initialization of FOAM invoked by the user.
Definition TFoam.cxx:321
virtual void Varedu(Double_t[5], Int_t &, Double_t &, Double_t &)
Internal method used by Initialize.
Definition TFoam.cxx:672
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:1180
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:1049
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:1207
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:520
virtual Double_t Eval(Double_t *)
Internal method.
Definition TFoam.cxx:1062
virtual void MakeAlpha()
Internal method used by Initialize.
Definition TFoam.cxx:845
virtual void ResetPseRan(TRandom *PseRan)
User may optionally reset random number generator using this method.
Definition TFoam.cxx:1006
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
! [fDim] Dynamic Mask for cell division
Definition TFoam.h:39
Definition TFoam.cxx:1519
virtual void SetRho(TFoamIntegrand *Rho)
User may use this method to set the distribution object.
Definition TFoam.cxx:1018
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:1029
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:1195
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:477
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:860
TFoam()
Default constructor for streamer, user should not use it.
Definition TFoam.cxx:145
virtual void Carver(Int_t &, Double_t &, Double_t &)
Internal method used by Initialize.
Definition TFoam.cxx:749
1-D histogram with a double per channel (see TH1 documentation)
Definition TH1.h:670
void Reset(Option_t *option="") override
Reset.
Definition TH1.cxx:10482
Sets the flag controlling the automatic add of histograms in memory.
Definition TH1.cxx:1294
virtual Int_t Fill(Double_t x)
Increment bin with abscissa X by 1.
Definition TH1.cxx:3344
Static function: cannot be inlined on Windows/NT.
Definition TH1.cxx:754
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.
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:973
virtual void Error(const char *method, const char *msgfmt,...) const
Issue error message.
Definition TObject.cxx:987
virtual void Info(const char *method, const char *msgfmt,...) const
Issue info message.
Definition TObject.cxx:961
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:595
Double_t Rndm() override
Machine independent random number generator.
Definition TRandom.cxx:559
Basic string class.
Definition TString.h:139
const char * Data() const
Definition TString.h:376
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
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
static uint64_t sum(uint64_t i)
Definition Factory.cxx:2345