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