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