Update of Single Muon analysis to be included in the official ESD->AOD analysis train...
[u/mrichter/AliRoot.git] / PWG3 / muon / AliAnalysisTaskSingleMu.cxx
CommitLineData
aad6618e 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
662e37fe 16//-----------------------------------------------------------------------------
17/// \class AliAnalysisTaskSingleMu
18/// Analysis task for single muons in the spectrometer.
9728bcfd 19/// The output is a list of histograms.
20/// The macro class can run on AOD or in the train with the ESD filter.
21/// If Monte Carlo information is present, some basics checks are performed.
662e37fe 22///
23/// \author Diego Stocco
24//-----------------------------------------------------------------------------
25
aad6618e 26//----------------------------------------------------------------------------
27// Implementation of the single muon analysis class
662e37fe 28// An example of usage can be found in the macro RunSingleMuAnalysisFromAOD.C.
aad6618e 29//----------------------------------------------------------------------------
30
31#define AliAnalysisTaskSingleMu_cxx
32
33// ROOT includes
aad6618e 34#include "TROOT.h"
9728bcfd 35#include "TH1.h"
36#include "TH2.h"
37#include "TH3.h"
38#include "TAxis.h"
39#include "TList.h"
662e37fe 40#include "TCanvas.h"
9728bcfd 41#include "TMath.h"
aad6618e 42
43// STEER includes
44#include "AliLog.h"
45
46#include "AliAODEvent.h"
47#include "AliAODTrack.h"
48#include "AliAODVertex.h"
aad6618e 49
9728bcfd 50#include "AliMCEvent.h"
51#include "AliMCParticle.h"
52#include "AliStack.h"
53
54
55#include "AliAnalysisTaskSE.h"
56//#include "AliAnalysisDataSlot.h"
57//#include "AliAnalysisManager.h"
aad6618e 58#include "AliAnalysisTaskSingleMu.h"
59
662e37fe 60/// \cond CLASSIMP
61ClassImp(AliAnalysisTaskSingleMu) // Class implementation in ROOT context
62/// \endcond
aad6618e 63
9728bcfd 64
aad6618e 65//________________________________________________________________________
66AliAnalysisTaskSingleMu::AliAnalysisTaskSingleMu(const char *name) :
9728bcfd 67 AliAnalysisTaskSE(name),
68 fUseMC(0),
69 fHistoList(0),
70 fHistoListMC(0)
aad6618e 71{
72 //
73 /// Constructor.
74 //
9728bcfd 75 DefineOutput(1, TList::Class());
76 DefineOutput(2, TList::Class());
aad6618e 77}
78
9728bcfd 79
80//________________________________________________________________________
81AliAnalysisTaskSingleMu::~AliAnalysisTaskSingleMu()
aad6618e 82{
83 //
9728bcfd 84 /// Destructor
aad6618e 85 //
9728bcfd 86 if ( fHistoList ) delete fHistoList;
87 if ( fHistoListMC ) delete fHistoListMC;
88}
aad6618e 89
9728bcfd 90
91//___________________________________________________________________________
92void AliAnalysisTaskSingleMu::SetFlagsFromHandler()
93{
94 //
95 /// Use the event handler information to correctly fill the analysis flags:
96 /// - check if Monte Carlo information is present
97 //
98
99 if ( MCEvent() ) {
100 fUseMC = kTRUE;
101 AliInfo("Monte Carlo information is present");
102 }
103 else {
104 AliInfo("No Monte Carlo information in run");
aad6618e 105 }
9728bcfd 106
aad6618e 107}
108
9728bcfd 109
aad6618e 110//___________________________________________________________________________
9728bcfd 111void AliAnalysisTaskSingleMu::UserCreateOutputObjects()
aad6618e 112{
113 //
114 /// Create output histograms
115 //
9728bcfd 116 AliInfo(Form(" CreateOutputObjects of task %s\n", GetName()));
aad6618e 117
9728bcfd 118 // initialize histogram lists
119 if(!fHistoList) fHistoList = new TList();
120 if(!fHistoListMC) fHistoListMC = new TList();
121
122 Int_t nPtBins = 60;
123 Float_t ptMin = 0., ptMax = 30.;
124 TString ptName("Pt"), ptTitle("p_{t}"), ptUnits("GeV/c");
125
126 Int_t nDcaBins = 100;
127 Float_t dcaMin = 0., dcaMax = 200.;
128 TString dcaName("DCA"), dcaTitle("DCA"), dcaUnits("cm");
662e37fe 129
9728bcfd 130 Int_t nVzBins = 60;
131 Float_t vzMin = -30, vzMax = 30;
132 TString vzName("Vz"), vzTitle("Vz"), vzUnits("cm");
133
134 Int_t nEtaBins = 25;
135 Float_t etaMin = -4.5, etaMax = -2.;
136 TString etaName("Eta"), etaTitle("#eta"), etaUnits("a.u.");
137
138 Int_t nRapidityBins = 25;
139 Float_t rapidityMin = -4.5, rapidityMax = -2.;
140 TString rapidityName("Rapidity"), rapidityTitle("rapidity"), rapidityUnits("a.u.");
141
142 TString trigName[kNtrigCuts];
143 trigName[kNoMatchTrig] = "NoMatch";
144 trigName[kAllPtTrig] = "AllPt";
145 trigName[kLowPtTrig] = "LowPt";
146 trigName[kHighPtTrig] = "HighPt";
147
148 TString srcName[kNtrackSources];
149 srcName[kCharmMu] = "Charm";
150 srcName[kBeautyMu] = "Beauty";
151 srcName[kPrimaryMu] = "Decay";
152 srcName[kSecondaryMu] = "Secondary";
153 srcName[kRecoHadron] = "Hadrons";
154 srcName[kUnknownPart] = "Unknown";
155
156 TH1F* histo1D = 0x0;
157 TH2F* histo2D = 0x0;
158 TString histoName, histoTitle;
159 Int_t histoIndex = 0;
160
161 // 1D histos
162 for (Int_t itrig=0; itrig < kNtrigCuts; itrig++) {
163 histoName = Form("%sDistrib%sTrig", ptName.Data(), trigName[itrig].Data());
164 histoTitle = Form("%s distribution. Trigger: %s", ptTitle.Data(), trigName[itrig].Data());
165 histo1D = new TH1F(histoName.Data(), histoTitle.Data(), nPtBins, ptMin, ptMax);
166 histo1D->GetXaxis()->SetTitle(Form("%s (%s)", ptTitle.Data(), ptUnits.Data()));
167 histoIndex = GetHistoIndex(kHistoPt, itrig);
168 fHistoList->AddAt(histo1D, histoIndex);
169 }
170
171 for (Int_t itrig=0; itrig < kNtrigCuts; itrig++) {
172 histoName = Form("%sDistrib%sTrig", dcaName.Data(), trigName[itrig].Data());
173 histoTitle = Form("%s distribution. Trigger: %s", dcaTitle.Data(), trigName[itrig].Data());
174 histo1D = new TH1F(histoName.Data(), histoTitle.Data(), nDcaBins, dcaMin, dcaMax);
175 histo1D->GetXaxis()->SetTitle(Form("%s (%s)", dcaTitle.Data(), dcaUnits.Data()));
176 histoIndex = GetHistoIndex(kHistoDCA, itrig);
177 fHistoList->AddAt(histo1D, histoIndex);
178 }
179
180 for (Int_t itrig=0; itrig < kNtrigCuts; itrig++) {
181 histoName = Form("%sDistrib%sTrig", vzName.Data(), trigName[itrig].Data());
182 histoTitle = Form("%s distribution. Trigger: %s", vzTitle.Data(), trigName[itrig].Data());
183 histo1D = new TH1F(histoName.Data(), histoTitle.Data(), nVzBins, vzMin, vzMax);
184 histo1D->GetXaxis()->SetTitle(Form("%s (%s)", vzTitle.Data(), vzUnits.Data()));
185 histoIndex = GetHistoIndex(kHistoVz, itrig);
186 fHistoList->AddAt(histo1D, histoIndex);
187 }
188
189 for (Int_t itrig=0; itrig < kNtrigCuts; itrig++) {
190 histoName = Form("%sDistrib%sTrig", etaName.Data(), trigName[itrig].Data());
191 histoTitle = Form("%s distribution. Trigger: %s", etaTitle.Data(), trigName[itrig].Data());
192 histo1D = new TH1F(histoName.Data(), histoTitle.Data(), nEtaBins, etaMin, etaMax);
193 histo1D->GetXaxis()->SetTitle(Form("%s (%s)", etaTitle.Data(), etaUnits.Data()));
194 histoIndex = GetHistoIndex(kHistoEta, itrig);
195 fHistoList->AddAt(histo1D, histoIndex);
196 }
197
198 for (Int_t itrig=0; itrig < kNtrigCuts; itrig++) {
199 histoName = Form("%sDistrib%sTrig", rapidityName.Data(), trigName[itrig].Data());
200 histoTitle = Form("%s distribution. Trigger: %s", rapidityTitle.Data(), trigName[itrig].Data());
201 histo1D = new TH1F(histoName.Data(), histoTitle.Data(), nRapidityBins, rapidityMin, rapidityMax);
202 histo1D->GetXaxis()->SetTitle(Form("%s (%s)", rapidityTitle.Data(), rapidityUnits.Data()));
203 histoIndex = GetHistoIndex(kHistoRapidity, itrig);
204 fHistoList->AddAt(histo1D, histoIndex);
205 }
662e37fe 206
9728bcfd 207 // 2D histos
208 for (Int_t itrig=0; itrig < kNtrigCuts; itrig++) {
209 histoName = Form("%s%sDistrib%sTrig", dcaName.Data(), ptName.Data(), trigName[itrig].Data());
210 histoTitle = Form("%s vs. %s distribution. Trigger: %s", dcaTitle.Data(), ptTitle.Data(), trigName[itrig].Data());
211 histo2D = new TH2F(histoName.Data(), histoTitle.Data(),
212 nPtBins, ptMin, ptMax,
213 nDcaBins, dcaMin, dcaMax);
214 histo2D->GetXaxis()->SetTitle(Form("%s (%s)", ptTitle.Data(), ptUnits.Data()));
215 histo2D->GetYaxis()->SetTitle(Form("%s (%s)", dcaTitle.Data(), dcaUnits.Data()));
216 histoIndex = GetHistoIndex(kHistoPtDCA, itrig);
217 fHistoList->AddAt(histo2D, histoIndex);
218 }
219
220 for (Int_t itrig=0; itrig < kNtrigCuts; itrig++) {
221 histoName = Form("%s%sDistrib%sTrig", vzName.Data(), ptName.Data(), trigName[itrig].Data());
222 histoTitle = Form("%s vs. %s distribution. Trigger: %s", vzTitle.Data(), ptTitle.Data(), trigName[itrig].Data());
223 histo2D = new TH2F(histoName.Data(), histoTitle.Data(),
224 nPtBins, ptMin, ptMax,
225 nVzBins, vzMin, vzMax);
226 histo2D->GetXaxis()->SetTitle(Form("%s (%s)", ptTitle.Data(), ptUnits.Data()));
227 histo2D->GetYaxis()->SetTitle(Form("%s (%s)", vzTitle.Data(), vzUnits.Data()));
228 histoIndex = GetHistoIndex(kHistoPtVz, itrig);
229 fHistoList->AddAt(histo2D, histoIndex);
230 }
231
232 for (Int_t itrig=0; itrig < kNtrigCuts; itrig++) {
233 histoName = Form("%s%sDistrib%sTrig", rapidityName.Data(), ptName.Data(), trigName[itrig].Data());
234 histoTitle = Form("%s vs. %s distribution. Trigger: %s", rapidityTitle.Data(), ptTitle.Data(), trigName[itrig].Data());
235 histo2D = new TH2F(histoName.Data(), histoTitle.Data(),
236 nPtBins, ptMin, ptMax,
237 nRapidityBins, rapidityMin, rapidityMax);
238 histo2D->GetXaxis()->SetTitle(Form("%s (%s)", ptTitle.Data(), ptUnits.Data()));
239 histo2D->GetYaxis()->SetTitle(Form("%s (%s)", rapidityTitle.Data(), rapidityUnits.Data()));
240 histoIndex = GetHistoIndex(kHistoPtRapidity, itrig);
241 fHistoList->AddAt(histo2D, histoIndex);
662e37fe 242 }
243
9728bcfd 244 // MC histos
245 for (Int_t itrig=0; itrig < kNtrigCuts; itrig++) {
246 for (Int_t isrc = 0; isrc < kNtrackSources; isrc++) {
247 histoName = Form("%sResolution%sTrig%s", ptName.Data(), trigName[itrig].Data(), srcName[isrc].Data());
248 histoTitle = Form("%s resolution. Trigger: %s (%s)", ptTitle.Data(), trigName[itrig].Data(), srcName[isrc].Data());
249 histo1D = new TH1F(histoName.Data(), histoTitle.Data(), 100, -5., 5.);
250 histo1D->GetXaxis()->SetTitle(Form("%s^{reco} - %s^{MC} (%s)", ptTitle.Data(), ptTitle.Data(), ptUnits.Data()));
251 histoIndex = GetHistoIndex(kHistoPtResolution, itrig, isrc);
252 fHistoListMC->AddAt(histo1D, histoIndex);
253 }
254 }
255
256 for (Int_t itrig=0; itrig < kNtrigCuts; itrig++) {
257 for (Int_t isrc = 0; isrc < kNtrackSources; isrc++) {
258 histoName = Form("%s%sDistrib%sTrig%s", dcaName.Data(), ptName.Data(),
259 trigName[itrig].Data(), srcName[isrc].Data());
260 histoTitle = Form("%s vs. %s distribution. Trigger: %s (%s)", dcaTitle.Data(), ptTitle.Data(),
261 trigName[itrig].Data(), srcName[isrc].Data());
262 histo2D = new TH2F(histoName.Data(), histoTitle.Data(),
263 nPtBins, ptMin, ptMax,
264 nDcaBins, dcaMin, dcaMax);
265 histo2D->GetXaxis()->SetTitle(Form("%s (%s)", ptTitle.Data(), ptUnits.Data()));
266 histo2D->GetYaxis()->SetTitle(Form("%s (%s)", dcaTitle.Data(), dcaUnits.Data()));
267 histoIndex = GetHistoIndex(kHistoPtDCAType, itrig, isrc);
268 fHistoListMC->AddAt(histo2D, histoIndex);
269 }
270 }
271
272 for (Int_t itrig=0; itrig < kNtrigCuts; itrig++) {
273 for (Int_t isrc = 0; isrc < kNtrackSources; isrc++) {
274 histoName = Form("%s%sDistrib%sTrig%s", vzName.Data(), ptName.Data(),
275 trigName[itrig].Data(), srcName[isrc].Data());
276 histoTitle = Form("%s vs. %s distribution. Trigger: %s (%s)", vzTitle.Data(), ptTitle.Data(),
277 trigName[itrig].Data(), srcName[isrc].Data());
278 histo2D = new TH2F(histoName.Data(), histoTitle.Data(),
279 nPtBins, ptMin, ptMax,
280 nVzBins, vzMin, vzMax);
281 histo2D->GetXaxis()->SetTitle(Form("%s (%s)", ptTitle.Data(), ptUnits.Data()));
282 histo2D->GetYaxis()->SetTitle(Form("%s (%s)", vzTitle.Data(), vzUnits.Data()));
283 histoIndex = GetHistoIndex(kHistoPtVzType, itrig, isrc);
284 fHistoListMC->AddAt(histo2D, histoIndex);
285 }
aad6618e 286 }
287}
288
289//________________________________________________________________________
9728bcfd 290void AliAnalysisTaskSingleMu::UserExec(Option_t * /*option*/)
aad6618e 291{
292 //
293 /// Main loop
294 /// Called for each event
295 //
296
9728bcfd 297 AliAODEvent* aodEvent = 0x0;
298 AliMCEvent* mcEvent = 0x0;
299
300 if ( Entry()==0 )
301 SetFlagsFromHandler();
aad6618e 302
9728bcfd 303 aodEvent = dynamic_cast<AliAODEvent*> (InputEvent());
304 if ( ! aodEvent ) aodEvent = AODEvent();
305 if ( ! aodEvent ) {
306 AliError ("AOD event not found. Nothing done!");
aad6618e 307 return;
308 }
309
9728bcfd 310 Int_t nTracks = aodEvent->GetNTracks();
311
312 if ( fUseMC ) mcEvent = MCEvent();
aad6618e 313
314 // Object declaration
9728bcfd 315 AliAODTrack *aodTrack = 0x0;
316 AliMCParticle* mcTrack = 0x0;
aad6618e 317
9728bcfd 318 const Float_t kDefaultFloat = -999.;
aad6618e 319
9728bcfd 320 Float_t pt = kDefaultFloat;
321 Float_t dca = kDefaultFloat;
322 Float_t xAtDca = kDefaultFloat;
323 Float_t yAtDca = kDefaultFloat;
324 Float_t eta = kDefaultFloat;
325 Float_t rapidity = kDefaultFloat;
326 Float_t vz = kDefaultFloat;
327 Int_t trackLabel = -1;
328 Int_t matchTrig = -1;
329
330 for (Int_t itrack = 0; itrack < nTracks; itrack++) {
aad6618e 331
9728bcfd 332 // Get variables
333 aodTrack = aodEvent->GetTrack(itrack);
334 if ( aodTrack->GetMostProbablePID() != AliAODTrack::kMuon ) continue;
335 pt = aodTrack->Pt();
336 xAtDca = aodTrack->XAtDCA();
337 yAtDca = aodTrack->YAtDCA();
338 dca = TMath::Sqrt( xAtDca * xAtDca + yAtDca * yAtDca );
339 eta = aodTrack->Eta();
340 rapidity = aodTrack->Y();
341 vz = aodTrack->Zv();
342 trackLabel = aodTrack->GetLabel();
343 matchTrig = aodTrack->GetMatchTrigger();
344
345 // Fill histograms
346 FillTriggerHistos(kHistoPt, matchTrig, -1, pt);
347 FillTriggerHistos(kHistoDCA, matchTrig, -1, dca);
348 FillTriggerHistos(kHistoVz, matchTrig, -1, vz);
349 FillTriggerHistos(kHistoEta, matchTrig, -1, eta);
350 FillTriggerHistos(kHistoRapidity, matchTrig, -1, rapidity);
351
352 FillTriggerHistos(kHistoPtDCA, matchTrig, -1, pt, dca);
353 FillTriggerHistos(kHistoPtVz, matchTrig, -1, pt, vz);
354 FillTriggerHistos(kHistoPtRapidity, matchTrig, -1, pt, rapidity);
355
356 // Monte Carlo part
357 if (! fUseMC ) continue;
358
359 Int_t motherType = kUnknownPart;
360
361 AliMCParticle* matchedMCTrack = 0x0;
362
363 Int_t nMCtracks = mcEvent->GetNumberOfTracks();
364 for(Int_t imctrack=0; imctrack<nMCtracks; imctrack++){
365 mcTrack = (AliMCParticle*)mcEvent->GetTrack(imctrack);
366 if ( trackLabel == mcTrack->GetLabel() ) {
367 matchedMCTrack = mcTrack;
368 break;
369 }
370 } // loop on MC tracks
371
372 if ( matchedMCTrack )
373 motherType = RecoTrackMother(matchedMCTrack, mcEvent);
374
375 AliDebug(1, Form("Found mother %i", motherType));
376
377 if ( motherType != kUnknownPart) {
378 FillTriggerHistos(kHistoPtResolution, matchTrig, motherType, pt - mcTrack->Pt());
379 }
380 FillTriggerHistos(kHistoPtDCAType, matchTrig, motherType, pt, dca);
381 FillTriggerHistos(kHistoPtVzType, matchTrig, motherType, pt, vz);
382
383 } // loop on tracks
384
aad6618e 385 // Post final data. It will be written to a file with option "RECREATE"
9728bcfd 386 PostData(1, fHistoList);
387 if ( fUseMC ) PostData(2, fHistoListMC);
aad6618e 388}
389
390//________________________________________________________________________
391void AliAnalysisTaskSingleMu::Terminate(Option_t *) {
392 //
393 /// Draw some histogram at the end.
394 //
395 if (!gROOT->IsBatch()) {
396 TCanvas *c1 = new TCanvas("c1","Vz vs Pt",10,10,310,310);
397 c1->SetFillColor(10); c1->SetHighLightColor(10);
9728bcfd 398 c1->SetLeftMargin(0.15); c1->SetBottomMargin(0.15);
399 Int_t histoIndex = GetHistoIndex(kHistoPtVz, kAllPtTrig);
400 ((TH2F*)fHistoList->At(histoIndex))->DrawCopy("COLZ");
aad6618e 401 }
402}
403
404//________________________________________________________________________
9728bcfd 405 void AliAnalysisTaskSingleMu::FillTriggerHistos(Int_t histoIndex, Int_t matchTrig, Int_t motherType,
406 Float_t var1, Float_t var2, Float_t var3)
aad6618e 407{
408 //
9728bcfd 409 /// Fill all histograms passing the trigger cuts
aad6618e 410 //
662e37fe 411
9728bcfd 412 Int_t nTrigs = TMath::Max(1, matchTrig);
413 TArrayI trigToFill(nTrigs);
414 trigToFill[0] = matchTrig;
415 for(Int_t itrig = 1; itrig < matchTrig; itrig++){
416 trigToFill[itrig] = itrig;
417 }
662e37fe 418
9728bcfd 419 TList* histoList = (motherType < 0 ) ? fHistoList : fHistoListMC;
420
421 TString className;
422 for(Int_t itrig = 0; itrig < nTrigs; itrig++){
423 Int_t currIndex = GetHistoIndex(histoIndex, trigToFill[itrig], motherType);
424 className = histoList->At(currIndex)->ClassName();
425 AliDebug(3, Form("matchTrig %i Fill %s", matchTrig, histoList->At(currIndex)->GetName()));
426 if (className.Contains("1"))
427 ((TH1F*)histoList->At(currIndex))->Fill(var1);
428 else if (className.Contains("2"))
429 ((TH2F*)histoList->At(currIndex))->Fill(var1, var2);
430 else if (className.Contains("3"))
431 ((TH3F*)histoList->At(currIndex))->Fill(var1, var2, var3);
432 else
433 AliWarning("Histogram not filled");
434 }
aad6618e 435}
436
aad6618e 437//________________________________________________________________________
9728bcfd 438Int_t AliAnalysisTaskSingleMu::RecoTrackMother(AliMCParticle* mcTrack, AliMCEvent* mcEvent)
aad6618e 439{
440 //
9728bcfd 441 /// Find track mother from kinematics
aad6618e 442 //
443
9728bcfd 444 Int_t recoPdg = mcTrack->PdgCode();
662e37fe 445
9728bcfd 446 Bool_t isMuon = (TMath::Abs(recoPdg) == 13) ? kTRUE : kFALSE;
aad6618e 447
9728bcfd 448 // Track is not a muon
449 if ( ! isMuon ) return kRecoHadron;
aad6618e 450
9728bcfd 451 Int_t imother = mcTrack->GetMother();
aad6618e 452
9728bcfd 453 if ( imother<0 )
454 return kPrimaryMu; // Drell-Yan Muon
aad6618e 455
9728bcfd 456 Int_t igrandma = imother;
662e37fe 457
9728bcfd 458 AliMCParticle* motherPart = (AliMCParticle*)mcEvent->GetTrack(imother);
459 Int_t motherPdg = motherPart->PdgCode();
662e37fe 460
9728bcfd 461 // Track is an heavy flavor muon
462 Int_t absPdg = TMath::Abs(motherPdg);
463 if(absPdg/100==5 || absPdg/1000==5) {
464 return kBeautyMu;
465 }
466 if(absPdg/100==4 || absPdg/1000==4){
467 Int_t newMother = -1;
468 igrandma = imother;
469 AliInfo("\nFound candidate B->C->mu. History:\n");
470 mcTrack->Print();
471 printf("- %2i ", igrandma);
472 motherPart->Print();
473 Int_t absGrandMotherPdg = TMath::Abs(motherPart->PdgCode());
474 while ( absGrandMotherPdg > 10 ) {
475 igrandma = ((AliMCParticle*)mcEvent->GetTrack(igrandma))->GetMother();
476 if( igrandma < 0 ) break;
477 printf("- %2i ", igrandma);
478 mcEvent->GetTrack(igrandma)->Print();
479 absGrandMotherPdg = TMath::Abs(((AliMCParticle*)mcEvent->GetTrack(igrandma))->PdgCode());
480 }
481
482 if (absGrandMotherPdg==5) newMother = kBeautyMu; // Charm from beauty
483 else if (absGrandMotherPdg==4) newMother = kCharmMu;
484
485 if(newMother<0) {
486 AliWarning("Mother not correctly found! Set to charm!\n");
487 newMother = kCharmMu;
488 }
489
490 return newMother;
491 }
492
493 Int_t nPrimaries = mcEvent->Stack()->GetNprimary();
662e37fe 494
9728bcfd 495 // Track is a bkg. muon
496 if (imother<nPrimaries) {
497 return kPrimaryMu;
498 }
499 else {
500 return kSecondaryMu;
501 }
502}
662e37fe 503
9728bcfd 504//________________________________________________________________________
505Int_t AliAnalysisTaskSingleMu::GetHistoIndex(Int_t histoTypeIndex, Int_t trigIndex, Int_t srcIndex)
506 {
507 if ( srcIndex < 0 ) return kNtrigCuts * histoTypeIndex + trigIndex;
508
509 return
510 kNtrackSources * kNtrigCuts * histoTypeIndex +
511 kNtrackSources * trigIndex +
512 srcIndex;
aad6618e 513}