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

#include <TParticle.h>
#include <TLorentzVector.h>
#include <TClonesArray.h>

#include "AliMCEventHandler.h"
#include "AliMCEvent.h"
#include "AliStack.h"
#include "AliAODMCParticle.h"
#include "AliESDEvent.h"
#include "AliESDMuonTrack.h"
#include "AliESDMuonCluster.h"
#include "AliAODTrack.h"
#include "AliAODMuonTrack.h"
#include "AliMCMuonTrack.h"

ClassImp(AliMCMuonTrack)

const Double_t AliMCMuonTrack::fgkSigmaCut = 10.;

//-----------------------------------------------------------------------------
AliMCMuonTrack::AliMCMuonTrack() :
AliAODMuonTrack(),
fIsFull(kFALSE),
fPGen(),
fTrackIndex(-1),
fTrackPDGCode(0),
fSource(-1),
fNParents(0),
fOscillation(kFALSE),
fWeight(0.)
{
  //
  // default constructor
  //
  for (Int_t i=0; i<fgkNParentsMax; i++) {
    fParentIndex[i] = -1;
    fParentPDGCode[i] = 0;
  }
  for (Int_t i=0; i<4; i++) {
    fQuarkIndex[i] = -1;
    fQuarkPDGCode[i] = 0;
  }
}

//-----------------------------------------------------------------------------
AliMCMuonTrack::AliMCMuonTrack(AliAODTrack *trkAOD, TClonesArray *mcClArr, Bool_t full) :
AliAODMuonTrack(trkAOD),
fIsFull(full),
fPGen(),
fTrackIndex(-1),
fTrackPDGCode(0),
fSource(-1),
fNParents(0),
fOscillation(kFALSE),
fWeight(0.)
{
  //
  // default constructor
  //
  AliAODMCParticle *pMC = this->FindTrackRef(trkAOD, mcClArr);
  if (pMC) this->SetMCInfo(pMC, mcClArr);
}

//-----------------------------------------------------------------------------
AliMCMuonTrack::AliMCMuonTrack(AliESDMuonTrack *trkESD, AliESDEvent *esd, AliMCEventHandler *mcH, Bool_t full) :
AliAODMuonTrack(trkESD),
fIsFull(full),
fPGen(),
fTrackIndex(-1),
fTrackPDGCode(0),
fSource(-1),
fNParents(0),
fOscillation(kFALSE),
fWeight(0.)
{
  //
  // default constructor
  //
  TParticle *pMC = this->FindTrackRef(trkESD, esd, mcH);  // do nothing when running on official train
  if (pMC) this->SetMCInfo(pMC, mcH);
  // MC infomation is not implement when running on the official train for they depend on the MUON module.
  // In the case of private running, one can get the information from ESD by uncommenting the definitions
  // of FindTrackRef() and CovESDtoMuonTrack()
  // and including the following head files in both AliMCMuonTrack.h and this file,
  // #include "AliMUONRecoCheck.h"
  // #include "AliMUONVTrackStore.h"
  // #include "AliMUONVClusterStore.h"
  // #include "AliMUONTrackParam.h"
  // #include "AliMUONTrack.h"
  // #include "AliMUONVCluster.h"
  // #include "AliMUONESDInterface.h"
}


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

//-----------------------------------------------------------------------------
AliAODMCParticle* AliMCMuonTrack::FindTrackRef(AliAODTrack *trkAOD, TClonesArray *mcClArr)
{
  AliAODMCParticle *pMC = 0;
  fTrackIndex = trkAOD->GetLabel();
  if (fTrackIndex>=0)
    pMC = (AliAODMCParticle*)mcClArr->At(fTrackIndex);
  return pMC;
}

