BeamSpotProducer/test/NtupleHelper.cc

0001 #define NtupleHelper_cc
0002 #include "NtupleHelper.h"
0003 #include <TH2.h>
0004 #include <TStyle.h>
0005 #include <TCanvas.h>
0006 #include <TMath.h>
0007
0008 ClassImp(NtupleHelper);
0009
0010 void NtupleHelper::Book() {
0011   hvxy = new TH2F("vxy", "x vers. y vertex", 100, -0.01, 0.01, 100, -0.01, 0.01);
0012   hvxyz = new TH3F("vxyz", "x y z vertex", 100, -0.01, 0.01, 200, -40., 40., 100, -0.01, 0.01);
0013   hdphi = new TH2F("dphi", "d vers phi", 100, -0., 6.283, 100, -0.2, 0.2);
0014   hd = new TH2F("d", "d vers. z ", 100, -40., 40., 100, -0.01, 0.01);
0015   hvxz = new TH2F("vxz", "x vers. z vertex", 100, -40., 40., 100, -0.01, 0.01);
0016   hsigma = new TH1F("sigma", "sigma d ", 100, 0., 0.01);
0017   hsx = new TH1F("sx", "sigma vers. z", 50, -40., 40.);
0018   hpt = new TH1F("pt", "pt", 100, 0., 20.);
0019   hsx->Sumw2();
0020   hsd = new TH1F("sd", "d0 vers. z", 50, -40., 40.);
0021   hsd->Sumw2();
0022   hsw = new TH1F("sw", "sigma weight vers. z", 50, -40., 40.);
0023   hsw->Sumw2();
0024 }
0025 zData NtupleHelper::Loop(int maxEvents) {
0026   //   In a ROOT session, you can do:
0027   //      Root > .L NtupleHelper.C
0028   //      Root > NtupleHelper t
0029   //      Root > t.GetEntry(12); // Fill t data members with entry number 12
0030   //      Root > t.Show();       // Show values of entry 12
0031   //      Root > t.Show(16);     // Read and show values of entry 16
0032   //      Root > t.Loop();       // Loop on all entries
0033   //
0034
0035   //     This is the loop skeleton where:
0036   //    jentry is the global entry number in the chain
0037   //    ientry is the entry number in the current Tree
0038   //  Note that the argument to GetEntry must be:
0039   //    jentry for TChain::GetEntry
0040   //    ientry for TTree::GetEntry and TBranch::GetEntry
0041   //
0042   //       To read only selected branches, Insert statements like:
0043   // METHOD1:
0044   //    fChain->SetBranchStatus("*",0);  // disable all branches
0045   //    fChain->SetBranchStatus("branchname",1);  // activate branchname
0046   // METHOD2: replace line
0047   //    fChain->GetEntry(jentry);       //read all branches
0048   //by  b_branchname->GetEntry(ientry); //read only this branch
0049   //   TH2F *hvxy = new TH2F("vxy","x vers. y vertex",100,-0.4,0.4,100,-0.4,0.4);
0050
0051   std::cout << "  loop over entries" << std::endl;
0052   std::cout << "  maximum number of entries: " << maxEvents << std::endl;
0053   zData zvector;
0054   if (fChain == 0)
0055     return zvector;
0056   zvector.erase(zvector.begin(), zvector.end());
0057   Long64_t nentries = fChain->GetEntriesFast();
0058
0059   std::cout << " total number of entries: " << fChain->GetEntries() << std::endl;
0060
0061   int theevent = 0;
0062   for (Long64_t jentry = 0; jentry < nentries; jentry++) {
0063     Long64_t ientry = LoadTree(jentry);
0064     if (ientry < 0)
0065       break;
0066     fChain->GetEntry(jentry);
0067     //      if (sigmaD <0.05&&pt>4.0)
0068     //  {
0069
0070     //if (z0>-20. && z0<20. && pt>4 ) {
0071     //if (pt>2 ) {
0072     //if (pt>5 && sigmaD<0.02) { //it was pt>5
0073     if (true) {
0074       if (maxEvents > 0 && maxEvents == theevent) {
0075         std::cout << " reached maximum number of tracks, continue" << std::endl;
0076         break;
0077       }
0078       hvxy->Fill(x, y);
0079       hvxyz->Fill(x, z0, y);
0080       hdphi->Fill(phi, d0);
0081       hvxz->Fill(z0, x);
0082       hsd->Fill(z0, fabs(d0));
0083       hd->Fill(z0, d0);
0084       hsx->Fill(z0, fabs(x));
0085       hpt->Fill(pt, 1.0);
0086       hsw->Fill(z0, 1.0);
0087       hsigma->Fill(sigmaD);
0088
0089       // for reco ntuples:
0090       if (pt > 1.2 && TMath::Abs(eta) < 2.4 && TMath::Abs(d0) < 5 && TMath::Abs(z0) < 60 && nPixelLayerMeas >= 3 &&
0091           nTotLayerMeas >= 11 && normchi2 < 2 && quality && algo) {
0092         //TMath::Abs(d0)<0.06 ) {
0093         //(chi2/ndof)<5 && TMath::Prob(chi2, (int)ndof)>0.02 && TMath::Abs(eta)<2.2 ) {
0094         zvector.push_back(data(z0, sigmaz0, d0, sigmaD, phi, pt, 1.));
0095         theevent++;
0096       }
0097     }
0098     // if (Cut(ientry) < 0) continue;
0099     //  }
0100   }
0101   std::cout << zvector.size() << " selected tracks read in.\n";
0102   hsx->Divide(hsw);
0103   hsd->Divide(hsw);
0104   return zvector;
0105 }