- Adding check and flagging for HG present
[u/mrichter/AliRoot.git] / MUON / AliMUONTrackLight.cxx
CommitLineData
55fd51b0 1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
b88403f3 3 * SigmaEffect_thetadegrees *
55fd51b0 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
3d1463c8 18//-----------------------------------------------------------------------------
55fd51b0 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
3d1463c8 27//
28// This class was prepared by INFN Cagliari, July 2006
29// (authors: H.Woehri, A.de Falco)
30//-----------------------------------------------------------------------------
55fd51b0 31
9bdda5f6 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
55fd51b0 39#include "AliMUONTrackLight.h"
40#include "AliMUONTrack.h"
41#include "AliMUONConstants.h"
bd5eb041 42#include "AliMUONVTrackStore.h"
933daf4c 43#include "AliMUONTrackParam.h"
55fd51b0 44
45#include "AliESDMuonTrack.h"
55fd51b0 46#include "AliStack.h"
f202486b 47#include "AliLog.h"
55fd51b0 48
49#include "TDatabasePDG.h"
55fd51b0 50#include "TParticle.h"
51#include "TString.h"
52
53ClassImp(AliMUONTrackLight)
54
55//===================================================================
56
57AliMUONTrackLight::AliMUONTrackLight()
58 : TObject(),
59 fPrec(),
60 fIsTriggered(kFALSE),
61 fCharge(-999),
b88403f3 62 fChi2(-1),
55fd51b0 63 fCentr(-1),
64 fPgen(),
65 fTrackPythiaLine(-999),
66 fTrackPDGCode(-999),
67 fOscillation(kFALSE),
b88403f3 68 fNParents(0),
69 fWeight(1)
55fd51b0 70{
71 /// default constructor
72 fPgen.SetPxPyPzE(0.,0.,0.,0.);
73 fPrec.SetPxPyPzE(0.,0.,0.,0.);
74 for (Int_t i=0; i<3; i++) fXYZ[i]=-999;
75 for (Int_t npar = 0; npar < fgkNParentsMax; npar++){
76 fParentPDGCode[npar] = -1;
77 fParentPythiaLine[npar] = -1;
78 }
79 for (Int_t i = 0; i < 4; i++){
80 fQuarkPDGCode[i] = -1;
81 fQuarkPythiaLine[i] = -1;
82 }
83}
84
85//============================================
86AliMUONTrackLight::AliMUONTrackLight(const AliMUONTrackLight &muonCopy)
87 : TObject(muonCopy),
88 fPrec(muonCopy.fPrec),
89 fIsTriggered(muonCopy.fIsTriggered),
90 fCharge(muonCopy.fCharge),
b88403f3 91 fChi2(muonCopy.fChi2),
55fd51b0 92 fCentr(muonCopy.fCentr),
93 fPgen(muonCopy.fPgen),
94 fTrackPythiaLine(muonCopy.fTrackPythiaLine),
95 fTrackPDGCode(muonCopy.fTrackPDGCode),
96 fOscillation(muonCopy.fOscillation),
b88403f3 97 fNParents(muonCopy.fNParents),
98 fWeight(muonCopy.fWeight)
55fd51b0 99{
100 /// copy constructor
101 for (Int_t i=0; i<3; i++) fXYZ[i]=muonCopy.fXYZ[i];
102 for (Int_t npar = 0; npar < fgkNParentsMax; npar++){
103 fParentPDGCode[npar] = muonCopy.fParentPDGCode[npar];
104 fParentPythiaLine[npar] = muonCopy.fParentPythiaLine[npar];
105 }
106 for (Int_t i = 0; i < 4; i++){
107 fQuarkPDGCode[i] = muonCopy.fQuarkPDGCode[i];
108 fQuarkPythiaLine[i] = muonCopy.fQuarkPythiaLine[i];
109 }
110}
111
112//============================================
113AliMUONTrackLight::AliMUONTrackLight(AliESDMuonTrack* muonTrack)
114 : TObject(),
115 fPrec(),
116 fIsTriggered(kFALSE),
117 fCharge(-999),
b88403f3 118 fChi2(-1),
55fd51b0 119 fCentr(-1),
120 fPgen(),
121 fTrackPythiaLine(-999),
122 fTrackPDGCode(-999),
123 fOscillation(kFALSE),
b88403f3 124 fNParents(0),
125 fWeight(1)
55fd51b0 126{
127 /// constructor
f202486b 128 fPgen.SetPxPyPzE(0.,0.,0.,0.);
129 for (Int_t npar = 0; npar < fgkNParentsMax; npar++){
130 fParentPDGCode[npar] = -1;
131 fParentPythiaLine[npar] = -1;
132 }
133 for (Int_t i = 0; i < 4; i++){
134 fQuarkPDGCode[i] = -1;
135 fQuarkPythiaLine[i] = -1;
136 }
137 FillFromESD(muonTrack);
55fd51b0 138}
139
140//============================================
b88403f3 141AliMUONTrackLight::~AliMUONTrackLight()
142{
143/// Destructor
144}
145
146//============================================
147
148void AliMUONTrackLight::FillFromAliMUONTrack(AliMUONTrack *trackReco,Double_t zvert){
149 /// this method sets the muon reconstructed momentum according to the value given by AliMUONTrack
f202486b 150 AliMUONTrackParam* trPar = trackReco->GetTrackParamAtVertex();
151 if (!trPar) {
152 AliError("The track must contain the parameters at vertex");
153 return;
61fed964 154 }
f202486b 155 this->SetCharge(Int_t(TMath::Sign(1.,trPar->GetInverseBendingMomentum())));
156 this->SetPxPyPz(trPar->Px(),trPar->Py(), trPar->Pz());
b88403f3 157 this->SetTriggered(trackReco->GetMatchTrigger());
158
f202486b 159 Double_t xyz[3] = { trPar->GetNonBendingCoor(),
160 trPar->GetBendingCoor(),
161 trPar->GetZ()};
61fed964 162 if (zvert!=-9999) xyz[2] = zvert;
b88403f3 163 this->SetVertex(xyz);
164}
165
166//============================================
167void AliMUONTrackLight::FillFromESD(AliESDMuonTrack* muonTrack,Double_t zvert){
55fd51b0 168 /// computes prec and charge from ESD track
169 Double_t mumass = TDatabasePDG::Instance()->GetParticle(13)->Mass();
170 Double_t thetaX = muonTrack->GetThetaX();
171 Double_t thetaY = muonTrack->GetThetaY();
172 Double_t tanthx = TMath::Tan(thetaX);
173 Double_t tanthy = TMath::Tan(thetaY);
174 Double_t pYZ = 1./TMath::Abs(muonTrack->GetInverseBendingMomentum());
175 Double_t pz = - pYZ / TMath::Sqrt(1.0 + tanthy * tanthy);
176 Double_t px = pz * tanthx;
177 Double_t py = pz * tanthy;
178 fCharge = Int_t(TMath::Sign(1.,muonTrack->GetInverseBendingMomentum()));
179 Double_t energy = TMath::Sqrt(mumass * mumass + px*px + py*py + pz*pz);
180 fPrec.SetPxPyPzE(px,py,pz,energy);
b88403f3 181 // get the position
182 fXYZ[0] = muonTrack->GetNonBendingCoor();
183 fXYZ[1] = muonTrack->GetBendingCoor();
184 if (zvert==-9999) fXYZ[2] = muonTrack->GetZ();
185 else fXYZ[2] = zvert;
186 // get the chi2 per d.o.f.
187 fChi2 = muonTrack->GetChi2()/ (2.0 * muonTrack->GetNHit() - 5);
188 fIsTriggered = muonTrack->GetMatchTrigger();
55fd51b0 189}
190
191//============================================
192void AliMUONTrackLight::SetPxPyPz(Double_t px, Double_t py, Double_t pz){
193 /// set the reconstructed 4-momentum, assuming the particle is a muon
194 Double_t mumass = TDatabasePDG::Instance()->GetParticle(13)->Mass();
195 Double_t energy = TMath::Sqrt(mumass * mumass + px*px + py*py + pz*pz);
196 fPrec.SetPxPyPzE(px,py,pz,energy);
197}
198
199//============================================
b88403f3 200void AliMUONTrackLight::FillMuonHistory(AliStack *stack, TParticle *part){
201 /// scans the muon history to determine parents pdg code and pythia line
55fd51b0 202 Int_t countP = -1;
203 Int_t parents[10], parLine[10];
204 Int_t lineM = part->GetFirstMother();//line in the Pythia output of the particle's mother
205
55fd51b0 206 TParticle *mother;
207 Int_t status=-1, pdg=-1;
208 while(lineM >= 0){
b88403f3 209
55fd51b0 210 mother = stack->Particle(lineM); //direct mother of rec. track
211 pdg = mother->GetPdgCode();//store PDG code of first mother
b88403f3 212 // break if a string, gluon, quark or diquark is found
213 if(pdg == 92 || pdg == 21 || TMath::Abs(pdg) < 10 || IsDiquark(pdg)) break;
55fd51b0 214 parents[++countP] = pdg;
215 parLine[countP] = lineM;
216 status = mother->GetStatusCode();//get its status code to check if oscillation occured
217 if(IsB0(parents[countP]) && status == 12) this->SetOscillation(kTRUE);
218 lineM = mother->GetFirstMother();
219 }
220 //store all the fragmented parents in an array:
221 for(int i = 0; i <= countP; i++){
222 this->SetParentPDGCode(i,parents[countP-i]);
223 this->SetParentPythiaLine(i,parLine[countP-i]);
224 }
55fd51b0 225 fNParents = countP+1;
226 countP = -1;
b88403f3 227
55fd51b0 228 //and store the lines of the string and further quarks in another array:
229 while(lineM >= 0){
230 mother = stack->Particle(lineM);
231 pdg = mother->GetPdgCode();
232 //now, get information before the fragmentation
233 this->SetQuarkPythiaLine(++countP, lineM);//store the line of the string in index 0
234 this->SetQuarkPDGCode(countP, pdg);//store the pdg of the quarks in index 1,2
235 lineM = mother->GetFirstMother();
236 }
b88403f3 237
55fd51b0 238 //check if in case of HF production, the string points to the correct end
239 //and correct it in case of need:
240 countP = 1;
b88403f3 241 for(int par = 0; par < 4; par++){
242 if(TMath::Abs(this->GetQuarkPDGCode(par)) < 6){
243 countP = par; //get the quark just before hadronisation
244 break;
245 }
246 }
55fd51b0 247 if(this->GetQuarkPythiaLine(countP) > -1 && (this->GetParentFlavour(0)==4 || this->GetParentFlavour(0)==5)){
248 if(this->GetParentFlavour(0) != TMath::Abs(this->GetQuarkPDGCode(countP))){
249
9bdda5f6 250 AliWarning(Form("quark flavour of parent and that of quark do not correspond: %d %d --> correcting\n",
251 this->GetParentFlavour(0), TMath::Abs(this->GetQuarkPDGCode(countP)))
252 );
253
7cdc832b 254 pdg = this->GetQuarkPDGCode(countP);
255 Int_t line = this->GetQuarkPythiaLine(countP);
55fd51b0 256 this->ResetQuarkInfo();
257 while(TMath::Abs(pdg) != this->GetParentFlavour(0)){//pdg of q,g in Pythia listing following the wrong string end
258 //must coincide with the flavour of the last fragmented mother
259
260 pdg = stack->Particle(++line)->GetPdgCode();
261 }
262 //now, we have the correct string end and the correct line
263 //continue to fill again all parents and corresponding lines
264 while(line >= 0){
265 mother = stack->Particle(line);//get again the mother
266 pdg = mother->GetPdgCode();
267 this->SetQuarkPythiaLine(countP, line);
268 this->SetQuarkPDGCode(countP++, pdg);
269 line = mother->GetFirstMother();
270 }
b88403f3 271 this->PrintInfo("h");
272 }//mismatch
55fd51b0 273 }
274}
b88403f3 275
55fd51b0 276//====================================
277void AliMUONTrackLight::ResetQuarkInfo(){
278 /// resets parton information
279 for(int pos = 1; pos < 4; pos++){//[0] is the string
280 this->SetQuarkPDGCode(pos,-1);
281 this->SetQuarkPythiaLine(pos,-1);
282 }
283}
284
285//====================================
286Bool_t AliMUONTrackLight::IsB0(Int_t intTest) const {
287 /// checks if the particle is a B0
288 Int_t bMes0[2] = {511,531};//flavour code of B0d and B0s
289 Bool_t answer = kFALSE;
290 for(int i = 0; i < 2; i++){
291 if(TMath::Abs(intTest) == bMes0[i]){
292 answer = kTRUE;
293 break;
294 }
295 }
296 return answer;
297}
298//====================================
299Bool_t AliMUONTrackLight::IsMotherAResonance(Int_t index) const {
300 /// checks if mother is a resonance
301 Int_t intTest = GetParentPDGCode(index);
302 // the resonance pdg code is built this way
303 // x00ffn where x=0,1,.. (1S,2S... states), f=quark flavour
304 Int_t id=intTest%100000;
305 return (!((id-id%10)%110));
306}
307//====================================
308Int_t AliMUONTrackLight::GetParentFlavour(Int_t idParent) const {
309 /// returns the flavour of parent idParent (idParent=0 is the oldest
310 /// hadronized parent)
311 Int_t pdg = GetParentPDGCode(idParent);
312 Int_t quark = TMath::Abs(pdg/100);
313 if(quark > 9) quark = quark/10;
314 return quark;
315}
316
317//====================================
57e2ad1a 318void AliMUONTrackLight::PrintInfo(const Option_t* opt){
55fd51b0 319 /// prints information about the track:
320 /// - "H" muon's decay history
321 /// - "K" muon kinematics
322 /// - "A" all variables
323 TString options(opt);
324 options.ToUpper();
325
326 if(options.Contains("H") || options.Contains("A")){ //muon decay history
327 char *name= new char[100];
328 TString pdg = "", line = "";
329 for(int i = 3; i >= 0; i--){
330 if(this->GetQuarkPythiaLine(i)>= 0){
331 sprintf(name, "%4d --> ", this->GetQuarkPythiaLine(i));
332 line += name;
333 sprintf(name, "%4d --> ", this->GetQuarkPDGCode(i));
334 pdg += name;
335 }
336 }
337 for(int i = 0; i < fNParents; i++){
338 if(this->GetParentPythiaLine(i)>= 0){
339 sprintf(name, "%7d --> ", this->GetParentPythiaLine(i));
340 line += name;
341 sprintf(name, "%7d --> ", this->GetParentPDGCode(i));
342 pdg += name;
343 }
344 }
345 sprintf(name, "%4d", this->GetTrackPythiaLine()); line += name;
346 sprintf(name, "%4d", this->GetTrackPDGCode()); pdg += name;
347
348 printf("\nmuon's decay history:\n");
349 printf(" PDG: %s\n", pdg.Data());
350 printf("line: %s\n", line.Data());
351 }
352 if(options.Contains("K") || options.Contains("A")){ //muon kinematic
353
354 Int_t charge = this->GetCharge();
355 Double_t *vtx = this->GetVertex();
356 TLorentzVector momRec = this->GetPRec();
357 TLorentzVector momGen = this->GetPGen();
358 printf("the track's charge is %d\n", charge);
359 printf("Primary vertex: Vx = %1.3f, Vy = %1.3f, Vz = %1.3f\n", vtx[0], vtx[1], vtx[2]);
360 printf("Generated: Px = %1.3f, Py = %1.3f, Pz = %1.3f\n", momGen.Px(), momGen.Py(), momGen.Pz());
361 printf("Reconstructed: Px = %1.3f, Py = %1.3f, Pz = %1.3f\n", momRec.Px(), momRec.Py(), momRec.Pz());
362 printf("Rec. variables: pT %1.3f, pseudo-rapidity %1.3f, theta %1.3f (%1.3f degree), phi %1.3f (%1.3f degree)\n",
363 momRec.Pt(), momRec.Eta(), momRec.Theta(), 180./TMath::Pi() * momRec.Theta(),
364 momRec.Phi(), 180./TMath::Pi() * momRec.Phi());
365 }
366}
b88403f3 367//====================================
368Bool_t AliMUONTrackLight::IsParentPionOrKaon(Int_t idparent){
369 /// checks if a muon comes from a pion or kaon or a particle that decays into one of these two
370 Int_t pdg = this->GetParentPDGCode(idparent);
371 if (TMath::Abs(pdg)==211 || //pi+
372 TMath::Abs(pdg)==321 || //K+
373 TMath::Abs(pdg)==213 || //rho+
374 TMath::Abs(pdg)==311 || //K0
375 TMath::Abs(pdg)==313 || //K*0
376 TMath::Abs(pdg)==323 //K*+
377 ) {
378 return kTRUE;
379 }
380 else return kFALSE;
381}
382//====================================
9ee1d6ff 383Bool_t AliMUONTrackLight::IsDiquark(Int_t pdg) const{
b88403f3 384 /// check if the provided pdg code corresponds to a diquark
385 pdg = TMath::Abs(pdg);
386 if((pdg > 1000) && (pdg%100 < 10)) return kTRUE;
387 else return kFALSE;
388}