GetESD.sh becomes runtest.sh
[u/mrichter/AliRoot.git] / ANALYSIS / AliD0toKpiAnalysis.cxx
CommitLineData
3a9a3487 1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
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 **************************************************************************/
15
16//----------------------------------------------------------------------------
17// Implementation of the D0toKpi reconstruction and analysis class
3a9a3487 18// Note: the two decay tracks are labelled: 0 (positive track)
19// 1 (negative track)
e6e6531b 20// An example of usage can be found in the macro AliD0toKpiTest.C
ef0182f7 21// Origin: A. Dainese andrea.dainese@lnl.infn.it
3a9a3487 22//----------------------------------------------------------------------------
3a9a3487 23#include <TKey.h>
24#include <TFile.h>
25#include <TTree.h>
26#include <TString.h>
27#include <TSystem.h>
28#include <TParticle.h>
29#include "AliESD.h"
ef0182f7 30#include "AliMC.h"
31#include "AliRun.h"
3a9a3487 32#include "AliRunLoader.h"
ef0182f7 33#include "AliVertexerTracks.h"
d681bb2d 34#include "AliESDVertex.h"
c7bafca9 35#include "AliESDv0.h"
3a9a3487 36#include "AliD0toKpi.h"
37#include "AliD0toKpiAnalysis.h"
59c6807d 38#include "AliLog.h"
3a9a3487 39
40typedef struct {
41 Int_t lab;
42 Int_t pdg;
43 Int_t mumlab;
44 Int_t mumpdg;
ef0182f7 45 Int_t mumprongs;
3a9a3487 46 Float_t Vx,Vy,Vz;
47 Float_t Px,Py,Pz;
48} REFTRACK;
49
50ClassImp(AliD0toKpiAnalysis)
51
52//----------------------------------------------------------------------------
53AliD0toKpiAnalysis::AliD0toKpiAnalysis() {
54 // Default constructor
ef0182f7 55
3a9a3487 56 SetPtCut();
57 Setd0Cut();
58 SetMassCut();
59 SetD0Cuts();
60 SetVertex1();
61 SetPID();
62 fVertexOnTheFly = kFALSE;
63 fSim = kFALSE;
64 fOnlySignal = kFALSE;
3a9a3487 65}
66//----------------------------------------------------------------------------
67AliD0toKpiAnalysis::~AliD0toKpiAnalysis() {}
68//----------------------------------------------------------------------------
69void AliD0toKpiAnalysis::ApplySelection(const Char_t *inName,const Char_t *outName) const {
70 // select candidates that pass fD0Cuts and write them to a new file
71
72 TFile *inFile = TFile::Open(inName);
73
74 TTree *treeD0in=(TTree*)inFile->Get("TreeD0");
75 AliD0toKpiAnalysis *inAnalysis = (AliD0toKpiAnalysis*)inFile->Get("D0toKpiAnalysis");
76 printf("+++\n+++ I N P U T S T A T U S:\n+++\n");
77 inAnalysis->PrintStatus();
78
79
80 AliD0toKpi *d = 0;
81 treeD0in->SetBranchAddress("D0toKpi",&d);
82 Int_t entries = (Int_t)treeD0in->GetEntries();
83
84 printf("+++\n+++ Number of D0 in input tree: %d\n+++\n",entries);
85
86 TTree *treeD0out = new TTree("TreeD0","Tree with selected D0 candidates");
87 treeD0out->Branch("D0toKpi","AliD0toKpi",&d,200000,0);
88
89
90 Int_t okD0=0,okD0bar=0;
91 Int_t nSel = 0;
92
93 for(Int_t i=0; i<entries; i++) {
94 // get event from tree
95 treeD0in->GetEvent(i);
96
97 if(fSim && fOnlySignal && !d->IsSignal()) continue;
98
99 // check if candidate passes selection (as D0 or D0bar)
100 if(d->Select(fD0Cuts,okD0,okD0bar)) {
101 nSel++;
102 treeD0out->Fill();
103 }
104
105 }
106
107 AliD0toKpiAnalysis *outAnalysis = (AliD0toKpiAnalysis*)inAnalysis->Clone("D0toKpiAnalysis");
108 outAnalysis->SetD0Cuts(fD0Cuts);
109 printf("------------------------------------------\n");
110 printf("+++\n+++ O U T P U T S T A T U S:\n+++\n");
111 outAnalysis->PrintStatus();
112
113 printf("+++\n+++ Number of D0 in output tree: %d\n+++\n",nSel);
114
115 TFile* outFile = new TFile(outName,"recreate");
116 treeD0out->Write();
117 outAnalysis->Write();
118 outFile->Close();
119
120 return;
121}
122//----------------------------------------------------------------------------
123Double_t AliD0toKpiAnalysis::CalculateTOFmass(Double_t mom,Double_t length,
124 Double_t time) const {
125 // calculated the mass from momentum, track length from vertex to TOF
126 // and time measured by the TOF
127 if(length==0.) return -1000.;
128 Double_t a = time*time/length/length;
129 if(a > 1.) {
130 a = TMath::Sqrt(a-1.);
131 } else {
132 a = -TMath::Sqrt(1.-a);
133 }
134
135 return mom*a;
136}
137//----------------------------------------------------------------------------
138void AliD0toKpiAnalysis::FindCandidates(Int_t evFirst,Int_t evLast,
139 const Char_t *outName) {
140 // Find D0 candidates and calculate parameters
141
3a9a3487 142
ef0182f7 143 TString esdName="AliESDs.root";
3a9a3487 144 if(gSystem->AccessPathName(esdName.Data(),kFileExists)) {
145 printf("AliD0toKpiAnalysis::FindCandidatesESD(): No ESDs file found!\n");
146 return;
147 }
148
149 TString outName1=outName;
ef0182f7 150 TString outName2="nTotEvents.dat";
3a9a3487 151
3a9a3487 152 Int_t nTotEv=0,nD0rec=0,nD0rec1ev=0;
153 Double_t dca;
154 Double_t v2[3],mom[6],d0[2];
3a9a3487 155 Int_t iTrkP,iTrkN,trkEntries;
156 Int_t nTrksP=0,nTrksN=0;
157 Int_t trkNum[2];
158 Double_t tofmass[2];
159 Int_t okD0=0,okD0bar=0;
c7bafca9 160 AliESDtrack *postrack = 0;
161 AliESDtrack *negtrack = 0;
3a9a3487 162
ef0182f7 163 // create the AliVertexerTracks object
3a9a3487 164 // (it will be used only if fVertexOnTheFly=kTrue)
ef0182f7 165 AliVertexerTracks *vertexer1 = new AliVertexerTracks;
166 if(fVertexOnTheFly) {
167 // open the mean vertex
168 TFile *invtx = new TFile("AliESDVertexMean.root");
169 AliESDVertex *initVertex = (AliESDVertex*)invtx->Get("vtxmean");
170 invtx->Close();
171 vertexer1->SetVtxStart(initVertex);
172 delete invtx;
173 }
3a9a3487 174 Int_t skipped[2];
175 Bool_t goodVtx1;
176
3a9a3487 177 // create tree for reconstructed D0s
178 AliD0toKpi *ioD0toKpi=0;
179 TTree *treeD0 = new TTree("TreeD0","Tree with D0 candidates");
180 treeD0->Branch("D0toKpi","AliD0toKpi",&ioD0toKpi,200000,0);
181
182 // open file with tracks
183 TFile *esdFile = TFile::Open(esdName.Data());
e6e6531b 184 AliESD* event = new AliESD;
185 TTree* tree = (TTree*) esdFile->Get("esdTree");
ef0182f7 186 if(!tree) {
e6e6531b 187 Error("FindCandidatesESD", "no ESD tree found");
188 return;
189 }
ef0182f7 190 tree->SetBranchAddress("ESD",&event);
3a9a3487 191
3a9a3487 192 // loop on events in file
ef0182f7 193 for(Int_t iEvent = evFirst; iEvent < tree->GetEntries(); iEvent++) {
194 if(iEvent > evLast) break;
e6e6531b 195 tree->GetEvent(iEvent);
31fd97b2 196 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.
3a9a3487 197 printf("--- Finding D0 -> Kpi in event %d\n",ev);
ef0182f7 198 // count the total number of events
199 nTotEv++;
3a9a3487 200
ef0182f7 201 trkEntries = (Int_t)event->GetNumberOfTracks();
3a9a3487 202 printf(" Number of tracks: %d\n",trkEntries);
ef0182f7 203 if(trkEntries<2) continue;
3a9a3487 204
ef0182f7 205 // retrieve primary vertex from file
206 if(!event->GetPrimaryVertex()) {
207 printf(" No vertex\n");
208 continue;
209 }
210 event->GetPrimaryVertex()->GetXYZ(fV1);
3a9a3487 211
212 // call function which applies sigle-track selection and
213 // separetes positives and negatives
214 TObjArray trksP(trkEntries/2);
215 Int_t *trkEntryP = new Int_t[trkEntries];
216 TObjArray trksN(trkEntries/2);
217 Int_t *trkEntryN = new Int_t[trkEntries];
218 TTree *trkTree = new TTree();
ef0182f7 219 SelectTracks(event,trksP,trkEntryP,nTrksP,
220 trksN,trkEntryN,nTrksN);
221
222 printf(" pos. tracks: %d neg .tracks: %d\n",nTrksP,nTrksN);
3a9a3487 223
3a9a3487 224
225 nD0rec1ev = 0;
226
ef0182f7 227 // LOOP ON POSITIVE TRACKS
3a9a3487 228 for(iTrkP=0; iTrkP<nTrksP; iTrkP++) {
ef0182f7 229 if(iTrkP%1==0) printf(" Processing positive track number %d of %d\n",iTrkP,nTrksP);
3a9a3487 230
231 // get track from track array
c7bafca9 232 postrack = (AliESDtrack*)trksP.UncheckedAt(iTrkP);
3a9a3487 233 trkNum[0] = trkEntryP[iTrkP];
234
ef0182f7 235 // LOOP ON NEGATIVE TRACKS
3a9a3487 236 for(iTrkN=0; iTrkN<nTrksN; iTrkN++) {
237
238 // get track from tracks array
c7bafca9 239 negtrack = (AliESDtrack*)trksN.UncheckedAt(iTrkN);
3a9a3487 240 trkNum[1] = trkEntryN[iTrkN];
241
3a9a3487 242 //
243 // ----------- DCA MINIMIZATION ------------------
244 //
c7bafca9 245 // find the DCA and propagate the tracks to the DCA
246 Double_t b=event->GetMagneticField();
247 AliESDtrack nt(*negtrack), pt(*postrack);
248 dca = nt.PropagateToDCA(&pt,b);
3a9a3487 249
c7bafca9 250 // define the AliESDv0 object
251 AliESDv0 vertex2(nt,trkNum[0],pt,trkNum[1]);
3a9a3487 252
253 // get position of the secondary vertex
c7bafca9 254 vertex2.GetXYZ(v2[0],v2[1],v2[2]);
255 vertex2.GetPPxPyPz(mom[0],mom[1],mom[2]);
256 vertex2.GetNPxPyPz(mom[3],mom[4],mom[5]);
257 // impact parameters of the tracks w.r.t. the primary vertex
c7bafca9 258 d0[0] = 10000.*pt.GetD(fV1[0],fV1[1],b);
259 d0[1] = -10000.*nt.GetD(fV1[0],fV1[1],b);
3a9a3487 260 goodVtx1 = kTRUE;
ef0182f7 261
3a9a3487 262 // no vertexing if DeltaMass > fMassCut
263 if(fVertexOnTheFly) {
264 goodVtx1 = kFALSE;
265 if(SelectInvMass(mom)) {
ef0182f7 266 // primary vertex from *other* tracks in the event
3a9a3487 267 skipped[0] = trkEntryP[iTrkP];
268 skipped[1] = trkEntryN[iTrkN];
269 vertexer1->SetSkipTracks(2,skipped);
d681bb2d 270 AliESDVertex *vertex1onfly =
ef0182f7 271 (AliESDVertex*)vertexer1->FindPrimaryVertex(event);
3a9a3487 272 if(vertex1onfly->GetNContributors()>0) goodVtx1 = kTRUE;
273 vertex1onfly->GetXYZ(fV1);
274 //vertex1onfly->PrintStatus();
275 delete vertex1onfly;
276 }
ef0182f7 277 }
278
3a9a3487 279
3a9a3487 280 // create the object AliD0toKpi
281 AliD0toKpi theD0(ev,trkNum,fV1,v2,dca,mom,d0);
3a9a3487 282 // select D0s
283 if(goodVtx1 && theD0.Select(fD0Cuts,okD0,okD0bar)) {
3a9a3487 284 // get PID info from ESD
285 AliESDtrack *t0 = (AliESDtrack*)event->GetTrack(trkNum[0]);
286 Double_t esdpid0[5];
287 t0->GetESDpid(esdpid0);
288 if(t0->GetStatus()&AliESDtrack::kTOFpid) {
289 tofmass[0] = CalculateTOFmass(t0->GetP(),
290 t0->GetIntegratedLength(),
291 t0->GetTOFsignal());
292 } else {
293 tofmass[0] = -1000.;
294 }
295 AliESDtrack *t1 = (AliESDtrack*)event->GetTrack(trkNum[1]);
296 Double_t esdpid1[5];
297 t1->GetESDpid(esdpid1);
298 if(t1->GetStatus()&AliESDtrack::kTOFpid) {
299 tofmass[1] = CalculateTOFmass(t1->GetP(),
300 t1->GetIntegratedLength(),
301 t1->GetTOFsignal());
302 } else {
303 tofmass[1] = -1000.;
304 }
305
306 theD0.SetPIDresponse(esdpid0,esdpid1);
307 theD0.SetTOFmasses(tofmass);
308
309 // fill the tree
310 ioD0toKpi=&theD0;
311 treeD0->Fill();
312
313 nD0rec++; nD0rec1ev++;
314 ioD0toKpi=0;
315 }
316
317 negtrack = 0;
318 } // loop on negative tracks
319 postrack = 0;
320 } // loop on positive tracks
ef0182f7 321
3a9a3487 322 delete [] trkEntryP;
323 delete [] trkEntryN;
324 delete trkTree;
3a9a3487 325
3a9a3487 326 printf(" Number of D0 candidates: %d\n",nD0rec1ev);
327 } // loop on events in file
328
329
330 printf("\n+++\n+++ Total number of events: %d\n+++\n",nTotEv);
331 printf("\n+++\n+++ Total number of D0 candidates: %d\n+++\n",nD0rec);
332
333 delete vertexer1;
3a9a3487 334
335 esdFile->Close();
336
337 // create a copy of this class to be written to output file
338 AliD0toKpiAnalysis *copy = (AliD0toKpiAnalysis*)this->Clone("D0toKpiAnalysis");
339
340 // add PDG codes to decay tracks in found candidates (in simulation mode)
341 // and store tree in the output file
342 if(!fSim) {
343 TFile *outroot = new TFile(outName1.Data(),"recreate");
344 treeD0->Write();
345 copy->Write();
346 outroot->Close();
347 delete outroot;
348 } else {
349 printf(" Now adding information from simulation (PDG codes) ...\n");
350 TTree *treeD0sim = new TTree("TreeD0","Tree with D0 candidates");
3a9a3487 351 SimulationInfo(treeD0,treeD0sim);
352 delete treeD0;
353 TFile *outroot = new TFile(outName1.Data(),"recreate");
354 treeD0sim->Write();
355 copy->Write();
356 outroot->Close();
357 delete outroot;
358 }
359
360 // write to a file the total number of events
361 FILE *outdat = fopen(outName2.Data(),"w");
362 fprintf(outdat,"%d\n",nTotEv);
363 fclose(outdat);
364
365 return;
366}
367//-----------------------------------------------------------------------------
368void AliD0toKpiAnalysis::PrintStatus() const {
369 // Print parameters being used
370
3a9a3487 371 printf("Preselections:\n");
372 printf(" fPtCut = %f GeV\n",fPtCut);
373 printf(" fd0Cut = %f micron\n",fd0Cut);
374 printf(" fMassCut = %f GeV\n",fMassCut);
375 if(fVertexOnTheFly) printf("Primary vertex on the fly\n");
376 if(fSim) {
377 printf("Simulation mode\n");
378 if(fOnlySignal) printf(" Only signal goes to file\n");
379 }
380 printf("Cuts on candidates:\n");
381 printf(" |M-MD0| [GeV] < %f\n",fD0Cuts[0]);
382 printf(" dca [micron] < %f\n",fD0Cuts[1]);
383 printf(" cosThetaStar < %f\n",fD0Cuts[2]);
384 printf(" pTK [GeV] > %f\n",fD0Cuts[3]);
385 printf(" pTpi [GeV] > %f\n",fD0Cuts[4]);
386 printf(" |d0K| [micron] < %f\n",fD0Cuts[5]);
387 printf(" |d0pi| [micron] < %f\n",fD0Cuts[6]);
388 printf(" d0d0 [micron^2] < %f\n",fD0Cuts[7]);
389 printf(" cosThetaPoint > %f\n",fD0Cuts[8]);
390
391 return;
392}
393//-----------------------------------------------------------------------------
394Bool_t AliD0toKpiAnalysis::SelectInvMass(const Double_t p[6]) const {
395 // Apply preselection in the invariant mass of the pair
396
397 Double_t mD0 = 1.8645;
398 Double_t mPi = 0.13957;
399 Double_t mKa = 0.49368;
400
401 Double_t energy[2];
402 Double_t mom2[2],momTot2;
403
404 mom2[0] = p[0]*p[0] + p[1]*p[1] + p[2]*p[2];
405 mom2[1] = p[3]*p[3] + p[4]*p[4] + p[5]*p[5];
406
407 momTot2 = (p[0]+p[3])*(p[0]+p[3])+
408 (p[1]+p[4])*(p[1]+p[4])+
409 (p[2]+p[5])*(p[2]+p[5]);
410
411 // D0 -> K- pi+
412 energy[1] = TMath::Sqrt(mKa*mKa+mom2[1]);
413 energy[0] = TMath::Sqrt(mPi*mPi+mom2[0]);
414
415 Double_t minvD0 = TMath::Sqrt((energy[0]+energy[1])*(energy[0]+energy[1])-momTot2);
416
417 // D0bar -> K+ pi-
418 energy[0] = TMath::Sqrt(mKa*mKa+mom2[0]);
419 energy[1] = TMath::Sqrt(mPi*mPi+mom2[1]);
420
421 Double_t minvD0bar = TMath::Sqrt((energy[0]+energy[1])*(energy[0]+energy[1])-momTot2);
422
423 if(TMath::Abs(minvD0-mD0) < fMassCut) return kTRUE;
424 if(TMath::Abs(minvD0bar-mD0) < fMassCut) return kTRUE;
425 return kFALSE;
426}
427//-----------------------------------------------------------------------------
ef0182f7 428void AliD0toKpiAnalysis::SelectTracks(AliESD *event,
3a9a3487 429 TObjArray &trksP,Int_t *trkEntryP,Int_t &nTrksP,
430 TObjArray &trksN,Int_t *trkEntryN,Int_t &nTrksN) const {
431 // Create two TObjArrays with positive and negative tracks and
432 // apply single-track preselection
433
434 nTrksP=0,nTrksN=0;
435
ef0182f7 436 Int_t entr = event->GetNumberOfTracks();
3a9a3487 437
438 // transfer ITS tracks from ESD to arrays and to a tree
439 for(Int_t i=0; i<entr; i++) {
3a9a3487 440
ef0182f7 441 AliESDtrack *esdtrack = event->GetTrack(i);
3a9a3487 442 UInt_t status = esdtrack->GetStatus();
443
ef0182f7 444 if(!(status&AliESDtrack::kITSin)) continue;
3a9a3487 445
3a9a3487 446 // single track selection
ef0182f7 447 if(!SingleTrkCuts(*esdtrack,event->GetMagneticField())) continue;
3a9a3487 448
c7bafca9 449 if(esdtrack->GetSign()<0) { // negative track
450 trksN.AddLast(esdtrack);
3a9a3487 451 trkEntryN[nTrksN] = i;
452 nTrksN++;
453 } else { // positive track
c7bafca9 454 trksP.AddLast(esdtrack);
3a9a3487 455 trkEntryP[nTrksP] = i;
456 nTrksP++;
457 }
458
459 } // loop on ESD tracks
460
461 return;
462}
463//-----------------------------------------------------------------------------
3a9a3487 464void AliD0toKpiAnalysis::SetD0Cuts(Double_t cut0,Double_t cut1,
465 Double_t cut2,Double_t cut3,Double_t cut4,
466 Double_t cut5,Double_t cut6,
467 Double_t cut7,Double_t cut8) {
468 // Set the cuts for D0 selection
469 fD0Cuts[0] = cut0;
470 fD0Cuts[1] = cut1;
471 fD0Cuts[2] = cut2;
472 fD0Cuts[3] = cut3;
473 fD0Cuts[4] = cut4;
474 fD0Cuts[5] = cut5;
475 fD0Cuts[6] = cut6;
476 fD0Cuts[7] = cut7;
477 fD0Cuts[8] = cut8;
478
479 return;
480}
481//-----------------------------------------------------------------------------
482void AliD0toKpiAnalysis::SetD0Cuts(const Double_t cuts[9]) {
483 // Set the cuts for D0 selection
484
485 for(Int_t i=0; i<9; i++) fD0Cuts[i] = cuts[i];
486
487 return;
488}
489//-----------------------------------------------------------------------------
c7bafca9 490Bool_t
491AliD0toKpiAnalysis::SingleTrkCuts(const AliESDtrack& trk, Double_t b) const {
3a9a3487 492 // Check if track passes some kinematical cuts
c7bafca9 493 // Magnetic field "b" (kG)
3a9a3487 494
c7bafca9 495 if(TMath::Abs(1./trk.GetParameter()[4]) < fPtCut)
3a9a3487 496 return kFALSE;
c7bafca9 497 if(TMath::Abs(10000.*trk.GetD(fV1[0],fV1[1],b)) < fd0Cut)
3a9a3487 498 return kFALSE;
499
500 return kTRUE;
501}
502//----------------------------------------------------------------------------
ef0182f7 503void AliD0toKpiAnalysis::MakeTracksRefFile(AliRun *gAlice,
504 Int_t evFirst,Int_t evLast) const {
3a9a3487 505 // Create a file with simulation info for the reconstructed tracks
506
ef0182f7 507 TFile *outFile = TFile::Open("D0TracksRefFile.root","recreate");
3a9a3487 508 TFile *esdFile = TFile::Open("AliESDs.root");
3a9a3487 509
ef0182f7 510 AliMC *mc = gAlice->GetMCApp();
3a9a3487 511
512 Int_t label;
513 TParticle *part;
514 TParticle *mumpart;
515 REFTRACK reftrk;
516
ef0182f7 517 AliESD* event = new AliESD;
518 TTree* tree = (TTree*) esdFile->Get("esdTree");
519 tree->SetBranchAddress("ESD",&event);
3a9a3487 520 // loop on events in file
ef0182f7 521 for(Int_t iEvent=evFirst; iEvent<tree->GetEntries(); iEvent++) {
522 if(iEvent>evLast) break;
523 tree->GetEvent(iEvent);
31fd97b2 524 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.
3a9a3487 525
ef0182f7 526 gAlice->GetEvent(ev);
3a9a3487 527
528 Int_t nentr=(Int_t)event->GetNumberOfTracks();
529
530 // Tree for true track parameters
531 char ttname[100];
532 sprintf(ttname,"Tree_Ref_%d",ev);
533 TTree *reftree = new TTree(ttname,"Tree with true track params");
534 reftree->Branch("rectracks",&reftrk,"lab/I:pdg:mumlab:mumpdg:Vx/F:Vy:Vz:Px:Py:Pz");
535
536 for(Int_t i=0; i<nentr; i++) {
537 AliESDtrack *esdtrack = (AliESDtrack*)event->GetTrack(i);
538 label = TMath::Abs(esdtrack->GetLabel());
539
ef0182f7 540 part = (TParticle*)mc->Particle(label);
3a9a3487 541 reftrk.lab = label;
542 reftrk.pdg = part->GetPdgCode();
543 reftrk.mumlab = part->GetFirstMother();
544 if(part->GetFirstMother()>=0) {
ef0182f7 545 mumpart = (TParticle*)gAlice->GetMCApp()->Particle(part->GetFirstMother());
3a9a3487 546 reftrk.mumpdg = mumpart->GetPdgCode();
ef0182f7 547 reftrk.mumprongs = mumpart->GetNDaughters();
3a9a3487 548 } else {
549 reftrk.mumpdg=-1;
ef0182f7 550 reftrk.mumprongs=-1;
3a9a3487 551 }
552 reftrk.Vx = part->Vx();
553 reftrk.Vy = part->Vy();
554 reftrk.Vz = part->Vz();
555 reftrk.Px = part->Px();
556 reftrk.Py = part->Py();
557 reftrk.Pz = part->Pz();
558
559 reftree->Fill();
560
561 } // loop on tracks
562
563 outFile->cd();
564 reftree->Write();
565
566 delete reftree;
3a9a3487 567 } // loop on events
568
569 esdFile->Close();
570 outFile->Close();
571
572 return;
573}
574//-----------------------------------------------------------------------------
575void AliD0toKpiAnalysis::SimulationInfo(TTree *treeD0in,TTree *treeD0out) const {
576 // add pdg codes to candidate decay tracks (for sim)
577
ef0182f7 578 TString refFileName("D0TracksRefFile.root");
3a9a3487 579 if(fSim && gSystem->AccessPathName(refFileName.Data(),kFileExists)) {
580 printf("AliD0toKpiAnalysis::SimulationInfo: no reference file found!\n");
581 return;
582 }
583 TFile *refFile = TFile::Open(refFileName.Data());
584
585 Char_t refTreeName[100];
586 Int_t event;
587 Int_t pdg[2],mumpdg[2],mumlab[2];
588 REFTRACK reftrk;
589
590 // read-in reference tree for event 0 (the only event for Pb-Pb)
591 sprintf(refTreeName,"Tree_Ref_%d",0);
592 TTree *refTree0 = (TTree*)refFile->Get(refTreeName);
593 refTree0->SetBranchAddress("rectracks",&reftrk);
594
595 AliD0toKpi *theD0 = 0;
596 treeD0in->SetBranchAddress("D0toKpi",&theD0);
597 treeD0out->Branch("D0toKpi","AliD0toKpi",&theD0,200000,0);
598
599 Int_t entries = (Int_t)treeD0in->GetEntries();
600
601 for(Int_t i=0; i<entries; i++) {
602 if(i%100==0) printf(" done %d candidates of %d\n",i,entries);
603
604 treeD0in->GetEvent(i);
605 event = theD0->EventNo();
606
607 if(event==0) { // always true for Pb-Pb (avoid to read-in tree every time)
608 refTree0->GetEvent(theD0->GetTrkNum(0));
609 pdg[0] = reftrk.pdg;
610 mumpdg[0] = reftrk.mumpdg;
611 mumlab[0] = reftrk.mumlab;
612 refTree0->GetEvent(theD0->GetTrkNum(1));
613 pdg[1] = reftrk.pdg;
614 mumpdg[1] = reftrk.mumpdg;
615 mumlab[1] = reftrk.mumlab;
616 } else {
617 sprintf(refTreeName,"Tree_Ref_%d",event);
618 TTree *refTree = (TTree*)refFile->Get(refTreeName);
619 refTree->SetBranchAddress("rectracks",&reftrk);
620 refTree->GetEvent(theD0->GetTrkNum(0));
621 pdg[0] = reftrk.pdg;
622 mumpdg[0] = reftrk.mumpdg;
623 mumlab[0] = reftrk.mumlab;
624 refTree->GetEvent(theD0->GetTrkNum(1));
625 pdg[1] = reftrk.pdg;
626 mumpdg[1] = reftrk.mumpdg;
627 mumlab[1] = reftrk.mumlab;
628 delete refTree;
629 }
630
631 theD0->SetPdgCodes(pdg);
632 theD0->SetMumPdgCodes(mumpdg);
633
ef0182f7 634 if(TMath::Abs(mumpdg[0])==421 &&
635 TMath::Abs(mumpdg[1])==421 &&
636 mumlab[0]==mumlab[1] &&
637 reftrk.mumprongs==2 &&
638 ((TMath::Abs(pdg[0])==211 && TMath::Abs(pdg[1])==321) ||
639 (TMath::Abs(pdg[0])==321 && TMath::Abs(pdg[1])==211))
640 ) theD0->SetSignal();
3a9a3487 641
642 if(!fOnlySignal || theD0->IsSignal()) treeD0out->Fill();
643
644 }
645
646 delete refTree0;
647
648 refFile->Close();
649
650 return;
651}
652
653
654
655
656
657