|
ATLAS Offline Software
|
Go to the documentation of this file.
6 #include "GaudiKernel/PhysicalConstants.h"
19 #include "RooRealVar.h"
20 #include "RooDataHist.h"
22 #include "RooBreitWigner.h"
23 #include "RooCBShape.h"
24 #include "RooFFTConvPdf.h"
25 #include "RooGlobalFunc.h"
26 #include "RooArgList.h"
27 #include "RooFitResult.h"
35 , m_triggerChainName(
"NoTrig")
90 m_varLabels[
"crtDiff"] =
"1/p_{T}(+) - 1/p_{T}(-) [GeV^{-1}]";
106 m_varRanges[
"crtDiff"] = std::make_pair(-0.15,0.15);
114 m_varRanges[
"crtDiff"] = std::make_pair(-0.03,0.03);
128 return StatusCode::SUCCESS;
157 m_stat =
new TH1F(
"statistics",
"Statistics",nCuts, 0.,
double(nCuts));
158 TString binLabels[] ={
"MuonIdTrk",
"TrkIsoPt40/Pt<.2",
"CombMuon",
"MCPgoodTrk",
"MuPtOK",
"eta<2.5",
"OppChargePair",
"InvMassOK"};
159 for(
int ilabel=0;ilabel<nCuts;ilabel++){
160 m_stat->GetXaxis()->SetBinLabel(ilabel+1,binLabels[ilabel]);
165 for (
const std::string& reg :
m_regions) {
169 m_varRanges[
"etaAll"] = std::make_pair(-1.05,1.05);
170 m_varRanges[
"etaPos"] = std::make_pair(-1.05,1.05);
171 m_varRanges[
"etaNeg"] = std::make_pair(-1.05,1.05);
172 }
else if (reg ==
"EAEA") {
176 }
else if (reg ==
"ECEC") {
177 m_varRanges[
"etaAll"] = std::make_pair(-2.5,-1.05);
178 m_varRanges[
"etaPos"] = std::make_pair(-2.5,-1.05);
179 m_varRanges[
"etaNeg"] = std::make_pair(-2.5,-1.05);
188 TString htitle =
hname +
"; Invmass[GeV/c^{2}]";
226 if (varD ==
"eta" || varD ==
"etaAll" || varD ==
"etaPos" || varD ==
"etaNeg" ){
228 }
else if (varD ==
"pt" || varD ==
"ptAll" || varD ==
"ptPos" || varD ==
"ptNeg" ){
230 }
else if (varD ==
"phi" || varD ==
"phiAll" || varD ==
"phiPos" || varD ==
"phiNeg" ){
240 return StatusCode::SUCCESS;
250 if(!muons.isValid()){
252 return StatusCode::FAILURE;
253 }
else ATH_MSG_DEBUG(
"Muon container successfully retrieved.");
259 for(
const auto*
muon : *muons ) {
262 if (!idTrk)
continue;
263 m_stat->Fill(
"MuonIdTrk",1);
264 double idTrkPt(0),ptSum(0);
267 return StatusCode::FAILURE;
270 idTrkPt = idTrk->
pt();
274 m_stat->Fill(
"TrkIsoPt40/Pt<.2",1);
275 if (!
muon->trackParticle(xAOD::Muon::CombinedTrackParticle))
continue;
276 m_stat->Fill(
"CombMuon",1);
279 m_stat->Fill(
"MCPgoodTrk",1);
281 if (idTrkPt <
m_ptCut*1000)
continue;
284 double idTrkEta = idTrk->
eta();
285 if (std::abs(idTrkEta)>2.5)
continue;
286 m_stat->Fill(
"eta<2.5",1);
292 int nMuons = goodMuons.size();
297 for (; mu1!=muEnd;++mu1){
298 const xAOD::TrackParticle *id1 = (*mu1)->trackParticle(xAOD::Muon::InnerDetectorTrackParticle);
300 for (; mu2!=muEnd; ++mu2){
303 double q1 = id1->
charge();
304 double q2 =
id2->charge();
305 if (q1*q2>0)
continue;
306 m_stat->Fill(
"OppChargePair",1);
316 if (invmass<m_minInvmass || invmass>
m_maxInvmass)
continue;
317 m_stat->Fill(
"InvMassOK",1);
321 double etaPos = idPos->
eta();
323 double etaNeg = idNeg->
eta();
326 double phiPos = idPos->
phi();
328 double phiNeg = idNeg->
phi();
338 double phiDiff = std::abs(phiPos - phiNeg);
345 std::string region =
"";
346 if ((std::abs(etaPos)<1.05 || std::abs(etaPos)==1.05) && (std::abs(etaNeg)<1.05 || std::abs(etaNeg)==1.05)) region=
"BB";
347 else if ((etaPos>1.05 && etaPos<2.5) && (etaNeg>1.05 && etaNeg<2.5)) region=
"EAEA";
348 else if ((etaPos<-1.05 && etaPos>-2.5) && (etaNeg<-1.05 && etaNeg>-2.5)) region=
"ECEC";
352 bool fillReg =
false;
356 if (!fillAll && !fillReg)
continue;
360 if (fillAll)
m_invmass[
"All"]->Fill(invmass);
361 if (fillReg)
m_invmass[region]->Fill(invmass);
365 std::map<std::string, TH2F*>* allVars =
nullptr;
372 for (
const std::pair<const std::string, TH2F*>&
p : *allVars) {
373 if (
p.first!=
"etaAll" &&
p.first!=
"phiAll" &&
p.first!=
"ptAll"){
376 }
else if (
p.first==
"etaAll"){
385 }
else if (
p.first==
"phiAll"){
394 }
else if (
p.first==
"ptAll"){
409 if (varD!=
"etaAll" && varD!=
"phiAll" && varD!=
"ptAll"){
412 }
else if (varD==
"etaAll"){
414 m_xDistr[
"All"][
"etaAll"]->Fill(etaPos);
415 m_xDistr[
"All"][
"etaAll"]->Fill(etaNeg);
418 m_xDistr[region][
"etaAll"]->Fill(etaPos);
419 m_xDistr[region][
"etaAll"]->Fill(etaNeg);
421 }
else if (varD==
"phiAll"){
423 m_xDistr[
"All"][
"phiAll"]->Fill(phiPos);
424 m_xDistr[
"All"][
"phiAll"]->Fill(phiNeg);
427 m_xDistr[region][
"phiAll"]->Fill(phiPos);
428 m_xDistr[region][
"phiAll"]->Fill(phiNeg);
430 }
else if (varD==
"ptAll"){
432 m_xDistr[
"All"][
"ptAll"]->Fill(ptPos);
433 m_xDistr[
"All"][
"ptAll"]->Fill(ptNeg);
436 m_xDistr[region][
"ptAll"]->Fill(ptPos);
437 m_xDistr[region][
"ptAll"]->Fill(ptNeg);
446 return StatusCode::SUCCESS;
455 for (
const std::string& reg :
m_regions) {
456 for (
const std::pair<const std::string, TH2F*>&
p :
m_2DinvmassVSx[reg]) {
458 std::vector<TH1F*> hout;
476 return StatusCode::SUCCESS;
482 TString
hname = hin->GetName();
484 TCanvas* ctemp =
new TCanvas(
"ctemp",
"ctemp",500,500);
486 hin->SetMarkerSize(1.2);
489 int nbins=hin->GetNbinsX();
491 std::ostringstream o; o<<
i;
492 TString projName =
hname + o.str();
493 TH1D* htemp = (TH1D*) (hin->ProjectionY(projName,
i+1,
i+1)->Clone());
494 htemp->SetTitle(projName);
496 if (htemp->GetEntries()>50){
497 double mean = 999., meanErr = 999.,
sigma = 999., sigmaErr = 999.,
chi2=0;
499 mean = htemp->GetMean();
500 sigma= htemp->GetRMS();
502 fn->SetParameters(
float(htemp->GetEntries())/10.,
mean,
sigma);
503 htemp->Fit(
"fn",
"RMLQN");
504 mean =
fn->GetParameter(1);
507 fn->SetParameters(
float(htemp->GetEntries())/10.,
mean,
sigma);
510 htemp->Fit(
"fn",
"RML");
511 ctemp->Print(psName);
512 }
else htemp->Fit(
"fn",
"RMLQN");
513 double frange = 2.4*
sigma;
514 double hrange = htemp->GetXaxis()->GetXmax()-htemp->GetXaxis()->GetXmin();
515 double ndf = frange/hrange*(htemp->GetNbinsX()) - 3;
518 mean =
fn->GetParameter(1);
519 meanErr =
fn->GetParError(1);
521 sigmaErr =
fn->GetParError(2);
526 RooRealVar
m(
"mass",
"dimuon invariant mass", 91.2, 71., 111.,
"GeV");
527 RooDataHist *
data =
nullptr;
528 data =
new RooDataHist(
"data",
"data",
m, htemp);
529 RooRealVar bwm0(
"bw_#mu",
"bw_#mu", 91.2, 85.2, 97.2) ;
530 RooRealVar bwsg(
"bw_#sigma",
"bw_#sigma", 2.4952) ;
531 RooBreitWigner bw(
"bw",
"bw",
m, bwm0, bwsg);
533 RooRealVar cbm0(
"cb_#mu",
"cb_#mu", 0 ) ;
534 RooRealVar cbsg(
"cb_#sigma",
"cb_#sigma", 3., 1., 10.) ;
535 RooRealVar cbal(
"cb_#alpha",
"cb_#alpha", 2.0) ;
536 RooRealVar cbn(
"cb_n",
"cb_n", 1., 0.05, 3.) ;
537 RooCBShape cb(
"cb",
"cb",
m, cbm0, cbsg, cbal, cbn);
540 RooFFTConvPdf bxc(
"bxc",
"BW (X) CB",
m, bw, cb) ;
541 bxc.fitTo(*
data, RooFit::PrintLevel(-1), RooFit::PrintEvalErrors(-1));
542 RooPlot* frame =
m.frame();
543 data->plotOn(frame, RooFit::MarkerSize(0.9));
545 bxc.plotOn (frame, RooFit::LineColor(kBlue));
548 ctemp->Print(psName);
550 mean = bwm0.getVal();
551 meanErr = bwm0.getError();
552 sigma = cbsg.getVal();
553 sigmaErr = cbsg.getError();
554 chi2 = frame->chiSquare();
560 hout.at(0)->SetBinContent(
i+1,
mean);
561 hout.at(0)->SetBinError(
i+1,meanErr);
562 hout.at(1)->SetBinContent(
i+1,
sigma);
563 hout.at(1)->SetBinError(
i+1,sigmaErr);
565 hout.at(0)->SetBinContent(
i+1,
mean);
566 hout.at(0)->SetBinError(
i+1,meanErr);
568 hout.at(0)->SetBinContent(
i+1,
sigma);
569 hout.at(0)->SetBinError(
i+1,sigmaErr);
577 if (hout.size()==2) hout.at(1)->Sumw2();
586 if (
sc.isFailure() ) {
605 if (nBLhits>0 || !(eBLhits)) countPass+=1;
606 if ((nhitsPIX+nPIXDS)>1) countPass+=1;
607 if ((nhitsSCT+nSCTDS)>5) countPass+=1;
608 if ((nPIXH+nSCTH)<3) countPass+=1;
611 int n = nTRTout + nTRThits;
612 if (std::abs(idTrk->
eta())<1.9){
613 if (
n>5 && nTRTout<(0.9*
n)) countPass+=1;
616 if (nTRTout<(0.9*
n)) countPass+=1;
621 ATH_MSG_WARNING(
"Trying to check trackquality but no xAOD::TrackParticle available!");
628 TLorentzVector particle1,particle2,
v;
631 v=particle1+particle2;
638 double px = id1->
p4().Px()+
id2->p4().Px();
639 double py = id1->
p4().Py()+
id2->p4().Py();
646 double px=id1->
p4().Px()+
id2->p4().Px();
647 double py = id1->
p4().Py()+
id2->p4().Py();
648 double pz = id1->
p4().Pz()+
id2->p4().Pz();
655 double px=id1->
p4().Px()+
id2->p4().Px();
656 double py = id1->
p4().Py()+
id2->p4().Py();
664 double qoverpt1=id1->
charge()/id1->
pt();
665 double qoverpt2=
id2->charge()/
id2->pt();
667 asym=std::abs(qoverpt1)-std::abs(qoverpt2);
virtual double pt() const override final
The transverse momentum ( ) of the particle.
@ numberOfPixelHoles
number of pixel layers on track with absence of hits [unit8_t].
virtual StatusCode procHistograms()
An inheriting class should either override this function or finalHists().
char data[hepevt_bytes_allocation_ATLAS]
std::map< std::string, std::pair< double, double > > m_varRanges
Const iterator class for DataVector/DataList.
virtual StatusCode fillHistograms()
An inheriting class should either override this function or fillHists().
void mean(std::vector< double > &bins, std::vector< double > &values, const std::vector< std::string > &files, const std::string &histname, const std::string &tplotname, const std::string &label="")
Scalar phi() const
phi method
double getEta(const xAOD::TrackParticle *id1, const xAOD::TrackParticle *id2) const
std::string find(const std::string &s)
return a remapped string
std::vector< std::string > m_varsVSmean
std::vector< std::string > m_varsVSwidth
@ VIEW_ELEMENTS
this data object is a view, it does not own its elmts
iterator end() noexcept
Return an iterator pointing past the end of the collection.
float charge() const
Returns the charge.
Scalar eta() const
pseudorapidity method
DataVector adapter that acts like it holds const pointers.
Gaudi::Details::PropertyBase & declareProperty(Gaudi::Property< T > &t)
virtual double eta() const override final
The pseudorapidity ( ) of the particle.
bool summaryValue(uint8_t &value, const SummaryType &information) const
Accessor for TrackSummary values.
double getPt(const xAOD::TrackParticle *id1, const xAOD::TrackParticle *id2) const
@ numberOfPixelHits
these are the pixel hits, including the b-layer [unit8_t].
@ expectInnermostPixelLayerHit
Do we expect a 0th-layer barrel hit for this track?
@ numberOfTRTHits
number of TRT hits [unit8_t].
def TH2F(name, title, nxbins, bins_par2, bins_par3, bins_par4, bins_par5=None, bins_par6=None, path='', **kwargs)
std::map< std::string, std::map< std::string, TH1F * > > m_xDistr
static Environment_t environment()
Returns the running environment of the monitoring application to help ManagedMonitorToolBase objects ...
std::map< std::string, TH1F * > m_invmass
virtual FourMom_t p4() const override final
The full 4-momentum of the particle.
static DataType_t dataType()
Returns the data type that the monitoring application is running over to help ManagedMonitorToolBase ...
std::map< std::string, std::map< std::string, TH2F * > > m_2DinvmassVSx
std::map< std::string, std::map< std::string, TH1F * > > m_invmassVSx
::StatusCode StatusCode
StatusCode definition for legacy code.
double chi2(TH1 *h0, TH1 *h1)
@ numberOfSCTHoles
number of SCT holes [unit8_t].
SG::ReadHandleKey< xAOD::MuonContainer > m_muonCollection
std::map< std::string, std::map< std::string, TH1F * > > m_widthVSx
std::string m_triggerChainName
value_type push_back(value_type pElem)
Add an element to the end of the collection.
double getPhi(const xAOD::TrackParticle *id1, const xAOD::TrackParticle *id2) const
std::map< std::string, double > m_varValues
double getCrtDiff(const xAOD::TrackParticle *id1, const xAOD::TrackParticle *id2) const
virtual StatusCode bookHistograms()
An inheriting class should either override this function or bookHists().
@ numberOfTRTOutliers
number of TRT outliers [unit8_t].
bool trackQuality(const xAOD::TrackParticle *idTrk)
std::vector< std::string > m_regions
virtual StatusCode initialize()
void iterativeGausFit(TH2F *hin, const std::vector< TH1F * > &hout, int mode)
DiMuMon(const std::string &type, const std::string &name, const IInterface *parent)
#define ATH_MSG_WARNING(x)
@ numberOfSCTDeadSensors
number of dead SCT sensors crossed [unit8_t].
DataVector adapter that acts like it holds const pointers.
void RegisterHisto(MonGroup &mon, T *histo)
@ numberOfSCTHits
number of hits in SCT [unit8_t].
@ numberOfPixelDeadSensors
number of dead pixel sensors crossed [unit8_t].
def TH1F(name, title, nxbins, bins_par2, bins_par3=None, path='', **kwargs)
std::vector< std::string > m_varsDistr
Class describing a TrackParticle.
Define macros for attributes used to control the static checker.
constexpr int pow(int base, int exp) noexcept
iterator begin() noexcept
Return an iterator pointing at the beginning of the collection.
double getInvmass(const xAOD::TrackParticle *track1, const xAOD::TrackParticle *track2, double Mass) const
virtual double phi() const override final
The azimuthal angle ( ) of the particle (has range to .)
@ numberOfInnermostPixelLayerHits
these are the hits in the 0th pixel barrel layer
std::map< std::string, std::string > m_varLabels