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