HcalHits/test/plot_HE.C

0001 // Commands executed in a GLOBAL scope, e.g. created hitograms aren't erased...
0002 void plot_HE(TString  inputfile="simevent_HE.root",
0003          TString outputfile="HE_histo.root",
0004          Int_t drawmode = 2,
0005              TString    reffile="../data/HE_ref.root"){
0006
0007   // Option to no-action(0)/draw(1)/save(2) (default = 0) histograms in gif.
0008   //int doDraw = 0;
0009   int doDraw = drawmode;
0010
0011   char * treename = "Events";        //The Title of Tree.
0012
0013   delete gROOT->GetListOfFiles()->FindObject(inputfile);
0014
0015   TFile * myf  = new TFile(inputfile);
0016
0017   TTree * tree = dynamic_cast<TTree*>(myf->Get("Events"));
0018   assert(tree != 0);
0019
0020   TBranch * branchLayer = tree->GetBranch("PHcalValidInfoLayer_g4SimHits_HcalInfoLayer_CaloTest.obj");
0021   assert(branchLayer != 0);
0022
0023   TBranch * branchNxN = tree->GetBranch("PHcalValidInfoNxN_g4SimHits_HcalInfoNxN_CaloTest.obj");
0024   assert(branchNxN != 0);
0025
0026   TBranch * branchJets = tree->GetBranch( "PHcalValidInfoJets_g4SimHits_HcalInfoJets_CaloTest.obj");  assert(branchJets != 0);
0027
0028   // Just number of entries (same for all branches)
0029   int  nent = branchLayer->GetEntries();
0030   cout << "Entries branchLayer : " << nent << endl;
0031   nent = branchJets->GetEntries();
0032   cout << "Entries branchJets  : " << nent << endl;
0033   nent = branchNxN->GetEntries();
0034   cout << "Entries branchNxN   : " << nent << endl;
0035
0036   // Variables from branches
0037   PHcalValidInfoJets infoJets;
0038   branchJets->SetAddress( &infoJets);
0039   PHcalValidInfoLayer infoLayer;
0040   branchLayer->SetAddress( &infoLayer);
0041   PHcalValidInfoNxN infoNxN;
0042   branchNxN->SetAddress( &infoNxN);
0043
0044   //***************************************************************************
0045   // Histo titles-labels
0046
0047   const int Nhist1     = 47, Nhist2 = 1;  // N simple and N combined histos
0048   const int Nhist1spec =  7;              // N special out of Nsimple total
0049   const int nLayersMAX = 20;
0050   const int nDepthsMAX =  5;
0051
0052   TH1F *h;                              // just a pointer
0053   TH1F *h1[Nhist1];
0054   TH1F *h1l[nLayersMAX];                // + all scint. layers separately
0055   TH1F *h1d[nDepthsMAX];                // + all depths
0056
0057   TH2F *h2[Nhist2];
0058   TH2F *h2g[5];         // +  eta-phi grid -related for all depthes
0059
0060   char *label1[Nhist1], *label2[Nhist2], *label1l[nLayersMAX ];
0061   char *label1d[nDepthsMAX], *label2g[5];
0062
0063
0064   // simple histos
0065   label1[0]  = &"rJetHits.gif";
0066   label1[1]  = &"tJetHits.gif";
0067   label1[2]  = &"eJetHits.gif";
0068   label1[3]  = &"ecalJet.gif";
0069   label1[4]  = &"hcalJet.gif";
0070   label1[5]  = &"hoJet.gif";
0071   label1[6]  = &"etotJet.gif";
0072   label1[7]  = &"detaJet.gif";
0073   label1[8]  = &"dphiJet.gif";
0074   label1[9]  = &"drJet.gif";
0075   label1[10] = &"jetE.gif";
0076   label1[11] = &"jetEta.gif";
0077   label1[12] = &"jetPhi.gif";
0078   label1[13] = &"dijetM.gif";
0079   label1[14] = &"ecalNxNr.gif";
0080   label1[15] = &"hcalNxNr.gif";
0081   label1[16] = &"hoNxNr.gif";
0082   label1[17] = &"etotNxNr.gif";
0083   label1[18] = &"ecalNxN.gif";
0084   label1[19] = &"hcalNxN.gif";
0085   label1[20] = &"hoNxN.gif";
0086   label1[21] = &"etotNxN.gif";
0087   label1[22] = &"layerHits.gif";
0088   label1[23] = &"etaHits.gif";
0089   label1[24] = &"phiHits.gif";
0090   label1[25] = &"eHits.gif";
0091   label1[26] = &"tHits.gif";
0092   label1[27] = &"idHits.gif";
0093   label1[28] = &"jitterHits.gif";
0094   label1[29] = &"eIxI.gif";
0095   label1[30] = &"tIxI.gif";
0096   label1[31] = &"eLayer.gif";
0097   label1[32] = &"eDepth.gif";
0098   label1[33] = &"eHO.gif";
0099   label1[34] = &"eHBHE.gif";
0100   label1[35] = &"elongHF.gif";
0101   label1[36] = &"eshortHF.gif";
0102   label1[37] = &"eEcalHF.gif";
0103   label1[38] = &"eHcalHF.gif";
0104
0105   // special
0106   label1[39] = &"NxN_trans_fraction.gif";
0107   label1[40] = &"tHist_50ns.gif";
0108   label1[41] = &"tHist_eweighted.gif";
0109   label1[42] = &"nHits_ECAL.gif";
0110   label1[43] = &"nHits_HCAL.gif";
0111   label1[44] = &"nHits.gif";
0112   label1[45] = &"longProf_eweighted.gif";
0113   label1[46] = &"E_hcal.gif";
0114
0115   label1l[0] = &"layer_0.gif";
0116   label1l[1] = &"layer_1.gif";
0117   label1l[2] = &"layer_2.gif";
0118   label1l[3] = &"layer_3.gif";
0119   label1l[4] = &"layer_4.gif";
0120   label1l[5] = &"layer_5.gif";
0121   label1l[6] = &"layer_6.gif";
0122   label1l[7] = &"layer_7.gif";
0123   label1l[8] = &"layer_8.gif";
0124   label1l[9] = &"layer_9.gif";
0125   label1l[10] = &"layer_10.gif";
0126   label1l[11] = &"layer_11.gif";
0127   label1l[12] = &"layer_12.gif";
0128   label1l[13] = &"layer_13.gif";
0129   label1l[14] = &"layer_14.gif";
0130   label1l[15] = &"layer_15.gif";
0131   label1l[16] = &"layer_16.gif";
0132   label1l[17] = &"layer_17.gif";
0133   label1l[18] = &"layer_18.gif";
0134   label1l[19] = &"layer_19.gif";
0135
0136   label1d[0] = &"depth_0.gif";
0137   label1d[1] = &"depth_1.gif";
0138   label1d[2] = &"depth_2.gif";
0139   label1d[3] = &"depth_3.gif";
0140   label1d[4] = &"depth_4.gif";
0141
0142   // more complicated histos and profiles
0143
0144   label2[0] = &"JetHCALvsECAL.gif";
0145
0146   label2g[0] = &"Eta-phi_grid_depth_0.gif";
0147   label2g[1] = &"Eta-phi_grid_depth_1.gif";
0148   label2g[2] = &"Eta-phi_grid_depth_2.gif";
0149   label2g[3] = &"Eta-phi_grid_depth_3.gif";
0150   label2g[4] = &"Eta-phi_grid_all_depths.gif";
0151
0152   // Some constants
0153
0154   const float fact = 117.0; // sampling factor which corresponds to those
0155                             // for layer = 0,1 in SimG4HcalValidation.cc
0156
0157   //***************************************************************************
0158   //...Book histograms
0159
0160   for (Int_t i = 0; i < Nhist1-Nhist1spec; i++) {
0161     char hname[3];
0162     sprintf(hname,"h%d",i);
0163
0164     if(i == 4 || i == 7 || i == 8 || i == 11 || i == 12 || i == 6) {
0165       if(i == 11) h1[i] = new TH1F(hname,label1[i],100,-5.,5.);
0166       if(i == 12) h1[i] = new TH1F(hname,label1[i],72,-3.1415926,3.1415926);
0167       if(i == 7 || i == 8) h1[i] = new TH1F(hname,label1[i],100,-0.1,0.1);
0168       if( i == 4)          h1[i] = new TH1F(hname,label1[i],50,0.,100.);
0169       if( i == 6)          h1[i] = new TH1F(hname,label1[i],50,0.,100.);
0170     }
0171     else {
0172       h1[i] = new TH1F(hname,label1[i],100,1.,0.);
0173     }
0174
0175   }
0176
0177   // Special : global timing < 50 ns
0178   h1[40] = new TH1F("h40",label1[40],50,0.,50.);
0179   // Special : timing in the cluster (7x7) enery-weighted
0180   h1[41] = new TH1F("h41",label1[41],30,0.,30.);
0181   // Special : number of ECAL&HCAL hits
0182   h1[42] = new TH1F("h42",label1[42],300,0.,3000.);
0183   h1[43] = new TH1F("h43",label1[43],300,0.,3000.);
0184   h1[44] = new TH1F("h44",label1[44],300,0.,3000.);
0185
0186   // Special : Longitudinal profile
0187   h1[45] = new TH1F("h45",label1[45],20,0.,20.);
0188   // Etot HCAL
0189   TH1F *h1[46] = new TH1F("h46",label1[46],50,0.,1.0);
0190
0191   for (int i = 0; i < Nhist1; i++) {
0192     if(i != 39)  h1[i]->Sumw2();
0193   }
0194
0195
0196   for (int i = 0; i < Nhist2; i++) {
0197     char hname[3];
0198     sprintf(hname,"D%d",i);
0199     h2[i] = new TH2F(hname,label2[i],100,0.,100.,100,0.,100.);
0200   }
0201   //  h[i]->Sumw2();         // to get errors properly calculated
0202
0203   // scint. layers
0204   for (int i = 0; i < nLayersMAX; i++) {
0205     char hname[4];
0206     sprintf(hname,"hl%d",i);
0207     h1l[i] = new TH1F(hname,label1l[i],40,0.,0.4);
0208   }
0209   // depths
0210   Float_t max[5] = {30000, 500, 500, 200, 200.};
0211   for (int i = 0; i < nDepthsMAX; i++) {
0212     char hname[3];
0213     sprintf(hname,"hd%d",i);
0214     h1d[i] = new TH1F(hname,label1d[i],100,0.,max[i]);
0215   }
0216
0217   // eta-phi grid (for muon samples)
0218   for (int i = 0; i < 5; i++) {
0219     char hname[3];
0220     sprintf(hname,"Dg%d",i);
0221     h2g[i] = new TH2F(hname,label2g[i],1000,-5.,5.,576,-3.1415927,3.1415927);
0222   }
0223
0224   //***************************************************************************
0225   //***************************************************************************
0226   //...Fetch the data and fill the histogram
0227
0228   // branches separately -
0229
0230   for (int i = 0; i<nent; i++) {
0231
0232     // cout << "Ev. " << i << endl;
0233
0234     // -- get entries
0235     branchLayer ->GetEntry(i);
0236     branchNxN   ->GetEntry(i);
0237     branchJets  ->GetEntry(i);
0238
0239     // -- Leading Jet
0240     int nJetHits =  infoJets.njethit();
0241     //cout << "nJetHits = " << nJetHits << endl;
0242
0243     std::vector<float> rJetHits(nJetHits);
0244     rJetHits = infoJets.jethitr();
0245     std::vector<float> tJetHits(nJetHits);
0246     tJetHits = infoJets.jethitt();
0247     std::vector<float> eJetHits(nJetHits);
0248     eJetHits = infoJets.jethite();
0249
0250     float ecalJet = infoJets.ecaljet();
0251     float hcalJet = infoJets.hcaljet();
0252     float   hoJet = infoJets.hojet();
0253     float etotJet = infoJets.etotjet();
0254
0255     float detaJet = infoJets.detajet();
0256     float dphiJet = infoJets.dphijet();
0257     float   drJet = infoJets.drjet();
0258     float  dijetM = infoJets.dijetm();
0259
0260
0261     for (int j = 0; j < nJetHits; j++) {
0262       h1[0]->Fill(rJetHits[j]);
0263       h1[1]->Fill(tJetHits[j]);
0264       h1[2]->Fill(eJetHits[j]);
0265     }
0266
0267     h1[3]->Fill(ecalJet); //
0268     h1[4]->Fill(hcalJet); //
0269
0270     h1[5]->Fill(hoJet);
0271     h1[6]->Fill(etotJet);
0272
0273     h2[0]->Fill(ecalJet,hcalJet); //
0274
0275     h1[7]->Fill(detaJet);
0276     h1[8]->Fill(dphiJet);
0277     h1[9]->Fill(drJet);
0278
0279     h1[13]->Fill(dijetM);
0280
0281     // All Jets
0282
0283     int                nJets  = infoJets.njet();
0284     std::vector<float> jetE  (nJets);
0285     jetE   = infoJets.jete();
0286     std::vector<float> jetEta(nJets);
0287     jetEta = infoJets.jeteta();
0288     std::vector<float> jetPhi(nJets);
0289     jetPhi = infoJets.jetphi();
0290
0291
0292     for (int j = 0; j < nJets; j++) {
0293       h1[10]->Fill(jetE[j]);
0294       h1[11]->Fill(jetEta[j]);
0295       h1[12]->Fill(jetPhi[j]);
0296     }
0297
0298
0299     // NxN quantities
0300     float ecalNxNr = infoNxN.ecalnxnr();
0301     float hcalNxNr = infoNxN.hcalnxnr();
0302     float   hoNxNr = infoNxN.honxnr();
0303     float etotNxNr = infoNxN.etotnxnr();
0304
0305     float ecalNxN  = infoNxN.ecalnxn();
0306     float hcalNxN  = infoNxN.hcalnxn();
0307     float   hoNxN  = infoNxN.honxn();
0308     float etotNxN  = infoNxN.etotnxn();
0309
0310     h1[14]->Fill(ecalNxNr);
0311     h1[15]->Fill(hcalNxNr);
0312     h1[16]->Fill(hoNxNr);
0313     h1[17]->Fill(etotNxNr);
0314
0315     h1[18]->Fill(ecalNxN);
0316     h1[19]->Fill(hcalNxN);
0317     h1[20]->Fill(hoNxN);
0318     h1[21]->Fill(etotNxN);
0319
0320
0321     // CaloHits from PHcalValidInfoLayer
0322
0323     int                    nHits = infoLayer.nHit();
0324     std::vector<float>    idHits(nHits);
0325     idHits = infoLayer.idHit();
0326     std::vector<float>   phiHits(nHits);
0327     phiHits = infoLayer.phiHit();
0328     std::vector<float>   etaHits(nHits);
0329     etaHits = infoLayer.etaHit();
0330     std::vector<float> layerHits(nHits);
0331     layerHits = infoLayer.layerHit();
0332     std::vector<float>     eHits(nHits);
0333     eHits = infoLayer.eHit();
0334     std::vector<float>     tHits(nHits);
0335     tHits = infoLayer.tHit();
0336
0337     int ne = 0, nh = 0;
0338     for (int j = 0; j < nHits; j++) {
0339       int layer = layerHits[j]-1;
0340       int id    = (int)idHits[j];
0341
0342       if(id >= 10) {ne++;}
0343       else {if (id < 5) nh++;}
0344
0345       //      cout << "Hit subdet = " << id  << "  lay = " << layer << endl;
0346
0347       h1[22]->Fill(Float_t(layer));
0348       h1[23]->Fill(etaHits[j]);
0349       h1[24]->Fill(phiHits[j]);
0350       h1[25]->Fill(eHits[j]);
0351       h1[26]->Fill(tHits[j]);
0352       h1[27]->Fill(idHits[j]);
0353       //      h1[28]->Fill(jitterHits[j]);   // no jitter anymore
0354
0355       h1[40]->Fill(tHits[j]);
0356       h1[41]->Fill(tHits[j],eHits[j]);
0357
0358       if(id < 6) { // HCAL only. Depth is needed, not layer !!!
0359     //if(layer == 0)               h2g[0]->Fill(etaHits[j],phiHits[j]);
0360     //if(layer == 1)               h2g[1]->Fill(etaHits[j],phiHits[j]);
0361     //if(layer == 2)               h2g[2]->Fill(etaHits[j],phiHits[j]);
0362     //if(layer == 3)               h2g[3]->Fill(etaHits[j],phiHits[j]);
0363     h2g[4]->Fill(etaHits[j],phiHits[j]);
0364       }
0365
0366
0367     }
0368
0369     h1[42]->Fill(Float_t(ne));
0370     h1[43]->Fill(Float_t(nh));
0371     h1[44]->Fill(Float_t(nHits));
0372
0373     // NxN  PHcalValidInfoNxN
0374     //    cout << " nIxI = " << nIxI << endl;
0375     int                    nIxI = infoNxN.nnxn();
0376     std::vector<float>    idIxI(nIxI);
0377     idIxI = infoNxN.idnxn();
0378     std::vector<float>     eIxI(nIxI);
0379     eIxI  = infoNxN.enxn();
0380     std::vector<float>     tIxI(nIxI);
0381     tIxI  = infoNxN.tnxn();
0382
0383     for (int j = 0; j < nIxI ; j++) {   // NB !!! j < nIxI
0384       h1[29]->Fill(eIxI[j]);
0385       h1[30]->Fill(tIxI[j]);
0386
0387       h1[39]->Fill(idIxI[j],eIxI[j]);  // transverse profile
0388
0389     }
0390
0391     // Layers and depths PHcalValidInfoLayer
0392
0393     std::vector<float> eLayer(nLayersMAX);
0394     eLayer = infoLayer.elayer();
0395     std::vector<float> eDepth(nDepthsMAX);
0396     eDepth = infoLayer.edepth();
0397
0398     float eTot = 0.;
0399
0400     for (int j = 0; j < nLayersMAX ; j++) {
0401       h1[31]->Fill(eLayer[j]);
0402       h1l[j]->Fill(eLayer[j]);
0403
0404     h1[45]->Fill((Float_t)(j),eLayer[j]);  // HCAL SimHits only
0405         eTot += eLayer[j];
0406     }
0407     for (int j = 0; j < nDepthsMAX; j++) {
0408       h1[32]->Fill(eDepth[j]);
0409       h1d[j]->Fill(eDepth[j]);
0410     }
0411
0412     h1[46]->Fill(eTot);
0413
0414     // The rest  PHcalValidInfoLayer
0415     float eHO      =  infoLayer.eho();
0416     float eHBHE    =  infoLayer.ehbhe();
0417     float elongHF  =  infoLayer.elonghf();
0418     float eshortHF =  infoLayer.eshorthf();
0419     float eEcalHF  =  infoLayer.eecalhf();
0420     float eHcalHF  =  infoLayer.ehcalhf();
0421
0422     h1[33]->Fill(eHO);
0423     h1[34]->Fill(eHBHE);
0424
0425     h1[35]->Fill(elongHF);
0426     h1[36]->Fill(eshortHF);
0427     h1[37]->Fill(eEcalHF);
0428     h1[38]->Fill(eHcalHF);
0429
0430   }
0431
0432   // cout << "After event cycle " << i << endl;
0433
0434
0435   //...Prepare the main canva
0436   TCanvas *myc = new TCanvas("myc","",800,600);
0437   gStyle->SetOptStat(1111);   // set stat         :0 - nothing
0438
0439   // Cycle for 1D distributions
0440   for (int ihist = 0; ihist < Nhist1 ; ihist++) {
0441     if(h1[ihist]->Integral() > 1.e-30 && h1[ihist]->Integral() < 1.e30 ) {
0442
0443       h1[ihist]->SetLineColor(45);
0444       h1[ihist]->SetLineWidth(2);
0445
0446       if(doDraw == 1) {
0447     h1[ihist]->Draw("h");
0448     myc->SaveAs(label1[ihist]);
0449       }
0450     }
0451   }
0452
0453   // Cycle for energy in all layers
0454   for (int ihist = 0; ihist < nLayersMAX; ihist++) {
0455     if(h1l[ihist]->Integral() > 1.e-30 && h1l[ihist]->Integral() < 1.e30 ) {
0456
0457       h1l[ihist]->SetLineColor(45);
0458       h1l[ihist]->SetLineWidth(2);
0459
0460
0461       if(doDraw == 1) {
0462     h1l[ihist]->Draw("h");
0463     myc->SaveAs(label1l[ihist]);
0464       }
0465     }
0466   }
0467
0468
0469   // Cycle for 2D distributions
0470   //  for (int ihist = 0; ihist < 1 ; ihist++) {
0471   for (int ihist = 0; ihist < Nhist2 ; ihist++) {
0472     if(h2[ihist]->Integral() > 1.e-30 && h2[ihist]->Integral() < 1.e30 ) {
0473
0474       h2[ihist]->SetMarkerColor(45);
0475       h2[ihist]->SetMarkerStyle(20);
0476       h2[ihist]->SetMarkerSize(0.7);  // marker size !
0477
0478       h2[ihist]->SetLineColor(45);
0479       h2[ihist]->SetLineWidth(2);
0480
0481       if(doDraw == 1) {
0482     h2[ihist]->Draw();
0483     myc->SaveAs(label2[ihist]);
0484       }
0485     }
0486   }
0487
0488
0489   // Cycle for eta-phi grids
0490   //  for (int ihist = 0; ihist < 5 ; ihist++) {
0491   for (int ihist = 4; ihist < 5 ; ihist++) {
0492     if(h2g[ihist]->Integral() > 1.e-30 && h2g[ihist]->Integral() < 1.e30 ) {
0493
0494       h2g[ihist]->SetMarkerColor(41);
0495       h2g[ihist]->SetMarkerStyle(20);
0496       h2g[ihist]->SetMarkerSize(0.2);
0497
0498       h2g[ihist]->SetLineColor(41);
0499       h2g[ihist]->SetLineWidth(2);
0500
0501       if(doDraw == 1) {
0502     h2g[ihist]->Draw();
0503     myc->SaveAs(label2g[ihist]);
0504       }
0505     }
0506   }
0507
0508
0509
0510
0511   // added by Julia Yarba
0512   //-----------------------
0513   // this is a temporary stuff that I've made
0514   // to create a reference ROOT histogram file
0515
0516
0517   if (doDraw == 2) {
0518     TFile OutFile(outputfile,"RECREATE") ;
0519     int ih = 0 ;
0520     for ( ih=0; ih<nLayersMAX; ih++ )
0521       {
0522     h1l[ih]->Write() ;
0523       }
0524     for ( ih=0; ih<Nhist1; ih++ )
0525       {
0526     h1[ih]->Write() ;
0527       }
0528
0529     OutFile.Write() ;
0530     OutFile.Close() ;
0531     cout << outputfile << " histogram file created" << endl ;
0532
0533     return ;
0534
0535   }
0536   /*
0537   return;
0538   */
0539
0540
0541    // now perform Chi2 test for histograms that hold
0542    // energy deposition in the Hcal layers 1-10, using
0543    // "reference" and "current" histograms
0544
0545
0546    // open up ref. ROOT file
0547    //
0548    TFile RefFile(reffile) ;
0549
0550    // service variables
0551    //
0552    TH1F* ref_hist = 0 ;
0553    int ih = 0 ;
0554
0555    // loop over layers 1-10
0556    //
0557    for ( ih=1; ih<11; ih++ )
0558    {
0559       // service - name of the ref histo
0560       //
0561       char ref_hname[4] ;
0562       sprintf( ref_hname, "hl%d", ih ) ;
0563
0564       // retrive ref.histos one by one
0565       //
0566       ref_hist = (TH1F*)RefFile.Get( ref_hname ) ;
0567
0568       // check if valid (no-NULL)
0569       //
0570       if ( ref_hist == NULL )
0571       {
0572          // print warning in case of trouble
0573      //
0574      cout << "No such ref. histogram" << *ref_hname << endl ;
0575       }
0576       else
0577       {
0578          // everything OK - perform Chi2 test
0579      //
0580     Double_t *res;
0581     Double_t pval = h1l[ih]->Chi2Test( ref_hist, "UU", res ) ;
0582
0583      // output Chi2 comparison results
0584      //
0585      cout << "[OVAL] : Edep in Layer " << ih << ", p-value= " << pval << endl ;
0586       }
0587    }
0588
0589
0590    // loop over specials : timing,  nhits(ECAL and HCAL)
0591    //
0592    for ( ih=40; ih<47; ih++ )
0593    {
0594       // service - name of the ref histo
0595       //
0596       char ref_hname[4] ;
0597       sprintf( ref_hname, "h%d", ih ) ;
0598
0599       // retrive ref.histos one by one
0600       //
0601       ref_hist = (TH1F*)RefFile.Get( ref_hname ) ;
0602
0603       // check if valid (no-NULL)
0604       //
0605       if ( ref_hist == NULL )
0606       {
0607          // print warning in case of trouble
0608      //
0609      cout << "No such ref. histogram" << *ref_hname << endl ;
0610       }
0611       else
0612       {
0613     // everything OK - perform Chi2 test
0614     //
0615     Double_t *res;
0616     Double_t pval = h1[ih]->Chi2Test( ref_hist, "UU", res) ;
0617
0618     // output Chi2 comparison results
0619     //
0620     cout << "[OVAL] : histo " << ih << ", p-value= " << pval << endl ;
0621       }
0622    }
0623
0624
0625    // close ref. ROOT file
0626    //
0627    RefFile.Close() ;
0628
0629   // at the end, close "current" ROOT tree file
0630   //
0631   myf->Close();
0632
0633
0634   return ;
0635 }