Rulechecker-complying update from P.Ganoti (pganoti@phys.uoa.gr)
[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 "TChain.h"
24 #include "TTree.h"
25 #include "TH2D.h"
26 #include "TParticle.h"
27 #include "TDatabasePDG.h"
28 #include "TParticlePDG.h"
29 #include "TF1.h"
30 #include "TList.h"
31 #include "TString.h"
32 #include "AliMCEventHandler.h"
33 #include "AliMCEvent.h"
34 #include "AliResonanceKink.h"
35 #include "AliESDkink.h"
36 #include "AliStack.h"
37 #include "AliESDtrack.h"
38 #include "AliESDEvent.h"
39 #include "AliExternalTrackParam.h"
40 #include "AliMCParticle.h"
41
42 ClassImp(AliResonanceKink)
43
44 //________________________________________________________________________
45 AliResonanceKink::AliResonanceKink() 
46   : 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), 
47   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), 
48 fMinTPCclusters(0),fMaxChi2PerTPCcluster(0), fMaxCov0(0), fMaxCov2(0), fMaxCov5(0) , fMaxCov9(0), fMaxCov14(0) //, fTPCrefitFlag(kFALSE)
49 , fInvmassPt(0), fInvmassPtTrue(0), fMCInvmassPt(0), fMCInvmassPtTrue(0)
50 {
51   // Constructor
52 }
53
54 //________________________________________________________________________
55 AliResonanceKink::AliResonanceKink(Int_t nbins, Float_t nlowx, Float_t nhighx, Int_t daughter1, Int_t daughter2, Int_t resonancePDG) 
56   : 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), 
57   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), 
58 fMaxDCAxy(0), fMaxDCAzaxis(0), fMinTPCclusters(0), fMaxChi2PerTPCcluster(0), fMaxCov0(0), fMaxCov2(0), fMaxCov5(0), fMaxCov9(0), fMaxCov14(0) //, fTPCrefitFlag(kFALSE)
59 , fInvmassPt(0), fInvmassPtTrue(0), fMCInvmassPt(0), fMCInvmassPtTrue(0)
60 {
61    // Constructor
62   
63    fOpeningAngle=new TH1D("fOpeningAngle"," ", 100,-1.0,1.0);
64
65    fInvariantMass=new TH1D("fInvariantMass"," ",fNbins,fLowX,fHighX);
66    fInvMassTrue=new TH1D("fInvMassTrue"," ",fNbins,fLowX,fHighX);
67    fPhiBothKinks=new TH1D("fPhiBothKinks"," ",fNbins,fLowX,fHighX);  // Applicable for phi(1020)
68
69    fRecPt=new TH1D("fRecPt"," ", 50,0.0,5.0);
70    fRecEta=new TH1D("fRecEta"," ", 44,-1.1,1.1);
71    fRecEtaPt=new TH2D("fRecEtaPt"," ", 50,0.0,5.0, 44,-1.1,1.1); 
72    fSimPt=new TH1D("fSimPt"," ", 50,0.0,5.0);
73    fSimEta=new TH1D("fSimEta"," ", 44,-1.1,1.1); 
74    fSimEtaPt=new TH2D("fSimEtaPt"," ", 50,0.0,5.0, 44,-1.1,1.1);
75    fSimPtKink=new TH1D("fSimPtKink"," ", 50,0.0,5.0);
76    fSimEtaKink=new TH1D("fSimEtaKink"," ", 44,-1.1,1.1);
77    fSimEtaPtKink=new TH2D("fSimEtaPtKink"," ", 50,0.0,5.0, 44,-1.1,1.1);                
78    fhdr=new TH1D("fhdr"," ", 100,0.0,5.0);  
79    fhdz=new TH1D("fhdz"," ", 100,0.0,5.0);
80    
81    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);
82    f1->SetParameter(0,0.493677);
83    f1->SetParameter(1,0.9127037);
84    f1->SetParameter(2,TMath::Pi());
85
86    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);
87    f2->SetParameter(0,0.13957018);
88    f2->SetParameter(1,0.2731374);
89    f2->SetParameter(2,TMath::Pi());
90    
91    fvtxz=new TH1D("fvtxz"," ", 100,-20.0,20.0);
92    
93    fInvmassPt=new TH2D("fInvmassPt"," ",fNbins,fLowX,fHighX,100,0.0,10.0);
94    fInvmassPtTrue=new TH2D("fInvmassPtTrue"," ",fNbins,fLowX,fHighX,100,0.0,10.0);  
95    fMCInvmassPt=new TH2D("fMCInvmassPt"," ",fNbins,fLowX,fHighX,100,0.0,10.0);
96    fMCInvmassPtTrue=new TH2D("fMCInvmassPtTrue"," ",fNbins,fLowX,fHighX,100,0.0,10.0); 
97 }
98
99 //________________________________________________________________________
100 AliResonanceKink:: ~AliResonanceKink()
101 {
102  //  Destructor
103  if(fOpeningAngle) delete fOpeningAngle;
104  if(fInvariantMass) delete fInvariantMass;
105  if(fInvMassTrue) delete fInvMassTrue;
106  if(fPhiBothKinks) delete fPhiBothKinks;
107  if(fRecPt) delete fRecPt;
108  if(fRecEta) delete fRecEta;
109  if(fRecEtaPt) delete fRecEtaPt;
110  if(fSimPt) delete fSimPt;
111  if(fSimEta) delete fSimEta;
112  if(fSimEtaPt) delete fSimEtaPt;
113  if(fSimPtKink) delete fSimPtKink;
114  if(fSimEtaKink) delete fSimEtaKink;
115  if(fSimEtaPtKink) delete fSimEtaPtKink;
116  if(fhdr) delete fhdr;
117  if(fhdz) delete fhdz;   
118  if(fvtxz) delete fvtxz;  
119  if(fInvmassPt) delete fInvmassPt; 
120  if(fInvmassPtTrue) delete fInvmassPtTrue; 
121  if(fMCInvmassPt) delete fMCInvmassPt; 
122  if(fMCInvmassPtTrue) delete fMCInvmassPtTrue;            
123 }
124 //________________________________________________________________________
125 TList* AliResonanceKink::GetHistogramList()
126 {
127   // Adding histograms to the list
128   fListOfHistos=new TList();
129  
130   fListOfHistos->Add(fOpeningAngle);
131   fListOfHistos->Add(fInvariantMass);
132   fListOfHistos->Add(fInvMassTrue);
133   fListOfHistos->Add(fPhiBothKinks);
134   fListOfHistos->Add(fRecPt);    
135   fListOfHistos->Add(fRecEta);   
136   fListOfHistos->Add(fRecEtaPt);    
137   fListOfHistos->Add(fSimPt);    
138   fListOfHistos->Add(fSimEta);   
139   fListOfHistos->Add(fSimEtaPt);     
140   fListOfHistos->Add(fSimPtKink);    
141   fListOfHistos->Add(fSimEtaKink);   
142   fListOfHistos->Add(fSimEtaPtKink);                                                           
143   fListOfHistos->Add(fhdr);
144   fListOfHistos->Add(fhdz);
145   fListOfHistos->Add(fvtxz);
146   fListOfHistos->Add(fInvmassPt);
147   fListOfHistos->Add(fInvmassPtTrue);
148   fListOfHistos->Add(fMCInvmassPt);
149   fListOfHistos->Add(fMCInvmassPtTrue);     
150      
151   return fListOfHistos;
152 }
153
154 //________________________________________________________________________
155 void AliResonanceKink::InitOutputHistograms(Int_t nbins, Float_t nlowx, Float_t nhighx)
156 {
157   //  Initialisation of the output histograms
158   fNbins=nbins; 
159   fLowX=nlowx; 
160   fHighX=nhighx;
161     
162   fOpeningAngle=new TH1D("fOpeningAngle"," ", 100,-1.0,1.0);
163
164   fInvariantMass=new TH1D("fInvariantMass"," ",fNbins,fLowX,fHighX);
165   fInvMassTrue=new TH1D("fInvMassTrue"," ",fNbins,fLowX,fHighX);
166   fPhiBothKinks=new TH1D("fPhiBothKinks"," ",fNbins,fLowX,fHighX);  // Applicable for phi(1020)
167
168   fRecPt=new TH1D("fRecPt"," ", 50,0.0,5.0);
169   fRecEta=new TH1D("fRecEta"," ", 44,-1.1,1.1);
170   fRecEtaPt=new TH2D("fRecEtaPt"," ", 50,0.0,5.0, 44,-1.1,1.1); 
171   fSimPt=new TH1D("fSimPt"," ", 50,0.0,5.0);
172   fSimEta=new TH1D("fSimEta"," ", 44,-1.1,1.1); 
173   fSimEtaPt=new TH2D("fSimEtaPt"," ", 50,0.0,5.0, 44,-1.1,1.1);
174   fSimPtKink=new TH1D("fSimPtKink"," ", 50,0.0,5.0);
175   fSimEtaKink=new TH1D("fSimEtaKink"," ", 44,-1.1,1.1);
176   fSimEtaPtKink=new TH2D("fSimEtaPtKink"," ", 50,0.0,5.0, 44,-1.1,1.1);                
177   fhdr=new TH1D("fhdr"," ", 100,0.0,5.0);  
178   fhdz=new TH1D("fhdz"," ", 100,0.0,5.0);
179    
180   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);
181   f1->SetParameter(0,0.493677);
182   f1->SetParameter(1,0.9127037);
183   f1->SetParameter(2,TMath::Pi());
184
185   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);
186   f2->SetParameter(0,0.13957018);
187   f2->SetParameter(1,0.2731374);
188   f2->SetParameter(2,TMath::Pi());
189    
190   fvtxz=new TH1D("fvtxz"," ", 100,-20.0,20.0);
191   
192   fInvmassPt=new TH2D("fInvmassPt"," ",fNbins,fLowX,fHighX,100,0.0,10.0);
193   fInvmassPtTrue=new TH2D("fInvmassPtTrue"," ",fNbins,fLowX,fHighX,100,0.0,10.0);  
194   fMCInvmassPt=new TH2D("fMCInvmassPt"," ",fNbins,fLowX,fHighX,100,0.0,10.0);
195   fMCInvmassPtTrue=new TH2D("fMCInvmassPtTrue"," ",fNbins,fLowX,fHighX,100,0.0,10.0);
196   
197 }
198   
199 //________________________________________________________________________
200 void AliResonanceKink::Analyse(AliESDEvent* esd, AliMCEvent* mcEvent) 
201 {
202   // Main loop
203   // Called for each event
204   Int_t resonancePDGcode, antiresonancePDGcode;
205   
206   if (fdaughter1pdg==kdaughterKaon)  {
207     resonancePDGcode=fresonancePDGcode;
208     antiresonancePDGcode=-fresonancePDGcode;
209   }
210   if (fdaughter1pdg!=kdaughterKaon)  {
211     resonancePDGcode=-fresonancePDGcode;
212     antiresonancePDGcode=fresonancePDGcode;
213   }  
214   if (fdaughter1pdg==fdaughter2pdg)  {
215     resonancePDGcode=fresonancePDGcode;
216     antiresonancePDGcode=fresonancePDGcode;
217   }  
218
219    Double_t daughter1pdgMass=TDatabasePDG::Instance()->GetParticle(fdaughter1pdg)->Mass();
220    Double_t daughter2pdgMass=TDatabasePDG::Instance()->GetParticle(fdaughter2pdg)->Mass();
221    
222   if (!esd) {
223     Printf("ERROR: fESD not available");
224     return;
225   }  
226
227     if (!mcEvent) {
228     Printf("ERROR: mcEvent not available");
229     return;
230   }  
231
232   AliStack* stack=mcEvent->Stack();
233
234   if(fAnalysisType == "MC") {
235
236   const AliESDVertex* vertex = GetEventVertex(esd);
237   if(!vertex) return;
238   Double_t vtx[3];
239   vertex->GetXYZ(vtx);
240   fvtxz->Fill(vtx[2]);
241
242   for (Int_t iMc = 0; iMc < stack->GetNprimary(); ++iMc)
243   {
244     TParticle* particle = stack->Particle(iMc);
245
246     if (!particle)
247     {
248       if (fDebug > 0) Printf("UNEXPECTED: particle with label %d not found in stack (mc loop)", iMc);
249       continue;
250     }
251
252        if(TMath::Abs(particle->GetPdgCode())==fresonancePDGcode) {
253        
254        if(particle->GetNDaughters()>2) continue;
255       
256        Int_t firstD=particle->GetFirstDaughter();
257        Int_t lastD=particle->GetLastDaughter();
258
259        if((firstD<0)||(firstD>=stack->GetNtrack())) continue;
260        if((lastD<0)||(lastD>=stack->GetNtrack())) continue;
261        
262        TParticle *daughterParticle1;
263        TParticle *daughterParticle2;
264        AliMCParticle *mcDaughter1;
265        AliMCParticle *mcDaughter2;       
266                
267         if(fdaughter1pdg==kdaughterKaon) {  
268           daughterParticle1=stack->Particle(firstD);
269           daughterParticle2=stack->Particle(lastD);
270           mcDaughter1= (AliMCParticle*) mcEvent->GetTrack(firstD);  
271           mcDaughter2= (AliMCParticle*) mcEvent->GetTrack(lastD);         
272         }
273         else {
274             if(fdaughter2pdg==kdaughterKaon) {
275               daughterParticle1=stack->Particle(lastD);
276               daughterParticle2=stack->Particle(firstD); 
277               mcDaughter1= (AliMCParticle*) mcEvent->GetTrack(lastD);
278               mcDaughter2= (AliMCParticle*) mcEvent->GetTrack(firstD); 
279              }    //to ensure that the first daughter is always the kaon
280             else continue;
281         }      
282  
283        TParticle* daughters1Daughter=0;
284        TParticle* daughters2Daughter=0;       
285        Int_t mcProcessDaughters1Daughter = -999;
286        Int_t mcProcessDaughters2Daughter = -999;       
287        AliMCParticle *mcDaughters1Daughter;
288        
289        if(mcDaughter1->Charge()==0) continue;
290        if(mcDaughter2->Charge()==0) continue;       //accept resonance decays in two charged daughters
291
292        Int_t nDecayKaonDaughter=-99; 
293        for(Int_t ia=0; ia<daughterParticle1->GetNDaughters(); ia++) {
294        if(((daughterParticle1->GetFirstDaughter()+ia)>0)&&((daughterParticle1->GetFirstDaughter()+ia)<stack->GetNtrack())) {
295            daughters1Daughter=stack->Particle(daughterParticle1->GetFirstDaughter()+ia);
296            mcProcessDaughters1Daughter=daughters1Daughter->GetUniqueID(); 
297            if(mcProcessDaughters1Daughter==4) {
298              nDecayKaonDaughter=daughterParticle1->GetFirstDaughter()+ia;
299              break;
300            }
301         }
302        }
303        
304        Int_t nProcessDaughter=-99; 
305        for(Int_t ib=0; ib<daughterParticle2->GetNDaughters(); ib++) {
306         if(((daughterParticle2->GetFirstDaughter()+ib)>0)&&((daughterParticle2->GetFirstDaughter()+ib)<stack->GetNtrack())) {
307            daughters2Daughter=stack->Particle(daughterParticle2->GetFirstDaughter()+ib);
308            mcProcessDaughters2Daughter=daughters2Daughter->GetUniqueID(); 
309            if((mcProcessDaughters2Daughter==4)||(mcProcessDaughters2Daughter==13)) {
310              nProcessDaughter=mcProcessDaughters2Daughter;     
311              break;
312            }  
313         }
314        }  
315        
316        Int_t numberOfCharged=0;
317        if((mcProcessDaughters1Daughter==4)&&(nDecayKaonDaughter>=0)) {
318          for(Int_t ic=nDecayKaonDaughter; ic<=daughterParticle1->GetLastDaughter(); ic++) {
319            if ((ic>=0)&&(ic<stack->GetNtrack())) mcDaughters1Daughter= dynamic_cast<AliMCParticle*>(mcEvent->GetTrack(ic));
320             else continue;
321             if(mcDaughters1Daughter->Charge()!=0) numberOfCharged=numberOfCharged+1;
322          }
323        }
324        
325          if ((particle->Pt()>0.25)&&(TMath::Abs(particle->Eta())<0.9)) {
326
327          fSimEta->Fill(particle->Eta());
328          fSimEtaPt->Fill(particle->Pt(), particle->Eta());
329          fSimPt->Fill(particle->Pt());
330          fMCInvmassPtTrue->Fill(particle->GetMass(), particle->Pt());
331          if(numberOfCharged>=2) continue; // leave out kaon decays to more than one charged daughter
332            
333
334          if((daughterParticle1->Pt()>0.25)&&(TMath::Abs(daughterParticle1->Eta())<0.9)&&(daughterParticle2->Pt()>0.25)&&(TMath::Abs(daughterParticle2->Eta())<0.9)) {
335            if((mcProcessDaughters1Daughter==4)&&(daughters1Daughter->R()>120.)&&(daughters1Daughter->R()<220.)&&( (nProcessDaughter<0)||((daughters2Daughter->R()>120.)&&(nProcessDaughter>0)))) { //below are the findable
336              fSimPtKink->Fill(particle->Pt());
337              fSimEtaKink->Fill(particle->Eta());
338              fSimEtaPtKink->Fill(particle->Pt(), particle->Eta());
339              fMCInvmassPt->Fill(particle->GetMass(), particle->Pt());
340            }
341          }
342        
343        } // for the generated spectra 
344              
345           }   //for the particular resonance
346   } //MC loop
347   
348   } // end fAnalysisType==MC
349   else 
350   
351   if(fAnalysisType == "ESD") {
352   const AliESDVertex* vertex = GetEventVertex(esd);
353   if(!vertex) return;
354   Double_t vtx[3];
355   vertex->GetXYZ(vtx);
356   fvtxz->Fill(vtx[2]);
357   Double_t ptrackpos[3], ptrackneg[3];
358   
359   TLorentzVector p4pos, anp4pos;
360   TLorentzVector p4neg, anp4neg;
361   TLorentzVector p4comb, anp4comb;
362   
363   for (Int_t iTracks = 0; iTracks < esd->GetNumberOfTracks(); iTracks++) {
364     AliESDtrack* trackpos = esd->GetTrack(iTracks);
365     if (!trackpos) {
366       if (fDebug > 0) Printf("ERROR: Could not receive track %d", iTracks);
367       continue;
368     }
369     if (trackpos->GetSign() < 0) continue;
370     
371     AliExternalTrackParam *tpcTrackpos = (AliExternalTrackParam *)trackpos->GetTPCInnerParam();
372     if(!tpcTrackpos) continue;
373     ptrackpos[0]=tpcTrackpos->Px();
374     ptrackpos[1]=tpcTrackpos->Py();   
375     ptrackpos[2]=tpcTrackpos->Pz();  
376     
377     Bool_t firstLevelAcceptPosTrack=IsAcceptedForKink(esd, vertex, trackpos);
378     if(firstLevelAcceptPosTrack==kFALSE) continue;
379     
380     TVector3 posTrackMom(ptrackpos[0],ptrackpos[1],ptrackpos[2]);
381         
382     TParticle * partpos = stack->Particle(TMath::Abs(trackpos->GetLabel()));
383     if (!partpos) continue;
384     Int_t pdgpos = partpos->GetPdgCode();
385     Int_t mumlabelpos=partpos->GetFirstMother();
386     mumlabelpos = TMath::Abs(mumlabelpos);
387     TParticle * mumpos=stack->Particle(mumlabelpos);
388     if (!mumpos) continue;
389     Int_t mumpdgpos = mumpos->GetPdgCode();
390     
391     Int_t indexKinkPos=trackpos->GetKinkIndex(0);
392     Bool_t posKaonKinkFlag=0;
393     if(indexKinkPos<0) posKaonKinkFlag=IsKink(esd, indexKinkPos, posTrackMom);
394     
395     if(posKaonKinkFlag==1) anp4pos.SetVectM(posTrackMom,daughter1pdgMass);
396     
397     if(indexKinkPos==0) {
398
399     Bool_t secondLevelAcceptPosTrack=IsAcceptedForTrack(esd, vertex, trackpos);
400     if(secondLevelAcceptPosTrack==kFALSE) continue;
401
402       p4pos.SetVectM(posTrackMom, daughter2pdgMass);
403     
404     }
405         
406       for (Int_t j=0; j<esd->GetNumberOfTracks(); j++) {
407         if(iTracks==j) continue;
408         AliESDtrack* trackneg=esd->GetTrack(j);
409         if (trackneg->GetSign() > 0) continue;
410         
411         AliExternalTrackParam *tpcTrackneg = (AliExternalTrackParam *)trackneg->GetTPCInnerParam();
412         if(!tpcTrackneg) continue;
413         ptrackneg[0]=tpcTrackneg->Px();
414         ptrackneg[1]=tpcTrackneg->Py();   
415         ptrackneg[2]=tpcTrackneg->Pz();  
416     
417         Bool_t firstLevelAcceptNegTrack=IsAcceptedForKink(esd, vertex, trackneg);
418         if(firstLevelAcceptNegTrack==kFALSE) continue;  
419
420         TVector3 negTrackMom(ptrackneg[0],ptrackneg[1],ptrackneg[2]);
421         
422         TParticle * partneg = stack->Particle(TMath::Abs(trackneg->GetLabel()));
423         if (!partneg) continue;
424         Int_t pdgneg = partneg->GetPdgCode();
425         Int_t mumlabelneg=partneg->GetFirstMother();
426         mumlabelneg = TMath::Abs(mumlabelneg);
427         TParticle * mumneg=stack->Particle(mumlabelneg);
428         if (!mumneg) continue;
429         Int_t mumpdgneg = mumneg->GetPdgCode();
430         
431         Int_t indexKinkNeg=trackneg->GetKinkIndex(0);
432         Bool_t negKaonKinkFlag=0;
433         if(indexKinkNeg<0) negKaonKinkFlag=IsKink(esd, indexKinkNeg, negTrackMom);
434         
435         if(negKaonKinkFlag==1) p4neg.SetVectM(negTrackMom,daughter1pdgMass);
436         
437         if(indexKinkNeg==0)  {
438  
439            Bool_t secondLevelAcceptNegTrack=IsAcceptedForTrack(esd, vertex, trackneg);
440            if(secondLevelAcceptNegTrack==kFALSE) continue;  
441           
442           anp4neg.SetVectM(negTrackMom, daughter2pdgMass);
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((posKaonKinkFlag==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           fInvmassPt->Fill(p4comb.M(), p4comb.Vect().Pt());
459           if ((mumpdgpos==(antiresonancePDGcode))&&(mumpdgneg==(antiresonancePDGcode))&&(mumlabelpos==mumlabelneg)
460           &&(pdgpos==fdaughter2pdg)&&(pdgneg==(-fdaughter1pdg))&&(TMath::Abs(trackpos->GetLabel())>=0)&&(TMath::Abs(trackneg->GetLabel())>=0)&&(mumlabelpos>=0)&&(mumlabelneg>=0)) {
461             fOpeningAngle->Fill(openingAngle);
462             fInvMassTrue->Fill(p4comb.M());
463             fInvmassPtTrue->Fill(p4comb.M(), p4comb.Vect().Pt());
464             if((TMath::Abs(p4pos.Vect().Eta())<1.1)&&(TMath::Abs(p4neg.Vect().Eta())<1.1)&&(p4comb.Vect().Eta()<1.1)) {
465               fRecPt->Fill(p4comb.Vect().Pt());
466               fRecEta->Fill(p4comb.Vect().Eta());
467               fRecEtaPt->Fill(p4comb.Vect().Perp(),p4comb.Vect().Eta());
468
469             }
470
471            }
472           
473         }
474         
475         if(posKaonKinkFlag==1) {
476           anp4comb=anp4pos;
477           anp4comb+=anp4neg;  
478           fInvariantMass->Fill(anp4comb.M());
479           fInvmassPt->Fill(anp4comb.M(), anp4comb.Vect().Pt());
480           if ((mumpdgpos==resonancePDGcode)&&(mumpdgneg==resonancePDGcode)&&(mumlabelpos==mumlabelneg)
481           &&(pdgpos==fdaughter1pdg)&&(pdgneg==(-fdaughter2pdg))&&(TMath::Abs(trackpos->GetLabel())>=0)&&(TMath::Abs(trackneg->GetLabel())>=0)&&(mumlabelpos>=0)  &&(mumlabelneg>=0)) {
482             fOpeningAngle->Fill(openingAngle);
483             fInvMassTrue->Fill(p4comb.M());
484             fInvmassPtTrue->Fill(anp4comb.M(), anp4comb.Vect().Pt());
485             if((TMath::Abs(anp4neg.Vect().Eta())<1.1)&&(TMath::Abs(anp4pos.Vect().Eta())<1.1)&&(anp4comb.Vect().Eta()<1.1)) {   
486              fRecPt->Fill(anp4comb.Vect().Pt());
487              fRecEta->Fill(anp4comb.Vect().Eta());
488              fRecEtaPt->Fill(anp4comb.Vect().Pt(), anp4comb.Vect().Eta());
489            }
490
491          }
492
493         }
494          
495       } //inner track loop 
496
497   } //outer track loop 
498   
499   } // end fAnalysisType == ESD
500   
501 }      
502
503 //____________________________________________________________________//
504 Float_t AliResonanceKink::GetSigmaToVertex(AliESDtrack* esdTrack) const {
505   // Calculates the number of sigma to the vertex.
506   
507   Float_t b[2];
508   Float_t bRes[2];
509   Float_t bCov[3];
510
511     esdTrack->GetImpactParametersTPC(b,bCov);
512   
513   if (bCov[0]<=0 || bCov[2]<=0) {
514     bCov[0]=0; bCov[2]=0;
515   }
516   
517   bRes[0] = TMath::Sqrt(bCov[0]);
518   bRes[1] = TMath::Sqrt(bCov[2]);
519   
520   if (bRes[0] == 0 || bRes[1] ==0) return -1;
521   
522   Float_t d = TMath::Sqrt(TMath::Power(b[0]/bRes[0],2) + TMath::Power(b[1]/bRes[1],2));
523   
524   if (TMath::Exp(-d * d / 2) < 1e-10) return 1000;
525   
526   d = TMath::ErfInverse(1 - TMath::Exp(-d * d / 2)) * TMath::Sqrt(2);
527   
528   return d;
529 }
530
531 //________________________________________________________________________
532 const AliESDVertex* AliResonanceKink::GetEventVertex(const AliESDEvent* esd) const
533 {
534   // Get the vertex 
535   
536   const AliESDVertex* vertex = esd->GetPrimaryVertex();
537
538   if((vertex->GetStatus()==kTRUE)&&(vertex->GetNContributors()>2)) return vertex;
539   else
540   { 
541      vertex = esd->GetPrimaryVertexSPD();
542      if((vertex->GetStatus()==kTRUE)&&(vertex->GetNContributors()>2)) return vertex;
543      else
544      return 0;
545   }
546 }
547
548 //________________________________________________________________________
549
550  Bool_t AliResonanceKink::IsAcceptedForKink(AliESDEvent *localesd,
551             const AliESDVertex *localvertex, AliESDtrack* localtrack) {
552    // Apply the selections for each kink
553
554   Double_t gPt = 0.0, gPx = 0.0, gPy = 0.0, gPz = 0.0;
555   Double_t dca[2] = {0.0,0.0}, cov[3] = {0.0,0.0,0.0};  //The impact parameters and their covariance.
556   Double_t dca3D = 0.0;
557   
558   AliExternalTrackParam *tpcTrack = (AliExternalTrackParam *)localtrack->GetTPCInnerParam();
559   if(!tpcTrack) {
560     gPt = 0.0; gPx = 0.0; gPy = 0.0; gPz = 0.0;
561     dca[0] = -100.; dca[1] = -100.; dca3D = -100.;
562     cov[0] = -100.; cov[1] = -100.; cov[2] = -100.;
563   }
564   else {
565     gPt = tpcTrack->Pt();
566     gPx = tpcTrack->Px();
567     gPy = tpcTrack->Py();
568     gPz = tpcTrack->Pz();
569     tpcTrack->PropagateToDCA(localvertex,
570                localesd->GetMagneticField(),100.,dca,cov);
571   }
572   
573   if(GetSigmaToVertex(localtrack) > fMaxNSigmaToVertex) {
574       if (fDebug > 1) Printf("IsAcceptedKink: Track rejected because it has a %lf sigmas to vertex TPC (max. requested: %lf)",   GetSigmaToVertex(localtrack),fMaxNSigmaToVertex);
575       return kFALSE;
576   }
577   
578   if(TMath::Abs(dca[0]) > fMaxDCAxy) {
579       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);
580       return kFALSE;
581   }
582     
583   if(TMath::Abs(dca[1]) > fMaxDCAzaxis) {
584       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);
585       return kFALSE;
586   }
587   
588   if(gPt < fMinPtTrackCut) {
589       if (fDebug > 1) Printf("IsAcceptedKink: Track rejected because it has a min value of pt of %lf (min. requested: %lf)", gPt, fMinPtTrackCut);
590       return kFALSE;
591   } 
592   
593   return kTRUE;
594 }
595
596 //________________________________________________________________________
597 Bool_t AliResonanceKink::IsAcceptedForTrack(AliESDEvent *localesd,                                                                                                                                          const AliESDVertex *localvertex, AliESDtrack *localtrack) {
598    // Apply the selections for each track
599
600   Double_t gPt = 0.0, gPx = 0.0, gPy = 0.0, gPz = 0.0;
601   Double_t dca[2] = {0.0,0.0}, cov[3] = {0.0,0.0,0.0};  //The impact parameters and their covariance.
602   Double_t dca3D = 0.0;
603   
604   AliExternalTrackParam *tpcTrack = (AliExternalTrackParam *)localtrack->GetTPCInnerParam();
605   if(!tpcTrack) {
606     gPt = 0.0; gPx = 0.0; gPy = 0.0; gPz = 0.0;
607     dca[0] = -100.; dca[1] = -100.; dca3D = -100.;
608     cov[0] = -100.; cov[1] = -100.; cov[2] = -100.;
609   }
610   else {
611     gPt = tpcTrack->Pt();
612     gPx = tpcTrack->Px();
613     gPy = tpcTrack->Py();
614     gPz = tpcTrack->Pz();
615     tpcTrack->PropagateToDCA(localvertex,
616                localesd->GetMagneticField(),100.,dca,cov);
617   }
618   
619   Int_t fcls[200];
620   Int_t nClustersTPC=localtrack->GetTPCclusters(fcls);
621   Float_t chi2perTPCcluster=-1.0;
622   if(nClustersTPC!=0) chi2perTPCcluster=(localtrack->GetTPCchi2())/Float_t(nClustersTPC);
623   
624   Double_t extCov[15];
625   localtrack->GetExternalCovariance(extCov);
626   
627   if((localtrack->GetStatus() & AliESDtrack::kTPCrefit) == 0) {
628       if (fDebug > 1) Printf("IsAccepted: Track rejected because of no refited in TPC");
629       return kFALSE;
630   } 
631
632   if(nClustersTPC < fMinTPCclusters) {
633       if (fDebug > 1) Printf("IsAccepted: Track rejected because it has a value of nclusters (TPC) of %ld (min. requested: %ld)", nClustersTPC, fMinTPCclusters);
634       return kFALSE;
635   } 
636   
637   if(chi2perTPCcluster > fMaxChi2PerTPCcluster) {
638       if (fDebug > 1) Printf("IsAccepted: Track rejected because it has a value of chi2perTPCcluster of %lf (max. requested: %lf)", chi2perTPCcluster, fMaxChi2PerTPCcluster);
639       return kFALSE;
640   } 
641
642   if(extCov[0] > fMaxCov0) {
643       if (fDebug > 1) Printf("IsAccepted: Track rejected because it has a value of cov[0] of %lf (max. requested: %lf)", cov[0], fMaxCov0);
644       return kFALSE;
645   }
646   
647   if(extCov[2] > fMaxCov2) {
648       if (fDebug > 1) Printf("IsAccepted: Track rejected because it has a value of cov[2] of %lf (max. requested: %lf)", cov[2], fMaxCov2);
649       return kFALSE;
650   }
651     
652   if(extCov[5] > fMaxCov5) {
653       if (fDebug > 1) Printf("IsAccepted: Track rejected because it has a value of cov[5] of %lf (max. requested: %lf)", cov[5], fMaxCov5);
654       return kFALSE;
655   }  
656   
657   if(extCov[9] > fMaxCov9) {
658       if (fDebug > 1) Printf("IsAccepted: Track rejected because it has a value of cov[9] of %lf (max. requested: %lf)", cov[9], fMaxCov9);
659       return kFALSE;
660   }  
661   
662   if(extCov[14] > fMaxCov14) {
663       if (fDebug > 1) Printf("IsAccepted: Track rejected because it has a value of cov[14] of %lf (max. requested: %lf)", cov[14], fMaxCov14);
664       return kFALSE;
665   } 
666  
667   return kTRUE;
668
669 }
670
671 //________________________________________________________________________
672 Bool_t AliResonanceKink::IsKink(AliESDEvent *localesd, Int_t kinkIndex, TVector3 trackMom) 
673 {
674    // Test some kinematical criteria for each kink
675
676          AliESDkink *kink=localesd->GetKink(TMath::Abs(kinkIndex)-1);
677          const TVector3 motherMfromKink(kink->GetMotherP());
678          const TVector3 daughterMKink(kink->GetDaughterP());
679          Float_t qt=kink->GetQt();
680
681          Double_t maxDecAngKmu=f1->Eval(motherMfromKink.Mag(),0.,0.,0.);
682          Double_t maxDecAngpimu=f2->Eval(motherMfromKink.Mag(),0.,0.,0.);
683
684          Float_t kinkAngle=TMath::RadToDeg()*kink->GetAngle(2);
685          
686          Float_t energyDaughterMu=TMath::Sqrt(daughterMKink.Mag()*daughterMKink.Mag()+0.105658*0.105658);
687          Float_t p1XM= motherMfromKink.Px();
688          Float_t p1YM= motherMfromKink.Py();
689          Float_t p1ZM= motherMfromKink.Pz();
690          Float_t p2XM= daughterMKink.Px();
691          Float_t p2YM= daughterMKink.Py();
692          Float_t p2ZM= daughterMKink.Pz();
693          Float_t p3Daughter=TMath::Sqrt(((p1XM-p2XM)*(p1XM-p2XM))+((p1YM-p2YM)*(p1YM-p2YM))+((p1ZM-p2ZM)*(p1ZM-p2ZM)));
694          Double_t invariantMassKmu= TMath::Sqrt((energyDaughterMu+p3Daughter)*(energyDaughterMu+p3Daughter)-motherMfromKink.Mag()*motherMfromKink.Mag());
695
696          if((kinkAngle>maxDecAngpimu)&&(qt>0.05)&&(qt<0.25)&&((kink->GetR()>110.)&&(kink->GetR()<230.))&&(TMath::Abs(trackMom.Eta())<1.1)&&(invariantMassKmu<0.6)) {
697
698            if(trackMom.Mag()<=1.1) {
699                 return kTRUE;
700            }
701            else 
702            if (kinkAngle<maxDecAngKmu) {
703                 return kTRUE;
704            }
705          }
706          return kFALSE;
707 }