]> git.uio.no Git - u/mrichter/AliRoot.git/blame_incremental - PWG3/muon/AliAnalysisTaskSingleMu.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame_incremental
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (non-incremental) | history | HEAD
cleanup of projection and slice functions
[u/mrichter/AliRoot.git] / PWG3 / muon / AliAnalysisTaskSingleMu.cxx
... / ...
CommitLineData
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//-----------------------------------------------------------------------------
17/// \class AliAnalysisTaskSingleMu
18/// Analysis task for single muons in the spectrometer.
19/// The output is a list of histograms.
20/// The macro class can run on AODs or ESDs.
21/// In the latter case a flag can be activated to produce a tree as output.
22/// If Monte Carlo information is present, some basics checks are performed.
23///
24/// \author Diego Stocco
25//-----------------------------------------------------------------------------
26
27//----------------------------------------------------------------------------
28// Implementation of the single muon analysis class
29// An example of usage can be found in the macro RunSingleMuAnalysisFromAOD.C.
30//----------------------------------------------------------------------------
31
32#define AliAnalysisTaskSingleMu_cxx
33
34// ROOT includes
35#include "TROOT.h"
36#include "TH1.h"
37#include "TH2.h"
38#include "TH3.h"
39#include "TAxis.h"
40#include "TList.h"
41#include "TCanvas.h"
42#include "TMath.h"
43#include "TTree.h"
44#include "TTimeStamp.h"
45#include "TMap.h"
46#include "TObjString.h"
47#include "TObjArray.h"
48#include "TIterator.h"
49#include "TParameter.h"
50#include "TMCProcess.h"
51
52// STEER includes
53#include "AliLog.h"
54
55#include "AliAODInputHandler.h"
56#include "AliAODEvent.h"
57#include "AliAODTrack.h"
58#include "AliAODVertex.h"
59
60#include "AliMCParticle.h"
61#include "AliStack.h"
62
63#include "AliESDInputHandler.h"
64#include "AliESDEvent.h"
65#include "AliESDMuonTrack.h"
66
67#include "AliMultiplicity.h"
68
69// ANALYSIS includes
70#include "AliAnalysisManager.h"
71#include "AliAnalysisTaskSE.h"
72#include "AliAnalysisDataSlot.h"
73#include "AliAnalysisDataContainer.h"
74
75// CORRFW includes
76#include "AliCFManager.h"
77#include "AliCFContainer.h"
78#include "AliCFGridSparse.h"
79#include "AliCFTrackKineCuts.h"
80#include "AliCFParticleGenCuts.h"
81#include "AliCFEventRecCuts.h"
82
83#include "AliAnalysisTaskSingleMu.h"
84
85/// \cond CLASSIMP
86ClassImp(AliAnalysisTaskSingleMu) // Class implementation in ROOT context
87/// \endcond
88
89
90//________________________________________________________________________
91AliAnalysisTaskSingleMu::AliAnalysisTaskSingleMu(const char *name, Int_t fillTreeScaleDown, Bool_t keepAll) :
92 AliAnalysisTaskSE(name),
93 fFillTreeScaleDown(fillTreeScaleDown),
94 fKeepAll(keepAll),
95 fkNvtxContribCut(1),
96 fTriggerClasses(0x0),
97 fCFManager(0),
98 fHistoList(0),
99 fHistoListMC(0),
100 fHistoListQA(0),
101 fTreeSingleMu(0),
102 fVarFloat(0),
103 fVarInt(0),
104 fVarChar(0),
105 fVarUInt(0),
106 fVarFloatMC(0),
107 fVarIntMC(0),
108 fAuxObjects(new TMap())
109{
110 //
111 /// Constructor.
112 //
113 if ( fFillTreeScaleDown <= 0 )
114 fKeepAll = kFALSE;
115
116 DefineOutput(1, AliCFContainer::Class());
117 DefineOutput(2, TList::Class());
118 DefineOutput(3, TList::Class());
119 DefineOutput(4, TList::Class());
120
121 if ( fFillTreeScaleDown > 0 )
122 DefineOutput(5, TTree::Class());
123
124 fAuxObjects->SetOwner();
125
126 SetTriggerClasses();
127}
128
129
130//________________________________________________________________________
131AliAnalysisTaskSingleMu::~AliAnalysisTaskSingleMu()
132{
133 //
134 /// Destructor
135 //
136
137 delete fTriggerClasses;
138 // For proof: do not delete output containers
139 if ( ! AliAnalysisManager::GetAnalysisManager()->IsProofMode() ) {
140 delete fCFManager->GetParticleContainer(); // The container is not deleted by framework
141 delete fCFManager;
142 delete fHistoList;
143 delete fHistoListMC;
144 delete fHistoListQA;
145 delete fTreeSingleMu;
146 }
147 delete fVarFloat;
148 delete fVarInt;
149 delete [] fVarChar;
150 delete fVarUInt;
151 delete fVarFloatMC;
152 delete fVarIntMC;
153 delete fAuxObjects;
154}
155
156
157//___________________________________________________________________________
158void AliAnalysisTaskSingleMu::NotifyRun()
159{
160 //
161 /// Use the event handler information to correctly fill the analysis flags:
162 /// - check if Monte Carlo information is present
163 //
164
165 if ( fMCEvent ) {
166 AliInfo("Monte Carlo information is present");
167 }
168 else {
169 AliInfo("No Monte Carlo information in run");
170 }
171}
172
173
174//___________________________________________________________________________
175void AliAnalysisTaskSingleMu::FinishTaskOutput()
176{
177 //
178 /// Perform histogram normalization after the last analyzed event
179 /// but before merging
180 //
181
182 // Set the correct run limits for the histograms
183 // vs run number to cope with the merging
184 Int_t histoIndex = -1;
185 Int_t indexPerRun[2] = {kHistoNeventsPerRun, kHistoNmuonsPerRun};
186 for ( Int_t ihisto=0; ihisto<2; ihisto++) {
187 histoIndex = GetHistoIndex(indexPerRun[ihisto]);
188 TH2F* histo2D = (TH2F*)fHistoList->At(histoIndex);
189 Double_t minX = 1e10, maxX = -1e10;
190 for (Int_t ibin=1; ibin<=histo2D->GetXaxis()->GetNbins(); ibin++){
191 TString runNum = histo2D->GetXaxis()->GetBinLabel(ibin);
192 minX = TMath::Min(runNum.Atof()-0.5, minX);
193 maxX = TMath::Max(runNum.Atof()+0.5, maxX);
194 }
195 histo2D->GetXaxis()->SetLimits(minX, maxX);
196 AliInfo(Form("Histogram %s run limits (%f, %f)",histo2D->GetName(), minX, maxX));
197 }
198}
199
200
201//___________________________________________________________________________
202void AliAnalysisTaskSingleMu::UserCreateOutputObjects()
203{
204 //
205 /// Create output histograms
206 //
207 AliInfo(Form(" CreateOutputObjects of task %s\n", GetName()));
208
209 // initialize histogram lists
210 if ( ! fHistoList ) fHistoList = new TList();
211 if ( ! fHistoListMC ) fHistoListMC = new TList();
212 if ( ! fHistoListQA ) fHistoListQA = new TList();
213
214 // Init variables
215 fVarFloat = new Float_t [kNvarFloat];
216 fVarInt = new Int_t [kNvarInt];
217 fVarChar = new Char_t *[kNvarChar];
218 fVarUInt = new UInt_t [kNvarUInt];
219 fVarFloatMC = new Float_t [kNvarFloatMC];
220 fVarIntMC = new Int_t [kNvarIntMC];
221
222 const Int_t charWidth[kNvarChar] = {255};
223 for(Int_t ivar=0; ivar<kNvarChar; ivar++){
224 fVarChar[ivar] = new Char_t [charWidth[ivar]];
225 }
226
227 Int_t nPtBins = 60;
228 Double_t ptMin = 0., ptMax = 30.;
229 TString ptName("Pt"), ptTitle("p_{t}"), ptUnits("GeV/c");
230
231 Int_t nRapidityBins = 25;
232 Double_t rapidityMin = -4.5, rapidityMax = -2.;
233 //TString rapidityName("Rapidity"), rapidityTitle("y"), rapidityUnits("");
234 TString rapidityName("Eta"), rapidityTitle("#eta"), rapidityUnits("");
235
236 Int_t nPhiBins = 36;
237 Double_t phiMin = 0.; Double_t phiMax = 2*TMath::Pi();
238 TString phiName("Phi"), phiTitle("#phi"), phiUnits("rad");
239
240 Int_t nDcaBins = 30;
241 Double_t dcaMin = 0., dcaMax = 300.;
242 TString dcaName("DCA"), dcaTitle("DCA"), dcaUnits("cm");
243
244 Int_t nVzBins = 40;
245 Double_t vzMin = -20., vzMax = 20.;
246 TString vzName("Vz"), vzTitle("Vz"), vzUnits("cm");
247
248 Int_t nThetaAbsEndBins = 4;
249 Double_t thetaAbsEndMin = -0.5, thetaAbsEndMax = 3.5;
250 TString thetaAbsEndName("ThetaAbsEnd"), thetaAbsEndTitle("#theta_{abs}"), thetaAbsEndUnits("a.u.");
251
252 Int_t nChargeBins = 2;
253 Double_t chargeMin = -2, chargeMax = 2.;
254 TString chargeName("Charge"), chargeTitle("charge"), chargeUnits("e");
255
256 Int_t nMatchTrigBins = 4;
257 Double_t matchTrigMin = -0.5, matchTrigMax = 3.5;
258 TString matchTrigName("MatchTrig"), matchTrigTitle("Trigger match"), matchTrigUnits("");
259
260 Int_t nTrigClassBins = fTriggerClasses->GetEntries();
261 Double_t trigClassMin = 0.5, trigClassMax = (Double_t)nTrigClassBins + 0.5;
262 TString trigClassName("TrigClass"), trigClassTitle("Fired trigger class"), trigClassUnits("");
263
264 Int_t nGoodVtxBins = 3;
265 Double_t goodVtxMin = -0.5, goodVtxMax = 2.5;
266 TString goodVtxName("GoodVtx"), goodVtxTitle("Vertex flags"), goodVtxUnits("");
267
268 Int_t nMotherTypeBins = kNtrackSources;
269 Double_t motherTypeMin = -0.5, motherTypeMax = (Double_t)kNtrackSources - 0.5;
270 TString motherType("MotherType"), motherTypeTitle("motherType"), motherTypeUnits("");
271
272 Int_t nPBins = 100;
273 Double_t pMin = 0., pMax = 400.;
274 TString pName("P"), pTitle("p"), pUnits("GeV/c");
275
276 TString trigName[kNtrigCuts];
277 trigName[kNoMatchTrig] = "NoMatch";
278 trigName[kAllPtTrig] = "AllPt";
279 trigName[kLowPtTrig] = "LowPt";
280 trigName[kHighPtTrig] = "HighPt";
281
282 TString srcName[kNtrackSources];
283 srcName[kCharmMu] = "Charm";
284 srcName[kBeautyMu] = "Beauty";
285 srcName[kPrimaryMu] = "Decay";
286 srcName[kSecondaryMu] = "Secondary";
287 srcName[kRecoHadron] = "Hadrons";
288 srcName[kUnknownPart] = "Unknown";
289
290 TH1F* histo1D = 0x0;
291 TH2F* histo2D = 0x0;
292 TString histoName, histoTitle;
293 Int_t histoIndex = 0;
294
295 // Multi-dimensional histo
296 Int_t nbins[kNvars] = {nPtBins, nRapidityBins, nPhiBins, nDcaBins, nVzBins, nThetaAbsEndBins, nChargeBins, nMatchTrigBins, nTrigClassBins, nGoodVtxBins, nMotherTypeBins, nPBins};
297 Double_t xmin[kNvars] = {ptMin, rapidityMin, phiMin, dcaMin, vzMin, thetaAbsEndMin, chargeMin, matchTrigMin, trigClassMin, goodVtxMin, motherTypeMin, pMin};
298 Double_t xmax[kNvars] = {ptMax, rapidityMax, phiMax, dcaMax, vzMax, thetaAbsEndMax, chargeMax, matchTrigMax, trigClassMax, goodVtxMax, motherTypeMax, pMax};
299 TString axisTitle[kNvars] = {ptTitle, rapidityTitle, phiTitle, dcaTitle, vzTitle, thetaAbsEndTitle, chargeTitle, matchTrigTitle, trigClassTitle, goodVtxTitle, motherTypeTitle, pTitle};
300 TString axisUnits[kNvars] = {ptUnits, rapidityUnits, phiUnits, dcaUnits, vzUnits, thetaAbsEndUnits, chargeUnits, matchTrigUnits, trigClassUnits, goodVtxUnits, motherTypeUnits, pUnits};
301
302 //TString stepTitle[kNsteps] = {"reconstructed", "in acceptance", "generated", "in acceptance (MC)"};
303 TString stepTitle[kNsteps] = {"reconstructed", "generated"};
304
305 // Create CF container
306
307 // The framework has problems if the name of the object
308 // and the one of container differ
309 // To be complaint, get the name from container and set it
310 TString containerName = GetOutputSlot(1)->GetContainer()->GetName();
311
312 AliCFContainer* container = new AliCFContainer(containerName.Data(),"container for tracks",kNsteps,kNvars,nbins);
313
314 for ( Int_t idim = 0; idim<kNvars; idim++){
315 histoTitle = Form("%s (%s)", axisTitle[idim].Data(), axisUnits[idim].Data());
316 histoTitle.ReplaceAll("()","");
317
318 container->SetVarTitle(idim, histoTitle.Data());
319 container->SetBinLimits(idim, xmin[idim], xmax[idim]);
320 }
321
322 for (Int_t istep=0; istep<kNsteps; istep++){
323 container->SetStepTitle(istep, stepTitle[istep].Data());
324 AliCFGridSparse* gridSparse = container->GetGrid(istep);
325
326 SetAxisLabel(gridSparse->GetAxis(kHvarTrigClass));
327 TAxis* isGoodVtxAxis = gridSparse->GetAxis(kHvarIsGoodVtx);
328 isGoodVtxAxis->SetBinLabel(1,Form("No vertex (contrib<%i)",fkNvtxContribCut));
329 isGoodVtxAxis->SetBinLabel(2,"Good vertex");
330 isGoodVtxAxis->SetBinLabel(3,"Pileup");
331 TAxis* motherTypeAxis = gridSparse->GetAxis(kHvarMotherType);
332 for (Int_t ibin=0; ibin<kNtrackSources; ibin++){
333 motherTypeAxis->SetBinLabel(ibin+1,srcName[ibin]);
334 }
335 }
336
337 // Create cuts
338
339 //Particle-Level cuts:
340 AliCFParticleGenCuts* mcGenCuts = new AliCFParticleGenCuts();
341 mcGenCuts->SetNameTitle("mcGenCuts","MC particle generation cuts");
342 mcGenCuts->SetRequirePdgCode(13,kTRUE);
343 mcGenCuts->SetQAOn(fHistoListQA);
344
345 /*
346 // MC kinematic cuts
347 AliCFTrackKineCuts *mcAccCuts = new AliCFTrackKineCuts();
348 mcAccCuts->SetNameTitle("mcAccCuts","MC-level acceptance cuts");
349 mcAccCuts->SetEtaRange(-4.,-2.5);
350 mcAccCuts->SetQAOn(fHistoListQA);
351
352 // Rec-Level kinematic cuts
353 AliCFTrackKineCuts *recAccCuts = new AliCFTrackKineCuts();
354 recAccCuts->SetNameTitle("recAccCuts","Reco-level acceptance cuts");
355 recAccCuts->SetEtaRange(-4.,-2.5);
356 recAccCuts->SetQAOn(fHistoListQA);
357 */
358
359 TObjArray* mcGenList = new TObjArray(0) ;
360 mcGenList->AddLast(mcGenCuts);
361
362 /*
363 TObjArray* mcAccList = new TObjArray(0);
364 mcAccList->AddLast(mcAccCuts);
365
366 TObjArray* recAccList = new TObjArray(0);
367 recAccList->AddLast(recAccCuts);
368 */
369
370 // Create CF manager
371 fCFManager = new AliCFManager() ;
372 fCFManager->SetParticleContainer(container);
373
374 // Add cuts
375 // Dummy event container
376 Int_t dummyBins[1] = {1};
377 AliCFContainer* evtContainer = new AliCFContainer("dummyContainer","dummy contaier for events",1,1,dummyBins);
378 fCFManager->SetEventContainer(evtContainer);
379 fCFManager->SetEventCutsList(0,0x0);
380
381 // Init empty cuts (avoid warnings in framework)
382 for (Int_t istep=0; istep<kNsteps; istep++) {
383 fCFManager->SetParticleCutsList(istep,0x0);
384 }
385 //fCFManager->SetParticleCutsList(kStepAcceptance,recAccList);
386 fCFManager->SetParticleCutsList(kStepGeneratedMC,mcGenList);
387 //fCFManager->SetParticleCutsList(kStepAcceptanceMC,mcAccList);
388
389 // Summary histos
390 histoName = "histoNeventsPerTrig";
391 histoTitle = "Number of events per trigger class";
392 histo2D = new TH2F(histoName.Data(), histoTitle.Data(), nTrigClassBins, trigClassMin, trigClassMax, 5, 0.5, 5.5);
393 histo2D->GetXaxis()->SetTitle("Fired class");
394 SetAxisLabel(histo2D->GetXaxis());
395 histo2D->GetYaxis()->SetBinLabel(1, "All events");
396 histo2D->GetYaxis()->SetBinLabel(2, "Good vtx events");
397 histo2D->GetYaxis()->SetBinLabel(3, "Pileup events");
398 histo2D->GetYaxis()->SetBinLabel(4, "All vertices");
399 histo2D->GetYaxis()->SetBinLabel(5, "Pileup vertices");
400 histoIndex = GetHistoIndex(kHistoNeventsPerTrig);
401 fHistoList->AddAt(histo2D, histoIndex);
402
403 histoName = "histoMuonMultiplicity";
404 histoTitle = "Muon track multiplicity";
405 histo2D = new TH2F(histoName.Data(), histoTitle.Data(), 15, -0.5, 15-0.5,
406 nTrigClassBins, trigClassMin, trigClassMax);
407 histo2D->GetXaxis()->SetTitle("# of muons");
408 SetAxisLabel(histo2D->GetYaxis());
409 histoIndex = GetHistoIndex(kHistoMuonMultiplicity);
410 fHistoList->AddAt(histo2D, histoIndex);
411
412 histoName = "histoEventVz";
413 histoTitle = "All events IP Vz distribution";
414 histo2D = new TH2F(histoName.Data(), histoTitle.Data(), nTrigClassBins, trigClassMin, trigClassMax, nVzBins, vzMin, vzMax);
415 histo2D->GetXaxis()->SetTitle("Fired class");
416 SetAxisLabel(histo2D->GetXaxis());
417 histoTitle = Form("%s (%s)", vzTitle.Data(), vzUnits.Data());
418 histo2D->GetYaxis()->SetTitle(histoTitle.Data());
419 histoIndex = GetHistoIndex(kHistoEventVz);
420 fHistoList->AddAt(histo2D, histoIndex);
421
422 Int_t hRunIndex[2] = {kHistoNeventsPerRun, kHistoNmuonsPerRun};
423 TString hRunName[2] = {"histoNeventsPerRun", "histoNmuonsPerRun"};
424 TString hRunTitle[2] = {"Number of events per run", "Number of muons per run"};
425 for ( Int_t ihisto=0; ihisto<2; ihisto++){
426 histoName = hRunName[ihisto];
427 histoTitle = hRunTitle[ihisto];
428 histo2D = new TH2F(histoName.Data(), histoTitle.Data(),
429 1, 0., 0.,
430 nTrigClassBins, trigClassMin, trigClassMax);
431 histo2D->GetXaxis()->SetTitle("Run number");
432 SetAxisLabel(histo2D->GetYaxis());
433 histo2D->Sumw2();
434 histoIndex = hRunIndex[ihisto];
435 fHistoList->AddAt(histo2D, histoIndex);
436 }
437
438 // MC histos summary
439 histoName = "histoCheckVz";
440 histoTitle = "Check IP Vz distribution";
441 histo2D = new TH2F(histoName.Data(), histoTitle.Data(),
442 nVzBins, vzMin, vzMax,
443 nVzBins, vzMin, vzMax);
444 histoTitle = Form("%s (%s)", vzTitle.Data(), vzUnits.Data());
445 histo2D->GetXaxis()->SetTitle(histoTitle.Data());
446 histoTitle = Form("%s MC (%s)", vzTitle.Data(), vzUnits.Data());
447 histo2D->GetYaxis()->SetTitle(histoTitle.Data());
448 histoIndex = GetHistoIndex(kHistoCheckVzMC);
449 fHistoListMC->AddAt(histo2D, histoIndex);
450
451 // MC histos
452 for (Int_t itrig=0; itrig < kNtrigCuts; itrig++) {
453 for (Int_t isrc = 0; isrc < kNtrackSources; isrc++) {
454 histoName = Form("%sResolution%sTrig%s", ptName.Data(), trigName[itrig].Data(), srcName[isrc].Data());
455 histoTitle = Form("%s resolution. Trigger: %s (%s)", ptTitle.Data(), trigName[itrig].Data(), srcName[isrc].Data());
456 histo1D = new TH1F(histoName.Data(), histoTitle.Data(), 100, -5., 5.);
457 histo1D->GetXaxis()->SetTitle(Form("%s^{reco} - %s^{MC} (%s)", ptTitle.Data(), ptTitle.Data(), ptUnits.Data()));
458 histoIndex = GetHistoIndex(kHistoPtResolutionMC, itrig, isrc);
459 fHistoListMC->AddAt(histo1D, histoIndex);
460 }
461 }
462
463 // Trees
464 if ( fFillTreeScaleDown > 0 ) {
465 TString leavesFloat[kNvarFloat] = {"Px", "Py", "Pz", "Pt",
466 "PxAtDCA", "PyAtDCA", "PzAtDCA", "PtAtDCA",
467 "PxUncorrected", "PyUncorrected", "PzUncorrected", "PtUncorrected",
468 "XUncorrected", "YUncorrected", "ZUncorrected",
469 "XatDCA", "YatDCA", "DCA",
470 "Eta", "Rapidity", "Charge", "RAtAbsEnd",
471 "IPVx", "IPVy", "IPVz"};
472 TString leavesInt[kNvarInt] = {"MatchTrig", "IsMuon", "IsGhost", "LoCircuit", "PassPhysicsSelection", "NVtxContrib", "NspdTracklets", "IsPileupVertex"};
473 TString leavesChar[kNvarChar] = {"FiredTrigClass"};
474 TString leavesUInt[kNvarUInt] = {"BunchCrossNum", "OrbitNum", "PeriodNum", "RunNum"};
475 TString leavesFloatMC[kNvarFloatMC] = {"PxMC", "PyMC", "PzMC", "PtMC", "EtaMC", "RapidityMC", "VxMC", "VyMC", "VzMC"};
476 TString leavesIntMC[kNvarIntMC] = {"Pdg", "MotherType"};
477
478 TString treeName = GetOutputSlot(5)->GetContainer()->GetName();
479 Bool_t hasMC = AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler();
480 TString treeTitle = ( hasMC ) ? "Single Mu" : "Single Mu MC";
481
482 OpenFile(5);
483 if ( ! fTreeSingleMu ) fTreeSingleMu = new TTree(treeName.Data(), treeTitle.Data());
484
485 for(Int_t itree=0; itree<2; itree++){
486 TParameter<Int_t>* par1 = new TParameter<Int_t>("fillTreeScaleDown",fFillTreeScaleDown);
487 TParameter<Int_t>* par2 = new TParameter<Int_t>("keepAllEvents",fKeepAll);
488 fTreeSingleMu->GetUserInfo()->Add(par1);
489 fTreeSingleMu->GetUserInfo()->Add(par2);
490 for(Int_t ivar=0; ivar<kNvarFloat; ivar++){
491 fTreeSingleMu->Branch(leavesFloat[ivar].Data(), &fVarFloat[ivar], Form("%s/F", leavesFloat[ivar].Data()));
492 }
493 for(Int_t ivar=0; ivar<kNvarInt; ivar++){
494 fTreeSingleMu->Branch(leavesInt[ivar].Data(), &fVarInt[ivar], Form("%s/I", leavesInt[ivar].Data()));
495 }
496 for(Int_t ivar=0; ivar<kNvarChar; ivar++){
497 TString addString = leavesChar[ivar] + "/C";
498 fTreeSingleMu->Branch(leavesChar[ivar].Data(), fVarChar[ivar], addString.Data());
499 }
500 for(Int_t ivar=0; ivar<kNvarUInt; ivar++){
501 fTreeSingleMu->Branch(leavesUInt[ivar].Data(), &fVarUInt[ivar], Form("%s/i", leavesUInt[ivar].Data()));
502 }
503 if ( hasMC ){
504 for(Int_t ivar=0; ivar<kNvarFloatMC; ivar++){
505 fTreeSingleMu->Branch(leavesFloatMC[ivar].Data(), &fVarFloatMC[ivar], Form("%s/F", leavesFloatMC[ivar].Data()));
506 }
507 for(Int_t ivar=0; ivar<kNvarIntMC; ivar++){
508 fTreeSingleMu->Branch(leavesIntMC[ivar].Data(), &fVarIntMC[ivar], Form("%s/I", leavesIntMC[ivar].Data()));
509 }
510 }
511 } // loop on trees
512 PostData(5, fTreeSingleMu);
513 } // fillNTuple
514
515 PostData(1, fCFManager->GetParticleContainer());
516 PostData(2, fHistoList);
517 PostData(3, fHistoListQA);
518 PostData(4, fHistoListMC);
519
520}
521
522//________________________________________________________________________
523void AliAnalysisTaskSingleMu::UserExec(Option_t * /*option*/)
524{
525 //
526 /// Main loop
527 /// Called for each event
528 //
529
530 AliESDEvent* esdEvent = 0x0;
531 AliAODEvent* aodEvent = 0x0;
532
533 esdEvent = dynamic_cast<AliESDEvent*> (InputEvent());
534 if ( ! esdEvent ){
535 aodEvent = dynamic_cast<AliAODEvent*> (InputEvent());
536 }
537
538 if ( ! aodEvent && ! esdEvent ) {
539 AliError ("AOD or ESD event not found. Nothing done!");
540 return;
541 }
542
543 if ( ! fMCEvent && InputEvent()->GetEventType() != 7 ) return; // Run only on physics events!
544
545 Bool_t fillCurrentEventTree = ( fFillTreeScaleDown == 0 ) ? kFALSE : ( Entry() % fFillTreeScaleDown == 0 );
546
547 Reset(kFALSE);
548
549 // Global event info
550 TString firedTrigClasses = ( esdEvent ) ? esdEvent->GetFiredTriggerClasses() : aodEvent->GetFiredTriggerClasses();
551 /*
552 if ( fMCEvent ) {
553 // Add the MB name (which is not there in simulation)
554 // CAVEAT: to be checked if we perform beam-gas simulations
555 //if ( ! firedTrigClasses.Contains("MB") ) printf("fired class %s\n", firedTrigClasses.Data()); // REMEMBER TO CUT
556 if ( ! firedTrigClasses.Contains("CINT") && ! firedTrigClasses.Contains("MB") ) {
557 TString trigName = ((TObjString*)fTriggerClasses->At(0))->GetString();
558 firedTrigClasses.Prepend(Form("%s ", trigName.Data()));
559 }
560 }
561 */
562
563 fVarFloat[kVarIPVz] = ( esdEvent ) ? esdEvent->GetPrimaryVertex()->GetZ() : aodEvent->GetPrimaryVertex()->GetZ();
564 fVarInt[kVarNVtxContrib] = ( esdEvent ) ? esdEvent->GetPrimaryVertex()->GetNContributors() : aodEvent->GetPrimaryVertex()->GetNContributors();
565 fVarInt[kVarNspdTracklets] = ( esdEvent ) ? esdEvent->GetMultiplicity()->GetNumberOfTracklets() : aodEvent->GetTracklets()->GetNumberOfTracklets();
566 fVarInt[kVarIsPileup] = ( esdEvent ) ? esdEvent->IsPileupFromSPD(3,0.8) : aodEvent->IsPileupFromSPD(3,0.8); // REMEMBER TO CHECK
567
568 fVarUInt[kVarRunNumber] = ( esdEvent ) ? esdEvent->GetRunNumber() : aodEvent->GetRunNumber();
569
570 Int_t isGoodVtxBin = ( fVarInt[kVarNVtxContrib] >= fkNvtxContribCut );
571 if ( fVarInt[kVarIsPileup] > 0 )
572 isGoodVtxBin = 2;
573
574 if ( fillCurrentEventTree ){
575 strncpy(fVarChar[kVarTrigMask], firedTrigClasses.Data(),255);
576
577 fVarInt[kVarPassPhysicsSelection] = ( esdEvent ) ? ((AliESDInputHandler*)AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler())->IsEventSelected() : ((AliAODInputHandler*)AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler())->IsEventSelected();
578
579 // Small workaround: in MC the bunch ID are not properly set and the timestamp is in seconds
580 // So fill bunchCrossing with the read timestamp
581 // fill the orbit and period number with a timestamp created while reading the run
582 TTimeStamp ts;
583 fVarUInt[kVarBunchCrossNumber] = ( fMCEvent ) ? (UInt_t)Entry() : esdEvent->GetBunchCrossNumber();
584 fVarUInt[kVarOrbitNumber] = ( fMCEvent ) ? ts.GetNanoSec() : esdEvent->GetOrbitNumber();
585 fVarUInt[kVarPeriodNumber] = ( fMCEvent ) ? ts.GetTime() : esdEvent->GetPeriodNumber();
586
587 fVarFloat[kVarIPVx] = ( esdEvent ) ? esdEvent->GetPrimaryVertex()->GetX() : aodEvent->GetPrimaryVertex()->GetX();
588 fVarFloat[kVarIPVy] = ( esdEvent ) ? esdEvent->GetPrimaryVertex()->GetY() : aodEvent->GetPrimaryVertex()->GetY();
589 }
590
591 // Object declaration
592 AliMCParticle* mcPart = 0x0;
593 AliVParticle* track = 0x0;
594
595 Double_t containerInput[kNvars];
596 Int_t histoIndex = -1;
597
598 fCFManager->SetRecEventInfo(InputEvent());
599
600 // Pure Monte Carlo part
601 if ( fMCEvent ) {
602 fCFManager->SetMCEventInfo (fMCEvent);
603 Int_t nMCtracks = fMCEvent->GetNumberOfTracks();
604 if ( nMCtracks > 0 ) {
605 containerInput[kHvarVz] = fMCEvent->Stack()->Particle(0)->Vz();
606 containerInput[kHvarIsGoodVtx] = 1;
607 histoIndex = GetHistoIndex(kHistoCheckVzMC);
608 ((TH2F*)fHistoListMC->At(histoIndex))->Fill(fVarFloat[kVarIPVz], fMCEvent->Stack()->Particle(0)->Vz());
609
610 }
611
612 for (Int_t ipart=0; ipart<nMCtracks; ipart++) {
613 mcPart = (AliMCParticle*)fMCEvent->GetTrack(ipart);
614
615 //check the MC-level cuts
616 if ( ! fCFManager->CheckParticleCuts(kStepGeneratedMC,mcPart) )
617 continue;
618
619 containerInput[kHvarPt] = mcPart->Pt();
620 //containerInput[kHvarY] = mcPart->Y();
621 containerInput[kHvarEta] = mcPart->Eta();
622 containerInput[kHvarPhi] = mcPart->Phi();
623 containerInput[kHvarDCA] = TMath::Sqrt(mcPart->Xv()*mcPart->Xv() +
624 mcPart->Yv()*mcPart->Yv());
625 containerInput[kHvarThetaZones] = GetBinThetaAbsEnd(TMath::Pi()-mcPart->Theta(),kTRUE);
626 containerInput[kHvarCharge] = mcPart->Charge()/3.;
627 containerInput[kHvarMatchTrig] = 1.;
628 containerInput[kHvarMotherType] = (Double_t)RecoTrackMother(mcPart);
629 containerInput[kHvarP] = mcPart->P();
630
631 for ( Int_t itrig=0; itrig<fTriggerClasses->GetEntries(); itrig++ ) {
632 TString trigName = ((TObjString*)fTriggerClasses->At(itrig))->GetString();
633 if ( ! firedTrigClasses.Contains(trigName.Data()) )
634 continue;
635 containerInput[kHvarTrigClass] = (Double_t)(itrig+1);
636 fCFManager->GetParticleContainer()->Fill(containerInput,kStepGeneratedMC);
637 // if ( fCFManager->CheckParticleCuts(kStepAcceptanceMC,mcPart) ) fCFManager->GetParticleContainer()->Fill(containerInput,kStepAcceptanceMC);
638 } // loop on trigger classes
639 } // loop on MC particles
640 } // is MC
641
642
643 Int_t trackLabel = -1;
644 Bool_t isGhost = kFALSE;
645 Int_t nGhosts = 0, nMuons = 0;
646
647 Int_t nTracks = ( esdEvent ) ? esdEvent->GetNumberOfMuonTracks() : aodEvent->GetNTracks();
648
649 for (Int_t itrack = 0; itrack < nTracks; itrack++) {
650
651 // Get variables
652 if ( esdEvent ){
653 // ESD only info
654
655 track = esdEvent->GetMuonTrack(itrack);
656 isGhost = ( ((AliESDMuonTrack*)track)->ContainTriggerData() && ! ((AliESDMuonTrack*)track)->ContainTrackerData() );
657
658 if ( fillCurrentEventTree ) {
659 if ( itrack > 0 ) Reset(kTRUE);
660 fVarFloat[kVarPxUncorrected] = ( isGhost ) ? -TMath::Tan(((AliESDMuonTrack*)track)->GetThetaXUncorrected()) : ((AliESDMuonTrack*)track)->PxUncorrected();
661 fVarFloat[kVarPyUncorrected] = ( isGhost ) ? -TMath::Tan(((AliESDMuonTrack*)track)->GetThetaYUncorrected()) : ((AliESDMuonTrack*)track)->PyUncorrected();
662 fVarFloat[kVarPzUncorrected] = ( isGhost ) ? -1 : ((AliESDMuonTrack*)track)->PzUncorrected();
663
664 fVarFloat[kVarPtUncorrected] =
665 TMath::Sqrt(fVarFloat[kVarPxUncorrected] * fVarFloat[kVarPxUncorrected] +
666 fVarFloat[kVarPyUncorrected] * fVarFloat[kVarPyUncorrected]);
667
668 fVarFloat[kVarXUncorrected] = ((AliESDMuonTrack*)track)->GetNonBendingCoorUncorrected();
669 fVarFloat[kVarYUncorrected] = ((AliESDMuonTrack*)track)->GetBendingCoorUncorrected();
670 fVarFloat[kVarZUncorrected] = ((AliESDMuonTrack*)track)->GetZUncorrected();
671
672 fVarInt[kVarLocalCircuit] = ((AliESDMuonTrack*)track)->LoCircuit();
673 }
674
675 if ( isGhost ) {
676 // If is ghost fill only a partial information
677 nGhosts++;
678 if ( fillCurrentEventTree ) {
679 fVarInt[kVarIsMuon] = 0;
680 fVarInt[kVarIsGhost] = nGhosts;
681 fVarInt[kVarMatchTrig] = ((AliESDMuonTrack*)track)->GetMatchTrigger();
682 fTreeSingleMu->Fill();
683 }
684 continue;
685 }
686 }
687 else {
688 track = aodEvent->GetTrack(itrack);
689 if ( ! ((AliAODTrack*)track)->IsMuonTrack() )
690 continue;
691 }
692
693 // Information for tracks in tracker
694 nMuons++;
695
696 fVarFloat[kVarPt] = track->Pt();
697 //fVarFloat[kVarRapidity] = track->Y();
698 fVarFloat[kVarEta] = track->Eta();
699 fVarFloat[kVarXatDCA] = (esdEvent ) ? ((AliESDMuonTrack*)track)->GetNonBendingCoorAtDCA() : ((AliAODTrack*)track)->XAtDCA();
700 fVarFloat[kVarYatDCA] = (esdEvent ) ? ((AliESDMuonTrack*)track)->GetBendingCoorAtDCA() : ((AliAODTrack*)track)->YAtDCA();
701 fVarFloat[kVarDCA] =
702 TMath::Sqrt( fVarFloat[kVarXatDCA] * fVarFloat[kVarXatDCA] +
703 fVarFloat[kVarYatDCA] * fVarFloat[kVarYatDCA] );
704 fVarFloat[kVarCharge] = (Float_t)track->Charge();
705 fVarFloat[kVarRAtAbsEnd] = (esdEvent ) ? ((AliESDMuonTrack*)track)->GetRAtAbsorberEnd() : ((AliAODTrack*)track)->GetRAtAbsorberEnd();
706 fVarInt[kVarMatchTrig] = (esdEvent ) ? ((AliESDMuonTrack*)track)->GetMatchTrigger() : ((AliAODTrack*)track)->GetMatchTrigger();
707
708 if ( fillCurrentEventTree ){
709 fVarInt[kVarIsMuon] = nMuons;
710 fVarInt[kVarIsGhost] = 0;
711 //fVarFloat[kVarEta] = track->Eta();
712 fVarFloat[kVarRapidity] = track->Y();
713 fVarFloat[kVarPx] = track->Px();
714 fVarFloat[kVarPy] = track->Py();
715 fVarFloat[kVarPz] = track->Pz();
716 fVarFloat[kVarPxAtDCA] = ( esdEvent ) ? ((AliESDMuonTrack*)track)->PxAtDCA() : ((AliAODTrack*)track)->PxAtDCA();
717 fVarFloat[kVarPyAtDCA] = ( esdEvent ) ? ((AliESDMuonTrack*)track)->PyAtDCA() : ((AliAODTrack*)track)->PyAtDCA();
718 fVarFloat[kVarPzAtDCA] = ( esdEvent ) ? ((AliESDMuonTrack*)track)->PzAtDCA() : ((AliAODTrack*)track)->PzAtDCA();
719 fVarFloat[kVarPtAtDCA] = TMath::Sqrt(fVarFloat[kVarPxAtDCA]*fVarFloat[kVarPxAtDCA] + fVarFloat[kVarPyAtDCA]*fVarFloat[kVarPyAtDCA]);
720 }
721
722 fVarIntMC[kVarMotherType] = kUnknownPart;
723
724 // Monte Carlo part
725 if ( fMCEvent ) {
726
727 trackLabel = track->GetLabel();
728
729 AliMCParticle* matchedMCTrack = 0x0;
730
731 Int_t nMCtracks = fMCEvent->GetNumberOfTracks();
732 for(Int_t imctrack=0; imctrack<nMCtracks; imctrack++){
733 mcPart = (AliMCParticle*)fMCEvent->GetTrack(imctrack);
734 if ( trackLabel == mcPart->GetLabel() ) {
735 matchedMCTrack = mcPart;
736 break;
737 }
738 } // loop on MC tracks
739
740 if ( matchedMCTrack ) {
741 fVarIntMC[kVarMotherType] = RecoTrackMother(matchedMCTrack);
742 fVarFloatMC[kVarPtMC] = matchedMCTrack->Pt();
743 FillTriggerHistos(kHistoPtResolutionMC, fVarInt[kVarMatchTrig], fVarIntMC[kVarMotherType], fVarFloat[kVarPt] - fVarFloatMC[kVarPtMC]);
744 if ( fillCurrentEventTree ){
745 fVarFloatMC[kVarPxMC] = matchedMCTrack->Px();
746 fVarFloatMC[kVarPyMC] = matchedMCTrack->Py();
747 fVarFloatMC[kVarPzMC] = matchedMCTrack->Pz();
748 fVarFloatMC[kVarEtaMC] = matchedMCTrack->Eta();
749 fVarFloatMC[kVarRapidityMC] = matchedMCTrack->Y();
750 fVarFloatMC[kVarVxMC] = matchedMCTrack->Xv();
751 fVarFloatMC[kVarVyMC] = matchedMCTrack->Yv();
752 fVarFloatMC[kVarVzMC] = matchedMCTrack->Zv();
753 fVarIntMC[kVarPdg] = matchedMCTrack->PdgCode();
754 }
755 }
756 AliDebug(1, Form("Found mother %i", fVarIntMC[kVarMotherType]));
757 } // if use MC
758
759 containerInput[kHvarPt] = fVarFloat[kVarPt];
760 //containerInput[kHvarY] = fVarFloat[kVarRapidity];
761 containerInput[kHvarEta] = fVarFloat[kVarEta];
762 containerInput[kHvarPhi] = track->Phi();
763 containerInput[kHvarDCA] = fVarFloat[kVarDCA];
764 containerInput[kHvarVz] = fVarFloat[kVarIPVz];
765 containerInput[kHvarThetaZones] = GetBinThetaAbsEnd(fVarFloat[kVarRAtAbsEnd]);
766 containerInput[kHvarCharge] = fVarFloat[kVarCharge];
767 containerInput[kHvarMatchTrig] = (Double_t)fVarInt[kVarMatchTrig];
768 containerInput[kHvarIsGoodVtx] = (Double_t)isGoodVtxBin;
769 containerInput[kHvarMotherType] = (Double_t)fVarIntMC[kVarMotherType];
770 containerInput[kHvarP] = track->P();
771
772 histoIndex = GetHistoIndex(kHistoNmuonsPerRun);
773 for ( Int_t itrig=0; itrig<fTriggerClasses->GetEntries(); itrig++ ) {
774 TString trigName = ((TObjString*)fTriggerClasses->At(itrig))->GetString();
775 if ( ! firedTrigClasses.Contains(trigName.Data()) )
776 continue;
777 containerInput[kHvarTrigClass] = (Double_t)(itrig+1);
778 fCFManager->GetParticleContainer()->Fill(containerInput,kStepReconstructed);
779 // if ( fCFManager->CheckParticleCuts(kStepAcceptance,track) ) fCFManager->GetParticleContainer()->Fill(containerInput,kStepAcceptance);
780 ((TH2F*)fHistoList->At(histoIndex))->Fill(Form("%u",fVarUInt[kVarRunNumber]), containerInput[kHvarTrigClass], 1.);
781 } // loop on trigger classes
782
783 if ( fillCurrentEventTree ) fTreeSingleMu->Fill();
784 } // loop on tracks
785
786 if ( fillCurrentEventTree && fKeepAll && ( ( nMuons + nGhosts ) == 0 ) ) {
787 // Fill event also if there is not muon (when explicitely required)
788 fTreeSingleMu->Fill();
789 }
790
791 for ( Int_t itrig=0; itrig<fTriggerClasses->GetEntries(); itrig++ ) {
792 TString trigName = ((TObjString*)fTriggerClasses->At(itrig))->GetString();
793 if ( ! firedTrigClasses.Contains(trigName.Data()) )
794 continue;
795 Double_t trigClassBin = (Double_t)(itrig+1);
796 histoIndex = GetHistoIndex(kHistoNeventsPerTrig);
797 ((TH2F*)fHistoList->At(histoIndex))->Fill(trigClassBin, 1., 1.); // All events
798 ((TH2F*)fHistoList->At(histoIndex))->Fill(trigClassBin, 4., (Double_t)(fVarInt[kVarIsPileup]+1)); // All vertices
799 if ( isGoodVtxBin == 1 )
800 ((TH2F*)fHistoList->At(histoIndex))->Fill(trigClassBin, 2., 1.); // Good vtx events
801 else if ( isGoodVtxBin == 2 ) {
802 ((TH2F*)fHistoList->At(histoIndex))->Fill(trigClassBin, 3., 1.); // Pileup events
803 ((TH2F*)fHistoList->At(histoIndex))->Fill(trigClassBin, 5., 1.); // Pileup vertices
804 }
805
806 histoIndex = GetHistoIndex(kHistoMuonMultiplicity);
807 ((TH1F*)fHistoList->At(histoIndex))->Fill(nMuons, trigClassBin);
808
809 if ( isGoodVtxBin == 1 ) {
810 histoIndex = GetHistoIndex(kHistoEventVz);
811 ((TH2F*)fHistoList->At(histoIndex))->Fill(trigClassBin,fVarFloat[kVarIPVz]);
812 }
813
814 histoIndex = GetHistoIndex(kHistoNeventsPerRun);
815 ((TH2F*)fHistoList->At(histoIndex))->Fill(Form("%u",fVarUInt[kVarRunNumber]), trigClassBin, 1.);
816 } // loop on trigger classes
817
818
819 // Post final data. It will be written to a file with option "RECREATE"
820 PostData(1, fCFManager->GetParticleContainer());
821 PostData(2, fHistoList);
822 PostData(3, fHistoListQA);
823 PostData(4, fHistoListMC);
824 if ( fFillTreeScaleDown > 0 )
825 PostData(5, fTreeSingleMu);
826}
827
828//________________________________________________________________________
829void AliAnalysisTaskSingleMu::Terminate(Option_t *) {
830 //
831 /// Draw some histogram at the end.
832 //
833
834 AliCFContainer* container = dynamic_cast<AliCFContainer*> (GetOutputData(1));
835 if ( ! container ) {
836 AliError("Cannot find container in file");
837 return;
838 }
839
840 //Int_t histoIndex = -1;
841
842 if ( ! gROOT->IsBatch() ) {
843 TString currName = GetName();
844 currName.Prepend("c1_");
845 TCanvas *c1_SingleMu = new TCanvas(currName.Data(),"Vz vs Pt",10,10,310,310);
846 c1_SingleMu->SetFillColor(10); c1_SingleMu->SetHighLightColor(10);
847 c1_SingleMu->SetLeftMargin(0.15); c1_SingleMu->SetBottomMargin(0.15);
848 TH2* histo = container->Project(kHvarPt,kHvarVz,kStepReconstructed);
849 currName = GetName();
850 currName.Prepend("hPtVz_");
851 histo->SetName(currName.Data());
852 histo->Draw("COLZ");
853 }
854
855}
856
857//________________________________________________________________________
858 void AliAnalysisTaskSingleMu::FillTriggerHistos(Int_t histoIndex, Int_t matchTrig, Int_t motherType,
859 Float_t var1, Float_t var2, Float_t var3)
860{
861 //
862 /// Fill all histograms passing the trigger cuts
863 //
864
865 Int_t nTrigs = TMath::Max(1, matchTrig);
866 TArrayI trigToFill(nTrigs);
867 trigToFill[0] = matchTrig;
868 for(Int_t itrig = 1; itrig < matchTrig; itrig++){
869 trigToFill[itrig] = itrig;
870 }
871
872 TList* histoList = (motherType < 0 ) ? fHistoList : fHistoListMC;
873
874 TString className;
875 for(Int_t itrig = 0; itrig < nTrigs; itrig++){
876 Int_t currIndex = GetHistoIndex(histoIndex, trigToFill[itrig], motherType);
877 className = histoList->At(currIndex)->ClassName();
878 AliDebug(3, Form("matchTrig %i Fill %s", matchTrig, histoList->At(currIndex)->GetName()));
879 if (className.Contains("1"))
880 ((TH1F*)histoList->At(currIndex))->Fill(var1);
881 else if (className.Contains("2"))
882 ((TH2F*)histoList->At(currIndex))->Fill(var1, var2);
883 else if (className.Contains("3"))
884 ((TH3F*)histoList->At(currIndex))->Fill(var1, var2, var3);
885 else
886 AliWarning("Histogram not filled");
887 }
888}
889
890//________________________________________________________________________
891Int_t AliAnalysisTaskSingleMu::RecoTrackMother(AliMCParticle* mcParticle)
892{
893 //
894 /// Find track mother from kinematics
895 //
896
897 Int_t recoPdg = mcParticle->PdgCode();
898
899 // Track is not a muon
900 if ( TMath::Abs(recoPdg) != 13 ) return kRecoHadron;
901
902 Int_t imother = mcParticle->GetMother();
903
904 Bool_t isFirstMotherHF = kFALSE;
905 Int_t step = 0;
906
907 while ( imother >= 0 ) {
908 TParticle* part = ((AliMCParticle*)fMCEvent->GetTrack(imother))->Particle();
909
910 if ( part->GetUniqueID() == kPHadronic )
911 return kSecondaryMu;
912
913 Int_t absPdg = TMath::Abs(part->GetPdgCode());
914
915 step++;
916 if ( step == 1 ) // only for first mother
917 isFirstMotherHF = ( ( absPdg >= 400 && absPdg < 600 ) ||
918 ( absPdg >= 4000 && absPdg < 6000 ) );
919
920 if ( absPdg < 4 )
921 return kPrimaryMu;
922
923 if ( isFirstMotherHF) {
924 if ( absPdg == 4 )
925 return kCharmMu;
926 else if ( absPdg == 5 )
927 return kBeautyMu;
928 }
929
930 imother = part->GetFirstMother();
931 }
932
933 return kPrimaryMu;
934}
935
936//________________________________________________________________________
937Int_t AliAnalysisTaskSingleMu::GetHistoIndex(Int_t histoTypeIndex, Int_t trigIndex, Int_t srcIndex)
938{
939 //
940 /// Get histogram index in the list
941 //
942
943 if ( srcIndex < 0 ) {
944 return histoTypeIndex;
945 }
946
947 return
948 kNsummaryHistosMC +
949 kNtrackSources * kNtrigCuts * histoTypeIndex +
950 kNtrackSources * trigIndex +
951 srcIndex;
952}
953
954//________________________________________________________________________
955Float_t AliAnalysisTaskSingleMu::GetBinThetaAbsEnd(Float_t RAtAbsEnd, Bool_t isTheta)
956{
957 //
958 /// Get bin of theta at absorber end region
959 //
960 Float_t thetaDeg = ( isTheta ) ? RAtAbsEnd : TMath::ATan( RAtAbsEnd / 505. );
961 thetaDeg *= TMath::RadToDeg();
962 if ( thetaDeg < 2. )
963 return 0.;
964 else if ( thetaDeg < 3. )
965 return 1.;
966 else if ( thetaDeg < 9. )
967 return 2.;
968
969 return 3.;
970}
971
972
973//________________________________________________________________________
974void AliAnalysisTaskSingleMu::Reset(Bool_t keepGlobal)
975{
976 //
977 /// Reset variables
978 //
979 Int_t lastVarFloat = ( keepGlobal ) ? kVarIPVx : kNvarFloat;
980 for(Int_t ivar=0; ivar<lastVarFloat; ivar++){
981 fVarFloat[ivar] = 0.;
982 }
983 Int_t lastVarInt = ( keepGlobal ) ? kVarPassPhysicsSelection : kNvarInt;
984 for(Int_t ivar=0; ivar<lastVarInt; ivar++){
985 fVarInt[ivar] = 0;
986 }
987 fVarInt[kVarMatchTrig] = -1;
988
989 if ( ! keepGlobal ){
990 for(Int_t ivar=0; ivar<kNvarChar; ivar++){
991 strncpy(fVarChar[ivar]," ",255);
992 }
993 for(Int_t ivar=0; ivar<kNvarUInt; ivar++){
994 fVarUInt[ivar] = 0;
995 }
996 }
997 if ( fMCEvent ){
998 for(Int_t ivar=0; ivar<kNvarFloatMC; ivar++){
999 fVarFloatMC[ivar] = 0.;
1000 }
1001 for(Int_t ivar=0; ivar<kNvarIntMC; ivar++){
1002 fVarIntMC[ivar] = -1;
1003 }
1004 }
1005}
1006
1007//________________________________________________________________________
1008void AliAnalysisTaskSingleMu::SetTriggerClasses(TString triggerClasses)
1009{
1010 /// Set trigger classes
1011 if ( fTriggerClasses )
1012 delete fTriggerClasses;
1013
1014 fTriggerClasses = triggerClasses.Tokenize(" ");
1015 fTriggerClasses->SetOwner();
1016
1017}
1018
1019//________________________________________________________________________
1020void AliAnalysisTaskSingleMu::SetAxisLabel(TAxis* axis)
1021{
1022 //
1023 /// Set the labels of trigger class axis
1024 //
1025 for ( Int_t itrig=0; itrig<fTriggerClasses->GetEntries(); itrig++ ) {
1026 TString trigName = ((TObjString*)fTriggerClasses->At(itrig))->GetString();
1027 axis->SetBinLabel(itrig+1,trigName.Data());
1028 }
1029 axis->SetTitle("Fired class");
1030}