X-Git-Url: http://git.uio.no/git/?p=u%2Fmrichter%2FAliRoot.git;a=blobdiff_plain;f=ANALYSIS%2FAliD0toKpiAnalysis.cxx;h=e9294d2dd6c455b52076f888f1ae50f32a04c159;hp=00f6e914d16f4b7ffbf7783aae2af22c68283cb3;hb=aef31c8d1dd03c07019169d65ccedcf7c0e3a7be;hpb=c7bafca9a7bac3b63264cb649d33a8172a486e94 diff --git a/ANALYSIS/AliD0toKpiAnalysis.cxx b/ANALYSIS/AliD0toKpiAnalysis.cxx index 00f6e914d16..e9294d2dd6c 100644 --- a/ANALYSIS/AliD0toKpiAnalysis.cxx +++ b/ANALYSIS/AliD0toKpiAnalysis.cxx @@ -18,7 +18,7 @@ // Note: the two decay tracks are labelled: 0 (positive track) // 1 (negative track) // An example of usage can be found in the macro AliD0toKpiTest.C -// Origin: A. Dainese andrea.dainese@pd.infn.it +// Origin: A. Dainese andrea.dainese@lnl.infn.it //---------------------------------------------------------------------------- #include #include @@ -27,9 +27,10 @@ #include #include #include "AliESD.h" -#include "AliStack.h" +#include "AliMC.h" +#include "AliRun.h" #include "AliRunLoader.h" -#include "AliITSVertexerTracks.h" +#include "AliVertexerTracks.h" #include "AliESDVertex.h" #include "AliESDv0.h" #include "AliD0toKpi.h" @@ -41,6 +42,7 @@ typedef struct { Int_t pdg; Int_t mumlab; Int_t mumpdg; + Int_t mumprongs; Float_t Vx,Vy,Vz; Float_t Px,Py,Pz; } REFTRACK; @@ -51,7 +53,6 @@ ClassImp(AliD0toKpiAnalysis) AliD0toKpiAnalysis::AliD0toKpiAnalysis() { // Default constructor - fBz=-9999; SetPtCut(); Setd0Cut(); SetMassCut(); @@ -133,261 +134,24 @@ Double_t AliD0toKpiAnalysis::CalculateTOFmass(Double_t mom,Double_t length, return mom*a; } -/* //---------------------------------------------------------------------------- void AliD0toKpiAnalysis::FindCandidates(Int_t evFirst,Int_t evLast, const Char_t *outName) { // Find D0 candidates and calculate parameters - if(fBz<-9000.) { - printf("AliD0toKpiAnalysis::FindCandidates(): Set B!\n"); - return; - } - - 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); - 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); - AliESDVertex *vertex1stored = (AliESDVertex*)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; iTrkPPropagateToDCA(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); - AliESDVertex *vertex1onfly = - (AliESDVertex*)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; - } - - 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"; 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]; @@ -396,27 +160,20 @@ void AliD0toKpiAnalysis::FindCandidatesESD(Int_t evFirst,Int_t evLast, 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); + 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"); @@ -426,31 +183,31 @@ void AliD0toKpiAnalysis::FindCandidatesESD(Int_t evFirst,Int_t evLast, TFile *esdFile = TFile::Open(esdName.Data()); AliESD* event = new AliESD; TTree* tree = (TTree*) esdFile->Get("esdTree"); - if (!tree) { + if(!tree) { Error("FindCandidatesESD", "no ESD tree found"); return; } - tree->SetBranchAddress("ESD", &event); + tree->SetBranchAddress("ESD",&event); // loop on events in file - for (Int_t iEvent = evFirst; iEvent < tree->GetEntries(); iEvent++) { - if (iEvent > evLast) break; + for(Int_t iEvent = evFirst; iEvent < tree->GetEntries(); iEvent++) { + if(iEvent > evLast) break; tree->GetEvent(iEvent); - Int_t ev = (Int_t)event->GetEventNumber(); + 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); - //AliESDVertex *vertex1stored = (AliESDVertex*)trkFile->Get(vtxName); - //vertex1stored->GetXYZ(fV1); - //delete vertex1stored; - event->GetVertex()->GetXYZ(fV1); - - 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 @@ -459,34 +216,29 @@ void AliD0toKpiAnalysis::FindCandidatesESD(Int_t evFirst,Int_t evLast, 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); - AliDebugClass(1,Form(" pos. tracks: %d neg .tracks: %d",nTrksP,nTrksN)); nD0rec1ev = 0; - // loop on positive tracks + // LOOP ON POSITIVE TRACKS for(iTrkP=0; iTrkPGetExternalParameters(x,par); - //ptP = 1./TMath::Abs(par[4]); - //alphaP = postrack->GetAlpha(); - //phiP = alphaP+TMath::ASin(par[2]); - postrack->GetPxPyPz(); - mom[0] = ptP*TMath::Cos(phiP); - mom[1] = ptP*TMath::Sin(phiP); - mom[2] = ptP*par[3]; - - 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]); + // primary vertex from *other* tracks in the event skipped[0] = trkEntryP[iTrkP]; skipped[1] = trkEntryN[iTrkN]; vertexer1->SetSkipTracks(2,skipped); AliESDVertex *vertex1onfly = - (AliESDVertex*)vertexer1->VertexOnTheFly(*trkTree); + (AliESDVertex*)vertexer1->FindPrimaryVertex(event); if(vertex1onfly->GetNContributors()>0) goodVtx1 = kTRUE; vertex1onfly->GetXYZ(fV1); //vertex1onfly->PrintStatus(); delete vertex1onfly; } - } + } + // 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]; @@ -587,14 +318,11 @@ void AliD0toKpiAnalysis::FindCandidatesESD(Int_t evFirst,Int_t evLast, } // 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 @@ -603,7 +331,6 @@ void AliD0toKpiAnalysis::FindCandidatesESD(Int_t evFirst,Int_t evLast, printf("\n+++\n+++ Total number of D0 candidates: %d\n+++\n",nD0rec); delete vertexer1; - //delete vertexer2; esdFile->Close(); @@ -621,7 +348,6 @@ void AliD0toKpiAnalysis::FindCandidatesESD(Int_t evFirst,Int_t evLast, } 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"); @@ -642,7 +368,6 @@ void AliD0toKpiAnalysis::FindCandidatesESD(Int_t evFirst,Int_t evLast, 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); @@ -699,46 +424,8 @@ Bool_t AliD0toKpiAnalysis::SelectInvMass(const Double_t p[6]) const { if(TMath::Abs(minvD0bar-mD0) < fMassCut) return kTRUE; 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; iGet1Pt()>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 @@ -746,18 +433,18 @@ void AliD0toKpiAnalysis::SelectTracksESD(AliESD &event, 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; iGetTrack(i); UInt_t status = esdtrack->GetStatus(); - if(!(status&AliESDtrack::kITSrefit)) continue; + if(!(status&AliESDtrack::kITSin)) continue; // single track selection - if(!SingleTrkCuts(*esdtrack,event.GetMagneticField())) continue; + if(!SingleTrkCuts(*esdtrack,event->GetMagneticField())) continue; if(esdtrack->GetSign()<0) { // negative track trksN.AddLast(esdtrack); @@ -774,54 +461,6 @@ void AliD0toKpiAnalysis::SelectTracksESD(AliESD &event, 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(); - - AliESDtrack *esdtrackfortree = 0; - trkTree->Branch("tracks","AliESDtrack",&esdtrackfortree,entr,0); - - - // transfer the tracks from ESD to arrays and to a tree - for(Int_t i=0; iGetStatus(); - - if(!(status&AliESDtrack::kITSrefit)) continue; - - // store track in the tree to be used for primary vertex finding - esdtrackfortree = new AliESDtrack(*esdtrack); - trkTree->Fill(); - //itstrackfortree = 0; - delete esdtrackfortree; - - // single track selection - if(!SingleTrkCuts(*esdtrack,event.GetMagneticField())) continue; - - if(esdtrack->GetSign()<0) { // negative track - trksN.AddLast(esdtrack); - trkEntryN[nTrksN] = i; - nTrksN++; - } else { // positive track - trksP.AddLast(esdtrack); - 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, @@ -860,110 +499,31 @@ AliD0toKpiAnalysis::SingleTrkCuts(const AliESDtrack& trk, Double_t b) const { 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; iGetEvent(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; iEventGetEntries(); 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(); @@ -977,15 +537,17 @@ void AliD0toKpiAnalysis::MakeTracksRefFileESD() const { 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(); @@ -1002,8 +564,6 @@ void AliD0toKpiAnalysis::MakeTracksRefFileESD() const { reftree->Write(); delete reftree; - delete event; - delete stack; } // loop on events esdFile->Close(); @@ -1015,7 +575,7 @@ void AliD0toKpiAnalysis::MakeTracksRefFileESD() const { 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; @@ -1071,8 +631,13 @@ void AliD0toKpiAnalysis::SimulationInfo(TTree *treeD0in,TTree *treeD0out) const 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();