memory consumption measures only memory used by preprocessor
[u/mrichter/AliRoot.git] / PWG3 / vertexingOld / AliD0toKpiAnalysis.cxx
CommitLineData
9af36e3a 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
18// Note: the two decay tracks are labelled: 0 (positive track)
19// 1 (negative track)
20// An example of usage can be found in the macro AliD0toKpiTest.C
21// Origin: A. Dainese andrea.dainese@lnl.infn.it
22//----------------------------------------------------------------------------
23#include <TKey.h>
24#include <TFile.h>
25#include <TTree.h>
26#include <TString.h>
27#include <TSystem.h>
28#include <TParticle.h>
5e6f195c 29#include "AliESDEvent.h"
9af36e3a 30#include "AliMC.h"
31#include "AliRun.h"
32#include "AliRunLoader.h"
33#include "AliVertexerTracks.h"
34#include "AliESDVertex.h"
35#include "AliESDv0.h"
36#include "AliD0toKpi.h"
37#include "AliD0toKpiAnalysis.h"
38#include "AliLog.h"
39
40typedef struct {
41 Int_t lab;
42 Int_t pdg;
43 Int_t mumlab;
44 Int_t mumpdg;
45 Int_t mumprongs;
46 Float_t Vx,Vy,Vz;
47 Float_t Px,Py,Pz;
48} REFTRACK;
49
50ClassImp(AliD0toKpiAnalysis)
51
52//----------------------------------------------------------------------------
dd1ac8a2 53AliD0toKpiAnalysis::AliD0toKpiAnalysis():
54fVertexOnTheFly(kFALSE),
55fSim(kFALSE),
56fOnlySignal(kFALSE),
57fPID("TOFparam_PbPb"),
58fPtCut(0.),
59fd0Cut(0.),
60fMassCut(1000.)
61 {
9af36e3a 62 // Default constructor
63
9af36e3a 64 SetD0Cuts();
65 SetVertex1();
9af36e3a 66}
67//----------------------------------------------------------------------------
68AliD0toKpiAnalysis::~AliD0toKpiAnalysis() {}
69//----------------------------------------------------------------------------
70void AliD0toKpiAnalysis::ApplySelection(const Char_t *inName,const Char_t *outName) const {
71 // select candidates that pass fD0Cuts and write them to a new file
72
73 TFile *inFile = TFile::Open(inName);
74
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();
79
80
81 AliD0toKpi *d = 0;
82 treeD0in->SetBranchAddress("D0toKpi",&d);
83 Int_t entries = (Int_t)treeD0in->GetEntries();
84
85 printf("+++\n+++ Number of D0 in input tree: %d\n+++\n",entries);
86
87 TTree *treeD0out = new TTree("TreeD0","Tree with selected D0 candidates");
88 treeD0out->Branch("D0toKpi","AliD0toKpi",&d,200000,0);
89
90
91 Int_t okD0=0,okD0bar=0;
92 Int_t nSel = 0;
93
94 for(Int_t i=0; i<entries; i++) {
95 // get event from tree
96 treeD0in->GetEvent(i);
97
98 if(fSim && fOnlySignal && !d->IsSignal()) continue;
99
100 // check if candidate passes selection (as D0 or D0bar)
101 if(d->Select(fD0Cuts,okD0,okD0bar)) {
102 nSel++;
103 treeD0out->Fill();
104 }
105
106 }
107
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();
113
114 printf("+++\n+++ Number of D0 in output tree: %d\n+++\n",nSel);
115
116 TFile* outFile = new TFile(outName,"recreate");
117 treeD0out->Write();
118 outAnalysis->Write();
119 outFile->Close();
120
121 return;
122}
123//----------------------------------------------------------------------------
124Double_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;
130 if(a > 1.) {
131 a = TMath::Sqrt(a-1.);
132 } else {
133 a = -TMath::Sqrt(1.-a);
134 }
135
136 return mom*a;
137}
138//----------------------------------------------------------------------------
139void AliD0toKpiAnalysis::FindCandidates(Int_t evFirst,Int_t evLast,
140 const Char_t *outName) {
141 // Find D0 candidates and calculate parameters
142
143
144 TString esdName="AliESDs.root";
145 if(gSystem->AccessPathName(esdName.Data(),kFileExists)) {
146 printf("AliD0toKpiAnalysis::FindCandidatesESD(): No ESDs file found!\n");
147 return;
148 }
149
150 TString outName1=outName;
151 TString outName2="nTotEvents.dat";
152
153 Int_t nTotEv=0,nD0rec=0,nD0rec1ev=0;
154 Double_t dca;
155 Double_t v2[3],mom[6],d0[2];
156 Int_t iTrkP,iTrkN,trkEntries;
157 Int_t nTrksP=0,nTrksN=0;
158 Int_t trkNum[2];
159 Double_t tofmass[2];
160 Int_t okD0=0,okD0bar=0;
161 AliESDtrack *postrack = 0;
162 AliESDtrack *negtrack = 0;
163
164 // create the AliVertexerTracks object
165 // (it will be used only if fVertexOnTheFly=kTrue)
dd1ac8a2 166 AliVertexerTracks *vertexer1 = new AliVertexerTracks();
9af36e3a 167 if(fVertexOnTheFly) {
168 // open the mean vertex
169 TFile *invtx = new TFile("AliESDVertexMean.root");
170 AliESDVertex *initVertex = (AliESDVertex*)invtx->Get("vtxmean");
171 invtx->Close();
172 vertexer1->SetVtxStart(initVertex);
173 delete invtx;
174 }
175 Int_t skipped[2];
176 Bool_t goodVtx1;
177
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);
182
183 // open file with tracks
184 TFile *esdFile = TFile::Open(esdName.Data());
5e6f195c 185 AliESDEvent* event = new AliESDEvent();
9af36e3a 186 TTree* tree = (TTree*) esdFile->Get("esdTree");
187 if(!tree) {
188 Error("FindCandidatesESD", "no ESD tree found");
189 return;
190 }
5e6f195c 191 event->ReadFromTree(tree);
9af36e3a 192
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
200 nTotEv++;
201
202 trkEntries = (Int_t)event->GetNumberOfTracks();
203 printf(" Number of tracks: %d\n",trkEntries);
204 if(trkEntries<2) continue;
205
206 // retrieve primary vertex from file
207 if(!event->GetPrimaryVertex()) {
208 printf(" No vertex\n");
209 continue;
210 }
211 event->GetPrimaryVertex()->GetXYZ(fV1);
212
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);
222
223 printf(" pos. tracks: %d neg .tracks: %d\n",nTrksP,nTrksN);
224
225
226 nD0rec1ev = 0;
227
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);
231
232 // get track from track array
233 postrack = (AliESDtrack*)trksP.UncheckedAt(iTrkP);
234 trkNum[0] = trkEntryP[iTrkP];
235
236 // LOOP ON NEGATIVE TRACKS
237 for(iTrkN=0; iTrkN<nTrksN; iTrkN++) {
238
239 // get track from tracks array
240 negtrack = (AliESDtrack*)trksN.UncheckedAt(iTrkN);
241 trkNum[1] = trkEntryN[iTrkN];
242
243 //
244 // ----------- DCA MINIMIZATION ------------------
245 //
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);
250
251 // define the AliESDv0 object
252 AliESDv0 vertex2(nt,trkNum[0],pt,trkNum[1]);
253
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);
261 goodVtx1 = kTRUE;
262
263 // no vertexing if DeltaMass > fMassCut
264 if(fVertexOnTheFly) {
265 goodVtx1 = kFALSE;
266 if(SelectInvMass(mom)) {
267 // primary vertex from *other* tracks in the event
dd1ac8a2 268 vertexer1->SetFieldkG(event->GetMagneticField());
9af36e3a 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();
277 delete vertex1onfly;
278 }
279 }
280
281
282 // create the object AliD0toKpi
283 AliD0toKpi theD0(ev,trkNum,fV1,v2,dca,mom,d0);
284 // select D0s
285 if(goodVtx1 && theD0.Select(fD0Cuts,okD0,okD0bar)) {
286 // get PID info from ESD
287 AliESDtrack *t0 = (AliESDtrack*)event->GetTrack(trkNum[0]);
288 Double_t esdpid0[5];
289 t0->GetESDpid(esdpid0);
290 if(t0->GetStatus()&AliESDtrack::kTOFpid) {
291 tofmass[0] = CalculateTOFmass(t0->GetP(),
292 t0->GetIntegratedLength(),
293 t0->GetTOFsignal());
294 } else {
295 tofmass[0] = -1000.;
296 }
297 AliESDtrack *t1 = (AliESDtrack*)event->GetTrack(trkNum[1]);
298 Double_t esdpid1[5];
299 t1->GetESDpid(esdpid1);
300 if(t1->GetStatus()&AliESDtrack::kTOFpid) {
301 tofmass[1] = CalculateTOFmass(t1->GetP(),
302 t1->GetIntegratedLength(),
303 t1->GetTOFsignal());
304 } else {
305 tofmass[1] = -1000.;
306 }
307
308 theD0.SetPIDresponse(esdpid0,esdpid1);
309 theD0.SetTOFmasses(tofmass);
310
311 // fill the tree
312 ioD0toKpi=&theD0;
313 treeD0->Fill();
314
315 nD0rec++; nD0rec1ev++;
316 ioD0toKpi=0;
317 }
318
319 negtrack = 0;
320 } // loop on negative tracks
321 postrack = 0;
322 } // loop on positive tracks
323
324 delete [] trkEntryP;
325 delete [] trkEntryN;
326 delete trkTree;
327
328 printf(" Number of D0 candidates: %d\n",nD0rec1ev);
329 } // loop on events in file
330
331
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);
334
335 delete vertexer1;
336
337 esdFile->Close();
338
339 // create a copy of this class to be written to output file
340 AliD0toKpiAnalysis *copy = (AliD0toKpiAnalysis*)this->Clone("D0toKpiAnalysis");
341
342 // add PDG codes to decay tracks in found candidates (in simulation mode)
343 // and store tree in the output file
344 if(!fSim) {
345 TFile *outroot = new TFile(outName1.Data(),"recreate");
346 treeD0->Write();
347 copy->Write();
348 outroot->Close();
349 delete outroot;
350 } else {
351 printf(" Now adding information from simulation (PDG codes) ...\n");
352 TTree *treeD0sim = new TTree("TreeD0","Tree with D0 candidates");
353 SimulationInfo(treeD0,treeD0sim);
354 delete treeD0;
355 TFile *outroot = new TFile(outName1.Data(),"recreate");
356 treeD0sim->Write();
357 copy->Write();
358 outroot->Close();
359 delete outroot;
360 }
361
362 // write to a file the total number of events
363 FILE *outdat = fopen(outName2.Data(),"w");
364 fprintf(outdat,"%d\n",nTotEv);
365 fclose(outdat);
366
367 return;
368}
369//-----------------------------------------------------------------------------
370void AliD0toKpiAnalysis::PrintStatus() const {
371 // Print parameters being used
372
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");
378 if(fSim) {
379 printf("Simulation mode\n");
380 if(fOnlySignal) printf(" Only signal goes to file\n");
381 }
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]);
392
393 return;
394}
395//-----------------------------------------------------------------------------
396Bool_t AliD0toKpiAnalysis::SelectInvMass(const Double_t p[6]) const {
397 // Apply preselection in the invariant mass of the pair
398
399 Double_t mD0 = 1.8645;
400 Double_t mPi = 0.13957;
401 Double_t mKa = 0.49368;
402
403 Double_t energy[2];
404 Double_t mom2[2],momTot2;
405
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];
408
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]);
412
413 // D0 -> K- pi+
414 energy[1] = TMath::Sqrt(mKa*mKa+mom2[1]);
415 energy[0] = TMath::Sqrt(mPi*mPi+mom2[0]);
416
417 Double_t minvD0 = TMath::Sqrt((energy[0]+energy[1])*(energy[0]+energy[1])-momTot2);
418
419 // D0bar -> K+ pi-
420 energy[0] = TMath::Sqrt(mKa*mKa+mom2[0]);
421 energy[1] = TMath::Sqrt(mPi*mPi+mom2[1]);
422
423 Double_t minvD0bar = TMath::Sqrt((energy[0]+energy[1])*(energy[0]+energy[1])-momTot2);
424
425 if(TMath::Abs(minvD0-mD0) < fMassCut) return kTRUE;
426 if(TMath::Abs(minvD0bar-mD0) < fMassCut) return kTRUE;
427 return kFALSE;
428}
429//-----------------------------------------------------------------------------
5e6f195c 430void AliD0toKpiAnalysis::SelectTracks(AliESDEvent *event,
9af36e3a 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
435
436 nTrksP=0,nTrksN=0;
437
438 Int_t entr = event->GetNumberOfTracks();
439
440 // transfer ITS tracks from ESD to arrays and to a tree
441 for(Int_t i=0; i<entr; i++) {
442
443 AliESDtrack *esdtrack = event->GetTrack(i);
444 UInt_t status = esdtrack->GetStatus();
445
446 if(!(status&AliESDtrack::kITSin)) continue;
447
448 // single track selection
449 if(!SingleTrkCuts(*esdtrack,event->GetMagneticField())) continue;
450
451 if(esdtrack->GetSign()<0) { // negative track
452 trksN.AddLast(esdtrack);
453 trkEntryN[nTrksN] = i;
454 nTrksN++;
455 } else { // positive track
456 trksP.AddLast(esdtrack);
457 trkEntryP[nTrksP] = i;
458 nTrksP++;
459 }
460
461 } // loop on ESD tracks
462
463 return;
464}
465//-----------------------------------------------------------------------------
466void 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
471 fD0Cuts[0] = cut0;
472 fD0Cuts[1] = cut1;
473 fD0Cuts[2] = cut2;
474 fD0Cuts[3] = cut3;
475 fD0Cuts[4] = cut4;
476 fD0Cuts[5] = cut5;
477 fD0Cuts[6] = cut6;
478 fD0Cuts[7] = cut7;
479 fD0Cuts[8] = cut8;
480
481 return;
482}
483//-----------------------------------------------------------------------------
484void AliD0toKpiAnalysis::SetD0Cuts(const Double_t cuts[9]) {
485 // Set the cuts for D0 selection
486
487 for(Int_t i=0; i<9; i++) fD0Cuts[i] = cuts[i];
488
489 return;
490}
491//-----------------------------------------------------------------------------
492Bool_t
493AliD0toKpiAnalysis::SingleTrkCuts(const AliESDtrack& trk, Double_t b) const {
494 // Check if track passes some kinematical cuts
495 // Magnetic field "b" (kG)
496
497 if(TMath::Abs(1./trk.GetParameter()[4]) < fPtCut)
498 return kFALSE;
499 if(TMath::Abs(10000.*trk.GetD(fV1[0],fV1[1],b)) < fd0Cut)
500 return kFALSE;
501
502 return kTRUE;
503}
504//----------------------------------------------------------------------------
71d22fe8 505void AliD0toKpiAnalysis::MakeTracksRefFile(AliRun *mygAlice,
9af36e3a 506 Int_t evFirst,Int_t evLast) const {
507 // Create a file with simulation info for the reconstructed tracks
508
509 TFile *outFile = TFile::Open("D0TracksRefFile.root","recreate");
510 TFile *esdFile = TFile::Open("AliESDs.root");
511
71d22fe8 512 AliMC *mc = mygAlice->GetMCApp();
9af36e3a 513
514 Int_t label;
515 TParticle *part;
516 TParticle *mumpart;
517 REFTRACK reftrk;
518
5e6f195c 519 AliESDEvent* event = new AliESDEvent();
9af36e3a 520 TTree* tree = (TTree*) esdFile->Get("esdTree");
5e6f195c 521 event->ReadFromTree(tree);
9af36e3a 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.
527
71d22fe8 528 mygAlice->GetEvent(ev);
9af36e3a 529
530 Int_t nentr=(Int_t)event->GetNumberOfTracks();
531
532 // Tree for true track parameters
533 char ttname[100];
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");
537
538 for(Int_t i=0; i<nentr; i++) {
539 AliESDtrack *esdtrack = (AliESDtrack*)event->GetTrack(i);
540 label = TMath::Abs(esdtrack->GetLabel());
541
542 part = (TParticle*)mc->Particle(label);
543 reftrk.lab = label;
544 reftrk.pdg = part->GetPdgCode();
545 reftrk.mumlab = part->GetFirstMother();
546 if(part->GetFirstMother()>=0) {
71d22fe8 547 mumpart = (TParticle*)mygAlice->GetMCApp()->Particle(part->GetFirstMother());
9af36e3a 548 reftrk.mumpdg = mumpart->GetPdgCode();
549 reftrk.mumprongs = mumpart->GetNDaughters();
550 } else {
551 reftrk.mumpdg=-1;
552 reftrk.mumprongs=-1;
553 }
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();
560
561 reftree->Fill();
562
563 } // loop on tracks
564
565 outFile->cd();
566 reftree->Write();
567
568 delete reftree;
569 } // loop on events
570
571 esdFile->Close();
572 outFile->Close();
573
574 return;
575}
576//-----------------------------------------------------------------------------
577void AliD0toKpiAnalysis::SimulationInfo(TTree *treeD0in,TTree *treeD0out) const {
578 // add pdg codes to candidate decay tracks (for sim)
579
580 TString refFileName("D0TracksRefFile.root");
581 if(fSim && gSystem->AccessPathName(refFileName.Data(),kFileExists)) {
582 printf("AliD0toKpiAnalysis::SimulationInfo: no reference file found!\n");
583 return;
584 }
585 TFile *refFile = TFile::Open(refFileName.Data());
586
587 Char_t refTreeName[100];
588 Int_t event;
589 Int_t pdg[2],mumpdg[2],mumlab[2];
590 REFTRACK reftrk;
591
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);
596
597 AliD0toKpi *theD0 = 0;
598 treeD0in->SetBranchAddress("D0toKpi",&theD0);
599 treeD0out->Branch("D0toKpi","AliD0toKpi",&theD0,200000,0);
600
601 Int_t entries = (Int_t)treeD0in->GetEntries();
602
603 for(Int_t i=0; i<entries; i++) {
604 if(i%100==0) printf(" done %d candidates of %d\n",i,entries);
605
606 treeD0in->GetEvent(i);
607 event = theD0->EventNo();
608
609 if(event==0) { // always true for Pb-Pb (avoid to read-in tree every time)
610 refTree0->GetEvent(theD0->GetTrkNum(0));
611 pdg[0] = reftrk.pdg;
612 mumpdg[0] = reftrk.mumpdg;
613 mumlab[0] = reftrk.mumlab;
614 refTree0->GetEvent(theD0->GetTrkNum(1));
615 pdg[1] = reftrk.pdg;
616 mumpdg[1] = reftrk.mumpdg;
617 mumlab[1] = reftrk.mumlab;
618 } else {
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));
623 pdg[0] = reftrk.pdg;
624 mumpdg[0] = reftrk.mumpdg;
625 mumlab[0] = reftrk.mumlab;
626 refTree->GetEvent(theD0->GetTrkNum(1));
627 pdg[1] = reftrk.pdg;
628 mumpdg[1] = reftrk.mumpdg;
629 mumlab[1] = reftrk.mumlab;
630 delete refTree;
631 }
632
633 theD0->SetPdgCodes(pdg);
634 theD0->SetMumPdgCodes(mumpdg);
635
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();
643
644 if(!fOnlySignal || theD0->IsSignal()) treeD0out->Fill();
645
646 }
647
648 delete refTree0;
649
650 refFile->Close();
651
652 return;
653}
654
655
656
657
658
659