29 #define AliAnalysisTaskMuonCuts_cxx
42 #include "TObjString.h"
43 #include "TObjArray.h"
50 #include "AliAODEvent.h"
51 #include "AliAODTrack.h"
52 #include "AliAODMCParticle.h"
53 #include "AliMCEvent.h"
54 #include "AliMCParticle.h"
56 #include "AliESDEvent.h"
57 #include "AliESDMuonTrack.h"
58 #include "AliInputEventHandler.h"
60 #include "AliAnalysisManager.h"
63 #include "AliMergeableCollection.h"
64 #include "AliMuonEventCuts.h"
65 #include "AliMuonTrackCuts.h"
66 #include "AliAnalysisMuonUtility.h"
67 #include "AliMuonAnalysisOutput.h"
87 AliVAnalysisMuon(name, cuts),
95 TString histoTypeKeys =
"hDCAxVsP hDCAyVsP hPdcaVsP hPDCAVsPCheck hDCAVsPCheck hChiProbVsP hSigmaVsPt hSigmaVsEta";
98 TString thetaAbsKeys =
"ThetaAbs23 ThetaAbs310";
122 for (
Int_t isrc=0; isrc<kNtrackSources; ++isrc ) {
124 histo =
new TH2F(histoName.Data(), histoName.Data(), 50, 0., 200., 100, -50., 50.);
125 histo->SetXTitle(
"p (GeV/c)");
126 histo->SetYTitle(
"DCA_{x} (cm)");
127 AddObjectToCollection(histo);
130 histo =
new TH2F(histoName.Data(), histoName.Data(), 50, 0., 200., 100, -50., 50.);
131 histo->SetXTitle(
"p (GeV/c)");
132 histo->SetYTitle(
"DCA_{y} (cm)");
133 AddObjectToCollection(histo);
136 histo =
new TH2F(histoName.Data(), histoName.Data(), 100, 0., 400., 50, 0., 500.);
137 histo->SetXTitle(
"p (GeV/c)");
138 histo->SetYTitle(
"p#timesDCA (cm #times GeV/c)");
139 AddObjectToCollection(histo);
142 histo =
new TH2F(histoName.Data(), histoName.Data(), 100, 0., 800., 100, 0., 5000.);
143 histo->SetXTitle(
"p (GeV/c)");
144 histo->SetYTitle(
"p#timesDCA (cm #times GeV/c)");
145 AddObjectToCollection(histo);
148 histo =
new TH2F(histoName.Data(), histoName.Data(), 100, 0., 800., 400, 0., 200.);
149 histo->SetXTitle(
"p (GeV/c)");
150 histo->SetYTitle(
"DCA (cm)");
151 AddObjectToCollection(histo);
154 histo =
new TH2F(histoName.Data(), histoName.Data(), 50, 0., 200., 50, 0., 1.);
155 histo->SetXTitle(
"p (GeV/c)");
156 histo->SetYTitle(
"Chisquare prob.");
157 AddObjectToCollection(histo);
161 histo->SetXTitle(
"p_{t} (GeV/c)");
162 histo->SetYTitle(
"#sigma_{p#timesDCA}");
164 histo->GetYaxis()->SetBinLabel(ibin+1,Form(
"%g",
fSigmaCuts[ibin]));
166 AddObjectToCollection(histo);
170 histo->SetXTitle(
"#eta");
171 histo->SetYTitle(
"#sigma_{p#timesDCA}");
173 histo->GetYaxis()->SetBinLabel(ibin+1,Form(
"%g",
fSigmaCuts[ibin]));
175 AddObjectToCollection(histo);
179 fMuonTrackCuts->Print();
200 AliVParticle* track = 0x0;
201 Int_t nTracks = AliAnalysisMuonUtility::GetNTracks(InputEvent());
202 for (
Int_t itrack = 0; itrack < nTracks; itrack++) {
203 track = AliAnalysisMuonUtility::GetTrack(itrack, InputEvent());
204 fMuonTrackCuts->CustomParam()->SetNSigmaPdca(1.e10);
205 if ( ! fMuonTrackCuts->IsSelected(track) )
continue;
207 Double_t thetaAbsEndDeg = AliAnalysisMuonUtility::GetThetaAbsDeg(track);
210 Int_t trackSrc = GetParticleType(track);
212 TVector3 dcaAtVz = fMuonTrackCuts->GetCorrectedDCA(track);
213 Double_t pTotMean = fMuonTrackCuts->GetAverageMomentum(track);
217 Double_t sigmaPdca = fMuonTrackCuts->IsThetaAbs23(track) ? fMuonTrackCuts->GetMuonTrackCutsParam().GetSigmaPdca23() : fMuonTrackCuts->GetMuonTrackCutsParam().GetSigmaPdca310();
220 Double_t chiProb = TMath::Prob(chi2, 2);
226 for (
Int_t itrig=0; itrig<selectTrigClasses.GetEntries(); ++itrig ) {
227 TString trigClassName = ((TObjString*)selectTrigClasses.At(itrig))->GetString();
229 ((
TH2*)GetMergeableObject(physSel, trigClassName, centrality, histoName))->Fill(pTot, dcaAtVz.X());
232 ((
TH2*)GetMergeableObject(physSel, trigClassName, centrality, histoName))->Fill(pTot, dcaAtVz.Y());
235 ((
TH2*)GetMergeableObject(physSel, trigClassName, centrality, histoName))->Fill(pTot, pDca);
238 ((
TH2*)GetMergeableObject(physSel, trigClassName, centrality, histoName))->Fill(pTot, pDca);
241 ((
TH2*)GetMergeableObject(physSel, trigClassName, centrality, histoName))->Fill(pTot, dca);
244 ((
TH2*)GetMergeableObject(physSel, trigClassName, centrality, histoName))->Fill(pTot, chiProb);
248 fMuonTrackCuts->CustomParam()->SetNSigmaPdca(
fSigmaCuts[isigma]);
249 if ( ! fMuonTrackCuts->IsSelected(track) )
continue;
250 for (
Int_t itrig=0; itrig<selectTrigClasses.GetEntries(); ++itrig ) {
251 TString trigClassName = ((TObjString*)selectTrigClasses.At(itrig))->GetString();
253 ((
TH2*)GetMergeableObject(physSel, trigClassName, centrality, histoName))->Fill(pt, isigma+1);
255 ((
TH2*)GetMergeableObject(physSel, trigClassName, centrality, histoName))->Fill(eta, isigma+1);
265 if ( ! sigmaCuts || nSigmaCuts < 0 ) {
272 Double_t cuts[] = {2., 3., 4., 5., 6., 7., 10., 15., 20., 25., 30., 1.e10};
273 Int_t nCuts =
sizeof(cuts)/
sizeof(cuts[0]);
288 AliVAnalysisMuon::Terminate(
"");
290 if ( ! fMergeableCollection )
return;
292 AliMuonAnalysisOutput muonOut(fOutputList);
294 TString physSel = fTerminateOptions->At(0)->GetName();
295 TString trigClassName = fTerminateOptions->At(1)->GetName();
296 TString centralityRange = fTerminateOptions->At(2)->GetName();
297 TString furtherOpt = fTerminateOptions->At(3)->GetName();
298 furtherOpt.ToUpper();
300 furtherOpt.ReplaceAll(
" ",
" ");
301 furtherOpt.ReplaceAll(
" =",
"=");
302 furtherOpt.ReplaceAll(
"= ",
"=");
303 TObjArray* optArray = furtherOpt.Tokenize(
" ");
305 for (
Int_t iopt=0; iopt<optArray->GetEntries(); ++iopt ) {
306 TString currOpt = optArray->At(iopt)->GetName();
307 if ( currOpt.Contains(
"REFSIGMA") ) {
308 currOpt.Remove(0,currOpt.Index(
"=")+1);
309 refSigmaCut = currOpt.Atof();
314 fMuonTrackCuts->Print(
"param");
316 Int_t srcColors[kNtrackSources] = {kBlack, kRed, kSpring, kTeal, kBlue, kViolet, kMagenta, kOrange, kGray};
329 TString histoName =
"", currName =
"", histoPattern =
"", drawOpt =
"";
330 currName = Form(
"%s_recoDCA", GetName());
331 can =
new TCanvas(currName.Data(),
"Reco DCA",igroup1*xshift,igroup2*yshift,600,600);
335 TString dcaName[2] = {
"DCAx",
"DCAy"};
337 printf(
"\nAverage reconstructed DCA:\n");
338 TF1* fitFuncMeanDca =
new TF1(
"fitFuncMeanDca",
"gausn",-20.,20.);
339 fitFuncMeanDca->SetParNames(
"Norm",
"Mean",
"Sigma");
341 for (
Int_t ihisto=0; ihisto<2; ++ihisto ) {
345 histo =
static_cast<TH2*
>(muonOut.GetSum(physSel, trigClassName, centralityRange, histoPattern));
346 if ( ! histo )
continue;
348 TH1* meanDcaVsP = histo->ProjectionX(Form(
"mean%sVsP_%s", dcaName[ihisto].
Data(),
fThetaAbsKeys->At(itheta)->GetName()));
350 meanDcaVsP->SetTitle(Form(
"Mean %s vs. p %s", dcaName[ihisto].
Data(),
fThetaAbsKeys->At(itheta)->GetName()));
351 meanDcaVsP->SetYTitle(Form(
"<%s> (cm)", dcaName[ihisto].
Data()));
352 meanDcaVsP->SetStats(kFALSE);
354 Int_t nPbins = histo->GetXaxis()->GetNbins();
358 TCanvas* meanDcaFitCan = 0x0;
364 for (
Int_t ibin=2; ibin<=nPbins; ++ibin ) {
365 currName = Form(
"hMean%s_%s_%s", dcaName[ihisto].
Data(), physSel.Data(), trigClassName.Data());
366 Int_t minBin = ( ibin == 0 ) ? 1 : ibin;
367 Int_t maxBin = ( ibin == 0 ) ? nPbins : ibin;
368 if ( ibin > 0 ) currName += Form(
"_pBin%i", ibin);
369 TH1* projHisto = histo->ProjectionY(currName.Data(), minBin, maxBin,
"e");
370 projHisto->SetTitle(Form(
"%s %s %g < p < %g (GeV/c)", dcaName[ihisto].
Data(),
fThetaAbsKeys->At(itheta)->GetName(), meanDcaVsP->GetXaxis()->GetBinLowEdge(minBin), meanDcaVsP->GetXaxis()->GetBinUpEdge(maxBin)));
372 if ( projHisto->GetEntries() == 0 )
continue;
373 if ( ipad % (nPadX*nPadY) == 0 ) {
374 currName = histo->GetName();
375 currName += Form(
"Fit_can_%i", ican++);
376 meanDcaFitCan =
new TCanvas(currName.Data(), currName.Data(),igroup1*xshift,igroup2*yshift,600,600);
377 meanDcaFitCan->Divide(nPadX,nPadY);
380 meanDcaFitCan->cd(++ipad);
382 if ( projHisto->Integral() > 50 ) {
383 fitFuncMeanDca->SetParameter(0, projHisto->Integral());
384 fitFuncMeanDca->SetParameter(1, projHisto->GetMean());
385 fitFuncMeanDca->SetParameter(2, projHisto->GetRMS());
387 fitDcaLim = TMath::Max(5., fitDcaLim);
388 projHisto->Fit(fitFuncMeanDca,
"RQ",
"e", -fitDcaLim, fitDcaLim);
389 Double_t chi2 = fitFuncMeanDca->GetChisquare();
390 Double_t ndf = fitFuncMeanDca->GetNDF();
391 if ( ndf <= 0.)
continue;
392 if ( chi2 / ndf > 100. )
continue;
393 Double_t meanDca = fitFuncMeanDca->GetParameter(1);
394 Double_t meanDcaErr = fitFuncMeanDca->GetParError(1);
395 meanDcaVsP->SetBinContent(ibin, meanDca);
396 meanDcaVsP->SetBinError(ibin, meanDcaErr);
399 else projHisto->Draw(
"e");
401 can->cd(2*itheta+ihisto+1);
405 meanDcaVsP->Fit(
"pol0",
"Q");
406 TF1* trendFit = (TF1*)meanDcaVsP->GetListOfFunctions()->FindObject(
"pol0");
408 averageDca[2*itheta+ihisto] = trendFit->GetParameter(0);
409 printf(
" %s: mean %s = %g +- %g (chi2/%i = %g)\n",
fThetaAbsKeys->At(itheta)->GetName(), dcaName[ihisto].Data(), trendFit->GetParameter(0), trendFit->GetParError(0), trendFit->GetNDF(), trendFit->GetChisquare()/((
Double_t)trendFit->GetNDF()));
424 printf(
"muonCuts->SetMeanDCA(%g, %g, 0.); // %s\n", averageDca[2*itheta], averageDca[2*itheta+1],
fThetaAbsKeys->At(itheta)->GetName());
430 Double_t nSigmaCut = fMuonTrackCuts->GetMuonTrackCutsParam().GetNSigmaPdca();
431 Double_t sigmaMeasCut[2] = { fMuonTrackCuts->GetMuonTrackCutsParam().GetSigmaPdca23(), fMuonTrackCuts->GetMuonTrackCutsParam().GetSigmaPdca310()};
432 Double_t relPResolution = fMuonTrackCuts->GetMuonTrackCutsParam().GetRelPResolution();
433 Double_t angleResolution = 535.*fMuonTrackCuts->GetMuonTrackCutsParam().GetSlopeResolution();
436 const Int_t kNCutFuncs = 2;
437 Int_t nShowFuncs = 1;
438 TString cutFormula =
"[1]*TMath::Sqrt( ( [0] / ( 1. - [1]*[2]*x / ( 1.+[1]*[2]*x ) ) ) * ( [0] / ( 1. - [1]*[2]*x / ( 1.+[1]*[2]*x ) ) ) + [3]*[3]*x*x)";
439 Double_t cutParam[kNCutFuncs][4] = {{sigmaMeasCut[0], nSigmaCut, relPResolution, angleResolution}, {sigmaMeasCut[0], nSigmaCut, 0., 0.32}};
440 Int_t cutColor[kNCutFuncs] = {kBlack, kRed};
443 can =
new TCanvas(
"pdcaSigmaFit",
"Sigma vs P fit",igroup1*xshift,igroup2*yshift,600,600);
446 TF1* fitFunc =
new TF1(
"fitFunc",
"x*gausn", 0., 400.);
447 fitFunc->SetParNames(
"Norm",
"Mean",
"Sigma");
448 gStyle->SetOptFit(1111);
451 printf(
"\nSigma p x DCA:\n");
454 TLegend* leg =
new TLegend(0.7, 0.7, 0.9, 0.9);
455 leg->SetBorderSize(1);
456 for (
Int_t isrc=0; isrc<kNtrackSources; ++isrc ) {
458 TH2* histo =
static_cast<TH2*
>(muonOut.GetSum(physSel, trigClassName, centralityRange, histoPattern));
459 if ( ! histo )
continue;
460 if ( histo->Integral() < 200 ) {
467 sigmaVsP->SetTitle(Form(
"#sigma_{p#timesDCA} vs. p %s %s",
fThetaAbsKeys->At(itheta)->GetName(), fSrcKeys->At(isrc)->GetName()));
468 sigmaVsP->SetYTitle(
"#sigma_{p#timesDCA} (cm #times GeV/c)");
469 sigmaVsP->SetStats(kFALSE);
471 Int_t nPbins = histo->GetXaxis()->GetNbins();
475 TCanvas* pdcaFitCan = 0x0;
481 for (
Int_t ibin=2; ibin<=nPbins; ++ibin ) {
482 currName = Form(
"hPDCA_%s_%s_%s_%s", fSrcKeys->At(isrc)->GetName(), physSel.Data(), trigClassName.Data(),
fThetaAbsKeys->At(itheta)->GetName());
483 Int_t minBin = ( ibin == 0 ) ? 1 : ibin;
484 Int_t maxBin = ( ibin == 0 ) ? nPbins : ibin;
485 if ( ibin > 0 ) currName += Form(
"_pBin%i", ibin);
486 TH1* projHisto = histo->ProjectionY(currName.Data(), minBin, maxBin,
"e");
487 if ( projHisto->GetEntries() == 0 )
continue;
488 projHisto->SetTitle(Form(
"P DCA %s %s %g < p < %g (GeV/c)",
fThetaAbsKeys->At(itheta)->GetName(), fSrcKeys->At(isrc)->GetName(), sigmaVsP->GetXaxis()->GetBinLowEdge(minBin), sigmaVsP->GetXaxis()->GetBinUpEdge(maxBin)));
489 if ( ipad % (nPadX*nPadY) == 0 ) {
490 currName = histo->GetName();
491 currName += Form(
"Fit_can_%i", ican++);
492 pdcaFitCan =
new TCanvas(currName.Data(), currName.Data(),igroup1*xshift,igroup2*yshift,600,600);
493 pdcaFitCan->Divide(nPadX,nPadY);
496 pdcaFitCan->cd(++ipad);
497 if ( projHisto->Integral() > 0.) projHisto->Scale(1./projHisto->Integral());
499 if ( projHisto->GetEntries() > 50 ) {
500 fitFunc->SetParameter(0, projHisto->Integral());
501 fitFunc->SetParameter(1, projHisto->GetMean());
502 fitFunc->SetParameter(2, projHisto->GetRMS());
503 projHisto->Fit(fitFunc,
"RQ",
"e", xMinFit[itheta], xMaxFit[itheta]);
504 Double_t chi2 = fitFunc->GetChisquare();
506 if ( ndf <= 0.)
continue;
507 if ( chi2 / ndf > 100. )
continue;
509 Double_t sigmaErr = fitFunc->GetParError(2);
510 sigmaVsP->SetBinContent(ibin, sigma);
511 sigmaVsP->SetBinError(ibin, sigmaErr);
514 else projHisto->Draw(
"e");
517 sigmaVsP->SetLineColor(srcColors[isrc]);
518 sigmaVsP->SetMarkerColor(srcColors[isrc]);
519 sigmaVsP->SetMarkerStyle(20+isrc);
521 if ( gPad->GetListOfPrimitives()->GetEntries() != 0 ) drawOpt +=
"same";
522 sigmaVsP->Draw(drawOpt.Data());
523 sigmaVsP->Fit(
"pol0",
"Q",drawOpt.Data());
524 TF1* trendFit = (TF1*)sigmaVsP->GetListOfFunctions()->FindObject(
"pol0");
526 averageSigmaPdca[kNtrackSources*itheta+isrc] = trendFit->GetParameter(0);
527 printf(
" %s %s: Sigma = %g +- %g (chi2/%i = %g)\n",
fThetaAbsKeys->At(itheta)->GetName(), fSrcKeys->At(isrc)->GetName(), trendFit->GetParameter(0), trendFit->GetParError(0), trendFit->GetNDF(), trendFit->GetChisquare()/((
Double_t)trendFit->GetNDF()));
529 leg->AddEntry(sigmaVsP, fSrcKeys->At(isrc)->GetName(),
"lp");
534 for (
Int_t icut=0; icut<nShowFuncs; ++icut ) {
535 currName = Form(
"sigmaCut_%s_%i",
fThetaAbsKeys->At(itheta)->GetName(), icut);
536 TF1* cutFunction =
new TF1(currName.Data(), cutFormula.Data(), pMinCut, pMaxCut);
537 cutParam[icut][0] = sigmaMeasCut[itheta];
538 cutParam[icut][1] = 1.;
539 cutFunction->SetParameters(cutParam[icut]);
540 cutFunction->SetLineColor(cutColor[icut]);
541 cutFunction->SetLineWidth(2);
542 cutFunction->Draw(
"same");
546 for (
Int_t isrc=0; isrc<kNtrackSources; ++isrc ) {
547 printf(
"muonCuts->SetSigmaPdca(%g, %g); // %s\n", averageSigmaPdca[isrc], averageSigmaPdca[kNtrackSources+isrc], fSrcKeys->At(isrc)->GetName());
552 for (
Int_t icheck=0; icheck<2; icheck++ ) {
555 for (
Int_t isrc=0; isrc<kNtrackSources; ++isrc ) {
557 TH2* histoCheck =
static_cast<TH2*
>(muonOut.GetSum(physSel, trigClassName, centralityRange, histoPattern));
558 if ( ! histoCheck )
continue;
559 currName = histoCheck->GetName();
560 currName.Append(
"_plotCut");
561 TCanvas* pdcaVsPcan =
new TCanvas(currName.Data(), currName.Data(), igroup1*xshift,(igroup2+icheck)*yshift,600,600);
562 pdcaVsPcan->SetRightMargin(0.12);
563 pdcaVsPcan->SetLogy();
564 pdcaVsPcan->SetLogz();
565 histoCheck->Draw(
"COLZ");
567 for (
Int_t icut=0; icut<nShowFuncs; ++icut ) {
568 currName = histoCheck->GetName();
569 currName.Append(Form(
"_cutFunc%i",icut));
570 TString currFormula = cutFormula;
571 if ( icheck == 1 ) currFormula.Append(
"/x");
572 TF1* cutFunction =
new TF1(currName.Data(),currFormula.Data(), pMinCut, pMaxCut);
573 cutParam[icut][0] = sigmaMeasCut[itheta];
574 cutParam[icut][1] = nSigmaCut;
575 cutFunction->SetParameters(cutParam[icut]);
576 cutFunction->SetLineWidth(2);
577 cutFunction->SetLineColor(cutColor[icut]);
578 cutFunction->Draw(
"same");
591 for (
Int_t isrc=0; isrc<kNtrackSources; ++isrc ) {
593 TH2* histo =
static_cast<TH2*
>(muonOut.GetSum(physSel, trigClassName, centralityRange, histoPattern));
594 if ( ! histo )
continue;
596 Int_t nPbins = histo->GetXaxis()->GetNbins();
600 TCanvas* chi2probCan = 0x0;
606 for (
Int_t ibin=0; ibin<=nPbins; ++ibin ) {
607 currName = Form(
"hChiProb_%s_%s_%s_%s", fSrcKeys->At(isrc)->GetName(), physSel.Data(), trigClassName.Data(),
fThetaAbsKeys->At(itheta)->GetName());
608 Int_t minBin = ( ibin == 0 ) ? 1 : ibin;
609 Int_t maxBin = ( ibin == 0 ) ? nPbins : ibin;
610 if ( ibin > 0 ) currName += Form(
"_pBin%i", ibin);
611 TH1* projHisto = histo->ProjectionY(currName.Data(), minBin, maxBin);
612 projHisto->SetTitle(Form(
"Chisquare prob %s %s %g < p < %g (GeV/c)",
fThetaAbsKeys->At(itheta)->GetName(), fSrcKeys->At(isrc)->GetName(), histo->GetXaxis()->GetBinLowEdge(minBin), histo->GetXaxis()->GetBinUpEdge(maxBin)));
613 if ( projHisto->GetEntries() == 0 )
continue;
614 if ( ipad % (nPadX*nPadY) == 0 ) {
615 currName = histo->GetName();
616 currName += Form(
"Fit_can_%i", ican++);
617 chi2probCan =
new TCanvas(currName.Data(), currName.Data(),igroup1*xshift,igroup2*yshift,600,600);
618 chi2probCan->Divide(nPadX,nPadY);
621 chi2probCan->cd(++ipad);
623 projHisto->Draw(
"e");
632 printf(
"\nReference sigma cut %g\n", refSigmaCut);
636 Int_t cutColors[13] = {kBlack, kRed, kBlue, kGreen, kCyan, kMagenta, kOrange, kViolet, kSpring, kGray, kPink, kAzure, kYellow};
637 Int_t cutStyles[13] = {20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 32, 33};
640 Bool_t useCustomSigma = furtherOpt.Contains(
"CUSTOMSIGMA");
645 for (
Int_t ihisto=0; ihisto<2; ++ihisto ) {
646 for (
Int_t isrc=0; isrc<kNtrackSources; ++isrc ) {
649 histoPattern =
GetHistoName(checkHistos[ihisto], itheta, isrc);
650 TH2* histo =
static_cast<TH2*
>(muonOut.GetSum(physSel, trigClassName, centralityRange, histoPattern));
651 if ( ! histo )
continue;
652 if ( histo->Integral() == 0. ) {
656 if ( orderCuts.GetSize() == 0 ) {
658 TAxis* axis = histo->GetYaxis();
659 Int_t nSigmaCuts = ( useCustomSigma ) ?
fSigmaCuts.GetSize() : axis->GetNbins();
660 orderCuts.Set(nSigmaCuts);
662 for (
Int_t isigma=0; isigma<axis->GetNbins(); ++isigma ) {
663 TString currLabel = axis->GetBinLabel(isigma+1);
664 Double_t currSigma = currLabel.Atof();
665 if ( useCustomSigma ) {
666 Bool_t foundMatch = kFALSE;
668 if ( TMath::Abs(currSigma -
fSigmaCuts[jsigma]) < 1e-4 ) {
673 if ( ! foundMatch )
continue;
675 Int_t currBin = ( TMath::Abs(currSigma - refSigmaCut) < 1e-4) ? 0 : ++countBin;
676 orderCuts[currBin] = isigma+1;
679 currName = histo->GetName();
680 currName.Append(
"_can");
681 can =
new TCanvas(currName.Data(), currName.Data(), igroup1*xshift,igroup2*yshift, 600, 600);
682 TLegend* leg =
new TLegend(0.6,0.6,0.9,0.9);
683 leg->SetBorderSize(1);
684 can->Divide(1,2,0,0);
686 gPad->SetPad(0., fracOfHeight, 0.99, 0.99);
688 gPad->SetTopMargin(0.03);
689 gPad->SetRightMargin(rightMargin);
693 gPad->SetPad(0., 0., 0.99, fracOfHeight);
694 gPad->SetRightMargin(rightMargin);
695 gPad->SetBottomMargin(0.08/fracOfHeight);
697 TH1* refCutHisto = 0x0;
699 for (
Int_t isigma=0; isigma<orderCuts.GetSize(); ++isigma ) {
700 currName = histo->GetName();
701 currName.Append(Form(
"_sigma%i", isigma));
702 Int_t currBin = orderCuts[isigma];
703 TH1* projectHisto = histo->ProjectionX(currName.Data(), currBin, currBin);
704 projectHisto->SetLineColor(cutColors[isigma]);
705 projectHisto->SetMarkerColor(cutColors[isigma]);
706 projectHisto->SetMarkerStyle(cutStyles[isigma]);
707 projectHisto->SetStats(kFALSE);
708 projectHisto->SetTitle(
"");
711 if ( gPad->GetListOfPrimitives()->GetEntries() != 0 ) drawOpt +=
"same";
712 projectHisto->Draw(drawOpt.Data());
713 TString currLabel = histo->GetYaxis()->GetBinLabel(currBin);
714 TString cutName = ( TMath::Abs(currLabel.Atof() ) > 1000. ) ?
"Total" : Form(
"%s #sigma cut", currLabel.Data());
715 leg->AddEntry(projectHisto, cutName.Data(),
"lp");
716 if ( ! refCutHisto ) {
717 refCutHisto = projectHisto;
718 legTitle = Form(
"(pass n #sigma) / (pass %s #sigma)", currLabel.Data());
721 currName.Append(
"_ratio");
722 TH1* histoRatio = (
TH1*)projectHisto->Clone(currName.Data());
724 TH1* auxNum = projectHisto;
725 TH1* auxDen = refCutHisto;
726 Bool_t mustInvert = kFALSE;
727 if ( auxNum->Integral()>auxDen->Integral() ) {
728 auxNum = refCutHisto;
729 auxDen = projectHisto;
732 histoRatio->Divide(auxNum,auxDen,1.,1.,
"B");
734 TH1* auxHistoOne =
static_cast<TH1*
>(histoRatio->Clone(
"auxHistoOne"));
735 for (
Int_t ibin=1; ibin<=auxHistoOne->GetXaxis()->GetNbins(); ibin++ ) {
736 auxHistoOne->SetBinContent(ibin,1.);
737 auxHistoOne->SetBinError(ibin,0.);
739 histoRatio->Divide(auxHistoOne,histoRatio);
742 histoRatio->SetTitle(
"");
743 histoRatio->GetYaxis()->SetTitle(legTitle.Data());
744 histoRatio->GetXaxis()->SetLabelSize(0.04/fracOfHeight);
745 histoRatio->GetXaxis()->SetTitleSize(0.035/fracOfHeight);
746 histoRatio->GetYaxis()->SetLabelSize(0.03/fracOfHeight);
747 histoRatio->GetYaxis()->SetTitleSize(0.03/fracOfHeight);
748 histoRatio->GetXaxis()->SetTitleOffset(1.);
749 histoRatio->GetYaxis()->SetTitleOffset(0.6);
750 histoRatio->GetYaxis()->SetRangeUser(0.5,1.5);
754 if ( gPad->GetListOfPrimitives()->GetEntries() != 0 ) drawOpt +=
"same";
755 histoRatio->Draw(drawOpt.Data());
766 Int_t nPtMins =
sizeof(
ptMin)/
sizeof(ptMin[0]);
767 for (
Int_t iptmin=0; iptmin<nPtMins; ++iptmin) {
768 for (
Int_t isrc=0; isrc<kNtrackSources; ++isrc ) {
773 for (
Int_t icent=1; icent<=GetCentralityClasses()->GetNbins(); ++icent ) {
774 TH2* histo =
static_cast<TH2*
>(muonOut.GetSum(physSel, trigClassName, GetCentralityClasses()->GetBinLabel(icent), histoName));
775 if ( ! histo )
continue;
776 Int_t ptMinBin = histo->GetXaxis()->FindBin(ptMin[iptmin]);
777 Int_t ptMaxBin = histo->GetXaxis()->GetNbins()+1;
778 currName = histo->GetName();
779 currName += Form(
"_%s_ptMin%i", GetCentralityClasses()->GetBinLabel(icent), TMath::Nint(ptMin[iptmin]));
780 TH1* projectHisto = histo->ProjectionY(currName.Data(), ptMinBin, ptMaxBin);
781 if ( projectHisto->GetMaximum() < 50. ) {
786 currName = Form(
"CutEffect_%s_ptMin%i", fSrcKeys->At(isrc)->GetName(), TMath::Nint(ptMin[iptmin]));
787 can =
new TCanvas(currName.Data(), currName.Data(), igroup1*xshift,igroup2*yshift, 600, 600);
791 leg =
new TLegend(0.6,0.6,0.9,0.9);
792 leg->SetBorderSize(1);
795 projectHisto->SetTitle(Form(
"Cut effect %s %s %g < p_{t} < %g (GeV/c)", fSrcKeys->At(isrc)->GetName(),
fThetaAbsKeys->At(itheta)->GetName(), histo->GetXaxis()->GetBinLowEdge(ptMinBin), histo->GetXaxis()->GetBinUpEdge(ptMaxBin)));
796 projectHisto->SetLineColor(cutColors[icent-1]);
797 projectHisto->SetMarkerColor(cutColors[icent-1]);
798 projectHisto->SetMarkerStyle(19+icent);
799 projectHisto->SetStats(0);
801 if ( gPad->GetListOfPrimitives()->GetEntries() != 0 ) drawOpt +=
"same";
802 projectHisto->Draw(drawOpt.Data());
803 leg->AddEntry(projectHisto, GetCentralityClasses()->GetBinLabel(icent),
"lp");
804 Double_t keptEvents = projectHisto->GetBinContent(orderCuts[0]);
805 Double_t totalEvents = projectHisto->GetBinContent(orderCuts[orderCuts.GetSize()-1]);
806 Double_t accepted = ( totalEvents > 0. ) ? keptEvents / totalEvents : 1.;
807 Double_t acceptedErr = ( totalEvents > 0. ) ? TMath::Sqrt(accepted*(1.-accepted)/totalEvents) : 1.;
808 printf(
"%12s %11s %6s (pt>%g) rejected evts : %6.2f +- %6.3f %%\n", fSrcKeys->At(isrc)->GetName(),
fThetaAbsKeys->At(itheta)->GetName(), GetCentralityClasses()->GetBinLabel(icent), ptMin[iptmin], (1.-accepted)*100., acceptedErr*100.);
811 if ( leg ) leg->Draw(
"same");
void MyUserCreateOutputObjects()
TObjArray * fThetaAbsKeys
Name of theta at absorber end.
AliAnalysisTaskMuonCuts()
p x DCA vs momentum (binning for fit)
virtual void Terminate(Option_t *option)
void ProcessEvent(TString physSel, const TObjArray &selectTrigClasses, TString centrality)
Chi square probability vs momentum.
TObjArray * fHistoTypeKeys
Base histogram name.
Bool_t Data(TH1F *h, Double_t *rangefit, Bool_t writefit, Double_t &sgn, Double_t &errsgn, Double_t &bkg, Double_t &errbkg, Double_t &sgnf, Double_t &errsgnf, Double_t &sigmafit, Int_t &status)
virtual ~AliAnalysisTaskMuonCuts()
eta distribution for different p x DCA sigma cuts
TArrayD fSigmaCuts
List of sigma cuts.
p x DCA vs momentum (check beam gas)
pt distribution for different p x DCA sigma cuts
void SetSigmaCuts(Int_t nSigmaCuts=-1, Double_t *sigmaCuts=0x0)
TString GetHistoName(Int_t histoTypeIndex, Int_t thetaAbsIndex, Int_t srcIndex)