ROOT  6.06/09
Reference Guide
TGeoMaterial.cxx
Go to the documentation of this file.
1 // @(#)root/geom:$Id$
2 // Author: Andrei Gheata 25/10/01
3 
4 /*************************************************************************
5  * Copyright (C) 1995-2000, Rene Brun and Fons Rademakers. *
6  * All rights reserved. *
7  * *
8  * For the licensing terms see $ROOTSYS/LICENSE. *
9  * For the list of contributors see $ROOTSYS/README/CREDITS. *
10  *************************************************************************/
11 
12 ////////////////////////////////////////////////////////////////////////////////
13 //
14 //Begin_Html
15 /*
16 <img src="gif/t_material.jpg">
17 */
18 //End_Html
19 #include "Riostream.h"
20 #include "TMath.h"
21 #include "TObjArray.h"
22 #include "TStyle.h"
23 #include "TList.h"
24 #include "TGeoManager.h"
25 #include "TGeoExtension.h"
26 #include "TGeoMaterial.h"
27 
28 // statics and globals
29 
31 
32 ////////////////////////////////////////////////////////////////////////////////
33 /// Default constructor
34 
36  :TNamed(), TAttFill(),
37  fIndex(0),
38  fA(0.),
39  fZ(0.),
40  fDensity(0.),
41  fRadLen(0.),
42  fIntLen(0.),
43  fTemperature(0.),
44  fPressure(0.),
45  fState(kMatStateUndefined),
46  fShader(NULL),
47  fCerenkov(NULL),
48  fElement(NULL),
49  fUserExtension(0),
50  fFWExtension(0)
51 {
52  SetUsed(kFALSE);
53  fIndex = -1;
54  fTemperature = STP_temperature;
55  fPressure = STP_pressure;
56  fState = kMatStateUndefined;
57 }
58 
59 ////////////////////////////////////////////////////////////////////////////////
60 /// constructor
61 
63  :TNamed(name, ""), TAttFill(),
64  fIndex(0),
65  fA(0.),
66  fZ(0.),
67  fDensity(0.),
68  fRadLen(0.),
69  fIntLen(0.),
70  fTemperature(0.),
71  fPressure(0.),
72  fState(kMatStateUndefined),
73  fShader(NULL),
74  fCerenkov(NULL),
75  fElement(NULL),
76  fUserExtension(0),
77  fFWExtension(0)
78 {
79  fName = fName.Strip();
80  SetUsed(kFALSE);
81  fIndex = -1;
85 
86  if (!gGeoManager) {
87  gGeoManager = new TGeoManager("Geometry", "default geometry");
88  }
89  gGeoManager->AddMaterial(this);
90 }
91 
92 ////////////////////////////////////////////////////////////////////////////////
93 /// constructor
94 
96  Double_t rho, Double_t radlen, Double_t intlen)
97  :TNamed(name, ""), TAttFill(),
98  fIndex(0),
99  fA(a),
100  fZ(z),
101  fDensity(rho),
102  fRadLen(0.),
103  fIntLen(0.),
104  fTemperature(0.),
105  fPressure(0.),
106  fState(kMatStateUndefined),
107  fShader(NULL),
108  fCerenkov(NULL),
109  fElement(NULL),
110  fUserExtension(0),
111  fFWExtension(0)
112 {
113  fName = fName.Strip();
114  SetUsed(kFALSE);
115  fIndex = -1;
116  fA = a;
117  fZ = z;
118  fDensity = rho;
122  SetRadLen(radlen, intlen);
123  if (!gGeoManager) {
124  gGeoManager = new TGeoManager("Geometry", "default geometry");
125  }
126  if (fZ - Int_t(fZ) > 1E-3)
127  Warning("ctor", "Material %s defined with fractional Z=%f", GetName(), fZ);
128  if (GetElement()) GetElement()->SetUsed();
129  gGeoManager->AddMaterial(this);
130 }
131 
132 ////////////////////////////////////////////////////////////////////////////////
133 /// Constructor with state, temperature and pressure.
134 
136  EGeoMaterialState state, Double_t temperature, Double_t pressure)
137  :TNamed(name, ""), TAttFill(),
138  fIndex(0),
139  fA(a),
140  fZ(z),
141  fDensity(rho),
142  fRadLen(0.),
143  fIntLen(0.),
144  fTemperature(temperature),
145  fPressure(pressure),
146  fState(state),
147  fShader(NULL),
148  fCerenkov(NULL),
149  fElement(NULL),
150  fUserExtension(0),
151  fFWExtension(0)
152 {
153  fName = fName.Strip();
154  SetUsed(kFALSE);
155  fIndex = -1;
156  SetRadLen(0,0);
157  if (!gGeoManager) {
158  gGeoManager = new TGeoManager("Geometry", "default geometry");
159  }
160  if (fZ - Int_t(fZ) > 1E-3)
161  Warning("ctor", "Material %s defined with fractional Z=%f", GetName(), fZ);
162  if (GetElement()) GetElement()->SetUsed();
163  gGeoManager->AddMaterial(this);
164 }
165 
166 ////////////////////////////////////////////////////////////////////////////////
167 /// constructor
168 
170  :TNamed(name, ""), TAttFill(),
171  fIndex(0),
172  fA(0.),
173  fZ(0.),
174  fDensity(rho),
175  fRadLen(0.),
176  fIntLen(0.),
177  fTemperature(0.),
178  fPressure(0.),
179  fState(kMatStateUndefined),
180  fShader(NULL),
181  fCerenkov(NULL),
182  fElement(elem),
183  fUserExtension(0),
184  fFWExtension(0)
185 {
186  fName = fName.Strip();
187  SetUsed(kFALSE);
188  fIndex = -1;
189  fA = elem->A();
190  fZ = elem->Z();
191  SetRadLen(0,0);
195  if (!gGeoManager) {
196  gGeoManager = new TGeoManager("Geometry", "default geometry");
197  }
198  if (fZ - Int_t(fZ) > 1E-3)
199  Warning("ctor", "Material %s defined with fractional Z=%f", GetName(), fZ);
200  if (GetElement()) GetElement()->SetUsed();
201  gGeoManager->AddMaterial(this);
202 }
203 
204 ////////////////////////////////////////////////////////////////////////////////
205 
207  TNamed(gm),
208  TAttFill(gm),
209  fIndex(gm.fIndex),
210  fA(gm.fA),
211  fZ(gm.fZ),
212  fDensity(gm.fDensity),
213  fRadLen(gm.fRadLen),
214  fIntLen(gm.fIntLen),
215  fTemperature(gm.fTemperature),
216  fPressure(gm.fPressure),
217  fState(gm.fState),
218  fShader(gm.fShader),
219  fCerenkov(gm.fCerenkov),
220  fElement(gm.fElement),
221  fUserExtension(gm.fUserExtension->Grab()),
222  fFWExtension(gm.fFWExtension->Grab())
223 
224 {
225  //copy constructor
226 }
227 
228 ////////////////////////////////////////////////////////////////////////////////
229 ///assignment operator
230 
232 {
233  if(this!=&gm) {
234  TNamed::operator=(gm);
236  fIndex=gm.fIndex;
237  fA=gm.fA;
238  fZ=gm.fZ;
239  fDensity=gm.fDensity;
240  fRadLen=gm.fRadLen;
241  fIntLen=gm.fIntLen;
243  fPressure=gm.fPressure;
244  fState=gm.fState;
245  fShader=gm.fShader;
246  fCerenkov=gm.fCerenkov;
247  fElement=gm.fElement;
250  }
251  return *this;
252 }
253 
254 ////////////////////////////////////////////////////////////////////////////////
255 /// Destructor
256 
258 {
261 }
262 
263 ////////////////////////////////////////////////////////////////////////////////
264 /// Connect user-defined extension to the material. The material "grabs" a copy, so
265 /// the original object can be released by the producer. Release the previously
266 /// connected extension if any.
267 ///==========================================================================
268 /// NOTE: This interface is intended for user extensions and is guaranteed not
269 /// to be used by TGeo
270 ///==========================================================================
271 
273 {
275  fUserExtension = 0;
276  if (ext) fUserExtension = ext->Grab();
277 }
278 
279 ////////////////////////////////////////////////////////////////////////////////
280 /// Connect framework defined extension to the material. The material "grabs" a copy,
281 /// so the original object can be released by the producer. Release the previously
282 /// connected extension if any.
283 ///==========================================================================
284 /// NOTE: This interface is intended for the use by TGeo and the users should
285 /// NOT connect extensions using this method
286 ///==========================================================================
287 
289 {
291  fFWExtension = 0;
292  if (ext) fFWExtension = ext->Grab();
293 }
294 
295 ////////////////////////////////////////////////////////////////////////////////
296 /// Get a copy of the user extension pointer. The user must call Release() on
297 /// the copy pointer once this pointer is not needed anymore (equivalent to
298 /// delete() after calling new())
299 
301 {
302  if (fUserExtension) return fUserExtension->Grab();
303  return 0;
304 }
305 
306 ////////////////////////////////////////////////////////////////////////////////
307 /// Get a copy of the framework extension pointer. The user must call Release() on
308 /// the copy pointer once this pointer is not needed anymore (equivalent to
309 /// delete() after calling new())
310 
312 {
313  if (fFWExtension) return fFWExtension->Grab();
314  return 0;
315 }
316 
317 ////////////////////////////////////////////////////////////////////////////////
318 /// Provide a pointer name containing uid.
319 
321 {
322  static TString name;
323  name = TString::Format("pMat%d", GetUniqueID());
324  return (char*)name.Data();
325 }
326 
327 ////////////////////////////////////////////////////////////////////////////////
328 /// Set radiation/absorbtion lengths. If the values are negative, their absolute value
329 /// is taken, otherwise radlen is recomputed using G3 formula.
330 
332 {
333  fRadLen = TMath::Abs(radlen);
334  fIntLen = TMath::Abs(intlen);
335  // Check for vacuum
336  if (fA<0.9 || fZ<0.9) {
337  if (radlen<-1e5 || intlen<-1e-5) {
338  Error("SetRadLen","Material %s: user values taken for vacuum: radlen=%g or intlen=%g - too small", GetName(),fRadLen, fIntLen);
339  return;
340  }
341  // Ignore positive values and take big numbers
342  if (radlen>=0) fRadLen = 1.E30;
343  if (intlen>=0) fIntLen = 1.E30;
344  return;
345  }
346  // compute radlen systematically with G3 formula for a valid material
347  if (radlen>=0) {
348  //taken grom Geant3 routine GSMATE
349  const Double_t alr2av=1.39621E-03, al183=5.20948;
351  (al183-TMath::Log(fZ)/3-TGeoMaterial::Coulomb(fZ)));
352  }
353  // Compute interaction length using the same formula as in GEANT4
354  if (intlen>=0) {
355  const Double_t cm = 1.;
356  const Double_t g = 6.2415e21; // [gram = 1E-3*joule*s*s/(m*m)]
357  const Double_t amu = 1.03642688246781065e-02; // [MeV/c^2]
358  const Double_t lambda0 = 35.*g/(cm*cm); // [g/cm^2]
359  Double_t nilinv = 0.0;
360  TGeoElement *elem = GetElement();
361  if (!elem) {
362  Fatal("SetRadLen", "Element not found for material %s", GetName());
363  return;
364  }
365  Double_t nbAtomsPerVolume = TMath::Na()*fDensity/elem->A();
366  nilinv += nbAtomsPerVolume*TMath::Power(elem->Neff(), 0.6666667);
367  nilinv *= amu/lambda0;
368  fIntLen = (nilinv<=0) ? TGeoShape::Big() : (1./nilinv);
369  }
370 }
371 
372 ////////////////////////////////////////////////////////////////////////////////
373 /// static function
374 /// Compute Coulomb correction for pair production and Brem
375 /// REFERENCE : EGS MANUAL SLAC 210 - UC32 - JUNE 78
376 /// FORMULA 2.7.17
377 
379 {
380  const Double_t alpha = 7.29927E-03;
381 
382  Double_t az = alpha*z;
383  Double_t az2 = az*az;
384  Double_t az4 = az2 * az2;
385  Double_t fp = ( 0.0083*az4 + 0.20206 + 1./(1.+az2) ) * az2;
386  Double_t fm = ( 0.0020*az4 + 0.0369 ) * az4;
387  return fp - fm;
388 }
389 
390 ////////////////////////////////////////////////////////////////////////////////
391 /// return true if the other material has the same physical properties
392 
394 {
395  if (other==this) return kTRUE;
396  if (other->IsMixture()) return kFALSE;
397  if (TMath::Abs(fA-other->GetA())>1E-3) return kFALSE;
398  if (TMath::Abs(fZ-other->GetZ())>1E-3) return kFALSE;
399  if (TMath::Abs(fDensity-other->GetDensity())>1E-6) return kFALSE;
400  if (GetCerenkovProperties() != other->GetCerenkovProperties()) return kFALSE;
401 // if (fRadLen != other->GetRadLen()) return kFALSE;
402 // if (fIntLen != other->GetIntLen()) return kFALSE;
403  return kTRUE;
404 }
405 
406 ////////////////////////////////////////////////////////////////////////////////
407 /// print characteristics of this material
408 
409 void TGeoMaterial::Print(const Option_t * /*option*/) const
410 {
411  printf("Material %s %s A=%g Z=%g rho=%g radlen=%g intlen=%g index=%i\n", GetName(), GetTitle(),
413 }
414 
415 ////////////////////////////////////////////////////////////////////////////////
416 /// Save a primitive as a C++ statement(s) on output stream "out".
417 
418 void TGeoMaterial::SavePrimitive(std::ostream &out, Option_t * /*option*/ /*= ""*/)
419 {
421  char *name = GetPointerName();
422  out << "// Material: " << GetName() << std::endl;
423  out << " a = " << fA << ";" << std::endl;
424  out << " z = " << fZ << ";" << std::endl;
425  out << " density = " << fDensity << ";" << std::endl;
426  out << " radl = " << fRadLen << ";" << std::endl;
427  out << " absl = " << fIntLen << ";" << std::endl;
428 
429  out << " " << name << " = new TGeoMaterial(\"" << GetName() << "\", a,z,density,radl,absl);" << std::endl;
430  out << " " << name << "->SetIndex(" << GetIndex() << ");" << std::endl;
432 }
433 
434 ////////////////////////////////////////////////////////////////////////////////
435 /// Get some default color related to this material.
436 
438 {
439  Int_t id = 1+ gGeoManager->GetListOfMaterials()->IndexOf(this);
440  return (2+id%6);
441 }
442 
443 ////////////////////////////////////////////////////////////////////////////////
444 /// Get a pointer to the element this material is made of.
445 
447 {
448  if (fElement) return fElement;
450  return table->GetElement(Int_t(fZ));
451 }
452 ////////////////////////////////////////////////////////////////////////////////
453 /// Single interface to get element properties.
454 
456 {
457  a = fA;
458  z = fZ;
459  w = 1.;
460 }
461 
462 ////////////////////////////////////////////////////////////////////////////////
463 /// Retreive material index in the list of materials
464 
466 {
467  if (fIndex>=0) return fIndex;
468  TList *matlist = gGeoManager->GetListOfMaterials();
469  fIndex = matlist->IndexOf(this);
470  return fIndex;
471 }
472 
473 ////////////////////////////////////////////////////////////////////////////////
474 /// Create the material representing the decay product of this material at a
475 /// given time. The precision represent the minimum cumulative branching ratio for
476 /// which decay products are still taken into account.
477 
479 {
480  TObjArray *pop = new TObjArray();
481  if (!fElement || !fElement->IsRadioNuclide()) return this;
482  FillMaterialEvolution(pop, precision);
483  Int_t ncomp = pop->GetEntriesFast();
484  if (!ncomp) return this;
485  TGeoElementRN *el;
486  Double_t *weight = new Double_t[ncomp];
487  Double_t amed = 0.;
488  Int_t i;
489  for (i=0; i<ncomp; i++) {
490  el = (TGeoElementRN *)pop->At(i);
491  weight[i] = el->Ratio()->Concentration(time) * el->A();
492  amed += weight[i];
493  }
494  Double_t rho = fDensity*amed/fA;
495  TGeoMixture *mix = 0;
496  Int_t ncomp1 = ncomp;
497  for (i=0; i<ncomp; i++) {
498  if ((weight[i]/amed)<precision) {
499  amed -= weight[i];
500  ncomp1--;
501  }
502  }
503  if (ncomp1<2) {
504  el = (TGeoElementRN *)pop->At(0);
505  delete [] weight;
506  delete pop;
507  if (ncomp1==1) return new TGeoMaterial(TString::Format("%s-evol",GetName()), el, rho);
508  return NULL;
509  }
510  mix = new TGeoMixture(TString::Format("%s-evol",GetName()), ncomp, rho);
511  for (i=0; i<ncomp; i++) {
512  weight[i] /= amed;
513  if (weight[i]<precision) continue;
514  el = (TGeoElementRN *)pop->At(i);
515  mix->AddElement(el, weight[i]);
516  }
517  delete [] weight;
518  delete pop;
519  return mix;
520 }
521 
522 ////////////////////////////////////////////////////////////////////////////////
523 /// Fills a user array with all the elements deriving from the possible
524 /// decay of the top element composing the mixture. Each element contained
525 /// by <population> may be a radionuclide having a Bateman solution attached.
526 /// The precision represent the minimum cumulative branching ratio for
527 /// which decay products are still taken into account.
528 /// To visualize the time evolution of each decay product one can use:
529 /// TGeoElement *elem = population->At(index);
530 /// TGeoElementRN *elemrn = 0;
531 /// if (elem->IsRadioNuclide()) elemrn = (TGeoElementRN*)elem;
532 /// One can get Ni/N1(t=0) at any moment of time. Ni is the number of atoms
533 /// of one of the decay products, N1(0) is the number of atoms of the top
534 /// element at t=0.
535 /// Double_t fraction_weight = elemrn->Ratio()->Concentration(time);
536 /// One can also display the time evolution of the fractional weigth:
537 /// elemrn->Ratio()->Draw(option);
538 
540 {
541  if (population->GetEntriesFast()) {
542  Error("FillMaterialEvolution", "Provide an empty array !");
543  return;
544  }
546  TGeoElement *elem;
547  TGeoElementRN *elemrn;
548  TIter next(table->GetElementsRN());
549  while ((elemrn=(TGeoElementRN*)next())) elemrn->ResetRatio();
550  elem = GetElement();
551  if (!elem) {
552  Fatal("FillMaterialEvolution", "Element not found for material %s", GetName());
553  return;
554  }
555  if (!elem->IsRadioNuclide()) {
556  population->Add(elem);
557  return;
558  }
559  elemrn = (TGeoElementRN*)elem;
560  elemrn->FillPopulation(population, precision);
561 }
562 
563 
564 /*************************************************************************
565  * TGeoMixture - mixtures of elements
566  *
567  *************************************************************************/
569 
570 ////////////////////////////////////////////////////////////////////////////////
571 /// Default constructor
572 
574 {
575  fNelements = 0;
576  fZmixture = 0;
577  fAmixture = 0;
578  fWeights = 0;
579  fNatoms = 0;
580  fElements = 0;
581 }
582 
583 ////////////////////////////////////////////////////////////////////////////////
584 /// constructor
585 
586 TGeoMixture::TGeoMixture(const char *name, Int_t /*nel*/, Double_t rho)
587  :TGeoMaterial(name)
588 {
589  fZmixture = 0;
590  fAmixture = 0;
591  fWeights = 0;
592  fNelements = 0;
593  fNatoms = 0;
594  fDensity = rho;
595  fElements = 0;
596  if (fDensity < 0) fDensity = 0.001;
597 }
598 
599 ////////////////////////////////////////////////////////////////////////////////
600 ///copy constructor
601 
603  TGeoMaterial(gm),
604  fNelements(gm.fNelements),
605  fZmixture(gm.fZmixture),
606  fAmixture(gm.fAmixture),
607  fWeights(gm.fWeights),
608  fNatoms(gm.fNatoms),
609  fElements(gm.fElements)
610 {
611 }
612 
613 ////////////////////////////////////////////////////////////////////////////////
614 ///assignment operator
615 
617 {
618  if(this!=&gm) {
621  fZmixture=gm.fZmixture;
622  fAmixture=gm.fAmixture;
623  fWeights=gm.fWeights;
624  fNatoms = gm.fNatoms;
625  fElements = gm.fElements;
626  }
627  return *this;
628 }
629 
630 ////////////////////////////////////////////////////////////////////////////////
631 /// Destructor
632 
634 {
635  if (fZmixture) delete[] fZmixture;
636  if (fAmixture) delete[] fAmixture;
637  if (fWeights) delete[] fWeights;
638  if (fNatoms) delete[] fNatoms;
639  if (fElements) delete fElements;
640 }
641 
642 ////////////////////////////////////////////////////////////////////////////////
643 /// Compute effective A/Z and radiation length
644 
646 {
647  const Double_t alr2av = 1.39621E-03 , al183 =5.20948;
648  const Double_t cm = 1.;
649  const Double_t g = 6.2415e21; // [gram = 1E-3*joule*s*s/(m*m)]
650  const Double_t amu = 1.03642688246781065e-02; // [MeV/c^2]
651  const Double_t lambda0 = 35.*g/(cm*cm); // [g/cm^2]
652  Double_t radinv = 0.0;
653  Double_t nilinv = 0.0;
654  Double_t nbAtomsPerVolume;
655  fA = 0;
656  fZ = 0;
657  for (Int_t j=0;j<fNelements;j++) {
658  if (fWeights[j] <= 0) continue;
659  fA += fWeights[j]*fAmixture[j];
660  fZ += fWeights[j]*fZmixture[j];
661  nbAtomsPerVolume = TMath::Na()*fDensity*fWeights[j]/GetElement(j)->A();
662  nilinv += nbAtomsPerVolume*TMath::Power(GetElement(j)->Neff(), 0.6666667);
663  Double_t zc = fZmixture[j];
664  Double_t alz = TMath::Log(zc)/3.;
665  Double_t xinv = zc*(zc+TGeoMaterial::ScreenFactor(zc))*
666  (al183-alz-TGeoMaterial::Coulomb(zc))/fAmixture[j];
667  radinv += xinv*fWeights[j];
668  }
669  radinv *= alr2av*fDensity;
670  if (radinv > 0) fRadLen = 1/radinv;
671  // Compute interaction length
672  nilinv *= amu/lambda0;
673  fIntLen = (nilinv<=0) ? TGeoShape::Big() : (1./nilinv);
674 }
675 
676 ////////////////////////////////////////////////////////////////////////////////
677 /// add an element to the mixture using fraction by weight
678 /// Check if the element is already defined
679 
681 {
683  if (z<1 || z>table->GetNelements()-1)
684  Fatal("AddElement", "Cannot add element having Z=%d to mixture %s", (Int_t)z, GetName());
685  Int_t i;
686  for (i=0; i<fNelements; i++) {
687  if (TMath::Abs(z-fZmixture[i])<1.e-6 && TMath::Abs(a-fAmixture[i])<1.e-6) {
688  fWeights[i] += weight;
690  return;
691  }
692  }
693  if (!fNelements) {
694  fZmixture = new Double_t[1];
695  fAmixture = new Double_t[1];
696  fWeights = new Double_t[1];
697  } else {
698  Int_t nelements = fNelements+1;
699  Double_t *zmixture = new Double_t[nelements];
700  Double_t *amixture = new Double_t[nelements];
701  Double_t *weights = new Double_t[nelements];
702  for (Int_t j=0; j<fNelements; j++) {
703  zmixture[j] = fZmixture[j];
704  amixture[j] = fAmixture[j];
705  weights[j] = fWeights[j];
706  }
707  delete [] fZmixture;
708  delete [] fAmixture;
709  delete [] fWeights;
710  fZmixture = zmixture;
711  fAmixture = amixture;
712  fWeights = weights;
713  }
714 
715  fNelements++;
716  i = fNelements - 1;
717  fZmixture[i] = z;
718  fAmixture[i] = a;
719  fWeights[i] = weight;
720  if (z - Int_t(z) > 1E-3)
721  Warning("DefineElement", "Mixture %s has element defined with fractional Z=%f", GetName(), z);
722  GetElement(i)->SetDefined();
723  table->GetElement((Int_t)z)->SetDefined();
724 
725  //compute equivalent radiation length (taken from Geant3/GSMIXT)
727 }
728 
729 ////////////////////////////////////////////////////////////////////////////////
730 /// Define one component of the mixture as an existing material/mixture.
731 
733 {
734  TGeoElement *elnew, *elem;
735  Double_t a,z;
736  if (!mat->IsMixture()) {
737  elem = mat->GetBaseElement();
738  if (elem) {
739  AddElement(elem, weight);
740  } else {
741  a = mat->GetA();
742  z = mat->GetZ();
743  AddElement(a, z, weight);
744  }
745  return;
746  }
747  // The material is a mixture.
748  TGeoMixture *mix = (TGeoMixture*)mat;
749  Double_t wnew;
750  Int_t nelem = mix->GetNelements();
751  Bool_t elfound;
752  Int_t i,j;
753  // loop the elements of the daughter mixture
754  for (i=0; i<nelem; i++) {
755  elfound = kFALSE;
756  elnew = mix->GetElement(i);
757  if (!elnew) continue;
758  // check if we have the element already defined in the parent mixture
759  for (j=0; j<fNelements; j++) {
760  if (fWeights[j]<=0) continue;
761  elem = GetElement(j);
762  if (elem == elnew) {
763  // element found, compute new weight
764  fWeights[j] += weight * (mix->GetWmixt())[i];
765  elfound = kTRUE;
766  break;
767  }
768  }
769  if (elfound) continue;
770  // element not found, define it
771  wnew = weight * (mix->GetWmixt())[i];
772  AddElement(elnew, wnew);
773  }
774 }
775 
776 ////////////////////////////////////////////////////////////////////////////////
777 /// add an element to the mixture using fraction by weight
778 
780 {
781  TGeoElement *elemold;
783  if (!fElements) fElements = new TObjArray(128);
784  Bool_t exist = kFALSE;
785  // If previous elements were defined by A/Z, add corresponding TGeoElements
786  for (Int_t i=0; i<fNelements; i++) {
787  elemold = (TGeoElement*)fElements->At(i);
788  if (!elemold) fElements->AddAt(elemold = table->GetElement((Int_t)fZmixture[i]), i);
789  if (elemold == elem) exist = kTRUE;
790  }
791  if (!exist) fElements->AddAtAndExpand(elem, fNelements);
792  AddElement(elem->A(), elem->Z(), weight);
793 }
794 
795 ////////////////////////////////////////////////////////////////////////////////
796 /// Add a mixture element by number of atoms in the chemical formula.
797 
799 {
800  Int_t i,j;
801  Double_t amol;
802  TGeoElement *elemold;
804  if (!fElements) fElements = new TObjArray(128);
805  // Check if the element is already defined
806  for (i=0; i<fNelements; i++) {
807  elemold = (TGeoElement*)fElements->At(i);
808  if (!elemold) fElements->AddAt(table->GetElement((Int_t)fZmixture[i]), i);
809  else if (elemold != elem) continue;
810  if ((elem==elemold) ||
811  (TMath::Abs(elem->Z()-fZmixture[i])<1.e-6 && TMath::Abs(elem->A()-fAmixture[i])<1.e-6)) {
812  fNatoms[i] += natoms;
813  amol = 0.;
814  for (j=0; j<fNelements; j++) amol += fAmixture[j]*fNatoms[j];
815  for (j=0; j<fNelements; j++) fWeights[j] = fNatoms[j]*fAmixture[j]/amol;
817  return;
818  }
819  }
820  // New element
821  if (!fNelements) {
822  fZmixture = new Double_t[1];
823  fAmixture = new Double_t[1];
824  fWeights = new Double_t[1];
825  fNatoms = new Int_t[1];
826  } else {
827  if (!fNatoms) {
828  Fatal("AddElement", "Cannot add element by natoms in mixture %s after defining elements by weight",
829  GetName());
830  return;
831  }
832  Int_t nelements = fNelements+1;
833  Double_t *zmixture = new Double_t[nelements];
834  Double_t *amixture = new Double_t[nelements];
835  Double_t *weights = new Double_t[nelements];
836  Int_t *nnatoms = new Int_t[nelements];
837  for (j=0; j<fNelements; j++) {
838  zmixture[j] = fZmixture[j];
839  amixture[j] = fAmixture[j];
840  weights[j] = fWeights[j];
841  nnatoms[j] = fNatoms[j];
842  }
843  delete [] fZmixture;
844  delete [] fAmixture;
845  delete [] fWeights;
846  delete [] fNatoms;
847  fZmixture = zmixture;
848  fAmixture = amixture;
849  fWeights = weights;
850  fNatoms = nnatoms;
851  }
852  fNelements++;
853  Int_t iel = fNelements-1;
854  fZmixture[iel] = elem->Z();
855  fAmixture[iel] = elem->A();
856  fNatoms[iel] = natoms;
857  fElements->AddAtAndExpand(elem, iel);
858  amol = 0.;
859  for (i=0; i<fNelements; i++) {
860  if (fNatoms[i]<=0) return;
861  amol += fAmixture[i]*fNatoms[i];
862  }
863  for (i=0; i<fNelements; i++) fWeights[i] = fNatoms[i]*fAmixture[i]/amol;
864  table->GetElement(elem->Z())->SetDefined();
866 }
867 
868 ////////////////////////////////////////////////////////////////////////////////
869 /// Define the mixture element at index iel by number of atoms in the chemical formula.
870 
871 void TGeoMixture::DefineElement(Int_t /*iel*/, Int_t z, Int_t natoms)
872 {
874  TGeoElement *elem = table->GetElement(z);
875  if (!elem) {
876  Fatal("DefineElement", "In mixture %s, element with Z=%i not found",GetName(),z);
877  return;
878  }
879  AddElement(elem, natoms);
880 }
881 
882 ////////////////////////////////////////////////////////////////////////////////
883 /// Retreive the pointer to the element corresponding to component I.
884 
886 {
887  if (i<0 || i>=fNelements) {
888  Error("GetElement", "Mixture %s has only %d elements", GetName(), fNelements);
889  return 0;
890  }
891  TGeoElement *elem = 0;
892  if (fElements) elem = (TGeoElement*)fElements->At(i);
893  if (elem) return elem;
895  return table->GetElement(Int_t(fZmixture[i]));
896 }
897 
898 ////////////////////////////////////////////////////////////////////////////////
899 /// Get specific activity (in Bq/gram) for the whole mixture (no argument) or
900 /// for a given component.
901 
903 {
904  if (i>=0 && i<fNelements) return fWeights[i]*GetElement(i)->GetSpecificActivity();
905  Double_t sa = 0;
906  for (Int_t iel=0; iel<fNelements; iel++) {
907  sa += fWeights[iel]*GetElement(iel)->GetSpecificActivity();
908  }
909  return sa;
910 }
911 
912 ////////////////////////////////////////////////////////////////////////////////
913 /// Return true if the other material has the same physical properties
914 
916 {
917  if (other->IsEqual(this)) return kTRUE;
918  if (!other->IsMixture()) return kFALSE;
919  TGeoMixture *mix = (TGeoMixture*)other;
920  if (!mix) return kFALSE;
921  if (fNelements != mix->GetNelements()) return kFALSE;
922  if (TMath::Abs(fA-other->GetA())>1E-3) return kFALSE;
923  if (TMath::Abs(fZ-other->GetZ())>1E-3) return kFALSE;
924  if (TMath::Abs(fDensity-other->GetDensity())>1E-6) return kFALSE;
925  if (GetCerenkovProperties() != other->GetCerenkovProperties()) return kFALSE;
926 // if (fRadLen != other->GetRadLen()) return kFALSE;
927 // if (fIntLen != other->GetIntLen()) return kFALSE;
928  for (Int_t i=0; i<fNelements; i++) {
929  if (TMath::Abs(fZmixture[i]-(mix->GetZmixt())[i])>1E-3) return kFALSE;
930  if (TMath::Abs(fAmixture[i]-(mix->GetAmixt())[i])>1E-3) return kFALSE;
931  if (TMath::Abs(fWeights[i]-(mix->GetWmixt())[i])>1E-3) return kFALSE;
932  }
933  return kTRUE;
934 }
935 
936 ////////////////////////////////////////////////////////////////////////////////
937 /// print characteristics of this material
938 
939 void TGeoMixture::Print(const Option_t * /*option*/) const
940 {
941  printf("Mixture %s %s Aeff=%g Zeff=%g rho=%g radlen=%g intlen=%g index=%i\n", GetName(), GetTitle(),
943  for (Int_t i=0; i<fNelements; i++) {
944  if (fNatoms) printf(" Element #%i : %s Z=%6.2f A=%6.2f w=%6.3f natoms=%d\n", i, GetElement(i)->GetName(),fZmixture[i],
945  fAmixture[i], fWeights[i], fNatoms[i]);
946  else printf(" Element #%i : %s Z=%6.2f A=%6.2f w=%6.3f\n", i, GetElement(i)->GetName(),fZmixture[i],
947  fAmixture[i], fWeights[i]);
948  }
949 }
950 
951 ////////////////////////////////////////////////////////////////////////////////
952 /// Save a primitive as a C++ statement(s) on output stream "out".
953 
954 void TGeoMixture::SavePrimitive(std::ostream &out, Option_t * /*option*/ /*= ""*/)
955 {
957  char *name = GetPointerName();
958  out << "// Mixture: " << GetName() << std::endl;
959  out << " nel = " << fNelements << ";" << std::endl;
960  out << " density = " << fDensity << ";" << std::endl;
961  out << " " << name << " = new TGeoMixture(\"" << GetName() << "\", nel,density);" << std::endl;
962  for (Int_t i=0; i<fNelements; i++) {
963  TGeoElement *el = GetElement(i);
964  out << " a = " << fAmixture[i] << "; z = "<< fZmixture[i] << "; w = " << fWeights[i] << "; // " << el->GetName() << std::endl;
965  out << " " << name << "->DefineElement(" << i << ",a,z,w);" << std::endl;
966  }
967  out << " " << name << "->SetIndex(" << GetIndex() << ");" << std::endl;
969 }
970 
971 ////////////////////////////////////////////////////////////////////////////////
972 /// Create the mixture representing the decay product of this material at a
973 /// given time. The precision represent the minimum cumulative branching ratio for
974 /// which decay products are still taken into account.
975 
977 {
978  TObjArray *pop = new TObjArray();
979  FillMaterialEvolution(pop, precision);
980  Int_t ncomp = pop->GetEntriesFast();
981  if (!ncomp) return this;
982  TGeoElement *elem;
983  TGeoElementRN *el;
984  Double_t *weight = new Double_t[ncomp];
985  Double_t amed = 0.;
986  Int_t i, j;
987  for (i=0; i<ncomp; i++) {
988  elem = (TGeoElement *)pop->At(i);
989  if (!elem->IsRadioNuclide()) {
990  j = fElements->IndexOf(elem);
991  weight[i] = fWeights[j]*fAmixture[0]/fWeights[0];
992  } else {
993  el = (TGeoElementRN*)elem;
994  weight[i] = el->Ratio()->Concentration(time) * el->A();
995  }
996  amed += weight[i];
997  }
998  Double_t rho = fDensity * fWeights[0] * amed/fAmixture[0];
999  TGeoMixture *mix = 0;
1000  Int_t ncomp1 = ncomp;
1001  for (i=0; i<ncomp; i++) {
1002  if ((weight[i]/amed)<precision) {
1003  amed -= weight[i];
1004  ncomp1--;
1005  }
1006  }
1007  if (ncomp1<2) {
1008  el = (TGeoElementRN *)pop->At(0);
1009  delete [] weight;
1010  delete pop;
1011  if (ncomp1==1) return new TGeoMaterial(TString::Format("%s-evol",GetName()), el, rho);
1012  return NULL;
1013  }
1014  mix = new TGeoMixture(TString::Format("%s-evol",GetName()), ncomp, rho);
1015  for (i=0; i<ncomp; i++) {
1016  weight[i] /= amed;
1017  if (weight[i]<precision) continue;
1018  el = (TGeoElementRN *)pop->At(i);
1019  mix->AddElement(el, weight[i]);
1020  }
1021  delete [] weight;
1022  delete pop;
1023  return mix;
1024 }
1025 
1026 ////////////////////////////////////////////////////////////////////////////////
1027 /// Fills a user array with all the elements deriving from the possible
1028 /// decay of the top elements composing the mixture. Each element contained
1029 /// by <population> may be a radionuclide having a Bateman solution attached.
1030 /// The precision represent the minimum cumulative branching ratio for
1031 /// which decay products are still taken into account.
1032 /// To visualize the time evolution of each decay product one can use:
1033 /// TGeoElement *elem = population->At(index);
1034 /// TGeoElementRN *elemrn = 0;
1035 /// if (elem->IsRadioNuclide()) elemrn = (TGeoElementRN*)elem;
1036 /// One can get Ni/N1(t=0) at any moment of time. Ni is the number of atoms
1037 /// of one of the decay products, N1(0) is the number of atoms of the first top
1038 /// element at t=0.
1039 /// Double_t fraction_weight = elemrn->Ratio()->Concentration(time);
1040 /// One can also display the time evolution of the fractional weigth:
1041 /// elemrn->Ratio()->Draw(option);
1042 
1044 {
1045  if (population->GetEntriesFast()) {
1046  Error("FillMaterialEvolution", "Provide an empty array !");
1047  return;
1048  }
1050  TGeoElement *elem;
1051  TGeoElementRN *elemrn;
1052  TIter next(table->GetElementsRN());
1053  while ((elemrn=(TGeoElementRN*)next())) elemrn->ResetRatio();
1054  Double_t factor;
1055  for (Int_t i=0; i<fNelements; i++) {
1056  elem = GetElement(i);
1057  if (!elem->IsRadioNuclide()) {
1058  population->Add(elem);
1059  continue;
1060  }
1061  elemrn = (TGeoElementRN*)elem;
1062  factor = fWeights[i]*fAmixture[0]/(fWeights[0]*fAmixture[i]);
1063  elemrn->FillPopulation(population, precision, factor);
1064  }
1065 }
1066 
1067 ////////////////////////////////////////////////////////////////////////////////
1068 /// static function
1069 /// Compute screening factor for pair production and Bremstrahlung
1070 /// REFERENCE : EGS MANUAL SLAC 210 - UC32 - JUNE 78
1071 /// FORMULA 2.7.22
1072 
1074 {
1075  const Double_t al183= 5.20948 , al1440 = 7.27239;
1076  Double_t alz = TMath::Log(z)/3.;
1077  Double_t factor = (al1440 - 2*alz) / (al183 - alz - TGeoMaterial::Coulomb(z));
1078  return factor;
1079 }
virtual const char * GetTitle() const
Returns title of object.
Definition: TNamed.h:52
Int_t * fNatoms
Definition: TGeoMaterial.h:166
void FillPopulation(TObjArray *population, Double_t precision=0.001, Double_t factor=1.)
Fills the input array with the set of RN elements resulting from the decay of this one...
An array of TObjects.
Definition: TObjArray.h:39
TList * GetListOfMaterials() const
Definition: TGeoManager.h:463
void DefineElement(Int_t iel, Double_t a, Double_t z, Double_t weight)
Definition: TGeoMaterial.h:212
Double_t * fWeights
Definition: TGeoMaterial.h:165
Double_t Log(Double_t x)
Definition: TMath.h:526
ClassImp(TSeqCollection) Int_t TSeqCollection TIter next(this)
Return index of object in collection.
TGeoExtension * fFWExtension
Transient user-defined extension to materials.
Definition: TGeoMaterial.h:77
Double_t Concentration(Double_t time) const
Find concentration of the element at a given time.
TObject * fCerenkov
Definition: TGeoMaterial.h:74
const char Option_t
Definition: RtypesCore.h:62
Int_t Z() const
Definition: TGeoElement.h:76
Int_t fNelements
Definition: TGeoMaterial.h:162
Double_t * GetZmixt() const
Definition: TGeoMaterial.h:196
TGeoElementTable * GetElementTable()
Returns material table. Creates it if not existing.
void SetRadLen(Double_t radlen, Double_t intlen=0.)
Set radiation/absorbtion lengths.
Basic string class.
Definition: TString.h:137
int Int_t
Definition: RtypesCore.h:41
bool Bool_t
Definition: RtypesCore.h:59
TArc * a
Definition: textangle.C:12
const Bool_t kFALSE
Definition: Rtypes.h:92
virtual Int_t GetDefaultColor() const
Get some default color related to this material.
virtual ~TGeoMaterial()
Destructor.
Int_t GetEntriesFast() const
Definition: TObjArray.h:66
virtual Bool_t IsEqual(const TObject *obj) const
Default equal comparison (objects are equal if they have the same address in memory).
Definition: TObject.cxx:527
void AverageProperties()
Compute effective A/Z and radiation length.
static const Double_t STP_pressure
Definition: TGeoMaterial.h:46
virtual TGeoElement * GetElement(Int_t i=0) const
Get a pointer to the element this material is made of.
Short_t Abs(Short_t d)
Definition: TMathBase.h:110
LongDouble_t Power(LongDouble_t x, LongDouble_t y)
Definition: TMath.h:501
void SetBit(UInt_t f, Bool_t set)
Set or unset the user status bits as specified in f.
Definition: TObject.cxx:732
virtual void FillMaterialEvolution(TObjArray *population, Double_t precision=0.001)
Fills a user array with all the elements deriving from the possible decay of the top elements composi...
const char * Data() const
Definition: TString.h:349
virtual Double_t GetSpecificActivity(Int_t i=-1) const
Get specific activity (in Bq/gram) for the whole mixture (no argument) or for a given component...
virtual void Fatal(const char *method, const char *msgfmt,...) const
Issue fatal error message.
Definition: TObject.cxx:946
Double_t A() const
Definition: TGeoElement.h:79
EGeoMaterialState fState
Definition: TGeoMaterial.h:72
Fill Area Attributes class.
Definition: TAttFill.h:32
TGeoBatemanSol * Ratio() const
Definition: TGeoElement.h:195
static TString Format(const char *fmt,...)
Static method which formats a string using a printf style format descriptor and return a TString...
Definition: TString.cxx:2334
Double_t * GetWmixt() const
Definition: TGeoMaterial.h:198
The TNamed class is the base class for all named ROOT classes.
Definition: TNamed.h:33
void SetUserExtension(TGeoExtension *ext)
Connect user-defined extension to the material.
ClassImp(TGeoMaterial) TGeoMaterial
Default constructor.
virtual void SavePrimitive(std::ostream &out, Option_t *option="")
Save a primitive as a C++ statement(s) on output stream "out".
virtual Bool_t IsMixture() const
Definition: TGeoMaterial.h:127
Double_t Neff() const
Returns effective number of nucleons.
virtual TGeoElement * GetElement(Int_t i=0) const
Retreive the pointer to the element corresponding to component I.
TGeoElement * GetBaseElement() const
Definition: TGeoMaterial.h:110
virtual void Error(const char *method, const char *msgfmt,...) const
Issue error message.
Definition: TObject.cxx:918
virtual void FillMaterialEvolution(TObjArray *population, Double_t precision=0.001)
Fills a user array with all the elements deriving from the possible decay of the top element composin...
virtual Double_t GetSpecificActivity() const
Definition: TGeoElement.h:85
char * out
Definition: TBase64.cxx:29
static Double_t Coulomb(Double_t z)
static function Compute Coulomb correction for pair production and Brem REFERENCE : EGS MANUAL SLAC 2...
virtual Double_t GetA() const
Definition: TGeoMaterial.h:103
A doubly linked list.
Definition: TList.h:47
Int_t GetIndex()
Retreive material index in the list of materials.
void AddElement(Double_t a, Double_t z, Double_t weight)
add an element to the mixture using fraction by weight Check if the element is already defined ...
Double_t fIntLen
Definition: TGeoMaterial.h:69
virtual Int_t GetNelements() const
Definition: TGeoMaterial.h:195
virtual Bool_t IsEq(const TGeoMaterial *other) const
Return true if the other material has the same physical properties.
virtual void Release() const =0
TNamed & operator=(const TNamed &rhs)
TNamed assignment operator.
Definition: TNamed.cxx:40
Int_t IndexOf(const TObject *obj) const
Definition: TObjArray.cxx:551
virtual void AddAtAndExpand(TObject *obj, Int_t idx)
Add object at position idx.
Definition: TObjArray.cxx:221
virtual TGeoMaterial * DecayMaterial(Double_t time, Double_t precision=0.001)
Create the material representing the decay product of this material at a given time.
Double_t fDensity
Definition: TGeoMaterial.h:67
char * GetPointerName() const
Provide a pointer name containing uid.
TGeoMixture & operator=(const TGeoMixture &)
assignment operator
TGeoExtension * GrabUserExtension() const
Get a copy of the user extension pointer.
Double_t fTemperature
Definition: TGeoMaterial.h:70
Bool_t TestBit(UInt_t f) const
Definition: TObject.h:173
TGeoElement * fElement
Definition: TGeoMaterial.h:75
Double_t E()
Definition: TMath.h:54
void SetFWExtension(TGeoExtension *ext)
Connect framework defined extension to the material.
virtual const char * GetName() const
Returns name of object.
Definition: TNamed.h:51
TSubString Strip(EStripType s=kTrailing, char c= ' ') const
Return a substring of self stripped at beginning and/or end.
Definition: TString.cxx:1069
TString fName
Definition: TNamed.h:36
virtual void AddAt(TObject *obj, Int_t idx)
Add object at position ids.
Definition: TObjArray.cxx:238
static const Double_t STP_temperature
Definition: TGeoMaterial.h:45
virtual TObject * GetCerenkovProperties() const
Definition: TGeoMaterial.h:115
Double_t * fZmixture
Definition: TGeoMaterial.h:163
TObjArray * GetElementsRN() const
Definition: TGeoElement.h:411
Double_t fZ
Definition: TGeoMaterial.h:66
static Double_t ScreenFactor(Double_t z)
static function Compute screening factor for pair production and Bremstrahlung REFERENCE : EGS MANUAL...
void SetUsed(Bool_t flag=kTRUE)
Definition: TGeoElement.h:92
virtual ~TGeoMixture()
Destructor.
R__EXTERN TGeoManager * gGeoManager
Definition: TGeoManager.h:556
double Double_t
Definition: RtypesCore.h:55
virtual Bool_t IsRadioNuclide() const
Definition: TGeoElement.h:88
ClassImp(TMCParticle) void TMCParticle printf(": p=(%7.3f,%7.3f,%9.3f) ;", fPx, fPy, fPz)
virtual void Print(const Option_t *option="") const
print characteristics of this material
virtual void GetElementProp(Double_t &a, Double_t &z, Double_t &w, Int_t i=0)
Single interface to get element properties.
virtual void SavePrimitive(std::ostream &out, Option_t *option="")
Save a primitive as a C++ statement(s) on output stream "out".
void SetDefined(Bool_t flag=kTRUE)
Definition: TGeoElement.h:91
Double_t * GetAmixt() const
Definition: TGeoMaterial.h:197
static Double_t Big()
Definition: TGeoShape.h:98
virtual UInt_t GetUniqueID() const
Return the unique object id.
Definition: TObject.cxx:433
Double_t Na()
Definition: TMath.h:104
#define name(a, b)
Definition: linkTestLib0.cpp:5
Binding & operator=(OUT(*fun)(void))
virtual Double_t GetZ() const
Definition: TGeoMaterial.h:104
Double_t fPressure
Definition: TGeoMaterial.h:71
TObjArray * fElements
Definition: TGeoMaterial.h:167
Double_t fRadLen
Definition: TGeoMaterial.h:68
virtual TGeoExtension * Grab()=0
virtual Double_t GetDensity() const
Definition: TGeoMaterial.h:106
void SetUsed(Bool_t flag=kTRUE)
Definition: TGeoMaterial.h:136
TGeoExtension * fUserExtension
Definition: TGeoMaterial.h:76
Double_t fA
Definition: TGeoMaterial.h:65
#define NULL
Definition: Rtypes.h:82
Int_t AddMaterial(const TGeoMaterial *material)
Add a material to the list. Returns index of the material in list.
std::vector< Double_t > fWeights
Definition: TKDE.cxx:52
TGeoMaterial & operator=(const TGeoMaterial &)
assignment operator
void Add(TObject *obj)
Definition: TObjArray.h:75
Int_t GetNelements() const
Definition: TGeoElement.h:415
virtual TGeoMaterial * DecayMaterial(Double_t time, Double_t precision=0.001)
Create the mixture representing the decay product of this material at a given time.
TGeoElement * GetElement(Int_t z)
Definition: TGeoElement.h:408
virtual Int_t IndexOf(const TObject *obj) const
TObject * At(Int_t idx) const
Definition: TObjArray.h:167
Double_t * fAmixture
Definition: TGeoMaterial.h:164
TObject * fShader
Definition: TGeoMaterial.h:73
const Bool_t kTRUE
Definition: Rtypes.h:91
virtual void Print(const Option_t *option="") const
print characteristics of this material
virtual Bool_t IsEq(const TGeoMaterial *other) const
return true if the other material has the same physical properties
TGeoExtension * GrabFWExtension() const
Get a copy of the framework extension pointer.
void ResetRatio()
Clears the existing ratio.
virtual void Warning(const char *method, const char *msgfmt,...) const
Issue warning message.
Definition: TObject.cxx:904