Merged Cvetans RowHoughTransformer, Anders latest developments in comp
[u/mrichter/AliRoot.git] / HLT / trigger / RunD0offline.C
CommitLineData
0bd0c1ef 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//sec. vertex
38double v2[3]={0,0,0};
39
40// sigle track cuts
41const Double_t kPtCut = 0.5; // GeV/c
42const Double_t kd0Cut = 50.; // micron
43const Double_t kd0CutHigh = 200.; // micron
44
45
46//cuts for combined tracks
47const Double_t cuts[7] = {0.005, // cuts[0] = lowest V0 cut (cm)
48 0.015, // cuts[1] = highest V0 cut (cm)
49 0.05, // cuts[2] = inv. mass cut (diferense) (Gev/c)
50 0.95, // cuts[3] = max cosine value for pointing angle
51 -5000, // cuts[4] = d0d0
52 0.8, // cuts[5] = costhetastar
53 0.5}; // cuts[6] = ptchild
54//cut for distance of closest aprach
55double cutDCA=0.01;
56
57// this function applies single track cuts
58Bool_t TrkCuts(const AliITStrackV2& trk);
59
60// this function creates TObjArrays with positive and negative tracks
61void SelectTracks(TTree& itsTree,
62 TObjArray& trksP,Int_t* itsEntryP,Int_t& nTrksP,
63 TObjArray& trksN,Int_t* itsEntryN,Int_t& nTrksN);
64
65//void GetPrimaryVertex(int i,Char_t* path="./");
66
67Int_t iTrkP,iTrkN,itsEntries;
68Char_t trksName[100];
69Int_t nTrksP=0,nTrksN=0;
70Int_t nD0=0;
71int ev=0;
72double mom[6];
73
74void RunD0offline(Int_t evFirst=0,Int_t evLast=1,Char_t* path="./") {
75
76 const Char_t *name="AliD0offline";
77 cerr<<'\n'<<name<<"...\n";
78 gBenchmark->Start(name);
79
80 AliKalmanTrack::SetConvConst(100/0.299792458/kBz);
81
82 // Open file with ITS tracks
83 Char_t fnameTrack[1024];
84 sprintf(fnameTrack,"%s/AliITStracksV2.root",path);
85 TFile* itstrks = TFile::Open(fnameTrack);
86
87 // tracks from ITS
88 sprintf(trksName,"TreeT_ITS_%d",ev);
89 TTree *itsTree=(TTree*)itstrks->Get(trksName);
90 if(!itsTree) continue;
91 itsEntries = (Int_t)itsTree->GetEntries();
92 printf("+++\n+++ Number of tracks in ITS: %d\n+++\n\n",itsEntries);
93
94 //Getting primary Vertex
95 GetPrimaryVertex(0,path);
96
97 // call function which applies sigle track selection and
98 // separetes positives and negatives
99 TObjArray trksP(itsEntries/2);
100 Int_t *itsEntryP = new Int_t[itsEntries];
101 TObjArray trksN(itsEntries/2);
102 Int_t *itsEntryN = new Int_t[itsEntries];
103 SelectTracks(*itsTree,trksP,itsEntryP,nTrksP,trksN,itsEntryN,nTrksN);
104
105 cout<<"#pos: "<<nTrksP<<endl;
106 cout<<"#neg: "<<nTrksN<<endl;
107
108 //the offline stuff
109 // define the cuts for vertexing
110 Double_t vtxcuts[]={33., // max. allowed chi2
111 0.0, // min. allowed negative daughter's impact param
112 0.0, // min. allowed positive daughter's impact param
113 1.0, // max. allowed DCA between the daughter tracks
114 -1.0, // min. allowed cosine of V0's pointing angle
115 0.0, // min. radius of the fiducial volume
116 2.9};// max. radius of the fiducial volume
117
118 // create the AliV0vertexer object
119 AliV0vertexer *vertexer2 = new AliV0vertexer(vtxcuts);
120
121 AliD0Trigger * D0 = new AliD0Trigger(cuts,kBz,primaryvertex);
122
123 double ptP,alphaP,phiP,ptN,alphaN,phiN,dca;
124
125 for(iTrkP=0; iTrkP<nTrksP; iTrkP++) {
126 postrack = (AliITStrackV2*)trksP.At(iTrkP);
127 for(iTrkN=0; iTrkN<nTrksN; iTrkN++) {
128 negtrack = (AliITStrackV2*)trksN.At(iTrkN);
129 D0.SetTracks(postrack,negtrack);
130 //
131 // ----------- DCA MINIMIZATION ------------------
132 //
133 // find the DCA and propagate the tracks to the DCA
134 double dca = vertexer2->PropagateToDCA(negtrack,postrack);
135
136 if(dca<cutDCA){
137 // define the AliV0vertex object
138 AliV0vertex *vertex2 = new AliV0vertex(*negtrack,*postrack);
139 // get position of the vertex
140 vertex2->GetXYZ(v2[0],v2[1],v2[2]);
141 delete vertex2;
142 if(D0.pTchild()){
143 if(D0.d0d0()){
144 if(D0.FindV0offline(v2)){
145
146 // momenta of the tracks at the vertex
147 ptP = 1./TMath::Abs(postrack->Get1Pt());
148 alphaP = postrack->GetAlpha();
149 phiP = alphaP+TMath::ASin(postrack->GetSnp());
150 mom[0] = ptP*TMath::Cos(phiP);
151 mom[1] = ptP*TMath::Sin(phiP);
152 mom[2] = ptP*postrack->GetTgl();
153
154 ptN = 1./TMath::Abs(negtrack->Get1Pt());
155 alphaN = negtrack->GetAlpha();
156 phiN = alphaN+TMath::ASin(negtrack->GetSnp());
157 mom[3] = ptN*TMath::Cos(phiN);
158 mom[4] = ptN*TMath::Sin(phiN);
159 mom[5] = ptN*negtrack->GetTgl();
160
161 D0.SetMomenta(mom);
162
163 if(D0.FindInvMass()){
164 if(D0.CosThetaStar()){
165 if(D0.PointingAngle()){
166 nD0++;
167 }
168 }
169 }
170 }
171 }
172 }
173 }
174 }
175 }
176 cout<<"#D0: "<<nD0<<endl;
177 gBenchmark->Stop(name);
178 gBenchmark->Show(name);
179}
180//___________________________________________________________________________
181void SelectTracks(TTree& itsTree,
182 TObjArray& trksP,Int_t* itsEntryP,Int_t& nTrksP,
183 TObjArray& trksN,Int_t* itsEntryN,Int_t& nTrksN) {
184//
185// this function creates two TObjArrays with positive and negative tracks
186//
187 nTrksP=0,nTrksN=0;
188
189
190 Int_t entr = (Int_t)itsTree.GetEntries();
191
192 // trasfer tracks from tree to arrays
193 for(Int_t i=0; i<entr; i++) {
194
195 AliITStrackV2 *itstrack = new AliITStrackV2;
196 itsTree.SetBranchAddress("tracks",&itstrack);
197
198 itsTree.GetEvent(i);
199
200 // single track selection
201 if(!TrkCuts(*itstrack)) { delete itstrack; continue; }
202
203 if(itstrack->Get1Pt()>0.) { // negative track
204 trksN.AddLast(itstrack);
205 itsEntryN[nTrksN] = i;
206 nTrksN++;
207 } else { // positive track
208 trksP.AddLast(itstrack);
209 itsEntryP[nTrksP] = i;
210 nTrksP++;
211 }
212
213 }
214
215 return;
216}
217//____________________________________________________________________________
218Bool_t TrkCuts(const AliITStrackV2& trk) {
219//
220// this function tells if track passes some kinematical cuts
221//
222 if(TMath::Abs(1./trk.Get1Pt()) < kPtCut) return kFALSE;
223 if(TMath::Abs(10000.*trk.GetD(primaryvertex[0],primaryvertex[1])) < kd0Cut) return kFALSE;
224 if(TMath::Abs(10000.*trk.GetD(primaryvertex[0],primaryvertex[1])) > kd0CutHigh) return kFALSE;
225
226 return kTRUE;
227}
228//____________________________________________________________________________
229void GetPrimaryVertex(int i,Char_t* path="./") {
230
231 int event=i;
232
233 Char_t falice[1024];
234 sprintf(falice,"%s/galice.root",path);
235 TFile * galice = new TFile(falice);
236
237 TDirectory * curdir;
238
239 Char_t vname[20];
240 galice->cd();
241
242 sprintf(vname,"Vertex_%d",event);
243 TArrayF o = 0;
244 o.Set(3);
245 AliHeader * header = 0;
246
247 TTree * treeE = (TTree*)gDirectory->Get("TE");
248 treeE->SetBranchAddress("Header",&header);
249 treeE->GetEntry(event);
250 AliGenEventHeader* genHeader = header->GenEventHeader();
251 if(genHeader){
252 genHeader->PrimaryVertex(o);
253 primaryvertex[0] = (Double_t)o[0];
254 primaryvertex[1] = (Double_t)o[1];
255 primaryvertex[2] = (Double_t)o[2];
256 }
257 else{
258 printf("Can't find Header");
259 }
260 delete header;
261 delete galice;
262}