//-----------------------------------------------------------------------------
TParticle* AliMCMuonTrack::FindTrackRef(AliESDMuonTrack *trkESD, AliESDEvent *esd, AliMCEventHandler *mcH)
{
  TParticle *pMCRef = 0;

  /*AliMUONTrack *trkMuon = CovESDtoMuonTrack(*trkESD);
  AliMUONRecoCheck rc(esd,mcH);
  AliMUONVTrackStore* trkRefArr = rc.TrackRefs(-1);
  TIter next(trkRefArr->CreateIterator());
  AliMUONTrack* trkRef = 0;
  while ((trkRef=static_cast<AliMUONTrack*>(next()))) {
    Bool_t trkCompArr[10];
    Int_t nMatchClusters = trkMuon->CompatibleTrack(trkRef, fgkSigmaCut, trkCompArr);
    Double_t matchClusterFrac = ((Double_t)nMatchClusters) / ((Double_t)trkMuon->GetNClusters());
    if ((trkCompArr[0] || trkCompArr[1] || trkCompArr[2] || trkCompArr[3]) &&
        (trkCompArr[6] || trkCompArr[7] || trkCompArr[8] || trkCompArr[9]) && matchClusterFrac>0.5) {
      fTrackIndex = trkRef->GetUniqueID();
      trkRefArr->Remove(*trkRef);
      break;
    }
  }
  if (fTrackIndex>=0) pMCRef = mcH->MCEvent()->Stack()->Particle(fTrackIndex);*/

  return pMCRef;
}

//-----------------------------------------------------------------------------
void AliMCMuonTrack::SetMCInfo(AliAODMCParticle *pMC, TClonesArray *mcClArr)
{
  fPGen.SetPxPyPzE(pMC->Px(), pMC->Py(), pMC->Pz(), pMC->E());
  fTrackPDGCode = pMC->GetPdgCode();
  if (TMath::Abs(fTrackPDGCode)!=13) {
    fSource = 4;
    return;
  } 

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

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

  this->FillHistoryParents(lineM, mcClArr);
  fSource = this->SelectHFMuon();
  return;
}

//-----------------------------------------------------------------------------
void AliMCMuonTrack::SetMCInfo(TParticle *pMC, AliMCEventHandler *mcH)
{
  fPGen.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;
  }

  AliStack *stack = mcH->MCEvent()->Stack();
  if (lineM>=stack->GetNprimary()) {
    fSource = 3;
    return;
  }

  this->FillHistoryParents(lineM, stack);
  fSource = this->SelectHFMuon();
  return;
}

//-----------------------------------------------------------------------------
void AliMCMuonTrack::FillHistoryParents(Int_t lineM, TClonesArray *mcClArr)
{
  Int_t countP=-1, pdg=0;
  Int_t parents[10], parLine[10];
  AliAODMCParticle *mother = 0;
  while(lineM>=0){
    mother = (AliAODMCParticle*)mcClArr->At(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->FillHistoryQuarks(lineM, mcClArr);
  return;
}

//-----------------------------------------------------------------------------
void AliMCMuonTrack::FillHistoryParents(Int_t lineM, AliStack *stack)
{
  Int_t countP=-1, pdg=0;
  Int_t parents[10], parLine[10];
  TParticle *mother = 0;
  while(lineM>=0){
    mother = stack->Particle(lineM);
    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->FillHistoryQuarks(lineM, stack);
  return;
}

//-----------------------------------------------------------------------------
void AliMCMuonTrack::FillHistoryQuarks(Int_t lineM, TClonesArray *mcClArr)
{
  // 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*)mcClArr->At(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->GetQuarkPDGCode(par))<6) { countP=par; break; }
  }
  if(this->GetQuarkIndex(countP)>-1 && (this->GetParentFlavour(0)==4 || this->GetParentFlavour(0)==5)) {
    if(this->GetParentFlavour(0)!=TMath::Abs(this->GetQuarkPDGCode(countP))) {
      AliWarning(Form("quark flavour of parent and that of quark do not correspond: %d %d --> correcting\n",
                 this->GetParentFlavour(0), TMath::Abs(this->GetQuarkPDGCode(countP))));

      pdg = this->GetQuarkPDGCode(countP);
      Int_t line = this->GetQuarkIndex(countP);
      this->ResetQuarkInfo();
      while(TMath::Abs(pdg)!=this->GetParentFlavour(0)) {
        pdg = ((AliAODMCParticle*)mcClArr->At(++line))->GetPdgCode();
      }
      while(line>=0){
        mother = (AliAODMCParticle*)mcClArr->At(line);
        pdg = mother->GetPdgCode();
        fQuarkIndex[countP] = line;
        fQuarkPDGCode[countP++] = pdg;
        line = mother->GetMother();
      }
    }
  }
  return;
}

//-----------------------------------------------------------------------------
void AliMCMuonTrack::FillHistoryQuarks(Int_t lineM, AliStack *stack)
{
  // method in $ALICE_ROOT/MUON/AliMUONTrackLight.cxx 
  if (lineM<0) return;
  Int_t countP=-1, pdg=0;
  TParticle *mother = 0;
  while(lineM>=0){
    mother = stack->Particle(lineM);
    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->GetQuarkPDGCode(par))<6) { countP=par; break; }
  }
  if(this->GetQuarkIndex(countP)>-1 && (this->GetParentFlavour(0)==4 || this->GetParentFlavour(0)==5)) {
    if(this->GetParentFlavour(0)!=TMath::Abs(this->GetQuarkPDGCode(countP))) {
      AliWarning(Form("quark flavour of parent and that of quark do not correspond: %d %d --> correcting\n",
                 this->GetParentFlavour(0), TMath::Abs(this->GetQuarkPDGCode(countP))));

      pdg = this->GetQuarkPDGCode(countP);
      Int_t line = this->GetQuarkIndex(countP);
      this->ResetQuarkInfo();
      while(TMath::Abs(pdg)!=this->GetParentFlavour(0)) {
        pdg = stack->Particle(++line)->GetPdgCode();
      }
      while(line>=0){
        mother = stack->Particle(line);
        pdg = mother->GetPdgCode();
        fQuarkIndex[countP] = line;
        fQuarkPDGCode[countP++] = pdg;
        line = mother->GetFirstMother();
      }
    }
  }
  return;
}

