| Geant4 Cross Reference |
1 struct SpeciesInfoAOS
2 {
3 SpeciesInfoAOS()
4 {
5 fNEvent = 0;
6 fNumber = 0;
7 fG = 0.;
8 fG2 = 0.;
9 }
10
11 SpeciesInfoAOS(const SpeciesInfoAOS& right) // Species A(B);
12 {
13 fNEvent = right.fNEvent;
14 fNumber = right.fNumber;
15 fG = right.fG;
16 fG2 = right.fG2;
17 fName = right.fName;
18 }
19
20 SpeciesInfoAOS& operator=(const SpeciesInfoAOS& right) // A = B
21 {
22 if(&right == this) return *this;
23 fNEvent = right.fNEvent;
24 fNumber = right.fNumber;
25 fG = right.fG;
26 fG2 = right.fG2;
27 fName = right.fName;
28 return *this;
29 }
30
31 Int_t fNEvent;
32 Int_t fNumber;
33 Double_t fG;
34 Double_t fG2;
35 string fName;
36 };
37
38 //------------------------------------------------------------------------
39
40 struct SpeciesInfoSOA
41 {
42 SpeciesInfoSOA()
43 {
44 fRelatErr = 0;
45 }
46
47 SpeciesInfoSOA(const SpeciesInfoSOA& right) :
48 fG(right.fG),
49 fGerr(right.fGerr),
50 fTime(right.fTime),
51 fRelatErr(right.fRelatErr),
52 fName(right.fName)
53 {}
54
55 SpeciesInfoSOA& operator=(const SpeciesInfoSOA& right)
56 {
57 if(this == &right) return *this;
58 fG = right.fG;
59 fGerr = right.fGerr;
60 fTime = right.fTime;
61 fRelatErr = right.fRelatErr;
62 fName = right.fName;
63 return *this;
64 }
65
66 std::vector<Double_t> fG;
67 std::vector<Double_t> fGerr;
68 std::vector<Double_t> fTime;
69 Double_t fRelatErr;
70 Double_t fMin, fMax;
71 string fName;
72 };
73
74 const char* filetypes[] = {
75 "PostScript", "*.ps",
76 "Encapsulated PostScript", "*.eps",
77 "PDF files", "*.pdf",
78 "Gif files", "*.gif",
79 "PNG files", "*.png",
80 "All files", "*",
81 0, 0
82 };
83
84 TGTab *gTab = nullptr;
85
86 void Save()
87 {
88 TGFileInfo fi;
89 fi.fFileTypes = filetypes;
90
91 new TGFileDialog(gClient->GetRoot(),gClient->GetRoot(),kFDSave,&fi);
92 gROOT->GetListOfCanvases()->At(gTab->GetCurrent())->SaveAs(fi.fFilename);
93 }
94
95
96 void plotG_time()
97 {
98 gROOT->SetStyle("Plain");
99 gStyle->SetPalette(1);
100 gStyle->SetCanvasBorderMode(0);
101 gStyle->SetFrameBorderMode(0);
102 gStyle->SetPadTickX(1);
103 gStyle->SetPadTickY(1);
104
105 TGMainFrame *main = new TGMainFrame(gClient->GetRoot(), 200, 200);
106 gTab = new TGTab(main, 200, 200);
107 Double_t timeA, sumG, sumG2;
108 Int_t speciesID, number, nEvent;
109 char speciesName[500];
110
111 TFile* file = new TFile;
112 file = TFile::Open("Species0.root");
113
114 TTree* tree = (TTree*)file->Get("species");
115 tree->SetBranchAddress("speciesID", &speciesID);
116 tree->SetBranchAddress("number", &number);
117 tree->SetBranchAddress("nEvent", &nEvent);
118 tree->SetBranchAddress("speciesName", &speciesName);
119 tree->SetBranchAddress("time", &timeA);
120 tree->SetBranchAddress("sumG", &sumG);
121 tree->SetBranchAddress("sumG2", &sumG2);
122
123 Long64_t nentries = tree->GetEntries();
124
125 if(nentries == 0) {
126 cout << "No entries found in the tree species contained in the file "
127 << file->GetPath() << endl;
128 exit(1);
129 }
130
131 std::map<int, std::map<double, SpeciesInfoAOS>> speciesTimeInfo;
132
133 for (Int_t j=0; j < nentries; j++) {
134 tree->GetEntry(j);
135
136 SpeciesInfoAOS& infoAOS = speciesTimeInfo[speciesID][timeA];
137
138 infoAOS.fNumber += number;
139 infoAOS.fG += sumG;
140 infoAOS.fG2 += sumG2;
141 infoAOS.fNEvent += nEvent;
142 infoAOS.fName = speciesName;
143 }
144
145 std::map<Int_t, SpeciesInfoSOA> speciesInfo;
146
147 auto it_SOA = speciesTimeInfo.begin();
148 auto end_SOA = speciesTimeInfo.end();
149
150 for (; it_SOA!=end_SOA;++it_SOA) {
151 const Int_t _speciesID = it_SOA->first;
152 SpeciesInfoSOA& info = speciesInfo[_speciesID];
153
154 auto it2 = it_SOA->second.begin();
155 auto end2 = it_SOA->second.end();
156
157 info.fName = it2->second.fName;
158 const size_t size2 = it_SOA->second.size();
159 info.fG.resize(size2);
160 info.fGerr.resize(size2);
161 info.fTime.resize(size2);
162
163 Int_t color = (2+_speciesID)%TColor::GetNumberOfColors();
164 if(color == 5 || color == 10 || color == 0) ++color;
165
166 for (int i2 = 0 ;it2!=end2;++it2, ++i2) {
167 SpeciesInfoAOS& infoAOS = it2->second;
168
169 Double_t _SumG2 = infoAOS.fG2;
170 Double_t _MeanG = infoAOS.fG/infoAOS.fNEvent;
171 Double_t _Gerr = (infoAOS.fNEvent > 1) ? sqrt((_SumG2/infoAOS.fNEvent - pow(_MeanG,2))
172 /(infoAOS.fNEvent-1) ) : 0.;
173
174 info.fG[i2] = _MeanG;
175 info.fGerr[i2] = _Gerr;
176 info.fTime[i2] = it2->first;
177
178 info.fRelatErr += _Gerr/(_MeanG + 1e-30);
179 }
180
181 TGraphErrors* gSpecies = new TGraphErrors(info.fG.size(),
182 info.fTime.data(),
183 info.fG.data(),
184 0,
185 info.fGerr.data());
186
187 TGCompositeFrame *tf = gTab->AddTab(info.fName.c_str());
188 TGCompositeFrame *frame = new TGCompositeFrame(tf, 60, 60,
189 kHorizontalFrame);
190
191 tf->AddFrame(frame, new TGLayoutHints(kLHintsExpandX|kLHintsExpandY,
192 10,10,10,2));
193
194 TRootEmbeddedCanvas *c1 = new TRootEmbeddedCanvas(info.fName.c_str(),
195 frame, 700, 500);
196 frame->AddFrame(c1, new TGLayoutHints(kLHintsExpandX|kLHintsExpandY,
197 10,10,10,2));
198 c1->GetCanvas()->SetLogx();
199
200 TGHorizontalFrame* hframe = new TGHorizontalFrame(tf, 200, 40);
201
202 TGTextButton* save = new TGTextButton(hframe, "&Save as ...",
203 "Save()");
204 hframe->AddFrame(save, new TGLayoutHints(kLHintsCenterX, 5, 5, 3, 4));
205
206 TGTextButton *exit = new TGTextButton(hframe, "&Exit ",
207 "gApplication->Terminate()");
208 hframe->AddFrame(exit, new TGLayoutHints(kLHintsCenterX, 5, 5, 3, 4));
209
210 tf->AddFrame(hframe, new TGLayoutHints(kLHintsCenterX, 2, 2, 2, 2));
211
212 gSpecies->SetTitle(info.fName.c_str());
213 gSpecies->SetMarkerStyle(20+_speciesID);
214 gSpecies->SetMarkerColor(color);
215 gSpecies->SetLineColor(color);
216 gSpecies->GetXaxis()->SetTitle("Time (ns)");
217 gSpecies->GetXaxis()->SetTitleOffset(1.1);
218 gSpecies->GetYaxis()->SetTitle("G value (molecules/100 eV)");
219 gSpecies->GetYaxis()->SetTitleOffset(1.2);
220 gSpecies->Draw("ALP");
221 }
222
223 main->AddFrame(gTab, new TGLayoutHints(kLHintsBottom | kLHintsExpandX |
224 kLHintsExpandY, 2, 2, 5, 1));
225
226 main->MapSubwindows();
227 main->Resize(); // resize to default size
228 main->MapWindow();
229
230 }
231