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>
25 #include <TClonesArray.h>
27 #include "AliMCEvent.h"
28 #include "AliMCEventHandler.h"
30 #include "AliAODMCParticle.h"
31 #include "AliESDMuonTrack.h"
32 #include "AliAODTrack.h"
33 #include "AliMuonInfoStoreRD.h"
34 #include "AliMuonInfoStoreMC.h"
38 ClassImp(AliMuonInfoStoreMC)
40 const TString AliMuonInfoStoreMC::fgkStdBranchName("MuonMC");
41 const Int_t AliMuonInfoStoreMC::fgkNSources = 7;
43 //-----------------------------------------------------------------------------
44 AliMuonInfoStoreMC::AliMuonInfoStoreMC() :
56 // default constructor
58 for (Int_t i=5; i--;) { fParentIndex[i] = -1; fParentPDGCode[i] = 0; }
59 for (Int_t i=4; i--;) { fQuarkIndex[i] = -1; fQuarkPDGCode[i] = 0; }
62 //-----------------------------------------------------------------------------
63 AliMuonInfoStoreMC::AliMuonInfoStoreMC(AliAODTrack *trkAOD, TClonesArray *mcClArr, Bool_t full) :
64 AliMuonInfoStoreRD(trkAOD),
75 // default constructor
77 for (Int_t i=5; i--;) { fParentIndex[i] = -1; fParentPDGCode[i] = 0; }
78 for (Int_t i=4; i--;) { fQuarkIndex[i] = -1; fQuarkPDGCode[i] = 0; }
80 AliAODMCParticle *pMC = this->FindTrackRef(trkAOD, mcClArr);
81 if (pMC) this->SetMCInfo(pMC, mcClArr);
84 //-----------------------------------------------------------------------------
85 AliMuonInfoStoreMC::AliMuonInfoStoreMC(AliESDMuonTrack *trkESD, AliMCEventHandler *mcH, Bool_t full) :
86 AliMuonInfoStoreRD(trkESD),
97 // default constructor
99 for (Int_t i=5; i--;) { fParentIndex[i] = -1; fParentPDGCode[i] = 0; }
100 for (Int_t i=4; i--;) { fQuarkIndex[i] = -1; fQuarkPDGCode[i] = 0; }
102 TParticle *pMC = this->FindTrackRef(trkESD, mcH);
103 if (pMC) this->SetMCInfo(pMC, mcH);
106 //-----------------------------------------------------------------------------
107 /*AliMuonInfoStoreMC::AliMuonInfoStoreMC(AliESDMuonTrack *trkESD, AliESDEvent *esd, AliMCEventHandler *mcH, Bool_t full) :
108 AliMuonInfoStoreRD(trkESD),
115 fOscillation(kFALSE),
118 #include "AliMUONRecoCheck.h"
119 #include "AliMUONVTrackStore.h"
120 #include "AliMUONTrack.h"
121 #include "AliMUONESDInterface.h"
123 // default constructor
125 for (Int_t i=5; i--;) { fParentIndex[i] = -1; fParentPDGCode[i] = 0; }
126 for (Int_t i=4; i--;) { fQuarkIndex[i] = -1; fQuarkPDGCode[i] = 0; }
128 TParticle *pMC = this->FindTrackRef(trkESD, esd, mcH);
129 if (pMC) this->SetMCInfo(pMC, mcH);
132 //-----------------------------------------------------------------------------
133 AliMuonInfoStoreMC::AliMuonInfoStoreMC(const AliMuonInfoStoreMC &src) :
134 AliMuonInfoStoreRD(src),
135 fIsFull(src.fIsFull),
136 fLorentzP(src.fLorentzP),
137 fTrackIndex(src.fTrackIndex),
138 fTrackPDGCode(src.fTrackPDGCode),
139 fSource(src.fSource),
140 fNParents(src.fNParents),
141 fOscillation(src.fOscillation),
147 for (Int_t i=5; i--;) {
148 fParentIndex[i] = src.fParentIndex[i];
149 fParentPDGCode[i] = src.fParentPDGCode[i];
151 for (Int_t i=4; i--;) {
152 fQuarkIndex[i] = src.fQuarkIndex[i];
153 fQuarkPDGCode[i] = src.fQuarkPDGCode[i];
157 //-----------------------------------------------------------------------------
158 AliMuonInfoStoreMC& AliMuonInfoStoreMC::operator=(const AliMuonInfoStoreMC &src)
161 // assignment constructor
163 if(&src==this) return *this;
164 AliMuonInfoStoreRD::operator=(src);
166 fIsFull = src.fIsFull;
167 fLorentzP = src.fLorentzP;
168 fTrackIndex = src.fTrackIndex;
169 fTrackPDGCode = src.fTrackPDGCode;
170 fSource = src.fSource;
171 fNParents = src.fNParents;
172 fOscillation = src.fOscillation;
173 fWeight = src.fWeight;
174 for (Int_t i=5; i--;) {
175 fParentIndex[i] = src.fParentIndex[i];
176 fParentPDGCode[i] = src.fParentPDGCode[i];
178 for (Int_t i=4; i--;) {
179 fQuarkIndex[i] = src.fQuarkIndex[i];
180 fQuarkPDGCode[i] = src.fQuarkPDGCode[i];
186 //-----------------------------------------------------------------------------
187 AliMuonInfoStoreMC::~AliMuonInfoStoreMC()
194 //-----------------------------------------------------------------------------
195 AliAODMCParticle* AliMuonInfoStoreMC::FindTrackRef(AliAODTrack* const trkAOD, TClonesArray* const mcClArr)
197 // find MC track ref with AOD base
199 AliAODMCParticle *pMC = 0;
200 fTrackIndex = trkAOD->GetLabel();
202 pMC = (AliAODMCParticle*)mcClArr->At(fTrackIndex);
206 //-----------------------------------------------------------------------------
207 TParticle* AliMuonInfoStoreMC::FindTrackRef(AliESDMuonTrack* const trkESD, AliMCEventHandler* const mcH)
209 // find MC track ref with ESD base
211 TParticle *pMCRef = 0;
212 fTrackIndex = trkESD->GetLabel();
213 if (fTrackIndex>=0) pMCRef = mcH->MCEvent()->Stack()->Particle(fTrackIndex);
217 //-----------------------------------------------------------------------------
218 /*TParticle* AliMuonInfoStoreMC::FindTrackRef(AliESDMuonTrack* const trkESD, AliESDEvent* const esd, AliMCEventHandler* const mcH)
220 // find MC track ref with ESD trackRef base
222 TParticle *pMCRef = 0;
223 AliMUONRecoCheck rc(esd,mcH);
224 AliMUONVTrackStore *trkRefArr = rc.TrackRefs(-1);
226 AliMUONTrack trkMuon;
227 AliMUONESDInterface::ESDToMUON(*trkESD, trkMuon, kFALSE);
229 Int_t nMatchClusters = 0;
230 AliMUONTrack *trkRef = rc.FindCompatibleTrack(trkMuon, *trkRefArr, nMatchClusters, kFALSE, 10.);
231 if (trkRef) fTrackIndex = trkRef->GetUniqueID();
232 if (fTrackIndex>=0) pMCRef = mcH->MCEvent()->Stack()->Particle(fTrackIndex);
236 //-----------------------------------------------------------------------------
237 void AliMuonInfoStoreMC::SetMCInfo(AliAODMCParticle *pMC, TClonesArray *mcClArr)
239 // fill track MC info with AOD base
241 fLorentzP.SetPxPyPzE(pMC->Px(), pMC->Py(), pMC->Pz(), pMC->E());
242 fTrackPDGCode = pMC->GetPdgCode();
243 if (TMath::Abs(fTrackPDGCode)!=13) {
248 Int_t lineM = pMC->GetMother();
254 Bool_t isPrimary = ((AliAODMCParticle*)mcClArr->At(lineM))->IsPrimary();
260 this->FillHistoryParents(lineM, mcClArr);
261 fSource = this->SelectHFMuon();
265 //-----------------------------------------------------------------------------
266 void AliMuonInfoStoreMC::SetMCInfo(TParticle *pMC, AliMCEventHandler* const mcH)
268 // fill track MC info with ESD base
270 fLorentzP.SetPxPyPzE(pMC->Px(), pMC->Py(), pMC->Pz(), pMC->Energy());
271 fTrackPDGCode = pMC->GetPdgCode();
272 if (TMath::Abs(fTrackPDGCode)!=13) {
277 Int_t lineM = pMC->GetFirstMother();
283 AliStack *stack = mcH->MCEvent()->Stack();
284 if (lineM>=stack->GetNprimary()) {
289 this->FillHistoryParents(lineM, stack);
290 fSource = this->SelectHFMuon();
294 //-----------------------------------------------------------------------------
295 void AliMuonInfoStoreMC::FillHistoryParents(Int_t lineM, TClonesArray *mcClArr)
297 // find track hadron parents with AOD base
299 Int_t countP=-1, pdg=0;
300 Int_t parents[10], parLine[10];
301 AliAODMCParticle *mother = 0;
303 mother = (AliAODMCParticle*)mcClArr->At(lineM);
304 pdg = mother->GetPdgCode();
305 if(pdg==92 || pdg==21 || TMath::Abs(pdg)<10 || IsDiquark(pdg)) break;
306 parents[++countP] = pdg;
307 parLine[countP] = lineM;
308 lineM = mother->GetMother();
310 for(Int_t i=0; i<=countP; i++){
311 fParentIndex[i] = parLine[countP-i];
312 fParentPDGCode[i] = parents[countP-i];
314 fNParents = countP + 1;
316 if (fIsFull && lineM>=0) this->FillHistoryQuarks(lineM, mcClArr);
320 //-----------------------------------------------------------------------------
321 void AliMuonInfoStoreMC::FillHistoryParents(Int_t lineM, AliStack *stack)
323 // find track hadron parents with ESD base
325 Int_t countP=-1, pdg=0;
326 Int_t parents[10], parLine[10];
327 TParticle *mother = 0;
329 mother = stack->Particle(lineM);
330 pdg = mother->GetPdgCode();
331 if(pdg==92 || pdg==21 || TMath::Abs(pdg)<10 || IsDiquark(pdg)) break;
332 parents[++countP] = pdg;
333 parLine[countP] = lineM;
334 lineM = mother->GetFirstMother();
336 for(Int_t i=0; i<=countP; i++){
337 fParentIndex[i] = parLine[countP-i];
338 fParentPDGCode[i] = parents[countP-i];
340 fNParents = countP + 1;
342 if (fIsFull && lineM>=0) this->FillHistoryQuarks(lineM, stack);
346 //-----------------------------------------------------------------------------
347 void AliMuonInfoStoreMC::FillHistoryQuarks(Int_t lineM, TClonesArray* const mcClArr)
349 // method in $ALICE_ROOT/MUON/AliMUONTrackLight.cxx
352 Int_t countP=-1, pdg=0;
353 AliAODMCParticle *mother = 0;
355 mother = (AliAODMCParticle*)mcClArr->At(lineM);
356 pdg = mother->GetPdgCode();
357 fQuarkIndex[++countP] = lineM;
358 fQuarkPDGCode[countP] = pdg;
359 lineM = mother->GetMother();
362 // for PYTHIA checking
364 for(Int_t par=0; par<4; par++) {
365 if(TMath::Abs(this->QuarkPDGCode(par))<6) { countP=par; break; }
367 if(this->QuarkIndex(countP)>-1 && (this->ParentFlavour(0)==4 || this->ParentFlavour(0)==5)) {
368 if(this->ParentFlavour(0)!=TMath::Abs(this->QuarkPDGCode(countP))) {
369 AliWarning(Form("quark flavour of parent and that of quark do not correspond: %d %d --> correcting\n",
370 this->ParentFlavour(0), TMath::Abs(this->QuarkPDGCode(countP))));
372 pdg = this->QuarkPDGCode(countP);
373 Int_t line = this->QuarkIndex(countP);
374 this->ResetQuarkInfo();
375 while(TMath::Abs(pdg)!=this->ParentFlavour(0)) {
376 pdg = ((AliAODMCParticle*)mcClArr->At(++line))->GetPdgCode();
379 mother = (AliAODMCParticle*)mcClArr->At(line);
380 pdg = mother->GetPdgCode();
381 fQuarkIndex[countP] = line;
382 fQuarkPDGCode[countP++] = pdg;
383 line = mother->GetMother();
390 //-----------------------------------------------------------------------------
391 void AliMuonInfoStoreMC::FillHistoryQuarks(Int_t lineM, AliStack* const stack)
393 // method in $ALICE_ROOT/MUON/AliMUONTrackLight.cxx
396 Int_t countP=-1, pdg=0;
397 TParticle *mother = 0;
399 mother = stack->Particle(lineM);
400 pdg = mother->GetPdgCode();
401 fQuarkIndex[++countP] = lineM;
402 fQuarkPDGCode[countP] = pdg;
403 lineM = mother->GetFirstMother();
406 // for PYTHIA checking
408 for(Int_t par=0; par<4; par++) {
409 if(TMath::Abs(this->QuarkPDGCode(par))<6) { countP=par; break; }
411 if(this->QuarkIndex(countP)>-1 && (this->ParentFlavour(0)==4 || this->ParentFlavour(0)==5)) {
412 if(this->ParentFlavour(0)!=TMath::Abs(this->QuarkPDGCode(countP))) {
413 AliWarning(Form("quark flavour of parent and that of quark do not correspond: %d %d --> correcting\n",
414 this->ParentFlavour(0), TMath::Abs(this->QuarkPDGCode(countP))));
416 pdg = this->QuarkPDGCode(countP);
417 Int_t line = this->QuarkIndex(countP);
418 this->ResetQuarkInfo();
419 while(TMath::Abs(pdg)!=this->ParentFlavour(0)) {
420 pdg = stack->Particle(++line)->GetPdgCode();
423 mother = stack->Particle(line);
424 pdg = mother->GetPdgCode();
425 fQuarkIndex[countP] = line;
426 fQuarkPDGCode[countP++] = pdg;
427 line = mother->GetFirstMother();
434 //-----------------------------------------------------------------------------
435 Int_t AliMuonInfoStoreMC::SelectHFMuon()
437 // set info of muon from HF
439 Int_t flv = ParentFlavour(0);
440 if (flv!=4 && flv!=5) return 2;
442 Bool_t isRes = kFALSE;
443 Int_t i=0, nparents=this->NParents();
444 while (i<nparents && !isRes) {
445 isRes = IsMotherAResonance(i++);
449 if (flv==5) return 0;
453 //-----------------------------------------------------------------------------
454 Bool_t AliMuonInfoStoreMC::IsDiquark(Int_t pdg)
456 // copy from $ALICE_ROOT/MUON/AliMUONTrackLight.cxx
457 pdg = TMath::Abs(pdg);
458 if(pdg>1000 && (pdg%100)<10) return kTRUE;
462 //-----------------------------------------------------------------------------
463 void AliMuonInfoStoreMC::ResetQuarkInfo()
465 // copy from $ALICE_ROOT/MUON/AliMUONTrackLight.cxx
466 for(Int_t pos=1; pos<4; pos++) {
467 fQuarkIndex[pos] = -1;
468 fQuarkPDGCode[pos] = 0;
473 //-----------------------------------------------------------------------------
474 Int_t AliMuonInfoStoreMC::ParentFlavour(Int_t i) const
476 // copy from $ALICE_ROOT/MUON/AliMUONTrackLight.cxx
477 Int_t pdg = ParentPDGCode(i);
478 pdg = TMath::Abs(pdg/100);
483 //-----------------------------------------------------------------------------
484 Bool_t AliMuonInfoStoreMC::IsMotherAResonance(Int_t i) const
486 // copy from $ALICE_ROOT/MUON/AliMUONTrackLight.cxx
487 Int_t pdg = ParentPDGCode(i);
489 return (!((id-id%10)%110));