Include Riostream.h
[u/mrichter/AliRoot.git] / MUON / MUONevaluation / 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 #include <cstdio>
54 #include <Riostream.h>
55
56 ClassImp(AliMUONTrackLight) 
57
58 //===================================================================
59
60 AliMUONTrackLight::AliMUONTrackLight() 
61   : TObject(), 
62     fPrec(), 
63     fIsTriggered(kFALSE),
64     fCharge(-999), 
65     fChi2(-1), 
66     fCentr(-1),
67     fPgen(), 
68     fTrackPythiaLine(-999),
69     fTrackPDGCode(-999),
70     fOscillation(kFALSE), 
71     fNParents(0),
72     fWeight(1)    
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 //============================================
89 AliMUONTrackLight::AliMUONTrackLight(const AliMUONTrackLight &muonCopy) 
90   : TObject(muonCopy), 
91     fPrec(muonCopy.fPrec), 
92     fIsTriggered(muonCopy.fIsTriggered),
93     fCharge(muonCopy.fCharge), 
94     fChi2(muonCopy.fChi2), 
95     fCentr(muonCopy.fCentr),
96     fPgen(muonCopy.fPgen), 
97     fTrackPythiaLine(muonCopy.fTrackPythiaLine),
98     fTrackPDGCode(muonCopy.fTrackPDGCode),
99     fOscillation(muonCopy.fOscillation), 
100     fNParents(muonCopy.fNParents),
101     fWeight(muonCopy.fWeight)
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 //============================================
116 AliMUONTrackLight::AliMUONTrackLight(AliESDMuonTrack* muonTrack)
117   : TObject(), 
118     fPrec(), 
119     fIsTriggered(kFALSE),
120     fCharge(-999), 
121     fChi2(-1),
122     fCentr(-1),
123     fPgen(), 
124     fTrackPythiaLine(-999),
125     fTrackPDGCode(-999),
126     fOscillation(kFALSE), 
127     fNParents(0),
128     fWeight(1)
129
130   /// constructor
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);
141 }
142
143 //============================================
144 AliMUONTrackLight::~AliMUONTrackLight()
145 {
146 /// Destructor
147
148
149 //============================================
150 AliMUONTrackLight& 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   }
180
181   return *this;
182 }    
183
184 //============================================
185
186 void AliMUONTrackLight::FillFromAliMUONTrack(AliMUONTrack *trackReco,Double_t zvert){
187   /// this method sets the muon reconstructed momentum according to the value given by AliMUONTrack
188   AliMUONTrackParam* trPar = trackReco->GetTrackParamAtVertex();
189   if (!trPar) {
190     AliError("The track must contain the parameters at vertex");
191     return;
192   }
193   this->SetCharge(Int_t(TMath::Sign(1.,trPar->GetInverseBendingMomentum())));
194   this->SetPxPyPz(trPar->Px(),trPar->Py(), trPar->Pz()); 
195   this->SetTriggered(trackReco->GetMatchTrigger()); 
196   
197   Double_t xyz[3] = { trPar->GetNonBendingCoor(), 
198                       trPar->GetBendingCoor(),
199                       trPar->GetZ()};
200   if (zvert!=-9999) xyz[2] = zvert;
201   this->SetVertex(xyz); 
202 }
203
204 //============================================
205 void AliMUONTrackLight::FillFromESD(AliESDMuonTrack* muonTrack,Double_t zvert){
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);
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();
227 }
228
229 //============================================
230 void 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 //============================================
238 void AliMUONTrackLight::FillMuonHistory(AliStack *stack, TParticle *part){
239   /// scans the muon history to determine parents pdg code and pythia line
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
244   TParticle *mother;
245   Int_t status=-1, pdg=-1;
246   while(lineM >= 0){
247     
248     mother = stack->Particle(lineM); //direct mother of rec. track
249     pdg = mother->GetPdgCode();//store PDG code of first mother
250     // break if a string, gluon, quark or diquark is found 
251     if(pdg == 92 || pdg == 21 || TMath::Abs(pdg) < 10 || IsDiquark(pdg)) break;
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   }
263   fNParents = countP+1;
264   countP = -1;
265
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   }
275   
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;
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   }
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
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       
292       pdg = this->GetQuarkPDGCode(countP);
293       Int_t line = this->GetQuarkPythiaLine(countP);
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       }
309       this->PrintInfo("h");
310     }//mismatch
311   }
312 }
313
314 //====================================
315 void 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 //====================================
324 Bool_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 //====================================
337 Bool_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 //====================================
346 Int_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 //====================================
356 void AliMUONTrackLight::PrintInfo(const Option_t* opt){
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){
369         snprintf(name, 100, "%4d --> ", this->GetQuarkPythiaLine(i));
370         line += name;
371         snprintf(name, 100, "%4d --> ", this->GetQuarkPDGCode(i));
372         pdg += name;
373       }
374     }
375     for(int i = 0; i < fNParents; i++){ 
376       if(this->GetParentPythiaLine(i)>= 0){
377         snprintf(name, 100, "%7d --> ", this->GetParentPythiaLine(i));
378         line += name;
379         snprintf(name, 100, "%7d --> ", this->GetParentPDGCode(i));
380         pdg += name;
381       }
382     }
383     snprintf(name, 100, "%4d", this->GetTrackPythiaLine()); line += name;
384     snprintf(name, 100, "%4d", this->GetTrackPDGCode()); pdg += name;
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 }
405 //====================================
406 Bool_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 //====================================
421 Bool_t AliMUONTrackLight::IsDiquark(Int_t pdg) const{
422   /// check if the provided pdg code corresponds to a diquark
423   pdg = TMath::Abs(pdg);
424   if((pdg > 1000) && (pdg%100 < 10)) return kTRUE;
425   else return kFALSE;
426 }
427 //====================================
428 void 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 //====================================
437 void 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 //====================================
446 void 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 //====================================
455 void 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 }