HcalHits/test/plot_HB.C

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