//-----------------------------------------------------------------------------
/*AliMUONTrack* AliMCMuonTrack::CovESDtoMuonTrack(AliESDMuonTrack &trkESD)
{
  // method in $ALICE_ROOT/PWG3/muondep/AliAnalysisTaskESDMCLabelAddition.cxx
  AliMUONTrack *trkMuon = new AliMUONTrack();
  if(!trkESD.ClustersStored()) return trkMuon;

  AliMUONTrackParam param;
  AliMUONESDInterface::GetParamAtFirstCluster(trkESD, param);
  AliMUONESDInterface::GetParamCov(trkESD, param);
  AliMUONVClusterStore *cStore = AliMUONESDInterface::NewClusterStore();
  AliMUONVCluster *cluster = cStore->CreateCluster(0,0,0);

  AliESDMuonCluster *esdCluster = (AliESDMuonCluster*)trkESD.GetClusters().First();
  while (esdCluster) {
    AliMUONESDInterface::ESDToMUON(*esdCluster, *cluster);
    param.SetZ(cluster->GetZ());
    trkMuon->AddTrackParamAtCluster(param, *cluster, kTRUE);
    esdCluster = (AliESDMuonCluster*)trkESD.GetClusters().After(esdCluster);
  }

  delete cluster; cluster = 0;
  delete cStore;  cStore = 0;
  return trkMuon;
}*/

