1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
16 //----------------------------------------------------------------------------
17 // Implementation of the D0toKpi reconstruction and analysis class
18 // Note: the two decay tracks are labelled: 0 (positive track)
20 // An example of usage can be found in the macro AliD0toKpiTest.C
21 // Origin: A. Dainese andrea.dainese@lnl.infn.it
22 //----------------------------------------------------------------------------
28 #include <TParticle.h>
29 #include "AliESDEvent.h"
32 #include "AliRunLoader.h"
33 #include "AliVertexerTracks.h"
34 #include "AliESDVertex.h"
36 #include "AliD0toKpi.h"
37 #include "AliD0toKpiAnalysis.h"
50 ClassImp(AliD0toKpiAnalysis)
52 //----------------------------------------------------------------------------
53 AliD0toKpiAnalysis::AliD0toKpiAnalysis():
54 fVertexOnTheFly(kFALSE),
57 fPID("TOFparam_PbPb"),
62 // Default constructor
67 //----------------------------------------------------------------------------
68 AliD0toKpiAnalysis::~AliD0toKpiAnalysis() {}
69 //----------------------------------------------------------------------------
70 void AliD0toKpiAnalysis::ApplySelection(const Char_t *inName,const Char_t *outName) const {
71 // select candidates that pass fD0Cuts and write them to a new file
73 TFile *inFile = TFile::Open(inName);
75 TTree *treeD0in=(TTree*)inFile->Get("TreeD0");
76 AliD0toKpiAnalysis *inAnalysis = (AliD0toKpiAnalysis*)inFile->Get("D0toKpiAnalysis");
77 printf("+++\n+++ I N P U T S T A T U S:\n+++\n");
78 inAnalysis->PrintStatus();
82 treeD0in->SetBranchAddress("D0toKpi",&d);
83 Int_t entries = (Int_t)treeD0in->GetEntries();
85 printf("+++\n+++ Number of D0 in input tree: %d\n+++\n",entries);
87 TTree *treeD0out = new TTree("TreeD0","Tree with selected D0 candidates");
88 treeD0out->Branch("D0toKpi","AliD0toKpi",&d,200000,0);
91 Int_t okD0=0,okD0bar=0;
94 for(Int_t i=0; i<entries; i++) {
95 // get event from tree
96 treeD0in->GetEvent(i);
98 if(fSim && fOnlySignal && !d->IsSignal()) continue;
100 // check if candidate passes selection (as D0 or D0bar)
101 if(d->Select(fD0Cuts,okD0,okD0bar)) {
108 AliD0toKpiAnalysis *outAnalysis = (AliD0toKpiAnalysis*)inAnalysis->Clone("D0toKpiAnalysis");
109 outAnalysis->SetD0Cuts(fD0Cuts);
110 printf("------------------------------------------\n");
111 printf("+++\n+++ O U T P U T S T A T U S:\n+++\n");
112 outAnalysis->PrintStatus();
114 printf("+++\n+++ Number of D0 in output tree: %d\n+++\n",nSel);
116 TFile* outFile = new TFile(outName,"recreate");
118 outAnalysis->Write();
123 //----------------------------------------------------------------------------
124 Double_t AliD0toKpiAnalysis::CalculateTOFmass(Double_t mom,Double_t length,
125 Double_t time) const {
126 // calculated the mass from momentum, track length from vertex to TOF
127 // and time measured by the TOF
128 if(length==0.) return -1000.;
129 Double_t a = time*time/length/length;
131 a = TMath::Sqrt(a-1.);
133 a = -TMath::Sqrt(1.-a);
138 //----------------------------------------------------------------------------
139 void AliD0toKpiAnalysis::FindCandidates(Int_t evFirst,Int_t evLast,
140 const Char_t *outName) {
141 // Find D0 candidates and calculate parameters
144 TString esdName="AliESDs.root";
145 if(gSystem->AccessPathName(esdName.Data(),kFileExists)) {
146 printf("AliD0toKpiAnalysis::FindCandidatesESD(): No ESDs file found!\n");
150 TString outName1=outName;
151 TString outName2="nTotEvents.dat";
153 Int_t nTotEv=0,nD0rec=0,nD0rec1ev=0;
155 Double_t v2[3],mom[6],d0[2];
156 Int_t iTrkP,iTrkN,trkEntries;
157 Int_t nTrksP=0,nTrksN=0;
160 Int_t okD0=0,okD0bar=0;
161 AliESDtrack *postrack = 0;
162 AliESDtrack *negtrack = 0;
164 // create the AliVertexerTracks object
165 // (it will be used only if fVertexOnTheFly=kTrue)
166 AliVertexerTracks *vertexer1 = new AliVertexerTracks();
167 if(fVertexOnTheFly) {
168 // open the mean vertex
169 TFile *invtx = new TFile("AliESDVertexMean.root");
170 AliESDVertex *initVertex = (AliESDVertex*)invtx->Get("vtxmean");
172 vertexer1->SetVtxStart(initVertex);
178 // create tree for reconstructed D0s
179 AliD0toKpi *ioD0toKpi=0;
180 TTree *treeD0 = new TTree("TreeD0","Tree with D0 candidates");
181 treeD0->Branch("D0toKpi","AliD0toKpi",&ioD0toKpi,200000,0);
183 // open file with tracks
184 TFile *esdFile = TFile::Open(esdName.Data());
185 AliESDEvent* event = new AliESDEvent();
186 TTree* tree = (TTree*) esdFile->Get("esdTree");
188 Error("FindCandidatesESD", "no ESD tree found");
191 event->ReadFromTree(tree);
193 // loop on events in file
194 for(Int_t iEvent = evFirst; iEvent < tree->GetEntries(); iEvent++) {
195 if(iEvent > evLast) break;
196 tree->GetEvent(iEvent);
197 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.
198 printf("--- Finding D0 -> Kpi in event %d\n",ev);
199 // count the total number of events
202 trkEntries = (Int_t)event->GetNumberOfTracks();
203 printf(" Number of tracks: %d\n",trkEntries);
204 if(trkEntries<2) continue;
206 // retrieve primary vertex from file
207 if(!event->GetPrimaryVertex()) {
208 printf(" No vertex\n");
211 event->GetPrimaryVertex()->GetXYZ(fV1);
213 // call function which applies sigle-track selection and
214 // separetes positives and negatives
215 TObjArray trksP(trkEntries/2);
216 Int_t *trkEntryP = new Int_t[trkEntries];
217 TObjArray trksN(trkEntries/2);
218 Int_t *trkEntryN = new Int_t[trkEntries];
219 TTree *trkTree = new TTree();
220 SelectTracks(event,trksP,trkEntryP,nTrksP,
221 trksN,trkEntryN,nTrksN);
223 printf(" pos. tracks: %d neg .tracks: %d\n",nTrksP,nTrksN);
228 // LOOP ON POSITIVE TRACKS
229 for(iTrkP=0; iTrkP<nTrksP; iTrkP++) {
230 if(iTrkP%1==0) printf(" Processing positive track number %d of %d\n",iTrkP,nTrksP);
232 // get track from track array
233 postrack = (AliESDtrack*)trksP.UncheckedAt(iTrkP);
234 trkNum[0] = trkEntryP[iTrkP];
236 // LOOP ON NEGATIVE TRACKS
237 for(iTrkN=0; iTrkN<nTrksN; iTrkN++) {
239 // get track from tracks array
240 negtrack = (AliESDtrack*)trksN.UncheckedAt(iTrkN);
241 trkNum[1] = trkEntryN[iTrkN];
244 // ----------- DCA MINIMIZATION ------------------
246 // find the DCA and propagate the tracks to the DCA
247 Double_t b=event->GetMagneticField();
248 AliESDtrack nt(*negtrack), pt(*postrack);
249 dca = nt.PropagateToDCA(&pt,b);
251 // define the AliESDv0 object
252 AliESDv0 vertex2(nt,trkNum[0],pt,trkNum[1]);
254 // get position of the secondary vertex
255 vertex2.GetXYZ(v2[0],v2[1],v2[2]);
256 vertex2.GetPPxPyPz(mom[0],mom[1],mom[2]);
257 vertex2.GetNPxPyPz(mom[3],mom[4],mom[5]);
258 // impact parameters of the tracks w.r.t. the primary vertex
259 d0[0] = 10000.*pt.GetD(fV1[0],fV1[1],b);
260 d0[1] = -10000.*nt.GetD(fV1[0],fV1[1],b);
263 // no vertexing if DeltaMass > fMassCut
264 if(fVertexOnTheFly) {
266 if(SelectInvMass(mom)) {
267 // primary vertex from *other* tracks in the event
268 vertexer1->SetFieldkG(event->GetMagneticField());
269 skipped[0] = trkEntryP[iTrkP];
270 skipped[1] = trkEntryN[iTrkN];
271 vertexer1->SetSkipTracks(2,skipped);
272 AliESDVertex *vertex1onfly =
273 (AliESDVertex*)vertexer1->FindPrimaryVertex(event);
274 if(vertex1onfly->GetNContributors()>0) goodVtx1 = kTRUE;
275 vertex1onfly->GetXYZ(fV1);
276 //vertex1onfly->PrintStatus();
282 // create the object AliD0toKpi
283 AliD0toKpi theD0(ev,trkNum,fV1,v2,dca,mom,d0);
285 if(goodVtx1 && theD0.Select(fD0Cuts,okD0,okD0bar)) {
286 // get PID info from ESD
287 AliESDtrack *t0 = (AliESDtrack*)event->GetTrack(trkNum[0]);
289 t0->GetESDpid(esdpid0);
290 if(t0->GetStatus()&AliESDtrack::kTOFpid) {
291 tofmass[0] = CalculateTOFmass(t0->GetP(),
292 t0->GetIntegratedLength(),
297 AliESDtrack *t1 = (AliESDtrack*)event->GetTrack(trkNum[1]);
299 t1->GetESDpid(esdpid1);
300 if(t1->GetStatus()&AliESDtrack::kTOFpid) {
301 tofmass[1] = CalculateTOFmass(t1->GetP(),
302 t1->GetIntegratedLength(),
308 theD0.SetPIDresponse(esdpid0,esdpid1);
309 theD0.SetTOFmasses(tofmass);
315 nD0rec++; nD0rec1ev++;
320 } // loop on negative tracks
322 } // loop on positive tracks
328 printf(" Number of D0 candidates: %d\n",nD0rec1ev);
329 } // loop on events in file
332 printf("\n+++\n+++ Total number of events: %d\n+++\n",nTotEv);
333 printf("\n+++\n+++ Total number of D0 candidates: %d\n+++\n",nD0rec);
339 // create a copy of this class to be written to output file
340 AliD0toKpiAnalysis *copy = (AliD0toKpiAnalysis*)this->Clone("D0toKpiAnalysis");
342 // add PDG codes to decay tracks in found candidates (in simulation mode)
343 // and store tree in the output file
345 TFile *outroot = new TFile(outName1.Data(),"recreate");
351 printf(" Now adding information from simulation (PDG codes) ...\n");
352 TTree *treeD0sim = new TTree("TreeD0","Tree with D0 candidates");
353 SimulationInfo(treeD0,treeD0sim);
355 TFile *outroot = new TFile(outName1.Data(),"recreate");
362 // write to a file the total number of events
363 FILE *outdat = fopen(outName2.Data(),"w");
364 fprintf(outdat,"%d\n",nTotEv);
369 //-----------------------------------------------------------------------------
370 void AliD0toKpiAnalysis::PrintStatus() const {
371 // Print parameters being used
373 printf("Preselections:\n");
374 printf(" fPtCut = %f GeV\n",fPtCut);
375 printf(" fd0Cut = %f micron\n",fd0Cut);
376 printf(" fMassCut = %f GeV\n",fMassCut);
377 if(fVertexOnTheFly) printf("Primary vertex on the fly\n");
379 printf("Simulation mode\n");
380 if(fOnlySignal) printf(" Only signal goes to file\n");
382 printf("Cuts on candidates:\n");
383 printf(" |M-MD0| [GeV] < %f\n",fD0Cuts[0]);
384 printf(" dca [micron] < %f\n",fD0Cuts[1]);
385 printf(" cosThetaStar < %f\n",fD0Cuts[2]);
386 printf(" pTK [GeV] > %f\n",fD0Cuts[3]);
387 printf(" pTpi [GeV] > %f\n",fD0Cuts[4]);
388 printf(" |d0K| [micron] < %f\n",fD0Cuts[5]);
389 printf(" |d0pi| [micron] < %f\n",fD0Cuts[6]);
390 printf(" d0d0 [micron^2] < %f\n",fD0Cuts[7]);
391 printf(" cosThetaPoint > %f\n",fD0Cuts[8]);
395 //-----------------------------------------------------------------------------
396 Bool_t AliD0toKpiAnalysis::SelectInvMass(const Double_t p[6]) const {
397 // Apply preselection in the invariant mass of the pair
399 Double_t mD0 = 1.8645;
400 Double_t mPi = 0.13957;
401 Double_t mKa = 0.49368;
404 Double_t mom2[2],momTot2;
406 mom2[0] = p[0]*p[0] + p[1]*p[1] + p[2]*p[2];
407 mom2[1] = p[3]*p[3] + p[4]*p[4] + p[5]*p[5];
409 momTot2 = (p[0]+p[3])*(p[0]+p[3])+
410 (p[1]+p[4])*(p[1]+p[4])+
411 (p[2]+p[5])*(p[2]+p[5]);
414 energy[1] = TMath::Sqrt(mKa*mKa+mom2[1]);
415 energy[0] = TMath::Sqrt(mPi*mPi+mom2[0]);
417 Double_t minvD0 = TMath::Sqrt((energy[0]+energy[1])*(energy[0]+energy[1])-momTot2);
420 energy[0] = TMath::Sqrt(mKa*mKa+mom2[0]);
421 energy[1] = TMath::Sqrt(mPi*mPi+mom2[1]);
423 Double_t minvD0bar = TMath::Sqrt((energy[0]+energy[1])*(energy[0]+energy[1])-momTot2);
425 if(TMath::Abs(minvD0-mD0) < fMassCut) return kTRUE;
426 if(TMath::Abs(minvD0bar-mD0) < fMassCut) return kTRUE;
429 //-----------------------------------------------------------------------------
430 void AliD0toKpiAnalysis::SelectTracks(AliESDEvent *event,
431 TObjArray &trksP,Int_t *trkEntryP,Int_t &nTrksP,
432 TObjArray &trksN,Int_t *trkEntryN,Int_t &nTrksN) const {
433 // Create two TObjArrays with positive and negative tracks and
434 // apply single-track preselection
438 Int_t entr = event->GetNumberOfTracks();
440 // transfer ITS tracks from ESD to arrays and to a tree
441 for(Int_t i=0; i<entr; i++) {
443 AliESDtrack *esdtrack = event->GetTrack(i);
444 UInt_t status = esdtrack->GetStatus();
446 if(!(status&AliESDtrack::kITSin)) continue;
448 // single track selection
449 if(!SingleTrkCuts(*esdtrack,event->GetMagneticField())) continue;
451 if(esdtrack->GetSign()<0) { // negative track
452 trksN.AddLast(esdtrack);
453 trkEntryN[nTrksN] = i;
455 } else { // positive track
456 trksP.AddLast(esdtrack);
457 trkEntryP[nTrksP] = i;
461 } // loop on ESD tracks
465 //-----------------------------------------------------------------------------
466 void AliD0toKpiAnalysis::SetD0Cuts(Double_t cut0,Double_t cut1,
467 Double_t cut2,Double_t cut3,Double_t cut4,
468 Double_t cut5,Double_t cut6,
469 Double_t cut7,Double_t cut8) {
470 // Set the cuts for D0 selection
483 //-----------------------------------------------------------------------------
484 void AliD0toKpiAnalysis::SetD0Cuts(const Double_t cuts[9]) {
485 // Set the cuts for D0 selection
487 for(Int_t i=0; i<9; i++) fD0Cuts[i] = cuts[i];
491 //-----------------------------------------------------------------------------
493 AliD0toKpiAnalysis::SingleTrkCuts(const AliESDtrack& trk, Double_t b) const {
494 // Check if track passes some kinematical cuts
495 // Magnetic field "b" (kG)
497 if(TMath::Abs(1./trk.GetParameter()[4]) < fPtCut)
499 if(TMath::Abs(10000.*trk.GetD(fV1[0],fV1[1],b)) < fd0Cut)
504 //----------------------------------------------------------------------------
505 void AliD0toKpiAnalysis::MakeTracksRefFile(AliRun *mygAlice,
506 Int_t evFirst,Int_t evLast) const {
507 // Create a file with simulation info for the reconstructed tracks
509 TFile *outFile = TFile::Open("D0TracksRefFile.root","recreate");
510 TFile *esdFile = TFile::Open("AliESDs.root");
512 AliMC *mc = mygAlice->GetMCApp();
519 AliESDEvent* event = new AliESDEvent();
520 TTree* tree = (TTree*) esdFile->Get("esdTree");
521 event->ReadFromTree(tree);
522 // loop on events in file
523 for(Int_t iEvent=evFirst; iEvent<tree->GetEntries(); iEvent++) {
524 if(iEvent>evLast) break;
525 tree->GetEvent(iEvent);
526 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.
528 mygAlice->GetEvent(ev);
530 Int_t nentr=(Int_t)event->GetNumberOfTracks();
532 // Tree for true track parameters
534 sprintf(ttname,"Tree_Ref_%d",ev);
535 TTree *reftree = new TTree(ttname,"Tree with true track params");
536 reftree->Branch("rectracks",&reftrk,"lab/I:pdg:mumlab:mumpdg:Vx/F:Vy:Vz:Px:Py:Pz");
538 for(Int_t i=0; i<nentr; i++) {
539 AliESDtrack *esdtrack = (AliESDtrack*)event->GetTrack(i);
540 label = TMath::Abs(esdtrack->GetLabel());
542 part = (TParticle*)mc->Particle(label);
544 reftrk.pdg = part->GetPdgCode();
545 reftrk.mumlab = part->GetFirstMother();
546 if(part->GetFirstMother()>=0) {
547 mumpart = (TParticle*)mygAlice->GetMCApp()->Particle(part->GetFirstMother());
548 reftrk.mumpdg = mumpart->GetPdgCode();
549 reftrk.mumprongs = mumpart->GetNDaughters();
554 reftrk.Vx = part->Vx();
555 reftrk.Vy = part->Vy();
556 reftrk.Vz = part->Vz();
557 reftrk.Px = part->Px();
558 reftrk.Py = part->Py();
559 reftrk.Pz = part->Pz();
576 //-----------------------------------------------------------------------------
577 void AliD0toKpiAnalysis::SimulationInfo(TTree *treeD0in,TTree *treeD0out) const {
578 // add pdg codes to candidate decay tracks (for sim)
580 TString refFileName("D0TracksRefFile.root");
581 if(fSim && gSystem->AccessPathName(refFileName.Data(),kFileExists)) {
582 printf("AliD0toKpiAnalysis::SimulationInfo: no reference file found!\n");
585 TFile *refFile = TFile::Open(refFileName.Data());
587 Char_t refTreeName[100];
589 Int_t pdg[2],mumpdg[2],mumlab[2];
592 // read-in reference tree for event 0 (the only event for Pb-Pb)
593 sprintf(refTreeName,"Tree_Ref_%d",0);
594 TTree *refTree0 = (TTree*)refFile->Get(refTreeName);
595 refTree0->SetBranchAddress("rectracks",&reftrk);
597 AliD0toKpi *theD0 = 0;
598 treeD0in->SetBranchAddress("D0toKpi",&theD0);
599 treeD0out->Branch("D0toKpi","AliD0toKpi",&theD0,200000,0);
601 Int_t entries = (Int_t)treeD0in->GetEntries();
603 for(Int_t i=0; i<entries; i++) {
604 if(i%100==0) printf(" done %d candidates of %d\n",i,entries);
606 treeD0in->GetEvent(i);
607 event = theD0->EventNo();
609 if(event==0) { // always true for Pb-Pb (avoid to read-in tree every time)
610 refTree0->GetEvent(theD0->GetTrkNum(0));
612 mumpdg[0] = reftrk.mumpdg;
613 mumlab[0] = reftrk.mumlab;
614 refTree0->GetEvent(theD0->GetTrkNum(1));
616 mumpdg[1] = reftrk.mumpdg;
617 mumlab[1] = reftrk.mumlab;
619 sprintf(refTreeName,"Tree_Ref_%d",event);
620 TTree *refTree = (TTree*)refFile->Get(refTreeName);
621 refTree->SetBranchAddress("rectracks",&reftrk);
622 refTree->GetEvent(theD0->GetTrkNum(0));
624 mumpdg[0] = reftrk.mumpdg;
625 mumlab[0] = reftrk.mumlab;
626 refTree->GetEvent(theD0->GetTrkNum(1));
628 mumpdg[1] = reftrk.mumpdg;
629 mumlab[1] = reftrk.mumlab;
633 theD0->SetPdgCodes(pdg);
634 theD0->SetMumPdgCodes(mumpdg);
636 if(TMath::Abs(mumpdg[0])==421 &&
637 TMath::Abs(mumpdg[1])==421 &&
638 mumlab[0]==mumlab[1] &&
639 reftrk.mumprongs==2 &&
640 ((TMath::Abs(pdg[0])==211 && TMath::Abs(pdg[1])==321) ||
641 (TMath::Abs(pdg[0])==321 && TMath::Abs(pdg[1])==211))
642 ) theD0->SetSignal();
644 if(!fOnlySignal || theD0->IsSignal()) treeD0out->Fill();