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