Update for Kink tasks:
[u/mrichter/AliRoot.git] / PWG2 / KINK / AliResonanceKinkLikeSign.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 AliResonanceKinkLikeSign
16// Example of an analysis task for producing a like-sign background for resonances having at least one
17// kaon-kink in their decay products.
7ccf0419 18// Background is calculated from a positive kaon kink and a negative track.
10eaad41 19//-----------------------------------------------------------------------------------------------------------------
20
be1a7181 21#include "AliESDEvent.h"
10eaad41 22#include "TH1D.h"
7ccf0419 23#include "TH2D.h"
be1a7181 24#include "TF1.h"
25#include "TDatabasePDG.h"
10eaad41 26#include "TLorentzVector.h"
be1a7181 27#include "AliESDkink.h"
10eaad41 28
29#include "AliResonanceKinkLikeSign.h"
10eaad41 30
31ClassImp(AliResonanceKinkLikeSign)
32
10eaad41 33//________________________________________________________________________
34AliResonanceKinkLikeSign::AliResonanceKinkLikeSign(const char *name)
be1a7181 35 : AliAnalysisTaskSE(name), fDebug(0), fListOfHistos(0), f1(0), f2(0), fPosKaonLikeSign(0), fLikeSignInvmassPt(0), fMaxNSigmaToVertex(0), fMinPtTrackCut(0), fMaxDCAxy(0), fMaxDCAzaxis(0),
36fMinTPCclusters(0),fMaxChi2PerTPCcluster(0), fMaxCov0(0), fMaxCov2(0), fMaxCov5(0) , fMaxCov9(0), fMaxCov14(0), fdaughter1pdg(0), fdaughter2pdg(0), fnbins(0), fnlowx(0), fnhighx(0)
10eaad41 37
38{
39 // Constructor
40
41 // Define input and output slots here
42 // Input slot #0 works with a TChain
10eaad41 43 DefineOutput(1, TList::Class());
44}
45
46//________________________________________________________________________
be1a7181 47void AliResonanceKinkLikeSign::UserCreateOutputObjects()
10eaad41 48{
49 // Create histograms
50 // Called once
51
52 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);
53 f1->SetParameter(0,0.493677);
54 f1->SetParameter(1,0.9127037);
55 f1->SetParameter(2,TMath::Pi());
56
57 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);
58 f2->SetParameter(0,0.13957018);
59 f2->SetParameter(1,0.2731374);
60 f2->SetParameter(2,TMath::Pi());
61
be1a7181 62 fPosKaonLikeSign=new TH1D("fPosKaonLikeSign"," ", fnbins, fnlowx, fnhighx);
63 fLikeSignInvmassPt=new TH2D("fLikeSignInvmassPt"," ", fnbins, fnlowx, fnhighx, 100,0.0,10.0);
64
10eaad41 65 fListOfHistos=new TList();
7ccf0419 66 fListOfHistos->Add(fPosKaonLikeSign);
67 fListOfHistos->Add(fLikeSignInvmassPt);
10eaad41 68
69}
70
71//________________________________________________________________________
be1a7181 72void AliResonanceKinkLikeSign::UserExec(Option_t *)
10eaad41 73{
74 // Main loop
75 // Called for each event
76
be1a7181 77 AliVEvent *event = InputEvent();
78 if (!event) {
79 Printf("ERROR: Could not retrieve event");
80 return;
81 }
10eaad41 82
be1a7181 83 AliESDEvent* esd = dynamic_cast<AliESDEvent*>(event);
84 if (!esd) {
85 Printf("ERROR: Could not retrieve esd");
86 return;
87 }
f2a16059 88
be1a7181 89 AliMCEvent* mcEvent = MCEvent();
90 if (!mcEvent) {
91 Printf("ERROR: Could not retrieve MC event");
92 return;
93 }
94
95 Double_t daughter1Mass, daughter2Mass;
96
97 if (fdaughter1pdg==kKPlus) {
98 daughter1Mass=TDatabasePDG::Instance()->GetParticle(fdaughter1pdg)->Mass();
99 daughter2Mass=TDatabasePDG::Instance()->GetParticle(fdaughter2pdg)->Mass();
100 }
101
102 if (fdaughter1pdg!=kKPlus) {
103 daughter1Mass=TDatabasePDG::Instance()->GetParticle(fdaughter2pdg)->Mass();
104 daughter2Mass=TDatabasePDG::Instance()->GetParticle(fdaughter1pdg)->Mass();
105 } // to ensure that daughter1Mass has always the kaon mass
106
107 if (!esd) {
108 Printf("ERROR: esd not available");
10eaad41 109 return;
110 }
111
be1a7181 112 const AliESDVertex* vertex = GetEventVertex(esd);
10eaad41 113 if(!vertex) return;
114
115 Double_t ptrackpos[3], ptrackneg[3];
116
117 TLorentzVector p4pos;
118 TLorentzVector p4neg;
119 TLorentzVector p4comb;
120
be1a7181 121 for (Int_t iTracks = 0; iTracks < esd->GetNumberOfTracks(); iTracks++) {
122
123 AliESDtrack* trackpos = esd->GetTrack(iTracks);
10eaad41 124 if (!trackpos) {
1ac5c71d 125 if (fDebug > 0) Printf("ERROR: Could not receive track %d", iTracks);
10eaad41 126 continue;
127 }
7ccf0419 128 if (trackpos->GetSign() < 0) continue;
10eaad41 129
be1a7181 130 AliExternalTrackParam *tpcTrackpos = (AliExternalTrackParam *)trackpos->GetTPCInnerParam();
131 if(!tpcTrackpos) continue;
132 ptrackpos[0]=tpcTrackpos->Px();
133 ptrackpos[1]=tpcTrackpos->Py();
134 ptrackpos[2]=tpcTrackpos->Pz();
10eaad41 135
be1a7181 136 Bool_t firstLevelAcceptPosTrack=IsAcceptedForKink(esd, vertex, trackpos);
137 if(firstLevelAcceptPosTrack==kFALSE) continue;
10eaad41 138
139 TVector3 posTrackMom(ptrackpos[0],ptrackpos[1],ptrackpos[2]);
10eaad41 140
be1a7181 141 Int_t indexKinkPos=trackpos->GetKinkIndex(0);
142 Bool_t posKaonKinkFlag=0;
143 if(indexKinkPos<0) posKaonKinkFlag=IsKink(esd, indexKinkPos, posTrackMom);
144 if(posKaonKinkFlag==0) continue;
145 if(posKaonKinkFlag==1) p4pos.SetVectM(posTrackMom, daughter1Mass);
7ccf0419 146
be1a7181 147 for (Int_t j=0; j<esd->GetNumberOfTracks(); j++) {
10eaad41 148 if(iTracks==j) continue;
be1a7181 149 AliESDtrack* trackneg=esd->GetTrack(j);
150 if (!trackneg) {
151 if (fDebug > 0) Printf("ERROR: Could not receive track %d", j);
152 continue;
153 }
7ccf0419 154 if (trackneg->GetSign() < 0) continue;
10eaad41 155
be1a7181 156 AliExternalTrackParam *tpcTrackneg = (AliExternalTrackParam *)trackneg->GetTPCInnerParam();
157 if(!tpcTrackneg) continue;
158 ptrackneg[0]=tpcTrackneg->Px();
159 ptrackneg[1]=tpcTrackneg->Py();
160 ptrackneg[2]=tpcTrackneg->Pz();
161
162 Bool_t firstLevelAcceptNegTrack=IsAcceptedForKink(esd, vertex, trackneg);
163 if(firstLevelAcceptNegTrack==kFALSE) continue;
10eaad41 164
be1a7181 165 TVector3 negTrackMom(ptrackneg[0],ptrackneg[1],ptrackneg[2]);
10eaad41 166
be1a7181 167 Bool_t secondLevelAcceptNegTrack=IsAcceptedForTrack(esd, vertex, trackneg);
168 if(secondLevelAcceptNegTrack==kFALSE) continue;
169
170 p4neg.SetVectM(negTrackMom, daughter2Mass);
10eaad41 171
172 p4comb=p4pos;
173 p4comb+=p4neg;
be1a7181 174
175 if(p4comb.Vect().Pt()<=0.25) continue;
176
177 if((TMath::Abs(p4pos.Vect().Eta())<0.9)&&(TMath::Abs(p4neg.Vect().Eta())<0.9)&&(p4comb.Vect().Eta()<0.9)) {
178
179 fPosKaonLikeSign->Fill(p4comb.M());
180 fLikeSignInvmassPt->Fill(p4comb.M(), p4comb.Vect().Pt());
181
182 }
10eaad41 183
184 } //inner track loop
185
186 } //outer track loop
187
188 PostData(1, fListOfHistos);
189}
190
be1a7181 191//________________________________________________________________________
192void AliResonanceKinkLikeSign::Terminate(Option_t *)
193{
194
195}
196
10eaad41 197//____________________________________________________________________//
198
199Float_t AliResonanceKinkLikeSign::GetSigmaToVertex(AliESDtrack* esdTrack) const {
200 // Calculates the number of sigma to the vertex.
201
202 Float_t b[2];
203 Float_t bRes[2];
204 Float_t bCov[3];
205
206 esdTrack->GetImpactParameters(b,bCov);
207
208 if (bCov[0]<=0 || bCov[2]<=0) {
209 //AliDebug(1, "Estimated b resolution lower or equal zero!");
210 bCov[0]=0; bCov[2]=0;
211 }
212 bRes[0] = TMath::Sqrt(bCov[0]);
213 bRes[1] = TMath::Sqrt(bCov[2]);
214
215 if (bRes[0] == 0 || bRes[1] ==0) return -1;
216
217 Float_t d = TMath::Sqrt(TMath::Power(b[0]/bRes[0],2) + TMath::Power(b[1]/bRes[1],2));
218
219 if (TMath::Exp(-d * d / 2) < 1e-10) return 1000;
220
221 d = TMath::ErfInverse(1 - TMath::Exp(-d * d / 2)) * TMath::Sqrt(2);
222
223 return d;
224}
225
226//____________________________________________________________________//
227
228const AliESDVertex* AliResonanceKinkLikeSign::GetEventVertex(const AliESDEvent* esd) const
229
230{
231 // Get the vertex
232
233 const AliESDVertex* vertex = esd->GetPrimaryVertex();
234
235 if((vertex->GetStatus()==kTRUE)&&(vertex->GetNContributors()>2)) return vertex;
236 else
237 {
238 vertex = esd->GetPrimaryVertexSPD();
239 if((vertex->GetStatus()==kTRUE)&&(vertex->GetNContributors()>2)) return vertex;
240 else
241 return 0;
242 }
243}
244
be1a7181 245//________________________________________________________________________
246
247 Bool_t AliResonanceKinkLikeSign::IsAcceptedForKink(AliESDEvent *localesd,
248 const AliESDVertex *localvertex, AliESDtrack* localtrack) {
249 // Apply the selections for each kink
250
251 Double_t gPt = 0.0, gPx = 0.0, gPy = 0.0, gPz = 0.0;
252 Double_t dca[2] = {0.0,0.0}, cov[3] = {0.0,0.0,0.0}; //The impact parameters and their covariance.
253 Double_t dca3D = 0.0;
254
255 AliExternalTrackParam *tpcTrack = (AliExternalTrackParam *)localtrack->GetTPCInnerParam();
256 if(!tpcTrack) {
257 gPt = 0.0; gPx = 0.0; gPy = 0.0; gPz = 0.0;
258 dca[0] = -100.; dca[1] = -100.; dca3D = -100.;
259 cov[0] = -100.; cov[1] = -100.; cov[2] = -100.;
260 }
261 else {
262 gPt = tpcTrack->Pt();
263 gPx = tpcTrack->Px();
264 gPy = tpcTrack->Py();
265 gPz = tpcTrack->Pz();
266 tpcTrack->PropagateToDCA(localvertex,
267 localesd->GetMagneticField(),100.,dca,cov);
268 }
269
270 if(GetSigmaToVertex(localtrack) > fMaxNSigmaToVertex) {
271 if (fDebug > 1) Printf("IsAcceptedKink: Track rejected because it has a %lf sigmas to vertex TPC (max. requested: %lf)", GetSigmaToVertex(localtrack),fMaxNSigmaToVertex);
272 return kFALSE;
273 }
274
275 if(TMath::Abs(dca[0]) > fMaxDCAxy) {
276 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);
277 return kFALSE;
278 }
279
280 if(TMath::Abs(dca[1]) > fMaxDCAzaxis) {
281 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);
282 return kFALSE;
283 }
284
285 if(gPt <= fMinPtTrackCut) {
286 if (fDebug > 1) Printf("IsAcceptedKink: Track rejected because it has a min value of pt of %lf (min. requested: %lf)", gPt, fMinPtTrackCut);
287 return kFALSE;
288 }
289
290 return kTRUE;
291}
292
293//________________________________________________________________________
294Bool_t AliResonanceKinkLikeSign::IsAcceptedForTrack(AliESDEvent *localesd, const AliESDVertex *localvertex, AliESDtrack *localtrack) {
295 // Apply the selections for each track
296
297 Double_t gPt = 0.0, gPx = 0.0, gPy = 0.0, gPz = 0.0;
298 Double_t dca[2] = {0.0,0.0}, cov[3] = {0.0,0.0,0.0}; //The impact parameters and their covariance.
299 Double_t dca3D = 0.0;
300
301 AliExternalTrackParam *tpcTrack = (AliExternalTrackParam *)localtrack->GetTPCInnerParam();
302 if(!tpcTrack) {
303 gPt = 0.0; gPx = 0.0; gPy = 0.0; gPz = 0.0;
304 dca[0] = -100.; dca[1] = -100.; dca3D = -100.;
305 cov[0] = -100.; cov[1] = -100.; cov[2] = -100.;
306 }
307 else {
308 gPt = tpcTrack->Pt();
309 gPx = tpcTrack->Px();
310 gPy = tpcTrack->Py();
311 gPz = tpcTrack->Pz();
312 tpcTrack->PropagateToDCA(localvertex,
313 localesd->GetMagneticField(),100.,dca,cov);
314 }
315
316 Int_t fcls[200];
317 Int_t nClustersTPC=localtrack->GetTPCclusters(fcls);
318 Float_t chi2perTPCcluster=-1.0;
319 if(nClustersTPC!=0) chi2perTPCcluster=(localtrack->GetTPCchi2())/Float_t(nClustersTPC);
320
321 Double_t extCov[15];
322 localtrack->GetExternalCovariance(extCov);
323
324 if((localtrack->GetStatus() & AliESDtrack::kTPCrefit) == 0) {
325 if (fDebug > 1) Printf("IsAccepted: Track rejected because of no refited in TPC");
326 return kFALSE;
327 }
328
329 if(nClustersTPC < fMinTPCclusters) {
330 if (fDebug > 1) Printf("IsAccepted: Track rejected because it has a value of nclusters (TPC) of %ld (min. requested: %ld)", nClustersTPC, fMinTPCclusters);
331 return kFALSE;
332 }
333
334 if(chi2perTPCcluster > fMaxChi2PerTPCcluster) {
335 if (fDebug > 1) Printf("IsAccepted: Track rejected because it has a value of chi2perTPCcluster of %lf (max. requested: %lf)", chi2perTPCcluster, fMaxChi2PerTPCcluster);
336 return kFALSE;
337 }
338
339 if(extCov[0] > fMaxCov0) {
340 if (fDebug > 1) Printf("IsAccepted: Track rejected because it has a value of cov[0] of %lf (max. requested: %lf)", cov[0], fMaxCov0);
341 return kFALSE;
342 }
343
344 if(extCov[2] > fMaxCov2) {
345 if (fDebug > 1) Printf("IsAccepted: Track rejected because it has a value of cov[2] of %lf (max. requested: %lf)", cov[2], fMaxCov2);
346 return kFALSE;
347 }
348
349 if(extCov[5] > fMaxCov5) {
350 if (fDebug > 1) Printf("IsAccepted: Track rejected because it has a value of cov[5] of %lf (max. requested: %lf)", cov[5], fMaxCov5);
351 return kFALSE;
352 }
353
354 if(extCov[9] > fMaxCov9) {
355 if (fDebug > 1) Printf("IsAccepted: Track rejected because it has a value of cov[9] of %lf (max. requested: %lf)", cov[9], fMaxCov9);
356 return kFALSE;
357 }
358
359 if(extCov[14] > fMaxCov14) {
360 if (fDebug > 1) Printf("IsAccepted: Track rejected because it has a value of cov[14] of %lf (max. requested: %lf)", cov[14], fMaxCov14);
361 return kFALSE;
362 }
363
364 return kTRUE;
365
366}
367
368//________________________________________________________________________
369Bool_t AliResonanceKinkLikeSign::IsKink(AliESDEvent *localesd, Int_t kinkIndex, TVector3 trackMom)
370{
371 // Test some kinematical criteria for each kink
372
373 AliESDkink *kink=localesd->GetKink(TMath::Abs(kinkIndex)-1);
374 const TVector3 motherMfromKink(kink->GetMotherP());
375 const TVector3 daughterMKink(kink->GetDaughterP());
376 Float_t qt=kink->GetQt();
377
378 Double_t maxDecAngKmu=f1->Eval(motherMfromKink.Mag(),0.,0.,0.);
379 Double_t maxDecAngpimu=f2->Eval(motherMfromKink.Mag(),0.,0.,0.);
380
381 Float_t kinkAngle=TMath::RadToDeg()*kink->GetAngle(2);
382
383 Float_t energyDaughterMu=TMath::Sqrt(daughterMKink.Mag()*daughterMKink.Mag()+0.105658*0.105658);
384 Float_t p1XM= motherMfromKink.Px();
385 Float_t p1YM= motherMfromKink.Py();
386 Float_t p1ZM= motherMfromKink.Pz();
387 Float_t p2XM= daughterMKink.Px();
388 Float_t p2YM= daughterMKink.Py();
389 Float_t p2ZM= daughterMKink.Pz();
390 Float_t p3Daughter=TMath::Sqrt(((p1XM-p2XM)*(p1XM-p2XM))+((p1YM-p2YM)*(p1YM-p2YM))+((p1ZM-p2ZM)*(p1ZM-p2ZM)));
391 Double_t invariantMassKmu= TMath::Sqrt((energyDaughterMu+p3Daughter)*(energyDaughterMu+p3Daughter)-motherMfromKink.Mag()*motherMfromKink.Mag());
392
393 if((kinkAngle>maxDecAngpimu)&&(qt>0.05)&&(qt<0.25)&&((kink->GetR()>120.)&&(kink->GetR()<220.))&&(TMath::Abs(trackMom.Eta())<0.9)&&(invariantMassKmu<0.6)) {
394
395 if(trackMom.Mag()<=1.1) {
396 return kTRUE;
397 }
398 else
399 if (kinkAngle<maxDecAngKmu) {
400 return kTRUE;
401 }
402 }
403 return kFALSE;
404}