HcalHits/test/plot_HF.C

0001 // Commands executed in a GLOBAL scope, e.g. created hitograms aren't erased...
0002 void plot_HF(TString  inputfile="simevent_HF.root",
0003          TString outputfile="HF_histo.root",
0004          Int_t drawmode = 2,
0005              TString    reffile="../data/HF_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     = 20, Nhist2 = 1;  // N simple and N combined histos
0048   const int Nhist1spec =  5;              // N special out of Nsimple total
0049
0050   TH1F *h;                              // just a pointer
0051   TH1F *h1[Nhist1];
0052   TH2F *h2g;                    // +  eta-phi grid -related for all depthes
0053
0054   char *label1[Nhist1];
0055   char *label2g;
0056
0057
0058   // simple histos
0059   label1[0]  = &"rJetHits.gif";
0060   label1[1]  = &"tJetHits.gif";
0061   label1[2]  = &"eJetHits.gif";
0062   label1[3]  = &"ecalJet.gif";
0063   label1[4]  = &"hcalJet.gif";
0064   label1[5]  = &"etotJet.gif";
0065   label1[6]  = &"jetE.gif";
0066   label1[7]  = &"jetEta.gif";
0067   label1[8]  = &"jetPhi.gif";
0068   label1[9]  = &"etaHits.gif";
0069   label1[10] = &"phiHits.gif";
0070   label1[11] = &"eHits.gif";
0071   label1[12] = &"tHits.gif";
0072   label1[13] = &"eEcalHF.gif";
0073   label1[14] = &"eHcalHF.gif";
0074
0075   // special
0076   label1[15] = &"tHits_60ns.gif";
0077   label1[16] = &"tHits_eweighted.gif";
0078   label1[17] = &"nHits_HF.gif";
0079   label1[18] = &"elongHF.gif";
0080   label1[19] = &"eshortHF.gif";
0081
0082   label2g    = &"Eta-phi_grid_all_depths.gif";
0083
0084
0085   //***************************************************************************
0086   //...Book histograms
0087
0088   for (Int_t i = 0; i < Nhist1-Nhist1spec; i++) {
0089     char hname[3];
0090     sprintf(hname,"h%d",i);
0091
0092     if(i == 4 || i == 7 || i == 8 ) {
0093       if(i == 7) h1[i] = new TH1F(hname,label1[i],100,-5.,5.);
0094       if(i == 8) h1[i] = new TH1F(hname,label1[i],72,-3.1415926,3.1415926);
0095       if(i == 4) h1[i] = new TH1F(hname,label1[i],30,0.,150.);
0096     }
0097     else {
0098       h1[i] = new TH1F(hname,label1[i],100,1.,0.);
0099     }
0100
0101   }
0102
0103   // Special : global timing < 60 ns
0104   h1[15] = new TH1F("h15",label1[15],60,0.,60.);
0105   // Special : timing in the cluster (7x7) enery-weighted
0106   h1[16] = new TH1F("h16",label1[16],60,0.,60.);
0107   // Special : number of HCAL hits
0108   h1[17] = new TH1F("h17",label1[17],20,0.,20.);
0109   // Special : signal in long fibers
0110   h1[18] = new TH1F("h18",label1[18],50,0.,50.);
0111   // Special : signal in short fibers
0112   h1[19] = new TH1F("h19",label1[19],50,0.,50.);
0113
0114   for (int i = 0;  i < Nhist1; i++) {
0115      h1[i]->Sumw2();
0116   }
0117
0118   // eta-phi grid (for muon samples)
0119   h2g = new TH2F("Grid",label2g,1000,-5.,5.,576,-3.1415927,3.1415927);
0120
0121   //***************************************************************************
0122   //***************************************************************************
0123   //...Fetch the data and fill the histogram
0124
0125   // branches separately -
0126
0127   for (int i = 0; i<nent; i++) {
0128
0129     // cout << "Ev. " << i << endl;
0130
0131     // -- get entries
0132     branchLayer ->GetEntry(i);
0133     branchNxN   ->GetEntry(i);
0134     branchJets  ->GetEntry(i);
0135
0136     // -- Leading Jet
0137     int nJetHits =  infoJets.njethit();
0138     //cout << "nJetHits = " << nJetHits << endl;
0139
0140     std::vector<float> rJetHits(nJetHits);
0141     rJetHits = infoJets.jethitr();
0142     std::vector<float> tJetHits(nJetHits);
0143     tJetHits = infoJets.jethitt();
0144     std::vector<float> eJetHits(nJetHits);
0145     eJetHits = infoJets.jethite();
0146
0147     float ecalJet = infoJets.ecaljet();
0148     float hcalJet = infoJets.hcaljet();
0149     float   hoJet = infoJets.hojet();
0150     float etotJet = infoJets.etotjet();
0151
0152     float detaJet = infoJets.detajet();
0153     float dphiJet = infoJets.dphijet();
0154     float   drJet = infoJets.drjet();
0155     float  dijetM = infoJets.dijetm();
0156
0157
0158     for (int j = 0; j < nJetHits; j++) {
0159       h1[0]->Fill(rJetHits[j]);
0160       h1[1]->Fill(tJetHits[j]);
0161       h1[2]->Fill(eJetHits[j]);
0162     }
0163
0164     h1[3]->Fill(ecalJet); //
0165     h1[4]->Fill(hcalJet); //
0166     h1[5]->Fill(etotJet);
0167
0168
0169     // All Jets
0170
0171     int                nJets  = infoJets.njet();
0172     std::vector<float> jetE(nJets);
0173     jetE  = infoJets.jete();
0174     std::vector<float> jetEta(nJets);
0175     jetEta = infoJets.jeteta();
0176     std::vector<float> jetPhi(nJets);
0177     jetPhi = infoJets.jetphi();
0178
0179     for (int j = 0; j < nJets; j++) {
0180       h1[6]->Fill(jetE[j]);
0181       h1[7]->Fill(jetEta[j]);
0182       h1[8]->Fill(jetPhi[j]);
0183     }
0184
0185
0186     // CaloHits from PHcalValidInfoLayer
0187
0188     int                    nHits = infoLayer.nHit();
0189     std::vector<float>    idHits (nHits);
0190     idHits = infoLayer.idHit();
0191     std::vector<float>   phiHits (nHits);
0192     phiHits = infoLayer.phiHit();
0193     std::vector<float>   etaHits (nHits);
0194     etaHits = infoLayer.etaHit();
0195     std::vector<float> layerHits (nHits);
0196     layerHits = infoLayer.layerHit();
0197     std::vector<float>     eHits (nHits);
0198     eHits = infoLayer.eHit();
0199     std::vector<float>     tHits (nHits);
0200     tHits  = infoLayer.tHit();
0201
0202     int ne = 0, nh = 0;
0203     for (int j = 0; j < nHits; j++) {
0204       int layer = layerHits[j]-1;
0205       int id    = (int)idHits[j];
0206
0207       if(id >= 10) {ne++;}
0208       else {nh++;}
0209
0210       //      cout << "Hit subdet = " << id  << "  lay = " << layer << endl;
0211
0212       h1[9]->Fill(etaHits[j]);
0213       h1[10]->Fill(phiHits[j]);
0214       h1[11]->Fill(eHits[j]);
0215       h1[12]->Fill(tHits[j]);
0216
0217       h1[15]->Fill(tHits[j]);
0218       h1[16]->Fill(tHits[j],eHits[j]);
0219
0220       if(id < 6) { // HCAL only. Depth is needed, not layer !!!
0221     h2g->Fill(etaHits[j],phiHits[j]);
0222       }
0223
0224
0225     }
0226
0227     h1[17]->Fill(Float_t(nh));
0228
0229
0230     // The rest  PHcalValidInfoLayer
0231     float elongHF  =  infoLayer.elonghf();
0232     float eshortHF =  infoLayer.eshorthf();
0233     float eEcalHF  =  infoLayer.eecalhf();
0234     float eHcalHF  =  infoLayer.ehcalhf();
0235
0236     h1[13]->Fill(eEcalHF);
0237     h1[14]->Fill(eHcalHF);
0238
0239     h1[18]->Fill(elongHF);
0240     h1[19]->Fill(eshortHF);
0241
0242   }
0243
0244   // cout << "After event cycle " << i << endl;
0245
0246   //...Prepare the main canva
0247   TCanvas *myc = new TCanvas("myc","",800,600);
0248   gStyle->SetOptStat(1111);   // set stat         :0 - nothing
0249
0250
0251   // Cycle for 1D distributions
0252   for (int ihist = 0; ihist < Nhist1 ; ihist++) {
0253     if(h1[ihist]->Integral() > 1.e-30 && h1[ihist]->Integral() < 1.e30 ) {
0254
0255       h1[ihist]->SetLineColor(45);
0256       h1[ihist]->SetLineWidth(2);
0257
0258       if(doDraw == 1) {
0259     h1[ihist]->Draw("h");
0260     myc->SaveAs(label1[ihist]);
0261       }
0262     }
0263   }
0264
0265
0266   // eta-phi grid
0267   if(h2g->Integral() > 1.e-30 && h2g->Integral() < 1.e30 ) {
0268
0269     h2g->SetMarkerColor(41);
0270     h2g->SetMarkerStyle(20);
0271     h2g->SetMarkerSize(0.2);
0272     h2g->SetLineColor(41);
0273     h2g->SetLineWidth(2);
0274
0275     if(doDraw == 1) {
0276       h2g->Draw();
0277       myc->SaveAs(label2g);
0278     }
0279   }
0280
0281
0282   // added by Julia Yarba
0283   //-----------------------
0284   // this is a temporary stuff that I've made
0285   // to create a reference ROOT histogram file
0286
0287   if (doDraw == 2) {
0288     TFile OutFile(outputfile,"RECREATE") ;
0289     int ih = 0 ;
0290
0291     for ( ih=0; ih<Nhist1; ih++ )
0292       {
0293     h1[ih]->Write() ;
0294       }
0295
0296     OutFile.Write() ;
0297     OutFile.Close() ;
0298     cout << outputfile << " histogram file created" << endl ;
0299
0300     return;
0301
0302   }
0303
0304   /*
0305   return;
0306   */
0307
0308
0309    // now perform Chi2 test for histograms using
0310    // "reference" and "current" histograms
0311
0312
0313    // open up ref. ROOT file
0314    //
0315    TFile RefFile(reffile) ;
0316
0317    // service variables
0318    //
0319    TH1F* ref_hist = 0 ;
0320    int ih = 0 ;
0321
0322
0323    // loop over specials : timing,  nhits, simhits-E
0324    //
0325    for ( ih=15; ih<20; ih++ )
0326    {
0327       // service - name of the ref histo
0328       //
0329       char ref_hname[4] ;
0330       sprintf( ref_hname, "h%d", ih ) ;
0331
0332       // retrive ref.histos one by one
0333       //
0334       ref_hist = (TH1F*)RefFile.Get( ref_hname ) ;
0335
0336       // check if valid (no-NULL)
0337       //
0338       if ( ref_hist == NULL )
0339       {
0340          // print warning in case of trouble
0341      //
0342      cout << "No such ref. histogram" << *ref_hname << endl ;
0343       }
0344       else
0345       {
0346          // everything OK - perform Chi2 test
0347      //
0348     Double_t *res;
0349     Double_t pval = h1[ih]->Chi2Test( ref_hist, "UU",res ) ;
0350
0351     // output Chi2 comparison results
0352     //
0353     cout << "[OVAL] : histo " << ih << ", p-value= " << pval << endl ;
0354       }
0355    }
0356
0357
0358    // close ref. ROOT file
0359    //
0360    RefFile.Close() ;
0361
0362   // at the end, close "current" ROOT tree file
0363   //
0364   myf->Close();
0365
0366
0367   return ;
0368 }