]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG2/KINK/AliResonanceKink.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
debug flag initialization added in AliResonanceKink
[u/mrichter/AliRoot.git] / PWG2 / KINK / AliResonanceKink.cxx
CommitLineData
10eaad41 1/**************************************************************************
2 * Author: Paraskevi Ganoti, University of Athens (pganoti@phys.uoa.gr) *
3 * Contributors are mentioned in the code where appropriate. *
4 * *
5 * Permission to use, copy, modify and distribute this software and its *
6 * documentation strictly for non-commercial purposes is hereby granted *
7 * without fee, provided that the above copyright notice appears in all *
8 * copies and that both the copyright notice and this permission notice *
9 * appear in the supporting documentation. The authors make no claims *
10 * about the suitability of this software for any purpose. It is *
11 * provided "as is" without express or implied warranty. *
12 **************************************************************************/
13
14//----------------------------------------------------------------------------------------------------------------
f9afc48d 15// class AliResonanceKink
10eaad41 16// Example of an analysis task for reconstructing resonances having at least one kaon-kink in their decay
17// products. It provides basic plots as well as plots helping to calculate the corrections.
18// Usage: To analyse a resonance having a kaon in its decay products, one has to modify the integer
f9afc48d 19// variables resonancePDG, fdaughter1pdg and fdaughter2pdg accordingly as well as daughter1pdgMass and daughter2pdgMass.
10eaad41 20// Also, depending on the analysis mode (ESD or MC), fAnalysisType in the constructor must also be changed
21//-----------------------------------------------------------------------------------------------------------------
22
23#include "TChain.h"
24#include "TTree.h"
25#include "TH2D.h"
26#include "TParticle.h"
92adf4f6 27#include "TDatabasePDG.h"
28#include "TParticlePDG.h"
10eaad41 29#include "TF1.h"
f9afc48d 30#include "TList.h"
31#include "TString.h"
10eaad41 32#include "AliMCEventHandler.h"
33#include "AliMCEvent.h"
f9afc48d 34#include "AliResonanceKink.h"
10eaad41 35#include "AliESDkink.h"
36#include "AliStack.h"
f9afc48d 37#include "AliESDtrack.h"
38#include "AliESDEvent.h"
92adf4f6 39#include "AliExternalTrackParam.h"
10eaad41 40
f9afc48d 41ClassImp(AliResonanceKink)
10eaad41 42
43//________________________________________________________________________
f9afc48d 44AliResonanceKink::AliResonanceKink()
a1cdcf46 45 : TObject(), fDebug(0), fListOfHistos(0), fOpeningAngle(0), fInvariantMass(0), fInvMassTrue(0), fPhiBothKinks(0), fRecPt(0), fRecEta(0), fRecEtaPt(0), fSimPt(0), fSimEta(0), fSimEtaPt(0), fSimPtKink(0), fSimEtaKink(0), fSimEtaPtKink(0),
92adf4f6 46 fhdr(0), fhdz(0), f1(0), f2(0), fAnalysisType(), fvtxz(0), fNbins(0), fLowX(0), fHighX(0), fdaughter1pdg(0), fdaughter2pdg(0), fresonancePDGcode(0), fMaxNSigmaToVertex(0), fMinPtTrackCut(0), fMaxDCAxy(0), fMaxDCAzaxis(0),
47fMinTPCclusters(0),fMaxChi2PerTPCcluster(0), fMaxCov0(0), fMaxCov2(0), fMaxCov5(0) , fMaxCov9(0), fMaxCov14(0) //, fTPCrefitFlag(kFALSE)
10eaad41 48
49{
50 // Constructor
51}
52
53//________________________________________________________________________
f9afc48d 54AliResonanceKink::AliResonanceKink(Int_t nbins, Float_t nlowx, Float_t nhighx, Int_t daughter1, Int_t daughter2, Int_t resonancePDG)
a1cdcf46 55 : TObject(), fDebug(0), fListOfHistos(0), fOpeningAngle(0), fInvariantMass(0), fInvMassTrue(0), fPhiBothKinks(0), fRecPt(0), fRecEta(0), fRecEtaPt(0), fSimPt(0), fSimEta(0), fSimEtaPt(0), fSimPtKink(0), fSimEtaKink(0), fSimEtaPtKink(0),
92adf4f6 56 fhdr(0), fhdz(0), f1(0), f2(0), fAnalysisType(), fvtxz(0), fNbins(nbins), fLowX(nlowx), fHighX(nhighx), fdaughter1pdg(daughter1), fdaughter2pdg(daughter2), fresonancePDGcode(resonancePDG), fMaxNSigmaToVertex(0), fMinPtTrackCut(0),
57fMaxDCAxy(0), fMaxDCAzaxis(0), fMinTPCclusters(0), fMaxChi2PerTPCcluster(0), fMaxCov0(0), fMaxCov2(0), fMaxCov5(0), fMaxCov9(0), fMaxCov14(0) //, fTPCrefitFlag(kFALSE)
10eaad41 58
10eaad41 59{
f9afc48d 60 // Constructor
10eaad41 61
10eaad41 62 fOpeningAngle=new TH1D("fOpeningAngle"," ", 100,-1.0,1.0);
63
f9afc48d 64 fInvariantMass=new TH1D("fInvariantMass"," ",fNbins,fLowX,fHighX);
65 fInvMassTrue=new TH1D("fInvMassTrue"," ",fNbins,fLowX,fHighX);
66 fPhiBothKinks=new TH1D("fPhiBothKinks"," ",fNbins,fLowX,fHighX); // Applicable for phi(1020)
10eaad41 67
68 fRecPt=new TH1D("fRecPt"," ", 50,0.0,5.0);
69 fRecEta=new TH1D("fRecEta"," ", 44,-1.1,1.1);
70 fRecEtaPt=new TH2D("fRecEtaPt"," ", 50,0.0,5.0, 44,-1.1,1.1);
71 fSimPt=new TH1D("fSimPt"," ", 50,0.0,5.0);
72 fSimEta=new TH1D("fSimEta"," ", 44,-1.1,1.1);
73 fSimEtaPt=new TH2D("fSimEtaPt"," ", 50,0.0,5.0, 44,-1.1,1.1);
74 fSimPtKink=new TH1D("fSimPtKink"," ", 50,0.0,5.0);
75 fSimEtaKink=new TH1D("fSimEtaKink"," ", 44,-1.1,1.1);
76 fSimEtaPtKink=new TH2D("fSimEtaPtKink"," ", 50,0.0,5.0, 44,-1.1,1.1);
77 fhdr=new TH1D("fhdr"," ", 100,0.0,5.0);
78 fhdz=new TH1D("fhdz"," ", 100,0.0,5.0);
f9afc48d 79
80 f1=new TF1("f1","((atan([0]*[1]*(1.0/(sqrt((x^2)*(1.0-([1]^2))-([0]^2)*([1]^2))))))*180.)/[2]",1.1,10.0);
81 f1->SetParameter(0,0.493677);
82 f1->SetParameter(1,0.9127037);
83 f1->SetParameter(2,TMath::Pi());
84
85 f2=new TF1("f2","((atan([0]*[1]*(1.0/(sqrt((x^2)*(1.0-([1]^2))-([0]^2)*([1]^2))))))*180.)/[2]",0.1,10.0);
86 f2->SetParameter(0,0.13957018);
87 f2->SetParameter(1,0.2731374);
88 f2->SetParameter(2,TMath::Pi());
89
10eaad41 90 fvtxz=new TH1D("fvtxz"," ", 100,-20.0,20.0);
92adf4f6 91
10eaad41 92}
93
94//________________________________________________________________________
f9afc48d 95AliResonanceKink:: ~AliResonanceKink()
10eaad41 96{
f9afc48d 97 // Destructor
98 if(fOpeningAngle) delete fOpeningAngle;
99 if(fInvariantMass) delete fInvariantMass;
100 if(fInvMassTrue) delete fInvMassTrue;
101 if(fPhiBothKinks) delete fPhiBothKinks;
102 if(fRecPt) delete fRecPt;
103 if(fRecEta) delete fRecEta;
104 if(fRecEtaPt) delete fRecEtaPt;
105 if(fSimPt) delete fSimPt;
106 if(fSimEta) delete fSimEta;
107 if(fSimEtaPt) delete fSimEtaPt;
108 if(fSimPtKink) delete fSimPtKink;
109 if(fSimEtaKink) delete fSimEtaKink;
110 if(fSimEtaPtKink) delete fSimEtaPtKink;
111 if(fhdr) delete fhdr;
112 if(fhdz) delete fhdz;
113 if(fvtxz) delete fvtxz;
114}
115//________________________________________________________________________
116TList* AliResonanceKink::GetHistogramList()
117{
118 // Adding histograms to the list
119 fListOfHistos=new TList();
10eaad41 120
f9afc48d 121 fListOfHistos->Add(fOpeningAngle);
122 fListOfHistos->Add(fInvariantMass);
123 fListOfHistos->Add(fInvMassTrue);
124 fListOfHistos->Add(fPhiBothKinks);
125 fListOfHistos->Add(fRecPt);
126 fListOfHistos->Add(fRecEta);
127 fListOfHistos->Add(fRecEtaPt);
128 fListOfHistos->Add(fSimPt);
129 fListOfHistos->Add(fSimEta);
130 fListOfHistos->Add(fSimEtaPt);
131 fListOfHistos->Add(fSimPtKink);
132 fListOfHistos->Add(fSimEtaKink);
133 fListOfHistos->Add(fSimEtaPtKink);
134 fListOfHistos->Add(fhdr);
135 fListOfHistos->Add(fhdz);
136 fListOfHistos->Add(fvtxz);
10eaad41 137
f9afc48d 138 return fListOfHistos;
139}
140
141//________________________________________________________________________
142void AliResonanceKink::InitOutputHistograms(Int_t nbins, Float_t nlowx, Float_t nhighx)
143{
144 // Initialisation of the output histograms
145 fNbins=nbins;
146 fLowX=nlowx;
147 fHighX=nhighx;
148
149 fOpeningAngle=new TH1D("fOpeningAngle"," ", 100,-1.0,1.0);
150
151 fInvariantMass=new TH1D("fInvariantMass"," ",fNbins,fLowX,fHighX);
152 fInvMassTrue=new TH1D("fInvMassTrue"," ",fNbins,fLowX,fHighX);
153 fPhiBothKinks=new TH1D("fPhiBothKinks"," ",fNbins,fLowX,fHighX); // Applicable for phi(1020)
154
155 fRecPt=new TH1D("fRecPt"," ", 50,0.0,5.0);
156 fRecEta=new TH1D("fRecEta"," ", 44,-1.1,1.1);
157 fRecEtaPt=new TH2D("fRecEtaPt"," ", 50,0.0,5.0, 44,-1.1,1.1);
158 fSimPt=new TH1D("fSimPt"," ", 50,0.0,5.0);
159 fSimEta=new TH1D("fSimEta"," ", 44,-1.1,1.1);
160 fSimEtaPt=new TH2D("fSimEtaPt"," ", 50,0.0,5.0, 44,-1.1,1.1);
161 fSimPtKink=new TH1D("fSimPtKink"," ", 50,0.0,5.0);
162 fSimEtaKink=new TH1D("fSimEtaKink"," ", 44,-1.1,1.1);
163 fSimEtaPtKink=new TH2D("fSimEtaPtKink"," ", 50,0.0,5.0, 44,-1.1,1.1);
164 fhdr=new TH1D("fhdr"," ", 100,0.0,5.0);
165 fhdz=new TH1D("fhdz"," ", 100,0.0,5.0);
10eaad41 166
f9afc48d 167 f1=new TF1("f1","((atan([0]*[1]*(1.0/(sqrt((x^2)*(1.0-([1]^2))-([0]^2)*([1]^2))))))*180.)/[2]",1.1,10.0);
168 f1->SetParameter(0,0.493677);
169 f1->SetParameter(1,0.9127037);
170 f1->SetParameter(2,TMath::Pi());
171
172 f2=new TF1("f2","((atan([0]*[1]*(1.0/(sqrt((x^2)*(1.0-([1]^2))-([0]^2)*([1]^2))))))*180.)/[2]",0.1,10.0);
173 f2->SetParameter(0,0.13957018);
174 f2->SetParameter(1,0.2731374);
175 f2->SetParameter(2,TMath::Pi());
f27d6e67 176
f9afc48d 177 fvtxz=new TH1D("fvtxz"," ", 100,-20.0,20.0);
178}
179
180//________________________________________________________________________
181void AliResonanceKink::Analyse(AliESDEvent* esd, AliMCEvent* mcEvent)
182{
183 // Main loop
184 // Called for each event
185 Int_t resonancePDGcode, antiresonancePDGcode;
186
187 if (fdaughter1pdg==kdaughterKaon) {
188 resonancePDGcode=fresonancePDGcode;
189 antiresonancePDGcode=-fresonancePDGcode;
10eaad41 190 }
f9afc48d 191 if (fdaughter1pdg!=kdaughterKaon) {
192 resonancePDGcode=-fresonancePDGcode;
193 antiresonancePDGcode=fresonancePDGcode;
194 }
195 if (fdaughter1pdg==fdaughter2pdg) {
196 resonancePDGcode=fresonancePDGcode;
197 antiresonancePDGcode=fresonancePDGcode;
198 }
10eaad41 199
f9afc48d 200 Double_t daughter1pdgMass=TDatabasePDG::Instance()->GetParticle(fdaughter1pdg)->Mass();
201 Double_t daughter2pdgMass=TDatabasePDG::Instance()->GetParticle(fdaughter2pdg)->Mass();
202
203 if (!esd) {
10eaad41 204 Printf("ERROR: fESD not available");
205 return;
206 }
207
f9afc48d 208 if (!mcEvent) {
209 Printf("ERROR: mcEvent not available");
210 return;
211 }
212
10eaad41 213 AliStack* stack=mcEvent->Stack();
214
215 if(fAnalysisType == "MC") {
10eaad41 216 for (Int_t iMc = 0; iMc < stack->GetNprimary(); ++iMc)
217 {
218 TParticle* particle = stack->Particle(iMc);
219
220 if (!particle)
221 {
1ac5c71d 222 if (fDebug > 0) Printf("UNEXPECTED: particle with label %d not found in stack (mc loop)", iMc);
10eaad41 223 continue;
224 }
225
f9afc48d 226 if(TMath::Abs(particle->GetPdgCode())==fresonancePDGcode) {
10eaad41 227 Int_t firstD=particle->GetFirstDaughter();
228 Int_t lastD=particle->GetLastDaughter();
229 TParticle *daughterParticle1=stack->Particle(firstD);
230 TParticle *daughterParticle2=stack->Particle(lastD);
231
232 TParticle* kaonFirstDaughter;
e40a0beb 233 Int_t mcProcessKaonFirstDaughter = -999;
10eaad41 234
235 for(Int_t ia=0; ia<daughterParticle1->GetNDaughters(); ia++){
236 if ((daughterParticle1->GetFirstDaughter()+ia)!=-1) {
237 kaonFirstDaughter=stack->Particle(daughterParticle1->GetFirstDaughter()+ia);
238 mcProcessKaonFirstDaughter=kaonFirstDaughter->GetUniqueID();
239 }
240 }
241
242 if((daughterParticle1->Pt()>0.25)&&(daughterParticle2->Pt()>0.25)&&(TMath::Abs(daughterParticle1->Eta())<1.1)&& (TMath::Abs(daughterParticle2->Eta())<1.1)&&(TMath::Abs(particle->Eta())<1.1)) {
243 fSimEta->Fill(particle->Eta());
244 fSimPt->Fill(particle->Pt());
245 fSimEtaPt->Fill(particle->Pt(), particle->Eta());
246 if(mcProcessKaonFirstDaughter==4) {
247 fSimPtKink->Fill(particle->Pt());
248 fSimEtaKink->Fill(particle->Eta());
249 fSimEtaPtKink->Fill(particle->Pt(), particle->Eta());
250 }
251 }
252 }
253 }
254
255 } // end fAnalysisType==MC
256 else
257
258 if(fAnalysisType == "ESD") {
f9afc48d 259 const AliESDVertex* vertex = GetEventVertex(esd);
10eaad41 260 if(!vertex) return;
261 Double_t vtx[3];
262 vertex->GetXYZ(vtx);
263 fvtxz->Fill(vtx[2]);
10eaad41 264 Double_t ptrackpos[3], ptrackneg[3];
265
266 TLorentzVector p4pos, anp4pos;
267 TLorentzVector p4neg, anp4neg;
268 TLorentzVector p4comb, anp4comb;
269
f9afc48d 270 for (Int_t iTracks = 0; iTracks < esd->GetNumberOfTracks(); iTracks++) {
271 AliESDtrack* trackpos = esd->GetTrack(iTracks);
10eaad41 272 if (!trackpos) {
1ac5c71d 273 if (fDebug > 0) Printf("ERROR: Could not receive track %d", iTracks);
10eaad41 274 continue;
275 }
276 if (trackpos->GetSign() < 0) continue;
277
92adf4f6 278 AliExternalTrackParam *tpcTrackpos = (AliExternalTrackParam *)trackpos->GetTPCInnerParam();
279 if(!tpcTrackpos) continue;
280 ptrackpos[0]=tpcTrackpos->Px();
281 ptrackpos[1]=tpcTrackpos->Py();
282 ptrackpos[2]=tpcTrackpos->Pz();
10eaad41 283
92adf4f6 284 Bool_t firstLevelAcceptPosTrack=IsAcceptedForKink(esd, vertex, trackpos);
285 if(firstLevelAcceptPosTrack==kFALSE) continue;
10eaad41 286
287 TVector3 posTrackMom(ptrackpos[0],ptrackpos[1],ptrackpos[2]);
92adf4f6 288
10eaad41 289 TParticle * partpos = stack->Particle(TMath::Abs(trackpos->GetLabel()));
290 if (!partpos) continue;
291 Int_t pdgpos = partpos->GetPdgCode();
292 Int_t mumlabelpos=partpos->GetFirstMother();
293 mumlabelpos = TMath::Abs(mumlabelpos);
294 TParticle * mumpos=stack->Particle(mumlabelpos);
295 if (!mumpos) continue;
296 Int_t mumpdgpos = mumpos->GetPdgCode();
297
298 Int_t indexKinkPos=trackpos->GetKinkIndex(0);
92adf4f6 299 Bool_t posKaonKinkFlag=0;
300 if(indexKinkPos<0) posKaonKinkFlag=IsKink(esd, indexKinkPos, posTrackMom);
10eaad41 301
92adf4f6 302 if(posKaonKinkFlag==1) anp4pos.SetVectM(posTrackMom,daughter1pdgMass);
10eaad41 303
304 if(indexKinkPos==0) {
92adf4f6 305
306 Bool_t secondLevelAcceptPosTrack=IsAcceptedForTrack(esd, vertex, trackpos);
307 if(secondLevelAcceptPosTrack==kFALSE) continue;
308
f9afc48d 309 p4pos.SetVectM(posTrackMom, daughter2pdgMass);
10eaad41 310
311 }
312
f9afc48d 313 for (Int_t j=0; j<esd->GetNumberOfTracks(); j++) {
10eaad41 314 if(iTracks==j) continue;
f9afc48d 315 AliESDtrack* trackneg=esd->GetTrack(j);
10eaad41 316 if (trackneg->GetSign() > 0) continue;
317
92adf4f6 318 AliExternalTrackParam *tpcTrackneg = (AliExternalTrackParam *)trackneg->GetTPCInnerParam();
319 if(!tpcTrackneg) continue;
320 ptrackneg[0]=tpcTrackneg->Px();
321 ptrackneg[1]=tpcTrackneg->Py();
322 ptrackneg[2]=tpcTrackneg->Pz();
10eaad41 323
92adf4f6 324 Bool_t firstLevelAcceptNegTrack=IsAcceptedForKink(esd, vertex, trackneg);
325 if(firstLevelAcceptNegTrack==kFALSE) continue;
10eaad41 326
327 TVector3 negTrackMom(ptrackneg[0],ptrackneg[1],ptrackneg[2]);
10eaad41 328
329 TParticle * partneg = stack->Particle(TMath::Abs(trackneg->GetLabel()));
330 if (!partneg) continue;
331 Int_t pdgneg = partneg->GetPdgCode();
332 Int_t mumlabelneg=partneg->GetFirstMother();
333 mumlabelneg = TMath::Abs(mumlabelneg);
334 TParticle * mumneg=stack->Particle(mumlabelneg);
335 if (!mumneg) continue;
336 Int_t mumpdgneg = mumneg->GetPdgCode();
337
338 Int_t indexKinkNeg=trackneg->GetKinkIndex(0);
92adf4f6 339 Bool_t negKaonKinkFlag=0;
340 if(indexKinkNeg<0) negKaonKinkFlag=IsKink(esd, indexKinkNeg, negTrackMom);
10eaad41 341
f9afc48d 342 if(negKaonKinkFlag==1) p4neg.SetVectM(negTrackMom,daughter1pdgMass);
10eaad41 343
344 if(indexKinkNeg==0) {
92adf4f6 345
346 Bool_t secondLevelAcceptNegTrack=IsAcceptedForTrack(esd, vertex, trackneg);
347 if(secondLevelAcceptNegTrack==kFALSE) continue;
348
f9afc48d 349 anp4neg.SetVectM(negTrackMom, daughter2pdgMass);
10eaad41 350
351 }
352
353 Double_t openingAngle=(ptrackpos[0]*ptrackneg[0]+ptrackpos[1]*ptrackneg[1]+ptrackpos[2]*ptrackneg[2])/(posTrackMom.Mag()*negTrackMom.Mag());
354
92adf4f6 355 if((posKaonKinkFlag==1)&&(negKaonKinkFlag==1)) {
10eaad41 356 p4comb=anp4pos;
357 p4comb+=p4neg;
358 if(openingAngle>0.6) fPhiBothKinks->Fill(p4comb.M());
359 }
360
361 if(negKaonKinkFlag==1) {
362 p4comb=p4pos;
363 p4comb+=p4neg;
364 fInvariantMass->Fill(p4comb.M());
f9afc48d 365 if ((mumpdgpos==(antiresonancePDGcode))&&(mumpdgneg==(antiresonancePDGcode))&&(mumlabelpos==mumlabelneg)
366 &&(pdgpos==fdaughter2pdg)&&(pdgneg==(-fdaughter1pdg))&&(TMath::Abs(trackpos->GetLabel())>=0)&&(TMath::Abs(trackneg->GetLabel())>=0)&&(mumlabelpos>=0)&&(mumlabelneg>=0)) {
10eaad41 367 fOpeningAngle->Fill(openingAngle);
368 fInvMassTrue->Fill(p4comb.M());
369 if((TMath::Abs(p4pos.Vect().Eta())<1.1)&&(TMath::Abs(p4neg.Vect().Eta())<1.1)&&(p4comb.Vect().Eta()<1.1)) {
370 fRecPt->Fill(p4comb.Vect().Pt());
371 fRecEta->Fill(p4comb.Vect().Eta());
372 fRecEtaPt->Fill(p4comb.Vect().Perp(),p4comb.Vect().Eta());
373
374 }
375
376 }
377
378 }
379
92adf4f6 380 if(posKaonKinkFlag==1) {
10eaad41 381 anp4comb=anp4pos;
382 anp4comb+=anp4neg;
383 fInvariantMass->Fill(anp4comb.M());
f9afc48d 384 if ((mumpdgpos==resonancePDGcode)&&(mumpdgneg==resonancePDGcode)&&(mumlabelpos==mumlabelneg)
385 &&(pdgpos==fdaughter1pdg)&&(pdgneg==(-fdaughter2pdg))&&(TMath::Abs(trackpos->GetLabel())>=0)&&(TMath::Abs(trackneg->GetLabel())>=0)&&(mumlabelpos>=0) &&(mumlabelneg>=0)) {
10eaad41 386 fOpeningAngle->Fill(openingAngle);
387 fInvMassTrue->Fill(p4comb.M());
388 if((TMath::Abs(anp4neg.Vect().Eta())<1.1)&&(TMath::Abs(anp4pos.Vect().Eta())<1.1)&&(anp4comb.Vect().Eta()<1.1)) {
389 fRecPt->Fill(anp4comb.Vect().Pt());
390 fRecEta->Fill(anp4comb.Vect().Eta());
391 fRecEtaPt->Fill(anp4comb.Vect().Pt(), anp4comb.Vect().Eta());
392 }
393
394 }
395
396 }
397
398 } //inner track loop
399
400 } //outer track loop
401
402 } // end fAnalysisType == ESD
403
10eaad41 404}
405
10eaad41 406//____________________________________________________________________//
f9afc48d 407Float_t AliResonanceKink::GetSigmaToVertex(AliESDtrack* esdTrack) const {
10eaad41 408 // Calculates the number of sigma to the vertex.
409
410 Float_t b[2];
411 Float_t bRes[2];
412 Float_t bCov[3];
413
92adf4f6 414 esdTrack->GetImpactParametersTPC(b,bCov);
10eaad41 415
416 if (bCov[0]<=0 || bCov[2]<=0) {
10eaad41 417 bCov[0]=0; bCov[2]=0;
418 }
f9afc48d 419
10eaad41 420 bRes[0] = TMath::Sqrt(bCov[0]);
421 bRes[1] = TMath::Sqrt(bCov[2]);
422
423 if (bRes[0] == 0 || bRes[1] ==0) return -1;
424
425 Float_t d = TMath::Sqrt(TMath::Power(b[0]/bRes[0],2) + TMath::Power(b[1]/bRes[1],2));
426
427 if (TMath::Exp(-d * d / 2) < 1e-10) return 1000;
428
429 d = TMath::ErfInverse(1 - TMath::Exp(-d * d / 2)) * TMath::Sqrt(2);
430
431 return d;
432}
433
f9afc48d 434//________________________________________________________________________
435const AliESDVertex* AliResonanceKink::GetEventVertex(const AliESDEvent* esd) const
10eaad41 436{
437 // Get the vertex
438
439 const AliESDVertex* vertex = esd->GetPrimaryVertex();
440
441 if((vertex->GetStatus()==kTRUE)&&(vertex->GetNContributors()>2)) return vertex;
442 else
443 {
444 vertex = esd->GetPrimaryVertexSPD();
445 if((vertex->GetStatus()==kTRUE)&&(vertex->GetNContributors()>2)) return vertex;
446 else
447 return 0;
448 }
449}
450
92adf4f6 451//________________________________________________________________________
452
453 Bool_t AliResonanceKink::IsAcceptedForKink(AliESDEvent *localesd,
454 const AliESDVertex *localvertex, AliESDtrack* localtrack) {
455 // Apply the selections for each kink
456
457 Double_t gPt = 0.0, gPx = 0.0, gPy = 0.0, gPz = 0.0;
458 Double_t dca[2] = {0.0,0.0}, cov[3] = {0.0,0.0,0.0}; //The impact parameters and their covariance.
459 Double_t dca3D = 0.0;
460
461 AliExternalTrackParam *tpcTrack = (AliExternalTrackParam *)localtrack->GetTPCInnerParam();
462 if(!tpcTrack) {
463 gPt = 0.0; gPx = 0.0; gPy = 0.0; gPz = 0.0;
464 dca[0] = -100.; dca[1] = -100.; dca3D = -100.;
465 cov[0] = -100.; cov[1] = -100.; cov[2] = -100.;
466 }
467 else {
468 gPt = tpcTrack->Pt();
469 gPx = tpcTrack->Px();
470 gPy = tpcTrack->Py();
471 gPz = tpcTrack->Pz();
472 tpcTrack->PropagateToDCA(localvertex,
473 localesd->GetMagneticField(),100.,dca,cov);
474 }
475
476 if(GetSigmaToVertex(localtrack) > fMaxNSigmaToVertex) {
1ac5c71d 477 if (fDebug > 1) Printf("IsAcceptedKink: Track rejected because it has a %lf sigmas to vertex TPC (max. requested: %lf)", GetSigmaToVertex(localtrack),fMaxNSigmaToVertex);
92adf4f6 478 return kFALSE;
479 }
480
481 if(TMath::Abs(dca[0]) > fMaxDCAxy) {
1ac5c71d 482 if (fDebug > 1) Printf("IsAcceptedKink: Track rejected because it has a value of dca(xy) (TPC) of %lf (max. requested: %lf)", TMath::Abs(dca[0]), fMaxDCAxy);
92adf4f6 483 return kFALSE;
484 }
485
486 if(TMath::Abs(dca[1]) > fMaxDCAzaxis) {
1ac5c71d 487 if (fDebug > 1) Printf("IsAcceptedKink: Track rejected because it has a value of dca(z) of %lf (max. requested: %lf)", TMath::Abs(dca[1]), fMaxDCAzaxis);
92adf4f6 488 return kFALSE;
489 }
490
491 if(gPt < fMinPtTrackCut) {
1ac5c71d 492 if (fDebug > 1) Printf("IsAcceptedKink: Track rejected because it has a min value of pt of %lf (min. requested: %lf)", gPt, fMinPtTrackCut);
92adf4f6 493 return kFALSE;
494 }
495
496 return kTRUE;
497}
498
499//________________________________________________________________________
500Bool_t AliResonanceKink::IsAcceptedForTrack(AliESDEvent *localesd, const AliESDVertex *localvertex, AliESDtrack *localtrack) {
501 // Apply the selections for each track
f9afc48d 502
92adf4f6 503 Double_t gPt = 0.0, gPx = 0.0, gPy = 0.0, gPz = 0.0;
504 Double_t dca[2] = {0.0,0.0}, cov[3] = {0.0,0.0,0.0}; //The impact parameters and their covariance.
505 Double_t dca3D = 0.0;
506
507 AliExternalTrackParam *tpcTrack = (AliExternalTrackParam *)localtrack->GetTPCInnerParam();
508 if(!tpcTrack) {
509 gPt = 0.0; gPx = 0.0; gPy = 0.0; gPz = 0.0;
510 dca[0] = -100.; dca[1] = -100.; dca3D = -100.;
511 cov[0] = -100.; cov[1] = -100.; cov[2] = -100.;
512 }
513 else {
514 gPt = tpcTrack->Pt();
515 gPx = tpcTrack->Px();
516 gPy = tpcTrack->Py();
517 gPz = tpcTrack->Pz();
518 tpcTrack->PropagateToDCA(localvertex,
519 localesd->GetMagneticField(),100.,dca,cov);
520 }
521
522 Int_t fcls[200];
523 Int_t nClustersTPC=localtrack->GetTPCclusters(fcls);
524 Float_t chi2perTPCcluster=-1.0;
525 if(nClustersTPC!=0) chi2perTPCcluster=(localtrack->GetTPCchi2())/Float_t(nClustersTPC);
526
527 Double_t extCov[15];
528 localtrack->GetExternalCovariance(extCov);
529
530 if((localtrack->GetStatus() & AliESDtrack::kTPCrefit) == 0) {
1ac5c71d 531 if (fDebug > 1) Printf("IsAccepted: Track rejected because of no refited in TPC");
92adf4f6 532 return kFALSE;
533 }
534
535 if(nClustersTPC < fMinTPCclusters) {
1ac5c71d 536 if (fDebug > 1) Printf("IsAccepted: Track rejected because it has a value of nclusters (TPC) of %ld (min. requested: %ld)", nClustersTPC, fMinTPCclusters);
92adf4f6 537 return kFALSE;
538 }
539
540 if(chi2perTPCcluster > fMaxChi2PerTPCcluster) {
1ac5c71d 541 if (fDebug > 1) Printf("IsAccepted: Track rejected because it has a value of chi2perTPCcluster of %lf (max. requested: %lf)", chi2perTPCcluster, fMaxChi2PerTPCcluster);
92adf4f6 542 return kFALSE;
543 }
544
545 if(extCov[0] > fMaxCov0) {
1ac5c71d 546 if (fDebug > 1) Printf("IsAccepted: Track rejected because it has a value of cov[0] of %lf (max. requested: %lf)", cov[0], fMaxCov0);
92adf4f6 547 return kFALSE;
548 }
549
550 if(extCov[2] > fMaxCov2) {
1ac5c71d 551 if (fDebug > 1) Printf("IsAccepted: Track rejected because it has a value of cov[2] of %lf (max. requested: %lf)", cov[2], fMaxCov2);
92adf4f6 552 return kFALSE;
553 }
554
555 if(extCov[5] > fMaxCov5) {
1ac5c71d 556 if (fDebug > 1) Printf("IsAccepted: Track rejected because it has a value of cov[5] of %lf (max. requested: %lf)", cov[5], fMaxCov5);
92adf4f6 557 return kFALSE;
558 }
559
560 if(extCov[9] > fMaxCov9) {
1ac5c71d 561 if (fDebug > 1) Printf("IsAccepted: Track rejected because it has a value of cov[9] of %lf (max. requested: %lf)", cov[9], fMaxCov9);
92adf4f6 562 return kFALSE;
563 }
564
565 if(extCov[14] > fMaxCov14) {
1ac5c71d 566 if (fDebug > 1) Printf("IsAccepted: Track rejected because it has a value of cov[14] of %lf (max. requested: %lf)", cov[14], fMaxCov14);
92adf4f6 567 return kFALSE;
568 }
f9afc48d 569
92adf4f6 570 return kTRUE;
571
572}
573
574//________________________________________________________________________
575Bool_t AliResonanceKink::IsKink(AliESDEvent *localesd, Int_t kinkIndex, TVector3 trackMom)
576{
577 // Test some kinematical criteria for each kink
578
579 AliESDkink *kink=localesd->GetKink(TMath::Abs(kinkIndex)-1);
580 const TVector3 motherMfromKink(kink->GetMotherP());
581 const TVector3 daughterMKink(kink->GetDaughterP());
582 Float_t qt=kink->GetQt();
583
584 Double_t maxDecAngKmu=f1->Eval(motherMfromKink.Mag(),0.,0.,0.);
585 Double_t maxDecAngpimu=f2->Eval(motherMfromKink.Mag(),0.,0.,0.);
586
587 Float_t kinkAngle=TMath::RadToDeg()*kink->GetAngle(2);
588
589 Float_t energyDaughterMu=TMath::Sqrt(daughterMKink.Mag()*daughterMKink.Mag()+0.105658*0.105658);
590 Float_t p1XM= motherMfromKink.Px();
591 Float_t p1YM= motherMfromKink.Py();
592 Float_t p1ZM= motherMfromKink.Pz();
593 Float_t p2XM= daughterMKink.Px();
594 Float_t p2YM= daughterMKink.Py();
595 Float_t p2ZM= daughterMKink.Pz();
596 Float_t p3Daughter=TMath::Sqrt(((p1XM-p2XM)*(p1XM-p2XM))+((p1YM-p2YM)*(p1YM-p2YM))+((p1ZM-p2ZM)*(p1ZM-p2ZM)));
597 Double_t invariantMassKmu= TMath::Sqrt((energyDaughterMu+p3Daughter)*(energyDaughterMu+p3Daughter)-motherMfromKink.Mag()*motherMfromKink.Mag());
598
599 if((kinkAngle>maxDecAngpimu)&&(qt>0.05)&&(qt<0.25)&&((kink->GetR()>110.)&&(kink->GetR()<230.))&&(TMath::Abs(trackMom.Eta())<1.1)&&(invariantMassKmu<0.6)) {
600
601 if(trackMom.Mag()<=1.1) {
602 return kTRUE;
603 }
604 else
605 if (kinkAngle<maxDecAngKmu) {
606 return kTRUE;
607 }
608 }
609 return kFALSE;
610}