]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG2/KINK/AliResonanceKinkPID.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Particle mass access updated: Evi Ganoti, University of Athens (pganoti@phys.uoa.gr)
[u/mrichter/AliRoot.git] / PWG2 / KINK / AliResonanceKinkPID.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//----------------------------------------------------------------------------------------------------------------
15// class AliResonanceKinkPID
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
19// variables resonancePDG, daughter1 and daughter2 accordingly as well as daughter1Mass and daughter2Mass.
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"
27#include <TVector3.h>
f27d6e67 28#include <TDatabasePDG.h>
29#include <TParticlePDG.h>
10eaad41 30#include "TF1.h"
31#include "AliAnalysisTaskSE.h"
32#include "AliAnalysisManager.h"
33
34#include "AliESDInputHandler.h"
35
36#include "AliMCEventHandler.h"
37#include "AliMCEvent.h"
38
39#include "AliResonanceKinkPID.h"
40#include "AliESDkink.h"
41#include "AliStack.h"
42
43ClassImp(AliResonanceKinkPID)
44
45//________________________________________________________________________
46AliResonanceKinkPID::AliResonanceKinkPID()
47 : AliAnalysisTaskSE(), fESD(0), fListOfHistos(0), fOpeningAngle(0), fInvariantMass(0), fInvMassTrue(0), fPhiBothKinks(0), fRecPt(0), fRecEta(0), fRecEtaPt(0),
48 fSimPt(0), fSimEta(0), fSimEtaPt(0), fSimPtKink(0), fSimEtaKink(0), fSimEtaPtKink(0),
49 fhdr(0), fhdz(0), f1(0), f2(0), fAnalysisType("ESD"), fvtxz(0)
50
51{
52 // Constructor
53}
54
55//________________________________________________________________________
56AliResonanceKinkPID::AliResonanceKinkPID(const char *name)
57 : AliAnalysisTaskSE(name), fESD(0), fListOfHistos(0), fOpeningAngle(0), fInvariantMass(0), fInvMassTrue(0), fPhiBothKinks(0), fRecPt(0), fRecEta(0), fRecEtaPt(0),
58 fSimPt(0), fSimEta(0), fSimEtaPt(0), fSimPtKink(0), fSimEtaKink(0), fSimEtaPtKink(0),
59 fhdr(0), fhdz(0), f1(0), f2(0), fAnalysisType("ESD"), fvtxz(0)
60
61{
62 // Constructor
63
64 // Define input and output slots here
65 // Input slot #0 works with a TChain
66 DefineInput(0, TChain::Class());
67 DefineOutput(1, TList::Class());
68}
69
70//________________________________________________________________________
71void AliResonanceKinkPID::ConnectInputData(Option_t *)
72{
73 // Connect ESD or AOD here
74 // Called once
10eaad41 75
76 TTree* tree = dynamic_cast<TTree*> (GetInputData(0));
77 if (!tree) {
78 Printf("ERROR: Could not read chain from input slot 0");
79 } else {
10eaad41 80
81 AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
82
83 if (!esdH) {
84 Printf("ERROR: Could not get ESDInputHandler");
85 } else
86 fESD = esdH->GetEvent();
87 }
88}
89
90//________________________________________________________________________
91void AliResonanceKinkPID::CreateOutputObjects()
92{
93 // Create histograms
94 // Called once
95
96 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);
97 f1->SetParameter(0,0.493677);
98 f1->SetParameter(1,0.9127037);
99 f1->SetParameter(2,TMath::Pi());
100
101 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);
102 f2->SetParameter(0,0.13957018);
103 f2->SetParameter(1,0.2731374);
104 f2->SetParameter(2,TMath::Pi());
105
106 // OpenFile(1); Uncomment for proof analysis
107 fOpeningAngle=new TH1D("fOpeningAngle"," ", 100,-1.0,1.0);
108
109 // for K*0(892)
110 fInvariantMass=new TH1D("fInvariantMass"," ", 60,0.6,1.2);
111 fInvMassTrue=new TH1D("fInvMassTrue"," ",60,0.6,1.2);
112 fPhiBothKinks=new TH1D("fPhiBothKinks"," ", 60,0.6,1.2); // Not applicable for K*0
113
114 //for phi(1020)
115 //fInvariantMass=new TH1D("fInvariantMass"," ", 70,0.99,1.088);
116 //fInvMassTrue=new TH1D("fInvMassTrue"," ", 70,0.99,1.088);
117 //fPhiBothKinks=new TH1D("fPhiBothKinks"," ", 70,0.99,1.088);
118
119 fRecPt=new TH1D("fRecPt"," ", 50,0.0,5.0);
120 fRecEta=new TH1D("fRecEta"," ", 44,-1.1,1.1);
121 fRecEtaPt=new TH2D("fRecEtaPt"," ", 50,0.0,5.0, 44,-1.1,1.1);
122 fSimPt=new TH1D("fSimPt"," ", 50,0.0,5.0);
123 fSimEta=new TH1D("fSimEta"," ", 44,-1.1,1.1);
124 fSimEtaPt=new TH2D("fSimEtaPt"," ", 50,0.0,5.0, 44,-1.1,1.1);
125 fSimPtKink=new TH1D("fSimPtKink"," ", 50,0.0,5.0);
126 fSimEtaKink=new TH1D("fSimEtaKink"," ", 44,-1.1,1.1);
127 fSimEtaPtKink=new TH2D("fSimEtaPtKink"," ", 50,0.0,5.0, 44,-1.1,1.1);
128 fhdr=new TH1D("fhdr"," ", 100,0.0,5.0);
129 fhdz=new TH1D("fhdz"," ", 100,0.0,5.0);
130 fvtxz=new TH1D("fvtxz"," ", 100,-20.0,20.0);
131
132 fListOfHistos=new TList();
133
134 fListOfHistos->Add(fOpeningAngle);
135 fListOfHistos->Add(fInvariantMass);
136 fListOfHistos->Add(fInvMassTrue);
137 fListOfHistos->Add(fPhiBothKinks);
138 fListOfHistos->Add(fRecPt);
139 fListOfHistos->Add(fRecEta);
140 fListOfHistos->Add(fRecEtaPt);
141 fListOfHistos->Add(fSimPt);
142 fListOfHistos->Add(fSimEta);
143 fListOfHistos->Add(fSimEtaPt);
144 fListOfHistos->Add(fSimPtKink);
145 fListOfHistos->Add(fSimEtaKink);
146 fListOfHistos->Add(fSimEtaPtKink);
147 fListOfHistos->Add(fhdr);
148 fListOfHistos->Add(fhdz);
149 fListOfHistos->Add(fvtxz);
150
151}
152
153//________________________________________________________________________
154void AliResonanceKinkPID::Exec(Option_t *)
155{
156 // Main loop
157 // Called for each event
158
159 enum ResonanceType {kPhi=333, kKstar0=313, kLambda1520=3124};
160 enum DaughterType {kdaughterPion=211, kdaughterKaon=321, kdaughterProton=2212};
161
162 Int_t resonancePDG=kKstar0;
163
164 Int_t daughter1=kdaughterKaon;
165 Int_t daughter2=kdaughterPion;
166
f27d6e67 167 Double_t daughter1Mass=TDatabasePDG::Instance()->GetParticle(kKPlus)->Mass();
168 Double_t daughter2Mass=TDatabasePDG::Instance()->GetParticle(kPiPlus)->Mass();
169
10eaad41 170 AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
171 if (!eventHandler) {
172 Printf("ERROR: Could not retrieve MC event handler");
173 return;
174 }
175
176 AliMCEvent* mcEvent = eventHandler->MCEvent();
177 if (!mcEvent) {
178 Printf("ERROR: Could not retrieve MC event");
179 return;
180 }
181
182 Printf("MC particles: %d", mcEvent->GetNumberOfTracks());
183
184 if (!fESD) {
185 Printf("ERROR: fESD not available");
186 return;
187 }
188
189 AliStack* stack=mcEvent->Stack();
190
191 if(fAnalysisType == "MC") {
192
193 for (Int_t iMc = 0; iMc < stack->GetNprimary(); ++iMc)
194 {
195 TParticle* particle = stack->Particle(iMc);
196
197 if (!particle)
198 {
199 Printf("UNEXPECTED: particle with label %d not found in stack (mc loop)", iMc);
200 continue;
201 }
202
203 if(TMath::Abs(particle->GetPdgCode())==resonancePDG) {
204 Int_t firstD=particle->GetFirstDaughter();
205 Int_t lastD=particle->GetLastDaughter();
206 TParticle *daughterParticle1=stack->Particle(firstD);
207 TParticle *daughterParticle2=stack->Particle(lastD);
208
209 TParticle* kaonFirstDaughter;
210 Int_t mcProcessKaonFirstDaughter;
211
212 for(Int_t ia=0; ia<daughterParticle1->GetNDaughters(); ia++){
213 if ((daughterParticle1->GetFirstDaughter()+ia)!=-1) {
214 kaonFirstDaughter=stack->Particle(daughterParticle1->GetFirstDaughter()+ia);
215 mcProcessKaonFirstDaughter=kaonFirstDaughter->GetUniqueID();
216 }
217 }
218
219 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)) {
220 fSimEta->Fill(particle->Eta());
221 fSimPt->Fill(particle->Pt());
222 fSimEtaPt->Fill(particle->Pt(), particle->Eta());
223 if(mcProcessKaonFirstDaughter==4) {
224 fSimPtKink->Fill(particle->Pt());
225 fSimEtaKink->Fill(particle->Eta());
226 fSimEtaPtKink->Fill(particle->Pt(), particle->Eta());
227 }
228 }
229 }
230 }
231
232 } // end fAnalysisType==MC
233 else
234
235 if(fAnalysisType == "ESD") {
236
237 const AliESDVertex* vertex = GetEventVertex(fESD);
238 if(!vertex) return;
239 Double_t vtx[3];
240 vertex->GetXYZ(vtx);
241 fvtxz->Fill(vtx[2]);
242
243 Double_t ptrackpos[3], ptrackneg[3];
244
245 TLorentzVector p4pos, anp4pos;
246 TLorentzVector p4neg, anp4neg;
247 TLorentzVector p4comb, anp4comb;
248
249
250 for (Int_t iTracks = 0; iTracks < fESD->GetNumberOfTracks(); iTracks++) {
251 AliESDtrack* trackpos = fESD->GetTrack(iTracks);
252 if (!trackpos) {
253 Printf("ERROR: Could not receive track %d", iTracks);
254 continue;
255 }
256 if (trackpos->GetSign() < 0) continue;
257
258 trackpos->GetPxPyPz(ptrackpos);
259
260 Float_t nSigmaToVertex = GetSigmaToVertex(trackpos);
261
262 Float_t bpos[2];
263 Float_t bCovpos[3];
264 trackpos->GetImpactParameters(bpos,bCovpos);
265
266 if (bCovpos[0]<=0 || bCovpos[2]<=0) {
267 Printf("Estimated b resolution lower or equal zero!");
268 bCovpos[0]=0; bCovpos[2]=0;
269 }
270
271 Float_t dcaToVertexXYpos = bpos[0];
272 Float_t dcaToVertexZpos = bpos[1];
273
274 fhdr->Fill(dcaToVertexXYpos);
275 fhdz->Fill(dcaToVertexZpos);
276
277 if(nSigmaToVertex>=4) continue;
278 if((dcaToVertexXYpos>3.0)||(dcaToVertexZpos>3.0)) continue;
279
280 TVector3 posTrackMom(ptrackpos[0],ptrackpos[1],ptrackpos[2]);
281
282 if(posTrackMom.Perp()<=0.25) continue;
283
284 TParticle * partpos = stack->Particle(TMath::Abs(trackpos->GetLabel()));
285 if (!partpos) continue;
286 Int_t pdgpos = partpos->GetPdgCode();
287 Int_t mumlabelpos=partpos->GetFirstMother();
288 mumlabelpos = TMath::Abs(mumlabelpos);
289 TParticle * mumpos=stack->Particle(mumlabelpos);
290 if (!mumpos) continue;
291 Int_t mumpdgpos = mumpos->GetPdgCode();
292
293 Int_t indexKinkPos=trackpos->GetKinkIndex(0);
294 Int_t kaonKinkFlag=0;
295 if(indexKinkPos<0){
296
297 AliESDkink *poskink=fESD->GetKink(TMath::Abs(indexKinkPos)-1);
298 const TVector3 motherMfromKinkPos(poskink->GetMotherP());
299 const TVector3 daughterMKinkPos(poskink->GetDaughterP());
300 Float_t posQt=poskink->GetQt();
301
302 Double_t maxDecAngKmuPos=f1->Eval(motherMfromKinkPos.Mag(),0.,0.,0.);
303 Double_t maxDecAngpimuPos=f2->Eval(motherMfromKinkPos.Mag(),0.,0.,0.);
304
305 Float_t kinkAnglePos=TMath::RadToDeg()*poskink->GetAngle(2);
306
307 Float_t energyDaughterMu=TMath::Sqrt(daughterMKinkPos.Mag()*daughterMKinkPos.Mag()+0.105658*0.105658);
308 Float_t p1XM= motherMfromKinkPos.Px();
309 Float_t p1YM= motherMfromKinkPos.Py();
310 Float_t p1ZM= motherMfromKinkPos.Pz();
311 Float_t p2XM= daughterMKinkPos.Px();
312 Float_t p2YM= daughterMKinkPos.Py();
313 Float_t p2ZM= daughterMKinkPos.Pz();
314 Float_t p3Daughter=TMath::Sqrt(((p1XM-p2XM)*(p1XM-p2XM))+((p1YM-p2YM)*(p1YM-p2YM))+((p1ZM-p2ZM)*(p1ZM-p2ZM)));
315 Double_t invariantMassKmuPos= TMath::Sqrt((energyDaughterMu+p3Daughter)*(energyDaughterMu+p3Daughter)-motherMfromKinkPos.Mag()*motherMfromKinkPos.Mag());
316
317 if((kinkAnglePos>maxDecAngpimuPos)&&(posQt>0.05)&&(posQt<0.25)&&((poskink->GetR()>110.)&&(poskink->GetR()<230.))&&(TMath::Abs(posTrackMom.Eta())<1.1)&&(invariantMassKmuPos<0.6)) {
318
319 if(posTrackMom.Mag()<=1.1) {
320 kaonKinkFlag=1;
321 }
322 else
323 if (kinkAnglePos<maxDecAngKmuPos) {
324 kaonKinkFlag=1;
325 }
326 }
327
328 } //End Kink Information
329
330 if(kaonKinkFlag==1) anp4pos.SetVectM(posTrackMom,daughter1Mass);
331
332 if(indexKinkPos==0) {
333 UInt_t status=trackpos->GetStatus();
334 if((status&AliESDtrack::kTPCrefit)==0) continue;
335 if(trackpos->GetTPCclusters(0)<50) continue;
336 if((trackpos->GetTPCchi2()/trackpos->GetTPCclusters(0))>3.5) continue;
337 Double_t extCovPos[15];
338 trackpos->GetExternalCovariance(extCovPos);
339 if(extCovPos[0]>2) continue;
340 if(extCovPos[2]>2) continue;
341 if(extCovPos[5]>0.5) continue;
342 if(extCovPos[9]>0.5) continue;
343 if(extCovPos[14]>2) continue;
344
345 p4pos.SetVectM(posTrackMom, daughter2Mass);
346
347 }
348
349 for (Int_t j=0; j<fESD->GetNumberOfTracks(); j++) {
350 if(iTracks==j) continue;
351 AliESDtrack* trackneg=fESD->GetTrack(j);
352 if (trackneg->GetSign() > 0) continue;
353
354 trackneg->GetPxPyPz(ptrackneg);
355 Float_t negSigmaToVertex = GetSigmaToVertex(trackneg);
356
357 Float_t bneg[2];
358 Float_t bCovneg[3];
359 trackneg->GetImpactParameters(bneg,bCovneg);
360 if (bCovneg[0]<=0 || bCovneg[2]<=0) {
361 Printf("Estimated b resolution lower or equal zero!");
362 bCovneg[0]=0; bCovneg[2]=0;
363 }
364
365 Float_t dcaToVertexXYneg = bneg[0];
366 Float_t dcaToVertexZneg = bneg[1];
367
368 fhdr->Fill(dcaToVertexXYneg);
369 fhdz->Fill(dcaToVertexZneg);
370
371 if(negSigmaToVertex>=4) continue;
372 if((dcaToVertexXYneg>3.0)||(dcaToVertexZneg>3.0)) continue;
373
374 TVector3 negTrackMom(ptrackneg[0],ptrackneg[1],ptrackneg[2]);
375
376 if(negTrackMom.Perp()<=0.25) continue;
377
378 TParticle * partneg = stack->Particle(TMath::Abs(trackneg->GetLabel()));
379 if (!partneg) continue;
380 Int_t pdgneg = partneg->GetPdgCode();
381 Int_t mumlabelneg=partneg->GetFirstMother();
382 mumlabelneg = TMath::Abs(mumlabelneg);
383 TParticle * mumneg=stack->Particle(mumlabelneg);
384 if (!mumneg) continue;
385 Int_t mumpdgneg = mumneg->GetPdgCode();
386
387 Int_t indexKinkNeg=trackneg->GetKinkIndex(0);
388 Int_t negKaonKinkFlag=0;
389 if(indexKinkNeg<0){
390
391 AliESDkink *negkink=fESD->GetKink(TMath::Abs(indexKinkNeg)-1);
392 const TVector3 motherMfromKinkNeg(negkink->GetMotherP());
393 const TVector3 daughterMKinkNeg(negkink->GetDaughterP());
394 Float_t negQt=negkink->GetQt();
395
396 Double_t maxDecAngKmuNeg=f1->Eval(motherMfromKinkNeg.Mag(),0.,0.,0.);
397 Double_t maxDecAngpimuNeg=f2->Eval(motherMfromKinkNeg.Mag(),0.,0.,0.);
398
399 Float_t kinkAngleNeg=TMath::RadToDeg()*negkink->GetAngle(2);
400
401 Float_t energyDaughterMuNeg=TMath::Sqrt(daughterMKinkNeg.Mag()*daughterMKinkNeg.Mag()+0.105658*0.105658);
402 Float_t p1XMNeg= motherMfromKinkNeg.Px();
403 Float_t p1YMNeg= motherMfromKinkNeg.Py();
404 Float_t p1ZMNeg= motherMfromKinkNeg.Pz();
405 Float_t p2XMNeg= daughterMKinkNeg.Px();
406 Float_t p2YMNeg= daughterMKinkNeg.Py();
407 Float_t p2ZMNeg= daughterMKinkNeg.Pz();
408 Float_t p3DaughterNeg=TMath::Sqrt(((p1XMNeg-p2XMNeg)*(p1XMNeg-p2XMNeg))+((p1YMNeg-p2YMNeg)*(p1YMNeg-p2YMNeg))+((p1ZMNeg-p2ZMNeg)*(p1ZMNeg-p2ZMNeg)));
409 Double_t invariantMassKmuNeg= TMath::Sqrt((energyDaughterMuNeg+p3DaughterNeg)*(energyDaughterMuNeg+p3DaughterNeg)-motherMfromKinkNeg.Mag()*motherMfromKinkNeg.Mag());
410
411 if((kinkAngleNeg>maxDecAngpimuNeg)&&(negQt>0.05)&&(negQt<0.25)&&((negkink->GetR()>110.)&&(negkink->GetR()<230.))&&(TMath::Abs(negTrackMom.Eta())<1.1)&&(invariantMassKmuNeg<0.6)) {
412
413 if(negTrackMom.Mag()<=1.1) {
414 negKaonKinkFlag=1;
415 }
416 else
417 if (kinkAngleNeg<maxDecAngKmuNeg) {
418 negKaonKinkFlag=1;
419 }
420 }
421
422 } //End Kink Information
423
424 if(negKaonKinkFlag==1) p4neg.SetVectM(negTrackMom,daughter1Mass);
425
426 if(indexKinkNeg==0) {
427 UInt_t statusneg=trackneg->GetStatus();
428
429 if((statusneg&AliESDtrack::kTPCrefit)==0) continue;
430
431 if(trackneg->GetTPCclusters(0)<50) continue;
432 if((trackneg->GetTPCchi2()/trackneg->GetTPCclusters(0))>3.5) continue;
433 Double_t extCovneg[15];
434 trackneg->GetExternalCovariance(extCovneg);
435 if(extCovneg[0]>2) continue;
436 if(extCovneg[2]>2) continue;
437 if(extCovneg[5]>0.5) continue;
438 if(extCovneg[9]>0.5) continue;
439 if(extCovneg[14]>2) continue;
440
441 anp4neg.SetVectM(negTrackMom, daughter2Mass);
442
443 }
444
445 Double_t openingAngle=(ptrackpos[0]*ptrackneg[0]+ptrackpos[1]*ptrackneg[1]+ptrackpos[2]*ptrackneg[2])/(posTrackMom.Mag()*negTrackMom.Mag());
446
447 if((kaonKinkFlag==1)&&(negKaonKinkFlag==1)) {
448 p4comb=anp4pos;
449 p4comb+=p4neg;
450 if(openingAngle>0.6) fPhiBothKinks->Fill(p4comb.M());
451 }
452
453 if(negKaonKinkFlag==1) {
454 p4comb=p4pos;
455 p4comb+=p4neg;
456 fInvariantMass->Fill(p4comb.M());
457 if ((mumpdgpos==(-resonancePDG))&&(mumpdgneg==(-resonancePDG))&&(mumlabelpos==mumlabelneg)
458 &&(pdgpos==daughter2)&&(pdgneg==(-daughter1))&&(TMath::Abs(trackpos->GetLabel())>=0)&&(TMath::Abs(trackneg->GetLabel())>=0)&&(mumlabelpos>=0)&&(mumlabelneg>=0)) {
459 fOpeningAngle->Fill(openingAngle);
460 fInvMassTrue->Fill(p4comb.M());
461 if((TMath::Abs(p4pos.Vect().Eta())<1.1)&&(TMath::Abs(p4neg.Vect().Eta())<1.1)&&(p4comb.Vect().Eta()<1.1)) {
462 fRecPt->Fill(p4comb.Vect().Pt());
463 fRecEta->Fill(p4comb.Vect().Eta());
464 fRecEtaPt->Fill(p4comb.Vect().Perp(),p4comb.Vect().Eta());
465
466 }
467
468 }
469
470 }
471
472 if(kaonKinkFlag==1) {
473 anp4comb=anp4pos;
474 anp4comb+=anp4neg;
475 fInvariantMass->Fill(anp4comb.M());
476 if ((mumpdgpos==resonancePDG)&&(mumpdgneg==resonancePDG)&&(mumlabelpos==mumlabelneg)
477 &&(pdgpos==daughter1)&&(pdgneg==(-daughter2))&&(TMath::Abs(trackpos->GetLabel())>=0)&&(TMath::Abs(trackneg->GetLabel())>=0)&&(mumlabelpos>=0) &&(mumlabelneg>=0)) {
478 fOpeningAngle->Fill(openingAngle);
479 fInvMassTrue->Fill(p4comb.M());
480 if((TMath::Abs(anp4neg.Vect().Eta())<1.1)&&(TMath::Abs(anp4pos.Vect().Eta())<1.1)&&(anp4comb.Vect().Eta()<1.1)) {
481 fRecPt->Fill(anp4comb.Vect().Pt());
482 fRecEta->Fill(anp4comb.Vect().Eta());
483 fRecEtaPt->Fill(anp4comb.Vect().Pt(), anp4comb.Vect().Eta());
484 }
485
486 }
487
488 }
489
490 } //inner track loop
491
492 } //outer track loop
493
494 } // end fAnalysisType == ESD
495
496 PostData(1, fListOfHistos);
497}
498
499//________________________________________________________________________
500void AliResonanceKinkPID::Terminate(Option_t *)
501{
502 // Draw result to the screen
503 // Called once at the end of the query
504
505// fHistPt = dynamic_cast<TH1F*> (GetOutputData(0));
506// if (!fHistPt) {
507// Printf("ERROR: fHistPt not available");
508// return;
509// }
510
511// TCanvas *c1 = new TCanvas("AliResonanceKinkPID","Pt MC",10,10,510,510);
512// c1->cd(1)->SetLogy();
513// fHistPt->DrawCopy("E");
514}
515
516//____________________________________________________________________//
517
518Float_t AliResonanceKinkPID::GetSigmaToVertex(AliESDtrack* esdTrack) const {
519 // Calculates the number of sigma to the vertex.
520
521 Float_t b[2];
522 Float_t bRes[2];
523 Float_t bCov[3];
524
525 esdTrack->GetImpactParameters(b,bCov);
526
527 if (bCov[0]<=0 || bCov[2]<=0) {
528 //AliDebug(1, "Estimated b resolution lower or equal zero!");
529 bCov[0]=0; bCov[2]=0;
530 }
531 bRes[0] = TMath::Sqrt(bCov[0]);
532 bRes[1] = TMath::Sqrt(bCov[2]);
533
534 if (bRes[0] == 0 || bRes[1] ==0) return -1;
535
536 Float_t d = TMath::Sqrt(TMath::Power(b[0]/bRes[0],2) + TMath::Power(b[1]/bRes[1],2));
537
538 if (TMath::Exp(-d * d / 2) < 1e-10) return 1000;
539
540 d = TMath::ErfInverse(1 - TMath::Exp(-d * d / 2)) * TMath::Sqrt(2);
541
542 return d;
543}
544
545//____________________________________________________________________//
546
547const AliESDVertex* AliResonanceKinkPID::GetEventVertex(const AliESDEvent* esd) const
548{
549 // Get the vertex
550
551 const AliESDVertex* vertex = esd->GetPrimaryVertex();
552
553 if((vertex->GetStatus()==kTRUE)&&(vertex->GetNContributors()>2)) return vertex;
554 else
555 {
556 vertex = esd->GetPrimaryVertexSPD();
557 if((vertex->GetStatus()==kTRUE)&&(vertex->GetNContributors()>2)) return vertex;
558 else
559 return 0;
560 }
561}
562