Logo ROOT   6.08/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 /// User may use this to set pointer to the global function (not descending
1090 /// from TFoamIntegrand) serving as a distribution for FOAM.
1091 /// It is useful for simple interactive applications.
1092 /// Note that persistency for FOAM object will not work in the case of such
1093 /// a distribution.
1094 
1095 void TFoam::SetRhoInt( void * /*fun*/)
1096 {
1097  Error("SetRhoInt", "Not used anymore.");
1098  /*
1099  const char *namefcn = gCling->Getp2f2funcname(fun); //name of integrand function
1100  if(namefcn) {
1101  fMethodCall=new TMethodCall();
1102  fMethodCall->InitWithPrototype(namefcn, "Int_t, Double_t *");
1103  }
1104  */
1105  fRho=0;
1106 }
1107 
1108 ////////////////////////////////////////////////////////////////////////////////
1109 /// Internal subprogram.
1110 /// Evaluates distribution to be generated.
1111 
1113 {
1114  Double_t result;
1115 
1116  if(!fRho) { //interactive mode
1117  Long_t paramArr[3];
1118  paramArr[0]=(Long_t)fDim;
1119  paramArr[1]=(Long_t)xRand;
1120  fMethodCall->SetParamPtrs(paramArr);
1121  fMethodCall->Execute(result);
1122  } else { //compiled mode
1123  result=fRho->Density(fDim,xRand);
1124  }
1125 
1126  return result;
1127 }
1128 
1129 ////////////////////////////////////////////////////////////////////////////////
1130 /// Internal subprogram.
1131 /// Return randomly chosen active cell with probability equal to its
1132 /// contribution into total driver integral using interpolation search.
1133 
1135 {
1136  Long_t lo, hi, hit;
1137  Double_t fhit, flo, fhi;
1138  Double_t random;
1139 
1140  random=fPseRan->Rndm();
1141  lo = 0; hi =fNoAct-1;
1142  flo = fPrimAcu[lo]; fhi=fPrimAcu[hi];
1143  while(lo+1<hi) {
1144  hit = lo + (Int_t)( (hi-lo)*(random-flo)/(fhi-flo)+0.5);
1145  if (hit<=lo)
1146  hit = lo+1;
1147  else if(hit>=hi)
1148  hit = hi-1;
1149  fhit=fPrimAcu[hit];
1150  if (fhit>random) {
1151  hi = hit;
1152  fhi = fhit;
1153  } else {
1154  lo = hit;
1155  flo = fhit;
1156  }
1157  }
1158  if (fPrimAcu[lo]>random)
1159  pCell = (TFoamCell *) fCellsAct->At(lo);
1160  else
1161  pCell = (TFoamCell *) fCellsAct->At(hi);
1162 } // TFoam::GenerCel2
1163 
1164 
1165 ////////////////////////////////////////////////////////////////////////////////
1166 /// User subprogram.
1167 /// It generates randomly point/vector according to user-defined distribution.
1168 /// Prior initialization with help of Initialize() is mandatory.
1169 /// Generated MC point/vector is available using GetMCvect and the MC weight with GetMCwt.
1170 /// MC point is generated with wt=1 or with variable weight, see OptRej switch.
1171 
1173 {
1174  Int_t j;
1175  Double_t wt,dx,mcwt;
1176  TFoamCell *rCell;
1177  //
1178  //********************** MC LOOP STARS HERE **********************
1179 ee0:
1180  GenerCel2(rCell); // choose randomly one cell
1181 
1182  MakeAlpha();
1183 
1184  TFoamVect cellPosi(fDim); TFoamVect cellSize(fDim);
1185  rCell->GetHcub(cellPosi,cellSize);
1186  for(j=0; j<fDim; j++)
1187  fMCvect[j]= cellPosi[j] +fAlpha[j]*cellSize[j];
1188  dx = rCell->GetVolume(); // Cartesian volume of the Cell
1189  // weight average normalized to PRIMARY integral over the cell
1190 
1191  wt=dx*Eval(fMCvect);
1192 
1193  mcwt = wt / rCell->GetPrim(); // PRIMARY controls normalization
1194  fNCalls++;
1195  fMCwt = mcwt;
1196  // accumulation of statistics for the main MC weight
1197  fSumWt += mcwt; // sum of Wt
1198  fSumWt2 += mcwt*mcwt; // sum of Wt**2
1199  fNevGen++; // sum of 1d0
1200  fWtMax = TMath::Max(fWtMax, mcwt); // maximum wt
1201  fWtMin = TMath::Min(fWtMin, mcwt); // minimum wt
1202  fMCMonit->Fill(mcwt);
1203  fHistWt->Fill(mcwt,1.0); // histogram
1204  //******* Optional rejection ******
1205  if(fOptRej == 1) {
1206  Double_t random;
1207  random=fPseRan->Rndm();
1208  if( fMaxWtRej*random > fMCwt) goto ee0; // Wt=1 events, internal rejection
1209  if( fMCwt<fMaxWtRej ) {
1210  fMCwt = 1.0; // normal Wt=1 event
1211  } else {
1212  fMCwt = fMCwt/fMaxWtRej; // weight for overweighted events! kept for debug
1213  fSumOve += fMCwt-fMaxWtRej; // contribution of overweighted
1214  }
1215  }
1216  //********************** MC LOOP ENDS HERE **********************
1217 } // MakeEvent
1218 
1219 ////////////////////////////////////////////////////////////////////////////////
1220 /// User may get generated MC point/vector with help of this method
1221 
1223 {
1224  for ( Int_t k=0 ; k<fDim ; k++) *(MCvect +k) = fMCvect[k];
1225 }//GetMCvect
1226 
1227 ////////////////////////////////////////////////////////////////////////////////
1228 /// User may get weight MC weight using this method
1229 
1231 {
1232  return(fMCwt);
1233 }
1234 ////////////////////////////////////////////////////////////////////////////////
1235 /// User may get weight MC weight using this method
1236 
1238 {
1239  mcwt=fMCwt;
1240 }
1241 
1242 ////////////////////////////////////////////////////////////////////////////////
1243 /// User subprogram which generates MC event and returns MC weight
1244 
1246 {
1247  MakeEvent();
1248  GetMCvect(MCvect);
1249  return(fMCwt);
1250 }//MCgenerate
1251 
1252 ////////////////////////////////////////////////////////////////////////////////
1253 /// User subprogram.
1254 /// It provides the value of the integral calculated from the averages of the MC run
1255 /// May be called after (or during) the MC run.
1256 
1257 void TFoam::GetIntegMC(Double_t &mcResult, Double_t &mcError)
1258 {
1259  Double_t mCerelat;
1260  mcResult = 0.0;
1261  mCerelat = 1.0;
1262  if (fNevGen>0) {
1263  mcResult = fPrime*fSumWt/fNevGen;
1264  mCerelat = sqrt( fSumWt2/(fSumWt*fSumWt) - 1/fNevGen);
1265  }
1266  mcError = mcResult *mCerelat;
1267 }//GetIntegMC
1268 
1269 ////////////////////////////////////////////////////////////////////////////////
1270 /// User subprogram.
1271 /// It returns NORMALIZATION integral to be combined with the average weights
1272 /// and content of the histograms in order to get proper absolute normalization
1273 /// of the integrand and distributions.
1274 /// It can be called after initialization, before or during the MC run.
1275 
1276 void TFoam::GetIntNorm(Double_t& IntNorm, Double_t& Errel )
1277 {
1278  if(fOptRej == 1) { // Wt=1 events, internal rejection
1279  Double_t intMC,errMC;
1280  GetIntegMC(intMC,errMC);
1281  IntNorm = intMC;
1282  Errel = errMC;
1283  } else { // Wted events, NO internal rejection
1284  IntNorm = fPrime;
1285  Errel = 0;
1286  }
1287 }//GetIntNorm
1288 
1289 ////////////////////////////////////////////////////////////////////////////////
1290 /// May be called optionally after the MC run.
1291 /// Returns various parameters of the MC weight for efficiency evaluation
1292 
1294 {
1295  Double_t mCeff, wtLim;
1296  fMCMonit->GetMCeff(eps, mCeff, wtLim);
1297  wtMax = wtLim;
1298  aveWt = fSumWt/fNevGen;
1299  sigma = sqrt( fSumWt2/fNevGen -aveWt*aveWt );
1300 }//GetmCeff
1301 
1302 ////////////////////////////////////////////////////////////////////////////////
1303 /// May be called optionally by the user after the MC run.
1304 /// It provides normalization and also prints some information/statistics on the MC run.
1305 
1306 void TFoam::Finalize(Double_t& IntNorm, Double_t& Errel)
1307 {
1308  GetIntNorm(IntNorm,Errel);
1309  Double_t mcResult,mcError;
1310  GetIntegMC(mcResult,mcError);
1311  Double_t mCerelat= mcError/mcResult;
1312  //
1313  if(fChat>0) {
1314  Double_t eps = 0.0005;
1315  Double_t mCeff, mcEf2, wtMax, aveWt, sigma;
1316  GetWtParams(eps, aveWt, wtMax, sigma);
1317  mCeff=0;
1318  if(wtMax>0.0) mCeff=aveWt/wtMax;
1319  mcEf2 = sigma/aveWt;
1320  Double_t driver = fCells[0]->GetDriv();
1321  //
1322  BXOPE;
1323  BXTXT("****************************************");
1324  BXTXT("****** TFoam::Finalize ******");
1325  BXTXT("****************************************");
1326  BX1I(" NevGen",fNevGen, "Number of generated events in the MC generation ");
1327  BX1I(" nCalls",fNCalls, "Total number of function calls ");
1328  BXTXT("----------------------------------------");
1329  BX1F(" AveWt",aveWt, "Average MC weight ");
1330  BX1F(" WtMin",fWtMin, "Minimum MC weight (absolute) ");
1331  BX1F(" WtMax",fWtMax, "Maximum MC weight (absolute) ");
1332  BXTXT("----------------------------------------");
1333  BX1F(" XPrime",fPrime, "Primary total integral, R_prime ");
1334  BX1F(" XDiver",driver, "Driver total integral, R_loss ");
1335  BXTXT("----------------------------------------");
1336  BX2F(" IntMC", mcResult, mcError, "Result of the MC Integral");
1337  BX1F(" mCerelat", mCerelat, "Relative error of the MC integral ");
1338  BX1F(" <w>/WtMax",mCeff, "MC efficiency, acceptance rate");
1339  BX1F(" Sigma/<w>",mcEf2, "MC efficiency, variance/ave_wt");
1340  BX1F(" WtMax",wtMax, "WtMax(esp= 0.0005) ");
1341  BX1F(" Sigma",sigma, "variance of MC weight ");
1342  if(fOptRej==1) {
1343  Double_t avOve=fSumOve/fSumWt;
1344  BX1F("<OveW>/<W>",avOve, "Contrib. of events wt>MaxWtRej");
1345  }
1346  BXCLO;
1347  }
1348 } // Finalize
1349 
1350 ////////////////////////////////////////////////////////////////////////////////
1351 /// This can be called before Initialize, after setting kDim
1352 /// It defines which variables are excluded in the process of the cell division.
1353 /// For example 'FoamX->SetInhiDiv(1, 1);' inhibits division of y-variable.
1354 /// The resulting map of cells in 2-dim. case will look as follows:
1355 ///
1356 /// <img src="gif/foam_Map2.gif">
1357 
1358 void TFoam::SetInhiDiv(Int_t iDim, Int_t InhiDiv)
1359 {
1360 
1361 
1362  if(fDim==0) Error("TFoam","SetInhiDiv: fDim=0 \n");
1363  if(fInhiDiv == 0) {
1364  fInhiDiv = new Int_t[ fDim ];
1365  for(Int_t i=0; i<fDim; i++) fInhiDiv[i]=0;
1366  }
1367  //
1368  if( ( 0<=iDim) && (iDim<fDim)) {
1369  fInhiDiv[iDim] = InhiDiv;
1370  } else
1371  Error("SetInhiDiv:","Wrong iDim \n");
1372 }//SetInhiDiv
1373 
1374 ////////////////////////////////////////////////////////////////////////////////
1375 /// This should be called before Initialize, after setting kDim
1376 /// It predefines values of the cell division for certain variable iDim.
1377 /// For example setting 3 predefined division lines using:
1378 /// xDiv[0]=0.30; xDiv[1]=0.40; xDiv[2]=0.65;
1379 /// FoamX->SetXdivPRD(0,3,xDiv);
1380 /// results in the following 2-dim. pattern of the cells:
1381 /// <img src="gif/foam_Map3.gif">
1382 
1383 void TFoam::SetXdivPRD(Int_t iDim, Int_t len, Double_t xDiv[])
1384 {
1385  Int_t i;
1386 
1387  if(fDim<=0) Error("SetXdivPRD", "fDim=0 \n");
1388  if( len<1 ) Error("SetXdivPRD", "len<1 \n");
1389  // allocate list of pointers, if it was not done before
1390  if(fXdivPRD == 0) {
1391  fXdivPRD = new TFoamVect*[fDim];
1392  for(i=0; i<fDim; i++) fXdivPRD[i]=0;
1393  }
1394  // set division list for direction iDim in H-cubic space!!!
1395  if( ( 0<=iDim) && (iDim<fDim)) {
1396  fOptPRD =1; // !!!!
1397  if( fXdivPRD[iDim] != 0)
1398  Error("SetXdivPRD", "Second allocation of XdivPRD not allowed \n");
1399  fXdivPRD[iDim] = new TFoamVect(len); // allocate list of division points
1400  for(i=0; i<len; i++) {
1401  (*fXdivPRD[iDim])[i]=xDiv[i]; // set list of division points
1402  }
1403  } else {
1404  Error("SetXdivPRD", "Wrong iDim \n");
1405  }
1406  // Printing predefined division points
1407  std::cout<<" SetXdivPRD, idim= "<<iDim<<" len= "<<len<<" "<<std::endl;
1408  for(i=0; i<len; i++) {
1409  if (iDim < fDim) std::cout<< (*fXdivPRD[iDim])[i] <<" ";
1410  }
1411  std::cout<<std::endl;
1412  for(i=0; i<len; i++) std::cout<< xDiv[i] <<" ";
1413  std::cout<<std::endl;
1414  //
1415 }//SetXdivPRD
1416 
1417 ////////////////////////////////////////////////////////////////////////////////
1418 /// User utility, miscellaneous and debug.
1419 /// Checks all pointers in the tree of cells. This is useful auto diagnostic.
1420 /// level=0, no printout, failures causes STOP
1421 /// level=1, printout, failures lead to WARNINGS only
1422 
1424 {
1425  Int_t errors, warnings;
1426  TFoamCell *cell;
1427  Long_t iCell;
1428 
1429  errors = 0; warnings = 0;
1430  if (level==1) std::cout << "///////////////////////////// FOAM_Checks /////////////////////////////////" << std::endl;
1431  for(iCell=1; iCell<=fLastCe; iCell++) {
1432  cell = fCells[iCell];
1433  // checking general rules
1434  if( ((cell->GetDau0()==0) && (cell->GetDau1()!=0) ) ||
1435  ((cell->GetDau1()==0) && (cell->GetDau0()!=0) ) ) {
1436  errors++;
1437  if (level==1) Error("CheckAll","ERROR: Cell's no %ld has only one daughter \n",iCell);
1438  }
1439  if( (cell->GetDau0()==0) && (cell->GetDau1()==0) && (cell->GetStat()==0) ) {
1440  errors++;
1441  if (level==1) Error("CheckAll","ERROR: Cell's no %ld has no daughter and is inactive \n",iCell);
1442  }
1443  if( (cell->GetDau0()!=0) && (cell->GetDau1()!=0) && (cell->GetStat()==1) ) {
1444  errors++;
1445  if (level==1) Error("CheckAll","ERROR: Cell's no %ld has two daughters and is active \n",iCell);
1446  }
1447 
1448  // checking parents
1449  if( (cell->GetPare())!=fCells[0] ) { // not child of the root
1450  if ( (cell != cell->GetPare()->GetDau0()) && (cell != cell->GetPare()->GetDau1()) ) {
1451  errors++;
1452  if (level==1) Error("CheckAll","ERROR: Cell's no %ld parent not pointing to this cell\n ",iCell);
1453  }
1454  }
1455 
1456  // checking daughters
1457  if(cell->GetDau0()!=0) {
1458  if(cell != (cell->GetDau0())->GetPare()) {
1459  errors++;
1460  if (level==1) Error("CheckAll","ERROR: Cell's no %ld daughter 0 not pointing to this cell \n",iCell);
1461  }
1462  }
1463  if(cell->GetDau1()!=0) {
1464  if(cell != (cell->GetDau1())->GetPare()) {
1465  errors++;
1466  if (level==1) Error("CheckAll","ERROR: Cell's no %ld daughter 1 not pointing to this cell \n",iCell);
1467  }
1468  }
1469  }// loop after cells;
1470 
1471  // Check for empty cells
1472  for(iCell=0; iCell<=fLastCe; iCell++) {
1473  cell = fCells[iCell];
1474  if( (cell->GetStat()==1) && (cell->GetDriv()==0) ) {
1475  warnings++;
1476  if(level==1) Warning("CheckAll", "Warning: Cell no. %ld is active but empty \n", iCell);
1477  }
1478  }
1479  // summary
1480  if(level==1){
1481  Info("CheckAll","Check has found %d errors and %d warnings \n",errors, warnings);
1482  }
1483  if(errors>0){
1484  Info("CheckAll","Check - found total %d errors \n",errors);
1485  }
1486 } // Check
1487 
1488 ////////////////////////////////////////////////////////////////////////////////
1489 /// Prints geometry of ALL cells of the FOAM
1490 
1492 {
1493  Long_t iCell;
1494 
1495  for(iCell=0; iCell<=fLastCe; iCell++) {
1496  std::cout<<"Cell["<<iCell<<"]={ ";
1497  //std::cout<<" "<< fCells[iCell]<<" "; // extra DEBUG
1498  std::cout<<std::endl;
1499  fCells[iCell]->Print("1");
1500  std::cout<<"}"<<std::endl;
1501  }
1502 }
1503 
1504 ////////////////////////////////////////////////////////////////////////////////
1505 /// Debugging tool which plots 2-dimensional cells as rectangles
1506 /// in C++ format readable for root
1507 
1509 {
1510  std::ofstream outfile(filename, std::ios::out);
1511  Double_t x1,y1,x2,y2,x,y;
1512  Long_t iCell;
1513  Double_t offs =0.1;
1514  Double_t lpag =1-2*offs;
1515  outfile<<"{" << std::endl;
1516  outfile<<"cMap = new TCanvas(\"Map1\",\" Cell Map \",600,600);"<<std::endl;
1517  //
1518  outfile<<"TBox*a=new TBox();"<<std::endl;
1519  outfile<<"a->SetFillStyle(0);"<<std::endl; // big frame
1520  outfile<<"a->SetLineWidth(4);"<<std::endl;
1521  outfile<<"a->SetLineColor(2);"<<std::endl;
1522  outfile<<"a->DrawBox("<<offs<<","<<offs<<","<<(offs+lpag)<<","<<(offs+lpag)<<");"<<std::endl;
1523  //
1524  outfile<<"TText*t=new TText();"<<std::endl; // text for numbering
1525  outfile<<"t->SetTextColor(4);"<<std::endl;
1526  if(fLastCe<51)
1527  outfile<<"t->SetTextSize(0.025);"<<std::endl; // text for numbering
1528  else if(fLastCe<251)
1529  outfile<<"t->SetTextSize(0.015);"<<std::endl;
1530  else
1531  outfile<<"t->SetTextSize(0.008);"<<std::endl;
1532  //
1533  outfile<<"TBox*b=new TBox();"<<std::endl; // single cell
1534  outfile <<"b->SetFillStyle(0);"<<std::endl;
1535  //
1536  if(fDim==2 && fLastCe<=2000) {
1537  TFoamVect cellPosi(fDim); TFoamVect cellSize(fDim);
1538  outfile << "// =========== Rectangular cells ==========="<< std::endl;
1539  for(iCell=1; iCell<=fLastCe; iCell++) {
1540  if( fCells[iCell]->GetStat() == 1) {
1541  fCells[iCell]->GetHcub(cellPosi,cellSize);
1542  x1 = offs+lpag*( cellPosi[0]); y1 = offs+lpag*( cellPosi[1]);
1543  x2 = offs+lpag*(cellPosi[0]+cellSize[0]); y2 = offs+lpag*(cellPosi[1]+cellSize[1]);
1544  // cell rectangle
1545  if(fLastCe<=2000)
1546  outfile<<"b->DrawBox("<<x1<<","<<y1<<","<<x2<<","<<y2<<");"<<std::endl;
1547  // cell number
1548  if(fLastCe<=250) {
1549  x = offs+lpag*(cellPosi[0]+0.5*cellSize[0]); y = offs+lpag*(cellPosi[1]+0.5*cellSize[1]);
1550  outfile<<"t->DrawText("<<x<<","<<y<<","<<"\""<<iCell<<"\""<<");"<<std::endl;
1551  }
1552  }
1553  }
1554  outfile<<"// ============== End Rectangles ==========="<< std::endl;
1555  }//kDim=2
1556  //
1557  //
1558  outfile << "}" << std::endl;
1559  outfile.close();
1560 }
1561 
1563 {
1564 // Void function for backward compatibility
1565 
1566  Info("LinkCells", "VOID function for backward compatibility \n");
1567  return;
1568 }
1569 
TRefArray * fCellsAct
Definition: TFoam.h:59
Int_t fNCells
Definition: TFoam.h:37
Int_t fNoAct
Lists of division values encoded in one vector per direction.
Definition: TFoam.h:53
virtual void CheckAll(Int_t)
User utility, miscellaneous and debug.
Definition: TFoam.cxx:1423
virtual void RootPlot2dim(Char_t *)
Debugging tool which plots 2-dimensional cells as rectangles in C++ format readable for root...
Definition: TFoam.cxx:1508
Int_t fLastCe
Definition: TFoam.h:54
Int_t GetStat() const
Definition: TFoamCell.h:70
virtual Int_t Fill(Double_t x)
Increment bin with abscissa X by 1.
Definition: TH1.cxx:3125
Double_t * fMCvect
Definition: TFoam.h:65
static long int sum(long int i)
Definition: Factory.cxx:1786
An array of TObjects.
Definition: TObjArray.h:39
Double_t fSumWt
Definition: TFoam.h:75
virtual void Info(const char *method, const char *msgfmt,...) const
Issue info message.
Definition: TObject.cxx:899
#define BX2F(name, numb, err, text)
Definition: TFoam.cxx:139
void Add(TObject *obj)
Definition: TRefArray.h:84
virtual void LinkCells(void)
Definition: TFoam.cxx:1562
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:329
virtual void GetWtParams(Double_t, Double_t &, Double_t &, Double_t &)
May be called optionally after the MC run.
Definition: TFoam.cxx:1293
TFoamCell ** fCells
Definition: TFoam.h:55
Int_t fOptPRD
[fDim] Flags for inhibiting cell division
Definition: TFoam.h:50
void SetDau1(TFoamCell *Daug)
Definition: TFoamCell.h:76
virtual void RndmArray(Int_t n, Float_t *array)
Return an array of n random numbers uniformly distributed in ]0,1].
Definition: TRandom.cxx:548
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:1245
Double_t fMCwt
Definition: TFoam.h:66
Double_t * fRvec
Definition: TFoam.h:67
Int_t fNBin
Definition: TFoam.h:44
Double_t fMaxWtRej
Definition: TFoam.h:58
Double_t fNevGen
Definition: TFoam.h:77
TString fVersion
Definition: TFoam.h:34
Basic string class.
Definition: TString.h:137
static Bool_t AddDirectoryStatus()
Static function: cannot be inlined on Windows/NT.
Definition: TH1.cxx:700
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:170
int Int_t
Definition: RtypesCore.h:41
virtual void MakeEvent()
User subprogram.
Definition: TFoam.cxx:1172
bool Bool_t
Definition: RtypesCore.h:59
TFoamCell * GetDau1() const
Definition: TFoamCell.h:74
Int_t fChat
Definition: TFoam.h:41
const Bool_t kFALSE
Definition: Rtypes.h:92
virtual ~TFoam()
Default destructor.
Definition: TFoam.cxx:263
Double_t fSumWt2
Definition: TFoam.h:75
TFoamIntegrand * fRho
Definition: TFoam.h:69
An array of references to TObjects.
Definition: TRefArray.h:43
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:51
virtual void GenerCel2(TFoamCell *&)
Internal subprogram.
Definition: TFoam.cxx:1134
Double_t fMCerror
Definition: TFoam.h:81
virtual void MakeActiveList()
Internal subrogram used by Initialize.
Definition: TFoam.cxx:994
Double_t fWtMin
Definition: TFoam.h:78
Double_t fMCresult
Definition: TFoam.h:80
Short_t Abs(Short_t d)
Definition: TMathBase.h:110
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:1218
const char * Name
Definition: TXMLSetup.cxx:67
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:1383
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:1230
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:70
This is the base class for the ROOT Random number generators.
Definition: TRandom.h:31
virtual void PrintCells()
Prints geometry of ALL cells of the FOAM.
Definition: TFoam.cxx:1491
Long_t fNEffev
Definition: TFoam.h:74
Int_t fOptRej
Definition: TFoam.h:42
TObject * At(Int_t idx) const
Definition: TRefArray.h:184
virtual void SetRhoInt(void *Rho) R__DEPRECATED(6
User may use this to set pointer to the global function (not descending from TFoamIntegrand) serving ...
Definition: TFoam.cxx:1095
#define BX1I(name, numb, text)
Definition: TFoam.cxx:135
Double_t Sqr(Double_t x) const
Definition: TFoam.h:153
const Double_t sigma
TObjArray * fHistEdg
Definition: TFoam.h:61
virtual void Reset(Option_t *option="")
Reset.
Definition: TH1.cxx:9386
void CalcVolume()
Calculates volume of the cell using size params which are calculated.
Definition: TFoamCell.cxx:174
TFoamCell * GetDau0() const
Definition: TFoamCell.h:73
Double_t GetPrim() const
Definition: TFoamCell.h:65
virtual void Finalize(Double_t &, Double_t &)
May be called optionally by the user after the MC run.
Definition: TFoam.cxx:1306
virtual Double_t Rndm()
Machine independent random number generator.
Definition: TRandom.cxx:512
void SetXdiv(Double_t Xdiv)
Definition: TFoamCell.h:57
#define BXTXT(text)
Definition: TFoam.cxx:133
Double_t * fPrimAcu
Definition: TFoam.h:60
virtual void MakeAlpha()
Internal subrogram used by Initialize.
Definition: TFoam.cxx:879
void SetDriv(Double_t Driv)
Definition: TFoamCell.h:67
virtual void GetMCvect(Double_t *)
User may get generated MC point/vector with help of this method.
Definition: TFoam.cxx:1222
Double_t * fAlpha
Definition: TFoam.h:83
Int_t fEvPerBin
Definition: TFoam.h:46
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:1358
virtual Double_t Density(Int_t ndim, Double_t *)=0
void SetIntg(Double_t Intg)
Definition: TFoamCell.h:66
Double_t GetIntg() const
Definition: TFoamCell.h:63
Double_t GetVolume() const
Definition: TFoamCell.h:62
virtual void Error(const char *method, const char *msgfmt,...) const
Issue error message.
Definition: TObject.cxx:925
virtual void Initialize()
Basic initialization of FOAM invoked by the user.
Definition: TFoam.cxx:357
Int_t * fMaskDiv
Definition: TFoam.h:48
Double_t GetDriv() const
Definition: TFoamCell.h:64
TFoamCell * GetPare() const
Definition: TFoamCell.h:72
TString fName
Definition: TFoam.h:33
tomato 1-D histogram with a double per channel (see TH1 documentation)}
Definition: TH1.h:618
void SetStat(Int_t Stat)
Definition: TFoamCell.h:71
Double_t fWtMax
Definition: TFoam.h:78
Int_t fDim
Definition: TFoam.h:36
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
long Long_t
Definition: RtypesCore.h:50
Int_t fRNmax
Definition: TFoam.h:38
Int_t GetBest() const
Definition: TFoamCell.h:55
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:279
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:76
Double_t y[n]
Definition: legend1.C:17
double func(double *x, double *p)
Definition: stressTF1.cxx:213
#define BXOPE
Definition: TFoam.cxx:130
virtual void SetPseRan(TRandom *PseRan)
Definition: TFoam.h:125
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:40
void SetSerial(Int_t Serial)
Definition: TFoamCell.h:77
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:1276
virtual void Explore(TFoamCell *Cell)
Internal subprogram used by Initialize.
Definition: TFoam.cxx:554
TString fDate
Definition: TFoam.h:35
TObjArray * fHistDbg
Definition: TFoam.h:62
void Execute(const char *, const char *, int *=0)
Execute method on this object with the given parameter string, e.g.
Definition: TMethodCall.h:68
Short_t Max(Short_t a, Short_t b)
Definition: TMathBase.h:202
Double_t fPrime
Definition: TFoam.h:79
void SetBest(Int_t Best)
Definition: TFoamCell.h:56
void SetPrim(Double_t Prim)
Definition: TFoamCell.h:68
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:301
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:57
float type_of_call hi(const int &, const int &)
double result[121]
virtual Double_t Eval(Double_t *)
Internal subprogram.
Definition: TFoam.cxx:1112
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:63
Int_t * fInhiDiv
[fDim] Dynamic Mask for cell division
Definition: TFoam.h:49
void SetDau0(TFoamCell *Daug)
Definition: TFoamCell.h:75
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:73
virtual void GetIntegMC(Double_t &, Double_t &)
User subprogram.
Definition: TFoam.cxx:1257
virtual void Warning(const char *method, const char *msgfmt,...) const
Issue warning message.
Definition: TObject.cxx:911
Definition: TFoam.h:29
Int_t fNSampl
Definition: TFoam.h:45
const char * Data() const
Definition: TString.h:349
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:71