From a9a2d814b50c77fd7f63fb379aa977c34132ed00 Mon Sep 17 00:00:00 2001 From: barbera Date: Tue, 29 Jan 2002 07:10:44 +0000 Subject: [PATCH] Updated code for tracking V2 (from Y. Belikov) --- ITS/AliCascadeComparison.C | 382 ++++++++++++++++++++++++++++++++ ITS/AliCascadeFindVertices.C | 42 ++++ ITS/AliCascadeVertex.cxx | 193 +++++++++++++++++ ITS/AliCascadeVertex.h | 60 ++++++ ITS/AliCascadeVertexer.cxx | 251 +++++++++++++++++++++ ITS/AliCascadeVertexer.h | 72 +++++++ ITS/AliITSComparisonV2.C | 2 +- ITS/AliITSFindClustersV2.C | 127 +++-------- ITS/AliITSrecoV2.h | 4 +- ITS/AliITStestV2.C | 14 +- ITS/AliITStrackV2.cxx | 407 +++++++++++------------------------ ITS/AliITStrackV2.h | 28 ++- ITS/AliITStrackerV2.cxx | 280 ++++++++++++------------ ITS/AliITStrackerV2.h | 4 +- ITS/AliV0Comparison.C | 17 +- ITS/AliV0FindVertices.C | 10 +- ITS/AliV0vertexer.cxx | 34 +-- ITS/ITSHitsToFastPoints.C | 13 +- ITS/ITSLinkDef.h | 2 + ITS/Makefile | 3 +- 20 files changed, 1381 insertions(+), 564 deletions(-) create mode 100644 ITS/AliCascadeComparison.C create mode 100644 ITS/AliCascadeFindVertices.C create mode 100644 ITS/AliCascadeVertex.cxx create mode 100644 ITS/AliCascadeVertex.h create mode 100644 ITS/AliCascadeVertexer.cxx create mode 100644 ITS/AliCascadeVertexer.h diff --git a/ITS/AliCascadeComparison.C b/ITS/AliCascadeComparison.C new file mode 100644 index 00000000000..92758ff743f --- /dev/null +++ b/ITS/AliCascadeComparison.C @@ -0,0 +1,382 @@ +/**************************************************************************** + * Very important, delicate and rather obscure macro. * + * * + * Creates list of "findable" cascades, * + * calculates efficiency, resolutions etc. * + * * + * Origin: Christian Kuhn, IReS, Strasbourg, christian.kuhn@ires.in2p3.fr * + ****************************************************************************/ + +#ifndef __CINT__ + #include + #include + + #include "TH1.h" + #include "TFile.h" + #include "TTree.h" + #include "TObjArray.h" + #include "TStyle.h" + #include "TCanvas.h" + #include "TLine.h" + #include "TText.h" + #include "TParticle.h" + #include "TStopwatch.h" + + #include "AliRun.h" + #include "AliPDG.h" + #include "AliCascadeVertex.h" +#endif + +struct GoodCascade { + Int_t nlab,plab; // V0's daughter labels + Int_t blab; // Bachelor label + Int_t code; + Float_t px,py,pz; + Float_t x,y,z; +}; +Int_t good_cascades(GoodCascade *gt, Int_t max); + +Int_t AliCascadeComparison(Int_t code=3312) { + //code= 3312; //kXiMinus + //code=-3312; //kXiPlusBar + //code= 3334; //kOmegaMinus + //code=-3334; //kOmegaPlusBar + + cerr<<"Doing comparison...\n"; + + const Double_t cascadeWindow=0.05, cascadeWidth=0.015; + Double_t cascadeMass=0.5; + switch (code) { + case kXiMinus: + case kXiPlusBar: cascadeMass=1.32131; break; + case kOmegaMinus: + case kOmegaPlusBar: cascadeMass=1.67245; break; + default: cerr<<"Invalid PDG code !\n"; return 1; + } + + TStopwatch timer; + + /*** Load reconstructed cascades ***/ + TFile *cf=TFile::Open("AliCascadeVertices.root"); + if (!cf->IsOpen()){cerr<<"Can't open AliCascadeVertices.root !\n";return 2;} + TObjArray carray(1000); + TTree *cTree=(TTree*)cf->Get("TreeCasc"); + TBranch *branch=cTree->GetBranch("cascades"); + Int_t nentr=(Int_t)cTree->GetEntries(); + for (Int_t i=0; iSetAddress(&iovertex); + cTree->GetEvent(i); + carray.AddLast(iovertex); + } + cf->Close(); + + /*** Check if the file with the "good" cascades exists ***/ + GoodCascade gc[100]; + Int_t ngood=0; + ifstream in("good_cascades"); + if (in) { + cerr<<"Reading good cascades...\n"; + while (in>>gc[ngood].nlab>>gc[ngood].plab>> + gc[ngood].blab>>gc[ngood].code>> + gc[ngood].px>>gc[ngood].py>>gc[ngood].pz>> + gc[ngood].x >>gc[ngood].y >>gc[ngood].z) { + ngood++; + cerr<SetXTitle("(mrad)"); hp->SetFillColor(2); + TH1F *hl=new TH1F("hl","Lambda Resolution",30,-30,30); + hl->SetXTitle("(mrad)"); hl->SetFillColor(1); hl->SetFillStyle(3013); + TH1F *hpt=new TH1F("hpt","Relative Pt Resolution",30,-10.,10.); + hpt->SetXTitle("(%)"); hpt->SetFillColor(2); + + TH1F *hx=new TH1F("hx","Position Resolution (X,Y)",30,-3.,3.); //x res. + hx->SetXTitle("(mm)"); hx->SetFillColor(6); + TH1F *hy=new TH1F("hy","Position Resolution (Y)",30,-3.,3.); //y res + hy->SetXTitle("(mm)"); hy->SetFillColor(1); hy->SetFillStyle(3013); + TH1F *hz=new TH1F("hz","Position Resolution (Z)",30,-3.,3.); //z res. + hz->SetXTitle("(mm)"); hz->SetFillColor(6); + + Double_t pmin=0.2, pmax=4.2; Int_t nchan=20; + TH1F *hgood=new TH1F("hgood","Good Cascades",nchan,pmin,pmax); + TH1F *hfound=new TH1F("hfound","Found Cascades",nchan,pmin,pmax); + TH1F *hfake=new TH1F("hfake","Fake Cascades",nchan,pmin,pmax); + TH1F *hg=new TH1F("hg","Efficiency for Good Cascades",nchan,pmin,pmax); + hg->SetLineColor(4); hg->SetLineWidth(2); + TH1F *hf=new TH1F("hf","Probability of Fake Cascades",nchan,pmin,pmax); + hf->SetFillColor(1); hf->SetFillStyle(3013); hf->SetLineWidth(2); + + Double_t mmin=cascadeMass-cascadeWindow, mmax=cascadeMass+cascadeWindow; + TH1F *cs =new TH1F("v0s","Cascade Effective Mass",40, mmin, mmax); + cs->SetXTitle("(GeV)"); cs->SetFillColor(6); + + Double_t pxg=0.,pyg=0.,ptg=0.; + Int_t nlab=-1, plab=-1, blab=-1; + Int_t i; + for (i=0; iGetV0labels(nlab,plab); + nlab=TMath::Abs(nlab); plab=TMath::Abs(plab); + blab=TMath::Abs(cascade->GetBachelorLabel()); + + cascade->ChangeMassHypothesis(code); + + Double_t mass=cascade->GetEffMass(); + cs->Fill(mass); + + if (TMath::Abs(mass-cascadeMass)>cascadeWidth) continue; + + Int_t j; + for (j=0; jGetPdgCode()) continue; + if (gc[j].nlab == nlab) + if (gc[j].plab == plab) + if (gc[j].blab == blab) break; + } + + Double_t px,py,pz; cascade->GetPxPyPz(px,py,pz); + Double_t pt=TMath::Sqrt(px*px+py*py); + + if (j==ngood) { + hfake->Fill(pt); + cerr<<"Fake cascade: ("<Fill((phi - phig)*1000.); + hl->Fill((lam - lamg)*1000.); + hpt->Fill((1/pt - 1/ptg)/(1/ptg)*100.); + + Double_t x,y,z; cascade->GetXYZ(x,y,z); + hx->Fill((x-gc[j].x)*10); + hy->Fill((y-gc[j].y)*10); + hz->Fill((z-gc[j].z)*10); + + hfound->Fill(ptg); + gc[j].nlab=-1; + + } + for (i=0; iFill(ptg); + nlab=gc[i].nlab; plab=gc[i].plab; blab=gc[i].blab; + if (nlab < 0) continue; + cerr<<"Cascade ("<GetEntries(); + Stat_t nf=hfound->GetEntries(); + + cerr<<"Number of found cascades: "<SetOptStat(111110); + gStyle->SetOptFit(1); + + TCanvas *c1=new TCanvas("c1","",0,0,580,610); + c1->Divide(2,2); + + c1->cd(1); + gPad->SetFillColor(42); gPad->SetFrameFillColor(10); + //hp->Fit("gaus"); + hp->Draw(); + hl->Draw("same"); c1->cd(); + + c1->cd(2); + gPad->SetFillColor(42); gPad->SetFrameFillColor(10); + //hpt->Fit("gaus"); c1->cd(); + hpt->Draw(); c1->cd(); + + c1->cd(3); + gPad->SetFillColor(42); gPad->SetFrameFillColor(10); + //hx->Fit("gaus"); + hx->Draw(); + hy->Draw("same"); c1->cd(); + + c1->cd(4); + gPad->SetFillColor(42); gPad->SetFrameFillColor(10); + //hz->Fit("gaus"); + hz->Draw(); + + + TCanvas *c2=new TCanvas("c2","",600,0,580,610); + c2->Divide(1,2); + + c2->cd(1); + gPad->SetFillColor(42); gPad->SetFrameFillColor(10); + hfound->Sumw2(); hgood->Sumw2(); hfake->Sumw2(); + hg->Divide(hfound,hgood,1,1.,"b"); + hf->Divide(hfake,hgood,1,1.,"b"); + hg->SetMaximum(1.4); + hg->SetYTitle("Cascade reconstruction efficiency"); + hg->SetXTitle("Pt (GeV/c)"); + hg->Draw(); + + TLine *line1 = new TLine(pmin,1.0,pmax,1.0); line1->SetLineStyle(4); + line1->Draw("same"); + TLine *line2 = new TLine(pmin,0.9,pmax,0.9); line2->SetLineStyle(4); + line2->Draw("same"); + + hf->SetFillColor(1); + hf->SetFillStyle(3013); + hf->SetLineColor(2); + hf->SetLineWidth(2); + hf->Draw("histsame"); + TText *text = new TText(0.461176,0.248448,"Fake cascades"); + text->SetTextSize(0.05); + text->Draw(); + text = new TText(0.453919,1.11408,"Good cascades"); + text->SetTextSize(0.05); + text->Draw(); + + + c2->cd(2); + gPad->SetFillColor(42); gPad->SetFrameFillColor(10); + cs->SetXTitle("(GeV/c)"); + cs->Fit("gaus","","",cascadeMass-cascadeWidth,cascadeMass+cascadeWidth); + Double_t max=cs->GetMaximum(); + TLine *line3 = + new TLine(cascadeMass-cascadeWidth,0.,cascadeMass-cascadeWidth,max); + line3->Draw("same"); + TLine *line4 = + new TLine(cascadeMass+cascadeWidth,0.,cascadeMass+cascadeWidth,max); + line4->Draw("same"); + + timer.Stop(); timer.Print(); + + return 0; +} + + + +Int_t good_cascades(GoodCascade *gc, Int_t max) { + Int_t nc=0; + /*** Get information about the cuts ***/ + Double_t r2min=0.9*0.9; + Double_t r2max=2.9*2.9; + + /*** Get labels of the "good" tracks ***/ + Double_t dd; Int_t id, label[15000], ngt=0; + ifstream in("good_tracks_its"); + if (!in) { + cerr<<"Can't open the file good_tracks_its \n"; + return nc; + } + while (in>>label[ngt]>>id>>dd>>dd>>dd>>dd>>dd>>dd) { + ngt++; + if (ngt>=15000) { + cerr<<"Too many good ITS tracks !\n"; + return nc; + } + } + if (!in.eof()) { + cerr<<"Read error (good_tracks_its) !\n"; + return nc; + } + + /*** Get an access to the kinematics ***/ + if (gAlice) {delete gAlice; gAlice=0;} + + TFile *file=TFile::Open("galice.root"); + if (!file->IsOpen()) {cerr<<"Can't open galice.root !\n"; exit(4);} + if (!(gAlice=(AliRun*)file->Get("gAlice"))) { + cerr<<"gAlice has not been found on galice.root !\n"; + exit(5); + } + + Int_t np=gAlice->GetEvent(0); + while (np--) { + cerr<Particle(np); + + /*** only these cascades are "good" ***/ + Int_t code=cp->GetPdgCode(); + if (code!=kXiMinus) if (code!=kXiPlusBar) + if (code!=kOmegaMinus) if (code!=kOmegaPlusBar) continue; + + /*** daughter tracks must be "good" ***/ + Int_t v0lab=cp->GetFirstDaughter(), blab=cp->GetLastDaughter(); + if (v0lab==blab) continue; + if (v0lab<0) continue; + if (blab<0) continue; + + TParticle *p0=gAlice->Particle(v0lab); + TParticle *bp=gAlice->Particle(blab); + if ((p0->GetPdgCode()!=kLambda0) && (p0->GetPdgCode()!=kLambda0Bar)) { + TParticle *p=p0; p0=bp; bp=p; + Int_t i=v0lab; v0lab=blab; blab=i; + if ((p0->GetPdgCode()!=kLambda0) && (p0->GetPdgCode()!=kLambda0Bar)) + continue; + } + + /** is the bachelor "good" ? **/ + Int_t i; + for (i=0; iGetFirstDaughter(), nlab=p0->GetLastDaughter(); + if (nlab==plab) continue; + if (nlab<0) continue; + if (plab<0) continue; + + for (i=0; iVx(), y=bp->Vy(), z=bp->Vz(), r2=x*x+y*y; //bachelor + if (r2r2max) continue; + TParticle *pp=gAlice->Particle(plab); + x=pp->Vx(); y=pp->Vy(); r2=x*x+y*y; //V0 + if (r2r2max) continue; + + if (gAlice->Particle(plab)->GetPDG()->Charge() < 0.) { + i=plab; plab=nlab; nlab=i; + } + + gc[nc].code=code; + gc[nc].plab=plab; gc[nc].nlab=nlab; gc[nc].blab=blab; + gc[nc].px=cp->Px(); gc[nc].py=cp->Py(); gc[nc].pz=cp->Pz(); + gc[nc].x=x; gc[nc].y=y; gc[nc].z=z; + nc++; + + } + + + return nc; +} diff --git a/ITS/AliCascadeFindVertices.C b/ITS/AliCascadeFindVertices.C new file mode 100644 index 00000000000..44eb584e23f --- /dev/null +++ b/ITS/AliCascadeFindVertices.C @@ -0,0 +1,42 @@ +#ifndef __CINT__ + #include + #include "AliCascadeVertexer.h" + #include "TFile.h" + #include "TStopwatch.h" +#endif + +Int_t AliCascadeFindVertices() { + cerr<<"Looking for cascade vertices...\n"; + + TFile *out=TFile::Open("AliCascadeVertices.root","new"); + if (!out->IsOpen()) { + cerr<<"Delete old AliCascadeVertices.root !\n"; return 1; + } + TFile *in=TFile::Open("AliITStracksV2.root"); + if (!in->IsOpen()) {cerr<<"Can't open AliITStracksV2.root !\n"; return 2;} + + TFile *file=TFile::Open("AliV0vertices.root"); + if (!file->IsOpen()) { + cerr<<"Can't open AliV0vertices.root !\n";return 3; + } + Double_t cuts[]={33., // max. allowed chi2 + 0.015,// min. allowed V0 impact parameter + 0.05, // window around the Lambda mass + 0.015,// min. allowed track impact parameter + 0.060,// max. allowed DCA between a V0 and a track + 0.997,// max. allowed cosine of the cascade pointing angle + 0.9, // min. radius of the fiducial volume + 2.9 // max. radius of the fiducial volume + }; + TStopwatch timer; + AliCascadeVertexer *vertexer=new AliCascadeVertexer(cuts); + Int_t rc=vertexer->V0sTracks2CascadeVertices(in,out); + delete vertexer; + timer.Stop(); timer.Print(); + + file->Close(); + in->Close(); + out->Close(); + + return rc; +} diff --git a/ITS/AliCascadeVertex.cxx b/ITS/AliCascadeVertex.cxx new file mode 100644 index 00000000000..d732940bbda --- /dev/null +++ b/ITS/AliCascadeVertex.cxx @@ -0,0 +1,193 @@ +/************************************************************************** + * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * + * * + * Author: The ALICE Off-line Project. * + * Contributors are mentioned in the code where appropriate. * + * * + * Permission to use, copy, modify and distribute this software and its * + * documentation strictly for non-commercial purposes is hereby granted * + * without fee, provided that the above copyright notice appears in all * + * copies and that both the copyright notice and this permission notice * + * appear in the supporting documentation. The authors make no claims * + * about the suitability of this software for any purpose. It is * + * provided "as is" without express or implied warranty. * + **************************************************************************/ + +//------------------------------------------------------------------------- +// Implementation of the cascade vertex class +// +// Origin: Christian Kuhn, IReS, Strasbourg, christian.kuhn@ires.in2p3.fr +//------------------------------------------------------------------------- +#include +#include + +#include "AliCascadeVertex.h" +#include "AliV0vertex.h" +#include "AliITStrackV2.h" + +ClassImp(AliCascadeVertex) + +AliCascadeVertex::AliCascadeVertex() : TObject() { + //-------------------------------------------------------------------- + // Default constructor (Xi-) + //-------------------------------------------------------------------- + fPdgCode=kXiMinus; + fEffMass=1.32131; + fChi2=1.e+33; + fPos[0]=fPos[1]=fPos[2]=0.; + fPosCov[0]=fPosCov[1]=fPosCov[2]=fPosCov[3]=fPosCov[4]=fPosCov[5]=0.; +} + + + +inline Double_t det(Double_t a00, Double_t a01, Double_t a10, Double_t a11){ + // determinant 2x2 + return a00*a11 - a01*a10; +} + +inline Double_t det (Double_t a00,Double_t a01,Double_t a02, + Double_t a10,Double_t a11,Double_t a12, + Double_t a20,Double_t a21,Double_t a22) { + // determinant 3x3 + return + a00*det(a11,a12,a21,a22)-a01*det(a10,a12,a20,a22)+a02*det(a10,a11,a20,a21); +} + + + +AliCascadeVertex::AliCascadeVertex(const AliV0vertex &v,const AliITStrackV2 &t) { + //-------------------------------------------------------------------- + // Main constructor + //-------------------------------------------------------------------- + fPdgCode=kXiMinus; + + fV0lab[0]=v.GetNlabel(); fV0lab[1]=v.GetPlabel(); + fBachLab=t.GetLabel(); + + //Trivial estimation of the vertex parameters + Double_t pt, phi, x, par[5]; + Double_t alpha, cs, sn; + + t.GetExternalParameters(x,par); alpha=t.GetAlpha(); + pt=1./TMath::Abs(par[4]); + phi=TMath::ASin(par[2]) + alpha; + + // momentum of the bachelor track + + Double_t px1=pt*TMath::Cos(phi), py1=pt*TMath::Sin(phi), pz1=pt*par[3]; + + cs=TMath::Cos(alpha); sn=TMath::Sin(alpha); + + Double_t x1=x*cs - par[0]*sn; // position of the bachelor at dca (bachelor,V0) + Double_t y1=x*sn + par[0]*cs; + Double_t z1=par[1]; + + Double_t x2,y2,z2; // position of the V0 + v.GetXYZ(x2,y2,z2); + + Double_t px2,py2,pz2; // momentum of V0 + v.GetPxPyPz(px2,py2,pz2); + + Double_t a2=((x1-x2)*px2+(y1-y2)*py2+(z1-z2)*pz2)/(px2*px2+py2*py2+pz2*pz2); + + Double_t xm=x2+a2*px2; + Double_t ym=y2+a2*py2; + Double_t zm=z2+a2*pz2; + + // position of the cascade decay + + fPos[0]=0.5*(x1+xm); fPos[1]=0.5*(y1+ym); fPos[2]=0.5*(z1+zm); + + + // momenta of the bachelor and the V0 + + fBachMom[0]=px1; fBachMom[1]=py1; fBachMom[2]=pz1; + fV0mom[0]=px2; fV0mom[1]=py2; fV0mom[2]=pz2; + + // invariant mass of the cascade (default is Ximinus) + + Double_t e1=TMath::Sqrt(0.13957*0.13957 + px1*px1 + py1*py1 + pz1*pz1); + Double_t e2=TMath::Sqrt(1.11568*1.11568 + px2*px2 + py2*py2 + pz2*pz2); + + fEffMass=TMath::Sqrt((e1+e2)*(e1+e2)- + (px1+px2)*(px1+px2)-(py1+py2)*(py1+py2)-(pz1+pz2)*(pz1+pz2)); + + fChi2=7.; + +} + +void AliCascadeVertex::ChangeMassHypothesis(Int_t code) { + //-------------------------------------------------------------------- + // This function changes the mass hypothesis for this cascade + //-------------------------------------------------------------------- + + // HOW TO DISTINGUISH BETWEEN A XIMINUS AND A XIPLUS ?????????? + // SAME QUESTION FOR (ANTI-)OMEGA'S (here) ... AND FOR (ANTI-)LAMBDAS (in AliV0vertex) ?? + // -> NEED ADDITIONAL CONDITION ON BACHELOR AND V0 PDGCODE !!!! BUT in the ANALYSIS MACROS !!! + + Double_t massBach, massV0; + + switch (code) { + + case kXiMinus: + massBach=0.13957; massV0=1.11568; break; + case kXiPlusBar: + massBach=0.13957; massV0=1.11568; break; + case kOmegaMinus: + massBach=0.49368; massV0=1.11568; break; + case kOmegaPlusBar: + massBach=0.49368; massV0=1.11568; break; + + default: + cerr<<"AliCascadeVertex::ChangeMassHypothesis: "; + cerr<<"invalide PDG code ! Assuming XiMinus's...\n"; + massBach=0.13957; massV0=1.11568; break; + } + + Double_t px1=fBachMom[0], py1=fBachMom[1], pz1=fBachMom[2]; + Double_t px2=fV0mom[0], py2=fV0mom[1], pz2=fV0mom[2]; + + Double_t e1=TMath::Sqrt(massBach*massBach + px1*px1 + py1*py1 + pz1*pz1); + Double_t e2=TMath::Sqrt(massV0*massV0 + px2*px2 + py2*py2 + pz2*pz2); + fEffMass=TMath::Sqrt((e1+e2)*(e1+e2)- + (px1+px2)*(px1+px2)-(py1+py2)*(py1+py2)-(pz1+pz2)*(pz1+pz2)); + + fPdgCode=code; +} + +void +AliCascadeVertex::GetPxPyPz(Double_t &px, Double_t &py, Double_t &pz) const { + //-------------------------------------------------------------------- + // This function returns the cascade momentum (global) + //-------------------------------------------------------------------- + px=fV0mom[0]+fBachMom[0]; + py=fV0mom[1]+fBachMom[1]; + pz=fV0mom[2]+fBachMom[2]; +} + +void AliCascadeVertex::GetXYZ(Double_t &x, Double_t &y, Double_t &z) const { + //-------------------------------------------------------------------- + // This function returns cascade position (global) + //-------------------------------------------------------------------- + x=fPos[0]; + y=fPos[1]; + z=fPos[2]; +} + +Double_t AliCascadeVertex::GetD(Double_t x0, Double_t y0, Double_t z0) const { + //-------------------------------------------------------------------- + // This function returns the cascade impact parameter + //-------------------------------------------------------------------- + + Double_t x=fPos[0],y=fPos[1],z=fPos[2]; + Double_t px=fV0mom[0]+fBachMom[0]; + Double_t py=fV0mom[1]+fBachMom[1]; + Double_t pz=fV0mom[2]+fBachMom[2]; + + Double_t dx=(y0-y)*pz - (z0-z)*py; + Double_t dy=(x0-x)*pz - (z0-z)*px; + Double_t dz=(x0-x)*py - (y0-y)*px; + Double_t d=TMath::Sqrt((dx*dx+dy*dy+dz*dz)/(px*px+py*py+pz*pz)); + + return d; +} diff --git a/ITS/AliCascadeVertex.h b/ITS/AliCascadeVertex.h new file mode 100644 index 00000000000..6979e738efb --- /dev/null +++ b/ITS/AliCascadeVertex.h @@ -0,0 +1,60 @@ +#ifndef ALICASCADEVERTEX_H +#define ALICASCADEVERTEX_H + +/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * + * See cxx source for full Copyright notice */ + +//------------------------------------------------------------------------- +// Cascade Vertex Class +// +// Origin: Christian Kuhn, IReS, Strasbourg, christian.kuhn@ires.in2p3.fr +//------------------------------------------------------------------------- + +#include +#include "AliPDG.h" + +class AliITStrackV2; +class AliV0vertex; + +#define kXiMinus 3312 +#define kXiPlusBar -3312 +#define kOmegaMinus 3334 +#define kOmegaPlusBar -3334 + +class AliCascadeVertex : public TObject { +public: + AliCascadeVertex(); + AliCascadeVertex(const AliV0vertex &vtx, const AliITStrackV2 &trk); + + void ChangeMassHypothesis(Int_t code=kXiMinus); + + Int_t GetPdgCode() const {return fPdgCode;} + Double_t GetEffMass() const {return fEffMass;} + Double_t GetChi2() const {return fChi2;} + void GetPxPyPz(Double_t &px, Double_t &py, Double_t &pz) const; + void GetXYZ(Double_t &x, Double_t &y, Double_t &z) const; + Double_t GetD(Double_t x0=0.,Double_t y0=0.,Double_t z0=0.) const; + void GetV0labels(Int_t &n, Int_t &p) const {n=fV0lab[0]; p=fV0lab[1];} + Int_t GetBachelorLabel() const {return fBachLab;} + +private: + Int_t fPdgCode; // reconstructed cascade type (PDG code) + Double_t fEffMass; // reconstructed cascade effective mass + Double_t fChi2; // chi2 value + Double_t fPos[3]; // cascade vertex position (global) + Double_t fPosCov[6]; // covariance matrix of the vertex position + + Int_t fV0lab[2]; // labels of the V0 daughter tracks + Double_t fV0mom[3]; // V0 momentum (global) + Double_t fV0momCov[6]; // covariance matrix of the V0 momentum. + + Int_t fBachLab; // label of the bachelor track + Double_t fBachMom[3]; // bachelor momentum (global) + Double_t fBachMomCov[6]; // covariance matrix of the bachelor momentum. + + ClassDef(AliCascadeVertex,1) // reconstructed cascade vertex +}; + +#endif + + diff --git a/ITS/AliCascadeVertexer.cxx b/ITS/AliCascadeVertexer.cxx new file mode 100644 index 00000000000..c692598d831 --- /dev/null +++ b/ITS/AliCascadeVertexer.cxx @@ -0,0 +1,251 @@ +/************************************************************************** + * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * + * * + * Author: The ALICE Off-line Project. * + * Contributors are mentioned in the code where appropriate. * + * * + * Permission to use, copy, modify and distribute this software and its * + * documentation strictly for non-commercial purposes is hereby granted * + * without fee, provided that the above copyright notice appears in all * + * copies and that both the copyright notice and this permission notice * + * appear in the supporting documentation. The authors make no claims * + * about the suitability of this software for any purpose. It is * + * provided "as is" without express or implied warranty. * + **************************************************************************/ + +//------------------------------------------------------------------------- +// Implementation of the cascade vertexer class +// +// Origin: Christian Kuhn, IReS, Strasbourg, christian.kuhn@ires.in2p3.fr +//------------------------------------------------------------------------- +#include +#include +#include +#include + +#include "AliCascadeVertex.h" +#include "AliCascadeVertexer.h" +#include "AliV0vertex.h" +#include "AliITStrackV2.h" + +ClassImp(AliCascadeVertexer) + +Int_t +AliCascadeVertexer::V0sTracks2CascadeVertices(const TFile *inp, TFile *out) { + + //-------------------------------------------------------------------- + //This function reconstructs cascade vertices + //-------------------------------------------------------------------- + TFile *in=(TFile*)inp; + TDirectory *savedir=gDirectory; + + // Tree for vertices(V0's) + + TTree *vtxTree=(TTree*)gDirectory->Get("TreeV"); + TBranch *branch=vtxTree->GetBranch("vertices"); + Int_t nentrV0=(Int_t)vtxTree->GetEntries(); + + TObjArray vtxV0(nentrV0); + + // fill TObjArray vtxV0 with vertices + + Int_t i, nV0=0; + for (i=0; iSetAddress(&ioVertex); + vtxTree->GetEvent(i); + nV0++; + vtxV0.AddLast(ioVertex); + + } + + + in->cd(); + + // Tree for tracks + + TTree *trkTree=(TTree*)in->Get("TreeT_ITS_0"); + branch=trkTree->GetBranch("tracks"); + Int_t nentr=(Int_t)trkTree->GetEntries(); + + TObjArray trks(nentr); + + // fill TObjArray trks with tracks + + Int_t ntrack=0; + + for (i=0; iSetAddress(&iotrack); + trkTree->GetEvent(i); + ntrack++; trks.AddLast(iotrack); + + } + + // create Tree for cascades + + out->cd(); + + TTree cascTree("TreeCasc","Tree with cascades"); + AliCascadeVertex *ioCascade=0; + cascTree.Branch("cascades","AliCascadeVertex",&ioCascade,32000,0); + + // loop on all vertices + + Double_t massLambda=1.11568; + Int_t ncasc=0; + + for (i=0; iChangeMassHypothesis(kLambda0); //I.B. + + if (V0ver->GetEffMass()GetEffMass()>massLambda+fMassWin) continue; + + if (V0ver->GetD(0,0,0)GetD())GetPdgCode()==kLambda0 && bachtrk->Get1Pt()<0.) continue; // condition on V0 label + if (V0ver->GetPdgCode()==kLambda0Bar && bachtrk->Get1Pt()>0.) continue; // + good sign for bachelor + + AliV0vertex V0(*V0ver), *pV0=&V0; + AliITStrackV2 bt(*bachtrk), *pbt=&bt; + + // calculation of the distance of closest approach between the V0 and the bachelor + + Double_t dca=PropagateToDCA(pV0,pbt); + if (dca > fDCAmax) continue; // cut on dca + + // construction of a cascade object + + AliCascadeVertex cascade(*pV0,*pbt); + if (cascade.GetChi2() > fChi2max) continue; + + // get cascade decay position (V0, bachelor) + + Double_t x,y,z; cascade.GetXYZ(x,y,z); + Double_t r2=x*x + y*y; + if (r2 > fRmax*fRmax) continue; // condition on fiducial zone + if (r2 < fRmin*fRmin) continue; + + // get cascade momentum + + Double_t px,py,pz; cascade.GetPxPyPz(px,py,pz); + Double_t p2=px*px+py*py+pz*pz; + Double_t cost=(x*px+y*py+z*pz)/TMath::Sqrt(p2*(r2+z*z)); + + if (costcd(); + + return 0; +} + +inline Double_t det(Double_t a00, Double_t a01, Double_t a10, Double_t a11){ + // determinant 2x2 + return a00*a11 - a01*a10; +} + +inline Double_t det (Double_t a00,Double_t a01,Double_t a02, + Double_t a10,Double_t a11,Double_t a12, + Double_t a20,Double_t a21,Double_t a22) { + // determinant 3x3 + return + a00*det(a11,a12,a21,a22)-a01*det(a10,a12,a20,a22)+a02*det(a10,a11,a20,a21); +} + + + + +Double_t +AliCascadeVertexer::PropagateToDCA(AliV0vertex *v, AliITStrackV2 *t) { + //-------------------------------------------------------------------- + // This function returns the DCA between the V0 and the track + //-------------------------------------------------------------------- + + Double_t phi, x, par[5]; + Double_t alpha, cs1, sn1; + + t->GetExternalParameters(x,par); alpha=t->GetAlpha(); + phi=TMath::ASin(par[2]) + alpha; + Double_t px1=TMath::Cos(phi), py1=TMath::Sin(phi), pz1=par[3]; + + + cs1=TMath::Cos(alpha); sn1=TMath::Sin(alpha); + Double_t x1=x*cs1 - par[0]*sn1; + Double_t y1=x*sn1 + par[0]*cs1; + Double_t z1=par[1]; + + + Double_t x2,y2,z2; // position and momentum of V0 + Double_t px2,py2,pz2; + + v->GetXYZ(x2,y2,z2); + v->GetPxPyPz(px2,py2,pz2); + +// calculation dca + + Double_t dd= det(x2-x1,y2-y1,z2-z1,px1,py1,pz1,px2,py2,pz2); + Double_t ax= det(py1,pz1,py2,pz2); + Double_t ay=-det(px1,pz1,px2,pz2); + Double_t az= det(px1,py1,px2,py2); + + Double_t dca=TMath::Abs(dd)/TMath::Sqrt(ax*ax + ay*ay + az*az); + +//points of the DCA + Double_t t1 = det(x2-x1,y2-y1,z2-z1,px2,py2,pz2,ax,ay,az)/ + det(px1,py1,pz1,px2,py2,pz2,ax,ay,az); + + x1 += px1*t1; y1 += py1*t1; //z1 += pz1*t1; + + + //propagate track to the points of DCA + + x1=x1*cs1 + y1*sn1; + if (!t->PropagateTo(x1,0.,0.)) { + cerr<<"AliV0vertexer::PropagateToDCA: propagation failed !\n"; + return 1.e+33; + } + + return dca; +} + + + + + + + + + + + diff --git a/ITS/AliCascadeVertexer.h b/ITS/AliCascadeVertexer.h new file mode 100644 index 00000000000..3dae6f3f884 --- /dev/null +++ b/ITS/AliCascadeVertexer.h @@ -0,0 +1,72 @@ +#ifndef ALICASCADEVERTEXER_H +#define ALICASCADEVERTEXER_H +/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * + * See cxx source for full Copyright notice */ + +//------------------------------------------------------------------ +// Cascade Vertexer Class +// +// Origin: Christian Kuhn, IReS, Strasbourg, christian.kuhn@ires.in2p3.fr +//------------------------------------------------------------------ + +#include "TObject.h" + +class TFile; +class AliITStrackV2; +class AliV0vertex; + +//_____________________________________________________________________________ +class AliCascadeVertexer : public TObject { +public: + AliCascadeVertexer(); + AliCascadeVertexer(const Double_t cuts[8]); + void SetCuts(const Double_t cuts[8]); + + Int_t V0sTracks2CascadeVertices(const TFile *in, TFile *out); + Double_t PropagateToDCA(AliV0vertex *vtx, AliITStrackV2 *trk); + + void GetCuts(Double_t cuts[8]) const; + +private: + Double_t fChi2max; // maximal allowed chi2 + Double_t fDV0min; // min. allowed V0 impact parameter + Double_t fMassWin; // window around the Lambda mass + Double_t fDBachMin; // min. allowed bachelor impact parameter + Double_t fDCAmax; // maximal allowed DCA between the V0 and the track + Double_t fCPAmax; // maximal allowed cosine of the cascade pointing angle + Double_t fRmin, fRmax;// max & min radii of the fiducial volume + + ClassDef(AliCascadeVertexer,1) // cascade verterxer +}; + +inline AliCascadeVertexer::AliCascadeVertexer() { + fChi2max=33.; + fDV0min=0.015; fMassWin=0.05; fDBachMin=0.015; + fDCAmax=0.01; fCPAmax=0.025; + fRmin=0.5; fRmax=2.5; +} + +inline AliCascadeVertexer::AliCascadeVertexer(const Double_t cuts[8]) { + fChi2max=cuts[0]; + fDV0min=cuts[1]; fMassWin=cuts[2]; fDBachMin=cuts[3]; + fDCAmax=cuts[4]; fCPAmax=cuts[5]; + fRmin=cuts[6]; fRmax=cuts[7]; +} + +inline void AliCascadeVertexer::SetCuts(const Double_t cuts[8]) { + fChi2max=cuts[0]; + fDV0min=cuts[1]; fMassWin=cuts[2]; fDBachMin=cuts[3]; + fDCAmax=cuts[4]; fCPAmax=cuts[5]; + fRmin=cuts[6]; fRmax=cuts[7]; +} + +inline void AliCascadeVertexer::GetCuts(Double_t cuts[8]) const { + cuts[0]=fChi2max; + cuts[1]=fDV0min; cuts[2]=fMassWin; cuts[3]=fDBachMin; + cuts[4]=fDCAmax; cuts[5]=fCPAmax; + cuts[6]=fRmin; cuts[7]=fRmax; +} + +#endif + + diff --git a/ITS/AliITSComparisonV2.C b/ITS/AliITSComparisonV2.C index 2f82cc2489e..04fe3422c9b 100644 --- a/ITS/AliITSComparisonV2.C +++ b/ITS/AliITSComparisonV2.C @@ -133,7 +133,7 @@ Int_t AliITSComparisonV2(Int_t event=0) { cerr<<"Track "<Propagate(track->GetAlpha(),3.,0.1/65.19*1.848,0.1*1.848); + track->PropagateTo(3.,0.0028,65.19); track->PropagateToVertex(); if (lab==tlab) hfound->Fill(ptg); diff --git a/ITS/AliITSFindClustersV2.C b/ITS/AliITSFindClustersV2.C index d915e30d956..28424495d84 100644 --- a/ITS/AliITSFindClustersV2.C +++ b/ITS/AliITSFindClustersV2.C @@ -6,108 +6,48 @@ #include "AliITSgeom.h" #include "AliITSRecPoint.h" #include "AliITSclusterV2.h" - #include "AliITSsimulationFastPoints.h" #include "TFile.h" #include "TTree.h" #include "TParticle.h" #endif -//Int_t AliITSFindClustersV2() { -Int_t AliITSFindClustersV2(char SlowOrFast='s') { - +Int_t AliITSFindClustersV2() { /**************************************************************** - * Just something to start with * + * This macro converts AliITSRecPoint(s) to AliITSclusterV2(s) * ****************************************************************/ - cerr<<"Looking for clusters...\n"; - cout<<"!!!! SlowOrFast = "< AliITSclusterV2 ////////////// - - // for the fast recpoints - if(SlowOrFast=='f') { - cout<<"22 !!!! SlowOrFast = "<IsOpen()) { - cerr<<"Can't open galice.root !\n"; - return 1; - } - - if (!(gAlice=(AliRun*)in->Get("gAlice"))) { - cerr<<"Can't find gAlice !\n"; - return 2; - } - - Int_t ev=0; - gAlice->GetEvent(ev); - - AliITS *ITS = (AliITS*)gAlice->GetModule("ITS"); - if (!ITS) { - cerr<<"Can't find the ITS !\n"; - return 3; - } - - gAlice->MakeTree("R"); ITS->MakeBranch("R",0); -//////////////// Taken from ITSHitsToFastPoints.C /////////////////////// - ITS->SetSimulationModel(0,new AliITSsimulationFastPoints()); - ITS->SetSimulationModel(1,new AliITSsimulationFastPoints()); - ITS->SetSimulationModel(2,new AliITSsimulationFastPoints()); - Int_t nsignal=25; - Int_t size=-1; - Int_t bgr_ev=Int_t(ev/nsignal); - ITS->HitsToFastRecPoints(ev,bgr_ev,size," ","All"," "); -////////////////////////////////////////////////////////////////////////// - - delete gAlice; gAlice=0; - in->Close(); - */ - } // end of fast recpoints - - - // for the slow points - - /*TFile */in=TFile::Open("galice.root"); + cerr<<"AliITSRecPoint(s) -> AliITSclusterV2(s)...\n"; if (gAlice) {delete gAlice; gAlice=0;} + TFile *in=TFile::Open("galice.root"); + if (!in->IsOpen()) { cerr<<"Can't open galice.root !\n"; return 1; } + if (!(gAlice=(AliRun*)in->Get("gAlice"))) { cerr<<"Can't find gAlice !\n"; - return 4; + return 2; } - gAlice->GetEvent(0); - /*AliITS */ITS = (AliITS*)gAlice->GetModule("ITS"); - if (!ITS) { - cerr<<"Can't find the ITS !\n"; - return 5; - } + AliITS *ITS = (AliITS*)gAlice->GetModule("ITS"); + if (!ITS) { cerr<<"Can't find the ITS !\n"; return 3; } AliITSgeom *geom=ITS->GetITSgeom(); TFile *out=TFile::Open("AliITSclustersV2.root","new"); if (!out->IsOpen()) { cerr<<"Delete old AliITSclustersV2.root !\n"; - return 6; + return 4; } geom->Write(); TClonesArray *clusters=new TClonesArray("AliITSclusterV2",10000); - //TTree *cTree=new TTree("cTree","ITS clusters"); TTree *cTree=new TTree("TreeC_ITS_0","ITS clusters"); cTree->Branch("Clusters",&clusters); TTree *pTree=gAlice->TreeR(); - if (!pTree) { - cerr<<"Can't get TreeR !\n"; - return 7; - } + if (!pTree) { cerr<<"Can't get TreeR !\n"; return 5; } TBranch *branch=pTree->GetBranch("ITSRecPoints"); - if (!branch) { - cerr<<"Can't get ITSRecPoints branch !\n"; - return 8; - } + if (!branch) { cerr<<"Can't get ITSRecPoints branch !\n"; return 6; } TClonesArray *points=new TClonesArray("AliITSRecPoint",10000); branch->SetAddress(&points); @@ -117,7 +57,7 @@ Int_t AliITSFindClustersV2(char SlowOrFast='s') { cerr<<"Number of entries: "<mip normalization factor for the SDD and SSD + Float_t lp[5]; Int_t lab[6]; //Why can't it be inside a loop ? for (Int_t i=0; iClear(); @@ -130,39 +70,36 @@ Int_t AliITSFindClustersV2(char SlowOrFast='s') { Int_t ndet=(lad-1)*geom->GetNdetectors(lay) + (det-1); nclusters+=ncl; - kmip=1.; - if(lay<5&&lay>2){kmip=280.;}; // b.b. - if(lay<7&&lay>4){kmip=38.;}; - //cout<<"i,ncl ="<mip normalization factor for the SDD and SSD + if(lay==4 || lay==3){kmip=280.;}; + if(lay==6 || lay==5){kmip=38.;}; for (Int_t j=0; jUncheckedAt(j); - Float_t lp[5]; + //Float_t lp[5]; lp[0]=-p->GetX()-yshift; if (lay==1) lp[0]=-lp[0]; lp[1]=p->GetZ()+zshift; lp[2]=p->GetSigmaX2(); lp[3]=p->GetSigmaZ2(); lp[4]=p->GetQ(); lp[4]/=kmip; - Int_t lab[6]; + //Int_t lab[6]; lab[0]=p->GetLabel(0);lab[1]=p->GetLabel(1);lab[2]=p->GetLabel(2); lab[3]=ndet; - Int_t label=lab[0]; - //if(label<=0) cout<<" !!!!!!!lab="<=0) { // b.b. - TParticle *part=(TParticle*)gAlice->Particle(label); - label=-3; - while (part->P() < 0.005) { - Int_t m=part->GetFirstMother(); - if (m<0) {cerr<<"Primary momentum: "<P()<Particle(label); - } - } - if (lab[1]<0) lab[1]=label; - else if (lab[2]<0) lab[2]=label; - else cerr<<"No empty labels !\n"; + Int_t label=lab[0]; + if (label>=0) { + TParticle *part=(TParticle*)gAlice->Particle(label); + label=-3; + while (part->P() < 0.005) { + Int_t m=part->GetFirstMother(); + if (m<0) {cerr<<"Primary momentum: "<P()<Particle(label); + } + if (lab[1]<0) lab[1]=label; + else if (lab[2]<0) lab[2]=label; + else cerr<<"No empty labels !\n"; + } new(cl[j]) AliITSclusterV2(lab,lp); } diff --git a/ITS/AliITSrecoV2.h b/ITS/AliITSrecoV2.h index ef117753a9e..16b91ce5b57 100644 --- a/ITS/AliITSrecoV2.h +++ b/ITS/AliITSrecoV2.h @@ -24,8 +24,8 @@ 1.44e-4, 1.44e-4, 7.84e-6, 7.84e-6, 0.006889, 0.006889 }; - const Double_t kChi2PerCluster=7.;//5.; - const Double_t kMaxChi2=15.;//12; + const Double_t kChi2PerCluster=5.;//7.; + const Double_t kMaxChi2=17.;//15.; const Double_t kMaxRoad=3.0; const Double_t kSigmaYV=0.005e+0; diff --git a/ITS/AliITStestV2.C b/ITS/AliITStestV2.C index 91160ce7dc7..afe92577aec 100644 --- a/ITS/AliITStestV2.C +++ b/ITS/AliITStestV2.C @@ -1,4 +1,4 @@ -Int_t AliITStestV2() { +Int_t AliITStestV2(Char_t SlowOrFast='s') { Int_t rc=0; if (gAlice) {delete gAlice; gAlice=0;} @@ -15,6 +15,18 @@ Int_t AliITStestV2() { delete gAlice; gAlice=0; in->Close(); + if (SlowOrFast=='f') { + cerr<<"Fast AliITSRecPoint(s) !\n"; + gROOT->LoadMacro("$(ALICE_ROOT)/ITS/ITSHitsToFastPoints.C"); + ITSHitsToFastPoints(); + } else { + cerr<<"Slow AliITSRecPoint(s) !\n"; + gROOT->LoadMacro("$(ALICE_ROOT)/ITS/AliITSHits2Digits.C"); + AliITSHits2Digits(); + gROOT->LoadMacro("$(ALICE_ROOT)/ITS/AliITSFindClusters.C"); + AliITSFindClusters(); + } + gROOT->LoadMacro("$(ALICE_ROOT)/ITS/AliITSFindClustersV2.C"); if (rc=AliITSFindClustersV2()) return rc; diff --git a/ITS/AliITStrackV2.cxx b/ITS/AliITStrackV2.cxx index 4c3609088ec..80f9a16a137 100644 --- a/ITS/AliITStrackV2.cxx +++ b/ITS/AliITStrackV2.cxx @@ -17,6 +17,7 @@ // Implementation of the ITS track class // // Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch +// dEdx analysis by: Boris Batyunya, JINR, Boris.Batiounia@cern.ch //------------------------------------------------------------------------- #include @@ -88,43 +89,26 @@ AliITStrackV2::AliITStrackV2(const AliITStrackV2& t) : AliKalmanTrack(t) { fC40=t.fC40; fC41=t.fC41; fC42=t.fC42; fC43=t.fC43; fC44=t.fC44; Int_t n=GetNumberOfClusters(); - //for (Int_t i=0; iGet1Pt()); - Double_t c =TMath::Abs(Get1Pt()); + //Double_t co=TMath::Abs(t->Get1Pt()); + //Double_t c =TMath::Abs(Get1Pt()); + Double_t co=t->GetSigmaY2()*t->GetSigmaZ2(); + Double_t c =GetSigmaY2()*GetSigmaZ2(); if (c>co) return 1; else if (cGet1Pt()*t->Get1Pt()); - Double_t bo2=po2/(po2 + t->GetMass()*t->GetMass()); - if (p2*b2>po2*bo2) return -1; - else if (p2*b2= 0.99999) { + if (TMath::Abs(f2) >= 0.9999) { Int_t n=GetNumberOfClusters(); if (n>kWARN) cerr<= 0.99999) { + if (TMath::Abs(f2) >= 0.9999) { Int_t n=GetNumberOfClusters(); if (n>kWARN) cerr<11.5) - //if (fP1*fP4<0) { - // if (n>kWARN) cerr<<"fP1*fP4="<=1) {if (n>kWARN) cerr<<"fP2="<kWARN) cerr<<"fC00="<kWARN) cerr<<"fC11="<kWARN) cerr<<"fC22="<kWARN) cerr<<"fC33="<kWARN) cerr<<"fC44="<kWARN) { cerr<<" bad determinant "<=0.9999){ + if (n>kWARN) cout<<"AliITStrackV2::Invariant : fP2="<9.) { + if (n>kWARN) cout<<"AliITStrackV2::Invariant : fC00="<9.) { + if (n>kWARN) cout<<"AliITStrackV2::Invariant : fC11="<1.) { + if (n>kWARN) cout<<"AliITStrackV2::Invariant : fC22="<1.) { + if (n>kWARN) cout<<"AliITStrackV2::Invariant : fC33="<6e-5) { + if (n>kWARN) cout<<"AliITStrackV2::Invariant : fC44="<= 0.99999) { + if (TMath::Abs(pp2) >= 0.9999) { Int_t n=GetNumberOfClusters(); if (n>kWARN) cerr<= 0.99999) { + if (TMath::Abs(f2) >= 0.9999) { Int_t n=GetNumberOfClusters(); if (n>kWARN) cerr< #include @@ -32,7 +33,7 @@ //#define DEBUG #ifdef DEBUG -Int_t LAB=7; +Int_t LAB=5126; #endif AliITStrackerV2::AliITSlayer AliITStrackerV2::fLayers[kMaxLayer]; // ITS layers @@ -177,7 +178,18 @@ Int_t AliITStrackerV2::Clusters2Tracks(const TFile *inp, TFile *out) { nentr=(Int_t)tpcTree->GetEntries(); for (Int_t i=0; iGetEvent(i); - itsTracks.AddLast(new AliITStrackV2(*itrack)); + AliITStrackV2 *t=new AliITStrackV2(*itrack); + if (TMath::Abs(t->GetD())>4) continue; + + t->PropagateTo(80.,0.0053); + if (TMath::Abs(t->GetY())>13.) t->CorrectForMaterial(0.03); + if (TMath::Abs(t->GetZ())<0.2) t->CorrectForMaterial(0.40); + t->PropagateTo(61.,0.0052); + Double_t xk=52.,x,y,z; t->GetGlobalXYZat(xk,x,y,z); + if (TMath::Abs(y)<7.77) t->PropagateTo(xk,0.19,24.); + t->PropagateTo(50.,0.001); + + itsTracks.AddLast(t); } delete tpcTree; //Thanks to Mariana Bondila delete itrack; @@ -198,11 +210,13 @@ Int_t AliITStrackerV2::Clusters2Tracks(const TFile *inp, TFile *out) { AliITStrackV2 *t=(AliITStrackV2*)itsTracks.UncheckedAt(i); if (t==0) continue; //this track has been already tracked Int_t tpcLabel=t->GetLabel(); //save the TPC track label + #ifdef DEBUG lbl=tpcLabel; if (TMath::Abs(tpcLabel)!=LAB) continue; cout<GetDetectorIndex()) { - idet=c->GetDetectorIndex(); + idet=c->GetDetectorIndex(); + if (idet != fTrackToFollow.GetDetectorIndex()) { const AliITSdetector &det=layer.GetDetector(idet); r=det.GetR(); phi=det.GetPhi(); - if (!fTrackToFollow.Propagate(phi,r,0.,0.)) throw ""; + if (!fTrackToFollow.Propagate(phi,r)) throw ""; + fTrackToFollow.SetDetectorIndex(idet); } Double_t chi2=fTrackToFollow.GetPredictedChi2(c); if (chi2GetDetectorIndex()) continue; + idet=cc->GetDetectorIndex(); + if (idet != fTrackToFollow.GetDetectorIndex()) continue; Double_t chi2=fTrackToFollow.GetPredictedChi2(cc); if (chi2 r+1) break; fTrackToFollow.SetDetectorIndex(idet); //Select possible prolongations and store the current track estimation track.~AliITStrackV2(); new(&track) AliITStrackV2(fTrackToFollow); - Double_t dz=3*TMath::Sqrt(track.GetSigmaZ2() + kSigmaZ2[fI]); + Double_t dz=3*TMath::Sqrt(track.GetSigmaZ2() + kSigmaZ2[i]); if (dz < 0.5*TMath::Abs(track.GetTgl())) dz=0.5*TMath::Abs(track.GetTgl()); if (dz > kMaxRoad) { //cerr<<"AliITStrackerV2::FollowProlongation: too broad road in Z !\n"; break; } - Double_t dy=4*TMath::Sqrt(track.GetSigmaY2() + kSigmaY2[fI]); + + if (TMath::Abs(fTrackToFollow.GetZ()-GetZ()) > r+dz) break; + + Double_t dy=4*TMath::Sqrt(track.GetSigmaY2() + kSigmaY2[i]); if (dy < 0.5*TMath::Abs(track.GetSnp())) dy=0.5*TMath::Abs(track.GetSnp()); if (dy > kMaxRoad) { //cerr<<"AliITStrackerV2::FollowProlongation: too broad road in Y !\n"; break; } + Double_t zmin=track.GetZ() - dz; Double_t zmax=track.GetZ() + dz; Double_t ymin=track.GetY() + r*phi - dy; Double_t ymax=track.GetY() + r*phi + dy; - layer.SelectClusters(zmin,zmax,ymin,ymax); + layer.SelectClusters(zmin,zmax,ymin,ymax); + fI--; //take another prolongation if (!TakeNextProlongation()) if (!tryAgain--) break; @@ -562,20 +525,16 @@ cout<GetLabel(0)!=TMath::Abs(lbl)) continue; -#endif - + //if (c->GetLabel(0)!=TMath::Abs(lbl)) continue; Int_t idet=c->GetDetectorIndex(); if (t.GetDetectorIndex()!=idet) { const AliITSdetector &det=layer.GetDetector(idet); - if (!t.Propagate(det.GetPhi(),det.GetR(),0.,0.)) { + if (!t.Propagate(det.GetPhi(),det.GetR())) { //cerr<<"AliITStrackerV2::TakeNextProlongation: " //"propagation failed !\n"; continue; } t.SetDetectorIndex(idet); + if (TMath::Abs(t.GetZ()-GetZ()) > layer.GetR()+dz) continue; #ifdef DEBUG cout<GetZ()) > dz) continue; if (TMath::Abs(t.GetY() - c->GetY()) > dy) continue; - //m[0]=fYV; m[1]=fZV; - //chi2=t.GetPredictedChi2(c,m,d/21.82*2.33); - chi2=t.GetPredictedChi2(c); - - if (chi2GetQ(), - fTrackToFollow.GetNumberOfClusters()-1); - if (constraint) fTrackToFollow.Improve(d/21.82*2.33,GetY(),GetZ()); + fTrackToFollow. + SetSampledEdx(c->GetQ(),fTrackToFollow.GetNumberOfClusters()-1); //b.b. + + { + Double_t d=layer.GetThickness(fTrackToFollow.GetY(),fTrackToFollow.GetZ()); + fTrackToFollow.CorrectForMaterial(d); + } + if (fConstraint[fPass]) { + Double_t d=GetEffectiveThickness(0,0); //Think of this !!!! + fTrackToFollow.Improve(d,GetY(),GetZ()); + } #ifdef DEBUG cout<<"accepted lab="<GetLabel(0)<<' ' <3.40) d+=0.0036; + if (TMath::Abs(y-2.50)<0.10) d+=(0.02-0.0036); + if (TMath::Abs(y+1.90)<0.10) d+=(0.02-0.0036); + for (Int_t i=0; i<12; i++) { + if (TMath::Abs(z-3.6*(i+0.5))<0.20) {d+=0.0036; break;} + if (TMath::Abs(z+3.6*(i+0.5))<0.20) {d+=0.0036; break;} + if (TMath::Abs(z-3.6*(i+0.929))<0.50) {d+=(0.02-0.0036); break;} + if (TMath::Abs(z+3.6*(i+0.104))<0.50) {d+=(0.02-0.0036); break;} + } + } else + if (373.40) d+=0.0036; + if (TMath::Abs(y-2.20)<0.10) d+=(0.02-0.0036); + if (TMath::Abs(y+2.20)<0.10) d+=(0.02-0.0036); + for (Int_t i=0; i<11; i++) { + if (TMath::Abs(z-3.6*i)<0.20) {d+=0.0036; break;} + if (TMath::Abs(z+3.6*i)<0.20) {d+=0.0036; break;} + if (TMath::Abs(z-3.6*(i+0.54))<0.50) {d+=(0.02-0.0036); break;} + if (TMath::Abs(z+3.6*(i+0.58))<0.50) {d+=(0.02-0.0036); break;} + } + } else + if (133.30) d+=0.0034; + if (TMath::Abs(y-2.10)<0.20) d+=0.0034*3; + if (TMath::Abs(y+2.10)<0.20) d+=0.0034*3; + for (Int_t i=0; i<4; i++) { + if (TMath::Abs(z-7.3*i)<0.60) {d+=0.0034; break;} + if (TMath::Abs(z+7.3*i)<0.60) {d+=0.0034; break;} + } + } else + if (60.45) d+=0.0064; + if (TMath::Abs(y-3.03)<0.10) d+=0.0192; + } else + if (30.55) d+=0.0064; + if (TMath::Abs(y+0.10)<0.10) d+=0.0192; + } + + d+=0.002; + return d; +} -Double_t AliITStrackerV2::GetEffectiveThickness(Double_t phi,Double_t z) const +Double_t AliITStrackerV2::GetEffectiveThickness(Double_t y,Double_t z) const { //-------------------------------------------------------------------- - //Returns the thickness between the current layer and the vertex + //Returns the thickness between the current layer and the vertex (units X0) //-------------------------------------------------------------------- - //-pi1) { Double_t xi=0.5*(fLayers[1].GetR()+fLayers[2].GetR()); - d+=0.034*xi*xi; + d+=0.011*xi*xi; } if (fI>3) { Double_t xi=0.5*(fLayers[3].GetR()+fLayers[4].GetR()); - d+=0.039*xi*xi; + d+=0.0053*xi*xi; } return d/(xn*xn); } @@ -817,7 +817,7 @@ Int_t AliITStrackerV2::AliITSlayer::InRoad() const { for (Int_t i=fI; iGetZ() > fZmax) break; - //if (c->IsUsed()) continue; + if (c->IsUsed()) continue; const AliITSdetector &det=GetDetector(c->GetDetectorIndex()); Double_t y=fR*det.GetPhi() + c->GetY(); diff --git a/ITS/AliITStrackerV2.h b/ITS/AliITStrackerV2.h index dc7f52eb6cc..b6ff0392ae1 100644 --- a/ITS/AliITStrackerV2.h +++ b/ITS/AliITStrackerV2.h @@ -53,7 +53,7 @@ public: AliITSclusterV2 *GetCluster(Int_t i) const {return fClusters[i];} AliITSdetector &GetDetector(Int_t n) const { return fDetectors[n]; } Int_t FindDetectorIndex(Double_t phi, Double_t z) const; - Double_t GetThickness(Double_t phi, Double_t z) const; + Double_t GetThickness(Double_t y, Double_t z) const; Int_t InRoad() const ; Int_t GetNumberOfClusters() const {return fN;} private: @@ -74,7 +74,7 @@ public: private: void CookLabel(AliKalmanTrack *t,Float_t wrong) const; - Double_t GetEffectiveThickness(Double_t phi, Double_t z) const; + Double_t GetEffectiveThickness(Double_t y, Double_t z) const; void FollowProlongation(); Int_t TakeNextProlongation(); void ResetBestTrack() { diff --git a/ITS/AliV0Comparison.C b/ITS/AliV0Comparison.C index 53bc2f31216..92bdc39b220 100644 --- a/ITS/AliV0Comparison.C +++ b/ITS/AliV0Comparison.C @@ -1,3 +1,12 @@ +/**************************************************************************** + * Very important, delicate and rather obscure macro. * + * * + * Creates list of "findable" V0s, * + * calculates efficiency, resolutions etc. * + * * + * Origin: I.Belikov, IReS, Strasbourg, Jouri.Belikov@cern.ch * + ****************************************************************************/ + #ifndef __CINT__ #include #include @@ -113,7 +122,7 @@ Int_t AliV0Comparison(Int_t code=310) { //Lambda=3122, LambdaBar=-3122 Double_t mmin=V0mass-V0window, mmax=V0mass+V0window; TH1F *v0s =new TH1F("v0s","V0s Effective Mass",40, mmin, mmax); - v0s->SetXTitle("(GeV/c)"); v0s->SetFillColor(6); + v0s->SetXTitle("(GeV)"); v0s->SetFillColor(6); Double_t pxg=0.,pyg=0.,ptg=0.; @@ -320,7 +329,11 @@ Int_t good_vertices(GoodVertex *gv, Int_t max) { Double_t x=p->Vx(), y=p->Vy(), z=p->Vz(), r2=x*x+y*y; if (r2r2max) continue; - + + if (gAlice->Particle(plab)->GetPDG()->Charge() < 0.) { + i=plab; plab=nlab; nlab=i; + } + gv[nv].code=code; gv[nv].plab=plab; gv[nv].nlab=nlab; gv[nv].px=p0->Px(); gv[nv].py=p0->Py(); gv[nv].pz=p0->Pz(); diff --git a/ITS/AliV0FindVertices.C b/ITS/AliV0FindVertices.C index fa10eb13773..6d1cf5c9f1e 100644 --- a/ITS/AliV0FindVertices.C +++ b/ITS/AliV0FindVertices.C @@ -14,11 +14,11 @@ Int_t AliV0FindVertices() { TFile *in=TFile::Open("AliITStracksV2.root"); if (!in->IsOpen()) {cerr<<"Can't open AliITStracksV2.root !\n"; return 2;} - Double_t cuts[]={33., // max. allowed chi2 - 0.015,// min. allowed negative daughter's impact parameter - 0.015,// min. allowed positive daughter's impact parameter - 0.060,// max. allowed DCA between the daughter tracks - 0.997,// max. allowed cosine of V0's pointing angle + Double_t cuts[]={33, // max. allowed chi2 + 0.02,// min. allowed negative daughter's impact parameter + 0.02,// min. allowed positive daughter's impact parameter + 0.080,// max. allowed DCA between the daughter tracks + 0.9985,// max. allowed cosine of V0's pointing angle 0.9, // min. radius of the fiducial volume 2.9 // max. radius of the fiducial volume }; diff --git a/ITS/AliV0vertexer.cxx b/ITS/AliV0vertexer.cxx index 077fb229047..435d25ed379 100644 --- a/ITS/AliV0vertexer.cxx +++ b/ITS/AliV0vertexer.cxx @@ -78,6 +78,7 @@ Int_t AliV0vertexer::Tracks2V0vertices(const TFile *inp, TFile *out) { for (i=0; iGetSigmaY2() + p->GetSigmaY2(); + Double_t dz2=n->GetSigmaZ2() + p->GetSigmaZ2(); + Double_t dx2=dy2; + + //dx2=dy2=dz2=1.; + Double_t p1[8]; External2Helix(n,p1); p1[6]=TMath::Sin(p1[2]); p1[7]=TMath::Cos(p1[2]); Double_t p2[8]; External2Helix(p,p2); @@ -183,19 +190,19 @@ Double_t AliV0vertexer::PropagateToDCA(AliITStrackV2 *n, AliITStrackV2 *p) { Evaluate(p2,t2,r2,g2,gg2); Double_t dx=r2[0]-r1[0], dy=r2[1]-r1[1], dz=r2[2]-r1[2]; - Double_t dm=dx*dx + dy*dy + dz*dz; + Double_t dm=dx*dx/dx2 + dy*dy/dy2 + dz*dz/dz2; Int_t max=27; while (max--) { - Double_t gt1=-(dx*g1[0] + dy*g1[1] + dz*g1[2]); - Double_t gt2=+(dx*g2[0] + dy*g2[1] + dz*g2[2]); - Double_t h11=g1[0]*g1[0] - dx*gg1[0] + - g1[1]*g1[1] - dy*gg1[1] + - g1[2]*g1[2] - dz*gg1[2]; - Double_t h22=g2[0]*g2[0] + dx*gg2[0] + - g2[1]*g2[1] + dy*gg2[1] + - g2[2]*g2[2] + dz*gg2[2]; - Double_t h12=-(g1[0]*g2[0] + g1[1]*g2[1] + g1[2]*g2[2]); + Double_t gt1=-(dx*g1[0]/dx2 + dy*g1[1]/dy2 + dz*g1[2]/dz2); + Double_t gt2=+(dx*g2[0]/dx2 + dy*g2[1]/dy2 + dz*g2[2]/dz2); + Double_t h11=(g1[0]*g1[0] - dx*gg1[0])/dx2 + + (g1[1]*g1[1] - dy*gg1[1])/dy2 + + (g1[2]*g1[2] - dz*gg1[2])/dz2; + Double_t h22=(g2[0]*g2[0] + dx*gg2[0])/dx2 + + (g2[1]*g2[1] + dy*gg2[1])/dy2 + + (g2[2]*g2[2] + dz*gg2[2])/dz2; + Double_t h12=-(g1[0]*g2[0]/dx2 + g1[1]*g2[1]/dy2 + g1[2]*g2[2]/dz2); Double_t det=h11*h22-h12*h12; @@ -215,7 +222,7 @@ Double_t AliV0vertexer::PropagateToDCA(AliITStrackV2 *n, AliITStrackV2 *p) { if (TMath::Abs(dt1)/(TMath::Abs(t1)+1.e-3) < 1.e-4) if (TMath::Abs(dt2)/(TMath::Abs(t2)+1.e-3) < 1.e-4) { - if ((gt1*gt1+gt2*gt2) > 1.e-4) + if ((gt1*gt1+gt2*gt2) > 1.e-4/dy2/dy2) cerr<<"AliV0vertexer::PropagateToDCA:" " stopped at not a stationary point !\n"; Double_t lmb=h11+h22; lmb=lmb-TMath::Sqrt(lmb*lmb-4*det); @@ -230,7 +237,7 @@ Double_t AliV0vertexer::PropagateToDCA(AliITStrackV2 *n, AliITStrackV2 *p) { Evaluate(p1,t1+dt1,r1,g1,gg1); Evaluate(p2,t2+dt2,r2,g2,gg2); dx=r2[0]-r1[0]; dy=r2[1]-r1[1]; dz=r2[2]-r1[2]; - dd=dx*dx + dy*dy + dz*dz; + dd=dx*dx/dx2 + dy*dy/dy2 + dz*dz/dz2; if (dd512) { @@ -263,7 +270,8 @@ Double_t AliV0vertexer::PropagateToDCA(AliITStrackV2 *n, AliITStrackV2 *p) { return 1.e+33; } - return TMath::Sqrt(dm); + //return TMath::Sqrt(dm); + return TMath::Sqrt(dx*dx + dy*dy + dz*dz); } diff --git a/ITS/ITSHitsToFastPoints.C b/ITS/ITSHitsToFastPoints.C index c91441f054e..11d07d051e7 100644 --- a/ITS/ITSHitsToFastPoints.C +++ b/ITS/ITSHitsToFastPoints.C @@ -40,12 +40,20 @@ void ITSHitsToFastPoints (Int_t evNumber1=0,Int_t evNumber2=0,Int_t nsignal=25, AliITS *ITS = (AliITS*) gAlice->GetModule("ITS"); if (!ITS) return; - // Set the simulation model + // Set the simulation model + + /* Bug !!! (I.Belikov) AliITSsimulationFastPoints *sim = new AliITSsimulationFastPoints(); for (Int_t i=0;i<3;i++) { ITS->SetSimulationModel(i,sim); } + */ + + for (Int_t i=0;i<3;i++) { + ITS->SetSimulationModel(i,new AliITSsimulationFastPoints()); + } + // // Event Loop @@ -71,8 +79,7 @@ void ITSHitsToFastPoints (Int_t evNumber1=0,Int_t evNumber2=0,Int_t nsignal=25, timer.Stop(); timer.Print(); } // event loop - //delete gAlice; - //gAlice=0; + delete gAlice; gAlice=0; file->Close(); } diff --git a/ITS/ITSLinkDef.h b/ITS/ITSLinkDef.h index 2ee407b44eb..4771f98460a 100644 --- a/ITS/ITSLinkDef.h +++ b/ITS/ITSLinkDef.h @@ -129,6 +129,8 @@ #pragma link C++ class AliITStrackerV2+; #pragma link C++ class AliV0vertex+; #pragma link C++ class AliV0vertexer+; +#pragma link C++ class AliCascadeVertex+; +#pragma link C++ class AliCascadeVertexer+; #pragma link C++ class AliITSVertex+; #endif diff --git a/ITS/Makefile b/ITS/Makefile index 2eeefa339e6..fd070188b66 100644 --- a/ITS/Makefile +++ b/ITS/Makefile @@ -44,7 +44,8 @@ SRCS = AliITS.cxx AliITSv1.cxx AliITSv3.cxx AliITSv5.cxx \ AliITSclusterV2.cxx AliITStrackV2.cxx AliITStrackerV2.cxx \ AliITSPid.cxx AliITStrackV2Pid.cxx \ AliITSVertex.cxx \ - AliV0vertex.cxx AliV0vertexer.cxx + AliV0vertex.cxx AliV0vertexer.cxx \ + AliCascadeVertex.cxx AliCascadeVertexer.cxx # AliITSAlignmentTrack.cxx AliITSAlignmentModule.cxx \ # Fortran sources -- 2.43.0