Merged HLT tag v1-2 with ALIROOT tag v3-09-Release.
[u/mrichter/AliRoot.git] / HLT / trigger / AliD0vtxFinderBkg_pp_VTX.C
CommitLineData
3e87ef69 1/****************************************************************************
2 * *
3 * This macro computes the position of the D0 decay vertex using *
4 * helix DCA minimization by J.Belikov *
5 * *
6 * Reconstructed D0 are stored in a tree as AliD0toKpi objects *
7 * *
8 ****************************************************************************/
9//#define __COMPILE__
10//#ifdef __COMPILE__
11#if !defined(__CINT__) || defined(__MAKECINT__)
12//-- --- standard headers-------------
13#include <iostream.h>
14//--------Root headers ---------------
15#include <TSystem.h>
16#include <TFile.h>
17#include <TString.h>
18#include <TStopwatch.h>
19#include <TObject.h>
20#include <TVector3.h>
21#include <TTree.h>
22#include <TParticle.h>
23#include <TArray.h>
24//----- AliRoot headers ---------------
25#include "alles.h"
26#include "AliRun.h"
27#include "AliKalmanTrack.h"
28#include "AliITStrackV2.h"
29#include "AliHeader.h"
30#include "AliGenEventHeader.h"
31#include "AliV0vertex.h"
32#include "AliV0vertexer.h"
33#include "AliITSVertex.h"
34#include "AliITSVertexer.h"
35#include "AliITSVertexerTracks.h"
36#include "AliD0toKpi.h"
37#endif
38//-------------------------------------
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} RECTRACK;
48
49// field (T)
50const Double_t kBz = 0.4;
51
52// primary vertex
53Double_t gv1[3] = {0.,0.,0,};
54
55// sigle track cuts
56//const Double_t kPtCut = 0.5; // GeV/c
57//const Double_t kd0Cut = 50.; // micron
58const Double_t kPtCut = 0.; // GeV/c
59const Double_t kd0Cut = 0.; // micron
60
61// cuts on D0 candidate (to be passed to function AliD0toKpi::Select())
62//
63// cuts[0] = inv. mass half width [GeV]
64// cuts[1] = dca [micron]
65// cuts[2] = cosThetaStar
66// cuts[3] = pTK [GeV/c]
67// cuts[4] = pTPi [GeV/c]
68// cuts[5] = d0K [micron] upper limit!
69// cuts[6] = d0Pi [micron] upper limit!
70// cuts[7] = d0d0 [micron^2]
71// cuts[8] = cosThetaPoint
72const Double_t kD0cuts[9] = {0.012,
73 500.,
74 0.8,
75 0.5,
76 0.5,
77 500.,
78 500.,
79 -5000.,
80 0.8};
81const Double_t kD0cutsAll[9] = {.1,
82 1.e5,
83 1.1,
84 0.,
85 0.,
86 1.e10,
87 1.e10,
88 1.e10,
89 -1.1};
90// mass cut for 1st level rejection
91const Double_t kMassCut = 0.1; // GeV
92
93// this function ckecks files existence
94Bool_t GetInputFiles();
95
96// 1st level rejection based on inv. mass
97Bool_t SelectInvMass(const Double_t p[6]);
98
99// this function applies single track cuts
100Bool_t TrkCuts(const AliITStrackV2& trk);
101
102// this function creates TObjArrays with positive and negative tracks
103void SelectTracks(TTree& itsTree,
104 TObjArray& trksP,Int_t* itsEntryP,Int_t& nTrksP,
105 TObjArray& trksN,Int_t* itsEntryN,Int_t& nTrksN);
106
107
108void AliD0vtxFinderBkg_pp_VTX(Int_t evFirst=0,Int_t evLast=4999) {
109
110 const Char_t *name="AliD0vtxFinderBkg_pp_VTX";
111 cerr<<'\n'<<name<<"...\n";
112 gBenchmark->Start(name);
113
114 AliKalmanTrack::SetConvConst(100/0.299792458/kBz);
115
116
117 // load library for D0 candidates class
118 gSystem->Load("libD0toKpi.so");
119
120 // check existence of input files
121 if(!GetInputFiles()) { cerr<<"No tracks file found"<<endl; return; }
122
123 TString outName("AliD0toKpiBkg.root");
124 TString outName2("nTotEv.dat");
125
126 Int_t ev;
127 Int_t nTotEv=0,nD0rec=0,nD0rec1ev=0;
128 Double_t dca;
129 Double_t v2[3],mom[6],d0[2];
130 Int_t pdg[2],mum[2];
131 Double_t alphaP,alphaN,ptP,ptN,phiP,phiN;
132 Int_t iTrkP,iTrkN,itsEntries;
133 Int_t nTrksP=0,nTrksN=0;
134 Int_t okD0=0,okD0bar=0;
135 Char_t trksName[100],refsName[100],vtxName[100];
136 RECTRACK rectrk;
137 AliITStrackV2 *postrack = 0;
138 AliITStrackV2 *negtrack = 0;
139
140 // create the AliITSVertexerTracks object
141 AliITSVertexerTracks *vertexer1 = new AliITSVertexerTracks;
142 //vertexer1->SetUseThrustFrame(1);
143 vertexer1->SetMinTracks(2);
144 vertexer1->SetDebug(0);
145 Int_t skipped[2];
146
147 // define the cuts for vertexing
148 Double_t vtxcuts[]={33., // max. allowed chi2
149 0.0, // min. allowed negative daughter's impact param
150 0.0, // min. allowed positive daughter's impact param
151 1.0, // max. allowed DCA between the daughter tracks
152 -1.0, // min. allowed cosine of V0's pointing angle
153 0.0, // min. radius of the fiducial volume
154 2.9};// max. radius of the fiducial volume
155
156 // create the AliV0vertexer object
157 AliV0vertexer *vertexer2 = new AliV0vertexer(vtxcuts);
158
159
160 // Create tree for reconstructed D0s
161 TTree D0Tree("TreeD0","Tree with D0 candidates");
162 AliD0toKpi *ioD0toKpi=0;
163 D0Tree.Branch("D0toKpi","AliD0toKpi",&ioD0toKpi,200000,0);
164
165 // Open file with ITS tracks
166 TFile* itstrks = TFile::Open("AliITStracksV2.root");
167
168 // Open file with ITS track references
169 TFile* itsrefs = TFile::Open("ITStracksRefFile.root");
170
171 // loop on events in file
172 for(ev=evFirst; ev<=evLast; ev++) {
173 printf(" --- Processing event %d ---\n",ev);
174 sprintf(vtxName,"VertexTracks_%d",ev);
175 sprintf(trksName,"TreeT_ITS_%d",ev);
176 sprintf(refsName,"Tree_Ref_%d",ev);
177
178
179 // tracks from ITS
180 TTree *itsTree=(TTree*)itstrks->Get(trksName);
181 if(!itsTree) continue;
182 itsEntries = (Int_t)itsTree->GetEntries();
183 printf("+++\n+++ Number of tracks in ITS: %d\n+++\n\n",itsEntries);
184
185 // tree from reference file
186 TTree *refTree=(TTree*)itsrefs->Get(refsName);
187 refTree->SetBranchAddress("rectracks",&rectrk);
188
189 // count the total number of events
190 nTotEv++;
191
192
193 // call function which applies sigle track selection and
194 // separetes positives and negatives
195 TObjArray trksP(itsEntries/2);
196 Int_t *itsEntryP = new Int_t[itsEntries];
197 TObjArray trksN(itsEntries/2);
198 Int_t *itsEntryN = new Int_t[itsEntries];
199 SelectTracks(*itsTree,trksP,itsEntryP,nTrksP,trksN,itsEntryN,nTrksN);
200
201 nD0rec1ev = 0;
202 // loop on positive tracks
203 for(iTrkP=0; iTrkP<nTrksP; iTrkP++) {
204 //if(iTrkP % 10==0) cerr<<"--- Processing positive track number: "<<iTrkP<<" of "<<nTrksP<<"\r";
205
206 // get track from track array
207 postrack = (AliITStrackV2*)trksP.At(iTrkP);
208
209 // get info on tracks PDG and mothers PDG from reference file
210 refTree->GetEvent(itsEntryP[iTrkP]);
211 pdg[0] = rectrk.pdg;
212 mum[0] = rectrk.mumpdg;
213
214 // loop on negative tracks
215 for(iTrkN=0; iTrkN<nTrksN; iTrkN++) {
216 // get track from track array
217 negtrack = (AliITStrackV2*)trksN.At(iTrkN);
218
219 // get info on tracks PDG and mothers PDG from reference file
220 refTree->GetEvent(itsEntryN[iTrkN]);
221 pdg[1] = rectrk.pdg;
222 mum[1] = rectrk.mumpdg;
223
224 //
225 // ----------- DCA MINIMIZATION ------------------
226 //
227 // find the DCA and propagate the tracks to the DCA
228 dca = vertexer2->PropagateToDCA(negtrack,postrack);
229
230 // define the AliV0vertex object
231 AliV0vertex *vertex2 = new AliV0vertex(*negtrack,*postrack);
232
233 // get position of the vertex
234 vertex2->GetXYZ(v2[0],v2[1],v2[2]);
235
236 delete vertex2;
237
238 // momenta of the tracks at the vertex
239 ptP = 1./TMath::Abs(postrack->Get1Pt());
240 alphaP = postrack->GetAlpha();
241 phiP = alphaP+TMath::ASin(postrack->GetSnp());
242 mom[0] = ptP*TMath::Cos(phiP);
243 mom[1] = ptP*TMath::Sin(phiP);
244 mom[2] = ptP*postrack->GetTgl();
245
246 ptN = 1./TMath::Abs(negtrack->Get1Pt());
247 alphaN = negtrack->GetAlpha();
248 phiN = alphaN+TMath::ASin(negtrack->GetSnp());
249 mom[3] = ptN*TMath::Cos(phiN);
250 mom[4] = ptN*TMath::Sin(phiN);
251 mom[5] = ptN*negtrack->GetTgl();
252
253
254 Bool_t goodVtx = kFALSE;
255 // no vertexing if DeltaMass > kMassCut
256 if(SelectInvMass(mom)) {
257 // primary vertex from other tracks in event
258 vertexer1->SetVtxStart(0.,0.);
259 skipped[0] = itsEntryP[iTrkP];
260 skipped[1] = itsEntryN[iTrkN];
261 vertexer1->SetSkipTracks(2,skipped);
262 AliITSVertex *vertex1 =
263 (AliITSVertex*)vertexer1->VertexOnTheFly(*itsTree);
264 if(vertex1->GetNContributors()>0) goodVtx = kTRUE;
265 vertex1->GetXYZ(gv1);
266 // impact parameters of the tracks w.r.t. the smeared primary vertex
267 d0[0] = 10000.*postrack->GetD(gv1[0],gv1[1]);
268 d0[1] = -10000.*negtrack->GetD(gv1[0],gv1[1]);
269 //vertex1->PrintStatus();
270 delete vertex1;
271 }
272
273 // create the object TD0rec and store it in the tree
274 AliD0toKpi theD0(kFALSE,ev,gv1,v2,dca,mom,d0,pdg,mum);
275
276
277 // select D0s
278 if(goodVtx && theD0.Select(kD0cutsAll,okD0,okD0bar)) {
279
280 // compute the weights
281 theD0.ComputeWgts();
282
283 // fill the tree
284 ioD0toKpi=&theD0;
285 D0Tree.Fill();
286
287 nD0rec++; nD0rec1ev++;
288 ioD0toKpi=0;
289 }
290
291 negtrack = 0;
292 } // loop on negative tracks
293 postrack = 0;
294 } // loop on positive tracks
295
296 delete [] itsEntryP;
297 delete [] itsEntryN;
298 delete itsTree;
299 delete refTree;
300
301 printf("\n+++\n+++ Number of D0 candidates: %d\n+++\n",nD0rec1ev);
302 printf("----------------------------------------------------------\n");
303
304 } // loop on events in file
305
306
307 printf("\n+++\n+++ Total number of events: %d\n+++\n",nTotEv);
308 printf("\n+++\n+++ Total number of D0 candidates: %d\n+++\n",nD0rec);
309
310 delete vertexer1;
311 delete vertexer2;
312
313 itstrks->Close();
314 itsrefs->Close();
315
316 // store tree in a file
317 TFile* outroot = new TFile(outName.Data(),"recreate");
318 D0Tree.Write();
319 outroot->Close();
320 delete outroot;
321
322 // write to a file the total number of events
323 FILE* outdat = fopen(outName2.Data(),"w");
324 fprintf(outdat,"%d\n",nTotEv);
325 fclose(outdat);
326
327
328 gBenchmark->Stop(name);
329 gBenchmark->Show(name);
330
331 return;
332}
333//___________________________________________________________________________
334Bool_t GetInputFiles() {
335 TString itsName("AliITStracksV2.root");
336 TString refName("ITStracksRefFile.root");
337
338 if(gSystem->AccessPathName(itsName.Data(),kFileExists) ||
339 gSystem->AccessPathName(refName.Data(),kFileExists)) return kFALSE;
340
341 return kTRUE;
342}
343//___________________________________________________________________________
344Bool_t SelectInvMass(const Double_t p[6]) {
345
346 Double_t mD0 = 1.8645;
347 Double_t mPi = 0.13957;
348 Double_t mKa = 0.49368;
349
350 Double_t energy[2];
351 Double_t mom2[2],momTot2;
352
353 mom2[0] = p[0]*p[0]+p[1]*p[1]+p[2]*p[2];
354 mom2[1] = p[3]*p[3]+p[4]*p[4]+p[5]*p[5];
355
356 momTot2 = (p[0]+p[3])*(p[0]+p[3])+
357 (p[1]+p[4])*(p[1]+p[4])+
358 (p[2]+p[5])*(p[2]+p[5]);
359
360 // D0 -> K- Pi+
361 energy[1] = TMath::Sqrt(mKa*mKa+mom2[1]);
362 energy[0] = TMath::Sqrt(mPi*mPi+mom2[0]);
363
364 Double_t minvD0 = TMath::Sqrt((energy[0]+energy[1])*(energy[0]+energy[1])-momTot2);
365
366
367 // D0bar -> K+ Pi-
368 energy[0] = TMath::Sqrt(mKa*mKa+mom2[0]);
369 energy[1] = TMath::Sqrt(mPi*mPi+mom2[1]);
370
371 Double_t minvD0bar = TMath::Sqrt((energy[0]+energy[1])*(energy[0]+energy[1])-momTot2);
372
373 if(TMath::Abs(minvD0-mD0) < kMassCut) return kTRUE;
374 if(TMath::Abs(minvD0bar-mD0) < kMassCut) return kTRUE;
375 return kFALSE;
376}
377//___________________________________________________________________________
378void SelectTracks(TTree& itsTree,
379 TObjArray& trksP,Int_t* itsEntryP,Int_t& nTrksP,
380 TObjArray& trksN,Int_t* itsEntryN,Int_t& nTrksN) {
381//
382// this function creates two TObjArrays with positive and negative tracks
383//
384 nTrksP=0,nTrksN=0;
385
386
387 Int_t entr = (Int_t)itsTree.GetEntries();
388
389 // trasfer tracks from tree to arrays
390 for(Int_t i=0; i<entr; i++) {
391
392 AliITStrackV2 *itstrack = new AliITStrackV2;
393 itsTree.SetBranchAddress("tracks",&itstrack);
394
395 itsTree.GetEvent(i);
396
397 // single track selection
398 if(!TrkCuts(*itstrack)) { delete itstrack; continue; }
399
400 if(itstrack->Get1Pt()>0.) { // negative track
401 trksN.AddLast(itstrack);
402 itsEntryN[nTrksN] = i;
403 nTrksN++;
404 } else { // positive track
405 trksP.AddLast(itstrack);
406 itsEntryP[nTrksP] = i;
407 nTrksP++;
408 }
409
410 }
411
412 return;
413}
414//____________________________________________________________________________
415Bool_t TrkCuts(const AliITStrackV2& trk) {
416//
417// this function tells if track passes some kinematical cuts
418//
419 if(TMath::Abs(1./trk.Get1Pt()) < kPtCut) return kFALSE;
420 if(TMath::Abs(10000.*trk.GetD(gv1[0],gv1[1])) < kd0Cut) return kFALSE;
421
422 return kTRUE;
423}
424
425
426
427
428
429
430
431