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