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"
47 //#include "TMCProcess.h"
50 #include "AliAODEvent.h"
51 #include "AliAODTrack.h"
52 #include "AliAODMCParticle.h"
53 #include "AliMCEvent.h"
54 #include "AliMCParticle.h"
55 //#include "AliStack.h"
56 #include "AliESDEvent.h"
57 #include "AliESDMuonTrack.h"
58 #include "AliInputEventHandler.h"
60 #include "AliAnalysisManager.h"
63 #include "AliMergeableCollection.h"
64 #include "AliMuonTrackCuts.h"
68 ClassImp(AliAnalysisTaskMuonCuts) // Class implementation in ROOT context
72 //________________________________________________________________________
73 AliAnalysisTaskMuonCuts::AliAnalysisTaskMuonCuts() :
82 //________________________________________________________________________
83 AliAnalysisTaskMuonCuts::AliAnalysisTaskMuonCuts(const char *name, const AliMuonTrackCuts& cuts) :
84 AliVAnalysisMuon(name, cuts),
92 TString histoTypeKeys = "hDCAxVsP hDCAyVsP hPdcaVsP hPDCAVsPCheck hDCAVsPCheck hChiProbVsP hSigmaVsPt hSigmaVsEta";
93 fHistoTypeKeys = histoTypeKeys.Tokenize(" ");
95 TString thetaAbsKeys = "ThetaAbs23 ThetaAbs310";
96 fThetaAbsKeys = thetaAbsKeys.Tokenize(" ");
102 //________________________________________________________________________
103 AliAnalysisTaskMuonCuts::~AliAnalysisTaskMuonCuts()
109 delete fHistoTypeKeys;
110 delete fThetaAbsKeys;
113 //___________________________________________________________________________
114 void AliAnalysisTaskMuonCuts::MyUserCreateOutputObjects()
117 TString histoName = "";
118 for ( Int_t itheta=0; itheta<kNthetaAbs; ++itheta ) {
119 for ( Int_t isrc=0; isrc<kNtrackSources; ++isrc ) {
120 histoName = GetHistoName(kDCAxVsP, itheta, isrc);
121 histo = new TH2F(histoName.Data(), histoName.Data(), 50, 0., 200., 100, -50., 50.);
122 histo->SetXTitle("p (GeV/c)");
123 histo->SetYTitle("DCA_{x} (cm)");
124 AddObjectToCollection(histo);
126 histoName = GetHistoName(kDCAyVsP, itheta, isrc);
127 histo = new TH2F(histoName.Data(), histoName.Data(), 50, 0., 200., 100, -50., 50.);
128 histo->SetXTitle("p (GeV/c)");
129 histo->SetYTitle("DCA_{y} (cm)");
130 AddObjectToCollection(histo);
132 histoName = GetHistoName(kPdcaVsP, itheta, isrc);
133 histo = new TH2F(histoName.Data(), histoName.Data(), 100, 0., 400., 50, 0., 500.);
134 histo->SetXTitle("p (GeV/c)");
135 histo->SetYTitle("p#timesDCA (cm #times GeV/c)");
136 AddObjectToCollection(histo);
138 histoName = GetHistoName(kPDCAVsPCheck, itheta, isrc);
139 histo = new TH2F(histoName.Data(), histoName.Data(), 100, 0., 800., 100, 0., 5000.);
140 histo->SetXTitle("p (GeV/c)");
141 histo->SetYTitle("p#timesDCA (cm #times GeV/c)");
142 AddObjectToCollection(histo);
144 histoName = GetHistoName(kDCAVsPCheck, itheta, isrc);
145 histo = new TH2F(histoName.Data(), histoName.Data(), 100, 0., 800., 100, 0., 200.);
146 histo->SetXTitle("p (GeV/c)");
147 histo->SetYTitle("DCA (cm)");
148 AddObjectToCollection(histo);
150 histoName = GetHistoName(kChiProbVsP, itheta, isrc);
151 histo = new TH2F(histoName.Data(), histoName.Data(), 50, 0., 200., 50, 0., 1.);
152 histo->SetXTitle("p (GeV/c)");
153 histo->SetYTitle("Chisquare prob.");
154 AddObjectToCollection(histo);
156 histoName = GetHistoName(kSigmaVsPt, itheta, isrc);
157 histo = new TH2F(histoName.Data(), histoName.Data(), 100, 0., 100., fSigmaCuts.GetSize(), 0.5, 0.5+(Double_t)fSigmaCuts.GetSize());
158 histo->SetXTitle("p_{t} (GeV/c)");
159 histo->SetYTitle("#sigma_{p#timesDCA}");
160 for ( Int_t ibin=0; ibin<fSigmaCuts.GetSize(); ++ibin ) {
161 histo->GetYaxis()->SetBinLabel(ibin+1,Form("%g", fSigmaCuts[ibin]));
163 AddObjectToCollection(histo);
165 histoName = GetHistoName(kSigmaVsEta, itheta, isrc);
166 histo = new TH2F(histoName.Data(), histoName.Data(), 25, -4.5, -2., fSigmaCuts.GetSize(), 0.5, 0.5+(Double_t)fSigmaCuts.GetSize());
167 histo->SetXTitle("#eta");
168 histo->SetYTitle("#sigma_{p#timesDCA}");
169 for ( Int_t ibin=0; ibin<fSigmaCuts.GetSize(); ++ibin ) {
170 histo->GetYaxis()->SetBinLabel(ibin+1,Form("%g", fSigmaCuts[ibin]));
172 AddObjectToCollection(histo);
173 } // loop on track sources
174 } // loop on theta abs
176 fMuonTrackCuts->Print();
180 //________________________________________________________________________
181 TString AliAnalysisTaskMuonCuts::GetHistoName(Int_t histoTypeIndex, Int_t thetaAbsIndex, Int_t srcIndex)
183 /// Get local histogram name
184 TString histoName = Form("%s%s%s", fHistoTypeKeys->At(histoTypeIndex)->GetName(), fThetaAbsKeys->At(thetaAbsIndex)->GetName(), fSrcKeys->At(srcIndex)->GetName());
189 //________________________________________________________________________
190 void AliAnalysisTaskMuonCuts::ProcessEvent(TString physSel, const TObjArray& selectTrigClasses, TString centrality)
196 if ( GetVertexSPD()->GetNContributors() < fMinNvtxContirbutors ) return;
198 TString histoName = "";
199 AliVParticle* track = 0x0;
200 Int_t nTracks = GetNTracks();
201 for (Int_t itrack = 0; itrack < nTracks; itrack++) {
202 track = GetTrack(itrack);
203 fMuonTrackCuts->SetNSigmaPdca(1.e10);
204 if ( ! fMuonTrackCuts->IsSelected(track) ) continue;
206 Double_t rAbsEnd = ( fAODEvent ) ? ((AliAODTrack*)track)->GetRAtAbsorberEnd(): ((AliESDMuonTrack*)track)->GetRAtAbsorberEnd();
207 Double_t thetaAbsEndDeg = TMath::ATan( rAbsEnd / 505. ) * TMath::RadToDeg();
208 Int_t thetaAbsBin = ( thetaAbsEndDeg < 3. ) ? kThetaAbs23 : kThetaAbs310;
210 Int_t trackSrc = GetParticleType(track);
212 TVector3 dcaAtVz = fMuonTrackCuts->GetCorrectedDCA(track);
213 Double_t pTotMean = fMuonTrackCuts->GetAverageMomentum(track);
214 Double_t dca = dcaAtVz.Mag();
215 Double_t pDca = pTotMean * dca;
217 Double_t chi2 = pDca / fMuonTrackCuts->GetSigmaPdca(rAbsEnd) ;
219 Double_t chiProb = TMath::Prob(chi2, 2);
221 Double_t pTot = track->P();
222 Double_t pt = track->Pt();
223 Double_t eta = track->Eta();
225 for ( Int_t itrig=0; itrig<selectTrigClasses.GetEntries(); ++itrig ) {
226 TString trigClassName = ((TObjString*)selectTrigClasses.At(itrig))->GetString();
227 histoName = GetHistoName(kDCAxVsP, thetaAbsBin, trackSrc);
228 ((TH2*)GetMergeableObject(physSel, trigClassName, centrality, histoName))->Fill(pTot, dcaAtVz.X());
230 histoName = GetHistoName(kDCAyVsP, thetaAbsBin, trackSrc);
231 ((TH2*)GetMergeableObject(physSel, trigClassName, centrality, histoName))->Fill(pTot, dcaAtVz.Y());
233 histoName = GetHistoName(kPdcaVsP, thetaAbsBin, trackSrc);
234 ((TH2*)GetMergeableObject(physSel, trigClassName, centrality, histoName))->Fill(pTot, pDca);
236 histoName = GetHistoName(kPDCAVsPCheck, thetaAbsBin, trackSrc);
237 ((TH2*)GetMergeableObject(physSel, trigClassName, centrality, histoName))->Fill(pTot, pDca);
239 histoName = GetHistoName(kDCAVsPCheck, thetaAbsBin, trackSrc);
240 ((TH2*)GetMergeableObject(physSel, trigClassName, centrality, histoName))->Fill(pTot, dca);
242 histoName = GetHistoName(kChiProbVsP, thetaAbsBin, trackSrc);
243 ((TH2*)GetMergeableObject(physSel, trigClassName, centrality, histoName))->Fill(pTot, chiProb);
244 } // loop on selected trigger classes
246 for ( Int_t isigma=0; isigma<fSigmaCuts.GetSize(); ++isigma) {
247 fMuonTrackCuts->SetNSigmaPdca(fSigmaCuts[isigma]);
248 if ( ! fMuonTrackCuts->IsSelected(track) ) continue;
249 for ( Int_t itrig=0; itrig<selectTrigClasses.GetEntries(); ++itrig ) {
250 TString trigClassName = ((TObjString*)selectTrigClasses.At(itrig))->GetString();
251 histoName = GetHistoName(kSigmaVsPt, thetaAbsBin, trackSrc);
252 ((TH2*)GetMergeableObject(physSel, trigClassName, centrality, histoName))->Fill(pt, isigma+1);
253 histoName = GetHistoName(kSigmaVsEta, thetaAbsBin, trackSrc);
254 ((TH2*)GetMergeableObject(physSel, trigClassName, centrality, histoName))->Fill(eta, isigma+1);
255 } // loop on selected trigger classes
260 //________________________________________________________________________
261 void AliAnalysisTaskMuonCuts::SetSigmaCuts(Int_t nSigmaCuts, Double_t* sigmaCuts)
263 /// Set number of sigmas
264 if ( ! sigmaCuts || nSigmaCuts < 0 ) {
265 // if ( defaultChiSquare ) {
266 // 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.};
267 // Int_t nCuts = sizeof(cuts)/sizeof(cuts[0]);
268 // fSigmaCuts.Set(nCuts, cuts);
271 Double_t cuts[] = {2., 3., 4., 5., 6., 7., 10., 15., 20., 25., 30., 1.e10};
272 Int_t nCuts = sizeof(cuts)/sizeof(cuts[0]);
273 fSigmaCuts.Set(nCuts, cuts);
277 fSigmaCuts.Set(nSigmaCuts, sigmaCuts);
281 //________________________________________________________________________
282 void AliAnalysisTaskMuonCuts::Terminate(Option_t *) {
284 /// Draw some histogram at the end.
287 AliVAnalysisMuon::Terminate("");
289 if ( ! fMergeableCollection ) return;
291 TString physSel = fTerminateOptions->At(0)->GetName();
292 TString trigClassName = fTerminateOptions->At(1)->GetName();
293 TString centralityRange = fTerminateOptions->At(2)->GetName();
294 TString furtherOpt = fTerminateOptions->At(3)->GetName();
295 furtherOpt.ToUpper();
297 furtherOpt.ReplaceAll(" ", " ");
298 furtherOpt.ReplaceAll(" =", "=");
299 furtherOpt.ReplaceAll("= ", "=");
300 TObjArray* optArray = furtherOpt.Tokenize(" ");
301 Double_t refSigmaCut = 6.;
302 for ( Int_t iopt=0; iopt<optArray->GetEntries(); ++iopt ) {
303 TString currOpt = optArray->At(iopt)->GetName();
304 if ( currOpt.Contains("REFSIGMA") ) {
305 currOpt.Remove(0,currOpt.Index("=")+1);
306 refSigmaCut = currOpt.Atof();
311 Int_t srcColors[kNtrackSources] = {kBlack, kRed, kGreen, kBlue, kViolet, 7, kOrange};
324 TString histoName = "", currName = "", histoPattern = "", drawOpt = "";
325 currName = Form("%s_recoDCA", GetName());
326 can = new TCanvas(currName.Data(),"Reco DCA",igroup1*xshift,igroup2*yshift,600,600);
329 Int_t recoDcaHisto[2] = {kDCAxVsP, kDCAyVsP};
330 TString dcaName[2] = {"DCAx", "DCAy"};
331 Double_t averageDca[4] = {0., 0.};
332 printf("\nAverage reconstructed DCA:\n");
333 TF1* fitFuncMeanDca = new TF1("fitFuncMeanDca","gausn",-20.,20.);
334 fitFuncMeanDca->SetParNames("Norm", "Mean", "Sigma");
335 for ( Int_t itheta=0; itheta<kNthetaAbs; ++itheta ) {
336 for ( Int_t ihisto=0; ihisto<2; ++ihisto ) {
339 histoPattern = Form("%s & %s", fHistoTypeKeys->At(recoDcaHisto[ihisto])->GetName(), fThetaAbsKeys->At(itheta)->GetName());
340 histo = (TH2*)GetSum(physSel, trigClassName, centralityRange, histoPattern);
341 if ( ! histo ) continue;
343 TH1* meanDcaVsP = histo->ProjectionX(Form("mean%sVsP_%s", dcaName[ihisto].Data(), fThetaAbsKeys->At(itheta)->GetName()));
345 meanDcaVsP->SetTitle(Form("Mean %s vs. p %s", dcaName[ihisto].Data(), fThetaAbsKeys->At(itheta)->GetName()));
346 meanDcaVsP->SetYTitle(Form("<%s> (cm)", dcaName[ihisto].Data()));
347 meanDcaVsP->SetStats(kFALSE);
349 Int_t nPbins = histo->GetXaxis()->GetNbins();
350 //Int_t nPadX = (Int_t)TMath::Sqrt(nPbins);
351 //Int_t nPadY = nPadX;
352 //if ( nPadX * nPadY < nPbins ) nPadX++;
353 TCanvas* meanDcaFitCan = 0x0;
359 for ( Int_t ibin=2; ibin<=nPbins; ++ibin ) {
360 currName = Form("hMean%s_%s_%s", dcaName[ihisto].Data(), physSel.Data(), trigClassName.Data());
361 Int_t minBin = ( ibin == 0 ) ? 1 : ibin;
362 Int_t maxBin = ( ibin == 0 ) ? nPbins : ibin;
363 if ( ibin > 0 ) currName += Form("_pBin%i", ibin);
364 TH1* projHisto = histo->ProjectionY(currName.Data(), minBin, maxBin, "e");
365 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)));
367 if ( projHisto->GetEntries() == 0 ) continue;
368 if ( ipad % (nPadX*nPadY) == 0 ) {
369 currName = histo->GetName();
370 currName += Form("Fit_can_%i", ican++);
371 meanDcaFitCan = new TCanvas(currName.Data(), currName.Data(),igroup1*xshift,igroup2*yshift,600,600);
372 meanDcaFitCan->Divide(nPadX,nPadY);
375 meanDcaFitCan->cd(++ipad);
377 if ( projHisto->Integral() > 50 ) {
378 fitFuncMeanDca->SetParameter(0, projHisto->Integral());
379 fitFuncMeanDca->SetParameter(1, projHisto->GetMean());
380 fitFuncMeanDca->SetParameter(2, projHisto->GetRMS());
381 Double_t fitDcaLim = ( ibin == 0 ) ? 40. : 40./((Double_t)ibin);
382 fitDcaLim = TMath::Max(5., fitDcaLim);
383 projHisto->Fit(fitFuncMeanDca, "RQ", "e", -fitDcaLim, fitDcaLim);
384 Double_t chi2 = fitFuncMeanDca->GetChisquare();
385 Double_t ndf = fitFuncMeanDca->GetNDF();
386 if ( ndf <= 0.) continue;
387 if ( chi2 / ndf > 100. ) continue;
388 Double_t meanDca = fitFuncMeanDca->GetParameter(1);
389 Double_t meanDcaErr = fitFuncMeanDca->GetParError(1);
390 meanDcaVsP->SetBinContent(ibin, meanDca);
391 meanDcaVsP->SetBinError(ibin, meanDcaErr);
392 //if ( ibin == 0 ) printf("%s meanDca = %g +- %g\n", fThetaAbsKeys->At(itheta)->GetName(), meanDca, meanDcaErr);
394 else projHisto->Draw("e");
395 } // loop on momentum bins
396 can->cd(2*itheta+ihisto+1);
397 //meanDcaVsP->SetLineColor(srcColors[isrc]);
398 //meanDcaVsP->SetMarkerColor(srcColors[isrc]);
399 //meanDcaVsP->SetMarkerStyle(20+isrc);
400 meanDcaVsP->Fit("pol0","Q");
401 TF1* trendFit = (TF1*)meanDcaVsP->GetListOfFunctions()->FindObject("pol0");
403 averageDca[2*itheta+ihisto] = trendFit->GetParameter(0);
404 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()));
407 //leg->AddEntry(meanDcaVsP, fSrcKeys->At(isrc)->GetName(), "lp");
414 //meanDca[ihisto] = histo->GetMean();
415 } // loop on histo type
416 } // loop on theta abs
418 for ( Int_t itheta=0; itheta<kNthetaAbs; ++itheta ) {
419 printf("muonCuts->SetMeanDCA(%g, %g, 0.); // %s\n", averageDca[2*itheta], averageDca[2*itheta+1], fThetaAbsKeys->At(itheta)->GetName());
425 Double_t nSigmaCut = fMuonTrackCuts->GetNSigmaPdca(); //6.;
426 Double_t sigmaMeasCut[2] = { fMuonTrackCuts->GetSigmaPdca(AliMuonTrackCuts::kThetaAbs23), fMuonTrackCuts->GetSigmaPdca(AliMuonTrackCuts::kThetaAbs23)}; //{99., 54.}; //{120., 63.};
427 Double_t relPResolution = fMuonTrackCuts->GetRelPResolution(); //4.5e-4; //6.e-3;//8.e-4;
428 Double_t angleResolution = 535.*fMuonTrackCuts->GetSlopeResolution(); //6.e-4;
429 Double_t pMinCut = 0.1;
430 Double_t pMaxCut = 800.;
431 const Int_t kNCutFuncs = 2;
432 Int_t nShowFuncs = 1;
433 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)";
434 Double_t cutParam[kNCutFuncs][4] = {{sigmaMeasCut[0], nSigmaCut, relPResolution, angleResolution}, {sigmaMeasCut[0], nSigmaCut, 0., 0.32}};
435 Int_t cutColor[kNCutFuncs] = {kBlack, kRed};
438 can = new TCanvas("pdcaSigmaFit","Sigma vs P fit",igroup1*xshift,igroup2*yshift,600,600);
441 TF1* fitFunc = new TF1("fitFunc", "x*gausn", 0., 400.);
442 fitFunc->SetParNames("Norm", "Mean", "Sigma");
443 gStyle->SetOptFit(1111);
444 Double_t xMinFit[2] = {0., 0.};
445 Double_t xMaxFit[2] = {320., 150.}; // {360., 180.};
446 printf("\nSigma p x DCA:\n");
447 Double_t averageSigmaPdca[kNtrackSources*kNthetaAbs] = {0.};
448 for ( Int_t itheta=0; itheta<kNthetaAbs; ++itheta ) {
449 TLegend* leg = new TLegend(0.7, 0.7, 0.9, 0.9);
450 leg->SetBorderSize(1);
451 for ( Int_t isrc=0; isrc<kNtrackSources; ++isrc ) {
452 histoPattern = GetHistoName(kPdcaVsP, itheta, isrc);
453 TH2* histo = (TH2*)GetSum(physSel, trigClassName, centralityRange, histoPattern);
454 if ( ! histo ) continue;
455 if ( histo->Integral() < 200 ) {
460 TH1* sigmaVsP = histo->ProjectionX(Form("sigma%s_%s_%s", fHistoTypeKeys->At(kPdcaVsP)->GetName(), fThetaAbsKeys->At(itheta)->GetName(), fSrcKeys->At(isrc)->GetName()));
462 sigmaVsP->SetTitle(Form("#sigma_{p#timesDCA} vs. p %s %s", fThetaAbsKeys->At(itheta)->GetName(), fSrcKeys->At(isrc)->GetName()));
463 sigmaVsP->SetYTitle("#sigma_{p#timesDCA} (cm #times GeV/c)");
464 sigmaVsP->SetStats(kFALSE);
466 Int_t nPbins = histo->GetXaxis()->GetNbins();
467 //Int_t nPadX = (Int_t)TMath::Sqrt(nPbins);
468 //Int_t nPadY = nPadX;
469 //if ( nPadX * nPadY < nPbins ) nPadX++;
470 TCanvas* pdcaFitCan = 0x0;
476 for ( Int_t ibin=2; ibin<=nPbins; ++ibin ) {
477 currName = Form("hPDCA_%s_%s_%s", fSrcKeys->At(isrc)->GetName(), physSel.Data(), trigClassName.Data());
478 Int_t minBin = ( ibin == 0 ) ? 1 : ibin;
479 Int_t maxBin = ( ibin == 0 ) ? nPbins : ibin;
480 if ( ibin > 0 ) currName += Form("_pBin%i", ibin);
481 TH1* projHisto = histo->ProjectionY(currName.Data(), minBin, maxBin, "e");
482 if ( projHisto->GetEntries() == 0 ) continue;
483 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)));
484 if ( ipad % (nPadX*nPadY) == 0 ) {
485 currName = histo->GetName();
486 currName += Form("Fit_can_%i", ican++);
487 pdcaFitCan = new TCanvas(currName.Data(), currName.Data(),igroup1*xshift,igroup2*yshift,600,600);
488 pdcaFitCan->Divide(nPadX,nPadY);
491 pdcaFitCan->cd(++ipad);
492 if ( projHisto->Integral() > 0.) projHisto->Scale(1./projHisto->Integral());
494 if ( projHisto->GetEntries() > 50 ) {
495 fitFunc->SetParameter(0, projHisto->Integral());
496 fitFunc->SetParameter(1, projHisto->GetMean());
497 fitFunc->SetParameter(2, projHisto->GetRMS());
498 projHisto->Fit(fitFunc, "RQ", "e", xMinFit[itheta], xMaxFit[itheta]);
499 Double_t chi2 = fitFunc->GetChisquare();
500 Double_t ndf = fitFunc->GetNDF();
501 if ( ndf <= 0.) continue;
502 if ( chi2 / ndf > 100. ) continue;
503 Double_t sigma = TMath::Abs(fitFunc->GetParameter(2));
504 Double_t sigmaErr = fitFunc->GetParError(2);
505 sigmaVsP->SetBinContent(ibin, sigma);
506 sigmaVsP->SetBinError(ibin, sigmaErr);
507 //if ( ibin == 0 ) printf("%s %s sigma = %g +- %g\n", fThetaAbsKeys->At(itheta)->GetName(), fSrcKeys->At(isrc)->GetName(), sigma, sigmaErr);
509 else projHisto->Draw("e");
510 } // loop on momentum bins
512 sigmaVsP->SetLineColor(srcColors[isrc]);
513 sigmaVsP->SetMarkerColor(srcColors[isrc]);
514 sigmaVsP->SetMarkerStyle(20+isrc);
516 if ( gPad->GetListOfPrimitives()->GetEntries() != 0 ) drawOpt += "same";
517 sigmaVsP->Draw(drawOpt.Data());
518 sigmaVsP->Fit("pol0","Q",drawOpt.Data());
519 TF1* trendFit = (TF1*)sigmaVsP->GetListOfFunctions()->FindObject("pol0");
521 averageSigmaPdca[kNtrackSources*itheta+isrc] = trendFit->GetParameter(0);
522 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()));
524 leg->AddEntry(sigmaVsP, fSrcKeys->At(isrc)->GetName(), "lp");
526 // Plot 2D function for check!
527 histoPattern = GetHistoName(kPDCAVsPCheck, itheta, isrc);
528 TH2* histoCheck = (TH2*)GetSum(physSel, trigClassName, centralityRange, histoPattern);
529 if ( ! histoCheck ) histoCheck = histo; // needed for old data
530 currName = histoCheck->GetName();
531 currName.Append("_plotCut");
532 TCanvas* pdcaVsPcan = new TCanvas(currName.Data(), currName.Data(), igroup1*xshift,(igroup2+1)*yshift,600,600);
533 pdcaVsPcan->SetLogz();
534 pdcaVsPcan->SetRightMargin(0.12);
535 histoCheck->Draw("COLZ");
537 for ( Int_t icut=0; icut<nShowFuncs; ++icut ) {
538 currName = Form("%s_cutFunc%i", histo->GetName(), icut);
539 TF1* cutFunction = new TF1(currName.Data(),cutFormula.Data(), pMinCut, pMaxCut);
540 cutParam[icut][0] = sigmaMeasCut[itheta];
541 cutParam[icut][1] = nSigmaCut;
542 cutFunction->SetParameters(cutParam[icut]);
543 cutFunction->SetLineWidth(2);
544 cutFunction->SetLineColor(cutColor[icut]);
545 cutFunction->Draw("same");
546 } // loop on cut func
551 for ( Int_t icut=0; icut<nShowFuncs; ++icut ) {
552 currName = Form("sigmaCut_%s_%i", fThetaAbsKeys->At(itheta)->GetName(), icut);
553 TF1* cutFunction = new TF1(currName.Data(), cutFormula.Data(), pMinCut, pMaxCut);
554 cutParam[icut][0] = sigmaMeasCut[itheta];
555 cutParam[icut][1] = 1.;
556 cutFunction->SetParameters(cutParam[icut]);
557 cutFunction->SetLineColor(cutColor[icut]);
558 cutFunction->SetLineWidth(2);
559 cutFunction->Draw("same");
561 } // loop on theta abs
563 for ( Int_t isrc=0; isrc<kNtrackSources; ++isrc ) {
564 printf("muonCuts->SetSigmaPdca(%g, %g); // %s\n", averageSigmaPdca[isrc], averageSigmaPdca[kNtrackSources+isrc], fSrcKeys->At(isrc)->GetName());
569 for ( Int_t itheta=0; itheta<kNthetaAbs; ++itheta ) {
570 for ( Int_t isrc=0; isrc<kNtrackSources; ++isrc ) {
571 histoPattern = GetHistoName(kDCAVsPCheck, itheta, isrc);
572 TH2* histoCheck = (TH2*)GetSum(physSel, trigClassName, centralityRange, histoPattern);
573 if ( ! histoCheck ) continue;
574 currName = histoCheck->GetName();
575 currName.Append("_plotCut");
576 TCanvas* pdcaVsPcan = new TCanvas(currName.Data(), currName.Data(), igroup1*xshift,(igroup2+1)*yshift,600,600);
577 pdcaVsPcan->SetRightMargin(0.12);
578 pdcaVsPcan->SetLogz();
579 histoCheck->Draw("COLZ");
581 for ( Int_t icut=0; icut<nShowFuncs; ++icut ) {
582 currName = histoCheck->GetName();
583 currName.Append(Form("_cutFunc%i",icut));
584 TString currFormula = cutFormula;
585 currFormula.Append("/x");
586 TF1* cutFunction = new TF1(currName.Data(),currFormula.Data(), pMinCut, pMaxCut);
587 cutParam[icut][0] = sigmaMeasCut[itheta];
588 cutParam[icut][1] = nSigmaCut;
589 cutFunction->SetParameters(cutParam[icut]);
590 cutFunction->SetLineWidth(2);
591 cutFunction->SetLineColor(cutColor[icut]);
592 cutFunction->Draw("same");
593 } // loop on cut functions
598 ///////////////////////////
599 // Check Chi square prob //
600 ///////////////////////////
603 for ( Int_t itheta=0; itheta<kNthetaAbs; ++itheta ) {
604 for ( Int_t isrc=0; isrc<kNtrackSources; ++isrc ) {
605 histoPattern = GetHistoName(kChiProbVsP, itheta, isrc);
606 TH2* histo = (TH2*)GetSum(physSel, trigClassName, centralityRange, histoPattern);
607 if ( ! histo ) continue;
609 Int_t nPbins = histo->GetXaxis()->GetNbins();
610 //Int_t nPadX = (Int_t)TMath::Sqrt(nPbins);
611 //Int_t nPadY = nPadX;
612 //if ( nPadX * nPadY < nPbins ) nPadX++;
613 TCanvas* chi2probCan = 0x0;
619 for ( Int_t ibin=0; ibin<=nPbins; ++ibin ) {
620 currName = Form("hChiProb_%s_%s_%s", fSrcKeys->At(isrc)->GetName(), physSel.Data(), trigClassName.Data());
621 Int_t minBin = ( ibin == 0 ) ? 1 : ibin;
622 Int_t maxBin = ( ibin == 0 ) ? nPbins : ibin;
623 if ( ibin > 0 ) currName += Form("_pBin%i", ibin);
624 TH1* projHisto = histo->ProjectionY(currName.Data(), minBin, maxBin);
625 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)));
626 if ( projHisto->GetEntries() == 0 ) continue;
627 if ( ipad % (nPadX*nPadY) == 0 ) {
628 currName = histo->GetName();
629 currName += Form("Fit_can_%i", ican++);
630 chi2probCan = new TCanvas(currName.Data(), currName.Data(),igroup1*xshift,igroup2*yshift,600,600);
631 chi2probCan->Divide(nPadX,nPadY);
634 chi2probCan->cd(++ipad);
636 projHisto->Draw("e");
637 } // loop on momentum bins
639 } // loop on theta abs
642 //////////////////////
643 // Check sigma cuts //
644 //////////////////////
645 printf("\nReference sigma cut %g\n", refSigmaCut);
647 Float_t fracOfHeight = 0.35;
648 Float_t rightMargin = 0.03;
649 Int_t cutColors[14] = {kBlack, kRed, kBlue, kGreen, kCyan, kMagenta, kOrange, kViolet, kSpring, kGray, kSpring, kAzure, kPink, kYellow};
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. ) {
668 if ( orderCuts.GetSize() == 0 ) {
670 TAxis* axis = histo->GetYaxis();
671 Int_t nSigmaCuts = ( useCustomSigma ) ? fSigmaCuts.GetSize() : axis->GetNbins();
672 orderCuts.Set(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<orderCuts.GetSize(); ++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[orderCuts.GetSize()-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[orderCuts.GetSize()-1]));
805 } // loop on centrality
806 if ( leg ) leg->Draw("same");
807 } // loop on theta abs