PWG3/muon/AliAnalysisTaskSingleMu.cxx

   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 AOD or in the train with the ESD filter.
  21 /// If Monte Carlo information is present, some basics checks are performed.
  22 ///
  23 /// \author Diego Stocco
  24 //-----------------------------------------------------------------------------
  25
  26 //----------------------------------------------------------------------------
  27 //    Implementation of the single muon analysis class
  28 // An example of usage can be found in the macro RunSingleMuAnalysisFromAOD.C.
  29 //----------------------------------------------------------------------------
  30
  31 #define AliAnalysisTaskSingleMu_cxx
  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 "TList.h"
  40 #include "TCanvas.h"
  41 #include "TMath.h"
  42
  43 // STEER includes
  44 #include "AliLog.h"
  45
  46 #include "AliAODEvent.h"
  47 #include "AliAODTrack.h"
  48 #include "AliAODVertex.h"
  49
  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"
  58 #include "AliAnalysisTaskSingleMu.h"
  59
  60 /// \cond CLASSIMP
  61 ClassImp(AliAnalysisTaskSingleMu) // Class implementation in ROOT context
  62 /// \endcond
  63
  64
  65 //________________________________________________________________________
  66 AliAnalysisTaskSingleMu::AliAnalysisTaskSingleMu(const char *name) :
  67   AliAnalysisTaskSE(name),
  68   fUseMC(0),
  69   fHistoList(0),
  70   fHistoListMC(0)
  71 {
  72   //
  73   /// Constructor.
  74   //
  75   DefineOutput(1, TList::Class());
  76   DefineOutput(2, TList::Class());
  77 }
  78
  79
  80 //________________________________________________________________________
  81 AliAnalysisTaskSingleMu::~AliAnalysisTaskSingleMu()
  82 {
  83   //
  84   /// Destructor
  85   //
  86   if ( fHistoList ) delete fHistoList;
  87   if ( fHistoListMC ) delete fHistoListMC;
  88 }
  89
  90
  91 //___________________________________________________________________________
  92 void 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");
 105   }
 106
 107 }
 108
 109
 110 //___________________________________________________________________________
 111 void AliAnalysisTaskSingleMu::UserCreateOutputObjects()
 112 {
 113   //
 114   /// Create output histograms
 115   //
 116   AliInfo(Form("   CreateOutputObjects of task %s\n", GetName()));
 117
 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");
 129
 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   }
 206
 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);
 242   }
 243
 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     }
 286   }
 287 }
 288
 289 //________________________________________________________________________
 290 void AliAnalysisTaskSingleMu::UserExec(Option_t * /*option*/)
 291 {
 292   //
 293   /// Main loop
 294   /// Called for each event
 295   //
 296
 297   AliAODEvent* aodEvent = 0x0;
 298   AliMCEvent*  mcEvent  = 0x0;
 299
 300   if ( Entry()==0 )
 301     SetFlagsFromHandler();
 302
 303   aodEvent = dynamic_cast<AliAODEvent*> (InputEvent());
 304   if ( ! aodEvent ) aodEvent = AODEvent();
 305   if ( ! aodEvent ) {
 306     AliError ("AOD event not found. Nothing done!");
 307     return;
 308   }
 309
 310   Int_t nTracks = aodEvent->GetNTracks();
 311
 312   if ( fUseMC ) mcEvent = MCEvent();
 313
 314   // Object declaration
 315   AliAODTrack *aodTrack = 0x0;
 316   AliMCParticle* mcTrack = 0x0;
 317
 318   const Float_t kDefaultFloat = -999.;
 319
 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++) {
 331
 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
 385   // Post final data. It will be written to a file with option "RECREATE"
 386   PostData(1, fHistoList);
 387   if ( fUseMC ) PostData(2, fHistoListMC);
 388 }
 389
 390 //________________________________________________________________________
 391 void 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);
 398     c1->SetLeftMargin(0.15); c1->SetBottomMargin(0.15);
 399     Int_t histoIndex = GetHistoIndex(kHistoPtVz, kAllPtTrig);
 400     ((TH2F*)fHistoList->At(histoIndex))->DrawCopy("COLZ");
 401   }
 402 }
 403
 404 //________________________________________________________________________
 405  void AliAnalysisTaskSingleMu::FillTriggerHistos(Int_t histoIndex, Int_t matchTrig, Int_t motherType,
 406                                                 Float_t var1, Float_t var2, Float_t var3)
 407 {
 408   //
 409   /// Fill all histograms passing the trigger cuts
 410   //
 411
 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   }
 418
 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   }
 435 }
 436
 437 //________________________________________________________________________
 438 Int_t AliAnalysisTaskSingleMu::RecoTrackMother(AliMCParticle* mcTrack, AliMCEvent* mcEvent)
 439 {
 440   //
 441   /// Find track mother from kinematics
 442   //
 443
 444   Int_t recoPdg = mcTrack->PdgCode();
 445
 446   Bool_t isMuon = (TMath::Abs(recoPdg) == 13) ? kTRUE : kFALSE;
 447
 448   // Track is not a muon
 449   if ( ! isMuon ) return kRecoHadron;
 450
 451   Int_t imother = mcTrack->GetMother();
 452
 453   if ( imother<0 )
 454     return kPrimaryMu; // Drell-Yan Muon
 455
 456   Int_t igrandma = imother;
 457
 458   AliMCParticle* motherPart = (AliMCParticle*)mcEvent->GetTrack(imother);
 459   Int_t motherPdg = motherPart->PdgCode();
 460
 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();
 494
 495   // Track is a bkg. muon
 496   if (imother<nPrimaries) {
 497     return kPrimaryMu;
 498   }
 499   else {
 500     return kSecondaryMu;
 501   }
 502 }
 503
 504 //________________________________________________________________________
 505 Int_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;
 513 }