readExfor.C
Go to the documentation of this file.
1 #include "Riostream.h"
2 
3 /**
4  * Example script for reading datafile.
5  * Usage:
6  * root -l RootScripts/iaeaBenchmark/readExfor.C
7  */
8 void readExfor() {
9  TString dir = gSystem->UnixPathName(gInterpreter->GetCurrentMacroName());
10  dir.ReplaceAll("basic.C","");
11  dir.ReplaceAll("/./","/");
12  ifstream in;
13  //in.open(Form("../data/fragmentEnergySpctra279mmWater0deg.dat",dir.Data()));
14  in.open(Form("experimentalData/fragmentEnergySpctra279mmWater0deg.dat",dir.Data()));
15  Float_t f1,f2,f3, f4,f5,f6;
16  Int_t nlines = 0;
17  TFile *f = new TFile("basic.root","RECREATE");
18  TH1F *h1 = new TH1F("h1","x distribution",100,-4,4);
19  TNtuple *ntuple = new TNtuple("ntuple","data from ascii file","Energy:He:B:H:Li:Be");
20 
21  Char_t DATAFLAG[4];
22  Int_t NDATA;
23  Char_t n1[6], n2[2], n3[2], n4[2], n5[2], n6[2],;
24  in >> DATAFLAG >> NDATA ; // Read EXFOR line: 'DATA 6'
25  in >> n1 >> n2 >> n3 >> n4 >> n5 >> n6; // Read column titles: 'Energy He B [...]'
26 
27  cout <<n1<<" "<<n2<<" "<<n3<<" "<<n4<<" "<<n5<<" "<<n6<<"\n";
28  while (1) {
29  in >> f1 >> f2 >> f3 >>f4 >> f5 >> f6;
30  if (!in.good()) break;
31  if (nlines < 500 ) printf("%i %0.2f %0.2f %0.2f %0.2f %0.2f \n",f1,f2,f3,f4,f5,f6);
32  //h1->Fill(f1);
33  ntuple->Fill(f1,f2,f3,f4,f5,f6);
34  nlines++;
35  }
36  ntuple->SetMarkerStyle(5);
37  ntuple->Draw("He:Energy","","l");
38  ntuple->Draw("B:Energy","","l,Same");
39  ntuple->Draw("H:Energy","","l,Same");
40  ntuple->Draw("Li:Energy","","l,Same");
41  ntuple->Draw("Be:Energy","","l,Same");
42  printf(" found %d points\n",nlines);
43  //Let's pull in the monte carlo analysis results
44 
45  TCanvas *mc = new TCanvas("c2", "Simulation");
46  TFile *MCData = TFile::Open("IAEA.root");
47  TH1F* MC_helium = (TH1F*)MCData->Get("heliumEnergyAfterPhantom");
48  TH1F* MC_hydrogen = (TH1F*)MCData->Get("hydrogenEnergyAfterPhantom");
49  //scale and plot
50  TNtuple *fragments = (TNtuple*) MCData->Get("fragmentNtuple");
51 
52  ScaleHelium = 1/(MC_helium->Integral());
53  ScaleHydrogen = 1/(MC_hydrogen->Integral());
54  //x should also be scaled to per nucleon
55 
56  MC_helium->Scale(ScaleHelium);
57  MC_helium->Draw("");
58  printf("Scaled helium by %.9f\n",ScaleHelium);
59 
60  MC_hydrogen->Scale(ScaleHydrogen);
61  MC_hydrogen->SetLineColor(kRed);
62  MC_hydrogen->Draw("Same");
63  printf("Scaled hydrogen by %.9f\n",ScaleHydrogen);
64 
65  TCanvas *fc = new TCanvas("fc", "Fragments");
66  fragments->Draw("energy");
67 
68  in.close();
69 
70  f->Write();
71 }
Float_t f4
enum BeamMode kRed
TTree * ntuple
Definition: macro.C:6
Float_t f2
nlines
Float_t f3
printf("%d Experimental points found\n", nlines)
Float_t f1
void readExfor()
Definition: readExfor.C:8
OStream cout
Definition: OStream.cxx:6
TH1F * h1
ifstream in
Definition: comparison.C:7
TDirectory * dir
Definition: macro.C:5