Removing the AliClusterTGeo class
[u/mrichter/AliRoot.git] / MUON / AliMUONTrackLight.cxx
CommitLineData
55fd51b0 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//This class was prepared by INFN Cagliari, July 2006
20//(authors: H.Woehri, A.de Falco)
21//
22//
23// Compact information for the muon generated tracks in the MUON arm
24// useful at the last stage of the analysis chain
25// provides a link between the reconstructed track and the generated particle
26// stores kinematical information at gen. and rec. level and
27// the decay history of the muon, allowing the identification of the
28// mother process
29//
30// To be used together with AliMUONPairLight
31//===================================================================
32
33#include "AliMUONTrackLight.h"
34#include "AliMUONTrack.h"
35#include "AliMUONConstants.h"
36
37#include "AliESDMuonTrack.h"
38#include "AliRunLoader.h"
39#include "AliStack.h"
40#include "AliHeader.h"
41
42#include "TDatabasePDG.h"
43#include "TClonesArray.h"
44#include "TParticle.h"
45#include "TString.h"
46
47ClassImp(AliMUONTrackLight)
48
49//===================================================================
50
51AliMUONTrackLight::AliMUONTrackLight()
52 : TObject(),
53 fPrec(),
54 fIsTriggered(kFALSE),
55 fCharge(-999),
56 fCentr(-1),
57 fPgen(),
58 fTrackPythiaLine(-999),
59 fTrackPDGCode(-999),
60 fOscillation(kFALSE),
61 fNParents(0)
62{
63 /// default constructor
64 fPgen.SetPxPyPzE(0.,0.,0.,0.);
65 fPrec.SetPxPyPzE(0.,0.,0.,0.);
66 for (Int_t i=0; i<3; i++) fXYZ[i]=-999;
67 for (Int_t npar = 0; npar < fgkNParentsMax; npar++){
68 fParentPDGCode[npar] = -1;
69 fParentPythiaLine[npar] = -1;
70 }
71 for (Int_t i = 0; i < 4; i++){
72 fQuarkPDGCode[i] = -1;
73 fQuarkPythiaLine[i] = -1;
74 }
75}
76
77//============================================
78AliMUONTrackLight::AliMUONTrackLight(const AliMUONTrackLight &muonCopy)
79 : TObject(muonCopy),
80 fPrec(muonCopy.fPrec),
81 fIsTriggered(muonCopy.fIsTriggered),
82 fCharge(muonCopy.fCharge),
83 fCentr(muonCopy.fCentr),
84 fPgen(muonCopy.fPgen),
85 fTrackPythiaLine(muonCopy.fTrackPythiaLine),
86 fTrackPDGCode(muonCopy.fTrackPDGCode),
87 fOscillation(muonCopy.fOscillation),
88 fNParents(muonCopy.fNParents)
89{
90 /// copy constructor
91 for (Int_t i=0; i<3; i++) fXYZ[i]=muonCopy.fXYZ[i];
92 for (Int_t npar = 0; npar < fgkNParentsMax; npar++){
93 fParentPDGCode[npar] = muonCopy.fParentPDGCode[npar];
94 fParentPythiaLine[npar] = muonCopy.fParentPythiaLine[npar];
95 }
96 for (Int_t i = 0; i < 4; i++){
97 fQuarkPDGCode[i] = muonCopy.fQuarkPDGCode[i];
98 fQuarkPythiaLine[i] = muonCopy.fQuarkPythiaLine[i];
99 }
100}
101
102//============================================
103AliMUONTrackLight::AliMUONTrackLight(AliESDMuonTrack* muonTrack)
104 : TObject(),
105 fPrec(),
106 fIsTriggered(kFALSE),
107 fCharge(-999),
108 fCentr(-1),
109 fPgen(),
110 fTrackPythiaLine(-999),
111 fTrackPDGCode(-999),
112 fOscillation(kFALSE),
113 fNParents(0)
114{
115 /// constructor
116 //AliMUONTrackLight();
117 ComputePRecAndChargeFromESD(muonTrack);
118}
119
120//============================================
121void AliMUONTrackLight::ComputePRecAndChargeFromESD(AliESDMuonTrack* muonTrack){
122 /// computes prec and charge from ESD track
123 Double_t mumass = TDatabasePDG::Instance()->GetParticle(13)->Mass();
124 Double_t thetaX = muonTrack->GetThetaX();
125 Double_t thetaY = muonTrack->GetThetaY();
126 Double_t tanthx = TMath::Tan(thetaX);
127 Double_t tanthy = TMath::Tan(thetaY);
128 Double_t pYZ = 1./TMath::Abs(muonTrack->GetInverseBendingMomentum());
129 Double_t pz = - pYZ / TMath::Sqrt(1.0 + tanthy * tanthy);
130 Double_t px = pz * tanthx;
131 Double_t py = pz * tanthy;
132 fCharge = Int_t(TMath::Sign(1.,muonTrack->GetInverseBendingMomentum()));
133 Double_t energy = TMath::Sqrt(mumass * mumass + px*px + py*py + pz*pz);
134 fPrec.SetPxPyPzE(px,py,pz,energy);
135}
136
137//============================================
138void AliMUONTrackLight::SetPxPyPz(Double_t px, Double_t py, Double_t pz){
139 /// set the reconstructed 4-momentum, assuming the particle is a muon
140 Double_t mumass = TDatabasePDG::Instance()->GetParticle(13)->Mass();
141 Double_t energy = TMath::Sqrt(mumass * mumass + px*px + py*py + pz*pz);
142 fPrec.SetPxPyPzE(px,py,pz,energy);
143}
144
145//============================================
146TParticle* AliMUONTrackLight::FindRefTrack(AliMUONTrack* trackReco, TClonesArray* trackRefArray, AliRunLoader *runLoader){
147 /// find the MC particle that corresponds to a given rec track
148 TParticle *part = 0;
149 const Double_t kSigma2Cut = 16; // 4 sigmas cut, kSigma2Cut = 4*4
150 Int_t nTrackRef = trackRefArray->GetEntriesFast();
151 Int_t compPart = 0;
152 for (Int_t iref = 0; iref < nTrackRef; iref++) {
153 AliMUONTrack *trackRef = (AliMUONTrack *)trackRefArray->At(iref);
154 // check if trackRef is compatible with trackReco:
155 //routine returns for each chamber a yes/no information if the
156 //hit of rec. track and hit of referenced track are compatible
157 Bool_t *compTrack = trackRef->CompatibleTrack(trackReco,kSigma2Cut);
158 Int_t iTrack = this->TrackCheck(compTrack); //returns number of validated conditions
159 if (iTrack==4) {
160 compPart++;
161 Int_t trackID = trackRef->GetTrackID();
162 this->SetTrackPythiaLine(trackID);
163 part = ((AliStack *)(((AliHeader *) runLoader->GetHeader())->Stack()))->Particle(trackID);
164 fTrackPDGCode = part->GetPdgCode();
165 }
166 }
167 if (compPart>1) {
168 printf ("<AliMUONTrackLight::FindRefTrack> ERROR: more than one particle compatible to the reconstructed track.\n");
169 Int_t i=0, j=1/i;
170 printf ("j=%d \n",j);
171 }
172 return part;
173}
174
175//============================================
176Int_t AliMUONTrackLight::TrackCheck(Bool_t *compTrack){
177 /// Apply reconstruction requirements
178 /// Return number of validated conditions
179 /// If all the tests are verified then TrackCheck = 4 (good track)
180 Int_t iTrack = 0;
181 Int_t hitsInLastStations = 0;
182
183 // apply reconstruction requirements
184 if (compTrack[0] || compTrack[1]) iTrack++; // at least one hit in st. 0
185 if (compTrack[2] || compTrack[3]) iTrack++; // at least one hit in st. 1
186 if (compTrack[4] || compTrack[5]) iTrack++; // at least one hit in st. 2
187 for (Int_t ch = 6; ch < AliMUONConstants::NTrackingCh(); ch++) {
188 if (compTrack[ch]) hitsInLastStations++;
189 }
190 if (hitsInLastStations > 2) iTrack++; // at least 3 hits in st. 3 & 4
191 return iTrack;
192}
193
194//============================================
195void AliMUONTrackLight::FillMuonHistory(AliRunLoader *runLoader, TParticle *part){
196 /// scans the muon history to determine parents pdg code and pythia line
197 Int_t countP = -1;
198 Int_t parents[10], parLine[10];
199 Int_t lineM = part->GetFirstMother();//line in the Pythia output of the particle's mother
200
201 AliStack *stack = runLoader->GetHeader()->Stack();
202 TParticle *mother;
203 Int_t status=-1, pdg=-1;
204 while(lineM >= 0){
205
206 mother = stack->Particle(lineM); //direct mother of rec. track
207 pdg = mother->GetPdgCode();//store PDG code of first mother
208 if(pdg == 92) break; // break if a string is found
209 parents[++countP] = pdg;
210 parLine[countP] = lineM;
211 status = mother->GetStatusCode();//get its status code to check if oscillation occured
212 if(IsB0(parents[countP]) && status == 12) this->SetOscillation(kTRUE);
213 lineM = mother->GetFirstMother();
214 }
215 //store all the fragmented parents in an array:
216 for(int i = 0; i <= countP; i++){
217 this->SetParentPDGCode(i,parents[countP-i]);
218 this->SetParentPythiaLine(i,parLine[countP-i]);
219 }
220
221 fNParents = countP+1;
222 countP = -1;
223 //and store the lines of the string and further quarks in another array:
224 while(lineM >= 0){
225 mother = stack->Particle(lineM);
226 pdg = mother->GetPdgCode();
227 //now, get information before the fragmentation
228 this->SetQuarkPythiaLine(++countP, lineM);//store the line of the string in index 0
229 this->SetQuarkPDGCode(countP, pdg);//store the pdg of the quarks in index 1,2
230 lineM = mother->GetFirstMother();
231 }
232
233 //check if in case of HF production, the string points to the correct end
234 //and correct it in case of need:
235 countP = 1;
236 if(this->GetQuarkPythiaLine(countP) > -1 && (this->GetParentFlavour(0)==4 || this->GetParentFlavour(0)==5)){
237 if(this->GetParentFlavour(0) != TMath::Abs(this->GetQuarkPDGCode(countP))){
238
239 Int_t pdg = this->GetQuarkPDGCode(1), line = this->GetQuarkPythiaLine(1);
240 this->ResetQuarkInfo();
241 while(TMath::Abs(pdg) != this->GetParentFlavour(0)){//pdg of q,g in Pythia listing following the wrong string end
242 //must coincide with the flavour of the last fragmented mother
243
244 pdg = stack->Particle(++line)->GetPdgCode();
245 }
246 //now, we have the correct string end and the correct line
247 //continue to fill again all parents and corresponding lines
248 while(line >= 0){
249 mother = stack->Particle(line);//get again the mother
250 pdg = mother->GetPdgCode();
251 this->SetQuarkPythiaLine(countP, line);
252 this->SetQuarkPDGCode(countP++, pdg);
253 line = mother->GetFirstMother();
254 }
255 }
256 }
257}
258//====================================
259void AliMUONTrackLight::ResetQuarkInfo(){
260 /// resets parton information
261 for(int pos = 1; pos < 4; pos++){//[0] is the string
262 this->SetQuarkPDGCode(pos,-1);
263 this->SetQuarkPythiaLine(pos,-1);
264 }
265}
266
267//====================================
268Bool_t AliMUONTrackLight::IsB0(Int_t intTest) const {
269 /// checks if the particle is a B0
270 Int_t bMes0[2] = {511,531};//flavour code of B0d and B0s
271 Bool_t answer = kFALSE;
272 for(int i = 0; i < 2; i++){
273 if(TMath::Abs(intTest) == bMes0[i]){
274 answer = kTRUE;
275 break;
276 }
277 }
278 return answer;
279}
280//====================================
281Bool_t AliMUONTrackLight::IsMotherAResonance(Int_t index) const {
282 /// checks if mother is a resonance
283 Int_t intTest = GetParentPDGCode(index);
284 // the resonance pdg code is built this way
285 // x00ffn where x=0,1,.. (1S,2S... states), f=quark flavour
286 Int_t id=intTest%100000;
287 return (!((id-id%10)%110));
288}
289//====================================
290Int_t AliMUONTrackLight::GetParentFlavour(Int_t idParent) const {
291 /// returns the flavour of parent idParent (idParent=0 is the oldest
292 /// hadronized parent)
293 Int_t pdg = GetParentPDGCode(idParent);
294 Int_t quark = TMath::Abs(pdg/100);
295 if(quark > 9) quark = quark/10;
296 return quark;
297}
298
299//====================================
300void AliMUONTrackLight::PrintInfo(Option_t* opt){
301 /// prints information about the track:
302 /// - "H" muon's decay history
303 /// - "K" muon kinematics
304 /// - "A" all variables
305 TString options(opt);
306 options.ToUpper();
307
308 if(options.Contains("H") || options.Contains("A")){ //muon decay history
309 char *name= new char[100];
310 TString pdg = "", line = "";
311 for(int i = 3; i >= 0; i--){
312 if(this->GetQuarkPythiaLine(i)>= 0){
313 sprintf(name, "%4d --> ", this->GetQuarkPythiaLine(i));
314 line += name;
315 sprintf(name, "%4d --> ", this->GetQuarkPDGCode(i));
316 pdg += name;
317 }
318 }
319 for(int i = 0; i < fNParents; i++){
320 if(this->GetParentPythiaLine(i)>= 0){
321 sprintf(name, "%7d --> ", this->GetParentPythiaLine(i));
322 line += name;
323 sprintf(name, "%7d --> ", this->GetParentPDGCode(i));
324 pdg += name;
325 }
326 }
327 sprintf(name, "%4d", this->GetTrackPythiaLine()); line += name;
328 sprintf(name, "%4d", this->GetTrackPDGCode()); pdg += name;
329
330 printf("\nmuon's decay history:\n");
331 printf(" PDG: %s\n", pdg.Data());
332 printf("line: %s\n", line.Data());
333 }
334 if(options.Contains("K") || options.Contains("A")){ //muon kinematic
335
336 Int_t charge = this->GetCharge();
337 Double_t *vtx = this->GetVertex();
338 TLorentzVector momRec = this->GetPRec();
339 TLorentzVector momGen = this->GetPGen();
340 printf("the track's charge is %d\n", charge);
341 printf("Primary vertex: Vx = %1.3f, Vy = %1.3f, Vz = %1.3f\n", vtx[0], vtx[1], vtx[2]);
342 printf("Generated: Px = %1.3f, Py = %1.3f, Pz = %1.3f\n", momGen.Px(), momGen.Py(), momGen.Pz());
343 printf("Reconstructed: Px = %1.3f, Py = %1.3f, Pz = %1.3f\n", momRec.Px(), momRec.Py(), momRec.Pz());
344 printf("Rec. variables: pT %1.3f, pseudo-rapidity %1.3f, theta %1.3f (%1.3f degree), phi %1.3f (%1.3f degree)\n",
345 momRec.Pt(), momRec.Eta(), momRec.Theta(), 180./TMath::Pi() * momRec.Theta(),
346 momRec.Phi(), 180./TMath::Pi() * momRec.Phi());
347 }
348}
349
350