const Double_t kBz = 0.4;
// primary vertex
-Double_t primaryvertex[3] = {0.,0.,0,};
+double primaryvertex[3]={0.,0.,0.};
//sec. vertex
double v2[3]={0,0,0};
// sigle track cuts
-const Double_t kPtCut = 0.5; // GeV/c
-const Double_t kd0Cut = 50.; // micron
-const Double_t kd0CutHigh = 200.; // micron
+const Double_t kPtCut = 0.5; // 0.5 GeV/c
+const Double_t kd0Cut = 50.; // 50 micron
+const Double_t kd0CutHigh = 400.; // 200 micron
//cuts for combined tracks
-const Double_t cuts[7] = {0.005, // cuts[0] = lowest V0 cut (cm)
- 0.015, // cuts[1] = highest V0 cut (cm)
- 0.05, // cuts[2] = inv. mass cut (diferense) (Gev/c)
- 0.95, // cuts[3] = max cosine value for pointing angle
- -5000, // cuts[4] = d0d0
- 0.8, // cuts[5] = costhetastar
- 0.5}; // cuts[6] = ptchild
+const Double_t cuts[7] = {0.005, // 0.005 cuts[0] = lowest V0 cut (cm)
+ 0.800, // 0.015 cuts[1] = highest V0 cut (cm)
+ 0.012, // 0.012 cuts[2] = inv. mass cut (diferense) (Gev/c)
+ 0.8, // 0.8 cuts[3] = min. cosine value for pointing angle
+ -5000, // -5000 cuts[4] = d0d0
+ 0, // 0.8 cuts[5] = costhetastar
+ 0.5}; // 0.5 cuts[6] = ptchild
//cut for distance of closest aprach
-double cutDCA=0.01;
+double cutDCA=0.05; //0.05
// this function applies single track cuts
Bool_t TrkCuts(const AliITStrackV2& trk);
//void GetPrimaryVertex(int i,Char_t* path="./");
+//void PtD0(Char_t* path="./");
+
Int_t iTrkP,iTrkN,itsEntries;
Char_t trksName[100];
Int_t nTrksP=0,nTrksN=0;
Int_t nD0=0;
int ev=0;
double mom[6];
+int event[10000];
+int index=0;
-void RunD0offline(Int_t evFirst=0,Int_t evLast=1,Char_t* path="./") {
+void RunD0offline(Char_t* path="./",bool h=false,bool PtD0=false) {
+ Char_t falice[1024];
+ sprintf(falice,"%s/galice.root",path);
+ TFile *f = new TFile(falice);
+ gAlice=(AliRun*)f->Get("gAlice");
+ int nEvent=gAlice->GetEventsPerRun();
+ //int nEvent=20;
+ delete gAlice;
const Char_t *name="AliD0offline";
cerr<<'\n'<<name<<"...\n";
+ cout<<"#Events: "<<nEvent<<endl;
gBenchmark->Start(name);
AliKalmanTrack::SetConvConst(100/0.299792458/kBz);
// Open file with ITS tracks
Char_t fnameTrack[1024];
+ //sprintf(fnameTrack,"%s/AliITStracksV2.root",path);
sprintf(fnameTrack,"%s/AliITStracksV2.root",path);
TFile* itstrks = TFile::Open(fnameTrack);
- // tracks from ITS
- sprintf(trksName,"TreeT_ITS_%d",ev);
- TTree *itsTree=(TTree*)itstrks->Get(trksName);
- if(!itsTree) continue;
- itsEntries = (Int_t)itsTree->GetEntries();
- printf("+++\n+++ Number of tracks in ITS: %d\n+++\n\n",itsEntries);
-
- //Getting primary Vertex
- GetPrimaryVertex(0,path);
-
- // call function which applies sigle track selection and
- // separetes positives and negatives
- TObjArray trksP(itsEntries/2);
- Int_t *itsEntryP = new Int_t[itsEntries];
- TObjArray trksN(itsEntries/2);
- Int_t *itsEntryN = new Int_t[itsEntries];
- SelectTracks(*itsTree,trksP,itsEntryP,nTrksP,trksN,itsEntryN,nTrksN);
-
- cout<<"#pos: "<<nTrksP<<endl;
- cout<<"#neg: "<<nTrksN<<endl;
-
//the offline stuff
// define the cuts for vertexing
Double_t vtxcuts[]={33., // max. allowed chi2
0.0, // min. allowed negative daughter's impact param
0.0, // min. allowed positive daughter's impact param
1.0, // max. allowed DCA between the daughter tracks
- -1.0, // min. allowed cosine of V0's pointing angle
+ -1.0, // min. allowed cosine of V0's pointing angle
0.0, // min. radius of the fiducial volume
2.9};// max. radius of the fiducial volume
-
- // create the AliV0vertexer object
- AliV0vertexer *vertexer2 = new AliV0vertexer(vtxcuts);
-
- AliD0Trigger * D0 = new AliD0Trigger(cuts,kBz,primaryvertex);
- double ptP,alphaP,phiP,ptN,alphaN,phiN,dca;
-
- for(iTrkP=0; iTrkP<nTrksP; iTrkP++) {
- postrack = (AliITStrackV2*)trksP.At(iTrkP);
- for(iTrkN=0; iTrkN<nTrksN; iTrkN++) {
- negtrack = (AliITStrackV2*)trksN.At(iTrkN);
- D0.SetTracks(postrack,negtrack);
- //
- // ----------- DCA MINIMIZATION ------------------
- //
- // find the DCA and propagate the tracks to the DCA
- double dca = vertexer2->PropagateToDCA(negtrack,postrack);
-
- if(dca<cutDCA){
- // define the AliV0vertex object
- AliV0vertex *vertex2 = new AliV0vertex(*negtrack,*postrack);
- // get position of the vertex
- vertex2->GetXYZ(v2[0],v2[1],v2[2]);
- delete vertex2;
- if(D0.pTchild()){
- if(D0.d0d0()){
- if(D0.FindV0offline(v2)){
-
- // momenta of the tracks at the vertex
- ptP = 1./TMath::Abs(postrack->Get1Pt());
- alphaP = postrack->GetAlpha();
- phiP = alphaP+TMath::ASin(postrack->GetSnp());
- mom[0] = ptP*TMath::Cos(phiP);
- mom[1] = ptP*TMath::Sin(phiP);
- mom[2] = ptP*postrack->GetTgl();
-
- ptN = 1./TMath::Abs(negtrack->Get1Pt());
- alphaN = negtrack->GetAlpha();
- phiN = alphaN+TMath::ASin(negtrack->GetSnp());
- mom[3] = ptN*TMath::Cos(phiN);
- mom[4] = ptN*TMath::Sin(phiN);
- mom[5] = ptN*negtrack->GetTgl();
-
- D0.SetMomenta(mom);
-
- if(D0.FindInvMass()){
- if(D0.CosThetaStar()){
- if(D0.PointingAngle()){
- nD0++;
- }
- }
+ TH1F *h1 = new TH1F("h1","Transvers momentun of reconstructed D0",100,0,10);
+ TH1F *h2 = new TH1F("h2","Transvers momentun of D0 with |Eta|<0.9",100,0,10);
+ TH1F *h3 = new TH1F("h3","Eta reconstructed of D0",100,-5,5);
+ TH1F *h4 = new TH1F("h4","Eta of D0",100,-5,5);
+
+ for(ev=0;ev<nEvent;ev++) {
+
+ cout<<"\n Event: "<<ev<<endl;
+
+ // tracks from ITS
+ sprintf(trksName,"TreeT_ITS_%d",ev);
+ TTree *itsTree=(TTree*)itstrks->Get(trksName);
+ if(!itsTree) continue;
+ itsEntries = (Int_t)itsTree->GetEntries();
+ printf("+++\n+++ Number of tracks in ITS: %d\n+++\n\n",itsEntries);
+
+ //Getting primary Vertex
+ GetPrimaryVertex(ev,path);
+
+ // call function which applies sigle track selection and
+ // separetes positives and negatives
+ TObjArray trksP(itsEntries/2);
+ Int_t *itsEntryP = new Int_t[itsEntries];
+ TObjArray trksN(itsEntries/2);
+ Int_t *itsEntryN = new Int_t[itsEntries];
+ SelectTracks(*itsTree,trksP,itsEntryP,nTrksP,trksN,itsEntryN,nTrksN);
+
+ cout<<"#pos: "<<nTrksP<<endl;
+ cout<<"#neg: "<<nTrksN<<endl;
+
+ // create the AliV0vertexer object
+ AliV0vertexer *vertexer2 = new AliV0vertexer(vtxcuts);
+
+ AliD0Trigger * D0 = new AliD0Trigger(cuts,kBz,primaryvertex);
+
+ double ptP,alphaP,phiP,ptN,alphaN,phiN,dca;
+
+ for(iTrkP=0; iTrkP<nTrksP; iTrkP++) {
+ postrack = (AliITStrackV2*)trksP.At(iTrkP);
+ for(iTrkN=0; iTrkN<nTrksN; iTrkN++) {
+ negtrack = (AliITStrackV2*)trksN.At(iTrkN);
+ D0->SetTracks(postrack,negtrack);
+ //
+ // ----------- DCA MINIMIZATION ------------------
+ //
+ // find the DCA and propagate the tracks to the DCA
+ double dca = vertexer2->PropagateToDCA(negtrack,postrack);
+
+ if(dca<cutDCA){
+ // define the AliV0vertex object
+ AliV0vertex *vertex2 = new AliV0vertex(*negtrack,*postrack);
+ // get position of the vertex
+ vertex2->GetXYZ(v2[0],v2[1],v2[2]);
+ delete vertex2;
+ //if(D0->FindV0offline(v2) && D0->d0d0()){
+ if(D0->d0d0()){
+ D0->SetV0(v2);
+ D0->FindMomentaOffline();
+ if(D0->FindInvMass() && D0->CosThetaStar() && D0->PointingAngle() && D0->pTchild()){
+ nD0++;
+ event[index]=ev; index++;
+ if(h){
+ h1->Fill(D0->Pt());
+ h3->Fill(D0->Eta());
}
- }
- }
+ }
+ }
}
}
- }
- }
- cout<<"#D0: "<<nD0<<endl;
+ }
+ }
+ cout<<"\n#D0: "<<nD0<<endl;
+ for(int i=0;i<nD0;i++)
+ cout<<event[i]<<endl;
gBenchmark->Stop(name);
gBenchmark->Show(name);
+
+ if(h){
+ if(!PtD0){
+ TCanvas *c = new TCanvas("c","c",700,1000);
+ c->Divide(1,2);
+ c->cd(1);
+ h1->Draw();
+ pt = new TPaveText(7.3,0.4,11,1.8,"br");
+ pt->SetFillColor(0);
+ pt->SetBorderSize(1);
+ pt->AddText("Cuts:");
+ Char_t st[1024];
+ sprintf(st,"First Pt: %g",kPtCut);
+ pt->AddText(st);
+ sprintf(st,"d0 low: %g",kd0Cut);
+ pt->AddText(st);
+ sprintf(st,"d0 high: %g",kd0CutHigh);
+ pt->AddText(st);
+ sprintf(st,"V0 low: %g",cuts[0]);
+ pt->AddText(st);
+ sprintf(st,"V0 high: %g",cuts[1]);
+ pt->AddText(st);
+ sprintf(st,"InvMass Diff: %g",cuts[2]);
+ pt->AddText(st);
+ sprintf(st,"cosPointingAngle: %g",cuts[3]);
+ pt->AddText(st);
+ sprintf(st,"d0d0: %g",cuts[4]);
+ pt->AddText(st);
+ sprintf(st,"cosThetaStar: %g",cuts[5]);
+ pt->AddText(st);
+ sprintf(st,"PtChild: %g",cuts[6]);
+ pt->AddText(st);
+ sprintf(st,"DCA: %g",cutDCA);
+ pt->AddText(st);
+ pt->Draw();
+ c_1->Modified();
+ c->cd();
+ c->cd(2);
+ h3->Draw();
+ }
+ if(PtD0){
+ PtD0(path,h2,h4);
+ TCanvas *c = new TCanvas("c","c",1000,700);
+ c->Divide(2,2);
+ c->cd(1);
+ h1->Draw();
+ pt = new TPaveText(7.3,0.4,11,1.8,"br");
+ pt->SetFillColor(0);
+ pt->SetBorderSize(1);
+ pt->AddText("Cuts:");
+ Char_t st[1024];
+ sprintf(st,"First Pt: %g",kPtCut);
+ pt->AddText(st);
+ sprintf(st,"d0 low: %g",kd0Cut);
+ pt->AddText(st);
+ sprintf(st,"d0 high: %g",kd0CutHigh);
+ pt->AddText(st);
+ sprintf(st,"V0 low: %g",cuts[0]);
+ pt->AddText(st);
+ sprintf(st,"V0 high: %g",cuts[1]);
+ pt->AddText(st);
+ sprintf(st,"InvMass Diff: %g",cuts[2]);
+ pt->AddText(st);
+ sprintf(st,"cosPointingAngle: %g",cuts[3]);
+ pt->AddText(st);
+ sprintf(st,"d0d0: %g",cuts[4]);
+ pt->AddText(st);
+ sprintf(st,"cosThetaStar: %g",cuts[5]);
+ pt->AddText(st);
+ sprintf(st,"PtChild: %g",cuts[6]);
+ pt->AddText(st);
+ sprintf(st,"DCA: %g",cutDCA);
+ pt->AddText(st);
+ pt->Draw();
+ c_1->Modified();
+ c->cd();
+ c->cd(2);
+ h3->Draw();
+ c->cd(3);
+ h2->Draw();
+ c->cd(4);
+ h4->Draw();
+ }
+ //TFile *outfile = new TFile("results.root","recreate");
+ //h1->Write();
+ //outfile->Close();
+
+ }
}
//___________________________________________________________________________
void SelectTracks(TTree& itsTree,
//
if(TMath::Abs(1./trk.Get1Pt()) < kPtCut) return kFALSE;
if(TMath::Abs(10000.*trk.GetD(primaryvertex[0],primaryvertex[1])) < kd0Cut) return kFALSE;
- if(TMath::Abs(10000.*trk.GetD(primaryvertex[0],primaryvertex[1])) > kd0CutHigh) return kFALSE;
+ //if(TMath::Abs(10000.*trk.GetD(primaryvertex[0],primaryvertex[1])) > kd0CutHigh) return kFALSE;
return kTRUE;
}
delete header;
delete galice;
}
+//____________________________________________________________________________
+PtD0(Char_t* path="./",TH1F * h1,TH1F * h4){
+
+ Char_t falice[1024];
+ sprintf(falice,"%s/galice.root",path);
+ TFile *f = new TFile(falice);
+ gAlice=(AliRun*)f->Get("gAlice");
+
+ TParticle *p;
+ int nd0=0;
+ int nkminus =0;
+ int npipluss = 0;
+ int nEvent=gAlice->GetEventsPerRun();
+
+ for (Int_t i = 0; i <nEvent; i++) {
+ cout<<"Event: "<<i<<endl;
+ gAlice->GetEvent(i);
+ Int_t nPart = gAlice->GetNtrack();
+ for (Int_t iPart = 0; iPart < nPart; iPart++) {
+ //cout<<"Particlenr.: "<<iPart<<endl;
+ p = (TParticle*)gAlice->Particle(iPart);
+ if (p->GetPdgCode()==421){
+ if(fabs(p->Eta())<0.9) h1->Fill(p.Pt());
+ h4->Fill(p.Eta());
+ nd0++;
+ }
+ if (p->GetPdgCode()==-321){
+ nkminus++;
+ }
+ if (p->GetPdgCode()==211){
+ npipluss++;
+ }
+ }
+ }
+ delete gAlice;
+}
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff