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