// Also, depending on the analysis mode (ESD or MC), fAnalysisType in the constructor must also be changed
//-----------------------------------------------------------------------------------------------------------------
-#include "TChain.h"
-#include "TTree.h"
#include "TH2D.h"
#include "TParticle.h"
-#include <TVector3.h>
-#include <TDatabasePDG.h>
-#include <TParticlePDG.h>
+#include "TDatabasePDG.h"
#include "TF1.h"
#include "TList.h"
-#include "TString.h"
-#include "AliMCEventHandler.h"
#include "AliMCEvent.h"
#include "AliResonanceKink.h"
#include "AliESDkink.h"
#include "AliStack.h"
#include "AliESDtrack.h"
#include "AliESDEvent.h"
+#include "AliExternalTrackParam.h"
+#include "AliMCParticle.h"
ClassImp(AliResonanceKink)
//________________________________________________________________________
AliResonanceKink::AliResonanceKink()
- : TObject(), 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),
- fhdr(0), fhdz(0), f1(0), f2(0), fAnalysisType(), fvtxz(0), fNbins(0), fLowX(0), fHighX(0), fdaughter1pdg(0), fdaughter2pdg(0), fresonancePDGcode(0)
-
+ : 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),
+ 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),
+fMinTPCclusters(0),fMaxChi2PerTPCcluster(0), fMaxCov0(0), fMaxCov2(0), fMaxCov5(0) , fMaxCov9(0), fMaxCov14(0) //, fTPCrefitFlag(kFALSE)
+, fInvmassPt(0), fInvmassPtTrue(0), fMCInvmassPt(0), fMCInvmassPtTrue(0)
{
// Constructor
}
//________________________________________________________________________
AliResonanceKink::AliResonanceKink(Int_t nbins, Float_t nlowx, Float_t nhighx, Int_t daughter1, Int_t daughter2, Int_t resonancePDG)
- : TObject(), 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),
- fhdr(0), fhdz(0), f1(0), f2(0), fAnalysisType(), fvtxz(0), fNbins(nbins), fLowX(nlowx), fHighX(nhighx), fdaughter1pdg(daughter1), fdaughter2pdg(daughter2), fresonancePDGcode(resonancePDG)
-
+ : 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),
+ 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),
+fMaxDCAxy(0), fMaxDCAzaxis(0), fMinTPCclusters(0), fMaxChi2PerTPCcluster(0), fMaxCov0(0), fMaxCov2(0), fMaxCov5(0), fMaxCov9(0), fMaxCov14(0) //, fTPCrefitFlag(kFALSE)
+, fInvmassPt(0), fInvmassPtTrue(0), fMCInvmassPt(0), fMCInvmassPtTrue(0)
{
// Constructor
fPhiBothKinks=new TH1D("fPhiBothKinks"," ",fNbins,fLowX,fHighX); // Applicable for phi(1020)
fRecPt=new TH1D("fRecPt"," ", 50,0.0,5.0);
- fRecEta=new TH1D("fRecEta"," ", 44,-1.1,1.1);
- fRecEtaPt=new TH2D("fRecEtaPt"," ", 50,0.0,5.0, 44,-1.1,1.1);
+ fRecEta=new TH1D("fRecEta"," ", 36,-0.9,0.9);
+ fRecEtaPt=new TH2D("fRecEtaPt"," ", 50,0.0,5.0, 36,-0.9,0.9);
fSimPt=new TH1D("fSimPt"," ", 50,0.0,5.0);
- fSimEta=new TH1D("fSimEta"," ", 44,-1.1,1.1);
- fSimEtaPt=new TH2D("fSimEtaPt"," ", 50,0.0,5.0, 44,-1.1,1.1);
+ fSimEta=new TH1D("fSimEta"," ", 36,-0.9,0.9);
+ fSimEtaPt=new TH2D("fSimEtaPt"," ", 50,0.0,5.0, 36,-0.9,0.9);
fSimPtKink=new TH1D("fSimPtKink"," ", 50,0.0,5.0);
- fSimEtaKink=new TH1D("fSimEtaKink"," ", 44,-1.1,1.1);
- fSimEtaPtKink=new TH2D("fSimEtaPtKink"," ", 50,0.0,5.0, 44,-1.1,1.1);
+ fSimEtaKink=new TH1D("fSimEtaKink"," ", 36,-0.9,0.9);
+ fSimEtaPtKink=new TH2D("fSimEtaPtKink"," ", 50,0.0,5.0, 36,-0.9,0.9);
fhdr=new TH1D("fhdr"," ", 100,0.0,5.0);
fhdz=new TH1D("fhdz"," ", 100,0.0,5.0);
f2->SetParameter(2,TMath::Pi());
fvtxz=new TH1D("fvtxz"," ", 100,-20.0,20.0);
-
+
+ fInvmassPt=new TH2D("fInvmassPt"," ",fNbins,fLowX,fHighX,100,0.0,10.0);
+ fInvmassPtTrue=new TH2D("fInvmassPtTrue"," ",fNbins,fLowX,fHighX,100,0.0,10.0);
+ fMCInvmassPt=new TH2D("fMCInvmassPt"," ",fNbins,fLowX,fHighX,100,0.0,10.0);
+ fMCInvmassPtTrue=new TH2D("fMCInvmassPtTrue"," ",fNbins,fLowX,fHighX,100,0.0,10.0);
}
//________________________________________________________________________
if(fSimEtaPtKink) delete fSimEtaPtKink;
if(fhdr) delete fhdr;
if(fhdz) delete fhdz;
- if(fvtxz) delete fvtxz;
+ if(fvtxz) delete fvtxz;
+ if(fInvmassPt) delete fInvmassPt;
+ if(fInvmassPtTrue) delete fInvmassPtTrue;
+ if(fMCInvmassPt) delete fMCInvmassPt;
+ if(fMCInvmassPtTrue) delete fMCInvmassPtTrue;
}
//________________________________________________________________________
TList* AliResonanceKink::GetHistogramList()
fListOfHistos->Add(fhdr);
fListOfHistos->Add(fhdz);
fListOfHistos->Add(fvtxz);
-
+ fListOfHistos->Add(fInvmassPt);
+ fListOfHistos->Add(fInvmassPtTrue);
+ fListOfHistos->Add(fMCInvmassPt);
+ fListOfHistos->Add(fMCInvmassPtTrue);
+
return fListOfHistos;
}
f2->SetParameter(2,TMath::Pi());
fvtxz=new TH1D("fvtxz"," ", 100,-20.0,20.0);
-}
+ fInvmassPt=new TH2D("fInvmassPt"," ",fNbins,fLowX,fHighX,100,0.0,10.0);
+ fInvmassPtTrue=new TH2D("fInvmassPtTrue"," ",fNbins,fLowX,fHighX,100,0.0,10.0);
+ fMCInvmassPt=new TH2D("fMCInvmassPt"," ",fNbins,fLowX,fHighX,100,0.0,10.0);
+ fMCInvmassPtTrue=new TH2D("fMCInvmassPtTrue"," ",fNbins,fLowX,fHighX,100,0.0,10.0);
+
+}
+
//________________________________________________________________________
void AliResonanceKink::Analyse(AliESDEvent* esd, AliMCEvent* mcEvent)
{
// Main loop
// Called for each event
Int_t resonancePDGcode, antiresonancePDGcode;
+ Double_t daughter1pdgMass, daughter2pdgMass;
- if (fdaughter1pdg==kdaughterKaon) {
+ if (fdaughter1pdg==kKPlus) {
resonancePDGcode=fresonancePDGcode;
antiresonancePDGcode=-fresonancePDGcode;
+ daughter1pdgMass=TDatabasePDG::Instance()->GetParticle(fdaughter1pdg)->Mass();
+ daughter2pdgMass=TDatabasePDG::Instance()->GetParticle(fdaughter2pdg)->Mass();
}
- if (fdaughter1pdg!=kdaughterKaon) {
+
+ if (fdaughter1pdg!=kKPlus) {
resonancePDGcode=-fresonancePDGcode;
antiresonancePDGcode=fresonancePDGcode;
- }
+ daughter1pdgMass=TDatabasePDG::Instance()->GetParticle(fdaughter2pdg)->Mass();
+ daughter2pdgMass=TDatabasePDG::Instance()->GetParticle(fdaughter1pdg)->Mass();
+ } // to ensure that daughter1pdgMass has always the kaon mass
+
if (fdaughter1pdg==fdaughter2pdg) {
resonancePDGcode=fresonancePDGcode;
antiresonancePDGcode=fresonancePDGcode;
}
- Double_t daughter1pdgMass=TDatabasePDG::Instance()->GetParticle(fdaughter1pdg)->Mass();
- Double_t daughter2pdgMass=TDatabasePDG::Instance()->GetParticle(fdaughter2pdg)->Mass();
-
if (!esd) {
Printf("ERROR: fESD not available");
return;
AliStack* stack=mcEvent->Stack();
if(fAnalysisType == "MC") {
+
+ const AliESDVertex* vertex = GetEventVertex(esd);
+ if(!vertex) return;
+ Double_t vtx[3];
+ vertex->GetXYZ(vtx);
+ fvtxz->Fill(vtx[2]);
+
for (Int_t iMc = 0; iMc < stack->GetNprimary(); ++iMc)
{
TParticle* particle = stack->Particle(iMc);
if (!particle)
{
- Printf("UNEXPECTED: particle with label %d not found in stack (mc loop)", iMc);
+ if (fDebug > 0) Printf("UNEXPECTED: particle with label %d not found in stack (mc loop)", iMc);
continue;
}
- if(TMath::Abs(particle->GetPdgCode())==fresonancePDGcode) {
+ if(TMath::Abs(particle->GetPdgCode())==fresonancePDGcode) {
+
+ if(particle->GetNDaughters()>2) continue;
+
Int_t firstD=particle->GetFirstDaughter();
Int_t lastD=particle->GetLastDaughter();
- TParticle *daughterParticle1=stack->Particle(firstD);
- TParticle *daughterParticle2=stack->Particle(lastD);
+
+ if((firstD<0)||(firstD>=stack->GetNtrack())) continue;
+ if((lastD<0)||(lastD>=stack->GetNtrack())) continue;
- TParticle* kaonFirstDaughter;
- Int_t mcProcessKaonFirstDaughter = -999;
+ TParticle *daughterParticle1 = 0;
+ TParticle *daughterParticle2 = 0;
+ AliMCParticle *mcDaughter1 = 0;
+ AliMCParticle *mcDaughter2 = 0;
+
+ if(fdaughter1pdg==kKPlus) {
+ daughterParticle1=stack->Particle(firstD);
+ daughterParticle2=stack->Particle(lastD);
+ mcDaughter1= (AliMCParticle*) mcEvent->GetTrack(firstD);
+ mcDaughter2= (AliMCParticle*) mcEvent->GetTrack(lastD);
+ }
+ else
+ if(fdaughter2pdg==kKPlus) {
+ daughterParticle1=stack->Particle(lastD);
+ daughterParticle2=stack->Particle(firstD);
+ mcDaughter1= (AliMCParticle*) mcEvent->GetTrack(lastD);
+ mcDaughter2= (AliMCParticle*) mcEvent->GetTrack(firstD);
+ } //to ensure that the first daughter is always the kaon
+
+ if(TMath::Abs(daughterParticle1->GetPdgCode())!=321) continue;
- for(Int_t ia=0; ia<daughterParticle1->GetNDaughters(); ia++){
- if ((daughterParticle1->GetFirstDaughter()+ia)!=-1) {
- kaonFirstDaughter=stack->Particle(daughterParticle1->GetFirstDaughter()+ia);
- mcProcessKaonFirstDaughter=kaonFirstDaughter->GetUniqueID();
+ TParticle * daughters1Daughter=0;
+ TParticle * daughters2Daughter=0;
+ Int_t mcProcessDaughters1Daughter = -999;
+ Int_t mcProcessDaughters2Daughter = -999;
+ AliMCParticle *mcDaughters1Daughter = 0;
+
+ if(mcDaughter1->Charge()==0) continue;
+ if(mcDaughter2->Charge()==0) continue; //accept resonance decays in two charged daughters
+
+ Int_t nDecayKaonDaughter=-99;
+ for(Int_t ia=0; ia<daughterParticle1->GetNDaughters(); ia++) {
+ if(((daughterParticle1->GetFirstDaughter()+ia)>0)&&((daughterParticle1->GetFirstDaughter()+ia)<stack->GetNtrack())) {
+ daughters1Daughter=stack->Particle(daughterParticle1->GetFirstDaughter()+ia);
+ mcProcessDaughters1Daughter=daughters1Daughter->GetUniqueID();
+ if(mcProcessDaughters1Daughter==4) {
+ nDecayKaonDaughter=daughterParticle1->GetFirstDaughter()+ia;
+ break;
+ }
}
- }
-
- if((daughterParticle1->Pt()>0.25)&&(daughterParticle2->Pt()>0.25)&&(TMath::Abs(daughterParticle1->Eta())<1.1)&& (TMath::Abs(daughterParticle2->Eta())<1.1)&&(TMath::Abs(particle->Eta())<1.1)) {
+ }
+
+ Int_t nProcessDaughter=-99;
+ for(Int_t ib=0; ib<daughterParticle2->GetNDaughters(); ib++) {
+ if(((daughterParticle2->GetFirstDaughter()+ib)>0)&&((daughterParticle2->GetFirstDaughter()+ib)<stack->GetNtrack())) {
+ daughters2Daughter=stack->Particle(daughterParticle2->GetFirstDaughter()+ib);
+ mcProcessDaughters2Daughter=daughters2Daughter->GetUniqueID();
+ if((mcProcessDaughters2Daughter==4)||(mcProcessDaughters2Daughter==13)) {
+ nProcessDaughter=mcProcessDaughters2Daughter;
+ break;
+ }
+ }
+ }
+
+ Int_t numberOfCharged=0;
+ if((mcProcessDaughters1Daughter==4)&&(nDecayKaonDaughter>=0)) {
+ for(Int_t ic=nDecayKaonDaughter; ic<=daughterParticle1->GetLastDaughter(); ic++) {
+ if ((ic>=0)&&(ic<stack->GetNtrack())) mcDaughters1Daughter= dynamic_cast<AliMCParticle*>(mcEvent->GetTrack(ic));
+ else continue;
+ if(mcDaughters1Daughter->Charge()!=0) numberOfCharged=numberOfCharged+1;
+ }
+ }
+
+ if(numberOfCharged>=2) continue; // leave out kaon decays to more than one charged daughter
+
+ if ((particle->Pt()>0.25)&&(TMath::Abs(particle->Eta())<0.9)) {
+
fSimEta->Fill(particle->Eta());
- fSimPt->Fill(particle->Pt());
fSimEtaPt->Fill(particle->Pt(), particle->Eta());
- if(mcProcessKaonFirstDaughter==4) {
- fSimPtKink->Fill(particle->Pt());
- fSimEtaKink->Fill(particle->Eta());
- fSimEtaPtKink->Fill(particle->Pt(), particle->Eta());
- }
- }
- }
- }
+ fSimPt->Fill(particle->Pt());
+ fMCInvmassPtTrue->Fill(particle->GetMass(), particle->Pt());
+
+ if((daughterParticle1->Pt()>0.25)&&(TMath::Abs(daughterParticle1->Eta())<0.9)&&(daughterParticle2->Pt()>0.25)&&(TMath::Abs(daughterParticle2->Eta())<0.9)) {
+ if((mcProcessDaughters1Daughter==4)&&(daughters1Daughter->R()>120.)&&(daughters1Daughter->R()<220.)&&( (nProcessDaughter<0)||((daughters2Daughter->R()>120.)&&(nProcessDaughter>0)))) { //below are the findable
+ fSimPtKink->Fill(particle->Pt());
+ fSimEtaKink->Fill(particle->Eta());
+ fSimEtaPtKink->Fill(particle->Pt(), particle->Eta());
+ fMCInvmassPt->Fill(particle->GetMass(), particle->Pt());
+ }
+ }
+
+ } // for the generated spectra
+
+ } //for the particular resonance
+ } //MC loop
} // end fAnalysisType==MC
else
for (Int_t iTracks = 0; iTracks < esd->GetNumberOfTracks(); iTracks++) {
AliESDtrack* trackpos = esd->GetTrack(iTracks);
if (!trackpos) {
- Printf("ERROR: Could not receive track %d", iTracks);
+ if (fDebug > 0) Printf("ERROR: Could not receive track %d", iTracks);
continue;
}
if (trackpos->GetSign() < 0) continue;
- trackpos->GetPxPyPz(ptrackpos);
+ AliExternalTrackParam *tpcTrackpos = (AliExternalTrackParam *)trackpos->GetTPCInnerParam();
+ if(!tpcTrackpos) continue;
+ ptrackpos[0]=tpcTrackpos->Px();
+ ptrackpos[1]=tpcTrackpos->Py();
+ ptrackpos[2]=tpcTrackpos->Pz();
- Float_t nSigmaToVertex = GetSigmaToVertex(trackpos);
-
- Float_t bpos[2];
- Float_t bCovpos[3];
- trackpos->GetImpactParameters(bpos,bCovpos);
-
- if (bCovpos[0]<=0 || bCovpos[2]<=0) {
- Printf("Estimated b resolution lower or equal zero!");
- bCovpos[0]=0; bCovpos[2]=0;
- }
-
- Float_t dcaToVertexXYpos = bpos[0];
- Float_t dcaToVertexZpos = bpos[1];
-
- fhdr->Fill(dcaToVertexXYpos);
- fhdz->Fill(dcaToVertexZpos);
-
- if(nSigmaToVertex>=4) continue;
- if((dcaToVertexXYpos>3.0)||(dcaToVertexZpos>3.0)) continue;
+ Bool_t firstLevelAcceptPosTrack=IsAcceptedForKink(esd, vertex, trackpos);
+ if(firstLevelAcceptPosTrack==kFALSE) continue;
TVector3 posTrackMom(ptrackpos[0],ptrackpos[1],ptrackpos[2]);
-
- if(posTrackMom.Perp()<=0.25) continue;
-
+
TParticle * partpos = stack->Particle(TMath::Abs(trackpos->GetLabel()));
if (!partpos) continue;
Int_t pdgpos = partpos->GetPdgCode();
Int_t mumpdgpos = mumpos->GetPdgCode();
Int_t indexKinkPos=trackpos->GetKinkIndex(0);
- Int_t kaonKinkFlag=0;
- if(indexKinkPos<0){
-
- AliESDkink *poskink=esd->GetKink(TMath::Abs(indexKinkPos)-1);
- const TVector3 motherMfromKinkPos(poskink->GetMotherP());
- const TVector3 daughterMKinkPos(poskink->GetDaughterP());
- Float_t posQt=poskink->GetQt();
-
- Double_t maxDecAngKmuPos=f1->Eval(motherMfromKinkPos.Mag(),0.,0.,0.);
- Double_t maxDecAngpimuPos=f2->Eval(motherMfromKinkPos.Mag(),0.,0.,0.);
-
- Float_t kinkAnglePos=TMath::RadToDeg()*poskink->GetAngle(2);
-
- Float_t energyDaughterMu=TMath::Sqrt(daughterMKinkPos.Mag()*daughterMKinkPos.Mag()+0.105658*0.105658);
- Float_t p1XM= motherMfromKinkPos.Px();
- Float_t p1YM= motherMfromKinkPos.Py();
- Float_t p1ZM= motherMfromKinkPos.Pz();
- Float_t p2XM= daughterMKinkPos.Px();
- Float_t p2YM= daughterMKinkPos.Py();
- Float_t p2ZM= daughterMKinkPos.Pz();
- Float_t p3Daughter=TMath::Sqrt(((p1XM-p2XM)*(p1XM-p2XM))+((p1YM-p2YM)*(p1YM-p2YM))+((p1ZM-p2ZM)*(p1ZM-p2ZM)));
- Double_t invariantMassKmuPos= TMath::Sqrt((energyDaughterMu+p3Daughter)*(energyDaughterMu+p3Daughter)-motherMfromKinkPos.Mag()*motherMfromKinkPos.Mag());
-
- if((kinkAnglePos>maxDecAngpimuPos)&&(posQt>0.05)&&(posQt<0.25)&&((poskink->GetR()>110.)&&(poskink->GetR()<230.))&&(TMath::Abs(posTrackMom.Eta())<1.1)&&(invariantMassKmuPos<0.6)) {
-
- if(posTrackMom.Mag()<=1.1) {
- kaonKinkFlag=1;
- }
- else
- if (kinkAnglePos<maxDecAngKmuPos) {
- kaonKinkFlag=1;
- }
- }
-
- } //End Kink Information
- if(kaonKinkFlag==1) anp4pos.SetVectM(posTrackMom,daughter1pdgMass);
+ if(indexKinkPos>0) continue;
+
+ Bool_t posKaonKinkFlag=0;
+
+ if(indexKinkPos<0) {
+ posKaonKinkFlag=IsKink(esd, indexKinkPos, posTrackMom);
+
+ if(posKaonKinkFlag==1) anp4pos.SetVectM(posTrackMom,daughter1pdgMass);
+ if(posKaonKinkFlag==0) continue;
+ }
if(indexKinkPos==0) {
- UInt_t status=trackpos->GetStatus();
- if((status&AliESDtrack::kTPCrefit)==0) continue;
- if(trackpos->GetTPCclusters(0)<50) continue;
- if((trackpos->GetTPCchi2()/trackpos->GetTPCclusters(0))>3.5) continue;
- Double_t extCovPos[15];
- trackpos->GetExternalCovariance(extCovPos);
- if(extCovPos[0]>2) continue;
- if(extCovPos[2]>2) continue;
- if(extCovPos[5]>0.5) continue;
- if(extCovPos[9]>0.5) continue;
- if(extCovPos[14]>2) continue;
-
+
+ Bool_t secondLevelAcceptPosTrack=IsAcceptedForTrack(esd, vertex, trackpos);
+ if(secondLevelAcceptPosTrack==kFALSE) continue;
+
p4pos.SetVectM(posTrackMom, daughter2pdgMass);
}
AliESDtrack* trackneg=esd->GetTrack(j);
if (trackneg->GetSign() > 0) continue;
- trackneg->GetPxPyPz(ptrackneg);
- Float_t negSigmaToVertex = GetSigmaToVertex(trackneg);
-
- Float_t bneg[2];
- Float_t bCovneg[3];
- trackneg->GetImpactParameters(bneg,bCovneg);
- if (bCovneg[0]<=0 || bCovneg[2]<=0) {
- Printf("Estimated b resolution lower or equal zero!");
- bCovneg[0]=0; bCovneg[2]=0;
- }
-
- Float_t dcaToVertexXYneg = bneg[0];
- Float_t dcaToVertexZneg = bneg[1];
+ AliExternalTrackParam *tpcTrackneg = (AliExternalTrackParam *)trackneg->GetTPCInnerParam();
+ if(!tpcTrackneg) continue;
+ ptrackneg[0]=tpcTrackneg->Px();
+ ptrackneg[1]=tpcTrackneg->Py();
+ ptrackneg[2]=tpcTrackneg->Pz();
- fhdr->Fill(dcaToVertexXYneg);
- fhdz->Fill(dcaToVertexZneg);
-
- if(negSigmaToVertex>=4) continue;
- if((dcaToVertexXYneg>3.0)||(dcaToVertexZneg>3.0)) continue;
+ Bool_t firstLevelAcceptNegTrack=IsAcceptedForKink(esd, vertex, trackneg);
+ if(firstLevelAcceptNegTrack==kFALSE) continue;
TVector3 negTrackMom(ptrackneg[0],ptrackneg[1],ptrackneg[2]);
-
- if(negTrackMom.Perp()<=0.25) continue;
TParticle * partneg = stack->Particle(TMath::Abs(trackneg->GetLabel()));
if (!partneg) continue;
Int_t mumpdgneg = mumneg->GetPdgCode();
Int_t indexKinkNeg=trackneg->GetKinkIndex(0);
- Int_t negKaonKinkFlag=0;
- if(indexKinkNeg<0){
-
- AliESDkink *negkink=esd->GetKink(TMath::Abs(indexKinkNeg)-1);
- const TVector3 motherMfromKinkNeg(negkink->GetMotherP());
- const TVector3 daughterMKinkNeg(negkink->GetDaughterP());
- Float_t negQt=negkink->GetQt();
-
- Double_t maxDecAngKmuNeg=f1->Eval(motherMfromKinkNeg.Mag(),0.,0.,0.);
- Double_t maxDecAngpimuNeg=f2->Eval(motherMfromKinkNeg.Mag(),0.,0.,0.);
-
- Float_t kinkAngleNeg=TMath::RadToDeg()*negkink->GetAngle(2);
-
- Float_t energyDaughterMuNeg=TMath::Sqrt(daughterMKinkNeg.Mag()*daughterMKinkNeg.Mag()+0.105658*0.105658);
- Float_t p1XMNeg= motherMfromKinkNeg.Px();
- Float_t p1YMNeg= motherMfromKinkNeg.Py();
- Float_t p1ZMNeg= motherMfromKinkNeg.Pz();
- Float_t p2XMNeg= daughterMKinkNeg.Px();
- Float_t p2YMNeg= daughterMKinkNeg.Py();
- Float_t p2ZMNeg= daughterMKinkNeg.Pz();
- Float_t p3DaughterNeg=TMath::Sqrt(((p1XMNeg-p2XMNeg)*(p1XMNeg-p2XMNeg))+((p1YMNeg-p2YMNeg)*(p1YMNeg-p2YMNeg))+((p1ZMNeg-p2ZMNeg)*(p1ZMNeg-p2ZMNeg)));
- Double_t invariantMassKmuNeg= TMath::Sqrt((energyDaughterMuNeg+p3DaughterNeg)*(energyDaughterMuNeg+p3DaughterNeg)-motherMfromKinkNeg.Mag()*motherMfromKinkNeg.Mag());
-
- if((kinkAngleNeg>maxDecAngpimuNeg)&&(negQt>0.05)&&(negQt<0.25)&&((negkink->GetR()>110.)&&(negkink->GetR()<230.))&&(TMath::Abs(negTrackMom.Eta())<1.1)&&(invariantMassKmuNeg<0.6)) {
-
- if(negTrackMom.Mag()<=1.1) {
- negKaonKinkFlag=1;
- }
- else
- if (kinkAngleNeg<maxDecAngKmuNeg) {
- negKaonKinkFlag=1;
- }
- }
-
- } //End Kink Information
- if(negKaonKinkFlag==1) p4neg.SetVectM(negTrackMom,daughter1pdgMass);
+ if(indexKinkNeg>0) continue;
+
+ Bool_t negKaonKinkFlag=0;
+
+ if(indexKinkNeg<0) {
+ negKaonKinkFlag=IsKink(esd, indexKinkNeg, negTrackMom);
+
+ if(negKaonKinkFlag==1) p4neg.SetVectM(negTrackMom,daughter1pdgMass);
+ if(negKaonKinkFlag==0) continue;
+ }
if(indexKinkNeg==0) {
- UInt_t statusneg=trackneg->GetStatus();
-
- if((statusneg&AliESDtrack::kTPCrefit)==0) continue;
-
- if(trackneg->GetTPCclusters(0)<50) continue;
- if((trackneg->GetTPCchi2()/trackneg->GetTPCclusters(0))>3.5) continue;
- Double_t extCovneg[15];
- trackneg->GetExternalCovariance(extCovneg);
- if(extCovneg[0]>2) continue;
- if(extCovneg[2]>2) continue;
- if(extCovneg[5]>0.5) continue;
- if(extCovneg[9]>0.5) continue;
- if(extCovneg[14]>2) continue;
-
+
+ Bool_t secondLevelAcceptNegTrack=IsAcceptedForTrack(esd, vertex, trackneg);
+ if(secondLevelAcceptNegTrack==kFALSE) continue;
+
anp4neg.SetVectM(negTrackMom, daughter2pdgMass);
}
Double_t openingAngle=(ptrackpos[0]*ptrackneg[0]+ptrackpos[1]*ptrackneg[1]+ptrackpos[2]*ptrackneg[2])/(posTrackMom.Mag()*negTrackMom.Mag());
- if((kaonKinkFlag==1)&&(negKaonKinkFlag==1)) {
- p4comb=anp4pos;
- p4comb+=p4neg;
- if(openingAngle>0.6) fPhiBothKinks->Fill(p4comb.M());
+ if((posKaonKinkFlag==1)&&(negKaonKinkFlag==1)) {
+ p4comb=anp4pos;
+ p4comb+=p4neg;
+ 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)) {
+ if(openingAngle>0.6) fPhiBothKinks->Fill(p4comb.M());
+ }
}
if(negKaonKinkFlag==1) {
p4comb=p4pos;
p4comb+=p4neg;
- fInvariantMass->Fill(p4comb.M());
- if ((mumpdgpos==(antiresonancePDGcode))&&(mumpdgneg==(antiresonancePDGcode))&&(mumlabelpos==mumlabelneg)
- &&(pdgpos==fdaughter2pdg)&&(pdgneg==(-fdaughter1pdg))&&(TMath::Abs(trackpos->GetLabel())>=0)&&(TMath::Abs(trackneg->GetLabel())>=0)&&(mumlabelpos>=0)&&(mumlabelneg>=0)) {
- fOpeningAngle->Fill(openingAngle);
- fInvMassTrue->Fill(p4comb.M());
- if((TMath::Abs(p4pos.Vect().Eta())<1.1)&&(TMath::Abs(p4neg.Vect().Eta())<1.1)&&(p4comb.Vect().Eta()<1.1)) {
+
+ if(p4comb.Vect().Pt()<=0.25) continue;
+
+ if((TMath::Abs(p4pos.Vect().Eta())<0.9)&&(TMath::Abs(p4neg.Vect().Eta())<0.9)&&(p4comb.Vect().Eta()<0.9)) {
+ fInvariantMass->Fill(p4comb.M());
+ fInvmassPt->Fill(p4comb.M(), p4comb.Vect().Pt());
+ if ((mumpdgpos==(antiresonancePDGcode))&&(mumpdgneg==(antiresonancePDGcode))&&(mumlabelpos==mumlabelneg)
+ &&(pdgpos==fdaughter2pdg)&&(pdgneg==(-fdaughter1pdg))&&(TMath::Abs(trackpos->GetLabel())>=0)&&(TMath::Abs(trackneg->GetLabel())>=0)&&(mumlabelpos>=0)&&(mumlabelneg>=0)) {
+ fOpeningAngle->Fill(openingAngle);
+ fInvMassTrue->Fill(p4comb.M());
+ fInvmassPtTrue->Fill(p4comb.M(), p4comb.Vect().Pt());
+
fRecPt->Fill(p4comb.Vect().Pt());
fRecEta->Fill(p4comb.Vect().Eta());
fRecEtaPt->Fill(p4comb.Vect().Perp(),p4comb.Vect().Eta());
}
- if(kaonKinkFlag==1) {
+ if(posKaonKinkFlag==1) {
anp4comb=anp4pos;
anp4comb+=anp4neg;
- fInvariantMass->Fill(anp4comb.M());
- if ((mumpdgpos==resonancePDGcode)&&(mumpdgneg==resonancePDGcode)&&(mumlabelpos==mumlabelneg)
- &&(pdgpos==fdaughter1pdg)&&(pdgneg==(-fdaughter2pdg))&&(TMath::Abs(trackpos->GetLabel())>=0)&&(TMath::Abs(trackneg->GetLabel())>=0)&&(mumlabelpos>=0) &&(mumlabelneg>=0)) {
- fOpeningAngle->Fill(openingAngle);
- fInvMassTrue->Fill(p4comb.M());
- if((TMath::Abs(anp4neg.Vect().Eta())<1.1)&&(TMath::Abs(anp4pos.Vect().Eta())<1.1)&&(anp4comb.Vect().Eta()<1.1)) {
- fRecPt->Fill(anp4comb.Vect().Pt());
- fRecEta->Fill(anp4comb.Vect().Eta());
- fRecEtaPt->Fill(anp4comb.Vect().Pt(), anp4comb.Vect().Eta());
+
+ if(anp4comb.Vect().Pt()<=0.25) continue;
+
+ if((TMath::Abs(anp4neg.Vect().Eta())<0.9)&&(TMath::Abs(anp4pos.Vect().Eta())<0.9)&&(anp4comb.Vect().Eta()<0.9)) {
+ fInvariantMass->Fill(anp4comb.M());
+ fInvmassPt->Fill(anp4comb.M(), anp4comb.Vect().Pt());
+ if ((mumpdgpos==resonancePDGcode)&&(mumpdgneg==resonancePDGcode)&&(mumlabelpos==mumlabelneg)
+ &&(pdgpos==fdaughter1pdg)&&(pdgneg==(-fdaughter2pdg))&&(TMath::Abs(trackpos->GetLabel())>=0)&&(TMath::Abs(trackneg->GetLabel())>=0)&&(mumlabelpos>=0) &&(mumlabelneg>=0)) {
+ fOpeningAngle->Fill(openingAngle);
+ fInvMassTrue->Fill(anp4comb.M());
+ fInvmassPtTrue->Fill(anp4comb.M(), anp4comb.Vect().Pt());
+
+ fRecPt->Fill(anp4comb.Vect().Pt());
+ fRecEta->Fill(anp4comb.Vect().Eta());
+ fRecEtaPt->Fill(anp4comb.Vect().Pt(), anp4comb.Vect().Eta());
}
}
Float_t bRes[2];
Float_t bCov[3];
- esdTrack->GetImpactParameters(b,bCov);
+ esdTrack->GetImpactParametersTPC(b,bCov);
if (bCov[0]<=0 || bCov[2]<=0) {
bCov[0]=0; bCov[2]=0;
}
}
+//________________________________________________________________________
+
+ Bool_t AliResonanceKink::IsAcceptedForKink(AliESDEvent *localesd,
+ const AliESDVertex *localvertex, AliESDtrack* localtrack) {
+ // Apply the selections for each kink
+
+ Double_t gPt = 0.0, gPx = 0.0, gPy = 0.0, gPz = 0.0;
+ Double_t dca[2] = {0.0,0.0}, cov[3] = {0.0,0.0,0.0}; //The impact parameters and their covariance.
+ Double_t dca3D = 0.0;
+
+ AliExternalTrackParam *tpcTrack = (AliExternalTrackParam *)localtrack->GetTPCInnerParam();
+ if(!tpcTrack) {
+ gPt = 0.0; gPx = 0.0; gPy = 0.0; gPz = 0.0;
+ dca[0] = -100.; dca[1] = -100.; dca3D = -100.;
+ cov[0] = -100.; cov[1] = -100.; cov[2] = -100.;
+ }
+ else {
+ gPt = tpcTrack->Pt();
+ gPx = tpcTrack->Px();
+ gPy = tpcTrack->Py();
+ gPz = tpcTrack->Pz();
+ tpcTrack->PropagateToDCA(localvertex,
+ localesd->GetMagneticField(),100.,dca,cov);
+ }
+
+ if(GetSigmaToVertex(localtrack) > fMaxNSigmaToVertex) {
+ if (fDebug > 1) Printf("IsAcceptedKink: Track rejected because it has a %lf sigmas to vertex TPC (max. requested: %lf)", GetSigmaToVertex(localtrack),fMaxNSigmaToVertex);
+ return kFALSE;
+ }
+
+ if(TMath::Abs(dca[0]) > fMaxDCAxy) {
+ 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);
+ return kFALSE;
+ }
+
+ if(TMath::Abs(dca[1]) > fMaxDCAzaxis) {
+ 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);
+ return kFALSE;
+ }
+
+ if(gPt <= fMinPtTrackCut) {
+ if (fDebug > 1) Printf("IsAcceptedKink: Track rejected because it has a min value of pt of %lf (min. requested: %lf)", gPt, fMinPtTrackCut);
+ return kFALSE;
+ }
+
+ return kTRUE;
+}
+
+//________________________________________________________________________
+Bool_t AliResonanceKink::IsAcceptedForTrack(AliESDEvent *localesd, const AliESDVertex *localvertex, AliESDtrack *localtrack) {
+ // Apply the selections for each track
+
+ Double_t gPt = 0.0, gPx = 0.0, gPy = 0.0, gPz = 0.0;
+ Double_t dca[2] = {0.0,0.0}, cov[3] = {0.0,0.0,0.0}; //The impact parameters and their covariance.
+ Double_t dca3D = 0.0;
+
+ AliExternalTrackParam *tpcTrack = (AliExternalTrackParam *)localtrack->GetTPCInnerParam();
+ if(!tpcTrack) {
+ gPt = 0.0; gPx = 0.0; gPy = 0.0; gPz = 0.0;
+ dca[0] = -100.; dca[1] = -100.; dca3D = -100.;
+ cov[0] = -100.; cov[1] = -100.; cov[2] = -100.;
+ }
+ else {
+ gPt = tpcTrack->Pt();
+ gPx = tpcTrack->Px();
+ gPy = tpcTrack->Py();
+ gPz = tpcTrack->Pz();
+ tpcTrack->PropagateToDCA(localvertex,
+ localesd->GetMagneticField(),100.,dca,cov);
+ }
+
+ Int_t fcls[200];
+ Int_t nClustersTPC=localtrack->GetTPCclusters(fcls);
+ Float_t chi2perTPCcluster=-1.0;
+ if(nClustersTPC!=0) chi2perTPCcluster=(localtrack->GetTPCchi2())/Float_t(nClustersTPC);
+
+ Double_t extCov[15];
+ localtrack->GetExternalCovariance(extCov);
+
+ if((localtrack->GetStatus() & AliESDtrack::kTPCrefit) == 0) {
+ if (fDebug > 1) Printf("IsAccepted: Track rejected because of no refited in TPC");
+ return kFALSE;
+ }
+ if(nClustersTPC < fMinTPCclusters) {
+ if (fDebug > 1) Printf("IsAccepted: Track rejected because it has a value of nclusters (TPC) of %ld (min. requested: %ld)", nClustersTPC, fMinTPCclusters);
+ return kFALSE;
+ }
+
+ if(chi2perTPCcluster > fMaxChi2PerTPCcluster) {
+ if (fDebug > 1) Printf("IsAccepted: Track rejected because it has a value of chi2perTPCcluster of %lf (max. requested: %lf)", chi2perTPCcluster, fMaxChi2PerTPCcluster);
+ return kFALSE;
+ }
+
+ if(extCov[0] > fMaxCov0) {
+ if (fDebug > 1) Printf("IsAccepted: Track rejected because it has a value of cov[0] of %lf (max. requested: %lf)", cov[0], fMaxCov0);
+ return kFALSE;
+ }
+
+ if(extCov[2] > fMaxCov2) {
+ if (fDebug > 1) Printf("IsAccepted: Track rejected because it has a value of cov[2] of %lf (max. requested: %lf)", cov[2], fMaxCov2);
+ return kFALSE;
+ }
+
+ if(extCov[5] > fMaxCov5) {
+ if (fDebug > 1) Printf("IsAccepted: Track rejected because it has a value of cov[5] of %lf (max. requested: %lf)", cov[5], fMaxCov5);
+ return kFALSE;
+ }
+
+ if(extCov[9] > fMaxCov9) {
+ if (fDebug > 1) Printf("IsAccepted: Track rejected because it has a value of cov[9] of %lf (max. requested: %lf)", cov[9], fMaxCov9);
+ return kFALSE;
+ }
+
+ if(extCov[14] > fMaxCov14) {
+ if (fDebug > 1) Printf("IsAccepted: Track rejected because it has a value of cov[14] of %lf (max. requested: %lf)", cov[14], fMaxCov14);
+ return kFALSE;
+ }
+ return kTRUE;
+
+}
+
+//________________________________________________________________________
+Bool_t AliResonanceKink::IsKink(AliESDEvent *localesd, Int_t kinkIndex, TVector3 trackMom)
+{
+ // Test some kinematical criteria for each kink
+
+ AliESDkink *kink=localesd->GetKink(TMath::Abs(kinkIndex)-1);
+ const TVector3 motherMfromKink(kink->GetMotherP());
+ const TVector3 daughterMKink(kink->GetDaughterP());
+ Float_t qt=kink->GetQt();
+
+ Double_t maxDecAngKmu=f1->Eval(motherMfromKink.Mag(),0.,0.,0.);
+ Double_t maxDecAngpimu=f2->Eval(motherMfromKink.Mag(),0.,0.,0.);
+
+ Float_t kinkAngle=TMath::RadToDeg()*kink->GetAngle(2);
+
+ Float_t energyDaughterMu=TMath::Sqrt(daughterMKink.Mag()*daughterMKink.Mag()+0.105658*0.105658);
+ Float_t p1XM= motherMfromKink.Px();
+ Float_t p1YM= motherMfromKink.Py();
+ Float_t p1ZM= motherMfromKink.Pz();
+ Float_t p2XM= daughterMKink.Px();
+ Float_t p2YM= daughterMKink.Py();
+ Float_t p2ZM= daughterMKink.Pz();
+ Float_t p3Daughter=TMath::Sqrt(((p1XM-p2XM)*(p1XM-p2XM))+((p1YM-p2YM)*(p1YM-p2YM))+((p1ZM-p2ZM)*(p1ZM-p2ZM)));
+ Double_t invariantMassKmu= TMath::Sqrt((energyDaughterMu+p3Daughter)*(energyDaughterMu+p3Daughter)-motherMfromKink.Mag()*motherMfromKink.Mag());
+
+ if((kinkAngle>maxDecAngpimu)&&(qt>0.05)&&(qt<0.25)&&((kink->GetR()>120.)&&(kink->GetR()<220.))&&(TMath::Abs(trackMom.Eta())<0.9)&&(invariantMassKmu<0.6)) {
+
+ if(trackMom.Mag()<=1.1) {
+ return kTRUE;
+ }
+ else
+ if (kinkAngle<maxDecAngKmu) {
+ return kTRUE;
+ }
+ }
+ return kFALSE;
+}
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff