//----------------------------------------------------------------------------
// Implementation of the D0toKpi reconstruction and analysis class
-//
// Note: the two decay tracks are labelled: 0 (positive track)
// 1 (negative track)
-//
-// Origin: A. Dainese andrea.dainese@pd.infn.it
+// An example of usage can be found in the macro AliD0toKpiTest.C
+// Origin: A. Dainese andrea.dainese@lnl.infn.it
//----------------------------------------------------------------------------
-#include <Riostream.h>
-#include <TH1.h>
-#include <TH2.h>
#include <TKey.h>
#include <TFile.h>
#include <TTree.h>
#include <TSystem.h>
#include <TParticle.h>
#include "AliESD.h"
-#include "AliStack.h"
+#include "AliMC.h"
+#include "AliRun.h"
#include "AliRunLoader.h"
-#include "AliITStrackV2.h"
-#include "AliITSVertexerTracks.h"
-#include "AliITSVertex.h"
-#include "AliV0vertexer.h"
-#include "AliV0vertex.h"
+#include "AliVertexerTracks.h"
+#include "AliESDVertex.h"
+#include "AliESDv0.h"
#include "AliD0toKpi.h"
#include "AliD0toKpiAnalysis.h"
+#include "AliLog.h"
typedef struct {
Int_t lab;
Int_t pdg;
Int_t mumlab;
Int_t mumpdg;
+ Int_t mumprongs;
Float_t Vx,Vy,Vz;
Float_t Px,Py,Pz;
} REFTRACK;
AliD0toKpiAnalysis::AliD0toKpiAnalysis() {
// Default constructor
- SetBz();
SetPtCut();
Setd0Cut();
SetMassCut();
fVertexOnTheFly = kFALSE;
fSim = kFALSE;
fOnlySignal = kFALSE;
- fDebug = kFALSE;
}
//----------------------------------------------------------------------------
AliD0toKpiAnalysis::~AliD0toKpiAnalysis() {}
const Char_t *outName) {
// Find D0 candidates and calculate parameters
- if(fBz<-9000.) {
- printf("AliD0toKpiAnalysis::FindCandidates(): Set B!\n");
- return;
- }
- AliKalmanTrack::SetConvConst(100/0.299792458/fBz);
-
- TString trkName("AliITStracksV2.root");
- if(gSystem->AccessPathName(trkName.Data(),kFileExists)) {
- printf("AliD0toKpiAnalysis::FindCandidates(): No tracks file found!\n");
- return;
- }
-
- TString outName1=outName;
- TString outName2("nTotEvents.dat");
-
- Int_t ev;
- Int_t nTotEv=0,nD0rec=0,nD0rec1ev=0;
- Double_t dca;
- Double_t v2[3],mom[6],d0[2];
- Double_t alphaP,alphaN,ptP,ptN,phiP,phiN;
- Int_t iTrkP,iTrkN,trkEntries;
- Int_t nTrksP=0,nTrksN=0;
- Int_t trkNum[2];
- Int_t okD0=0,okD0bar=0;
- Char_t trkTreeName[100],vtxName[100];
- AliITStrackV2 *postrack = 0;
- AliITStrackV2 *negtrack = 0;
-
- // create the AliITSVertexerTracks object
- // (it will be used only if fVertexOnTheFly=kTrue)
- AliITSVertexerTracks *vertexer1 = new AliITSVertexerTracks;
- vertexer1->SetMinTracks(2);
- vertexer1->SetDebug(0);
- Int_t skipped[2];
- Bool_t goodVtx1;
-
-
- // define the cuts for vertexing
- Double_t vtxcuts[]={50., // max. allowed chi2
- 0.0, // min. allowed negative daughter's impact param
- 0.0, // min. allowed positive daughter's impact param
- 1.0, // max. allowed DCA between the daughter tracks
- -1.0, // min. allowed cosine of pointing angle
- 0.0, // min. radius of the fiducial volume
- 2.9};// max. radius of the fiducial volume
-
- // create the AliV0vertexer object
- AliV0vertexer *vertexer2 = new AliV0vertexer(vtxcuts);
-
- // create tree for reconstructed D0s
- AliD0toKpi *ioD0toKpi=0;
- TTree *treeD0 = new TTree("TreeD0","Tree with D0 candidates");
- treeD0->Branch("D0toKpi","AliD0toKpi",&ioD0toKpi,200000,0);
-
- // open file with tracks
- TFile *trkFile = TFile::Open(trkName.Data());
-
- // loop on events in file
- for(ev=evFirst; ev<=evLast; ev++) {
- printf(" --- Looking for D0->Kpi in event %d ---\n",ev);
-
- // retrieve primary vertex from file
- sprintf(vtxName,"VertexTracks_%d",ev);
- AliITSVertex *vertex1stored = (AliITSVertex*)trkFile->Get(vtxName);
- vertex1stored->GetXYZ(fV1);
- delete vertex1stored;
-
- // retrieve tracks from file
- sprintf(trkTreeName,"TreeT_ITS_%d",ev);
- TTree *trkTree=(TTree*)trkFile->Get(trkTreeName);
- if(!trkTree) {
- printf("AliD0toKpiAnalysis::FindCandidates():\n tracks tree not found for evet %d\n",ev);
- continue;
- }
- trkEntries = (Int_t)trkTree->GetEntries();
- printf(" Number of tracks: %d\n",trkEntries);
-
- // count the total number of events
- nTotEv++;
-
- // call function which applies sigle-track selection and
- // separetes positives and negatives
- TObjArray trksP(trkEntries/2);
- Int_t *trkEntryP = new Int_t[trkEntries];
- TObjArray trksN(trkEntries/2);
- Int_t *trkEntryN = new Int_t[trkEntries];
- SelectTracks(*trkTree,trksP,trkEntryP,nTrksP,trksN,trkEntryN,nTrksN);
-
- nD0rec1ev = 0;
-
- // loop on positive tracks
- for(iTrkP=0; iTrkP<nTrksP; iTrkP++) {
- if(iTrkP%10==0) printf(" Processing positive track number %d of %d\n",iTrkP,nTrksP);
-
- // get track from track array
- postrack = (AliITStrackV2*)trksP.At(iTrkP);
- trkNum[0] = trkEntryP[iTrkP];
-
- // loop on negative tracks
- for(iTrkN=0; iTrkN<nTrksN; iTrkN++) {
- // get track from tracks array
- negtrack = (AliITStrackV2*)trksN.At(iTrkN);
- trkNum[1] = trkEntryN[iTrkN];
- AliITStrackV2 nt(*negtrack), pt(*postrack), *pnt=&nt, *ppt=&pt;
-
- //
- // ----------- DCA MINIMIZATION ------------------
- //
- // find the DCA and propagate the tracks to the DCA
- dca = vertexer2->PropagateToDCA(pnt,ppt);
-
- // define the AliV0vertex object
- AliV0vertex *vertex2 = new AliV0vertex(*pnt,*ppt);
-
- // get position of the secondary vertex
- vertex2->GetXYZ(v2[0],v2[1],v2[2]);
-
- delete vertex2;
-
- // momenta of the tracks at the vertex
- ptP = 1./TMath::Abs(ppt->Get1Pt());
- alphaP = ppt->GetAlpha();
- phiP = alphaP+TMath::ASin(ppt->GetSnp());
- mom[0] = ptP*TMath::Cos(phiP);
- mom[1] = ptP*TMath::Sin(phiP);
- mom[2] = ptP*ppt->GetTgl();
-
- ptN = 1./TMath::Abs(pnt->Get1Pt());
- alphaN = pnt->GetAlpha();
- phiN = alphaN+TMath::ASin(pnt->GetSnp());
- mom[3] = ptN*TMath::Cos(phiN);
- mom[4] = ptN*TMath::Sin(phiN);
- mom[5] = ptN*pnt->GetTgl();
-
- goodVtx1 = kTRUE;
- // no vertexing if DeltaMass > fMassCut
- if(fVertexOnTheFly) {
- goodVtx1 = kFALSE;
- if(SelectInvMass(mom)) {
- // primary vertex from *other* tracks in event
- vertexer1->SetVtxStart(fV1[0],fV1[1]);
- skipped[0] = trkEntryP[iTrkP];
- skipped[1] = trkEntryN[iTrkN];
- vertexer1->SetSkipTracks(2,skipped);
- AliITSVertex *vertex1onfly =
- (AliITSVertex*)vertexer1->VertexOnTheFly(*trkTree);
- if(vertex1onfly->GetNContributors()>0) goodVtx1 = kTRUE;
- vertex1onfly->GetXYZ(fV1);
- //vertex1onfly->PrintStatus();
- delete vertex1onfly;
- }
- }
-
- // impact parameters of the tracks w.r.t. the primary vertex
- d0[0] = 10000.*ppt->GetD(fV1[0],fV1[1]);
- d0[1] = -10000.*pnt->GetD(fV1[0],fV1[1]);
-
- // create the object AliD0toKpi
- AliD0toKpi theD0(ev,trkNum,fV1,v2,dca,mom,d0);
-
- // select D0s
- if(goodVtx1 && theD0.Select(fD0Cuts,okD0,okD0bar)) {
-
- // fill the tree
- ioD0toKpi=&theD0;
- treeD0->Fill();
-
- nD0rec++; nD0rec1ev++;
- ioD0toKpi=0;
- }
-
- negtrack = 0;
- } // loop on negative tracks
- postrack = 0;
- } // loop on positive tracks
-
- trksP.Delete();
- trksN.Delete();
- delete [] trkEntryP;
- delete [] trkEntryN;
- delete trkTree;
-
- printf(" Number of D0 candidates: %d\n",nD0rec1ev);
- } // loop on events in file
-
-
- printf("\n+++\n+++ Total number of events: %d\n+++\n",nTotEv);
- printf("\n+++\n+++ Total number of D0 candidates: %d\n+++\n",nD0rec);
-
- delete vertexer1;
- delete vertexer2;
-
- trkFile->Close();
-
- // create a copy of this class to be written to output file
- AliD0toKpiAnalysis *copy = (AliD0toKpiAnalysis*)this->Clone("D0toKpiAnalysis");
- copy->PrintStatus();
-
- // add PDG codes to decay tracks in found candidates (in simulation mode)
- // and store tree in the output file
- if(!fSim) {
- TFile *outroot = new TFile(outName1.Data(),"recreate");
- treeD0->Write();
- copy->Write();
- outroot->Close();
- delete outroot;
- } else {
- printf(" Now adding information from simulation (PDG codes) ...\n");
- TTree *treeD0sim = new TTree("TreeD0","Tree with D0 candidates");
- MakeTracksRefFile(evFirst,evLast);
- SimulationInfo(treeD0,treeD0sim);
- delete treeD0;
- TFile *outroot = new TFile(outName1.Data(),"recreate");
- treeD0sim->Write();
- copy->Write();
- outroot->Close();
- delete outroot;
- }
-
- // write to a file the total number of events
- FILE *outdat = fopen(outName2.Data(),"w");
- fprintf(outdat,"%d\n",nTotEv);
- fclose(outdat);
-
- return;
-}
-//----------------------------------------------------------------------------
-void AliD0toKpiAnalysis::FindCandidatesESD(Int_t evFirst,Int_t evLast,
- const Char_t *outName) {
- // Find D0 candidates and calculate parameters
-
- if(fBz<-9000.) {
- printf("AliD0toKpiAnalysis::FindCandidatesESD(): Set B!\n");
- return;
- }
- AliKalmanTrack::SetConvConst(100/0.299792458/fBz);
-
- TString esdName("AliESDs.root");
+ TString esdName="AliESDs.root";
if(gSystem->AccessPathName(esdName.Data(),kFileExists)) {
printf("AliD0toKpiAnalysis::FindCandidatesESD(): No ESDs file found!\n");
return;
}
TString outName1=outName;
- TString outName2("nTotEvents.dat");
+ TString outName2="nTotEvents.dat";
- Double_t covV1[6];
Int_t nTotEv=0,nD0rec=0,nD0rec1ev=0;
Double_t dca;
Double_t v2[3],mom[6],d0[2];
- Double_t alphaP,alphaN,ptP,ptN,phiP,phiN;
Int_t iTrkP,iTrkN,trkEntries;
Int_t nTrksP=0,nTrksN=0;
Int_t trkNum[2];
Double_t tofmass[2];
Int_t okD0=0,okD0bar=0;
- AliITStrackV2 *postrack = 0;
- AliITStrackV2 *negtrack = 0;
+ AliESDtrack *postrack = 0;
+ AliESDtrack *negtrack = 0;
- // create the AliITSVertexerTracks object
+ // create the AliVertexerTracks object
// (it will be used only if fVertexOnTheFly=kTrue)
- AliITSVertexerTracks *vertexer1 = new AliITSVertexerTracks;
- vertexer1->SetMinTracks(2);
- vertexer1->SetDebug(0);
+ AliVertexerTracks *vertexer1 = new AliVertexerTracks;
+ if(fVertexOnTheFly) {
+ // open the mean vertex
+ TFile *invtx = new TFile("AliESDVertexMean.root");
+ AliESDVertex *initVertex = (AliESDVertex*)invtx->Get("vtxmean");
+ invtx->Close();
+ vertexer1->SetVtxStart(initVertex);
+ delete invtx;
+ }
Int_t skipped[2];
Bool_t goodVtx1;
-
- // define the cuts for vertexing
- Double_t vtxcuts[]={50., // max. allowed chi2
- 0.0, // min. allowed negative daughter's impact param
- 0.0, // min. allowed positive daughter's impact param
- 1.0, // max. allowed DCA between the daughter tracks
- -1.0, // min. allowed cosine of pointing angle
- 0.0, // min. radius of the fiducial volume
- 2.9};// max. radius of the fiducial volume
-
- // create the AliV0vertexer object
- AliV0vertexer *vertexer2 = new AliV0vertexer(vtxcuts);
-
// create tree for reconstructed D0s
AliD0toKpi *ioD0toKpi=0;
TTree *treeD0 = new TTree("TreeD0","Tree with D0 candidates");
// open file with tracks
TFile *esdFile = TFile::Open(esdName.Data());
+ AliESD* event = new AliESD;
+ TTree* tree = (TTree*) esdFile->Get("esdTree");
+ if(!tree) {
+ Error("FindCandidatesESD", "no ESD tree found");
+ return;
+ }
+ tree->SetBranchAddress("ESD",&event);
- TKey *key=0;
- TIter next(esdFile->GetListOfKeys());
// loop on events in file
- while ((key=(TKey*)next())!=0) {
- AliESD *event=(AliESD*)key->ReadObj();
- Int_t ev = (Int_t)event->GetEventNumber();
- if(ev<evFirst || ev>evLast) { delete event; continue; }
+ for(Int_t iEvent = evFirst; iEvent < tree->GetEntries(); iEvent++) {
+ if(iEvent > evLast) break;
+ tree->GetEvent(iEvent);
+ Int_t ev = (Int_t)event->GetEventNumberInFile(); // This is most likely NOT the event number you'd like to use. It has nothing to do with the 'real' event number.
printf("--- Finding D0 -> Kpi in event %d\n",ev);
+ // count the total number of events
+ nTotEv++;
- // retrieve primary vertex from file
- //sprintf(vtxName,"Vertex_%d",ev);
- //AliITSVertex *vertex1stored = (AliITSVertex*)trkFile->Get(vtxName);
- //vertex1stored->GetXYZ(fV1);
- //delete vertex1stored;
- event->GetVertex(fV1,covV1);
-
- trkEntries = event->GetNumberOfTracks();
+ trkEntries = (Int_t)event->GetNumberOfTracks();
printf(" Number of tracks: %d\n",trkEntries);
+ if(trkEntries<2) continue;
- // count the total number of events
- nTotEv++;
+ // retrieve primary vertex from file
+ if(!event->GetPrimaryVertex()) {
+ printf(" No vertex\n");
+ continue;
+ }
+ event->GetPrimaryVertex()->GetXYZ(fV1);
// call function which applies sigle-track selection and
// separetes positives and negatives
TObjArray trksN(trkEntries/2);
Int_t *trkEntryN = new Int_t[trkEntries];
TTree *trkTree = new TTree();
- if(fVertexOnTheFly) {
- SelectTracksESDvtx(*event,trkTree,trksP,trkEntryP,nTrksP,
- trksN,trkEntryN,nTrksN);
- } else {
- SelectTracksESD(*event,trksP,trkEntryP,nTrksP,
- trksN,trkEntryN,nTrksN);
- }
+ SelectTracks(event,trksP,trkEntryP,nTrksP,
+ trksN,trkEntryN,nTrksN);
+
+ printf(" pos. tracks: %d neg .tracks: %d\n",nTrksP,nTrksN);
- if(fDebug) printf(" pos. tracks: %d neg .tracks: %d\n",nTrksP,nTrksN);
nD0rec1ev = 0;
- // loop on positive tracks
+ // LOOP ON POSITIVE TRACKS
for(iTrkP=0; iTrkP<nTrksP; iTrkP++) {
- if((iTrkP%10==0) || fDebug) printf(" Processing positive track number %d of %d\n",iTrkP,nTrksP);
+ if(iTrkP%1==0) printf(" Processing positive track number %d of %d\n",iTrkP,nTrksP);
// get track from track array
- postrack = (AliITStrackV2*)trksP.UncheckedAt(iTrkP);
+ postrack = (AliESDtrack*)trksP.UncheckedAt(iTrkP);
trkNum[0] = trkEntryP[iTrkP];
- // loop on negative tracks
+ // LOOP ON NEGATIVE TRACKS
for(iTrkN=0; iTrkN<nTrksN; iTrkN++) {
// get track from tracks array
- negtrack = (AliITStrackV2*)trksN.UncheckedAt(iTrkN);
+ negtrack = (AliESDtrack*)trksN.UncheckedAt(iTrkN);
trkNum[1] = trkEntryN[iTrkN];
- AliITStrackV2 nt(*negtrack), pt(*postrack), *pnt=&nt, *ppt=&pt;
-
//
// ----------- DCA MINIMIZATION ------------------
//
- // find the DCA and propagate the tracks to the DCA
- dca = vertexer2->PropagateToDCA(pnt,ppt);
+ // find the DCA and propagate the tracks to the DCA
+ Double_t b=event->GetMagneticField();
+ AliESDtrack nt(*negtrack), pt(*postrack);
+ dca = nt.PropagateToDCA(&pt,b);
- // define the AliV0vertex object
- AliV0vertex *vertex2 = new AliV0vertex(*pnt,*ppt);
+ // define the AliESDv0 object
+ AliESDv0 vertex2(nt,trkNum[0],pt,trkNum[1]);
// get position of the secondary vertex
- vertex2->GetXYZ(v2[0],v2[1],v2[2]);
-
- delete vertex2;
-
- // momenta of the tracks at the vertex
- ptP = 1./TMath::Abs(ppt->Get1Pt());
- alphaP = ppt->GetAlpha();
- phiP = alphaP+TMath::ASin(ppt->GetSnp());
- mom[0] = ptP*TMath::Cos(phiP);
- mom[1] = ptP*TMath::Sin(phiP);
- mom[2] = ptP*ppt->GetTgl();
-
- ptN = 1./TMath::Abs(pnt->Get1Pt());
- alphaN = pnt->GetAlpha();
- phiN = alphaN+TMath::ASin(pnt->GetSnp());
- mom[3] = ptN*TMath::Cos(phiN);
- mom[4] = ptN*TMath::Sin(phiN);
- mom[5] = ptN*pnt->GetTgl();
-
+ vertex2.GetXYZ(v2[0],v2[1],v2[2]);
+ vertex2.GetPPxPyPz(mom[0],mom[1],mom[2]);
+ vertex2.GetNPxPyPz(mom[3],mom[4],mom[5]);
+ // impact parameters of the tracks w.r.t. the primary vertex
+ d0[0] = 10000.*pt.GetD(fV1[0],fV1[1],b);
+ d0[1] = -10000.*nt.GetD(fV1[0],fV1[1],b);
goodVtx1 = kTRUE;
+
// no vertexing if DeltaMass > fMassCut
if(fVertexOnTheFly) {
goodVtx1 = kFALSE;
if(SelectInvMass(mom)) {
- // primary vertex from *other* tracks in event
- vertexer1->SetVtxStart(fV1[0],fV1[1]);
+ // primary vertex from *other* tracks in the event
skipped[0] = trkEntryP[iTrkP];
skipped[1] = trkEntryN[iTrkN];
vertexer1->SetSkipTracks(2,skipped);
- AliITSVertex *vertex1onfly =
- (AliITSVertex*)vertexer1->VertexOnTheFly(*trkTree);
+ AliESDVertex *vertex1onfly =
+ (AliESDVertex*)vertexer1->FindPrimaryVertex(event);
if(vertex1onfly->GetNContributors()>0) goodVtx1 = kTRUE;
vertex1onfly->GetXYZ(fV1);
//vertex1onfly->PrintStatus();
delete vertex1onfly;
}
- }
-
- // impact parameters of the tracks w.r.t. the primary vertex
- d0[0] = 10000.*ppt->GetD(fV1[0],fV1[1]);
- d0[1] = -10000.*pnt->GetD(fV1[0],fV1[1]);
+ }
+
// create the object AliD0toKpi
AliD0toKpi theD0(ev,trkNum,fV1,v2,dca,mom,d0);
-
// select D0s
if(goodVtx1 && theD0.Select(fD0Cuts,okD0,okD0bar)) {
-
// get PID info from ESD
AliESDtrack *t0 = (AliESDtrack*)event->GetTrack(trkNum[0]);
Double_t esdpid0[5];
} // loop on negative tracks
postrack = 0;
} // loop on positive tracks
-
- trksP.Delete();
- trksN.Delete();
+
delete [] trkEntryP;
delete [] trkEntryN;
delete trkTree;
- delete event;
-
printf(" Number of D0 candidates: %d\n",nD0rec1ev);
} // loop on events in file
printf("\n+++\n+++ Total number of D0 candidates: %d\n+++\n",nD0rec);
delete vertexer1;
- delete vertexer2;
esdFile->Close();
} else {
printf(" Now adding information from simulation (PDG codes) ...\n");
TTree *treeD0sim = new TTree("TreeD0","Tree with D0 candidates");
- MakeTracksRefFileESD();
SimulationInfo(treeD0,treeD0sim);
delete treeD0;
TFile *outroot = new TFile(outName1.Data(),"recreate");
void AliD0toKpiAnalysis::PrintStatus() const {
// Print parameters being used
- printf(" fBz = %f T\n",fBz);
printf("Preselections:\n");
printf(" fPtCut = %f GeV\n",fPtCut);
printf(" fd0Cut = %f micron\n",fd0Cut);
return kFALSE;
}
//-----------------------------------------------------------------------------
-void AliD0toKpiAnalysis::SelectTracks(TTree &trkTree,
- TObjArray &trksP,Int_t *trkEntryP,Int_t &nTrksP,
- TObjArray &trksN,Int_t *trkEntryN,Int_t &nTrksN) const {
- // Create two TObjArrays with positive and negative tracks and
- // apply single-track preselection
-
- nTrksP=0,nTrksN=0;
-
- Int_t entr = (Int_t)trkTree.GetEntries();
-
- // trasfer tracks from tree to arrays
- for(Int_t i=0; i<entr; i++) {
-
- AliITStrackV2 *track = new AliITStrackV2;
- trkTree.SetBranchAddress("tracks",&track);
-
- trkTree.GetEvent(i);
-
- // single track selection
- if(!SingleTrkCuts(*track)) { delete track; continue; }
-
- if(track->Get1Pt()>0.) { // negative track
- trksN.AddLast(track);
- trkEntryN[nTrksN] = i;
- nTrksN++;
- } else { // positive track
- trksP.AddLast(track);
- trkEntryP[nTrksP] = i;
- nTrksP++;
- }
-
- }
-
- return;
-}
-//-----------------------------------------------------------------------------
-void AliD0toKpiAnalysis::SelectTracksESD(AliESD &event,
+void AliD0toKpiAnalysis::SelectTracks(AliESD *event,
TObjArray &trksP,Int_t *trkEntryP,Int_t &nTrksP,
TObjArray &trksN,Int_t *trkEntryN,Int_t &nTrksN) const {
// Create two TObjArrays with positive and negative tracks and
nTrksP=0,nTrksN=0;
- Int_t entr = event.GetNumberOfTracks();
+ Int_t entr = event->GetNumberOfTracks();
// transfer ITS tracks from ESD to arrays and to a tree
for(Int_t i=0; i<entr; i++) {
- //if(fDebug) printf(" SelectTracksESD: %d/%d\n",i,entr);
- AliESDtrack *esdtrack = (AliESDtrack*)event.GetTrack(i);
+ AliESDtrack *esdtrack = event->GetTrack(i);
UInt_t status = esdtrack->GetStatus();
- if(!(status&AliESDtrack::kITSrefit)) continue;
-
- AliITStrackV2 *itstrack = new AliITStrackV2(*esdtrack);
+ if(!(status&AliESDtrack::kITSin)) continue;
// single track selection
- if(!SingleTrkCuts(*itstrack))
- { delete itstrack; continue; }
+ if(!SingleTrkCuts(*esdtrack,event->GetMagneticField())) continue;
- if(itstrack->Get1Pt()>0.) { // negative track
- trksN.AddLast(itstrack);
+ if(esdtrack->GetSign()<0) { // negative track
+ trksN.AddLast(esdtrack);
trkEntryN[nTrksN] = i;
nTrksN++;
} else { // positive track
- trksP.AddLast(itstrack);
+ trksP.AddLast(esdtrack);
trkEntryP[nTrksP] = i;
nTrksP++;
}
return;
}
//-----------------------------------------------------------------------------
-void AliD0toKpiAnalysis::SelectTracksESDvtx(AliESD &event,TTree *trkTree,
- TObjArray &trksP,Int_t *trkEntryP,Int_t &nTrksP,
- TObjArray &trksN,Int_t *trkEntryN,Int_t &nTrksN) const {
- // Create two TObjArrays with positive and negative tracks and
- // apply single-track preselection
-
- nTrksP=0,nTrksN=0;
-
- Int_t entr = event.GetNumberOfTracks();
-
- AliITStrackV2 *itstrackfortree = 0;
- trkTree->Branch("tracks","AliITStrackV2",&itstrackfortree,entr,0);
-
-
- // transfer ITS tracks from ESD to arrays and to a tree
- for(Int_t i=0; i<entr; i++) {
-
- AliESDtrack *esdtrack = (AliESDtrack*)event.GetTrack(i);
- UInt_t status = esdtrack->GetStatus();
-
- if(!(status&AliESDtrack::kITSrefit)) continue;
-
- AliITStrackV2 *itstrack = new AliITStrackV2(*esdtrack);
-
- // store track in the tree to be used for primary vertex finding
- itstrackfortree = new AliITStrackV2(*esdtrack);
- trkTree->Fill();
- itstrackfortree = 0;
-
- // single track selection
- if(!SingleTrkCuts(*itstrack))
- { delete itstrack; continue; }
-
- if(itstrack->Get1Pt()>0.) { // negative track
- trksN.AddLast(itstrack);
- trkEntryN[nTrksN] = i;
- nTrksN++;
- } else { // positive track
- trksP.AddLast(itstrack);
- trkEntryP[nTrksP] = i;
- nTrksP++;
- }
-
- } // loop on esd tracks
-
- delete itstrackfortree;
-
- return;
-}
-//-----------------------------------------------------------------------------
void AliD0toKpiAnalysis::SetD0Cuts(Double_t cut0,Double_t cut1,
Double_t cut2,Double_t cut3,Double_t cut4,
Double_t cut5,Double_t cut6,
return;
}
//-----------------------------------------------------------------------------
-Bool_t AliD0toKpiAnalysis::SingleTrkCuts(const AliITStrackV2& trk) const {
+Bool_t
+AliD0toKpiAnalysis::SingleTrkCuts(const AliESDtrack& trk, Double_t b) const {
// Check if track passes some kinematical cuts
+ // Magnetic field "b" (kG)
- if(TMath::Abs(1./trk.Get1Pt()) < fPtCut)
+ if(TMath::Abs(1./trk.GetParameter()[4]) < fPtCut)
return kFALSE;
- if(TMath::Abs(10000.*trk.GetD(fV1[0],fV1[1])) < fd0Cut)
+ if(TMath::Abs(10000.*trk.GetD(fV1[0],fV1[1],b)) < fd0Cut)
return kFALSE;
return kTRUE;
}
//----------------------------------------------------------------------------
-void AliD0toKpiAnalysis::MakeTracksRefFile(Int_t evFirst,Int_t evLast)
- const {
- // Create a file with simulation info for the reconstructed tracks
-
- TFile *out = TFile::Open("ITStracksRefFile.root","recreate");
- TFile *trk = TFile::Open("AliITStracksV2.root");
- AliRunLoader *rl = AliRunLoader::Open("galice.root");
-
- // load kinematics
- rl->LoadKinematics();
-
- Int_t label;
- TParticle *part;
- TParticle *mumpart;
- REFTRACK reftrk;
-
- for(Int_t ev=evFirst; ev<=evLast; ev++){
- rl->GetEvent(ev);
- AliStack *stack = rl->Stack();
-
- trk->cd();
-
- // Tree with tracks
- char tname[100];
- sprintf(tname,"TreeT_ITS_%d",ev);
-
- TTree *tracktree=(TTree*)trk->Get(tname);
- if(!tracktree) continue;
- AliITStrackV2 *track = new AliITStrackV2;
- tracktree->SetBranchAddress("tracks",&track);
- Int_t nentr=(Int_t)tracktree->GetEntries();
-
- // Tree for true track parameters
- char ttname[100];
- sprintf(ttname,"Tree_Ref_%d",ev);
- TTree *reftree = new TTree(ttname,"Tree with true track params");
- reftree->Branch("rectracks",&reftrk,"lab/I:pdg:mumlab:mumpdg:Vx/F:Vy:Vz:Px:Py:Pz");
-
- for(Int_t i=0; i<nentr; i++) {
- tracktree->GetEvent(i);
- label = TMath::Abs(track->GetLabel());
-
- part = (TParticle*)stack->Particle(label);
- reftrk.lab = label;
- reftrk.pdg = part->GetPdgCode();
- reftrk.mumlab = part->GetFirstMother();
- if(part->GetFirstMother()>=0) {
- mumpart = (TParticle*)stack->Particle(part->GetFirstMother());
- reftrk.mumpdg = mumpart->GetPdgCode();
- } else {
- reftrk.mumpdg=-1;
- }
- reftrk.Vx = part->Vx();
- reftrk.Vy = part->Vy();
- reftrk.Vz = part->Vz();
- reftrk.Px = part->Px();
- reftrk.Py = part->Py();
- reftrk.Pz = part->Pz();
-
- reftree->Fill();
- } // loop on tracks
-
- out->cd();
- reftree->Write();
-
- delete track;
- delete reftree;
- delete tracktree;
- delete stack;
- } // loop on events
-
- trk->Close();
- out->Close();
-
- return;
-}
-//----------------------------------------------------------------------------
-void AliD0toKpiAnalysis::MakeTracksRefFileESD() const {
+void AliD0toKpiAnalysis::MakeTracksRefFile(AliRun *gAlice,
+ Int_t evFirst,Int_t evLast) const {
// Create a file with simulation info for the reconstructed tracks
- TFile *outFile = TFile::Open("ITStracksRefFile.root","recreate");
+ TFile *outFile = TFile::Open("D0TracksRefFile.root","recreate");
TFile *esdFile = TFile::Open("AliESDs.root");
- AliRunLoader *rl = AliRunLoader::Open("galice.root");
- // load kinematics
- rl->LoadKinematics();
+ AliMC *mc = gAlice->GetMCApp();
Int_t label;
TParticle *part;
TParticle *mumpart;
REFTRACK reftrk;
- TKey *key=0;
- TIter next(esdFile->GetListOfKeys());
+ AliESD* event = new AliESD;
+ TTree* tree = (TTree*) esdFile->Get("esdTree");
+ tree->SetBranchAddress("ESD",&event);
// loop on events in file
- while ((key=(TKey*)next())!=0) {
- AliESD *event=(AliESD*)key->ReadObj();
- Int_t ev = (Int_t)event->GetEventNumber();
+ for(Int_t iEvent=evFirst; iEvent<tree->GetEntries(); iEvent++) {
+ if(iEvent>evLast) break;
+ tree->GetEvent(iEvent);
+ Int_t ev = (Int_t)event->GetEventNumberInFile(); // This is most likely NOT the event number you'd like to use. It has nothing to do with the 'real' event number.
- rl->GetEvent(ev);
- AliStack *stack = rl->Stack();
+ gAlice->GetEvent(ev);
Int_t nentr=(Int_t)event->GetNumberOfTracks();
AliESDtrack *esdtrack = (AliESDtrack*)event->GetTrack(i);
label = TMath::Abs(esdtrack->GetLabel());
- part = (TParticle*)stack->Particle(label);
+ part = (TParticle*)mc->Particle(label);
reftrk.lab = label;
reftrk.pdg = part->GetPdgCode();
reftrk.mumlab = part->GetFirstMother();
if(part->GetFirstMother()>=0) {
- mumpart = (TParticle*)stack->Particle(part->GetFirstMother());
+ mumpart = (TParticle*)gAlice->GetMCApp()->Particle(part->GetFirstMother());
reftrk.mumpdg = mumpart->GetPdgCode();
+ reftrk.mumprongs = mumpart->GetNDaughters();
} else {
reftrk.mumpdg=-1;
+ reftrk.mumprongs=-1;
}
reftrk.Vx = part->Vx();
reftrk.Vy = part->Vy();
reftree->Write();
delete reftree;
- delete event;
- delete stack;
} // loop on events
esdFile->Close();
void AliD0toKpiAnalysis::SimulationInfo(TTree *treeD0in,TTree *treeD0out) const {
// add pdg codes to candidate decay tracks (for sim)
- TString refFileName("ITStracksRefFile.root");
+ TString refFileName("D0TracksRefFile.root");
if(fSim && gSystem->AccessPathName(refFileName.Data(),kFileExists)) {
printf("AliD0toKpiAnalysis::SimulationInfo: no reference file found!\n");
return;
theD0->SetPdgCodes(pdg);
theD0->SetMumPdgCodes(mumpdg);
- if(TMath::Abs(mumpdg[0])==421 && TMath::Abs(mumpdg[1])==421
- && mumlab[0]==mumlab[1]) theD0->SetSignal();
+ if(TMath::Abs(mumpdg[0])==421 &&
+ TMath::Abs(mumpdg[1])==421 &&
+ mumlab[0]==mumlab[1] &&
+ reftrk.mumprongs==2 &&
+ ((TMath::Abs(pdg[0])==211 && TMath::Abs(pdg[1])==321) ||
+ (TMath::Abs(pdg[0])==321 && TMath::Abs(pdg[1])==211))
+ ) theD0->SetSignal();
if(!fOnlySignal || theD0->IsSignal()) treeD0out->Fill();
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff