CMS Tutorial @ PKU -- From Zero to Expert

Upsilon Analysis Exercises

ZhenHu, FangzhouZhang, Sijin Qian

Goal of this twiki

This Twiki page presents a tutorial prepared for new CMS students at Peking University. We assume the attendees to be familiar with basic C++ and ROOT.

Acknowledgments

CMS Group at Purdue University, Fermilab LPC, EJTerm

YuZheng, ZoltanGecse, NunoLeonardo, IanShipsey

Introduction

In these Exercises we will measure the Upsilon resonance cross section .

The differential cross section is given by

The determination of the l.h.s is obtained from the r.h.s measured inputs:

• acceptance A (see Exercise 2),
• efficiency ε (see Exercise 3),
• fitted yield N (see Exercise 4),
• and luminosity L.

The implementation presented is based on the analysis of dimuon events, collected with the HLT_Mu3 trigger. The basic objects are tracker muons.

Analysis Steps

Exercise 1 - Setup the working area

Release & Software

Log to a computing node (Currently we have node02~node15):

ssh node02

Create working directory, set cmssw release:

cp /md3000/scratch01/sw/cmsset_default.sh ./
source cmsset_default.sh

mkdir UpsilonAna 
cd UpsilonAna
scramv1 list | grep CMSSW
cmsrel CMSSW_3_1_3

cd CMSSW_3_1_3/src
cmsenv

Set up CVS

export CVS_RSH=/usr/bin/ssh
export CVSROOT=:pserver:anonymous@cmscvs.cern.ch:/cvs_server/repositories/CMSSW
cvs login

If you see an error message like Unknown host cmscvs.cern.ch on some nodes, please change to another node. You'll need to use CVS to check out code and configuration files to your area so that you can modify and use them. ( Notice The frequently-used CVS package for Jpsi and Upsilon analysis is the official Onia2MuMu package. But in this tutorial, we will use a much more condensed version PKUTutorial instead .)

Checkout and compile required software tags:

cvs co -d zhenhu/PKUTutorial UserCode/zhenhu/PKUTutorial

addpkg MuonAnalysis/MuonAssociators                     V01-01-00
addpkg PhysicsTools/TagAndProbe                         V01-08-01
addpkg MuonAnalysis/TagAndProbe                         V06-02-01
addpkg CondCore/PhysicsToolsPlugins                     V00-00-03-01
addpkg CondFormats/BTauObjects                          V01-03-00
addpkg CondFormats/DataRecord                           V06-07-03-01
addpkg CondFormats/PhysicsToolsObjects                  V01-02-07-01
addpkg RecoBTag/PerformanceDB                           V00-01-00
addpkg RecoBTag/Records                                 V00-01-00
addpkg RecoMuon/Records                                 V00-01-03

cd PhysicsTools/TagAndProbe
patch -p0 < ../../zhenhu/PKUTutorial/test/tnp-v010801-p5.patch

cd ../..
addpkg GeneratorInterface/Pythia6Interface           V01-04-39
cd GeneratorInterface/Pythia6Interface
patch -p0 < ../../zhenhu/PKUTutorial/test/pythia6PtGun_Rapidity.patch

cd ../../PhysicsTools/TagAndProbe
patch -p0 < ../../zhenhu/PKUTutorial/test/tripleGaussianLineShape.patch

cd ../..
scramv1 b -j4

Consider the test directory as the default working location

cd zhenhu/PKUTutorial/test

Exercise 2 - Calculation of Acceptance

The CMS detector does not cover the whole eta range therefore a fraction of Upsilons decaying to muons will not be detected. Also tracks from muons with low pT will not be identified as muons as the low pT muons cannot reach the muon system. Therefore there is a need to estimate what fraction of Upsilons can be detected assuming 100% detector efficiency. This fraction is called geometric and kinematic acceptance. In this analysis we accept muons in the pseudorapidity range |η| < 2.1. We make this choice to match the eta range of the single muon trigger. Within this range we assume full detector coverage. We also require both muons to have pT > 3.5 GeV to have sufficient muon identification and triggering efficiencies. The acceptance is a function of Upsilon pT, rapidity and polarization.

To estimate the acceptance of Upsilons we rely on Monte Carlo simulations. In contrast, the efficiencies of muon identification and triggering can be estimated from the data with the Tag and Probe method, described below in Exercise 3. However, the efficiency of a muon being reconstructed with a track in the silicon tracking system needs MC simulations and therefore it is combined with the acceptance. Thus, out final definition of acceptance is the fraction of generated Upsilons decaying to muons that are reconstructed with silicon tracks of opposite charge with pT > 3.5 GeV and |η| < 2.1 and also the invariant mass is in the mass window of the Upsilon. Formally:

