]> 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
one more change (from P.Ganoti pganoti@phys.uoa.gr)
[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"
7ccf0419 40#include "AliMCParticle.h"
10eaad41 41
f9afc48d 42ClassImp(AliResonanceKink)
10eaad41 43
44//________________________________________________________________________
f9afc48d 45AliResonanceKink::AliResonanceKink()
a1cdcf46 46 : 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 47 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),
48fMinTPCclusters(0),fMaxChi2PerTPCcluster(0), fMaxCov0(0), fMaxCov2(0), fMaxCov5(0) , fMaxCov9(0), fMaxCov14(0) //, fTPCrefitFlag(kFALSE)
7ccf0419 49, fInvmassPt(0), fInvmassPtTrue(0), fMCInvmassPt(0), fMCInvmassPtTrue(0)
10eaad41 50{
51 // Constructor
52}
53
54//________________________________________________________________________
f9afc48d 55AliResonanceKink::AliResonanceKink(Int_t nbins, Float_t nlowx, Float_t nhighx, Int_t daughter1, Int_t daughter2, Int_t resonancePDG)
a1cdcf46 56 : 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 57 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),
58fMaxDCAxy(0), fMaxDCAzaxis(0), fMinTPCclusters(0), fMaxChi2PerTPCcluster(0), fMaxCov0(0), fMaxCov2(0), fMaxCov5(0), fMaxCov9(0), fMaxCov14(0) //, fTPCrefitFlag(kFALSE)
7ccf0419 59, fInvmassPt(0), fInvmassPtTrue(0), fMCInvmassPt(0), fMCInvmassPtTrue(0)
10eaad41 60{
f9afc48d 61 // Constructor
10eaad41 62
10eaad41 63 fOpeningAngle=new TH1D("fOpeningAngle"," ", 100,-1.0,1.0);
64
f9afc48d 65 fInvariantMass=new TH1D("fInvariantMass"," ",fNbins,fLowX,fHighX);
66 fInvMassTrue=new TH1D("fInvMassTrue"," ",fNbins,fLowX,fHighX);
67 fPhiBothKinks=new TH1D("fPhiBothKinks"," ",fNbins,fLowX,fHighX); // Applicable for phi(1020)
10eaad41 68
69 fRecPt=new TH1D("fRecPt"," ", 50,0.0,5.0);
70 fRecEta=new TH1D("fRecEta"," ", 44,-1.1,1.1);
71 fRecEtaPt=new TH2D("fRecEtaPt"," ", 50,0.0,5.0, 44,-1.1,1.1);
72 fSimPt=new TH1D("fSimPt"," ", 50,0.0,5.0);
73 fSimEta=new TH1D("fSimEta"," ", 44,-1.1,1.1);
74 fSimEtaPt=new TH2D("fSimEtaPt"," ", 50,0.0,5.0, 44,-1.1,1.1);
75 fSimPtKink=new TH1D("fSimPtKink"," ", 50,0.0,5.0);
76 fSimEtaKink=new TH1D("fSimEtaKink"," ", 44,-1.1,1.1);
77 fSimEtaPtKink=new TH2D("fSimEtaPtKink"," ", 50,0.0,5.0, 44,-1.1,1.1);
78 fhdr=new TH1D("fhdr"," ", 100,0.0,5.0);
79 fhdz=new TH1D("fhdz"," ", 100,0.0,5.0);
f9afc48d 80
81 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);
82 f1->SetParameter(0,0.493677);
83 f1->SetParameter(1,0.9127037);
84 f1->SetParameter(2,TMath::Pi());
85
86 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);
87 f2->SetParameter(0,0.13957018);
88 f2->SetParameter(1,0.2731374);
89 f2->SetParameter(2,TMath::Pi());
90
10eaad41 91 fvtxz=new TH1D("fvtxz"," ", 100,-20.0,20.0);
92adf4f6 92
7ccf0419 93 fInvmassPt=new TH2D("fInvmassPt"," ",fNbins,fLowX,fHighX,100,0.0,10.0);
94 fInvmassPtTrue=new TH2D("fInvmassPtTrue"," ",fNbins,fLowX,fHighX,100,0.0,10.0);
95 fMCInvmassPt=new TH2D("fMCInvmassPt"," ",fNbins,fLowX,fHighX,100,0.0,10.0);
96 fMCInvmassPtTrue=new TH2D("fMCInvmassPtTrue"," ",fNbins,fLowX,fHighX,100,0.0,10.0);
10eaad41 97}
98
99//________________________________________________________________________
f9afc48d 100AliResonanceKink:: ~AliResonanceKink()
10eaad41 101{
f9afc48d 102 // Destructor
103 if(fOpeningAngle) delete fOpeningAngle;
104 if(fInvariantMass) delete fInvariantMass;
105 if(fInvMassTrue) delete fInvMassTrue;
106 if(fPhiBothKinks) delete fPhiBothKinks;
107 if(fRecPt) delete fRecPt;
108 if(fRecEta) delete fRecEta;
109 if(fRecEtaPt) delete fRecEtaPt;
110 if(fSimPt) delete fSimPt;
111 if(fSimEta) delete fSimEta;
112 if(fSimEtaPt) delete fSimEtaPt;
113 if(fSimPtKink) delete fSimPtKink;
114 if(fSimEtaKink) delete fSimEtaKink;
115 if(fSimEtaPtKink) delete fSimEtaPtKink;
116 if(fhdr) delete fhdr;
117 if(fhdz) delete fhdz;
7ccf0419 118 if(fvtxz) delete fvtxz;
119 if(fInvmassPt) delete fInvmassPt;
120 if(fInvmassPtTrue) delete fInvmassPtTrue;
121 if(fMCInvmassPt) delete fMCInvmassPt;
122 if(fMCInvmassPtTrue) delete fMCInvmassPtTrue;
f9afc48d 123}
124//________________________________________________________________________
125TList* AliResonanceKink::GetHistogramList()
126{
127 // Adding histograms to the list
128 fListOfHistos=new TList();
10eaad41 129
f9afc48d 130 fListOfHistos->Add(fOpeningAngle);
131 fListOfHistos->Add(fInvariantMass);
132 fListOfHistos->Add(fInvMassTrue);
133 fListOfHistos->Add(fPhiBothKinks);
134 fListOfHistos->Add(fRecPt);
135 fListOfHistos->Add(fRecEta);
136 fListOfHistos->Add(fRecEtaPt);
137 fListOfHistos->Add(fSimPt);
138 fListOfHistos->Add(fSimEta);
139 fListOfHistos->Add(fSimEtaPt);
140 fListOfHistos->Add(fSimPtKink);
141 fListOfHistos->Add(fSimEtaKink);
142 fListOfHistos->Add(fSimEtaPtKink);
143 fListOfHistos->Add(fhdr);
144 fListOfHistos->Add(fhdz);
145 fListOfHistos->Add(fvtxz);
7ccf0419 146 fListOfHistos->Add(fInvmassPt);
147 fListOfHistos->Add(fInvmassPtTrue);
148 fListOfHistos->Add(fMCInvmassPt);
149 fListOfHistos->Add(fMCInvmassPtTrue);
150
f9afc48d 151 return fListOfHistos;
152}
153
154//________________________________________________________________________
155void AliResonanceKink::InitOutputHistograms(Int_t nbins, Float_t nlowx, Float_t nhighx)
156{
157 // Initialisation of the output histograms
158 fNbins=nbins;
159 fLowX=nlowx;
160 fHighX=nhighx;
161
162 fOpeningAngle=new TH1D("fOpeningAngle"," ", 100,-1.0,1.0);
163
164 fInvariantMass=new TH1D("fInvariantMass"," ",fNbins,fLowX,fHighX);
165 fInvMassTrue=new TH1D("fInvMassTrue"," ",fNbins,fLowX,fHighX);
166 fPhiBothKinks=new TH1D("fPhiBothKinks"," ",fNbins,fLowX,fHighX); // Applicable for phi(1020)
167
168 fRecPt=new TH1D("fRecPt"," ", 50,0.0,5.0);
169 fRecEta=new TH1D("fRecEta"," ", 44,-1.1,1.1);
170 fRecEtaPt=new TH2D("fRecEtaPt"," ", 50,0.0,5.0, 44,-1.1,1.1);
171 fSimPt=new TH1D("fSimPt"," ", 50,0.0,5.0);
172 fSimEta=new TH1D("fSimEta"," ", 44,-1.1,1.1);
173 fSimEtaPt=new TH2D("fSimEtaPt"," ", 50,0.0,5.0, 44,-1.1,1.1);
174 fSimPtKink=new TH1D("fSimPtKink"," ", 50,0.0,5.0);
175 fSimEtaKink=new TH1D("fSimEtaKink"," ", 44,-1.1,1.1);
176 fSimEtaPtKink=new TH2D("fSimEtaPtKink"," ", 50,0.0,5.0, 44,-1.1,1.1);
177 fhdr=new TH1D("fhdr"," ", 100,0.0,5.0);
178 fhdz=new TH1D("fhdz"," ", 100,0.0,5.0);
10eaad41 179
f9afc48d 180 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);
181 f1->SetParameter(0,0.493677);
182 f1->SetParameter(1,0.9127037);
183 f1->SetParameter(2,TMath::Pi());
184
185 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);
186 f2->SetParameter(0,0.13957018);
187 f2->SetParameter(1,0.2731374);
188 f2->SetParameter(2,TMath::Pi());
f27d6e67 189
f9afc48d 190 fvtxz=new TH1D("fvtxz"," ", 100,-20.0,20.0);
7ccf0419 191
192 fInvmassPt=new TH2D("fInvmassPt"," ",fNbins,fLowX,fHighX,100,0.0,10.0);
193 fInvmassPtTrue=new TH2D("fInvmassPtTrue"," ",fNbins,fLowX,fHighX,100,0.0,10.0);
194 fMCInvmassPt=new TH2D("fMCInvmassPt"," ",fNbins,fLowX,fHighX,100,0.0,10.0);
195 fMCInvmassPtTrue=new TH2D("fMCInvmassPtTrue"," ",fNbins,fLowX,fHighX,100,0.0,10.0);
196
f9afc48d 197}
198
199//________________________________________________________________________
200void AliResonanceKink::Analyse(AliESDEvent* esd, AliMCEvent* mcEvent)
201{
202 // Main loop
203 // Called for each event
204 Int_t resonancePDGcode, antiresonancePDGcode;
205
206 if (fdaughter1pdg==kdaughterKaon) {
207 resonancePDGcode=fresonancePDGcode;
208 antiresonancePDGcode=-fresonancePDGcode;
10eaad41 209 }
f9afc48d 210 if (fdaughter1pdg!=kdaughterKaon) {
211 resonancePDGcode=-fresonancePDGcode;
212 antiresonancePDGcode=fresonancePDGcode;
213 }
214 if (fdaughter1pdg==fdaughter2pdg) {
215 resonancePDGcode=fresonancePDGcode;
216 antiresonancePDGcode=fresonancePDGcode;
217 }
10eaad41 218
f9afc48d 219 Double_t daughter1pdgMass=TDatabasePDG::Instance()->GetParticle(fdaughter1pdg)->Mass();
220 Double_t daughter2pdgMass=TDatabasePDG::Instance()->GetParticle(fdaughter2pdg)->Mass();
221
222 if (!esd) {
10eaad41 223 Printf("ERROR: fESD not available");
224 return;
225 }
226
f9afc48d 227 if (!mcEvent) {
228 Printf("ERROR: mcEvent not available");
229 return;
230 }
231
10eaad41 232 AliStack* stack=mcEvent->Stack();
233
234 if(fAnalysisType == "MC") {
7ccf0419 235
236 const AliESDVertex* vertex = GetEventVertex(esd);
237 if(!vertex) return;
238 Double_t vtx[3];
239 vertex->GetXYZ(vtx);
240 fvtxz->Fill(vtx[2]);
241
10eaad41 242 for (Int_t iMc = 0; iMc < stack->GetNprimary(); ++iMc)
243 {
244 TParticle* particle = stack->Particle(iMc);
245
246 if (!particle)
247 {
1ac5c71d 248 if (fDebug > 0) Printf("UNEXPECTED: particle with label %d not found in stack (mc loop)", iMc);
10eaad41 249 continue;
250 }
251
7ccf0419 252 if(TMath::Abs(particle->GetPdgCode())==fresonancePDGcode) {
04c3c355 253
254 if(particle->GetNDaughters()>2) continue;
7ccf0419 255
10eaad41 256 Int_t firstD=particle->GetFirstDaughter();
257 Int_t lastD=particle->GetLastDaughter();
7ccf0419 258
259 if((firstD<0)||(firstD>=stack->GetNtrack())) continue;
260 if((lastD<0)||(lastD>=stack->GetNtrack())) continue;
10eaad41 261
7ccf0419 262 TParticle *daughterParticle1;
263 TParticle *daughterParticle2;
264 AliMCParticle *mcDaughter1;
265 AliMCParticle *mcDaughter2;
266
267 if(fdaughter1pdg==kdaughterKaon) {
268 daughterParticle1=stack->Particle(firstD);
269 daughterParticle2=stack->Particle(lastD);
270 mcDaughter1= (AliMCParticle*) mcEvent->GetTrack(firstD);
271 mcDaughter2= (AliMCParticle*) mcEvent->GetTrack(lastD);
10eaad41 272 }
79b67623 273 else
7ccf0419 274 if(fdaughter2pdg==kdaughterKaon) {
275 daughterParticle1=stack->Particle(lastD);
276 daughterParticle2=stack->Particle(firstD);
277 mcDaughter1= (AliMCParticle*) mcEvent->GetTrack(lastD);
278 mcDaughter2= (AliMCParticle*) mcEvent->GetTrack(firstD);
279 } //to ensure that the first daughter is always the kaon
79b67623 280
281 if(TMath::Abs(daughterParticle1->GetPdgCode())!=321) continue;
282
7ccf0419 283 TParticle* daughters1Daughter=0;
284 TParticle* daughters2Daughter=0;
285 Int_t mcProcessDaughters1Daughter = -999;
286 Int_t mcProcessDaughters2Daughter = -999;
287 AliMCParticle *mcDaughters1Daughter;
288
289 if(mcDaughter1->Charge()==0) continue;
290 if(mcDaughter2->Charge()==0) continue; //accept resonance decays in two charged daughters
291
292 Int_t nDecayKaonDaughter=-99;
293 for(Int_t ia=0; ia<daughterParticle1->GetNDaughters(); ia++) {
294 if(((daughterParticle1->GetFirstDaughter()+ia)>0)&&((daughterParticle1->GetFirstDaughter()+ia)<stack->GetNtrack())) {
295 daughters1Daughter=stack->Particle(daughterParticle1->GetFirstDaughter()+ia);
296 mcProcessDaughters1Daughter=daughters1Daughter->GetUniqueID();
297 if(mcProcessDaughters1Daughter==4) {
298 nDecayKaonDaughter=daughterParticle1->GetFirstDaughter()+ia;
299 break;
300 }
301 }
302 }
303
304 Int_t nProcessDaughter=-99;
305 for(Int_t ib=0; ib<daughterParticle2->GetNDaughters(); ib++) {
306 if(((daughterParticle2->GetFirstDaughter()+ib)>0)&&((daughterParticle2->GetFirstDaughter()+ib)<stack->GetNtrack())) {
307 daughters2Daughter=stack->Particle(daughterParticle2->GetFirstDaughter()+ib);
308 mcProcessDaughters2Daughter=daughters2Daughter->GetUniqueID();
309 if((mcProcessDaughters2Daughter==4)||(mcProcessDaughters2Daughter==13)) {
310 nProcessDaughter=mcProcessDaughters2Daughter;
311 break;
312 }
313 }
314 }
315
316 Int_t numberOfCharged=0;
317 if((mcProcessDaughters1Daughter==4)&&(nDecayKaonDaughter>=0)) {
318 for(Int_t ic=nDecayKaonDaughter; ic<=daughterParticle1->GetLastDaughter(); ic++) {
319 if ((ic>=0)&&(ic<stack->GetNtrack())) mcDaughters1Daughter= dynamic_cast<AliMCParticle*>(mcEvent->GetTrack(ic));
320 else continue;
321 if(mcDaughters1Daughter->Charge()!=0) numberOfCharged=numberOfCharged+1;
322 }
323 }
324
325 if ((particle->Pt()>0.25)&&(TMath::Abs(particle->Eta())<0.9)) {
326
10eaad41 327 fSimEta->Fill(particle->Eta());
10eaad41 328 fSimEtaPt->Fill(particle->Pt(), particle->Eta());
7ccf0419 329 fSimPt->Fill(particle->Pt());
330 fMCInvmassPtTrue->Fill(particle->GetMass(), particle->Pt());
331 if(numberOfCharged>=2) continue; // leave out kaon decays to more than one charged daughter
332
333
334 if((daughterParticle1->Pt()>0.25)&&(TMath::Abs(daughterParticle1->Eta())<0.9)&&(daughterParticle2->Pt()>0.25)&&(TMath::Abs(daughterParticle2->Eta())<0.9)) {
335 if((mcProcessDaughters1Daughter==4)&&(daughters1Daughter->R()>120.)&&(daughters1Daughter->R()<220.)&&( (nProcessDaughter<0)||((daughters2Daughter->R()>120.)&&(nProcessDaughter>0)))) { //below are the findable
336 fSimPtKink->Fill(particle->Pt());
337 fSimEtaKink->Fill(particle->Eta());
338 fSimEtaPtKink->Fill(particle->Pt(), particle->Eta());
339 fMCInvmassPt->Fill(particle->GetMass(), particle->Pt());
340 }
341 }
342
343 } // for the generated spectra
344
345 } //for the particular resonance
346 } //MC loop
10eaad41 347
348 } // end fAnalysisType==MC
349 else
350
351 if(fAnalysisType == "ESD") {
f9afc48d 352 const AliESDVertex* vertex = GetEventVertex(esd);
10eaad41 353 if(!vertex) return;
354 Double_t vtx[3];
355 vertex->GetXYZ(vtx);
356 fvtxz->Fill(vtx[2]);
10eaad41 357 Double_t ptrackpos[3], ptrackneg[3];
358
359 TLorentzVector p4pos, anp4pos;
360 TLorentzVector p4neg, anp4neg;
361 TLorentzVector p4comb, anp4comb;
362
f9afc48d 363 for (Int_t iTracks = 0; iTracks < esd->GetNumberOfTracks(); iTracks++) {
364 AliESDtrack* trackpos = esd->GetTrack(iTracks);
10eaad41 365 if (!trackpos) {
1ac5c71d 366 if (fDebug > 0) Printf("ERROR: Could not receive track %d", iTracks);
10eaad41 367 continue;
368 }
369 if (trackpos->GetSign() < 0) continue;
370
92adf4f6 371 AliExternalTrackParam *tpcTrackpos = (AliExternalTrackParam *)trackpos->GetTPCInnerParam();
372 if(!tpcTrackpos) continue;
373 ptrackpos[0]=tpcTrackpos->Px();
374 ptrackpos[1]=tpcTrackpos->Py();
375 ptrackpos[2]=tpcTrackpos->Pz();
10eaad41 376
92adf4f6 377 Bool_t firstLevelAcceptPosTrack=IsAcceptedForKink(esd, vertex, trackpos);
378 if(firstLevelAcceptPosTrack==kFALSE) continue;
10eaad41 379
380 TVector3 posTrackMom(ptrackpos[0],ptrackpos[1],ptrackpos[2]);
92adf4f6 381
10eaad41 382 TParticle * partpos = stack->Particle(TMath::Abs(trackpos->GetLabel()));
383 if (!partpos) continue;
384 Int_t pdgpos = partpos->GetPdgCode();
385 Int_t mumlabelpos=partpos->GetFirstMother();
386 mumlabelpos = TMath::Abs(mumlabelpos);
387 TParticle * mumpos=stack->Particle(mumlabelpos);
388 if (!mumpos) continue;
389 Int_t mumpdgpos = mumpos->GetPdgCode();
390
391 Int_t indexKinkPos=trackpos->GetKinkIndex(0);
92adf4f6 392 Bool_t posKaonKinkFlag=0;
393 if(indexKinkPos<0) posKaonKinkFlag=IsKink(esd, indexKinkPos, posTrackMom);
10eaad41 394
92adf4f6 395 if(posKaonKinkFlag==1) anp4pos.SetVectM(posTrackMom,daughter1pdgMass);
10eaad41 396
397 if(indexKinkPos==0) {
92adf4f6 398
399 Bool_t secondLevelAcceptPosTrack=IsAcceptedForTrack(esd, vertex, trackpos);
400 if(secondLevelAcceptPosTrack==kFALSE) continue;
401
f9afc48d 402 p4pos.SetVectM(posTrackMom, daughter2pdgMass);
10eaad41 403
404 }
405
f9afc48d 406 for (Int_t j=0; j<esd->GetNumberOfTracks(); j++) {
10eaad41 407 if(iTracks==j) continue;
f9afc48d 408 AliESDtrack* trackneg=esd->GetTrack(j);
10eaad41 409 if (trackneg->GetSign() > 0) continue;
410
92adf4f6 411 AliExternalTrackParam *tpcTrackneg = (AliExternalTrackParam *)trackneg->GetTPCInnerParam();
412 if(!tpcTrackneg) continue;
413 ptrackneg[0]=tpcTrackneg->Px();
414 ptrackneg[1]=tpcTrackneg->Py();
415 ptrackneg[2]=tpcTrackneg->Pz();
10eaad41 416
92adf4f6 417 Bool_t firstLevelAcceptNegTrack=IsAcceptedForKink(esd, vertex, trackneg);
418 if(firstLevelAcceptNegTrack==kFALSE) continue;
10eaad41 419
420 TVector3 negTrackMom(ptrackneg[0],ptrackneg[1],ptrackneg[2]);
10eaad41 421
422 TParticle * partneg = stack->Particle(TMath::Abs(trackneg->GetLabel()));
423 if (!partneg) continue;
424 Int_t pdgneg = partneg->GetPdgCode();
425 Int_t mumlabelneg=partneg->GetFirstMother();
426 mumlabelneg = TMath::Abs(mumlabelneg);
427 TParticle * mumneg=stack->Particle(mumlabelneg);
428 if (!mumneg) continue;
429 Int_t mumpdgneg = mumneg->GetPdgCode();
430
431 Int_t indexKinkNeg=trackneg->GetKinkIndex(0);
92adf4f6 432 Bool_t negKaonKinkFlag=0;
433 if(indexKinkNeg<0) negKaonKinkFlag=IsKink(esd, indexKinkNeg, negTrackMom);
10eaad41 434
f9afc48d 435 if(negKaonKinkFlag==1) p4neg.SetVectM(negTrackMom,daughter1pdgMass);
10eaad41 436
437 if(indexKinkNeg==0) {
92adf4f6 438
439 Bool_t secondLevelAcceptNegTrack=IsAcceptedForTrack(esd, vertex, trackneg);
440 if(secondLevelAcceptNegTrack==kFALSE) continue;
441
f9afc48d 442 anp4neg.SetVectM(negTrackMom, daughter2pdgMass);
10eaad41 443
444 }
445
446 Double_t openingAngle=(ptrackpos[0]*ptrackneg[0]+ptrackpos[1]*ptrackneg[1]+ptrackpos[2]*ptrackneg[2])/(posTrackMom.Mag()*negTrackMom.Mag());
447
92adf4f6 448 if((posKaonKinkFlag==1)&&(negKaonKinkFlag==1)) {
10eaad41 449 p4comb=anp4pos;
450 p4comb+=p4neg;
451 if(openingAngle>0.6) fPhiBothKinks->Fill(p4comb.M());
452 }
453
454 if(negKaonKinkFlag==1) {
455 p4comb=p4pos;
456 p4comb+=p4neg;
457 fInvariantMass->Fill(p4comb.M());
7ccf0419 458 fInvmassPt->Fill(p4comb.M(), p4comb.Vect().Pt());
f9afc48d 459 if ((mumpdgpos==(antiresonancePDGcode))&&(mumpdgneg==(antiresonancePDGcode))&&(mumlabelpos==mumlabelneg)
460 &&(pdgpos==fdaughter2pdg)&&(pdgneg==(-fdaughter1pdg))&&(TMath::Abs(trackpos->GetLabel())>=0)&&(TMath::Abs(trackneg->GetLabel())>=0)&&(mumlabelpos>=0)&&(mumlabelneg>=0)) {
10eaad41 461 fOpeningAngle->Fill(openingAngle);
462 fInvMassTrue->Fill(p4comb.M());
7ccf0419 463 fInvmassPtTrue->Fill(p4comb.M(), p4comb.Vect().Pt());
10eaad41 464 if((TMath::Abs(p4pos.Vect().Eta())<1.1)&&(TMath::Abs(p4neg.Vect().Eta())<1.1)&&(p4comb.Vect().Eta()<1.1)) {
465 fRecPt->Fill(p4comb.Vect().Pt());
466 fRecEta->Fill(p4comb.Vect().Eta());
467 fRecEtaPt->Fill(p4comb.Vect().Perp(),p4comb.Vect().Eta());
468
469 }
470
471 }
472
473 }
474
92adf4f6 475 if(posKaonKinkFlag==1) {
10eaad41 476 anp4comb=anp4pos;
477 anp4comb+=anp4neg;
478 fInvariantMass->Fill(anp4comb.M());
7ccf0419 479 fInvmassPt->Fill(anp4comb.M(), anp4comb.Vect().Pt());
f9afc48d 480 if ((mumpdgpos==resonancePDGcode)&&(mumpdgneg==resonancePDGcode)&&(mumlabelpos==mumlabelneg)
481 &&(pdgpos==fdaughter1pdg)&&(pdgneg==(-fdaughter2pdg))&&(TMath::Abs(trackpos->GetLabel())>=0)&&(TMath::Abs(trackneg->GetLabel())>=0)&&(mumlabelpos>=0) &&(mumlabelneg>=0)) {
10eaad41 482 fOpeningAngle->Fill(openingAngle);
483 fInvMassTrue->Fill(p4comb.M());
7ccf0419 484 fInvmassPtTrue->Fill(anp4comb.M(), anp4comb.Vect().Pt());
10eaad41 485 if((TMath::Abs(anp4neg.Vect().Eta())<1.1)&&(TMath::Abs(anp4pos.Vect().Eta())<1.1)&&(anp4comb.Vect().Eta()<1.1)) {
486 fRecPt->Fill(anp4comb.Vect().Pt());
487 fRecEta->Fill(anp4comb.Vect().Eta());
488 fRecEtaPt->Fill(anp4comb.Vect().Pt(), anp4comb.Vect().Eta());
489 }
490
491 }
492
493 }
494
495 } //inner track loop
496
497 } //outer track loop
498
499 } // end fAnalysisType == ESD
500
10eaad41 501}
502
10eaad41 503//____________________________________________________________________//
f9afc48d 504Float_t AliResonanceKink::GetSigmaToVertex(AliESDtrack* esdTrack) const {
10eaad41 505 // Calculates the number of sigma to the vertex.
506
507 Float_t b[2];
508 Float_t bRes[2];
509 Float_t bCov[3];
510
92adf4f6 511 esdTrack->GetImpactParametersTPC(b,bCov);
10eaad41 512
513 if (bCov[0]<=0 || bCov[2]<=0) {
10eaad41 514 bCov[0]=0; bCov[2]=0;
515 }
f9afc48d 516
10eaad41 517 bRes[0] = TMath::Sqrt(bCov[0]);
518 bRes[1] = TMath::Sqrt(bCov[2]);
519
520 if (bRes[0] == 0 || bRes[1] ==0) return -1;
521
522 Float_t d = TMath::Sqrt(TMath::Power(b[0]/bRes[0],2) + TMath::Power(b[1]/bRes[1],2));
523
524 if (TMath::Exp(-d * d / 2) < 1e-10) return 1000;
525
526 d = TMath::ErfInverse(1 - TMath::Exp(-d * d / 2)) * TMath::Sqrt(2);
527
528 return d;
529}
530
f9afc48d 531//________________________________________________________________________
532const AliESDVertex* AliResonanceKink::GetEventVertex(const AliESDEvent* esd) const
10eaad41 533{
534 // Get the vertex
535
536 const AliESDVertex* vertex = esd->GetPrimaryVertex();
537
538 if((vertex->GetStatus()==kTRUE)&&(vertex->GetNContributors()>2)) return vertex;
539 else
540 {
541 vertex = esd->GetPrimaryVertexSPD();
542 if((vertex->GetStatus()==kTRUE)&&(vertex->GetNContributors()>2)) return vertex;
543 else
544 return 0;
545 }
546}
547
92adf4f6 548//________________________________________________________________________
549
550 Bool_t AliResonanceKink::IsAcceptedForKink(AliESDEvent *localesd,
551 const AliESDVertex *localvertex, AliESDtrack* localtrack) {
552 // Apply the selections for each kink
553
554 Double_t gPt = 0.0, gPx = 0.0, gPy = 0.0, gPz = 0.0;
555 Double_t dca[2] = {0.0,0.0}, cov[3] = {0.0,0.0,0.0}; //The impact parameters and their covariance.
556 Double_t dca3D = 0.0;
557
558 AliExternalTrackParam *tpcTrack = (AliExternalTrackParam *)localtrack->GetTPCInnerParam();
559 if(!tpcTrack) {
560 gPt = 0.0; gPx = 0.0; gPy = 0.0; gPz = 0.0;
561 dca[0] = -100.; dca[1] = -100.; dca3D = -100.;
562 cov[0] = -100.; cov[1] = -100.; cov[2] = -100.;
563 }
564 else {
565 gPt = tpcTrack->Pt();
566 gPx = tpcTrack->Px();
567 gPy = tpcTrack->Py();
568 gPz = tpcTrack->Pz();
569 tpcTrack->PropagateToDCA(localvertex,
570 localesd->GetMagneticField(),100.,dca,cov);
571 }
572
573 if(GetSigmaToVertex(localtrack) > fMaxNSigmaToVertex) {
1ac5c71d 574 if (fDebug > 1) Printf("IsAcceptedKink: Track rejected because it has a %lf sigmas to vertex TPC (max. requested: %lf)", GetSigmaToVertex(localtrack),fMaxNSigmaToVertex);
92adf4f6 575 return kFALSE;
576 }
577
578 if(TMath::Abs(dca[0]) > fMaxDCAxy) {
1ac5c71d 579 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 580 return kFALSE;
581 }
582
583 if(TMath::Abs(dca[1]) > fMaxDCAzaxis) {
1ac5c71d 584 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 585 return kFALSE;
586 }
587
588 if(gPt < fMinPtTrackCut) {
1ac5c71d 589 if (fDebug > 1) Printf("IsAcceptedKink: Track rejected because it has a min value of pt of %lf (min. requested: %lf)", gPt, fMinPtTrackCut);
92adf4f6 590 return kFALSE;
591 }
592
593 return kTRUE;
594}
595
596//________________________________________________________________________
597Bool_t AliResonanceKink::IsAcceptedForTrack(AliESDEvent *localesd, const AliESDVertex *localvertex, AliESDtrack *localtrack) {
598 // Apply the selections for each track
f9afc48d 599
92adf4f6 600 Double_t gPt = 0.0, gPx = 0.0, gPy = 0.0, gPz = 0.0;
601 Double_t dca[2] = {0.0,0.0}, cov[3] = {0.0,0.0,0.0}; //The impact parameters and their covariance.
602 Double_t dca3D = 0.0;
603
604 AliExternalTrackParam *tpcTrack = (AliExternalTrackParam *)localtrack->GetTPCInnerParam();
605 if(!tpcTrack) {
606 gPt = 0.0; gPx = 0.0; gPy = 0.0; gPz = 0.0;
607 dca[0] = -100.; dca[1] = -100.; dca3D = -100.;
608 cov[0] = -100.; cov[1] = -100.; cov[2] = -100.;
609 }
610 else {
611 gPt = tpcTrack->Pt();
612 gPx = tpcTrack->Px();
613 gPy = tpcTrack->Py();
614 gPz = tpcTrack->Pz();
615 tpcTrack->PropagateToDCA(localvertex,
616 localesd->GetMagneticField(),100.,dca,cov);
617 }
618
619 Int_t fcls[200];
620 Int_t nClustersTPC=localtrack->GetTPCclusters(fcls);
621 Float_t chi2perTPCcluster=-1.0;
622 if(nClustersTPC!=0) chi2perTPCcluster=(localtrack->GetTPCchi2())/Float_t(nClustersTPC);
623
624 Double_t extCov[15];
625 localtrack->GetExternalCovariance(extCov);
626
627 if((localtrack->GetStatus() & AliESDtrack::kTPCrefit) == 0) {
1ac5c71d 628 if (fDebug > 1) Printf("IsAccepted: Track rejected because of no refited in TPC");
92adf4f6 629 return kFALSE;
630 }
631
632 if(nClustersTPC < fMinTPCclusters) {
1ac5c71d 633 if (fDebug > 1) Printf("IsAccepted: Track rejected because it has a value of nclusters (TPC) of %ld (min. requested: %ld)", nClustersTPC, fMinTPCclusters);
92adf4f6 634 return kFALSE;
635 }
636
637 if(chi2perTPCcluster > fMaxChi2PerTPCcluster) {
1ac5c71d 638 if (fDebug > 1) Printf("IsAccepted: Track rejected because it has a value of chi2perTPCcluster of %lf (max. requested: %lf)", chi2perTPCcluster, fMaxChi2PerTPCcluster);
92adf4f6 639 return kFALSE;
640 }
641
642 if(extCov[0] > fMaxCov0) {
1ac5c71d 643 if (fDebug > 1) Printf("IsAccepted: Track rejected because it has a value of cov[0] of %lf (max. requested: %lf)", cov[0], fMaxCov0);
92adf4f6 644 return kFALSE;
645 }
646
647 if(extCov[2] > fMaxCov2) {
1ac5c71d 648 if (fDebug > 1) Printf("IsAccepted: Track rejected because it has a value of cov[2] of %lf (max. requested: %lf)", cov[2], fMaxCov2);
92adf4f6 649 return kFALSE;
650 }
651
652 if(extCov[5] > fMaxCov5) {
1ac5c71d 653 if (fDebug > 1) Printf("IsAccepted: Track rejected because it has a value of cov[5] of %lf (max. requested: %lf)", cov[5], fMaxCov5);
92adf4f6 654 return kFALSE;
655 }
656
657 if(extCov[9] > fMaxCov9) {
1ac5c71d 658 if (fDebug > 1) Printf("IsAccepted: Track rejected because it has a value of cov[9] of %lf (max. requested: %lf)", cov[9], fMaxCov9);
92adf4f6 659 return kFALSE;
660 }
661
662 if(extCov[14] > fMaxCov14) {
1ac5c71d 663 if (fDebug > 1) Printf("IsAccepted: Track rejected because it has a value of cov[14] of %lf (max. requested: %lf)", cov[14], fMaxCov14);
92adf4f6 664 return kFALSE;
665 }
f9afc48d 666
92adf4f6 667 return kTRUE;
668
669}
670
671//________________________________________________________________________
672Bool_t AliResonanceKink::IsKink(AliESDEvent *localesd, Int_t kinkIndex, TVector3 trackMom)
673{
674 // Test some kinematical criteria for each kink
675
676 AliESDkink *kink=localesd->GetKink(TMath::Abs(kinkIndex)-1);
677 const TVector3 motherMfromKink(kink->GetMotherP());
678 const TVector3 daughterMKink(kink->GetDaughterP());
679 Float_t qt=kink->GetQt();
680
681 Double_t maxDecAngKmu=f1->Eval(motherMfromKink.Mag(),0.,0.,0.);
682 Double_t maxDecAngpimu=f2->Eval(motherMfromKink.Mag(),0.,0.,0.);
683
684 Float_t kinkAngle=TMath::RadToDeg()*kink->GetAngle(2);
685
686 Float_t energyDaughterMu=TMath::Sqrt(daughterMKink.Mag()*daughterMKink.Mag()+0.105658*0.105658);
687 Float_t p1XM= motherMfromKink.Px();
688 Float_t p1YM= motherMfromKink.Py();
689 Float_t p1ZM= motherMfromKink.Pz();
690 Float_t p2XM= daughterMKink.Px();
691 Float_t p2YM= daughterMKink.Py();
692 Float_t p2ZM= daughterMKink.Pz();
693 Float_t p3Daughter=TMath::Sqrt(((p1XM-p2XM)*(p1XM-p2XM))+((p1YM-p2YM)*(p1YM-p2YM))+((p1ZM-p2ZM)*(p1ZM-p2ZM)));
694 Double_t invariantMassKmu= TMath::Sqrt((energyDaughterMu+p3Daughter)*(energyDaughterMu+p3Daughter)-motherMfromKink.Mag()*motherMfromKink.Mag());
695
696 if((kinkAngle>maxDecAngpimu)&&(qt>0.05)&&(qt<0.25)&&((kink->GetR()>110.)&&(kink->GetR()<230.))&&(TMath::Abs(trackMom.Eta())<1.1)&&(invariantMassKmu<0.6)) {
697
698 if(trackMom.Mag()<=1.1) {
699 return kTRUE;
700 }
701 else
702 if (kinkAngle<maxDecAngKmu) {
703 return kTRUE;
704 }
705 }
706 return kFALSE;
707}