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 **************************************************************************/
18 /////////////////////////////////////////////////////////////
20 // class used to extract and store info of MC particle
22 // Author: X-M. Zhang, zhang@clermont.in2p3.fr
23 // zhangxm@iopp.ccnu.edu.cn
24 /////////////////////////////////////////////////////////////
26 #include <TParticle.h>
28 #include "AliMCEvent.h"
29 #include "AliAODMCParticle.h"
30 #include "AliESDMuonTrack.h"
31 #include "AliAODTrack.h"
32 #include "AliMuonInfoStoreRD.h"
33 #include "AliMuonInfoStoreMC.h"
34 #include "AliMCEvent.h"
36 ClassImp(AliMuonInfoStoreMC)
38 const TString AliMuonInfoStoreMC::fgkStdBranchName("MuonMC");
39 const Int_t AliMuonInfoStoreMC::fgkSourcesN = 6;
41 //-----------------------------------------------------------------------------
42 AliMuonInfoStoreMC::AliMuonInfoStoreMC() :
54 // default constructor
56 for (Int_t i=5; i--;) { fParentIndex[i] = -1; fParentPDGCode[i] = 0; }
57 for (Int_t i=4; i--;) { fQuarkIndex[i] = -1; fQuarkPDGCode[i] = 0; }
60 //-----------------------------------------------------------------------------
61 AliMuonInfoStoreMC::AliMuonInfoStoreMC(AliAODTrack *trkAOD, AliMCEvent *mcEvent, UInt_t selMask, Bool_t full) :
62 AliMuonInfoStoreRD(trkAOD,selMask),
73 // default constructor
75 for (Int_t i=5; i--;) { fParentIndex[i] = -1; fParentPDGCode[i] = 0; }
76 for (Int_t i=4; i--;) { fQuarkIndex[i] = -1; fQuarkPDGCode[i] = 0; }
78 this->SetMCInfoAOD(mcEvent, trkAOD->GetLabel());
81 //-----------------------------------------------------------------------------
82 AliMuonInfoStoreMC::AliMuonInfoStoreMC(AliESDMuonTrack *trkESD, AliMCEvent *mcEvent, UInt_t selMask, Bool_t full) :
83 AliMuonInfoStoreRD(trkESD,selMask),
94 // default constructor
96 for (Int_t i=5; i--;) { fParentIndex[i] = -1; fParentPDGCode[i] = 0; }
97 for (Int_t i=4; i--;) { fQuarkIndex[i] = -1; fQuarkPDGCode[i] = 0; }
99 this->SetMCInfoESD(mcEvent, trkESD->GetLabel());
102 //-----------------------------------------------------------------------------
103 AliMuonInfoStoreMC::AliMuonInfoStoreMC(const AliMuonInfoStoreMC &src) :
104 AliMuonInfoStoreRD(src),
105 fIsFull(src.fIsFull),
106 fLorentzP(src.fLorentzP),
107 fTrackIndex(src.fTrackIndex),
108 fTrackPDGCode(src.fTrackPDGCode),
109 fSource(src.fSource),
110 fNParents(src.fNParents),
111 fOscillation(src.fOscillation),
117 for (Int_t i=5; i--;) {
118 fParentIndex[i] = src.fParentIndex[i];
119 fParentPDGCode[i] = src.fParentPDGCode[i];
121 for (Int_t i=4; i--;) {
122 fQuarkIndex[i] = src.fQuarkIndex[i];
123 fQuarkPDGCode[i] = src.fQuarkPDGCode[i];
127 //-----------------------------------------------------------------------------
128 AliMuonInfoStoreMC& AliMuonInfoStoreMC::operator=(const AliMuonInfoStoreMC &src)
131 // assignment constructor
133 if(&src==this) return *this;
134 AliMuonInfoStoreRD::operator=(src);
136 fIsFull = src.fIsFull;
137 fLorentzP = src.fLorentzP;
138 fTrackIndex = src.fTrackIndex;
139 fTrackPDGCode = src.fTrackPDGCode;
140 fSource = src.fSource;
141 fNParents = src.fNParents;
142 fOscillation = src.fOscillation;
143 fWeight = src.fWeight;
144 for (Int_t i=5; i--;) {
145 fParentIndex[i] = src.fParentIndex[i];
146 fParentPDGCode[i] = src.fParentPDGCode[i];
148 for (Int_t i=4; i--;) {
149 fQuarkIndex[i] = src.fQuarkIndex[i];
150 fQuarkPDGCode[i] = src.fQuarkPDGCode[i];
156 //-----------------------------------------------------------------------------
157 AliMuonInfoStoreMC::~AliMuonInfoStoreMC()
164 //-----------------------------------------------------------------------------
165 void AliMuonInfoStoreMC::SetMCInfoAOD(AliMCEvent *mcEvent, Int_t label)
167 // fill track MC info with AOD base
169 if (fTrackIndex<0) { fSource=5; return; }
171 AliAODMCParticle *pMC = (AliAODMCParticle*)mcEvent->GetTrack(fTrackIndex);
172 fLorentzP.SetPxPyPzE(pMC->Px(), pMC->Py(), pMC->Pz(), pMC->E());
174 fTrackPDGCode = pMC->PdgCode();
175 if (TMath::Abs(fTrackPDGCode)!=13) { fSource=4; return; }
177 Int_t lineM = pMC->GetMother();
178 if (lineM<0) { fSource=2; return; }
180 Bool_t isPrimary = ((AliAODMCParticle*)mcEvent->GetTrack(lineM))->IsPrimary();
181 if (!isPrimary) { fSource=3; return; }
183 Int_t countP=-1, pdg=0;
184 Int_t parents[10], parLine[10];
185 AliAODMCParticle *mother = 0;
187 mother = (AliAODMCParticle*)mcEvent->GetTrack(lineM);
188 pdg = mother->GetPdgCode();
189 if(pdg==92 || pdg==21 || TMath::Abs(pdg)<10 || IsDiquark(pdg)) break;
190 parents[++countP] = pdg;
191 parLine[countP] = lineM;
192 lineM = mother->GetMother();
194 for(Int_t i=0; i<=countP; i++){
195 fParentIndex[i] = parLine[countP-i];
196 fParentPDGCode[i] = parents[countP-i];
198 fNParents = countP + 1;
200 if (fIsFull && lineM>=0) this->FillHistoryQuarksAOD(mcEvent, lineM);
202 fSource = this->SelectHFMuon();
206 //-----------------------------------------------------------------------------
207 void AliMuonInfoStoreMC::SetMCInfoESD(AliMCEvent *mcEvent, Int_t label)
209 // fill track MC info with ESD base
211 if (fTrackIndex<0) { fSource=5; return; }
213 TParticle *pMC = ((AliMCParticle*)mcEvent->GetTrack(fTrackIndex))->Particle();
214 fLorentzP.SetPxPyPzE(pMC->Px(), pMC->Py(), pMC->Pz(), pMC->Energy());
216 fTrackPDGCode = pMC->GetPdgCode();
217 if (TMath::Abs(fTrackPDGCode)!=13) { fSource=4; return; }
219 Int_t lineM = pMC->GetFirstMother();
220 if (lineM<0) { fSource=2; return; }
222 if (lineM>=mcEvent->Stack()->GetNprimary()) { fSource=3; return; }
224 Int_t countP=-1, pdg=0;
225 Int_t parents[10], parLine[10];
226 TParticle *mother = 0;
228 mother = ((AliMCParticle*)mcEvent->GetTrack(lineM))->Particle();
229 pdg = mother->GetPdgCode();
230 if(pdg==92 || pdg==21 || TMath::Abs(pdg)<10 || IsDiquark(pdg)) break;
231 parents[++countP] = pdg;
232 parLine[countP] = lineM;
233 lineM = mother->GetFirstMother();
235 for(Int_t i=0; i<=countP; i++){
236 fParentIndex[i] = parLine[countP-i];
237 fParentPDGCode[i] = parents[countP-i];
239 fNParents = countP + 1;
241 if (fIsFull && lineM>=0) this->FillHistoryQuarksESD(mcEvent, lineM);
243 fSource = this->SelectHFMuon();
247 //-----------------------------------------------------------------------------
248 void AliMuonInfoStoreMC::FillHistoryQuarksAOD(AliMCEvent* const mcEvent, Int_t lineM)
250 // method in $ALICE_ROOT/MUON/AliMUONTrackLight.cxx
253 Int_t countP=-1, pdg=0;
254 AliAODMCParticle *mother = 0;
256 mother = (AliAODMCParticle*)mcEvent->GetTrack(lineM);
257 pdg = mother->GetPdgCode();
258 fQuarkIndex[++countP] = lineM;
259 fQuarkPDGCode[countP] = pdg;
260 lineM = mother->GetMother();
263 // for PYTHIA checking
265 for(Int_t par=0; par<4; par++) {
266 if(TMath::Abs(this->QuarkPDGCode(par))<6) { countP=par; break; }
268 if(this->QuarkIndex(countP)>-1 && (this->ParentFlavour(0)==4 || this->ParentFlavour(0)==5)) {
269 if(this->ParentFlavour(0)!=TMath::Abs(this->QuarkPDGCode(countP))) {
270 AliWarning(Form("quark flavour of parent and that of quark do not correspond: %d %d --> correcting\n",
271 this->ParentFlavour(0), TMath::Abs(this->QuarkPDGCode(countP))));
273 pdg = this->QuarkPDGCode(countP);
274 Int_t line = this->QuarkIndex(countP);
275 this->ResetQuarkInfo();
276 while(TMath::Abs(pdg)!=this->ParentFlavour(0)) {
277 pdg = ((AliAODMCParticle*)mcEvent->GetTrack(++line))->GetPdgCode();
280 mother = (AliAODMCParticle*)mcEvent->GetTrack(line);
281 pdg = mother->GetPdgCode();
282 fQuarkIndex[countP] = line;
283 fQuarkPDGCode[countP++] = pdg;
284 line = mother->GetMother();
291 //-----------------------------------------------------------------------------
292 void AliMuonInfoStoreMC::FillHistoryQuarksESD(AliMCEvent* const mcEvent, Int_t lineM)
294 // method in $ALICE_ROOT/MUON/AliMUONTrackLight.cxx
297 Int_t countP=-1, pdg=0;
298 TParticle *mother = 0;
300 mother = ((AliMCParticle*)mcEvent->GetTrack(lineM))->Particle();
301 pdg = mother->GetPdgCode();
302 fQuarkIndex[++countP] = lineM;
303 fQuarkPDGCode[countP] = pdg;
304 lineM = mother->GetFirstMother();
307 // for PYTHIA checking
309 for(Int_t par=0; par<4; par++) {
310 if(TMath::Abs(this->QuarkPDGCode(par))<6) { countP=par; break; }
312 if(this->QuarkIndex(countP)>-1 && (this->ParentFlavour(0)==4 || this->ParentFlavour(0)==5)) {
313 if(this->ParentFlavour(0)!=TMath::Abs(this->QuarkPDGCode(countP))) {
314 AliWarning(Form("quark flavour of parent and that of quark do not correspond: %d %d --> correcting\n",
315 this->ParentFlavour(0), TMath::Abs(this->QuarkPDGCode(countP))));
317 pdg = this->QuarkPDGCode(countP);
318 Int_t line = this->QuarkIndex(countP);
319 this->ResetQuarkInfo();
320 while(TMath::Abs(pdg)!=this->ParentFlavour(0)) {
321 pdg = ((AliMCParticle*)mcEvent->GetTrack(++lineM))->Particle()->GetPdgCode();
324 mother = ((AliMCParticle*)mcEvent->GetTrack(++lineM))->Particle();
325 pdg = mother->GetPdgCode();
326 fQuarkIndex[countP] = line;
327 fQuarkPDGCode[countP++] = pdg;
328 line = mother->GetFirstMother();
335 //-----------------------------------------------------------------------------
336 Int_t AliMuonInfoStoreMC::SelectHFMuon()
338 // set info of muon from HF
340 Int_t flv = ParentFlavour(0);
341 if (flv!=4 && flv!=5) return 2;
343 Bool_t isRes = kFALSE;
344 Int_t i=0, nparents=this->ParentsN();
345 while (i<nparents && !isRes) {
346 isRes = IsMotherAResonance(i++);
350 if (flv==5) return 0;
354 //-----------------------------------------------------------------------------
355 Bool_t AliMuonInfoStoreMC::IsDiquark(Int_t pdg)
357 // copy from $ALICE_ROOT/MUON/AliMUONTrackLight.cxx
358 pdg = TMath::Abs(pdg);
359 if(pdg>1000 && (pdg%100)<10) return kTRUE;
363 //-----------------------------------------------------------------------------
364 void AliMuonInfoStoreMC::ResetQuarkInfo()
366 // copy from $ALICE_ROOT/MUON/AliMUONTrackLight.cxx
367 for(Int_t pos=1; pos<4; pos++) {
368 fQuarkIndex[pos] = -1;
369 fQuarkPDGCode[pos] = 0;
374 //-----------------------------------------------------------------------------
375 Int_t AliMuonInfoStoreMC::ParentFlavour(Int_t i) const
377 // copy from $ALICE_ROOT/MUON/AliMUONTrackLight.cxx
378 Int_t pdg = ParentPDGCode(i);
379 pdg = TMath::Abs(pdg/100);
384 //-----------------------------------------------------------------------------
385 Bool_t AliMuonInfoStoreMC::IsMotherAResonance(Int_t i) const
387 // copy from $ALICE_ROOT/MUON/AliMUONTrackLight.cxx
388 Int_t pdg = ParentPDGCode(i);
390 return (!((id-id%10)%110));