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