Merged HLT tag v1-2 with ALIROOT tag v3-09-Release.
[u/mrichter/AliRoot.git] / HLT / trigger / AliD0vtxFinderSgn_pp_VTX.C
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
40 typedef 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)
50 const Double_t kBz = 0.4;
51
52 // primary vertex
53 Double_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
58 const Double_t kPtCut = 0.;  // GeV/c
59 const 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
72 const 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};
81 const 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
91 const Double_t kMassCut = 0.1; // GeV
92
93 // this function ckecks files existence
94 Bool_t GetInputFiles();
95
96 // 1st level rejection based on inv. mass
97 Bool_t SelectInvMass(const Double_t p[6]);
98
99 // this function creates TObjArrays with positive and negative tracks
100 void   SelectTracks(TTree& itsTree,
101                     TObjArray& trksP,Int_t* itsEntryP,Int_t& nTrksP,
102                     TObjArray& trksN,Int_t* itsEntryN,Int_t& nTrksN);
103
104 // this function applies single track cuts
105 Bool_t TrkCuts(const AliITStrackV2& trk);
106
107 void   AliD0vtxFinderSgn_pp_VTX(Int_t evFirst=0,Int_t evLast=999) {
108
109   const Char_t *name="AliD0vtxFinderSgn_pp_VTX";
110   cerr<<'\n'<<name<<"...\n";
111   gBenchmark->Start(name);  
112
113   // load library for D0 candidates class
114   gSystem->Load("libD0toKpi.so");
115
116   // check existence of input files 
117   if(!GetInputFiles()) { cerr<<"No tracks file found"<<endl; return; }
118
119   TString outNameSgn("AliD0toKpiSgn.root");
120   TString outNameBkg("AliD0toKpiBkgS.root");
121
122   AliKalmanTrack::SetConvConst(100/0.299792458/kBz);
123
124   Bool_t isSignal;
125   Int_t ev;
126   Int_t nTotEv=0,nD0recSgn=0,nD0recBkgS=0,nD0recBkg=0,nD0rec1ev=0;
127   Double_t dca;
128   Double_t v2[3],mom[6],d0[2];
129   Int_t pdg[2],mum[2],mumlab[2];
130   Double_t alphaP,alphaN,ptP,ptN,phiP,phiN;
131   Int_t l,k,iTrkP,iTrkN,itsEntries;
132   Int_t mesonD0,mesonD0label[10];
133   Int_t nTrksP=0,nTrksN=0;
134   Int_t okD0=0,okD0bar=0;
135   Char_t vtxName[100],trksName[100],refsName[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->SetMinTracks(2);
143   vertexer1->SetDebug(0);
144   Int_t skipped[2];
145
146   // define the cuts for vertexing
147   Double_t vtxcuts[]={33., // max. allowed chi2
148                       0.0, // min. allowed negative daughter's impact param 
149                       0.0, // min. allowed positive daughter's impact param 
150                       1.0, // max. allowed DCA between the daughter tracks
151                      -1.0, // min. allowed cosine of V0's pointing angle
152                       0.0, // min. radius of the fiducial volume
153                       2.9};// max. radius of the fiducial volume
154   
155   // create the AliV0vertexer object
156   AliV0vertexer *vertexer2 = new AliV0vertexer(vtxcuts);
157
158
159   // Create trees for reconstructed D0s
160   AliD0toKpi *ioD0toKpi=0;
161   // Signal
162   TTree D0TreeSgn("TreeD0","Tree with D0 candidates");
163   D0TreeSgn.Branch("D0toKpi","AliD0toKpi",&ioD0toKpi,200000,0);
164   // Background
165   TTree D0TreeBkg("TreeD0","Tree with D0 candidates");
166   D0TreeBkg.Branch("D0toKpi","AliD0toKpi",&ioD0toKpi,200000,0);
167
168   // Open file with ITS tracks
169   TFile* itstrks = TFile::Open("AliITStracksV2.root");
170
171   // Open file with ITS track references
172   TFile* itsrefs = TFile::Open("ITStracksRefFile.root");
173
174
175
176   // loop on events in file
177   for(ev=evFirst; ev<=evLast; ev++) {
178     printf(" --- Processing event  %d ---\n",ev);
179     sprintf(trksName,"TreeT_ITS_%d",ev);
180     sprintf(refsName,"Tree_Ref_%d",ev);
181
182
183
184     // tracks from ITS
185     TTree *itsTree=(TTree*)itstrks->Get(trksName);
186     if(!itsTree) continue;
187     itsEntries = (Int_t)itsTree->GetEntries();
188     printf("+++\n+++ Number of tracks in ITS: %d\n+++\n\n",itsEntries);
189     
190     // tree from reference file
191     TTree *refTree=(TTree*)itsrefs->Get(refsName);
192     refTree->SetBranchAddress("rectracks",&rectrk);
193     
194     Int_t refEntries = (Int_t)refTree->GetEntries();
195     k=0;
196     for(l=0; l<10; l++) mesonD0label[l]=-1;
197     for(l=0; l<refEntries; l++) {
198       refTree->GetEvent(l);
199       if(TMath::Abs(rectrk.mumpdg)!=421) continue;
200       //cerr<<" pdg: "<<rectrk.pdg<<"  mumpdg: "<<rectrk.mumpdg<<" mumlabel: "<<rectrk.mumlab<<endl;
201       mesonD0label[k] = rectrk.mumlab;
202       k++;
203     }
204
205     mesonD0 = mesonD0label[0];
206     for(l=0; l<10; l++) {
207       if(mesonD0label[l]==-1) continue;
208       for(k=0; k<10; k++) {
209         if(k==l) continue;
210         if(mesonD0label[k] == mesonD0label[l]) mesonD0 = mesonD0label[k]; 
211       }
212     }
213     //cerr<<" mesonD0: "<<mesonD0<<endl;
214
215     Double_t *brWgt = new Double_t[refEntries];
216     for(l=0; l<refEntries; l++) {
217       refTree->GetEvent(l);
218       // normally tracks have weight = 1
219       brWgt[l] = 1.;
220       // decay products of non-good D0 weighted with D0->Kpi B.R.
221       if(TMath::Abs(rectrk.mumpdg)==421 && rectrk.mumlab!=mesonD0) 
222         brWgt[l]=0.038;
223       // decay products of D+ weighted with D+->Kpipi B.R.
224       if(TMath::Abs(rectrk.mumpdg)==411) 
225         brWgt[l]=0.09;
226     }
227
228     // count the total number of events
229     nTotEv++;
230
231     // call function which applies sigle track selection and
232     // separetes positives and negatives
233     TObjArray trksP(itsEntries/2);
234     Int_t *itsEntryP = new Int_t[itsEntries];
235     TObjArray trksN(itsEntries/2);
236     Int_t *itsEntryN = new Int_t[itsEntries];
237     SelectTracks(*itsTree,trksP,itsEntryP,nTrksP,trksN,itsEntryN,nTrksN);
238       
239     nD0rec1ev = 0;
240     // loop on positive tracks
241     for(iTrkP=0; iTrkP<nTrksP; iTrkP++) {
242       //if(iTrkP % 10==0) cerr<<"--- Processing positive track number: "<<iTrkP<<" of "<<nTrksP<<"\r";
243           
244       // get track from track array
245       postrack = (AliITStrackV2*)trksP.At(iTrkP);
246
247       // get info on tracks PDG and mothers PDG from reference file
248       refTree->GetEvent(itsEntryP[iTrkP]);
249       pdg[0] = rectrk.pdg;
250       mum[0] = rectrk.mumpdg;
251       mumlab[0] = rectrk.mumlab;
252
253       // loop on negative tracks 
254       for(iTrkN=0; iTrkN<nTrksN; iTrkN++) {
255         // get track from track array
256         negtrack = (AliITStrackV2*)trksN.At(iTrkN);
257       
258         // get info on tracks PDG and mothers PDG from reference file
259         refTree->GetEvent(itsEntryN[iTrkN]);
260         pdg[1] = rectrk.pdg;
261         mum[1] = rectrk.mumpdg;
262         mumlab[1] = rectrk.mumlab;
263  
264         //
265         // ----------- DCA MINIMIZATION ------------------
266         //
267         // find the DCA and propagate the tracks to the DCA 
268         dca = vertexer2->PropagateToDCA(negtrack,postrack);
269
270         // define the AliV0vertex object
271         AliV0vertex *vertex2 = new AliV0vertex(*negtrack,*postrack);
272           
273         // get position of the vertex
274         vertex2->GetXYZ(v2[0],v2[1],v2[2]);
275
276         delete vertex2;
277           
278         // momenta of the tracks at the vertex
279         ptP = 1./TMath::Abs(postrack->Get1Pt());
280         alphaP = postrack->GetAlpha();
281         phiP = alphaP+TMath::ASin(postrack->GetSnp());
282         mom[0] = ptP*TMath::Cos(phiP); 
283         mom[1] = ptP*TMath::Sin(phiP);
284         mom[2] = ptP*postrack->GetTgl();
285           
286         ptN = 1./TMath::Abs(negtrack->Get1Pt());
287         alphaN = negtrack->GetAlpha();
288         phiN = alphaN+TMath::ASin(negtrack->GetSnp());
289         mom[3] = ptN*TMath::Cos(phiN); 
290         mom[4] = ptN*TMath::Sin(phiN);
291         mom[5] = ptN*negtrack->GetTgl();
292
293         Bool_t goodVtx = kFALSE;
294         // no vertexing if DeltaMass > kMassCut 
295         if(SelectInvMass(mom)) {
296           // primary vertex from other tracks in event
297           vertexer1->SetVtxStart(0.,0.);
298           skipped[0] = itsEntryP[iTrkP];
299           skipped[1] = itsEntryN[iTrkN];
300           vertexer1->SetSkipTracks(2,skipped);
301           AliITSVertex *vertex1 = 
302             (AliITSVertex*)vertexer1->VertexOnTheFly(*itsTree); 
303           if(vertex1->GetNContributors()>0) goodVtx = kTRUE;
304           vertex1->GetXYZ(gv1);
305           // impact parameters of the tracks w.r.t. the primary vertex
306           d0[0] =  10000.*postrack->GetD(gv1[0],gv1[1]);
307           d0[1] = -10000.*negtrack->GetD(gv1[0],gv1[1]);
308           //vertex1->PrintStatus();
309           delete vertex1;
310         }
311
312         isSignal = kFALSE;
313         if(mumlab[0]==mesonD0 && mumlab[1]==mesonD0) isSignal = kTRUE;    
314
315
316         // create the object TD0rec and store it in the tree
317         AliD0toKpi theD0(isSignal,ev,gv1,v2,dca,mom,d0,pdg,mum);
318
319         // select D0s
320         if(goodVtx && theD0.Select(kD0cutsAll,okD0,okD0bar)) {
321               
322           // compute the weights
323           theD0.ComputeWgts();
324           theD0.CorrectWgt4BR(brWgt[itsEntryN[iTrkN]]*brWgt[itsEntryP[iTrkP]]);
325
326           // fill the tree
327           ioD0toKpi=&theD0;
328           if(isSignal) { 
329             D0TreeSgn.Fill(); 
330             nD0recSgn++; 
331           } else { 
332             D0TreeBkg.Fill();
333             if(TMath::Abs(mum[0])==421 || TMath::Abs(mum[0])==411 ||
334                TMath::Abs(mum[1])==421 || TMath::Abs(mum[1])==411) {
335               nD0recBkgS++;
336             } else {
337               nD0recBkg++;
338             }  
339           }
340
341           nD0rec1ev++;
342
343           ioD0toKpi=0;
344         } // end if select
345
346         negtrack = 0;
347       } // loop on negative tracks
348       postrack = 0;
349     }   // loop on positive tracks
350
351     delete [] itsEntryP;
352     delete [] itsEntryN;
353     delete [] brWgt;
354     delete itsTree;
355     delete refTree;
356
357     printf("\n+++\n+++ Number of D0 candidates: %d\n+++\n",nD0rec1ev);
358     printf("----------------------------------------------------------\n");
359
360   }    // loop on events in file
361
362
363   printf("\n+++\n+++ Total number of events: %d\n+++\n",nTotEv);
364   printf("\n+++\n+++ Total number of D0 candidates:\n
365                  +++    Sgn:   %d\n
366                  +++    BkgS:  %d\n
367                  +++    Bkg:   %d\n+++\n",nD0recSgn,nD0recBkgS,nD0recBkg);
368
369   delete vertexer2;
370
371   itstrks->Close();
372   itsrefs->Close();
373
374   // store trees in files
375   TFile* outrootSgn = new TFile(outNameSgn.Data(),"recreate");
376   D0TreeSgn.Write();
377   outrootSgn->Close();
378   delete outrootSgn;
379   TFile* outrootBkg = new TFile(outNameBkg.Data(),"recreate");
380   D0TreeBkg.Write();
381   outrootBkg->Close();
382   delete outrootBkg;
383
384   gBenchmark->Stop(name);
385   gBenchmark->Show(name);
386
387   return;
388 }
389 //___________________________________________________________________________
390 Bool_t   GetInputFiles() {
391   TString itsName("AliITStracksV2.root");
392   TString refName("ITStracksRefFile.root");
393
394   if(gSystem->AccessPathName(itsName.Data(),kFileExists) ||
395      gSystem->AccessPathName(refName.Data(),kFileExists)) return kFALSE;
396
397   return kTRUE;
398 }
399 //___________________________________________________________________________
400 Bool_t SelectInvMass(const Double_t p[6]) {
401
402   Double_t mD0 = 1.8645;
403   Double_t mPi = 0.13957;
404   Double_t mKa = 0.49368;
405
406   Double_t energy[2];
407   Double_t mom2[2],momTot2;
408
409   mom2[0] = p[0]*p[0]+p[1]*p[1]+p[2]*p[2];
410   mom2[1] = p[3]*p[3]+p[4]*p[4]+p[5]*p[5];
411
412   momTot2 = (p[0]+p[3])*(p[0]+p[3])+
413             (p[1]+p[4])*(p[1]+p[4])+
414             (p[2]+p[5])*(p[2]+p[5]);
415
416   // D0 -> K- Pi+
417   energy[1] = TMath::Sqrt(mKa*mKa+mom2[1]);
418   energy[0] = TMath::Sqrt(mPi*mPi+mom2[0]);
419
420   Double_t minvD0 = TMath::Sqrt((energy[0]+energy[1])*(energy[0]+energy[1])-momTot2);
421     
422
423   // D0bar -> K+ Pi-
424   energy[0] = TMath::Sqrt(mKa*mKa+mom2[0]);
425   energy[1] = TMath::Sqrt(mPi*mPi+mom2[1]);
426
427   Double_t minvD0bar = TMath::Sqrt((energy[0]+energy[1])*(energy[0]+energy[1])-momTot2);
428
429   if(TMath::Abs(minvD0-mD0) < kMassCut)    return kTRUE;
430   if(TMath::Abs(minvD0bar-mD0) < kMassCut) return kTRUE;
431   return kFALSE;
432 }
433 //___________________________________________________________________________
434 void   SelectTracks(TTree& itsTree,
435                     TObjArray& trksP,Int_t* itsEntryP,Int_t& nTrksP,
436                     TObjArray& trksN,Int_t* itsEntryN,Int_t& nTrksN) {
437 //
438 // this function creates two TObjArrays with positive and negative tracks
439 //
440   nTrksP=0,nTrksN=0;
441
442  
443   Int_t entr = (Int_t)itsTree.GetEntries();
444
445   // trasfer tracks from tree to arrays
446   for(Int_t i=0; i<entr; i++) {
447
448     AliITStrackV2 *itstrack = new AliITStrackV2; 
449     itsTree.SetBranchAddress("tracks",&itstrack);
450
451     itsTree.GetEvent(i);
452
453     // single track selection
454     if(!TrkCuts(*itstrack)) { delete itstrack; continue; }
455
456     if(itstrack->Get1Pt()>0.) { // negative track
457       trksN.AddLast(itstrack);
458       itsEntryN[nTrksN] = i;
459       nTrksN++;
460     } else {                    // positive track
461       trksP.AddLast(itstrack);
462       itsEntryP[nTrksP] = i;
463       nTrksP++;
464     }
465
466   }
467
468   return;
469 }
470 //____________________________________________________________________________
471 Bool_t TrkCuts(const AliITStrackV2& trk) {
472 // 
473 // this function tells if track passes some kinematical cuts  
474 //
475   if(TMath::Abs(1./trk.Get1Pt()) < kPtCut)                return kFALSE;
476   if(TMath::Abs(10000.*trk.GetD(gv1[0],gv1[1])) < kd0Cut) return kFALSE;
477
478   return kTRUE;
479 }
480 //____________________________________________________________________________
481
482
483
484
485
486
487
488
489
490
491