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