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