[u/mrichter/AliRoot.git] / PWG3 / muon / AliAnalysisTaskSingleMu.cxx

/**************************************************************************
 * Copyright(c) 1998-2007, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/

//-----------------------------------------------------------------------------
/// \class AliAnalysisTaskSingleMu
/// Analysis task for single muons in the spectrometer.
/// The output is a list of histograms.
/// The macro class can run on AOD or in the train with the ESD filter.
/// If Monte Carlo information is present, some basics checks are performed.
///
/// \author Diego Stocco
//-----------------------------------------------------------------------------

//----------------------------------------------------------------------------
//    Implementation of the single muon analysis class
// An example of usage can be found in the macro RunSingleMuAnalysisFromAOD.C.
//----------------------------------------------------------------------------

#define AliAnalysisTaskSingleMu_cxx

// ROOT includes
#include "TROOT.h"
#include "TH1.h"
#include "TH2.h"
#include "TH3.h"
#include "TAxis.h"
#include "TList.h"
#include "TCanvas.h"
#include "TMath.h"

// STEER includes
#include "AliLog.h"

#include "AliAODEvent.h"
#include "AliAODTrack.h"
#include "AliAODVertex.h"

#include "AliMCEvent.h"
#include "AliMCParticle.h"
#include "AliStack.h"


#include "AliAnalysisTaskSE.h"
//#include "AliAnalysisDataSlot.h"
//#include "AliAnalysisManager.h"
#include "AliAnalysisTaskSingleMu.h"

/// \cond CLASSIMP
ClassImp(AliAnalysisTaskSingleMu) // Class implementation in ROOT context
/// \endcond


//________________________________________________________________________
AliAnalysisTaskSingleMu::AliAnalysisTaskSingleMu(const char *name) :
  AliAnalysisTaskSE(name),
  fUseMC(0),
  fHistoList(0),
  fHistoListMC(0)
{
  //
  /// Constructor.
  //
  DefineOutput(1, TList::Class());
  DefineOutput(2, TList::Class());
}


//________________________________________________________________________
AliAnalysisTaskSingleMu::~AliAnalysisTaskSingleMu()
{
  //
  /// Destructor
  //
  if ( fHistoList ) delete fHistoList;
  if ( fHistoListMC ) delete fHistoListMC;
}


//___________________________________________________________________________
void AliAnalysisTaskSingleMu::SetFlagsFromHandler()
{
  //
  /// Use the event handler information to correctly fill the analysis flags:
  /// - check if Monte Carlo information is present
  //  
  
  if ( MCEvent() ) {
    fUseMC = kTRUE;
    AliInfo("Monte Carlo information is present");
  }
  else {
    AliInfo("No Monte Carlo information in run");
  }

}


//___________________________________________________________________________
void AliAnalysisTaskSingleMu::UserCreateOutputObjects() 
{
  //
  /// Create output histograms
  //
  AliInfo(Form("   CreateOutputObjects of task %s\n", GetName()));

  // initialize histogram lists
  if(!fHistoList) fHistoList = new TList();
  if(!fHistoListMC) fHistoListMC = new TList();
  
  Int_t nPtBins = 60;
  Float_t ptMin = 0., ptMax = 30.;
  TString ptName("Pt"), ptTitle("p_{t}"), ptUnits("GeV/c");
  
  Int_t nDcaBins = 100;
  Float_t dcaMin = 0., dcaMax = 200.;
  TString dcaName("DCA"), dcaTitle("DCA"), dcaUnits("cm");

  Int_t nVzBins = 60;
  Float_t vzMin = -30, vzMax = 30;
  TString vzName("Vz"), vzTitle("Vz"), vzUnits("cm");
  
  Int_t nEtaBins = 25;
  Float_t etaMin = -4.5, etaMax = -2.;
  TString etaName("Eta"), etaTitle("#eta"), etaUnits("a.u.");
  
  Int_t nRapidityBins = 25;
  Float_t rapidityMin = -4.5, rapidityMax = -2.;
  TString rapidityName("Rapidity"), rapidityTitle("rapidity"), rapidityUnits("a.u.");

  TString trigName[kNtrigCuts];
  trigName[kNoMatchTrig] = "NoMatch";
  trigName[kAllPtTrig]   = "AllPt";
  trigName[kLowPtTrig]   = "LowPt";
  trigName[kHighPtTrig]  = "HighPt";

  TString srcName[kNtrackSources];
  srcName[kCharmMu]     = "Charm";
  srcName[kBeautyMu]    = "Beauty";
  srcName[kPrimaryMu]   = "Decay";
  srcName[kSecondaryMu] = "Secondary";
  srcName[kRecoHadron]  = "Hadrons";
  srcName[kUnknownPart] = "Unknown";

  TH1F* histo1D = 0x0;
  TH2F* histo2D = 0x0;
  TString histoName, histoTitle;
  Int_t histoIndex = 0;

  // 1D histos
  for (Int_t itrig=0; itrig < kNtrigCuts; itrig++) {
    histoName = Form("%sDistrib%sTrig", ptName.Data(), trigName[itrig].Data());
    histoTitle = Form("%s distribution. Trigger: %s", ptTitle.Data(), trigName[itrig].Data());
    histo1D = new TH1F(histoName.Data(), histoTitle.Data(), nPtBins, ptMin, ptMax);
    histo1D->GetXaxis()->SetTitle(Form("%s (%s)", ptTitle.Data(), ptUnits.Data()));
    histoIndex = GetHistoIndex(kHistoPt, itrig);
    fHistoList->AddAt(histo1D, histoIndex);
  }
  
  for (Int_t itrig=0; itrig < kNtrigCuts; itrig++) {
    histoName = Form("%sDistrib%sTrig", dcaName.Data(), trigName[itrig].Data());
    histoTitle = Form("%s distribution. Trigger: %s", dcaTitle.Data(), trigName[itrig].Data());
    histo1D = new TH1F(histoName.Data(), histoTitle.Data(), nDcaBins, dcaMin, dcaMax);
    histo1D->GetXaxis()->SetTitle(Form("%s (%s)", dcaTitle.Data(), dcaUnits.Data()));
    histoIndex = GetHistoIndex(kHistoDCA, itrig);
    fHistoList->AddAt(histo1D, histoIndex);
  }
  
  for (Int_t itrig=0; itrig < kNtrigCuts; itrig++) {
    histoName = Form("%sDistrib%sTrig", vzName.Data(), trigName[itrig].Data());
    histoTitle = Form("%s distribution. Trigger: %s", vzTitle.Data(), trigName[itrig].Data());
    histo1D = new TH1F(histoName.Data(), histoTitle.Data(), nVzBins, vzMin, vzMax);
    histo1D->GetXaxis()->SetTitle(Form("%s (%s)", vzTitle.Data(), vzUnits.Data()));
    histoIndex = GetHistoIndex(kHistoVz, itrig);
    fHistoList->AddAt(histo1D, histoIndex);
  }
  
  for (Int_t itrig=0; itrig < kNtrigCuts; itrig++) {
    histoName = Form("%sDistrib%sTrig", etaName.Data(), trigName[itrig].Data());
    histoTitle = Form("%s distribution. Trigger: %s", etaTitle.Data(), trigName[itrig].Data());
    histo1D = new TH1F(histoName.Data(), histoTitle.Data(), nEtaBins, etaMin, etaMax);
    histo1D->GetXaxis()->SetTitle(Form("%s (%s)", etaTitle.Data(), etaUnits.Data()));
    histoIndex = GetHistoIndex(kHistoEta, itrig);
    fHistoList->AddAt(histo1D, histoIndex);
  }
  
  for (Int_t itrig=0; itrig < kNtrigCuts; itrig++) {
    histoName = Form("%sDistrib%sTrig", rapidityName.Data(), trigName[itrig].Data());
    histoTitle = Form("%s distribution. Trigger: %s", rapidityTitle.Data(), trigName[itrig].Data());
    histo1D = new TH1F(histoName.Data(), histoTitle.Data(), nRapidityBins, rapidityMin, rapidityMax);
    histo1D->GetXaxis()->SetTitle(Form("%s (%s)", rapidityTitle.Data(), rapidityUnits.Data()));
    histoIndex = GetHistoIndex(kHistoRapidity, itrig);
    fHistoList->AddAt(histo1D, histoIndex);
  }

  // 2D histos
  for (Int_t itrig=0; itrig < kNtrigCuts; itrig++) {
    histoName = Form("%s%sDistrib%sTrig", dcaName.Data(), ptName.Data(), trigName[itrig].Data());
    histoTitle = Form("%s vs. %s distribution. Trigger: %s", dcaTitle.Data(), ptTitle.Data(), trigName[itrig].Data());
    histo2D = new TH2F(histoName.Data(), histoTitle.Data(),
		       nPtBins, ptMin, ptMax,
		       nDcaBins, dcaMin, dcaMax);
    histo2D->GetXaxis()->SetTitle(Form("%s (%s)", ptTitle.Data(), ptUnits.Data()));
    histo2D->GetYaxis()->SetTitle(Form("%s (%s)", dcaTitle.Data(), dcaUnits.Data()));
    histoIndex = GetHistoIndex(kHistoPtDCA, itrig);
    fHistoList->AddAt(histo2D, histoIndex);
  }
  
  for (Int_t itrig=0; itrig < kNtrigCuts; itrig++) {
    histoName = Form("%s%sDistrib%sTrig", vzName.Data(), ptName.Data(), trigName[itrig].Data());
    histoTitle = Form("%s vs. %s distribution. Trigger: %s", vzTitle.Data(), ptTitle.Data(), trigName[itrig].Data());
    histo2D = new TH2F(histoName.Data(), histoTitle.Data(),
		       nPtBins, ptMin, ptMax,
		       nVzBins, vzMin, vzMax);
    histo2D->GetXaxis()->SetTitle(Form("%s (%s)", ptTitle.Data(), ptUnits.Data()));
    histo2D->GetYaxis()->SetTitle(Form("%s (%s)", vzTitle.Data(), vzUnits.Data()));
    histoIndex = GetHistoIndex(kHistoPtVz, itrig);
    fHistoList->AddAt(histo2D, histoIndex);
  }
  
  for (Int_t itrig=0; itrig < kNtrigCuts; itrig++) {
    histoName = Form("%s%sDistrib%sTrig", rapidityName.Data(), ptName.Data(), trigName[itrig].Data());
    histoTitle = Form("%s vs. %s distribution. Trigger: %s", rapidityTitle.Data(), ptTitle.Data(), trigName[itrig].Data());
    histo2D = new TH2F(histoName.Data(), histoTitle.Data(),
		     nPtBins, ptMin, ptMax,
		     nRapidityBins, rapidityMin, rapidityMax);
    histo2D->GetXaxis()->SetTitle(Form("%s (%s)", ptTitle.Data(), ptUnits.Data()));
    histo2D->GetYaxis()->SetTitle(Form("%s (%s)", rapidityTitle.Data(), rapidityUnits.Data()));
    histoIndex = GetHistoIndex(kHistoPtRapidity, itrig);
    fHistoList->AddAt(histo2D, histoIndex);
  }

  // MC histos
  for (Int_t itrig=0; itrig < kNtrigCuts; itrig++) {
    for (Int_t isrc = 0; isrc < kNtrackSources; isrc++) {
      histoName = Form("%sResolution%sTrig%s", ptName.Data(), trigName[itrig].Data(), srcName[isrc].Data());
      histoTitle = Form("%s resolution. Trigger: %s (%s)", ptTitle.Data(), trigName[itrig].Data(), srcName[isrc].Data());
      histo1D = new TH1F(histoName.Data(), histoTitle.Data(), 100, -5., 5.);
      histo1D->GetXaxis()->SetTitle(Form("%s^{reco} - %s^{MC} (%s)", ptTitle.Data(), ptTitle.Data(), ptUnits.Data()));
      histoIndex = GetHistoIndex(kHistoPtResolution, itrig, isrc);
      fHistoListMC->AddAt(histo1D, histoIndex);
    }
  }

  for (Int_t itrig=0; itrig < kNtrigCuts; itrig++) {
    for (Int_t isrc = 0; isrc < kNtrackSources; isrc++) {
      histoName = Form("%s%sDistrib%sTrig%s", dcaName.Data(), ptName.Data(),
		       trigName[itrig].Data(), srcName[isrc].Data());
      histoTitle = Form("%s vs. %s distribution. Trigger: %s (%s)", dcaTitle.Data(), ptTitle.Data(),
			trigName[itrig].Data(), srcName[isrc].Data());
      histo2D = new TH2F(histoName.Data(), histoTitle.Data(),
			 nPtBins, ptMin, ptMax,
			 nDcaBins, dcaMin, dcaMax);
      histo2D->GetXaxis()->SetTitle(Form("%s (%s)", ptTitle.Data(), ptUnits.Data()));
      histo2D->GetYaxis()->SetTitle(Form("%s (%s)", dcaTitle.Data(), dcaUnits.Data()));
      histoIndex = GetHistoIndex(kHistoPtDCAType, itrig, isrc);
      fHistoListMC->AddAt(histo2D, histoIndex);
    }
  }
  
  for (Int_t itrig=0; itrig < kNtrigCuts; itrig++) {
    for (Int_t isrc = 0; isrc < kNtrackSources; isrc++) {
      histoName = Form("%s%sDistrib%sTrig%s", vzName.Data(), ptName.Data(), 
		       trigName[itrig].Data(), srcName[isrc].Data());
      histoTitle = Form("%s vs. %s distribution. Trigger: %s (%s)", vzTitle.Data(), ptTitle.Data(),
			trigName[itrig].Data(), srcName[isrc].Data());
      histo2D = new TH2F(histoName.Data(), histoTitle.Data(),
			 nPtBins, ptMin, ptMax,
			 nVzBins, vzMin, vzMax);
      histo2D->GetXaxis()->SetTitle(Form("%s (%s)", ptTitle.Data(), ptUnits.Data()));
      histo2D->GetYaxis()->SetTitle(Form("%s (%s)", vzTitle.Data(), vzUnits.Data()));
      histoIndex = GetHistoIndex(kHistoPtVzType, itrig, isrc);
      fHistoListMC->AddAt(histo2D, histoIndex);
    }
  }
}

//________________________________________________________________________
void AliAnalysisTaskSingleMu::UserExec(Option_t * /*option*/) 
{
  //
  /// Main loop
  /// Called for each event
  //

  AliAODEvent* aodEvent = 0x0;
  AliMCEvent*  mcEvent  = 0x0;

  if ( Entry()==0 ) 
    SetFlagsFromHandler();

  aodEvent = dynamic_cast<AliAODEvent*> (InputEvent());
  if ( ! aodEvent ) aodEvent = AODEvent();
  if ( ! aodEvent ) {
    AliError ("AOD event not found. Nothing done!");
    return;
  }

  Int_t nTracks = aodEvent->GetNTracks();

  if ( fUseMC ) mcEvent = MCEvent();

  // Object declaration
  AliAODTrack *aodTrack = 0x0;
  AliMCParticle* mcTrack = 0x0;

  const Float_t kDefaultFloat = -999.;

  Float_t pt       = kDefaultFloat;
  Float_t dca      = kDefaultFloat;
  Float_t xAtDca   = kDefaultFloat;
  Float_t yAtDca   = kDefaultFloat;
  Float_t eta      = kDefaultFloat;
  Float_t rapidity = kDefaultFloat;
  Float_t vz       = kDefaultFloat;
  Int_t trackLabel = -1;
  Int_t matchTrig = -1;

  for (Int_t itrack = 0; itrack < nTracks; itrack++) {

    // Get variables
    aodTrack = aodEvent->GetTrack(itrack);
    if ( aodTrack->GetMostProbablePID() != AliAODTrack::kMuon ) continue;
    pt = aodTrack->Pt();
    xAtDca = aodTrack->XAtDCA();
    yAtDca = aodTrack->YAtDCA();
    dca = TMath::Sqrt( xAtDca * xAtDca + yAtDca * yAtDca );
    eta = aodTrack->Eta();
    rapidity = aodTrack->Y();
    vz = aodTrack->Zv();
    trackLabel = aodTrack->GetLabel();
    matchTrig = aodTrack->GetMatchTrigger();

    // Fill histograms
    FillTriggerHistos(kHistoPt,       matchTrig, -1, pt);
    FillTriggerHistos(kHistoDCA,      matchTrig, -1, dca);
    FillTriggerHistos(kHistoVz,       matchTrig, -1, vz);
    FillTriggerHistos(kHistoEta,      matchTrig, -1, eta);
    FillTriggerHistos(kHistoRapidity, matchTrig, -1, rapidity);

    FillTriggerHistos(kHistoPtDCA,      matchTrig, -1, pt, dca);
    FillTriggerHistos(kHistoPtVz,       matchTrig, -1, pt, vz);
    FillTriggerHistos(kHistoPtRapidity, matchTrig, -1, pt, rapidity);

    // Monte Carlo part
    if (! fUseMC ) continue;

    Int_t motherType = kUnknownPart;

    AliMCParticle* matchedMCTrack = 0x0;

    Int_t nMCtracks = mcEvent->GetNumberOfTracks();
    for(Int_t imctrack=0; imctrack<nMCtracks; imctrack++){
      mcTrack = (AliMCParticle*)mcEvent->GetTrack(imctrack);
      if ( trackLabel == mcTrack->GetLabel() ) {
	matchedMCTrack = mcTrack;
	break;
      }
    } // loop on MC tracks

    if ( matchedMCTrack )
      motherType = RecoTrackMother(matchedMCTrack, mcEvent);

    AliDebug(1, Form("Found mother %i", motherType));

    if ( motherType != kUnknownPart) {
      FillTriggerHistos(kHistoPtResolution, matchTrig, motherType, pt - mcTrack->Pt());
    }
    FillTriggerHistos(kHistoPtDCAType, matchTrig, motherType, pt, dca);
    FillTriggerHistos(kHistoPtVzType,  matchTrig, motherType, pt, vz);

  } // loop on tracks
  
  // Post final data. It will be written to a file with option "RECREATE"
  PostData(1, fHistoList);
  if ( fUseMC ) PostData(2, fHistoListMC);
}

