AliAODEvent::GetHeader now return AliVHeader
[u/mrichter/AliRoot.git] / PWG / muon / AliUtilityMuonAncestor.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-2007, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
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 purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15
16 /* $Id: AliUtilityMuonAncestor.cxx 47782 2011-02-24 18:37:31Z martinez $ */
17
18 //-----------------------------------------------------------------------------
19 /// \class AliUtilityMuonAncestor
20 /// Static utilities to get the muon ancestor in MC
21 ///
22 /// \author Diego Stocco
23 //-----------------------------------------------------------------------------
24
25 #include "AliUtilityMuonAncestor.h"
26
27 // ROOT includes
28 #include "TDatabasePDG.h"
29 #include "TParticlePDG.h"
30 #include "TPDGCode.h"
31 #include "TMCProcess.h"
32 #include "TMath.h"
33
34 // STEER includes
35 #include "AliMCEvent.h"
36 #include "AliVParticle.h"
37 #include "AliMCParticle.h"
38 #include "AliAODMCParticle.h"
39 #include "AliLog.h"
40
41 // PWGmuon includes
42 #include "AliAnalysisMuonUtility.h"
43
44 /// \cond CLASSIMP
45 ClassImp(AliUtilityMuonAncestor) // Class implementation in ROOT context
46 /// \endcond
47
48 //_________________________________________________________
49 AliUtilityMuonAncestor::AliUtilityMuonAncestor() :
50 TObject(),
51 fPx(0.),
52 fPy(0.),
53 fPz(0.),
54 fMask(0),
55 fAncestor(-999)
56 {
57   /// Default constructor
58 }
59
60
61 //_________________________________________________________
62 AliUtilityMuonAncestor::~AliUtilityMuonAncestor()
63 {
64   /// Default destructor
65 }
66
67 //_________________________________________________________
68 AliUtilityMuonAncestor::AliUtilityMuonAncestor(const AliUtilityMuonAncestor& obj) :
69 TObject(obj),
70 fPx(obj.fPx),
71 fPy(obj.fPy),
72 fPz(obj.fPz),
73 fMask(obj.fMask),
74 fAncestor(obj.fAncestor)
75 {
76   /// Copy constructor
77 }
78
79 //_________________________________________________________
80 AliUtilityMuonAncestor& AliUtilityMuonAncestor::operator=(const AliUtilityMuonAncestor& obj)
81 {
82   /// Copy operator
83   if ( this != &obj ) {
84     TObject::operator=(obj);
85     fPx = obj.fPx;
86     fPy = obj.fPy;
87     fPz = obj.fPz;
88     fMask = obj.fMask;
89     fAncestor = obj.fAncestor;
90   }
91   return *this;
92 }
93
94 //_________________________________________________________
95 Bool_t AliUtilityMuonAncestor::BuildAncestor ( AliVParticle* track, const AliMCEvent* mcEvent )
96 {
97   /// Build ancestor
98   
99   // If track is the same as the one in memory, do not re-build the track ancestor
100   if ( track->Px() == fPx && track->Py() == fPy && track->Pz() == fPz ) return kTRUE;
101   fPx = track->Px();
102   fPy = track->Py();
103   fPz = track->Pz();
104   fMask = 0;
105   fAncestor = -999;
106   if ( ! mcEvent ) return kFALSE;
107   
108   AliVParticle* mcParticle = 0x0;
109   
110   if ( AliAnalysisMuonUtility::IsMCTrack(track) ) mcParticle = track;
111   else {
112     Int_t trackLabel = track->GetLabel();
113     if ( trackLabel < 0 ) return kFALSE;
114     mcParticle = mcEvent->GetTrack(trackLabel);
115   }
116   
117   // Track is MC (or matches a MC track)
118   SETBIT(fMask,kIsID);
119   
120   Int_t recoPdg = mcParticle->PdgCode();
121   
122   // Track is not a muon
123   if ( TMath::Abs(recoPdg) == 13 ) SETBIT(fMask,kIsMuon);
124   
125   Int_t imother = AliAnalysisMuonUtility::GetMotherIndex(mcParticle);
126
127   while ( imother >= 0 ) {
128     AliVParticle* part = mcEvent->GetTrack(imother);
129     
130     Int_t absPdg = TMath::Abs(part->PdgCode());
131     
132     // This is a quark
133     if ( absPdg < 10 ) return kTRUE;
134     
135     fAncestor = imother;
136     Bool_t isPrimary = AliAnalysisMuonUtility::IsPrimary(part, mcEvent);
137     
138     if ( isPrimary ) {
139       Int_t mpdg = absPdg%100000;
140       if ( mpdg >= 100 && mpdg < 10000 ) {
141         Int_t flv  = Int_t ( mpdg / TMath::Power(10, Int_t(TMath::Log10(mpdg) )));
142         if ( flv < 4 ) SETBIT(fMask,kHasLightParent);
143         else if ( flv >= 6 ) continue;
144         else {
145           TParticlePDG* partPdg = TDatabasePDG::Instance()->GetParticle(part->PdgCode());
146           if ( partPdg && ! partPdg->AntiParticle() ) SETBIT(fMask,kHasQuarkoniumParent);
147           else if ( flv == 4 ) SETBIT(fMask,kHasCharmParent);
148           else SETBIT(fMask,kHasBeautyParent);
149         }
150       } // absPdg within 100 and 10000
151     } // is primary
152     else {
153       UInt_t mcProcess = AliAnalysisMuonUtility::GetMCProcess(part);
154       if ( mcProcess == kPHadronic ||
155            ( mcProcess == 0 && part->Zv() < -90. ) ) {
156         // The MC process is not well computed in the AODs of old MC productions
157         // In this case, declare the particle as "secondary" if it is produced inside the front absorber
158         SETBIT(fMask,kIsSecondary);
159 //        return kTRUE;
160       }
161     } // is secondary
162     
163     imother = AliAnalysisMuonUtility::GetMotherIndex(part);
164     
165   } // loop on mothers
166   return kTRUE;
167 }
168
169 //_________________________________________________________
170 Bool_t AliUtilityMuonAncestor::CheckAncestor ( AliVParticle* track, const AliMCEvent* mcEvent, Int_t ancestorPdg, Bool_t matchAbsPdg )
171 {
172   /// Check ancestor
173   Int_t pdgCode = GetAncestorPdg(track, mcEvent);
174   if ( matchAbsPdg ) {
175     pdgCode = TMath::Abs(pdgCode);
176     ancestorPdg = TMath::Abs(ancestorPdg);
177   }
178   return ( pdgCode == ancestorPdg );
179 }
180
181 //_________________________________________________________
182 Int_t AliUtilityMuonAncestor::GetAncestor ( AliVParticle* track, const AliMCEvent* mcEvent )
183 {
184   /// Return ancestor (compute it if necessary)
185   BuildAncestor(track,mcEvent);
186   return fAncestor;
187 }
188
189 //_________________________________________________________
190 Int_t AliUtilityMuonAncestor::GetAncestorPdg ( AliVParticle* track, const AliMCEvent* mcEvent )
191 {
192   /// Return ancestor Pdg
193   if ( BuildAncestor(track,mcEvent) ) {
194     if ( fAncestor >= 0 ) return mcEvent->GetTrack(fAncestor)->PdgCode();
195   }
196   return 0;
197 }
198
199 //_________________________________________________________
200 Long64_t AliUtilityMuonAncestor::GetMask ( AliVParticle* track, const AliMCEvent* mcEvent )
201 {
202   /// Return mask
203   BuildAncestor(track,mcEvent);
204   return fMask;
205 }
206
207
208 //_________________________________________________________
209 Bool_t AliUtilityMuonAncestor::IsBeautyMu ( AliVParticle* track, const AliMCEvent* mcEvent )
210 {
211   /// Muon from beauty decays
212   Long64_t mask = GetMask(track,mcEvent);
213   return ( TESTBIT(mask,kIsID) && TESTBIT(mask,kIsMuon) && TESTBIT(mask,kHasBeautyParent) & ! TESTBIT(mask,kHasLightParent) );
214 }
215
216 //_________________________________________________________
217 Bool_t AliUtilityMuonAncestor::IsBJpsiMu ( AliVParticle* track, const AliMCEvent* mcEvent )
218 {
219   /// Muon B->J/psi decays
220   if ( IsBeautyMu(track,mcEvent) ) {
221     Int_t imother = AliAnalysisMuonUtility::GetMotherIndex(track);
222     if ( imother >= 0 ) return ( mcEvent->GetTrack(imother)->PdgCode() == 443 );
223   }
224   return kFALSE;
225 }
226
227 //_________________________________________________________
228 Bool_t AliUtilityMuonAncestor::IsCharmMu ( AliVParticle* track, const AliMCEvent* mcEvent )
229 {
230   /// Muon from charm decays
231   Long64_t mask = GetMask(track,mcEvent);
232   return ( TESTBIT(mask,kIsID) && TESTBIT(mask,kIsMuon) && TESTBIT(mask,kHasCharmParent) && ! TESTBIT(mask,kHasBeautyParent) && ! TESTBIT(mask,kHasLightParent) );
233 }
234
235 //_________________________________________________________
236 Bool_t AliUtilityMuonAncestor::IsDecayMu ( AliVParticle* track, const AliMCEvent* mcEvent )
237 {
238   /// Muon from light hadron decays
239   Long64_t mask = GetMask(track,mcEvent);
240   return ( TESTBIT(mask,kIsID) && TESTBIT(mask,kIsMuon) && TESTBIT(mask,kHasLightParent) && ! TESTBIT(mask,kIsSecondary) );
241 }
242
243 //_________________________________________________________
244 Bool_t AliUtilityMuonAncestor::IsHadron ( AliVParticle* track, const AliMCEvent* mcEvent )
245 {
246   /// Reconstructed hadron
247   Long64_t mask = GetMask(track,mcEvent);
248   return ( TESTBIT(mask,kIsID) && ! TESTBIT(mask,kIsMuon) );
249 }
250
251 //_________________________________________________________
252 Bool_t AliUtilityMuonAncestor::IsMuon ( AliVParticle* track, const AliMCEvent* mcEvent )
253 {
254   /// Track is muon
255   Long64_t mask = GetMask(track,mcEvent);
256   return ( TESTBIT(mask,kIsID) && TESTBIT(mask,kIsMuon) );
257 }
258
259 //_________________________________________________________
260 Bool_t AliUtilityMuonAncestor::IsQuarkoniumMu ( AliVParticle* track, const AliMCEvent* mcEvent )
261 {
262   /// Mu from quarkonium decay
263   Long64_t mask = GetMask(track,mcEvent);
264   return ( TESTBIT(mask,kIsID) && TESTBIT(mask,kIsMuon) && TESTBIT(mask,kHasQuarkoniumParent) );
265 }
266
267 //_________________________________________________________
268 Bool_t AliUtilityMuonAncestor::IsSecondaryMu ( AliVParticle* track, const AliMCEvent* mcEvent )
269 {
270   /// Muon from secondary decays in absorber
271   Long64_t mask = GetMask(track,mcEvent);
272   return ( TESTBIT(mask,kIsID) && TESTBIT(mask,kIsMuon) && TESTBIT(mask,kIsSecondary) );
273 }
274
275 //_________________________________________________________
276 Bool_t AliUtilityMuonAncestor::IsUnidentified ( AliVParticle* track, const AliMCEvent* mcEvent )
277 {
278   /// Unidentified muon
279   Long64_t mask = GetMask(track,mcEvent);
280   return ( ! TESTBIT(mask,kIsID) );
281 }
282
283 //_________________________________________________________
284 Bool_t AliUtilityMuonAncestor::IsWBosonMu ( AliVParticle* track, const AliMCEvent* mcEvent )
285 {
286   /// Muon from W boson decays
287   return ( IsMuon(track,mcEvent) && CheckAncestor(track,mcEvent,24) );
288 }
289           
290 //_________________________________________________________
291 Bool_t AliUtilityMuonAncestor::IsZBosonMu ( AliVParticle* track, const AliMCEvent* mcEvent )
292 {
293   /// Muon from Z boson decays
294   return ( IsMuon(track,mcEvent) && CheckAncestor(track,mcEvent,kZ0) );
295 }