git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Added Gautes changes from Bergen.
[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>
68772666 6#include <fstream.h>
0bd0c1ef 7//--------Root headers ---------------
8#include <TSystem.h>
9#include <TFile.h>
10#include <TString.h>
11#include <TStopwatch.h>
12#include <TObject.h>
13#include <TVector3.h>
14#include <TTree.h>
68772666 15#include <TObjArray>
0bd0c1ef 16#include <TParticle.h>
17#include <TArray.h>
18//----- AliRoot headers ---------------
19#include "alles.h"
20#include "AliRun.h"
21#include "AliKalmanTrack.h"
22#include "AliITStrackV2.h"
23#include "AliHeader.h"
24#include "AliGenEventHeader.h"
25#include "AliV0vertex.h"
26#include "AliV0vertexer.h"
27#include "AliITSVertex.h"
28#include "AliITSVertexer.h"
29#include "AliITSVertexerTracks.h"
30#include "AliD0Trigger.h"
31#endif
32//-------------------------------------
68772666 33typedef struct {
34 Int_t lab;
35 Int_t pdg;
36 Int_t mumlab;
37 Int_t mumpdg;
38 Float_t Vx,Vy,Vz;
39 Float_t Px,Py,Pz;
40} RECTRACK;
41
0bd0c1ef 42// field (T)
43const Double_t kBz = 0.4;
44
45// primary vertex
de3c3890 46double primaryvertex[3]={0.,0.,0.};
0bd0c1ef 47
48//sec. vertex
49double v2[3]={0,0,0};
50
51// sigle track cuts
de3c3890 52const Double_t kPtCut = 0.5; // 0.5 GeV/c
53const Double_t kd0Cut = 50.; // 50 micron
54const Double_t kd0CutHigh = 400.; // 200 micron
0bd0c1ef 55
56
57//cuts for combined tracks
de3c3890 58const Double_t cuts[7] = {0.005, // 0.005 cuts[0] = lowest V0 cut (cm)
59 0.800, // 0.015 cuts[1] = highest V0 cut (cm)
60 0.012, // 0.012 cuts[2] = inv. mass cut (diferense) (Gev/c)
61 0.8, // 0.8 cuts[3] = min. cosine value for pointing angle
62 -5000, // -5000 cuts[4] = d0d0
68772666 63 0, // 0.8 cuts[5] = costhetastar
de3c3890 64 0.5}; // 0.5 cuts[6] = ptchild
0bd0c1ef 65//cut for distance of closest aprach
de3c3890 66double cutDCA=0.05; //0.05
0bd0c1ef 67
68// this function applies single track cuts
69Bool_t TrkCuts(const AliITStrackV2& trk);
70
71// this function creates TObjArrays with positive and negative tracks
72void SelectTracks(TTree& itsTree,
73 TObjArray& trksP,Int_t* itsEntryP,Int_t& nTrksP,
74 TObjArray& trksN,Int_t* itsEntryN,Int_t& nTrksN);
75
76//void GetPrimaryVertex(int i,Char_t* path="./");
77
de3c3890 78//void PtD0(Char_t* path="./");
79
0bd0c1ef 80Int_t iTrkP,iTrkN,itsEntries;
68772666 81Char_t trksName[100],refsName[100];
0bd0c1ef 82Int_t nTrksP=0,nTrksN=0;
83Int_t nD0=0;
84int ev=0;
85double mom[6];
de3c3890 86int event[10000];
87int index=0;
68772666 88Bool_t isSignal;
89Int_t nTotEv=0,nD0recSgn=0,nD0recBkgS=0,nD0recBkg=0,nD0rec1ev=0;
90Int_t pdg[2],mum[2],mumlab[2];
91RECTRACK rectrk;
0bd0c1ef 92
de3c3890 93void RunD0offline(Char_t* path="./",bool h=false,bool PtD0=false) {
68772666 94
0bd0c1ef 95 AliKalmanTrack::SetConvConst(100/0.299792458/kBz);
68772666 96
0bd0c1ef 97 // Open file with ITS tracks
98 Char_t fnameTrack[1024];
de3c3890 99 //sprintf(fnameTrack,"%s/AliITStracksV2.root",path);
0bd0c1ef 100 sprintf(fnameTrack,"%s/AliITStracksV2.root",path);
101 TFile* itstrks = TFile::Open(fnameTrack);
68772666 102
103 Char_t refFile[1024];
104 //sprintf(fnameTrack,"%s/ITStracksRefFile.root",path);
105 sprintf(refFile,"%s/ITStracksRefFile.root",path);
106 TFile* itsrefs = TFile::Open(refFile);
107
0bd0c1ef 108
0bd0c1ef 109 //the offline stuff
110 // define the cuts for vertexing
111 Double_t vtxcuts[]={33., // max. allowed chi2
112 0.0, // min. allowed negative daughter's impact param
113 0.0, // min. allowed positive daughter's impact param
114 1.0, // max. allowed DCA between the daughter tracks
de3c3890 115 -1.0, // min. allowed cosine of V0's pointing angle
0bd0c1ef 116 0.0, // min. radius of the fiducial volume
117 2.9};// max. radius of the fiducial volume
68772666 118
de3c3890 119 TH1F *h1 = new TH1F("h1","Transvers momentun of reconstructed D0",100,0,10);
120 TH1F *h2 = new TH1F("h2","Transvers momentun of D0 with |Eta|<0.9",100,0,10);
121 TH1F *h3 = new TH1F("h3","Eta reconstructed of D0",100,-5,5);
122 TH1F *h4 = new TH1F("h4","Eta of D0",100,-5,5);
68772666 123
124 Char_t falice[1024];
125 sprintf(falice,"%s/galice.root",path);
126 TFile *f = new TFile(falice);
127 gAlice=(AliRun*)f->Get("gAlice");
128 int nEvent=gAlice->GetEventsPerRun();
129 //int nEvent=20;
130 cout<<"#Events: "<<nEvent<<endl;
131 delete gAlice;
de3c3890 132
68772666 133 const Char_t *name="AliD0offline";
134 cerr<<'\n'<<name<<"...\n";
135 gBenchmark->Start(name);
136
de3c3890 137 for(ev=0;ev<nEvent;ev++) {
138
139 cout<<"\n Event: "<<ev<<endl;
68772666 140
de3c3890 141 // tracks from ITS
142 sprintf(trksName,"TreeT_ITS_%d",ev);
143 TTree *itsTree=(TTree*)itstrks->Get(trksName);
144 if(!itsTree) continue;
145 itsEntries = (Int_t)itsTree->GetEntries();
146 printf("+++\n+++ Number of tracks in ITS: %d\n+++\n\n",itsEntries);
147
68772666 148 // tree from reference file
149
150 sprintf(refsName,"Tree_Ref_%d",ev);
151 TTree *refTree=(TTree*)itsrefs->Get(refsName);
152 refTree->SetBranchAddress("rectracks",&rectrk);
153
de3c3890 154 //Getting primary Vertex
155 GetPrimaryVertex(ev,path);
156
68772666 157 // count the total number of events
158 nTotEv++;
159
de3c3890 160 // call function which applies sigle track selection and
161 // separetes positives and negatives
162 TObjArray trksP(itsEntries/2);
163 Int_t *itsEntryP = new Int_t[itsEntries];
164 TObjArray trksN(itsEntries/2);
165 Int_t *itsEntryN = new Int_t[itsEntries];
166 SelectTracks(*itsTree,trksP,itsEntryP,nTrksP,trksN,itsEntryN,nTrksN);
167
168 cout<<"#pos: "<<nTrksP<<endl;
169 cout<<"#neg: "<<nTrksN<<endl;
170
171 // create the AliV0vertexer object
172 AliV0vertexer *vertexer2 = new AliV0vertexer(vtxcuts);
173
174 AliD0Trigger * D0 = new AliD0Trigger(cuts,kBz,primaryvertex);
175
176 double ptP,alphaP,phiP,ptN,alphaN,phiN,dca;
177
178 for(iTrkP=0; iTrkP<nTrksP; iTrkP++) {
179 postrack = (AliITStrackV2*)trksP.At(iTrkP);
68772666 180
181 // get info on tracks PDG and mothers PDG from reference file
182 refTree->GetEvent(itsEntryP[iTrkP]);
183 pdg[0] = rectrk.pdg;
184 mum[0] = rectrk.mumpdg;
185 mumlab[0] = rectrk.mumlab;
186
de3c3890 187 for(iTrkN=0; iTrkN<nTrksN; iTrkN++) {
188 negtrack = (AliITStrackV2*)trksN.At(iTrkN);
68772666 189
190 // get info on tracks PDG and mothers PDG from reference file
191 refTree->GetEvent(itsEntryN[iTrkN]);
192 pdg[1] = rectrk.pdg;
193 mum[1] = rectrk.mumpdg;
194 mumlab[1] = rectrk.mumlab;
195
de3c3890 196 D0->SetTracks(postrack,negtrack);
68772666 197
de3c3890 198 // ----------- DCA MINIMIZATION ------------------
199 //
200 // find the DCA and propagate the tracks to the DCA
201 double dca = vertexer2->PropagateToDCA(negtrack,postrack);
202
203 if(dca<cutDCA){
204 // define the AliV0vertex object
205 AliV0vertex *vertex2 = new AliV0vertex(*negtrack,*postrack);
206 // get position of the vertex
207 vertex2->GetXYZ(v2[0],v2[1],v2[2]);
208 delete vertex2;
209 //if(D0->FindV0offline(v2) && D0->d0d0()){
210 if(D0->d0d0()){
211 D0->SetV0(v2);
212 D0->FindMomentaOffline();
213 if(D0->FindInvMass() && D0->CosThetaStar() && D0->PointingAngle() && D0->pTchild()){
214 nD0++;
215 event[index]=ev; index++;
216 if(h){
217 h1->Fill(D0->Pt());
218 h3->Fill(D0->Eta());
0bd0c1ef 219 }
68772666 220
221 if(mumlab[0]==mumlab[1] && TMath::Abs(mum[0])==421 && TMath::Abs(mum[1])==421) {
222 nD0recSgn++;
223 }
224 else {
225 if(TMath::Abs(mum[0])==421 || TMath::Abs(mum[0])==411 ||
226 TMath::Abs(mum[1])==421 || TMath::Abs(mum[1])==411) {
227 nD0recBkgS++;
228 } else {
229 nD0recBkg++;
230 }
231 }
de3c3890 232 }
233 }
0bd0c1ef 234 }
235 }
de3c3890 236 }
68772666 237
238 //delete D0;
239 //delete itstrks;
240 //delete itsTree;
241 //delete trksP;
242 //delete itsEntryP;
243 //delete trksN;
244 //delete itsEntryN;
de3c3890 245 }
68772666 246
247 cout<<"\nMy #D0: "<<nD0<<"\n"<<endl;
248 printf("\n+++\n+++ Total number of events: %d\n+++\n",nTotEv);
249 printf("\n+++\n+++ Total number of D0 candidates:\n
250 +++ Sgn: %d\n
251 +++ BkgS: %d\n
252 +++ Bkg: %d\n+++\n",nD0recSgn,nD0recBkgS,nD0recBkg);
253
0bd0c1ef 254 gBenchmark->Stop(name);
255 gBenchmark->Show(name);
68772666 256
de3c3890 257 if(h){
258 if(!PtD0){
259 TCanvas *c = new TCanvas("c","c",700,1000);
260 c->Divide(1,2);
261 c->cd(1);
262 h1->Draw();
263 pt = new TPaveText(7.3,0.4,11,1.8,"br");
264 pt->SetFillColor(0);
265 pt->SetBorderSize(1);
266 pt->AddText("Cuts:");
267 Char_t st[1024];
268 sprintf(st,"First Pt: %g",kPtCut);
269 pt->AddText(st);
270 sprintf(st,"d0 low: %g",kd0Cut);
271 pt->AddText(st);
272 sprintf(st,"d0 high: %g",kd0CutHigh);
273 pt->AddText(st);
274 sprintf(st,"V0 low: %g",cuts[0]);
275 pt->AddText(st);
276 sprintf(st,"V0 high: %g",cuts[1]);
277 pt->AddText(st);
278 sprintf(st,"InvMass Diff: %g",cuts[2]);
279 pt->AddText(st);
280 sprintf(st,"cosPointingAngle: %g",cuts[3]);
281 pt->AddText(st);
282 sprintf(st,"d0d0: %g",cuts[4]);
283 pt->AddText(st);
284 sprintf(st,"cosThetaStar: %g",cuts[5]);
285 pt->AddText(st);
286 sprintf(st,"PtChild: %g",cuts[6]);
287 pt->AddText(st);
288 sprintf(st,"DCA: %g",cutDCA);
289 pt->AddText(st);
290 pt->Draw();
291 c_1->Modified();
292 c->cd();
293 c->cd(2);
294 h3->Draw();
295 }
296 if(PtD0){
297 PtD0(path,h2,h4);
298 TCanvas *c = new TCanvas("c","c",1000,700);
299 c->Divide(2,2);
300 c->cd(1);
301 h1->Draw();
302 pt = new TPaveText(7.3,0.4,11,1.8,"br");
303 pt->SetFillColor(0);
304 pt->SetBorderSize(1);
305 pt->AddText("Cuts:");
306 Char_t st[1024];
307 sprintf(st,"First Pt: %g",kPtCut);
308 pt->AddText(st);
309 sprintf(st,"d0 low: %g",kd0Cut);
310 pt->AddText(st);
311 sprintf(st,"d0 high: %g",kd0CutHigh);
312 pt->AddText(st);
313 sprintf(st,"V0 low: %g",cuts[0]);
314 pt->AddText(st);
315 sprintf(st,"V0 high: %g",cuts[1]);
316 pt->AddText(st);
317 sprintf(st,"InvMass Diff: %g",cuts[2]);
318 pt->AddText(st);
68772666 319 sprintf(st,"cosPointAng: %g",cuts[3]);
de3c3890 320 pt->AddText(st);
321 sprintf(st,"d0d0: %g",cuts[4]);
322 pt->AddText(st);
68772666 323 sprintf(st,"cosTheta*: %g",cuts[5]);
de3c3890 324 pt->AddText(st);
325 sprintf(st,"PtChild: %g",cuts[6]);
326 pt->AddText(st);
327 sprintf(st,"DCA: %g",cutDCA);
328 pt->AddText(st);
329 pt->Draw();
330 c_1->Modified();
331 c->cd();
332 c->cd(2);
333 h3->Draw();
334 c->cd(3);
335 h2->Draw();
336 c->cd(4);
337 h4->Draw();
68772666 338 }
339 }
340 if(h){
341 if(!PtD0){
342 Char_t outName[1024];
343 sprintf(outName,"%s/ReconstructedD0.root",path);
344 TFile* outroot = new TFile(outName,"recreate");
345 h1->Write();
346 h3->Write();
347 outroot->Close();
348 delete outroot;
de3c3890 349 }
68772666 350 }
351 if(PtD0){
352 Char_t outName[1024];
353 sprintf(outName,"%s/ReconstructedD0.root",path);
354 TFile* outroot = new TFile(outName,"recreate");
355 h1->Write();
356 h2->Write();
357 h3->Write();
358 h4->Write();
359 outroot->Close();
360 delete outroot;
361 }
362 Char_t foutName[1024];
363 sprintf(foutName,"%s/Cuts",path);
364 ofstream fout(foutName);
365 Char_t st2[1024];
366 sprintf(st2,"First Pt: %g",kPtCut);
367 fout<<st2<<endl;
368 sprintf(st2,"d0 low: %g",kd0Cut);
369 fout<<st2<<endl;
370 sprintf(st2,"d0 high: %g",kd0CutHigh);
371 fout<<st2<<endl;
372 sprintf(st2,"V0 low: %g",cuts[0]);
373 fout<<st2<<endl;
374 sprintf(st2,"V0 high: %g",cuts[1]);
375 fout<<st2<<endl;
376 sprintf(st2,"InvMass Diff: %g",cuts[2]);
377 fout<<st2<<endl;
378 sprintf(st2,"cosPointAng: %g",cuts[3]);
379 fout<<st2<<endl;
380 sprintf(st2,"d0d0: %g",cuts[4]);
381 fout<<st2<<endl;
382 sprintf(st2,"cosTheta*: %g",cuts[5]);
383 fout<<st2<<endl;
384 sprintf(st2,"PtChild: %g",cuts[6]);
385 fout<<st2<<endl;
386 sprintf(st2,"DCA: %g",cutDCA);
387 fout<<st2<<endl;
388 fout.close();
389
390 Char_t fName[1024];
391 sprintf(fName,"%s/Events",path);
392 ofstream fevent(fName);
393 for(int i=0;i<nD0;i++){fevent<<event[i]<<endl;}
394 fevent.close();
0bd0c1ef 395}
396//___________________________________________________________________________
397void SelectTracks(TTree& itsTree,
398 TObjArray& trksP,Int_t* itsEntryP,Int_t& nTrksP,
399 TObjArray& trksN,Int_t* itsEntryN,Int_t& nTrksN) {
400//
401// this function creates two TObjArrays with positive and negative tracks
402//
403 nTrksP=0,nTrksN=0;
404
405
406 Int_t entr = (Int_t)itsTree.GetEntries();
407
408 // trasfer tracks from tree to arrays
409 for(Int_t i=0; i<entr; i++) {
410
411 AliITStrackV2 *itstrack = new AliITStrackV2;
412 itsTree.SetBranchAddress("tracks",&itstrack);
413
414 itsTree.GetEvent(i);
415
416 // single track selection
417 if(!TrkCuts(*itstrack)) { delete itstrack; continue; }
418
419 if(itstrack->Get1Pt()>0.) { // negative track
420 trksN.AddLast(itstrack);
421 itsEntryN[nTrksN] = i;
422 nTrksN++;
423 } else { // positive track
424 trksP.AddLast(itstrack);
425 itsEntryP[nTrksP] = i;
426 nTrksP++;
427 }
428
429 }
430
431 return;
432}
433//____________________________________________________________________________
434Bool_t TrkCuts(const AliITStrackV2& trk) {
435//
436// this function tells if track passes some kinematical cuts
437//
438 if(TMath::Abs(1./trk.Get1Pt()) < kPtCut) return kFALSE;
439 if(TMath::Abs(10000.*trk.GetD(primaryvertex[0],primaryvertex[1])) < kd0Cut) return kFALSE;
de3c3890 440 //if(TMath::Abs(10000.*trk.GetD(primaryvertex[0],primaryvertex[1])) > kd0CutHigh) return kFALSE;
0bd0c1ef 441
442 return kTRUE;
443}
444//____________________________________________________________________________
445void GetPrimaryVertex(int i,Char_t* path="./") {
446
447 int event=i;
448
449 Char_t falice[1024];
450 sprintf(falice,"%s/galice.root",path);
451 TFile * galice = new TFile(falice);
452
453 TDirectory * curdir;
454
455 Char_t vname[20];
456 galice->cd();
457
458 sprintf(vname,"Vertex_%d",event);
459 TArrayF o = 0;
460 o.Set(3);
461 AliHeader * header = 0;
462
463 TTree * treeE = (TTree*)gDirectory->Get("TE");
464 treeE->SetBranchAddress("Header",&header);
465 treeE->GetEntry(event);
466 AliGenEventHeader* genHeader = header->GenEventHeader();
467 if(genHeader){
468 genHeader->PrimaryVertex(o);
469 primaryvertex[0] = (Double_t)o[0];
470 primaryvertex[1] = (Double_t)o[1];
471 primaryvertex[2] = (Double_t)o[2];
472 }
473 else{
474 printf("Can't find Header");
475 }
476 delete header;
477 delete galice;
478}
de3c3890 479//____________________________________________________________________________
480PtD0(Char_t* path="./",TH1F * h1,TH1F * h4){
481
482 Char_t falice[1024];
483 sprintf(falice,"%s/galice.root",path);
484 TFile *f = new TFile(falice);
485 gAlice=(AliRun*)f->Get("gAlice");
486
487 TParticle *p;
488 int nd0=0;
489 int nkminus =0;
490 int npipluss = 0;
491 int nEvent=gAlice->GetEventsPerRun();
68772666 492
de3c3890 493 for (Int_t i = 0; i <nEvent; i++) {
494 cout<<"Event: "<<i<<endl;
495 gAlice->GetEvent(i);
496 Int_t nPart = gAlice->GetNtrack();
497 for (Int_t iPart = 0; iPart < nPart; iPart++) {
498 //cout<<"Particlenr.: "<<iPart<<endl;
499 p = (TParticle*)gAlice->Particle(iPart);
500 if (p->GetPdgCode()==421){
501 if(fabs(p->Eta())<0.9) h1->Fill(p.Pt());
502 h4->Fill(p.Eta());
503 nd0++;
504 }
505 if (p->GetPdgCode()==-321){
506 nkminus++;
507 }
508 if (p->GetPdgCode()==211){
509 npipluss++;
510 }
511 }
512 }
513 delete gAlice;
514}