//________________________________________________________________________
void AliAnalysisTaskSingleMu::Terminate(Option_t *) {
  //
  /// Draw some histogram at the end.
  //
  if (!gROOT->IsBatch()) {
    fHistoList = dynamic_cast<TList*> (GetOutputData(1));
    TCanvas *c1_SingleMu = new TCanvas("c1_SingleMu","Vz vs Pt",10,10,310,310);
    c1_SingleMu->SetFillColor(10); c1_SingleMu->SetHighLightColor(10);
    c1_SingleMu->SetLeftMargin(0.15); c1_SingleMu->SetBottomMargin(0.15);
    Int_t histoIndex = GetHistoIndex(kHistoPtVz, kAllPtTrig);
    ((TH2F*)fHistoList->At(histoIndex))->DrawCopy("COLZ");
  }
}

//________________________________________________________________________
 void AliAnalysisTaskSingleMu::FillTriggerHistos(Int_t histoIndex, Int_t matchTrig, Int_t motherType,
						Float_t var1, Float_t var2, Float_t var3)
{
  //
  /// Fill all histograms passing the trigger cuts
  //

  Int_t nTrigs = TMath::Max(1, matchTrig);
  TArrayI trigToFill(nTrigs);
  trigToFill[0] = matchTrig;
  for(Int_t itrig = 1; itrig < matchTrig; itrig++){
    trigToFill[itrig] = itrig;
  }

  TList* histoList = (motherType < 0 ) ? fHistoList : fHistoListMC;

  TString className;
  for(Int_t itrig = 0; itrig < nTrigs; itrig++){
    Int_t currIndex = GetHistoIndex(histoIndex, trigToFill[itrig], motherType);
    className = histoList->At(currIndex)->ClassName();
    AliDebug(3, Form("matchTrig %i  Fill %s", matchTrig, histoList->At(currIndex)->GetName()));
    if (className.Contains("1"))
      ((TH1F*)histoList->At(currIndex))->Fill(var1);
    else if (className.Contains("2"))
      ((TH2F*)histoList->At(currIndex))->Fill(var1, var2);
    else if (className.Contains("3"))
      ((TH3F*)histoList->At(currIndex))->Fill(var1, var2, var3);
    else
      AliWarning("Histogram not filled");
  }
}

