]> git.uio.no Git - u/mrichter/AliRoot.git/blame - ANALYSIS/AliD0toKpiAnalysis.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Do not convert bool to pointer
[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
3a9a3487 21// Origin: A. Dainese andrea.dainese@pd.infn.it
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"
30#include "AliStack.h"
31#include "AliRunLoader.h"
32#include "AliITStrackV2.h"
33#include "AliITSVertexerTracks.h"
d681bb2d 34#include "AliESDVertex.h"
3a9a3487 35#include "AliV0vertexer.h"
36#include "AliV0vertex.h"
37#include "AliD0toKpi.h"
38#include "AliD0toKpiAnalysis.h"
59c6807d 39#include "AliLog.h"
3a9a3487 40
41typedef struct {
42 Int_t lab;
43 Int_t pdg;
44 Int_t mumlab;
45 Int_t mumpdg;
46 Float_t Vx,Vy,Vz;
47 Float_t Px,Py,Pz;
48} REFTRACK;
49
50ClassImp(AliD0toKpiAnalysis)
51
52//----------------------------------------------------------------------------
53AliD0toKpiAnalysis::AliD0toKpiAnalysis() {
54 // Default constructor
55
c84a5e9e 56 fBz=-9999;
3a9a3487 57 SetPtCut();
58 Setd0Cut();
59 SetMassCut();
60 SetD0Cuts();
61 SetVertex1();
62 SetPID();
63 fVertexOnTheFly = kFALSE;
64 fSim = kFALSE;
65 fOnlySignal = kFALSE;
3a9a3487 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 if(fBz<-9000.) {
144 printf("AliD0toKpiAnalysis::FindCandidates(): Set B!\n");
145 return;
146 }
3a9a3487 147
148 TString trkName("AliITStracksV2.root");
149 if(gSystem->AccessPathName(trkName.Data(),kFileExists)) {
150 printf("AliD0toKpiAnalysis::FindCandidates(): No tracks file found!\n");
151 return;
152 }
153
154 TString outName1=outName;
155 TString outName2("nTotEvents.dat");
156
157 Int_t ev;
158 Int_t nTotEv=0,nD0rec=0,nD0rec1ev=0;
159 Double_t dca;
160 Double_t v2[3],mom[6],d0[2];
161 Double_t alphaP,alphaN,ptP,ptN,phiP,phiN;
162 Int_t iTrkP,iTrkN,trkEntries;
163 Int_t nTrksP=0,nTrksN=0;
164 Int_t trkNum[2];
165 Int_t okD0=0,okD0bar=0;
166 Char_t trkTreeName[100],vtxName[100];
167 AliITStrackV2 *postrack = 0;
168 AliITStrackV2 *negtrack = 0;
169
170 // create the AliITSVertexerTracks object
171 // (it will be used only if fVertexOnTheFly=kTrue)
172 AliITSVertexerTracks *vertexer1 = new AliITSVertexerTracks;
173 vertexer1->SetMinTracks(2);
3a9a3487 174 Int_t skipped[2];
175 Bool_t goodVtx1;
176
177
178 // define the cuts for vertexing
179 Double_t vtxcuts[]={50., // max. allowed chi2
180 0.0, // min. allowed negative daughter's impact param
181 0.0, // min. allowed positive daughter's impact param
182 1.0, // max. allowed DCA between the daughter tracks
183 -1.0, // min. allowed cosine of pointing angle
184 0.0, // min. radius of the fiducial volume
185 2.9};// max. radius of the fiducial volume
186
187 // create the AliV0vertexer object
188 AliV0vertexer *vertexer2 = new AliV0vertexer(vtxcuts);
189
190 // create tree for reconstructed D0s
191 AliD0toKpi *ioD0toKpi=0;
192 TTree *treeD0 = new TTree("TreeD0","Tree with D0 candidates");
193 treeD0->Branch("D0toKpi","AliD0toKpi",&ioD0toKpi,200000,0);
194
195 // open file with tracks
196 TFile *trkFile = TFile::Open(trkName.Data());
197
198 // loop on events in file
199 for(ev=evFirst; ev<=evLast; ev++) {
200 printf(" --- Looking for D0->Kpi in event %d ---\n",ev);
201
202 // retrieve primary vertex from file
203 sprintf(vtxName,"VertexTracks_%d",ev);
d681bb2d 204 AliESDVertex *vertex1stored = (AliESDVertex*)trkFile->Get(vtxName);
3a9a3487 205 vertex1stored->GetXYZ(fV1);
206 delete vertex1stored;
207
208 // retrieve tracks from file
209 sprintf(trkTreeName,"TreeT_ITS_%d",ev);
210 TTree *trkTree=(TTree*)trkFile->Get(trkTreeName);
211 if(!trkTree) {
212 printf("AliD0toKpiAnalysis::FindCandidates():\n tracks tree not found for evet %d\n",ev);
213 continue;
214 }
215 trkEntries = (Int_t)trkTree->GetEntries();
216 printf(" Number of tracks: %d\n",trkEntries);
217
218 // count the total number of events
219 nTotEv++;
220
221 // call function which applies sigle-track selection and
222 // separetes positives and negatives
223 TObjArray trksP(trkEntries/2);
224 Int_t *trkEntryP = new Int_t[trkEntries];
225 TObjArray trksN(trkEntries/2);
226 Int_t *trkEntryN = new Int_t[trkEntries];
227 SelectTracks(*trkTree,trksP,trkEntryP,nTrksP,trksN,trkEntryN,nTrksN);
228
229 nD0rec1ev = 0;
230
231 // loop on positive tracks
232 for(iTrkP=0; iTrkP<nTrksP; iTrkP++) {
233 if(iTrkP%10==0) printf(" Processing positive track number %d of %d\n",iTrkP,nTrksP);
234
235 // get track from track array
236 postrack = (AliITStrackV2*)trksP.At(iTrkP);
237 trkNum[0] = trkEntryP[iTrkP];
238
239 // loop on negative tracks
240 for(iTrkN=0; iTrkN<nTrksN; iTrkN++) {
241 // get track from tracks array
242 negtrack = (AliITStrackV2*)trksN.At(iTrkN);
243 trkNum[1] = trkEntryN[iTrkN];
244
245 AliITStrackV2 nt(*negtrack), pt(*postrack), *pnt=&nt, *ppt=&pt;
246
247 //
248 // ----------- DCA MINIMIZATION ------------------
249 //
250 // find the DCA and propagate the tracks to the DCA
251 dca = vertexer2->PropagateToDCA(pnt,ppt);
252
253 // define the AliV0vertex object
254 AliV0vertex *vertex2 = new AliV0vertex(*pnt,*ppt);
255
256 // get position of the secondary vertex
257 vertex2->GetXYZ(v2[0],v2[1],v2[2]);
258
259 delete vertex2;
260
261 // momenta of the tracks at the vertex
262 ptP = 1./TMath::Abs(ppt->Get1Pt());
263 alphaP = ppt->GetAlpha();
264 phiP = alphaP+TMath::ASin(ppt->GetSnp());
265 mom[0] = ptP*TMath::Cos(phiP);
266 mom[1] = ptP*TMath::Sin(phiP);
267 mom[2] = ptP*ppt->GetTgl();
268
269 ptN = 1./TMath::Abs(pnt->Get1Pt());
270 alphaN = pnt->GetAlpha();
271 phiN = alphaN+TMath::ASin(pnt->GetSnp());
272 mom[3] = ptN*TMath::Cos(phiN);
273 mom[4] = ptN*TMath::Sin(phiN);
274 mom[5] = ptN*pnt->GetTgl();
275
276 goodVtx1 = kTRUE;
277 // no vertexing if DeltaMass > fMassCut
278 if(fVertexOnTheFly) {
279 goodVtx1 = kFALSE;
280 if(SelectInvMass(mom)) {
281 // primary vertex from *other* tracks in event
282 vertexer1->SetVtxStart(fV1[0],fV1[1]);
283 skipped[0] = trkEntryP[iTrkP];
284 skipped[1] = trkEntryN[iTrkN];
285 vertexer1->SetSkipTracks(2,skipped);
d681bb2d 286 AliESDVertex *vertex1onfly =
287 (AliESDVertex*)vertexer1->VertexOnTheFly(*trkTree);
3a9a3487 288 if(vertex1onfly->GetNContributors()>0) goodVtx1 = kTRUE;
289 vertex1onfly->GetXYZ(fV1);
290 //vertex1onfly->PrintStatus();
291 delete vertex1onfly;
292 }
293 }
294
295 // impact parameters of the tracks w.r.t. the primary vertex
296 d0[0] = 10000.*ppt->GetD(fV1[0],fV1[1]);
297 d0[1] = -10000.*pnt->GetD(fV1[0],fV1[1]);
298
299 // create the object AliD0toKpi
300 AliD0toKpi theD0(ev,trkNum,fV1,v2,dca,mom,d0);
301
302 // select D0s
303 if(goodVtx1 && theD0.Select(fD0Cuts,okD0,okD0bar)) {
304
305 // fill the tree
306 ioD0toKpi=&theD0;
307 treeD0->Fill();
308
309 nD0rec++; nD0rec1ev++;
310 ioD0toKpi=0;
311 }
312
313 negtrack = 0;
314 } // loop on negative tracks
315 postrack = 0;
316 } // loop on positive tracks
317
318 trksP.Delete();
319 trksN.Delete();
320 delete [] trkEntryP;
321 delete [] trkEntryN;
322 delete trkTree;
323
324 printf(" Number of D0 candidates: %d\n",nD0rec1ev);
325 } // loop on events in file
326
327
328 printf("\n+++\n+++ Total number of events: %d\n+++\n",nTotEv);
329 printf("\n+++\n+++ Total number of D0 candidates: %d\n+++\n",nD0rec);
330
331 delete vertexer1;
332 delete vertexer2;
333
334 trkFile->Close();
335
336 // create a copy of this class to be written to output file
337 AliD0toKpiAnalysis *copy = (AliD0toKpiAnalysis*)this->Clone("D0toKpiAnalysis");
338 copy->PrintStatus();
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");
351 MakeTracksRefFile(evFirst,evLast);
352 SimulationInfo(treeD0,treeD0sim);
353 delete treeD0;
354 TFile *outroot = new TFile(outName1.Data(),"recreate");
355 treeD0sim->Write();
356 copy->Write();
357 outroot->Close();
358 delete outroot;
359 }
360
361 // write to a file the total number of events
362 FILE *outdat = fopen(outName2.Data(),"w");
363 fprintf(outdat,"%d\n",nTotEv);
364 fclose(outdat);
365
366 return;
367}
368//----------------------------------------------------------------------------
369void AliD0toKpiAnalysis::FindCandidatesESD(Int_t evFirst,Int_t evLast,
370 const Char_t *outName) {
371 // Find D0 candidates and calculate parameters
372
373 if(fBz<-9000.) {
374 printf("AliD0toKpiAnalysis::FindCandidatesESD(): Set B!\n");
375 return;
376 }
3a9a3487 377
378 TString esdName("AliESDs.root");
379 if(gSystem->AccessPathName(esdName.Data(),kFileExists)) {
380 printf("AliD0toKpiAnalysis::FindCandidatesESD(): No ESDs file found!\n");
381 return;
382 }
383
384 TString outName1=outName;
385 TString outName2("nTotEvents.dat");
386
3a9a3487 387 Int_t nTotEv=0,nD0rec=0,nD0rec1ev=0;
388 Double_t dca;
389 Double_t v2[3],mom[6],d0[2];
390 Double_t alphaP,alphaN,ptP,ptN,phiP,phiN;
391 Int_t iTrkP,iTrkN,trkEntries;
392 Int_t nTrksP=0,nTrksN=0;
393 Int_t trkNum[2];
394 Double_t tofmass[2];
395 Int_t okD0=0,okD0bar=0;
396 AliITStrackV2 *postrack = 0;
397 AliITStrackV2 *negtrack = 0;
398
399 // create the AliITSVertexerTracks object
400 // (it will be used only if fVertexOnTheFly=kTrue)
401 AliITSVertexerTracks *vertexer1 = new AliITSVertexerTracks;
402 vertexer1->SetMinTracks(2);
3a9a3487 403 Int_t skipped[2];
404 Bool_t goodVtx1;
405
406
407 // define the cuts for vertexing
408 Double_t vtxcuts[]={50., // max. allowed chi2
409 0.0, // min. allowed negative daughter's impact param
410 0.0, // min. allowed positive daughter's impact param
411 1.0, // max. allowed DCA between the daughter tracks
412 -1.0, // min. allowed cosine of pointing angle
413 0.0, // min. radius of the fiducial volume
414 2.9};// max. radius of the fiducial volume
415
416 // create the AliV0vertexer object
417 AliV0vertexer *vertexer2 = new AliV0vertexer(vtxcuts);
418
419 // create tree for reconstructed D0s
420 AliD0toKpi *ioD0toKpi=0;
421 TTree *treeD0 = new TTree("TreeD0","Tree with D0 candidates");
422 treeD0->Branch("D0toKpi","AliD0toKpi",&ioD0toKpi,200000,0);
423
424 // open file with tracks
425 TFile *esdFile = TFile::Open(esdName.Data());
e6e6531b 426 AliESD* event = new AliESD;
427 TTree* tree = (TTree*) esdFile->Get("esdTree");
428 if (!tree) {
429 Error("FindCandidatesESD", "no ESD tree found");
430 return;
431 }
432 tree->SetBranchAddress("ESD", &event);
3a9a3487 433
3a9a3487 434 // loop on events in file
e6e6531b 435 for (Int_t iEvent = evFirst; iEvent < tree->GetEntries(); iEvent++) {
436 if (iEvent > evLast) break;
437 tree->GetEvent(iEvent);
3a9a3487 438 Int_t ev = (Int_t)event->GetEventNumber();
3a9a3487 439 printf("--- Finding D0 -> Kpi in event %d\n",ev);
440
441 // retrieve primary vertex from file
442 //sprintf(vtxName,"Vertex_%d",ev);
d681bb2d 443 //AliESDVertex *vertex1stored = (AliESDVertex*)trkFile->Get(vtxName);
3a9a3487 444 //vertex1stored->GetXYZ(fV1);
445 //delete vertex1stored;
2257f27e 446 event->GetVertex()->GetXYZ(fV1);
3a9a3487 447
448 trkEntries = event->GetNumberOfTracks();
449 printf(" Number of tracks: %d\n",trkEntries);
450
451 // count the total number of events
452 nTotEv++;
453
454 // call function which applies sigle-track selection and
455 // separetes positives and negatives
456 TObjArray trksP(trkEntries/2);
457 Int_t *trkEntryP = new Int_t[trkEntries];
458 TObjArray trksN(trkEntries/2);
459 Int_t *trkEntryN = new Int_t[trkEntries];
460 TTree *trkTree = new TTree();
461 if(fVertexOnTheFly) {
462 SelectTracksESDvtx(*event,trkTree,trksP,trkEntryP,nTrksP,
463 trksN,trkEntryN,nTrksN);
464 } else {
465 SelectTracksESD(*event,trksP,trkEntryP,nTrksP,
466 trksN,trkEntryN,nTrksN);
467 }
468
59c6807d 469 AliDebugClass(1,Form(" pos. tracks: %d neg .tracks: %d",nTrksP,nTrksN));
3a9a3487 470
471 nD0rec1ev = 0;
472
473 // loop on positive tracks
474 for(iTrkP=0; iTrkP<nTrksP; iTrkP++) {
59c6807d 475 if(iTrkP%10==0) AliDebugClass(1,Form(" Processing positive track number %d of %d",iTrkP,nTrksP));
3a9a3487 476
477 // get track from track array
478 postrack = (AliITStrackV2*)trksP.UncheckedAt(iTrkP);
479 trkNum[0] = trkEntryP[iTrkP];
480
481 // loop on negative tracks
482 for(iTrkN=0; iTrkN<nTrksN; iTrkN++) {
483
484 // get track from tracks array
485 negtrack = (AliITStrackV2*)trksN.UncheckedAt(iTrkN);
486 trkNum[1] = trkEntryN[iTrkN];
487
488 AliITStrackV2 nt(*negtrack), pt(*postrack), *pnt=&nt, *ppt=&pt;
489
490 //
491 // ----------- DCA MINIMIZATION ------------------
492 //
493 // find the DCA and propagate the tracks to the DCA
494 dca = vertexer2->PropagateToDCA(pnt,ppt);
495
496 // define the AliV0vertex object
497 AliV0vertex *vertex2 = new AliV0vertex(*pnt,*ppt);
498
499 // get position of the secondary vertex
500 vertex2->GetXYZ(v2[0],v2[1],v2[2]);
501
502 delete vertex2;
503
504 // momenta of the tracks at the vertex
505 ptP = 1./TMath::Abs(ppt->Get1Pt());
506 alphaP = ppt->GetAlpha();
507 phiP = alphaP+TMath::ASin(ppt->GetSnp());
508 mom[0] = ptP*TMath::Cos(phiP);
509 mom[1] = ptP*TMath::Sin(phiP);
510 mom[2] = ptP*ppt->GetTgl();
511
512 ptN = 1./TMath::Abs(pnt->Get1Pt());
513 alphaN = pnt->GetAlpha();
514 phiN = alphaN+TMath::ASin(pnt->GetSnp());
515 mom[3] = ptN*TMath::Cos(phiN);
516 mom[4] = ptN*TMath::Sin(phiN);
517 mom[5] = ptN*pnt->GetTgl();
518
519 goodVtx1 = kTRUE;
520 // no vertexing if DeltaMass > fMassCut
521 if(fVertexOnTheFly) {
522 goodVtx1 = kFALSE;
523 if(SelectInvMass(mom)) {
524 // primary vertex from *other* tracks in event
525 vertexer1->SetVtxStart(fV1[0],fV1[1]);
526 skipped[0] = trkEntryP[iTrkP];
527 skipped[1] = trkEntryN[iTrkN];
528 vertexer1->SetSkipTracks(2,skipped);
d681bb2d 529 AliESDVertex *vertex1onfly =
530 (AliESDVertex*)vertexer1->VertexOnTheFly(*trkTree);
3a9a3487 531 if(vertex1onfly->GetNContributors()>0) goodVtx1 = kTRUE;
532 vertex1onfly->GetXYZ(fV1);
533 //vertex1onfly->PrintStatus();
534 delete vertex1onfly;
535 }
536 }
537
538 // impact parameters of the tracks w.r.t. the primary vertex
539 d0[0] = 10000.*ppt->GetD(fV1[0],fV1[1]);
540 d0[1] = -10000.*pnt->GetD(fV1[0],fV1[1]);
541
542 // create the object AliD0toKpi
543 AliD0toKpi theD0(ev,trkNum,fV1,v2,dca,mom,d0);
544
545 // select D0s
546 if(goodVtx1 && theD0.Select(fD0Cuts,okD0,okD0bar)) {
547
548 // get PID info from ESD
549 AliESDtrack *t0 = (AliESDtrack*)event->GetTrack(trkNum[0]);
550 Double_t esdpid0[5];
551 t0->GetESDpid(esdpid0);
552 if(t0->GetStatus()&AliESDtrack::kTOFpid) {
553 tofmass[0] = CalculateTOFmass(t0->GetP(),
554 t0->GetIntegratedLength(),
555 t0->GetTOFsignal());
556 } else {
557 tofmass[0] = -1000.;
558 }
559 AliESDtrack *t1 = (AliESDtrack*)event->GetTrack(trkNum[1]);
560 Double_t esdpid1[5];
561 t1->GetESDpid(esdpid1);
562 if(t1->GetStatus()&AliESDtrack::kTOFpid) {
563 tofmass[1] = CalculateTOFmass(t1->GetP(),
564 t1->GetIntegratedLength(),
565 t1->GetTOFsignal());
566 } else {
567 tofmass[1] = -1000.;
568 }
569
570 theD0.SetPIDresponse(esdpid0,esdpid1);
571 theD0.SetTOFmasses(tofmass);
572
573 // fill the tree
574 ioD0toKpi=&theD0;
575 treeD0->Fill();
576
577 nD0rec++; nD0rec1ev++;
578 ioD0toKpi=0;
579 }
580
581 negtrack = 0;
582 } // loop on negative tracks
583 postrack = 0;
584 } // loop on positive tracks
585
586 trksP.Delete();
587 trksN.Delete();
588 delete [] trkEntryP;
589 delete [] trkEntryN;
590 delete trkTree;
3a9a3487 591
592
593 printf(" Number of D0 candidates: %d\n",nD0rec1ev);
594 } // loop on events in file
595
596
597 printf("\n+++\n+++ Total number of events: %d\n+++\n",nTotEv);
598 printf("\n+++\n+++ Total number of D0 candidates: %d\n+++\n",nD0rec);
599
600 delete vertexer1;
601 delete vertexer2;
602
603 esdFile->Close();
604
605 // create a copy of this class to be written to output file
606 AliD0toKpiAnalysis *copy = (AliD0toKpiAnalysis*)this->Clone("D0toKpiAnalysis");
607
608 // add PDG codes to decay tracks in found candidates (in simulation mode)
609 // and store tree in the output file
610 if(!fSim) {
611 TFile *outroot = new TFile(outName1.Data(),"recreate");
612 treeD0->Write();
613 copy->Write();
614 outroot->Close();
615 delete outroot;
616 } else {
617 printf(" Now adding information from simulation (PDG codes) ...\n");
618 TTree *treeD0sim = new TTree("TreeD0","Tree with D0 candidates");
619 MakeTracksRefFileESD();
620 SimulationInfo(treeD0,treeD0sim);
621 delete treeD0;
622 TFile *outroot = new TFile(outName1.Data(),"recreate");
623 treeD0sim->Write();
624 copy->Write();
625 outroot->Close();
626 delete outroot;
627 }
628
629 // write to a file the total number of events
630 FILE *outdat = fopen(outName2.Data(),"w");
631 fprintf(outdat,"%d\n",nTotEv);
632 fclose(outdat);
633
634 return;
635}
636//-----------------------------------------------------------------------------
637void AliD0toKpiAnalysis::PrintStatus() const {
638 // Print parameters being used
639
640 printf(" fBz = %f T\n",fBz);
641 printf("Preselections:\n");
642 printf(" fPtCut = %f GeV\n",fPtCut);
643 printf(" fd0Cut = %f micron\n",fd0Cut);
644 printf(" fMassCut = %f GeV\n",fMassCut);
645 if(fVertexOnTheFly) printf("Primary vertex on the fly\n");
646 if(fSim) {
647 printf("Simulation mode\n");
648 if(fOnlySignal) printf(" Only signal goes to file\n");
649 }
650 printf("Cuts on candidates:\n");
651 printf(" |M-MD0| [GeV] < %f\n",fD0Cuts[0]);
652 printf(" dca [micron] < %f\n",fD0Cuts[1]);
653 printf(" cosThetaStar < %f\n",fD0Cuts[2]);
654 printf(" pTK [GeV] > %f\n",fD0Cuts[3]);
655 printf(" pTpi [GeV] > %f\n",fD0Cuts[4]);
656 printf(" |d0K| [micron] < %f\n",fD0Cuts[5]);
657 printf(" |d0pi| [micron] < %f\n",fD0Cuts[6]);
658 printf(" d0d0 [micron^2] < %f\n",fD0Cuts[7]);
659 printf(" cosThetaPoint > %f\n",fD0Cuts[8]);
660
661 return;
662}
663//-----------------------------------------------------------------------------
664Bool_t AliD0toKpiAnalysis::SelectInvMass(const Double_t p[6]) const {
665 // Apply preselection in the invariant mass of the pair
666
667 Double_t mD0 = 1.8645;
668 Double_t mPi = 0.13957;
669 Double_t mKa = 0.49368;
670
671 Double_t energy[2];
672 Double_t mom2[2],momTot2;
673
674 mom2[0] = p[0]*p[0] + p[1]*p[1] + p[2]*p[2];
675 mom2[1] = p[3]*p[3] + p[4]*p[4] + p[5]*p[5];
676
677 momTot2 = (p[0]+p[3])*(p[0]+p[3])+
678 (p[1]+p[4])*(p[1]+p[4])+
679 (p[2]+p[5])*(p[2]+p[5]);
680
681 // D0 -> K- pi+
682 energy[1] = TMath::Sqrt(mKa*mKa+mom2[1]);
683 energy[0] = TMath::Sqrt(mPi*mPi+mom2[0]);
684
685 Double_t minvD0 = TMath::Sqrt((energy[0]+energy[1])*(energy[0]+energy[1])-momTot2);
686
687 // D0bar -> K+ pi-
688 energy[0] = TMath::Sqrt(mKa*mKa+mom2[0]);
689 energy[1] = TMath::Sqrt(mPi*mPi+mom2[1]);
690
691 Double_t minvD0bar = TMath::Sqrt((energy[0]+energy[1])*(energy[0]+energy[1])-momTot2);
692
693 if(TMath::Abs(minvD0-mD0) < fMassCut) return kTRUE;
694 if(TMath::Abs(minvD0bar-mD0) < fMassCut) return kTRUE;
695 return kFALSE;
696}
697//-----------------------------------------------------------------------------
698void AliD0toKpiAnalysis::SelectTracks(TTree &trkTree,
699 TObjArray &trksP,Int_t *trkEntryP,Int_t &nTrksP,
700 TObjArray &trksN,Int_t *trkEntryN,Int_t &nTrksN) const {
701 // Create two TObjArrays with positive and negative tracks and
702 // apply single-track preselection
703
704 nTrksP=0,nTrksN=0;
705
706 Int_t entr = (Int_t)trkTree.GetEntries();
707
708 // trasfer tracks from tree to arrays
709 for(Int_t i=0; i<entr; i++) {
710
711 AliITStrackV2 *track = new AliITStrackV2;
712 trkTree.SetBranchAddress("tracks",&track);
713
714 trkTree.GetEvent(i);
715
716 // single track selection
717 if(!SingleTrkCuts(*track)) { delete track; continue; }
718
719 if(track->Get1Pt()>0.) { // negative track
720 trksN.AddLast(track);
721 trkEntryN[nTrksN] = i;
722 nTrksN++;
723 } else { // positive track
724 trksP.AddLast(track);
725 trkEntryP[nTrksP] = i;
726 nTrksP++;
727 }
728
729 }
730
731 return;
732}
733//-----------------------------------------------------------------------------
734void AliD0toKpiAnalysis::SelectTracksESD(AliESD &event,
735 TObjArray &trksP,Int_t *trkEntryP,Int_t &nTrksP,
736 TObjArray &trksN,Int_t *trkEntryN,Int_t &nTrksN) const {
737 // Create two TObjArrays with positive and negative tracks and
738 // apply single-track preselection
739
740 nTrksP=0,nTrksN=0;
741
742 Int_t entr = event.GetNumberOfTracks();
743
744 // transfer ITS tracks from ESD to arrays and to a tree
745 for(Int_t i=0; i<entr; i++) {
3a9a3487 746
747 AliESDtrack *esdtrack = (AliESDtrack*)event.GetTrack(i);
748 UInt_t status = esdtrack->GetStatus();
749
750 if(!(status&AliESDtrack::kITSrefit)) continue;
751
752 AliITStrackV2 *itstrack = new AliITStrackV2(*esdtrack);
753
754 // single track selection
755 if(!SingleTrkCuts(*itstrack))
756 { delete itstrack; continue; }
757
758 if(itstrack->Get1Pt()>0.) { // negative track
759 trksN.AddLast(itstrack);
760 trkEntryN[nTrksN] = i;
761 nTrksN++;
762 } else { // positive track
763 trksP.AddLast(itstrack);
764 trkEntryP[nTrksP] = i;
765 nTrksP++;
766 }
767
768 } // loop on ESD tracks
769
770 return;
771}
772//-----------------------------------------------------------------------------
773void AliD0toKpiAnalysis::SelectTracksESDvtx(AliESD &event,TTree *trkTree,
774 TObjArray &trksP,Int_t *trkEntryP,Int_t &nTrksP,
775 TObjArray &trksN,Int_t *trkEntryN,Int_t &nTrksN) const {
776 // Create two TObjArrays with positive and negative tracks and
777 // apply single-track preselection
778
779 nTrksP=0,nTrksN=0;
780
781 Int_t entr = event.GetNumberOfTracks();
782
783 AliITStrackV2 *itstrackfortree = 0;
784 trkTree->Branch("tracks","AliITStrackV2",&itstrackfortree,entr,0);
785
786
787 // transfer ITS tracks from ESD to arrays and to a tree
788 for(Int_t i=0; i<entr; i++) {
789
790 AliESDtrack *esdtrack = (AliESDtrack*)event.GetTrack(i);
791 UInt_t status = esdtrack->GetStatus();
792
793 if(!(status&AliESDtrack::kITSrefit)) continue;
794
795 AliITStrackV2 *itstrack = new AliITStrackV2(*esdtrack);
796
797 // store track in the tree to be used for primary vertex finding
798 itstrackfortree = new AliITStrackV2(*esdtrack);
799 trkTree->Fill();
800 itstrackfortree = 0;
801
802 // single track selection
803 if(!SingleTrkCuts(*itstrack))
804 { delete itstrack; continue; }
805
806 if(itstrack->Get1Pt()>0.) { // negative track
807 trksN.AddLast(itstrack);
808 trkEntryN[nTrksN] = i;
809 nTrksN++;
810 } else { // positive track
811 trksP.AddLast(itstrack);
812 trkEntryP[nTrksP] = i;
813 nTrksP++;
814 }
815
816 } // loop on esd tracks
817
818 delete itstrackfortree;
819
820 return;
821}
822//-----------------------------------------------------------------------------
823void AliD0toKpiAnalysis::SetD0Cuts(Double_t cut0,Double_t cut1,
824 Double_t cut2,Double_t cut3,Double_t cut4,
825 Double_t cut5,Double_t cut6,
826 Double_t cut7,Double_t cut8) {
827 // Set the cuts for D0 selection
828 fD0Cuts[0] = cut0;
829 fD0Cuts[1] = cut1;
830 fD0Cuts[2] = cut2;
831 fD0Cuts[3] = cut3;
832 fD0Cuts[4] = cut4;
833 fD0Cuts[5] = cut5;
834 fD0Cuts[6] = cut6;
835 fD0Cuts[7] = cut7;
836 fD0Cuts[8] = cut8;
837
838 return;
839}
840//-----------------------------------------------------------------------------
841void AliD0toKpiAnalysis::SetD0Cuts(const Double_t cuts[9]) {
842 // Set the cuts for D0 selection
843
844 for(Int_t i=0; i<9; i++) fD0Cuts[i] = cuts[i];
845
846 return;
847}
848//-----------------------------------------------------------------------------
849Bool_t AliD0toKpiAnalysis::SingleTrkCuts(const AliITStrackV2& trk) const {
850 // Check if track passes some kinematical cuts
851
852 if(TMath::Abs(1./trk.Get1Pt()) < fPtCut)
853 return kFALSE;
854 if(TMath::Abs(10000.*trk.GetD(fV1[0],fV1[1])) < fd0Cut)
855 return kFALSE;
856
857 return kTRUE;
858}
859//----------------------------------------------------------------------------
860void AliD0toKpiAnalysis::MakeTracksRefFile(Int_t evFirst,Int_t evLast)
861 const {
862 // Create a file with simulation info for the reconstructed tracks
863
864 TFile *out = TFile::Open("ITStracksRefFile.root","recreate");
865 TFile *trk = TFile::Open("AliITStracksV2.root");
866 AliRunLoader *rl = AliRunLoader::Open("galice.root");
867
868 // load kinematics
869 rl->LoadKinematics();
870
871 Int_t label;
872 TParticle *part;
873 TParticle *mumpart;
874 REFTRACK reftrk;
875
876 for(Int_t ev=evFirst; ev<=evLast; ev++){
877 rl->GetEvent(ev);
878 AliStack *stack = rl->Stack();
879
880 trk->cd();
881
882 // Tree with tracks
883 char tname[100];
884 sprintf(tname,"TreeT_ITS_%d",ev);
885
886 TTree *tracktree=(TTree*)trk->Get(tname);
887 if(!tracktree) continue;
888 AliITStrackV2 *track = new AliITStrackV2;
889 tracktree->SetBranchAddress("tracks",&track);
890 Int_t nentr=(Int_t)tracktree->GetEntries();
891
892 // Tree for true track parameters
893 char ttname[100];
894 sprintf(ttname,"Tree_Ref_%d",ev);
895 TTree *reftree = new TTree(ttname,"Tree with true track params");
896 reftree->Branch("rectracks",&reftrk,"lab/I:pdg:mumlab:mumpdg:Vx/F:Vy:Vz:Px:Py:Pz");
897
898 for(Int_t i=0; i<nentr; i++) {
899 tracktree->GetEvent(i);
900 label = TMath::Abs(track->GetLabel());
901
902 part = (TParticle*)stack->Particle(label);
903 reftrk.lab = label;
904 reftrk.pdg = part->GetPdgCode();
905 reftrk.mumlab = part->GetFirstMother();
906 if(part->GetFirstMother()>=0) {
907 mumpart = (TParticle*)stack->Particle(part->GetFirstMother());
908 reftrk.mumpdg = mumpart->GetPdgCode();
909 } else {
910 reftrk.mumpdg=-1;
911 }
912 reftrk.Vx = part->Vx();
913 reftrk.Vy = part->Vy();
914 reftrk.Vz = part->Vz();
915 reftrk.Px = part->Px();
916 reftrk.Py = part->Py();
917 reftrk.Pz = part->Pz();
918
919 reftree->Fill();
920 } // loop on tracks
921
922 out->cd();
923 reftree->Write();
924
925 delete track;
926 delete reftree;
927 delete tracktree;
928 delete stack;
929 } // loop on events
930
931 trk->Close();
932 out->Close();
933
934 return;
935}
936//----------------------------------------------------------------------------
937void AliD0toKpiAnalysis::MakeTracksRefFileESD() const {
938 // Create a file with simulation info for the reconstructed tracks
939
940 TFile *outFile = TFile::Open("ITStracksRefFile.root","recreate");
941 TFile *esdFile = TFile::Open("AliESDs.root");
942 AliRunLoader *rl = AliRunLoader::Open("galice.root");
943
944 // load kinematics
945 rl->LoadKinematics();
946
947 Int_t label;
948 TParticle *part;
949 TParticle *mumpart;
950 REFTRACK reftrk;
951
952 TKey *key=0;
953 TIter next(esdFile->GetListOfKeys());
954 // loop on events in file
955 while ((key=(TKey*)next())!=0) {
956 AliESD *event=(AliESD*)key->ReadObj();
957 Int_t ev = (Int_t)event->GetEventNumber();
958
959 rl->GetEvent(ev);
960 AliStack *stack = rl->Stack();
961
962 Int_t nentr=(Int_t)event->GetNumberOfTracks();
963
964 // Tree for true track parameters
965 char ttname[100];
966 sprintf(ttname,"Tree_Ref_%d",ev);
967 TTree *reftree = new TTree(ttname,"Tree with true track params");
968 reftree->Branch("rectracks",&reftrk,"lab/I:pdg:mumlab:mumpdg:Vx/F:Vy:Vz:Px:Py:Pz");
969
970 for(Int_t i=0; i<nentr; i++) {
971 AliESDtrack *esdtrack = (AliESDtrack*)event->GetTrack(i);
972 label = TMath::Abs(esdtrack->GetLabel());
973
974 part = (TParticle*)stack->Particle(label);
975 reftrk.lab = label;
976 reftrk.pdg = part->GetPdgCode();
977 reftrk.mumlab = part->GetFirstMother();
978 if(part->GetFirstMother()>=0) {
979 mumpart = (TParticle*)stack->Particle(part->GetFirstMother());
980 reftrk.mumpdg = mumpart->GetPdgCode();
981 } else {
982 reftrk.mumpdg=-1;
983 }
984 reftrk.Vx = part->Vx();
985 reftrk.Vy = part->Vy();
986 reftrk.Vz = part->Vz();
987 reftrk.Px = part->Px();
988 reftrk.Py = part->Py();
989 reftrk.Pz = part->Pz();
990
991 reftree->Fill();
992
993 } // loop on tracks
994
995 outFile->cd();
996 reftree->Write();
997
998 delete reftree;
999 delete event;
1000 delete stack;
1001 } // loop on events
1002
1003 esdFile->Close();
1004 outFile->Close();
1005
1006 return;
1007}
1008//-----------------------------------------------------------------------------
1009void AliD0toKpiAnalysis::SimulationInfo(TTree *treeD0in,TTree *treeD0out) const {
1010 // add pdg codes to candidate decay tracks (for sim)
1011
1012 TString refFileName("ITStracksRefFile.root");
1013 if(fSim && gSystem->AccessPathName(refFileName.Data(),kFileExists)) {
1014 printf("AliD0toKpiAnalysis::SimulationInfo: no reference file found!\n");
1015 return;
1016 }
1017 TFile *refFile = TFile::Open(refFileName.Data());
1018
1019 Char_t refTreeName[100];
1020 Int_t event;
1021 Int_t pdg[2],mumpdg[2],mumlab[2];
1022 REFTRACK reftrk;
1023
1024 // read-in reference tree for event 0 (the only event for Pb-Pb)
1025 sprintf(refTreeName,"Tree_Ref_%d",0);
1026 TTree *refTree0 = (TTree*)refFile->Get(refTreeName);
1027 refTree0->SetBranchAddress("rectracks",&reftrk);
1028
1029 AliD0toKpi *theD0 = 0;
1030 treeD0in->SetBranchAddress("D0toKpi",&theD0);
1031 treeD0out->Branch("D0toKpi","AliD0toKpi",&theD0,200000,0);
1032
1033 Int_t entries = (Int_t)treeD0in->GetEntries();
1034
1035 for(Int_t i=0; i<entries; i++) {
1036 if(i%100==0) printf(" done %d candidates of %d\n",i,entries);
1037
1038 treeD0in->GetEvent(i);
1039 event = theD0->EventNo();
1040
1041 if(event==0) { // always true for Pb-Pb (avoid to read-in tree every time)
1042 refTree0->GetEvent(theD0->GetTrkNum(0));
1043 pdg[0] = reftrk.pdg;
1044 mumpdg[0] = reftrk.mumpdg;
1045 mumlab[0] = reftrk.mumlab;
1046 refTree0->GetEvent(theD0->GetTrkNum(1));
1047 pdg[1] = reftrk.pdg;
1048 mumpdg[1] = reftrk.mumpdg;
1049 mumlab[1] = reftrk.mumlab;
1050 } else {
1051 sprintf(refTreeName,"Tree_Ref_%d",event);
1052 TTree *refTree = (TTree*)refFile->Get(refTreeName);
1053 refTree->SetBranchAddress("rectracks",&reftrk);
1054 refTree->GetEvent(theD0->GetTrkNum(0));
1055 pdg[0] = reftrk.pdg;
1056 mumpdg[0] = reftrk.mumpdg;
1057 mumlab[0] = reftrk.mumlab;
1058 refTree->GetEvent(theD0->GetTrkNum(1));
1059 pdg[1] = reftrk.pdg;
1060 mumpdg[1] = reftrk.mumpdg;
1061 mumlab[1] = reftrk.mumlab;
1062 delete refTree;
1063 }
1064
1065 theD0->SetPdgCodes(pdg);
1066 theD0->SetMumPdgCodes(mumpdg);
1067
1068 if(TMath::Abs(mumpdg[0])==421 && TMath::Abs(mumpdg[1])==421
1069 && mumlab[0]==mumlab[1]) theD0->SetSignal();
1070
1071 if(!fOnlySignal || theD0->IsSignal()) treeD0out->Fill();
1072
1073 }
1074
1075 delete refTree0;
1076
1077 refFile->Close();
1078
1079 return;
1080}
1081
1082
1083
1084
1085
1086