added PMD raw digits reading class
[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
de3c3890 35double primaryvertex[3]={0.,0.,0.};
0bd0c1ef 36
37//sec. vertex
38double v2[3]={0,0,0};
39
40// sigle track cuts
de3c3890 41const Double_t kPtCut = 0.5; // 0.5 GeV/c
42const Double_t kd0Cut = 50.; // 50 micron
43const Double_t kd0CutHigh = 400.; // 200 micron
0bd0c1ef 44
45
46//cuts for combined tracks
de3c3890 47const Double_t cuts[7] = {0.005, // 0.005 cuts[0] = lowest V0 cut (cm)
48 0.800, // 0.015 cuts[1] = highest V0 cut (cm)
49 0.012, // 0.012 cuts[2] = inv. mass cut (diferense) (Gev/c)
50 0.8, // 0.8 cuts[3] = min. cosine value for pointing angle
51 -5000, // -5000 cuts[4] = d0d0
52 0, // 0.8 cuts[5] = costhetastar
53 0.5}; // 0.5 cuts[6] = ptchild
0bd0c1ef 54//cut for distance of closest aprach
de3c3890 55double cutDCA=0.05; //0.05
0bd0c1ef 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
de3c3890 67//void PtD0(Char_t* path="./");
68
0bd0c1ef 69Int_t iTrkP,iTrkN,itsEntries;
70Char_t trksName[100];
71Int_t nTrksP=0,nTrksN=0;
72Int_t nD0=0;
73int ev=0;
74double mom[6];
de3c3890 75int event[10000];
76int index=0;
0bd0c1ef 77
de3c3890 78void RunD0offline(Char_t* path="./",bool h=false,bool PtD0=false) {
79 Char_t falice[1024];
80 sprintf(falice,"%s/galice.root",path);
81 TFile *f = new TFile(falice);
82 gAlice=(AliRun*)f->Get("gAlice");
83 int nEvent=gAlice->GetEventsPerRun();
84 //int nEvent=20;
85 delete gAlice;
0bd0c1ef 86
87 const Char_t *name="AliD0offline";
88 cerr<<'\n'<<name<<"...\n";
de3c3890 89 cout<<"#Events: "<<nEvent<<endl;
0bd0c1ef 90 gBenchmark->Start(name);
91
92 AliKalmanTrack::SetConvConst(100/0.299792458/kBz);
93
94 // Open file with ITS tracks
95 Char_t fnameTrack[1024];
de3c3890 96 //sprintf(fnameTrack,"%s/AliITStracksV2.root",path);
0bd0c1ef 97 sprintf(fnameTrack,"%s/AliITStracksV2.root",path);
98 TFile* itstrks = TFile::Open(fnameTrack);
99
0bd0c1ef 100 //the offline stuff
101 // define the cuts for vertexing
102 Double_t vtxcuts[]={33., // max. allowed chi2
103 0.0, // min. allowed negative daughter's impact param
104 0.0, // min. allowed positive daughter's impact param
105 1.0, // max. allowed DCA between the daughter tracks
de3c3890 106 -1.0, // min. allowed cosine of V0's pointing angle
0bd0c1ef 107 0.0, // min. radius of the fiducial volume
108 2.9};// max. radius of the fiducial volume
0bd0c1ef 109
de3c3890 110 TH1F *h1 = new TH1F("h1","Transvers momentun of reconstructed D0",100,0,10);
111 TH1F *h2 = new TH1F("h2","Transvers momentun of D0 with |Eta|<0.9",100,0,10);
112 TH1F *h3 = new TH1F("h3","Eta reconstructed of D0",100,-5,5);
113 TH1F *h4 = new TH1F("h4","Eta of D0",100,-5,5);
114
115 for(ev=0;ev<nEvent;ev++) {
116
117 cout<<"\n Event: "<<ev<<endl;
118
119 // tracks from ITS
120 sprintf(trksName,"TreeT_ITS_%d",ev);
121 TTree *itsTree=(TTree*)itstrks->Get(trksName);
122 if(!itsTree) continue;
123 itsEntries = (Int_t)itsTree->GetEntries();
124 printf("+++\n+++ Number of tracks in ITS: %d\n+++\n\n",itsEntries);
125
126 //Getting primary Vertex
127 GetPrimaryVertex(ev,path);
128
129 // call function which applies sigle track selection and
130 // separetes positives and negatives
131 TObjArray trksP(itsEntries/2);
132 Int_t *itsEntryP = new Int_t[itsEntries];
133 TObjArray trksN(itsEntries/2);
134 Int_t *itsEntryN = new Int_t[itsEntries];
135 SelectTracks(*itsTree,trksP,itsEntryP,nTrksP,trksN,itsEntryN,nTrksN);
136
137 cout<<"#pos: "<<nTrksP<<endl;
138 cout<<"#neg: "<<nTrksN<<endl;
139
140 // create the AliV0vertexer object
141 AliV0vertexer *vertexer2 = new AliV0vertexer(vtxcuts);
142
143 AliD0Trigger * D0 = new AliD0Trigger(cuts,kBz,primaryvertex);
144
145 double ptP,alphaP,phiP,ptN,alphaN,phiN,dca;
146
147 for(iTrkP=0; iTrkP<nTrksP; iTrkP++) {
148 postrack = (AliITStrackV2*)trksP.At(iTrkP);
149 for(iTrkN=0; iTrkN<nTrksN; iTrkN++) {
150 negtrack = (AliITStrackV2*)trksN.At(iTrkN);
151 D0->SetTracks(postrack,negtrack);
152 //
153 // ----------- DCA MINIMIZATION ------------------
154 //
155 // find the DCA and propagate the tracks to the DCA
156 double dca = vertexer2->PropagateToDCA(negtrack,postrack);
157
158 if(dca<cutDCA){
159 // define the AliV0vertex object
160 AliV0vertex *vertex2 = new AliV0vertex(*negtrack,*postrack);
161 // get position of the vertex
162 vertex2->GetXYZ(v2[0],v2[1],v2[2]);
163 delete vertex2;
164 //if(D0->FindV0offline(v2) && D0->d0d0()){
165 if(D0->d0d0()){
166 D0->SetV0(v2);
167 D0->FindMomentaOffline();
168 if(D0->FindInvMass() && D0->CosThetaStar() && D0->PointingAngle() && D0->pTchild()){
169 nD0++;
170 event[index]=ev; index++;
171 if(h){
172 h1->Fill(D0->Pt());
173 h3->Fill(D0->Eta());
0bd0c1ef 174 }
de3c3890 175 }
176 }
0bd0c1ef 177 }
178 }
de3c3890 179 }
180 }
181 cout<<"\n#D0: "<<nD0<<endl;
182 for(int i=0;i<nD0;i++)
183 cout<<event[i]<<endl;
0bd0c1ef 184 gBenchmark->Stop(name);
185 gBenchmark->Show(name);
de3c3890 186
187 if(h){
188 if(!PtD0){
189 TCanvas *c = new TCanvas("c","c",700,1000);
190 c->Divide(1,2);
191 c->cd(1);
192 h1->Draw();
193 pt = new TPaveText(7.3,0.4,11,1.8,"br");
194 pt->SetFillColor(0);
195 pt->SetBorderSize(1);
196 pt->AddText("Cuts:");
197 Char_t st[1024];
198 sprintf(st,"First Pt: %g",kPtCut);
199 pt->AddText(st);
200 sprintf(st,"d0 low: %g",kd0Cut);
201 pt->AddText(st);
202 sprintf(st,"d0 high: %g",kd0CutHigh);
203 pt->AddText(st);
204 sprintf(st,"V0 low: %g",cuts[0]);
205 pt->AddText(st);
206 sprintf(st,"V0 high: %g",cuts[1]);
207 pt->AddText(st);
208 sprintf(st,"InvMass Diff: %g",cuts[2]);
209 pt->AddText(st);
210 sprintf(st,"cosPointingAngle: %g",cuts[3]);
211 pt->AddText(st);
212 sprintf(st,"d0d0: %g",cuts[4]);
213 pt->AddText(st);
214 sprintf(st,"cosThetaStar: %g",cuts[5]);
215 pt->AddText(st);
216 sprintf(st,"PtChild: %g",cuts[6]);
217 pt->AddText(st);
218 sprintf(st,"DCA: %g",cutDCA);
219 pt->AddText(st);
220 pt->Draw();
221 c_1->Modified();
222 c->cd();
223 c->cd(2);
224 h3->Draw();
225 }
226 if(PtD0){
227 PtD0(path,h2,h4);
228 TCanvas *c = new TCanvas("c","c",1000,700);
229 c->Divide(2,2);
230 c->cd(1);
231 h1->Draw();
232 pt = new TPaveText(7.3,0.4,11,1.8,"br");
233 pt->SetFillColor(0);
234 pt->SetBorderSize(1);
235 pt->AddText("Cuts:");
236 Char_t st[1024];
237 sprintf(st,"First Pt: %g",kPtCut);
238 pt->AddText(st);
239 sprintf(st,"d0 low: %g",kd0Cut);
240 pt->AddText(st);
241 sprintf(st,"d0 high: %g",kd0CutHigh);
242 pt->AddText(st);
243 sprintf(st,"V0 low: %g",cuts[0]);
244 pt->AddText(st);
245 sprintf(st,"V0 high: %g",cuts[1]);
246 pt->AddText(st);
247 sprintf(st,"InvMass Diff: %g",cuts[2]);
248 pt->AddText(st);
249 sprintf(st,"cosPointingAngle: %g",cuts[3]);
250 pt->AddText(st);
251 sprintf(st,"d0d0: %g",cuts[4]);
252 pt->AddText(st);
253 sprintf(st,"cosThetaStar: %g",cuts[5]);
254 pt->AddText(st);
255 sprintf(st,"PtChild: %g",cuts[6]);
256 pt->AddText(st);
257 sprintf(st,"DCA: %g",cutDCA);
258 pt->AddText(st);
259 pt->Draw();
260 c_1->Modified();
261 c->cd();
262 c->cd(2);
263 h3->Draw();
264 c->cd(3);
265 h2->Draw();
266 c->cd(4);
267 h4->Draw();
268 }
269 //TFile *outfile = new TFile("results.root","recreate");
270 //h1->Write();
271 //outfile->Close();
272
273 }
0bd0c1ef 274}
275//___________________________________________________________________________
276void SelectTracks(TTree& itsTree,
277 TObjArray& trksP,Int_t* itsEntryP,Int_t& nTrksP,
278 TObjArray& trksN,Int_t* itsEntryN,Int_t& nTrksN) {
279//
280// this function creates two TObjArrays with positive and negative tracks
281//
282 nTrksP=0,nTrksN=0;
283
284
285 Int_t entr = (Int_t)itsTree.GetEntries();
286
287 // trasfer tracks from tree to arrays
288 for(Int_t i=0; i<entr; i++) {
289
290 AliITStrackV2 *itstrack = new AliITStrackV2;
291 itsTree.SetBranchAddress("tracks",&itstrack);
292
293 itsTree.GetEvent(i);
294
295 // single track selection
296 if(!TrkCuts(*itstrack)) { delete itstrack; continue; }
297
298 if(itstrack->Get1Pt()>0.) { // negative track
299 trksN.AddLast(itstrack);
300 itsEntryN[nTrksN] = i;
301 nTrksN++;
302 } else { // positive track
303 trksP.AddLast(itstrack);
304 itsEntryP[nTrksP] = i;
305 nTrksP++;
306 }
307
308 }
309
310 return;
311}
312//____________________________________________________________________________
313Bool_t TrkCuts(const AliITStrackV2& trk) {
314//
315// this function tells if track passes some kinematical cuts
316//
317 if(TMath::Abs(1./trk.Get1Pt()) < kPtCut) return kFALSE;
318 if(TMath::Abs(10000.*trk.GetD(primaryvertex[0],primaryvertex[1])) < kd0Cut) return kFALSE;
de3c3890 319 //if(TMath::Abs(10000.*trk.GetD(primaryvertex[0],primaryvertex[1])) > kd0CutHigh) return kFALSE;
0bd0c1ef 320
321 return kTRUE;
322}
323//____________________________________________________________________________
324void GetPrimaryVertex(int i,Char_t* path="./") {
325
326 int event=i;
327
328 Char_t falice[1024];
329 sprintf(falice,"%s/galice.root",path);
330 TFile * galice = new TFile(falice);
331
332 TDirectory * curdir;
333
334 Char_t vname[20];
335 galice->cd();
336
337 sprintf(vname,"Vertex_%d",event);
338 TArrayF o = 0;
339 o.Set(3);
340 AliHeader * header = 0;
341
342 TTree * treeE = (TTree*)gDirectory->Get("TE");
343 treeE->SetBranchAddress("Header",&header);
344 treeE->GetEntry(event);
345 AliGenEventHeader* genHeader = header->GenEventHeader();
346 if(genHeader){
347 genHeader->PrimaryVertex(o);
348 primaryvertex[0] = (Double_t)o[0];
349 primaryvertex[1] = (Double_t)o[1];
350 primaryvertex[2] = (Double_t)o[2];
351 }
352 else{
353 printf("Can't find Header");
354 }
355 delete header;
356 delete galice;
357}
de3c3890 358//____________________________________________________________________________
359PtD0(Char_t* path="./",TH1F * h1,TH1F * h4){
360
361 Char_t falice[1024];
362 sprintf(falice,"%s/galice.root",path);
363 TFile *f = new TFile(falice);
364 gAlice=(AliRun*)f->Get("gAlice");
365
366 TParticle *p;
367 int nd0=0;
368 int nkminus =0;
369 int npipluss = 0;
370 int nEvent=gAlice->GetEventsPerRun();
371
372 for (Int_t i = 0; i <nEvent; i++) {
373 cout<<"Event: "<<i<<endl;
374 gAlice->GetEvent(i);
375 Int_t nPart = gAlice->GetNtrack();
376 for (Int_t iPart = 0; iPart < nPart; iPart++) {
377 //cout<<"Particlenr.: "<<iPart<<endl;
378 p = (TParticle*)gAlice->Particle(iPart);
379 if (p->GetPdgCode()==421){
380 if(fabs(p->Eta())<0.9) h1->Fill(p.Pt());
381 h4->Fill(p.Eta());
382 nd0++;
383 }
384 if (p->GetPdgCode()==-321){
385 nkminus++;
386 }
387 if (p->GetPdgCode()==211){
388 npipluss++;
389 }
390 }
391 }
392 delete gAlice;
393}