AliPhysics  vAN-20150409 (2c6aa91)
 All Classes Namespaces Files Functions Variables Enumerations Enumerator Macros
AliAnalysisTaskMuonCuts.cxx
Go to the documentation of this file.
1 /**************************************************************************
2  * Copyright(c) 1998-2007, ALICE Experiment at CERN, All rights reserved. *
3  * *
4  * Author: The ALICE Off-line Project. *
5  * Contributors are mentioned in the code where appropriate. *
6  * *
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  **************************************************************************/
15 
16 /* $Id: AliAnalysisTaskMuonCuts.cxx 47782 2011-02-24 18:37:31Z martinez $ */
17 
28 
29 #define AliAnalysisTaskMuonCuts_cxx
30 
31 #include "AliAnalysisTaskMuonCuts.h"
32 
33 // ROOT includes
34 #include "TROOT.h"
35 #include "TH1.h"
36 #include "TH2.h"
37 #include "TH3.h"
38 #include "TAxis.h"
39 #include "TCanvas.h"
40 #include "TLegend.h"
41 #include "TMath.h"
42 #include "TObjString.h"
43 #include "TObjArray.h"
44 #include "TF1.h"
45 #include "TStyle.h"
46 #include "TArrayI.h"
47 //#include "TMCProcess.h"
48 
49 // STEER includes
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"
59 
60 #include "AliAnalysisManager.h"
61 
62 // PWG includes
63 #include "AliMergeableCollection.h"
64 #include "AliMuonEventCuts.h"
65 #include "AliMuonTrackCuts.h"
66 #include "AliAnalysisMuonUtility.h"
67 
68 
70 ClassImp(AliAnalysisTaskMuonCuts) // Class implementation in ROOT context
72 
73 
74 //________________________________________________________________________
75 AliAnalysisTaskMuonCuts::AliAnalysisTaskMuonCuts() :
76  AliVAnalysisMuon(),
77  fHistoTypeKeys(0x0),
78  fThetaAbsKeys(0x0),
79  fSigmaCuts(TArrayD())
80 {
82 }
83 
84 //________________________________________________________________________
85 AliAnalysisTaskMuonCuts::AliAnalysisTaskMuonCuts(const char *name, const AliMuonTrackCuts& cuts) :
86  AliVAnalysisMuon(name, cuts),
87  fHistoTypeKeys(0x0),
88  fThetaAbsKeys(0x0),
89  fSigmaCuts(TArrayD())
90 {
91  //
93  //
94  TString histoTypeKeys = "hDCAxVsP hDCAyVsP hPdcaVsP hPDCAVsPCheck hDCAVsPCheck hChiProbVsP hSigmaVsPt hSigmaVsEta";
95  fHistoTypeKeys = histoTypeKeys.Tokenize(" ");
96 
97  TString thetaAbsKeys = "ThetaAbs23 ThetaAbs310";
98  fThetaAbsKeys = thetaAbsKeys.Tokenize(" ");
99 
100  SetSigmaCuts();
101 }
102 
103 
104 //________________________________________________________________________
105 AliAnalysisTaskMuonCuts::~AliAnalysisTaskMuonCuts()
106 {
107  //
109  //
110 
111  delete fHistoTypeKeys;
112  delete fThetaAbsKeys;
113 }
114 
115 //___________________________________________________________________________
116 void AliAnalysisTaskMuonCuts::MyUserCreateOutputObjects()
117 {
118  TH1* histo = 0x0;
119  TString histoName = "";
120  for ( Int_t itheta=0; itheta<kNthetaAbs; ++itheta ) {
121  for ( Int_t isrc=0; isrc<kNtrackSources; ++isrc ) {
122  histoName = GetHistoName(kDCAxVsP, itheta, isrc);
123  histo = new TH2F(histoName.Data(), histoName.Data(), 50, 0., 200., 100, -50., 50.);
124  histo->SetXTitle("p (GeV/c)");
125  histo->SetYTitle("DCA_{x} (cm)");
126  AddObjectToCollection(histo);
127 
128  histoName = GetHistoName(kDCAyVsP, itheta, isrc);
129  histo = new TH2F(histoName.Data(), histoName.Data(), 50, 0., 200., 100, -50., 50.);
130  histo->SetXTitle("p (GeV/c)");
131  histo->SetYTitle("DCA_{y} (cm)");
132  AddObjectToCollection(histo);
133 
134  histoName = GetHistoName(kPdcaVsP, itheta, isrc);
135  histo = new TH2F(histoName.Data(), histoName.Data(), 100, 0., 400., 50, 0., 500.);
136  histo->SetXTitle("p (GeV/c)");
137  histo->SetYTitle("p#timesDCA (cm #times GeV/c)");
138  AddObjectToCollection(histo);
139 
140  histoName = GetHistoName(kPDCAVsPCheck, itheta, isrc);
141  histo = new TH2F(histoName.Data(), histoName.Data(), 100, 0., 800., 100, 0., 5000.);
142  histo->SetXTitle("p (GeV/c)");
143  histo->SetYTitle("p#timesDCA (cm #times GeV/c)");
144  AddObjectToCollection(histo);
145 
146  histoName = GetHistoName(kDCAVsPCheck, itheta, isrc);
147  histo = new TH2F(histoName.Data(), histoName.Data(), 100, 0., 800., 400, 0., 200.);
148  histo->SetXTitle("p (GeV/c)");
149  histo->SetYTitle("DCA (cm)");
150  AddObjectToCollection(histo);
151 
152  histoName = GetHistoName(kChiProbVsP, itheta, isrc);
153  histo = new TH2F(histoName.Data(), histoName.Data(), 50, 0., 200., 50, 0., 1.);
154  histo->SetXTitle("p (GeV/c)");
155  histo->SetYTitle("Chisquare prob.");
156  AddObjectToCollection(histo);
157 
158  histoName = GetHistoName(kSigmaVsPt, itheta, isrc);
159  histo = new TH2F(histoName.Data(), histoName.Data(), 100, 0., 100., fSigmaCuts.GetSize(), 0.5, 0.5+(Double_t)fSigmaCuts.GetSize());
160  histo->SetXTitle("p_{t} (GeV/c)");
161  histo->SetYTitle("#sigma_{p#timesDCA}");
162  for ( Int_t ibin=0; ibin<fSigmaCuts.GetSize(); ++ibin ) {
163  histo->GetYaxis()->SetBinLabel(ibin+1,Form("%g", fSigmaCuts[ibin]));
164  }
165  AddObjectToCollection(histo);
166 
167  histoName = GetHistoName(kSigmaVsEta, itheta, isrc);
168  histo = new TH2F(histoName.Data(), histoName.Data(), 25, -4.5, -2., fSigmaCuts.GetSize(), 0.5, 0.5+(Double_t)fSigmaCuts.GetSize());
169  histo->SetXTitle("#eta");
170  histo->SetYTitle("#sigma_{p#timesDCA}");
171  for ( Int_t ibin=0; ibin<fSigmaCuts.GetSize(); ++ibin ) {
172  histo->GetYaxis()->SetBinLabel(ibin+1,Form("%g", fSigmaCuts[ibin]));
173  }
174  AddObjectToCollection(histo);
175  } // loop on track sources
176  } // loop on theta abs
177 
178  fMuonTrackCuts->Print();
179 
180 }
181 
182 //________________________________________________________________________
183 TString AliAnalysisTaskMuonCuts::GetHistoName(Int_t histoTypeIndex, Int_t thetaAbsIndex, Int_t srcIndex)
184 {
186  TString histoName = Form("%s%s%s", fHistoTypeKeys->At(histoTypeIndex)->GetName(), fThetaAbsKeys->At(thetaAbsIndex)->GetName(), fSrcKeys->At(srcIndex)->GetName());
187 
188  return histoName;
189 }
190 
191 //________________________________________________________________________
192 void AliAnalysisTaskMuonCuts::ProcessEvent(TString physSel, const TObjArray& selectTrigClasses, TString centrality)
193 {
194  //
196  //
197 
198  TString histoName = "";
199  AliVParticle* track = 0x0;
200  Int_t nTracks = AliAnalysisMuonUtility::GetNTracks(InputEvent());
201  for (Int_t itrack = 0; itrack < nTracks; itrack++) {
202  track = AliAnalysisMuonUtility::GetTrack(itrack, InputEvent());
203  fMuonTrackCuts->CustomParam()->SetNSigmaPdca(1.e10);
204  if ( ! fMuonTrackCuts->IsSelected(track) ) continue;
205 
206  Double_t thetaAbsEndDeg = AliAnalysisMuonUtility::GetThetaAbsDeg(track);
207  Int_t thetaAbsBin = ( thetaAbsEndDeg < 3. ) ? kThetaAbs23 : kThetaAbs310;
208 
209  Int_t trackSrc = GetParticleType(track);
210 
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;
215 
216  Double_t sigmaPdca = fMuonTrackCuts->IsThetaAbs23(track) ? fMuonTrackCuts->GetMuonTrackCutsParam().GetSigmaPdca23() : fMuonTrackCuts->GetMuonTrackCutsParam().GetSigmaPdca310();
217  Double_t chi2 = pDca / sigmaPdca;
218  chi2 *= chi2;
219  Double_t chiProb = TMath::Prob(chi2, 2);
220 
221  Double_t pTot = track->P();
222  Double_t pt = track->Pt();
223  Double_t eta = track->Eta();
224 
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());
229 
230  histoName = GetHistoName(kDCAyVsP, thetaAbsBin, trackSrc);
231  ((TH2*)GetMergeableObject(physSel, trigClassName, centrality, histoName))->Fill(pTot, dcaAtVz.Y());
232 
233  histoName = GetHistoName(kPdcaVsP, thetaAbsBin, trackSrc);
234  ((TH2*)GetMergeableObject(physSel, trigClassName, centrality, histoName))->Fill(pTot, pDca);
235 
236  histoName = GetHistoName(kPDCAVsPCheck, thetaAbsBin, trackSrc);
237  ((TH2*)GetMergeableObject(physSel, trigClassName, centrality, histoName))->Fill(pTot, pDca);
238 
239  histoName = GetHistoName(kDCAVsPCheck, thetaAbsBin, trackSrc);
240  ((TH2*)GetMergeableObject(physSel, trigClassName, centrality, histoName))->Fill(pTot, dca);
241 
242  histoName = GetHistoName(kChiProbVsP, thetaAbsBin, trackSrc);
243  ((TH2*)GetMergeableObject(physSel, trigClassName, centrality, histoName))->Fill(pTot, chiProb);
244  } // loop on selected trigger classes
245 
246  for ( Int_t isigma=0; isigma<fSigmaCuts.GetSize(); ++isigma) {
247  fMuonTrackCuts->CustomParam()->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
256  } // loop on sigmas
257  }
258 }
259 
260 //________________________________________________________________________
261 void AliAnalysisTaskMuonCuts::SetSigmaCuts(Int_t nSigmaCuts, Double_t* sigmaCuts)
262 {
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);
269 // }
270 // else {
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);
274  // }
275  }
276  else {
277  fSigmaCuts.Set(nSigmaCuts, sigmaCuts);
278  }
279 }
280 
281 //________________________________________________________________________
282 void AliAnalysisTaskMuonCuts::Terminate(Option_t *) {
283  //
285  //
286 
287  AliVAnalysisMuon::Terminate("");
288 
289  if ( ! fMergeableCollection ) return;
290 
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();
296 
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();
307  }
308  }
309  delete optArray;
310 
311  fMuonTrackCuts->Print("param");
312 
313  Int_t srcColors[kNtrackSources] = {kBlack, kRed, kSpring, kTeal, kBlue, kViolet, kMagenta, kOrange, kGray};
314 
315  TCanvas* can = 0x0;
316  Int_t xshift = 100;
317  Int_t yshift = 100;
318  Int_t igroup1 = -1;
319  Int_t igroup2 = 0;
320 
322  // Reco DCA //
324  igroup1++;
325  igroup2 = 0;
326  TString histoName = "", currName = "", histoPattern = "", drawOpt = "";
327  currName = Form("%s_recoDCA", GetName());
328  can = new TCanvas(currName.Data(),"Reco DCA",igroup1*xshift,igroup2*yshift,600,600);
329  can->Divide(2,2);
330  igroup2++;
331  Int_t recoDcaHisto[2] = {kDCAxVsP, kDCAyVsP};
332  TString dcaName[2] = {"DCAx", "DCAy"};
333  Double_t averageDca[4] = {0., 0.};
334  printf("\nAverage reconstructed DCA:\n");
335  TF1* fitFuncMeanDca = new TF1("fitFuncMeanDca","gausn",-20.,20.);
336  fitFuncMeanDca->SetParNames("Norm", "Mean", "Sigma");
337  for ( Int_t itheta=0; itheta<kNthetaAbs; ++itheta ) {
338  for ( Int_t ihisto=0; ihisto<2; ++ihisto ) {
339  TH2* histo = 0x0;
340  histoPattern = "";
341  histoPattern = Form("%s%s*", fHistoTypeKeys->At(recoDcaHisto[ihisto])->GetName(), fThetaAbsKeys->At(itheta)->GetName());
342  histo = (TH2*)GetSum(physSel, trigClassName, centralityRange, histoPattern);
343  if ( ! histo ) continue;
344 
345  TH1* meanDcaVsP = histo->ProjectionX(Form("mean%sVsP_%s", dcaName[ihisto].Data(), fThetaAbsKeys->At(itheta)->GetName()));
346  meanDcaVsP->Reset();
347  meanDcaVsP->SetTitle(Form("Mean %s vs. p %s", dcaName[ihisto].Data(), fThetaAbsKeys->At(itheta)->GetName()));
348  meanDcaVsP->SetYTitle(Form("<%s> (cm)", dcaName[ihisto].Data()));
349  meanDcaVsP->SetStats(kFALSE);
350 
351  Int_t nPbins = histo->GetXaxis()->GetNbins();
352  //Int_t nPadX = (Int_t)TMath::Sqrt(nPbins);
353  //Int_t nPadY = nPadX;
354  //if ( nPadX * nPadY < nPbins ) nPadX++;
355  TCanvas* meanDcaFitCan = 0x0;
356  Int_t nPadX = 5;
357  Int_t nPadY = 5;
358  Int_t ipad = 0;
359  Int_t ican = 0;
360 
361  for ( Int_t ibin=2; ibin<=nPbins; ++ibin ) {
362  currName = Form("hMean%s_%s_%s", dcaName[ihisto].Data(), physSel.Data(), trigClassName.Data());
363  Int_t minBin = ( ibin == 0 ) ? 1 : ibin;
364  Int_t maxBin = ( ibin == 0 ) ? nPbins : ibin;
365  if ( ibin > 0 ) currName += Form("_pBin%i", ibin);
366  TH1* projHisto = histo->ProjectionY(currName.Data(), minBin, maxBin, "e");
367  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)));
368 
369  if ( projHisto->GetEntries() == 0 ) continue;
370  if ( ipad % (nPadX*nPadY) == 0 ) {
371  currName = histo->GetName();
372  currName += Form("Fit_can_%i", ican++);
373  meanDcaFitCan = new TCanvas(currName.Data(), currName.Data(),igroup1*xshift,igroup2*yshift,600,600);
374  meanDcaFitCan->Divide(nPadX,nPadY);
375  ipad = 0;
376  }
377  meanDcaFitCan->cd(++ipad);
378  gPad->SetLogy();
379  if ( projHisto->Integral() > 50 ) {
380  fitFuncMeanDca->SetParameter(0, projHisto->Integral());
381  fitFuncMeanDca->SetParameter(1, projHisto->GetMean());
382  fitFuncMeanDca->SetParameter(2, projHisto->GetRMS());
383  Double_t fitDcaLim = ( ibin == 0 ) ? 40. : 40./((Double_t)ibin);
384  fitDcaLim = TMath::Max(5., fitDcaLim);
385  projHisto->Fit(fitFuncMeanDca, "RQ", "e", -fitDcaLim, fitDcaLim);
386  Double_t chi2 = fitFuncMeanDca->GetChisquare();
387  Double_t ndf = fitFuncMeanDca->GetNDF();
388  if ( ndf <= 0.) continue;
389  if ( chi2 / ndf > 100. ) continue;
390  Double_t meanDca = fitFuncMeanDca->GetParameter(1);
391  Double_t meanDcaErr = fitFuncMeanDca->GetParError(1);
392  meanDcaVsP->SetBinContent(ibin, meanDca);
393  meanDcaVsP->SetBinError(ibin, meanDcaErr);
394  //if ( ibin == 0 ) printf("%s meanDca = %g +- %g\n", fThetaAbsKeys->At(itheta)->GetName(), meanDca, meanDcaErr);
395  }
396  else projHisto->Draw("e");
397  } // loop on momentum bins
398  can->cd(2*itheta+ihisto+1);
399  //meanDcaVsP->SetLineColor(srcColors[isrc]);
400  //meanDcaVsP->SetMarkerColor(srcColors[isrc]);
401  //meanDcaVsP->SetMarkerStyle(20+isrc);
402  meanDcaVsP->Fit("pol0","Q");
403  TF1* trendFit = (TF1*)meanDcaVsP->GetListOfFunctions()->FindObject("pol0");
404  if ( trendFit ) {
405  averageDca[2*itheta+ihisto] = trendFit->GetParameter(0);
406  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  }
408  //drawOpt = "esame";
409  //leg->AddEntry(meanDcaVsP, fSrcKeys->At(isrc)->GetName(), "lp");
410  //} // loop on src
411  //can->cd(itheta+1);
412  //leg->Draw("same");
413 
414  //can->cd(ipad++);
415  //histo->Draw();
416  //meanDca[ihisto] = histo->GetMean();
417  } // loop on histo type
418  } // loop on theta abs
419 
420  for ( Int_t itheta=0; itheta<kNthetaAbs; ++itheta ) {
421  printf("muonCuts->SetMeanDCA(%g, %g, 0.); // %s\n", averageDca[2*itheta], averageDca[2*itheta+1], fThetaAbsKeys->At(itheta)->GetName());
422  }
423 
425  // Fit pDCA //
427  Double_t nSigmaCut = fMuonTrackCuts->GetMuonTrackCutsParam().GetNSigmaPdca(); //6.;
428  Double_t sigmaMeasCut[2] = { fMuonTrackCuts->GetMuonTrackCutsParam().GetSigmaPdca23(), fMuonTrackCuts->GetMuonTrackCutsParam().GetSigmaPdca310()}; //{99., 54.}; //{120., 63.};
429  Double_t relPResolution = fMuonTrackCuts->GetMuonTrackCutsParam().GetRelPResolution(); //4.5e-4; //6.e-3;//8.e-4;
430  Double_t angleResolution = 535.*fMuonTrackCuts->GetMuonTrackCutsParam().GetSlopeResolution(); //6.e-4;
431  Double_t pMinCut = 0.1;
432  Double_t pMaxCut = 800.;
433  const Int_t kNCutFuncs = 2;
434  Int_t nShowFuncs = 1;
435  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)";
436  Double_t cutParam[kNCutFuncs][4] = {{sigmaMeasCut[0], nSigmaCut, relPResolution, angleResolution}, {sigmaMeasCut[0], nSigmaCut, 0., 0.32}};
437  Int_t cutColor[kNCutFuncs] = {kBlack, kRed};
438  igroup1++;
439  igroup2 = 0;
440  can = new TCanvas("pdcaSigmaFit","Sigma vs P fit",igroup1*xshift,igroup2*yshift,600,600);
441  can->Divide(2,1);
442  igroup2++;
443  TF1* fitFunc = new TF1("fitFunc", "x*gausn", 0., 400.);
444  fitFunc->SetParNames("Norm", "Mean", "Sigma");
445  gStyle->SetOptFit(1111);
446  Double_t xMinFit[2] = {0., 0.};
447  Double_t xMaxFit[2] = {260., 150.}; //{320., 150.}; // {360., 180.};
448  printf("\nSigma p x DCA:\n");
449  Double_t averageSigmaPdca[kNtrackSources*kNthetaAbs] = {0.};
450  for ( Int_t itheta=0; itheta<kNthetaAbs; ++itheta ) {
451  TLegend* leg = new TLegend(0.7, 0.7, 0.9, 0.9);
452  leg->SetBorderSize(1);
453  for ( Int_t isrc=0; isrc<kNtrackSources; ++isrc ) {
454  histoPattern = GetHistoName(kPdcaVsP, itheta, isrc);
455  TH2* histo = (TH2*)GetSum(physSel, trigClassName, centralityRange, histoPattern);
456  if ( ! histo ) continue;
457  if ( histo->Integral() < 200 ) {
458  delete histo;
459  continue;
460  }
461 
462  TH1* sigmaVsP = histo->ProjectionX(Form("sigma%s_%s_%s", fHistoTypeKeys->At(kPdcaVsP)->GetName(), fThetaAbsKeys->At(itheta)->GetName(), fSrcKeys->At(isrc)->GetName()));
463  sigmaVsP->Reset();
464  sigmaVsP->SetTitle(Form("#sigma_{p#timesDCA} vs. p %s %s", fThetaAbsKeys->At(itheta)->GetName(), fSrcKeys->At(isrc)->GetName()));
465  sigmaVsP->SetYTitle("#sigma_{p#timesDCA} (cm #times GeV/c)");
466  sigmaVsP->SetStats(kFALSE);
467 
468  Int_t nPbins = histo->GetXaxis()->GetNbins();
469  //Int_t nPadX = (Int_t)TMath::Sqrt(nPbins);
470  //Int_t nPadY = nPadX;
471  //if ( nPadX * nPadY < nPbins ) nPadX++;
472  TCanvas* pdcaFitCan = 0x0;
473  Int_t nPadX = 5;
474  Int_t nPadY = 5;
475  Int_t ipad = 0;
476  Int_t ican = 0;
477 
478  for ( Int_t ibin=2; ibin<=nPbins; ++ibin ) {
479  currName = Form("hPDCA_%s_%s_%s_%s", fSrcKeys->At(isrc)->GetName(), physSel.Data(), trigClassName.Data(), fThetaAbsKeys->At(itheta)->GetName());
480  Int_t minBin = ( ibin == 0 ) ? 1 : ibin;
481  Int_t maxBin = ( ibin == 0 ) ? nPbins : ibin;
482  if ( ibin > 0 ) currName += Form("_pBin%i", ibin);
483  TH1* projHisto = histo->ProjectionY(currName.Data(), minBin, maxBin, "e");
484  if ( projHisto->GetEntries() == 0 ) continue;
485  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)));
486  if ( ipad % (nPadX*nPadY) == 0 ) {
487  currName = histo->GetName();
488  currName += Form("Fit_can_%i", ican++);
489  pdcaFitCan = new TCanvas(currName.Data(), currName.Data(),igroup1*xshift,igroup2*yshift,600,600);
490  pdcaFitCan->Divide(nPadX,nPadY);
491  ipad = 0;
492  }
493  pdcaFitCan->cd(++ipad);
494  if ( projHisto->Integral() > 0.) projHisto->Scale(1./projHisto->Integral());
495  gPad->SetLogy();
496  if ( projHisto->GetEntries() > 50 ) {
497  fitFunc->SetParameter(0, projHisto->Integral());
498  fitFunc->SetParameter(1, projHisto->GetMean());
499  fitFunc->SetParameter(2, projHisto->GetRMS());
500  projHisto->Fit(fitFunc, "RQ", "e", xMinFit[itheta], xMaxFit[itheta]);
501  Double_t chi2 = fitFunc->GetChisquare();
502  Double_t ndf = fitFunc->GetNDF();
503  if ( ndf <= 0.) continue;
504  if ( chi2 / ndf > 100. ) continue;
505  Double_t sigma = TMath::Abs(fitFunc->GetParameter(2));
506  Double_t sigmaErr = fitFunc->GetParError(2);
507  sigmaVsP->SetBinContent(ibin, sigma);
508  sigmaVsP->SetBinError(ibin, sigmaErr);
509  //if ( ibin == 0 ) printf("%s %s sigma = %g +- %g\n", fThetaAbsKeys->At(itheta)->GetName(), fSrcKeys->At(isrc)->GetName(), sigma, sigmaErr);
510  }
511  else projHisto->Draw("e");
512  } // loop on momentum bins
513  can->cd(itheta+1);
514  sigmaVsP->SetLineColor(srcColors[isrc]);
515  sigmaVsP->SetMarkerColor(srcColors[isrc]);
516  sigmaVsP->SetMarkerStyle(20+isrc);
517  drawOpt = "e";
518  if ( gPad->GetListOfPrimitives()->GetEntries() != 0 ) drawOpt += "same";
519  sigmaVsP->Draw(drawOpt.Data());
520  sigmaVsP->Fit("pol0","Q",drawOpt.Data());
521  TF1* trendFit = (TF1*)sigmaVsP->GetListOfFunctions()->FindObject("pol0");
522  if ( trendFit ) {
523  averageSigmaPdca[kNtrackSources*itheta+isrc] = trendFit->GetParameter(0);
524  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()));
525  }
526  leg->AddEntry(sigmaVsP, fSrcKeys->At(isrc)->GetName(), "lp");
527  } // loop on src
528  can->cd(itheta+1);
529  leg->Draw("same");
530 
531  for ( Int_t icut=0; icut<nShowFuncs; ++icut ) {
532  currName = Form("sigmaCut_%s_%i", fThetaAbsKeys->At(itheta)->GetName(), icut);
533  TF1* cutFunction = new TF1(currName.Data(), cutFormula.Data(), pMinCut, pMaxCut);
534  cutParam[icut][0] = sigmaMeasCut[itheta];
535  cutParam[icut][1] = 1.;
536  cutFunction->SetParameters(cutParam[icut]);
537  cutFunction->SetLineColor(cutColor[icut]);
538  cutFunction->SetLineWidth(2);
539  cutFunction->Draw("same");
540  }
541  } // loop on theta abs
542 
543  for ( Int_t isrc=0; isrc<kNtrackSources; ++isrc ) {
544  printf("muonCuts->SetSigmaPdca(%g, %g); // %s\n", averageSigmaPdca[isrc], averageSigmaPdca[kNtrackSources+isrc], fSrcKeys->At(isrc)->GetName());
545  }
546  printf("\n");
547 
548  igroup2++;
549  for ( Int_t icheck=0; icheck<2; icheck++ ) {
550  Int_t hDraw = ( icheck == 0 ) ? kPDCAVsPCheck : kDCAVsPCheck;
551  for ( Int_t itheta=0; itheta<kNthetaAbs; ++itheta ) {
552  for ( Int_t isrc=0; isrc<kNtrackSources; ++isrc ) {
553  histoPattern = GetHistoName(hDraw, itheta, isrc);
554  TH2* histoCheck = (TH2*)GetSum(physSel, trigClassName, centralityRange, histoPattern);
555  if ( ! histoCheck ) continue;
556  currName = histoCheck->GetName();
557  currName.Append("_plotCut");
558  TCanvas* pdcaVsPcan = new TCanvas(currName.Data(), currName.Data(), igroup1*xshift,(igroup2+icheck)*yshift,600,600);
559  pdcaVsPcan->SetRightMargin(0.12);
560  pdcaVsPcan->SetLogy();
561  pdcaVsPcan->SetLogz();
562  histoCheck->Draw("COLZ");
563 
564  for ( Int_t icut=0; icut<nShowFuncs; ++icut ) {
565  currName = histoCheck->GetName();
566  currName.Append(Form("_cutFunc%i",icut));
567  TString currFormula = cutFormula;
568  if ( icheck == 1 ) currFormula.Append("/x");
569  TF1* cutFunction = new TF1(currName.Data(),currFormula.Data(), pMinCut, pMaxCut);
570  cutParam[icut][0] = sigmaMeasCut[itheta];
571  cutParam[icut][1] = nSigmaCut;
572  cutFunction->SetParameters(cutParam[icut]);
573  cutFunction->SetLineWidth(2);
574  cutFunction->SetLineColor(cutColor[icut]);
575  cutFunction->Draw("same");
576  } // loop on cut functions
577  } // loop on src
578  } //loop on theta
579  } // loop on check
580 
581 
583  // Check Chi square prob //
585  igroup1++;
586  igroup2 = 0;
587  for ( Int_t itheta=0; itheta<kNthetaAbs; ++itheta ) {
588  for ( Int_t isrc=0; isrc<kNtrackSources; ++isrc ) {
589  histoPattern = GetHistoName(kChiProbVsP, itheta, isrc);
590  TH2* histo = (TH2*)GetSum(physSel, trigClassName, centralityRange, histoPattern);
591  if ( ! histo ) continue;
592 
593  Int_t nPbins = histo->GetXaxis()->GetNbins();
594  //Int_t nPadX = (Int_t)TMath::Sqrt(nPbins);
595  //Int_t nPadY = nPadX;
596  //if ( nPadX * nPadY < nPbins ) nPadX++;
597  TCanvas* chi2probCan = 0x0;
598  Int_t nPadX = 5;
599  Int_t nPadY = 5;
600  Int_t ipad = 0;
601  Int_t ican = 0;
602 
603  for ( Int_t ibin=0; ibin<=nPbins; ++ibin ) {
604  currName = Form("hChiProb_%s_%s_%s_%s", fSrcKeys->At(isrc)->GetName(), physSel.Data(), trigClassName.Data(), fThetaAbsKeys->At(itheta)->GetName());
605  Int_t minBin = ( ibin == 0 ) ? 1 : ibin;
606  Int_t maxBin = ( ibin == 0 ) ? nPbins : ibin;
607  if ( ibin > 0 ) currName += Form("_pBin%i", ibin);
608  TH1* projHisto = histo->ProjectionY(currName.Data(), minBin, maxBin);
609  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)));
610  if ( projHisto->GetEntries() == 0 ) continue;
611  if ( ipad % (nPadX*nPadY) == 0 ) {
612  currName = histo->GetName();
613  currName += Form("Fit_can_%i", ican++);
614  chi2probCan = new TCanvas(currName.Data(), currName.Data(),igroup1*xshift,igroup2*yshift,600,600);
615  chi2probCan->Divide(nPadX,nPadY);
616  ipad = 0;
617  }
618  chi2probCan->cd(++ipad);
619  gPad->SetLogy();
620  projHisto->Draw("e");
621  } // loop on momentum bins
622  } // loop on src
623  } // loop on theta abs
624 
625 
627  // Check sigma cuts //
629  printf("\nReference sigma cut %g\n", refSigmaCut);
630 
631  Float_t fracOfHeight = 0.35;
632  Float_t rightMargin = 0.03;
633  Int_t cutColors[14] = {kBlack, kRed, kBlue, kGreen, kCyan, kMagenta, kOrange, kViolet, kSpring, kGray, kSpring, kAzure, kPink, kYellow};
634  TArrayI orderCuts;
635  Int_t checkHistos[2] = {kSigmaVsPt, kSigmaVsEta};
636  Bool_t useCustomSigma = furtherOpt.Contains("CUSTOMSIGMA");
637 
638  igroup1++;
639  igroup2 = 0;
640 
641  for ( Int_t ihisto=0; ihisto<2; ++ihisto ) {
642  for ( Int_t isrc=0; isrc<kNtrackSources; ++isrc ) {
643  can = 0x0;
644  for ( Int_t itheta=0; itheta<kNthetaAbs; ++itheta ) {
645  histoPattern = GetHistoName(checkHistos[ihisto], itheta, isrc);
646  TH2* histo = (TH2*)GetSum(physSel, trigClassName, centralityRange, histoPattern);
647  if ( ! histo ) continue;
648  if ( histo->Integral() == 0. ) {
649  delete histo;
650  continue;
651  }
652  if ( orderCuts.GetSize() == 0 ) {
653  // Re-order axis
654  TAxis* axis = histo->GetYaxis();
655  Int_t nSigmaCuts = ( useCustomSigma ) ? fSigmaCuts.GetSize() : axis->GetNbins();
656  orderCuts.Set(nSigmaCuts);
657  Int_t countBin = 0;
658  for ( Int_t isigma=0; isigma<axis->GetNbins(); ++isigma ) {
659  TString currLabel = axis->GetBinLabel(isigma+1);
660  Double_t currSigma = currLabel.Atof();
661  if ( useCustomSigma ) {
662  Bool_t foundMatch = kFALSE;
663  for ( Int_t jsigma=0; jsigma<fSigmaCuts.GetSize(); ++jsigma ) {
664  if ( TMath::Abs(currSigma - fSigmaCuts[jsigma]) < 1e-4 ) {
665  foundMatch = kTRUE;
666  break;
667  }
668  }
669  if ( ! foundMatch ) continue;
670  }
671  Int_t currBin = ( TMath::Abs(currSigma - refSigmaCut) < 1e-4) ? 0 : ++countBin;
672  orderCuts[currBin] = isigma+1;
673  }
674  }
675  currName = histo->GetName();
676  currName.Append("_can");
677  can = new TCanvas(currName.Data(), currName.Data(), igroup1*xshift,igroup2*yshift, 600, 600);
678  TLegend* leg = new TLegend(0.6,0.6,0.9,0.9);
679  leg->SetBorderSize(1);
680  can->Divide(1,2,0,0);
681  can->cd(1);
682  gPad->SetPad(0., fracOfHeight, 0.99, 0.99);
683  //gPad->SetTopMargin(0.08);
684  gPad->SetTopMargin(0.03);
685  gPad->SetRightMargin(rightMargin);
686  gPad->SetLogy();
687 
688  can->cd(2);
689  gPad->SetPad(0., 0., 0.99, fracOfHeight);
690  gPad->SetRightMargin(rightMargin);
691  gPad->SetBottomMargin(0.08/fracOfHeight);
692 
693  TH1* refCutHisto = 0x0;
694  TString legTitle = "";
695  for ( Int_t isigma=0; isigma<orderCuts.GetSize(); ++isigma ) {
696  currName = histo->GetName();
697  currName.Append(Form("_sigma%i", isigma));
698  Int_t currBin = orderCuts[isigma];
699  TH1* projectHisto = histo->ProjectionX(currName.Data(), currBin, currBin);
700  projectHisto->SetLineColor(cutColors[isigma]);
701  projectHisto->SetMarkerColor(cutColors[isigma]);
702  projectHisto->SetMarkerStyle(20+isigma);
703  projectHisto->SetStats(kFALSE);
704  projectHisto->SetTitle("");
705  can->cd(1);
706  drawOpt = "e";
707  if ( gPad->GetListOfPrimitives()->GetEntries() != 0 ) drawOpt += "same";
708  projectHisto->Draw(drawOpt.Data());
709  TString currLabel = histo->GetYaxis()->GetBinLabel(currBin);
710  TString cutName = ( TMath::Abs(currLabel.Atof() ) > 1000. ) ? "Total" : Form("%s #sigma cut", currLabel.Data());
711  leg->AddEntry(projectHisto, cutName.Data(), "lp");
712  if ( ! refCutHisto ) {
713  refCutHisto = projectHisto;
714  legTitle = Form("(pass n #sigma) / (pass %s #sigma)", currLabel.Data());
715  continue;
716  }
717  currName.Append("_ratio");
718  TH1* histoRatio = (TH1*)projectHisto->Clone(currName.Data());
719  histoRatio->Sumw2();
720  TH1* auxNum = projectHisto;
721  TH1* auxDen = refCutHisto;
722  Bool_t mustInvert = kFALSE;
723  if ( auxNum->Integral()>auxDen->Integral() ) {
724  auxNum = refCutHisto;
725  auxDen = projectHisto;
726  mustInvert = kTRUE;
727  }
728  histoRatio->Divide(auxNum,auxDen,1.,1.,"B");
729  if ( mustInvert ) {
730  TH1* auxHistoOne = static_cast<TH1*>(histoRatio->Clone("auxHistoOne"));
731  for ( Int_t ibin=1; ibin<=auxHistoOne->GetXaxis()->GetNbins(); ibin++ ) {
732  auxHistoOne->SetBinContent(ibin,1.);
733  auxHistoOne->SetBinError(ibin,0.);
734  }
735  histoRatio->Divide(auxHistoOne,histoRatio);
736  delete auxHistoOne;
737  }
738  histoRatio->SetTitle("");
739  histoRatio->GetYaxis()->SetTitle(legTitle.Data());
740  histoRatio->GetXaxis()->SetLabelSize(0.04/fracOfHeight);
741  histoRatio->GetXaxis()->SetTitleSize(0.035/fracOfHeight);
742  histoRatio->GetYaxis()->SetLabelSize(0.03/fracOfHeight);
743  histoRatio->GetYaxis()->SetTitleSize(0.03/fracOfHeight);
744  histoRatio->GetXaxis()->SetTitleOffset(1.);
745  histoRatio->GetYaxis()->SetTitleOffset(0.6);
746  histoRatio->GetYaxis()->SetRangeUser(0.5,1.5);
747 
748  can->cd(2);
749  drawOpt = "e";
750  if ( gPad->GetListOfPrimitives()->GetEntries() != 0 ) drawOpt += "same";
751  histoRatio->Draw(drawOpt.Data());
752  }// loop on sigma cuts
753  can->cd(1);
754  leg->Draw("same");
755  } // loop on theta abs
756  } // loop on src
757  } // loop on histo type
758 
759  igroup1++;
760  igroup2=0;
761  Double_t ptMin[] = {0., 2., 4., 15., 60.};
762  Int_t nPtMins = sizeof(ptMin)/sizeof(ptMin[0]);
763  for ( Int_t iptmin=0; iptmin<nPtMins; ++iptmin) {
764  for ( Int_t isrc=0; isrc<kNtrackSources; ++isrc ) {
765  can = 0x0;
766  for ( Int_t itheta=0; itheta<kNthetaAbs; ++itheta ) {
767  TLegend* leg = 0x0;
768  histoName = GetHistoName(kSigmaVsPt, itheta, isrc);
769  for ( Int_t icent=1; icent<=GetCentralityClasses()->GetNbins(); ++icent ) {
770  TH2* histo = (TH2*)GetSum(physSel, trigClassName, GetCentralityClasses()->GetBinLabel(icent), histoName);
771  if ( ! histo ) continue;
772  Int_t ptMinBin = histo->GetXaxis()->FindBin(ptMin[iptmin]);
773  Int_t ptMaxBin = histo->GetXaxis()->GetNbins()+1;
774  currName = histo->GetName();
775  currName += Form("_%s_ptMin%i", GetCentralityClasses()->GetBinLabel(icent), TMath::Nint(ptMin[iptmin]));
776  TH1* projectHisto = histo->ProjectionY(currName.Data(), ptMinBin, ptMaxBin);
777  if ( projectHisto->GetMaximum() < 50. ) {
778  delete projectHisto;
779  continue;
780  }
781  if ( ! can ) {
782  currName = Form("CutEffect_%s_ptMin%i", fSrcKeys->At(isrc)->GetName(), TMath::Nint(ptMin[iptmin]));
783  can = new TCanvas(currName.Data(), currName.Data(), igroup1*xshift,igroup2*yshift, 600, 600);
784  can->Divide(2,1);
785  }
786  if ( ! leg ) {
787  leg = new TLegend(0.6,0.6,0.9,0.9);
788  leg->SetBorderSize(1);
789  }
790  can->cd(itheta+1);
791  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)));
792  projectHisto->SetLineColor(cutColors[icent-1]);
793  projectHisto->SetMarkerColor(cutColors[icent-1]);
794  projectHisto->SetMarkerStyle(19+icent);
795  projectHisto->SetStats(0);
796  drawOpt = "e";
797  if ( gPad->GetListOfPrimitives()->GetEntries() != 0 ) drawOpt += "same";
798  projectHisto->Draw(drawOpt.Data());
799  leg->AddEntry(projectHisto, GetCentralityClasses()->GetBinLabel(icent), "lp");
800  Double_t keptEvents = projectHisto->GetBinContent(orderCuts[0]);
801  Double_t totalEvents = projectHisto->GetBinContent(orderCuts[orderCuts.GetSize()-1]);
802  Double_t accepted = ( totalEvents > 0. ) ? keptEvents / totalEvents : 1.;
803  Double_t acceptedErr = ( totalEvents > 0. ) ? TMath::Sqrt(accepted*(1.-accepted)/totalEvents) : 1.;
804  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.);
805  //printf(" rejected %g total %g (%s vs %s)\n",totalEvents-keptEvents,totalEvents,projectHisto->GetXaxis()->GetBinLabel(orderCuts[0]),projectHisto->GetXaxis()->GetBinLabel(orderCuts[orderCuts.GetSize()-1]));
806  } // loop on centrality
807  if ( leg ) leg->Draw("same");
808  } // loop on theta abs
809  } // loop on sources
810  } // loop on pt min
811 }
ClassImp
ClassImp(AliAnalysisTaskTriggerRates) AliAnalysisTaskTriggerRates
Definition: AliAnalysisTaskTriggerRates.cxx:34
ptMin
Double_t ptMin
Definition: MuonTrackingEfficiency.C:52
AliAnalysisTaskMuonCuts::kPDCAVsPCheck
p x DCA vs momentum (check beam gas)
Definition: AliAnalysisTaskMuonCuts.h:49
AliAnalysisTaskMuonCuts::fThetaAbsKeys
TObjArray * fThetaAbsKeys
Name of theta at absorber end.
Definition: AliAnalysisTaskMuonCuts.h:60
AliAnalysisTaskMuonCuts::kChiProbVsP
Chi square probability vs momentum.
Definition: AliAnalysisTaskMuonCuts.h:51
AliAnalysisTaskMuonCuts::kSigmaVsEta
eta distribution for different p x DCA sigma cuts
Definition: AliAnalysisTaskMuonCuts.h:53
AliAnalysisTaskMuonCuts::Terminate
virtual void Terminate(Option_t *option)
Definition: AliAnalysisTaskMuonCuts.cxx:282
AliAnalysisTaskMuonCuts::ProcessEvent
void ProcessEvent(TString physSel, const TObjArray &selectTrigClasses, TString centrality)
Definition: AliAnalysisTaskMuonCuts.cxx:192
AliAnalysisTaskMuonCuts::fHistoTypeKeys
TObjArray * fHistoTypeKeys
Base histogram name.
Definition: AliAnalysisTaskMuonCuts.h:59
AliAnalysisTaskMuonCuts::kSigmaVsPt
pt distribution for different p x DCA sigma cuts
Definition: AliAnalysisTaskMuonCuts.h:52
AliAnalysisTaskMuonCuts.h
AliAnalysisTaskMuonCuts::fSigmaCuts
TArrayD fSigmaCuts
List of sigma cuts.
Definition: AliAnalysisTaskMuonCuts.h:61
AliAnalysisTaskMuonCuts::kPdcaVsP
p x DCA vs momentum (binning for fit)
Definition: AliAnalysisTaskMuonCuts.h:48
AliAnalysisTaskMuonCuts::SetSigmaCuts
void SetSigmaCuts(Int_t nSigmaCuts=-1, Double_t *sigmaCuts=0x0)
Definition: AliAnalysisTaskMuonCuts.cxx:261
AliAnalysisTaskMuonCuts::GetHistoName
TString GetHistoName(Int_t histoTypeIndex, Int_t thetaAbsIndex, Int_t srcIndex)
Definition: AliAnalysisTaskMuonCuts.cxx:183