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