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