dna/scavenger/plotG.C

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