From 92adf4f6b7345f0f4dc83298341687a9eab1516e Mon Sep 17 00:00:00 2001 From: bhippoly Date: Tue, 28 Jul 2009 15:53:16 +0000 Subject: [PATCH] Update for Kink tasks from Evi Ganoti --- PWG2/KINK/AliAnalysisKinkESDMC.cxx | 245 ++++++++------- PWG2/KINK/AliAnalysisKinkESDMC.h | 84 ++--- PWG2/KINK/AliAnalysisTaskKinkResonance.cxx | 79 ++--- PWG2/KINK/AliAnalysisTaskKinkResonance.h | 16 +- PWG2/KINK/AliResonanceKink.cxx | 348 ++++++++++++--------- PWG2/KINK/AliResonanceKink.h | 118 +++++-- PWG2/KINK/macros/AddTaskKinkResonance.C | 11 + 7 files changed, 507 insertions(+), 394 deletions(-) diff --git a/PWG2/KINK/AliAnalysisKinkESDMC.cxx b/PWG2/KINK/AliAnalysisKinkESDMC.cxx index 6207ea00767..a75f752caaf 100644 --- a/PWG2/KINK/AliAnalysisKinkESDMC.cxx +++ b/PWG2/KINK/AliAnalysisKinkESDMC.cxx @@ -29,51 +29,29 @@ #include #include "TF1.h" -#include "AliAnalysisTaskSE.h" +#include "AliAnalysisTask.h" #include "AliAnalysisManager.h" +#include "AliVEvent.h" #include "AliESDEvent.h" -#include "AliESDInputHandler.h" - -#include "AliMCEventHandler.h" #include "AliMCEvent.h" - #include "AliAnalysisKinkESDMC.h" - #include "AliStack.h" #include "AliESDpid.h" #include "AliESDkink.h" ClassImp(AliAnalysisKinkESDMC) -//________________________________________________________________________ -AliAnalysisKinkESDMC::AliAnalysisKinkESDMC() - : AliAnalysisTaskSE(), fESD(0),fHistPtESD(),fHistPt(0),fHistQtAll(0),fHistQt1(0),fHistQt2(0) - , fHistPtKaon(0),fHistPtKPDG(0),fHistEta(0),fHistEtaK(0),fptKMC(0),fMultiplMC(0),fESDMult(0),fgenpt(0),frad(0), - fKinkKaon(0), fKinkKaonBg(0), fM1kaon(0), fgenPtEtR(0),fPtKink(0), fptKink(0), - fcodeH(0), fdcodeH(0), fAngMomK(0),fAngMomPi(0), - fRpr(0),fZpr(0), - fAngMomKC(0), - fZvXv(0),fZvYv(0), fXvYv(0), fPtPrKink(0), f1(0), f2(0), - fListOfHistos(0) - -{ - // Constructor - -} - //________________________________________________________________________ AliAnalysisKinkESDMC::AliAnalysisKinkESDMC(const char *name) - : AliAnalysisTaskSE(name), fESD(0),fHistPtESD(0),fHistPt(0),fHistQtAll(0),fHistQt1(0),fHistQt2(0) + : AliAnalysisTaskSE(name), fHistPtESD(0),fHistPt(0),fHistQtAll(0),fHistQt1(0),fHistQt2(0) , fHistPtKaon(0),fHistPtKPDG(0),fHistEta(0),fHistEtaK(0),fptKMC(0),fMultiplMC(0),fESDMult(0),fgenpt(0),frad(0), fKinkKaon(0), fKinkKaonBg(0), fM1kaon(0), fgenPtEtR(0),fPtKink(0), fptKink(0), - fcodeH(0), fdcodeH(0), fAngMomK(0),fAngMomPi(0), + fcodeH(0), fdcodeH(0), fAngMomK(0),fAngMomPi(0), fAngMomKC(0), fMultESDK(0), fMultMCK(0), fRpr(0),fZpr(0), - fAngMomKC(0), fZvXv(0),fZvYv(0), fXvYv(0), fPtPrKink(0), f1(0), f2(0), fListOfHistos(0) - // , faod(0), fSignalTree(0) { // Constructor @@ -85,27 +63,7 @@ AliAnalysisKinkESDMC::AliAnalysisKinkESDMC(const char *name) } //________________________________________________________________________ -void AliAnalysisKinkESDMC::ConnectInputData(Option_t *) -{ - // Connect ESD or AOD here - // Called once - - TTree* tree = dynamic_cast (GetInputData(0)); - if (!tree) { - Printf("ERROR: Could not read chain from input slot 0"); - } else { - - AliESDInputHandler *esdH = dynamic_cast (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()); - - if (!esdH) { - Printf("ERROR: Could not get ESDInputHandler"); - } else - fESD = esdH->GetEvent(); - } -} - -//________________________________________________________________________ -void AliAnalysisKinkESDMC::CreateOutputObjects() +void AliAnalysisKinkESDMC::UserCreateOutputObjects() { // Create histograms // Called once @@ -123,41 +81,42 @@ void AliAnalysisKinkESDMC::CreateOutputObjects() //Open file 1= CAF //OpenFile(1); - fHistPtESD = new TH1F("fHistPtESD", "P_{T} distribution", 60, 0.0, 6.0); + fHistPtESD = new TH1F("fHistPtESD", "P_{T} distribution",100, 0.0,10.0); fHistPtESD->GetXaxis()->SetTitle("P_{T} (GeV/c)"); fHistPtESD->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)"); fHistPtESD->SetMarkerStyle(kFullCircle); - fHistPt = new TH1F("fHistPt", "P_{T} distribution", 60, 0.0, 6.0); + fHistPt = new TH1F("fHistPt", "P_{T} distribution",100, 0.0,10.0); fHistQtAll = new TH1F("fHistQtAll", "Q_{T} distr All Kinks ",100, 0.0,.250); fHistQt1= new TH1F("fHistQt1", "Q_{T} distribution",100, 0.0,.300); fHistQt2= new TH1F("fHistQt2", "Q_{T} distribution",100, 0.0,.300); - fHistPtKaon = new TH1F("fHistPtKaon", "P_{T}Kaon distribution", 60, 0.0, 6.0); - fHistPtKPDG = new TH1F("fHistPtKPDG", "P_{T}Kaon distribution", 60, 0.0, 6.0); + fHistPtKaon = new TH1F("fHistPtKaon", "P_{T}Kaon distribution",100, 0.0,10.0); + fHistPtKPDG = new TH1F("fHistPtKPDG", "P_{T}Kaon distribution",100, 0.0,10.0); fHistEta= new TH1F("fHistEta", "Eta distribution", 26,-1.3, 1.3); fHistEtaK= new TH1F("fHistEtaK", "EtaK distribution", 26,-1.3, 1.3); - fptKMC= new TH1F("fptKMC", "P_{T}Kaon generated", 60, 0.0, 6.0); - fMultiplMC= new TH1F("fMultiplMC", "charge multiplicity MC", 60, 0.5,120.5); + fptKMC= new TH1F("fptKMC", "P_{T}Kaon generated",100, 0.0,10.0); + fMultiplMC= new TH1F("fMultiplMC", "charge multiplicity MC",60, 0.5,120.5); fESDMult= new TH1F("fESDMult", "charge multipliESD", 60, 0.5,120.5); - fgenpt= new TH1F("fgenpt", "genpt K distribution", 60, 0.0, 6.0); + fgenpt= new TH1F("fgenpt", "genpt K distribution",100, 0.0,10.0); frad= new TH1F("frad", "radius K generated",100,0.,400.); - fKinkKaon= new TH1F("fKinkKaon", "P_{T}Kaon kinks identi", 60, 0.0, 6.0); - fKinkKaonBg= new TH1F("fKinkKaonBg", "P_{T}Kaon kinks backgr", 60, 0.0, 6.0); + fKinkKaon= new TH1F("fKinkKaon", "P_{T}Kaon kinks identi",100, 0.0,10.0); + fKinkKaonBg= new TH1F("fKinkKaonBg", "P_{T}Kaon kinks backgr",100, 0.0,10.0); fM1kaon= new TH1F("fM1kaon","Invar m(kaon) from kink->mu+netrino decay",80,0.0, 0.8); - fgenPtEtR= new TH1F("fgenPtEtR", "P_{T}Kaon distribution", 60, 0.0, 6.0); - fPtKink= new TH1F("fPtKink", "P_{T}Kaon Kink bution", 60, 0.0, 6.0); - fptKink= new TH1F("fptKink", "P_{T}Kaon Kink bution", 60, 0.0, 6.0); + fgenPtEtR= new TH1F("fgenPtEtR", "P_{T}Kaon distribution",100, 0.0,10.0); + fPtKink= new TH1F("fPtKink", "P_{T}Kaon Kink bution",100, 0.0,10.0); + fptKink= new TH1F("fptKink", "P_{T}Kaon Kink bution",100, 0.0,10.0); fcodeH = new TH2F("fcodeH", "code vrs dcode dist. kinks,K",100,0.,2500.,100,0.,2500.); fdcodeH = new TH2F("fdcodeH", "code vrs dcode dist. kinks,K",100,0.,2500.,100,0.,2500.); fAngMomK= new TH2F("fAngMomK","Decay angle vrs Mother Mom,K",100,0.0,5.0,80,0.,80.); fAngMomPi= new TH2F("fAngMomPi","Decay angle vrs Mother Mom,Pi",100,0.0,5.0,80,0.,80.); - + fAngMomKC= new TH2F("fAngMomKC","Decay angle vrs Mother Mom,K",100,0.0,5.0,80,0.,80.); + fMultESDK=new TH1F("fMultESDK", "charge multipliESD kaons", 60, 0.5,240.5); + fMultMCK=new TH1F("fMultMCK", "charge multipli MC kaons",100, 0.5,600.5); fRpr = new TH1D("fRpr", "rad distribution PID pr",50,0.0, 2.5); fZpr = new TH1D("fZpr", "z distribution PID pr ",60,-15.,15.); - fAngMomKC= new TH2F("fAngMomKC","Decay angle vrs Mother Mom,K",100,0.0,5.0,80,0.,80.); fZvXv= new TH2F("fZvXv","Xv-Zv main vtx",60,-0.5,0.5,60, -15., 15.0); fZvYv= new TH2F("fZvYv","Yv-Zv main vtx",60,-0.5,0.5, 60, -15., 15.); fXvYv= new TH2F("fXvYv","Xv-Yv main vtx", 60,-1.5,1.5, 60, -1.5, 1.5); - fPtPrKink=new TH1F("fPtPrKink","pt of ESD kaonKink tracks", 60, 0.0, 6.0); + fPtPrKink=new TH1F("fPtPrKink","pt of ESD kaonKink tracks",100, 0.0,10.0); fListOfHistos=new TList(); @@ -185,9 +144,11 @@ void AliAnalysisKinkESDMC::CreateOutputObjects() fListOfHistos->Add(fdcodeH); fListOfHistos->Add(fAngMomK); fListOfHistos->Add(fAngMomPi); + fListOfHistos->Add(fAngMomKC); + fListOfHistos->Add(fMultESDK); + fListOfHistos->Add(fMultMCK); fListOfHistos->Add(fRpr); fListOfHistos->Add(fZpr); - fListOfHistos->Add(fAngMomKC); fListOfHistos->Add(fZvXv); fListOfHistos->Add(fZvYv); fListOfHistos->Add(fXvYv); @@ -195,21 +156,27 @@ void AliAnalysisKinkESDMC::CreateOutputObjects() } //________________________________________________________________________ -void AliAnalysisKinkESDMC::Exec(Option_t *) +void AliAnalysisKinkESDMC::UserExec(Option_t *) { // Main loop // Called for each event // Process MC truth, therefore we receive the AliAnalysisManager and ask it for the AliMCEventHandler // This handler can return the current MC event - - AliMCEventHandler* eventHandler = dynamic_cast (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler()); - if (!eventHandler) { - Printf("ERROR: Could not retrieve MC event handler"); - return; + + AliVEvent *event = InputEvent(); + if (!event) { + Printf("ERROR: Could not retrieve event"); + return; + } + + AliESDEvent* esd = dynamic_cast(event); + if (!esd) { + Printf("ERROR: Could not retrieve esd"); + return; } - AliMCEvent* mcEvent = eventHandler->MCEvent(); + AliMCEvent* mcEvent = MCEvent(); if (!mcEvent) { Printf("ERROR: Could not retrieve MC event"); return; @@ -217,34 +184,34 @@ void AliAnalysisKinkESDMC::Exec(Option_t *) Printf("MC particles: %d", mcEvent->GetNumberOfTracks()); - if (!fESD) { - Printf("ERROR: fESD not available"); - return; - } AliStack* stack=mcEvent->Stack(); //primary tracks in MC Int_t nPrim = stack->GetNprimary(); // - const AliESDVertex *vertex=GetEventVertex(fESD); + const AliESDVertex *vertex=GetEventVertex(esd); if(!vertex) return; // Double_t vpos[3]; vertex->GetXYZ(vpos); fZpr->Fill(vpos[2]); - fZvXv->Fill( vpos[1], vpos[2]); - fZvYv->Fill( vpos[0], vpos[2]); +// fZvXv->Fill( vpos[1], vpos[2]); +// fZvYv->Fill( vpos[0], vpos[2]); Double_t vtrack[3], ptrack[3]; + // Int_t nESDTracK = 0; - Int_t nGoodTracks = fESD->GetNumberOfTracks(); + Int_t nGoodTracks = esd->GetNumberOfTracks(); fESDMult->Fill(nGoodTracks); + // // track loop - for (Int_t iTracks = 0; iTracks < fESD->GetNumberOfTracks(); iTracks++) { - AliESDtrack* track = fESD->GetTrack(iTracks); + for (Int_t iTracks = 0; iTracks < esd->GetNumberOfTracks(); iTracks++) { +// loop only on primary tracks + + AliESDtrack* track = esd->GetTrack(iTracks); if (!track) { Printf("ERROR: Could not receive track %d", iTracks); continue; @@ -253,26 +220,39 @@ void AliAnalysisKinkESDMC::Exec(Option_t *) fHistPt->Fill(track->Pt()); + + Int_t label = track->GetLabel(); + label = TMath::Abs(label); + TParticle * part = stack->Particle(label); + if (!part) continue; +// loop only on Primary tracks +// gia oles tis troxies 5/7/2009 if (label > nPrim) continue; +// if (label > nPrim) continue; /// primary tracks only + // ta nominal cauts tis apomakrynoun ols if (label < nPrim) continue; // gia tis secondary , 5/7/2009 + // pt cut at 300 MeV if ( (track->Pt())<.300)continue; - UInt_t status=track->GetStatus(); +// UInt_t status=track->GetStatus(); // if((status&AliESDtrack::kITSrefit)==0) continue; - if((status&AliESDtrack::kTPCrefit)==0) continue; - if((track->GetTPCchi2()/track->GetTPCclusters(0))>3.5) continue; + // if((status&AliESDtrack::kTPCrefit)==0) continue; + // ayksanei to ratio BG/real if((track->GetTPCchi2()/track->GetTPCclusters(0))>4.5) continue; + if((track->GetTPCchi2()/track->GetTPCclusters(0))>4.5) continue; Double_t extCovPos[15]; track->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; + // 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; track->GetXYZ(vtrack); fXvYv->Fill(vtrack[0],vtrack[1]); + fZvYv->Fill(vtrack[0],vtrack[2]); + fZvXv->Fill(vtrack[1],vtrack[2]); // track momentum track->GetPxPyPz(ptrack); @@ -298,18 +278,17 @@ void AliAnalysisKinkESDMC::Exec(Option_t *) Float_t dcaToVertexZpos = bpos[1]; fRpr->Fill(dcaToVertexZpos); - +// test for secondary kinks , 5/7/2009 if(nSigmaToVertex>=4) continue; - if((dcaToVertexXYpos>3.0)||(dcaToVertexZpos>3.0)) continue; + // if((dcaToVertexXYpos>3.0)||(dcaToVertexZpos>3.0)) continue; //arxiko + //if((dcaToVertexXYpos>0.4)||(dcaToVertexZpos>3.0)) continue;// me Zpos > 5 clen the secondaries + if((dcaToVertexXYpos>2.4)||(dcaToVertexZpos>5.0)) continue;// me Zpos > 5 clen the secondaries + // if((dcaToVertexXYpos>3.0)||(dcaToVertexZpos>3.0)) continue; // me ayto krataei 4 fakes apo ta 5 + // if((dcaToVertexXYpos>4.0)||(dcaToVertexZpos>5.0)) continue; - - Int_t label = track->GetLabel(); - label = TMath::Abs(label); - TParticle * part = stack->Particle(label); - if (!part) continue; // cut on eta - if( (TMath::Abs(trackEta )) > 1.1 ) continue; + if( (TMath::Abs(trackEta )) > 0.9 ) continue; fHistPtESD->Fill(track->Pt()); // Add Kink analysis @@ -321,7 +300,7 @@ void AliAnalysisKinkESDMC::Exec(Option_t *) // select kink class - AliESDkink *kink=fESD->GetKink(TMath::Abs(indexKinkPos)-1); + AliESDkink *kink=esd->GetKink(TMath::Abs(indexKinkPos)-1); // Int_t eSDfLabel1=kink->GetLabel(0); @@ -336,24 +315,34 @@ void AliAnalysisKinkESDMC::Exec(Option_t *) const TVector3 daughterMKink(kink->GetDaughterP()); Float_t qT=kink->GetQt(); - fHistEta->Fill(trackEta) ; // Eta distr + //fHistEta->Fill(trackEta) ; // Eta distr fHistQtAll->Fill(qT) ; // Qt distr fptKink->Fill(motherMfromKink.Pt()); /// pt from kink + Float_t kinkAngle=TMath::RadToDeg()*kink->GetAngle(2); + // fake kinks, small Qt and small kink angle + if(( kinkAngle<1.)&&(TMath::Abs(code1)==321)&&(TMath::Abs(dcode1)==13)) fHistQt1 ->Fill(qT) ; // Qt distr if ( qT<0.012) fcodeH->Fill( TMath::Abs(code1), TMath::Abs(dcode1) ); if(qT<0.012) continue; // remove fake kinks + if( (kinkAngle<1.) ) continue; - Float_t kinkAngle=TMath::RadToDeg()*kink->GetAngle(2); - //remove the double taracks - if( (kinkAngle<1.) ) continue; + //if( (kinkAngle<1.) ) continue; // - if( (qT>0.05)&& ( qT<0.25) ) fHistQt2->Fill(qT); // candidate kaon kinks - if (TMath::Abs(code1)==321) fHistPtKPDG->Fill(track->Pt()); // ALL KAONS (pdg) inside kink sample + if((kink->GetR()>120.)&&(kink->GetR()<220.)&&(TMath::Abs(trackEta)<0.9)&&(labelFill(track->Pt()); // ALL KAONS (pdg) inside ESD kink sample + fHistEta->Fill(trackEta) ; // Eta distr of PDG kink ESD kaons + fMultESDK->Fill(nGoodTracks); + } + } + if( (qT>0.05)&& ( qT<0.25) ) fHistQt2->Fill(qT); // candidate kaon kinks // maximum decay angle at a given mother momentum Double_t maxDecAngKmu=f1->Eval(motherMfromKink.Mag(),0.,0.,0.); @@ -376,23 +365,29 @@ void AliAnalysisKinkESDMC::Exec(Option_t *) fM1kaon->Fill(invariantMassKmu); // kaon selection from kinks -if((kinkAngle>maxDecAngpimu)&&(qT>0.05)&&(qT<0.25)&&((kink->GetR()>110.)&&(kink->GetR()<230.))&&(TMath::Abs(trackEta)<1.1)&&(invariantMassKmu<0.6)) { - fHistEtaK->Fill(trackEta); +if((kinkAngle>maxDecAngpimu)&&(qT>0.05)&&(qT<0.25)&&((kink->GetR()>120.)&&(kink->GetR()<220.))&&(TMath::Abs(trackEta)<0.9)&&(invariantMassKmu<0.6)) { + // fHistEtaK->Fill(trackEta); - if(TMath::Abs(code1)==321) fAngMomKC->Fill(motherMfromKink.Mag(), kinkAngle); // real kaons + // if(TMath::Abs(code1)==321) fAngMomKC->Fill(motherMfromKink.Mag(), kinkAngle); // real kaons - fHistQt1 ->Fill(qT) ; // Qt distr +// fHistQt1 ->Fill(qT) ; // Qt distr if ( (kinkAngle>maxDecAngKmu) && ( motherMfromKink.Mag()> 1.1 ) ) continue; // maximum angle selection + if(TMath::Abs(code1)==321) fAngMomKC->Fill(motherMfromKink.Mag(), kinkAngle); // real kaons fHistPtKaon->Fill(track->Pt()); //all PID kink-kaon // background inside the identified kaons, e.g KK , Kp if((TMath::Abs(code1)==321)&&(( TMath::Abs(dcode1)) >=(TMath::Abs(code1)) )) fdcodeH->Fill( TMath::Abs(code1), TMath::Abs(dcode1)); - - if((TMath::Abs(code1)==321) ) { +// kaons from k to mun and k to pipi and to e decay + if( ( (TMath::Abs(code1)==321)&&(TMath::Abs(dcode1)==13))|| + ( (TMath::Abs(code1)==321)&&(TMath::Abs(dcode1)==11)) || + ( (TMath::Abs(code1)==321)&&(TMath::Abs(dcode1)==211)) ) { + // if((TMath::Abs(code1)==321) ) { if ( label< nPrim ) fPtPrKink->Fill(track->Pt()); - fKinkKaon->Fill(track->Pt()); } + fKinkKaon->Fill(track->Pt()); + fHistEtaK->Fill(trackEta); + } else { fKinkKaonBg->Fill(track->Pt()); } // primary and secondary @@ -409,8 +404,9 @@ if((kinkAngle>maxDecAngpimu)&&(qT>0.05)&&(qT<0.25)&&((kink->GetR()>110.)&&(kink- // variable to count tracks Int_t nMCTracks = 0; + Int_t nMCKinkKs = 0; -for (Int_t iMc = 0; iMc < nPrim; ++iMc) +for (Int_t iMc = 0; iMc < nPrim; iMc++) { // AliDebug(AliLog::kDebug+1, Form("MC Loop: Processing particle %d.", iMc)); @@ -425,7 +421,7 @@ for (Int_t iMc = 0; iMc < nPrim; ++iMc) // Float_t eta = particle->Eta(); - if (TMath::Abs(eta) > 1.1) + if (TMath::Abs(eta) > 0.9) continue; Double_t ptK = particle->Pt(); @@ -433,10 +429,14 @@ for (Int_t iMc = 0; iMc < nPrim; ++iMc) if( ptK <0.300) continue; Float_t code = particle->GetPdgCode(); - if (code==321 || code==-321){ + Int_t mcProcess=-1011; + if ((code==321)||(code==-321)){ + fptKMC->Fill(ptK); + + // fMultiplMC->Fill(nPrim); Int_t firstD=particle->GetFirstDaughter(); Int_t lastD=particle->GetLastDaughter(); @@ -445,27 +445,33 @@ for (Int_t iMc = 0; iMc < nPrim; ++iMc) TParticle* daughter1=stack->Particle(k+1); Float_t dcode = daughter1->GetPdgCode(); + mcProcess=daughter1->GetUniqueID(); + if (mcProcess==4) { - frad->Fill(daughter1->R()); + // frad->Fill(daughter1->R()); - if (((daughter1->R())>110)&&((daughter1->R())<230) ){ + if (((daughter1->R())>120)&&((daughter1->R())<220) ){ if (( ( code==321 )&& ( dcode ==-13 ))|| (( code == -321 )&& ( dcode == 13))) fgenPtEtR->Fill( ptK );//to muon if (( ( code==321 )&& ( dcode ==211 ))|| (( code == -321 )&& ( dcode ==-211))) fgenPtEtR->Fill( ptK );//to pion - // if (( ( code==321 )&& ( dcode ==-11 ))|| (( code == -321 )&& ( dcode ==11))) fgenPtEtR->Fill( ptK );// to electr + if (( ( code==321 )&& ( dcode ==-11 ))|| (( code == -321 )&& ( dcode ==11))) fgenPtEtR->Fill( ptK );// to electr + //fMultMCK->Fill(nPrim); + frad->Fill(daughter1->R()); + nMCKinkKs++; } // next only to mu decay if (((code==321)&&(dcode==-13))||((code==-321)&&(dcode==13))){ - if (((daughter1->R())>110)&&((daughter1->R())<230)){ + if (((daughter1->R())>120)&&((daughter1->R())<220)){ fgenpt->Fill(ptK); } } - } + }// decay + }// daughters - } + } /// kaons @@ -474,6 +480,7 @@ for (Int_t iMc = 0; iMc < nPrim; ++iMc) // printf("hello mc"); fMultiplMC->Fill(nMCTracks); + if( nMCKinkKs>0) fMultMCK->Fill(nPrim); PostData(1, fListOfHistos); diff --git a/PWG2/KINK/AliAnalysisKinkESDMC.h b/PWG2/KINK/AliAnalysisKinkESDMC.h index b04799d2333..6f36a711716 100644 --- a/PWG2/KINK/AliAnalysisKinkESDMC.h +++ b/PWG2/KINK/AliAnalysisKinkESDMC.h @@ -13,68 +13,68 @@ // mspyrop@phys.uoa.gr //----------------------------------------------------------------- -class AliESDVertex; class AliESDEvent; +class AliESDVertex; class AliESDtrack; -class TTree; class TF1; class TH1F; class TH2F; class TH1D; class TH2D; +class TList; #include "AliAnalysisTaskSE.h" class AliAnalysisKinkESDMC : public AliAnalysisTaskSE { public: - AliAnalysisKinkESDMC(); - AliAnalysisKinkESDMC(const char *name); + // AliAnalysisKinkESDMC(); + AliAnalysisKinkESDMC(const char *name = "AliAnalysisKinkESDMC"); virtual ~AliAnalysisKinkESDMC() {} - - virtual void ConnectInputData(Option_t *); - virtual void CreateOutputObjects(); - virtual void Exec(Option_t *option); + + virtual void UserCreateOutputObjects(); + virtual void UserExec(Option_t *option); virtual void Terminate(Option_t *); - + Float_t GetSigmaToVertex(AliESDtrack* esdTrack) const; const AliESDVertex *GetEventVertex(const AliESDEvent* esd) const; private: - AliESDEvent *fESD; //!ESD object - TH1F *fHistPtESD; //!Pt spectrum of all ESD inside eta, Pt cuts - TH1F *fHistPt; //!Pt spectrum of all ESD tracks - TH1F *fHistQtAll; //!Qt spectrum of all kinks - TH1F *fHistQt1; //!Qt spectrum of Kaon selected sample - TH1F *fHistQt2; //!Qt spectrum in Qt region of kaons - TH1F *fHistPtKaon; //!Pt Kaon spectrum of clean sample - TH1F *fHistPtKPDG; //!Pt Kaon spectrum , confirmed by PDG,inside kaon Qt region - TH1F *fHistEta; //!Eta spectrum of all kinks - TH1F *fHistEtaK; //!Eta spectrum of kaons selected by kink topology - TH1F *fptKMC; //!Pt Kaon spectrum MC, inside eta and pt cuts - TH1F *fMultiplMC; //!charge multipl MC - TH1F *fESDMult; //!ESD charged mult - TH1F *fgenpt; //!Pt Kaon-Kink->mu spectrum , MC, inside eta, Pt, radius cuts - TH1F *frad; //!radius of kinks, MC , inside the eta nad Pt cuts - TH1F *fKinkKaon; //!Pt of PDG Kaons inside the selcted ones by the KInk topology - TH1F *fKinkKaonBg; //!Pt of the BG inside the kink-Kaon identified spectrum - TH1F *fM1kaon; //!inv mass of kink-tracks taken as kaons decaying to mu + neutrino - TH1F *fgenPtEtR; //!MC Pt spectrum of kaons decaying to muon+neutrino and pi +pi, inside eta,Pt,Rad cuts - TH1F *fPtKink; //!Pt spectrum of all kinks from track bank - TH1F *fptKink; //!Pt spectrum of all kinks from kink bank - TH2F *fcodeH ; //!PDG code(mother) vrs PDG dcode(daughter) of kinks with Qt <0.12 (fake) - TH2F *fdcodeH ; //!Kinks, code vrs dcode of BG,if mother code is 321 and daughter code > - TH2F *fAngMomK; //! Decay angle vrs Mother Mom for pdg kaons - TH2F *fAngMomPi; //! Decay angle vrs Mother Mom for pdg pions - TH1D *fRpr; //! Radius of VTX at Y , X plane - TH1D *fZpr; //! Z distrio of main vertex - TH2F *fAngMomKC; //! Decay angle vrs Mother Mom for pdg kaons, inside the selected sample - TH2F *fZvXv; //! two dime of Z vrs X of vtx main - TH2F *fZvYv; //! two dime of Z vrs Y of vtx main - TH2F *fXvYv; //! two dime of X vrs Y of main tracks vtx main - TH1F *fPtPrKink; //! pt of Primary PDG kaons inside the selected ones by the kink topology + TH1F *fHistPtESD; //Pt spectrum of all ESD inside eta, Pt cuts + TH1F *fHistPt; //Pt spectrum of all ESD tracks + TH1F *fHistQtAll; //Qt spectrum of all kinks + TH1F *fHistQt1; //Qt spectrum of Kaon selected sample + TH1F *fHistQt2; //Qt spectrum in Qt region of kaons + TH1F *fHistPtKaon; //Pt Kaon spectrum of clean sample + TH1F *fHistPtKPDG; //Pt Kaon spectrum , confirmed by PDG,inside kaon Qt region + TH1F *fHistEta; //Eta spectrum of all kinks + TH1F *fHistEtaK; //Eta spectrum of kaons selected by kink topology + TH1F *fptKMC; //Pt Kaon spectrum MC, inside eta and pt cuts + TH1F *fMultiplMC; //charge multipl MC + TH1F *fESDMult; //ESD charged mult + TH1F *fgenpt; //Pt Kaon-Kink->mu spectrum , MC, inside eta, Pt, radius cuts + TH1F *frad; //radius of kinks, MC , inside the eta nad Pt cuts + TH1F *fKinkKaon; //Pt of PDG Kaons inside the selcted ones by the KInk topology + TH1F *fKinkKaonBg; //Pt of the BG inside the kink-Kaon identified spectrum + TH1F *fM1kaon; //inv mass of kink-tracks taken as kaons decaying to mu + neutrino + TH1F *fgenPtEtR; //MC Pt spectrum of kaons decaying to muon+neutrino and pi +pi, inside eta,Pt,Rad cuts + TH1F *fPtKink; //Pt spectrum of all kinks from track bank + TH1F *fptKink; //Pt spectrum of all kinks from kink bank + TH2F *fcodeH ; //PDG code(mother) vrs PDG dcode(daughter) of kinks with Qt <0.12 (fake) + TH2F *fdcodeH ; //inks, code vrs dcode of BG,if mother code is 321 and daughter code > + TH2F *fAngMomK; // Decay angle vrs Mother Mom for pdg kaons + TH2F *fAngMomPi; // Decay angle vrs Mother Mom for pdg pions + TH2F *fAngMomKC; //Decay angle vrs Mother Mom for pdg kaons, inside the selected sample + TH1F *fMultESDK; //ESD charged mult + TH1F *fMultMCK; //MC K charged mult + TH1D *fRpr; // Radius of VTX at Y , X plane + TH1D *fZpr; //Z distrio of main vertex + TH2F *fZvXv; //two dime of Z vrs X of vtx main + TH2F *fZvYv; // two dime of Z vrs Y of vtx main + TH2F *fXvYv; // two dime of X vrs Y of main tracks vtx main + TH1F *fPtPrKink; // pt of Primary PDG kaons inside the selected ones by the kink topology TF1 *f1; TF1 *f2; - TList *fListOfHistos; //! list of histos + TList *fListOfHistos; // list of histos AliAnalysisKinkESDMC(const AliAnalysisKinkESDMC&); // not implemented AliAnalysisKinkESDMC& operator=(const AliAnalysisKinkESDMC&); // not implemented diff --git a/PWG2/KINK/AliAnalysisTaskKinkResonance.cxx b/PWG2/KINK/AliAnalysisTaskKinkResonance.cxx index c71c2e192f5..082331295e0 100644 --- a/PWG2/KINK/AliAnalysisTaskKinkResonance.cxx +++ b/PWG2/KINK/AliAnalysisTaskKinkResonance.cxx @@ -14,94 +14,69 @@ //---------------------------------------------------------------------------------------------------------------- // class AliAnalysisTaskKinkResonance // Example of an analysis task for reconstructing resonances having at least one kaon-kink in their decay -// products. It provides basic plots as well as plots helping to calculate the corrections. -// Usage: To analyse a resonance having a kaon in its decay products, one has to modify the integer -// variables resonancePDG, daughter1 and daughter2 accordingly as well as daughter1Mass and daughter2Mass. -// Also, depending on the analysis mode (ESD or MC), fAnalysisType in the constructor must also be changed +// products. //----------------------------------------------------------------------------------------------------------------- -#include "TCanvas.h" #include "TChain.h" #include "TTree.h" #include "TH2D.h" -#include "AliAnalysisTaskSE.h" #include "AliAnalysisManager.h" -#include "AliESDInputHandler.h" -#include "AliMCEventHandler.h" #include "AliMCEvent.h" +#include "AliVEvent.h" +#include "AliESDEvent.h" #include "AliResonanceKink.h" #include "AliAnalysisTaskKinkResonance.h" ClassImp(AliAnalysisTaskKinkResonance) -//_____________________________________________________________________________ -AliAnalysisTaskKinkResonance::AliAnalysisTaskKinkResonance() - : AliAnalysisTaskSE(), fESD(0), fmcEventH(0), fList(0), fKinkResonance(0) - -{ - // Constructor -} - //______________________________________________________________________________ -AliAnalysisTaskKinkResonance::AliAnalysisTaskKinkResonance(const char *name) - : AliAnalysisTaskSE(name), fESD(0), fmcEventH(0), fList(0), fKinkResonance(0) +AliAnalysisTaskKinkResonance::AliAnalysisTaskKinkResonance(const char *dname) + : AliAnalysisTaskSE(dname), fList(0), fKinkResonance(0) { // Constructor // Define input and output slots here // Input slot #0 works with a TChain - DefineInput(0, TChain::Class()); + //DefineInput(0, TChain::Class()); DefineOutput(1, TList::Class()); } //________________________________________________________________________ -void AliAnalysisTaskKinkResonance::ConnectInputData(Option_t *) -{ - // Connect ESD or AOD here - // Called once - - TTree* tree = dynamic_cast (GetInputData(0)); - if (!tree) { - Printf("ERROR: Could not read chain from input slot 0"); - } else { - - AliESDInputHandler *esdH = dynamic_cast (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()); - - if (!esdH) { - Printf("ERROR: Could not get ESDInputHandler"); - } else - fESD = esdH->GetEvent(); - - AliMCEventHandler* eventHandler = dynamic_cast (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler()); - - if (!eventHandler) { - Printf("ERROR: Could not retrieve MC event handler"); - return; - } - - fmcEventH = eventHandler->MCEvent(); - } -} - -//________________________________________________________________________ -void AliAnalysisTaskKinkResonance::CreateOutputObjects() +void AliAnalysisTaskKinkResonance::UserCreateOutputObjects() { // Create histograms // Called once - fList=new TList(); fList= fKinkResonance->GetHistogramList(); } //________________________________________________________________________ -void AliAnalysisTaskKinkResonance::Exec(Option_t *) +void AliAnalysisTaskKinkResonance::UserExec(Option_t *) { // Main loop // Called for each event - fKinkResonance->Analyse(fESD, fmcEventH); + AliVEvent *event = InputEvent(); + if (!event) { + Printf("ERROR: Could not retrieve event"); + return; + } + + AliESDEvent* esd = dynamic_cast(event); + if (!esd) { + Printf("ERROR: Could not retrieve esd"); + return; + } + + AliMCEvent* mcEvent = MCEvent(); + if (!mcEvent) { + Printf("ERROR: Could not retrieve MC event"); + return; + } + + fKinkResonance->Analyse(esd, mcEvent); PostData(1, fList); diff --git a/PWG2/KINK/AliAnalysisTaskKinkResonance.h b/PWG2/KINK/AliAnalysisTaskKinkResonance.h index d62482dc554..e33a85fe554 100644 --- a/PWG2/KINK/AliAnalysisTaskKinkResonance.h +++ b/PWG2/KINK/AliAnalysisTaskKinkResonance.h @@ -11,29 +11,25 @@ //------------------------------------------------------------------------------ class TList; -class AliESDEvent; -class AliMCEvent; class AliResonanceKink; class TH1D; +#include "AliAnalysisTaskSE.h" + class AliAnalysisTaskKinkResonance : public AliAnalysisTaskSE { public: - AliAnalysisTaskKinkResonance(); - AliAnalysisTaskKinkResonance(const char *name); + AliAnalysisTaskKinkResonance(const char *dname = "AliAnalysisTaskKinkResonance"); virtual ~AliAnalysisTaskKinkResonance() {} - virtual void ConnectInputData(Option_t *); - virtual void CreateOutputObjects(); - virtual void Exec(Option_t *option); + virtual void UserCreateOutputObjects(); + virtual void UserExec(Option_t *option); virtual void Terminate(Option_t *); void SetAnalysisKinkObject(AliResonanceKink * const kinkResonance) { fKinkResonance=kinkResonance;} private: - AliESDEvent *fESD; //! ESD object - AliMCEvent *fmcEventH; - TList *fList; //! List + TList *fList; // List AliResonanceKink *fKinkResonance; AliAnalysisTaskKinkResonance(const AliAnalysisTaskKinkResonance&); // not implemented diff --git a/PWG2/KINK/AliResonanceKink.cxx b/PWG2/KINK/AliResonanceKink.cxx index 6f8cf49bc47..c155308c8c7 100644 --- a/PWG2/KINK/AliResonanceKink.cxx +++ b/PWG2/KINK/AliResonanceKink.cxx @@ -24,9 +24,8 @@ #include "TTree.h" #include "TH2D.h" #include "TParticle.h" -#include -#include -#include +#include "TDatabasePDG.h" +#include "TParticlePDG.h" #include "TF1.h" #include "TList.h" #include "TString.h" @@ -37,13 +36,15 @@ #include "AliStack.h" #include "AliESDtrack.h" #include "AliESDEvent.h" +#include "AliExternalTrackParam.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) + 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) { // Constructor @@ -52,7 +53,8 @@ AliResonanceKink::AliResonanceKink() //________________________________________________________________________ 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) + 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) { // Constructor @@ -86,7 +88,7 @@ AliResonanceKink::AliResonanceKink(Int_t nbins, Float_t nlowx, Float_t nhighx, I f2->SetParameter(2,TMath::Pi()); fvtxz=new TH1D("fvtxz"," ", 100,-20.0,20.0); - + } //________________________________________________________________________ @@ -273,32 +275,17 @@ void AliResonanceKink::Analyse(AliESDEvent* esd, AliMCEvent* mcEvent) } if (trackpos->GetSign() < 0) continue; - trackpos->GetPxPyPz(ptrackpos); - - 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]; + AliExternalTrackParam *tpcTrackpos = (AliExternalTrackParam *)trackpos->GetTPCInnerParam(); + if(!tpcTrackpos) continue; + ptrackpos[0]=tpcTrackpos->Px(); + ptrackpos[1]=tpcTrackpos->Py(); + ptrackpos[2]=tpcTrackpos->Pz(); - 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(); @@ -309,57 +296,16 @@ void AliResonanceKink::Analyse(AliESDEvent* esd, AliMCEvent* mcEvent) 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 (kinkAnglePosGetStatus(); - 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); } @@ -369,29 +315,16 @@ void AliResonanceKink::Analyse(AliESDEvent* esd, AliMCEvent* mcEvent) 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; @@ -403,66 +336,23 @@ void AliResonanceKink::Analyse(AliESDEvent* esd, AliMCEvent* mcEvent) 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 (kinkAngleNegGetStatus(); - - 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)) { + if((posKaonKinkFlag==1)&&(negKaonKinkFlag==1)) { p4comb=anp4pos; p4comb+=p4neg; if(openingAngle>0.6) fPhiBothKinks->Fill(p4comb.M()); @@ -487,7 +377,7 @@ void AliResonanceKink::Analyse(AliESDEvent* esd, AliMCEvent* mcEvent) } - if(kaonKinkFlag==1) { + if(posKaonKinkFlag==1) { anp4comb=anp4pos; anp4comb+=anp4neg; fInvariantMass->Fill(anp4comb.M()); @@ -521,7 +411,7 @@ Float_t AliResonanceKink::GetSigmaToVertex(AliESDtrack* esdTrack) const { 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; @@ -558,5 +448,165 @@ const AliESDVertex* AliResonanceKink::GetEventVertex(const AliESDEvent* esd) con } } +//________________________________________________________________________ + + 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) { + 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) { + 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) { + 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) { + 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) { + Printf("IsAccepted: Track rejected because of no refited in TPC"); + return kFALSE; + } + + if(nClustersTPC < fMinTPCclusters) { + Printf("IsAccepted: Track rejected because it has a value of nclusters (TPC) of %ld (min. requested: %ld)", nClustersTPC, fMinTPCclusters); + return kFALSE; + } + + if(chi2perTPCcluster > fMaxChi2PerTPCcluster) { + Printf("IsAccepted: Track rejected because it has a value of chi2perTPCcluster of %lf (max. requested: %lf)", chi2perTPCcluster, fMaxChi2PerTPCcluster); + return kFALSE; + } + + if(extCov[0] > fMaxCov0) { + 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) { + 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) { + 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) { + 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) { + 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()>110.)&&(kink->GetR()<230.))&&(TMath::Abs(trackMom.Eta())<1.1)&&(invariantMassKmu<0.6)) { + + if(trackMom.Mag()<=1.1) { + return kTRUE; + } + else + if (kinkAngleGetMCtruthEventHandler()) { kinkResonanceObject->InitOutputHistograms(60, 0.6, 1.2); kinkResonanceObject->SetPDGCodes(kKPlus, kPiPlus, AliResonanceKink::kKstar0); kinkResonanceObject->SetAnalysisType("ESD"); // "ESD" or "MC" + kinkResonanceObject->SetMaxNsigmaToVertex(4.0); + kinkResonanceObject->SetMaxDCAxy(3.0); + kinkResonanceObject->SetMaxDCAzaxis(3.0); + kinkResonanceObject->SetPtTrackCut(0.25); + kinkResonanceObject->SetMinTPCclusters(50); + kinkResonanceObject->SetMaxChi2PerTPCcluster(3.5); + kinkResonanceObject->SetMaxCov0(2.0); + kinkResonanceObject->SetMaxCov2(2.0); + kinkResonanceObject->SetMaxCov5(0.5); + kinkResonanceObject->SetMaxCov9(0.5); + kinkResonanceObject->SetMaxCov14(2.0); // Create and configure the task AliAnalysisTaskKinkResonance *taskresonanceKstar = new AliAnalysisTaskKinkResonance("TaskResKinkPID"); -- 2.43.0