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. *
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 **************************************************************************/
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
26 // To be used together with AliMUONPairLight
28 // This class was prepared by INFN Cagliari, July 2006
29 // (authors: H.Woehri, A.de Falco)
30 //-----------------------------------------------------------------------------
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>
39 #include "AliMUONTrackLight.h"
40 #include "AliMUONTrack.h"
41 #include "AliMUONConstants.h"
42 #include "AliMUONVTrackStore.h"
43 #include "AliMUONTrackParam.h"
45 #include "AliESDMuonTrack.h"
49 #include "TDatabasePDG.h"
50 #include "TParticle.h"
55 ClassImp(AliMUONTrackLight)
57 //===================================================================
59 AliMUONTrackLight::AliMUONTrackLight()
67 fTrackPythiaLine(-999),
73 /// default constructor
74 fPgen.SetPxPyPzE(0.,0.,0.,0.);
75 fPrec.SetPxPyPzE(0.,0.,0.,0.);
76 for (Int_t i=0; i<3; i++) fXYZ[i]=-999;
77 for (Int_t npar = 0; npar < fgkNParentsMax; npar++){
78 fParentPDGCode[npar] = -1;
79 fParentPythiaLine[npar] = -1;
81 for (Int_t i = 0; i < 4; i++){
82 fQuarkPDGCode[i] = -1;
83 fQuarkPythiaLine[i] = -1;
87 //============================================
88 AliMUONTrackLight::AliMUONTrackLight(const AliMUONTrackLight &muonCopy)
90 fPrec(muonCopy.fPrec),
91 fIsTriggered(muonCopy.fIsTriggered),
92 fCharge(muonCopy.fCharge),
93 fChi2(muonCopy.fChi2),
94 fCentr(muonCopy.fCentr),
95 fPgen(muonCopy.fPgen),
96 fTrackPythiaLine(muonCopy.fTrackPythiaLine),
97 fTrackPDGCode(muonCopy.fTrackPDGCode),
98 fOscillation(muonCopy.fOscillation),
99 fNParents(muonCopy.fNParents),
100 fWeight(muonCopy.fWeight)
103 for (Int_t i=0; i<3; i++) fXYZ[i]=muonCopy.fXYZ[i];
104 for (Int_t npar = 0; npar < fgkNParentsMax; npar++){
105 fParentPDGCode[npar] = muonCopy.fParentPDGCode[npar];
106 fParentPythiaLine[npar] = muonCopy.fParentPythiaLine[npar];
108 for (Int_t i = 0; i < 4; i++){
109 fQuarkPDGCode[i] = muonCopy.fQuarkPDGCode[i];
110 fQuarkPythiaLine[i] = muonCopy.fQuarkPythiaLine[i];
114 //============================================
115 AliMUONTrackLight::AliMUONTrackLight(AliESDMuonTrack* muonTrack)
118 fIsTriggered(kFALSE),
123 fTrackPythiaLine(-999),
125 fOscillation(kFALSE),
130 fPgen.SetPxPyPzE(0.,0.,0.,0.);
131 for (Int_t npar = 0; npar < fgkNParentsMax; npar++){
132 fParentPDGCode[npar] = -1;
133 fParentPythiaLine[npar] = -1;
135 for (Int_t i = 0; i < 4; i++){
136 fQuarkPDGCode[i] = -1;
137 fQuarkPythiaLine[i] = -1;
139 FillFromESD(muonTrack);
142 //============================================
143 AliMUONTrackLight::~AliMUONTrackLight()
148 //============================================
149 AliMUONTrackLight& AliMUONTrackLight::operator=(const AliMUONTrackLight& muonCopy)
151 // check assignment to self
152 if (this == &muonCopy) return *this;
154 // base class assignment
155 TObject::operator=(muonCopy);
157 // assignment operator
158 fPrec = muonCopy.fPrec;
159 fIsTriggered = muonCopy.fIsTriggered;
160 fCharge = muonCopy.fCharge;
161 fChi2 = muonCopy.fChi2;
162 fCentr = muonCopy.fCentr;
163 fPgen = muonCopy.fPgen;
164 fTrackPythiaLine = muonCopy.fTrackPythiaLine;
165 fTrackPDGCode = muonCopy.fTrackPDGCode;
166 fOscillation = muonCopy.fOscillation;
167 fNParents = muonCopy.fNParents;
168 fWeight = muonCopy.fWeight;
170 for (Int_t i=0; i<3; i++) fXYZ[i]=muonCopy.fXYZ[i];
171 for (Int_t npar = 0; npar < fgkNParentsMax; npar++){
172 fParentPDGCode[npar] = muonCopy.fParentPDGCode[npar];
173 fParentPythiaLine[npar] = muonCopy.fParentPythiaLine[npar];
175 for (Int_t i = 0; i < 4; i++){
176 fQuarkPDGCode[i] = muonCopy.fQuarkPDGCode[i];
177 fQuarkPythiaLine[i] = muonCopy.fQuarkPythiaLine[i];
181 //============================================
183 void AliMUONTrackLight::FillFromAliMUONTrack(AliMUONTrack *trackReco,Double_t zvert){
184 /// this method sets the muon reconstructed momentum according to the value given by AliMUONTrack
185 AliMUONTrackParam* trPar = trackReco->GetTrackParamAtVertex();
187 AliError("The track must contain the parameters at vertex");
190 this->SetCharge(Int_t(TMath::Sign(1.,trPar->GetInverseBendingMomentum())));
191 this->SetPxPyPz(trPar->Px(),trPar->Py(), trPar->Pz());
192 this->SetTriggered(trackReco->GetMatchTrigger());
194 Double_t xyz[3] = { trPar->GetNonBendingCoor(),
195 trPar->GetBendingCoor(),
197 if (zvert!=-9999) xyz[2] = zvert;
198 this->SetVertex(xyz);
201 //============================================
202 void AliMUONTrackLight::FillFromESD(AliESDMuonTrack* muonTrack,Double_t zvert){
203 /// computes prec and charge from ESD track
204 Double_t mumass = TDatabasePDG::Instance()->GetParticle(13)->Mass();
205 Double_t thetaX = muonTrack->GetThetaX();
206 Double_t thetaY = muonTrack->GetThetaY();
207 Double_t tanthx = TMath::Tan(thetaX);
208 Double_t tanthy = TMath::Tan(thetaY);
209 Double_t pYZ = 1./TMath::Abs(muonTrack->GetInverseBendingMomentum());
210 Double_t pz = - pYZ / TMath::Sqrt(1.0 + tanthy * tanthy);
211 Double_t px = pz * tanthx;
212 Double_t py = pz * tanthy;
213 fCharge = Int_t(TMath::Sign(1.,muonTrack->GetInverseBendingMomentum()));
214 Double_t energy = TMath::Sqrt(mumass * mumass + px*px + py*py + pz*pz);
215 fPrec.SetPxPyPzE(px,py,pz,energy);
217 fXYZ[0] = muonTrack->GetNonBendingCoor();
218 fXYZ[1] = muonTrack->GetBendingCoor();
219 if (zvert==-9999) fXYZ[2] = muonTrack->GetZ();
220 else fXYZ[2] = zvert;
221 // get the chi2 per d.o.f.
222 fChi2 = muonTrack->GetChi2()/ (2.0 * muonTrack->GetNHit() - 5);
223 fIsTriggered = muonTrack->GetMatchTrigger();
226 //============================================
227 void AliMUONTrackLight::SetPxPyPz(Double_t px, Double_t py, Double_t pz){
228 /// set the reconstructed 4-momentum, assuming the particle is a muon
229 Double_t mumass = TDatabasePDG::Instance()->GetParticle(13)->Mass();
230 Double_t energy = TMath::Sqrt(mumass * mumass + px*px + py*py + pz*pz);
231 fPrec.SetPxPyPzE(px,py,pz,energy);
234 //============================================
235 void AliMUONTrackLight::FillMuonHistory(AliStack *stack, TParticle *part){
236 /// scans the muon history to determine parents pdg code and pythia line
238 Int_t parents[10], parLine[10];
239 Int_t lineM = part->GetFirstMother();//line in the Pythia output of the particle's mother
242 Int_t status=-1, pdg=-1;
245 mother = stack->Particle(lineM); //direct mother of rec. track
246 pdg = mother->GetPdgCode();//store PDG code of first mother
247 // break if a string, gluon, quark or diquark is found
248 if(pdg == 92 || pdg == 21 || TMath::Abs(pdg) < 10 || IsDiquark(pdg)) break;
249 parents[++countP] = pdg;
250 parLine[countP] = lineM;
251 status = mother->GetStatusCode();//get its status code to check if oscillation occured
252 if(IsB0(parents[countP]) && status == 12) this->SetOscillation(kTRUE);
253 lineM = mother->GetFirstMother();
255 //store all the fragmented parents in an array:
256 for(int i = 0; i <= countP; i++){
257 this->SetParentPDGCode(i,parents[countP-i]);
258 this->SetParentPythiaLine(i,parLine[countP-i]);
260 fNParents = countP+1;
263 //and store the lines of the string and further quarks in another array:
265 mother = stack->Particle(lineM);
266 pdg = mother->GetPdgCode();
267 //now, get information before the fragmentation
268 this->SetQuarkPythiaLine(++countP, lineM);//store the line of the string in index 0
269 this->SetQuarkPDGCode(countP, pdg);//store the pdg of the quarks in index 1,2
270 lineM = mother->GetFirstMother();
273 //check if in case of HF production, the string points to the correct end
274 //and correct it in case of need:
276 for(int par = 0; par < 4; par++){
277 if(TMath::Abs(this->GetQuarkPDGCode(par)) < 6){
278 countP = par; //get the quark just before hadronisation
282 if(this->GetQuarkPythiaLine(countP) > -1 && (this->GetParentFlavour(0)==4 || this->GetParentFlavour(0)==5)){
283 if(this->GetParentFlavour(0) != TMath::Abs(this->GetQuarkPDGCode(countP))){
285 AliWarning(Form("quark flavour of parent and that of quark do not correspond: %d %d --> correcting\n",
286 this->GetParentFlavour(0), TMath::Abs(this->GetQuarkPDGCode(countP)))
289 pdg = this->GetQuarkPDGCode(countP);
290 Int_t line = this->GetQuarkPythiaLine(countP);
291 this->ResetQuarkInfo();
292 while(TMath::Abs(pdg) != this->GetParentFlavour(0)){//pdg of q,g in Pythia listing following the wrong string end
293 //must coincide with the flavour of the last fragmented mother
295 pdg = stack->Particle(++line)->GetPdgCode();
297 //now, we have the correct string end and the correct line
298 //continue to fill again all parents and corresponding lines
300 mother = stack->Particle(line);//get again the mother
301 pdg = mother->GetPdgCode();
302 this->SetQuarkPythiaLine(countP, line);
303 this->SetQuarkPDGCode(countP++, pdg);
304 line = mother->GetFirstMother();
306 this->PrintInfo("h");
311 //====================================
312 void AliMUONTrackLight::ResetQuarkInfo(){
313 /// resets parton information
314 for(int pos = 1; pos < 4; pos++){//[0] is the string
315 this->SetQuarkPDGCode(pos,-1);
316 this->SetQuarkPythiaLine(pos,-1);
320 //====================================
321 Bool_t AliMUONTrackLight::IsB0(Int_t intTest) const {
322 /// checks if the particle is a B0
323 Int_t bMes0[2] = {511,531};//flavour code of B0d and B0s
324 Bool_t answer = kFALSE;
325 for(int i = 0; i < 2; i++){
326 if(TMath::Abs(intTest) == bMes0[i]){
333 //====================================
334 Bool_t AliMUONTrackLight::IsMotherAResonance(Int_t index) const {
335 /// checks if mother is a resonance
336 Int_t intTest = GetParentPDGCode(index);
337 // the resonance pdg code is built this way
338 // x00ffn where x=0,1,.. (1S,2S... states), f=quark flavour
339 Int_t id=intTest%100000;
340 return (!((id-id%10)%110));
342 //====================================
343 Int_t AliMUONTrackLight::GetParentFlavour(Int_t idParent) const {
344 /// returns the flavour of parent idParent (idParent=0 is the oldest
345 /// hadronized parent)
346 Int_t pdg = GetParentPDGCode(idParent);
347 Int_t quark = TMath::Abs(pdg/100);
348 if(quark > 9) quark = quark/10;
352 //====================================
353 void AliMUONTrackLight::PrintInfo(const Option_t* opt){
354 /// prints information about the track:
355 /// - "H" muon's decay history
356 /// - "K" muon kinematics
357 /// - "A" all variables
358 TString options(opt);
361 if(options.Contains("H") || options.Contains("A")){ //muon decay history
362 char *name= new char[100];
363 TString pdg = "", line = "";
364 for(int i = 3; i >= 0; i--){
365 if(this->GetQuarkPythiaLine(i)>= 0){
366 snprintf(name, 100, "%4d --> ", this->GetQuarkPythiaLine(i));
368 snprintf(name, 100, "%4d --> ", this->GetQuarkPDGCode(i));
372 for(int i = 0; i < fNParents; i++){
373 if(this->GetParentPythiaLine(i)>= 0){
374 snprintf(name, 100, "%7d --> ", this->GetParentPythiaLine(i));
376 snprintf(name, 100, "%7d --> ", this->GetParentPDGCode(i));
380 snprintf(name, 100, "%4d", this->GetTrackPythiaLine()); line += name;
381 snprintf(name, 100, "%4d", this->GetTrackPDGCode()); pdg += name;
383 printf("\nmuon's decay history:\n");
384 printf(" PDG: %s\n", pdg.Data());
385 printf("line: %s\n", line.Data());
387 if(options.Contains("K") || options.Contains("A")){ //muon kinematic
389 Int_t charge = this->GetCharge();
390 Double_t *vtx = this->GetVertex();
391 TLorentzVector momRec = this->GetPRec();
392 TLorentzVector momGen = this->GetPGen();
393 printf("the track's charge is %d\n", charge);
394 printf("Primary vertex: Vx = %1.3f, Vy = %1.3f, Vz = %1.3f\n", vtx[0], vtx[1], vtx[2]);
395 printf("Generated: Px = %1.3f, Py = %1.3f, Pz = %1.3f\n", momGen.Px(), momGen.Py(), momGen.Pz());
396 printf("Reconstructed: Px = %1.3f, Py = %1.3f, Pz = %1.3f\n", momRec.Px(), momRec.Py(), momRec.Pz());
397 printf("Rec. variables: pT %1.3f, pseudo-rapidity %1.3f, theta %1.3f (%1.3f degree), phi %1.3f (%1.3f degree)\n",
398 momRec.Pt(), momRec.Eta(), momRec.Theta(), 180./TMath::Pi() * momRec.Theta(),
399 momRec.Phi(), 180./TMath::Pi() * momRec.Phi());
402 //====================================
403 Bool_t AliMUONTrackLight::IsParentPionOrKaon(Int_t idparent){
404 /// checks if a muon comes from a pion or kaon or a particle that decays into one of these two
405 Int_t pdg = this->GetParentPDGCode(idparent);
406 if (TMath::Abs(pdg)==211 || //pi+
407 TMath::Abs(pdg)==321 || //K+
408 TMath::Abs(pdg)==213 || //rho+
409 TMath::Abs(pdg)==311 || //K0
410 TMath::Abs(pdg)==313 || //K*0
411 TMath::Abs(pdg)==323 //K*+
417 //====================================
418 Bool_t AliMUONTrackLight::IsDiquark(Int_t pdg) const{
419 /// check if the provided pdg code corresponds to a diquark
420 pdg = TMath::Abs(pdg);
421 if((pdg > 1000) && (pdg%100 < 10)) return kTRUE;