//________________________________________________________________________
Int_t AliAnalysisTaskSingleMu::RecoTrackMother(AliMCParticle* mcTrack, AliMCEvent* mcEvent)
{
  //
  /// Find track mother from kinematics
  //

  Int_t recoPdg = mcTrack->PdgCode();

  Bool_t isMuon = (TMath::Abs(recoPdg) == 13) ? kTRUE : kFALSE;

  // Track is not a muon
  if ( ! isMuon ) return kRecoHadron;

  Int_t imother = mcTrack->GetMother();

  if ( imother<0 ) 
    return kPrimaryMu; // Drell-Yan Muon

  Int_t igrandma = imother;

  AliMCParticle* motherPart = (AliMCParticle*)mcEvent->GetTrack(imother);
  Int_t motherPdg = motherPart->PdgCode();

  // Track is an heavy flavor muon
  Int_t absPdg = TMath::Abs(motherPdg);
  if(absPdg/100==5 || absPdg/1000==5) {
    return kBeautyMu;
  }
  if(absPdg/100==4 || absPdg/1000==4){
    Int_t newMother = -1;
    igrandma = imother;
    AliInfo("\nFound candidate B->C->mu. History:\n");
    mcTrack->Print();
    printf("- %2i   ", igrandma);
    motherPart->Print();
    Int_t absGrandMotherPdg = TMath::Abs(motherPart->PdgCode());
    while ( absGrandMotherPdg > 10 ) {
      igrandma = ((AliMCParticle*)mcEvent->GetTrack(igrandma))->GetMother();
      if( igrandma < 0 ) break;
      printf("- %2i   ", igrandma);
      mcEvent->GetTrack(igrandma)->Print();
      absGrandMotherPdg = TMath::Abs(((AliMCParticle*)mcEvent->GetTrack(igrandma))->PdgCode());
    }

    if (absGrandMotherPdg==5) newMother = kBeautyMu; // Charm from beauty
    else if (absGrandMotherPdg==4) newMother = kCharmMu;

    if(newMother<0) {
      AliWarning("Mother not correctly found! Set to charm!\n");
      newMother = kCharmMu;
    }

    return newMother;
  }

  Int_t nPrimaries = mcEvent->Stack()->GetNprimary();

  // Track is a bkg. muon
  if (imother<nPrimaries) {
    return kPrimaryMu;
  }
  else {
    return kSecondaryMu;
  }
}

