Improve plotted histograms in Terminate (Diego)
[u/mrichter/AliRoot.git] / PWGPP / MUON / lite / AliAnalysisTaskMuonCuts.cxx
CommitLineData
35ac82fa 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
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
26///
27/// \author Diego Stocco
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"
9e6be1fd 46#include "TArrayI.h"
35ac82fa 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"
ebba9ef0 64#include "AliMuonEventCuts.h"
35ac82fa 65#include "AliMuonTrackCuts.h"
ebba9ef0 66#include "AliAnalysisMuonUtility.h"
35ac82fa 67
68
69/// \cond CLASSIMP
70ClassImp(AliAnalysisTaskMuonCuts) // Class implementation in ROOT context
71/// \endcond
72
73
74//________________________________________________________________________
75AliAnalysisTaskMuonCuts::AliAnalysisTaskMuonCuts() :
76 AliVAnalysisMuon(),
77 fHistoTypeKeys(0x0),
78 fThetaAbsKeys(0x0),
79 fSigmaCuts(TArrayD())
80{
81 /// Default ctor.
82}
83
84//________________________________________________________________________
85AliAnalysisTaskMuonCuts::AliAnalysisTaskMuonCuts(const char *name, const AliMuonTrackCuts& cuts) :
86 AliVAnalysisMuon(name, cuts),
87 fHistoTypeKeys(0x0),
88 fThetaAbsKeys(0x0),
89 fSigmaCuts(TArrayD())
90{
91 //
92 /// Constructor.
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//________________________________________________________________________
105AliAnalysisTaskMuonCuts::~AliAnalysisTaskMuonCuts()
106{
107 //
108 /// Destructor
109 //
110
111 delete fHistoTypeKeys;
112 delete fThetaAbsKeys;
113}
114
115//___________________________________________________________________________
116void 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);
064a0d13 147 histo = new TH2F(histoName.Data(), histoName.Data(), 100, 0., 800., 400, 0., 200.);
35ac82fa 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//________________________________________________________________________
183TString AliAnalysisTaskMuonCuts::GetHistoName(Int_t histoTypeIndex, Int_t thetaAbsIndex, Int_t srcIndex)
184{
185 /// Get local histogram name
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//________________________________________________________________________
192void AliAnalysisTaskMuonCuts::ProcessEvent(TString physSel, const TObjArray& selectTrigClasses, TString centrality)
193{
194 //
195 /// Fill histogram
196 //
197
35ac82fa 198 TString histoName = "";
199 AliVParticle* track = 0x0;
ebba9ef0 200 Int_t nTracks = AliAnalysisMuonUtility::GetNTracks(InputEvent());
35ac82fa 201 for (Int_t itrack = 0; itrack < nTracks; itrack++) {
ebba9ef0 202 track = AliAnalysisMuonUtility::GetTrack(itrack, InputEvent());
68e9b988 203 fMuonTrackCuts->CustomParam()->SetNSigmaPdca(1.e10);
35ac82fa 204 if ( ! fMuonTrackCuts->IsSelected(track) ) continue;
205
ebba9ef0 206 Double_t thetaAbsEndDeg = AliAnalysisMuonUtility::GetThetaAbsDeg(track);
35ac82fa 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
68e9b988 216 Double_t sigmaPdca = fMuonTrackCuts->IsThetaAbs23(track) ? fMuonTrackCuts->GetMuonTrackCutsParam().GetSigmaPdca23() : fMuonTrackCuts->GetMuonTrackCutsParam().GetSigmaPdca310();
217 Double_t chi2 = pDca / sigmaPdca;
35ac82fa 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) {
68e9b988 247 fMuonTrackCuts->CustomParam()->SetNSigmaPdca(fSigmaCuts[isigma]);
35ac82fa 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//________________________________________________________________________
261void AliAnalysisTaskMuonCuts::SetSigmaCuts(Int_t nSigmaCuts, Double_t* sigmaCuts)
262{
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);
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//________________________________________________________________________
282void AliAnalysisTaskMuonCuts::Terminate(Option_t *) {
283 //
284 /// Draw some histogram at the end.
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;
064a0d13 310
311 fMuonTrackCuts->Print("param");
35ac82fa 312
064a0d13 313 Int_t srcColors[kNtrackSources] = {kBlack, kRed, kSpring, kTeal, kBlue, kViolet, kMagenta, kOrange, kGray};
35ac82fa 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
321 //////////////
322 // Reco DCA //
323 //////////////
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 = "";
ebba9ef0 341 histoPattern = Form("%s%s*", fHistoTypeKeys->At(recoDcaHisto[ihisto])->GetName(), fThetaAbsKeys->At(itheta)->GetName());
35ac82fa 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
424 //////////////
425 // Fit pDCA //
426 //////////////
68e9b988 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;
35ac82fa 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.};
064a0d13 447 Double_t xMaxFit[2] = {260., 150.}; //{320., 150.}; // {360., 180.};
35ac82fa 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 ) {
c38835b3 479 currName = Form("hPDCA_%s_%s_%s_%s", fSrcKeys->At(isrc)->GetName(), physSel.Data(), trigClassName.Data(), fThetaAbsKeys->At(itheta)->GetName());
35ac82fa 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");
35ac82fa 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++;
064a0d13 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;
35ac82fa 556 currName = histoCheck->GetName();
064a0d13 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
35ac82fa 580
581
582 ///////////////////////////
583 // Check Chi square prob //
584 ///////////////////////////
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 ) {
c38835b3 604 currName = Form("hChiProb_%s_%s_%s_%s", fSrcKeys->At(isrc)->GetName(), physSel.Data(), trigClassName.Data(), fThetaAbsKeys->At(itheta)->GetName());
35ac82fa 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
626 //////////////////////
627 // Check sigma cuts //
628 //////////////////////
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};
9e6be1fd 634 TArrayI orderCuts;
35ac82fa 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 }
9e6be1fd 652 if ( orderCuts.GetSize() == 0 ) {
35ac82fa 653 // Re-order axis
654 TAxis* axis = histo->GetYaxis();
9e6be1fd 655 Int_t nSigmaCuts = ( useCustomSigma ) ? fSigmaCuts.GetSize() : axis->GetNbins();
656 orderCuts.Set(nSigmaCuts);
35ac82fa 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 = "";
9e6be1fd 695 for ( Int_t isigma=0; isigma<orderCuts.GetSize(); ++isigma ) {
35ac82fa 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();
064a0d13 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 }
35ac82fa 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;
064a0d13 761 Double_t ptMin[] = {0., 2., 4., 15., 60.};
35ac82fa 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);
ebba9ef0 769 for ( Int_t icent=1; icent<=GetCentralityClasses()->GetNbins(); ++icent ) {
770 TH2* histo = (TH2*)GetSum(physSel, trigClassName, GetCentralityClasses()->GetBinLabel(icent), histoName);
35ac82fa 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();
ebba9ef0 775 currName += Form("_%s_ptMin%i", GetCentralityClasses()->GetBinLabel(icent), TMath::Nint(ptMin[iptmin]));
35ac82fa 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());
ebba9ef0 799 leg->AddEntry(projectHisto, GetCentralityClasses()->GetBinLabel(icent), "lp");
35ac82fa 800 Double_t keptEvents = projectHisto->GetBinContent(orderCuts[0]);
9e6be1fd 801 Double_t totalEvents = projectHisto->GetBinContent(orderCuts[orderCuts.GetSize()-1]);
35ac82fa 802 Double_t accepted = ( totalEvents > 0. ) ? keptEvents / totalEvents : 1.;
803 Double_t acceptedErr = ( totalEvents > 0. ) ? TMath::Sqrt(accepted*(1.-accepted)/totalEvents) : 1.;
ebba9ef0 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.);
9e6be1fd 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]));
35ac82fa 806 } // loop on centrality
807 if ( leg ) leg->Draw("same");
808 } // loop on theta abs
809 } // loop on sources
810 } // loop on pt min
35ac82fa 811}