Include Riostream.h
[u/mrichter/AliRoot.git] / MUON / MUONevaluation / 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
b785c4bd 53#include <cstdio>
ebf19635 54#include <Riostream.h>
b785c4bd 55
55fd51b0 56ClassImp(AliMUONTrackLight)
57
58//===================================================================
59
60AliMUONTrackLight::AliMUONTrackLight()
61 : TObject(),
62 fPrec(),
63 fIsTriggered(kFALSE),
64 fCharge(-999),
b88403f3 65 fChi2(-1),
55fd51b0 66 fCentr(-1),
67 fPgen(),
68 fTrackPythiaLine(-999),
69 fTrackPDGCode(-999),
70 fOscillation(kFALSE),
b88403f3 71 fNParents(0),
72 fWeight(1)
55fd51b0 73{
74 /// default constructor
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//============================================
89AliMUONTrackLight::AliMUONTrackLight(const AliMUONTrackLight &muonCopy)
90 : TObject(muonCopy),
91 fPrec(muonCopy.fPrec),
92 fIsTriggered(muonCopy.fIsTriggered),
93 fCharge(muonCopy.fCharge),
b88403f3 94 fChi2(muonCopy.fChi2),
55fd51b0 95 fCentr(muonCopy.fCentr),
96 fPgen(muonCopy.fPgen),
97 fTrackPythiaLine(muonCopy.fTrackPythiaLine),
98 fTrackPDGCode(muonCopy.fTrackPDGCode),
99 fOscillation(muonCopy.fOscillation),
b88403f3 100 fNParents(muonCopy.fNParents),
101 fWeight(muonCopy.fWeight)
55fd51b0 102{
103 /// copy constructor
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//============================================
116AliMUONTrackLight::AliMUONTrackLight(AliESDMuonTrack* muonTrack)
117 : TObject(),
118 fPrec(),
119 fIsTriggered(kFALSE),
120 fCharge(-999),
b88403f3 121 fChi2(-1),
55fd51b0 122 fCentr(-1),
123 fPgen(),
124 fTrackPythiaLine(-999),
125 fTrackPDGCode(-999),
126 fOscillation(kFALSE),
b88403f3 127 fNParents(0),
128 fWeight(1)
55fd51b0 129{
130 /// constructor
f202486b 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);
55fd51b0 141}
142
143//============================================
b88403f3 144AliMUONTrackLight::~AliMUONTrackLight()
145{
146/// Destructor
147}
148
149//============================================
e733e3dd 150AliMUONTrackLight& AliMUONTrackLight::operator=(const AliMUONTrackLight& muonCopy)
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;
165 fTrackPythiaLine = muonCopy.fTrackPythiaLine;
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 }
72421d27 180
181 return *this;
e733e3dd 182}
183
184//============================================
b88403f3 185
186void AliMUONTrackLight::FillFromAliMUONTrack(AliMUONTrack *trackReco,Double_t zvert){
187 /// this method sets the muon reconstructed momentum according to the value given by AliMUONTrack
f202486b 188 AliMUONTrackParam* trPar = trackReco->GetTrackParamAtVertex();
189 if (!trPar) {
190 AliError("The track must contain the parameters at vertex");
191 return;
61fed964 192 }
f202486b 193 this->SetCharge(Int_t(TMath::Sign(1.,trPar->GetInverseBendingMomentum())));
194 this->SetPxPyPz(trPar->Px(),trPar->Py(), trPar->Pz());
b88403f3 195 this->SetTriggered(trackReco->GetMatchTrigger());
196
f202486b 197 Double_t xyz[3] = { trPar->GetNonBendingCoor(),
198 trPar->GetBendingCoor(),
199 trPar->GetZ()};
61fed964 200 if (zvert!=-9999) xyz[2] = zvert;
b88403f3 201 this->SetVertex(xyz);
202}
203
204//============================================
205void AliMUONTrackLight::FillFromESD(AliESDMuonTrack* muonTrack,Double_t zvert){
55fd51b0 206 /// computes prec and charge from ESD track
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);
b88403f3 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();
55fd51b0 227}
228
229//============================================
230void AliMUONTrackLight::SetPxPyPz(Double_t px, Double_t py, Double_t pz){
231 /// set the reconstructed 4-momentum, assuming the particle is a muon
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//============================================
b88403f3 238void AliMUONTrackLight::FillMuonHistory(AliStack *stack, TParticle *part){
239 /// scans the muon history to determine parents pdg code and pythia line
55fd51b0 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
55fd51b0 244 TParticle *mother;
245 Int_t status=-1, pdg=-1;
246 while(lineM >= 0){
b88403f3 247
55fd51b0 248 mother = stack->Particle(lineM); //direct mother of rec. track
249 pdg = mother->GetPdgCode();//store PDG code of first mother
b88403f3 250 // break if a string, gluon, quark or diquark is found
251 if(pdg == 92 || pdg == 21 || TMath::Abs(pdg) < 10 || IsDiquark(pdg)) break;
55fd51b0 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 }
55fd51b0 263 fNParents = countP+1;
264 countP = -1;
b88403f3 265
55fd51b0 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 }
b88403f3 275
55fd51b0 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;
b88403f3 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 }
55fd51b0 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
9bdda5f6 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
7cdc832b 292 pdg = this->GetQuarkPDGCode(countP);
293 Int_t line = this->GetQuarkPythiaLine(countP);
55fd51b0 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 }
b88403f3 309 this->PrintInfo("h");
310 }//mismatch
55fd51b0 311 }
312}
b88403f3 313
55fd51b0 314//====================================
315void AliMUONTrackLight::ResetQuarkInfo(){
316 /// resets parton information
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//====================================
324Bool_t AliMUONTrackLight::IsB0(Int_t intTest) const {
325 /// checks if the particle is a B0
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//====================================
337Bool_t AliMUONTrackLight::IsMotherAResonance(Int_t index) const {
338 /// checks if mother is a resonance
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//====================================
346Int_t AliMUONTrackLight::GetParentFlavour(Int_t idParent) const {
347 /// returns the flavour of parent idParent (idParent=0 is the oldest
348 /// hadronized parent)
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//====================================
57e2ad1a 356void AliMUONTrackLight::PrintInfo(const Option_t* opt){
55fd51b0 357 /// prints information about the track:
358 /// - "H" muon's decay history
359 /// - "K" muon kinematics
360 /// - "A" all variables
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){
b785c4bd 369 snprintf(name, 100, "%4d --> ", this->GetQuarkPythiaLine(i));
55fd51b0 370 line += name;
b785c4bd 371 snprintf(name, 100, "%4d --> ", this->GetQuarkPDGCode(i));
55fd51b0 372 pdg += name;
373 }
374 }
375 for(int i = 0; i < fNParents; i++){
376 if(this->GetParentPythiaLine(i)>= 0){
b785c4bd 377 snprintf(name, 100, "%7d --> ", this->GetParentPythiaLine(i));
55fd51b0 378 line += name;
b785c4bd 379 snprintf(name, 100, "%7d --> ", this->GetParentPDGCode(i));
55fd51b0 380 pdg += name;
381 }
382 }
b785c4bd 383 snprintf(name, 100, "%4d", this->GetTrackPythiaLine()); line += name;
384 snprintf(name, 100, "%4d", this->GetTrackPDGCode()); pdg += name;
55fd51b0 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}
b88403f3 405//====================================
406Bool_t AliMUONTrackLight::IsParentPionOrKaon(Int_t idparent){
407 /// checks if a muon comes from a pion or kaon or a particle that decays into one of these two
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//====================================
9ee1d6ff 421Bool_t AliMUONTrackLight::IsDiquark(Int_t pdg) const{
4adf1b4c 422 /// check if the provided pdg code corresponds to a diquark
b88403f3 423 pdg = TMath::Abs(pdg);
424 if((pdg > 1000) && (pdg%100 < 10)) return kTRUE;
425 else return kFALSE;
426}
4adf1b4c 427//====================================
428void AliMUONTrackLight::SetParentPDGCode(Int_t index, Int_t pdg) {
429 /// Set hadronised parents and grandparents
430 if ( index < fgkNParentsMax ) {
431 fParentPDGCode[index] = pdg;
432 } else {
433 AliErrorStream() << "Index outside the array size." << std::endl;
434 }
435}
436//====================================
437void AliMUONTrackLight::SetParentPythiaLine(Int_t index, Int_t line) {
438 /// Set line of Pythia output for hadronised parents & grandparents
439 if ( index < fgkNParentsMax ) {
440 fParentPythiaLine[index] = line;
441 } else {
442 AliErrorStream() << "Index outside the array size." << std::endl;
443 }
444}
445//====================================
446void AliMUONTrackLight::SetQuarkPDGCode(Int_t index, Int_t pdg){
447 /// Set pdg of the string [0], quarks/gluons [1,2], sometimes proton [3]
448 if ( index < 4 ) {
449 fQuarkPDGCode[index] = pdg;
450 } else {
451 AliErrorStream() << "Index outside the array size." << std::endl;
452 }
453}
454//====================================
455void AliMUONTrackLight::SetQuarkPythiaLine(Int_t index, Int_t line){
456 /// Set line of Pythia output for string [0] and quarks [1,2], sometimes proton [3]
457 if ( index < 4 ) {
458 fQuarkPythiaLine[index] = line;
459 } else {
460 AliErrorStream() << "Index outside the array size." << std::endl;
461 }
462}