//________________________________________________________________________
Int_t AliAnalysisTaskSingleMu::GetHistoIndex(Int_t histoTypeIndex, Int_t trigIndex, Int_t srcIndex)
 {
   if ( srcIndex < 0 ) return kNtrigCuts * histoTypeIndex + trigIndex;

   return 
     kNtrackSources * kNtrigCuts * histoTypeIndex  + 
     kNtrackSources * trigIndex  + 
     srcIndex;
}
Commit	Line	Data
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	//----------------------------------------------------------------------------
aad6618e	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"
aad6618e	59
662e37fe	60	/// \cond CLASSIMP
	61	ClassImp(AliAnalysisTaskSingleMu) // Class implementation in ROOT context
	62	/// \endcond
aad6618e	63
9728bcfd	64
aad6618e	65	//________________________________________________________________________
aad6618e	66	AliAnalysisTaskSingleMu::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());
9728bcfd	76	DefineOutput(2, TList::Class());
aad6618e	77	}
aad6618e	78
9728bcfd	79
	80	//________________________________________________________________________
	81	AliAnalysisTaskSingleMu::~AliAnalysisTaskSingleMu()
aad6618e	82	{
aad6618e	83	//
9728bcfd	84	/// Destructor
aad6618e	85	//
9728bcfd	86	if ( fHistoList ) delete fHistoList;
	87	if ( fHistoListMC ) delete fHistoListMC;
	88	}
