ALIROOT-5433 Transition to CDHv3 in HLT
[u/mrichter/AliRoot.git] / HLT / trigger / RunD0Trigger.C
CommitLineData
5a31e9df 1//#define __COMPILE__
2//#ifdef __COMPILE__
3#if !defined(__CINT__) || defined(__MAKECINT__)
4//-- --- standard headers-------------
5#include <iostream.h>
6//--------Root headers ---------------
7#include <TSystem.h>
8#include <TFile.h>
9#include <TString.h>
10#include <TStopwatch.h>
11#include <TObject.h>
12#include <TVector3.h>
13#include <TTree.h>
14#include <TParticle.h>
15#include <TArray.h>
16//----- AliRoot headers ---------------
17#include "alles.h"
18#include "AliRun.h"
19#include "AliKalmanTrack.h"
20#include "AliITStrackV2.h"
21#include "AliHeader.h"
22#include "AliGenEventHeader.h"
23#include "AliV0vertex.h"
24#include "AliV0vertexer.h"
25#include "AliITSVertex.h"
26#include "AliITSVertexer.h"
27#include "AliITSVertexerTracks.h"
28#include "AliD0Trigger.h"
29#endif
30//-------------------------------------
31// field (T)
32const Double_t kBz = 0.4;
33
34// primary vertex
35Double_t primaryvertex[3] = {0.,0.,0,};
36
37// sigle track cuts
0bd0c1ef 38const Double_t kPtCut = 0.5; // GeV/c
39const Double_t kd0Cut = 50.; // micron
5a31e9df 40const Double_t kd0CutHigh = 200.; // micron
41
42
43//cuts for combined tracks
0bd0c1ef 44const Double_t cuts[5] = {0.005, // cuts[0] = lowest V0 cut (cm)
45 200, // cuts[1] = highest V0 cut (cm) //0.02
46 0.1, // cuts[2] = inv. mass cut (diferense) (Gev/c)
47 0.95, // cuts[3] = max cosine value for pointing angle
48 -5000}; // cuts[4] = max cosine value for pointing angle
5a31e9df 49
50Int_t iTrkP,iTrkN,itsEntries;
51Char_t trksName[100];
52Int_t nTrksP=0,nTrksN=0;
53Int_t nD0=0;
54int ev=0;
55
0bd0c1ef 56void RunD0Trigger(Int_t evFirst=0,Int_t evLast=1,Char_t* path="./") {
5a31e9df 57
58 const Char_t *name="AliD0Trigger";
59 cerr<<'\n'<<name<<"...\n";
60 gBenchmark->Start(name);
61
62 AliKalmanTrack::SetConvConst(100/0.299792458/kBz);
63
64 // Open file with ITS tracks
65 Char_t fnameTrack[1024];
66 sprintf(fnameTrack,"%s/AliITStracksV2.root",path);
67 TFile* itstrks = TFile::Open(fnameTrack);
68
69 // tracks from ITS
70 sprintf(trksName,"TreeT_ITS_%d",ev);
71 TTree *itsTree=(TTree*)itstrks->Get(trksName);
72 if(!itsTree) continue;
73 itsEntries = (Int_t)itsTree->GetEntries();
74 printf("+++\n+++ Number of tracks in ITS: %d\n+++\n\n",itsEntries);
75
76 //Getting primary Vertex
77 GetPrimaryVertex(0,path);
78
79 // call function which applies sigle track selection and
80 // separetes positives and negatives
81 TObjArray trksP(itsEntries/2);
82 Int_t *itsEntryP = new Int_t[itsEntries];
83 TObjArray trksN(itsEntries/2);
84 Int_t *itsEntryN = new Int_t[itsEntries];
85 SelectTracks(*itsTree,trksP,itsEntryP,nTrksP,trksN,itsEntryN,nTrksN);
86
87 cout<<"#pos: "<<nTrksP<<endl;
88 cout<<"#neg: "<<nTrksN<<endl;
89
90 AliD0Trigger * D0 = new AliD0Trigger(cuts,kBz,primaryvertex);
91
92 for(iTrkP=0; iTrkP<nTrksP; iTrkP++) {
93 postrack = (AliITStrackV2*)trksP.At(iTrkP);
94 for(iTrkN=0; iTrkN<nTrksN; iTrkN++) {
95 negtrack = (AliITStrackV2*)trksN.At(iTrkN);
96 D0.SetTracks(postrack,negtrack);
97 if(D0.FindV0()){
98 D0.FindMomentaAtVertex();
99 if(D0.FindInvMass()){
100 if(D0.PointingAngle()){
101 nD0++;
102 }
103 }
104 }
105 }
106 }
107 cout<<"#D0: "<<nD0<<endl;
108 gBenchmark->Stop(name);
109 gBenchmark->Show(name);
110}
111//___________________________________________________________________________
112void SelectTracks(TTree& itsTree,
113 TObjArray& trksP,Int_t* itsEntryP,Int_t& nTrksP,
114 TObjArray& trksN,Int_t* itsEntryN,Int_t& nTrksN) {
115//
116// this function creates two TObjArrays with positive and negative tracks
117//
118 nTrksP=0,nTrksN=0;
119
120
121 Int_t entr = (Int_t)itsTree.GetEntries();
122
123 // trasfer tracks from tree to arrays
124 for(Int_t i=0; i<entr; i++) {
125
126 AliITStrackV2 *itstrack = new AliITStrackV2;
127 itsTree.SetBranchAddress("tracks",&itstrack);
128
129 itsTree.GetEvent(i);
130
131 // single track selection
132 if(!TrkCuts(*itstrack)) { delete itstrack; continue; }
133
134 if(itstrack->Get1Pt()>0.) { // negative track
135 trksN.AddLast(itstrack);
136 itsEntryN[nTrksN] = i;
137 nTrksN++;
138 } else { // positive track
139 trksP.AddLast(itstrack);
140 itsEntryP[nTrksP] = i;
141 nTrksP++;
142 }
143
144 }
145
146 return;
147}
148//____________________________________________________________________________
149Bool_t TrkCuts(const AliITStrackV2& trk) {
150//
151// this function tells if track passes some kinematical cuts
152//
153 if(TMath::Abs(1./trk.Get1Pt()) < kPtCut) return kFALSE;
154 if(TMath::Abs(10000.*trk.GetD(primaryvertex[0],primaryvertex[1])) < kd0Cut) return kFALSE;
155 if(TMath::Abs(10000.*trk.GetD(primaryvertex[0],primaryvertex[1])) > kd0CutHigh) return kFALSE;
156
157 return kTRUE;
158}
159//____________________________________________________________________________
160void GetPrimaryVertex(int i,Char_t* path="./") {
161
162 int event=i;
163
164 Char_t falice[1024];
165 sprintf(falice,"%s/galice.root",path);
166 TFile * galice = new TFile(falice);
167
168 TDirectory * curdir;
169
170 Char_t vname[20];
171 galice->cd();
172
173 sprintf(vname,"Vertex_%d",event);
174 TArrayF o = 0;
175 o.Set(3);
176 AliHeader * header = 0;
177
178 TTree * treeE = (TTree*)gDirectory->Get("TE");
179 treeE->SetBranchAddress("Header",&header);
180 treeE->GetEntry(event);
181 AliGenEventHeader* genHeader = header->GenEventHeader();
182 if(genHeader){
183 genHeader->PrimaryVertex(o);
184 primaryvertex[0] = (Double_t)o[0];
185 primaryvertex[1] = (Double_t)o[1];
186 primaryvertex[2] = (Double_t)o[2];
187 }
188 else{
189 printf("Can't find Header");
190 }
191 delete header;
192 delete galice;
193}