GammaWare  Head Version for release 0.9
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
FirstGEM.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 
29 #include <iostream>
30 using namespace std;
31 #endif
32 
33 using namespace Gw;
34 
35 #define MAX_STRING_LENGTH 50
36 
49 void testBaseGEM(const char *pathtodata, Int_t nbcas = 100)
50 {
51  TStopwatch watch;
52 
53  watch.Start();
54 
55  // Experimental conditions (from PACE like)
56  const Int_t nb_residu = 6; // 151DY, 152DY, 153DY, 151TB, 149GD, 148GD
57 
58  BaseGEM nucleus[nb_residu]; // gem related
59 
60  TString filename[nb_residu]; // database related
61  TString nucleusname[nb_residu]; // database related
62 
63  Float_t sigma[nb_residu] = { 29.7, 47.1, 19.3, 19.3, 0.4, 1.31 };
64 
65  filename[0] = pathtodata; filename[0] += "DY_ENSDF.ens"; nucleusname[0] = "151DY";
66  filename[1] = pathtodata; filename[1] += "152DY.ags"; nucleusname[1] = "152DY";
67  filename[2] = pathtodata; filename[2] += "DY_ENSDF.ens"; nucleusname[2] = "153DY";
68  filename[3] = pathtodata; filename[3] += "TB_ENSDF.ens"; nucleusname[3] = "151TB";
69  filename[4] = pathtodata; filename[4] += "GD_ENSDF.ens"; nucleusname[4] = "149GD";
70  filename[5] = pathtodata; filename[5] += "148GD.ags"; nucleusname[5] = "148GD";
71 
72  // same prototype function to init all nucleus
73  Int_t ok = 0;
74  for (int i = 0; i < nb_residu; i++ )
75  { ok += nucleus[i].Import(filename[i],nucleusname[i]); }
76 
77  if ( ok != nb_residu ) {
78  cout << " A Level Scheme is missing !! " << endl; return;
79  }
80  // create the RandObj to select the nb_residu nuclei.
81  RandObj channel;
82  for (int i = 0; i < nb_residu; i++ ) channel.Add(nucleus+i,sigma[i]);
83 
84  // some spectra to be filled
85  TH1F *K = new TH1F("BaseGEM_K","K",100,0,100);
86  TH1F *projtot = new TH1F("BaseGEM","BaseGEM",8192,0,4096);
87 
88  watch.Stop();
89  cout << "Initialisation time " << endl; watch.Print();
90  // start the Monte-Carlo
91 
92  watch.Reset(); watch.Start();
93  Cascade cas; BaseGEM *gem; GammaLink *gam; Int_t which;
94  for (Int_t i = 0; i < nbcas; i++ ) {
95 
96  gem = (BaseGEM *)channel.Rand(which); // select one nucleus
97  gem->DoCascade(cas); // determine a cascade on this nucleus
98 
99  K->Fill(cas.GetSize()) ;
100  for (Int_t j = 0; j < cas.GetSize(); j++ ) {
101  gam = (GammaLink *)cas.At(j); projtot->Fill(gam->GetEnergy().Get());
102  } // j
103  } //i
104  watch.Stop();
105  cout << "Monte-Carlo time " << endl; watch.Print();
106 
107  projtot->Draw();
108 }
109 
124 void testGEM(const char *pathtodata, Int_t nbcas = 100)
125 {
126  TStopwatch watch;
127 
128  watch.Start();
129 
130  // Experimental conditions (from PACE like)
131  const Int_t nb_residu = 6; // 151DY, 152DY, 153DY, 151TB, 149GD, 148GD
132 
133  GEM nucleus[nb_residu]; // gem related
134 
135  TString filename[nb_residu]; // database related
136  TString nucleusname[nb_residu]; // database related
137 
138  Float_t sigma[nb_residu] = { 29.7, 47.1, 19.3, 19.3, 0.4, 1.31 };
139 
140  filename[0] = pathtodata; filename[0] += "DY_ENSDF.ens"; nucleusname[0] = "151DY";
141  filename[1] = pathtodata; filename[1] += "152DY.ags"; nucleusname[1] = "152DY";
142  filename[2] = pathtodata; filename[2] += "DY_ENSDF.ens"; nucleusname[2] = "153DY";
143  filename[3] = pathtodata; filename[3] += "TB_ENSDF.ens"; nucleusname[3] = "151TB";
144  filename[4] = pathtodata; filename[4] += "GD_ENSDF.ens"; nucleusname[4] = "149GD";
145  filename[5] = pathtodata; filename[5] += "148GD.ags"; nucleusname[5] = "148GD";
146 
147  // same prototype function to init all nucleus
148  Int_t ok = 0;
149  for (int i = 0; i < nb_residu; i++ ) {
150  ok += nucleus[i].Import(filename[i],nucleusname[i]);
151  // entry point 80 MeV, spin 60
152  nucleus[i].SetFeeding("TestFeeding"); nucleus[i].SetExI(80000,60);
153  }
154 
155  if ( ok != nb_residu ) {
156  cout << " A Level Scheme is missing !! " << endl; return;
157  }
158  // create the RandObj to select the nb_residu nuclei.
159  RandObj channel;
160  for (int i = 0; i < nb_residu; i++ ) channel.Add(nucleus+i,sigma[i]);
161 
162  // some spectra to be filled
163  TH1F *K = new TH1F("GEM_K","K",100,0,100);
164  TH1F *projtot = new TH1F("GEM","GEM",8192,0,4096);
165 
166  watch.Stop();
167  cout << "Initialisation time " << endl; watch.Print();
168  // start the Monte-Carlo
169 
170  watch.Reset(); watch.Start();
171  Cascade cas; BaseGEM *gem; GammaLink *gam; Int_t which;
172  for (Int_t i = 0; i < nbcas; i++ ) {
173 
174  gem = (BaseGEM *)channel.Rand(which); // select one nucleus
175  gem->DoCascade(cas); // determine a cascade on this nucleus
176 
177  K->Fill(cas.GetSize()) ;
178  for (Int_t j = 0; j < cas.GetSize(); j++ ) {
179  gam = (GammaLink *)cas.At(j); projtot->Fill(gam->GetEnergy().Get());
180  } // j
181  } //i
182  watch.Stop();
183  cout << "Monte-Carlo time " << endl; watch.Print();
184 
185  projtot->Draw();
186 }
187 
206 void testGammaFilter1(const char *pathtodata, Int_t nbcas = 100, const char *filtername = "AGATA_FILTER.root")
207 {
208  TStopwatch watch;
209 
210  watch.Start();
211 
212  // Experimental conditions (from PACE like)
213  const Int_t nb_residu = 6; // 151DY, 152DY, 153DY, 151TB, 149GD, 148GD
214 
215  BaseGEM nucleus[nb_residu]; // gem related
216 
217  TString filename[nb_residu]; // database related
218  TString nucleusname[nb_residu]; // database related
219 
220  Float_t sigma[nb_residu] = { 29.7, 47.1, 19.3, 19.3, 0.4, 1.31 };
221 
222  filename[0] = pathtodata; filename[0] += "DY_ENSDF.ens"; nucleusname[0] = "151DY";
223  filename[1] = pathtodata; filename[1] += "152DY.ags"; nucleusname[1] = "152DY";
224  filename[2] = pathtodata; filename[2] += "DY_ENSDF.ens"; nucleusname[2] = "153DY";
225  filename[3] = pathtodata; filename[3] += "TB_ENSDF.ens"; nucleusname[3] = "151TB";
226  filename[4] = pathtodata; filename[4] += "GD_ENSDF.ens"; nucleusname[4] = "149GD";
227  filename[5] = pathtodata; filename[5] += "148GD.ags"; nucleusname[5] = "148GD";
228 
229  // same prototype function to init all nucleus
230  Int_t ok = 0;
231  for (int i = 0; i < nb_residu; i++ )
232  { ok += nucleus[i].Import(filename[i],nucleusname[i]); }
233 
234  if ( ok != nb_residu ) {
235  cout << " A Level Scheme is missing !! " << endl; return;
236  }
237  // create the RandObj to select the nb_residu nuclei.
238  RandObj channel;
239  for (int i = 0; i < nb_residu; i++ ) channel.Add(nucleus+i,sigma[i]);
240 
241  // initiate the GammaFilter
242  Float_t e[1000]; // up to 1000 energies .. should be ok for the moment
243  GammaFilter filter; filter.InitFilter(filtername);
244  if ( filter.IsEffective() == false ) return ;
245 
246  // some spectra to be filled
247  TH1F *K = new TH1F("GammaFilter1_K","K",100,0,100);
248  TH1F *projtot = new TH1F("GammaFilter1","GammaFilter1",8192,0,4096);
249 
250  watch.Stop();
251  cout << "Initialisation time " << endl; watch.Print();
252  // start the Monte-Carlo
253 
254  watch.Reset(); watch.Start();
255  Cascade cas; BaseGEM *gem; GammaLink *gam; Int_t which;
256  for (Int_t i = 0; i < nbcas; i++ ) {
257 
258  gem = (BaseGEM *)channel.Rand(which); // select one nucleus
259  gem->DoCascade(cas); // determine a cascade on this nucleus
260 
261  // apply experimental filter on this cascade.
262  Int_t nbgamma = filter.ApplyE(cas,e);
263 
264  K->Fill(nbgamma) ;
265  for (Int_t j = 0; j < nbgamma; j++ ) {
266  projtot->Fill(e[j]);
267  } // j
268  } //i
269  watch.Stop();
270  cout << "Monte-Carlo time " << endl; watch.Print();
271 
272  projtot->Draw();
273 }
274 
294 void testGammaFilter2(const char *pathtodata, Int_t nbcas = 100, const char *filtername = "AGATA_FILTER.root")
295 {
296  TStopwatch watch;
297 
298  watch.Start();
299 
300  // Experimental conditions (from PACE like)
301  const Int_t nb_residu = 6; // 151DY, 152DY, 153DY, 151TB, 149GD, 148GD
302 
303  GEM nucleus[nb_residu]; // gem related
304 
305  TString filename[nb_residu]; // database related
306  TString nucleusname[nb_residu]; // database related
307 
308  Float_t sigma[nb_residu] = { 29.7, 47.1, 19.3, 19.3, 0.4, 1.31 };
309 
310  filename[0] = pathtodata; filename[0] += "DY_ENSDF.ens"; nucleusname[0] = "151DY";
311  filename[1] = pathtodata; filename[1] += "152DY.ags"; nucleusname[1] = "152DY";
312  filename[2] = pathtodata; filename[2] += "DY_ENSDF.ens"; nucleusname[2] = "153DY";
313  filename[3] = pathtodata; filename[3] += "TB_ENSDF.ens"; nucleusname[3] = "151TB";
314  filename[4] = pathtodata; filename[4] += "GD_ENSDF.ens"; nucleusname[4] = "149GD";
315  filename[5] = pathtodata; filename[5] += "148GD.ags"; nucleusname[5] = "148GD";
316 
317  // same prototype function to init all nucleus
318  Int_t ok = 0;
319  for (int i = 0; i < nb_residu; i++ ) {
320  ok += nucleus[i].Import(filename[i],nucleusname[i]);
321  // entry point 80 MeV, spin 60
322  nucleus[i].SetFeeding("TestFeeding"); nucleus[i].SetExI(80000,60);
323  }
324 
325  if ( ok != nb_residu ) {
326  cout << " A Level Scheme is missing !! " << endl; return;
327  }
328  // create the RandObj to select the nb_residu nuclei.
329  RandObj channel;
330  for (int i = 0; i < nb_residu; i++ ) channel.Add(nucleus+i,sigma[i]);
331 
332  // initiate the GammaFilter
333  Int_t nbgamma; Float_t e[1000]; // up to 1000 energies ..should be ok for the moment
334  GammaFilter filter; filter.InitFilter(filtername);
335  if ( filter.IsEffective() == false ) return ;
336 
337  // some spectra to be filled
338  TH1F *K = new TH1F("GammaFilter2_K","K",100,0,100);
339  TH1F *projtot = new TH1F("GammaFilter2","GammaFilter2",8192,0,4096);
340 
341  // open a root file to store the produced events in a ROOT TTree
342  TFile fileout("SimuDY.root","recreate");
343 
344  TTree *treeout;
345  treeout = (TTree *)fileout.Get("SimuDY");
346  if ( treeout == NULL ) {
347  treeout = new TTree("SimuDY","FromGammaFilter2");
348  treeout->Branch("mult",&nbgamma,"nbgamma/I"); treeout->Branch("e",e,"e[nbgamma]/F");
349  }
350 // else {
351 // treeout->SetBranchAddress("nbgamma",&nbgamma);
352 //
353 // }
354 
355  watch.Stop();
356  cout << "Initialisation time " << endl; watch.Print();
357  // start the Monte-Carlo
358 
359  watch.Reset(); watch.Start();
360  Cascade cas; BaseGEM *gem; GammaLink *gam; Int_t which;
361  for (Int_t i = 0; i < nbcas; i++ ) {
362 
363  gem = (BaseGEM *)channel.Rand(which); // select one nucleus
364  gem->DoCascade(cas); // determine a cascade on this nucleus
365 
366  // apply experimental filter on this cascade.
367  nbgamma = filter.ApplyE(cas,e);
368 
369  K->Fill(nbgamma) ;
370  for (Int_t j = 0; j < nbgamma; j++ ) {
371  projtot->Fill(e[j]);
372  } // j
373 
374  treeout->Fill();
375  } //i
376 
377  treeout->Write();
378 
379  watch.Stop();
380  cout << "Monte-Carlo time " << endl; watch.Print();
381 
382  projtot->Draw();
383 }
384 
385 
Gw::RandObj::Rand
virtual TObject * Rand()
it returns a pointer to a randomly selected object from the collection
Definition: RandObj.cpp:93
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
Gw::GEM
Class to get randomly cascades (with lateral feeding) of gammas on the basis of a level scheme...
Definition: GEM.h:55
Gw::GammaFilter::InitFilter
Bool_t InitFilter(const char *)
to modifiy the efficency curve.
Definition: GammaFilter.cpp:64
Gw::BaseGEM
Class to get randomly cascades of gammas on the basis of a level scheme.
Definition: BaseGEM.h:107
Gw::RandObj
A class to select randomly an object in a TObjArray of objects.
Definition: RandObj.h:55
Gw::GEM::SetFeeding
Link * SetFeeding(const char *feedingtype="TestFeeding")
Definition: GEM.cpp:47
testGEM
void testGEM(const char *pathtodata, Int_t nbcas=100)
Function to generate discret gamma-rays between several level schemes ...
Definition: FirstGEM.C:124
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)
Gw::GammaFilter
Standard filter for gamma-rays (EUROBALL-like appproach)
Definition: GammaFilter.h:41
testGammaFilter1
void testGammaFilter1(const char *pathtodata, Int_t nbcas=100, const char *filtername="AGATA_FILTER.root")
Function to generate tracked like events.
Definition: FirstGEM.C:206
testGammaFilter2
void testGammaFilter2(const char *pathtodata, Int_t nbcas=100, const char *filtername="AGATA_FILTER.root")
Function to generate tracked like events and store them in a ROOT TTree.
Definition: FirstGEM.C:294
Gw::BaseGEM::DoCascade
virtual Int_t DoCascade(TSeqCollection &cascade, Option_t *opt="")
cascade simulation
Definition: BaseGEM.cpp:523
testBaseGEM
void testBaseGEM(const char *pathtodata, Int_t nbcas=100)
Function to generate discret gamma-rays between several level schemes.
Definition: FirstGEM.C:49
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
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