Update for Kink tasks:
[u/mrichter/AliRoot.git] / PWG2 / KINK / AliResonanceKink.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 AliResonanceKink
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, fdaughter1pdg and fdaughter2pdg accordingly as well as daughter1pdgMass  and daughter2pdgMass.
20 //        Also, depending on the analysis mode (ESD or MC), fAnalysisType in the constructor must also be changed 
21 //-----------------------------------------------------------------------------------------------------------------
22
23 #include "TH2D.h"
24 #include "TParticle.h"
25 #include "TDatabasePDG.h"
26 #include "TF1.h"
27 #include "TList.h"
28 #include "AliMCEvent.h"
29 #include "AliResonanceKink.h"
30 #include "AliESDkink.h"
31 #include "AliStack.h"
32 #include "AliESDtrack.h"
33 #include "AliESDEvent.h"
34 #include "AliExternalTrackParam.h"
35 #include "AliMCParticle.h"
36
37 ClassImp(AliResonanceKink)
38
39 //________________________________________________________________________
40 AliResonanceKink::AliResonanceKink() 
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), 
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), 
43 fMinTPCclusters(0),fMaxChi2PerTPCcluster(0), fMaxCov0(0), fMaxCov2(0), fMaxCov5(0) , fMaxCov9(0), fMaxCov14(0) //, fTPCrefitFlag(kFALSE)
44 , fInvmassPt(0), fInvmassPtTrue(0), fMCInvmassPt(0), fMCInvmassPtTrue(0)
45 {
46   // Constructor
47 }
48
49 //________________________________________________________________________
50 AliResonanceKink::AliResonanceKink(Int_t nbins, Float_t nlowx, Float_t nhighx, Int_t daughter1, Int_t daughter2, Int_t resonancePDG) 
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), 
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), 
53 fMaxDCAxy(0), fMaxDCAzaxis(0), fMinTPCclusters(0), fMaxChi2PerTPCcluster(0), fMaxCov0(0), fMaxCov2(0), fMaxCov5(0), fMaxCov9(0), fMaxCov14(0) //, fTPCrefitFlag(kFALSE)
54 , fInvmassPt(0), fInvmassPtTrue(0), fMCInvmassPt(0), fMCInvmassPtTrue(0)
55 {
56    // Constructor
57   
58    fOpeningAngle=new TH1D("fOpeningAngle"," ", 100,-1.0,1.0);
59
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)
63
64    fRecPt=new TH1D("fRecPt"," ", 50,0.0,5.0);
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); 
67    fSimPt=new TH1D("fSimPt"," ", 50,0.0,5.0);
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);
70    fSimPtKink=new TH1D("fSimPtKink"," ", 50,0.0,5.0);
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);                
73    fhdr=new TH1D("fhdr"," ", 100,0.0,5.0);  
74    fhdz=new TH1D("fhdz"," ", 100,0.0,5.0);
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    
86    fvtxz=new TH1D("fvtxz"," ", 100,-20.0,20.0);
87    
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); 
92 }
93
94 //________________________________________________________________________
95 AliResonanceKink:: ~AliResonanceKink()
96 {
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;   
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;            
118 }
119 //________________________________________________________________________
120 TList* AliResonanceKink::GetHistogramList()
121 {
122   // Adding histograms to the list
123   fListOfHistos=new TList();
124  
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);
141   fListOfHistos->Add(fInvmassPt);
142   fListOfHistos->Add(fInvmassPtTrue);
143   fListOfHistos->Add(fMCInvmassPt);
144   fListOfHistos->Add(fMCInvmassPtTrue);     
145      
146   return fListOfHistos;
147 }
148
149 //________________________________________________________________________
150 void 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);
174    
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());
184    
185   fvtxz=new TH1D("fvtxz"," ", 100,-20.0,20.0);
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   
192 }
193
194 //________________________________________________________________________
195 void AliResonanceKink::Analyse(AliESDEvent* esd, AliMCEvent* mcEvent) 
196 {
197   // Main loop
198   // Called for each event
199   Int_t resonancePDGcode, antiresonancePDGcode;
200   Double_t daughter1pdgMass, daughter2pdgMass;
201   
202   if (fdaughter1pdg==kKPlus)  {
203     resonancePDGcode=fresonancePDGcode;
204     antiresonancePDGcode=-fresonancePDGcode;
205     daughter1pdgMass=TDatabasePDG::Instance()->GetParticle(fdaughter1pdg)->Mass();
206     daughter2pdgMass=TDatabasePDG::Instance()->GetParticle(fdaughter2pdg)->Mass(); 
207   }
208   
209   if (fdaughter1pdg!=kKPlus)  {
210     resonancePDGcode=-fresonancePDGcode;
211     antiresonancePDGcode=fresonancePDGcode;
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   
216   if (fdaughter1pdg==fdaughter2pdg)  {
217     resonancePDGcode=fresonancePDGcode;
218     antiresonancePDGcode=fresonancePDGcode;
219   }  
220
221   if (!esd) {
222     Printf("ERROR: fESD not available");
223     return;
224   }  
225
226     if (!mcEvent) {
227     Printf("ERROR: mcEvent not available");
228     return;
229   }  
230
231   AliStack* stack=mcEvent->Stack();
232
233   if(fAnalysisType == "MC") {
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
241   for (Int_t iMc = 0; iMc < stack->GetNprimary(); ++iMc)
242   {
243     TParticle* particle = stack->Particle(iMc);
244
245     if (!particle)
246     {
247       if (fDebug > 0) Printf("UNEXPECTED: particle with label %d not found in stack (mc loop)", iMc);
248       continue;
249     }
250
251        if(TMath::Abs(particle->GetPdgCode())==fresonancePDGcode) {
252        
253        if(particle->GetNDaughters()>2) continue;
254       
255        Int_t firstD=particle->GetFirstDaughter();
256        Int_t lastD=particle->GetLastDaughter();
257
258        if((firstD<0)||(firstD>=stack->GetNtrack())) continue;
259        if((lastD<0)||(lastD>=stack->GetNtrack())) continue;
260        
261        TParticle *daughterParticle1 = 0;
262        TParticle *daughterParticle2 = 0;
263        AliMCParticle   *mcDaughter1 = 0;
264        AliMCParticle   *mcDaughter2 = 0;       
265                
266         if(fdaughter1pdg==kKPlus) {  
267           daughterParticle1=stack->Particle(firstD);
268           daughterParticle2=stack->Particle(lastD);
269           mcDaughter1= (AliMCParticle*) mcEvent->GetTrack(firstD);  
270           mcDaughter2= (AliMCParticle*) mcEvent->GetTrack(lastD);         
271         }
272         else 
273             if(fdaughter2pdg==kKPlus) {
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
279              
280        if(TMath::Abs(daughterParticle1->GetPdgCode())!=321) continue;
281        
282        TParticle      * daughters1Daughter=0;
283        TParticle      * daughters2Daughter=0;       
284        Int_t mcProcessDaughters1Daughter = -999;
285        Int_t mcProcessDaughters2Daughter = -999;       
286        AliMCParticle *mcDaughters1Daughter = 0;
287        
288        if(mcDaughter1->Charge()==0) continue;
289        if(mcDaughter2->Charge()==0) continue;       //accept resonance decays in two charged daughters
290
291        Int_t nDecayKaonDaughter=-99; 
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        
303        Int_t nProcessDaughter=-99; 
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        
315        Int_t numberOfCharged=0;
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        
324          if(numberOfCharged>=2) continue; // leave out kaon decays to more than one charged daughter
325          
326          if ((particle->Pt()>0.25)&&(TMath::Abs(particle->Eta())<0.9)) {
327
328          fSimEta->Fill(particle->Eta());
329          fSimEtaPt->Fill(particle->Pt(), particle->Eta());
330          fSimPt->Fill(particle->Pt());
331          fMCInvmassPtTrue->Fill(particle->GetMass(), particle->Pt());      
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
346   
347   } // end fAnalysisType==MC
348   else 
349   
350   if(fAnalysisType == "ESD") {
351   const AliESDVertex* vertex = GetEventVertex(esd);
352   if(!vertex) return;
353   Double_t vtx[3];
354   vertex->GetXYZ(vtx);
355   fvtxz->Fill(vtx[2]);
356   Double_t ptrackpos[3], ptrackneg[3];
357   
358   TLorentzVector p4pos, anp4pos;
359   TLorentzVector p4neg, anp4neg;
360   TLorentzVector p4comb, anp4comb;
361   
362   for (Int_t iTracks = 0; iTracks < esd->GetNumberOfTracks(); iTracks++) {
363     AliESDtrack* trackpos = esd->GetTrack(iTracks);
364     if (!trackpos) {
365       if (fDebug > 0) Printf("ERROR: Could not receive track %d", iTracks);
366       continue;
367     }
368     if (trackpos->GetSign() < 0) continue;
369     
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();  
375     
376     Bool_t firstLevelAcceptPosTrack=IsAcceptedForKink(esd, vertex, trackpos);
377     if(firstLevelAcceptPosTrack==kFALSE) continue;
378     
379     TVector3 posTrackMom(ptrackpos[0],ptrackpos[1],ptrackpos[2]);
380         
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);
391     
392     if(indexKinkPos>0) continue;
393     
394     Bool_t posKaonKinkFlag=0;
395     
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     }
402     
403     if(indexKinkPos==0) {
404
405     Bool_t secondLevelAcceptPosTrack=IsAcceptedForTrack(esd, vertex, trackpos);
406     if(secondLevelAcceptPosTrack==kFALSE) continue;
407
408       p4pos.SetVectM(posTrackMom, daughter2pdgMass);
409     
410     }
411         
412       for (Int_t j=0; j<esd->GetNumberOfTracks(); j++) {
413         if(iTracks==j) continue;
414         AliESDtrack* trackneg=esd->GetTrack(j);
415         if (trackneg->GetSign() > 0) continue;
416         
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();  
422     
423         Bool_t firstLevelAcceptNegTrack=IsAcceptedForKink(esd, vertex, trackneg);
424         if(firstLevelAcceptNegTrack==kFALSE) continue;  
425
426         TVector3 negTrackMom(ptrackneg[0],ptrackneg[1],ptrackneg[2]);
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);
438         
439         if(indexKinkNeg>0) continue;
440         
441         Bool_t negKaonKinkFlag=0;
442         
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         }
449         
450         if(indexKinkNeg==0)  {
451  
452            Bool_t secondLevelAcceptNegTrack=IsAcceptedForTrack(esd, vertex, trackneg);
453            if(secondLevelAcceptNegTrack==kFALSE) continue;  
454           
455           anp4neg.SetVectM(negTrackMom, daughter2pdgMass);
456         
457         }
458         
459         Double_t openingAngle=(ptrackpos[0]*ptrackneg[0]+ptrackpos[1]*ptrackneg[1]+ptrackpos[2]*ptrackneg[2])/(posTrackMom.Mag()*negTrackMom.Mag());
460
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           }
467         }
468                 
469         if(negKaonKinkFlag==1) {
470           p4comb=p4pos;
471           p4comb+=p4neg;
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
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         
494         if(posKaonKinkFlag==1) {
495           anp4comb=anp4pos;
496           anp4comb+=anp4neg;  
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());
512            }
513
514          }
515
516         }
517          
518       } //inner track loop 
519
520   } //outer track loop 
521   
522   } // end fAnalysisType == ESD
523   
524 }      
525
526 //____________________________________________________________________//
527 Float_t AliResonanceKink::GetSigmaToVertex(AliESDtrack* esdTrack) const {
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
534     esdTrack->GetImpactParametersTPC(b,bCov);
535   
536   if (bCov[0]<=0 || bCov[2]<=0) {
537     bCov[0]=0; bCov[2]=0;
538   }
539   
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
554 //________________________________________________________________________
555 const AliESDVertex* AliResonanceKink::GetEventVertex(const AliESDEvent* esd) const
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
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) {
597       if (fDebug > 1) Printf("IsAcceptedKink: Track rejected because it has a %lf sigmas to vertex TPC (max. requested: %lf)",   GetSigmaToVertex(localtrack),fMaxNSigmaToVertex);
598       return kFALSE;
599   }
600   
601   if(TMath::Abs(dca[0]) > fMaxDCAxy) {
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);
603       return kFALSE;
604   }
605     
606   if(TMath::Abs(dca[1]) > fMaxDCAzaxis) {
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);
608       return kFALSE;
609   }
610   
611   if(gPt <= fMinPtTrackCut) {
612       if (fDebug > 1) Printf("IsAcceptedKink: Track rejected because it has a min value of pt of %lf (min. requested: %lf)", gPt, fMinPtTrackCut);
613       return kFALSE;
614   } 
615   
616   return kTRUE;
617 }
618
619 //________________________________________________________________________
620 Bool_t AliResonanceKink::IsAcceptedForTrack(AliESDEvent *localesd,                                                                                                                                          const AliESDVertex *localvertex, AliESDtrack *localtrack) {
621    // Apply the selections for each track
622
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) {
651       if (fDebug > 1) Printf("IsAccepted: Track rejected because of no refited in TPC");
652       return kFALSE;
653   } 
654
655   if(nClustersTPC < fMinTPCclusters) {
656       if (fDebug > 1) Printf("IsAccepted: Track rejected because it has a value of nclusters (TPC) of %ld (min. requested: %ld)", nClustersTPC, fMinTPCclusters);
657       return kFALSE;
658   } 
659   
660   if(chi2perTPCcluster > fMaxChi2PerTPCcluster) {
661       if (fDebug > 1) Printf("IsAccepted: Track rejected because it has a value of chi2perTPCcluster of %lf (max. requested: %lf)", chi2perTPCcluster, fMaxChi2PerTPCcluster);
662       return kFALSE;
663   } 
664
665   if(extCov[0] > fMaxCov0) {
666       if (fDebug > 1) Printf("IsAccepted: Track rejected because it has a value of cov[0] of %lf (max. requested: %lf)", cov[0], fMaxCov0);
667       return kFALSE;
668   }
669   
670   if(extCov[2] > fMaxCov2) {
671       if (fDebug > 1) Printf("IsAccepted: Track rejected because it has a value of cov[2] of %lf (max. requested: %lf)", cov[2], fMaxCov2);
672       return kFALSE;
673   }
674     
675   if(extCov[5] > fMaxCov5) {
676       if (fDebug > 1) Printf("IsAccepted: Track rejected because it has a value of cov[5] of %lf (max. requested: %lf)", cov[5], fMaxCov5);
677       return kFALSE;
678   }  
679   
680   if(extCov[9] > fMaxCov9) {
681       if (fDebug > 1) Printf("IsAccepted: Track rejected because it has a value of cov[9] of %lf (max. requested: %lf)", cov[9], fMaxCov9);
682       return kFALSE;
683   }  
684   
685   if(extCov[14] > fMaxCov14) {
686       if (fDebug > 1) Printf("IsAccepted: Track rejected because it has a value of cov[14] of %lf (max. requested: %lf)", cov[14], fMaxCov14);
687       return kFALSE;
688   } 
689  
690   return kTRUE;
691
692 }
693
694 //________________________________________________________________________
695 Bool_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
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)) {
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 }