Extended JTerm Upsilon Analysis Exercises

Goal of this twiki

This Twiki page presents one of the long exercises prepared for EJTERM, Commissioning and Analysis of Early Data with CMS, at LPC, FNAL, Jan. 5-9, 2010 (see the full programme). We assume the attendees to be familiar with basic C++ and Python and to have performed the pre-EJTERM exercises: see the EJTERM front page, https://www.physics.purdue.edu/particle/ejterm/. Important note for attendees of EJTERM: please post problems/questions/comments on the EJTERM Upsilon Logbook.

Introduction

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

The differential cross section is given by

   \begin{displaymath} \frac{d\sigma (pp \to \Upsilon(nS))}{dp_T} \times Br(\Upsilon(nS) \to \mu^{+}\mu^{-}) = \frac{N^{fit}_{\Upsilon(nS)}(p_T; {\cal{A}}, \epsilon)}{\int {\cal{L}} dt \cdot \Delta p_T}   \end{displaymath}

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

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

Along with executing the various steps involved in a physics analysis, the attendee will have become familiar with:

  • determining detector acceptances
  • data based methods for measuring reconstruction efficiencies
  • mc generation and sample manipulations, reconstruction, fitting
  • using various standard CMS software tools and packages
  • discover 'realistic' features in data emulated samples, and account for their effects

Analysis Steps

Exercise 1: Setup the working area

Release & Software

Create working directory, set cmssw release:

mkdir UpsilonAna 
cd UpsilonAna
source /uscmst1/prod/sw/cms/cshrc uaf 
cmsrel CMSSW_3_1_3
cd CMSSW_3_1_3/src
cmsenv

Checkout and compile required software tags:

cmscvsroot CMSSW

cvs co -r EJTERMjan2010 -d YZheng/EJtermUpsAna UserCode/YZheng/EJtermUpsAna

addpkg MuonAnalysis/MuonAssociators                     V01-01-00
addpkg PhysicsTools/TagAndProbe                         V01-08-01
cvs co MuonAnalysis/TagAndProbe
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 < ../../YZheng/EJtermUpsAna/test/tnp-v010801-p5.patch
cd ../..
addpkg GeneratorInterface/Pythia6Interface
cd GeneratorInterface/Pythia6Interface
patch -p0 < ../../YZheng/EJtermUpsAna/test/pythia6PtGun_Rapidity.patch
cd ../..
cd PhysicsTools/TagAndProbe
patch -p0 < ../../YZheng/EJtermUpsAna/test/tripleGaussianLineShape.patch
cd ../..

scramv1 b -j4
this takes a few minutes

Consider the test directory as the default working location

cd YZheng/EJtermUpsAna/test

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 cmslpc: /uscms_data/d2/zgecse/EJTerm/skim_1.87pb-1.root

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
    • /Upsilon1S/Summer09-MC_31X_V3_7TeV-v1/GEN-SIM-RECO
    • /Upsilon1S/Summer09-MC_31X_V3_7TeV-v1/GEN-SIM-RECO
    • /ppMuMuX/Summer09-MC_31X_V3_7TeV-v1/GEN-SIM-RECO

Take Upsilon1S for example. Go to the DBS webpage, search the dataset "/Upsilon1S/Summer09-MC_31X_V3_7TeV-v1/GEN-SIM-RAW". First make a note of the total number of events in this dataset: "contains 1166738 events". Then click on "Conf.files". You can see the configuration file used for generating Upsilon1S to dimuon samples. Inspect the contents and notice the lines below:

    filterEfficiency = cms.untracked.double(0.139),
The "filterEfficiency" has to do with the filter that is used in the generation (like cuts on muon Pt and Eta). It is defined as the ratio between events passing the filter and the total number of generated events by PYTHIA. For example if the filterEfficiency is 0.01, it means of every 100 events PYTHIA generates, only 1 events is of interest for us, and 99 events are thrown away.
    comEnergy = cms.double(7000.0),
The center-of-mass energy is 7 TeV.
    crossSection = cms.untracked.double(99900.0),
