1 /**************************************************************************
2 * Copyright(c) 1998-2007, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
16 /* $Id: AliAnalysisTaskMuonCuts.cxx 47782 2011-02-24 18:37:31Z martinez $ */
18 /// Task for tuning cuts for single muons in the spectrometer.
19 /// The output is a list of histograms.
20 /// Task need to be run in steps:
21 /// Step 1: find <DCAx> and <DCAy> in MC
22 /// Step 2: find signal sigma_pDCA in MC
23 /// Step 3: find <DCAx> and <DCAy> in data
24 /// Step 4: apply sigma_pDCA cut on simulations and check the
25 /// resulting kinematic variables
27 /// \author Diego Stocco
29 #define AliAnalysisTaskMuonCuts_cxx
31 #include "AliAnalysisTaskMuonCuts.h"
42 #include "TObjString.h"
43 #include "TObjArray.h"
46 //#include "TMCProcess.h"
49 #include "AliAODEvent.h"
50 #include "AliAODTrack.h"
51 #include "AliAODMCParticle.h"
52 #include "AliMCEvent.h"
53 #include "AliMCParticle.h"
54 //#include "AliStack.h"
55 #include "AliESDEvent.h"
56 #include "AliESDMuonTrack.h"
57 #include "AliInputEventHandler.h"
59 #include "AliAnalysisManager.h"
62 #include "AliMergeableCollection.h"
63 #include "AliMuonTrackCuts.h"
67 ClassImp(AliAnalysisTaskMuonCuts) // Class implementation in ROOT context
71 //________________________________________________________________________
72 AliAnalysisTaskMuonCuts::AliAnalysisTaskMuonCuts() :
81 //________________________________________________________________________
82 AliAnalysisTaskMuonCuts::AliAnalysisTaskMuonCuts(const char *name, const AliMuonTrackCuts& cuts) :
83 AliVAnalysisMuon(name, cuts),
91 TString histoTypeKeys = "hDCAxVsP hDCAyVsP hPdcaVsP hPDCAVsPCheck hDCAVsPCheck hChiProbVsP hSigmaVsPt hSigmaVsEta";
92 fHistoTypeKeys = histoTypeKeys.Tokenize(" ");
94 TString thetaAbsKeys = "ThetaAbs23 ThetaAbs310";
95 fThetaAbsKeys = thetaAbsKeys.Tokenize(" ");
101 //________________________________________________________________________
102 AliAnalysisTaskMuonCuts::~AliAnalysisTaskMuonCuts()
108 delete fHistoTypeKeys;
109 delete fThetaAbsKeys;
112 //___________________________________________________________________________
113 void AliAnalysisTaskMuonCuts::MyUserCreateOutputObjects()
116 TString histoName = "";
117 for ( Int_t itheta=0; itheta<kNthetaAbs; ++itheta ) {
118 for ( Int_t isrc=0; isrc<kNtrackSources; ++isrc ) {
119 histoName = GetHistoName(kDCAxVsP, itheta, isrc);
120 histo = new TH2F(histoName.Data(), histoName.Data(), 50, 0., 200., 100, -50., 50.);
121 histo->SetXTitle("p (GeV/c)");
122 histo->SetYTitle("DCA_{x} (cm)");
123 AddObjectToCollection(histo);
125 histoName = GetHistoName(kDCAyVsP, itheta, isrc);
126 histo = new TH2F(histoName.Data(), histoName.Data(), 50, 0., 200., 100, -50., 50.);
127 histo->SetXTitle("p (GeV/c)");
128 histo->SetYTitle("DCA_{y} (cm)");
129 AddObjectToCollection(histo);
131 histoName = GetHistoName(kPdcaVsP, itheta, isrc);
132 histo = new TH2F(histoName.Data(), histoName.Data(), 100, 0., 400., 50, 0., 500.);
133 histo->SetXTitle("p (GeV/c)");
134 histo->SetYTitle("p#timesDCA (cm #times GeV/c)");
135 AddObjectToCollection(histo);
137 histoName = GetHistoName(kPDCAVsPCheck, itheta, isrc);
138 histo = new TH2F(histoName.Data(), histoName.Data(), 100, 0., 800., 100, 0., 5000.);
139 histo->SetXTitle("p (GeV/c)");
140 histo->SetYTitle("p#timesDCA (cm #times GeV/c)");
141 AddObjectToCollection(histo);
143 histoName = GetHistoName(kDCAVsPCheck, itheta, isrc);
144 histo = new TH2F(histoName.Data(), histoName.Data(), 100, 0., 800., 100, 0., 200.);
145 histo->SetXTitle("p (GeV/c)");
146 histo->SetYTitle("DCA (cm)");
147 AddObjectToCollection(histo);
149 histoName = GetHistoName(kChiProbVsP, itheta, isrc);
150 histo = new TH2F(histoName.Data(), histoName.Data(), 50, 0., 200., 50, 0., 1.);
151 histo->SetXTitle("p (GeV/c)");
152 histo->SetYTitle("Chisquare prob.");
153 AddObjectToCollection(histo);
155 histoName = GetHistoName(kSigmaVsPt, itheta, isrc);
156 histo = new TH2F(histoName.Data(), histoName.Data(), 100, 0., 100., fSigmaCuts.GetSize(), 0.5, 0.5+(Double_t)fSigmaCuts.GetSize());
157 histo->SetXTitle("p_{t} (GeV/c)");
158 histo->SetYTitle("#sigma_{p#timesDCA}");
159 for ( Int_t ibin=0; ibin<fSigmaCuts.GetSize(); ++ibin ) {
160 histo->GetYaxis()->SetBinLabel(ibin+1,Form("%g", fSigmaCuts[ibin]));
162 AddObjectToCollection(histo);
164 histoName = GetHistoName(kSigmaVsEta, itheta, isrc);
165 histo = new TH2F(histoName.Data(), histoName.Data(), 25, -4.5, -2., fSigmaCuts.GetSize(), 0.5, 0.5+(Double_t)fSigmaCuts.GetSize());
166 histo->SetXTitle("#eta");
167 histo->SetYTitle("#sigma_{p#timesDCA}");
168 for ( Int_t ibin=0; ibin<fSigmaCuts.GetSize(); ++ibin ) {
169 histo->GetYaxis()->SetBinLabel(ibin+1,Form("%g", fSigmaCuts[ibin]));
171 AddObjectToCollection(histo);
172 } // loop on track sources
173 } // loop on theta abs
175 fMuonTrackCuts->Print();
179 //________________________________________________________________________
180 TString AliAnalysisTaskMuonCuts::GetHistoName(Int_t histoTypeIndex, Int_t thetaAbsIndex, Int_t srcIndex)
182 /// Get local histogram name
183 TString histoName = Form("%s%s%s", fHistoTypeKeys->At(histoTypeIndex)->GetName(), fThetaAbsKeys->At(thetaAbsIndex)->GetName(), fSrcKeys->At(srcIndex)->GetName());
188 //________________________________________________________________________
189 void AliAnalysisTaskMuonCuts::ProcessEvent(TString physSel, const TObjArray& selectTrigClasses, TString centrality)
195 if ( GetVertexSPD()->GetNContributors() < fMinNvtxContirbutors ) return;
197 TString histoName = "";
198 AliVParticle* track = 0x0;
199 Int_t nTracks = GetNTracks();
200 for (Int_t itrack = 0; itrack < nTracks; itrack++) {
201 track = GetTrack(itrack);
202 fMuonTrackCuts->SetNSigmaPdca(1.e10);
203 if ( ! fMuonTrackCuts->IsSelected(track) ) continue;
205 Double_t rAbsEnd = ( fAODEvent ) ? ((AliAODTrack*)track)->GetRAtAbsorberEnd(): ((AliESDMuonTrack*)track)->GetRAtAbsorberEnd();
206 Double_t thetaAbsEndDeg = TMath::ATan( rAbsEnd / 505. ) * TMath::RadToDeg();
207 Int_t thetaAbsBin = ( thetaAbsEndDeg < 3. ) ? kThetaAbs23 : kThetaAbs310;
209 Int_t trackSrc = GetParticleType(track);
211 TVector3 dcaAtVz = fMuonTrackCuts->GetCorrectedDCA(track);
212 Double_t pTotMean = fMuonTrackCuts->GetAverageMomentum(track);
213 Double_t dca = dcaAtVz.Mag();
214 Double_t pDca = pTotMean * dca;
216 Double_t chi2 = pDca / fMuonTrackCuts->GetSigmaPdca(rAbsEnd) ;
218 Double_t chiProb = TMath::Prob(chi2, 2);
220 Double_t pTot = track->P();
221 Double_t pt = track->Pt();
222 Double_t eta = track->Eta();
224 for ( Int_t itrig=0; itrig<selectTrigClasses.GetEntries(); ++itrig ) {
225 TString trigClassName = ((TObjString*)selectTrigClasses.At(itrig))->GetString();
226 histoName = GetHistoName(kDCAxVsP, thetaAbsBin, trackSrc);
227 ((TH2*)GetMergeableObject(physSel, trigClassName, centrality, histoName))->Fill(pTot, dcaAtVz.X());
229 histoName = GetHistoName(kDCAyVsP, thetaAbsBin, trackSrc);
230 ((TH2*)GetMergeableObject(physSel, trigClassName, centrality, histoName))->Fill(pTot, dcaAtVz.Y());
232 histoName = GetHistoName(kPdcaVsP, thetaAbsBin, trackSrc);
233 ((TH2*)GetMergeableObject(physSel, trigClassName, centrality, histoName))->Fill(pTot, pDca);
235 histoName = GetHistoName(kPDCAVsPCheck, thetaAbsBin, trackSrc);
236 ((TH2*)GetMergeableObject(physSel, trigClassName, centrality, histoName))->Fill(pTot, pDca);
238 histoName = GetHistoName(kDCAVsPCheck, thetaAbsBin, trackSrc);
239 ((TH2*)GetMergeableObject(physSel, trigClassName, centrality, histoName))->Fill(pTot, dca);
241 histoName = GetHistoName(kChiProbVsP, thetaAbsBin, trackSrc);
242 ((TH2*)GetMergeableObject(physSel, trigClassName, centrality, histoName))->Fill(pTot, chiProb);
243 } // loop on selected trigger classes
245 for ( Int_t isigma=0; isigma<fSigmaCuts.GetSize(); ++isigma) {
246 fMuonTrackCuts->SetNSigmaPdca(fSigmaCuts[isigma]);
247 if ( ! fMuonTrackCuts->IsSelected(track) ) continue;
248 for ( Int_t itrig=0; itrig<selectTrigClasses.GetEntries(); ++itrig ) {
249 TString trigClassName = ((TObjString*)selectTrigClasses.At(itrig))->GetString();
250 histoName = GetHistoName(kSigmaVsPt, thetaAbsBin, trackSrc);
251 ((TH2*)GetMergeableObject(physSel, trigClassName, centrality, histoName))->Fill(pt, isigma+1);
252 histoName = GetHistoName(kSigmaVsEta, thetaAbsBin, trackSrc);
253 ((TH2*)GetMergeableObject(physSel, trigClassName, centrality, histoName))->Fill(eta, isigma+1);
254 } // loop on selected trigger classes
259 //________________________________________________________________________
260 void AliAnalysisTaskMuonCuts::SetSigmaCuts(Int_t nSigmaCuts, Double_t* sigmaCuts)
262 /// Set number of sigmas
263 if ( ! sigmaCuts || nSigmaCuts < 0 ) {
264 // if ( defaultChiSquare ) {
265 // Double_t cuts[] = {0.3, 0.1, 0.05, 0.03, 0.01, 0.005, 0.003, 0.001, 0.0005, 0.0003, 0.0001, 0.};
266 // Int_t nCuts = sizeof(cuts)/sizeof(cuts[0]);
267 // fSigmaCuts.Set(nCuts, cuts);
270 Double_t cuts[] = {2., 3., 4., 5., 6., 7., 10., 15., 20., 25., 30., 1.e10};
271 Int_t nCuts = sizeof(cuts)/sizeof(cuts[0]);
272 fSigmaCuts.Set(nCuts, cuts);
276 fSigmaCuts.Set(nSigmaCuts, sigmaCuts);
280 //________________________________________________________________________
281 void AliAnalysisTaskMuonCuts::Terminate(Option_t *) {
283 /// Draw some histogram at the end.
286 AliVAnalysisMuon::Terminate("");
288 if ( ! fMergeableCollection ) return;
290 TString physSel = fTerminateOptions->At(0)->GetName();
291 TString trigClassName = fTerminateOptions->At(1)->GetName();
292 TString centralityRange = fTerminateOptions->At(2)->GetName();
293 TString furtherOpt = fTerminateOptions->At(3)->GetName();
294 furtherOpt.ToUpper();
296 furtherOpt.ReplaceAll(" ", " ");
297 furtherOpt.ReplaceAll(" =", "=");
298 furtherOpt.ReplaceAll("= ", "=");
299 TObjArray* optArray = furtherOpt.Tokenize(" ");
300 Double_t refSigmaCut = 6.;
301 for ( Int_t iopt=0; iopt<optArray->GetEntries(); ++iopt ) {
302 TString currOpt = optArray->At(iopt)->GetName();
303 if ( currOpt.Contains("REFSIGMA") ) {
304 currOpt.Remove(0,currOpt.Index("=")+1);
305 refSigmaCut = currOpt.Atof();
310 Int_t srcColors[kNtrackSources] = {kBlack, kRed, kGreen, kBlue, kViolet, 7, kOrange};
323 TString histoName = "", currName = "", histoPattern = "", drawOpt = "";
324 currName = Form("%s_recoDCA", GetName());
325 can = new TCanvas(currName.Data(),"Reco DCA",igroup1*xshift,igroup2*yshift,600,600);
328 Int_t recoDcaHisto[2] = {kDCAxVsP, kDCAyVsP};
329 TString dcaName[2] = {"DCAx", "DCAy"};
330 Double_t averageDca[4] = {0., 0.};
331 printf("\nAverage reconstructed DCA:\n");
332 TF1* fitFuncMeanDca = new TF1("fitFuncMeanDca","gausn",-20.,20.);
333 fitFuncMeanDca->SetParNames("Norm", "Mean", "Sigma");
334 for ( Int_t itheta=0; itheta<kNthetaAbs; ++itheta ) {
335 for ( Int_t ihisto=0; ihisto<2; ++ihisto ) {
338 histoPattern = Form("%s & %s", fHistoTypeKeys->At(recoDcaHisto[ihisto])->GetName(), fThetaAbsKeys->At(itheta)->GetName());
339 histo = (TH2*)GetSum(physSel, trigClassName, centralityRange, histoPattern);
340 if ( ! histo ) continue;
342 TH1* meanDcaVsP = histo->ProjectionX(Form("mean%sVsP_%s", dcaName[ihisto].Data(), fThetaAbsKeys->At(itheta)->GetName()));
344 meanDcaVsP->SetTitle(Form("Mean %s vs. p %s", dcaName[ihisto].Data(), fThetaAbsKeys->At(itheta)->GetName()));
345 meanDcaVsP->SetYTitle(Form("<%s> (cm)", dcaName[ihisto].Data()));
346 meanDcaVsP->SetStats(kFALSE);
348 Int_t nPbins = histo->GetXaxis()->GetNbins();
349 //Int_t nPadX = (Int_t)TMath::Sqrt(nPbins);
350 //Int_t nPadY = nPadX;
351 //if ( nPadX * nPadY < nPbins ) nPadX++;
352 TCanvas* meanDcaFitCan = 0x0;
358 for ( Int_t ibin=2; ibin<=nPbins; ++ibin ) {
359 currName = Form("hMean%s_%s_%s", dcaName[ihisto].Data(), physSel.Data(), trigClassName.Data());
360 Int_t minBin = ( ibin == 0 ) ? 1 : ibin;
361 Int_t maxBin = ( ibin == 0 ) ? nPbins : ibin;
362 if ( ibin > 0 ) currName += Form("_pBin%i", ibin);
363 TH1* projHisto = histo->ProjectionY(currName.Data(), minBin, maxBin, "e");
364 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)));
366 if ( projHisto->GetEntries() == 0 ) continue;
367 if ( ipad % (nPadX*nPadY) == 0 ) {
368 currName = histo->GetName();
369 currName += Form("Fit_can_%i", ican++);
370 meanDcaFitCan = new TCanvas(currName.Data(), currName.Data(),igroup1*xshift,igroup2*yshift,600,600);
371 meanDcaFitCan->Divide(nPadX,nPadY);
374 meanDcaFitCan->cd(++ipad);
376 if ( projHisto->Integral() > 50 ) {
377 fitFuncMeanDca->SetParameter(0, projHisto->Integral());
378 fitFuncMeanDca->SetParameter(1, projHisto->GetMean());
379 fitFuncMeanDca->SetParameter(2, projHisto->GetRMS());
380 Double_t fitDcaLim = ( ibin == 0 ) ? 40. : 40./((Double_t)ibin);
381 fitDcaLim = TMath::Max(5., fitDcaLim);
382 projHisto->Fit(fitFuncMeanDca, "RQ", "e", -fitDcaLim, fitDcaLim);
383 Double_t chi2 = fitFuncMeanDca->GetChisquare();
384 Double_t ndf = fitFuncMeanDca->GetNDF();
385 if ( ndf <= 0.) continue;
386 if ( chi2 / ndf > 100. ) continue;
387 Double_t meanDca = fitFuncMeanDca->GetParameter(1);
388 Double_t meanDcaErr = fitFuncMeanDca->GetParError(1);
389 meanDcaVsP->SetBinContent(ibin, meanDca);
390 meanDcaVsP->SetBinError(ibin, meanDcaErr);
391 //if ( ibin == 0 ) printf("%s meanDca = %g +- %g\n", fThetaAbsKeys->At(itheta)->GetName(), meanDca, meanDcaErr);
393 else projHisto->Draw("e");
394 } // loop on momentum bins
395 can->cd(2*itheta+ihisto+1);
396 //meanDcaVsP->SetLineColor(srcColors[isrc]);
397 //meanDcaVsP->SetMarkerColor(srcColors[isrc]);
398 //meanDcaVsP->SetMarkerStyle(20+isrc);
399 meanDcaVsP->Fit("pol0","Q");
400 TF1* trendFit = (TF1*)meanDcaVsP->GetListOfFunctions()->FindObject("pol0");
402 averageDca[2*itheta+ihisto] = trendFit->GetParameter(0);
403 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()));
406 //leg->AddEntry(meanDcaVsP, fSrcKeys->At(isrc)->GetName(), "lp");
413 //meanDca[ihisto] = histo->GetMean();
414 } // loop on histo type
415 } // loop on theta abs
417 for ( Int_t itheta=0; itheta<kNthetaAbs; ++itheta ) {
418 printf("muonCuts->SetMeanDCA(%g, %g, 0.); // %s\n", averageDca[2*itheta], averageDca[2*itheta+1], fThetaAbsKeys->At(itheta)->GetName());
424 Double_t nSigmaCut = fMuonTrackCuts->GetNSigmaPdca(); //6.;
425 Double_t sigmaMeasCut[2] = { fMuonTrackCuts->GetSigmaPdca(AliMuonTrackCuts::kThetaAbs23), fMuonTrackCuts->GetSigmaPdca(AliMuonTrackCuts::kThetaAbs23)}; //{99., 54.}; //{120., 63.};
426 Double_t relPResolution = fMuonTrackCuts->GetRelPResolution(); //4.5e-4; //6.e-3;//8.e-4;
427 Double_t angleResolution = 535.*fMuonTrackCuts->GetSlopeResolution(); //6.e-4;
428 Double_t pMinCut = 0.1;
429 Double_t pMaxCut = 800.;
430 const Int_t kNCutFuncs = 2;
431 Int_t nShowFuncs = 1;
432 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)";
433 Double_t cutParam[kNCutFuncs][4] = {{sigmaMeasCut[0], nSigmaCut, relPResolution, angleResolution}, {sigmaMeasCut[0], nSigmaCut, 0., 0.32}};
434 Int_t cutColor[kNCutFuncs] = {kBlack, kRed};
437 can = new TCanvas("pdcaSigmaFit","Sigma vs P fit",igroup1*xshift,igroup2*yshift,600,600);
440 TF1* fitFunc = new TF1("fitFunc", "x*gausn", 0., 400.);
441 fitFunc->SetParNames("Norm", "Mean", "Sigma");
442 gStyle->SetOptFit(1111);
443 Double_t xMinFit[2] = {0., 0.};
444 Double_t xMaxFit[2] = {320., 150.}; // {360., 180.};
445 printf("\nSigma p x DCA:\n");
446 Double_t averageSigmaPdca[kNtrackSources*kNthetaAbs] = {0.};
447 for ( Int_t itheta=0; itheta<kNthetaAbs; ++itheta ) {
448 TLegend* leg = new TLegend(0.7, 0.7, 0.9, 0.9);
449 leg->SetBorderSize(1);
450 for ( Int_t isrc=0; isrc<kNtrackSources; ++isrc ) {
451 histoPattern = GetHistoName(kPdcaVsP, itheta, isrc);
452 TH2* histo = (TH2*)GetSum(physSel, trigClassName, centralityRange, histoPattern);
453 if ( ! histo ) continue;
454 if ( histo->Integral() < 200 ) {
459 TH1* sigmaVsP = histo->ProjectionX(Form("sigma%s_%s_%s", fHistoTypeKeys->At(kPdcaVsP)->GetName(), fThetaAbsKeys->At(itheta)->GetName(), fSrcKeys->At(isrc)->GetName()));
461 sigmaVsP->SetTitle(Form("#sigma_{p#timesDCA} vs. p %s %s", fThetaAbsKeys->At(itheta)->GetName(), fSrcKeys->At(isrc)->GetName()));
462 sigmaVsP->SetYTitle("#sigma_{p#timesDCA} (cm #times GeV/c)");
463 sigmaVsP->SetStats(kFALSE);
465 Int_t nPbins = histo->GetXaxis()->GetNbins();
466 //Int_t nPadX = (Int_t)TMath::Sqrt(nPbins);
467 //Int_t nPadY = nPadX;
468 //if ( nPadX * nPadY < nPbins ) nPadX++;
469 TCanvas* pdcaFitCan = 0x0;
475 for ( Int_t ibin=2; ibin<=nPbins; ++ibin ) {
476 currName = Form("hPDCA_%s_%s_%s", fSrcKeys->At(isrc)->GetName(), physSel.Data(), trigClassName.Data());
477 Int_t minBin = ( ibin == 0 ) ? 1 : ibin;
478 Int_t maxBin = ( ibin == 0 ) ? nPbins : ibin;
479 if ( ibin > 0 ) currName += Form("_pBin%i", ibin);
480 TH1* projHisto = histo->ProjectionY(currName.Data(), minBin, maxBin, "e");
481 if ( projHisto->GetEntries() == 0 ) continue;
482 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)));
483 if ( ipad % (nPadX*nPadY) == 0 ) {
484 currName = histo->GetName();
485 currName += Form("Fit_can_%i", ican++);
486 pdcaFitCan = new TCanvas(currName.Data(), currName.Data(),igroup1*xshift,igroup2*yshift,600,600);
487 pdcaFitCan->Divide(nPadX,nPadY);
490 pdcaFitCan->cd(++ipad);
491 if ( projHisto->Integral() > 0.) projHisto->Scale(1./projHisto->Integral());
493 if ( projHisto->GetEntries() > 50 ) {
494 fitFunc->SetParameter(0, projHisto->Integral());
495 fitFunc->SetParameter(1, projHisto->GetMean());
496 fitFunc->SetParameter(2, projHisto->GetRMS());
497 projHisto->Fit(fitFunc, "RQ", "e", xMinFit[itheta], xMaxFit[itheta]);
498 Double_t chi2 = fitFunc->GetChisquare();
499 Double_t ndf = fitFunc->GetNDF();
500 if ( ndf <= 0.) continue;
501 if ( chi2 / ndf > 100. ) continue;
502 Double_t sigma = TMath::Abs(fitFunc->GetParameter(2));
503 Double_t sigmaErr = fitFunc->GetParError(2);
504 sigmaVsP->SetBinContent(ibin, sigma);
505 sigmaVsP->SetBinError(ibin, sigmaErr);
506 //if ( ibin == 0 ) printf("%s %s sigma = %g +- %g\n", fThetaAbsKeys->At(itheta)->GetName(), fSrcKeys->At(isrc)->GetName(), sigma, sigmaErr);
508 else projHisto->Draw("e");
509 } // loop on momentum bins
511 sigmaVsP->SetLineColor(srcColors[isrc]);
512 sigmaVsP->SetMarkerColor(srcColors[isrc]);
513 sigmaVsP->SetMarkerStyle(20+isrc);
515 if ( gPad->GetListOfPrimitives()->GetEntries() != 0 ) drawOpt += "same";
516 sigmaVsP->Draw(drawOpt.Data());
517 sigmaVsP->Fit("pol0","Q",drawOpt.Data());
518 TF1* trendFit = (TF1*)sigmaVsP->GetListOfFunctions()->FindObject("pol0");
520 averageSigmaPdca[kNtrackSources*itheta+isrc] = trendFit->GetParameter(0);
521 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()));
523 leg->AddEntry(sigmaVsP, fSrcKeys->At(isrc)->GetName(), "lp");
525 // Plot 2D function for check!
526 histoPattern = GetHistoName(kPDCAVsPCheck, itheta, isrc);
527 TH2* histoCheck = (TH2*)GetSum(physSel, trigClassName, centralityRange, histoPattern);
528 if ( ! histoCheck ) histoCheck = histo; // needed for old data
529 currName = histoCheck->GetName();
530 currName.Append("_plotCut");
531 TCanvas* pdcaVsPcan = new TCanvas(currName.Data(), currName.Data(), igroup1*xshift,(igroup2+1)*yshift,600,600);
532 pdcaVsPcan->SetLogz();
533 pdcaVsPcan->SetRightMargin(0.12);
534 histoCheck->Draw("COLZ");
536 for ( Int_t icut=0; icut<nShowFuncs; ++icut ) {
537 currName = Form("%s_cutFunc%i", histo->GetName(), icut);
538 TF1* cutFunction = new TF1(currName.Data(),cutFormula.Data(), pMinCut, pMaxCut);
539 cutParam[icut][0] = sigmaMeasCut[itheta];
540 cutParam[icut][1] = nSigmaCut;
541 cutFunction->SetParameters(cutParam[icut]);
542 cutFunction->SetLineWidth(2);
543 cutFunction->SetLineColor(cutColor[icut]);
544 cutFunction->Draw("same");
545 } // loop on cut func
550 for ( Int_t icut=0; icut<nShowFuncs; ++icut ) {
551 currName = Form("sigmaCut_%s_%i", fThetaAbsKeys->At(itheta)->GetName(), icut);
552 TF1* cutFunction = new TF1(currName.Data(), cutFormula.Data(), pMinCut, pMaxCut);
553 cutParam[icut][0] = sigmaMeasCut[itheta];
554 cutParam[icut][1] = 1.;
555 cutFunction->SetParameters(cutParam[icut]);
556 cutFunction->SetLineColor(cutColor[icut]);
557 cutFunction->SetLineWidth(2);
558 cutFunction->Draw("same");
560 } // loop on theta abs
562 for ( Int_t isrc=0; isrc<kNtrackSources; ++isrc ) {
563 printf("muonCuts->SetSigmaPdca(%g, %g); // %s\n", averageSigmaPdca[isrc], averageSigmaPdca[kNtrackSources+isrc], fSrcKeys->At(isrc)->GetName());
568 for ( Int_t itheta=0; itheta<kNthetaAbs; ++itheta ) {
569 for ( Int_t isrc=0; isrc<kNtrackSources; ++isrc ) {
570 histoPattern = GetHistoName(kDCAVsPCheck, itheta, isrc);
571 TH2* histoCheck = (TH2*)GetSum(physSel, trigClassName, centralityRange, histoPattern);
572 if ( ! histoCheck ) continue;
573 currName = histoCheck->GetName();
574 currName.Append("_plotCut");
575 TCanvas* pdcaVsPcan = new TCanvas(currName.Data(), currName.Data(), igroup1*xshift,(igroup2+1)*yshift,600,600);
576 pdcaVsPcan->SetRightMargin(0.12);
577 pdcaVsPcan->SetLogz();
578 histoCheck->Draw("COLZ");
580 for ( Int_t icut=0; icut<nShowFuncs; ++icut ) {
581 currName = histoCheck->GetName();
582 currName.Append(Form("_cutFunc%i",icut));
583 TString currFormula = cutFormula;
584 currFormula.Append("/x");
585 TF1* cutFunction = new TF1(currName.Data(),currFormula.Data(), pMinCut, pMaxCut);
586 cutParam[icut][0] = sigmaMeasCut[itheta];
587 cutParam[icut][1] = nSigmaCut;
588 cutFunction->SetParameters(cutParam[icut]);
589 cutFunction->SetLineWidth(2);
590 cutFunction->SetLineColor(cutColor[icut]);
591 cutFunction->Draw("same");
592 } // loop on cut functions
597 ///////////////////////////
598 // Check Chi square prob //
599 ///////////////////////////
602 for ( Int_t itheta=0; itheta<kNthetaAbs; ++itheta ) {
603 for ( Int_t isrc=0; isrc<kNtrackSources; ++isrc ) {
604 histoPattern = GetHistoName(kChiProbVsP, itheta, isrc);
605 TH2* histo = (TH2*)GetSum(physSel, trigClassName, centralityRange, histoPattern);
606 if ( ! histo ) continue;
608 Int_t nPbins = histo->GetXaxis()->GetNbins();
609 //Int_t nPadX = (Int_t)TMath::Sqrt(nPbins);
610 //Int_t nPadY = nPadX;
611 //if ( nPadX * nPadY < nPbins ) nPadX++;
612 TCanvas* chi2probCan = 0x0;
618 for ( Int_t ibin=0; ibin<=nPbins; ++ibin ) {
619 currName = Form("hChiProb_%s_%s_%s", fSrcKeys->At(isrc)->GetName(), physSel.Data(), trigClassName.Data());
620 Int_t minBin = ( ibin == 0 ) ? 1 : ibin;
621 Int_t maxBin = ( ibin == 0 ) ? nPbins : ibin;
622 if ( ibin > 0 ) currName += Form("_pBin%i", ibin);
623 TH1* projHisto = histo->ProjectionY(currName.Data(), minBin, maxBin);
624 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)));
625 if ( projHisto->GetEntries() == 0 ) continue;
626 if ( ipad % (nPadX*nPadY) == 0 ) {
627 currName = histo->GetName();
628 currName += Form("Fit_can_%i", ican++);
629 chi2probCan = new TCanvas(currName.Data(), currName.Data(),igroup1*xshift,igroup2*yshift,600,600);
630 chi2probCan->Divide(nPadX,nPadY);
633 chi2probCan->cd(++ipad);
635 projHisto->Draw("e");
636 } // loop on momentum bins
638 } // loop on theta abs
641 //////////////////////
642 // Check sigma cuts //
643 //////////////////////
644 printf("\nReference sigma cut %g\n", refSigmaCut);
646 Float_t fracOfHeight = 0.35;
647 Float_t rightMargin = 0.03;
648 Int_t cutColors[14] = {kBlack, kRed, kBlue, kGreen, kCyan, kMagenta, kOrange, kViolet, kSpring, kGray, kSpring, kAzure, kPink, kYellow};
649 Int_t* orderCuts = 0x0;
650 Int_t nSigmaCuts = 0;
651 Int_t checkHistos[2] = {kSigmaVsPt, kSigmaVsEta};
652 Bool_t useCustomSigma = furtherOpt.Contains("CUSTOMSIGMA");
657 for ( Int_t ihisto=0; ihisto<2; ++ihisto ) {
658 for ( Int_t isrc=0; isrc<kNtrackSources; ++isrc ) {
660 for ( Int_t itheta=0; itheta<kNthetaAbs; ++itheta ) {
661 histoPattern = GetHistoName(checkHistos[ihisto], itheta, isrc);
662 TH2* histo = (TH2*)GetSum(physSel, trigClassName, centralityRange, histoPattern);
663 if ( ! histo ) continue;
664 if ( histo->Integral() == 0. ) {
670 TAxis* axis = histo->GetYaxis();
671 nSigmaCuts = ( useCustomSigma ) ? fSigmaCuts.GetSize() : axis->GetNbins();
672 orderCuts = new Int_t[nSigmaCuts];
674 for ( Int_t isigma=0; isigma<axis->GetNbins(); ++isigma ) {
675 TString currLabel = axis->GetBinLabel(isigma+1);
676 Double_t currSigma = currLabel.Atof();
677 if ( useCustomSigma ) {
678 Bool_t foundMatch = kFALSE;
679 for ( Int_t jsigma=0; jsigma<fSigmaCuts.GetSize(); ++jsigma ) {
680 if ( TMath::Abs(currSigma - fSigmaCuts[jsigma]) < 1e-4 ) {
685 if ( ! foundMatch ) continue;
687 Int_t currBin = ( TMath::Abs(currSigma - refSigmaCut) < 1e-4) ? 0 : ++countBin;
688 orderCuts[currBin] = isigma+1;
691 currName = histo->GetName();
692 currName.Append("_can");
693 can = new TCanvas(currName.Data(), currName.Data(), igroup1*xshift,igroup2*yshift, 600, 600);
694 TLegend* leg = new TLegend(0.6,0.6,0.9,0.9);
695 leg->SetBorderSize(1);
696 can->Divide(1,2,0,0);
698 gPad->SetPad(0., fracOfHeight, 0.99, 0.99);
699 //gPad->SetTopMargin(0.08);
700 gPad->SetTopMargin(0.03);
701 gPad->SetRightMargin(rightMargin);
705 gPad->SetPad(0., 0., 0.99, fracOfHeight);
706 gPad->SetRightMargin(rightMargin);
707 gPad->SetBottomMargin(0.08/fracOfHeight);
709 TH1* refCutHisto = 0x0;
710 TString legTitle = "";
711 for ( Int_t isigma=0; isigma<nSigmaCuts; ++isigma ) {
712 currName = histo->GetName();
713 currName.Append(Form("_sigma%i", isigma));
714 Int_t currBin = orderCuts[isigma];
715 TH1* projectHisto = histo->ProjectionX(currName.Data(), currBin, currBin);
716 projectHisto->SetLineColor(cutColors[isigma]);
717 projectHisto->SetMarkerColor(cutColors[isigma]);
718 projectHisto->SetMarkerStyle(20+isigma);
719 projectHisto->SetStats(kFALSE);
720 projectHisto->SetTitle("");
723 if ( gPad->GetListOfPrimitives()->GetEntries() != 0 ) drawOpt += "same";
724 projectHisto->Draw(drawOpt.Data());
725 TString currLabel = histo->GetYaxis()->GetBinLabel(currBin);
726 TString cutName = ( TMath::Abs(currLabel.Atof() ) > 1000. ) ? "Total" : Form("%s #sigma cut", currLabel.Data());
727 leg->AddEntry(projectHisto, cutName.Data(), "lp");
728 if ( ! refCutHisto ) {
729 refCutHisto = projectHisto;
730 legTitle = Form("(pass n #sigma) / (pass %s #sigma)", currLabel.Data());
733 currName.Append("_ratio");
734 TH1* histoRatio = (TH1*)projectHisto->Clone(currName.Data());
736 histoRatio->Divide(refCutHisto);
737 histoRatio->SetTitle("");
738 histoRatio->GetYaxis()->SetTitle(legTitle.Data());
739 histoRatio->GetXaxis()->SetLabelSize(0.04/fracOfHeight);
740 histoRatio->GetXaxis()->SetTitleSize(0.035/fracOfHeight);
741 histoRatio->GetYaxis()->SetLabelSize(0.03/fracOfHeight);
742 histoRatio->GetYaxis()->SetTitleSize(0.03/fracOfHeight);
743 histoRatio->GetXaxis()->SetTitleOffset(1.);
744 histoRatio->GetYaxis()->SetTitleOffset(0.6);
745 histoRatio->GetYaxis()->SetRangeUser(0.5,1.5);
749 if ( gPad->GetListOfPrimitives()->GetEntries() != 0 ) drawOpt += "same";
750 histoRatio->Draw(drawOpt.Data());
751 }// loop on sigma cuts
754 } // loop on theta abs
756 } // loop on histo type
760 Double_t ptMin[] = {0., 2., 4., 15.};
761 Int_t nPtMins = sizeof(ptMin)/sizeof(ptMin[0]);
762 for ( Int_t iptmin=0; iptmin<nPtMins; ++iptmin) {
763 for ( Int_t isrc=0; isrc<kNtrackSources; ++isrc ) {
765 for ( Int_t itheta=0; itheta<kNthetaAbs; ++itheta ) {
767 histoName = GetHistoName(kSigmaVsPt, itheta, isrc);
768 for ( Int_t icent=1; icent<=fCentralityClasses->GetNbins(); ++icent ) {
769 TH2* histo = (TH2*)GetSum(physSel, trigClassName, fCentralityClasses->GetBinLabel(icent), histoName);
770 if ( ! histo ) continue;
771 Int_t ptMinBin = histo->GetXaxis()->FindBin(ptMin[iptmin]);
772 Int_t ptMaxBin = histo->GetXaxis()->GetNbins()+1;
773 currName = histo->GetName();
774 currName += Form("_%s_ptMin%i", fCentralityClasses->GetBinLabel(icent), TMath::Nint(ptMin[iptmin]));
775 TH1* projectHisto = histo->ProjectionY(currName.Data(), ptMinBin, ptMaxBin);
776 if ( projectHisto->GetMaximum() < 50. ) {
781 currName = Form("CutEffect_%s_ptMin%i", fSrcKeys->At(isrc)->GetName(), TMath::Nint(ptMin[iptmin]));
782 can = new TCanvas(currName.Data(), currName.Data(), igroup1*xshift,igroup2*yshift, 600, 600);
786 leg = new TLegend(0.6,0.6,0.9,0.9);
787 leg->SetBorderSize(1);
790 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)));
791 projectHisto->SetLineColor(cutColors[icent-1]);
792 projectHisto->SetMarkerColor(cutColors[icent-1]);
793 projectHisto->SetMarkerStyle(19+icent);
794 projectHisto->SetStats(0);
796 if ( gPad->GetListOfPrimitives()->GetEntries() != 0 ) drawOpt += "same";
797 projectHisto->Draw(drawOpt.Data());
798 leg->AddEntry(projectHisto, fCentralityClasses->GetBinLabel(icent), "lp");
799 Double_t keptEvents = projectHisto->GetBinContent(orderCuts[0]);
800 Double_t totalEvents = projectHisto->GetBinContent(orderCuts[nSigmaCuts-1]);
801 Double_t accepted = ( totalEvents > 0. ) ? keptEvents / totalEvents : 1.;
802 Double_t acceptedErr = ( totalEvents > 0. ) ? TMath::Sqrt(accepted*(1.-accepted)/totalEvents) : 1.;
803 printf("%12s %11s %6s (pt>%g) rejected evts : %6.2f +- %6.3f %%\n", fSrcKeys->At(isrc)->GetName(), fThetaAbsKeys->At(itheta)->GetName(), fCentralityClasses->GetBinLabel(icent), ptMin[iptmin], (1.-accepted)*100., acceptedErr*100.);
804 //printf(" rejected %g total %g (%s vs %s)\n",totalEvents-keptEvents,totalEvents,projectHisto->GetXaxis()->GetBinLabel(orderCuts[0]),projectHisto->GetXaxis()->GetBinLabel(orderCuts[nSigmaCuts-1]));
805 } // loop on centrality
806 if ( leg ) leg->Draw("same");
807 } // loop on theta abs