Step 1 - Generate an "Upsilon Gun" Data Sample

When generating Upsilon events with PYTHIA, the pT spectrum of Upsilons is steeply falling and the statistics of events at high pT is low. In addition, the presence of the underlying event prolongs the simulation of events significantly. Hence we use an "Upsilon gun" that produces events that contain only the Upsilon, distributed flat in pT and Rapidity, and the daughter muons with zero polarization. Then all events are full-simulated and reconstructed. The upsilonGun_cfg.py job configuration in the test directory does the above mentioned steps. (For another generator examples, please see Upsilon1S, and H200ZZ4L.)

upsilonGun_cfg.py contains the parameters of the "Pythia6PtGun" module:

process.generator = cms.EDProducer("Pythia6PtGun",
    PGunParameters = cms.PSet(
        MinPhi = cms.double(-3.14159265359),
        MinPt = cms.double(0.0),
        ParticleID = cms.vint32(553),
        MaxRapidity = cms.double(2.0),
        MaxPhi = cms.double(3.14159265359),
        MinRapidity = cms.double(-2.0),
        AddAntiParticle = cms.bool(False),
        MaxPt = cms.double(20.0)
    ),
...

pythiaUpsilonDecays = cms.vstring(
    'MDME(1034,1)= 0 ! 0.025200    ! e- e+',
    'MDME(1035,1)= 1 ! 0.024800    ! mu- mu+',
    'MDME(1036,1)= 0 ! 0.026700    ! tau- tau+',
    'MDME(1037,1)= 0 ! 0.015000    ! d dbar',
    'MDME(1038,1)= 0 ! 0.045000    ! u ubar',
    'MDME(1039,1)= 0 ! 0.015000    ! s sbar',
    'MDME(1040,1)= 0 ! 0.045000    ! c cbar',
    'MDME(1041,1)= 0 ! 0.774300    ! g g g',
    'MDME(1042,1)= 0 ! 0.029000    ! gamma g'),

process.schedule = cms.Schedule(
  process.generation_step,
  process.simulation_step,
  process.digitisation_step,process.L1simulation_step,process.digi2raw_step,process.raw2digi_step,
  process.reconstruction_step)

It would require significant space to save every object in the output file. Therefore, to save space, an analyzer has been prepared DimuonTree.cc in the ../src directory that dumps only the necessary variables into a TTree. This basic code will also allow you to understand how to access information from a data sample root file. Fifteen parameters are saved into a TTree, including the generated momenta of Upsilons and muons along with the momenta of reconstructed tracks.

Adjust the number of required events to 100 events:

process.maxEvents = cms.untracked.PSet(
    input = cms.untracked.int32(100)
)

cmsRun upsilonGun_cfg.py

root -l upsilonGun.root
root [1] TBrowser b

root [2] UpsTree->Scan("*");

Notice: The generated values of the Y and muons depends on the random number generator in Pythia and so differ from the above example. However, the structure of the tree is the same.

If everything looks fine, please log on the CERN computing server (lxplus.cern.ch), or log on our T3 site at Peking University (grid08.phy.pku.edu.cn). Then you can submit a CRAB job to produce larger statistics.

Notice: You have to create the same CMSSW working area before submitting CRAB jobs.

Show-CRAB-instruction Hide-CRAB-instruction

voms-proxy-init -voms cms

source /home/cmssoft/CRAB/CRAB_2_6_6/crab.sh

crab -create -submit -cfg upsilonGun.crab

crab -status -c upsilonGun

crab -getoutput -c upsilonGun

The CRAB results will be in the upsilonGun/res/ directory. To merge them, use the standard ROOT tool

hadd -f upsilonGun.root upsilonGun/res/upsilonGun_*.root

The resulting data sample is saved in upsilonGun.root (overwriting the previous version).

Step 2 - Calculate Acceptance

Since the acceptance depends on the pT and Rapidity of the Upsilon, we calculate it in bins of pT and Rapidity. The polarization dependence will be studied as a systematic uncertainty. We are still in the directory "zhenhu/PKUTutorial/test". There is a ROOT macro named "acceptance.C". It defines a two dimensional histogram for the generated Upsilon pT and Rapidity and it is filled for each event of the upsilonGun.root sample generated in the previous exercise.

TH2F* genPtRap = new TH2F("genUps","",4,0,20,4,-2,2);
for(int i=0; i < t->GetEntries(); i++){
  genPtRap->Fill( genUpsPt[0], genUpsRapidity[0] );
  ...
}

TH2F* recoPtRap = (TH2F*)genPtRap->Clone("recoUps");
for(int i=0; i < t->GetEntries(); i++){
  if ( recoMuCharge[tr1]*recoMuCharge[tr2] == -1 && recoMuPt[tr1] >= 3.5 && recoMuPt[tr2] >= 3.5 && fabs(recoMuEta[tr1]) < 2.1 && fabs(recoMuEta[tr2]) < 2.1 ){
    ...
    if( minDeltaM < 1.0 ){
      recoPtRap->Fill( recoUpsPt, recoUpsRapidity );
    }
  }
}

acc->Divide(recoPtRap, genPtRap, 1, 1, "B");

make acceptance #compiles acceptance.C
./acceptance upsilonGun.root acceptance.root

root -l acceptance.root
root [1] acceptance->Print("range")

Show result... Hide result...

To measure the acceptance with better precision a larger sample of generated events is required. To save time, we have prepared a tree with 10M events following step 1. It is located at the SE of our T3 Site (grid09.phy.pku.edu.cn). You can transfer it to your local directory and give it a new name upsilonGun_new.root:

scp PKUTutorial@grid09.phy.pku.edu.cn:/md1000_1/user/zhenhu/Tutorial/upsilonGun.root ./upsilonGun_new.root

Notice: If you can not use the scp transfer, try to logout the node you are working on, and then log in hepfarm01 or hepfarm02.

Now the binning of the acceptance can be refined. Please adjust it to pT (20,0,20) and Rapidity (40,-2,2) in the acceptance.C and run:

make acceptance
./acceptance upsilonGun_new.root acceptance.root

root -b -q acceptancePlot.C

Show result... Hide result...

Quiz1: Please write a piece of C++ code to plot the muon pT distribution, muon eta distribution, Upsilon pT distribution, and Upsilon eta distribution at generator level.

Quiz2: Please write a piece of C++ code to plot the muon pT distribution, muon eta distribution, Upsilon pT distribution, and Upsilon eta distribution at reconstruction level.

Quiz3: In what pT range can the Upsilon cross section be measured if we accept muons only with pT > 10 GeV?

Exercise 3 - Using Tag and Probe to measure efficiencies

Tag and Probe methodology

Tag and Probe (T&P) is a data driven method for efficiency calculations. It utilizes resonances to identify probe objects. Assuming no efficiency correlation between the two decay products of the resonance it identifies one of the muons, called a Tag, using tight identification criteria and measures the efficiency of the other muon, referred to as a Probe. For detailed explanation of Tag and Probe analysis of Upsilon, please read CMS AN2009_119.

Efficiency components:

• Tracking efficiency (assuming it to be 1)
• Tracker muon ID efficiency
• Trigger efficiency

Step 1 - Skim the Data Samples

A MC data sample equivalent to 1.87/pb integrated luminosity (455492 events) has been produced to use in the exercise, and is located at grid09.phy.pku.edu.cn:/md1000_1/user/zhenhu/Tutorial/skim_1.87pb-1.root . Please scp transfer it to your working directory.

Next it's explained how this was put together, starting from officially produced samples.

• List of signal and background central monte carlo datasets on DBS
  • /Upsilon1S/Summer09-MC_31X_V3_7TeV-v1/GEN-SIM-RECO
  • /Upsilon2S/Summer09-MC_31X_V3_7TeV-v1/GEN-SIM-RECO
  • /Upsilon3S/Summer09-MC_31X_V3_7TeV-v1/GEN-SIM-RECO
  • /ppMuMuX/Summer09-MC_31X_V3_7TeV-v1/GEN-SIM-RECO

Take Upsilon1S for example.

Show Examples Hide Examples

    filterEfficiency = cms.untracked.double(0.139),

    comEnergy = cms.double(7000.0),

    crossSection = cms.untracked.double(99900.0),

'MDME(1035,1)=1   ! 0.024800    mu- mu+',

crab -create -submit -cfg skimUpsilon1S_7TeV.crab
crab -status -c skimUpsilon1S_7TeV

crab -getoutput -c skimUpsilon1S_7TeV
cd skimUpsilon1S_7TeV/res

Similarly, we can prepare the same skimmed samples for Upsilon 2S and 3S along with the background sample ppMuMuX. And then we can merge them together with the script merge_cfg.py in order to make a date-like soup. As it will take a while to produce them, we have prepared it as pointed to above.

Step 2 - Produce tag-and-probe Ntuples

We will use tracker muons in this analysis. First, we need to select the good tracker muons that satisfy the offline selection cuts (with good track quality in geometric and kinematic acceptance).

We prepared this file "MuonSelectorUpsilon.cc" to select good tracker muons.

Please move this module code:

cp ../src/MuonSelectorUpsilon.cc ../../../MuonAnalysis/TagAndProbe/plugins/

Open the file MuonSelectorUpsilon.cc

vi  ../../../MuonAnalysis/TagAndProbe/plugins/MuonSelectorUpsilon.cc

/*
/// Register this as a CMSSW Plugin
#include "FWCore/Framework/interface/MakerMacros.h"
DEFINE_FWK_MODULE(MuonSelectorUpsilon);
*/

Then compile:

cd ../../..
scramv1 b

Make sure you are still in the /test directory. We prepared this configuration file "tnp_from_skim.py" to produce the Tag and Probe Ntuple.

cmsRun tnp_from_skim.py

Two root files will be produced. (For your convenience, a backup of these two files are also stored in grid09.phy.pku.edu.cn:/md1000_1/user/zhenhu/Tutorial/). Each of them will contain a TTree, with invariant mass, pt, eta and passing/failing information for measuring tracker muon ID or trigger efficiency

• histo_output_MuFromTkUpsTkM.root is for the tracker muon identification efficiency
• histo_output_HltFromUpsTkM.root is for the trigger efficiency

For example:

root -l histo_output_MuFromTkUpsTkM.root
root [1] fitter_tree->Scan("");

Step 3 - Compute the efficiences

We supposed that you are still in /test . Start fitting:

cmsRun fit_tnp_Id_cfg.py
cmsRun fit_tnp_trigger_cfg.py

Once the jobs finish, there will be a root file containing fitting results produced for each tag and probe Ntuple:

process.fitMuFromTkUpsTkM = RunFit.clone(
    ReadFromFiles = [ 'histo_output_MuFromTkUpsTkM.root' ],
    FitFileName   =     'fit_result_MuFromTkUpsTkM.root'  ,
)

process.fitHltFromJPsiTkM = RunFit.clone(
    ReadFromFiles = [ 'histo_output_HltFromUpsTkM.root' ],
    FitFileName   =     'fit_result_HltFromUpsTkM.root'  ,
    Var2BinBoundaries   = cms.untracked.vdouble( -2.1,-1.2,0.0,1.2,2.1),
)

Each fit_result file contains fitting plots in each pt/eta bin along with efficiencies via Tag and Probe or MC truth.

For example:

root -l fit_result_MuFromTkUpsTkM.root
root [1] TBrowser b

Click on the root file name:

root [2] fit_canvas_pt_1_eta_1->Draw(" ");

You will see the mass fitting plots of "all probes", "passing probes", "failing probes" for tracker muon ID efficiency in pt (4.5, 6), eta (-1.2,0):

Show result... Hide result...

root [3] gStyle->SetPalette(1);
root [4] gStyle->SetOptStat(0);
root [5] fit_eff_pt_eta->Draw("colz");
root [6] fit_eff_pt_eta->DrawCopy("text same");
root [7] fit_canvas_pt_1_eta_1->SetLogx();
root [8] fit_eff_pt_eta->GetXaxis()->SetMoreLogLabels();

The tracker muon ID efficiency from Tag and Probe in each eta and pt bin is displayed:

Show result... Hide result...

You can aslo examine the MC truth efficiency "truth_eff_pt_eta" by yourself. The commands are similar to the above.

Show Trigger-Bias-Study Hide Trigger-Bias-Study

process.triggerMuons = cms.EDFilter("PATMuonRefSelector",
    src = cms.InputTag("patMuons"),
    cut = cms.string( PASS_HLT ),
    filter = cms.bool(True),
)

process.tagMuons = cms.EDFilter("PATMuonRefSelector",
    src = cms.InputTag("patMuons"),
    cut = cms.string("isGlobalMuon && pt > 3 && abs(eta) < 2.4 "), 
    filter = cms.bool(True),
)

Exercise 4 - Upsilon Selection and Yield Tree

The purpose of this exercise is to create a ROOT TTree with the Upsilon candidate information to be used for fitting the signal yield. The input is the skimmed data sample that contains the PAT objects, produced in the previous exercise. The output TTree has to contain the mass, pt and rapidity of the Upsilon candidates, as well as the pt and eta of the daughter muons for efficiency corrections. Technically this task consist of two parts: A CMSSW python configuration upsilonYield_cfg.py file that defines the selection of Upsilon candidates. And there a C++ module YieldTree.cc that dumps the variables of the candidates to a TTree.

Step 1 - Produce the data tree from the skim

Change the Event number to -1 and run the prepared configuration:

  cmsRun upsilonYield_cfg.py

Inspect the output root file:

root -l upsilonYield.root
root [1] upsilonYield->cd()
root [2] upsilonYield->Scan()

Show result... Hide result...

Step 2 - Apply selection

• number of valid silicon hits > 10
• normalized chi2 of silicon fit < 5
• TMOneStationTight tracker muon selection

process.goodMuonsOld = cms.EDFilter("PATMuonRefSelector",
    src = cms.InputTag("patMuons"),
    cut = cms.string("track.isNonnull && pt > 3.5 && abs(eta) < 2.1"),
)
process.goodMuons = cms.EDProducer("MuonSelectorUpsilon",
    src = cms.InputTag("goodMuonsOld"),
    selectGlobalMuons = cms.bool(False)
)

Exercise 5 - Cross Section

We have now the various ingredients in place for estimating the cross section.

• L denotes the luminosity of the sample.
• Δp_T is the width of the dimuon pT bin
• N is the corrected yield

The corrected yield is extracted from the fit to the dimuon invariant mass distribution. In the fit, each event carries as weight which corresponds to the inverse of the product of the Υ-acceptance and μ-efficiencies. The latter corrections are determined from the acceptance and efficiency maps determined previously.

Step 1 - Acceptance and efficiency corrections

In the first step, the event weights are computed, and appended to the data yield tree. The ROOT macro makeWeights.C will take the yield file, upsilonYield.root , as well as the efficiency and acceptance files as input and produce an output file upsilonYieldWeighted.root that contains the same upsilonYield tree but with the weight column included in addition. Please verify these files are available in the current working directory, or otherwise gather them there.

upsilonYield.root
acceptance.root
fit_result_MuFromTkUpsTkM.root
fit_result_HltFromUpsTm.root

To run it execute:

root -b -q 'makeWeights.C("upsilonYield.root", "acceptance.root", "fit_result_MuFromTkUpsTkM.root", "fit_result_HltFromUpsTm.root", "upsilonYieldWeighted.root")'

Show result... Hide result...

Step 2 - Fit the upsilon mass spectra

The baseline code for the cross section fitting part of the exercise is implemented in oniafitter.C. The various steps are implemented as separate functions therein, eg buildPdf, readData, fitYield .

Fitting is performed with the standard minuit/root/roofit tools. The program may be executed as follows.

root -l -b -q oniafitter.C >& log &

root -l 
root [0]  .L oniafitter.C+
root [0]  oniafitter()

The model description for the signal and background mass distributions need to be specified. A simple model is provided, which assigns a second order polynomial pdf for the background, while the signal peaks are each described by a double gaussian (along with tail polynomial to account for final state radiation).

Quiz: Inspect the provided pdf implementations

show hide

 RooAbsPdf  *pdf  = new RooAddPdf("pdf","total signal+background pdf", 
				   RooArgList(*sigPdf,*bkgPdf),RooArgList(*nsig,*nbkg));

Quiz: Run the raw fit, with no weights, over the full sample.

solution hide

  /// fit raw yield
  fitYield(*ws,"total");
  plotMass(figs_ + "/mass_raw.gif",*ws);

Show result... Hide result...

Suggested exercise: Inspect/quantify the goodness of the fit, and make adaptions to the model to possibly improve the description of the data.

Step 3 - Extract the corrected yield

The efficiency and acceptance corrections are embedded in the event weights, which were computed and added to the input root three previously. The corrected yield is extracted by repeating the fit to the weighted dataset.

Quiz: Run the weighted fit over the full sample.

solution hide

 /// fit corrected yield 
  ((RooDataSet*)ws->data("data"))->setWeightVar(*ws->var("weight"));
  fitYield(*ws,"total");
  plotMass(figs_ + "/mass_wei.gif",*ws);

Show result... Hide result...

Quiz: Compute the average weight. Compare the obtained weighted and raw yields as a cross-check.

Step 4 - Extract the cross section

Total, σ

All left to do is to retrieve the measured yields from the saved fit results, and normalize by the sample's luminosity (1.87 pb^-1).

Quiz: Compute the total cross section.

solution hide

 /// total cross section
  xsectTot(*ws);

Differential, dσ/dp_T

For extracting the differential cross section, the dataset is divided in bins of dimuon transverse momentum (this quantity is part of the input root tree), and the fitting procedure used above is applied to the individual pT-bin sub-samples.

Quiz: Compute the differential total cross section.

solution hide

 /// cross section vs pT
  xsectDif(*ws,true);

Quiz: Compute the total and differential cross section for the individual 1S, 2S, 3S resonances.

hint hide

Your comments and feedback

-- ZhenHu - 13-Jan-2010