AliRoot Core  a565103 (a565103)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AliMUONTrackLight.cxx
Go to the documentation of this file.
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  * SigmaEffect_thetadegrees *
4  * Author: The ALICE Off-line Project. *
5  * Contributors are mentioned in the code where appropriate. *
6  * *
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 purpeateose. It is *
13  * provided "as is" without express or implied warranty. *
14  **************************************************************************/
15 
16 /* $Id$ */
17 
18 //-----------------------------------------------------------------------------
19 // Compact information for the muon generated tracks in the MUON arm
20 // useful at the last stage of the analysis chain
21 // provides a link between the reconstructed track and the generated particle
22 // stores kinematical information at gen. and rec. level and
23 // the decay history of the muon, allowing the identification of the
24 // mother process
25 //
26 // To be used together with AliMUONPairLight
27 //
28 // This class was prepared by INFN Cagliari, July 2006
29 // (authors: H.Woehri, A.de Falco)
30 //-----------------------------------------------------------------------------
31 
32 // 13 Nov 2007:
33 // Added a temporary fix to FindRefTrack to be able to handle reconstructed tracks
34 // generated from ESD muon track information. The problem is that the ESD data at
35 // the moment only contains the first hit on chamber 1. Hopefully in the near future
36 // this will be fixed and all hit information will be available.
37 // - Artur Szostak <artursz@iafrica.com>
38 
39 #include "AliMUONTrackLight.h"
40 #include "AliMUONTrack.h"
41 #include "AliMUONConstants.h"
42 #include "AliMUONVTrackStore.h"
43 #include "AliMUONTrackParam.h"
44 
45 #include "AliESDMuonTrack.h"
46 #include "AliStack.h"
47 #include "AliLog.h"
48 
49 #include "TDatabasePDG.h"
50 #include "TParticle.h"
51 #include "TString.h"
52 
53 #include <cstdio>
54 #include <Riostream.h>
55 
57 
58 //===================================================================
59 
61  : TObject(),
62  fPrec(),
63  fIsTriggered(kFALSE),
64  fCharge(-999),
65  fChi2(-1),
66  fCentr(-1),
67  fPgen(),
68  fTrackPythiaLine(-999),
69  fTrackPDGCode(-999),
70  fOscillation(kFALSE),
71  fNParents(0),
72  fWeight(1)
73 {
75  fPgen.SetPxPyPzE(0.,0.,0.,0.);
76  fPrec.SetPxPyPzE(0.,0.,0.,0.);
77  for (Int_t i=0; i<3; i++) fXYZ[i]=-999;
78  for (Int_t npar = 0; npar < fgkNParentsMax; npar++){
79  fParentPDGCode[npar] = -1;
80  fParentPythiaLine[npar] = -1;
81  }
82  for (Int_t i = 0; i < 4; i++){
83  fQuarkPDGCode[i] = -1;
84  fQuarkPythiaLine[i] = -1;
85  }
86 }
87 
88 //============================================
90  : TObject(muonCopy),
91  fPrec(muonCopy.fPrec),
92  fIsTriggered(muonCopy.fIsTriggered),
93  fCharge(muonCopy.fCharge),
94  fChi2(muonCopy.fChi2),
95  fCentr(muonCopy.fCentr),
96  fPgen(muonCopy.fPgen),
97  fTrackPythiaLine(muonCopy.fTrackPythiaLine),
98  fTrackPDGCode(muonCopy.fTrackPDGCode),
99  fOscillation(muonCopy.fOscillation),
100  fNParents(muonCopy.fNParents),
101  fWeight(muonCopy.fWeight)
102 {
104  for (Int_t i=0; i<3; i++) fXYZ[i]=muonCopy.fXYZ[i];
105  for (Int_t npar = 0; npar < fgkNParentsMax; npar++){
106  fParentPDGCode[npar] = muonCopy.fParentPDGCode[npar];
107  fParentPythiaLine[npar] = muonCopy.fParentPythiaLine[npar];
108  }
109  for (Int_t i = 0; i < 4; i++){
110  fQuarkPDGCode[i] = muonCopy.fQuarkPDGCode[i];
111  fQuarkPythiaLine[i] = muonCopy.fQuarkPythiaLine[i];
112  }
113 }
114 
115 //============================================
116 AliMUONTrackLight::AliMUONTrackLight(AliESDMuonTrack* muonTrack)
117  : TObject(),
118  fPrec(),
119  fIsTriggered(kFALSE),
120  fCharge(-999),
121  fChi2(-1),
122  fCentr(-1),
123  fPgen(),
124  fTrackPythiaLine(-999),
125  fTrackPDGCode(-999),
126  fOscillation(kFALSE),
127  fNParents(0),
128  fWeight(1)
129 {
131  fPgen.SetPxPyPzE(0.,0.,0.,0.);
132  for (Int_t npar = 0; npar < fgkNParentsMax; npar++){
133  fParentPDGCode[npar] = -1;
134  fParentPythiaLine[npar] = -1;
135  }
136  for (Int_t i = 0; i < 4; i++){
137  fQuarkPDGCode[i] = -1;
138  fQuarkPythiaLine[i] = -1;
139  }
140  FillFromESD(muonTrack);
141 }
142 
143 //============================================
145 {
147 }
148 
149 //============================================
151 {
152  // check assignment to self
153  if (this == &muonCopy) return *this;
154 
155  // base class assignment
156  TObject::operator=(muonCopy);
157 
158  // assignment operator
159  fPrec = muonCopy.fPrec;
160  fIsTriggered = muonCopy.fIsTriggered;
161  fCharge = muonCopy.fCharge;
162  fChi2 = muonCopy.fChi2;
163  fCentr = muonCopy.fCentr;
164  fPgen = muonCopy.fPgen;
166  fTrackPDGCode = muonCopy.fTrackPDGCode;
167  fOscillation = muonCopy.fOscillation;
168  fNParents = muonCopy.fNParents;
169  fWeight = muonCopy.fWeight;
170 
171  for (Int_t i=0; i<3; i++) fXYZ[i]=muonCopy.fXYZ[i];
172  for (Int_t npar = 0; npar < fgkNParentsMax; npar++){
173  fParentPDGCode[npar] = muonCopy.fParentPDGCode[npar];
174  fParentPythiaLine[npar] = muonCopy.fParentPythiaLine[npar];
175  }
176  for (Int_t i = 0; i < 4; i++){
177  fQuarkPDGCode[i] = muonCopy.fQuarkPDGCode[i];
178  fQuarkPythiaLine[i] = muonCopy.fQuarkPythiaLine[i];
179  }
180 
181  return *this;
182 }
183 
184 //============================================
185 
188  AliMUONTrackParam* trPar = trackReco->GetTrackParamAtVertex();
189  if (!trPar) {
190  AliError("The track must contain the parameters at vertex");
191  return;
192  }
193  this->SetCharge(Int_t(TMath::Sign(1.,trPar->GetInverseBendingMomentum())));
194  this->SetPxPyPz(trPar->Px(),trPar->Py(), trPar->Pz());
195  this->SetTriggered(trackReco->GetMatchTrigger());
196 
197  Double_t xyz[3] = { trPar->GetNonBendingCoor(),
198  trPar->GetBendingCoor(),
199  trPar->GetZ()};
200  if (zvert!=-9999) xyz[2] = zvert;
201  this->SetVertex(xyz);
202 }
203 
204 //============================================
205 void AliMUONTrackLight::FillFromESD(AliESDMuonTrack* muonTrack,Double_t zvert){
207  Double_t mumass = TDatabasePDG::Instance()->GetParticle(13)->Mass();
208  Double_t thetaX = muonTrack->GetThetaX();
209  Double_t thetaY = muonTrack->GetThetaY();
210  Double_t tanthx = TMath::Tan(thetaX);
211  Double_t tanthy = TMath::Tan(thetaY);
212  Double_t pYZ = 1./TMath::Abs(muonTrack->GetInverseBendingMomentum());
213  Double_t pz = - pYZ / TMath::Sqrt(1.0 + tanthy * tanthy);
214  Double_t px = pz * tanthx;
215  Double_t py = pz * tanthy;
216  fCharge = Int_t(TMath::Sign(1.,muonTrack->GetInverseBendingMomentum()));
217  Double_t energy = TMath::Sqrt(mumass * mumass + px*px + py*py + pz*pz);
218  fPrec.SetPxPyPzE(px,py,pz,energy);
219  // get the position
220  fXYZ[0] = muonTrack->GetNonBendingCoor();
221  fXYZ[1] = muonTrack->GetBendingCoor();
222  if (zvert==-9999) fXYZ[2] = muonTrack->GetZ();
223  else fXYZ[2] = zvert;
224  // get the chi2 per d.o.f.
225  fChi2 = muonTrack->GetChi2()/ (2.0 * muonTrack->GetNHit() - 5);
226  fIsTriggered = muonTrack->GetMatchTrigger();
227 }
228 
229 //============================================
230 void AliMUONTrackLight::SetPxPyPz(Double_t px, Double_t py, Double_t pz){
232  Double_t mumass = TDatabasePDG::Instance()->GetParticle(13)->Mass();
233  Double_t energy = TMath::Sqrt(mumass * mumass + px*px + py*py + pz*pz);
234  fPrec.SetPxPyPzE(px,py,pz,energy);
235 }
236 
237 //============================================
238 void AliMUONTrackLight::FillMuonHistory(AliStack *stack, TParticle *part){
240  Int_t countP = -1;
241  Int_t parents[10], parLine[10];
242  Int_t lineM = part->GetFirstMother();//line in the Pythia output of the particle's mother
243 
244  TParticle *mother;
245  Int_t status=-1, pdg=-1;
246  while(lineM >= 0){
247 
248  mother = stack->Particle(lineM); //direct mother of rec. track
249  pdg = mother->GetPdgCode();//store PDG code of first mother
250  // break if a string, gluon, quark or diquark is found
251  if(pdg == 92 || pdg == 21 || TMath::Abs(pdg) < 10 || IsDiquark(pdg)) break;
252  parents[++countP] = pdg;
253  parLine[countP] = lineM;
254  status = mother->GetStatusCode();//get its status code to check if oscillation occured
255  if(IsB0(parents[countP]) && status == 12) this->SetOscillation(kTRUE);
256  lineM = mother->GetFirstMother();
257  }
258  //store all the fragmented parents in an array:
259  for(int i = 0; i <= countP; i++){
260  this->SetParentPDGCode(i,parents[countP-i]);
261  this->SetParentPythiaLine(i,parLine[countP-i]);
262  }
263  fNParents = countP+1;
264  countP = -1;
265 
266  //and store the lines of the string and further quarks in another array:
267  while(lineM >= 0){
268  mother = stack->Particle(lineM);
269  pdg = mother->GetPdgCode();
270  //now, get information before the fragmentation
271  this->SetQuarkPythiaLine(++countP, lineM);//store the line of the string in index 0
272  this->SetQuarkPDGCode(countP, pdg);//store the pdg of the quarks in index 1,2
273  lineM = mother->GetFirstMother();
274  }
275 
276  //check if in case of HF production, the string points to the correct end
277  //and correct it in case of need:
278  countP = 1;
279  for(int par = 0; par < 4; par++){
280  if(TMath::Abs(this->GetQuarkPDGCode(par)) < 6){
281  countP = par; //get the quark just before hadronisation
282  break;
283  }
284  }
285  if(this->GetQuarkPythiaLine(countP) > -1 && (this->GetParentFlavour(0)==4 || this->GetParentFlavour(0)==5)){
286  if(this->GetParentFlavour(0) != TMath::Abs(this->GetQuarkPDGCode(countP))){
287 
288  AliWarning(Form("quark flavour of parent and that of quark do not correspond: %d %d --> correcting\n",
289  this->GetParentFlavour(0), TMath::Abs(this->GetQuarkPDGCode(countP)))
290  );
291 
292  pdg = this->GetQuarkPDGCode(countP);
293  Int_t line = this->GetQuarkPythiaLine(countP);
294  this->ResetQuarkInfo();
295  while(TMath::Abs(pdg) != this->GetParentFlavour(0)){//pdg of q,g in Pythia listing following the wrong string end
296  //must coincide with the flavour of the last fragmented mother
297 
298  pdg = stack->Particle(++line)->GetPdgCode();
299  }
300  //now, we have the correct string end and the correct line
301  //continue to fill again all parents and corresponding lines
302  while(line >= 0){
303  mother = stack->Particle(line);//get again the mother
304  pdg = mother->GetPdgCode();
305  this->SetQuarkPythiaLine(countP, line);
306  this->SetQuarkPDGCode(countP++, pdg);
307  line = mother->GetFirstMother();
308  }
309  this->PrintInfo("h");
310  }//mismatch
311  }
312 }
313 
314 //====================================
317  for(int pos = 1; pos < 4; pos++){//[0] is the string
318  this->SetQuarkPDGCode(pos,-1);
319  this->SetQuarkPythiaLine(pos,-1);
320  }
321 }
322 
323 //====================================
324 Bool_t AliMUONTrackLight::IsB0(Int_t intTest) const {
326  Int_t bMes0[2] = {511,531};//flavour code of B0d and B0s
327  Bool_t answer = kFALSE;
328  for(int i = 0; i < 2; i++){
329  if(TMath::Abs(intTest) == bMes0[i]){
330  answer = kTRUE;
331  break;
332  }
333  }
334  return answer;
335 }
336 //====================================
337 Bool_t AliMUONTrackLight::IsMotherAResonance(Int_t index) const {
339  Int_t intTest = GetParentPDGCode(index);
340  // the resonance pdg code is built this way
341  // x00ffn where x=0,1,.. (1S,2S... states), f=quark flavour
342  Int_t id=intTest%100000;
343  return (!((id-id%10)%110));
344 }
345 //====================================
346 Int_t AliMUONTrackLight::GetParentFlavour(Int_t idParent) const {
349  Int_t pdg = GetParentPDGCode(idParent);
350  Int_t quark = TMath::Abs(pdg/100);
351  if(quark > 9) quark = quark/10;
352  return quark;
353 }
354 
355 //====================================
356 void AliMUONTrackLight::PrintInfo(const Option_t* opt){
361  TString options(opt);
362  options.ToUpper();
363 
364  if(options.Contains("H") || options.Contains("A")){ //muon decay history
365  char *name= new char[100];
366  TString pdg = "", line = "";
367  for(int i = 3; i >= 0; i--){
368  if(this->GetQuarkPythiaLine(i)>= 0){
369  snprintf(name, 100, "%4d --> ", this->GetQuarkPythiaLine(i));
370  line += name;
371  snprintf(name, 100, "%4d --> ", this->GetQuarkPDGCode(i));
372  pdg += name;
373  }
374  }
375  for(int i = 0; i < fNParents; i++){
376  if(this->GetParentPythiaLine(i)>= 0){
377  snprintf(name, 100, "%7d --> ", this->GetParentPythiaLine(i));
378  line += name;
379  snprintf(name, 100, "%7d --> ", this->GetParentPDGCode(i));
380  pdg += name;
381  }
382  }
383  snprintf(name, 100, "%4d", this->GetTrackPythiaLine()); line += name;
384  snprintf(name, 100, "%4d", this->GetTrackPDGCode()); pdg += name;
385 
386  printf("\nmuon's decay history:\n");
387  printf(" PDG: %s\n", pdg.Data());
388  printf("line: %s\n", line.Data());
389  }
390  if(options.Contains("K") || options.Contains("A")){ //muon kinematic
391 
392  Int_t charge = this->GetCharge();
393  Double_t *vtx = this->GetVertex();
394  TLorentzVector momRec = this->GetPRec();
395  TLorentzVector momGen = this->GetPGen();
396  printf("the track's charge is %d\n", charge);
397  printf("Primary vertex: Vx = %1.3f, Vy = %1.3f, Vz = %1.3f\n", vtx[0], vtx[1], vtx[2]);
398  printf("Generated: Px = %1.3f, Py = %1.3f, Pz = %1.3f\n", momGen.Px(), momGen.Py(), momGen.Pz());
399  printf("Reconstructed: Px = %1.3f, Py = %1.3f, Pz = %1.3f\n", momRec.Px(), momRec.Py(), momRec.Pz());
400  printf("Rec. variables: pT %1.3f, pseudo-rapidity %1.3f, theta %1.3f (%1.3f degree), phi %1.3f (%1.3f degree)\n",
401  momRec.Pt(), momRec.Eta(), momRec.Theta(), 180./TMath::Pi() * momRec.Theta(),
402  momRec.Phi(), 180./TMath::Pi() * momRec.Phi());
403  }
404 }
405 //====================================
408  Int_t pdg = this->GetParentPDGCode(idparent);
409  if (TMath::Abs(pdg)==211 || //pi+
410  TMath::Abs(pdg)==321 || //K+
411  TMath::Abs(pdg)==213 || //rho+
412  TMath::Abs(pdg)==311 || //K0
413  TMath::Abs(pdg)==313 || //K*0
414  TMath::Abs(pdg)==323 //K*+
415  ) {
416  return kTRUE;
417  }
418  else return kFALSE;
419 }
420 //====================================
421 Bool_t AliMUONTrackLight::IsDiquark(Int_t pdg) const{
423  pdg = TMath::Abs(pdg);
424  if((pdg > 1000) && (pdg%100 < 10)) return kTRUE;
425  else return kFALSE;
426 }
427 //====================================
428 void AliMUONTrackLight::SetParentPDGCode(Int_t index, Int_t pdg) {
430  if ( index < fgkNParentsMax ) {
431  fParentPDGCode[index] = pdg;
432  } else {
433  AliErrorStream() << "Index outside the array size." << std::endl;
434  }
435 }
436 //====================================
437 void AliMUONTrackLight::SetParentPythiaLine(Int_t index, Int_t line) {
439  if ( index < fgkNParentsMax ) {
440  fParentPythiaLine[index] = line;
441  } else {
442  AliErrorStream() << "Index outside the array size." << std::endl;
443  }
444 }
445 //====================================
446 void AliMUONTrackLight::SetQuarkPDGCode(Int_t index, Int_t pdg){
448  if ( index < 4 ) {
449  fQuarkPDGCode[index] = pdg;
450  } else {
451  AliErrorStream() << "Index outside the array size." << std::endl;
452  }
453 }
454 //====================================
455 void AliMUONTrackLight::SetQuarkPythiaLine(Int_t index, Int_t line){
457  if ( index < 4 ) {
458  fQuarkPythiaLine[index] = line;
459  } else {
460  AliErrorStream() << "Index outside the array size." << std::endl;
461  }
462 }
printf("Chi2/npoints = %f\n", TMath::Sqrt(chi2/npoints))
void FillFromAliMUONTrack(AliMUONTrack *trackReco, Double_t zvert=-9999)
Int_t fTrackPDGCode
pdg code of the rec. track (in general will be a muon)
TLorentzVector fPgen
4-momentum of the generated particle
void SetQuarkPDGCode(Int_t index, Int_t pdg)
Set pdg of the string [0], quarks/gluons [1,2], sometimes proton [3].
Int_t GetQuarkPythiaLine(Int_t index=0) const
Return line of Pythia output for string [0] and quarks [1,2], sometimes proton [3].
Int_t GetMatchTrigger(void) const
return 1,2,3 if track matches with trigger track, 0 if not
Definition: AliMUONTrack.h:79
Bool_t IsB0(Int_t intTest) const
Double_t GetBendingCoor() const
return bending coordinate (cm)
TLorentzVector GetPRec() const
Return reconstructed 4-momentum.
AliMUONTrackLight & operator=(const AliMUONTrackLight &)
Bool_t fOscillation
flag for oscillation
ClassImp(AliMUONTrackLight) AliMUONTrackLight
Bool_t IsDiquark(Int_t pdg) const
for(Int_t isec=54;isec< 71;isec+=1)
Definition: AnalyzeLaser.C:32
Double_t GetZ() const
return Z coordinate (cm)
void SetTriggered(Bool_t isTriggered)
Set flag for trigger.
Int_t fNParents
acually filled no. of fragmented parents
Track parameters in ALICE dimuon spectrometer.
Int_t GetTrackPythiaLine() const
Return line of Pythia output for string [0] and quarks [1,2], sometimes proton [3].
Double_t fChi2
chi2 / ndf in the MUON track fit
void SetVertex(Double_t *xyz)
Set primary vertex position from the ITS.
Double_t GetInverseBendingMomentum() const
return inverse bending momentum (GeV/c ** -1) times the charge (assumed forward motion) ...
Int_t fQuarkPDGCode[4]
pdg of the string [0], quarks/gluons [1,2], sometimes proton [3]
Int_t fParentPythiaLine[fgkNParentsMax]
line of Pythia output for hadronised parents & grandparents
Int_t fParentPDGCode[fgkNParentsMax]
hadronised parents and grandparents
Int_t GetQuarkPDGCode(Int_t index=0) const
Return pdg of the string [0], quarks/gluons [1,2], sometimes proton [3].
void SetParentPDGCode(Int_t index, Int_t pdg)
Set hadronised parents and grandparents.
void SetCharge(Int_t charge)
Set muon charge.
Double_t Py() const
Double_t * GetVertex()
Return primary vertex position from the ITS.
Int_t fTrackPythiaLine
line of kin. stack where rec. track (in general, the muon) is stored
Bool_t IsMotherAResonance(Int_t index=0) const
void SetOscillation(Bool_t oscillation)
Set flag for oscillation.
Bool_t IsParentPionOrKaon(Int_t idParent=0)
Float_t fCentr
centrality
Int_t GetParentFlavour(Int_t idParent=0) const
TLorentzVector fPrec
reconstructed 4-momentum
Double_t Px() const
virtual void PrintInfo(const Option_t *opt)
Compact information for the muon generated tracks.
void SetPxPyPz(Double_t px, Double_t py, Double_t pz)
void FillFromESD(AliESDMuonTrack *muonTrack, Double_t zvert=-9999)
Double_t GetNonBendingCoor() const
return non bending coordinate (cm)
Int_t GetCharge() const
Return muon charge.
Int_t fQuarkPythiaLine[4]
line of Pythia output for string [0] and quarks [1,2], sometimes proton [3]
void SetParentPythiaLine(Int_t index, Int_t line)
Set line of Pythia output for hadronised parents & grandparents.
Double_t Pz() const
TLorentzVector GetPGen() const
Return 4-momentum of the generated particle.
Int_t fCharge
muon charge
Reconstructed track in ALICE dimuon spectrometer.
Definition: AliMUONTrack.h:24
Double_t fXYZ[3]
primary vertex position from the ITS
Int_t GetParentPDGCode(Int_t index=0) const
Return hadronised parents and grandparents.
void FillMuonHistory(AliStack *stack, TParticle *part)
Int_t GetParentPythiaLine(Int_t index=0) const
Return line of Pythia output for hadronised parents & grandparents.
Int_t GetTrackPDGCode() const
Return pdg code of the rec. track (in general will be a muon)
void SetQuarkPythiaLine(Int_t index, Int_t line)
Set line of Pythia output for string [0] and quarks [1,2], sometimes proton [3].
Bool_t fIsTriggered
flag for trigger
static const Int_t fgkNParentsMax
maximum number of parents
AliMUONTrackParam * GetTrackParamAtVertex() const
return pointer to track parameters at vertex (can be 0x0)
Definition: AliMUONTrack.h:98