Keep e+e- from pair production of primary gammas.
[u/mrichter/AliRoot.git] / PWG3 / muon / AliMCMuonTrack.cxx
CommitLineData
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
16ClassImp(AliMCMuonTrack)
17
18const Double_t AliMCMuonTrack::fgkSigmaCut = 10.;
19
20//-----------------------------------------------------------------------------
21AliMCMuonTrack::AliMCMuonTrack() :
22AliAODMuonTrack(),
23fIsFull(kFALSE),
24fPGen(),
25fTrackIndex(-1),
26fTrackPDGCode(0),
27fSource(-1),
28fNParents(0),
29fOscillation(kFALSE),
30fWeight(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//-----------------------------------------------------------------------------
46AliMCMuonTrack::AliMCMuonTrack(AliAODTrack *trkAOD, TClonesArray *mcClArr, Bool_t full) :
47AliAODMuonTrack(trkAOD),
48fIsFull(full),
49fPGen(),
50fTrackIndex(-1),
51fTrackPDGCode(0),
52fSource(-1),
53fNParents(0),
54fOscillation(kFALSE),
55fWeight(0.)
56{
57 //
58 // default constructor
59 //
60 AliAODMCParticle *pMC = this->FindTrackRef(trkAOD, mcClArr);
61 if (pMC) this->SetMCInfo(pMC, mcClArr);
62}
63
64//-----------------------------------------------------------------------------
65AliMCMuonTrack::AliMCMuonTrack(AliESDMuonTrack *trkESD, AliESDEvent *esd, AliMCEventHandler *mcH, Bool_t full) :
66AliAODMuonTrack(trkESD),
67fIsFull(full),
68fPGen(),
69fTrackIndex(-1),
70fTrackPDGCode(0),
71fSource(-1),
72fNParents(0),
73fOscillation(kFALSE),
74fWeight(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//-----------------------------------------------------------------------------
97AliMCMuonTrack::~AliMCMuonTrack()
98{
99 //
100 // destructor
101 //
102}
103
104//-----------------------------------------------------------------------------
105AliAODMCParticle* 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//-----------------------------------------------------------------------------
115TParticle* 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//-----------------------------------------------------------------------------
141void 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//-----------------------------------------------------------------------------
169void 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//-----------------------------------------------------------------------------
196void 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//-----------------------------------------------------------------------------
220void 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//-----------------------------------------------------------------------------
244void 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//-----------------------------------------------------------------------------
287void 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//-----------------------------------------------------------------------------
356Int_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//-----------------------------------------------------------------------------
373Bool_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//-----------------------------------------------------------------------------
382void 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//-----------------------------------------------------------------------------
392Int_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//-----------------------------------------------------------------------------
401Bool_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}