GlobalHits/test/MakeValidation.C

0001 ///////////////////////////////////////////////////////////////////////////////
0002 // Macro to compare histograms from the GlobalHitsProducer for validation
0003 //
0004 // root -b -q MakeValidation.C
0005 ///////////////////////////////////////////////////////////////////////////////
0006 #include "TFile.h"
0007 #include "TTree.h"
0008 #include "TBranch.h"
0009 #include "TString.h"
0010
0011 void MakeValidation(TString sfilename = "GlobalHitsHistograms.root",
0012             TString rfilename = "GlobalHitsHistograms-reference.root",
0013             TString filename = "GlobalHitsHistogramsCompare")
0014 {
0015   gROOT->Reset();
0016
0017   // zoom in on x-axis of plots
0018   bool zoomx = false;
0019
0020   // create new histograms for comparison
0021   //gROOT->LoadMacro("MakeHistograms.C");
0022   //MakeHistograms();
0023
0024   // setup names
0025   //TString sfilename = "GlobalHitsHistograms.root";
0026   //TString rfilename = "GlobalHitsHistograms-reference.root";
0027
0028   // clear memory of file names
0029   delete gROOT->GetListOfFiles()->FindObject(sfilename);
0030   delete gROOT->GetListOfFiles()->FindObject(rfilename);
0031
0032   // open reference file
0033   TFile *sfile = new TFile(sfilename);
0034   TFile *rfile = new TFile(rfilename);
0035
0036   // create canvas
0037   Int_t cWidth = 928, cHeight = 1218;
0038   //TCanvas *myCanvas = new TCanvas("globalhits","globalhits",cWidth,cHeight);
0039   TCanvas *myCanvas = new TCanvas("globalhits","globalhits");
0040   //myCanvas->Size(21.59, 27.94);
0041
0042   // open output ps file
0043   //TString filename = "GlobalHitsHistogramsCompare";
0044   TString psfile = filename+".ps";
0045   TString psfileopen = filename+".ps[";
0046   TString psfileclose = filename+".ps]";
0047   myCanvas->Print(psfileopen);
0048
0049   // create label
0050   TLatex *label = new TLatex();
0051   label->SetNDC();
0052   label->SetTextSize(0.03);
0053   label->SetTextAlign(22);
0054   TString labeltitle;
0055
0056   // ceate text
0057   TText* te = new TText();
0058   te->SetTextSize(0.075);
0059
0060   // create attributes
0061   Int_t rcolor = kBlue;
0062   Int_t rtype = kDotted;
0063   Int_t stype = kSolid;
0064   Int_t scolor = kRed;
0065   Int_t linewidth = 2;
0066
0067   vector<Int_t> histnames;
0068
0069   vector<string> mchistname;
0070   mchistname.push_back("hMCRGP1");
0071   mchistname.push_back("hMCRGP1");
0072   histnames.push_back(0);
0073
0074   vector<string> vtxhistname;
0075   vtxhistname.push_back("hMCG4Vtx1");
0076   vtxhistname.push_back("hMCG4Vtx2");
0077   vtxhistname.push_back("hGeantVtxX1");
0078   vtxhistname.push_back("hGeantVtxX2");
0079   vtxhistname.push_back("hGeantVtxY1");
0080   vtxhistname.push_back("hGeantVtxY2");
0081   vtxhistname.push_back("hGeantVtxZ1");
0082   vtxhistname.push_back("hGeantVtxZ2");
0083   histnames.push_back(1);
0084
0085   vector<string> trkhistname;
0086   trkhistname.push_back("hMCG4Trk1");
0087   trkhistname.push_back("hMCG4Trk2");
0088   trkhistname.push_back("hGeantTrkPt");
0089   trkhistname.push_back("hGeantTrkE");
0090   histnames.push_back(2);
0091
0092   vector<string> ecalhistname;
0093   ecalhistname.push_back("hCaloEcal1");
0094   ecalhistname.push_back("hCaloEcal2");
0095   ecalhistname.push_back("hCaloEcalE1");
0096   ecalhistname.push_back("hCaloEcalE2");
0097   ecalhistname.push_back("hCaloEcalToF1");
0098   ecalhistname.push_back("hCaloEcalToF2");
0099   ecalhistname.push_back("hCaloEcalPhi");
0100   ecalhistname.push_back("hCaloEcalEta");
0101   histnames.push_back(3);
0102
0103   vector<string> preshhistname;
0104   preshhistname.push_back("hCaloPreSh1");
0105   preshhistname.push_back("hCaloPreSh2");
0106   preshhistname.push_back("hCaloPreShE1");
0107   preshhistname.push_back("hCaloPreShE2");
0108   preshhistname.push_back("hCaloPreShToF1");
0109   preshhistname.push_back("hCaloPreShToF2");
0110   preshhistname.push_back("hCaloPreShPhi");
0111   preshhistname.push_back("hCaloPreShEta");
0112   histnames.push_back(4);
0113
0114   vector<string> hcalhistname;
0115   hcalhistname.push_back("hCaloHcal1");
0116   hcalhistname.push_back("hCaloHcal2");
0117   hcalhistname.push_back("hCaloHcalE1");
0118   hcalhistname.push_back("hCaloHcalE2");
0119   hcalhistname.push_back("hCaloHcalToF1");
0120   hcalhistname.push_back("hCaloHcalToF2");
0121   hcalhistname.push_back("hCaloHcalPhi");
0122   hcalhistname.push_back("hCaloHcalEta");
0123   histnames.push_back(5);
0124
0125   vector<string> pxlhistname;
0126   pxlhistname.push_back("hTrackerPx1");
0127   pxlhistname.push_back("hTrackerPx2");
0128   pxlhistname.push_back("hTrackerPxPhi");
0129   pxlhistname.push_back("hTrackerPxEta");
0130   pxlhistname.push_back("hTrackerPxBToF");
0131   pxlhistname.push_back("hTrackerPxBR");
0132   pxlhistname.push_back("hTrackerPxFToF");
0133   pxlhistname.push_back("hTrackerPxFZ");
0134   histnames.push_back(6);
0135
0136   vector<string> sihistname;
0137   sihistname.push_back("hTrackerSi1");
0138   sihistname.push_back("hTrackerSi2");
0139   sihistname.push_back("hTrackerSiPhi");
0140   sihistname.push_back("hTrackerSiEta");
0141   sihistname.push_back("hTrackerSiBToF");
0142   sihistname.push_back("hTrackerSiBR");
0143   sihistname.push_back("hTrackerSiFToF");
0144   sihistname.push_back("hTrackerSiFZ");
0145   histnames.push_back(7);
0146
0147   vector<string> muonhistname;
0148   muonhistname.push_back("hMuon1");
0149   muonhistname.push_back("hMuon2");
0150   muonhistname.push_back("hMuonPhi");
0151   muonhistname.push_back("hMuonEta");
0152   histnames.push_back(8);
0153
0154   vector<string> cschistname;
0155   cschistname.push_back("hMuonCscToF1");
0156   cschistname.push_back("hMuonCscToF2");
0157   cschistname.push_back("hMuonCscZ");
0158   histnames.push_back(9);
0159
0160   vector<string> dthistname;
0161   dthistname.push_back("hMuonDtToF1");
0162   dthistname.push_back("hMuonDtToF2");
0163   dthistname.push_back("hMuonDtR");
0164   histnames.push_back(10);
0165
0166   vector<string> rpchistname;
0167   rpchistname.push_back("hMuonRpcFToF1");
0168   rpchistname.push_back("hMuonRpcFToF2");
0169   rpchistname.push_back("hMuonRpcFZ");
0170   rpchistname.push_back("hMuonRpcBToF1");
0171   rpchistname.push_back("hMuonRpcBToF2");
0172   rpchistname.push_back("hMuonRpcBR");
0173   histnames.push_back(11);
0174
0175   //loop through histograms to prepare output
0176   for (Int_t i = 0; i < histnames.size(); ++i) {
0177
0178     vector<string> names;
0179
0180     // setup canvas depending on group of plots
0181     TCanvas *Canvas;
0182     if (i == 0) {
0183       names = mchistname;
0184       //Canvas = new TCanvas("MCRGP","MCRGP",cWidth,cHeight);
0185       Canvas = new TCanvas("MCRGP","MCRGP");
0186       //Canvas->Size(21.59, 27.94);
0187       Canvas->Divide(1,2);
0188       myCanvas = Canvas;
0189       myCanvas->cd(0);
0190       //label->DrawLatex(0.5,1.00,"Monte Carlo RawGenPart");
0191     }
0192     if (i == 1) {
0193       names = vtxhistname;
0194       //Canvas = new TCanvas("G4Vtx","G4Vtx",cWidth,cHeight);
0195       Canvas = new TCanvas("G4Vtx","G4Vtx");
0196       //Canvas->Size(21.59, 27.94);
0197       Canvas->Divide(2,4);
0198       myCanvas = Canvas;
0199       myCanvas->cd(0);
0200       //label->DrawLatex(0.5,1.00,"Geant4 Vertices");
0201     }
0202     if (i == 2) {
0203       names = trkhistname;
0204       //Canvas = new TCanvas("G4Trk","G4Trk",cWidth,cWidth);
0205       Canvas = new TCanvas("G4Trk","G4Trk");
0206       //Canvas->Size(21.59, 27.94);
0207       Canvas->Divide(2,2);
0208       myCanvas = Canvas;
0209       myCanvas->cd(0);
0210       //label->DrawLatex(0.5,1.00,"Geant4 Tracks");
0211     }
0212     if (i == 3) {
0213       names = ecalhistname;
0214       //Canvas = new TCanvas("ECalHits","ECalHits",cWidth,cHeight);
0215       Canvas = new TCanvas("ECalHits","ECalHits");
0216       //Canvas->Size(21.59, 27.94);
0217       Canvas->Divide(2,4);
0218       myCanvas = Canvas;
0219       myCanvas->cd(0);
0220       //label->DrawLatex(0.5,1.00,"ECal Information");
0221     }
0222     if (i == 4) {
0223       names = preshhistname;
0224       //Canvas = new TCanvas("PreShHits","PreShHits",cWidth,cHeight);
0225       Canvas = new TCanvas("PreShHits","PreShHits");
0226       //Canvas->Size(21.59, 27.94);
0227       Canvas->Divide(2,4);
0228       myCanvas = Canvas;
0229       myCanvas->cd(0);
0230       //label->DrawLatex(0.5,1.00,"PreShower Information");
0231     }
0232     if (i == 5) {
0233       names = hcalhistname;
0234       //Canvas = new TCanvas("HCalHits","HCalHits",cWidth,cHeight);
0235       Canvas = new TCanvas("HCalHits","HCalHits");
0236       //Canvas->Size(21.59, 27.94);
0237       Canvas->Divide(2,4);
0238       myCanvas = Canvas;
0239       myCanvas->cd(0);
0240       //label->DrawLatex(0.5,1.00,"HCal Information");
0241     }
0242     if (i == 6) {
0243       names = pxlhistname;
0244       //Canvas = new TCanvas("PixelHits","PixelHits",cWidth,cHeight);
0245       Canvas = new TCanvas("PixelHits","PixelHits");
0246       //Canvas->Size(21.59, 27.94);
0247       Canvas->Divide(2,4);
0248       myCanvas = Canvas;
0249       myCanvas->cd(0);
0250       //label->DrawLatex(0.5,1.00,"Pixel Information");
0251     }
0252     if (i == 7) {
0253       names = sihistname;
0254       //Canvas = new TCanvas("StripHits","StripHits",cWidth,cHeight);
0255       Canvas = new TCanvas("StripHits","StripHits");
0256       //Canvas->Size(21.59, 27.94);
0257       Canvas->Divide(2,4);
0258       myCanvas = Canvas;
0259       myCanvas->cd(0);
0260       //label->DrawLatex(0.5,1.00,"Strip Information");
0261     }
0262     if (i == 8) {
0263       names = muonhistname;
0264       //Canvas = new TCanvas("MuonHits","MuonHits",cWidth,cWidth);
0265       Canvas = new TCanvas("MuonHits","MuonHits");
0266       //Canvas->Size(21.59, 27.94);
0267       Canvas->Divide(2,2);
0268       myCanvas = Canvas;
0269       myCanvas->cd(0);
0270       //label->DrawLatex(0.5,1.00,"Muon Information");
0271     }
0272     if (i == 9) {
0273       names = cschistname;
0274       //Canvas = new TCanvas("MuonCscHits","MuonCscHits",cWidth,cWidth);
0275       Canvas = new TCanvas("MuonCscHits","MuonCscHits");
0276       //Canvas->Size(21.59, 27.94);
0277       Canvas->Divide(2,2);
0278       myCanvas = Canvas;
0279       myCanvas->cd(0);
0280       //label->DrawLatex(0.5,1.00,"Muon CSC Information");
0281     }
0282     if (i == 10) {
0283       names = dthistname;
0284       //Canvas = new TCanvas("MuonDtHits","MuonDtHits",cWidth,cWidth);
0285       Canvas = new TCanvas("MuonDtHits","MuonDtHits");
0286       //Canvas->Size(21.59, 27.94);
0287       Canvas->Divide(2,2);
0288       myCanvas = Canvas;
0289       myCanvas->cd(0);
0290       //label->DrawLatex(0.5,1.00,"Muon DT Information");
0291     }
0292     if (i == 11) {
0293       names = rpchistname;
0294       //Canvas = new TCanvas("MuonRpcHits","MuonRpcHits",cWidth,cWidth);
0295       Canvas = new TCanvas("MuonRpcHits","MuonRpcHits");
0296       //Canvas->Size(21.59, 27.94);
0297       Canvas->Divide(2,3);
0298       myCanvas = Canvas;
0299       myCanvas->cd(0);
0300       //label->DrawLatex(0.5,1.00,"Muon RPC Information");
0301     }
0302
0303     // loop through plots
0304     for (Int_t j = 0; j < names.size(); ++j) {
0305
0306       TH1F *sh;
0307       TH1F *rh;
0308
0309       // set axis info for the histograms
0310       if (i == 0) {
0311     TString hpath = "DQMData/GlobalHitsV/MCGeant/"+names[j];
0312     sh = (TH1F*)sfile->Get(hpath);
0313     rh = (TH1F*)rfile->Get(hpath);
0314       }
0315       if (i == 1) {
0316     TString hpath = "DQMData/GlobalHitsV/MCGeant/"+names[j];
0317     sh = (TH1F*)sfile->Get(hpath);
0318     rh = (TH1F*)rfile->Get(hpath);
0319       }
0320       if (i == 2) {
0321     TString hpath = "DQMData/GlobalHitsV/MCGeant/"+names[j];
0322     sh = (TH1F*)srcfile->Get(hpath);
0323     rh = (TH1F*)rfile->Get(hpath);
0324       }
0325       if (i == 3 || i == 4 || i == 5) {
0326     if (i == 3 || i == 4) {
0327       TString hpath = "DQMData/GlobalHitsV/ECals/"+names[j];
0328       sh = (TH1F*)sfile->Get(hpath);
0329       rh = (TH1F*)rfile->Get(hpath);
0330     }
0331     if (i == 5) {
0332       TString hpath = "DQMData/GlobalHitsV/HCals/"+names[j];
0333       sh = (TH1F*)sfile->Get(hpath);
0334       rh = (TH1F*)rfile->Get(hpath);
0335     }
0336       }
0337       if (i == 6) {
0338     TString hpath = "DQMData/GlobalHitsV/SiPixels/"+names[j];
0339     sh = (TH1F*)sfile->Get(hpath);
0340     rh = (TH1F*)rfile->Get(hpath);
0341       }
0342       if (i == 7) {
0343     TString hpath = "DQMData/GlobalHitsV/SiStrips/"+names[j];
0344     sh = (TH1F*)sfile->Get(hpath);
0345     rh = (TH1F*)rfile->Get(hpath);
0346       }
0347       if (i == 8) {
0348     TString hpath = "DQMData/GlobalHitsV/Muons/"+names[j];
0349     sh = (TH1F*)sfile->Get(hpath);
0350     rh = (TH1F*)rfile->Get(hpath);
0351       }
0352       if (i == 9) {
0353     TString hpath = "DQMData/GlobalHitsV/Muons/"+names[j];
0354     sh = (TH1F*)sfile->Get(hpath);
0355     rh = (TH1F*)rfile->Get(hpath);
0356       }
0357       if (i == 10) {
0358     TString hpath = "DQMData/GlobalHitsV/Muons/"+names[j];
0359     sh = (TH1F*)sfile->Get(hpath);
0360     rh = (TH1F*)rfile->Get(hpath);
0361       }
0362       if (i == 11) {
0363     TString hpath = "DQMData/GlobalHitsV/Muon/"+names[j];
0364     sh = (TH1F*)sfile->Get(hpath);
0365     rh = (TH1F*)rfile->Get(hpath);
0366       }
0367
0368       // extract plot from both files
0369       //TH1F *sh = (TH1F*)sfile->Get(names[j].c_str());
0370       sh->SetLineColor(scolor);
0371       sh->SetLineWidth(linewidth);
0372       sh->SetLineStyle(stype);
0373       Double_t smax = sh->GetMaximum();
0374       //TH1F *rh = (TH1F*)rfile->Get(names[j].c_str());
0375       rh->SetLineColor(rcolor);
0376       rh->SetLineWidth(linewidth);
0377       rh->SetLineStyle(rtype);
0378       Double_t rmax = rh->GetMaximum();
0379
0380       // find the probability value for the two plots being from the same
0381       // source
0382       //double pv = rh->Chi2Test(sh,"OU");
0383       double pv = rh->KolmogorovTest(sh);
0384       std::strstream buf;
0385       std::string value;
0386       buf << "PV=" << pv <<std::endl;
0387       buf >> value;
0388
0389       // set maximum of y axis
0390       Double_t max = smax;
0391       if (rmax > smax) max = rmax;
0392       max *= 1.10;
0393       rh->SetMaximum(max);
0394
0395       // zoom in on x axis
0396       if (zoomx) {
0397     Int_t nbins = rh->GetNbinsX();
0398     Int_t rlbin, slbin;
0399     Int_t rfbin = -1, sfbin = -1;
0400     for (Int_t k = 1; k < nbins; ++k) {
0401       Int_t rbincont = rh->GetBinContent(k);
0402       if (rbincont != 0) {
0403         rlbin = k;
0404         if (rfbin == -1) rfbin = k;
0405       }
0406       Int_t sbincont = sh->GetBinContent(k);
0407       if (sbincont != 0) {
0408         slbin = k;
0409         if (sfbin == -1) sfbin = k;
0410       }
0411     }
0412     Int_t binmin = rfbin, binmax = rlbin+1;
0413     if (sfbin < binmin) binmin = sfbin;
0414     if (slbin > binmax) binmax = slbin+1;
0415     rh->SetAxisRange(rh->GetBinLowEdge(binmin),rh->GetBinLowEdge(binmax));
0416       }
0417
0418       // make plots
0419       myCanvas->cd(j+1);
0420       //gPad->SetLogy();
0421       sh->Draw();
0422       rh->Draw("sames");
0423
0424       te->DrawTextNDC(0.15,0.8, value.c_str());
0425       std::cout << "[OVAL] " << rh->GetName()
0426         << " PV = " << pv << std::endl;
0427     }
0428     myCanvas->Print(psfile);
0429
0430   } // end loop through histnames
0431
0432   // close output ps file
0433   myCanvas->Print(psfileclose);
0434
0435   // close root files
0436   rfile->Close();
0437   sfile->Close();
0438
0439   //convert to pdf
0440   TString cmnd;
0441   cmnd = "ps2pdf "+psfile+" "+filename+".pdf";
0442   gSystem->Exec(cmnd);
0443   cmnd = "rm "+psfile;
0444   gSystem->Exec(cmnd);
0445
0446   return;
0447 }