The number in "crossSection" is a value in [pb] and is the product of two factors: 1) The production cross-section "sigma" with which the particle is produced. 2) The forced branching ratio, "BR", leading to the desired decay mode in which the particle should be looked for:
'MDME(1035,1)=1   ! 0.024800    mu- mu+',

The integrated luminosity is the number of events before filter divided by the production cross-section. In this case, it should be (1166738/0.139)/99900 = 84.02 /pb . Similarly, we calculated the intergrated luminosity of the other datasets. And we see the background sample "ppMuMuX" has the smallest integrated luminosity 1.87/pb. So to make things more realistic, we decide to use 1.87/pb signal samples to make the level of signal and background compatible. To do this, we just need to calcuate how many events we need for the signals. Take Upsilon 1S for example again. We already measured the integrated luminosity of Upsilon 1S dataset is 84.02/pb corresponding to 11667380 events. So for 1.87/pb, we should have 1166738*(1.87/84.02) = 25967 events. For Upsilon 2S and 3S, we just repeated this calculation.

Skim the Data Samples

Inspect the python configuration file: skimUpsilon_cfg.py

Submit CRAB jobs since the dataset is not located at FNAL.

First initiate your proxy:

voms-proxy-init -voms cms

Set up CRAB:

source /uscmst1/prod/grid/CRAB/crab.csh

Inspect the CRAB configuration file test/skimUpsilon.crab . To submit the job:

crab -create -submit -cfg skimUpsilon.crab
crab -status -c skimUpsilon
When the jobs are done:
crab -getoutput -c skimUpsilon
cd skimUpsilon/res

The produced skimmed Upsilon1S data sample, named "skimUpsilon_1.root" is used in the next step. Similarly, we can prepare the same skimmed samples for Upsilon 2S and 3S along with the background sample ppMuMuX. As it will take a while to produce them, we prepared the data-like soup pointed to above.

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. 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:

   \begin{displaymath}   {\mathcal A}\left(p_T^\Upsilon,y^\Upsilon \right) = \frac{N^{\rm reco}\left( \left. p_T^\Upsilon,y^\Upsilon \right| p_T^{\rm track}>3.5{\rm GeV},|\eta^{\rm track}|<2.1 \right)} {N^{\rm gen}\left({p'}_T^\Upsilon,y'^\Upsilon \right)}   \end{displaymath}

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. It 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)
    ),
...
The ID of Upsilon(1S) is 553. We are planning to measure the Upsilon cross section in the Rapidity range (-2,2) and pT range (0,20) GeV. Adjust the kinematic phase-space of the Upsilon to somewhat larger range as due to resolution effects an Upsilon generated outside the measured range can be reconstructed within. Please ensure that the Upsilons decay only to muons by changing the branching fraction as shown below:
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'),
The schedule clearly shows the executed steps: generation, simulation,...,reconstruction:
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)
)
Test the configuration by running
cmsRun upsilonGun_cfg.py
To check the output, open the resulting upsilonGun.root file:
root -l upsilonGun.root
root [1] TBrowser b
Click on "ROOT Files", and then "upsilonGun.root" and then "UpsTree", you will find 15 branches inside. Click on each branch, to see the distribution of each parameter. Or you can scan the contents of the tree:
root [2] UpsTree->Scan("*");
It will display all the branches in the tree like below:

scanout.png

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 submit a CRAB job to produce larger statistics. Modify the upsilonGun_cfg.crab configuration to submit 100 jobs producing 500 events each. Then to submit execute:

crab -create -submit -cfg upsilonGun.crab
crab -status -c upsilonGun
Each event takes about 2 seconds, so the jobs will run for about 20 minutes. In the meantime one could start working on Step 2. When all jobs are "Done", retrieve the output with
crab -getoutput -c upsilonGun
The 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 "YZheng/EJtermUpsAna/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] );
  ...
}
There is another 2D histogram, which is a copy of the above, that collects events in which the Upsilon is reconstructed with two tracks of opposite charge, pT > 3.5 GeV, |η| < 2.1 and invariant mass within 1. GeV of the PDG value 9.46 GeV.
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 );
    }
  }
}
The ratio of the two histograms is the acceptance by definition.
acc->Divide(recoPtRap, genPtRap, 1, 1, "B");
To run the macro:
make acceptance #compiles acceptance.C
./acceptance upsilonGun.root acceptance.root
To print out the result:
root -l acceptance.root
root [1] acceptance->Print("range")
acceptancePrint.png

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 here: /uscms_data/d2/zgecse/EJTerm/upsilonGun.root. Now the binning of 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 /uscms_data/d2/zgecse/EJTerm/upsilonGun.root acceptance.root
The result can be plotted with acceptancePlot.C that saves an acceptance.png figure.
root -b -q acceptancePlot.C
acceptance.png

