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   Int_t speciesID, number, nEvent;
0104   char speciesName[500];
0105
0106
0107   TFile *file = new TFile;
0108   file = TFile::Open("Dose_0.root");
0109
0110   TDirectoryFile *dir = dynamic_cast<TDirectoryFile *>(file->Get("ntuple"));
0111   TTree *tree = (TTree *) dir->Get("0.010000");
0112   tree->SetBranchAddress("speciesID", &speciesID);
0113   tree->SetBranchAddress("number", &number);
0114   tree->SetBranchAddress("nEvent", &nEvent);
0115   tree->SetBranchAddress("speciesName", &speciesName);
0116   tree->SetBranchAddress("time", &timeA);
0117   tree->SetBranchAddress("sumG", &sumG);
0118   tree->SetBranchAddress("sumG2", &sumG2);
0119
0120   Long64_t nentries = tree->GetEntries();
0121
0122   if (nentries == 0) {
0123     cout << "No entries found in the tree species contained in the file "
0124          << file->GetPath() << endl;
0125     exit(1);
0126   }
0127
0128   std::map<int, std::map<double, SpeciesInfoAOS>> speciesTimeInfo;
0129
0130   for (Int_t j = 0; j < nentries; j++) {
0131     tree->GetEntry(j);
0132
0133     SpeciesInfoAOS &infoAOS = speciesTimeInfo[speciesID][timeA];
0134
0135     infoAOS.fNumber += number;
0136     infoAOS.fG += sumG;
0137     infoAOS.fG2 += sumG2;
0138     infoAOS.fNEvent += nEvent;
0139     infoAOS.fName = speciesName;
0140   }
0141
0142   std::map <Int_t, SpeciesInfoSOA> speciesInfo;
0143
0144   auto it_SOA = speciesTimeInfo.begin();
0145   auto end_SOA = speciesTimeInfo.end();
0146
0147
0148   for (; it_SOA != end_SOA; ++it_SOA) {
0149     const Int_t _speciesID = it_SOA->first;
0150     SpeciesInfoSOA &info = speciesInfo[_speciesID];
0151
0152     auto it2 = it_SOA->second.begin();
0153     auto end2 = it_SOA->second.end();
0154
0155     info.fName = it2->second.fName;
0156     const size_t size2 = it_SOA->second.size();
0157     info.fG.resize(size2);
0158     info.fGerr.resize(size2);
0159     info.fTime.resize(size2);
0160
0161     Int_t color = (2 + _speciesID) % TColor::GetNumberOfColors();
0162     if (color == 5 || color == 10 || color == 0) ++color;
0163
0164     for (int i2 = 0; it2 != end2; ++it2, ++i2) {
0165       SpeciesInfoAOS &infoAOS = it2->second;
0166
0167       Double_t _SumG2 = infoAOS.fG2;
0168       Double_t _MeanG = infoAOS.fG / infoAOS.fNEvent;
0169       Double_t _Gerr = sqrt((_SumG2 / infoAOS.fNEvent - pow(_MeanG, 2))
0170                             / (infoAOS.fNEvent - 1));
0171
0172       info.fG[i2] = _MeanG;
0173       info.fGerr[i2] = _Gerr;
0174       info.fTime[i2] = it2->first;
0175
0176       info.fRelatErr += _Gerr / (_MeanG + 1e-30);
0177       if(info.fG[i2] != 0)
0178       {
0179         std::cout<<info.fTime[i2]<<" "<<info.fG[i2]<<" "<<info.fGerr[i2]<<" "<<info.fName<<std::endl;
0180       }
0181     }
0182
0183     TGraphErrors *gSpecies = new TGraphErrors(info.fG.size(),
0184                                               info.fTime.data(),
0185                                               info.fG.data(),
0186                                               0,
0187                                               info.fGerr.data());
0188
0189     TGCompositeFrame *tf = fTab->AddTab(info.fName.c_str());
0190     TGCompositeFrame *frame = new TGCompositeFrame(tf, 60, 60,
0191                                                    kHorizontalFrame);
0192
0193     tf->AddFrame(frame, new TGLayoutHints(kLHintsExpandX | kLHintsExpandY,
0194                                           10, 10, 10, 2));
0195
0196     TRootEmbeddedCanvas *c1 = new TRootEmbeddedCanvas(info.fName.c_str(),
0197                                                       frame, 700, 500);
0198     frame->AddFrame(c1, new TGLayoutHints(kLHintsExpandX | kLHintsExpandY,
0199                                           10, 10, 10, 2));
0200     c1->GetCanvas()->SetLogx();
0201
0202     TGHorizontalFrame *hframe = new TGHorizontalFrame(tf, 200, 40);
0203
0204     TGTextButton *save = new TGTextButton(hframe, "&Save as ...",
0205                                           "Save()");
0206     hframe->AddFrame(save, new TGLayoutHints(kLHintsCenterX, 5, 5, 3, 4));
0207
0208     TGTextButton *exit = new TGTextButton(hframe, "&Exit ",
0209                                           "gApplication->Terminate()");
0210     hframe->AddFrame(exit, new TGLayoutHints(kLHintsCenterX, 5, 5, 3, 4));
0211
0212     tf->AddFrame(hframe, new TGLayoutHints(kLHintsCenterX, 2, 2, 2, 2));
0213
0214     gSpecies->SetTitle(info.fName.c_str());
0215     gSpecies->SetMarkerStyle(20 + _speciesID);
0216     gSpecies->SetMarkerColor(color);
0217     gSpecies->SetLineColor(color);
0218     gSpecies->GetXaxis()->SetTitle("Time (ns)");
0219     gSpecies->GetXaxis()->SetTitleOffset(1.1);
0220     gSpecies->GetYaxis()->SetTitle("G value (molecules/100 eV)");
0221     gSpecies->GetYaxis()->SetTitleOffset(1.2);
0222     gSpecies->Draw("ALP");
0223   }
0224
0225   main->AddFrame(fTab, new TGLayoutHints(kLHintsBottom | kLHintsExpandX |
0226                                          kLHintsExpandY, 2, 2, 5, 1));
0227
0228   main->MapSubwindows();
0229   main->Resize();   // resize to default size
0230   main->MapWindow();
0231
0232 }
0233