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