aad6618e	89
9728bcfd	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");
aad6618e	105	}
9728bcfd	106
aad6618e	107	}
aad6618e	108
9728bcfd	109
aad6618e	110	//___________________________________________________________________________
9728bcfd	111	void 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	}
662e37fe	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	290	void 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
aad6618e	314	// Object declaration
9728bcfd	315	AliAODTrack *aodTrack = 0x0;
9728bcfd	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);
9728bcfd	387	if ( fUseMC ) PostData(2, fHistoListMC);
aad6618e	388	}
	389
	390	//________________________________________________________________________
	391	void AliAnalysisTaskSingleMu::Terminate(Option_t *) {
	392	//
	393	/// Draw some histogram at the end.
	394	//
	395	if (!gROOT->IsBatch()) {
93fd3322	396	fHistoList = dynamic_cast<TList*> (GetOutputData(1));
	397	TCanvas *c1_SingleMu = new TCanvas("c1_SingleMu","Vz vs Pt",10,10,310,310);
	398	c1_SingleMu->SetFillColor(10); c1_SingleMu->SetHighLightColor(10);
	399	c1_SingleMu->SetLeftMargin(0.15); c1_SingleMu->SetBottomMargin(0.15);
9728bcfd	400	Int_t histoIndex = GetHistoIndex(kHistoPtVz, kAllPtTrig);
9728bcfd	401	((TH2F*)fHistoList->At(histoIndex))->DrawCopy("COLZ");
aad6618e	402	}
	403	}
	404
	405	//________________________________________________________________________
