StRoot: macros/qa/Kolmogorov.C Source File

00001 // $ID:$
00002 // Compare with Kolmogorov test set of QA Histogram into *.hist.root files
00003 //________________________________________________________________________________
00004 #if !defined(__CINT__) || defined(__MAKECINT__)
00005 #include "Riostream.h"
00006 #include <stdio.h>
00007 #include "TROOT.h"
00008 #include "TSystem.h"
00009 #include "TMath.h"
00010 #include "TH1.h"
00011 #include "TH2.h"
00012 #include "TH3.h"
00013 #include "TStyle.h"
00014 #include "TF1.h"
00015 #include "TProfile.h"
00016 #include "TTree.h"
00017 #include "TChain.h"
00018 #include "TFile.h"
00019 #include "TNtuple.h"
00020 #include "TCanvas.h"
00021 #include "TFileSet.h"
00022 #include "TDataSetIter.h"
00023 #include "TDataSet.h"
00024 #include "TDataSetIter.h"
00025 #include "TClassTable.h"
00026 //#include "DeDxTree.C"
00027 #include "TMinuit.h"
00028 #include "TSpectrum.h"
00029 #include "TString.h"
00030 #include "TLine.h"
00031 #include "TText.h"
00032 #include "TList.h"
00033 #include "TPolyMarker.h"
00034 #include "TKey.h"
00035 #include "TLegend.h"
00036 #include "StTree.h"
00037 #endif
00038 TCanvas *c1 = 0;
00039 TFile *F12[2] = {0,0};
00040 //________________________________________________________________________________
00041 void Kolmogorov(const Char_t *file1 = "", const Char_t *file2 = "", 
00042                 const Char_t *histName = "", 
00043                 Bool_t allProc = kTRUE,
00044                 Bool_t makePNG = kTRUE,
00045                 Double_t probCut = 1.e-2) {
00046   TString File1(file1);
00047   TString File2(file2);
00048   if (File1 == "" && File2 == "") {
00049     cout << "Usage: " << endl;
00050     cout << "root.exe \'Kolmogorov.C+("
00051          << "\"/star/rcf/test/dev/trs_sl302/Tue/year_2005/cucu200_minbias/rcf1216_05_200evts.hist.root\","
00052          << "\"/star/rcf/test/dev/trs_sl302/Wed/year_2005/cucu200_minbias/rcf1216_05_200evts.hist.root\")\'" << endl;
00053     return;
00054   }
00055   if ( gClassTable->GetID("StIOEvent") < 0) { 
00056     gSystem->Load("St_base");
00057     gSystem->Load("StUtilities");
00058   }
00059   TString HistName(histName);
00060   cout << "Run Kolmogorov Test for files:" << File1 << " and " << File2;
00061   if (HistName != "") cout << " for Histogram: " << HistName;
00062   else                cout << " for all histograms in file";
00063   if (allProc)        cout << " Process them non stop";
00064   else                cout << " Start Dialog after test fails";
00065   cout << endl;
00066   if (makePNG)        cout << " Create png files with result of failed comparision" << endl;
00067   cout << "Probability cut is " << probCut << endl;
00068   Int_t NF = 0;
00069   if (File1 != "")  F12[0] = TFile::Open(file1);
00070   if (File2 != "")  F12[1] = TFile::Open(file2);
00071   TSeqCollection *files = gROOT->GetListOfFiles();
00072   if (! files) return;
00073   TIter nextF(files);
00074   TFile *f = 0;
00075   TString FileN[2];
00076   while ((f = (TFile *) nextF())) {
00077     if (NF == 2) {cout << "more than " << NF << " files. Skip the rest." << endl; break;}
00078     F12[NF] = f;
00079     FileN[NF] = f->GetName();
00080     FileN[NF].ReplaceAll(".hist.root","");
00081     if (NF == 1) {FileN[NF].ReplaceAll(FileN[0].Data(),""); FileN[NF].ReplaceAll("/star/rcf/test/","");}
00082     TListIter nextkey( f->GetListOfKeys() );
00083     TKey *key = 0;
00084     while ((key = (TKey*) nextkey())) {
00085       TString Name(key->GetName());
00086       Int_t nh  = 0;
00087       if (Name.BeginsWith("histBranch")) {
00088         StIOEvent *io = (StIOEvent *) f->Get(key->GetName());
00089         TList *makerList = (TList *) io->fObj;
00090         if (makerList) {
00091           TListIter nextMaker(makerList);
00092           TObjectSet *o = 0;
00093           while ((o = (TObjectSet *) nextMaker())) {
00094             TList *histList = (TList *) o->GetObject();
00095             if (histList) {
00096               TH1 *h = 0;
00097               TListIter nextH(histList);
00098               while ((h = (TH1 *) nextH())) {
00099                 if (h->InheritsFrom("TH1")) {
00100                   h->SetDirectory(f);
00101                   nh++;
00102                 }
00103               }
00104             }
00105           }
00106         }
00107         cout << "Source file " << NF++ << ": " << f->GetName() << " with " << nh << " Histograms" <<  endl;
00108       }
00109     }
00110   }
00111   if (NF != 2) {cout << "Problem to get 2 files" << endl; return;}
00112   TListIter next(F12[0]->GetList());
00113   TObject *obj = 0;
00114   TH1 *h, *h1, *h2;
00115   if (! gROOT->IsBatch()) c1 = new TCanvas("c1","c1",800,500);
00116   Int_t nhists = 0;
00117   while ((obj = next()) && nhists <10000 ) {
00118     nhists++; cout << obj->GetName() << " hist " << nhists << endl;
00119     Double_t prob = 0;
00120     Double_t probN= 0;
00121     Int_t ifMult = 0;
00122     Bool_t makepng = makePNG;
00123     cout << "obj " << obj->GetName() << "/" << obj->GetTitle() << endl;
00124     if ( obj->IsA()->InheritsFrom( "TH1" )) {
00125       if (obj->IsA()->InheritsFrom( "StMultiH1F" )) ifMult = 1;
00126       if (ifMult) {
00127         h1 = new TH2F(* ((TH2F *) obj));
00128       } else 
00129         h1 = (TH1 *) obj;
00130       Stat_t e1 = h1->GetEntries();
00131       if (HistName != "" && HistName != TString(h1->GetName())) continue;
00132       cout << "Found histogram " << h1->GetName() << " in " << F12[0]->GetName() << " with " << e1 << " Entries" << endl;
00133       h = (TH1 *) F12[1]->Get(h1->GetName());
00134       if (! h) {cout << "Can't find histogram " << h1->GetName() << " in " << F12[1]->GetName() << endl; continue;}
00135       if (ifMult) {
00136         h2 = new TH2F(* ((TH2F *)h));
00137       } else 
00138         h2 = (TH1 *) h;
00139       Stat_t e2 = h2->GetEntries();
00140       cout << "Found histogram " << h2->GetName() << " in " << F12[1]->GetName() << " with " << e2 << " Entries" << endl;
00141       if (e1 <= 0.0 || e2 <= 0.0) {
00142         if (e1 + e2 > 0.0) 
00143           cout << "==== Incomparibale Histograms" << h1->GetName() << " e1: " << e1 << "\te2: " << e2 << endl;
00144         continue;
00145       }
00146       if ( h1->Integral() <= 0.0) { cout << "\tis empty" << endl; continue;}
00147       cout << "Found histogram " << h2->GetName() << " in " << F12[1]->GetName() << endl;
00148       
00149       h1->SetXTitle(FileN[0]);
00150       h2->SetXTitle(FileN[1]);
00151       prob = h1->KolmogorovTest(h2,"UOD");
00152       probN = h1->KolmogorovTest(h2,"UOND");
00153       TString pngName(h1->GetName());
00154       Int_t Fail = prob < probCut || probN < probCut;
00155       if (! Fail) makepng = kFALSE;
00156       if (prob < probCut) pngName += "FailS";
00157       if (probN < probCut) pngName += "FailN";
00158       h2->SetLineColor(2);
00159       h2->SetMarkerColor(2);
00160       //      TLegend *leg = new TLegend(0.7,0.6,1.1,0.8,"");
00161       TLegend leg(0.7,0.6,1.1,0.8,"");
00162       TString Title("_file0 (new)");//F12[0]->GetName());
00163       leg.AddEntry(h1,Title.Data());
00164       Title = Form("prob = %f",prob);
00165       leg.AddEntry(h1,Title.Data());
00166       Title = "_file1 (ref)"; //F12[1]->GetName();
00167       leg.AddEntry(h2,Title.Data());
00168       Title = Form("probN = %f",probN);
00169       leg.AddEntry(h2,Title.Data());
00170       if (! gROOT->IsBatch()) {
00171         c1->Clear();
00172 #if 0
00173         TString cName("");
00174         cName += h1->GetName();
00175         c1->SetTitle(cName.Data());
00176         c1->SetName(cName.Data());
00177 #endif
00178         c1->Divide(2,1);
00179         
00180         if (  h1->InheritsFrom( "TH2" ) ) {// && !obj->IsA()->InheritsFrom( "StMultiH1F" )) {
00181           c1->cd(1);
00182           h1->Draw();
00183           leg.Draw();
00184           c1->cd(2);
00185           h2->Draw();
00186         } else if (h1->InheritsFrom( "TH1" )) {
00187           c1->cd(2);
00188           h2->Draw();
00189           c1->cd(1);
00190           h1->Draw();
00191           //      h2->SetNormFactor(e1);
00192           h2->Draw("same");
00193           leg.Draw();
00194         } 
00195         c1->Update();
00196         pngName += ".png";
00197         pngName.ReplaceAll(" ","");
00198         if (makepng)
00199           c1->SaveAs(pngName);
00200         //        TVirtualX::Instance()->WritePixmap(c1->GetCanvasID(),-1,-1,pngName.Data());
00201       }
00202       if (prob < probCut || probN < probCut) {
00203         cout << "========  KolmogorovTest fails for " <<  h1->GetName() << " with  prob " 
00204              << prob << " ===================" << endl; 
00205         if (! allProc) {
00206           Int_t i;
00207           cout << "Type in a number (<0 break, >=0 go to the next histogram <";
00208           cin >> i; 
00209           if (i < 0) break;
00210         }
00211       }
00212     }
00213     if (ifMult) {delete h1; delete h2;}
00214   }
00215 }