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