Quiz: 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.

Efficiency components:

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

Step 1 - 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/

Uncomment the last three lines in this file to make it as a CMSSW Plugin now:

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

Then compile:

cd ../../..
scramv1 b
cd -

We prepared this configuration file "tnp_from_skim.py" to produce the Tag and Probe Ntuple.

Modify it by changing the source root file name /uscms_data/d2/zgecse/EJTerm/skim_1.87pb-1.root

cmsRun tnp_from_skim.py
it will take ~10 minutes to finish

Two root files will be produced. 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("");

tnptree.png

Step 2 - 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
This part will take ~15min to finish

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):

fitplot_3signal.jpg

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:

IdEffPtEta.png

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

Studies on Trigger Bias:

Eventually, Tag and Probe method will be performed on real data to get efficiencies. In Monte Carlos data samples, we don't have any trigger bias in them. That is to say, we have both the events which pass trigger and the events which are not be able to pass trigger. And we are able to know whether an event passing trigger or not. However, in real data, the events will pass trigger first and then come to us. Therefore, we always have a subset of all the events. Then the question is, will the trigger efficiency measured from real data (with trigger bias) be consistent with the one measured from Monte Carlos data samples (without trigger bias)? Let's do this studies step by step as below.

1) Produce data samples with trigger bias We don't want to waste time to make the skimmed soup again. Instead, we will change the configuration file "tnp_from_skim.py" to force Tag and Probe only performed on those events which pass trigger.

Add the following lines in the "tnp_from_skim.py" just after the definition of "PASS_HLT":

process.triggerMuons = cms.EDFilter("PATMuonRefSelector",
    src = cms.InputTag("patMuons"),
    cut = cms.string( PASS_HLT ),
    filter = cms.bool(True),
)
This filter simply filters events by checking whether the event contains at least one pat muon which can pass the "HLT_Mu3" trigger path. If the size of the so-called "triggerMuons" is not zero, then this event must pass the "HLT_Mu3" trigger. Otherwise, it must not. Only the events passing "HLT_Mu3" will be keeped after this filter. Don't forget to add this process beyond all of the other processes.

2) Repeat Exercise 3 with the new "tnp_from_skim.py". Notice, when fitting the histograms, you only need to run the "fit_tnp_trigger_cfg.py" because only the trigger efficiencies may be affected due to the trigger bias. Also, please save the output files into different names. You can add "_bias" as an extension of each file name.

QUIZ: Are these differences between the tag and probe efficiencies measured from non-bias samples and bias samples?

If the answer to the above quiz is "yes", please continue doing the next step.

3) Open the "tnp_from_skim.py" again and check the tagMuons requirements:

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

QUIZ: Modify the above lines to require tag muons passing trigger.

Once you've done the quiz, repeat 2) with the new "tnp_from_skim.py". Save output files into new names.

QUIZ: Compare the efficiencies measured from 3) to the previous two tset of trigger efficienies. What will you conclude?

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()

upsilonYield.png

Step 2 - Apply selection

QUIZ: Modify test/upsilonYield_cfg.py to select muons with additional quality cuts:
  • number of valid silicon hits > 10
  • normalized chi2 of silicon fit < 5
  • TMOneStationTight tracker muon selection
Hint: upsilonYield_cfg.py should contain the following:
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)
)
Rerun the job to obtain the new tree.

