Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/examples/extended/medical/dna/chem6/plotG_time.C

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

Diff markup

Differences between /examples/extended/medical/dna/chem6/plotG_time.C (Version 11.3.0) and /examples/extended/medical/dna/chem6/plotG_time.C (Version 11.1.2)


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