1 /**************************************************************************
2 * Copyright(c) 1998-2006, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
16 /////////////////////////////////////////////////////////////
18 // class used to extract and store info of MC particle
20 // Author: X-M. Zhang, zhang@clermont.in2p3.fr
21 // zhangxm@iopp.ccnu.edu.cn
22 /////////////////////////////////////////////////////////////
24 #include <TParticle.h>
26 #include "AliMCEvent.h"
27 #include "AliAODMCParticle.h"
28 #include "AliESDMuonTrack.h"
29 #include "AliAODTrack.h"
30 #include "AliMuonInfoStoreRD.h"
31 #include "AliMuonInfoStoreMC.h"
32 #include "AliMCEvent.h"
34 ClassImp(AliMuonInfoStoreMC)
36 const TString AliMuonInfoStoreMC::fgkStdBranchName("MuonMC");
37 const Int_t AliMuonInfoStoreMC::fgkSourcesN = 6;
39 //-----------------------------------------------------------------------------
40 AliMuonInfoStoreMC::AliMuonInfoStoreMC() :
52 // default constructor
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; }
58 //-----------------------------------------------------------------------------
59 AliMuonInfoStoreMC::AliMuonInfoStoreMC(AliAODTrack *trkAOD, AliMCEvent *mcEvent, Bool_t full) :
60 AliMuonInfoStoreRD(trkAOD),
71 // default constructor
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; }
76 this->SetMCInfoAOD(mcEvent, trkAOD->GetLabel());
79 //-----------------------------------------------------------------------------
80 AliMuonInfoStoreMC::AliMuonInfoStoreMC(AliESDMuonTrack *trkESD, AliMCEvent *mcEvent, Bool_t full) :
81 AliMuonInfoStoreRD(trkESD),
92 // default constructor
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; }
97 this->SetMCInfoESD(mcEvent, trkESD->GetLabel());
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),
115 for (Int_t i=5; i--;) {
116 fParentIndex[i] = src.fParentIndex[i];
117 fParentPDGCode[i] = src.fParentPDGCode[i];
119 for (Int_t i=4; i--;) {
120 fQuarkIndex[i] = src.fQuarkIndex[i];
121 fQuarkPDGCode[i] = src.fQuarkPDGCode[i];
125 //-----------------------------------------------------------------------------
126 AliMuonInfoStoreMC& AliMuonInfoStoreMC::operator=(const AliMuonInfoStoreMC &src)
129 // assignment constructor
131 if(&src==this) return *this;
132 AliMuonInfoStoreRD::operator=(src);
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];
146 for (Int_t i=4; i--;) {
147 fQuarkIndex[i] = src.fQuarkIndex[i];
148 fQuarkPDGCode[i] = src.fQuarkPDGCode[i];
154 //-----------------------------------------------------------------------------
155 AliMuonInfoStoreMC::~AliMuonInfoStoreMC()
162 //-----------------------------------------------------------------------------
163 void AliMuonInfoStoreMC::SetMCInfoAOD(AliMCEvent *mcEvent, Int_t label)
165 // fill track MC info with AOD base
167 if (fTrackIndex<0) { fSource=5; return; }
169 AliAODMCParticle *pMC = (AliAODMCParticle*)mcEvent->GetTrack(fTrackIndex);
170 fLorentzP.SetPxPyPzE(pMC->Px(), pMC->Py(), pMC->Pz(), pMC->E());
172 fTrackPDGCode = pMC->PdgCode();
173 if (TMath::Abs(fTrackPDGCode)!=13) { fSource=4; return; }
175 Int_t lineM = pMC->GetMother();
176 if (lineM<0) { fSource=2; return; }
178 Bool_t isPrimary = ((AliAODMCParticle*)mcEvent->GetTrack(lineM))->IsPrimary();
179 if (!isPrimary) { fSource=3; return; }
181 Int_t countP=-1, pdg=0;
182 Int_t parents[10], parLine[10];
183 AliAODMCParticle *mother = 0;
185 mother = (AliAODMCParticle*)mcEvent->GetTrack(lineM);
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();
192 for(Int_t i=0; i<=countP; i++){
193 fParentIndex[i] = parLine[countP-i];
194 fParentPDGCode[i] = parents[countP-i];
196 fNParents = countP + 1;
198 if (fIsFull && lineM>=0) this->FillHistoryQuarksAOD(mcEvent, lineM);
200 fSource = this->SelectHFMuon();
204 //-----------------------------------------------------------------------------
205 void AliMuonInfoStoreMC::SetMCInfoESD(AliMCEvent *mcEvent, Int_t label)
207 // fill track MC info with ESD base
209 if (fTrackIndex<0) { fSource=5; return; }
211 TParticle *pMC = ((AliMCParticle*)mcEvent->GetTrack(fTrackIndex))->Particle();
212 fLorentzP.SetPxPyPzE(pMC->Px(), pMC->Py(), pMC->Pz(), pMC->Energy());
214 fTrackPDGCode = pMC->GetPdgCode();
215 if (TMath::Abs(fTrackPDGCode)!=13) { fSource=4; return; }
217 Int_t lineM = pMC->GetFirstMother();
218 if (lineM<0) { fSource=2; return; }
220 if (lineM>=mcEvent->Stack()->GetNprimary()) { fSource=3; return; }
222 Int_t countP=-1, pdg=0;
223 Int_t parents[10], parLine[10];
224 TParticle *mother = 0;
226 mother = ((AliMCParticle*)mcEvent->GetTrack(lineM))->Particle();
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();
233 for(Int_t i=0; i<=countP; i++){
234 fParentIndex[i] = parLine[countP-i];
235 fParentPDGCode[i] = parents[countP-i];
237 fNParents = countP + 1;
239 if (fIsFull && lineM>=0) this->FillHistoryQuarksESD(mcEvent, lineM);
241 fSource = this->SelectHFMuon();
245 //-----------------------------------------------------------------------------
246 void AliMuonInfoStoreMC::FillHistoryQuarksAOD(AliMCEvent* const mcEvent, Int_t lineM)
248 // method in $ALICE_ROOT/MUON/AliMUONTrackLight.cxx
251 Int_t countP=-1, pdg=0;
252 AliAODMCParticle *mother = 0;
254 mother = (AliAODMCParticle*)mcEvent->GetTrack(lineM);
255 pdg = mother->GetPdgCode();
256 fQuarkIndex[++countP] = lineM;
257 fQuarkPDGCode[countP] = pdg;
258 lineM = mother->GetMother();
261 // for PYTHIA checking
263 for(Int_t par=0; par<4; par++) {
264 if(TMath::Abs(this->QuarkPDGCode(par))<6) { countP=par; break; }
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))));
271 pdg = this->QuarkPDGCode(countP);
272 Int_t line = this->QuarkIndex(countP);
273 this->ResetQuarkInfo();
274 while(TMath::Abs(pdg)!=this->ParentFlavour(0)) {
275 pdg = ((AliAODMCParticle*)mcEvent->GetTrack(++line))->GetPdgCode();
278 mother = (AliAODMCParticle*)mcEvent->GetTrack(line);
279 pdg = mother->GetPdgCode();
280 fQuarkIndex[countP] = line;
281 fQuarkPDGCode[countP++] = pdg;
282 line = mother->GetMother();
289 //-----------------------------------------------------------------------------
290 void AliMuonInfoStoreMC::FillHistoryQuarksESD(AliMCEvent* const mcEvent, Int_t lineM)
292 // method in $ALICE_ROOT/MUON/AliMUONTrackLight.cxx
295 Int_t countP=-1, pdg=0;
296 TParticle *mother = 0;
298 mother = ((AliMCParticle*)mcEvent->GetTrack(lineM))->Particle();
299 pdg = mother->GetPdgCode();
300 fQuarkIndex[++countP] = lineM;
301 fQuarkPDGCode[countP] = pdg;
302 lineM = mother->GetFirstMother();
305 // for PYTHIA checking
307 for(Int_t par=0; par<4; par++) {
308 if(TMath::Abs(this->QuarkPDGCode(par))<6) { countP=par; break; }
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))));
315 pdg = this->QuarkPDGCode(countP);
316 Int_t line = this->QuarkIndex(countP);
317 this->ResetQuarkInfo();
318 while(TMath::Abs(pdg)!=this->ParentFlavour(0)) {
319 pdg = ((AliMCParticle*)mcEvent->GetTrack(++lineM))->Particle()->GetPdgCode();
322 mother = ((AliMCParticle*)mcEvent->GetTrack(++lineM))->Particle();
323 pdg = mother->GetPdgCode();
324 fQuarkIndex[countP] = line;
325 fQuarkPDGCode[countP++] = pdg;
326 line = mother->GetFirstMother();
333 //-----------------------------------------------------------------------------
334 Int_t AliMuonInfoStoreMC::SelectHFMuon()
336 // set info of muon from HF
338 Int_t flv = ParentFlavour(0);
339 if (flv!=4 && flv!=5) return 2;
341 Bool_t isRes = kFALSE;
342 Int_t i=0, nparents=this->ParentsN();
343 while (i<nparents && !isRes) {
344 isRes = IsMotherAResonance(i++);
348 if (flv==5) return 0;
352 //-----------------------------------------------------------------------------
353 Bool_t AliMuonInfoStoreMC::IsDiquark(Int_t pdg)
355 // copy from $ALICE_ROOT/MUON/AliMUONTrackLight.cxx
356 pdg = TMath::Abs(pdg);
357 if(pdg>1000 && (pdg%100)<10) return kTRUE;
361 //-----------------------------------------------------------------------------
362 void AliMuonInfoStoreMC::ResetQuarkInfo()
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;
372 //-----------------------------------------------------------------------------
373 Int_t AliMuonInfoStoreMC::ParentFlavour(Int_t i) const
375 // copy from $ALICE_ROOT/MUON/AliMUONTrackLight.cxx
376 Int_t pdg = ParentPDGCode(i);
377 pdg = TMath::Abs(pdg/100);
382 //-----------------------------------------------------------------------------
383 Bool_t AliMuonInfoStoreMC::IsMotherAResonance(Int_t i) const
385 // copy from $ALICE_ROOT/MUON/AliMUONTrackLight.cxx
386 Int_t pdg = ParentPDGCode(i);
388 return (!((id-id%10)%110));