Coding conventions
[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 "AliMUONTrackExtrap.h"
44 #include "AliMUONTrackParam.h"
45
46 #include "AliESDMuonTrack.h"
47 #include "AliRunLoader.h"
48 #include "AliStack.h"
49 #include "AliHeader.h"
50
51 #include "TDatabasePDG.h"
52 #include "TParticle.h"
53 #include "TString.h"
54
55 ClassImp(AliMUONTrackLight) 
56
57 //===================================================================
58
59 AliMUONTrackLight::AliMUONTrackLight() 
60   : TObject(), 
61     fPrec(), 
62     fIsTriggered(kFALSE),
63     fCharge(-999), 
64     fChi2(-1), 
65     fCentr(-1),
66     fPgen(), 
67     fTrackPythiaLine(-999),
68     fTrackPDGCode(-999),
69     fOscillation(kFALSE), 
70     fNParents(0),
71     fWeight(1)    
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 //============================================
88 AliMUONTrackLight::AliMUONTrackLight(const AliMUONTrackLight &muonCopy) 
89   : TObject(muonCopy), 
90     fPrec(muonCopy.fPrec), 
91     fIsTriggered(muonCopy.fIsTriggered),
92     fCharge(muonCopy.fCharge), 
93     fChi2(muonCopy.fChi2), 
94     fCentr(muonCopy.fCentr),
95     fPgen(muonCopy.fPgen), 
96     fTrackPythiaLine(muonCopy.fTrackPythiaLine),
97     fTrackPDGCode(muonCopy.fTrackPDGCode),
98     fOscillation(muonCopy.fOscillation), 
99     fNParents(muonCopy.fNParents),
100     fWeight(muonCopy.fWeight)
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 //============================================
115 AliMUONTrackLight::AliMUONTrackLight(AliESDMuonTrack* muonTrack)
116   : TObject(), 
117     fPrec(), 
118     fIsTriggered(kFALSE),
119     fCharge(-999), 
120     fChi2(-1),
121     fCentr(-1),
122     fPgen(), 
123     fTrackPythiaLine(-999),
124     fTrackPDGCode(-999),
125     fOscillation(kFALSE), 
126     fNParents(0),
127     fWeight(1)
128
129   /// constructor
130   //AliMUONTrackLight(); 
131   FillFromESD(muonTrack); 
132 }
133
134 //============================================
135 AliMUONTrackLight::~AliMUONTrackLight()
136 {
137 /// Destructor
138
139
140 //============================================
141
142 void AliMUONTrackLight::FillFromAliMUONTrack(AliMUONTrack *trackReco,Double_t zvert){
143   /// this method sets the muon reconstructed momentum according to the value given by AliMUONTrack
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   }
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(),
156                       trPar.GetZ()};
157   if (zvert!=-9999) xyz[2] = zvert;
158   this->SetVertex(xyz); 
159 }
160
161 //============================================
162 void AliMUONTrackLight::FillFromESD(AliESDMuonTrack* muonTrack,Double_t zvert){
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);
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();
184 }
185
186 //============================================
187 void 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 //============================================
195 TParticle* AliMUONTrackLight::FindRefTrack(
196                 AliMUONTrack* trackReco, AliMUONVTrackStore* trackRefArray,
197                 AliStack* stack)
198 {
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.
205   
206   TParticle *part = 0; 
207   const Double_t kSigmaCut = 4;  // 4 sigmas cut
208   Int_t compPart = 0;
209   Bool_t compTrack;
210   
211   TIter next(trackRefArray->CreateIterator());
212   AliMUONTrack* trackRef;
213   while ( (trackRef = static_cast<AliMUONTrack*>(next())) ) {
214     // check if trackRef is compatible with trackReco:
215     compTrack = kFALSE;
216     
217     if (trackReco->GetNClusters() > 1) {
218       
219       // check cluster by cluster if trackReco contain info at each cluster
220       Bool_t compTrackArray[10];
221       trackRef->CompatibleTrack(trackReco,kSigmaCut,compTrackArray);
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) {
236       compPart++;
237       Int_t trackID = trackRef->GetTrackID();
238       this->SetTrackPythiaLine(trackID);
239       part = stack->Particle(trackID);
240       fTrackPDGCode = part->GetPdgCode();
241     }
242   }
243   if (compPart>1) {
244     AliFatal("More than one particle compatible with the reconstructed track.");
245   } 
246   return part;
247 }
248
249 //============================================
250 Int_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)
254   
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 //============================================
270 void AliMUONTrackLight::FillMuonHistory(AliStack *stack, TParticle *part){
271   /// scans the muon history to determine parents pdg code and pythia line
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
276   TParticle *mother;
277   Int_t status=-1, pdg=-1;
278   while(lineM >= 0){
279     
280     mother = stack->Particle(lineM); //direct mother of rec. track
281     pdg = mother->GetPdgCode();//store PDG code of first mother
282     // break if a string, gluon, quark or diquark is found 
283     if(pdg == 92 || pdg == 21 || TMath::Abs(pdg) < 10 || IsDiquark(pdg)) break;
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   }
295   fNParents = countP+1;
296   countP = -1;
297
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   }
307   
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;
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   }
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
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       
324       pdg = this->GetQuarkPDGCode(countP);
325       Int_t line = this->GetQuarkPythiaLine(countP);
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       }
341       this->PrintInfo("h");
342     }//mismatch
343   }
344 }
345
346 //====================================
347 void 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 //====================================
356 Bool_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 //====================================
369 Bool_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 //====================================
378 Int_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 //====================================
388 void 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 }
437 //====================================
438 Bool_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 //====================================
453 Bool_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 }