Back to home page

Project CMSSW displayed by LXR

 
 

    


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

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 }