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