GammaWare  Head Version for release 0.9
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
LStoSToGS.C
Go to the documentation of this file.
1 
19 #ifndef __CINT__
20 #include <Link.h>
21 #include <XGammaLink.h>
22 #include <GEM.h>
23 #include <GammaFilter.h>
24 
25 #include <TStopwatch.h>
26 #include <TH1.h>
27 #include <TTree.h>
28 #include <TFile.h>
29 #include <TDatabasePDG.h>
30 #include <TVector3.h>
31 
32 #include <SToGS_BaseROOTEvents.h>
33 
34 #include <iostream>
35 using namespace std;
36 #endif
37 
38 #include <Random.h>
39 
40 using namespace Gw;
41 
42 
52 void LStoSToGS1(const char *pathtoLS, Int_t nbcas = 100, const char *rfilename = "PrimariesFromGW.root")
53 {
54  // Random engine is TRandom3 based on clock ... just remove clock to have identical random sequences
55  Random::Instance()->SetCurrent("large_clock");
56 
57  TStopwatch watch;
58  watch.Start();
59 
60  // Inittialisation: experimental conditions
61  const Int_t nb_residu = 6; // number of residus
62  //
63  TString filename[nb_residu], nucleusname[nb_residu]; // database related: to get Level Scheme
64  //
65  Float_t sigma[nb_residu] = { 29.7, 47.1, 19.3, 19.3, 0.4, 1.31 }; // cross section to select randomly one residu
66  Int_t neutrons[nb_residu] = { 5, 4, 3, 4, 3, 4 }; // number of neutron emitted
67  Int_t protons[nb_residu] = { 0, 0, 0, 1, 0, 0 }; // number of proton emitted
68  Int_t alphas [nb_residu] = { 0, 0, 0, 0, 1, 1 }; // number of alpha emitted
69 
70  BaseGEM nucleus[nb_residu]; // gem related
71 
72  filename[0] = pathtoLS; filename[0] += "DY_ENSDF.ens"; nucleusname[0] = "151DY";
73  filename[1] = pathtoLS; filename[1] += "152DY.ags"; nucleusname[1] = "152DY";
74  filename[2] = pathtoLS; filename[2] += "DY_ENSDF.ens"; nucleusname[2] = "153DY";
75  filename[3] = pathtoLS; filename[3] += "TB_ENSDF.ens"; nucleusname[3] = "151TB";
76  filename[4] = pathtoLS; filename[4] += "GD_ENSDF.ens"; nucleusname[4] = "149GD";
77  filename[5] = pathtoLS; filename[5] += "148GD.ags"; nucleusname[5] = "148GD";
78 
79  // same prototype function to init all nucleus
80  Int_t ok = 0;
81  for (int i = 0; i < nb_residu; i++ ) {
82  ok += nucleus[i].Import(filename[i],nucleusname[i]);
83  }
84 
85  if ( ok != nb_residu ) {
86  cout << " A Level Scheme is missing !! " << endl;
87  return;
88  }
89  // create the RandObj to select amongst the nb_residu nuclei.
90  RandObj channel;
91  for (int i = 0; i < nb_residu; i++ )
92  channel.Add(nucleus+i,sigma[i]);
93 
94  // ROOT file
95  TFile rf(rfilename,"UPDATE"); TTree *tree_primary = new TTree("GWLS_to_SToGS","Primaries from GW");
96  //
97  SBRPEvent pr_event;
98  //
99  tree_primary->Branch("Pr.",&pr_event);
100 
101  // some spectra to be filled
102  TH1F *hnbgamma = new TH1F("NbGamma","NbGamma",100,0,100); TH1F *hegamma = new TH1F("GammaSpectrum","Gamma Spectrum",8192,0,4096);
103  Int_t gPDG = TDatabasePDG::Instance()->GetParticle("gamma")->PdgCode();
104  TH1F *hnbprotons = new TH1F("NbProtons","NbProtons",100,0,100); TH1F *heprotons = new TH1F("ProtonSpectrum","Proton Spectrum",8192,0,4096);
105  Int_t pPDG = TDatabasePDG::Instance()->GetParticle("proton")->PdgCode();
106  TH1F *hnbneutrons = new TH1F("NbNeutrons","NbNeutrons",100,0,100); TH1F *heneutrons = new TH1F("NeutronSpectrum","Neutron Spectrum",8192,0,4096);
107  Int_t nPDG = TDatabasePDG::Instance()->GetParticle("neutron")->PdgCode();
108  // PDG: Nuclear codes are given as 10-digit numbers ±10LZZZAAAI
109  // PDG: L = the total number of strange quarks - 0 if not hypernuclei
110  // PDG:I gives the isomer level, with I = 0 corresponding to the ground state and I > 0 to excitations
111  TH1F *hnbalphas = new TH1F("NbAlphas","NbAlphas",100,0,100); TH1F *healphas = new TH1F("AlphaSpectrum","Alpha Spectrum",8192,0,4096);
112  Int_t aPDG = 1000020040;
113 
114  watch.Stop();
115  cout << "Initialisation time " << endl; watch.Print();
116 
117  // start the Monte-Carlo
118  watch.Reset();
119  watch.Start();
120 
121  Cascade cas; BaseGEM *gem; GammaLink *gam; Int_t which; Double_t ESUM = 0; Int_t MULTTOT = 0; TVector3 angles;
122  for (Int_t i = 0; i < nbcas; i++ ) {
123 
124  ESUM = 0.0;
125  MULTTOT = 0;
126  //
127  gem = (BaseGEM *)channel.Rand(which); pr_event.Clear(""); // select one nucleus
128 
129  // GAMMA-rays from LS
130  gem->DoCascade(cas);
131  //
132  MULTTOT += cas.GetSize();
133  for (Int_t j = 0; j < MULTTOT; j++ ) {
134  gam = (GammaLink *)cas.At(j);
135 
136  SBRPHit *stogshit = pr_event.AddHit();
137 
138  // now copy gam content into stogshit
139  stogshit->fE = gam->GetEnergy().Get();
140  stogshit->fX = 0.; // X position
141  stogshit->fY = 0.; // Y position
142  stogshit->fZ = 0.; // Z position
143 
144  Random::Instance()->Current()->Sphere(stogshit->fPX,stogshit->fPY,stogshit->fPZ,1.0);
145 
146  stogshit->fT = 0.; //
147 
148  stogshit->fPDG = gPDG;
149  stogshit->fFlag = j;
150  stogshit->fUID = 0;
151 
152  hegamma->Fill(stogshit->fE); ESUM += stogshit->fE;
153  } // j
154  hnbgamma->Fill(MULTTOT) ;
155 
156  // Neutrons
157  MULTTOT += neutrons[which];
158  for (Int_t j = 0; j < neutrons[which]; j++) {
159  SBRPHit *stogshit = pr_event.AddHit();
160 
161  // now copy gam content into stogshit
162  stogshit->fE = 1000.0;
163  stogshit->fX = 0.; // X position
164  stogshit->fY = 0.; // Y position
165  stogshit->fZ = 0.; // Z position
166 
167  Random::Instance()->Current()->Sphere(stogshit->fPX,stogshit->fPY,stogshit->fPZ,1.0);
168 
169  stogshit->fT = 0.; //
170 
171  stogshit->fPDG = nPDG;
172  stogshit->fFlag = j;
173  stogshit->fUID = 0;
174 
175  heneutrons->Fill(stogshit->fE);; ESUM += stogshit->fE;
176  }
177  hnbneutrons->Fill(neutrons[which]) ;
178 
179  // Protons
180  MULTTOT += protons[which];
181  for (Int_t j = 0; j < protons[which]; j++) {
182  SBRPHit *stogshit = pr_event.AddHit();
183 
184  // now copy gam content into stogshit
185  stogshit->fE = 2000.0;
186  stogshit->fX = 0.; // X position
187  stogshit->fY = 0.; // Y position
188  stogshit->fZ = 0.; // Z position
189 
190  Random::Instance()->Current()->Sphere(stogshit->fPX,stogshit->fPY,stogshit->fPZ,1.0);
191 
192  stogshit->fT = 0.; //
193 
194  stogshit->fPDG = pPDG;
195  stogshit->fFlag = j;
196  stogshit->fUID = 0;
197 
198  heprotons->Fill(stogshit->fE);; ESUM += stogshit->fE;
199  }
200  hnbprotons->Fill(protons[which]) ;
201 
202  // alphas
203  MULTTOT += alphas[which];
204  for (Int_t j = 0; j < alphas[which]; j++) {
205  SBRPHit *stogshit = pr_event.AddHit();
206 
207  // now copy gam content into stogshit
208  stogshit->fE = 3000;
209  stogshit->fX = 0.; // X position
210  stogshit->fY = 0.; // Y position
211  stogshit->fZ = 0.; // Z position
212 
213  Random::Instance()->Current()->Sphere(stogshit->fPX,stogshit->fPY,stogshit->fPZ,1.0);
214 
215  stogshit->fT = 0.; //
216 
217  stogshit->fPDG = aPDG;
218  stogshit->fFlag = j;
219  stogshit->fUID = 0;
220 
221  healphas->Fill(stogshit->fE);; ESUM += stogshit->fE;
222  }
223  hnbalphas->Fill(alphas[which]) ;
224 
225  pr_event.SetEMult(ESUM,MULTTOT);
226  tree_primary->Fill();
227 
228  } //i
229 
230  rf.Write();
231  watch.Stop();
232  cout << "Monte-Carlo time " << endl; watch.Print();
233 
234 }
235 
250 void testGEM(const char *pathtoLS, Int_t nbcas = 100)
251 {
252  TStopwatch watch;
253 
254  watch.Start();
255 
256  // Experimental conditions (from PACE like)
257  const Int_t nb_residu = 6; // 151DY, 152DY, 153DY, 151TB, 149GD, 148GD
258 
259  GEM nucleus[nb_residu]; // gem related
260 
261  TString filename[nb_residu]; // database related
262  TString nucleusname[nb_residu]; // database related
263 
264  Float_t sigma[nb_residu] = { 29.7, 47.1, 19.3, 19.3, 0.4, 1.31 };
265 
266  filename[0] = pathtoLS; filename[0] += "DY_ENSDF.ens"; nucleusname[0] = "151DY";
267  filename[1] = pathtoLS; filename[1] += "152DY.ags"; nucleusname[1] = "152DY";
268  filename[2] = pathtoLS; filename[2] += "DY_ENSDF.ens"; nucleusname[2] = "153DY";
269  filename[3] = pathtoLS; filename[3] += "TB_ENSDF.ens"; nucleusname[3] = "151TB";
270  filename[4] = pathtoLS; filename[4] += "GD_ENSDF.ens"; nucleusname[4] = "149GD";
271  filename[5] = pathtoLS; filename[5] += "148GD.ags"; nucleusname[5] = "148GD";
272 
273  // same prototype function to init all nucleus
274  Int_t ok = 0;
275  for (int i = 0; i < nb_residu; i++ ) {
276  ok += nucleus[i].Import(filename[i],nucleusname[i]);
277  // entry point 80 MeV, spin 60
278  nucleus[i].SetFeeding("TestFeeding"); nucleus[i].SetExI(80000,60);
279  }
280 
281  if ( ok != nb_residu ) {
282  cout << " A Level Scheme is missing !! " << endl; return;
283  }
284  // create the RandObj to select the nb_residu nuclei.
285  RandObj channel;
286  for (int i = 0; i < nb_residu; i++ ) channel.Add(nucleus+i,sigma[i]);
287 
288  // some spectra to be filled
289  TH1F *K = new TH1F("GEM_K","K",100,0,100);
290  TH1F *projtot = new TH1F("GEM","GEM",8192,0,4096);
291 
292  watch.Stop();
293  cout << "Initialisation time " << endl; watch.Print();
294  // start the Monte-Carlo
295 
296  watch.Reset(); watch.Start();
297  Cascade cas; BaseGEM *gem; GammaLink *gam; Int_t which;
298  for (Int_t i = 0; i < nbcas; i++ ) {
299 
300  gem = (BaseGEM *)channel.Rand(which); // select one nucleus
301  gem->DoCascade(cas); // determine a cascade on this nucleus
302 
303  K->Fill(cas.GetSize()) ;
304  for (Int_t j = 0; j < cas.GetSize(); j++ ) {
305  gam = (GammaLink *)cas.At(j); projtot->Fill(gam->GetEnergy().Get());
306  } // j
307  } //i
308  watch.Stop();
309  cout << "Monte-Carlo time " << endl; watch.Print();
310 
311  projtot->Draw();
312 }
313 
332 void testGammaFilter1(const char *pathtoLS, Int_t nbcas = 100, const char *filtername = "AGATA_FILTER.root")
333 {
334  TStopwatch watch;
335 
336  watch.Start();
337 
338  // Experimental conditions (from PACE like)
339  const Int_t nb_residu = 6; // 151DY, 152DY, 153DY, 151TB, 149GD, 148GD
340 
341  BaseGEM nucleus[nb_residu]; // gem related
342 
343  TString filename[nb_residu]; // database related
344  TString nucleusname[nb_residu]; // database related
345 
346  Float_t sigma[nb_residu] = { 29.7, 47.1, 19.3, 19.3, 0.4, 1.31 };
347 
348  filename[0] = pathtoLS; filename[0] += "DY_ENSDF.ens"; nucleusname[0] = "151DY";
349  filename[1] = pathtoLS; filename[1] += "152DY.ags"; nucleusname[1] = "152DY";
350  filename[2] = pathtoLS; filename[2] += "DY_ENSDF.ens"; nucleusname[2] = "153DY";
351  filename[3] = pathtoLS; filename[3] += "TB_ENSDF.ens"; nucleusname[3] = "151TB";
352  filename[4] = pathtoLS; filename[4] += "GD_ENSDF.ens"; nucleusname[4] = "149GD";
353  filename[5] = pathtoLS; filename[5] += "148GD.ags"; nucleusname[5] = "148GD";
354 
355  // same prototype function to init all nucleus
356  Int_t ok = 0;
357  for (int i = 0; i < nb_residu; i++ )
358  { ok += nucleus[i].Import(filename[i],nucleusname[i]); }
359 
360  if ( ok != nb_residu ) {
361  cout << " A Level Scheme is missing !! " << endl; return;
362  }
363  // create the RandObj to select the nb_residu nuclei.
364  RandObj channel;
365  for (int i = 0; i < nb_residu; i++ ) channel.Add(nucleus+i,sigma[i]);
366 
367  // initiate the GammaFilter
368  Float_t e[1000]; // up to 1000 energies .. should be ok for the moment
369  GammaFilter filter; filter.InitFilter(filtername);
370  if ( filter.IsEffective() == false ) return ;
371 
372  // some spectra to be filled
373  TH1F *K = new TH1F("GammaFilter1_K","K",100,0,100);
374  TH1F *projtot = new TH1F("GammaFilter1","GammaFilter1",8192,0,4096);
375 
376  watch.Stop();
377  cout << "Initialisation time " << endl; watch.Print();
378  // start the Monte-Carlo
379 
380  watch.Reset(); watch.Start();
381  Cascade cas; BaseGEM *gem; GammaLink *gam; Int_t which;
382  for (Int_t i = 0; i < nbcas; i++ ) {
383 
384  gem = (BaseGEM *)channel.Rand(which); // select one nucleus
385  gem->DoCascade(cas); // determine a cascade on this nucleus
386 
387  // apply experimental filter on this cascade.
388  Int_t nbgamma = filter.ApplyE(cas,e);
389 
390  K->Fill(nbgamma) ;
391  for (Int_t j = 0; j < nbgamma; j++ ) {
392  projtot->Fill(e[j]);
393  } // j
394  } //i
395  watch.Stop();
396  cout << "Monte-Carlo time " << endl; watch.Print();
397 
398  projtot->Draw();
399 }
400 
420 void testGammaFilter2(const char *pathtoLS, Int_t nbcas = 100, const char *filtername = "AGATA_FILTER.root")
421 {
422  TStopwatch watch;
423 
424  watch.Start();
425 
426  // Experimental conditions (from PACE like)
427  const Int_t nb_residu = 6; // 151DY, 152DY, 153DY, 151TB, 149GD, 148GD
428 
429  GEM nucleus[nb_residu]; // gem related
430 
431  TString filename[nb_residu]; // database related
432  TString nucleusname[nb_residu]; // database related
433 
434  Float_t sigma[nb_residu] = { 29.7, 47.1, 19.3, 19.3, 0.4, 1.31 };
435 
436  filename[0] = pathtoLS; filename[0] += "DY_ENSDF.ens"; nucleusname[0] = "151DY";
437  filename[1] = pathtoLS; filename[1] += "152DY.ags"; nucleusname[1] = "152DY";
438  filename[2] = pathtoLS; filename[2] += "DY_ENSDF.ens"; nucleusname[2] = "153DY";
439  filename[3] = pathtoLS; filename[3] += "TB_ENSDF.ens"; nucleusname[3] = "151TB";
440  filename[4] = pathtoLS; filename[4] += "GD_ENSDF.ens"; nucleusname[4] = "149GD";
441  filename[5] = pathtoLS; filename[5] += "148GD.ags"; nucleusname[5] = "148GD";
442 
443  // same prototype function to init all nucleus
444  Int_t ok = 0;
445  for (int i = 0; i < nb_residu; i++ ) {
446  ok += nucleus[i].Import(filename[i],nucleusname[i]);
447  // entry point 80 MeV, spin 60
448  nucleus[i].SetFeeding("TestFeeding"); nucleus[i].SetExI(80000,60);
449  }
450 
451  if ( ok != nb_residu ) {
452  cout << " A Level Scheme is missing !! " << endl; return;
453  }
454  // create the RandObj to select the nb_residu nuclei.
455  RandObj channel;
456  for (int i = 0; i < nb_residu; i++ ) channel.Add(nucleus+i,sigma[i]);
457 
458  // initiate the GammaFilter
459  Int_t nbgamma; Float_t e[1000]; // up to 1000 energies ..should be ok for the moment
460  GammaFilter filter; filter.InitFilter(filtername);
461  if ( filter.IsEffective() == false ) return ;
462 
463  // some spectra to be filled
464  TH1F *K = new TH1F("GammaFilter2_K","K",100,0,100);
465  TH1F *projtot = new TH1F("GammaFilter2","GammaFilter2",8192,0,4096);
466 
467  // open a root file to store the produced events in a ROOT TTree
468  TFile fileout("SimuDY.root","recreate");
469 
470  TTree *treeout;
471  treeout = (TTree *)fileout.Get("SimuDY");
472  if ( treeout == NULL ) {
473  treeout = new TTree("SimuDY","FromGammaFilter2");
474  treeout->Branch("mult",&nbgamma,"nbgamma/I"); treeout->Branch("e",e,"e[nbgamma]/F");
475  }
476  // else {
477  // treeout->SetBranchAddress("nbgamma",&nbgamma);
478  //
479  // }
480 
481  watch.Stop();
482  cout << "Initialisation time " << endl; watch.Print();
483  // start the Monte-Carlo
484 
485  watch.Reset(); watch.Start();
486  Cascade cas; BaseGEM *gem; GammaLink *gam; Int_t which;
487  for (Int_t i = 0; i < nbcas; i++ ) {
488 
489  gem = (BaseGEM *)channel.Rand(which); // select one nucleus
490  gem->DoCascade(cas); // determine a cascade on this nucleus
491 
492  // apply experimental filter on this cascade.
493  nbgamma = filter.ApplyE(cas,e);
494 
495  K->Fill(nbgamma) ;
496  for (Int_t j = 0; j < nbgamma; j++ ) {
497  projtot->Fill(e[j]);
498  } // j
499 
500  treeout->Fill();
501  } //i
502 
503  treeout->Write();
504 
505  watch.Stop();
506  cout << "Monte-Carlo time " << endl; watch.Print();
507 
508  projtot->Draw();
509 }
510 
511 
Gw::RandObj::Rand
virtual TObject * Rand()
it returns a pointer to a randomly selected object from the collection
Definition: RandObj.cpp:93
LStoSToGS1
void LStoSToGS1(const char *pathtoLS, Int_t nbcas=100, const char *rfilename="PrimariesFromGW.root")
Function to generate discret gamma-rays between several level schemes.
Definition: LStoSToGS.C:52
Gw::GEM::SetExI
void SetExI(Float_t E, Float_t I)
starting point in the ExI plan from which starts the gamma-ray emission
Definition: GEM.cpp:59
Gw::Cascade
A cascade is a list of links.
Definition: Cascade.h:51
Gw::GammaFilter::ApplyE
virtual Int_t ApplyE(const TSeqCollection &, TBuffer &)
Apply this filter to the list of gamma-rays.
Definition: GammaFilter.cpp:117
SBRPHit
SToGS Base Root Polarized Hit.
Definition: SToGS_BaseROOTEvents.h:116
SBRHit::fZ
Double32_t fZ
Definition: SToGS_BaseROOTEvents.h:86
Gw::GEM
Class to get randomly cascades (with lateral feeding) of gammas on the basis of a level scheme...
Definition: GEM.h:55
SBRHit::fE
Double32_t fE
Definition: SToGS_BaseROOTEvents.h:80
SBRPHit::fPZ
Double32_t fPZ
Definition: SToGS_BaseROOTEvents.h:122
testGammaFilter1
void testGammaFilter1(const char *pathtoLS, Int_t nbcas=100, const char *filtername="AGATA_FILTER.root")
Function to generate tracked like events.
Definition: LStoSToGS.C:332
Gw::GammaFilter::InitFilter
Bool_t InitFilter(const char *)
to modifiy the efficency curve.
Definition: GammaFilter.cpp:64
SBRPEvent::AddHit
SBRPHit * AddHit()
add a hit to the current event
Definition: SToGS_BaseROOTEvents.C:204
Gw::BaseGEM
Class to get randomly cascades of gammas on the basis of a level scheme.
Definition: BaseGEM.h:107
SBRHit::fX
Double32_t fX
Definition: SToGS_BaseROOTEvents.h:84
SBRPHit::fPX
Double32_t fPX
Definition: SToGS_BaseROOTEvents.h:120
SToGS_BaseROOTEvents.h
Gw::RandObj
A class to select randomly an object in a TObjArray of objects.
Definition: RandObj.h:55
SBRHit::fT
Double32_t fT
Definition: SToGS_BaseROOTEvents.h:87
Gw::GEM::SetFeeding
Link * SetFeeding(const char *feedingtype="TestFeeding")
Definition: GEM.cpp:47
Gw::Data::Get
virtual Data_T Get() const
get the value, can be overloaded
Definition: Data.h:70
Gw::RandObj::Add
virtual void Add(TObject *, Float_t)
it adds an object with a weight to this collection.
Definition: RandObj.cpp:49
GEM.h
header file for a GEM (gamma-rays generator)
SBRPHit::fPY
Double32_t fPY
Definition: SToGS_BaseROOTEvents.h:121
Gw::GammaFilter
Standard filter for gamma-rays (EUROBALL-like appproach)
Definition: GammaFilter.h:41
Random.h
header file for Random
testGEM
void testGEM(const char *pathtoLS, Int_t nbcas=100)
Function to generate discret gamma-rays between several level schemes ...
Definition: LStoSToGS.C:250
SBRPEvent::Clear
virtual void Clear(Option_t *opt)
clear the collection of hits, set H, K to 0
Definition: SToGS_BaseROOTEvents.C:196
Gw::BaseGEM::DoCascade
virtual Int_t DoCascade(TSeqCollection &cascade, Option_t *opt="")
cascade simulation
Definition: BaseGEM.cpp:523
endl
ADF::LogMessage & endl(ADF::LogMessage &log)
Definition: ADFLogMessage.cpp:265
Gw::GammaFilter::IsEffective
Bool_t IsEffective()
check out if it has been correctly initiated
Definition: GammaFilter.cpp:95
testGammaFilter2
void testGammaFilter2(const char *pathtoLS, Int_t nbcas=100, const char *filtername="AGATA_FILTER.root")
Function to generate tracked like events and store them in a ROOT TTree.
Definition: LStoSToGS.C:420
Gw::BaseGEM::Import
virtual Int_t Import(const Char_t *, Option_t *opt="152DY")
Load the level scheme and init the Monte-Carlo.
Definition: BaseGEM.cpp:194
SBRPEvent::SetEMult
void SetEMult(Double_t h, Int_t k)
Definition: SToGS_BaseROOTEvents.h:274
SBRHit::fY
Double32_t fY
Definition: SToGS_BaseROOTEvents.h:85