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.
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 "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
121class FoamIntegrandFunction : public TFoamIntegrand {
122
123public:
124
125 typedef Double_t (*FunctionPtr)(Int_t, Double_t*);
126
127 FoamIntegrandFunction(FunctionPtr func) : fFunc(func) {}
128
129 virtual ~FoamIntegrandFunction() {}
130
131 // evaluate the density using the provided function pointer
132 Double_t Density (Int_t nDim, Double_t * x) {
133 return fFunc(nDim,x);
134 }
135
136private:
137
138 FunctionPtr fFunc;
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,binTot;
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 binTot = 0.0;
784 for(iBin=0;iBin<fNBin;iBin++){
785 carvTot = carvTot + (binMax-bins[iBin]); // Total Carve (more stable)
786 binTot +=bins[iBin];
787 }
788 // primTot = binMax*fNBin;
789 //std::cout <<"Carver: CarvTot "<<CarvTot<< " primTot "<<primTot<<std::endl;
790 jLow =0;
791 jUp =fNBin-1;
792 carvOne = gVlow;
793 Double_t yLevel = gVlow;
794 for(iBin=0; iBin<fNBin;iBin++) {
795 theBin = bins[iBin];
796 //----- walk to the left and find first bin > theBin
797 iLow = iBin;
798 for(j=iBin; j>-1; j-- ) {
799 if(theBin< bins[j]) break;
800 iLow = j;
801 }
802 //iLow = iBin;
803 //if(iLow>0) while( (theBin >= bins[iLow-1])&&(iLow >0) ){iLow--;} // horror!!!
804 //------ walk to the right and find first bin > theBin
805 iUp = iBin;
806 for(j=iBin; j<fNBin; j++) {
807 if(theBin< bins[j]) break;
808 iUp = j;
809 }
810 //iUp = iBin;
811 //if(iUp<fNBin-1) while( (theBin >= bins[iUp+1])&&( iUp<fNBin-1 ) ){iUp++;} // horror!!!
812 //
813 carve = (iUp-iLow+1)*(binMax-theBin);
814 if( carve > carvOne) {
815 carvOne = carve;
816 jLow = iLow;
817 jUp = iUp;
818 yLevel = theBin;
819 }
820 }//iBin
821 if( carvTot > carvMax) {
822 carvMax = carvTot;
823 //primMax = primTot;
824 //std::cout <<"Carver: primMax "<<primMax<<std::endl;
825 kBest = kProj; // Best edge
826 xBest = ((Double_t)(jLow))/fNBin;
827 yBest = ((Double_t)(jUp+1))/fNBin;
828 if(jLow == 0 ) xBest = yBest;
829 if(jUp == fNBin-1) yBest = xBest;
830 // division ratio in units of 1/fNBin, testing
831 // jDivi = jLow;
832 // if(jLow == 0 ) jDivi=jUp+1;
833 }
834 //====== extra histograms for debug purposes
835 //std::cout<<"kProj= "<<kProj<<" jLow= "<<jLow<<" jUp= "<<jUp<<std::endl;
836 for(iBin=0; iBin<fNBin; iBin++)
837 ((TH1D *)(*fHistDbg)[kProj])->SetBinContent(iBin+1,binMax);
838 for(iBin=jLow; iBin<jUp+1; iBin++)
839 ((TH1D *)(*fHistDbg)[kProj])->SetBinContent(iBin+1,yLevel);
840 }//kProj
841 if( (kBest >= fDim) || (kBest<0) ) Error("Carver", "Something wrong with kBest - kBest = %d dim = %d\n",kBest,fDim);
842 delete [] bins;
843} //TFoam::Carver
844
845////////////////////////////////////////////////////////////////////////////////
846/// Internal method used by Initialize.
847/// Provides random vector Alpha 0< Alpha(i) < 1
848
850{
851 Int_t k;
852 if(fDim<1) return;
853
854 // simply generate and load kDim uniform random numbers
855 fPseRan->RndmArray(fDim,fRvec); // kDim random numbers needed
856 for(k=0; k<fDim; k++) fAlpha[k] = fRvec[k];
857} //MakeAlpha
858
859
860////////////////////////////////////////////////////////////////////////////////
861/// Internal method used by Initialize.
862/// It grow new cells by the binary division process.
863
865{
866 Long_t iCell;
867 TFoamCell* newCell;
868
869 while ( (fLastCe+2) < fNCells ) { // this condition also checked inside Divide
870 iCell = PeekMax(); // peek up cell with maximum driver integral
871 if( (iCell<0) || (iCell>fLastCe) ) {
872 Error("Grow", "Wrong iCell \n");
873 break;
874 }
875 newCell = fCells[iCell];
876
877 if(fLastCe !=0) {
878 Int_t kEcho=10;
879 if(fLastCe>=10000) kEcho=100;
880 if( (fLastCe%kEcho)==0) {
881 if (fChat>0) {
882 if(fDim<10)
883 std::cout<<fDim<<std::flush;
884 else
885 std::cout<<"."<<std::flush;
886 if( (fLastCe%(100*kEcho))==0) std::cout<<"|"<<fLastCe<<std::endl<<std::flush;
887 }
888 }
889 }
890 if( Divide( newCell )==0) break; // and divide it into two
891 }
892 if (fChat>0) {
893 std::cout<<std::endl<<std::flush;
894 }
895 CheckAll(0); // set arg=1 for more info
896}
897
898////////////////////////////////////////////////////////////////////////////////
899/// Internal method used by Initialize.
900/// It finds cell with maximal driver integral for the purpose of the division.
901
903{
904 Long_t i;
905 Long_t iCell = -1;
906 Double_t drivMax, driv;
907
908 drivMax = gVlow;
909 for(i=0; i<=fLastCe; i++) {//without root
910 if( fCells[i]->GetStat() == 1 ) {
911 driv = TMath::Abs( fCells[i]->GetDriv());
912 //std::cout<<"PeekMax: Driv = "<<driv<<std::endl;
913 if(driv > drivMax) {
914 drivMax = driv;
915 iCell = i;
916 }
917 }
918 }
919 // std::cout << 'TFoam_PeekMax: iCell=' << iCell << std::endl;
920 if (iCell == -1)
921 std::cout << "STOP in TFoam::PeekMax: not found iCell=" << iCell << std::endl;
922 return(iCell);
923} // TFoam_PeekMax
924
925////////////////////////////////////////////////////////////////////////////////
926/// Internal method used by Initialize.
927///
928/// It divides cell iCell into two daughter cells.
929/// The iCell is retained and tagged as inactive, daughter cells are appended
930/// at the end of the buffer.
931/// New vertex is added to list of vertices.
932/// List of active cells is updated, iCell removed, two daughters added
933/// and their properties set with help of MC sampling (TFoam_Explore)
934/// Returns Code RC=-1 of buffer limit is reached, fLastCe=fnBuf.
935
937{
938 Int_t kBest;
939
940 if(fLastCe+1 >= fNCells) Error("Divide", "Buffer limit is reached, fLastCe=fnBuf \n");
941
942 cell->SetStat(0); // reset cell as inactive
943 fNoAct--;
944
945 kBest = cell->GetBest();
946 if( kBest<0 || kBest>=fDim ) Error("Divide", "Wrong kBest \n");
947
948 //////////////////////////////////////////////////////////////////
949 // define two daughter cells (active) //
950 //////////////////////////////////////////////////////////////////
951
952 Int_t d1 = CellFill(1, cell);
953 Int_t d2 = CellFill(1, cell);
954 cell->SetDau0((fCells[d1]));
955 cell->SetDau1((fCells[d2]));
956 Explore( (fCells[d1]) );
957 Explore( (fCells[d2]) );
958 return 1;
959} // TFoam_Divide
960
961
962////////////////////////////////////////////////////////////////////////////////
963/// Internal method used by Initialize.
964///
965/// It finds out number of active cells fNoAct,
966/// creates list of active cell fCellsAct and primary cumulative fPrimAcu.
967/// They are used during the MC generation to choose randomly an active cell.
968
970{
971 Long_t iCell;
973 // flush previous result
974 if(fPrimAcu != 0) delete [] fPrimAcu;
975 fCellsAct.clear();
976 fCellsAct.reserve(fNoAct);
977
978 // Count Active cells and find total Primary
979 // Fill-in tables of active cells
980
981 fPrime = 0.0;
982 for(iCell=0; iCell<=fLastCe; iCell++) {
983 if (fCells[iCell]->GetStat()==1) {
984 fPrime += fCells[iCell]->GetPrim();
985 fCellsAct.push_back(iCell);
986 }
987 }
988
989 if(fNoAct != static_cast<Int_t>(fCellsAct.size())) Error("MakeActiveList", "Wrong fNoAct \n" );
990 if(fPrime == 0.) Error("MakeActiveList", "Integrand function is zero \n" );
991
992 fPrimAcu = new Double_t[fNoAct]; // cumulative primary for MC generation
993 if( fPrimAcu==0 ) Error("MakeActiveList", "Cant allocate fCellsAct or fPrimAcu \n");
994
995 sum =0.0;
996 for(iCell=0; iCell<fNoAct; iCell++) {
997 sum = sum + ( fCells[fCellsAct[iCell]] )->GetPrim()/fPrime;
998 fPrimAcu[iCell]=sum;
999 }
1000
1001} //MakeActiveList
1002
1003////////////////////////////////////////////////////////////////////////////////
1004/// User may optionally reset random number generator using this method.
1005///
1006/// Usually it is done when FOAM object is restored from the disk.
1007/// IMPORTANT: this method deletes existing random number generator registered in the FOAM object.
1008/// In particular such an object is created by the streamer during the disk-read operation.
1009
1011{
1012 if(fPseRan) {
1013 Info("ResetPseRan", "Resetting random number generator \n");
1014 delete fPseRan;
1015 }
1016 SetPseRan(PseRan);
1017}
1018
1019////////////////////////////////////////////////////////////////////////////////
1020/// User may use this method to set the distribution object
1021
1023{
1024 if (fun)
1025 fRho=fun;
1026 else
1027 Error("SetRho", "Bad function \n" );
1028}
1029////////////////////////////////////////////////////////////////////////////////
1030/// User may use this method to set the distribution object as a global function pointer
1031/// (and not as an interpreted function).
1032
1034{
1035 // This is needed for both AClic and Cling
1036 if (fun) {
1037 // delete function object if it has been created here in SetRho
1038 if (fRho && dynamic_cast<FoamIntegrandFunction*>(fRho) ) delete fRho;
1039 fRho= new FoamIntegrandFunction(fun);
1040 } else
1041 Error("SetRho", "Bad function \n" );
1042}
1043
1044////////////////////////////////////////////////////////////////////////////////
1045/// User may optionally reset the distribution using this method.
1046///
1047/// Usually it is done when FOAM object is restored from the disk.
1048/// IMPORTANT: this method deletes existing distribution object registered in the FOAM object.
1049/// In particular such an object is created by the streamer diring the disk-read operation.
1050/// This method is used only in very special cases, because the distribution in most cases
1051/// should be "owned" by the FOAM object and should not be replaced by another one after initialization.
1052
1054{
1055 if(fRho) {
1056 Info("ResetRho", "!!! Resetting distribution function !!!\n");
1057 delete fRho;
1058 }
1059 SetRho(fun);
1060}
1061
1062////////////////////////////////////////////////////////////////////////////////
1063/// Internal method.
1064/// Evaluates distribution to be generated.
1065
1067{
1068 Double_t result;
1069
1070 if(!fRho) { //interactive mode
1071 Long_t paramArr[3];
1072 paramArr[0]=(Long_t)fDim;
1073 paramArr[1]=(Long_t)xRand;
1074 fMethodCall->SetParamPtrs(paramArr);
1075 fMethodCall->Execute(result);
1076 } else { //compiled mode
1077 result=fRho->Density(fDim,xRand);
1078 }
1079
1080 return result;
1081}
1082
1083////////////////////////////////////////////////////////////////////////////////
1084/// Internal method.
1085/// Return randomly chosen active cell with probability equal to its
1086/// contribution into total driver integral using interpolation search.
1087
1089{
1090 Long_t lo, hi, hit;
1091 Double_t fhit, flo, fhi;
1092 Double_t random;
1093
1094 random=fPseRan->Rndm();
1095 lo = 0; hi =fNoAct-1;
1096 flo = fPrimAcu[lo]; fhi=fPrimAcu[hi];
1097 while(lo+1<hi) {
1098 hit = lo + (Int_t)( (hi-lo)*(random-flo)/(fhi-flo)+0.5);
1099 if (hit<=lo)
1100 hit = lo+1;
1101 else if(hit>=hi)
1102 hit = hi-1;
1103 fhit=fPrimAcu[hit];
1104 if (fhit>random) {
1105 hi = hit;
1106 fhi = fhit;
1107 } else {
1108 lo = hit;
1109 flo = fhit;
1110 }
1111 }
1112 if (fPrimAcu[lo]>random)
1113 pCell = fCells[fCellsAct[lo]];
1114 else
1115 pCell = fCells[fCellsAct[hi]];
1116} // TFoam::GenerCel2
1117
1118
1119////////////////////////////////////////////////////////////////////////////////
1120/// User method.
1121/// It generates randomly point/vector according to user-defined distribution.
1122/// Prior initialization with help of Initialize() is mandatory.
1123/// Generated MC point/vector is available using GetMCvect and the MC weight with GetMCwt.
1124/// MC point is generated with wt=1 or with variable weight, see OptRej switch.
1125
1127{
1128 Int_t j;
1129 Double_t wt,dx,mcwt;
1130 TFoamCell *rCell;
1131 //
1132 //********************** MC LOOP STARS HERE **********************
1133ee0:
1134 GenerCel2(rCell); // choose randomly one cell
1135
1136 MakeAlpha();
1137
1138 TFoamVect cellPosi(fDim); TFoamVect cellSize(fDim);
1139 rCell->GetHcub(cellPosi,cellSize);
1140 for(j=0; j<fDim; j++)
1141 fMCvect[j]= cellPosi[j] +fAlpha[j]*cellSize[j];
1142 dx = rCell->GetVolume(); // Cartesian volume of the Cell
1143 // weight average normalized to PRIMARY integral over the cell
1144
1145 wt=dx*Eval(fMCvect);
1146
1147 mcwt = wt / rCell->GetPrim(); // PRIMARY controls normalization
1148 fNCalls++;
1149 fMCwt = mcwt;
1150 // accumulation of statistics for the main MC weight
1151 fSumWt += mcwt; // sum of Wt
1152 fSumWt2 += mcwt*mcwt; // sum of Wt**2
1153 fNevGen++; // sum of 1d0
1154 fWtMax = TMath::Max(fWtMax, mcwt); // maximum wt
1155 fWtMin = TMath::Min(fWtMin, mcwt); // minimum wt
1156 fMCMonit->Fill(mcwt);
1157 fHistWt->Fill(mcwt,1.0); // histogram
1158 //******* Optional rejection ******
1159 if(fOptRej == 1) {
1160 Double_t random;
1161 random=fPseRan->Rndm();
1162 if( fMaxWtRej*random > fMCwt) goto ee0; // Wt=1 events, internal rejection
1163 if( fMCwt<fMaxWtRej ) {
1164 fMCwt = 1.0; // normal Wt=1 event
1165 } else {
1166 fMCwt = fMCwt/fMaxWtRej; // weight for overweighted events! kept for debug
1167 fSumOve += fMCwt-fMaxWtRej; // contribution of overweighted
1168 }
1169 }
1170 //********************** MC LOOP ENDS HERE **********************
1171} // MakeEvent
1172
1173////////////////////////////////////////////////////////////////////////////////
1174/// User may get generated MC point/vector with help of this method
1175
1177{
1178 for ( Int_t k=0 ; k<fDim ; k++) *(MCvect +k) = fMCvect[k];
1179}
1180
1181////////////////////////////////////////////////////////////////////////////////
1182/// User may get weight MC weight using this method
1183
1185{
1186 return(fMCwt);
1187}
1188////////////////////////////////////////////////////////////////////////////////
1189/// User may get weight MC weight using this method
1190
1192{
1193 mcwt=fMCwt;
1194}
1195
1196////////////////////////////////////////////////////////////////////////////////
1197/// User method which generates MC event and returns MC weight
1198
1200{
1201 MakeEvent();
1202 GetMCvect(MCvect);
1203 return(fMCwt);
1204}
1205
1206////////////////////////////////////////////////////////////////////////////////
1207/// User method.
1208/// It provides the value of the integral calculated from the averages of the MC run
1209/// May be called after (or during) the MC run.
1210
1211void TFoam::GetIntegMC(Double_t &mcResult, Double_t &mcError)
1212{
1213 Double_t mCerelat;
1214 mcResult = 0.0;
1215 mCerelat = 1.0;
1216 if (fNevGen>0) {
1217 mcResult = fPrime*fSumWt/fNevGen;
1218 mCerelat = sqrt( fSumWt2/(fSumWt*fSumWt) - 1/fNevGen);
1219 }
1220 mcError = mcResult *mCerelat;
1221}
1222
1223////////////////////////////////////////////////////////////////////////////////
1224/// User method.
1225/// It returns NORMALIZATION integral to be combined with the average weights
1226/// and content of the histograms in order to get proper absolute normalization
1227/// of the integrand and distributions.
1228/// It can be called after initialization, before or during the MC run.
1229
1230void TFoam::GetIntNorm(Double_t& IntNorm, Double_t& Errel )
1231{
1232 if(fOptRej == 1) { // Wt=1 events, internal rejection
1233 Double_t intMC,errMC;
1234 GetIntegMC(intMC,errMC);
1235 IntNorm = intMC;
1236 Errel = errMC;
1237 } else { // Wted events, NO internal rejection
1238 IntNorm = fPrime;
1239 Errel = 0;
1240 }
1241}
1242
1243////////////////////////////////////////////////////////////////////////////////
1244/// May be called optionally after the MC run.
1245/// Returns various parameters of the MC weight for efficiency evaluation
1246
1248{
1249 Double_t mCeff, wtLim;
1250 fMCMonit->GetMCeff(eps, mCeff, wtLim);
1251 wtMax = wtLim;
1252 aveWt = fSumWt/fNevGen;
1253 sigma = sqrt( fSumWt2/fNevGen -aveWt*aveWt );
1254}
1255
1256////////////////////////////////////////////////////////////////////////////////
1257/// May be called optionally by the user after the MC run.
1258/// It provides normalization and also prints some information/statistics on the MC run.
1259
1260void TFoam::Finalize(Double_t& IntNorm, Double_t& Errel)
1261{
1262 GetIntNorm(IntNorm,Errel);
1263 Double_t mcResult,mcError;
1264 GetIntegMC(mcResult,mcError);
1265 Double_t mCerelat= mcError/mcResult;
1266 //
1267 if(fChat>0) {
1268 Double_t eps = 0.0005;
1269 Double_t mCeff, mcEf2, wtMax, aveWt, sigma;
1270 GetWtParams(eps, aveWt, wtMax, sigma);
1271 mCeff=0;
1272 if(wtMax>0.0) mCeff=aveWt/wtMax;
1273 mcEf2 = sigma/aveWt;
1274 Double_t driver = fCells[0]->GetDriv();
1275 //
1276 BXOPE;
1277 BXTXT("****************************************");
1278 BXTXT("****** TFoam::Finalize ******");
1279 BXTXT("****************************************");
1280 BX1I(" NevGen",fNevGen, "Number of generated events in the MC generation ");
1281 BX1I(" nCalls",fNCalls, "Total number of function calls ");
1282 BXTXT("----------------------------------------");
1283 BX1F(" AveWt",aveWt, "Average MC weight ");
1284 BX1F(" WtMin",fWtMin, "Minimum MC weight (absolute) ");
1285 BX1F(" WtMax",fWtMax, "Maximum MC weight (absolute) ");
1286 BXTXT("----------------------------------------");
1287 BX1F(" XPrime",fPrime, "Primary total integral, R_prime ");
1288 BX1F(" XDiver",driver, "Driver total integral, R_loss ");
1289 BXTXT("----------------------------------------");
1290 BX2F(" IntMC", mcResult, mcError, "Result of the MC Integral");
1291 BX1F(" mCerelat", mCerelat, "Relative error of the MC integral ");
1292 BX1F(" <w>/WtMax",mCeff, "MC efficiency, acceptance rate");
1293 BX1F(" Sigma/<w>",mcEf2, "MC efficiency, variance/ave_wt");
1294 BX1F(" WtMax",wtMax, "WtMax(esp= 0.0005) ");
1295 BX1F(" Sigma",sigma, "variance of MC weight ");
1296 if(fOptRej==1) {
1297 Double_t avOve=fSumOve/fSumWt;
1298 BX1F("<OveW>/<W>",avOve, "Contrib. of events wt>MaxWtRej");
1299 }
1300 BXCLO;
1301 }
1302} // Finalize
1303
1304////////////////////////////////////////////////////////////////////////////////
1305/// This can be called before Initialize, after setting kDim
1306/// It defines which variables are excluded in the process of the cell division.
1307/// For example 'FoamX->SetInhiDiv(1, 1);' inhibits division of y-variable.
1308/// The resulting map of cells in 2-dim. case will look as follows:
1309///
1310/// \image html foam_Map2.png width=400
1311
1312void TFoam::SetInhiDiv(Int_t iDim, Int_t InhiDiv)
1313{
1314
1315
1316 if(fDim==0) Error("TFoam","SetInhiDiv: fDim=0 \n");
1317 if(fInhiDiv == 0) {
1318 fInhiDiv = new Int_t[ fDim ];
1319 for(Int_t i=0; i<fDim; i++) fInhiDiv[i]=0;
1320 }
1321 //
1322 if( ( 0<=iDim) && (iDim<fDim)) {
1323 fInhiDiv[iDim] = InhiDiv;
1324 } else
1325 Error("SetInhiDiv:","Wrong iDim \n");
1326}
1327
1328////////////////////////////////////////////////////////////////////////////////
1329/// This should be called before Initialize, after setting kDim
1330/// It predefines values of the cell division for certain variable iDim.
1331/// For example setting 3 predefined division lines using:
1332///
1333/// ~~~ {.cpp}
1334/// xDiv[0]=0.30; xDiv[1]=0.40; xDiv[2]=0.65;
1335/// FoamX->SetXdivPRD(0,3,xDiv);
1336/// ~~~
1337///
1338/// results in the following 2-dim. pattern of the cells:
1339///
1340/// \image html foam_Map3.png width=400
1341
1342void TFoam::SetXdivPRD(Int_t iDim, Int_t len, Double_t xDiv[])
1343{
1344 Int_t i;
1345
1346 if(fDim<=0) Error("SetXdivPRD", "fDim=0 \n");
1347 if( len<1 ) Error("SetXdivPRD", "len<1 \n");
1348 // allocate list of pointers, if it was not done before
1349 if(fXdivPRD == 0) {
1350 fXdivPRD = new TFoamVect*[fDim];
1351 for(i=0; i<fDim; i++) fXdivPRD[i]=0;
1352 }
1353 // set division list for direction iDim in H-cubic space!!!
1354 if( ( 0<=iDim) && (iDim<fDim)) {
1355 fOptPRD =1; // !!!!
1356 if( fXdivPRD[iDim] != 0)
1357 Error("SetXdivPRD", "Second allocation of XdivPRD not allowed \n");
1358 fXdivPRD[iDim] = new TFoamVect(len); // allocate list of division points
1359 for(i=0; i<len; i++) {
1360 (*fXdivPRD[iDim])[i]=xDiv[i]; // set list of division points
1361 }
1362 } else {
1363 Error("SetXdivPRD", "Wrong iDim \n");
1364 }
1365 // Printing predefined division points
1366 std::cout<<" SetXdivPRD, idim= "<<iDim<<" len= "<<len<<" "<<std::endl;
1367 for(i=0; i<len; i++) {
1368 if (iDim < fDim) std::cout<< (*fXdivPRD[iDim])[i] <<" ";
1369 }
1370 std::cout<<std::endl;
1371 for(i=0; i<len; i++) std::cout<< xDiv[i] <<" ";
1372 std::cout<<std::endl;
1373 //
1374}
1375
1376////////////////////////////////////////////////////////////////////////////////
1377/// User utility, miscellaneous and debug.
1378/// Checks all pointers in the tree of cells. This is useful auto diagnostic.
1379/// - level=0, no printout, failures causes STOP
1380/// - level=1, printout, failures lead to WARNINGS only
1381
1383{
1384 Int_t errors, warnings;
1385 TFoamCell *cell;
1386 Long_t iCell;
1387
1388 errors = 0; warnings = 0;
1389 if (level==1) std::cout << "///////////////////////////// FOAM_Checks /////////////////////////////////" << std::endl;
1390 for(iCell=1; iCell<=fLastCe; iCell++) {
1391 cell = fCells[iCell];
1392 // checking general rules
1393 if( ((cell->GetDau0()==0) && (cell->GetDau1()!=0) ) ||
1394 ((cell->GetDau1()==0) && (cell->GetDau0()!=0) ) ) {
1395 errors++;
1396 if (level==1) Error("CheckAll","ERROR: Cell's no %ld has only one daughter \n",iCell);
1397 }
1398 if( (cell->GetDau0()==0) && (cell->GetDau1()==0) && (cell->GetStat()==0) ) {
1399 errors++;
1400 if (level==1) Error("CheckAll","ERROR: Cell's no %ld has no daughter and is inactive \n",iCell);
1401 }
1402 if( (cell->GetDau0()!=0) && (cell->GetDau1()!=0) && (cell->GetStat()==1) ) {
1403 errors++;
1404 if (level==1) Error("CheckAll","ERROR: Cell's no %ld has two daughters and is active \n",iCell);
1405 }
1406
1407 // checking parents
1408 if( (cell->GetPare())!=fCells[0] ) { // not child of the root
1409 if ( (cell != cell->GetPare()->GetDau0()) && (cell != cell->GetPare()->GetDau1()) ) {
1410 errors++;
1411 if (level==1) Error("CheckAll","ERROR: Cell's no %ld parent not pointing to this cell\n ",iCell);
1412 }
1413 }
1414
1415 // checking daughters
1416 if(cell->GetDau0()!=0) {
1417 if(cell != (cell->GetDau0())->GetPare()) {
1418 errors++;
1419 if (level==1) Error("CheckAll","ERROR: Cell's no %ld daughter 0 not pointing to this cell \n",iCell);
1420 }
1421 }
1422 if(cell->GetDau1()!=0) {
1423 if(cell != (cell->GetDau1())->GetPare()) {
1424 errors++;
1425 if (level==1) Error("CheckAll","ERROR: Cell's no %ld daughter 1 not pointing to this cell \n",iCell);
1426 }
1427 }
1428 }// loop after cells;
1429
1430 // Check for empty cells
1431 for(iCell=0; iCell<=fLastCe; iCell++) {
1432 cell = fCells[iCell];
1433 if( (cell->GetStat()==1) && (cell->GetDriv()==0) ) {
1434 warnings++;
1435 if(level==1) Warning("CheckAll", "Warning: Cell no. %ld is active but empty \n", iCell);
1436 }
1437 }
1438 // summary
1439 if(level==1){
1440 Info("CheckAll","Check has found %d errors and %d warnings \n",errors, warnings);
1441 }
1442 if(errors>0){
1443 Info("CheckAll","Check - found total %d errors \n",errors);
1444 }
1445} // Check
1446
1447////////////////////////////////////////////////////////////////////////////////
1448/// Prints geometry of ALL cells of the FOAM
1449
1451{
1452 Long_t iCell;
1453
1454 for(iCell=0; iCell<=fLastCe; iCell++) {
1455 std::cout<<"Cell["<<iCell<<"]={ ";
1456 //std::cout<<" "<< fCells[iCell]<<" "; // extra DEBUG
1457 std::cout<<std::endl;
1458 fCells[iCell]->Print("1");
1459 std::cout<<"}"<<std::endl;
1460 }
1461}
1462
1463////////////////////////////////////////////////////////////////////////////////
1464/// Debugging tool which plots 2-dimensional cells as rectangles
1465/// in C++ format readable for root
1466
1468{
1469 std::ofstream outfile(filename, std::ios::out);
1470 Double_t x1,y1,x2,y2,x,y;
1471 Long_t iCell;
1472 Double_t offs =0.1;
1473 Double_t lpag =1-2*offs;
1474 outfile<<"{" << std::endl;
1475 outfile<<"cMap = new TCanvas(\"Map1\",\" Cell Map \",600,600);"<<std::endl;
1476 //
1477 outfile<<"TBox*a=new TBox();"<<std::endl;
1478 outfile<<"a->SetFillStyle(0);"<<std::endl; // big frame
1479 outfile<<"a->SetLineWidth(4);"<<std::endl;
1480 outfile<<"a->SetLineColor(2);"<<std::endl;
1481 outfile<<"a->DrawBox("<<offs<<","<<offs<<","<<(offs+lpag)<<","<<(offs+lpag)<<");"<<std::endl;
1482 //
1483 outfile<<"TText*t=new TText();"<<std::endl; // text for numbering
1484 outfile<<"t->SetTextColor(4);"<<std::endl;
1485 if(fLastCe<51)
1486 outfile<<"t->SetTextSize(0.025);"<<std::endl; // text for numbering
1487 else if(fLastCe<251)
1488 outfile<<"t->SetTextSize(0.015);"<<std::endl;
1489 else
1490 outfile<<"t->SetTextSize(0.008);"<<std::endl;
1491 //
1492 outfile<<"TBox*b=new TBox();"<<std::endl; // single cell
1493 outfile <<"b->SetFillStyle(0);"<<std::endl;
1494 //
1495 if(fDim==2 && fLastCe<=2000) {
1496 TFoamVect cellPosi(fDim); TFoamVect cellSize(fDim);
1497 outfile << "// =========== Rectangular cells ==========="<< std::endl;
1498 for(iCell=1; iCell<=fLastCe; iCell++) {
1499 if( fCells[iCell]->GetStat() == 1) {
1500 fCells[iCell]->GetHcub(cellPosi,cellSize);
1501 x1 = offs+lpag*( cellPosi[0]); y1 = offs+lpag*( cellPosi[1]);
1502 x2 = offs+lpag*(cellPosi[0]+cellSize[0]); y2 = offs+lpag*(cellPosi[1]+cellSize[1]);
1503 // cell rectangle
1504 if(fLastCe<=2000)
1505 outfile<<"b->DrawBox("<<x1<<","<<y1<<","<<x2<<","<<y2<<");"<<std::endl;
1506 // cell number
1507 if(fLastCe<=250) {
1508 x = offs+lpag*(cellPosi[0]+0.5*cellSize[0]); y = offs+lpag*(cellPosi[1]+0.5*cellSize[1]);
1509 outfile<<"t->DrawText("<<x<<","<<y<<","<<"\""<<iCell<<"\""<<");"<<std::endl;
1510 }
1511 }
1512 }
1513 outfile<<"// ============== End Rectangles ==========="<< std::endl;
1514 }//kDim=2
1515 //
1516 //
1517 outfile << "}" << std::endl;
1518 outfile.close();
1519}
1520
1522{
1523// Void function for backward compatibility
1524
1525 Info("LinkCells", "VOID function for backward compatibility \n");
1526 return;
1527}
1528
#define e(i)
Definition: RSha256.hxx:103
static const double x2[5]
static const double x1[5]
int Int_t
Definition: RtypesCore.h:45
char Char_t
Definition: RtypesCore.h:37
const Bool_t kFALSE
Definition: RtypesCore.h:92
long Long_t
Definition: RtypesCore.h:54
double Double_t
Definition: RtypesCore.h:59
#define ClassImp(name)
Definition: Rtypes.h:364
#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
float type_of_call hi(const int &, const int &)
double sqrt(double)
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
void Print(Option_t *option) const
Printout of the cell geometry parameters for the debug purpose.
Definition: TFoamCell.cxx:185
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
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:969
virtual void CheckAll(Int_t)
User utility, miscellaneous and debug.
Definition: TFoam.cxx:1382
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:936
virtual void PrintCells()
Prints geometry of ALL cells of the FOAM.
Definition: TFoam.cxx:1450
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:1176
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:1126
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
virtual Long_t PeekMax()
Internal method used by Initialize.
Definition: TFoam.cxx:902
virtual void GenerCel2(TFoamCell *&)
Internal method.
Definition: TFoam.cxx:1088
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:1230
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:1184
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:1053
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:1211
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:1467
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:1066
virtual void MakeAlpha()
Internal method used by Initialize.
Definition: TFoam.cxx:849
virtual void ResetPseRan(TRandom *PseRan)
User may optionally reset random number generator using this method.
Definition: TFoam.cxx:1010
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:1342
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:1312
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
virtual ~TFoam()
Default destructor.
Definition: TFoam.cxx:230
Int_t * fMaskDiv
! [fDim] Dynamic Mask for cell division
Definition: TFoam.h:39
virtual void LinkCells(void)
Definition: TFoam.cxx:1521
virtual void SetRho(TFoamIntegrand *Rho)
User may use this method to set the distribution object.
Definition: TFoam.cxx:1022
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:1033
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:1260
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:1199
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:1247
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:864
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:618
virtual void Reset(Option_t *option="")
Reset.
Definition: TH1.cxx:10126
static void AddDirectory(Bool_t add=kTRUE)
Sets the flag controlling the automatic add of histograms in memory.
Definition: TH1.cxx:1282
virtual Int_t Fill(Double_t x)
Increment bin with abscissa X by 1.
Definition: TH1.cxx:3350
static Bool_t AddDirectoryStatus()
Static function: cannot be inlined on Windows/NT.
Definition: TH1.cxx:750
void Execute(const char *, const char *, int *=0)
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:37
virtual void Delete(Option_t *option="")
Remove all objects from the array AND delete all heap based objects.
Definition: TObjArray.cxx:356
Mother of all ROOT objects.
Definition: TObject.h:37
virtual void Warning(const char *method, const char *msgfmt,...) const
Issue warning message.
Definition: TObject.cxx:879
virtual void Error(const char *method, const char *msgfmt,...) const
Issue error message.
Definition: TObject.cxx:893
virtual void Info(const char *method, const char *msgfmt,...) const
Issue info message.
Definition: TObject.cxx:867
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
virtual Double_t Rndm()
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
Short_t Max(Short_t a, Short_t b)
Definition: TMathBase.h:212
Short_t Min(Short_t a, Short_t b)
Definition: TMathBase.h:180
Short_t Abs(Short_t d)
Definition: TMathBase.h:120
const char * Name
Definition: TXMLSetup.cxx:67
static uint64_t sum(uint64_t i)
Definition: Factory.cxx:2345