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