Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 14:33:00

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 }