Exercise 5: Cross Section

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

   \begin{displaymath} \sigma_{\rm{tot}} \times Br(\Upsilon \to \mu\mu) =  \frac{N_{\rm{fit}}^{\rm corr.}}{L}   \end{displaymath}

   \begin{displaymath} \frac{d\sigma}{dp_T} \times Br(\Upsilon \to \mu\mu) = \frac{N_{\rm{fit}}^{\rm corr.}(p_T)}{L \cdot \Delta p_T}   \end{displaymath}

  • L denotes the luminosity of the sample.
  • ΔpT 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")'

upsilonYieldWeighted.png

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. For executing the fits, it is recommended to point to a more recent version of root

source /uscms_data/d2/leonardo/root-v5.23/bin/thisroot.csh
The program may be executed as follows.
root -l -b -q oniafitter.C >& log &
or
root -l 
root [0]  .L oniafitter.C+
root [1]  oniafitter()
You'll need to edit the code to select the various steps you like to be run.

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 shape 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

void buildPdf(RooWorkspace& w) {
...
RooAbsPdf  *pdf   = new RooAddPdf ("pdf","total signal+background pdf", 
				    RooArgList(*sig1S,*sig2S,*sig3S,*bkgPdf),
				    RooArgList(*nsig1,*nsig2,*nsig3,*nbkgd));
   w.import(*pdf);
}
Verify the pdf components, and how model parameters are constrained and shared.

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

un-comment code:
  /// fit raw yield
  fitYield(*ws,"total");
  plotMass(figs_ + "/mass_raw.gif",*ws);

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.

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

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

Step 3 - 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.

un-comment code:
 /// total cross section
  xsectTot(*ws);

Differential, dσ/dpT

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.

un-comment code:
 /// cross section vs pT
  xsectDif(*ws,true);
note: set argument to false if fits do not need to be repeated

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

Extract the yields from the individual peaks, and repeat the procedure above.

Exercise 6: Analysis of a Realistic Data Sample

A sample with an induced detector or trigger anomaly will be provided.

The task of this exercise involves discovering the resulting effect on the would-be-data sample (hint: check sample characteristics, by inspecting observable distributions of interest), and implement the required corrections accounting for this, before repeating the cross section measurement.

Quiz: perform the cross section measurement on the presumable-data-like sample provided.

Exercise 7: Evaluate the Systematic Uncertainties

There are several sources of systematic uncertainties that contribute to the measured value of the Y(nS) cross section.

The luminosity is potentially the dominant source. This will be provided by CMS luminosity monitoring group. It is expected to be of order 10% during the initial running of the LHC. This folds in directly into the cross section uncertainty. Other effects need to be propagated from the source to the measured values, this is exemplified next.

Source 1 - Unknown Polarization

The Upsilon polarization is currently uncertain. The MC sample used to determine the acceptance was generated with zero polarization, α=0.

The polarization affects the angular distribution of the decay muons.

   \begin{displaymath} I(\theta^{*}) \sim 1 + \alpha \cos^{2} \theta^{*}   \end{displaymath}

θ* is the angle between the muon momentum in the Upsilon rest frame and the Upsilon momentum in the lab frame.

Quiz: Investigate the effect of non-zero polarization on the acceptance, by re-weighting the sample according to the above expression. Overlay the acceptance 1D curves for α=0, α=1 and α=-1. Use the obtained acceptance 2D maps to verify the effect on the final cross section result.

Hint: the acceptance.C code already contains the definition of cosTheta variable. Use it to set the weight variable.

Quiz: Repeat the above exercise, by taking the latest Upsilon polarization results by CDF. For the scope of the exercise, consider assigning the verified deviation as the systematic uncertainty.

Note: if interested, further details on this topic are found in cern.ch/cms-quarkonia/Polarization

Source 2 - Tag and Probe bias

The effect of potential biases introduced by the tag and probe method implementation should be considered.

Think: Why does Tag and Probe method have a bias?

Quiz: Compare the efficiency results obtained with the tag-n-probe method applied to a MC sample, to the MC-based determination. Estimate the associated systematic uncertainty arising from this source.

Source 3 - Statistical Uncertainties

The acceptance and efficiency calculations have associated statistical uncertainties (from finite sample sizes).

