6f646b5373b108eb4015fd88e245c4316cb8d87b
[u/mrichter/AliRoot.git] / PWG3 / muon / AliMCMuonTrack.cxx
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 }