3 #if !defined(__CINT__) || defined(__MAKECINT__)
4 //-- --- standard headers-------------
7 //--------Root headers ---------------
11 #include <TStopwatch.h>
16 #include <TParticle.h>
18 //----- AliRoot headers ---------------
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"
32 //-------------------------------------
43 const Double_t kBz = 0.4;
46 double primaryvertex[3]={0.,0.,0.};
52 const Double_t kPtCut = 0.5; // 0.5 GeV/c
53 const Double_t kd0Cut = 50.; // 50 micron
54 const Double_t kd0CutHigh = 400.; // 200 micron
57 //cuts for combined tracks
58 const 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
63 0, // 0.8 cuts[5] = costhetastar
64 0.5}; // 0.5 cuts[6] = ptchild
65 //cut for distance of closest aprach
66 double cutDCA=0.05; //0.05
68 // this function applies single track cuts
69 Bool_t TrkCuts(const AliITStrackV2& trk);
71 // this function creates TObjArrays with positive and negative tracks
72 void SelectTracks(TTree& itsTree,
73 TObjArray& trksP,Int_t* itsEntryP,Int_t& nTrksP,
74 TObjArray& trksN,Int_t* itsEntryN,Int_t& nTrksN);
76 //void GetPrimaryVertex(int i,Char_t* path="./");
78 //void PtD0(Char_t* path="./");
80 Int_t iTrkP,iTrkN,itsEntries;
81 Char_t trksName[100],refsName[100];
82 Int_t nTrksP=0,nTrksN=0;
89 Int_t nTotEv=0,nD0recSgn=0,nD0recBkgS=0,nD0recBkg=0,nD0rec1ev=0;
90 Int_t pdg[2],mum[2],mumlab[2];
93 void RunD0offline(Char_t* path="./",bool h=false,bool PtD0=false) {
95 AliKalmanTrack::SetConvConst(100/0.299792458/kBz);
97 // Open file with ITS tracks
98 Char_t fnameTrack[1024];
99 //sprintf(fnameTrack,"%s/AliITStracksV2.root",path);
100 sprintf(fnameTrack,"%s/AliITStracksV2.root",path);
101 TFile* itstrks = TFile::Open(fnameTrack);
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);
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
115 -1.0, // min. allowed cosine of V0's pointing angle
116 0.0, // min. radius of the fiducial volume
117 2.9};// max. radius of the fiducial volume
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);
125 sprintf(falice,"%s/galice.root",path);
126 TFile *f = new TFile(falice);
127 gAlice=(AliRun*)f->Get("gAlice");
128 int nEvent=gAlice->GetEventsPerRun();
130 cout<<"#Events: "<<nEvent<<endl;
133 const Char_t *name="AliD0offline";
134 cerr<<'\n'<<name<<"...\n";
135 gBenchmark->Start(name);
137 for(ev=0;ev<nEvent;ev++) {
139 cout<<"\n Event: "<<ev<<endl;
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);
148 // tree from reference file
150 sprintf(refsName,"Tree_Ref_%d",ev);
151 TTree *refTree=(TTree*)itsrefs->Get(refsName);
152 refTree->SetBranchAddress("rectracks",&rectrk);
154 //Getting primary Vertex
155 GetPrimaryVertex(ev,path);
157 // count the total number of events
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);
168 cout<<"#pos: "<<nTrksP<<endl;
169 cout<<"#neg: "<<nTrksN<<endl;
171 // create the AliV0vertexer object
172 AliV0vertexer *vertexer2 = new AliV0vertexer(vtxcuts);
174 AliD0Trigger * D0 = new AliD0Trigger(cuts,kBz,primaryvertex);
176 double ptP,alphaP,phiP,ptN,alphaN,phiN,dca;
178 for(iTrkP=0; iTrkP<nTrksP; iTrkP++) {
179 postrack = (AliITStrackV2*)trksP.At(iTrkP);
181 // get info on tracks PDG and mothers PDG from reference file
182 refTree->GetEvent(itsEntryP[iTrkP]);
184 mum[0] = rectrk.mumpdg;
185 mumlab[0] = rectrk.mumlab;
187 for(iTrkN=0; iTrkN<nTrksN; iTrkN++) {
188 negtrack = (AliITStrackV2*)trksN.At(iTrkN);
190 // get info on tracks PDG and mothers PDG from reference file
191 refTree->GetEvent(itsEntryN[iTrkN]);
193 mum[1] = rectrk.mumpdg;
194 mumlab[1] = rectrk.mumlab;
196 D0->SetTracks(postrack,negtrack);
198 // ----------- DCA MINIMIZATION ------------------
200 // find the DCA and propagate the tracks to the DCA
201 double dca = vertexer2->PropagateToDCA(negtrack,postrack);
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]);
209 //if(D0->FindV0offline(v2) && D0->d0d0()){
212 D0->FindMomentaOffline();
213 if(D0->FindInvMass() && D0->CosThetaStar() && D0->PointingAngle() && D0->pTchild()){
215 event[index]=ev; index++;
221 if(mumlab[0]==mumlab[1] && TMath::Abs(mum[0])==421 && TMath::Abs(mum[1])==421) {
225 if(TMath::Abs(mum[0])==421 || TMath::Abs(mum[0])==411 ||
226 TMath::Abs(mum[1])==421 || TMath::Abs(mum[1])==411) {
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
252 +++ Bkg: %d\n+++\n",nD0recSgn,nD0recBkgS,nD0recBkg);
254 gBenchmark->Stop(name);
255 gBenchmark->Show(name);
259 TCanvas *c = new TCanvas("c","c",700,1000);
263 pt = new TPaveText(7.3,0.4,11,1.8,"br");
265 pt->SetBorderSize(1);
266 pt->AddText("Cuts:");
268 sprintf(st,"First Pt: %g",kPtCut);
270 sprintf(st,"d0 low: %g",kd0Cut);
272 sprintf(st,"d0 high: %g",kd0CutHigh);
274 sprintf(st,"V0 low: %g",cuts[0]);
276 sprintf(st,"V0 high: %g",cuts[1]);
278 sprintf(st,"InvMass Diff: %g",cuts[2]);
280 sprintf(st,"cosPointingAngle: %g",cuts[3]);
282 sprintf(st,"d0d0: %g",cuts[4]);
284 sprintf(st,"cosThetaStar: %g",cuts[5]);
286 sprintf(st,"PtChild: %g",cuts[6]);
288 sprintf(st,"DCA: %g",cutDCA);
298 TCanvas *c = new TCanvas("c","c",1000,700);
302 pt = new TPaveText(7.3,0.4,11,1.8,"br");
304 pt->SetBorderSize(1);
305 pt->AddText("Cuts:");
307 sprintf(st,"First Pt: %g",kPtCut);
309 sprintf(st,"d0 low: %g",kd0Cut);
311 sprintf(st,"d0 high: %g",kd0CutHigh);
313 sprintf(st,"V0 low: %g",cuts[0]);
315 sprintf(st,"V0 high: %g",cuts[1]);
317 sprintf(st,"InvMass Diff: %g",cuts[2]);
319 sprintf(st,"cosPointAng: %g",cuts[3]);
321 sprintf(st,"d0d0: %g",cuts[4]);
323 sprintf(st,"cosTheta*: %g",cuts[5]);
325 sprintf(st,"PtChild: %g",cuts[6]);
327 sprintf(st,"DCA: %g",cutDCA);
342 Char_t outName[1024];
343 sprintf(outName,"%s/ReconstructedD0.root",path);
344 TFile* outroot = new TFile(outName,"recreate");
352 Char_t outName[1024];
353 sprintf(outName,"%s/ReconstructedD0.root",path);
354 TFile* outroot = new TFile(outName,"recreate");
362 Char_t foutName[1024];
363 sprintf(foutName,"%s/Cuts",path);
364 ofstream fout(foutName);
366 sprintf(st2,"First Pt: %g",kPtCut);
368 sprintf(st2,"d0 low: %g",kd0Cut);
370 sprintf(st2,"d0 high: %g",kd0CutHigh);
372 sprintf(st2,"V0 low: %g",cuts[0]);
374 sprintf(st2,"V0 high: %g",cuts[1]);
376 sprintf(st2,"InvMass Diff: %g",cuts[2]);
378 sprintf(st2,"cosPointAng: %g",cuts[3]);
380 sprintf(st2,"d0d0: %g",cuts[4]);
382 sprintf(st2,"cosTheta*: %g",cuts[5]);
384 sprintf(st2,"PtChild: %g",cuts[6]);
386 sprintf(st2,"DCA: %g",cutDCA);
391 sprintf(fName,"%s/Events",path);
392 ofstream fevent(fName);
393 for(int i=0;i<nD0;i++){fevent<<event[i]<<endl;}
396 //___________________________________________________________________________
397 void SelectTracks(TTree& itsTree,
398 TObjArray& trksP,Int_t* itsEntryP,Int_t& nTrksP,
399 TObjArray& trksN,Int_t* itsEntryN,Int_t& nTrksN) {
401 // this function creates two TObjArrays with positive and negative tracks
406 Int_t entr = (Int_t)itsTree.GetEntries();
408 // trasfer tracks from tree to arrays
409 for(Int_t i=0; i<entr; i++) {
411 AliITStrackV2 *itstrack = new AliITStrackV2;
412 itsTree.SetBranchAddress("tracks",&itstrack);
416 // single track selection
417 if(!TrkCuts(*itstrack)) { delete itstrack; continue; }
419 if(itstrack->Get1Pt()>0.) { // negative track
420 trksN.AddLast(itstrack);
421 itsEntryN[nTrksN] = i;
423 } else { // positive track
424 trksP.AddLast(itstrack);
425 itsEntryP[nTrksP] = i;
433 //____________________________________________________________________________
434 Bool_t TrkCuts(const AliITStrackV2& trk) {
436 // this function tells if track passes some kinematical cuts
438 if(TMath::Abs(1./trk.Get1Pt()) < kPtCut) return kFALSE;
439 if(TMath::Abs(10000.*trk.GetD(primaryvertex[0],primaryvertex[1])) < kd0Cut) return kFALSE;
440 //if(TMath::Abs(10000.*trk.GetD(primaryvertex[0],primaryvertex[1])) > kd0CutHigh) return kFALSE;
444 //____________________________________________________________________________
445 void GetPrimaryVertex(int i,Char_t* path="./") {
450 sprintf(falice,"%s/galice.root",path);
451 TFile * galice = new TFile(falice);
458 sprintf(vname,"Vertex_%d",event);
461 AliHeader * header = 0;
463 TTree * treeE = (TTree*)gDirectory->Get("TE");
464 treeE->SetBranchAddress("Header",&header);
465 treeE->GetEntry(event);
466 AliGenEventHeader* genHeader = header->GenEventHeader();
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];
474 printf("Can't find Header");
479 //____________________________________________________________________________
480 PtD0(Char_t* path="./",TH1F * h1,TH1F * h4){
483 sprintf(falice,"%s/galice.root",path);
484 TFile *f = new TFile(falice);
485 gAlice=(AliRun*)f->Get("gAlice");
491 int nEvent=gAlice->GetEventsPerRun();
493 for (Int_t i = 0; i <nEvent; i++) {
494 cout<<"Event: "<<i<<endl;
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());
505 if (p->GetPdgCode()==-321){
508 if (p->GetPdgCode()==211){