Quiz: Propagate the effect of these, as a systematic uncertainty to the cross section measurement.

Source 4 - Choice of PDF

The model adopted to describe the signal and background components may induce biases in the extracted yields. Example causes for this uncertainty: average resolutions, radiative tail, intrinsic fit biases.

Quiz: Provide a simple estimate of an uncertainty of this type by modifying the fit model, and re-fit the yield.

Note: The uncertainties obtained from weighted fits are underestimated (weights>1 effectively inflates the statistics) in the roofit version appended to the cmssw release; this has been corrected in more recent root versions. For obtaining the correct statistical error on the yield, it's needed to setup a more recent root version as follows:

source /uscms_data/d2/leonardo/root-v5.24/bin/thisroot.csh
and make the following change in the oniafitter's fitYield() method (pass the option SumW2Error , see details)
void fitYield(RooWorkspace& w, RooAbsData* data, TString name) {
...
  RooFitResult* fitres = w.pdf("pdf")->fitTo(*data, Save(), Extended(kTRUE), SumW2Error(data->isWeighted()));//,Range(mmin_,mmax_));
...
}

Discussion of the Results

Additional Materials

Links to quarkonium related CMS pages and presentations

Responsible: IanShipsey YuZheng ZoltanGecse NunoLeonardo

Topic attachments
I Attachment History Action Size Date Who Comment
PNGpng IdEffPtEta.png r1 manage 52.3 K 2010-01-05 - 20:52 YuZheng  
PNGpng acceptance.png r1 manage 7.3 K 2009-12-02 - 21:27 YuZheng acceptance
PNGpng acceptance2.png r1 manage 10.2 K 2010-01-05 - 22:16 ZoltanGecse  
PNGpng acceptancePrint.png r1 manage 21.3 K 2010-01-05 - 22:40 ZoltanGecse  
PNGpng fitplot.png r2 r1 manage 64.7 K 2009-12-04 - 06:00 YuZheng fit
JPEGjpg fitplot_3signal.jpg r1 manage 72.9 K 2010-01-05 - 20:45 YuZheng  
GIFgif mass_raw.gif r1 manage 9.0 K 2009-12-04 - 13:51 NunoLeonardo fit to the upsilon raw mass distribution
GIFgif mass_wei.gif r1 manage 9.7 K 2009-12-04 - 13:52 NunoLeonardo fit to the upsilon mass distribution (weighted)
PNGpng scan.png r1 manage 79.0 K 2009-12-01 - 22:54 YuZheng output of scan
PNGpng scanout.png r1 manage 79.0 K 2009-12-01 - 22:57 YuZheng content
Texttxt skimUps_cfg.py.txt r1 manage 20.3 K 2009-11-30 - 18:33 YuZheng skim upsilon
PDFpdf talk.pdf r1 manage 1602.0 K 2010-01-09 - 17:57 ZoltanGecse  
PNGpng tnptree.png r1 manage 26.1 K 2009-12-02 - 23:22 YuZheng tnp tree
PNGpng triggerEffPt.png r2 r1 manage 40.0 K 2010-01-05 - 20:50 YuZheng  
PNGpng upsilonYield.png r1 manage 74.8 K 2009-12-02 - 20:31 ZoltanGecse  
PNGpng upsilonYieldWeighted.png r1 manage 37.5 K 2009-12-04 - 00:13 ZoltanGecse  
PNGpng upsilonpt.png r2 r1 manage 13.5 K 2009-11-18 - 22:45 YuZheng upsilon_pt_w/wo_cuts
Edit | Attach | Watch | Print version | History: r94 < r93 < r92 < r91 < r90 | Backlinks | Raw View | Raw edit | More topic actions...
Topic revision: r93 - 2010-08-18 - NunoLeonardo
 
    • Cern Search Icon Cern Search
    • TWiki Search Icon TWiki Search
    • Google Search Icon Google Search

    CMSPublic All webs login

This site is powered by the TWiki collaboration platform Powered by PerlCopyright & 2008-2020 by the contributing authors. All material on this collaboration platform is the property of the contributing authors.
Ideas, requests, problems regarding TWiki? Send feedback