9728bcfd	406	void AliAnalysisTaskSingleMu::FillTriggerHistos(Int_t histoIndex, Int_t matchTrig, Int_t motherType,
9728bcfd	407	Float_t var1, Float_t var2, Float_t var3)
aad6618e	408	{
aad6618e	409	//
9728bcfd	410	/// Fill all histograms passing the trigger cuts
aad6618e	411	//
662e37fe	412
9728bcfd	413	Int_t nTrigs = TMath::Max(1, matchTrig);
	414	TArrayI trigToFill(nTrigs);
	415	trigToFill[0] = matchTrig;
	416	for(Int_t itrig = 1; itrig < matchTrig; itrig++){
	417	trigToFill[itrig] = itrig;
	418	}
662e37fe	419
9728bcfd	420	TList* histoList = (motherType < 0 ) ? fHistoList : fHistoListMC;
	421
	422	TString className;
	423	for(Int_t itrig = 0; itrig < nTrigs; itrig++){
	424	Int_t currIndex = GetHistoIndex(histoIndex, trigToFill[itrig], motherType);
	425	className = histoList->At(currIndex)->ClassName();
	426	AliDebug(3, Form("matchTrig %i Fill %s", matchTrig, histoList->At(currIndex)->GetName()));
	427	if (className.Contains("1"))
	428	((TH1F*)histoList->At(currIndex))->Fill(var1);
	429	else if (className.Contains("2"))
	430	((TH2F*)histoList->At(currIndex))->Fill(var1, var2);
	431	else if (className.Contains("3"))
	432	((TH3F*)histoList->At(currIndex))->Fill(var1, var2, var3);
	433	else
	434	AliWarning("Histogram not filled");
	435	}
aad6618e	436	}
aad6618e	437
aad6618e	438	//________________________________________________________________________
9728bcfd	439	Int_t AliAnalysisTaskSingleMu::RecoTrackMother(AliMCParticle* mcTrack, AliMCEvent* mcEvent)
aad6618e	440	{
aad6618e	441	//
9728bcfd	442	/// Find track mother from kinematics
aad6618e	443	//
aad6618e	444
9728bcfd	445	Int_t recoPdg = mcTrack->PdgCode();
662e37fe	446
9728bcfd	447	Bool_t isMuon = (TMath::Abs(recoPdg) == 13) ? kTRUE : kFALSE;
aad6618e	448
9728bcfd	449	// Track is not a muon
9728bcfd	450	if ( ! isMuon ) return kRecoHadron;
aad6618e	451
9728bcfd	452	Int_t imother = mcTrack->GetMother();
aad6618e	453
9728bcfd	454	if ( imother<0 )
9728bcfd	455	return kPrimaryMu; // Drell-Yan Muon
aad6618e	456
9728bcfd	457	Int_t igrandma = imother;
662e37fe	458
9728bcfd	459	AliMCParticle* motherPart = (AliMCParticle*)mcEvent->GetTrack(imother);
9728bcfd	460	Int_t motherPdg = motherPart->PdgCode();
662e37fe	461
9728bcfd	462	// Track is an heavy flavor muon
	463	Int_t absPdg = TMath::Abs(motherPdg);
	464	if(absPdg/100==5 \|\| absPdg/1000==5) {
	465	return kBeautyMu;
	466	}
	467	if(absPdg/100==4 \|\| absPdg/1000==4){
	468	Int_t newMother = -1;
	469	igrandma = imother;
	470	AliInfo("\nFound candidate B->C->mu. History:\n");
	471	mcTrack->Print();
	472	printf("- %2i ", igrandma);
	473	motherPart->Print();
	474	Int_t absGrandMotherPdg = TMath::Abs(motherPart->PdgCode());
	475	while ( absGrandMotherPdg > 10 ) {
	476	igrandma = ((AliMCParticle*)mcEvent->GetTrack(igrandma))->GetMother();
	477	if( igrandma < 0 ) break;
	478	printf("- %2i ", igrandma);
	479	mcEvent->GetTrack(igrandma)->Print();
	480	absGrandMotherPdg = TMath::Abs(((AliMCParticle*)mcEvent->GetTrack(igrandma))->PdgCode());
	481	}
	482
	483	if (absGrandMotherPdg==5) newMother = kBeautyMu; // Charm from beauty
	484	else if (absGrandMotherPdg==4) newMother = kCharmMu;
	485
	486	if(newMother<0) {
	487	AliWarning("Mother not correctly found! Set to charm!\n");
	488	newMother = kCharmMu;
	489	}
	490
	491	return newMother;
	492	}
	493
	494	Int_t nPrimaries = mcEvent->Stack()->GetNprimary();
662e37fe	495
9728bcfd	496	// Track is a bkg. muon
	497	if (imother<nPrimaries) {
	498	return kPrimaryMu;
	499	}
	500	else {
	501	return kSecondaryMu;
	502	}
	503	}
662e37fe	504
9728bcfd	505	//________________________________________________________________________
	506	Int_t AliAnalysisTaskSingleMu::GetHistoIndex(Int_t histoTypeIndex, Int_t trigIndex, Int_t srcIndex)
	507	{
	508	if ( srcIndex < 0 ) return kNtrigCuts * histoTypeIndex + trigIndex;
	509
	510	return
	511	kNtrackSources * kNtrigCuts * histoTypeIndex +
	512	kNtrackSources * trigIndex +
	513	srcIndex;
aad6618e	514	}