GammaWare  Head Version for release 0.9
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
FirstPax.C
Go to the documentation of this file.
1 #ifndef __CINT__
2 
3 #ifndef ROOT_TObjArray
4 #include "TObjArray.h"
5 #endif
6 #ifndef ROOT_TH1
7 #include "TH1.h"
8 #endif
9 #ifndef ROOT_TH2
10 #include "TH2.h"
11 #endif
12 #ifndef ROOT_TMultiGraph
13 #include "TMultiGraph.h"
14 #endif
15 #ifndef ROOT_TGraphErrors
16 #include "TGraphErrors.h"
17 #endif
18 #ifndef ROOT_TLegend
19 #include "TLegend.h"
20 #endif
21 #ifndef ROOT_TLeaf
22 #include "TLeaf.h"
23 #endif
24 #ifndef ROOT_TLine
25 #include "TLine.h"
26 #endif
27 #ifndef ROOT_TPad
28 #include "TPad.h"
29 #endif
30 #ifndef ROOT_TMarker
31 #include "TMarker.h"
32 #endif
33 #ifndef GW_EVAPLMOF_H
34 #include "EvapLMOF.h"
35 #endif
36 
37 #endif
38 
39 #include "EvapLMOF.h"
40 
41 const Int_t MAXZ = 120; const Int_t MAXN = 200; // upper limit for the number of protons and neutrons for one nucleus
42 
43 using namespace Gw;
44 
52 void Scan(const char *fname)
53 {
54  EvapLMOF pax;
55  pax.AddFile(fname); pax.VERBOSE = 1; pax << idonly; pax.Scan(0,100); pax.ShowCurrentConditions();
56 }
57 
67 TH2D *Production2D(const char *hname, const char *fname, const char *pathtofile = "./",bool fis = false)
68 {
69  // allocate 2D plot
70  TH2D *h = new TH2D(hname,"Cross Section",MAXN,0,MAXN,MAXZ,0,MAXZ);
71 
72  // open chain of files
73  EvapLMOF pax; pax.AddFile(fname); pax.SetPath(pathtofile); pax << idonly;
74 
75  while( pax.NextEvent() ) { // loop on events in the generated pax file
76  if ( pax.GetId().FIS == fis ) {
77  h->Fill(pax.GetId().Z,pax.GetId().N);
78  }
79  } // while reader
80 
81  h->Fill(pax.GetGlobal().ZC,pax.GetGlobal().AC-pax.GetGlobal().ZC,pax.GetGlobal().NCASC);
82 
83  // now scale to get the cross sections
84  h->Scale(pax.GetGlobal().SIGMA/pax.GetGlobal().NCASC);
85 
86  pax.ShowCurrentConditions();
87 
88  h->Draw("col"); return h;
89 }
90 
103 TMultiGraph *Excitation(const char *filebasename, const char *pathtofile = "./", Int_t istart = 1, Int_t iend = 10, Float_t threshold = 0.01)
104 {
105  // style for each graph ... colmarq_Z color for an isotope, stymarq_N style for isotone
106  // firstZ firstN indicate the starting point in (Z,N) .. could be the compound nucleus
107  // so the indice #i in colmarq_Z (stymarq_N) corresponds to the nucleux with firstZ-i (firstN-i) protons (neutrons)
108  // Fission have a different style as well as nuclei 'out of range'
109  const Short_t NBSTYLE = 11;
110  Short_t firstZ = 64;
111  Short_t firstN = 88;
112  const Short_t colmarq_Z[NBSTYLE] = { 2, 3, 4, 6, 7, 8, 9,20,30,40,50}; const Short_t colmarq_FISSION = 1;
113  const Short_t stymarq_N[NBSTYLE] = {20,21,22,23,24,25,26,27,28,29,30}; const Short_t stymarq_FISSION = 5;
114 
115  const Int_t nbin = 2000, min = 0, max = 200; // histograms definitions
116 
117  // array of histograms that are transformed to TGraph at the end.
118  TObjArray ar(MAXN*MAXZ); ar.SetOwner(); TH1F *h; TString tmp;
119 
120  // open chain of files that are in the pathtofile directory
121  EvapLMOF pax; pax.SetPath(pathtofile); pax << idonly;
122 
123  // add files to the pax reader
124  for (Int_t i = istart; i < iend; i++ ) {
125  tmp = filebasename;
126  if ( i < 10 )
127  tmp += "00";
128  if ( 10 <= i && i < 100 )
129  tmp += "0";
130 
131  tmp += i; tmp += ".pax";
132  if ( pax.VERBOSE ) cout << "Adding " << tmp.Data() << endl; pax.AddFile(tmp.Data());
133  }
134 
135  while( pax.NextEvent() ) { // loop on events in the generated pax file
136  if ( pax.GetId().FIS == kTRUE ) { // slot 0 of the TObjArray to store fission like events
137  if ( ar[0] == NULL ) {
138  tmp ="Fission"; h = new TH1F(tmp.Data(),filebasename,nbin,min,max); ar.AddAt(h,0);
139 
140  // set attributes
141  h->SetMarkerColor(colmarq_FISSION); h->SetMarkerStyle(stymarq_FISSION); h->SetLineColor(colmarq_FISSION);
142  }
143  h = (TH1F *)(ar.At(0));
144  if ( h ) h->Fill(pax.GetGlobal().ENER,pax.GetGlobal().SIGMA/pax.GetGlobal().NCASC);
145  }
146  else { // residus
147  if ( pax.GetId().N > 0 && pax.GetId().Z > 0 ) {
148  if ( ar[pax.GetId().Z * MAXN + pax.GetId().N] == NULL ) {
149  tmp = pax.GetId().Z; tmp += "-"; tmp += pax.GetId().N;
150  h = new TH1F(tmp.Data(),filebasename,nbin,min,max); ar.AddAt(h,pax.GetId().Z * MAXN + pax.GetId().N);
151 
152  // set histo style color and marker
153  if ( (firstN - pax.GetId().N) > -1 && (firstN - pax.GetId().N) < NBSTYLE
154  && (firstZ - pax.GetId().Z) > -1 && (firstZ - pax.GetId().Z) < NBSTYLE ) {
155  h->SetMarkerColor(colmarq_Z[firstZ - pax.GetId().Z]);
156  h->SetMarkerStyle(stymarq_N[firstN - pax.GetId().N]);
157  h->SetLineColor(colmarq_Z[firstZ - pax.GetId().Z]);
158  }
159  else {
160  h->SetMarkerColor(14);
161  h->SetMarkerStyle(2);
162  h->SetLineColor(14);
163  h->SetLineStyle(2);
164  }
165  }
166  h = (TH1F *)(ar.At(pax.GetId().Z * MAXN + pax.GetId().N));
167  if ( h ) h->Fill(pax.GetGlobal().ENER,pax.GetGlobal().SIGMA/pax.GetGlobal().NCASC);
168  }
169  }
170  } // while reader
171  pax.ShowCurrentConditions();
172 
173  // to put all histograms in graphs
174  TMultiGraph *stack = new TMultiGraph(filebasename,filebasename);
175 
176  TLegend *leg = new TLegend(0.8,0.8,0.99,0.99);
177 
178  Short_t ihisto = 0;
179  if ( ar.GetEntries() > 0 ) { // at least one histo exists
180  for (Int_t i = 0; i < MAXN*MAXZ; i++ ) {
181  if ( ar[i] ) { // an histogram exists at this slot
182 
183  h = (TH1F *)ar[i]; if ( h->GetMaximum() < threshold ) continue;
184 
185  // new graph added to the collection
186  TGraphErrors *g = new TGraphErrors();
187  stack->Add(g);
188 
189  // to set different style for each graph and the legend
190  g->SetMarkerColor(h->GetMarkerColor()); g->SetMarkerStyle(h->GetMarkerStyle());
191  g->SetLineColor(h->GetLineColor());
192  g->SetLineStyle(h->GetLineStyle());
193 
194  g->SetNameTitle(h->GetName(),h->GetTitle()); // to get the histogram that is used to build the TGraph
195 
196  TString tmp;
197  if ( i == 0 ) tmp = "Fission";
198  else {
199  tmp ="X ^{"; tmp += (i / MAXN + i % MAXN) ; tmp +="}_{"; tmp += i / MAXN; tmp +="}";
200  }
201  leg->AddEntry(g,tmp.Data(),"L P");
202 
203  // fill the graph
204  Int_t ipoint = 0;
205  for (Int_t ibin = 1; ibin < nbin+1; ibin++ ) { // keep bincontent > threshold
206  if ( h->GetBinContent(ibin) > threshold ) {
207  g->SetPoint(ipoint,h->GetBinCenter(ibin),h->GetBinContent(ibin));
208  g->SetPointError(ipoint,h->GetBinWidth(ibin),h->GetBinError(ibin));
209  ipoint++;
210  }
211  }
212  ihisto++;
213  }
214  }
215  stack->Draw("ALP");
216  stack->GetXaxis()->SetTitle("Excitation energy (MeV)"); stack->GetYaxis()->SetTitle("Cross section (mbarns)");
217  leg->Draw();
218  }
219  return stack;
220 }
221 
233 void Distribution(TH1F &histo, const char *filebasename, const char *pathtofile = "./", Int_t istart = 1, Int_t iend = 10)
234 {
235  TString tmp;
236 
237  // open chain of files that are in the pathtofile directory
238  EvapLMOF pax; pax.SetPath(pathtofile); pax << full; // pax << idonly; pax.VERBOSE = 1;
239 
240  // add files to the pax reader
241  for (Int_t i = istart; i < iend; i++ ) {
242  tmp = filebasename;
243  if ( i < 10 )
244  tmp += "00";
245  if ( 10 <= i && i < 100 )
246  tmp += "0";
247 
248  tmp += i; tmp += ".pax";
249  if ( pax.VERBOSE ) cout << "Adding " << tmp.Data() << endl; pax.AddFile(tmp.Data());
250  }
251 
252  TLeaf *x = pax.GetLeaf(histo.GetTitle()), *w = pax.GetLeaf("WT");
253  if ( x == NULL )
254  { cout << " Unknown distribution " << tmp.Data() << endl; return; }
255 
256  // fill the histogram and draw it
257  while( pax.NextEvent() ) { // loop on events in the generated pax file
258  histo.Fill(x->GetValue(),w->GetValue());
259  }
260  pax.ShowCurrentConditions(); histo.Draw(); delete x; delete w;
261 }
262 
267 TGraph *Evolution1D(const char *var, const char *filebasename, const char *pathtofile = "./", Int_t istart = 1, Int_t iend = 10, Int_t which = 0)
268 {
269  TString tmp;
270 
271  // open chain of files that are in the pathtofile directory
272  EvapLMOF pax; pax.SetPath(pathtofile); pax << full;
273 
274  // add files to the pax reader
275  for (Int_t i = istart; i < iend; i++ ) {
276  tmp = filebasename;
277  if ( i < 10 )
278  tmp += "00";
279  if ( 10 <= i && i < 100 )
280  tmp += "0";
281 
282  tmp += i; tmp += ".pax";
283  if ( pax.VERBOSE ) cout << "Adding " << tmp.Data() << endl; pax.AddFile(tmp.Data());
284  }
285 
286  // to get the variable to be inspected
287  TLeaf *x = pax.GetLeaf(var);
288  if ( x == NULL )
289  { cout << " Unknown variable " << tmp.Data() << endl; return NULL; }
290 
291  // fill the graph and draw it
292  TGraph *gr = new TGraph();
293  Int_t i = 0;
294  while( pax.NextEvent() ) { // loop on events in the generated pax file
295  /* if ( x->GetValue(which) > 1 ) pax.ShowCurrentEvent() ; */ gr->SetPoint(i,i,x->GetValue()); i++;
296  }
297  pax.ShowCurrentConditions(); gr->Draw("A*"); delete x; return gr;
298 }
299 
300 
309 void Distribution2D(TH2 &histo, const char *filebasename, const char *pathtofile = "./", Int_t istart = 1, Int_t iend = 10)
310 {
311  TString tmp;
312 
313  // open chain of files that are in the pathtofile directory
314  EvapLMOF pax; pax.SetPath(pathtofile); pax << full; // pax << idonly; pax.VERBOSE = 1;
315 
316  // add files to the pax reader
317  for (Int_t i = istart; i < iend; i++ ) {
318  tmp = filebasename;
319  if ( i < 10 )
320  tmp += "00";
321  if ( 10 <= i && i < 100 )
322  tmp += "0";
323 
324  tmp += i; tmp += ".pax";
325  if ( pax.VERBOSE ) cout << "Adding " << tmp.Data() << endl; pax.AddFile(tmp.Data());
326  }
327 
328  // to get x, y from the histogram's name
329  TLeaf *x = NULL, *y = NULL, *w = pax.GetLeaf("WT");
330 
331  tmp = histo.GetTitle(); TObjArray *arr = tmp.Tokenize(":");
332  if ( arr && arr->GetEntries() == 2 ) { // look for what is to be drawn
333  x = pax.GetLeaf(arr->At(0)->GetName());
334  y = pax.GetLeaf(arr->At(1)->GetName());
335 
336  delete arr;
337  }
338 
339  if ( x == NULL || y == NULL )
340  { cout << " Unknown distribution " << tmp.Data() << endl; return; }
341 
342  // fill the histogram and draw it
343  while( pax.NextEvent() ) { // loop on events in the generated pax file
344 // if ( pax.GetId().FIS == false ) histo.Fill(x->GetValue(),y->GetValue(),w->GetValue());
345  histo.Fill(x->GetValue(),y->GetValue(),w->GetValue());
346  }
347  pax.ShowCurrentConditions(); histo.Draw(); delete x; delete y; delete w;
348 }
349 
353 void EPxTHETA(TH2 &histo, const char *filebasename, const char *pathtofile = "./", Int_t istart = 1, Int_t iend = 10, Int_t mode = 1)
354 {
355  TString tmp;
356 
357  // open chain of files that are in the pathtofile directory
358  EvapLMOF pax; pax.SetPath(pathtofile); pax << full; // pax << idonly; pax.VERBOSE = 1;
359 
360  // add files to the pax reader
361  for (Int_t i = istart; i < iend; i++ ) {
362  tmp = filebasename;
363  if ( i < 10 )
364  tmp += "00";
365  if ( 10 <= i && i < 100 )
366  tmp += "0";
367 
368  tmp += i; tmp += ".pax";
369  if ( pax.VERBOSE ) cout << "Adding " << tmp.Data() << endl; pax.AddFile(tmp.Data());
370  }
371 
372  // to get x, y from the histogram's name
373  TLeaf *x = pax.GetLeaf("THETA"), *y = pax.GetLeaf("EEJ"), *w = pax.GetLeaf("WT");
374 
375 
376  // fill the histogram and draw it
377  while( pax.NextEvent() ) { // loop on events in the generated pax file
378  for( Short_t i = 0; i < pax.GetId().NUM ; i++){
379  if ( pax.GetData().MODE[i] == mode ) histo.Fill(x->GetValue(i)*180/TMath::Pi(),y->GetValue(i),w->GetValue());
380  }
381  }
382  pax.ShowCurrentConditions(); histo.Draw(); delete x; delete y; delete w;
383 }
384 
390 void ShowCascades(const char *filebasename, const char *pathtofile = "./", Int_t istart = 1, Int_t iend = 10)
391 {
392  TString tmp; TLine line; line.SetLineWidth(2); TMarker mark; mark.SetMarkerStyle(23);
393 
394  Short_t colors[11] = {0,1,30,28,2,4,6,8,9,10,11};
395  Short_t styles[11] = {0,2,1,1,1,1,1,8,9,10,11};
396 
397  if ( gPad == NULL ) return;
398 
399  // open chain of files that are in the pathtofile directory
400  EvapLMOF pax; pax.SetPath(pathtofile); pax << full; // pax << idonly; pax.VERBOSE = 1;
401 
402  // add files to the pax reader
403  for (Int_t i = istart; i < iend; i++ ) {
404  tmp = filebasename;
405  if ( i < 10 )
406  tmp += "00";
407  if ( 10 <= i && i < 100 )
408  tmp += "0";
409 
410  tmp += i; tmp += ".pax";
411  if ( pax.VERBOSE ) cout << "Adding " << tmp.Data() << endl; pax.AddFile(tmp.Data());
412  }
413 
414  Int_t refresh = 0;
415  while( pax.NextEvent() ) { // loop on events in the generated pax file
416 
417  if ( pax.GetId().FIS == kTRUE ) continue;
418 
419  if ( refresh % 3 == 0 )
420  { refresh = 0 ; gPad->Clear(); gPad->DrawFrame(0,0,TMath::Max(1.1*pax.GetId().L+5,10.),1.1*pax.GetGlobal().ENER); }
421  refresh++;
422 
423  mark.DrawMarker(pax.GetId().JE,pax.GetId().EE); // entry point as given by the id
424 
425  for( Int_t i = 0; i < pax.GetId().NUM ; i++){ // to draw the full cascade
426  line.SetLineColor(colors[pax.GetData().MODE[i]]);
427  line.SetLineStyle(styles[pax.GetData().MODE[i]]);
428  if ( i == pax.GetId().NUM - 1 )
429  line.DrawLine(pax.GetData().JI[i],pax.GetData().EI[i],0,0);
430  else
431  line.DrawLine(pax.GetData().JI[i],pax.GetData().EI[i],pax.GetKin().JF[i],pax.GetData().EI[i+1]);
432  }
433  gPad->Update(); TObject *obj = gPad->WaitPrimitive();
434 
435  // pax.ShowCurrentEvent();
436  if ( gROOT->IsInterrupted() ) { cout << "loop stopped " << endl; break; }
437  }
438 }
439 
443 void Distribution(const char *filebasename, const char *pathtofile = "./", Int_t istart = 1, Int_t iend = 10)
444 {
445  TH1F *histo = new TH1F("","",200,-100,100);
446 
447  TString tmp;
448 
449  // open chain of files that are in the pathtofile directory
450  EvapLMOF pax; pax.SetPath(pathtofile); pax << full; // pax << idonly; pax.VERBOSE = 1;
451 
452  // add files to the pax reader
453  for (Int_t i = istart; i < iend; i++ ) {
454  tmp = filebasename;
455  if ( i < 10 )
456  tmp += "00";
457  if ( 10 <= i && i < 100 )
458  tmp += "0";
459 
460  tmp += i; tmp += ".pax";
461  if ( pax.VERBOSE ) cout << "Adding " << tmp.Data() << endl; pax.AddFile(tmp.Data());
462  }
463 
464  // fill the histogram and draw it
465  while( pax.NextEvent() ) { // loop on events in the generated pax file
466  // if ( pax.GetId().L-pax.GetData().JI[0] > -1 ) {
467  histo->Fill(pax.GetId().L-pax.GetData().JI[0]);
468  // pax.ShowCurrentEvent();}
469  }
470  pax.ShowCurrentConditions(); histo->Draw();
471 }
472 
473 
Gw::EvapLMOF::AddFile
void AddFile(const char *fname)
to add an input (pax) file
Definition: EvapLMOF.h:219
Gw::EvapLMOF
Interface to read/decode EvapOR List Mode Output Files (so-called .pax files)
Definition: EvapLMOF.h:90
Scan
void Scan(const char *fname)
Function to scan a pax file and display the content on the standard ouput.
Definition: FirstPax.C:52
__LoadWatchers__::tmp
TString tmp
Definition: GANILLoadWatchers.C:47
Excitation
TMultiGraph * Excitation(const char *filebasename, const char *pathtofile="./", Int_t istart=1, Int_t iend=10, Float_t threshold=0.01)
Function to plot Excitation function.
Definition: FirstPax.C:103
ShowCascades
void ShowCascades(const char *filebasename, const char *pathtofile="./", Int_t istart=1, Int_t iend=10)
Function to show cascades in the E x I plan.
Definition: FirstPax.C:390
Gw::full
EvapLMOF & full(EvapLMOF &reader)
Definition: EvapLMOF.h:295
histo
TH1F * histo[MaxValue]
Definition: ReadDaqAlone.C:31
Gw::EvapLMOF::ShowCurrentConditions
void ShowCurrentConditions(std::ostream &out=std::cout)
The current status is shown on the standard ouptut.
Definition: EvapLMOF.cpp:242
Gw::EvapLMOF::GetLeaf
TLeaf * GetLeaf(const char *leafname="no")
To get a ROOT leaf corresponding to one of the variable.
Definition: EvapLMOF.cpp:633
Evolution1D
TGraph * Evolution1D(const char *var, const char *filebasename, const char *pathtofile="./", Int_t istart=1, Int_t iend=10, Int_t which=0)
Function to draw a single distribution.
Definition: FirstPax.C:267
Gw::EvapLMOF::SetPath
void SetPath(const char *path="./")
set the directory where to find the input files
Definition: EvapLMOF.h:222
Gw::EvapLMOF::Scan
void Scan(UInt_t start=0, UInt_t end=10)
To scan events and show them on the standard output.
Definition: EvapLMOF.cpp:615
Gw::EvapLMOF::GetData
const DATA & GetData()
To get the current data structure.
Definition: EvapLMOF.h:258
Distribution
void Distribution(TH1F &histo, const char *filebasename, const char *pathtofile="./", Int_t istart=1, Int_t iend=10)
Function to draw a single distribution.
Definition: FirstPax.C:233
Gw::EvapLMOF::GetGlobal
const GLOBAL & GetGlobal()
To get the current global structure.
Definition: EvapLMOF.h:254
MAXN
const Int_t MAXN
Definition: FirstPax.C:41
MAXZ
const Int_t MAXZ
Definition: FirstPax.C:41
Distribution2D
void Distribution2D(TH2 &histo, const char *filebasename, const char *pathtofile="./", Int_t istart=1, Int_t iend=10)
Function to draw a 2D distribution.
Definition: FirstPax.C:309
Gw::EvapLMOF::GetId
const ID & GetId()
To get the current id structure.
Definition: EvapLMOF.h:256
endl
ADF::LogMessage & endl(ADF::LogMessage &log)
Definition: ADFLogMessage.cpp:265
Production2D
TH2D * Production2D(const char *hname, const char *fname, const char *pathtofile="./", bool fis=false)
Function to plot a 2D plot of all the produced residus.
Definition: FirstPax.C:67
Gw::EvapLMOF::VERBOSE
static Short_t VERBOSE
Set verbosity level.
Definition: EvapLMOF.h:97
Gw::EvapLMOF::GetKin
const KINMAT & GetKin()
To get the current kinematic structure.
Definition: EvapLMOF.h:260
Gw::idonly
EvapLMOF & idonly(EvapLMOF &reader)
Definition: EvapLMOF.h:294
Gw::EvapLMOF::NextEvent
Bool_t NextEvent(UInt_t next=0)
To load the next event.
Definition: EvapLMOF.cpp:572
EPxTHETA
void EPxTHETA(TH2 &histo, const char *filebasename, const char *pathtofile="./", Int_t istart=1, Int_t iend=10, Int_t mode=1)
Function to draw the Energy of the particule as a function of the theta angle.
Definition: FirstPax.C:353