dna/UHDR/plotG_time.C

0001 // from chem6 example
0002 struct SpeciesInfoAOS {
0003   SpeciesInfoAOS() {
0004     fNEvent = 0;
0005     fNumber = 0;
0006     fG = 0.;
0007     fG2 = 0.;
0008   }
0009
0010   SpeciesInfoAOS(const SpeciesInfoAOS &right) // Species A(B);
0011   {
0012     fNEvent = right.fNEvent;
0013     fNumber = right.fNumber;
0014     fG = right.fG;
0015     fG2 = right.fG2;
0016     fName = right.fName;
0017   }
0018
0019   SpeciesInfoAOS &operator=(const SpeciesInfoAOS &right) // A = B
0020   {
0021     if (&right == this) return *this;
0022     fNEvent = right.fNEvent;
0023     fNumber = right.fNumber;
0024     fG = right.fG;
0025     fG2 = right.fG2;
0026     fName = right.fName;
0027     return *this;
0028   }
0029
0030   Int_t fNEvent;
0031   Int_t fNumber;
0032   Double_t fG;
0033   Double_t fG2;
0034   string fName;
0035 };
0036
0037 //------------------------------------------------------------------------
0038
0039 struct SpeciesInfoSOA {
0040   SpeciesInfoSOA() {
0041     fRelatErr = 0;
0042   }
0043
0044   SpeciesInfoSOA(const SpeciesInfoSOA &right) :
0045       fG(right.fG),
0046       fGerr(right.fGerr),
0047       fTime(right.fTime),
0048       fRelatErr(right.fRelatErr),
0049       fName(right.fName) {}
0050
0051   SpeciesInfoSOA &operator=(const SpeciesInfoSOA &right) {
0052     if (this == &right) return *this;
0053     fG = right.fG;
0054     fGerr = right.fGerr;
0055     fTime = right.fTime;
0056     fRelatErr = right.fRelatErr;
0057     fName = right.fName;
0058     return *this;
0059   }
0060
0061   std::vector <Double_t> fG;
0062   std::vector <Double_t> fGerr;
0063   std::vector <Double_t> fTime;
0064   Double_t fRelatErr;
0065   Double_t fMin, fMax;
0066   string fName;
0067 };
0068
0069 const char *filetypes[] = {
0070     "PostScript", "*.ps",
0071     "Encapsulated PostScript", "*.eps",
0072     "PDF files", "*.pdf",
0073     "Gif files", "*.gif",
0074     "PNG files", "*.png",
0075     "All files", "*",
0076     0, 0
0077 };
0078
0079 TGTab *gTab = nullptr;
0080
0081 void Save()
0082 {
0083   TGFileInfo fi;
0084   fi.fFileTypes = filetypes;
0085
0086   new TGFileDialog(gClient->GetRoot(),gClient->GetRoot(),kFDSave,&fi);
0087   gROOT->GetListOfCanvases()->At(gTab->GetCurrent())->SaveAs(fi.fFilename);
0088 }
0089
0090 void plotG_time() {
0091
0092   gROOT->SetStyle("Plain");
0093   gStyle->SetPalette(1);
0094   gStyle->SetCanvasBorderMode(0);
0095   gStyle->SetFrameBorderMode(0);
0096   gStyle->SetPadTickX(1);
0097   gStyle->SetPadTickY(1);
0098
0099   TGMainFrame *main = new TGMainFrame(gClient->GetRoot(), 200, 200);
0100   TGTab *fTab = new TGTab(main, 200, 200);
0101
0102   Double_t timeA, sumG, sumG2;
0103   Double_t sumDose, sumDose2;
0104   Double_t sumDoseRate, sumDoseRate2;
0105   Int_t speciesID, number, nEvent;
0106   char speciesName[500];
0107
0108
0109   TFile *file = new TFile;
0110   file = TFile::Open("Dose_0.root");
0111
0112   TDirectoryFile *dir = dynamic_cast<TDirectoryFile *>(file->Get("ntuple"));
0113   TTree *tree = (TTree *) dir->Get("Gvalue");
0114   tree->SetBranchAddress("speciesID", &speciesID);
0115   tree->SetBranchAddress("number", &number);
0116   tree->SetBranchAddress("nEvent", &nEvent);
0117   tree->SetBranchAddress("speciesName", &speciesName);
0118   tree->SetBranchAddress("time", &timeA);
0119   tree->SetBranchAddress("sumG", &sumG);
0120   tree->SetBranchAddress("sumG2", &sumG2);
0121   tree->SetBranchAddress("TotalDose", &sumDose);
0122   tree->SetBranchAddress("TotalDose2", &sumDose2);
0123   tree->SetBranchAddress("TotalDoseRate", &sumDoseRate);
0124   tree->SetBranchAddress("TotalDoseRate2", &sumDoseRate2);
0125
0126   Long64_t nentries = tree->GetEntries();
0127
0128   if (nentries == 0) {
0129     cout << "No entries found in the tree species contained in the file "
0130          << file->GetPath() << endl;
0131     exit(1);
0132   }
0133   tree->GetEntry(1);
0134   Double_t errDose = sqrt((sumDose2 / nEvent - pow(sumDose/nEvent, 2))
0135                           / (nEvent - 1));
0136   std::cout<<"Average dose (in Gy) : "<<sumDose/nEvent<<" +- "<<errDose<<std::endl;
0137
0138   Double_t errDoseRate = sqrt((sumDoseRate2 / nEvent - pow(sumDoseRate / nEvent, 2))
0139                           / (nEvent - 1));
0140   std::cout<<"Average dose rate (in Gy/s) : "<<sumDoseRate/nEvent<<" +- "<<errDoseRate<<std::endl;
0141
0142
0143   std::map<int, std::map<double, SpeciesInfoAOS>> speciesTimeInfo;
0144
0145   for (Int_t j = 0; j < nentries; j++) {
0146     tree->GetEntry(j);
0147
0148     SpeciesInfoAOS &infoAOS = speciesTimeInfo[speciesID][timeA];
0149
0150     infoAOS.fNumber += number;
0151     infoAOS.fG += sumG;
0152     infoAOS.fG2 += sumG2;
0153     infoAOS.fNEvent += nEvent;
0154     infoAOS.fName = speciesName;
0155   }
0156
0157   std::map <Int_t, SpeciesInfoSOA> speciesInfo;
0158
0159   auto it_SOA = speciesTimeInfo.begin();
0160   auto end_SOA = speciesTimeInfo.end();
0161
0162
0163   for (; it_SOA != end_SOA; ++it_SOA) {
0164     const Int_t _speciesID = it_SOA->first;
0165     SpeciesInfoSOA &info = speciesInfo[_speciesID];
0166
0167     auto it2 = it_SOA->second.begin();
0168     auto end2 = it_SOA->second.end();
0169
0170     info.fName = it2->second.fName;
0171     const size_t size2 = it_SOA->second.size();
0172     info.fG.resize(size2);
0173     info.fGerr.resize(size2);
0174     info.fTime.resize(size2);
0175
0176     Int_t color = (2 + _speciesID) % TColor::GetNumberOfColors();
0177     if (color == 5 || color == 10 || color == 0) ++color;
0178
0179     for (int i2 = 0; it2 != end2; ++it2, ++i2) {
0180       SpeciesInfoAOS &infoAOS = it2->second;
0181
0182       Double_t _SumG2 = infoAOS.fG2;
0183       Double_t _MeanG = infoAOS.fG / infoAOS.fNEvent;
0184       Double_t _Gerr = sqrt((_SumG2 / infoAOS.fNEvent - pow(_MeanG, 2))
0185                             / (infoAOS.fNEvent - 1));
0186
0187       info.fG[i2] = _MeanG;
0188       info.fGerr[i2] = _Gerr;
0189       info.fTime[i2] = it2->first;
0190
0191       info.fRelatErr += _Gerr / (_MeanG + 1e-30);
0192       if(info.fG[i2] != 0)
0193       {
0194         std::cout<<info.fTime[i2]<<" "<<info.fG[i2]<<" "<<info.fGerr[i2]<<" "<<info.fName<<std::endl;
0195       }
0196     }
0197
0198     TGraphErrors *gSpecies = new TGraphErrors(info.fG.size(),
0199                                               info.fTime.data(),
0200                                               info.fG.data(),
0201                                               0,
0202                                               info.fGerr.data());
0203
0204     TGCompositeFrame *tf = fTab->AddTab(info.fName.c_str());
0205     TGCompositeFrame *frame = new TGCompositeFrame(tf, 60, 60,
0206                                                    kHorizontalFrame);
0207
0208     tf->AddFrame(frame, new TGLayoutHints(kLHintsExpandX | kLHintsExpandY,
0209                                           10, 10, 10, 2));
0210
0211     TRootEmbeddedCanvas *c1 = new TRootEmbeddedCanvas(info.fName.c_str(),
0212                                                       frame, 700, 500);
0213     frame->AddFrame(c1, new TGLayoutHints(kLHintsExpandX | kLHintsExpandY,
0214                                           10, 10, 10, 2));
0215     c1->GetCanvas()->SetLogx();
0216
0217     TGHorizontalFrame *hframe = new TGHorizontalFrame(tf, 200, 40);
0218
0219     TGTextButton *save = new TGTextButton(hframe, "&Save as ...",
0220                                           "Save()");
0221     hframe->AddFrame(save, new TGLayoutHints(kLHintsCenterX, 5, 5, 3, 4));
0222
0223     TGTextButton *exit = new TGTextButton(hframe, "&Exit ",
0224                                           "gApplication->Terminate()");
0225     hframe->AddFrame(exit, new TGLayoutHints(kLHintsCenterX, 5, 5, 3, 4));
0226
0227     tf->AddFrame(hframe, new TGLayoutHints(kLHintsCenterX, 2, 2, 2, 2));
0228
0229     gSpecies->SetTitle(info.fName.c_str());
0230     gSpecies->SetMarkerStyle(20 + _speciesID);
0231     gSpecies->SetMarkerColor(color);
0232     gSpecies->SetLineColor(color);
0233     gSpecies->GetXaxis()->SetTitle("Time (ns)");
0234     gSpecies->GetXaxis()->SetTitleOffset(1.1);
0235     gSpecies->GetYaxis()->SetTitle("G value (molecules/100 eV)");
0236     gSpecies->GetYaxis()->SetTitleOffset(1.2);
0237     gSpecies->Draw("ALP");
0238   }
0239
0240   main->AddFrame(fTab, new TGLayoutHints(kLHintsBottom | kLHintsExpandX |
0241                                          kLHintsExpandY, 2, 2, 5, 1));
0242
0243   main->MapSubwindows();
0244   main->Resize();   // resize to default size
0245   main->MapWindow();
0246
0247 }
0248