//-----------------------------------------------------------------------------
Int_t AliMCMuonTrack::SelectHFMuon()
{
  Int_t flv = GetParentFlavour(0);
  if (flv!=4 && flv!=5) return 2;

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

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

//-----------------------------------------------------------------------------
Bool_t AliMCMuonTrack::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 AliMCMuonTrack::ResetQuarkInfo()
{
  for(Int_t pos=1; pos<4; pos++) {
    fQuarkIndex[pos] = -1;
    fQuarkPDGCode[pos] = 0;
  }
  return;
}

//-----------------------------------------------------------------------------
Int_t AliMCMuonTrack::GetParentFlavour(Int_t i) const
{
  Int_t pdg = GetParentPDGCode(i);
  pdg = TMath::Abs(pdg/100);
  if(pdg>9) pdg /= 10;
  return pdg;
}

//-----------------------------------------------------------------------------
Bool_t AliMCMuonTrack::IsMotherAResonance(Int_t i) const
{
  // copy from $ALICE_ROOT/MUON/AliMUONTrackLight.cxx
  Int_t pdg = GetParentPDGCode(i);
  Int_t id=pdg%100000;
  return (!((id-id%10)%110));
}
Commit	Line	Data
4292b3b6	1	#include <TParticle.h>
	2	#include <TLorentzVector.h>
	3	#include <TClonesArray.h>
	4
	5	#include "AliMCEventHandler.h"
	6	#include "AliMCEvent.h"
	7	#include "AliStack.h"
	8	#include "AliAODMCParticle.h"
	9	#include "AliESDEvent.h"
	10	#include "AliESDMuonTrack.h"
	11	#include "AliESDMuonCluster.h"
	12	#include "AliAODTrack.h"
	13	#include "AliAODMuonTrack.h"
	14	#include "AliMCMuonTrack.h"
	15
	16	ClassImp(AliMCMuonTrack)
	17
	18	const Double_t AliMCMuonTrack::fgkSigmaCut = 10.;
	19
	20	//-----------------------------------------------------------------------------
	21	AliMCMuonTrack::AliMCMuonTrack() :
	22	AliAODMuonTrack(),
	23	fIsFull(kFALSE),
	24	fPGen(),
	25	fTrackIndex(-1),
	26	fTrackPDGCode(0),
	27	fSource(-1),
	28	fNParents(0),
	29	fOscillation(kFALSE),
	30	fWeight(0.)
	31	{
	32	//
	33	// default constructor
	34	//
	35	for (Int_t i=0; i<fgkNParentsMax; i++) {
	36	fParentIndex[i] = -1;
	37	fParentPDGCode[i] = 0;
	38	}
	39	for (Int_t i=0; i<4; i++) {
	40	fQuarkIndex[i] = -1;
	41	fQuarkPDGCode[i] = 0;
	42	}
	43	}
	44
	45	//-----------------------------------------------------------------------------
	46	AliMCMuonTrack::AliMCMuonTrack(AliAODTrack trkAOD, TClonesArray mcClArr, Bool_t full) :
	47	AliAODMuonTrack(trkAOD),
	48	fIsFull(full),
	49	fPGen(),
	50	fTrackIndex(-1),
	51	fTrackPDGCode(0),
	52	fSource(-1),
	53	fNParents(0),
	54	fOscillation(kFALSE),
	55	fWeight(0.)
	56	{
	57	//
	58	// default constructor
	59	//
	60	AliAODMCParticle *pMC = this->FindTrackRef(trkAOD, mcClArr);
	61	if (pMC) this->SetMCInfo(pMC, mcClArr);
	62	}
	63
	64	//-----------------------------------------------------------------------------
65	AliMCMuonTrack::AliMCMuonTrack(AliESDMuonTrack trkESD, AliESDEvent esd, AliMCEventHandler *mcH, Bool_t full) :
66	AliAODMuonTrack(trkESD),
67	fIsFull(full),
68	fPGen(),
69	fTrackIndex(-1),
70	fTrackPDGCode(0),
71	fSource(-1),
72	fNParents(0),
73	fOscillation(kFALSE),
74	fWeight(0.)
75	{
76	//
77	// default constructor
78	//
79	TParticle *pMC = this->FindTrackRef(trkESD, esd, mcH); // do nothing when running on official train
80	if (pMC) this->SetMCInfo(pMC, mcH);
81	// MC infomation is not implement when running on the official train for they depend on the MUON module.
82	// In the case of private running, one can get the information from ESD by uncommenting the definitions
83	// of FindTrackRef() and CovESDtoMuonTrack()
84	// and including the following head files in both AliMCMuonTrack.h and this file,
85	// #include "AliMUONRecoCheck.h"
86	// #include "AliMUONVTrackStore.h"
87	// #include "AliMUONVClusterStore.h"
88	// #include "AliMUONTrackParam.h"
89	// #include "AliMUONTrack.h"
90	// #include "AliMUONVCluster.h"
91	// #include "AliMUONESDInterface.h"
92	}
93
94
95
96	//-----------------------------------------------------------------------------
97	AliMCMuonTrack::~AliMCMuonTrack()
98	{
99	//
100	// destructor
101	//
102	}
103
104	//-----------------------------------------------------------------------------
105	AliAODMCParticle* AliMCMuonTrack::FindTrackRef(AliAODTrack trkAOD, TClonesArray mcClArr)
106	{
107	AliAODMCParticle *pMC = 0;
108	fTrackIndex = trkAOD->GetLabel();
109	if (fTrackIndex>=0)
110	pMC = (AliAODMCParticle*)mcClArr->At(fTrackIndex);
111	return pMC;
112	}
113
114	//-----------------------------------------------------------------------------
115	TParticle* AliMCMuonTrack::FindTrackRef(AliESDMuonTrack trkESD, AliESDEvent esd, AliMCEventHandler *mcH)
116	{
117	TParticle *pMCRef = 0;
118
119	/AliMUONTrack trkMuon = CovESDtoMuonTrack(*trkESD);
120	AliMUONRecoCheck rc(esd,mcH);
121	AliMUONVTrackStore* trkRefArr = rc.TrackRefs(-1);
122	TIter next(trkRefArr->CreateIterator());
123	AliMUONTrack* trkRef = 0;
124	while ((trkRef=static_cast<AliMUONTrack*>(next()))) {
125	Bool_t trkCompArr[10];
126	Int_t nMatchClusters = trkMuon->CompatibleTrack(trkRef, fgkSigmaCut, trkCompArr);
127	Double_t matchClusterFrac = ((Double_t)nMatchClusters) / ((Double_t)trkMuon->GetNClusters());
128	if ((trkCompArr[0] \|\| trkCompArr[1] \|\| trkCompArr[2] \|\| trkCompArr[3]) &&
129	(trkCompArr[6] \|\| trkCompArr[7] \|\| trkCompArr[8] \|\| trkCompArr[9]) && matchClusterFrac>0.5) {
130	fTrackIndex = trkRef->GetUniqueID();
131	trkRefArr->Remove(*trkRef);
132	break;
133	}
134	}
135	if (fTrackIndex>=0) pMCRef = mcH->MCEvent()->Stack()->Particle(fTrackIndex);*/
136
137	return pMCRef;
138	}
139
140	//-----------------------------------------------------------------------------
141	void AliMCMuonTrack::SetMCInfo(AliAODMCParticle pMC, TClonesArray mcClArr)
142	{
143	fPGen.SetPxPyPzE(pMC->Px(), pMC->Py(), pMC->Pz(), pMC->E());
144	fTrackPDGCode = pMC->GetPdgCode();
145	if (TMath::Abs(fTrackPDGCode)!=13) {
146	fSource = 4;
147	return;
148	}
149
150	Int_t lineM = pMC->GetMother();
151	if (lineM<0) {
152	fSource = 2;
153	return;
154	}
155
156	Bool_t isPrimary =
157	((AliAODMCParticle*)mcClArr->At(lineM))->IsPrimary();
158	if (!isPrimary) {
159	fSource = 3;
160	return;
161	}
162
163	this->FillHistoryParents(lineM, mcClArr);
164	fSource = this->SelectHFMuon();
165	return;
166	}
167
168	//-----------------------------------------------------------------------------
169	void AliMCMuonTrack::SetMCInfo(TParticle pMC, AliMCEventHandler mcH)
170	{
171	fPGen.SetPxPyPzE(pMC->Px(), pMC->Py(), pMC->Pz(), pMC->Energy());
172	fTrackPDGCode = pMC->GetPdgCode();
173	if (TMath::Abs(fTrackPDGCode)!=13) {
174	fSource = 4;
175	return;
176	}
177
178	Int_t lineM = pMC->GetFirstMother();
179	if (lineM<0) {
180	fSource = 2;
181	return;
182	}
183
184	AliStack *stack = mcH->MCEvent()->Stack();
185	if (lineM>=stack->GetNprimary()) {
186	fSource = 3;
187	return;
188	}
189
190	this->FillHistoryParents(lineM, stack);
191	fSource = this->SelectHFMuon();
192	return;
193	}
194
195	//-----------------------------------------------------------------------------
196	void AliMCMuonTrack::FillHistoryParents(Int_t lineM, TClonesArray *mcClArr)
197	{
198	Int_t countP=-1, pdg=0;
199	Int_t parents[10], parLine[10];
200	AliAODMCParticle *mother = 0;
201	while(lineM>=0){
202	mother = (AliAODMCParticle*)mcClArr->At(lineM);
203	pdg = mother->GetPdgCode();
204	if(pdg==92 \|\| pdg==21 \|\| TMath::Abs(pdg)<10 \|\| IsDiquark(pdg)) break;
205	parents[++countP] = pdg;
206	parLine[countP] = lineM;
207	lineM = mother->GetMother();
208	}
209	for(Int_t i=0; i<=countP; i++){
210	fParentIndex[i] = parLine[countP-i];
211	fParentPDGCode[i] = parents[countP-i];
212	}
213	fNParents = countP + 1;
214
215	if (fIsFull && lineM>=0) this->FillHistoryQuarks(lineM, mcClArr);
216	return;
217	}
218
219	//-----------------------------------------------------------------------------
220	void AliMCMuonTrack::FillHistoryParents(Int_t lineM, AliStack *stack)
221	{
222	Int_t countP=-1, pdg=0;
223	Int_t parents[10], parLine[10];
224	TParticle *mother = 0;
225	while(lineM>=0){
226	mother = stack->Particle(lineM);
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
239	if (fIsFull && lineM>=0) this->FillHistoryQuarks(lineM, stack);
240	return;
241	}
242
243	//-----------------------------------------------------------------------------
244	void AliMCMuonTrack::FillHistoryQuarks(Int_t lineM, TClonesArray *mcClArr)
245	{
246	// method in $ALICE_ROOT/MUON/AliMUONTrackLight.cxx
247	if (lineM<0) return;
248	Int_t countP=-1, pdg=0;
249	AliAODMCParticle *mother = 0;
250	while(lineM>=0){
251	mother = (AliAODMCParticle*)mcClArr->At(lineM);
252	pdg = mother->GetPdgCode();
253	fQuarkIndex[++countP] = lineM;
254	fQuarkPDGCode[countP] = pdg;
255	lineM = mother->GetMother();
256	}
257
258	// for PYTHIA checking
259	countP = 1;
260	for(Int_t par=0; par<4; par++) {
261	if(TMath::Abs(this->GetQuarkPDGCode(par))<6) { countP=par; break; }
262	}
263	if(this->GetQuarkIndex(countP)>-1 && (this->GetParentFlavour(0)==4 \|\| this->GetParentFlavour(0)==5)) {
264	if(this->GetParentFlavour(0)!=TMath::Abs(this->GetQuarkPDGCode(countP))) {
265	AliWarning(Form("quark flavour of parent and that of quark do not correspond: %d %d --> correcting\n",
266	this->GetParentFlavour(0), TMath::Abs(this->GetQuarkPDGCode(countP))));
267
268	pdg = this->GetQuarkPDGCode(countP);
269	Int_t line = this->GetQuarkIndex(countP);
270	this->ResetQuarkInfo();
271	while(TMath::Abs(pdg)!=this->GetParentFlavour(0)) {
272	pdg = ((AliAODMCParticle*)mcClArr->At(++line))->GetPdgCode();
273	}
274	while(line>=0){
275	mother = (AliAODMCParticle*)mcClArr->At(line);
276	pdg = mother->GetPdgCode();
277	fQuarkIndex[countP] = line;
278	fQuarkPDGCode[countP++] = pdg;
279	line = mother->GetMother();
280	}
281	}
282	}
283	return;
284	}
285
286	//-----------------------------------------------------------------------------
287	void AliMCMuonTrack::FillHistoryQuarks(Int_t lineM, AliStack *stack)
288	{
289	// method in $ALICE_ROOT/MUON/AliMUONTrackLight.cxx
290	if (lineM<0) return;
291	Int_t countP=-1, pdg=0;
292	TParticle *mother = 0;
293	while(lineM>=0){
294	mother = stack->Particle(lineM);
295	pdg = mother->GetPdgCode();
296	fQuarkIndex[++countP] = lineM;
297	fQuarkPDGCode[countP] = pdg;
298	lineM = mother->GetFirstMother();
299	}
300
301	// for PYTHIA checking
302	countP = 1;
303	for(Int_t par=0; par<4; par++) {
304	if(TMath::Abs(this->GetQuarkPDGCode(par))<6) { countP=par; break; }
305	}
306	if(this->GetQuarkIndex(countP)>-1 && (this->GetParentFlavour(0)==4 \|\| this->GetParentFlavour(0)==5)) {
307	if(this->GetParentFlavour(0)!=TMath::Abs(this->GetQuarkPDGCode(countP))) {
308	AliWarning(Form("quark flavour of parent and that of quark do not correspond: %d %d --> correcting\n",
309	this->GetParentFlavour(0), TMath::Abs(this->GetQuarkPDGCode(countP))));
310
311	pdg = this->GetQuarkPDGCode(countP);
312	Int_t line = this->GetQuarkIndex(countP);
313	this->ResetQuarkInfo();
314	while(TMath::Abs(pdg)!=this->GetParentFlavour(0)) {
315	pdg = stack->Particle(++line)->GetPdgCode();
316	}
317	while(line>=0){
318	mother = stack->Particle(line);
319	pdg = mother->GetPdgCode();
320	fQuarkIndex[countP] = line;
321	fQuarkPDGCode[countP++] = pdg;
322	line = mother->GetFirstMother();
323	}
324	}
325	}
326	return;
327	}
328
329	//-----------------------------------------------------------------------------
330	/AliMUONTrack AliMCMuonTrack::CovESDtoMuonTrack(AliESDMuonTrack &trkESD)
331	{
332	// method in $ALICE_ROOT/PWG3/muondep/AliAnalysisTaskESDMCLabelAddition.cxx
333	AliMUONTrack *trkMuon = new AliMUONTrack();
334	if(!trkESD.ClustersStored()) return trkMuon;
335
336	AliMUONTrackParam param;
337	AliMUONESDInterface::GetParamAtFirstCluster(trkESD, param);
338	AliMUONESDInterface::GetParamCov(trkESD, param);
339	AliMUONVClusterStore *cStore = AliMUONESDInterface::NewClusterStore();
340	AliMUONVCluster *cluster = cStore->CreateCluster(0,0,0);
341
342	AliESDMuonCluster esdCluster = (AliESDMuonCluster)trkESD.GetClusters().First();
343	while (esdCluster) {
344	AliMUONESDInterface::ESDToMUON(esdCluster, cluster);
345	param.SetZ(cluster->GetZ());
346	trkMuon->AddTrackParamAtCluster(param, *cluster, kTRUE);
347	esdCluster = (AliESDMuonCluster*)trkESD.GetClusters().After(esdCluster);
348	}
349
350	delete cluster; cluster = 0;
351	delete cStore; cStore = 0;
352	return trkMuon;
353	}*/
354
355	//-----------------------------------------------------------------------------
356	Int_t AliMCMuonTrack::SelectHFMuon()
357	{
358	Int_t flv = GetParentFlavour(0);
359	if (flv!=4 && flv!=5) return 2;
360
361	Bool_t isRes = kFALSE;
362	Int_t i=0, nparents=this->GetNParents();
363	while (i<nparents && !isRes) {
364	isRes = IsMotherAResonance(i++);
365	}
366
367	if (isRes) return 2;
368	if (flv==5) return 0;
369	else return 1;
370	}
371
372	//-----------------------------------------------------------------------------
373	Bool_t AliMCMuonTrack::IsDiquark(Int_t pdg)
374	{
375	// copy from $ALICE_ROOT/MUON/AliMUONTrackLight.cxx
376	pdg = TMath::Abs(pdg);
377	if(pdg>1000 && (pdg%100)<10) return kTRUE;
378	else return kFALSE;
379	}
380
381	//-----------------------------------------------------------------------------
382	void AliMCMuonTrack::ResetQuarkInfo()
383	{
384	for(Int_t pos=1; pos<4; pos++) {
385	fQuarkIndex[pos] = -1;
386	fQuarkPDGCode[pos] = 0;
387	}
388	return;
389	}
390
391	//-----------------------------------------------------------------------------
392	Int_t AliMCMuonTrack::GetParentFlavour(Int_t i) const
393	{
394	Int_t pdg = GetParentPDGCode(i);
395	pdg = TMath::Abs(pdg/100);
396	if(pdg>9) pdg /= 10;
397	return pdg;
398	}
399
400	//-----------------------------------------------------------------------------
401	Bool_t AliMCMuonTrack::IsMotherAResonance(Int_t i) const
402	{
403	// copy from $ALICE_ROOT/MUON/AliMUONTrackLight.cxx
404	Int_t pdg = GetParentPDGCode(i);
405	Int_t id=pdg%100000;
406	return (!((id-id%10)%110));
407	}