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

/**************************************************************************
 * Copyright(c) 1998-2006, 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 used to extract and store info of MC particle
//
// Author: X-M. Zhang, zhang@clermont.in2p3.fr
//                     zhangxm@iopp.ccnu.edu.cn
/////////////////////////////////////////////////////////////

#include <TParticle.h>

#include "AliMCEvent.h"
#include "AliAODMCParticle.h"
#include "AliESDMuonTrack.h"
#include "AliAODTrack.h"
#include "AliMuonInfoStoreRD.h"
#include "AliMuonInfoStoreMC.h"
#include "AliMCEvent.h"

ClassImp(AliMuonInfoStoreMC)

const TString AliMuonInfoStoreMC::fgkStdBranchName("MuonMC");
const Int_t   AliMuonInfoStoreMC::fgkSourcesN = 6;

//-----------------------------------------------------------------------------
AliMuonInfoStoreMC::AliMuonInfoStoreMC() :
AliMuonInfoStoreRD(),
fIsFull(kFALSE),
fLorentzP(),
fTrackIndex(-1),
fTrackPDGCode(0),
fSource(-1),
fNParents(0),
fOscillation(kFALSE),
fWeight(0.)
{
  //
  // default constructor
  //
  for (Int_t i=5; i--;) { fParentIndex[i] = -1; fParentPDGCode[i] = 0; }
  for (Int_t i=4; i--;) { fQuarkIndex[i]  = -1; fQuarkPDGCode[i]  = 0; }
}

//-----------------------------------------------------------------------------
AliMuonInfoStoreMC::AliMuonInfoStoreMC(AliAODTrack *trkAOD, AliMCEvent *mcEvent, Bool_t full) :
AliMuonInfoStoreRD(trkAOD),
fIsFull(full),
fLorentzP(),
fTrackIndex(-1),
fTrackPDGCode(0),
fSource(-1),
fNParents(0),
fOscillation(kFALSE),
fWeight(0.)
{
  //
  // default constructor
  //
  for (Int_t i=5; i--;) { fParentIndex[i] = -1; fParentPDGCode[i] = 0; }
  for (Int_t i=4; i--;) { fQuarkIndex[i]  = -1; fQuarkPDGCode[i]  = 0; }

  this->SetMCInfoAOD(mcEvent, trkAOD->GetLabel());
}

//-----------------------------------------------------------------------------
AliMuonInfoStoreMC::AliMuonInfoStoreMC(AliESDMuonTrack *trkESD, AliMCEvent *mcEvent, Bool_t full) :
AliMuonInfoStoreRD(trkESD),
fIsFull(full),
fLorentzP(),
fTrackIndex(-1),
fTrackPDGCode(0),
fSource(-1),
fNParents(0),
fOscillation(kFALSE),
fWeight(0.)
{
  //
  // default constructor
  //
  for (Int_t i=5; i--;) { fParentIndex[i] = -1; fParentPDGCode[i] = 0; }
  for (Int_t i=4; i--;) { fQuarkIndex[i]  = -1; fQuarkPDGCode[i]  = 0; }

  this->SetMCInfoESD(mcEvent, trkESD->GetLabel());
}

//-----------------------------------------------------------------------------
AliMuonInfoStoreMC::AliMuonInfoStoreMC(const AliMuonInfoStoreMC &src) :
AliMuonInfoStoreRD(src),
fIsFull(src.fIsFull),
fLorentzP(src.fLorentzP),
fTrackIndex(src.fTrackIndex),
fTrackPDGCode(src.fTrackPDGCode),
fSource(src.fSource),
fNParents(src.fNParents),
fOscillation(src.fOscillation),
fWeight(src.fWeight)
{
  //
  // copy constructor
  //
  for (Int_t i=5; i--;) {
    fParentIndex[i]   = src.fParentIndex[i];
    fParentPDGCode[i] = src.fParentPDGCode[i];
  }
  for (Int_t i=4; i--;) {
    fQuarkIndex[i]    = src.fQuarkIndex[i];
    fQuarkPDGCode[i]  = src.fQuarkPDGCode[i];
  }
}

//-----------------------------------------------------------------------------
AliMuonInfoStoreMC& AliMuonInfoStoreMC::operator=(const AliMuonInfoStoreMC &src)
{
  //
  // assignment constructor
  //
  if(&src==this) return *this;
  AliMuonInfoStoreRD::operator=(src);

  fIsFull       = src.fIsFull;
  fLorentzP     = src.fLorentzP;
  fTrackIndex   = src.fTrackIndex;
  fTrackPDGCode = src.fTrackPDGCode;
  fSource       = src.fSource;
  fNParents     = src.fNParents;
  fOscillation  = src.fOscillation;
  fWeight       = src.fWeight;
  for (Int_t i=5; i--;) {
    fParentIndex[i]   = src.fParentIndex[i];
    fParentPDGCode[i] = src.fParentPDGCode[i];
  }
  for (Int_t i=4; i--;) {
    fQuarkIndex[i]    = src.fQuarkIndex[i];
    fQuarkPDGCode[i]  = src.fQuarkPDGCode[i];
  }

  return *this;
}

//-----------------------------------------------------------------------------
AliMuonInfoStoreMC::~AliMuonInfoStoreMC()
{
  //
  // destructor
  //
}

//-----------------------------------------------------------------------------
void AliMuonInfoStoreMC::SetMCInfoAOD(AliMCEvent *mcEvent, Int_t label)
{
  // fill track MC info with AOD base
  fTrackIndex = label;
  if (fTrackIndex<0) { fSource=5; return; }

  AliAODMCParticle *pMC = (AliAODMCParticle*)mcEvent->GetTrack(fTrackIndex);
  fLorentzP.SetPxPyPzE(pMC->Px(), pMC->Py(), pMC->Pz(), pMC->E());

  fTrackPDGCode = pMC->PdgCode();
  if (TMath::Abs(fTrackPDGCode)!=13) { fSource=4; return; } 

  Int_t lineM = pMC->GetMother();
  if (lineM<0) { fSource=2; return; }

  Bool_t isPrimary = ((AliAODMCParticle*)mcEvent->GetTrack(lineM))->IsPrimary();
  if (!isPrimary) { fSource=3; return; }

  Int_t countP=-1, pdg=0;
  Int_t parents[10], parLine[10];
  AliAODMCParticle *mother = 0;
  while(lineM>=0){
    mother = (AliAODMCParticle*)mcEvent->GetTrack(lineM);
    pdg = mother->GetPdgCode();
    if(pdg==92 || pdg==21 || TMath::Abs(pdg)<10 || IsDiquark(pdg)) break;
    parents[++countP] = pdg;
    parLine[countP] = lineM;
    lineM = mother->GetMother();
  }
  for(Int_t i=0; i<=countP; i++){
    fParentIndex[i] = parLine[countP-i];
    fParentPDGCode[i] = parents[countP-i];
  }
  fNParents = countP + 1;

  if (fIsFull && lineM>=0) this->FillHistoryQuarksAOD(mcEvent, lineM);

  fSource = this->SelectHFMuon();
  return;
}

//-----------------------------------------------------------------------------
void AliMuonInfoStoreMC::SetMCInfoESD(AliMCEvent *mcEvent, Int_t label)
{
  // fill track MC info with ESD base
  fTrackIndex = label;
  if (fTrackIndex<0) { fSource=5; return; }

  TParticle *pMC = ((AliMCParticle*)mcEvent->GetTrack(fTrackIndex))->Particle();
  fLorentzP.SetPxPyPzE(pMC->Px(), pMC->Py(), pMC->Pz(), pMC->Energy());

  fTrackPDGCode = pMC->GetPdgCode();
  if (TMath::Abs(fTrackPDGCode)!=13) { fSource=4; return; }

  Int_t lineM = pMC->GetFirstMother();
  if (lineM<0) { fSource=2; return; }

  if (lineM>=mcEvent->Stack()->GetNprimary()) { fSource=3; return; }

  Int_t countP=-1, pdg=0;
  Int_t parents[10], parLine[10];
  TParticle *mother = 0;
  while(lineM>=0){
    mother = ((AliMCParticle*)mcEvent->GetTrack(lineM))->Particle();
    pdg = mother->GetPdgCode();
    if(pdg==92 || pdg==21 || TMath::Abs(pdg)<10 || IsDiquark(pdg)) break;
    parents[++countP] = pdg;
    parLine[countP] = lineM;
    lineM = mother->GetFirstMother();
  }
  for(Int_t i=0; i<=countP; i++){
    fParentIndex[i] = parLine[countP-i];
    fParentPDGCode[i] = parents[countP-i];
  }
  fNParents = countP + 1;

  if (fIsFull && lineM>=0) this->FillHistoryQuarksESD(mcEvent, lineM);

  fSource = this->SelectHFMuon();
  return;
}

//-----------------------------------------------------------------------------
void AliMuonInfoStoreMC::FillHistoryQuarksAOD(AliMCEvent* const mcEvent, Int_t lineM)
{
  // method in $ALICE_ROOT/MUON/AliMUONTrackLight.cxx 

  if (lineM<0) return;
  Int_t countP=-1, pdg=0;
  AliAODMCParticle *mother = 0;
  while(lineM>=0){
    mother = (AliAODMCParticle*)mcEvent->GetTrack(lineM);
    pdg = mother->GetPdgCode();
    fQuarkIndex[++countP] = lineM;
    fQuarkPDGCode[countP] = pdg;
    lineM = mother->GetMother();
  }

  // for PYTHIA checking
  countP = 1;
  for(Int_t par=0; par<4; par++) {
    if(TMath::Abs(this->QuarkPDGCode(par))<6) { countP=par; break; }
  }
  if(this->QuarkIndex(countP)>-1 && (this->ParentFlavour(0)==4 || this->ParentFlavour(0)==5)) {
    if(this->ParentFlavour(0)!=TMath::Abs(this->QuarkPDGCode(countP))) {
      AliWarning(Form("quark flavour of parent and that of quark do not correspond: %d %d --> correcting\n",
                 this->ParentFlavour(0), TMath::Abs(this->QuarkPDGCode(countP))));

      pdg = this->QuarkPDGCode(countP);
      Int_t line = this->QuarkIndex(countP);
      this->ResetQuarkInfo();
      while(TMath::Abs(pdg)!=this->ParentFlavour(0)) {
        pdg = ((AliAODMCParticle*)mcEvent->GetTrack(++line))->GetPdgCode();
      }
      while(line>=0){
        mother = (AliAODMCParticle*)mcEvent->GetTrack(line);
        pdg = mother->GetPdgCode();
        fQuarkIndex[countP] = line;
        fQuarkPDGCode[countP++] = pdg;
        line = mother->GetMother();
      }
    }
  }
  return;
}

//-----------------------------------------------------------------------------
void AliMuonInfoStoreMC::FillHistoryQuarksESD(AliMCEvent* const mcEvent, Int_t lineM)
{
  // method in $ALICE_ROOT/MUON/AliMUONTrackLight.cxx 

  if (lineM<0) return;
  Int_t countP=-1, pdg=0;
  TParticle *mother = 0;
  while(lineM>=0){
    mother = ((AliMCParticle*)mcEvent->GetTrack(lineM))->Particle();
    pdg = mother->GetPdgCode();
    fQuarkIndex[++countP] = lineM;
    fQuarkPDGCode[countP] = pdg;
    lineM = mother->GetFirstMother();
  }

  // for PYTHIA checking
  countP = 1;
  for(Int_t par=0; par<4; par++) {
    if(TMath::Abs(this->QuarkPDGCode(par))<6) { countP=par; break; }
  }
  if(this->QuarkIndex(countP)>-1 && (this->ParentFlavour(0)==4 || this->ParentFlavour(0)==5)) {
    if(this->ParentFlavour(0)!=TMath::Abs(this->QuarkPDGCode(countP))) {
      AliWarning(Form("quark flavour of parent and that of quark do not correspond: %d %d --> correcting\n",
                 this->ParentFlavour(0), TMath::Abs(this->QuarkPDGCode(countP))));

      pdg = this->QuarkPDGCode(countP);
      Int_t line = this->QuarkIndex(countP);
      this->ResetQuarkInfo();
      while(TMath::Abs(pdg)!=this->ParentFlavour(0)) {
        pdg = ((AliMCParticle*)mcEvent->GetTrack(++lineM))->Particle()->GetPdgCode();
      }
      while(line>=0){
        mother = ((AliMCParticle*)mcEvent->GetTrack(++lineM))->Particle();
        pdg = mother->GetPdgCode();
        fQuarkIndex[countP] = line;
        fQuarkPDGCode[countP++] = pdg;
        line = mother->GetFirstMother();
      }
    }
  }
  return;
}

//-----------------------------------------------------------------------------
Int_t AliMuonInfoStoreMC::SelectHFMuon()
{
  // set info of muon from HF

  Int_t flv = ParentFlavour(0);
  if (flv!=4 && flv!=5) return 2;

  Bool_t isRes = kFALSE;
  Int_t i=0, nparents=this->ParentsN();
  while (i<nparents && !isRes) {
    isRes = IsMotherAResonance(i++);
  }

  if (isRes) return 2;
  if (flv==5) return 0;
  else return 1;
}

//-----------------------------------------------------------------------------
Bool_t AliMuonInfoStoreMC::IsDiquark(Int_t pdg)
{
  // copy from $ALICE_ROOT/MUON/AliMUONTrackLight.cxx
  pdg = TMath::Abs(pdg);
  if(pdg>1000 && (pdg%100)<10) return kTRUE;
  else return kFALSE;
}

//-----------------------------------------------------------------------------
void AliMuonInfoStoreMC::ResetQuarkInfo()
{
  // copy from $ALICE_ROOT/MUON/AliMUONTrackLight.cxx
  for(Int_t pos=1; pos<4; pos++) {
    fQuarkIndex[pos] = -1;
    fQuarkPDGCode[pos] = 0;
  }
  return;
}

//-----------------------------------------------------------------------------
Int_t AliMuonInfoStoreMC::ParentFlavour(Int_t i) const
{
  // copy from $ALICE_ROOT/MUON/AliMUONTrackLight.cxx
  Int_t pdg = ParentPDGCode(i);
  pdg = TMath::Abs(pdg/100);
  if(pdg>9) pdg /= 10;
  return pdg;
}

//-----------------------------------------------------------------------------
Bool_t AliMuonInfoStoreMC::IsMotherAResonance(Int_t i) const
{
  // copy from $ALICE_ROOT/MUON/AliMUONTrackLight.cxx
  Int_t pdg = ParentPDGCode(i);
  Int_t id=pdg%100000;
  return (!((id-id%10)%110));
}
Commit	Line	Data
fd1d0cb9	1	/**************************************************************************
	2	* Copyright(c) 1998-2006, 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	//
	18	// class used to extract and store info of MC particle
	19	//
	20	// Author: X-M. Zhang, zhang@clermont.in2p3.fr
	21	// zhangxm@iopp.ccnu.edu.cn
	22	/////////////////////////////////////////////////////////////
	23
	24	#include <TParticle.h>
fd1d0cb9	25
fd1d0cb9	26	#include "AliMCEvent.h"
fd1d0cb9	27	#include "AliAODMCParticle.h"
	28	#include "AliESDMuonTrack.h"
	29	#include "AliAODTrack.h"
	30	#include "AliMuonInfoStoreRD.h"
	31	#include "AliMuonInfoStoreMC.h"
1195bb6f	32	#include "AliMCEvent.h"
fd1d0cb9	33
	34	ClassImp(AliMuonInfoStoreMC)
	35
	36	const TString AliMuonInfoStoreMC::fgkStdBranchName("MuonMC");
1195bb6f	37	const Int_t AliMuonInfoStoreMC::fgkSourcesN = 6;
fd1d0cb9	38
	39	//-----------------------------------------------------------------------------
	40	AliMuonInfoStoreMC::AliMuonInfoStoreMC() :
	41	AliMuonInfoStoreRD(),
	42	fIsFull(kFALSE),
	43	fLorentzP(),
	44	fTrackIndex(-1),
	45	fTrackPDGCode(0),
	46	fSource(-1),
	47	fNParents(0),
	48	fOscillation(kFALSE),
	49	fWeight(0.)
	50	{
	51	//
	52	// default constructor
	53	//
	54	for (Int_t i=5; i--;) { fParentIndex[i] = -1; fParentPDGCode[i] = 0; }
	55	for (Int_t i=4; i--;) { fQuarkIndex[i] = -1; fQuarkPDGCode[i] = 0; }
	56	}
	57
	58	//-----------------------------------------------------------------------------
1195bb6f	59	AliMuonInfoStoreMC::AliMuonInfoStoreMC(AliAODTrack trkAOD, AliMCEvent mcEvent, Bool_t full) :
fd1d0cb9	60	AliMuonInfoStoreRD(trkAOD),
	61	fIsFull(full),
	62	fLorentzP(),
	63	fTrackIndex(-1),
	64	fTrackPDGCode(0),
	65	fSource(-1),
	66	fNParents(0),
	67	fOscillation(kFALSE),
	68	fWeight(0.)
	69	{
	70	//
	71	// default constructor
	72	//
	73	for (Int_t i=5; i--;) { fParentIndex[i] = -1; fParentPDGCode[i] = 0; }
	74	for (Int_t i=4; i--;) { fQuarkIndex[i] = -1; fQuarkPDGCode[i] = 0; }
	75
1195bb6f	76	this->SetMCInfoAOD(mcEvent, trkAOD->GetLabel());
fd1d0cb9	77	}
	78
	79	//-----------------------------------------------------------------------------
1195bb6f	80	AliMuonInfoStoreMC::AliMuonInfoStoreMC(AliESDMuonTrack trkESD, AliMCEvent mcEvent, Bool_t full) :
fd1d0cb9	81	AliMuonInfoStoreRD(trkESD),
	82	fIsFull(full),
	83	fLorentzP(),
	84	fTrackIndex(-1),
	85	fTrackPDGCode(0),
	86	fSource(-1),
	87	fNParents(0),
	88	fOscillation(kFALSE),
	89	fWeight(0.)
	90	{
	91	//
	92	// default constructor
	93	//
	94	for (Int_t i=5; i--;) { fParentIndex[i] = -1; fParentPDGCode[i] = 0; }
	95	for (Int_t i=4; i--;) { fQuarkIndex[i] = -1; fQuarkPDGCode[i] = 0; }
	96
1195bb6f	97	this->SetMCInfoESD(mcEvent, trkESD->GetLabel());
fd1d0cb9	98	}
fd1d0cb9	99
fd1d0cb9	100	//-----------------------------------------------------------------------------
	101	AliMuonInfoStoreMC::AliMuonInfoStoreMC(const AliMuonInfoStoreMC &src) :
	102	AliMuonInfoStoreRD(src),
	103	fIsFull(src.fIsFull),
	104	fLorentzP(src.fLorentzP),
	105	fTrackIndex(src.fTrackIndex),
	106	fTrackPDGCode(src.fTrackPDGCode),
	107	fSource(src.fSource),
	108	fNParents(src.fNParents),
	109	fOscillation(src.fOscillation),
	110	fWeight(src.fWeight)
	111	{
	112	//
	113	// copy constructor
	114	//
	115	for (Int_t i=5; i--;) {
	116	fParentIndex[i] = src.fParentIndex[i];
	117	fParentPDGCode[i] = src.fParentPDGCode[i];
	118	}
	119	for (Int_t i=4; i--;) {
	120	fQuarkIndex[i] = src.fQuarkIndex[i];
	121	fQuarkPDGCode[i] = src.fQuarkPDGCode[i];
	122	}
	123	}
	124
	125	//-----------------------------------------------------------------------------
	126	AliMuonInfoStoreMC& AliMuonInfoStoreMC::operator=(const AliMuonInfoStoreMC &src)
	127	{
	128	//
	129	// assignment constructor
	130	//
	131	if(&src==this) return *this;
	132	AliMuonInfoStoreRD::operator=(src);
	133
	134	fIsFull = src.fIsFull;
	135	fLorentzP = src.fLorentzP;
	136	fTrackIndex = src.fTrackIndex;
	137	fTrackPDGCode = src.fTrackPDGCode;
	138	fSource = src.fSource;
	139	fNParents = src.fNParents;
	140	fOscillation = src.fOscillation;
	141	fWeight = src.fWeight;
	142	for (Int_t i=5; i--;) {
	143	fParentIndex[i] = src.fParentIndex[i];
	144	fParentPDGCode[i] = src.fParentPDGCode[i];
	145	}
	146	for (Int_t i=4; i--;) {
	147	fQuarkIndex[i] = src.fQuarkIndex[i];
	148	fQuarkPDGCode[i] = src.fQuarkPDGCode[i];
	149	}
	150
	151	return *this;
	152	}
	153
	154	//-----------------------------------------------------------------------------
	155	AliMuonInfoStoreMC::~AliMuonInfoStoreMC()
	156	{
	157	//
	158	// destructor
	159	//
	160	}
	161
	162	//-----------------------------------------------------------------------------
1195bb6f	163	void AliMuonInfoStoreMC::SetMCInfoAOD(AliMCEvent *mcEvent, Int_t label)
fd1d0cb9	164	{
fd1d0cb9	165	// fill track MC info with AOD base
1195bb6f	166	fTrackIndex = label;
1195bb6f	167	if (fTrackIndex<0) { fSource=5; return; }
fd1d0cb9	168
1195bb6f	169	AliAODMCParticle pMC = (AliAODMCParticle)mcEvent->GetTrack(fTrackIndex);
fd1d0cb9	170	fLorentzP.SetPxPyPzE(pMC->Px(), pMC->Py(), pMC->Pz(), pMC->E());
fd1d0cb9	171
1195bb6f	172	fTrackPDGCode = pMC->PdgCode();
1195bb6f	173	if (TMath::Abs(fTrackPDGCode)!=13) { fSource=4; return; }
fd1d0cb9	174
1195bb6f	175	Int_t lineM = pMC->GetMother();
1195bb6f	176	if (lineM<0) { fSource=2; return; }
fd1d0cb9	177
1195bb6f	178	Bool_t isPrimary = ((AliAODMCParticle*)mcEvent->GetTrack(lineM))->IsPrimary();
1195bb6f	179	if (!isPrimary) { fSource=3; return; }
fd1d0cb9	180
	181	Int_t countP=-1, pdg=0;
	182	Int_t parents[10], parLine[10];
	183	AliAODMCParticle *mother = 0;
	184	while(lineM>=0){
1195bb6f	185	mother = (AliAODMCParticle*)mcEvent->GetTrack(lineM);
fd1d0cb9	186	pdg = mother->GetPdgCode();
	187	if(pdg==92 \|\| pdg==21 \|\| TMath::Abs(pdg)<10 \|\| IsDiquark(pdg)) break;
	188	parents[++countP] = pdg;
	189	parLine[countP] = lineM;
	190	lineM = mother->GetMother();
	191	}
	192	for(Int_t i=0; i<=countP; i++){
	193	fParentIndex[i] = parLine[countP-i];
	194	fParentPDGCode[i] = parents[countP-i];
	195	}
	196	fNParents = countP + 1;
	197
1195bb6f	198	if (fIsFull && lineM>=0) this->FillHistoryQuarksAOD(mcEvent, lineM);
	199
	200	fSource = this->SelectHFMuon();
fd1d0cb9	201	return;
	202	}
	203
	204	//-----------------------------------------------------------------------------
1195bb6f	205	void AliMuonInfoStoreMC::SetMCInfoESD(AliMCEvent *mcEvent, Int_t label)
fd1d0cb9	206	{
1195bb6f	207	// fill track MC info with ESD base
	208	fTrackIndex = label;
	209	if (fTrackIndex<0) { fSource=5; return; }
	210
	211	TParticle pMC = ((AliMCParticle)mcEvent->GetTrack(fTrackIndex))->Particle();
	212	fLorentzP.SetPxPyPzE(pMC->Px(), pMC->Py(), pMC->Pz(), pMC->Energy());
	213
	214	fTrackPDGCode = pMC->GetPdgCode();
	215	if (TMath::Abs(fTrackPDGCode)!=13) { fSource=4; return; }
	216
	217	Int_t lineM = pMC->GetFirstMother();
	218	if (lineM<0) { fSource=2; return; }
	219
	220	if (lineM>=mcEvent->Stack()->GetNprimary()) { fSource=3; return; }
fd1d0cb9	221
	222	Int_t countP=-1, pdg=0;
	223	Int_t parents[10], parLine[10];
	224	TParticle *mother = 0;
	225	while(lineM>=0){
1195bb6f	226	mother = ((AliMCParticle*)mcEvent->GetTrack(lineM))->Particle();
fd1d0cb9	227	pdg = mother->GetPdgCode();
	228	if(pdg==92 \|\| pdg==21 \|\| TMath::Abs(pdg)<10 \|\| IsDiquark(pdg)) break;
	229	parents[++countP] = pdg;
	230	parLine[countP] = lineM;
	231	lineM = mother->GetFirstMother();
	232	}
	233	for(Int_t i=0; i<=countP; i++){
	234	fParentIndex[i] = parLine[countP-i];
	235	fParentPDGCode[i] = parents[countP-i];
	236	}
	237	fNParents = countP + 1;
	238
1195bb6f	239	if (fIsFull && lineM>=0) this->FillHistoryQuarksESD(mcEvent, lineM);
	240
	241	fSource = this->SelectHFMuon();
fd1d0cb9	242	return;
	243	}
	244
	245	//-----------------------------------------------------------------------------
1195bb6f	246	void AliMuonInfoStoreMC::FillHistoryQuarksAOD(AliMCEvent* const mcEvent, Int_t lineM)
fd1d0cb9	247	{
	248	// method in $ALICE_ROOT/MUON/AliMUONTrackLight.cxx
	249
	250	if (lineM<0) return;
	251	Int_t countP=-1, pdg=0;
	252	AliAODMCParticle *mother = 0;
	253	while(lineM>=0){
1195bb6f	254	mother = (AliAODMCParticle*)mcEvent->GetTrack(lineM);
fd1d0cb9	255	pdg = mother->GetPdgCode();
	256	fQuarkIndex[++countP] = lineM;
	257	fQuarkPDGCode[countP] = pdg;
	258	lineM = mother->GetMother();
	259	}
	260
	261	// for PYTHIA checking
	262	countP = 1;
	263	for(Int_t par=0; par<4; par++) {
	264	if(TMath::Abs(this->QuarkPDGCode(par))<6) { countP=par; break; }
	265	}
	266	if(this->QuarkIndex(countP)>-1 && (this->ParentFlavour(0)==4 \|\| this->ParentFlavour(0)==5)) {
	267	if(this->ParentFlavour(0)!=TMath::Abs(this->QuarkPDGCode(countP))) {
	268	AliWarning(Form("quark flavour of parent and that of quark do not correspond: %d %d --> correcting\n",
	269	this->ParentFlavour(0), TMath::Abs(this->QuarkPDGCode(countP))));
	270
	271	pdg = this->QuarkPDGCode(countP);
	272	Int_t line = this->QuarkIndex(countP);
	273	this->ResetQuarkInfo();
	274	while(TMath::Abs(pdg)!=this->ParentFlavour(0)) {
1195bb6f	275	pdg = ((AliAODMCParticle*)mcEvent->GetTrack(++line))->GetPdgCode();
fd1d0cb9	276	}
fd1d0cb9	277	while(line>=0){
1195bb6f	278	mother = (AliAODMCParticle*)mcEvent->GetTrack(line);
fd1d0cb9	279	pdg = mother->GetPdgCode();
	280	fQuarkIndex[countP] = line;
	281	fQuarkPDGCode[countP++] = pdg;
	282	line = mother->GetMother();
	283	}
	284	}
	285	}
	286	return;
	287	}
	288
	289	//-----------------------------------------------------------------------------
1195bb6f	290	void AliMuonInfoStoreMC::FillHistoryQuarksESD(AliMCEvent* const mcEvent, Int_t lineM)
fd1d0cb9	291	{
	292	// method in $ALICE_ROOT/MUON/AliMUONTrackLight.cxx
	293
	294	if (lineM<0) return;
	295	Int_t countP=-1, pdg=0;
	296	TParticle *mother = 0;
	297	while(lineM>=0){
1195bb6f	298	mother = ((AliMCParticle*)mcEvent->GetTrack(lineM))->Particle();
fd1d0cb9	299	pdg = mother->GetPdgCode();
	300	fQuarkIndex[++countP] = lineM;
	301	fQuarkPDGCode[countP] = pdg;
	302	lineM = mother->GetFirstMother();
	303	}
	304
	305	// for PYTHIA checking
	306	countP = 1;
	307	for(Int_t par=0; par<4; par++) {
	308	if(TMath::Abs(this->QuarkPDGCode(par))<6) { countP=par; break; }
	309	}
	310	if(this->QuarkIndex(countP)>-1 && (this->ParentFlavour(0)==4 \|\| this->ParentFlavour(0)==5)) {
	311	if(this->ParentFlavour(0)!=TMath::Abs(this->QuarkPDGCode(countP))) {
	312	AliWarning(Form("quark flavour of parent and that of quark do not correspond: %d %d --> correcting\n",
	313	this->ParentFlavour(0), TMath::Abs(this->QuarkPDGCode(countP))));
	314
	315	pdg = this->QuarkPDGCode(countP);
	316	Int_t line = this->QuarkIndex(countP);
	317	this->ResetQuarkInfo();
	318	while(TMath::Abs(pdg)!=this->ParentFlavour(0)) {
1195bb6f	319	pdg = ((AliMCParticle*)mcEvent->GetTrack(++lineM))->Particle()->GetPdgCode();
fd1d0cb9	320	}
fd1d0cb9	321	while(line>=0){
1195bb6f	322	mother = ((AliMCParticle*)mcEvent->GetTrack(++lineM))->Particle();
fd1d0cb9	323	pdg = mother->GetPdgCode();
	324	fQuarkIndex[countP] = line;
	325	fQuarkPDGCode[countP++] = pdg;
	326	line = mother->GetFirstMother();
	327	}
	328	}
	329	}
	330	return;
	331	}
	332
	333	//-----------------------------------------------------------------------------
	334	Int_t AliMuonInfoStoreMC::SelectHFMuon()
	335	{
	336	// set info of muon from HF
	337
	338	Int_t flv = ParentFlavour(0);
	339	if (flv!=4 && flv!=5) return 2;
	340
	341	Bool_t isRes = kFALSE;
1195bb6f	342	Int_t i=0, nparents=this->ParentsN();
fd1d0cb9	343	while (i<nparents && !isRes) {
	344	isRes = IsMotherAResonance(i++);
	345	}
	346
	347	if (isRes) return 2;
	348	if (flv==5) return 0;
	349	else return 1;
	350	}
	351
	352	//-----------------------------------------------------------------------------
	353	Bool_t AliMuonInfoStoreMC::IsDiquark(Int_t pdg)
	354	{
	355	// copy from $ALICE_ROOT/MUON/AliMUONTrackLight.cxx
	356	pdg = TMath::Abs(pdg);
	357	if(pdg>1000 && (pdg%100)<10) return kTRUE;
	358	else return kFALSE;
	359	}
	360
	361	//-----------------------------------------------------------------------------
	362	void AliMuonInfoStoreMC::ResetQuarkInfo()
	363	{
	364	// copy from $ALICE_ROOT/MUON/AliMUONTrackLight.cxx
	365	for(Int_t pos=1; pos<4; pos++) {
	366	fQuarkIndex[pos] = -1;
	367	fQuarkPDGCode[pos] = 0;
	368	}
	369	return;
	370	}
	371
	372	//-----------------------------------------------------------------------------
	373	Int_t AliMuonInfoStoreMC::ParentFlavour(Int_t i) const
	374	{
	375	// copy from $ALICE_ROOT/MUON/AliMUONTrackLight.cxx
	376	Int_t pdg = ParentPDGCode(i);
	377	pdg = TMath::Abs(pdg/100);
	378	if(pdg>9) pdg /= 10;
	379	return pdg;
	380	}
	381
	382	//-----------------------------------------------------------------------------
	383	Bool_t AliMuonInfoStoreMC::IsMotherAResonance(Int_t i) const
	384	{
	385	// copy from $ALICE_ROOT/MUON/AliMUONTrackLight.cxx
	386	Int_t pdg = ParentPDGCode(i);
	387	Int_t id=pdg%100000;
	388	return (!((id-id%10)%110));
	389	}