Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 14:32:59

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 }